/* 

  Pareto (type I) distribution in Picat.

  https://en.wikipedia.org/wiki/Pareto_distribution
  Mathematica: ParetoDistribution[k, alpha]

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

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

main => go.

/*

  var : V1
  Probabilities (truncated):
  63818.943695483379997: 0.0001000000000000
  58370.201776736954344: 0.0001000000000000
  56831.090769014452235: 0.0001000000000000
  50765.364110406611871: 0.0001000000000000
  .........
  6820.510760699130515: 0.0001000000000000
  6820.449828006448115: 0.0001000000000000
  6820.401925808168926: 0.0001000000000000
  6820.270248104660823: 0.0001000000000000
  mean = 9065.54
  HPD intervals:
  HPD interval (0.84): 6820.27024810466082..10748.85514708295705
  Percentiles:
  (0.001 6821.6522447665375)
  (0.01 6840.7197854048291)
  (0.025 6865.1718472313441)
  (0.05 6908.1262835353928)
  (0.1 7000.9348004027279)
  (0.25 7328.3966802419536)
  (0.5 8112.5780611019673)
  (0.75 9603.513970077569)
  (0.84 10749.126251258307)
  (0.9 12031.25280855036)
  (0.95 14321.477964480857)
  (0.975 16856.72651381086)
  (0.99 21378.50136119435)
  (0.999 37658.774212631273)
  (0.9999 58370.746650924979)
  (0.99999 63274.123991034474)

  var : prob1
  Probabilities:
  true: 0.9002000000000000
  false: 0.0998000000000000
  mean = [true = 0.9002,false = 0.0998]
  HPD intervals:
  show_hpd_intervals: data is not numeric
  Percentiles:
  (0.001 false)
  (0.01 false)
  (0.025 false)
  (0.05 false)
  (0.1 true)
  (0.25 true)
  (0.5 true)
  (0.75 true)
  (0.84 true)
  (0.9 true)
  (0.95 true)
  (0.975 true)
  (0.99 true)
  (0.999 true)
  (0.9999 true)
  (0.99999 true)

   Theoretical for V1:
   Picat> X=pareto1_dist_mean(6820,4)  
   X = 9093.33333333333394


   Theoretical for Prob1:
   Picat> X=1-pareto1_dist_cdf(6820,4,7000)
   X = 0.901042629637651

   Quantiles:
   Picat> X=pareto1_dist_quantile(6820,4,0.99)   
   X = 21566.733642348342983

   Picat> X=pareto1_dist_quantile(6820,4,0.99999)
   X = 121278.65576479252195


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_hpd_intervals,hpd_intervals=[0.84],
                           show_percentiles]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  X = 6820,
  Alpha = 4,
  V1 = pareto1_dist(X,Alpha),
  Prob1 = check(V1 >= 7000),
  
  add("V1",V1),
  add("prob1",Prob1).