/* Earthquake Extreme Value Theory using Gumbel distribution in Picat. Example and data from Mathematica GumbelDistribution: yearly maximum magnitude of earthquakes in the United States in the past 200 years: """ mYdist = EstimatedDistribution(maxYearly, GumbelDistribution(Alpha, Beta)) -> GumbelDistribution(7.03496, 0.745586) ... Find the probability of the annual maximum earthquake having a magnitude of at least 6: Probability(x >= 6, x -> mYdist) -> 0.779155 Simulate the magnitudes of the annual maximum earthquake for 30 years: maxMag = RandomVariate(mYdist, 30) -> (6.60781, 5.93618, 7.73324, 6.22724, 3.95527, 6.11688, 7.49837, 6.53918, 8.12725, 5.62827, 7.77249, 6.60604, 7.20895, 6.6224, 7.78966, 6.77433, 6.95503, 7.08528, 6.34501, 6.95825, 6.48127, 7.54685, 7.37537, 6.6997, 7.07498, 7.34062, 6.74814, 7.24383, 5.45639, 6.85268) Max@% -> 8.12725 """ Cf my Gamble model gamble_gumbel_earthquake.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. /* Data: [min = 4.2,mean = 6.61987,median = 6.7,max = 8.4,variance = 0.777751,stdev = 0.881902] Histogram (total 151) 4.200: 1 #### (0.007 / 0.007) 4.305: 0 (0.000 / 0.007) 4.410: 1 #### (0.007 / 0.013) 4.515: 2 ######## (0.013 / 0.026) 4.620: 3 ########### (0.020 / 0.046) 4.725: 0 (0.000 / 0.046) 4.830: 2 ######## (0.013 / 0.060) 4.935: 1 #### (0.007 / 0.066) 5.040: 1 #### (0.007 / 0.073) 5.145: 1 #### (0.007 / 0.079) 5.250: 2 ######## (0.013 / 0.093) 5.355: 0 (0.000 / 0.093) 5.460: 5 ################### (0.033 / 0.126) 5.565: 1 #### (0.007 / 0.132) 5.670: 1 #### (0.007 / 0.139) 5.775: 3 ########### (0.020 / 0.159) 5.880: 6 ####################### (0.040 / 0.199) 5.985: 9 ################################## (0.060 / 0.258) 6.090: 3 ########### (0.020 / 0.278) 6.195: 7 ########################## (0.046 / 0.325) 6.300: 5 ################### (0.033 / 0.358) 6.405: 1 #### (0.007 / 0.364) 6.510: 6 ####################### (0.040 / 0.404) 6.615: 1 #### (0.007 / 0.411) 6.720: 16 ############################################################ (0.106 / 0.517) 6.825: 4 ############### (0.026 / 0.543) 6.930: 5 ################### (0.033 / 0.576) 7.035: 10 ###################################### (0.066 / 0.642) 7.140: 12 ############################################# (0.079 / 0.722) 7.245: 8 ############################## (0.053 / 0.775) 7.350: 12 ############################################# (0.079 / 0.854) 7.455: 1 #### (0.007 / 0.861) 7.560: 4 ############### (0.026 / 0.887) 7.665: 5 ################### (0.033 / 0.921) 7.770: 3 ########### (0.020 / 0.940) 7.875: 4 ############### (0.026 / 0.967) 7.980: 1 #### (0.007 / 0.974) 8.085: 1 #### (0.007 / 0.980) 8.190: 1 #### (0.007 / 0.987) 8.295: 2 ######## (0.013 / 1.000) Parameter recovery: min_accepted_samples = 100 var : mu Probabilities (truncated): 7.497517925058419: 0.0100000000000000 7.477460438997932: 0.0100000000000000 7.468260511179835: 0.0100000000000000 7.46628186148715: 0.0100000000000000 ......... 6.563004215835451: 0.0100000000000000 6.553511703673153: 0.0100000000000000 6.530408974326853: 0.0100000000000000 6.486204620002434: 0.0100000000000000 mean = 6.9904 var : sigma Probabilities (truncated): 1.171256485940542: 0.0100000000000000 1.169427997045884: 0.0100000000000000 1.121369409897071: 0.0100000000000000 1.079707984384013: 0.0100000000000000 ......... 0.369432526812624: 0.0100000000000000 0.355782155113193: 0.0100000000000000 0.349372490471868: 0.0100000000000000 0.342439697283525: 0.0100000000000000 mean = 0.712491 var : p Probabilities: true: 0.7700000000000000 false: 0.2300000000000000 mean = 0.77 var : post Probabilities (truncated): 8.619461332855753: 0.0100000000000000 8.087081857857973: 0.0100000000000000 8.037239707989007: 0.0100000000000000 7.854395772268506: 0.0100000000000000 ......... 3.876574243306138: 0.0100000000000000 3.861388184641182: 0.0100000000000000 3.623479374730913: 0.0100000000000000 3.051415881773203: 0.0100000000000000 mean = 6.47751 var : max Probabilities (truncated): 9.291442611730488: 0.0100000000000000 9.239551878446843: 0.0100000000000000 9.100978482909309: 0.0100000000000000 8.98285885647635: 0.0100000000000000 ......... 7.341344923089266: 0.0100000000000000 7.328026299732573: 0.0100000000000000 7.119568173640898: 0.0100000000000000 7.070199808252523: 0.0100000000000000 mean = 8.1805 Post: [min = 3.05142,mean = 6.47751,median = 6.66334,max = 8.61946,variance = 1.02809,stdev = 1.01395] Histogram (total 100) 3.051: 1 ###### (0.010 / 0.010) 3.191: 0 (0.000 / 0.010) 3.330: 0 (0.000 / 0.010) 3.469: 0 (0.000 / 0.010) 3.608: 1 ###### (0.010 / 0.020) 3.747: 0 (0.000 / 0.020) 3.887: 3 ################## (0.030 / 0.050) 4.026: 0 (0.000 / 0.050) 4.165: 0 (0.000 / 0.050) 4.304: 1 ###### (0.010 / 0.060) 4.443: 0 (0.000 / 0.060) 4.583: 2 ############ (0.020 / 0.080) 4.722: 0 (0.000 / 0.080) 4.861: 1 ###### (0.010 / 0.090) 5.000: 0 (0.000 / 0.090) 5.139: 0 (0.000 / 0.090) 5.279: 4 ######################## (0.040 / 0.130) 5.418: 2 ############ (0.020 / 0.150) 5.557: 2 ############ (0.020 / 0.170) 5.696: 2 ############ (0.020 / 0.190) 5.835: 1 ###### (0.010 / 0.200) 5.975: 4 ######################## (0.040 / 0.240) 6.114: 3 ################## (0.030 / 0.270) 6.253: 6 #################################### (0.060 / 0.330) 6.392: 6 #################################### (0.060 / 0.390) 6.531: 5 ############################## (0.050 / 0.440) 6.671: 10 ############################################################ (0.100 / 0.540) 6.810: 10 ############################################################ (0.100 / 0.640) 6.949: 8 ################################################ (0.080 / 0.720) 7.088: 5 ############################## (0.050 / 0.770) 7.227: 4 ######################## (0.040 / 0.810) 7.367: 7 ########################################## (0.070 / 0.880) 7.506: 3 ################## (0.030 / 0.910) 7.645: 2 ############ (0.020 / 0.930) 7.784: 3 ################## (0.030 / 0.960) 7.923: 1 ###### (0.010 / 0.970) 8.063: 2 ############ (0.020 / 0.990) 8.202: 0 (0.000 / 0.990) 8.341: 0 (0.000 / 0.990) 8.480: 1 ###### (0.010 / 1.000) Simulate the magnitudes of the annual maximum earthquake for 30 years: [5.89749,6.81716,7.11757,5.93636,6.10053,6.46253,6.88279,5.89685,4.60106,5.23501,5.67748,6.05007,7.35823,7.38467,6.66702,7.60522,4.0078,7.99416,5.07711,6.24894,6.74539,6.38831,7.09759,6.54674,7.28041,7.08274,7.5861,7.06819,6.3332,7.11866] */ go ?=> % This is from Mathematica GumbelDistribution: % The yearly maximum value of earthquakes in USA the past 200 year Data = [6.0,4.5,5.3,6.5,4.4,6.0,7.2,7.4,5.0,5.5,4.6,4.5,6.8,7.0, 4.6,6.0,4.2,4.9,6.5,4.6,6.0,5.5,7.6,6.1,6.3,5.6,5.9,4.8, 5.9,6.3,5.8,5.1,7.9,6.1,6.0,7.0,7.3,6.7,5.8,4.8,6.2,5.8, 6.0,6.2,6.0,5.9,6.2,6.7,6.3,5.9,6.0,6.3,5.5,6.7,5.5,5.9, 5.9,6.0,6.4,7.6,8.0,7.7,7.1,7.0,7.0,7.3,7.4,7.8,7.4,7.0,7.4, 7.0,7.1,7.3,7.2,6.7,7.7,7.7,7.9,6.8,6.1,5.5,5.2,7.3,7.2, 5.7,6.7,7.0,7.1,6.8,7.8,6.5,6.2,7.2,6.9,7.1,7.1,6.2,7.3, 8.3,6.2,7.4,6.7,6.9,7.4,7.1,6.7,7.3,7.2,7.5,6.9,6.7,7.1, 7.2,7.1,7.2,7.0,6.8,8.1,7.9,7.7,6.9,6.7,6.7,6.7,8.4,8.2, 7.0,6.7,7.1,6.5,6.7,7.1,7.6,6.7,6.5,7.6,6.3,6.7,6.7,7.1, 6.9,7.0,6.2,7.2,6.6,6.5,7.9,7.8,7.7,7.1], println("Data:"), show_simple_stats(Data), show_histogram(Data), nl, println("Parameter recovery:"), reset_store, run_model(1000,$model(Data),[show_probs_trunc,mean, % show_percentiles, % show_hpd_intervals,hpd_intervals=[0.94], % show_histogram, min_accepted_samples=100,show_accepted_samples=false ]), println("Post:"), Store = get_store(), Post = Store.get("post"), show_simple_stats(Post), show_histogram(Post), nl, Mu = Store.get("mu").mean, Sigma = Store.get("sigma").mean, println("\nSimulate the magnitudes of the annual maximum earthquake for 30 years:"), println(gumbel_dist_n(Mu,Sigma,30)), nl, % show_store_lengths,nl, % fail, nl. go => true. model(Data) => Mu = normal_dist(Data.mean,Data.stdev), Sigma = uniform(0,10), Y = gumbel_dist_n(Mu,Sigma,Data.len), observe_abc(Data,Y,1/2), Post = gumbel_dist(Mu,Sigma), P = check(Post >= 6), % Find the max value Max = Y.max, if observed_ok then add("mu",Mu), add("sigma",Sigma), add("p",P), add("post",Post), add("max",Max) end.