/* 

  Fair coin in Picat.

  From Statistics101 (Resample Stats)
  http://www.statistics101.net/QuickReference.pdf
  Page 46
  """
  Say that you flipped the coin ten times and got seven heads. Is the coin fair? 
  Your research hypothesis is that the coin is unfair. Therefore, the null hypothesis, 
  which is the opposite of the research hypothesis, and which you are hoping to
  contradict, is that the coin is fair. You then choose a probability criterion which you consider so
  low that if the computed probability were equal to or less than the criterion then it would be
  reasonable to say that the sample could not have come from the benchmark universe. Say you
  choose 0.05. Then you would determine the probability of seven heads out of ten throws from
  that universe and compare it to your criterion.
  ...
  The result comes in around 0.12. Since this is greater than your criterion of 0.05, 
  you must conclude that seven heads out of ten is not sufficiently unusual for a fair 
  coin, so you cannot reject the null hypothesis.
  """

  The probability of getting 7 heads in 10 tosses with a fair coin (p) is about 11.7%,
  so this does not indicate an unfair coin.

  However, if we got 0, 1, 9, or 10 head this would be surprising for a fair coin.
  Note that taken these cases together, i.e. the probability of getting 
  either 0, 1, 9,or 10 coins (p2) is about 2.15% which is not _very_ uncommon.
  

  Picat> X=binomial_dist_quantile(10,1/2,0.01)
  X = 1

  Picat> X=binomial_dist_quantile(10,1/2,0.99)
  X = 9

  Picat> pdf_all($binomial_dist(10,1/2),0.000001,0.999999).sort_down(2).printf_list
   5 0.24609375 
   6 0.205078125 
   4 0.205078125 
   7 0.1171875 
   3 0.1171875 
   8 0.0439453125 
   2 0.0439453125 
   9 0.009765625 
   1 0.009765625 
  10 0.0009765625 
   0 0.0009765625 

  Cf my Gamble model gamble_fair_coin2.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):
  5: 0.2461000000000000
  6: 0.2032000000000000
  4: 0.2016000000000000
  7: 0.1231000000000000
  .........
  9: 0.0103000000000000
  1: 0.0095000000000000
  0: 0.0006000000000000
  10: 0.0004000000000000
  mean = 5.03
  Percentiles:
  (0.001 1)
  (0.01 1)
  (0.025 2)
  (0.05 2)
  (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 9)
  (0.9999 10)
  (0.99999 10)
  HPD intervals:
  HPD interval (0.84): 3.00000000000000..7.00000000000000

  var : p
  Probabilities:
  false: 0.8769000000000000
  true: 0.1231000000000000
  mean = [false = 0.8769,true = 0.1231]
  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 true)
  (0.95 true)
  (0.975 true)
  (0.99 true)
  (0.999 true)
  (0.9999 true)
  (0.99999 true)
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : p2
  Probabilities:
  false: 0.9792000000000000
  true: 0.0208000000000000
  mean = [false = 0.9792,true = 0.0208]
  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 false)
  (0.99 true)
  (0.999 true)
  (0.9999 true)
  (0.99999 true)
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : surprise
  Probabilities (truncated):
  0.0: 0.2461000000000000
  2.041241452319315: 0.1494000000000000
  -2.041241452319315: 0.1466000000000000
  4.364357804719847: 0.1198000000000000
  .........
  13.333333333333334: 0.0053000000000000
  13.333333333333332: 0.0050000000000000
  4.364357804719848: 0.0033000000000000
  0: 0.0010000000000000
  mean = 0.0764451
  Percentiles:
  (0.001 -13.333333333333332)
  (0.01 -7.5)
  (0.025 -7.4999999999999982)
  (0.05 -7.4999999999999982)
  (0.1 -4.3643578047198481)
  (0.25 -2.0412414523193148)
  (0.5 0)
  (0.75 2.0412414523193148)
  (0.84 4.3643578047198472)
  (0.9 4.3643578047198472)
  (0.95 7.4999999999999982)
  (0.975 7.5)
  (0.99 13.333333333333332)
  (0.999 13.333333333333334)
  (0.9999 13.333333333333334)
  (0.99999 13.333333333333334)
  HPD intervals:
  HPD interval (0.84): -4.36435780471985..4.36435780471985


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_percentiles,show_hpd_intervals,hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  Sample = bern_n(1/2,10),
  S = Sample.sum,
  Stdev = Sample.stdev,
  P = check(S == 7),
  P2 = check( (S <= 1 ; S >= 9)),
  Surprise = cond(Stdev != 0.0, (S - 5) / Stdev, 0),
  
  add("s",S),
  add("p",P),
  add("p2",P2),
  add("surprise",Surprise).