/* Banach's match box problem in Picat. Banach's match box problem in Racket Gamble. From Gunnar Blom, Lars Holst, Dennis Sandell: "Problems and Snapshots from the World of Probability" Page 9f, Problem 1.7 Number of walks until no shoes """ A person has, in each of his two pockets, a box with n matches. Now and then he talkes a match from a randomly chosen box until he finds the selected box empty. Find the expectation of the number, R, of remaining matches in the other box. """ Cf my Gamble model gamble_banachs_match_box_problem.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. /* n = 10 theoretical = 2.7020575410170005 var : len Probabilities (truncated): 20: 0.1850000000000000 19: 0.1791000000000000 18: 0.1643000000000000 17: 0.1403000000000000 ......... 13: 0.0355000000000000 12: 0.0149000000000000 11: 0.0058000000000000 10: 0.0011000000000000 mean = 17.3391 HPD intervals: HPD interval (0.84): 15.00000000000000..20.00000000000000 Histogram: 10: 11 (0.001 / 0.001) 11: 58 ## (0.006 / 0.007) 12: 149 ##### (0.015 / 0.022) 13: 355 ############ (0.035 / 0.057) 14: 558 ################## (0.056 / 0.113) 15: 938 ############################## (0.094 / 0.207) 16: 1244 ######################################## (0.124 / 0.331) 17: 1403 ############################################## (0.140 / 0.472) 18: 1643 ##################################################### (0.164 / 0.636) 19: 1791 ########################################################## (0.179 / 0.815) 20: 1850 ############################################################ (0.185 / 1.000) var : num left in other box Probabilities (truncated): 0: 0.1850000000000000 1: 0.1791000000000000 2: 0.1643000000000000 3: 0.1403000000000000 ......... 7: 0.0355000000000000 8: 0.0149000000000000 9: 0.0058000000000000 10: 0.0011000000000000 mean = 2.6609 HPD intervals: HPD interval (0.84): 0.00000000000000..5.00000000000000 Histogram: 0: 1850 ############################################################ (0.185 / 0.185) 1: 1791 ########################################################## (0.179 / 0.364) 2: 1643 ##################################################### (0.164 / 0.528) 3: 1403 ############################################## (0.140 / 0.669) 4: 1244 ######################################## (0.124 / 0.793) 5: 938 ############################## (0.094 / 0.887) 6: 558 ################## (0.056 / 0.943) 7: 355 ############ (0.035 / 0.978) 8: 149 ##### (0.015 / 0.993) 9: 58 ## (0.006 / 0.999) 10: 11 (0.001 / 1.000) n = 50 theoretical = 7.0386869500888691 var : len Probabilities (truncated): 100: 0.0796000000000000 97: 0.0776000000000000 98: 0.0766000000000000 96: 0.0749000000000000 ......... 71: 0.0003000000000000 72: 0.0001000000000000 70: 0.0001000000000000 69: 0.0001000000000000 mean = 92.9519 HPD intervals: HPD interval (0.84): 87.00000000000000..100.00000000000000 Histogram: 69: 1 (0.000 / 0.000) 70: 1 (0.000 / 0.000) 71: 3 (0.000 / 0.001) 72: 1 (0.000 / 0.001) 73: 6 (0.001 / 0.001) 74: 5 (0.001 / 0.002) 75: 12 # (0.001 / 0.003) 76: 20 ## (0.002 / 0.005) 77: 33 ## (0.003 / 0.008) 78: 54 #### (0.005 / 0.014) 79: 69 ##### (0.007 / 0.020) 80: 65 ##### (0.006 / 0.027) 81: 93 ####### (0.009 / 0.036) 82: 102 ######## (0.010 / 0.046) 83: 141 ########### (0.014 / 0.061) 84: 177 ############# (0.018 / 0.078) 85: 228 ################# (0.023 / 0.101) 86: 282 ##################### (0.028 / 0.129) 87: 335 ######################### (0.034 / 0.163) 88: 364 ########################### (0.036 / 0.199) 89: 429 ################################ (0.043 / 0.242) 90: 466 ################################### (0.047 / 0.289) 91: 575 ########################################### (0.058 / 0.346) 92: 625 ############################################### (0.062 / 0.409) 93: 647 ################################################# (0.065 / 0.473) 94: 728 ####################################################### (0.073 / 0.546) 95: 714 ###################################################### (0.071 / 0.618) 96: 749 ######################################################## (0.075 / 0.692) 97: 776 ########################################################## (0.078 / 0.770) 98: 766 ########################################################## (0.077 / 0.847) 99: 737 ######################################################## (0.074 / 0.920) 100: 796 ############################################################ (0.080 / 1.000) var : num left in other box Probabilities (truncated): 0: 0.0796000000000000 3: 0.0776000000000000 2: 0.0766000000000000 4: 0.0749000000000000 ......... 29: 0.0003000000000000 31: 0.0001000000000000 30: 0.0001000000000000 28: 0.0001000000000000 mean = 7.0481 HPD intervals: HPD interval (0.84): 0.00000000000000..13.00000000000000 Histogram: 0: 796 ############################################################ (0.080 / 0.080) 1: 737 ######################################################## (0.074 / 0.153) 2: 766 ########################################################## (0.077 / 0.230) 3: 776 ########################################################## (0.078 / 0.307) 4: 749 ######################################################## (0.075 / 0.382) 5: 714 ###################################################### (0.071 / 0.454) 6: 728 ####################################################### (0.073 / 0.527) 7: 647 ################################################# (0.065 / 0.591) 8: 625 ############################################### (0.062 / 0.654) 9: 575 ########################################### (0.058 / 0.711) 10: 466 ################################### (0.047 / 0.758) 11: 429 ################################ (0.043 / 0.801) 12: 364 ########################### (0.036 / 0.837) 13: 335 ######################### (0.034 / 0.871) 14: 282 ##################### (0.028 / 0.899) 15: 228 ################# (0.023 / 0.922) 16: 177 ############# (0.018 / 0.939) 17: 141 ########### (0.014 / 0.954) 18: 102 ######## (0.010 / 0.964) 19: 93 ####### (0.009 / 0.973) 20: 65 ##### (0.006 / 0.979) 21: 69 ##### (0.007 / 0.986) 22: 54 #### (0.005 / 0.992) 23: 33 ## (0.003 / 0.995) 24: 20 ## (0.002 / 0.997) 25: 12 # (0.001 / 0.998) 26: 5 (0.001 / 0.999) 27: 6 (0.001 / 0.999) 28: 1 (0.000 / 0.999) 29: 3 (0.000 / 1.000) 30: 1 (0.000 / 1.000) 31: 1 (0.000 / 1.000) n = 100 theoretical = 10.3261058897212070 var : len Probabilities (truncated): 200: 0.0603000000000000 199: 0.0578000000000000 198: 0.0575000000000000 196: 0.0575000000000000 ......... 157: 0.0004000000000000 153: 0.0001000000000000 152: 0.0001000000000000 151: 0.0001000000000000 mean = 189.739 HPD intervals: HPD interval (0.84): 181.00000000000000..200.00000000000000 Histogram: 151.000: 1 (0.000 / 0.000) 152.225: 1 (0.000 / 0.000) 153.450: 1 (0.000 / 0.000) 154.675: 0 (0.000 / 0.000) 155.900: 0 (0.000 / 0.000) 157.125: 4 (0.000 / 0.001) 158.350: 0 (0.000 / 0.001) 159.575: 10 # (0.001 / 0.002) 160.800: 6 (0.001 / 0.002) 162.025: 7 (0.001 / 0.003) 163.250: 6 (0.001 / 0.004) 164.475: 28 # (0.003 / 0.006) 165.700: 15 # (0.002 / 0.008) 166.925: 19 # (0.002 / 0.010) 168.150: 30 ## (0.003 / 0.013) 169.375: 36 ## (0.004 / 0.016) 170.600: 110 ###### (0.011 / 0.027) 171.825: 61 ### (0.006 / 0.033) 173.050: 85 #### (0.009 / 0.042) 174.275: 89 ##### (0.009 / 0.051) 175.500: 204 ########## (0.020 / 0.071) 176.725: 145 ####### (0.015 / 0.086) 177.950: 154 ######## (0.015 / 0.101) 179.175: 178 ######### (0.018 / 0.119) 180.400: 423 ##################### (0.042 / 0.161) 181.625: 237 ############ (0.024 / 0.185) 182.850: 263 ############# (0.026 / 0.211) 184.075: 292 ############### (0.029 / 0.240) 185.300: 313 ################ (0.031 / 0.272) 186.525: 717 #################################### (0.072 / 0.343) 187.750: 342 ################# (0.034 / 0.378) 188.975: 394 #################### (0.039 / 0.417) 190.200: 454 ####################### (0.045 / 0.462) 191.425: 950 ################################################ (0.095 / 0.557) 192.650: 505 ########################## (0.051 / 0.608) 193.875: 524 ########################### (0.052 / 0.660) 195.100: 500 ######################### (0.050 / 0.710) 196.325: 575 ############################# (0.058 / 0.768) 197.550: 1140 ########################################################## (0.114 / 0.882) 198.775: 1181 ############################################################ (0.118 / 1.000) var : num left in other box Probabilities (truncated): 0: 0.0603000000000000 1: 0.0578000000000000 4: 0.0575000000000000 2: 0.0575000000000000 ......... 43: 0.0004000000000000 49: 0.0001000000000000 48: 0.0001000000000000 47: 0.0001000000000000 mean = 10.2614 HPD intervals: HPD interval (0.84): 0.00000000000000..19.00000000000000 Histogram: 0.000: 603 ################################ (0.060 / 0.060) 1.225: 578 ############################## (0.058 / 0.118) 2.450: 1140 ############################################################ (0.114 / 0.232) 3.675: 575 ############################## (0.058 / 0.290) 4.900: 500 ########################## (0.050 / 0.340) 6.125: 524 ############################ (0.052 / 0.392) 7.350: 505 ########################### (0.051 / 0.443) 8.575: 950 ################################################## (0.095 / 0.537) 9.800: 454 ######################## (0.045 / 0.583) 11.025: 394 ##################### (0.039 / 0.622) 12.250: 342 ################## (0.034 / 0.656) 13.475: 717 ###################################### (0.072 / 0.728) 14.700: 313 ################ (0.031 / 0.759) 15.925: 292 ############### (0.029 / 0.789) 17.150: 263 ############## (0.026 / 0.815) 18.375: 237 ############ (0.024 / 0.839) 19.600: 423 ###################### (0.042 / 0.881) 20.825: 178 ######### (0.018 / 0.899) 22.050: 154 ######## (0.015 / 0.914) 23.275: 145 ######## (0.015 / 0.929) 24.500: 204 ########### (0.020 / 0.949) 25.725: 89 ##### (0.009 / 0.958) 26.950: 85 #### (0.009 / 0.966) 28.175: 61 ### (0.006 / 0.973) 29.400: 110 ###### (0.011 / 0.984) 30.625: 36 ## (0.004 / 0.987) 31.850: 30 ## (0.003 / 0.990) 33.075: 19 # (0.002 / 0.992) 34.300: 15 # (0.002 / 0.994) 35.525: 28 # (0.003 / 0.996) 36.750: 6 (0.001 / 0.997) 37.975: 7 (0.001 / 0.998) 39.200: 6 (0.001 / 0.998) 40.425: 10 # (0.001 / 0.999) 41.650: 0 (0.000 / 0.999) 42.875: 4 (0.000 / 1.000) 44.100: 0 (0.000 / 1.000) 45.325: 0 (0.000 / 1.000) 46.550: 1 (0.000 / 1.000) 47.775: 2 (0.000 / 1.000) */ % The Stirling approach (from the book page 11): theoretical(N) = to_fstring("%.16f",2* sqrt(N/pi) - 1 + 3/(4*sqrt(N*pi))). go ?=> Ns = [10,50,100,1000], foreach(NN in Ns) println(NN=theoretical(NN)), end, nl, member(N,Ns.delete(1000)), println(n=N), println(theoretical=theoretical(N)), reset_store, run_model(10_000,$model(N),[show_probs_trunc,mean,show_hpd_intervals,show_histogram, hpd_intervals=[0.84] %,random_seed=0 ]), nl, fail, nl. go => true. select_match_box(A,Left,Right) = Res => Pick = condt(flip(),l,r), if (Pick == l,Left == 0) ; (Pick == r, Right == 0) then Res = [A,cond(Pick == l,Right,Left)] else if Pick == l then Res = select_match_box(A ++ [Pick],Left-1,Right) else Res = select_match_box(A ++ [Pick],Left,Right-1) end end. model(N) => Res = select_match_box([],N,N), [A,NumLeftInOtherBox] = Res, % Freq=A.get_freq, % add("a",A), add("len",A.len), add("num left in other box",NumLeftInOtherBox). % add("freq",Freq).