/* 

  Logistic distribution (parameter recovering) in Picat.



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

/*
  Recover paramer for logistic_dist(3,2)

  Num accepted samples: 1000  Total samples: 744281  (0.001%)
  var : mu
  Probabilities (truncated):
  5.082916042340415: 0.0010000000000000
  4.866902653531592: 0.0010000000000000
  4.824978851166079: 0.0010000000000000
  4.662576096440933: 0.0010000000000000
  .........
  2.52567454358827: 0.0010000000000000
  2.489688509371918: 0.0010000000000000
  2.432859503865642: 0.0010000000000000
  2.26718247042372: 0.0010000000000000
  mean = 3.58611

  var : s
  Probabilities (truncated):
  2.572457609964748: 0.0010000000000000
  2.545323326506337: 0.0010000000000000
  2.539887452749483: 0.0010000000000000
  2.524290323501588: 0.0010000000000000
  .........
  1.448890562843946: 0.0010000000000000
  1.440243200138325: 0.0010000000000000
  1.414347608301019: 0.0010000000000000
  1.388094342028766: 0.0010000000000000
  mean = 1.95265

  var : post
  Probabilities (truncated):
  17.355392376230892: 0.0010000000000000
  15.925722448558258: 0.0010000000000000
  15.175394324288625: 0.0010000000000000
  14.69376610581139: 0.0010000000000000
  .........
  -6.949520075130648: 0.0010000000000000
  -8.179190198360644: 0.0010000000000000
  -12.818360735451639: 0.0010000000000000
  -13.134062925846006: 0.0010000000000000
  mean = 3.64959


*/
go ?=>
  Data = logistic_dist_n(3,2,100),
  println(Data),
  println([min=Data.min,mean=Data.mean,max=Data.max,stdev=Data.stdev]),
  reset_store,
  run_model(10_000,$model(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.
go => true.
  
  
model(Data) =>
  Mu = uniform(0,20),
  S  = uniform(0,10),

  X = logistic_dist_n(Mu,S,Data.len),

  observe_abc(Data,X,1/10),

  Post = logistic_dist(Mu,S),

  if observed_ok then
    add("mu",Mu),
    add("s",S),
    add("post",Post),    
  end.