/* 

  Coin toss in Picat.

  From https://people.duke.edu/~ccc14/sta-663/PyStan.html

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

/*

  var : p
  Probabilities (truncated):
  0.76078018558566: 0.00234741784037559 (1 / 426)
  0.748142167818657: 0.00234741784037559 (1 / 426)
  0.742222787487951: 0.00234741784037559 (1 / 426)
  0.736101589472614: 0.00234741784037559 (1 / 426)
  .........
  0.478827697776306: 0.00234741784037559 (1 / 426)
  0.472261720139761: 0.00234741784037559 (1 / 426)
  0.471444745534792: 0.00234741784037559 (1 / 426)
  0.452778173935529: 0.00234741784037559 (1 / 426)

  mean = 0.605242
  HPD intervals:
  HPD interval (0.5): 0.57714416623694..0.63512975544558
  HPD interval (0.84): 0.54897554368472..0.68508296970814
  HPD interval (0.9): 0.52164492350616..0.68238951504908
  HPD interval (0.95): 0.49952265236452..0.69174702640968
  HPD interval (0.99): 0.47144474553479..0.73610158947261
  HPD interval (0.999): 0.45277817393553..0.76078018558566
  HPD interval (0.9999): 0.45277817393553..0.76078018558566
  HPD interval (0.99999): 0.45277817393553..0.76078018558566


*/
go ?=>

  reset_store,
  run_model(30_000,$model,[show_probs_rat_trunc,mean,show_hpd_intervals]),
  nl.
go => true.

model() =>
  N = 100, % Number of tosses
  P = beta_dist(2,2),
  Y = binomial_dist(N,P),

  observe(Y==61),
  if observed_ok()  then
    add("p",P)
  end.