/* Generalized extreme value dist in Picat. This is from Wikipedia https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution GEV(mu, sigma, xi) See ppl_distributions.pi and ppl_distributions_test.pi for more on this. Cf - ppl_min_stable_dist.pi - ppl_max_stable_dist.pi - my Gamble model gamble_generalized_extreme_value_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. /* PDF: 2,2,0,1 = 0.158521 2,2,0,-1 = 0.0253536 2,2,2,1 = 0.270671 2,2,1,-1 = 0 2,2,-1,1 = 0.111565 2,2,-1,-1 = 0.0410425 CDF: 2,2,0,1 = 0.192296 2,2,0,-1 = 0.0113143 2,2,1,1 = 0.135335 2,2,1,-1 = 0 2,2-1,1 = 0.22313 2,2,-1,-1 = 0.082085 10,5,10,100 = 0.551778 Quantile: 2,2,0,0.01 = -1.05436 2,2,0,0.99 = 11.2003 2,2,0,0.9999 = 20.4206 2,2,1,0.01 = 0.434294 2,2,1,0,99 = 198.998 2,2,1,0,9999 = 19999.0 2,2,-1,0.1 = -5.21034 2,2-1,0.99 = 3.9799 2,2,-1,0.9999 = 3.9998 Mean: 2,2,0 theoretical = 3.15443 sample_mean = 3.16804 2,2,0.99 theoretical = 200.854 sample_mean = 35.503 2,2,1 theoretical = 1.0e+100 sample_mean = 24.5846 2,2,-1 theoretical = 2.0 sample_mean = 1.99658 Histogram for generalized_extreme_value_dist(10,3,0.1) Histogram (total 10000) 4.085: 14 (0.001 / 0.001) 5.519: 201 ####### (0.020 / 0.021) 6.952: 850 ############################ (0.085 / 0.107) 8.385: 1484 ################################################ (0.148 / 0.255) 9.818: 1839 ############################################################ (0.184 / 0.439) 11.251: 1565 ################################################### (0.157 / 0.595) 12.684: 1168 ###################################### (0.117 / 0.712) 14.118: 874 ############################# (0.087 / 0.800) 15.551: 608 #################### (0.061 / 0.860) 16.984: 435 ############## (0.043 / 0.904) 18.417: 292 ########## (0.029 / 0.933) 19.850: 191 ###### (0.019 / 0.952) 21.283: 135 #### (0.013 / 0.966) 22.717: 84 ### (0.008 / 0.974) 24.150: 67 ## (0.007 / 0.981) 25.583: 51 ## (0.005 / 0.986) 27.016: 43 # (0.004 / 0.990) 28.449: 21 # (0.002 / 0.992) 29.882: 20 # (0.002 / 0.994) 31.316: 12 (0.001 / 0.995) 32.749: 10 (0.001 / 0.996) 34.182: 11 (0.001 / 0.997) 35.615: 5 (0.001 / 0.998) 37.048: 4 (0.000 / 0.998) 38.482: 5 (0.001 / 0.999) 39.915: 3 (0.000 / 0.999) 41.348: 4 (0.000 / 1.000) 42.781: 0 (0.000 / 1.000) 44.214: 1 (0.000 / 1.000) 45.647: 1 (0.000 / 1.000) 47.081: 0 (0.000 / 1.000) 48.514: 0 (0.000 / 1.000) 49.947: 0 (0.000 / 1.000) 51.380: 0 (0.000 / 1.000) 52.813: 0 (0.000 / 1.000) 54.246: 0 (0.000 / 1.000) 55.680: 1 (0.000 / 1.000) 57.113: 0 (0.000 / 1.000) 58.546: 0 (0.000 / 1.000) 59.979: 1 (0.000 / 1.000) */ go ?=> _ = random2(), println("PDF:"), println("2,2,0,1"=generalized_extreme_value_dist_pdf(2,2,0,1)), % 0.158521 println("2,2,0,-1"=generalized_extreme_value_dist_pdf(2,2,0,-1)), % 0.02535356 ok println("2,2,2,1"=generalized_extreme_value_dist_pdf(2,2,1, 1)), % 0.2706706 ok println("2,2,1,-1"=generalized_extreme_value_dist_pdf(2,2,1,-1)), % 0 ok <-- println("2,2,-1,1"=generalized_extreme_value_dist_pdf(2,2,-1,1)), % 0.1115651 ok println("2,2,-1,-1"=generalized_extreme_value_dist_pdf(2,2,-1,-1)), % 0.0410425 ok nl, println("CDF:"), println("2,2,0,1"=generalized_extreme_value_dist_cdf(2,2,0,1)), % 0.1922956 ok println("2,2,0,-1"=generalized_extreme_value_dist_cdf(2,2,0,-1)), % 0.01131429 ok println("2,2,1,1"=generalized_extreme_value_dist_cdf(2,2,1,1)), % 0.1353353 ok println("2,2,1,-1"=generalized_extreme_value_dist_cdf(2,2,1,-1)), % 0 ok println("2,2-1,1"=generalized_extreme_value_dist_cdf(2,2,-1,1)), % 0.2231302 ok println("2,2,-1,-1"=generalized_extreme_value_dist_cdf(2,2,-1,-1)), % 0.082085 ok println("10,5,10,100"=generalized_extreme_value_dist_cdf(10,5,10,100)), % 0.551778 nl, println("Quantile:"), println("2,2,0,0.01"=generalized_extreme_value_dist_quantile(2,2,0,0.01)), % -1.054359 ok println("2,2,0,0.99"=generalized_extreme_value_dist_quantile(2,2,0,0.99)), % 11.2003 ok println("2,2,0,0.9999"=generalized_extreme_value_dist_quantile(2,2,0,0.9999)), println("2,2,1,0.01"=generalized_extreme_value_dist_quantile(2,2,1,0.01)), % 0.4342945 ok println("2,2,1,0,99"=generalized_extreme_value_dist_quantile(2,2,1,0.99)), % 198.9983 ol println("2,2,1,0,9999"=generalized_extreme_value_dist_quantile(2,2,1,0.9999)), % 198.9983 ok println("2,2,-1,0.1"=generalized_extreme_value_dist_quantile(2,2,-1,0.01)), % -5.21034 ok println("2,2-1,0.99"=generalized_extreme_value_dist_quantile(2,2,-1,0.99)), % 3.979899 ok println("2,2,-1,0.9999"=generalized_extreme_value_dist_quantile(2,2,-1,0.9999)), nl, println("Mean:"), println("2,2,0"), println(theoretical=generalized_extreme_value_dist_mean(2,2,0)), println(sample_mean=generalized_extreme_value_dist_n(2,2,0,100000).mean), nl, println("2,2,0.99"), println(theoretical=generalized_extreme_value_dist_mean(2,2,0.99)), println(sample_mean=generalized_extreme_value_dist_n(2,2,0.99,100000).mean), nl, println("2,2,1"), println(theoretical=generalized_extreme_value_dist_mean(2,2,1)), println(sample_mean=generalized_extreme_value_dist_n(2,2,1,100000).mean), nl, println("2,2,-1"), println(theoretical=generalized_extreme_value_dist_mean(2,2,-1)), println(sample_mean=generalized_extreme_value_dist_n(2,2,-1,100000).mean), nl, println("Histogram for generalized_extreme_value_dist(10,3,0.1)"), generalized_extreme_value_dist_n(10,3,0.1,10000).show_histogram, nl. go => true. /* var : d Probabilities (truncated): 62.291933155257404: 0.0001000000000000 52.896883278542767: 0.0001000000000000 49.732144056362429: 0.0001000000000000 44.274454670294674: 0.0001000000000000 43.581166391819963: 0.0001000000000000 42.60749537657918: 0.0001000000000000 42.158608300659992: 0.0001000000000000 41.947062496254119: 0.0001000000000000 41.331034191216162: 0.0001000000000000 39.247543940185224: 0.0001000000000000 ......... 4.617127702701773: 0.0001000000000000 4.605055169940218: 0.0001000000000000 4.572215470731838: 0.0001000000000000 4.553332768713202: 0.0001000000000000 4.531223710605813: 0.0001000000000000 4.477480474153186: 0.0001000000000000 4.23957511288611: 0.0001000000000000 4.166666575381338: 0.0001000000000000 4.066680344166731: 0.0001000000000000 3.748486969085178: 0.0001000000000000 mean = 12.0426 HPD intervals: HPD interval (0.84): 6.09894397482339..16.08065565361940 HPD interval (0.9): 5.94364378434305..18.03195174343821 HPD interval (0.99): 4.53122371060581..27.21809085236705 HPD interval (0.999): 3.74848696908518..38.74092541550811 HPD interval (0.9999): 3.74848696908518..52.89688327854277 HPD interval (0.99999): 3.74848696908518..62.29193315525740 */ go2 ?=> reset_store, run_model(10_000,$model2,[show_probs_trunc,mean, % show_percentiles, truncate_size=10, show_hpd_intervals,hpd_intervals=[0.84,0.9,0.99,0.999,0.9999,0.99999] % show_histogram, % min_accepted_samples=1000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go2 => true. model2() => Mu = 10, Sigma = 3, Xi = 0.1, D = generalized_extreme_value_dist(Mu,Sigma,Xi), add("d",D). /* Recover parameters 100 max values from 1000 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] [min = 139.124,mean = 148.425,median = 147.765,max = 168.628,variance = 26.1524,stdev = 5.11394] Histogram (total 100) 139.124: 2 ########### (0.020 / 0.020) 139.862: 0 (0.000 / 0.020) 140.599: 1 ##### (0.010 / 0.030) 141.337: 4 ###################### (0.040 / 0.070) 142.074: 5 ########################### (0.050 / 0.120) 142.812: 1 ##### (0.010 / 0.130) 143.550: 7 ###################################### (0.070 / 0.200) 144.287: 1 ##### (0.010 / 0.210) 145.025: 6 ################################# (0.060 / 0.270) 145.762: 9 ################################################# (0.090 / 0.360) 146.500: 11 ############################################################ (0.110 / 0.470) 147.238: 1 ##### (0.010 / 0.480) 147.975: 7 ###################################### (0.070 / 0.550) 148.713: 6 ################################# (0.060 / 0.610) 149.450: 6 ################################# (0.060 / 0.670) 150.188: 3 ################ (0.030 / 0.700) 150.926: 5 ########################### (0.050 / 0.750) 151.663: 1 ##### (0.010 / 0.760) 152.401: 4 ###################### (0.040 / 0.800) 153.139: 5 ########################### (0.050 / 0.850) 153.876: 3 ################ (0.030 / 0.880) 154.614: 0 (0.000 / 0.880) 155.351: 2 ########### (0.020 / 0.900) 156.089: 2 ########### (0.020 / 0.920) 156.827: 3 ################ (0.030 / 0.950) 157.564: 3 ################ (0.030 / 0.980) 158.302: 0 (0.000 / 0.980) 159.039: 0 (0.000 / 0.980) 159.777: 0 (0.000 / 0.980) 160.515: 0 (0.000 / 0.980) 161.252: 0 (0.000 / 0.980) 161.990: 1 ##### (0.010 / 0.990) 162.727: 0 (0.000 / 0.990) 163.465: 0 (0.000 / 0.990) 164.203: 0 (0.000 / 0.990) 164.940: 0 (0.000 / 0.990) 165.678: 0 (0.000 / 0.990) 166.416: 0 (0.000 / 0.990) 167.153: 0 (0.000 / 0.990) 167.891: 1 ##### (0.010 / 1.000) min_accepted_samples = 100 var : mu Probabilities (truncated): 149.371044734302245: 0.0100000000000000 148.692060507541129: 0.0100000000000000 148.67775304384034: 0.0100000000000000 148.492446461793776: 0.0100000000000000 148.40707968408546: 0.0100000000000000 148.396961502609713: 0.0100000000000000 148.229231619752682: 0.0100000000000000 148.19554383726171: 0.0100000000000000 148.109953202416136: 0.0100000000000000 148.074019991893664: 0.0100000000000000 ......... 145.99111392105857: 0.0100000000000000 145.890338223781242: 0.0100000000000000 145.873749521823299: 0.0100000000000000 145.847092983688867: 0.0100000000000000 145.798854914522821: 0.0100000000000000 145.703350275977641: 0.0100000000000000 145.674388619155764: 0.0100000000000000 145.627993579889562: 0.0100000000000000 145.559039782028549: 0.0100000000000000 145.478742746649345: 0.0100000000000000 mean = 147.007 HPD intervals: HPD interval (0.84): 145.87374952182330..148.07401999189366 HPD interval (0.9): 145.62799357988956..148.10995320241614 HPD interval (0.99): 145.47874274664935..148.69206050754113 HPD interval (0.999): 145.47874274664935..149.37104473430225 HPD interval (0.9999): 145.47874274664935..149.37104473430225 HPD interval (0.99999): 145.47874274664935..149.37104473430225 var : sigma Probabilities (truncated): 6.860010468801488: 0.0100000000000000 6.717581286429233: 0.0100000000000000 6.634329616387529: 0.0100000000000000 6.605968646056004: 0.0100000000000000 6.547066116028962: 0.0100000000000000 6.446615469849956: 0.0100000000000000 6.416063423462242: 0.0100000000000000 6.383391733459845: 0.0100000000000000 6.37351978867013: 0.0100000000000000 6.345198567186109: 0.0100000000000000 ......... 2.560521523729209: 0.0100000000000000 2.253592806986343: 0.0100000000000000 2.218740695258016: 0.0100000000000000 2.111797464132215: 0.0100000000000000 1.768423859853495: 0.0100000000000000 1.725724787323608: 0.0100000000000000 1.422583778119918: 0.0100000000000000 1.21516509038171: 0.0100000000000000 0.973026818117605: 0.0100000000000000 0.772509659068896: 0.0100000000000000 mean = 4.82829 HPD intervals: HPD interval (0.84): 3.29708654587021..6.54706611602896 HPD interval (0.9): 2.80311256777640..6.86001046880149 HPD interval (0.99): 0.97302681811760..6.86001046880149 HPD interval (0.999): 0.77250965906890..6.86001046880149 HPD interval (0.9999): 0.77250965906890..6.86001046880149 HPD interval (0.99999): 0.77250965906890..6.86001046880149 var : xi Probabilities (truncated): 0.998067158273453: 0.0100000000000000 0.652586054826428: 0.0100000000000000 0.645551212898247: 0.0100000000000000 0.604356387445869: 0.0100000000000000 0.523834964038727: 0.0100000000000000 0.509092181692409: 0.0100000000000000 0.423090683027678: 0.0100000000000000 0.402441032883917: 0.0100000000000000 0.37070964247487: 0.0100000000000000 0.327880111209992: 0.0100000000000000 ......... -0.840283108800781: 0.0100000000000000 -0.875151826010622: 0.0100000000000000 -0.8783015161186: 0.0100000000000000 -0.883483478745205: 0.0100000000000000 -0.884227202219994: 0.0100000000000000 -0.891821089150301: 0.0100000000000000 -0.920209294147887: 0.0100000000000000 -0.945072983366006: 0.0100000000000000 -0.953790769890785: 0.0100000000000000 -0.965038130043558: 0.0100000000000000 mean = -0.326927 HPD intervals: HPD interval (0.84): -0.96503813004356..0.07755458963921 HPD interval (0.9): -0.96503813004356..0.24573496228351 HPD interval (0.99): -0.96503813004356..0.65258605482643 HPD interval (0.999): -0.96503813004356..0.99806715827345 HPD interval (0.9999): -0.96503813004356..0.99806715827345 HPD interval (0.99999): -0.96503813004356..0.99806715827345 var : post Probabilities (truncated): 204.522066843925586: 0.0100000000000000 159.627825244437645: 0.0100000000000000 156.763982472568102: 0.0100000000000000 156.673198765303709: 0.0100000000000000 155.861808723863987: 0.0100000000000000 155.741182250659307: 0.0100000000000000 155.64384800273001: 0.0100000000000000 155.397440587728852: 0.0100000000000000 155.314530741149241: 0.0100000000000000 154.746689182072771: 0.0100000000000000 ......... 141.8638456871426: 0.0100000000000000 141.549931489141471: 0.0100000000000000 141.53635236276881: 0.0100000000000000 141.270339633306691: 0.0100000000000000 141.048109300131586: 0.0100000000000000 140.9180010455608: 0.0100000000000000 140.619264597736873: 0.0100000000000000 139.911690924914012: 0.0100000000000000 136.937143235785214: 0.0100000000000000 135.386040768918434: 0.0100000000000000 mean = 148.647 HPD intervals: HPD interval (0.84): 140.91800104556080..153.90909265066151 HPD interval (0.9): 141.27033963330669..155.86180872386399 HPD interval (0.99): 135.38604076891843..159.62782524443764 HPD interval (0.999): 135.38604076891843..204.52206684392559 HPD interval (0.9999): 135.38604076891843..204.52206684392559 HPD interval (0.99999): 135.38604076891843..204.52206684392559 var : p Probabilities: false: 0.9900000000000000 true: 0.0100000000000000 mean = 0.01 HPD intervals: show_hpd_intervals: data is not numeric [min = 135.386,mean = 148.647,median = 148.114,max = 204.522,variance = 54.4297,stdev = 7.37765] Histogram (total 100) 135.386: 1 #### (0.010 / 0.010) 137.114: 1 #### (0.010 / 0.020) 138.843: 0 (0.000 / 0.020) 140.571: 5 ##################### (0.050 / 0.070) 142.300: 10 ########################################### (0.100 / 0.170) 144.028: 10 ########################################### (0.100 / 0.270) 145.756: 14 ############################################################ (0.140 / 0.410) 147.485: 11 ############################################### (0.110 / 0.520) 149.213: 11 ############################################### (0.110 / 0.630) 150.942: 12 ################################################### (0.120 / 0.750) 152.670: 10 ########################################### (0.100 / 0.850) 154.398: 6 ########################## (0.060 / 0.910) 156.127: 7 ############################## (0.070 / 0.980) 157.855: 0 (0.000 / 0.980) 159.584: 1 #### (0.010 / 0.990) 161.312: 0 (0.000 / 0.990) 163.040: 0 (0.000 / 0.990) 164.769: 0 (0.000 / 0.990) 166.497: 0 (0.000 / 0.990) 168.226: 0 (0.000 / 0.990) 169.954: 0 (0.000 / 0.990) 171.682: 0 (0.000 / 0.990) 173.411: 0 (0.000 / 0.990) 175.139: 0 (0.000 / 0.990) 176.868: 0 (0.000 / 0.990) 178.596: 0 (0.000 / 0.990) 180.324: 0 (0.000 / 0.990) 182.053: 0 (0.000 / 0.990) 183.781: 0 (0.000 / 0.990) 185.510: 0 (0.000 / 0.990) 187.238: 0 (0.000 / 0.990) 188.966: 0 (0.000 / 0.990) 190.695: 0 (0.000 / 0.990) 192.423: 0 (0.000 / 0.990) 194.152: 0 (0.000 / 0.990) 195.880: 0 (0.000 / 0.990) 197.608: 0 (0.000 / 0.990) 199.337: 0 (0.000 / 0.990) 201.065: 0 (0.000 / 0.990) 202.794: 1 #### (0.010 / 1.000) Mathematica: FindDistributionParameters[data,MaxStableDistribution[mu, sigma, xi]] -> {xi -> -0.042961, sigma -> 4.21246, mu -> 146.162} */ go3 ?=> Data = [ normal_dist_n(100,15,1000).max : _ in 1..100], println(data=Data), show_simple_stats(Data), show_histogram(Data), reset_store, run_model(10_000,$model3(Data),[show_probs_trunc,mean, % show_percentiles, truncate_size=10, show_hpd_intervals,hpd_intervals=[0.84,0.9,0.99,0.999,0.9999,0.99999], % show_histogram, min_accepted_samples=100,show_accepted_samples=true ]), Post = get_store().get("post"), show_simple_stats(Post), show_histogram(Post), nl, % show_store_lengths,nl, % fail, nl. go3 => true. model3(Data) => Mu = normal_dist(Data.mean,Data.stdev), Sigma = uniform(0,20), Xi = uniform(-1,1), Y = generalized_extreme_value_dist_n(Mu,Sigma,Xi,Data.len), observe_abc(Data,Y,1/10), Post = generalized_extreme_value_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.