/* 

  Mean and stdev for 10 flips of a fair coin in Picat.

  From http://www.statistics101.net/statistics101web_000007.htm
  """
  From CliffsQuickReview Statistics, p. 47, example 7 
  What is the mean and standard deviation for a 
  binomial probability distribution for 10 flips 
  of a fair coin? 
  Calculated result using binomial formula:  
  mean = 5, standard deviation = 1.58  
  """

  Cf my Gamble model gamble_mean_and_stdev_for_10_flips_of_a_fair_coin.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 : s
  Probabilities:
  5: 0.2436000000000000
  4: 0.2069000000000000
  6: 0.2065000000000000
  3: 0.1161000000000000
  7: 0.1157000000000000
  8: 0.0503000000000000
  2: 0.0408000000000000
  9: 0.0096000000000000
  1: 0.0089000000000000
  0: 0.0009000000000000
  10: 0.0007000000000000
  mean = 5.0291
  stdev = 1.57976
  HPD intervals:
  HPD interval (0.84): 3.00000000000000..7.00000000000000
  HPD interval (0.9): 2.00000000000000..7.00000000000000
  HPD interval (0.95): 2.00000000000000..8.00000000000000
  HPD interval (0.99): 1.00000000000000..9.00000000000000


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs,mean,stdev,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84,0.9,0.95,0.99]
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  S = resample(10,[0,1]).sum,
  
  add("s",S).


/* 
   Using binomial(10,1/2)

   var : s
   Probabilities:
   5: 0.2428000000000000
   6: 0.2114000000000000
   4: 0.2051000000000000
   7: 0.1173000000000000
   3: 0.1133000000000000
   8: 0.0465000000000000
   2: 0.0412000000000000
   1: 0.0107000000000000
   9: 0.0099000000000000
   10: 0.0009000000000000
   0: 0.0009000000000000
   mean = 5.027
   stdev = 1.58179
   HPD intervals:
   HPD interval (0.84): 3.00000000000000..7.00000000000000
   HPD interval (0.9): 2.00000000000000..7.00000000000000
   HPD interval (0.95): 2.00000000000000..8.00000000000000
   HPD interval (0.99): 1.00000000000000..9.00000000000000


   Picat> X=binomial_dist_mean(10,1/2)          
   X = 5.0

   Picat> X=binomial_dist_variance(10,1/2)
   X = 2.5

   Picat> X=sqrt(binomial_dist_variance(10,1/2))
   X = 1.58113883008419


*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs,mean,stdev,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84,0.9,0.95,0.99]
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  S = binomial_dist(10,1/2),
  
  add("s",S).