/* 

  Zip dist (max values) in Picat.

  Compare with
  Picat> quantile_all($zipf_dist(10, 0.3)).printf_list
  0.000001  1 
  0.00001  1 
  0.001  1 
  0.01  1 
  0.025  1 
  0.05  1 
  0.25  2 
  0.5  5 
  0.75  7 
  0.84  8 
  0.975 10 
  0.99 10 
  0.999 10 
  0.99999 10 
  0.999999 10 

  Picat> X = zipf_dist_mean(10, 0.3)                  
  X = 4.899692026871479


  Cf my Gamble model gamble_zipf_dist_max_values.rkt

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import ppl_distributions, ppl_utils.
import util.
% import ordset.

main => go.

/*

  num_samples = 1
  var : z max
  Probabilities:
  1: 0.1478000000000000
  2: 0.1303000000000000
  3: 0.1139000000000000
  4: 0.0976000000000000
  5: 0.0952000000000000
  6: 0.0934000000000000
  7: 0.0858000000000000
  8: 0.0848000000000000
  9: 0.0789000000000000
  10: 0.0723000000000000
  mean = 4.889


  num_samples = 3
  var : z max
  Probabilities:
  10: 0.2137000000000000
  9: 0.1896000000000000
  8: 0.1549000000000000
  7: 0.1300000000000000
  6: 0.1083000000000000
  5: 0.0817000000000000
  4: 0.0605000000000000
  3: 0.0375000000000000
  2: 0.0196000000000000
  1: 0.0042000000000000
  mean = 7.4488


  num_samples = 10
  var : z max
  Probabilities:
  10: 0.5504000000000000
  9: 0.2720000000000000
  8: 0.1132000000000000
  7: 0.0451000000000000
  6: 0.0157000000000000
  5: 0.0029000000000000
  4: 0.0006000000000000
  3: 0.0001000000000000
  mean = 9.2847


  num_samples = 20
  var : z max
  Probabilities:
  10: 0.7947000000000000
  9: 0.1711000000000000
  8: 0.0294000000000000
  7: 0.0045000000000000
  6: 0.0003000000000000
  mean = 9.7554


  num_samples = 40
  var : z max
  Probabilities:
  10: 0.9650000000000000
  9: 0.0342000000000000
  8: 0.0008000000000000
  mean = 9.9642

*/
go ?=>
  member(NumSamples,[1,3,10,20,40]),
  println(num_samples=NumSamples),
  reset_store,
  run_model(10_000,$model(NumSamples),[show_probs_trunc,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(NumSamples) =>
  N = 10,
  S = 0.3,
  Z = zipf_dist_n(N,S,NumSamples),
  ZMax = Z.max,
  
  add("z max",ZMax).