/* 

  Extreme value distribution test in Picat.

  From Mathematica ExtremeValueDistribution
  """
  ExtremeValueDistribution can be used to model monthly maximum wind speeds. Recorded monthly 
  maxima of wind speeds in km/h for Boston, MA, from January 1950 till December 2009:
   ...
  Fit the distribution to the data:
  edist = EstimatedDistribution(maxWinds, ExtremeValueDistribution(a, b))
  -> ExtremeValueDistribution(47.2374, 9.1497))

  ... 

  Find the probability of monthly maximum wind speed exceeding 60 mph:
  Probability(x > Quantity(60, "mph"), x ~ edist)
  -> 0.00454842
  """

  Note that the original data was in km/h so let's keep that unit:
  """
  Probability(x > Quantity(60, "km per hour"), x ~ edist)
  -> 0.219536
  """

  Cf my Gamble model gamble_extreme_value_test.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 : d
  Probabilities (truncated):
  126.016038095973357: 0.0001000000000000
  123.920251128210651: 0.0001000000000000
  121.178751253953095: 0.0001000000000000
  120.819548297824554: 0.0001000000000000
  .........
  26.356197938898255: 0.0001000000000000
  25.961885175907025: 0.0001000000000000
  25.952295216077886: 0.0001000000000000
  25.884473976707714: 0.0001000000000000
  mean = 52.5041
  Percentiles:
  (0.001 29.56036182716468)
  (0.01 33.535675866798435)
  (0.025 35.165396951693033)
  (0.05 36.988344782623017)
  (0.1 39.457953655825264)
  (0.25 44.163461491879531)
  (0.5 50.607774444859409)
  (0.75 58.692087788330745)
  (0.84 63.311569119045345)
  (0.9 67.619545850838136)
  (0.95 74.178005834416211)
  (0.975 81.433893957874247)
  (0.99 89.353354365385144)
  (0.999 111.25421586959706)
  (0.9999 123.92046070690594)
  (0.99999 125.80648035706929)

  var : p
  Probabilities:
  false: 0.7776000000000000
  true: 0.2224000000000000
  mean = [false = 0.7776,true = 0.2224]
  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 true)
  (0.9 true)
  (0.95 true)
  (0.975 true)
  (0.99 true)
  (0.999 true)
  (0.9999 true)
  (0.99999 true)



*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_percentiles]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  A = 47.2374,
  B = 9.1497,
  D = extreme_value_dist(A,B),
  P = check(D > 60),
  
  add("d",D),
  add("p",P).