/* 

  Poisson distribution in Picat.

  Please see ppl_distributions.pi and ppl_distributions_test.pi 
  for more on this.

  Cf my Gamble model gamble_poisson_dist.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 : d
  Probabilities:
  3: 0.1995500000000000
  4: 0.1988000000000000
  5: 0.1527000000000000
  2: 0.1491000000000000
  6: 0.1005000000000000
  1: 0.0717500000000000
  7: 0.0582500000000000
  8: 0.0302500000000000
  0: 0.0194000000000000
  9: 0.0129000000000000
  10: 0.0043000000000000
  11: 0.0016500000000000
  12: 0.0006000000000000
  13: 0.0002000000000000
  14: 0.0000500000000000
  mean = 3.9678


  theoretical_mean = 4

  Theretical PDF:
  0 0.018315638888734 
  1 0.073262555554937 
  2 0.146525111109873 
  3 0.195366814813165 
  4 0.195366814813165 
  5 0.156293451850531 
  6 0.104195634567021 
  7 0.059540362609726 
  8 0.029770181304863 
  9 0.01323119169105 
  10 0.00529247667642 
  11 0.001924536973244 
  12 0.000641512324415 
  13 0.000197388407512 
  14 0.000056396687861 
  15 0.000015039116763 


*/
go ?=>
  reset_store,
  run_model(20_000,$model,[show_probs,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  println(theoretical_mean=poisson_dist_mean(4)),
  nl,
  println("Theretical PDF:"),
  pdf_all($poisson_dist(4),0.01,0.99999).printf_list,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  Lambda = 4,
  
  D = poisson_dist(Lambda),
  
  add("d",D).

/*
  Parameter recover 100 samples of poisson_dist(4)

  Num accepted samples: 1000  Total samples: 23415  (0.043%)
  var : lambda
  Probabilities (truncated):
  6.565701950325492: 0.0010000000000000
  6.436958083015381: 0.0010000000000000
  6.401441737050861: 0.0010000000000000
  6.358992387884758: 0.0010000000000000
  .........
  1.797274959044193: 0.0010000000000000
  1.760845123557768: 0.0010000000000000
  1.749620875399383: 0.0010000000000000
  1.62950554729882: 0.0010000000000000
  mean = 3.98422
  HPD intervals:
  HPD interval (0.94): 1.95055225227054..5.92848645988781



*/
go2 ?=>
  Data = poisson_dist_n(4,100),
  reset_store,
  run_model(10_000,$model2(Data),[show_probs_trunc,mean,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2(Data) =>
  Lambda = uniform(0.01,100),
  
  X = poisson_dist_n(Lambda,Data.len),

  observe_abc(Data,X),

  if observed_ok then
    add("lambda",Lambda)
  end.