/* 8 schools in Picat. Cf my Gamble model gamble_8_schools.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. /* Forcing 1000 accepted samples Num accepted samples: 995 Total samples: 1300 (0.765%) Num accepted samples: 996 Total samples: 1301 (0.766%) Num accepted samples: 997 Total samples: 1302 (0.766%) Num accepted samples: 998 Total samples: 1303 (0.766%) Num accepted samples: 999 Total samples: 1304 (0.766%) Num accepted samples: 1000 Total samples: 1305 (0.766%) var : tau Probabilities (truncated): 15.265756971936828: 0.0010000000000000 14.185762316017307: 0.0010000000000000 14.040961066126583: 0.0010000000000000 13.724223869226909: 0.0010000000000000 ......... 6.056032114296944: 0.0010000000000000 5.737246332767662: 0.0010000000000000 5.428844572682983: 0.0010000000000000 4.349959113983455: 0.0010000000000000 mean = 9.83543 HPD intervals: HPD interval (0.84): 7.99064707948034..12.17097490862742 Histogram (total 1000) 4.350: 1 # (0.001 / 0.001) 4.623: 0 (0.000 / 0.001) 4.896: 0 (0.000 / 0.001) 5.169: 0 (0.000 / 0.001) 5.442: 1 # (0.001 / 0.002) 5.714: 1 # (0.001 / 0.003) 5.987: 2 # (0.002 / 0.005) 6.260: 3 ## (0.003 / 0.008) 6.533: 6 #### (0.006 / 0.014) 6.806: 9 ###### (0.009 / 0.023) 7.079: 13 ######### (0.013 / 0.036) 7.352: 14 ########## (0.014 / 0.050) 7.625: 29 #################### (0.029 / 0.079) 7.898: 28 ################### (0.028 / 0.107) 8.170: 44 ############################## (0.044 / 0.151) 8.443: 49 ################################# (0.049 / 0.200) 8.716: 61 ########################################## (0.061 / 0.261) 8.989: 61 ########################################## (0.061 / 0.322) 9.262: 69 ############################################### (0.069 / 0.391) 9.535: 88 ############################################################ (0.088 / 0.479) 9.808: 61 ########################################## (0.061 / 0.540) 10.081: 68 ############################################## (0.068 / 0.608) 10.354: 68 ############################################## (0.068 / 0.676) 10.627: 59 ######################################## (0.059 / 0.735) 10.899: 49 ################################# (0.049 / 0.784) 11.172: 52 ################################### (0.052 / 0.836) 11.445: 37 ######################### (0.037 / 0.873) 11.718: 29 #################### (0.029 / 0.902) 11.991: 31 ##################### (0.031 / 0.933) 12.264: 24 ################ (0.024 / 0.957) 12.537: 13 ######### (0.013 / 0.970) 12.810: 7 ##### (0.007 / 0.977) 13.083: 9 ###### (0.009 / 0.986) 13.355: 6 #### (0.006 / 0.992) 13.628: 5 ### (0.005 / 0.997) 13.901: 0 (0.000 / 0.997) 14.174: 2 # (0.002 / 0.999) 14.447: 0 (0.000 / 0.999) 14.720: 0 (0.000 / 0.999) 14.993: 1 # (0.001 / 1.000) var : mu Probabilities (truncated): 14.333919603486269: 0.0010000000000000 14.074402228424436: 0.0010000000000000 13.732774108627201: 0.0010000000000000 13.504448205596358: 0.0010000000000000 ......... 6.241419943131531: 0.0010000000000000 6.01005238286889: 0.0010000000000000 5.783640757179354: 0.0010000000000000 5.739027455921886: 0.0010000000000000 mean = 9.91927 HPD intervals: HPD interval (0.84): 8.16085591041359..11.91388514782323 Histogram (total 1000) 5.739: 2 ## (0.002 / 0.002) 5.954: 1 # (0.001 / 0.003) 6.169: 1 # (0.001 / 0.004) 6.384: 2 ## (0.002 / 0.006) 6.599: 2 ## (0.002 / 0.008) 6.813: 3 ### (0.003 / 0.011) 7.028: 5 #### (0.005 / 0.016) 7.243: 13 ########### (0.013 / 0.029) 7.458: 17 ############### (0.017 / 0.046) 7.673: 12 ########## (0.012 / 0.058) 7.888: 13 ########### (0.013 / 0.071) 8.103: 25 ###################### (0.025 / 0.096) 8.317: 33 ############################# (0.033 / 0.129) 8.532: 34 ############################## (0.034 / 0.163) 8.747: 48 ########################################## (0.048 / 0.211) 8.962: 58 ################################################## (0.058 / 0.269) 9.177: 56 ################################################# (0.056 / 0.325) 9.392: 69 ############################################################ (0.069 / 0.394) 9.607: 68 ########################################################### (0.068 / 0.462) 9.822: 56 ################################################# (0.056 / 0.518) 10.036: 52 ############################################# (0.052 / 0.570) 10.251: 56 ################################################# (0.056 / 0.626) 10.466: 54 ############################################### (0.054 / 0.680) 10.681: 48 ########################################## (0.048 / 0.728) 10.896: 46 ######################################## (0.046 / 0.774) 11.111: 44 ###################################### (0.044 / 0.818) 11.326: 51 ############################################ (0.051 / 0.869) 11.541: 24 ##################### (0.024 / 0.893) 11.755: 24 ##################### (0.024 / 0.917) 11.970: 27 ####################### (0.027 / 0.944) 12.185: 14 ############ (0.014 / 0.958) 12.400: 19 ################# (0.019 / 0.977) 12.615: 7 ###### (0.007 / 0.984) 12.830: 5 #### (0.005 / 0.989) 13.045: 5 #### (0.005 / 0.994) 13.260: 1 # (0.001 / 0.995) 13.474: 2 ## (0.002 / 0.997) 13.689: 1 # (0.001 / 0.998) 13.904: 0 (0.000 / 0.998) 14.119: 2 ## (0.002 / 1.000) var : eta[1] Probabilities (truncated): 3.577068890504524: 0.0010000000000000 2.700787001954348: 0.0010000000000000 2.697346015770439: 0.0010000000000000 2.57947466793492: 0.0010000000000000 ......... -2.541534595655675: 0.0010000000000000 -2.710624066064451: 0.0010000000000000 -3.421055806150165: 0.0010000000000000 -3.477373484611852: 0.0010000000000000 mean = -0.00988128 HPD intervals: HPD interval (0.84): -1.31818717224933..1.35811603239546 Histogram (total 1000) -3.477: 2 ## (0.002 / 0.002) -3.301: 0 (0.000 / 0.002) -3.125: 0 (0.000 / 0.002) -2.948: 0 (0.000 / 0.002) -2.772: 1 # (0.001 / 0.003) -2.596: 2 ## (0.002 / 0.005) -2.419: 5 #### (0.005 / 0.010) -2.243: 3 ## (0.003 / 0.013) -2.066: 6 ##### (0.006 / 0.019) -1.890: 13 ########### (0.013 / 0.032) -1.714: 13 ########### (0.013 / 0.045) -1.537: 21 ################# (0.021 / 0.066) -1.361: 26 ##################### (0.026 / 0.092) -1.185: 32 ########################## (0.032 / 0.124) -1.008: 31 ######################### (0.031 / 0.155) -0.832: 68 ####################################################### (0.068 / 0.223) -0.656: 60 ################################################# (0.060 / 0.283) -0.479: 74 ############################################################ (0.074 / 0.357) -0.303: 71 ########################################################## (0.071 / 0.428) -0.127: 72 ########################################################## (0.072 / 0.500) 0.050: 67 ###################################################### (0.067 / 0.567) 0.226: 66 ###################################################### (0.066 / 0.633) 0.403: 65 ##################################################### (0.065 / 0.698) 0.579: 72 ########################################################## (0.072 / 0.770) 0.755: 51 ######################################### (0.051 / 0.821) 0.932: 29 ######################## (0.029 / 0.850) 1.108: 36 ############################# (0.036 / 0.886) 1.284: 36 ############################# (0.036 / 0.922) 1.461: 24 ################### (0.024 / 0.946) 1.637: 19 ############### (0.019 / 0.965) 1.813: 13 ########### (0.013 / 0.978) 1.990: 3 ## (0.003 / 0.981) 2.166: 9 ####### (0.009 / 0.990) 2.343: 3 ## (0.003 / 0.993) 2.519: 4 ### (0.004 / 0.997) 2.695: 2 ## (0.002 / 0.999) 2.872: 0 (0.000 / 0.999) 3.048: 0 (0.000 / 0.999) 3.224: 0 (0.000 / 0.999) 3.401: 1 # (0.001 / 1.000) var : eta mean Probabilities (truncated): 0.929310988757471: 0.0010000000000000 0.914821258863032: 0.0010000000000000 0.868061758897515: 0.0010000000000000 0.862973545593311: 0.0010000000000000 ......... -0.880244436778075: 0.0010000000000000 -0.889137524316275: 0.0010000000000000 -1.043052867319072: 0.0010000000000000 -1.09571170981086: 0.0010000000000000 mean = -0.0119344 HPD intervals: HPD interval (0.84): -0.47961330971879..0.38678503749665 Histogram (total 1000) -1.096: 1 # (0.001 / 0.001) -1.045: 1 # (0.001 / 0.002) -0.994: 0 (0.000 / 0.002) -0.944: 0 (0.000 / 0.002) -0.893: 2 ## (0.002 / 0.004) -0.843: 2 ## (0.002 / 0.006) -0.792: 2 ## (0.002 / 0.008) -0.741: 6 ##### (0.006 / 0.014) -0.691: 4 ### (0.004 / 0.018) -0.640: 10 ######## (0.010 / 0.028) -0.589: 11 ######### (0.011 / 0.039) -0.539: 12 ########## (0.012 / 0.051) -0.488: 17 ############## (0.017 / 0.068) -0.438: 28 ###################### (0.028 / 0.096) -0.387: 37 ############################## (0.037 / 0.133) -0.336: 43 ################################## (0.043 / 0.176) -0.286: 45 #################################### (0.045 / 0.221) -0.235: 54 ########################################### (0.054 / 0.275) -0.184: 60 ################################################ (0.060 / 0.335) -0.134: 60 ################################################ (0.060 / 0.395) -0.083: 57 ############################################## (0.057 / 0.452) -0.033: 52 ########################################## (0.052 / 0.504) 0.018: 60 ################################################ (0.060 / 0.564) 0.069: 75 ############################################################ (0.075 / 0.639) 0.119: 52 ########################################## (0.052 / 0.691) 0.170: 51 ######################################### (0.051 / 0.742) 0.221: 41 ################################# (0.041 / 0.783) 0.271: 48 ###################################### (0.048 / 0.831) 0.322: 44 ################################### (0.044 / 0.875) 0.372: 30 ######################## (0.030 / 0.905) 0.423: 18 ############## (0.018 / 0.923) 0.474: 21 ################# (0.021 / 0.944) 0.524: 8 ###### (0.008 / 0.952) 0.575: 19 ############### (0.019 / 0.971) 0.626: 9 ####### (0.009 / 0.980) 0.676: 3 ## (0.003 / 0.983) 0.727: 8 ###### (0.008 / 0.991) 0.777: 1 # (0.001 / 0.992) 0.828: 2 ## (0.002 / 0.994) 0.879: 6 ##### (0.006 / 1.000) var : theta[1] Probabilities (truncated): 42.495347124123256: 0.0010000000000000 40.177924495551991: 0.0010000000000000 38.143360214274949: 0.0010000000000000 37.969383153463788: 0.0010000000000000 ......... -15.054196655408308: 0.0010000000000000 -17.11111409976035: 0.0010000000000000 -20.451920628419153: 0.0010000000000000 -20.959489971144716: 0.0010000000000000 mean = 9.84895 HPD intervals: HPD interval (0.84): -3.46200916053440..23.13879862383114 Histogram (total 1000) -20.959: 2 # (0.002 / 0.002) -19.373: 0 (0.000 / 0.002) -17.787: 1 # (0.001 / 0.003) -16.200: 0 (0.000 / 0.003) -14.614: 5 #### (0.005 / 0.008) -13.028: 4 ### (0.004 / 0.012) -11.441: 3 ## (0.003 / 0.015) -9.855: 6 #### (0.006 / 0.021) -8.269: 16 ############ (0.016 / 0.037) -6.682: 16 ############ (0.016 / 0.053) -5.096: 18 ############# (0.018 / 0.071) -3.509: 21 ################ (0.021 / 0.092) -1.923: 31 ####################### (0.031 / 0.123) -0.337: 37 ########################### (0.037 / 0.160) 1.250: 37 ########################### (0.037 / 0.197) 2.836: 58 ########################################### (0.058 / 0.255) 4.422: 62 ############################################## (0.062 / 0.317) 6.009: 62 ############################################## (0.062 / 0.379) 7.595: 62 ############################################## (0.062 / 0.441) 9.182: 81 ############################################################ (0.081 / 0.522) 10.768: 52 ####################################### (0.052 / 0.574) 12.354: 63 ############################################### (0.063 / 0.637) 13.941: 72 ##################################################### (0.072 / 0.709) 15.527: 55 ######################################### (0.055 / 0.764) 17.113: 46 ################################## (0.046 / 0.810) 18.700: 35 ########################## (0.035 / 0.845) 20.286: 35 ########################## (0.035 / 0.880) 21.873: 24 ################## (0.024 / 0.904) 23.459: 32 ######################## (0.032 / 0.936) 25.045: 16 ############ (0.016 / 0.952) 26.632: 7 ##### (0.007 / 0.959) 28.218: 13 ########## (0.013 / 0.972) 29.804: 9 ####### (0.009 / 0.981) 31.391: 4 ### (0.004 / 0.985) 32.977: 3 ## (0.003 / 0.988) 34.563: 3 ## (0.003 / 0.991) 36.150: 4 ### (0.004 / 0.995) 37.736: 3 ## (0.003 / 0.998) 39.323: 0 (0.000 / 0.998) 40.909: 2 # (0.002 / 1.000) var : theta mean Probabilities (truncated): 20.854066353556306: 0.0010000000000000 20.582185653682828: 0.0010000000000000 19.52912596900314: 0.0010000000000000 19.289168092022571: 0.0010000000000000 ......... 0.61140104104705: 0.0010000000000000 0.562628346052259: 0.0010000000000000 -0.302481731088621: 0.0010000000000000 -0.364162622956823: 0.0010000000000000 mean = 9.8358 HPD intervals: HPD interval (0.84): 5.38056581671003..14.42246170376671 Histogram (total 1000) -0.364: 2 ## (0.002 / 0.002) 0.166: 0 (0.000 / 0.002) 0.697: 4 ### (0.004 / 0.006) 1.227: 3 ### (0.003 / 0.009) 1.758: 1 # (0.001 / 0.010) 2.288: 5 #### (0.005 / 0.015) 2.819: 4 ### (0.004 / 0.019) 3.349: 5 #### (0.005 / 0.024) 3.879: 9 ######## (0.009 / 0.033) 4.410: 17 ############### (0.017 / 0.050) 4.940: 21 ################## (0.021 / 0.071) 5.471: 27 ####################### (0.027 / 0.098) 6.001: 38 ################################# (0.038 / 0.136) 6.532: 40 ################################### (0.040 / 0.176) 7.062: 54 ############################################### (0.054 / 0.230) 7.593: 47 ######################################### (0.047 / 0.277) 8.123: 69 ############################################################ (0.069 / 0.346) 8.654: 61 ##################################################### (0.061 / 0.407) 9.184: 61 ##################################################### (0.061 / 0.468) 9.714: 66 ######################################################### (0.066 / 0.534) 10.245: 59 ################################################### (0.059 / 0.593) 10.775: 50 ########################################### (0.050 / 0.643) 11.306: 68 ########################################################### (0.068 / 0.711) 11.836: 38 ################################# (0.038 / 0.749) 12.367: 52 ############################################# (0.052 / 0.801) 12.897: 35 ############################## (0.035 / 0.836) 13.428: 37 ################################ (0.037 / 0.873) 13.958: 32 ############################ (0.032 / 0.905) 14.489: 26 ####################### (0.026 / 0.931) 15.019: 20 ################# (0.020 / 0.951) 15.550: 9 ######## (0.009 / 0.960) 16.080: 9 ######## (0.009 / 0.969) 16.610: 8 ####### (0.008 / 0.977) 17.141: 8 ####### (0.008 / 0.985) 17.671: 5 #### (0.005 / 0.990) 18.202: 1 # (0.001 / 0.991) 18.732: 2 ## (0.002 / 0.993) 19.263: 4 ### (0.004 / 0.997) 19.793: 1 # (0.001 / 0.998) 20.324: 2 ## (0.002 / 1.000) */ go ?=> reset_store, run_model(100_000,$model,[show_probs_trunc,mean, show_hpd_intervals,hpd_intervals=[0.84], show_histogram, min_accepted_samples=1000,show_accepted_samples=true]), nl, % fail, nl. go => true. model() => Ys = [28.0,8.0,-3.0,7.0,-1.0,1.0,18.0,12.0], Sigmas = [15.0,10.0,16.0,11.0,9.0,11.0,10.0,18.0], N = Ys.len, Mu = normal_dist(10,sqrt(2)), % Mean is 10 Tau = normal_dist(10,sqrt(2)), % (It should be neat with a truncated_normal_dist) observe(Tau > 1), Eta = [normal_dist(0,1) : _ in 1..N], EtaMean = Eta.mean, Theta = [Mu + Tau * Eta[I] : I in 1..N], ThetaMean = Theta.mean, Y = [normal_dist(Theta[I], Sigmas[I]) : I in 1..N], % likelihood % This variant works in Gamble, but is too slow in Picat PPL % foreach(I in 1..N) % observe(abs(Ys[I]-Y[I]) <= 5) % end, % This approach is much faster Mean = Ys.mean, Stdev = Ys.stdev, YMean = Y.mean, YStdev = Y.stdev, observe( abs(Mean-YMean) < Stdev ), observe( abs(Stdev-YStdev) < Stdev ), if observed_ok then add("tau",Tau), add("mu",Mu), add("eta[1]",Eta[1]), add("eta mean",EtaMean), add("theta[1]",Theta[1]), add("theta mean",ThetaMean), end.