/* 

  Log normal distribition in Picat.

  See more in ppl_distributions.pi and ppl_distributions_test.pi

  Cf my Gamble model gamble_log_normal_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 (truncated):
  44.552124304365236: 0.0001000000000000
  39.643417171204156: 0.0001000000000000
  36.986110183013594: 0.0001000000000000
  36.857983614293609: 0.0001000000000000
  .........
  0.028361022214494: 0.0001000000000000
  0.027138891224294: 0.0001000000000000
  0.017881964874886: 0.0001000000000000
  0.012377024527369: 0.0001000000000000
  mean = 1.64822
  HPD intervals:
  HPD interval (0.84): 0.06779310158896..2.68707467938546

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean ,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84]
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  D = lognormal_dist(0,1),
  add("d",D).


/*

  Simple parameter recovery. 

  Not very bad...

  Num accepted samples: 100  Total samples: 63669  (0.002%)

  var : mu
  Probabilities (truncated):
  0.591430636398229: 0.0100000000000000
  0.58986401212861: 0.0100000000000000
  0.558204395025132: 0.0100000000000000
  0.55303592726264: 0.0100000000000000
  .........
  0.019673230135661: 0.0100000000000000
  0.014925869188703: 0.0100000000000000
  0.012605395173936: 0.0100000000000000
  0.00899540260853: 0.0100000000000000
  mean = 0.265274
  HPD intervals:
  HPD interval (0.84): 0.01967323013566..0.45142157955627

  var : sigma
  Probabilities (truncated):
  1.504707677059205: 0.0100000000000000
  1.415654565866876: 0.0100000000000000
  1.369353864048307: 0.0100000000000000
  1.223371080692565: 0.0100000000000000
  .........
  0.651920801332184: 0.0100000000000000
  0.650396822323276: 0.0100000000000000
  0.636552772781138: 0.0100000000000000
  0.632888637777832: 0.0100000000000000
  mean = 0.930499
  HPD intervals:
  HPD interval (0.84): 0.63288863777783..1.06657702525453

  var : post
  Probabilities (truncated):
  10.287474903283449: 0.0100000000000000
  8.890338983898724: 0.0100000000000000
  8.251583599661419: 0.0100000000000000
  4.494437017137865: 0.0100000000000000
  .........
  0.269698628631155: 0.0100000000000000
  0.252143603615873: 0.0100000000000000
  0.172346119071202: 0.0100000000000000
  0.130656494535122: 0.0100000000000000
  mean = 1.73019
  HPD intervals:
  HPD interval (0.84): 0.13065649453512..2.88870135715500


  Mathematica: 
  FindDistributionParameters[data, LogNormalDistribution[mu, sigma]] 
  -> {mu -> -0.142795, sigma -> 1.04496}



*/
go2 ?=>
  Data = lognormal_dist_n(0,1,100),
  println(data=Data),
  show_simple_stats(Data),
  reset_store,
  run_model(10_000,$model2(Data),[show_probs_trunc,mean,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84],
                           % show_histogram,
                           min_accepted_samples=100,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2(Data) =>
  Mu = uniform(0,10),
  Sigma = uniform(0,10),

  Y = lognormal_dist_n(Mu,Sigma,Data.len),
  observe_abc(Data,Y,1/3),

  Post = lognormal_dist(Mu,Sigma),

  if observed_ok then
    add("mu",Mu),
    add("sigma",Sigma),
    add("post",Post),
  end.