/* 

  Laplace distribution in Picat.

  From Handbook on probability distributions
  page 72ff

  """
  Let U be a uniform variate. Then the algorithm is
  * V = U − 1/2
  * X = m + sigma sign(V ) log(1 − 2|V |)
  return X
  """

  See ppl_distributions.pi and ppl_distributions_test.pi for more on this.

  Example from Mathematica's LaplaceDistribution
  """
  The difference of flood stages between river stations A and B in a year has 
  been estimated to follow a Laplace distribution with mu==10 feet and b==3.4 feet. 
  Find the probability that the difference is greater than 15 feet: 
  D = LaplaceDistribution[10, 3.4]

  Probability[Abs[u] > 15,  u e D]
  -> 
  0.115215
  """

  Some exact calculations:
  Picat> X=laplace_dist_mean(10,3.4)    
  X = 10

  Picat> X=1-laplace_dist_cdf(10,3.4,15)+laplace_dist_cdf(10,3.4,-15)
  X = 0.115215489916249


  Cf my Gamble model gamble_laplace_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):
  38.472402154593304: 0.0001000000000000
  35.070469307367851: 0.0001000000000000
  34.640604757623748: 0.0001000000000000
  33.904177209639727: 0.0001000000000000
  .........
  -14.253876088436119: 0.0001000000000000
  -17.326734457709325: 0.0001000000000000
  -17.902115617306077: 0.0001000000000000
  -19.983557836785142: 0.0001000000000000
  mean = 10.0732
  Percentiles:
  (0.001 -9.5254770896875929)
  (0.01 -2.9064160033629989)
  (0.025 0.062853761713077924)
  (0.05 2.2528791494664748)
  (0.1 4.7077966570929917)
  (0.25 7.7071094149528738)
  (0.5 10.078140841433592)
  (0.75 12.350609122495367)
  (0.84 13.90838744629569)
  (0.9 15.524798945657368)
  (0.95 17.919111158893685)
  (0.975 20.119800621312898)
  (0.99 23.384345610996593)
  (0.999 30.369129475837404)
  (0.9999 35.070809500650171)
  (0.99999 38.132242889203326)
  HPD intervals:
  HPD interval (0.84): 3.88538137399407..16.23145811408844

  var : p
  Probabilities:
  false: 0.8818000000000000
  true: 0.1182000000000000
  mean = [false = 0.8818,true = 0.1182]
  Percentiles:
  (0.001 false)
  (0.01 false)
  (0.025 false)
  (0.05 false)
  (0.1 false)
  (0.25 false)
  (0.5 false)
  (0.75 false)
  (0.84 false)
  (0.9 true)
  (0.95 true)
  (0.975 true)
  (0.99 true)
  (0.999 true)
  (0.9999 true)
  (0.99999 true)
  HPD intervals:
  show_hpd_intervals: data is not numeric

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  Mu = 10,
  Sigma = 3.4,
  D = laplace_dist(Mu,Sigma),

  P = check(abs(D) > 15),
  
  add("d",D),
  add("p",P).


/*
  Parameter recovery laplace(10,3.4)

  Num accepted samples: 100  Total samples: 280143  (0.000%)

  var : mu
  Probabilities (truncated):
  13.312547334149736: 0.0100000000000000
  12.84461641350976: 0.0100000000000000
  12.528418056912916: 0.0100000000000000
  12.353862688110146: 0.0100000000000000
  .........
  9.475632016302846: 0.0100000000000000
  9.470903179361907: 0.0100000000000000
  9.349480368825365: 0.0100000000000000
  8.583851069483837: 0.0100000000000000
  mean = 10.7733
  Percentiles:
  (0.001 8.6596483701186475)
  (0.01 9.3418240758319495)
  (0.025 9.4731493769088537)
  (0.05 9.4931194649511568)
  (0.1 9.6111427525110269)
  (0.25 10.261946409131376)
  (0.5 10.814504633105596)
  (0.75 11.260008339891213)
  (0.84 11.497713410992972)
  (0.9 11.684559495972731)
  (0.95 12.148961961338744)
  (0.975 12.445504256731599)
  (0.99 12.849295722716162)
  (0.999 13.266222173006376)
  (0.9999 13.3079148180354)
  (0.99999 13.312084082538302)
  HPD intervals:
  HPD interval (0.94): 9.34948036882536..12.14236030920519

  var : sigma
  Probabilities (truncated):
  6.497090759918601: 0.0100000000000000
  5.795355236993802: 0.0100000000000000
  5.521518553384356: 0.0100000000000000
  5.340142904473534: 0.0100000000000000
  .........
  2.719064914956253: 0.0100000000000000
  2.609460066356445: 0.0100000000000000
  2.602519943659436: 0.0100000000000000
  2.196091097870884: 0.0100000000000000
  mean = 3.79015
  Percentiles:
  (0.001 2.2363275536039504)
  (0.01 2.5984556552015503)
  (0.025 2.6615223694413541)
  (0.05 2.7819687960585435)
  (0.1 2.9229196714809724)
  (0.25 3.2023345298144661)
  (0.5 3.7787458457885057)
  (0.75 4.1580782221434998)
  (0.84 4.4663521619822601)
  (0.9 4.6311574404272982)
  (0.95 5.1362558110320276)
  (0.975 5.4353651201517144)
  (0.99 5.8023725922230538)
  (0.999 6.4276189431490431)
  (0.9999 6.4901435782416446)
  (0.99999 6.4963960417509039)
  HPD intervals:
  HPD interval (0.94): 2.60946006635645..5.13766806811917


  Compare with Mathematica:
  FindDistributionParameters[data, LaplaceDistribution[mu, sigma]] 
  -> {mu -> 10.1308, sigma -> 3.30752}

*/
go2 ?=>
  Data = laplace_dist_n(10,3.4,40),
  println(data=Data),
  reset_store,
  run_model(10_000,$model2(Data),[show_probs_trunc,mean,show_percentiles,
                                  show_hpd_intervals,hpd_intervals=[0.94],
                                  min_accepted_samples=100,show_accepted_samples=true
                                  ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2(Data) =>
  Mu = uniform(0,100),
  Sigma = uniform(0,20),
  
  X = laplace_dist_n(Mu,Sigma,Data.len),

  observe_abc(Data,X,1/10),

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