/* 

  HPD intervals in Picat.

  HPD (Highest Posterior Density).

  For comparison, also shown is the - generally wider - exact credible 
  intervals  for a probability distribution (i.e. not the data).

  Cf my Gamble model gamble_credible_interval.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.


/*

  dist = poisson_dist(10)
  Compare with exact (and the wider) credible Intervals for the distribution:
  Credible intervals (exact):
  0.5: 8..12
  0.84: 6..15
  0.9: 5..15
  0.94: 5..16
  0.95: 4..17
  0.99: 3..19
  0.999: 2..22
  0.9999: 1..24
  0.99999: 0..27

  var : p
  Probabilities:
  9: 0.1285000000000000
  10: 0.1249000000000000
  8: 0.1149000000000000
  11: 0.1074000000000000
  7: 0.0944000000000000
  12: 0.0933000000000000
  13: 0.0699000000000000
  6: 0.0629000000000000
  14: 0.0555000000000000
  5: 0.0362000000000000
  15: 0.0317000000000000
  16: 0.0219000000000000
  4: 0.0201000000000000
  17: 0.0138000000000000
  3: 0.0087000000000000
  18: 0.0073000000000000
  19: 0.0033000000000000
  20: 0.0017000000000000
  2: 0.0015000000000000
  1: 0.0006000000000000
  23: 0.0005000000000000
  21: 0.0005000000000000
  22: 0.0004000000000000
  24: 0.0001000000000000
  mean = 9.9625
  Percentiles:
  (0.001 2)
  (0.01 3)
  (0.025 4)
  (0.05 5)
  (0.1 6)
  (0.25 8)
  (0.5 10)
  (0.75 12)
  (0.84 13)
  (0.9 14)
  (0.95 15)
  (0.975 17)
  (0.99 18)
  (0.999 21.001000000000204)
  (0.9999 23.000099999999293)
  (0.99999 23.900010000001203)
  HPD intervals:
  HPD interval (0.5): 6.00000000000000..10.00000000000000
  HPD interval (0.84): 6.00000000000000..14.00000000000000
  HPD interval (0.9): 4.00000000000000..14.00000000000000
  HPD interval (0.95): 4.00000000000000..16.00000000000000
  HPD interval (0.99): 3.00000000000000..18.00000000000000
  HPD interval (0.999): 1.00000000000000..21.00000000000000
  HPD interval (0.9999): 1.00000000000000..23.00000000000000
  HPD interval (0.99999): 1.00000000000000..24.00000000000000

  dist = binomial_dist(10,1 / 2)
  Compare with exact (and the wider) credible Intervals for the distribution:
  Credible intervals (exact):
  0.5: 4..6
  0.84: 3..7
  0.9: 2..8
  0.94: 2..8
  0.95: 2..8
  0.99: 1..9
  0.999: 0..10
  0.9999: 0..10
  0.99999: 0..10


  var : p
  Probabilities:
  5: 0.2549000000000000
  6: 0.2020000000000000
  4: 0.1979000000000000
  3: 0.1191000000000000
  7: 0.1178000000000000
  2: 0.0435000000000000
  8: 0.0423000000000000
  9: 0.0112000000000000
  1: 0.0088000000000000
  0: 0.0018000000000000
  10: 0.0007000000000000
  mean = 5.002
  Percentiles:
  (0.001 0)
  (0.01 1)
  (0.025 2)
  (0.05 2)
  (0.1 3)
  (0.25 4)
  (0.5 5)
  (0.75 6)
  (0.84 7)
  (0.9 7)
  (0.95 8)
  (0.975 8)
  (0.99 9)
  (0.999 9)
  (0.9999 10)
  (0.99999 10)
  HPD intervals:
  HPD interval (0.5): 3.00000000000000..5.00000000000000
  HPD interval (0.84): 3.00000000000000..7.00000000000000
  HPD interval (0.9): 2.00000000000000..7.00000000000000
  HPD interval (0.95): 2.00000000000000..8.00000000000000
  HPD interval (0.99): 1.00000000000000..9.00000000000000
  HPD interval (0.999): 0.00000000000000..9.00000000000000
  HPD interval (0.9999): 0.00000000000000..10.00000000000000
  HPD interval (0.99999): 0.00000000000000..10.00000000000000


  dist = normal_dist(100,15)
  Compare with exact (and the wider) credible Intervals for the distribution:
  Credible intervals (exact):
  0.5: 89.882653747058782 .. 110.117346252941218
  0.84: 78.923926595355525 .. 121.076073404644475
  0.9: 75.327195595727929 .. 124.672804404272114
  0.94: 71.788095877731266 .. 128.211904122268720
  0.95: 70.600540231899288 .. 129.399459768100712
  0.99: 61.362560446766857 .. 138.637439553233094
  0.999: 50.642099027621875 .. 149.357900972377649
  0.9999: 41.641121703774118 .. 158.358878296229932
  0.99999: 33.742398798071591 .. 166.257601202072436

  var : p
  Probabilities (truncated):
  158.610359134664918: 0.0001000000000000
  158.134249012264092: 0.0001000000000000
  157.873100904502053: 0.0001000000000000
  153.477922536420692: 0.0001000000000000
  .........
  51.807670128362297: 0.0001000000000000
  51.532781546631888: 0.0001000000000000
  50.621582277142011: 0.0001000000000000
  49.518638095060837: 0.0001000000000000
  mean = 100.106
  Percentiles:
  (0.001 55.760500359643657)
  (0.01 65.332873519774722)
  (0.025 70.80189606177467)
  (0.05 75.428173045179975)
  (0.1 80.996693456516553)
  (0.25 89.751752848664125)
  (0.5 100.12005035012571)
  (0.75 110.21677191029204)
  (0.84 115.19124519577575)
  (0.9 119.59488559886002)
  (0.95 125.00342386502955)
  (0.975 129.58817289345365)
  (0.99 135.21310166553471)
  (0.999 148.02955903875446)
  (0.9999 158.13429662327599)
  (0.99999 158.56275288352663)
  HPD intervals:
  HPD interval (0.5): 89.56394742037000..110.00701385233180
  HPD interval (0.84): 79.30730498008219..121.43089605041530
  HPD interval (0.9): 75.71521858091357..125.19909243241395
  HPD interval (0.95): 71.14681140157236..129.73869469446538
  HPD interval (0.99): 61.36552118100384..138.55944778211261
  HPD interval (0.999): 53.20785669963710..149.73929813502934
  HPD interval (0.9999): 50.62158227714201..158.61035913466492
  HPD interval (0.99999): 49.51863809506084..158.61035913466492


*/
go ?=>
  member([Dist,Type], [ [$poisson_dist(10),discrete],
                        [$binomial_dist(10,1/2),discrete],
                        [$normal_dist(100,15),continuous]
                       ]),
  println(dist=Dist),
  println("Compare with exact (and the wider) credible Intervals for the distribution:"),
  % Qs = quantile_qs(),  
  Qs = [0.5,0.84,0.90,0.94,0.95,0.99,0.999,0.9999,0.99999],
  show_credible_intervals_exact(Dist,Qs),
  reset_store,
  if Type == discrete then
    run_model(10_000,$model(Dist),[show_probs,mean,show_percentiles,
                                   show_hpd_intervals,hpd_intervals=Qs                                
                                   ])
  else
    run_model(10_000,$model(Dist),[show_probs_trunc,mean,show_percentiles,
                                   show_hpd_intervals,hpd_intervals=Qs])  
  end,
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.


model(Dist) =>
  P = apply(Dist),
  add("p",P).