/* 

  Multinomial balls in Picat.

  From Mathematica (Multinomial)
  """
  Distribute 5 balls among 3 containers, picking each container with equal probability. 
  Find the probability that no container is empty:

    Probability(x) > 0 && x2 > 0 && x3 > 0, (x1, x2, x3) ->) MultinomialDistribution(5, (1/3, 1/3, 1/3)))
    -> 
    50/81  = 0.61728395061728395062
  """

  Cf my Gamble model gamble_multinomial_balls.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 : p not empty
  Probabilities:
  true: 0.6269000000000000
  false: 0.3731000000000000

  var : balls
  Probabilities:
  [1,2,2]: 0.1257000000000000
  [2,2,1]: 0.1233000000000000
  [2,1,2]: 0.1222000000000000
  [3,1,1]: 0.0858000000000000
  [1,3,1]: 0.0853000000000000
  [1,1,3]: 0.0846000000000000
  [0,3,2]: 0.0420000000000000
  [2,3,0]: 0.0419000000000000
  [3,2,0]: 0.0401000000000000
  [2,0,3]: 0.0401000000000000
  [0,2,3]: 0.0396000000000000
  [3,0,2]: 0.0380000000000000
  [1,4,0]: 0.0226000000000000
  [4,0,1]: 0.0209000000000000
  [1,0,4]: 0.0197000000000000
  [4,1,0]: 0.0190000000000000
  [0,4,1]: 0.0187000000000000
  [0,1,4]: 0.0185000000000000
  [0,5,0]: 0.0044000000000000
  [0,0,5]: 0.0040000000000000
  [5,0,0]: 0.0036000000000000

  Theoretical:
  [2,2,1]: 0.12345679012346
  [2,1,2]: 0.12345679012346
  [1,2,2]: 0.12345679012346
  [3,1,1]: 0.08230452674897
  [1,3,1]: 0.08230452674897
  [1,1,3]: 0.08230452674897
  [3,2,0]: 0.04115226337449
  [3,0,2]: 0.04115226337449
  [2,3,0]: 0.04115226337449
  [2,0,3]: 0.04115226337449
  [0,3,2]: 0.04115226337449
  [0,2,3]: 0.04115226337449
  [4,1,0]: 0.02057613168724
  [4,0,1]: 0.02057613168724
  [1,4,0]: 0.02057613168724
  [1,0,4]: 0.02057613168724
  [0,4,1]: 0.02057613168724
  [0,1,4]: 0.02057613168724
  [5,0,0]: 0.00411522633745
  [0,5,0]: 0.00411522633745
  [0,0,5]: 0.00411522633745


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  println("Theoretical:"),
  PDFs = [Bs=multinomial_dist_pdf(5, [1/3,1/3,1/3],Bs) :
                               Bs in [ [B1,B2,B3] : B1 in 0..5,B2 in 0..5, B3 in 0..5, B1+B2+B3 == 5]],
  foreach(B=PDF in PDFs.sort_down(2))
    printf("%w: %0.14f\n",B,PDF)
  end,

  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 5,
  Ps = [1,1,1],
  Balls = multinomial_dist(N,Ps),
  [B1,B2,B3] = Balls,

  % Prob that not container is empty
  PNotEmpty = check( (B1>0,B2>0,B3>0)),
  
  add("p not empty",PNotEmpty),
  add("balls",Balls).