/* GEV for maximum wind speeds 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 """ The original data was in km/h but the last query is in mph, so let's convert 60mph to km/h: """ UnitConvert[Quantity[60, ("Miles")/("Hours")], "km/h"] // N -> Quantity[96.5606, ("Kilometers")/("Hours")] Probability[x > QuantityMagnitude@%, x \[Distributed] edist] 0.00454842 """ Note: Here we use generalized_extreme_value_dist/3 instead. Cf my Gamble model gamble_generalized_extreme_value_maximum_wind_speed.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. /* [len = 72,min = 55.44,mean = 73.56,median = 72.36,max = 138.96,variance = 192.254,stdev = 13.8656] Histogram (total 72) 55.440: 3 ############## (0.042 / 0.042) 58.423: 7 ################################ (0.097 / 0.139) 61.406: 2 ######### (0.028 / 0.167) 64.389: 13 ############################################################ (0.181 / 0.347) 67.371: 9 ########################################## (0.125 / 0.472) 70.354: 1 ##### (0.014 / 0.486) 73.337: 9 ########################################## (0.125 / 0.611) 76.320: 9 ########################################## (0.125 / 0.736) 79.303: 2 ######### (0.028 / 0.764) 82.286: 4 ################## (0.056 / 0.819) 85.269: 3 ############## (0.042 / 0.861) 88.251: 3 ############## (0.042 / 0.903) 91.234: 0 (0.000 / 0.903) 94.217: 1 ##### (0.014 / 0.917) 97.200: 3 ############## (0.042 / 0.958) 100.183: 1 ##### (0.014 / 0.972) 103.166: 0 (0.000 / 0.972) 106.149: 1 ##### (0.014 / 0.986) 109.131: 0 (0.000 / 0.986) 112.114: 0 (0.000 / 0.986) 115.097: 0 (0.000 / 0.986) 118.080: 0 (0.000 / 0.986) 121.063: 0 (0.000 / 0.986) 124.046: 0 (0.000 / 0.986) 127.029: 0 (0.000 / 0.986) 130.011: 0 (0.000 / 0.986) 132.994: 0 (0.000 / 0.986) 135.977: 1 ##### (0.014 / 1.000) HPD intervals: HPD interval (0.84): 55.44000000000000..85.31999999999999 HPD interval (0.9): 55.44000000000000..88.92000000000000 HPD interval (0.99): 55.44000000000000..138.96000000000001 HPD interval (0.999): 55.44000000000000..138.96000000000001 HPD interval (0.9999): 55.44000000000000..138.96000000000001 HPD interval (0.99999): 55.44000000000000..138.96000000000001 methods = [median,iqr,mad,q10,q90,kurtosis] explain = Data appears strongly skewed and heavy-tailed. Using robust shape statistics. min_accepted_samples = 1000 var : mu Probabilities (truncated): 79.988881806192097: 0.0010000000000000 79.666861697004208: 0.0010000000000000 79.544045552477485: 0.0010000000000000 79.475753928559882: 0.0010000000000000 ......... 53.794345225341431: 0.0010000000000000 53.785597635972657: 0.0010000000000000 53.781464141860567: 0.0010000000000000 53.36076012095441: 0.0010000000000000 mean = 68.5333 HPD intervals: HPD interval (0.84): 61.70152173329794..77.38578969488222 HPD interval (0.9): 58.76943581946869..77.03932592045813 HPD interval (0.99): 54.93360312869081..79.20853753252395 HPD interval (0.999): 53.78146414186057..79.98888180619210 HPD interval (0.9999): 53.36076012095441..79.98888180619210 HPD interval (0.99999): 53.36076012095441..79.98888180619210 var : sigma Probabilities (truncated): 9.992676028047072: 0.0010000000000000 9.986632699140642: 0.0010000000000000 9.980842531649788: 0.0010000000000000 9.980226960023971: 0.0010000000000000 ......... 0.621740641361913: 0.0010000000000000 0.544566288844015: 0.0010000000000000 0.535759008739031: 0.0010000000000000 0.531018432477032: 0.0010000000000000 mean = 5.81061 HPD intervals: HPD interval (0.84): 2.26565337379726..9.10608009859271 HPD interval (0.9): 2.12172720680094..9.60123783424554 HPD interval (0.99): 0.97898673777421..9.99267602804707 HPD interval (0.999): 0.53101843247703..9.98663269914064 HPD interval (0.9999): 0.53101843247703..9.99267602804707 HPD interval (0.99999): 0.53101843247703..9.99267602804707 var : xi Probabilities (truncated): 0.998325931838865: 0.0010000000000000 0.998274889308156: 0.0010000000000000 0.996427352538531: 0.0010000000000000 0.990325803863968: 0.0010000000000000 ......... 0.002141065430893: 0.0010000000000000 0.00139811076289: 0.0010000000000000 0.001296414062984: 0.0010000000000000 0.000707593746813: 0.0010000000000000 mean = 0.35403 HPD intervals: HPD interval (0.84): 0.00070759374681..0.60511904377729 HPD interval (0.9): 0.00129641406298..0.70430315318718 HPD interval (0.99): 0.00070759374681..0.96566610176380 HPD interval (0.999): 0.00129641406298..0.99832593183887 HPD interval (0.9999): 0.00070759374681..0.99832593183887 HPD interval (0.99999): 0.00070759374681..0.99832593183887 var : post Probabilities (truncated): 256.21380990367777: 0.0010000000000000 236.585616554346103: 0.0010000000000000 215.00475109088643: 0.0010000000000000 214.808012202518796: 0.0010000000000000 ......... 49.21094374017558: 0.0010000000000000 49.140074096972008: 0.0010000000000000 47.977468364177653: 0.0010000000000000 46.665078276642241: 0.0010000000000000 mean = 74.1215 HPD intervals: HPD interval (0.84): 58.35382213592970..87.60811926655106 HPD interval (0.9): 54.46045653175763..90.67760783086450 HPD interval (0.99): 47.97746836417765..148.16452658204486 HPD interval (0.999): 46.66507827664224..236.58561655434610 HPD interval (0.9999): 46.66507827664224..256.21380990367777 HPD interval (0.99999): 46.66507827664224..256.21380990367777 var : p Probabilities: false: 0.9510000000000000 true: 0.0490000000000000 mean = 0.049 HPD intervals: show_hpd_intervals: data is not numeric var : p2 Probabilities: false: 0.9890000000000000 true: 0.0110000000000000 mean = 0.011 HPD intervals: show_hpd_intervals: data is not numeric Post: [len = 1000,min = 46.6651,mean = 74.1215,median = 71.1463,max = 256.214,variance = 314.125,stdev = 17.7236] Histogram (total 1000) 46.665: 4 # (0.004 / 0.004) 51.904: 22 ##### (0.022 / 0.026) 57.143: 65 ################ (0.065 / 0.091) 62.381: 122 ############################## (0.122 / 0.213) 67.620: 243 ############################################################ (0.243 / 0.456) 72.859: 229 ######################################################### (0.229 / 0.685) 78.097: 125 ############################### (0.125 / 0.810) 83.336: 68 ################# (0.068 / 0.878) 88.575: 46 ########### (0.046 / 0.924) 93.814: 26 ###### (0.026 / 0.950) 99.052: 13 ### (0.013 / 0.963) 104.291: 6 # (0.006 / 0.969) 109.530: 6 # (0.006 / 0.975) 114.768: 3 # (0.003 / 0.978) 120.007: 4 # (0.004 / 0.982) 125.246: 3 # (0.003 / 0.985) 130.485: 3 # (0.003 / 0.988) 135.723: 1 (0.001 / 0.989) 140.962: 0 (0.000 / 0.989) 146.201: 2 (0.002 / 0.991) 151.439: 0 (0.000 / 0.991) 156.678: 1 (0.001 / 0.992) 161.917: 1 (0.001 / 0.993) 167.156: 0 (0.000 / 0.993) 172.394: 0 (0.000 / 0.993) 177.633: 0 (0.000 / 0.993) 182.872: 0 (0.000 / 0.993) 188.110: 0 (0.000 / 0.993) 193.349: 1 (0.001 / 0.994) 198.588: 0 (0.000 / 0.994) 203.827: 1 (0.001 / 0.995) 209.065: 0 (0.000 / 0.995) 214.304: 3 # (0.003 / 0.998) 219.543: 0 (0.000 / 0.998) 224.782: 0 (0.000 / 0.998) 230.020: 0 (0.000 / 0.998) 235.259: 1 (0.001 / 0.999) 240.498: 0 (0.000 / 0.999) 245.736: 0 (0.000 / 0.999) 250.975: 1 (0.001 / 1.000) HPD intervals: HPD interval (0.84): 58.35382213592970..87.60811926655106 HPD interval (0.9): 54.46045653175763..90.67760783086450 HPD interval (0.99): 47.97746836417765..148.16452658204486 HPD interval (0.999): 46.66507827664224..236.58561655434610 HPD interval (0.9999): 46.66507827664224..256.21380990367777 HPD interval (0.99999): 46.66507827664224..256.21380990367777 */ go ?=> % Annual maximum wind speed in Boston 1938..2009 % From Mathematica (FrechetDistribution, see above) Data = [72.36,64.8,72.36,68.4,77.76,55.44,96.48,96.48,66.6,72.36, 61.2,55.8,96.48,74.16,77.76,88.92,138.96,77.76,85.32,66.6, 70.56,63,81.72,66.6,74.16,77.76,64.8,57.6,66.6,77.76,63, 64.8,74.16,74.16,79.56,77.76,68.76,59.4,86.76,84.96,92.88, 81.72,79.56,81.72,64.8,105.48,64.8,85.32,64.8,63,72.36, 100.08,59.4,75.96,81.72,87.12,77.76,68.76,74.16,68.4,75.96, 63,63,64.8,57.6,61.2,59.4,66.6,57.24,63,55.44,59.4], show_simple_stats(Data), show_histogram(Data), show_hpd_intervals(Data,[0.84,0.9,0.99,0.999,0.9999,0.99999]), [Methods,Explain] = recommend_abc_stats_explained(Data), println(methods=Methods), println(explain=Explain), reset_store, run_model(10_000,$model(Data),[show_probs_trunc,mean, % show_percentiles, 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 ]), println("Post:"), Post = get_store().get("post"), show_simple_stats(Post), show_histogram(Post), show_hpd_intervals(Post,[0.84,0.9,0.99,0.999,0.9999,0.99999]), nl, % show_store_lengths,nl, % fail, nl. go => true. model(Data) => Mu = normal_dist(Data.mean,Data.stdev), Sigma = uniform(0,10), Xi = uniform(0,1), Y = generalized_extreme_value_dist_n(Mu,Sigma,Xi,Data.len), % observe_abc(Data,Y,1/10), observe_abc(Data,Y,1,[median,iqr,mad,q10,q90,kurtosis]), Post = generalized_extreme_value_dist(Mu,Sigma,Xi), P = check(Post >= 96.5606), P2 = check(Post >= Data.max), if observed_ok then add("mu",Mu), add("sigma",Sigma), add("xi",Xi), add("post",Post), add("p",P), add("p2",P2), end.