/* Bayesian Linear regression in Picat. This is a port of my Racket/Gamble model gamble_bayesian_linear_regression.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. main => go. /* var : b Probabilities (truncated): 3.066889684964614: 0.0028409090909091 2.590439185100948: 0.0028409090909091 2.453108343658046: 0.0028409090909091 2.318288390271639: 0.0028409090909091 ......... -2.110609110212689: 0.0028409090909091 -2.248233409230633: 0.0028409090909091 -2.769346240912291: 0.0028409090909091 -3.648259134321716: 0.0028409090909091 mean = -0.225598 Percentiles: (0.001 -3.3397607087350076) (0.01 -2.1019009667521602) (0.025 -1.3342496838504587) (0.05 -1.1215808140992582) (0.1 -0.94017402858987142) (0.25 -0.66029457493296573) (0.5 -0.31660322149685183) (0.75 0.14924261533457958) (0.84 0.36770377768003565) (0.9 0.62832300007095832) (0.95 0.99263422057932671) (0.975 1.2783546228404217) (0.99 2.070048706878588) (0.999 2.8996555595124676) (0.9999 3.050166272419399) (0.99999 3.0652173437100894) HPD intervals: HPD interval (0.5): -0.72454423762679..0.01179705214185 HPD interval (0.84): -1.12600140685322..0.50378023886626 HPD interval (0.9): -1.14597587947581..0.89982210711531 HPD interval (0.95): -1.22007613657002..1.30991837952125 HPD interval (0.99): -2.24823340923063..2.59043918510095 HPD interval (0.999): -3.64825913432172..3.06688968496461 HPD interval (0.9999): -3.64825913432172..3.06688968496461 HPD interval (0.99999): -3.64825913432172..3.06688968496461 Histogram: -3.648: 1 # (0.003 / 0.003) -3.480: 0 (0.000 / 0.003) -3.313: 0 (0.000 / 0.003) -3.145: 0 (0.000 / 0.003) -2.977: 0 (0.000 / 0.003) -2.809: 1 # (0.003 / 0.006) -2.641: 0 (0.000 / 0.006) -2.473: 0 (0.000 / 0.006) -2.305: 1 # (0.003 / 0.009) -2.137: 2 ## (0.006 / 0.014) -1.969: 0 (0.000 / 0.014) -1.802: 0 (0.000 / 0.014) -1.634: 1 # (0.003 / 0.017) -1.466: 2 ## (0.006 / 0.023) -1.298: 3 #### (0.009 / 0.031) -1.130: 13 ################ (0.037 / 0.068) -0.962: 22 ########################### (0.062 / 0.131) -0.794: 29 #################################### (0.082 / 0.213) -0.626: 39 ################################################ (0.111 / 0.324) -0.459: 49 ############################################################ (0.139 / 0.463) -0.291: 40 ################################################# (0.114 / 0.577) -0.123: 31 ###################################### (0.088 / 0.665) 0.045: 26 ################################ (0.074 / 0.739) 0.213: 23 ############################ (0.065 / 0.804) 0.381: 22 ########################### (0.062 / 0.866) 0.549: 13 ################ (0.037 / 0.903) 0.717: 7 ######### (0.020 / 0.923) 0.884: 8 ########## (0.023 / 0.946) 1.052: 5 ###### (0.014 / 0.960) 1.220: 5 ###### (0.014 / 0.974) 1.388: 4 ##### (0.011 / 0.986) 1.556: 0 (0.000 / 0.986) 1.724: 0 (0.000 / 0.986) 1.892: 1 # (0.003 / 0.989) 2.060: 0 (0.000 / 0.989) 2.227: 0 (0.000 / 0.989) 2.395: 2 ## (0.006 / 0.994) 2.563: 1 # (0.003 / 0.997) 2.731: 0 (0.000 / 0.997) 2.899: 1 # (0.003 / 1.000) var : c Probabilities: 0: 0.4439750000000000 1: 0.3847500000000000 2: 0.1317000000000000 3: 0.0307750000000000 4: 0.0088000000000000 mean = 0.775675 Percentiles: (0.001 0) (0.01 0) (0.025 0) (0.05 0) (0.1 0) (0.25 0) (0.5 1) (0.75 1) (0.84 2) (0.9 2) (0.95 2) (0.975 3) (0.99 3) (0.999 4) (0.9999 4) (0.99999 4) HPD intervals: HPD interval (0.5): 0.00000000000000..1.00000000000000 HPD interval (0.84): 0.00000000000000..2.00000000000000 HPD interval (0.9): 0.00000000000000..2.00000000000000 HPD interval (0.95): 0.00000000000000..2.00000000000000 HPD interval (0.99): 0.00000000000000..3.00000000000000 HPD interval (0.999): 0.00000000000000..4.00000000000000 HPD interval (0.9999): 0.00000000000000..4.00000000000000 HPD interval (0.99999): 0.00000000000000..4.00000000000000 Histogram: 0: 17759 ############################################################ (0.444 / 0.444) 1: 15390 #################################################### (0.385 / 0.829) 2: 5268 ################## (0.132 / 0.960) 3: 1231 #### (0.031 / 0.991) 4: 352 # (0.009 / 1.000) var : m Probabilities (truncated): 3.584300216864791: 0.0028409090909091 3.291515139218108: 0.0028409090909091 3.223448205365786: 0.0028409090909091 3.204078215345814: 0.0028409090909091 ......... 0.670058399697084: 0.0028409090909091 0.579633597561929: 0.0028409090909091 0.24141567123122: 0.0028409090909091 -0.346450221490876: 0.0028409090909091 mean = 1.98486 Percentiles: (0.001 -0.14010929314542034) (0.01 0.69603064183869434) (0.025 0.97893076034891913) (0.05 1.3682593312539222) (0.1 1.5610920785438347) (0.25 1.7833632578715681) (0.5 1.9936113672152884) (0.75 2.2047824569147787) (0.84 2.3301925996961019) (0.9 2.4545959039727601) (0.95 2.552064313546492) (0.975 2.6438856336545271) (0.99 3.1944388764589462) (0.999 3.4815326546108056) (0.9999 3.5740234606393928) (0.99999 3.5832725412422497) HPD intervals: HPD interval (0.5): 1.75729560799170..2.16369087969785 HPD interval (0.84): 1.53023135919329..2.48091855411771 HPD interval (0.9): 1.53023135919329..2.64462630628488 HPD interval (0.95): 1.12708937015546..2.65659141675684 HPD interval (0.99): 0.57963359756193..3.29151513921811 HPD interval (0.999): -0.34645022149088..3.58430021686479 HPD interval (0.9999): -0.34645022149088..3.58430021686479 HPD interval (0.99999): -0.34645022149088..3.58430021686479 Histogram: -0.346: 1 # (0.003 / 0.003) -0.248: 0 (0.000 / 0.003) -0.150: 0 (0.000 / 0.003) -0.052: 0 (0.000 / 0.003) 0.047: 0 (0.000 / 0.003) 0.145: 0 (0.000 / 0.003) 0.243: 1 # (0.003 / 0.006) 0.341: 0 (0.000 / 0.006) 0.440: 0 (0.000 / 0.006) 0.538: 1 # (0.003 / 0.009) 0.636: 1 # (0.003 / 0.011) 0.735: 1 # (0.003 / 0.014) 0.833: 1 # (0.003 / 0.017) 0.931: 3 #### (0.009 / 0.026) 1.029: 1 # (0.003 / 0.028) 1.128: 2 ### (0.006 / 0.034) 1.226: 2 ### (0.006 / 0.040) 1.324: 4 ##### (0.011 / 0.051) 1.422: 4 ##### (0.011 / 0.063) 1.521: 15 #################### (0.043 / 0.105) 1.619: 21 ########################### (0.060 / 0.165) 1.717: 21 ########################### (0.060 / 0.224) 1.815: 40 #################################################### (0.114 / 0.338) 1.914: 46 ############################################################ (0.131 / 0.469) 2.012: 37 ################################################ (0.105 / 0.574) 2.110: 42 ####################################################### (0.119 / 0.693) 2.209: 35 ############################################## (0.099 / 0.793) 2.307: 21 ########################### (0.060 / 0.852) 2.405: 16 ##################### (0.045 / 0.898) 2.503: 18 ####################### (0.051 / 0.949) 2.602: 10 ############# (0.028 / 0.977) 2.700: 2 ### (0.006 / 0.983) 2.798: 1 # (0.003 / 0.986) 2.896: 0 (0.000 / 0.986) 2.995: 0 (0.000 / 0.986) 3.093: 0 (0.000 / 0.986) 3.191: 3 #### (0.009 / 0.994) 3.289: 1 # (0.003 / 0.997) 3.388: 0 (0.000 / 0.997) 3.486: 1 # (0.003 / 1.000) var : sigma Probabilities (truncated): 2.799426808865843: 0.0028409090909091 2.741399025167206: 0.0028409090909091 2.39325275886671: 0.0028409090909091 1.996057090686863: 0.0028409090909091 ......... 0.009012481431094: 0.0028409090909091 0.005202190963521: 0.0028409090909091 0.00118857218553: 0.0028409090909091 0.000878108199589: 0.0028409090909091 mean = 0.517874 Percentiles: (0.001 0.00098708105865449479) (0.01 0.0092675803464122401) (0.025 0.02033877767014863) (0.05 0.037436297031799751) (0.1 0.071702395898694385) (0.25 0.19794970022115738) (0.5 0.40146315884087347) (0.75 0.72179743221722048) (0.84 0.93307650096885486) (0.9 1.1302275893431502) (0.95 1.3902338505924352) (0.975 1.5340113983226997) (0.99 1.8920835661956725) (0.999 2.7790590567876214) (0.9999 2.7973900336580204) (0.99999 2.7992231313450602) HPD intervals: HPD interval (0.5): 0.06087862944222..0.45898117408151 HPD interval (0.84): 0.00087810819959..0.93401744465880 HPD interval (0.9): 0.00087810819959..1.13155860404723 HPD interval (0.95): 0.00087810819959..1.40011360112504 HPD interval (0.99): 0.00087810819959..1.99605709068686 HPD interval (0.999): 0.00087810819959..2.79942680886584 HPD interval (0.9999): 0.00087810819959..2.79942680886584 HPD interval (0.99999): 0.00087810819959..2.79942680886584 Histogram: 0.001: 18 ############################### (0.051 / 0.051) 0.071: 35 ############################################################ (0.099 / 0.151) 0.141: 29 ################################################## (0.082 / 0.233) 0.211: 29 ################################################## (0.082 / 0.315) 0.281: 35 ############################################################ (0.099 / 0.415) 0.351: 26 ############################################# (0.074 / 0.489) 0.421: 28 ################################################ (0.080 / 0.568) 0.491: 18 ############################### (0.051 / 0.619) 0.561: 16 ########################### (0.045 / 0.665) 0.631: 13 ###################### (0.037 / 0.702) 0.701: 24 ######################################### (0.068 / 0.770) 0.770: 12 ##################### (0.034 / 0.804) 0.840: 5 ######### (0.014 / 0.818) 0.910: 8 ############## (0.023 / 0.841) 0.980: 9 ############### (0.026 / 0.866) 1.050: 9 ############### (0.026 / 0.892) 1.120: 5 ######### (0.014 / 0.906) 1.190: 6 ########## (0.017 / 0.923) 1.260: 4 ####### (0.011 / 0.935) 1.330: 4 ####### (0.011 / 0.946) 1.400: 4 ####### (0.011 / 0.957) 1.470: 4 ####### (0.011 / 0.969) 1.540: 2 ### (0.006 / 0.974) 1.610: 1 ## (0.003 / 0.977) 1.680: 3 ##### (0.009 / 0.986) 1.750: 0 (0.000 / 0.986) 1.820: 1 ## (0.003 / 0.989) 1.890: 0 (0.000 / 0.989) 1.960: 0 (0.000 / 0.989) 2.030: 1 ## (0.003 / 0.991) 2.100: 0 (0.000 / 0.991) 2.170: 0 (0.000 / 0.991) 2.240: 0 (0.000 / 0.991) 2.310: 0 (0.000 / 0.991) 2.380: 1 ## (0.003 / 0.994) 2.450: 0 (0.000 / 0.994) 2.520: 0 (0.000 / 0.994) 2.590: 0 (0.000 / 0.994) 2.659: 0 (0.000 / 0.994) 2.729: 2 ### (0.006 / 1.000) var : y4 Probabilities (truncated): 8.808831800928637: 0.0028409090909091 8.689391739947862: 0.0028409090909091 8.181990335282649: 0.0028409090909091 8.12546163501238: 0.0028409090909091 ......... 2.573213883313889: 0.0028409090909091 1.579227794218249: 0.0028409090909091 1.20316883382541: 0.0028409090909091 0.435318039806809: 0.0028409090909091 mean = 5.70581 Percentiles: (0.001 0.70483366850733775) (0.01 2.8656485276688102) (0.025 3.5536721340879334) (0.05 4.086949413741241) (0.1 4.6763381221592688) (0.25 5.2020569272704389) (0.5 5.7004693115907346) (0.75 6.2376339039692583) (0.84 6.582692369831566) (0.9 6.8597131953832369) (0.95 7.1909749935355274) (0.975 7.6874862500085879) (0.99 8.0275370599286813) (0.999 8.7669083395243845) (0.9999 8.8046394547882123) (0.99999 8.808412566314594) HPD intervals: HPD interval (0.5): 5.16886068371922..6.18747746870571 HPD interval (0.84): 4.67305501405337..6.97360528242417 HPD interval (0.9): 4.01598104889658..7.04674844838890 HPD interval (0.95): 3.72089560427775..7.78122671787557 HPD interval (0.99): 2.57321388331389..8.80883180092864 HPD interval (0.999): 0.43531803980681..8.80883180092864 HPD interval (0.9999): 0.43531803980681..8.80883180092864 HPD interval (0.99999): 0.43531803980681..8.80883180092864 Histogram: 0.435: 1 # (0.003 / 0.003) 0.645: 0 (0.000 / 0.003) 0.854: 0 (0.000 / 0.003) 1.063: 0 (0.000 / 0.003) 1.273: 1 # (0.003 / 0.006) 1.482: 1 # (0.003 / 0.009) 1.691: 0 (0.000 / 0.009) 1.901: 0 (0.000 / 0.009) 2.110: 0 (0.000 / 0.009) 2.319: 0 (0.000 / 0.009) 2.529: 1 # (0.003 / 0.011) 2.738: 0 (0.000 / 0.011) 2.947: 0 (0.000 / 0.011) 3.157: 1 # (0.003 / 0.014) 3.366: 3 #### (0.009 / 0.023) 3.575: 3 #### (0.009 / 0.031) 3.785: 3 #### (0.009 / 0.040) 3.994: 6 ######## (0.017 / 0.057) 4.203: 5 ####### (0.014 / 0.071) 4.413: 4 ##### (0.011 / 0.082) 4.622: 9 ############ (0.026 / 0.108) 4.831: 14 ################### (0.040 / 0.148) 5.041: 28 ###################################### (0.080 / 0.227) 5.250: 28 ###################################### (0.080 / 0.307) 5.459: 44 ############################################################ (0.125 / 0.432) 5.669: 34 ############################################## (0.097 / 0.528) 5.878: 35 ################################################ (0.099 / 0.628) 6.087: 38 #################################################### (0.108 / 0.736) 6.297: 22 ############################## (0.062 / 0.798) 6.506: 18 ######################### (0.051 / 0.849) 6.715: 17 ####################### (0.048 / 0.898) 6.925: 15 #################### (0.043 / 0.940) 7.134: 6 ######## (0.017 / 0.957) 7.343: 1 # (0.003 / 0.960) 7.553: 4 ##### (0.011 / 0.972) 7.762: 4 ##### (0.011 / 0.983) 7.971: 2 ### (0.006 / 0.989) 8.181: 2 ### (0.006 / 0.994) 8.390: 0 (0.000 / 0.994) 8.599: 2 ### (0.006 / 1.000) */ go ?=> Ys = [0,1,4,6], run_model(40_000,$model(Ys),[show_probs_trunc,mean,show_percentiles,show_hpd_intervals,show_histogram]), % run_model(30_000,$model(Ys),[mean]), % fail, nl. go => true. model(Ys) => Len = Ys.len, M = normal_dist(0,2), B = normal_dist(0,2), Sigma = gamma_dist(1,1), Y = [normal_dist(M * (I-1) + B,Sigma) : I in 1..Len], C = 0, foreach(I in 1..Len) if abs(Ys[I]-Y[I]) <= 1.0 then C := C + 1 end end, observe(C == Len), Y4 = normal_dist(M * (4-1) + B,Sigma), add("c",C), if observed_ok then add("m",M), add("b",B), add("sigma",Sigma), add("y4",Y4), end.