/* 

  Frechet distribution in Picat.

  From Mathematica (FrechetDistribution)
  """
  Quantile(FrechetDistribution(alpha, beta), x)
  -> 
  beta*(-Log(x))^(-1/alpha)    0 < x < 1
  0                             x <= 0
  Infinity                      True
  """
      
  Example (op.cit)
  """
  According to a study, the annual maximal tephra (solid material) volume in volcanic eruptions 
  follows a FrechetDistribution with shape parameter 0.71 and scale parameter 6.3, given in cubic 
  kilometers:
  D = FrechetDistribution(0.71,6.3)
  ...
  Find the probability that the annual maximal tephra volume is greater than 30 cubic kilometers:
  Probability(x > 30,  Element(x, D))
  -> 0.281219
  ""

  Picat> X=1-frechet_dist_cdf(0.71,6.3,30)
  X = 0.281219288418227

  Picat> X=frechet_dist_quantile(0.71,6.3,0.5)
  X = 10.55674523617512


  For more on the Frechet distribution, see
  ppl_distributions.pi and ppl_distributions_test.pi

  Cf my Gamble model gamble_frechet_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):
  11.175199720963459: 0.0002000000000000
  1132496.705504813697189: 0.0001000000000000
  829856.47760959668085: 0.0001000000000000
  410755.327565882180352: 0.0001000000000000
  .........
  0.364476064967467: 0.0001000000000000
  0.352702131206569: 0.0001000000000000
  0.340943013909871: 0.0001000000000000
  0.340850951682057: 0.0001000000000000
  mean = 630.602
  HPD intervals:
  HPD interval (0.84): 0.36447606496747..71.06928968297143

  var : p
  Probabilities:
  false: 0.7272000000000000
  true: 0.2728000000000000
  mean = [false = 0.7272,true = 0.2728]
  HPD intervals:
  show_hpd_intervals: data is not numeric

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,
                           show_hpd_intervals,hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  Alpha = 0.71,
  Beta = 6.3,
  D = frechet_dist(Alpha,Beta),
  P = check(D > 30),
  
  add("d",D),
  add("p",P).