/* 

  Unbiased_die in Picat.

  From https://www.reddit.com/r/askmath/comments/1hr5wxs/check_whether_the_die_is_unbiased_with_hypothesis/
  """
  Here is a problem of hypothesis which took me almost 2 hours to complete because 
  i was confused as the level of significance wasn't given but somewhere i find out we 
  can simply get it by calculating 1-(confidence interval).

  Can somebody check whether the solution given in image 2 is correct or not. 
  Plus isn't the integral given wrong in the image 1 as the exponential should 
  be e-(x^2/2) dx so i assume that's a printing mistake.

  [The problem in the first image is
     15. A die was thrown 9000 times and of these 3230 yielded a 3 or 4. Is this consistent
         with the hypotheses that the die was unbiased? [0.95]
  ]

  """


  Cf my Gamble model gamble_unbiased_die.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 (truncated):
  2987: 0.0104000000000000
  3020: 0.0102000000000000
  3007: 0.0100000000000000
  2995: 0.0100000000000000
  .........
  2840: 0.0001000000000000
  2839: 0.0001000000000000
  2837: 0.0001000000000000
  2831: 0.0001000000000000
  mean = 2999.67
  [len = 10000,min = 2831,mean = 2999.67,median = 2999.0,max = 3163,variance = 2015.51,stdev = 44.8944]

  var : p
  Probabilities:
  false: 1.0000000000000000
  mean = 0.0
  show_simple_stats: no numeric data

  var : p2
  Probabilities:
  false: 0.9370000000000001
  true: 0.0630000000000000
  mean = 0.063
  show_simple_stats: no numeric data

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,
                           show_simple_stats % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.84,0.9,0.94,0.99,0.99999],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  Throws = random_integer1_n(6,9000),
  S = [cond( (T == 3 ; T == 4),1,0) : T in Throws].sum,

  Target = 3230,
  P = check(S == Target),
  P2 = check((S >= Target * 0.95, S <= Target * 1.05)),
  
  add("s",S),
  add("p",P),
  add("p2",P2).