/* 

  Bugs book example 3.3.2 in Picat.

  Example 3.3.2 Surgery (continued): beta-binomial analysis using BUGS

  Cf my Gamble model gamble_bugs_book_3_3_2.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 : p crit
  Probabilities:
  false: 0.5809334657398213
  true: 0.4190665342601788
  mean = [false = 0.580933,true = 0.419067]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : y pred
  Probabilities (truncated):
  1: 0.3041211519364449
  0: 0.2768123138033763
  2: 0.2043197616683217
  3: 0.1226415094339623
  .........
  7: 0.0037239324726912
  8: 0.0014895729890765
  9: 0.0004965243296922
  10: 0.0002482621648461
  mean = 1.51117
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..3.00000000000000

  var : theta
  Probabilities (truncated):
  0.30100614472645: 0.0002482621648461
  0.287229490072522: 0.0002482621648461
  0.259438490754922: 0.0002482621648461
  0.258046091942835: 0.0002482621648461
  .........
  0.003930707878165: 0.0002482621648461
  0.003825830817417: 0.0002482621648461
  0.002953661229263: 0.0002482621648461
  0.002185704328694: 0.0002482621648461
  mean = 0.0742754
  HPD intervals:
  HPD interval (0.84): 0.01730018004156..0.12063067019741

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_hpd_intervals,hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  
  A = 3,
  B = 27,
  N = 10,
  NPred = 20,
  NCrit = 2,

  Theta = beta_dist(A,B),
  
  Y = binomial_dist(N, Theta),
  YPred = binomial_dist(NPred, Theta),
  PCrit = check( (YPred + 0.5) > NCrit), %  toReal(y_pred - n_crit) + 0.5 > 0.0;

  observe(Y == 0),
    
  if observed_ok then
    add("p crit",PCrit),
    add("y pred",YPred),
    add("theta",Theta),        
  end.