/* 

  Max stable distribution in Picat.

  From Mathematica MaxStableDistribution

  This is also called Generalized Extreme Value distribution 
  (for maximum values).

  Cf my Gamble model gamble_max_stable_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 : d
  Probabilities (truncated):
  68495.996019391619484: 0.0001000000000000
  8649.559211176370809: 0.0001000000000000
  7990.078222660015854: 0.0001000000000000
  4847.668887637149965: 0.0001000000000000
  .........
  0.106898646273663: 0.0001000000000000
  0.102210362951394: 0.0001000000000000
  0.10140909994393: 0.0001000000000000
  0.10105672708781: 0.0001000000000000
  mean = 20.2922

  var : p
  Probabilities:
  true: 0.6324000000000000
  false: 0.3676000000000000
  mean = [true = 0.6324,false = 0.3676]

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  Mu = 2,
  Sigma = 2,
  Xi = 0.9,

  D = max_stable_dist(Mu,Sigma,Xi),
  P = check(D >= Mu),
  add("d",D),
  add("p",P).


/*
  Parameter recovery for some different distributions
  100 generated values from 1000 max values

  * normal_dist(100,15) 

    data = [155.376,141.58,149.338,157.625,146.455,143.3,157.919,146.002,153.301,152.368,149.343,141.632,150.357,156.6,156.912,156.01,149.093,156.894,144.694,146.576,142.044,168.628,143.707,153.847,145.425,146.375,146.214,144.844,145.694,147.871,151.576,149.516,146.165,152.988,161.728,150.949,145.892,146.54,143.382,152.206,151.004,140.624,148.342,143.682,153.354,145.529,142.769,147.662,151.275,145.053,148.437,145.443,153.865,139.409,148.885,148.399,148.978,146.699,142.097,150.717,145.676,145.009,153.394,146.247,149.333,155.858,148.792,142.044,139.124,157.862,149.912,142.168,145.489,146.924,147.867,146.698,145.326,142.409,141.346,146.557,141.66,143.409,155.559,143.714,145.261,147.614,144.09,152.592,146.431,148.117,149.211,143.746,153.198,150.268,147.94,145.949,151.279,148.911,153.793,152.501]

    Num accepted samples: 100  Total samples: 55180  (0.002%)

    var : mu
    Probabilities (truncated):
    150.784666199971383: 0.0100000000000000
    150.773994376290375: 0.0100000000000000
    150.379236459050048: 0.0100000000000000
    150.370964068801754: 0.0100000000000000
    .........
    145.540390118918964: 0.0100000000000000
    145.530773790903197: 0.0100000000000000
    145.351951250202148: 0.0100000000000000
    145.201707758644943: 0.0100000000000000
    mean = 148.012
    HPD intervals:
    HPD interval (0.84): 145.74753196996940..149.77278852121441
    HPD interval (0.999): 145.20170775864494..150.78466619997138

    var : sigma
    Probabilities (truncated):
    7.681276596934198: 0.0100000000000000
    7.004477655051499: 0.0100000000000000
    6.899923052126505: 0.0100000000000000
    6.646660820882144: 0.0100000000000000
    .........
    0.926419059246042: 0.0100000000000000
    0.893320227457825: 0.0100000000000000
    0.566900614913041: 0.0100000000000000
    0.0504046026852: 0.0100000000000000
    mean = 4.23855
    HPD intervals:
    HPD interval (0.84): 2.39016848727603..6.46665619987373
    HPD interval (0.999): 0.05040460268520..7.68127659693420

    var : xi
    Probabilities (truncated):
    1.359120382163264: 0.0100000000000000
    0.831995380032805: 0.0100000000000000
    0.424410194821847: 0.0100000000000000
    0.418644046140203: 0.0100000000000000
    .........
    -2.715636033432389: 0.0100000000000000
    -2.780946736121991: 0.0100000000000000
    -2.88373482780705: 0.0100000000000000
    -2.89111817250546: 0.0100000000000000
    mean = -0.973859
    HPD intervals:
    HPD interval (0.84): -1.92241579150894..0.39664368768998
    HPD interval (0.999): -2.89111817250546..1.35912038216326

    var : post
    Probabilities (truncated):
    160.242157150257384: 0.0100000000000000
    157.144434783436594: 0.0100000000000000
    156.419535072074893: 0.0100000000000000
    155.329122256481838: 0.0100000000000000
    .........
    139.506715394343331: 0.0100000000000000
    138.28056131617771: 0.0100000000000000
    135.535272660382105: 0.0100000000000000
    113.13583077870328: 0.0100000000000000
    mean = 148.812
    HPD intervals:
    HPD interval (0.84): 144.88540400944333..154.39956138475594
    HPD interval (0.999): 113.13583077870328..160.24215715025738

    Mathematica for this data:
    FindDistributionParameters[data, MaxStableDistribution[mu, sigma, xi]]
    -> {xi -> -0.042961, sigma -> 4.21246, mu -> 146.162}

  * binomial_dist_n(100,0.5,1000)
    data = [64,67,67,68,64,68,64,66,66,64,63,66,67,65,65,64,65,65,65,66,66,64,67,66,66,67,67,69,69,67,68,66,66,63,67,66,65,62,67,63,64,67,65,66,66,67,65,66,66,68,68,67,66,64,63,67,65,64,63,63,64,65,65,64,66,65,68,66,65,63,65,67,67,64,67,64,67,63,68,66,64,66,64,66,67,65,65,68,68,66,65,65,65,71,66,66,68,64,67,63]

    min_accepted_samples = 100

    var : mu
    Probabilities (truncated):
    66.317567416342541: 0.0100000000000000
    66.308501635343262: 0.0100000000000000
    66.281455526318126: 0.0100000000000000
    66.237092139141865: 0.0100000000000000
    .........
    64.753265316406399: 0.0100000000000000
    64.667491567634585: 0.0100000000000000
    64.593703391101243: 0.0100000000000000
    64.543224229188709: 0.0100000000000000
    mean = 65.526
    HPD intervals:
    HPD interval (0.84): 64.95936583734380..66.20738462696869
    HPD interval (0.999): 64.54322422918871..66.31756741634254

    var : sigma
    Probabilities (truncated):
    2.410377069567506: 0.0100000000000000
    2.243373590634844: 0.0100000000000000
    2.232286782112106: 0.0100000000000000
    2.202022598219115: 0.0100000000000000
    .........
    0.273775095247559: 0.0100000000000000
    0.130604701177499: 0.0100000000000000
    0.105830282022166: 0.0100000000000000
    0.029393052695968: 0.0100000000000000
    mean = 1.35437
    HPD intervals:
    HPD interval (0.84): 0.52915893519724..2.03298613523738
    HPD interval (0.999): 0.02939305269597..2.41037706956751

    var : xi
    Probabilities (truncated):
    1.499923779861966: 0.0100000000000000
    1.139363806759643: 0.0100000000000000
    0.958473644200002: 0.0100000000000000
    0.904045744754395: 0.0100000000000000
    .........
    -2.492045142451322: 0.0100000000000000
    -2.496890159089533: 0.0100000000000000
    -2.529429760542433: 0.0100000000000000
    -2.938032444537632: 0.0100000000000000
    mean = -0.908073
    HPD intervals:
    HPD interval (0.84): -1.90498641454847..0.41775143491931
    HPD interval (0.999): -2.93803244453763..1.49992377986197

    var : post
    Probabilities (truncated):
    67.237483615938373: 0.0100000000000000
    67.092665059401114: 0.0100000000000000
    67.063740201181076: 0.0100000000000000
    67.035867033032233: 0.0100000000000000
    .........
    61.847206913141285: 0.0100000000000000
    58.161567934765657: 0.0100000000000000
    55.387663242298032: 0.0100000000000000
    54.866370788435646: 0.0100000000000000
    mean = 65.3422
    HPD intervals:
    HPD interval (0.84): 64.15271554148748..67.06374020118108
    HPD interval (0.999): 54.86637078843565..67.23748361593837


    Mathematica for this data:
    FindDistributionParameters[data, MaxStableDistribution[mu, sigma, xi]]
    -> {xi -> -0.192996, sigma -> 1.57653, mu -> 65.0293}

  * exponential_dist_n(1/3,1000)

    data = [20.3686,29.9039,23.4305,21.3033,23.3937,22.1753,22.671,28.2861,21.1869,20.6018,19.5325,17.7353,21.0796,22.3865,20.0765,19.1643,26.9191,19.622,20.5386,30.148,17.8932,25.3397,20.829,23.6383,24.5573,18.7436,22.7233,31.6542,21.808,24.4283,18.9566,22.235,19.353,20.4063,21.0214,39.7966,25.3183,18.9478,17.9771,20.4333,23.268,23.5584,31.4024,18.9491,17.3201,21.7407,22.092,27.4217,17.9928,25.9311,18.319,23.5066,19.9457,23.685,33.3652,26.9951,19.1964,23.908,19.468,19.3625,18.1883,16.9218,23.8251,18.1796,16.8194,23.1363,25.8771,17.3123,21.6287,33.0263,23.0487,20.6824,19.3702,23.187,20.2264,20.5285,26.6023,21.8571,20.8598,18.2995,19.3248,21.1916,21.2578,24.9312,18.2481,20.5757,19.7764,19.8163,18.4614,19.7039,20.4325,22.4178,27.8575,20.5657,22.4555,25.5838,20.021,19.8057,24.3259,24.2055]
    [min = 16.8194,mean = 22.2455,median = 21.1893,max = 39.7966,variance = 15.8185,stdev = 3.97725]

    min_accepted_samples = 100

    var : mu
    Probabilities (truncated):
    23.934169924046635: 0.0100000000000000
    23.933831693327743: 0.0100000000000000
    23.884821824687407: 0.0100000000000000
    23.85344556554266: 0.0100000000000000
    .........
    19.90787769191882: 0.0100000000000000
    19.812185474535884: 0.0100000000000000
    19.790986194983269: 0.0100000000000000
    19.784245723612454: 0.0100000000000000
    mean = 21.987
    HPD intervals:
    HPD interval (0.84): 20.45590111865587..23.54606459201455
    HPD interval (0.999): 19.78424572361245..23.93416992404664

    var : sigma
    Probabilities (truncated):
    6.181442433121354: 0.0100000000000000
    5.547073541929514: 0.0100000000000000
    5.519084076173177: 0.0100000000000000
    5.478874363693816: 0.0100000000000000
    .........
    1.04671970989868: 0.0100000000000000
    0.988802882371844: 0.0100000000000000
    0.962240444944352: 0.0100000000000000
    0.266359718640502: 0.0100000000000000
    mean = 3.39536
    HPD intervals:
    HPD interval (0.84): 1.48543346742421..4.85949205460934
    HPD interval (0.999): 0.26635971864050..6.18144243312135

    var : xi
    Probabilities (truncated):
    0.859240144425649: 0.0100000000000000
    0.715203142592312: 0.0100000000000000
    0.534213931082848: 0.0100000000000000
    0.332519015917796: 0.0100000000000000
    .........
    -2.654742934580307: 0.0100000000000000
    -2.666559803144336: 0.0100000000000000
    -2.735627802897072: 0.0100000000000000
    -2.754427741167335: 0.0100000000000000
    mean = -1.05462
    HPD intervals:
    HPD interval (0.84): -2.28822734266903..0.13580304623386
    HPD interval (0.999): -2.75442774116733..0.85924014442565

    var : post
    Probabilities (truncated):
    35.00884864884739: 0.0100000000000000
    31.797010459007382: 0.0100000000000000
    29.69882002545512: 0.0100000000000000
    28.231046689603794: 0.0100000000000000
    .........
    13.296175427706954: 0.0100000000000000
    13.162445264265413: 0.0100000000000000
    12.607874521195278: 0.0100000000000000
    11.944655343571814: 0.0100000000000000
    mean = 22.1397
    HPD intervals:
    HPD interval (0.84): 16.85337074917767..26.49446598421405
    HPD interval (0.999): 11.94465534357181..35.00884864884739

    [min = 11.9447,mean = 22.1397,median = 22.9527,max = 35.0088,variance = 15.4815,stdev = 3.93466]

    Mathematica for this data:
    FindDistributionParameters[data, MaxStableDistribution[mu, sigma, xi]]
    -> {xi -> 0.154302, sigma -> 2.55254, mu -> 20.3325}

*/
go2 ?=>
  Data=[normal_dist_n(100,15,1000).max : _ in 1..100],
  % Data=[binomial_dist_n(100,0.5,1000).max : _ in 1..100],
  % Data=[exponential_dist_n(1/3,1000).max : _ in 1..100],
  % Data=[gamma_dist_n(10,10,1000).max : _ in 1..4],
  println(data=Data),
  show_simple_stats(Data),
  reset_store,
  run_model(1_000_000,$model2(Data),[show_probs_trunc,mean,
                                     show_hpd_intervals,hpd_intervals=[0.84,0.999],
                                     min_accepted_samples=100,show_accepted_samples=false % true
                                     ]),
  Post = get_store().get("post"),
  show_simple_stats(Post),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go2 => true.
  
  
model2(Data) =>
  Mu    = normal_dist(Data.avg,Data.stdev),
  Sigma = uniform(0,20),
  Xi    = uniform(-3,3),

  Y = max_stable_dist_n(Mu,Sigma,Xi,Data.len),
  observe_abc(Data,Y,1/10),
  
  Post = max_stable_dist(Mu,Sigma,Xi),
  % P = check(Post >= Data.max),
  
  if observed_ok then
    add("mu",Mu),
    add("sigma",Sigma),
    add("xi",Xi),
    add("post",Post),
    % add("p",P),    
  end.