/* 

  Bernoilli test in Picat.

  Testing how well we can recover probabilities for bernoulli

  This is a port of my Gamble model gamble_bernoulli_test.rkt
  though it has a smaller sample (10 values instead of 120).

  Here we use 10 values with prob 0.3
  Picat> Ys=[bernoulli_dist(0.3) : _ in 1..10]
  Ys = [0,1,0,0,0,0,0,1,0,0]

  

  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.

/*

  It did recover the prob parameter fairly well.

  var : num_same
  Probabilities (truncated):
  5: 0.2054500000000000
  4: 0.1857500000000000
  6: 0.1855000000000000
  3: 0.1290500000000000
  .........
  9: 0.0177000000000000
  1: 0.0160000000000000
  0: 0.0021000000000000
  10: 0.0019000000000000

  mean = 5.01675

  var : prior
  Probabilities (truncated):
  0.891564359153024: 0.0263157894736842
  0.795149171391658: 0.0263157894736842
  0.725779756732261: 0.0263157894736842
  0.720484483181122: 0.0263157894736842
  .........
  0.282926108345996: 0.0263157894736842
  0.241769428115255: 0.0263157894736842
  0.161760836285106: 0.0263157894736842
  0.120804841814513: 0.0263157894736842

  mean = 0.495028

  var : prob
  Probabilities (truncated):
  0.586519264226658: 0.0263157894736842
  0.513038714464109: 0.0263157894736842
  0.489595867297162: 0.0263157894736842
  0.476674781545027: 0.0263157894736842
  .........
  0.132281745038396: 0.0263157894736842
  0.129021943920673: 0.0263157894736842
  0.128581713541266: 0.0263157894736842
  0.118263103797334: 0.0263157894736842

  mean = 0.305155

*/
go ?=>
  % Generated by 
  % Picat> Ys=[bernoulli_dist(0.3) : _ in 1..10]
  Ys = [0,1,0,0,0,0,0,1,0,0],
  reset_store,
  run_model(20_000,$model(Ys),[show_probs_trunc,mean]),
  fail,  
  nl.
go => true.
  

model(Ys) =>
  Len = Ys.len,
  Prob = beta_dist(4,4),  % The probability
  Prior = beta_dist(4,4), % The prior for comparison

  Y = [bernoulli_dist(Prob) : _ in 1..Len],

  % Observe the values in Ys
  C = [cond(Y[I] == Ys[I],1,0) : I  in 1..Len].sum,
  observe(C == Len),
  add("num_same",C),
  if observed_ok then
    add("prob",Prob),
    add("prior",Prior)
  end.