/* 

  BUGS book, 3.5.1 in Picat.

  Example 3.5.1 Heart transplants: learning from data

  Cf my WebPPL model bugs_book_3_5_1.wppl

  WebPPL:
  expectation:
  [ [ 'Is', 3.0248484423091933 ],
    [ 'pT', 0.7564845131501468 ],
    [ 'surv_t', 5.024848442309195 ],
    [ 'theta', 0.1726172625722976 ] ]

  BUGS book:
  Output:
           Mean      SD  Naive SE Time-series SE
  Is     3.1442 2.28250 0.0080699      0.0085034
  pT     0.7500 0.12014 0.0004248      0.0005856
  surv.t 5.1442 2.28250 0.0080699      0.0085034
  theta  0.1667 0.05899 0.0002086      0.0002086
 

  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 : Is
  Probabilities (truncated):
  10.42085849410797: 0.0100000000000000
  9.062668029043822: 0.0100000000000000
  8.771809863651031: 0.0100000000000000
  8.118882067018143: 0.0100000000000000
  .........
  0.156838271301281: 0.0100000000000000
  0.137775989062067: 0.0100000000000000
  0.113719415183676: 0.0100000000000000
  -0.546773870833125: 0.0100000000000000
  mean = 3.34637
  [len = 100,min = -0.546774,mean = 3.34637,median = 2.89538,max = 10.4209,variance = 5.27879,stdev = 2.29756]
  HPD intervals:
  HPD interval (0.84): 0.11371941518368..5.90814412133397
  HPD interval (0.9): 0.11371941518368..6.64740879651361
  HPD interval (0.94): 0.13777598906207..7.51347613239676
  HPD interval (0.99): -0.54677387083313..9.06266802904382
  HPD interval (0.99999): -0.54677387083313..10.42085849410797

  var : Pt
  Probabilities (truncated):
  0.973809913719916: 0.0100000000000000
  0.951466801553716: 0.0100000000000000
  0.938170963403848: 0.0100000000000000
  0.931638691542502: 0.0100000000000000
  .........
  0.424981740501235: 0.0100000000000000
  0.42201121264231: 0.0100000000000000
  0.417404146593718: 0.0100000000000000
  0.374692314478891: 0.0100000000000000
  mean = 0.730436
  [len = 100,min = 0.374692,mean = 0.730436,median = 0.745311,max = 0.97381,variance = 0.0164709,stdev = 0.128339]
  HPD intervals:
  HPD interval (0.84): 0.55321540662703..0.89184548095420
  HPD interval (0.9): 0.51400546613802..0.90587008507264
  HPD interval (0.94): 0.51400546613802..0.95146680155372
  HPD interval (0.99): 0.41740414659372..0.97380991371992
  HPD interval (0.99999): 0.37469231447889..0.97380991371992

  var : Surv t
  Probabilities (truncated):
  12.42085849410797: 0.0100000000000000
  11.062668029043822: 0.0100000000000000
  10.771809863651031: 0.0100000000000000
  10.118882067018143: 0.0100000000000000
  .........
  2.156838271301281: 0.0100000000000000
  2.137775989062067: 0.0100000000000000
  2.113719415183676: 0.0100000000000000
  1.453226129166875: 0.0100000000000000
  mean = 5.34637
  [len = 100,min = 1.45323,mean = 5.34637,median = 4.89538,max = 12.4209,variance = 5.27879,stdev = 2.29756]
  HPD intervals:
  HPD interval (0.84): 2.11371941518368..7.90814412133397
  HPD interval (0.9): 2.11371941518368..8.64740879651361
  HPD interval (0.94): 2.13777598906207..9.51347613239676
  HPD interval (0.99): 1.45322612916687..11.06266802904382
  HPD interval (0.99999): 1.45322612916687..12.42085849410797

  var : theta
  Probabilities (truncated):
  0.441672538617947: 0.0100000000000000
  0.379428660994163: 0.0100000000000000
  0.376282886304965: 0.0100000000000000
  0.31835344559027: 0.0100000000000000
  .........
  0.071988610505216: 0.0100000000000000
  0.070116002223914: 0.0100000000000000
  0.067509992253977: 0.0100000000000000
  0.065934097668466: 0.0100000000000000
  mean = 0.159796
  [len = 100,min = 0.0659341,mean = 0.159796,median = 0.150471,max = 0.441673,variance = 0.00465593,stdev = 0.0682344]
  HPD intervals:
  HPD interval (0.84): 0.06750999225398..0.21130647016857
  HPD interval (0.9): 0.06593409766847..0.24469006262517
  HPD interval (0.94): 0.06593409766847..0.26930956500188
  HPD interval (0.99): 0.06593409766847..0.37942866099416
  HPD interval (0.99999): 0.06593409766847..0.44167253861795

*/
go ?=>
  SPs = [2.0,3.0,4.0,4.0,6.0,7.0,10.0,12.0], %  data	
  reset_store,
  run_model(10_000,$model(SPs),[show_probs_trunc,mean,
                           show_simple_stats,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84,0.9,0.94,0.99,0.99999],
                           % show_histogram,
                           min_accepted_samples=100 % ,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(Data) =>
  % Fixed integer YT = 8,
  
  NT = 10,

  Theta = gamma_dist(0.1,0.1),
  PT = uniform(0,1),
  YT = binomial_dist(NT,PT),

  SurvT = PT/Theta,
  Is = SurvT - 2.0,

  observe(YT == 8),

  SP = exponential_dist_n(Theta,NT),
  observe_abc(Data,SP),

  if observed_ok then
    add("Is",Is),
    add("Pt",PT), 
    add("Surv t",SurvT),
    add("theta",Theta),    
  end.