/* 

  Sum distribution in Picat.


  Distribution of summing discrete uniform numbers to a specific sum.

  Inspired by Fletcher Thompson "What’s the Probability Ten Dice Add Up To 12?"
  https://medium.com/puzzle-sphere/whats-the-probability-ten-dice-add-up-to-12-83f637205505

  Probability of getting the sum 12 when throwing 10 fair dice:
  Picat> P = sum_prob_dist_pdf(1,6,10,12), X = 1/P

  Picat>  P = sum_prob_dist_pdf(1,6,10,12), X = 1/P
  P = 0.000000909599443
  X = 1099385.018181818304583

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


/*
  Probability of getting the sum 12 when throwing 10 fair dice:
  p:0.000000909599 (about 1 in 1099385.02)

*/
go ?=>
  println("Probability of getting the sum 12 when throwing 10 fair dice:"),
  P = sum_prob_dist_pdf(1,6,10,12),
  printf("p:%.12f (about 1 in %.2f)\n",P,1/P),
  nl.
go => true.

/*
  Probability that throwing 4 dice give the sum 12

  [n = 4,s = 12]
  var : ss
  Probabilities:
  14: 0.1157000000000000
  13: 0.1031000000000000
  15: 0.1002000000000000
  12: 0.0980000000000000
  16: 0.0959000000000000
  11: 0.0835000000000000
  17: 0.0799000000000000
  18: 0.0640000000000000
  10: 0.0623000000000000
  19: 0.0427000000000000
  9: 0.0404000000000000
  8: 0.0310000000000000
  20: 0.0293000000000000
  7: 0.0161000000000000
  21: 0.0158000000000000
  22: 0.0086000000000000
  6: 0.0066000000000000
  23: 0.0031000000000000
  5: 0.0024000000000000
  24: 0.0008000000000000
  4: 0.0006000000000000
  mean = 14.0124
  stdev = 3.44091

  var : p
  Probabilities:
  false: 0.9020000000000000
  true: 0.0980000000000000
  mean = 0.098
  stdev = 0.297315

  theoretical_mean = 14.0
  theoretical_stdev = 3.41565
  theoretical_prob = 0.0964506

  All probs:
  14 0.112654320987654 
  15 0.108024691358025 
  13 0.108024691358025 
  16 0.096450617283951 
  12 0.096450617283951 
  17 0.080246913580247 
  11 0.080246913580247 
  18 0.061728395061728 
  10 0.061728395061728 
  19 0.04320987654321 
  9 0.04320987654321 
  20 0.027006172839506 
  8 0.027006172839506 
  21 0.015432098765432 
  7 0.015432098765432 
  22 0.007716049382716 
  6 0.007716049382716 
  23 0.003086419753086 
  5 0.003086419753086 
  24 0.000771604938272 

*/
go2 ?=>
  N = 4,
  S = 12,
  println([n=N,s=S]),
  reset_store,
  run_model(10_000,$model2(N,S),[show_probs,mean,stdev % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  println(theoretical_mean=sum_prob_dist_mean(1,6,N)),
  println(theoretical_stdev=sum_prob_dist_variance(1,6,N).sqrt),
  println(theoretical_prob=sum_prob_dist_pdf(1,6,N,S)),
  println("All probs:"),
  pdf_all($sum_prob_dist(1,6,4),0.001,0.9999999).sort_down(2).printf_list,
  
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2(N,S) =>

  X = random_integer1_n(6,N),
  Ss = X.sum,
  P = check(Ss = S),

  add("ss",Ss),
  add("p",P).