/* Heights in Picat. From "Probabilistic programming", John Winn and Tom Minka, Machine Learning Summer School (ProbabilisticProgramming_slides.pdf) """ Suppose we take a woman at random and a man at random from the UK population The woman turns out to be taller than the man. - What is the probability of this event? - What is the posterior distribution over the woman’s height? - What is the posterior distribution over the man’s height? .... double heightMan = random((Normal 177,64)); double heightWoman = random((Normal 164,64)); Bernoulli dist = infer(heightWoman > heightMan); constrain(heightWoman > heightMan); Gaussian distWoman = infer(heightWoman); Gaussian distMan = infer(heightMan); """ I have split this into two questions: - model with obs #f: What is the probability? - model with obs #t: Posterior distributions (we observe that the woman is taller than the man) Cf my Gamble model gamble_heights.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. /* Compiling:: /home/hakank/picat/me/run_program.pi run_program.pi compiled in 2 milliseconds loading... program = ppl_heights.pi Compiling:: ppl_heights.pi ppl_heights.pi compiled in 3 milliseconds loading... loading.../home/hakank/picat/me/ppl_distributions.qi loading.../home/hakank/picat/me/ppl_utils.qi goal = go obs = false var : height female Probabilities (truncated): 196.073713969843396: 0.0001000000000000 192.954183121136396: 0.0001000000000000 190.92321567948639: 0.0001000000000000 190.742908449887068: 0.0001000000000000 ......... 136.809963612525365: 0.0001000000000000 134.206277063382231: 0.0001000000000000 133.535795110348204: 0.0001000000000000 131.92504599987646: 0.0001000000000000 mean = 163.997 HPD intervals: HPD interval (0.84): 152.13877247385190..174.54027293715163 HPD interval (0.95): 147.67209885761076..178.64531347452481 HPD interval (0.99): 143.20560822208449..183.92675237115640 Histogram (total 10000) 131.925: 1 (0.000 / 0.000) 133.529: 2 (0.000 / 0.000) 135.132: 0 (0.000 / 0.000) 136.736: 2 (0.000 / 0.001) 138.340: 6 (0.001 / 0.001) 139.944: 6 (0.001 / 0.002) 141.547: 19 # (0.002 / 0.004) 143.151: 25 ## (0.003 / 0.006) 144.755: 28 ## (0.003 / 0.009) 146.358: 69 ##### (0.007 / 0.016) 147.962: 127 ######### (0.013 / 0.029) 149.566: 146 ########### (0.015 / 0.043) 151.170: 221 ################ (0.022 / 0.065) 152.773: 297 ###################### (0.030 / 0.095) 154.377: 387 ############################ (0.039 / 0.134) 155.981: 518 ###################################### (0.052 / 0.185) 157.585: 553 ######################################### (0.055 / 0.241) 159.188: 682 ################################################## (0.068 / 0.309) 160.792: 732 ###################################################### (0.073 / 0.382) 162.396: 815 ############################################################ (0.082 / 0.464) 163.999: 780 ######################################################### (0.078 / 0.542) 165.603: 753 ####################################################### (0.075 / 0.617) 167.207: 779 ######################################################### (0.078 / 0.695) 168.811: 613 ############################################# (0.061 / 0.756) 170.414: 605 ############################################# (0.060 / 0.817) 172.018: 457 ################################## (0.046 / 0.862) 173.622: 417 ############################### (0.042 / 0.904) 175.225: 293 ###################### (0.029 / 0.933) 176.829: 249 ################## (0.025 / 0.958) 178.433: 152 ########### (0.015 / 0.973) 180.037: 95 ####### (0.009 / 0.983) 181.640: 74 ##### (0.007 / 0.990) 183.244: 46 ### (0.005 / 0.995) 184.848: 16 # (0.002 / 0.996) 186.451: 20 # (0.002 / 0.998) 188.055: 8 # (0.001 / 0.999) 189.659: 3 (0.000 / 1.000) 191.263: 2 (0.000 / 1.000) 192.866: 1 (0.000 / 1.000) 194.470: 1 (0.000 / 1.000) var : height male Probabilities (truncated): 200.585507640480643: 0.0001000000000000 197.509913073356046: 0.0001000000000000 197.25201789867296: 0.0001000000000000 197.157674132415707: 0.0001000000000000 ......... 145.382165421712699: 0.0001000000000000 144.910615022722197: 0.0001000000000000 143.45905620874089: 0.0001000000000000 135.924403078982238: 0.0001000000000000 mean = 172.066 HPD intervals: HPD interval (0.84): 161.26602378748996..183.45757330512276 HPD interval (0.95): 156.47387418362140..187.46060059297639 HPD interval (0.99): 151.88107314035838..191.98410535853321 Histogram (total 10000) 135.924: 1 (0.000 / 0.000) 137.541: 0 (0.000 / 0.000) 139.157: 0 (0.000 / 0.000) 140.774: 0 (0.000 / 0.000) 142.391: 0 (0.000 / 0.000) 144.007: 1 (0.000 / 0.000) 145.624: 3 (0.000 / 0.001) 147.240: 7 (0.001 / 0.001) 148.857: 15 # (0.002 / 0.003) 150.473: 15 # (0.002 / 0.004) 152.090: 31 ## (0.003 / 0.007) 153.706: 56 #### (0.006 / 0.013) 155.323: 78 ###### (0.008 / 0.021) 156.939: 143 ########## (0.014 / 0.035) 158.556: 194 ############## (0.019 / 0.054) 160.172: 271 ################### (0.027 / 0.082) 161.789: 359 ########################## (0.036 / 0.117) 163.405: 472 ################################## (0.047 / 0.165) 165.022: 520 ##################################### (0.052 / 0.217) 166.638: 654 ############################################### (0.065 / 0.282) 168.255: 708 ################################################### (0.071 / 0.353) 169.871: 768 ####################################################### (0.077 / 0.430) 171.488: 841 ############################################################ (0.084 / 0.514) 173.105: 795 ######################################################### (0.080 / 0.593) 174.721: 743 ##################################################### (0.074 / 0.667) 176.338: 661 ############################################### (0.066 / 0.734) 177.954: 658 ############################################### (0.066 / 0.799) 179.571: 551 ####################################### (0.055 / 0.855) 181.187: 416 ############################## (0.042 / 0.896) 182.804: 324 ####################### (0.032 / 0.928) 184.420: 221 ################ (0.022 / 0.951) 186.037: 170 ############ (0.017 / 0.968) 187.653: 124 ######### (0.012 / 0.980) 189.270: 87 ###### (0.009 / 0.989) 190.886: 51 #### (0.005 / 0.994) 192.503: 31 ## (0.003 / 0.997) 194.119: 12 # (0.001 / 0.998) 195.736: 9 # (0.001 / 0.999) 197.352: 9 # (0.001 / 1.000) 198.969: 1 (0.000 / 1.000) var : diff Probabilities (truncated): 41.453928000025911: 0.0001000000000000 34.254099826019399: 0.0001000000000000 34.035056714259497: 0.0001000000000000 32.605472216670279: 0.0001000000000000 ......... -42.768204935120792: 0.0001000000000000 -46.069352393226978: 0.0001000000000000 -50.032288714260346: 0.0001000000000000 -58.738654180635336: 0.0001000000000000 mean = -8.06898 HPD intervals: HPD interval (0.84): -24.51683351853256..7.16161060066540 HPD interval (0.95): -30.16951892328865..13.84130681585020 HPD interval (0.99): -36.98677579241419..20.99865198905604 Histogram (total 10000) -58.739: 1 (0.000 / 0.000) -56.234: 0 (0.000 / 0.000) -53.729: 0 (0.000 / 0.000) -51.224: 1 (0.000 / 0.000) -48.719: 0 (0.000 / 0.000) -46.215: 1 (0.000 / 0.000) -43.710: 3 (0.000 / 0.001) -41.205: 18 # (0.002 / 0.002) -38.700: 20 # (0.002 / 0.004) -36.195: 44 ### (0.004 / 0.009) -33.691: 68 ##### (0.007 / 0.016) -31.186: 92 ###### (0.009 / 0.025) -28.681: 173 ############ (0.017 / 0.042) -26.176: 218 ############### (0.022 / 0.064) -23.671: 375 ######################### (0.037 / 0.101) -21.166: 475 ################################ (0.048 / 0.149) -18.662: 609 ######################################### (0.061 / 0.210) -16.157: 643 ############################################ (0.064 / 0.274) -13.652: 795 ###################################################### (0.080 / 0.354) -11.147: 845 ######################################################### (0.085 / 0.438) -8.642: 883 ############################################################ (0.088 / 0.526) -6.138: 871 ########################################################### (0.087 / 0.614) -3.633: 802 ###################################################### (0.080 / 0.694) -1.128: 751 ################################################### (0.075 / 0.769) 1.377: 612 ########################################## (0.061 / 0.830) 3.882: 474 ################################ (0.047 / 0.877) 6.387: 355 ######################## (0.035 / 0.913) 8.891: 330 ###################### (0.033 / 0.946) 11.396: 192 ############# (0.019 / 0.965) 13.901: 141 ########## (0.014 / 0.979) 16.406: 84 ###### (0.008 / 0.988) 18.911: 52 #### (0.005 / 0.993) 21.415: 40 ### (0.004 / 0.997) 23.920: 15 # (0.002 / 0.998) 26.425: 2 (0.000 / 0.998) 28.930: 9 # (0.001 / 0.999) 31.435: 3 (0.000 / 1.000) 33.939: 2 (0.000 / 1.000) 36.444: 0 (0.000 / 1.000) 38.949: 1 (0.000 / 1.000) var : p Probabilities: false: 0.7662000000000000 true: 0.2338000000000000 mean = [false = 0.7662,true = 0.2338] HPD intervals: show_hpd_intervals: data is not numeric Histogram (total 10000) false : 7662 ############################################################ (0.766 / 0.766) true : 2338 ################## (0.234 / 1.000) obs = true var : height female Probabilities (truncated): 192.789746690744067: 0.0004387889425186 191.963119354026929: 0.0004387889425186 191.765551904903504: 0.0004387889425186 191.700120199701587: 0.0004387889425186 ......... 153.838830435273621: 0.0004387889425186 153.108060580743341: 0.0004387889425186 151.759868400811683: 0.0004387889425186 150.920646600363767: 0.0004387889425186 mean = 171.209 HPD intervals: HPD interval (0.84): 162.48540075301108..180.85482185163323 HPD interval (0.95): 158.56633496588259..183.85766767818180 HPD interval (0.99): 155.37113293625862..187.31351508103873 Histogram (total 2279) 150.921: 1 (0.000 / 0.000) 151.967: 1 (0.000 / 0.001) 153.014: 1 (0.000 / 0.001) 154.061: 6 ## (0.003 / 0.004) 155.108: 4 # (0.002 / 0.006) 156.154: 8 ### (0.004 / 0.009) 157.201: 22 ######## (0.010 / 0.019) 158.248: 22 ######## (0.010 / 0.029) 159.294: 34 ############ (0.015 / 0.043) 160.341: 34 ############ (0.015 / 0.058) 161.388: 52 ################## (0.023 / 0.081) 162.435: 55 ################### (0.024 / 0.105) 163.481: 56 ################### (0.025 / 0.130) 164.528: 72 ######################### (0.032 / 0.161) 165.575: 114 ####################################### (0.050 / 0.211) 166.622: 109 ###################################### (0.048 / 0.259) 167.668: 121 ########################################## (0.053 / 0.312) 168.715: 142 ################################################# (0.062 / 0.375) 169.762: 174 ############################################################ (0.076 / 0.451) 170.808: 146 ################################################## (0.064 / 0.515) 171.855: 148 ################################################### (0.065 / 0.580) 172.902: 125 ########################################### (0.055 / 0.635) 173.949: 136 ############################################### (0.060 / 0.695) 174.995: 127 ############################################ (0.056 / 0.750) 176.042: 116 ######################################## (0.051 / 0.801) 177.089: 91 ############################### (0.040 / 0.841) 178.136: 76 ########################## (0.033 / 0.875) 179.182: 63 ###################### (0.028 / 0.902) 180.229: 62 ##################### (0.027 / 0.929) 181.276: 41 ############## (0.018 / 0.947) 182.322: 34 ############ (0.015 / 0.962) 183.369: 27 ######### (0.012 / 0.974) 184.416: 18 ###### (0.008 / 0.982) 185.463: 15 ##### (0.007 / 0.989) 186.509: 10 ### (0.004 / 0.993) 187.556: 5 ## (0.002 / 0.995) 188.603: 4 # (0.002 / 0.997) 189.650: 3 # (0.001 / 0.998) 190.696: 0 (0.000 / 0.998) 191.743: 4 # (0.002 / 1.000) var : height male Probabilities (truncated): 182.994528616014605: 0.0004387889425186 182.99016820029351: 0.0004387889425186 182.313025850475555: 0.0004387889425186 181.6623022114251: 0.0004387889425186 ......... 145.900986044353715: 0.0004387889425186 145.057009958680084: 0.0004387889425186 143.612649842687205: 0.0004387889425186 142.765884213856936: 0.0004387889425186 mean = 164.447 HPD intervals: HPD interval (0.84): 154.14457323673565..172.25301193341483 HPD interval (0.95): 152.13955422970105..176.95394630530288 HPD interval (0.99): 147.42630310208759..179.69599006967857 Histogram (total 2279) 142.766: 1 (0.000 / 0.000) 143.772: 1 (0.000 / 0.001) 144.777: 1 (0.000 / 0.001) 145.783: 1 (0.000 / 0.002) 146.789: 7 ### (0.003 / 0.005) 147.794: 6 ### (0.003 / 0.007) 148.800: 9 #### (0.004 / 0.011) 149.806: 12 ##### (0.005 / 0.017) 150.812: 10 #### (0.004 / 0.021) 151.817: 18 ######## (0.008 / 0.029) 152.823: 27 ############ (0.012 / 0.041) 153.829: 39 ################# (0.017 / 0.058) 154.834: 51 ###################### (0.022 / 0.080) 155.840: 64 ############################ (0.028 / 0.108) 156.846: 81 ################################### (0.036 / 0.144) 157.852: 74 ################################ (0.032 / 0.176) 158.857: 75 ################################# (0.033 / 0.209) 159.863: 112 ################################################# (0.049 / 0.258) 160.869: 121 ##################################################### (0.053 / 0.312) 161.874: 136 ############################################################ (0.060 / 0.371) 162.880: 131 ######################################################### (0.057 / 0.429) 163.886: 137 ############################################################ (0.060 / 0.489) 164.892: 135 ########################################################### (0.059 / 0.548) 165.897: 136 ############################################################ (0.060 / 0.608) 166.903: 135 ########################################################### (0.059 / 0.667) 167.909: 136 ############################################################ (0.060 / 0.727) 168.915: 127 ######################################################## (0.056 / 0.782) 169.920: 111 ################################################# (0.049 / 0.831) 170.926: 81 ################################### (0.036 / 0.867) 171.932: 67 ############################# (0.029 / 0.896) 172.937: 52 ####################### (0.023 / 0.919) 173.943: 51 ###################### (0.022 / 0.941) 174.949: 47 ##################### (0.021 / 0.962) 175.955: 21 ######### (0.009 / 0.971) 176.960: 19 ######## (0.008 / 0.979) 177.966: 18 ######## (0.008 / 0.987) 178.972: 14 ###### (0.006 / 0.993) 179.977: 6 ### (0.003 / 0.996) 180.983: 4 ## (0.002 / 0.998) 181.989: 5 ## (0.002 / 1.000) var : diff Probabilities (truncated): 34.004368152505435: 0.0004387889425186 30.906093443632727: 0.0004387889425186 28.841708234131858: 0.0004387889425186 27.808862767270455: 0.0004387889425186 ......... 0.013438749506349: 0.0004387889425186 0.009540542306098: 0.0004387889425186 0.009451062347068: 0.0004387889425186 0.005792465202319: 0.0004387889425186 mean = 6.76249 HPD intervals: HPD interval (0.84): 0.00945106234707..12.35103303631669 HPD interval (0.95): 0.00579246520232..17.46854901300600 HPD interval (0.99): 0.00579246520232..23.41438454974877 Histogram (total 2279) 0.006: 116 ################################## (0.051 / 0.051) 0.856: 205 ############################################################ (0.090 / 0.141) 1.706: 194 ######################################################### (0.085 / 0.226) 2.556: 170 ################################################## (0.075 / 0.301) 3.406: 172 ################################################## (0.075 / 0.376) 4.256: 185 ###################################################### (0.081 / 0.457) 5.106: 144 ########################################## (0.063 / 0.520) 5.956: 128 ##################################### (0.056 / 0.577) 6.806: 98 ############################# (0.043 / 0.620) 7.655: 100 ############################# (0.044 / 0.663) 8.505: 80 ####################### (0.035 / 0.699) 9.355: 87 ######################### (0.038 / 0.737) 10.205: 88 ########################## (0.039 / 0.775) 11.055: 73 ##################### (0.032 / 0.807) 11.905: 72 ##################### (0.032 / 0.839) 12.755: 65 ################### (0.029 / 0.867) 13.605: 39 ########### (0.017 / 0.885) 14.455: 39 ########### (0.017 / 0.902) 15.305: 47 ############## (0.021 / 0.922) 16.155: 32 ######### (0.014 / 0.936) 17.005: 29 ######## (0.013 / 0.949) 17.855: 22 ###### (0.010 / 0.959) 18.705: 23 ####### (0.010 / 0.969) 19.555: 13 #### (0.006 / 0.975) 20.405: 8 ## (0.004 / 0.978) 21.255: 10 ### (0.004 / 0.982) 22.105: 6 ## (0.003 / 0.985) 22.955: 11 ### (0.005 / 0.990) 23.805: 8 ## (0.004 / 0.993) 24.655: 4 # (0.002 / 0.995) 25.505: 3 # (0.001 / 0.996) 26.355: 4 # (0.002 / 0.998) 27.205: 0 (0.000 / 0.998) 28.055: 1 (0.000 / 0.999) 28.905: 1 (0.000 / 0.999) 29.755: 0 (0.000 / 0.999) 30.605: 1 (0.000 / 1.000) 31.454: 0 (0.000 / 1.000) 32.304: 0 (0.000 / 1.000) 33.154: 1 (0.000 / 1.000) var : p Probabilities: true: 1.0000000000000000 mean = [true = 1.0] HPD intervals: show_hpd_intervals: data is not numeric Histogram (total 2279) true : 2279 ############################################################ (1.000 / 1.000) */ go ?=> member(Obs,[false,true]), println(obs=Obs), reset_store, run_model(10_000,$model(Obs),[show_probs_trunc,mean, show_hpd_intervals,hpd_intervals=[0.84,0.95,0.99], show_histogram ]), nl, % show_store_lengths,nl, fail, nl. go => true. model(Obs) => HeightMale = normal_dist(172,sqrt(64)), HeightFemale = normal_dist(164,sqrt(64)), Diff = HeightFemale - HeightMale, P = check(HeightFemale > HeightMale), if Obs == true then observe(HeightFemale > HeightMale) end, if Obs == false ; observed_ok then add("height female",HeightFemale), add("height male",HeightMale), add("diff",Diff), add("p",P), end.