/* Extreme value distribution - maximum wind speed 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 """ Cf - ppl_extreme_value_dist.pi - my Gamble model gamble_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 = 720,mean = 52.5545,stdev = 12.042] var : a Probabilities (truncated): 62.874051810226895: 0.0010000000000000 61.929552305336934: 0.0010000000000000 61.857019309257609: 0.0010000000000000 61.781833972439074: 0.0010000000000000 ......... 42.968414255373482: 0.0010000000000000 42.914218188183312: 0.0010000000000000 42.697416496892103: 0.0010000000000000 42.411156075892052: 0.0010000000000000 mean = 48.1618 HPD intervals: HPD interval (0.94): 42.98947097420529..57.11837922526808 var : b Probabilities (truncated): 12.702494388372724: 0.0010000000000000 12.570774885687891: 0.0010000000000000 12.549804444641955: 0.0010000000000000 12.439136836618793: 0.0010000000000000 ......... -11.976856395460882: 0.0010000000000000 -12.042615911799889: 0.0010000000000000 -12.102116876646225: 0.0010000000000000 -12.18359017555678: 0.0010000000000000 mean = 7.90005 HPD intervals: HPD interval (0.94): -8.01079291919938..12.57077488568789 var : post Probabilities (truncated): 111.13291453837968: 0.0010000000000000 101.188786286639811: 0.0010000000000000 100.125990004461428: 0.0010000000000000 99.487054522583691: 0.0010000000000000 ......... 25.56190256868549: 0.0010000000000000 25.241102873723104: 0.0010000000000000 23.037283933947528: 0.0010000000000000 22.94556226016217: 0.0010000000000000 mean = 53.0027 HPD intervals: HPD interval (0.94): 32.75924297772318..76.34705027057161 var : p Probabilities: false: 0.9950000000000000 true: 0.0050000000000000 mean = [false = 0.995,true = 0.005] HPD intervals: show_hpd_intervals: data is not numeric */ go ?=> % The data from Mathematica ExtremeValueDistribution: Data = [66.6,55.8,55.8,51.84,50.04,50.04,42.48,38.88,66.6,44.64, 96.48,55.8,48.24,72.36,55.8,53.64,53.64,44.64,42.48,48.24, 44.64,44.64,74.16,50.04,55.8,77.76,55.8,59.4,51.84,48.24, 37.08,42.48,55.8,44.64,53.64,59.4,63,64.8,55.8,55.8,51.84, 46.44,46.44,55.8,48.24,57.6,88.92,55.8,55.8,50.04,59.4, 53.64,61.2,51.84,44.64,138.96,113.04,51.84,68.76,74.16, 77.76,61.2,70.56,70.56,51.84,40.68,44.64,63,55.8,57.6, 72.36,55.8,81.72,70.56,79.56,77.76,63,46.44,51.84,37.08, 44.64,55.8,85.32,77.76,66.6,61.2,66.6,63,63,51.84,55.8, 42.48,46.44,42.48,55.8,55.8,63,64.8,63,53.64,44.64,51.84, 48.24,55.8,38.88,50.04,70.56,55.8,59.4,55.8,61.2,48.24, 38.88,38.88,40.68,33.48,35.28,42.48,57.6,63,61.2,64.8, 74.16,53.64,44.64,50.04,46.44,37.08,81.72,44.64,46.44,59.4, 55.8,63,64.8,63,44.64,42.48,40.68,42.48,66.6,55.8,61.2, 61.2,61.2,66.6,70.56,66.6,55.8,61.2,44.64,53.64,63,53.64, 57.6,74.16,70.56,55.8,53.64,77.76,51.84,50.04,46.44,53.64, 51.84,55.8,70.56,57.6,59.4,64.8,50.04,48.24,51.84,46.44, 42.48,40.68,46.44,48.24,53.64,61.2,57.6,53.64,44.64,40.68, 46.44,38.88,37.08,46.44,48.24,44.64,55.8,40.68,61.2,48.24, 55.8,38.88,46.44,33.48,44.64,33.48,44.64,46.44,61.2,66.6, 50.04,59.4,46.44,59.4,77.76,44.64,37.08,44.64,42.48,48.24, 44.64,48.24,55.8,55.8,59.4,42.48,51.84,37.08,35.28,38.88, 46.44,38.88,63,51.84,59.4,55.8,64.8,59.4,46.44,44.64,44.64, 42.48,42.48,46.44,46.44,55.8,55.8,72.36,57.6,70.56,37.08, 74.16,42.48,42.48,40.68,50.04,59.4,64.8,63,64.8,74.16,55.8, 46.44,33.48,35.28,55.8,37.08,38.88,48.24,51.84,64.8,64.8, 48.24,37.08,46.44,48.24,33.48,79.56,51.84,50.04,66.6,61.2, 63,53.64,61.2,61.2,42.48,37.08,29.52,40.68,77.76,55.8,55.8, 51.84,63,68.76,55.8,59.04,55.44,46.44,37.08,44.64,51.84, 40.68,48.24,61.2,50.04,59.4,51.84,57.6,40.68,46.44,38.88, 42.48,37.08,40.68,44.64,46.44,66.6,59.4,64.8,86.76,50.04, 38.88,47.88,55.8,77.76,55.8,51.84,57.6,64.8,55.8,66.6, 50.04,66.6,74.16,37.08,84.96,40.68,53.64,59.4,51.84,74.16, 81.72,48.24,44.64,50.04,38.88,92.88,46.08,42.48,37.08,44.64, 92.52,64.8,51.84,46.44,55.8,59.4,38.88,35.28,42.48,48.24, 51.84,55.8,81.72,63,42.48,59.4,79.56,57.6,37.08,42.48, 46.08,48.24,66.6,48.24,57.6,48.24,74.16,50.04,51.84,81.72, 40.68,37.08,38.88,40.68,40.68,38.88,48.24,59.4,51.84,50.04, 64.8,44.28,48.24,40.68,46.44,35.28,48.24,55.8,53.64,51.84, 64.8,105.48,59.4,40.68,40.68,46.44,48.24,38.88,46.44,51.84, 64.8,40.68,53.64,64.8,44.64,44.64,38.88,40.68,35.28,38.88, 46.44,46.44,37.08,44.64,51.84,40.68,46.44,44.64,46.44,40.68, 37.08,85.32,40.68,44.64,44.64,48.24,38.88,51.84,44.64,37.08, 40.68,37.08,40.68,42.48,44.64,61.2,64.8,55.8,46.44,59.4, 59.4,46.44,40.68,37.08,44.64,46.44,63,50.04,48.24,50.04, 59.4,48.24,46.44,46.44,72.36,40.68,46.44,44.64,66.6,55.44, 48.24,53.64,53.64,55.44,100.08,50.04,48.24,59.4,85.32,55.44, 46.44,61.2,51.84,53.64,53.64,59.4,51.84,51.84,42.48,35.28, 38.88,40.68,59.4,57.6,51.84,50.04,59.4,59.4,61.2,46.44, 53.64,38.88,75.96,48.24,59.4,46.44,50.04,64.8,50.04,57.6, 44.64,53.64,48.24,38.88,40.68,38.88,46.44,55.8,81.72,50.04, 49.68,87.12,48.24,49.68,48.24,55.44,40.68,37.08,51.84,53.64, 53.64,70.2,50.04,55.8,51.84,44.64,46.08,38.88,42.48,47.88, 38.88,66.6,77.76,46.44,68.76,53.64,68.4,48.24,35.28,37.08, 44.64,47.88,53.64,64.8,50.04,74.16,53.64,46.44,55.8,51.84, 35.28,46.44,29.52,53.64,70.2,51.84,59.4,50.04,50.04,63, 61.2,50.04,44.28,38.88,46.44,40.68,51.84,68.4,48.24,46.44, 59.4,61.2,59.4,50.04,46.44,38.88,38.88,75.96,50.04,50.04, 51.84,57.6,42.48,61.2,46.44,44.64,42.48,57.24,35.28,51.84, 63,63,51.84,59.4,53.64,51.84,59.4,51.84,53.64,43.2,38.88, 38.88,50.04,38.88,63,40.68,55.8,64.8,44.64,51.84,59.4, 42.48,42.48,44.64,50.04,48.24,40.68,53.64,46.44,48.24,51.84, 55.8,44.28,38.88,42.48,55.8,50.04,53.64,57.6,51.84,51.84, 51.84,42.48,37.08,37.08,37.08,35.28,38.88,59.4,61.2,59.4, 46.44,40.68,48.24,46.44,38.88,42.48,35.28,33.48,46.44,37.08, 59.4,55.44,66.6,50.04,59.4,46.44,61.2,40.68,42.48,59.4, 59.4,66.6,44.64,51.84,57.24,57.24,44.28,50.04,51.84,55.44, 37.08,38.88,40.68,55.44,50.04,51.84,46.44,61.2,55.44,63, 46.44,38.88,37.08,42.48,38.88,51.84,44.28,57.24,51.84,51.84, 51.84,50.04,48.24,50.04,42.48,37.08,53.64,53.64,53.64,55.44, 51.84,51.84,42.48,46.44,51.84,48.24,46.44,37.08,46.44,46.44, 48.24,59.4], println([len=Data.len,mean=Data.mean,stdev=Data.stdev]), reset_store, run_model(10_000,$model(Data),[show_probs_trunc,mean, % show_percentiles,show_histogram, show_hpd_intervals,hpd_intervals=[0.94], min_accepted_samples=10_000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go => true. model(Data) => A = normal_dist(Data.mean,10), B = normal_dist(Data.stdev, 10), Y = extreme_value_dist_n(A,B,Data.len), Post = extreme_value_dist(A,B), P = check(Post > 96.5606), % 60mph -> 96.5606km/h observe_abc(Data,Y,1), if observed_ok then add("a",A), add("b",B), add("post",Post), add("p",P), end. /* This is the data that's used in my Gamble model and here we got about the same results. [len = 72,mean = 73.56,stdev = 13.8656] var : a Probabilities (truncated): 86.236432308256653: 0.0001000000000000 85.902854851150465: 0.0001000000000000 85.777356498032347: 0.0001000000000000 85.381663788746522: 0.0001000000000000 ......... 59.693867650269482: 0.0001000000000000 59.671700035950664: 0.0001000000000000 59.439915655417778: 0.0001000000000000 58.522769135852968: 0.0001000000000000 mean = 68.1968 HPD intervals: HPD interval (0.94): 61.41830538613792..73.43623166002418 var : b Probabilities (truncated): 20.071377567876944: 0.0001000000000000 19.22307322564162: 0.0001000000000000 19.036916833364444: 0.0001000000000000 18.8566662855422: 0.0001000000000000 ......... -15.345150917606784: 0.0001000000000000 -15.715566662448243: 0.0001000000000000 -15.861173638410431: 0.0001000000000000 -16.316688139488335: 0.0001000000000000 mean = 10.2366 HPD intervals: HPD interval (0.94): 6.91500282057589..16.53853248813058 var : post Probabilities (truncated): 180.591407682469622: 0.0001000000000000 179.790088342425747: 0.0001000000000000 168.031537930471814: 0.0001000000000000 167.442432732749836: 0.0001000000000000 ......... 29.623228080154746: 0.0001000000000000 27.839605222268531: 0.0001000000000000 20.350272703813175: 0.0001000000000000 17.530597188222615: 0.0001000000000000 mean = 74.0833 HPD intervals: HPD interval (0.94): 49.25347331413262..102.65512248439688 var : p Probabilities: false: 0.9270000000000000 true: 0.0730000000000000 mean = [false = 0.927,true = 0.073] HPD intervals: show_hpd_intervals: data is not numeric */ go2 ?=> % 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], println([len=Data.len,mean=Data.mean,stdev=Data.stdev]), reset_store, run_model(10_000,$model2(Data),[show_probs_trunc,mean, % show_percentiles,show_histogram, show_hpd_intervals,hpd_intervals=[0.94], min_accepted_samples=1000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go2 => true. model2(Data) => A = normal_dist(Data.mean,10), B = normal_dist(Data.stdev, 10), Y = extreme_value_dist_n(A,B,Data.len), Post = extreme_value_dist(A,B), P = check(Post > 96.5606), % 60mph -> 96.5606km/h observe_abc(Data,Y,1/4), if observed_ok then add("a",A), add("b",B), add("post",Post), add("p",P), end.