/* Min stable distribution in Picat. From Mathematica's MinStableDistribution This is also called Generalized Extreme Value distribution (for minimum values). Cf - ppl_max_stable_dist.pi - my Gamble model gamble_min_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. % import ordset. main => go. /* var : d Probabilities (truncated): 2.964564976566995: 0.0001000000000000 2.963386214769604: 0.0001000000000000 2.951024855077758: 0.0001000000000000 2.948487277897518: 0.0001000000000000 ......... -2254.440318373098762: 0.0001000000000000 -2610.07741754912513: 0.0001000000000000 -2854.791577610351851: 0.0001000000000000 -3691.067519791808536: 0.0001000000000000 mean = -3.9555 var : p Probabilities: true: 0.6344000000000000 false: 0.3656000000000000 mean = 0.6344 */ go ?=> reset_store, run_model(10_000,$model,[show_probs_trunc,mean % , % show_percentiles, % show_hpd_intervals,hpd_intervals=[0.94], % show_histogram, % min_accepted_samples=1000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go => true. model() => Mu = 2, Sigma = 1, Xi = 0.91, D = min_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) ata = [52.2457,50.3331,50.5745,59.5358,44.6794,45.9373,43.8911,47.8911,46.6448,52.1717,57.6365,51.5661,50.2668,55.1991,61.4971,55.142,52.6328,36.1755,45.3101,58.4281,59.77,51.1347,49.5661,54.5427,55.7496,48.0242,54.704,52.4439,55.6391,57.5156,49.8339,56.578,55.0912,55.0812,57.7574,54.9812,52.288,59.4496,54.5976,58.1715,50.0853,40.4693,47.0353,51.3869,50.9317,53.0083,48.5121,50.8732,48.9543,57.1818,36.5173,56.3951,48.6537,44.6968,45.1373,50.4733,43.0206,55.0862,48.1573,48.3745,46.9145,58.2011,58.5196,51.7453,50.7022,57.3354,42.6152,57.7422,53.8274,58.2036,48.3229,56.2557,39.6128,47.8297,50.9057,46.2272,47.7191,47.1141,53.2916,54.0832,46.1916,47.5187,55.6883,50.3734,34.8906,58.9863,45.457,51.7422,44.4938,43.6683,49.1532,54.3772,55.9552,45.5637,53.7519,53.16,44.8064,53.1535,57.3142,49.8931] [min = 34.8906,mean = 51.0894,median = 51.2608,max = 61.4971,variance = 30.3505,stdev = 5.50913] min_accepted_samples = 100 var : mu Probabilities (truncated): 54.456561899874472: 0.0100000000000000 54.276167088468227: 0.0100000000000000 54.143891812826894: 0.0100000000000000 53.928154889768798: 0.0100000000000000 ......... 49.008150087937466: 0.0100000000000000 48.950085618231803: 0.0100000000000000 48.895011255474394: 0.0100000000000000 48.745750047497737: 0.0100000000000000 mean = 51.4034 HPD intervals: HPD interval (0.84): 49.79626369440624..53.53525390001045 HPD interval (0.999): 48.74575004749774..54.45656189987447 var : sigma Probabilities (truncated): 7.491220304505537: 0.0100000000000000 6.969778270912253: 0.0100000000000000 6.91008071736902: 0.0100000000000000 6.721770095975032: 0.0100000000000000 ......... 0.974777732498375: 0.0100000000000000 0.767940376311513: 0.0100000000000000 0.260928114997655: 0.0100000000000000 0.086371973197149: 0.0100000000000000 mean = 4.5685 HPD intervals: HPD interval (0.84): 2.28717372859231..6.72177009597503 HPD interval (0.999): 0.08637197319715..7.49122030450554 var : xi Probabilities (truncated): 1.477427698894138: 0.0100000000000000 1.378552980897274: 0.0100000000000000 0.658955652107929: 0.0100000000000000 0.633978019298044: 0.0100000000000000 ......... -2.339606864070337: 0.0100000000000000 -2.641911255028989: 0.0100000000000000 -2.666657273502395: 0.0100000000000000 -2.694217792569761: 0.0100000000000000 mean = -0.929357 HPD intervals: HPD interval (0.84): -2.26476547087765..0.09826078968973 HPD interval (0.999): -2.69421779256976..1.47742769889414 var : post Probabilities (truncated): 63.924292080061456: 0.0100000000000000 61.246901587544855: 0.0100000000000000 61.089334533148062: 0.0100000000000000 59.439670440778769: 0.0100000000000000 ......... 40.162805430743617: 0.0100000000000000 35.256090407433312: 0.0100000000000000 35.0295459695255: 0.0100000000000000 9.897029195662192: 0.0100000000000000 mean = 51.4586 HPD intervals: HPD interval (0.84): 45.22638101451440..58.29452507531460 HPD interval (0.999): 9.89702919566219..63.92429208006146 [min = 9.89703,mean = 51.4586,median = 52.3555,max = 63.9243,variance = 43.3449,stdev = 6.58369] Mathematica for this data: FindDistributionParameters[data, MinStableDistribution[mu, sigma, xi]] FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. -> {mu -> 53.4555, xi -> -0.0964524, sigma -> 4.78664} * binomial_dist_n(100,0.5,1000) data = [34,35,35,34,35,35,33,33,35,35,35,35,34,35,34,35,35,37,37,35,31,35,35,36,35,34,32,36,35,36,36,33,35,32,34,35,33,34,33,32,35,33,35,31,35,30,36,33,34,35,34,34,35,32,34,33,33,35,33,34,34,34,34,36,33,33,34,35,29,34,35,36,33,33,33,35,34,34,33,33,33,36,33,34,31,37,36,34,36,36,35,34,36,35,35,36,32,34,35,33] [min = 29,mean = 34.18,median = 34.0,max = 37,variance = 2.1476,stdev = 1.46547] min_accepted_samples = 100 var : mu Probabilities (truncated): 35.084634496567425: 0.0100000000000000 35.023618047128473: 0.0100000000000000 35.018805281424576: 0.0100000000000000 35.000568939438907: 0.0100000000000000 ......... 33.809208410876707: 0.0100000000000000 33.752251019176349: 0.0100000000000000 33.700641268367164: 0.0100000000000000 33.627171429831712: 0.0100000000000000 mean = 34.3533 HPD intervals: HPD interval (0.84): 33.88110631996670..34.86884194372266 HPD interval (0.999): 33.62717142983171..35.08463449656742 var : sigma Probabilities (truncated): 1.949369982792702: 0.0100000000000000 1.917647645770408: 0.0100000000000000 1.875646823027938: 0.0100000000000000 1.845168975109779: 0.0100000000000000 ......... 0.290682027251777: 0.0100000000000000 0.263271117705512: 0.0100000000000000 0.029178047566292: 0.0100000000000000 0.018193461009391: 0.0100000000000000 mean = 1.22581 HPD intervals: HPD interval (0.84): 0.60270564658693..1.72903229562986 HPD interval (0.999): 0.01819346100939..1.94936998279270 var : xi Probabilities (truncated): 2.079901523459657: 0.0100000000000000 1.249220471945228: 0.0100000000000000 0.745297181301423: 0.0100000000000000 0.665374941036745: 0.0100000000000000 ......... -2.260395948430708: 0.0100000000000000 -2.390236730403377: 0.0100000000000000 -2.576896810707122: 0.0100000000000000 -2.844944276774742: 0.0100000000000000 mean = -0.814822 HPD intervals: HPD interval (0.84): -2.12365121632984..0.16374155840079 HPD interval (0.999): -2.84494427677474..2.07990152345966 var : post Probabilities (truncated): 39.724178858574973: 0.0100000000000000 37.28428442908308: 0.0100000000000000 37.247441680001479: 0.0100000000000000 36.882553534525485: 0.0100000000000000 ......... 31.090352533366353: 0.0100000000000000 28.439170769379505: 0.0100000000000000 26.652078830234778: 0.0100000000000000 21.501609618654594: 0.0100000000000000 mean = 34.4104 HPD intervals: HPD interval (0.84): 33.20372877040590..36.39220313845838 HPD interval (0.999): 21.50160961865459..39.72417885857497 [min = 21.5016,mean = 34.4104,median = 34.5187,max = 39.7242,variance = 4.21212,stdev = 2.05234] Mathematica for this data: FindDistributionParameters[data, MinStableDistribution[mu, sigma, xi]] -> {mu -> 34.8038, xi -> -0.0918462, sigma -> 1.27101} */ go2 ?=> Data=[normal_dist_n(100,15,1000).min : _ in 1..100], % Data=[binomial_dist_n(100,0.5,1000).min : _ in 1..100], % Data=[exponential_dist_n(1/3,1000).min : _ in 1..100], % Data=[gamma_dist_n(10,10,1000).min : _ 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 = min_stable_dist_n(Mu,Sigma,Xi,Data.len), observe_abc(Data,Y,1/10), Post = max_stable_dist(Mu,Sigma,Xi), if observed_ok then add("mu",Mu), add("sigma",Sigma), add("xi",Xi), add("post",Post), end.