/* Extreme value distribution in Picat. From https://www.randomservices.org/random/special/ExtremeValue.html """ The distribution is also known as the standard Gumbel distribution in honor of Emil Gumbel. As we will show below in (13), it arises as the limit of the maximum of independent random variables, each with the standard exponential distribution (when this maximum is appropriately centered). This fact is the main reason that the distribution is special, and is the reason for the name. For the remainder of this discussion, suppose that random variable has the standard Gumbel distribution. ... The quantile function G^-1 of V given by G^-1(p) = -ln(-ln(p)) p in 0..1 """ This is implemented as extreme_value_dist1. This seems to be about the same as Mathematica's ExtremeValueDistribution(0,1) From Mathematica ExtremeValueDistribution """ ExtremeValueDistribution(alpha,beta) represents an extreme value distribution with location parameter alpha and scale parameter beta --- The extreme value distribution gives the asymptotic distribution of the maximum value in a sample from a distribution such as the normal distribution. ... Quantile(ExtremeValueDistribution(a,b), x) -> a-b Log(-Log(x)) 0 < x < 1 Infinity x <= 0 if 0 <= x <= 1 Infinity True """ Cf my Gamble model gamble_extreme_value_distribution.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. main => go. /* var : d Probabilities (truncated): 9.05305559649878: 0.0001000000000000 7.647161225743283: 0.0001000000000000 7.629292328770541: 0.0001000000000000 7.600934078208708: 0.0001000000000000 ......... -2.17645556155996: 0.0001000000000000 -2.198562444626482: 0.0001000000000000 -2.366825154653017: 0.0001000000000000 -2.446137853283312: 0.0001000000000000 mean = 0.562238 HPD intervals: HPD interval (0.5): -0.64897424795701..0.80800771281626 HPD interval (0.84): -1.11386684403736..2.13665301674635 HPD interval (0.9): -1.43078300380912..2.44089865500852 HPD interval (0.95): -1.65192364533403..3.08928293378513 HPD interval (0.99): -1.86795888383507..4.72736936205418 HPD interval (0.999): -2.19856244462648..6.78075456631161 HPD interval (0.9999): -2.44613785328331..7.64716122574328 HPD interval (0.99999): -2.44613785328331..9.05305559649878 var : p Probabilities: true: 0.6298000000000000 false: 0.3702000000000000 mean = [true = 0.6298,false = 0.3702] HPD intervals: show_hpd_intervals: data is not numeric */ go ?=> reset_store, run_model(10_000,$model,[show_probs_trunc,mean,show_hpd_intervals]), nl, % show_store_lengths, % fail, nl. go => true. model() => D = extreme_value_dist(), % extreme_value_dist(0,1) P = check(D > 0), add("d",D), add("p",P). /* var : d Probabilities (truncated): 278.69228266709689: 0.0001000000000000 272.468040051222147: 0.0001000000000000 253.675471396821848: 0.0001000000000000 252.97571514290064: 0.0001000000000000 ......... 68.209043757513328: 0.0001000000000000 67.675577548020868: 0.0001000000000000 67.40646236271391: 0.0001000000000000 63.874646322336318: 0.0001000000000000 mean = 108.634 HPD intervals: HPD interval (0.5): 90.75889994407336..112.77123303588375 HPD interval (0.84): 81.06878589408197..129.02578177678339 HPD interval (0.9): 80.34262289515361..137.37527134384587 HPD interval (0.95): 76.87848537069999..146.61628900252916 HPD interval (0.99): 73.57836202064360..174.42701398983178 HPD interval (0.999): 63.87464632233632..214.11495828897830 HPD interval (0.9999): 63.87464632233632..272.46804005122215 HPD interval (0.99999): 63.87464632233632..278.69228266709689 Histogram: 63.875: 1 (0.000 / 0.000) 69.245: 16 # (0.002 / 0.002) 74.616: 89 #### (0.009 / 0.011) 79.986: 319 ############### (0.032 / 0.042) 85.356: 655 ############################## (0.066 / 0.108) 90.727: 1051 ################################################ (0.105 / 0.213) 96.097: 1248 ######################################################### (0.125 / 0.338) 101.468: 1305 ############################################################ (0.131 / 0.468) 106.838: 1223 ######################################################## (0.122 / 0.591) 112.209: 1029 ############################################### (0.103 / 0.694) 117.579: 818 ###################################### (0.082 / 0.775) 122.949: 613 ############################ (0.061 / 0.837) 128.320: 460 ##################### (0.046 / 0.883) 133.690: 353 ################ (0.035 / 0.918) 139.061: 233 ########### (0.023 / 0.941) 144.431: 186 ######### (0.019 / 0.960) 149.802: 120 ###### (0.012 / 0.972) 155.172: 76 ### (0.008 / 0.980) 160.543: 45 ## (0.004 / 0.984) 165.913: 42 ## (0.004 / 0.988) 171.283: 39 ## (0.004 / 0.992) 176.654: 19 # (0.002 / 0.994) 182.024: 10 (0.001 / 0.995) 187.395: 14 # (0.001 / 0.996) 192.765: 10 (0.001 / 0.997) 198.136: 5 (0.001 / 0.998) 203.506: 3 (0.000 / 0.998) 208.877: 5 (0.001 / 0.999) 214.247: 3 (0.000 / 0.999) 219.617: 2 (0.000 / 0.999) 224.988: 1 (0.000 / 0.999) 230.358: 1 (0.000 / 0.999) 235.729: 0 (0.000 / 0.999) 241.099: 1 (0.000 / 0.999) 246.470: 1 (0.000 / 1.000) 251.840: 2 (0.000 / 1.000) 257.211: 0 (0.000 / 1.000) 262.581: 0 (0.000 / 1.000) 267.951: 0 (0.000 / 1.000) 273.322: 2 (0.000 / 1.000) var : p Probabilities: true: 0.6350000000000000 false: 0.3650000000000000 mean = [true = 0.635,false = 0.365] HPD intervals: show_hpd_intervals: data is not numeric Histogram: false : 3650 ################################## (0.365 / 0.365) true : 6350 ############################################################ (0.635 / 1.000) */ go2 ?=> reset_store, run_model(10_000,$model2,[show_probs_trunc,mean,show_hpd_intervals,show_histogram]), nl, % show_store_lengths, % fail, nl. go2 => true. model2() => D = extreme_value_dist(100,15), P = check(D > 100), add("d",D), add("p",P). /* Trying to recover the parameters Picat> X=extreme_value_dist_n(10,5,10) X = [15.281819132444104,10.028270196182174,14.014207692800706,9.901431868321195,22.504941544790793,35.934094047899976,10.887493952334752,5.223940870946778,22.304607209505008,0.61166403311217] This is too hard for PPL... For just a few values it works a little better: X = [15.281819132444104,10.028270196182174,14.014207692800706], var : alpha Probabilities (truncated): 18.835902250761119: 0.0151515151515152 17.759030637219098: 0.0151515151515152 17.647883723325972: 0.0151515151515152 17.532151200590725: 0.0151515151515152 ......... 8.083275700026785: 0.0151515151515152 6.530470730052549: 0.0151515151515152 6.092193362346009: 0.0151515151515152 0.532628987232516: 0.0151515151515152 mean = 12.0814 HPD intervals: HPD interval (0.5): 10.77726959752723..14.34114283618571 HPD interval (0.84): 8.08327570002679..15.66007757357325 HPD interval (0.9): 8.08327570002679..17.53215120059073 HPD interval (0.95): 8.08327570002679..18.83590225076112 HPD interval (0.99): 0.53262898723252..18.83590225076112 HPD interval (0.999): 0.53262898723252..18.83590225076112 HPD interval (0.9999): 0.53262898723252..18.83590225076112 HPD interval (0.99999): 0.53262898723252..18.83590225076112 var : beta Probabilities (truncated): 8.908524857325723: 0.0151515151515152 8.855256367873054: 0.0151515151515152 8.088239812333249: 0.0151515151515152 7.65414126573789: 0.0151515151515152 ......... 1.656144210908629: 0.0151515151515152 1.604891005719495: 0.0151515151515152 1.154426597596345: 0.0151515151515152 1.122260480710427: 0.0151515151515152 mean = 4.30281 HPD intervals: HPD interval (0.5): 2.60837067505735..5.03178365297233 HPD interval (0.84): 1.60489100571950..6.43583061473250 HPD interval (0.9): 1.12226048071043..6.75197821890562 HPD interval (0.95): 1.12226048071043..7.65414126573789 HPD interval (0.99): 1.12226048071043..8.90852485732572 HPD interval (0.999): 1.12226048071043..8.90852485732572 HPD interval (0.9999): 1.12226048071043..8.90852485732572 HPD interval (0.99999): 1.12226048071043..8.90852485732572 var : post Probabilities (truncated): 59.578763172989923: 0.0151515151515152 34.327655197029642: 0.0151515151515152 26.970524502106827: 0.0151515151515152 25.647462809837371: 0.0151515151515152 ......... 3.852857913261569: 0.0151515151515152 3.168561512480159: 0.0151515151515152 1.530042821027592: 0.0151515151515152 -2.00468169208193: 0.0151515151515152 mean = 14.4123 HPD intervals: HPD interval (0.5): 7.10247673410081..14.56975245137271 HPD interval (0.84): 6.49538184875579..21.79813399774052 HPD interval (0.9): 6.49538184875579..26.97052450210683 HPD interval (0.95): 1.53004282102759..26.97052450210683 HPD interval (0.99): -2.00468169208193..59.57876317298992 HPD interval (0.999): -2.00468169208193..59.57876317298992 HPD interval (0.9999): -2.00468169208193..59.57876317298992 HPD interval (0.99999): -2.00468169208193..59.57876317298992 */ go3 ?=> reset_store, run_model(100_000,$model3,[show_probs_trunc,mean,show_hpd_intervals]), nl, % show_store_lengths, % fail, nl. go3 => true. model3() => % X = [15.281819132444104,10.028270196182174,14.014207692800706,9.901431868321195,22.504941544790793,35.934094047899976,10.887493952334752,5.223940870946778,22.304607209505008,0.61166403311217], X = [15.281819132444104,10.028270196182174,14.014207692800706], Alpha = uniform(0,20), Beta = uniform(0,10), foreach(I in 1..X.len) observe( abs(extreme_value_dist(Alpha,Beta)-X[I]) <= 1.0) % observe( abs( normal_dist(extreme_value_dist(Alpha,Beta),10)-X[I]) <= 1.0) end, Post = extreme_value_dist(Alpha,Beta), if observed_ok then % println(Alpha=Beta=Post), add("alpha",Alpha), add("beta",Beta), add("post",Post) end.