/* 

  The Euro coin problem in Picat.

  From Think Bayes, page 33ff
  """
  A statistical statement appeared in "The Guardian" on Friday January 4, 2002:
      When spun on edge 250 times, a Belgian one-euro coin
      came up heads 140 times and tails 110. 'It looks very
      suspicious to me,' said Barry Blight, a statistics lecturer
      at the London School of Economics. 'If the coin were
      unbiased, the chance of getting a result as extreme as
      that would be less than 7%.'

  But do these data give evidence that the coin is biased rather than fair?
  """

  Below are two models, the first uses bernoulli_dist, the 
  second uses binomial dist. 
  (There's no discernible speed difference between these model.)

  This is a port of my Gamble model gamble_euro_coin_problem.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.

main => go.

/*
  Model 1:
  var : prob
  Probabilities (truncated):
  0.650123651405981: 0.0053763440860215
  0.646109778660077: 0.0053763440860215
  0.632276532616487: 0.0053763440860215
  0.622607824737545: 0.0053763440860215
  .........
  0.491183652664872: 0.0053763440860215
  0.487510000851834: 0.0053763440860215
  0.479311539132756: 0.0053763440860215
  0.475704321896695: 0.0053763440860215
  mean = 0.556003

  var : prob < 0.5
  Probabilities:
  0: 0.9516129032258065
  1: 0.0483870967741935
  mean = 0.0483871

  var : prob > 0.5
  Probabilities:
  1: 0.9516129032258065
  0: 0.0483870967741935
  mean = 0.951613

  Model 2:
  var : prob
  Probabilities (truncated):
  0.641318741227119: 0.0060606060606061
  0.633225348530738: 0.0060606060606061
  0.623409036170211: 0.0060606060606061
  0.61681555712329: 0.0060606060606061
  .........
  0.491349550077013: 0.0060606060606061
  0.489896305948883: 0.0060606060606061
  0.480348323693898: 0.0060606060606061
  0.478730313870223: 0.0060606060606061
  mean = 0.559071

  var : prob < 0.5
  Probabilities:
  0: 0.9636363636363636
  1: 0.0363636363636364
  mean = 0.0363636

  var : prob > 0.5
  Probabilities:
  1: 0.9636363636363636
  0: 0.0363636363636364
  mean = 0.963636

*/
go ?=>
  reset_store,
  println("Model 1:"),
  run_model(30_000,$model1,[show_probs_trunc,mean]),
  nl,
  reset_store,
  println("Model 2:"),  
  run_model(30_000,$model2,[show_probs_trunc,mean]),
  nl.
go => true.

%
% Using bernoulli dist
%
model1() =>
  N = 250,

  % Probability of throwing head
  Prob = beta_dist(2,2), 

  Coin = bernoulli_dist_n(Prob,N),
  Sum250 = Coin.sum,

  observe(Sum250 == 140),
  if observed_ok then
    add("prob",Prob),  
    add("prob > 0.5",cond(Prob > 0.5,1,0)),
    add("prob < 0.5",cond(Prob < 0.5,1,0))
  end.

%
% Using binomial dist
%
model2() =>
  N = 250,

  % Probability of throwing head
  Prob = beta_dist(2,2), 

  Heads = binomial_dist(N,Prob),

  observe(Heads == 140),
  
  if observed_ok then
    add("prob",Prob),  
    add("prob > 0.5",cond(Prob > 0.5,1,0)),
    add("prob < 0.5",cond(Prob < 0.5,1,0))
  end.