/* 

  Binomial coin in Picat.

  From https://reference.wolfram.com/language/ref/BinomialDistribution.html
  """
  Compute the probability that there are between 60 and 80 heads in 100 coin flips.

  Probability[60 <= x <= 80, x -> heads[100]]
  -> 0.028444

  ... and an unfair coin (p:0.6)
  Probability[60 <= x <= 80, x -> uheads[100]]
  -> 0.543289
  
  """

  Picat>  X=(1-binomial_dist_cdf(100,1/2,59))- (1-binomial_dist_cdf(100,1/2,80))
  X = 0.028443966685352

  Picat> X=(1-binomial_dist_cdf(100,0.6,59))- (1-binomial_dist_cdf(100,0.6,80))
  X = 0.543288604557029


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

/*

  p = 0.5
  var : b
  Probabilities (truncated):
  49: 0.0815000000000000
  51: 0.0806000000000000
  50: 0.0724000000000000
  48: 0.0710000000000000
  .........
  33: 0.0005000000000000
  68: 0.0002000000000000
  32: 0.0002000000000000
  66: 0.0001000000000000
  mean = 49.9393

  var : prob
  Probabilities:
  false: 0.9729000000000000
  true: 0.0271000000000000
  mean = [false = 0.9729,true = 0.0271]


  p = 0.6
  var : b
  Probabilities (truncated):
  59: 0.0850000000000000
  60: 0.0844000000000000
  62: 0.0785000000000000
  61: 0.0756000000000000
  .........
  44: 0.0004000000000000
  76: 0.0001000000000000
  40: 0.0001000000000000
  38: 0.0001000000000000
  mean = 59.9693

  var : prob
  Probabilities:
  true: 0.5424000000000000
  false: 0.4576000000000000
  mean = [true = 0.5424,false = 0.4576]

*/
go ?=>
  member(P,[1/2,0.6]),
  println(p=P),
  reset_store,
  run_model(10_000,$model(P),[show_probs_trunc,mean % ,
                           % show_percentiles,show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(P) =>

  B = binomial_dist(100,P),
  
  Prob = check( (60 <= B, B <= 80) ),
  
  add("b",B),
  add("prob",Prob).