/* 

  8 boys and 2 girls in Picat.

  A friend of mine has 10 grandchildren: 8 boys and 2 girls. 
  How unusual/surprising it this?

  Let's assume that the birth of girls is slighly more common than
  the birth of boys: 0.51 vs 0.49.

  The probability of 8 boys of 10 children is about 3.9%.
  * p-binomial: the exact result from (binomial 10 0.49): 0.038897483585189484


  Picat> X=binomial_dist_pdf(10,0.49,8)                               
  X = 0.038897483585189

  Picat> pdf_all($binomial_dist(10,0.49)).sort_down(2).printf_list    
   5 0.245601956092531 
   4 0.213022104774134 
   6 0.196642089028334 
   3 0.126695362606191 
   7 0.107960362603791 
   2 0.04944997571109 
   8 0.038897483585189 
   1 0.011437409348143 


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

/*

  [p_boy = 0.49,p_girl = 0.51,n = 10,obs_num_boys = 8]
  p_exact = 0.0388975
  var : num boys
  Probabilities:
  5: 0.2431000000000000
  4: 0.2165000000000000
  6: 0.1981000000000000
  3: 0.1272000000000000
  7: 0.1077000000000000
  2: 0.0475000000000000
  8: 0.0347000000000000
  1: 0.0142000000000000
  9: 0.0087000000000000
  0: 0.0014000000000000
  10: 0.0009000000000000
  mean = 4.8797
  Percentiles:
  (0.001 0)
  (0.01 1)
  (0.025 2)
  (0.05 2)
  (0.1 3)
  (0.25 4)
  (0.5 5)
  (0.75 6)
  (0.84 6)
  (0.9 7)
  (0.95 7)
  (0.975 8)
  (0.99 8)
  (0.999 9)
  (0.9999 10)
  (0.99999 10)
  HPD intervals:
  HPD interval (0.94): 2.00000000000000..7.00000000000000

  var : num girls
  Probabilities:
  5: 0.2431000000000000
  6: 0.2165000000000000
  4: 0.1981000000000000
  7: 0.1272000000000000
  3: 0.1077000000000000
  8: 0.0475000000000000
  2: 0.0347000000000000
  9: 0.0142000000000000
  1: 0.0087000000000000
  10: 0.0014000000000000
  0: 0.0009000000000000
  mean = 5.1203
  Percentiles:
  (0.001 1)
  (0.01 2)
  (0.025 2)
  (0.05 3)
  (0.1 3)
  (0.25 4)
  (0.5 5)
  (0.75 6)
  (0.84 7)
  (0.9 7)
  (0.95 8)
  (0.975 8)
  (0.99 9)
  (0.999 10)
  (0.9999 10)
  (0.99999 10)
  HPD intervals:
  HPD interval (0.94): 3.00000000000000..8.00000000000000

  var : p
  Probabilities:
  false: 0.9653000000000000
  true: 0.0347000000000000
  mean = [false = 0.9653,true = 0.0347]
  Percentiles:
  (0.001 false)
  (0.01 false)
  (0.025 false)
  (0.05 false)
  (0.1 false)
  (0.25 false)
  (0.5 false)
  (0.75 false)
  (0.84 false)
  (0.9 false)
  (0.95 false)
  (0.975 true)
  (0.99 true)
  (0.999 true)
  (0.9999 true)
  (0.99999 true)
  HPD intervals:
  show_hpd_intervals: data is not numeric



*/
go ?=>
  PBoy = 0.49,
  PGirl = 1-PBoy,
  N = 10,
  ObsNumBoys = 8,
  println([p_boy=PBoy,p_girl=PGirl,n=N,obs_num_boys=ObsNumBoys]),
  println(p_exact=binomial_dist_pdf(N,PBoy,ObsNumBoys)),
  reset_store,
  run_model(10_000,$model(PBoy,PGirl,N,ObsNumBoys),
                          [show_probs,mean,show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.94]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(PBoy,PGirl,N,ObsNumBoys) =>
  Children = categorical_n([PBoy,PGirl],[boy,girl],N),

  NumBoys = count_occurrences(boy,Children),
  NumGirls = count_occurrences(girl,Children),

  P = check(NumBoys == ObsNumBoys),
  
  add("num boys",NumBoys),
  add("num girls",NumGirls),
  add("p",P).