/* 

  Min, median, and max of Poisson distributions in Picat.

  Comparing with the order statistics (in go2/0).

  Cf my Gamble model gamble_poisson_min_max.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.

/*

  var : min val
  Probabilities:
  5: 0.2526000000000000
  6: 0.2519000000000000
  4: 0.1688000000000000
  7: 0.1554000000000000
  3: 0.0768000000000000
  8: 0.0512000000000000
  2: 0.0239000000000000
  9: 0.0115000000000000
  1: 0.0059000000000000
  10: 0.0012000000000000
  0: 0.0008000000000000
  mean = 5.3466
  HPD intervals:
  HPD interval (0.94): 3.00000000000000..8.00000000000000

  var : median val
  Probabilities:
  10: 0.3236000000000000
  9: 0.2639000000000000
  11: 0.2102000000000000
  8: 0.1022000000000000
  12: 0.0668000000000000
  7: 0.0175000000000000
  13: 0.0119000000000000
  6: 0.0021000000000000
  14: 0.0017000000000000
  5: 0.0001000000000000
  mean = 9.8566
  HPD intervals:
  HPD interval (0.94): 8.00000000000000..12.00000000000000

  var : max val
  Probabilities:
  15: 0.1995000000000000
  14: 0.1857000000000000
  16: 0.1592000000000000
  13: 0.1220000000000000
  17: 0.1170000000000000
  18: 0.0663000000000000
  12: 0.0597000000000000
  19: 0.0373000000000000
  20: 0.0192000000000000
  11: 0.0151000000000000
  21: 0.0090000000000000
  22: 0.0045000000000000
  10: 0.0028000000000000
  23: 0.0015000000000000
  24: 0.0005000000000000
  25: 0.0004000000000000
  29: 0.0001000000000000
  26: 0.0001000000000000
  8: 0.0001000000000000
  mean = 15.2619
  HPD intervals:
  HPD interval (0.94): 12.00000000000000..19.00000000000000

  var : diff
  Probabilities:
  10: 0.1592000000000000
  9: 0.1581000000000000
  8: 0.1417000000000000
  11: 0.1281000000000000
  12: 0.0975000000000000
  7: 0.0895000000000000
  13: 0.0661000000000000
  6: 0.0516000000000000
  14: 0.0430000000000000
  15: 0.0207000000000000
  5: 0.0174000000000000
  16: 0.0110000000000000
  17: 0.0060000000000000
  4: 0.0044000000000000
  18: 0.0023000000000000
  19: 0.0013000000000000
  20: 0.0008000000000000
  3: 0.0007000000000000
  21: 0.0005000000000000
  25: 0.0001000000000000
  mean = 9.9153
  HPD intervals:
  HPD interval (0.94): 5.00000000000000..14.00000000000000

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs,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() =>
  N = 11,
  Lambda = 10,
  Ps = poisson_dist_n(Lambda,N),
  MinVal = Ps.min,
  MedianVal = Ps.median,
  MaxVal = Ps.max,
  Diff = MaxVal - MinVal,
  
  add("min val",MinVal),
  add("median val",MedianVal),
  add("max val",MaxVal),
  add("diff",Diff).


/*
  Order statistics. 
  Compare with the value from the model.

  Min:
   5 0.255524089385398 
   6 0.250120433933389 
   4 0.170614069430921 
   7 0.150925210713132 
   3 0.077954165059379 
   8 0.053157637325407 
   2 0.024565317689993 
   9 0.010471941570619 
   1 0.004980409174759 
  10 0.001121211905279 
   0 0.000499285879377 
  11 0.000064235438485 
  12 0.0000019602184 
  13 0.000000031991941 
  14 0.000000000282154 
  15 0.000000000001364 
  16 0.000000000000003 
  19 0.000000000000001 
  23 0.000000000000001 
  27 0.000000000000001 
  25 0.0 
  20 0.0 
  28 0.0 
  21 0.0 
  30 -0.0 
  18 -0.0 
  29 -0.0 
  17 -0.0 
  26 -0.0 
  24 -0.000000000000001 
  22 -0.000000000000001 

  Median:
  10 0.327345923378689 
   9 0.266131927209918 
  11 0.203276490889527 
   8 0.10258101969298 
  12 0.067621626459721 
   7 0.017487206521785 
  13 0.012758881967567 
  14 0.00144237919701 
   6 0.001215345604459 
  15 0.000102847765482 
   5 0.000031090407005 
  16 0.000004845215829 
   4 0.000000254470249 
  17 0.000000157000919 
  18 0.000000003618322 
   3 0.000000000538626 
  19 0.000000000060934 
  20 0.000000000000766 
   2 0.000000000000206 
  21 0.000000000000007 
  23 0.000000000000001 
  27 0.000000000000001 
  25 0.0 
  28 0.0 
   1 0.0 
   0 0.0 
  30 -0.0 
  29 -0.0 
  26 -0.0 
  22 -0.000000000000001 
  24 -0.000000000000001 

  Max:
  15 0.193732158407237 
  14 0.181945284926558 
  16 0.162517045008593 
  13 0.125037171321823 
  17 0.11402303408936 
  18 0.070037376784896 
  12 0.057642237512694 
  19 0.038923168061375 
  20 0.020015444967688 
  11 0.016147745882951 
  21 0.009663495831037 
  22 0.004420990131765 
  10 0.002460858850023 
  23 0.001927748472905 
  24 0.000804228386795 
  25 0.000321856174384 
   9 0.000180139308496 
  26 0.000123815923081 
  27 0.000045861286614 
  28 0.000016379494341 
  29 0.000005648158049 
   8 0.000005490996768 
  30 0.000001882725795 
   7 0.000000058898339 
   6 0.000000000181248 
   5 0.000000000000124 
   4 0.0 
   3 0.0 
   2 0.0 
   1 0.0 
   0 0.0 


*/
go2 ?=>
  println("Order statistics:"),
  PDF = $poisson_dist_pdf(10),
  CDF = $poisson_dist_cdf(10),

  println("Min:"),
  Min = [V=order_statistics_with_replacement_discrete_pdf(PDF,CDF,11,1,V) : V in 0..30].sort_down(2),
  Min.printf_list,
  nl,
  println("Median:"),
  Median = [V=order_statistics_with_replacement_discrete_pdf(PDF,CDF,11,6,V) : V in 0..30].sort_down(2),
  Median.printf_list,
  nl,

  println("Max:"),
  Max = [V=order_statistics_with_replacement_discrete_pdf(PDF,CDF,11,11,V) : V in 0..30].sort_down(2),
  Max.printf_list,
  nl.
go2 => true.