/* Geometric cereal box (PPL) in Picat. From Mathematica (GeometricDistribution) """ A cereal box contains one out of a set of n different plastic animals. The animals are equally likely to occur, independently of what animals are in other boxes. Simulate the animal collection process, assuming there are 10 animals for 25 boxes: RandomVariate[DiscreteUniformDistribution[{1, 10}], 25] -> {8,5,4,2,2,9,2,5,2,3,4,5,1,10,1,6,8,8,6,10,8,10,4,7,4} After k unique animals have been collected, the number of boxes needed to find a new unique animal among the remaining n-k follows a geometric distribution with parameter 1-k/n. Find the expected number of boxes needed to get a new unique animal: e = Expectation[x + 1, x -> GeometricDistribution[1 - k/n]] -> -(n/(k - n)) Number of boxes before next unique animal: Block[{n = 5}, DiscretePlot[e, {k, 1, n - 1}, ExtentSize -> 1/2]] Find the expected number of boxes needed to collect 6 unique animals: Block[{n = 5}, Sum[e, {k, 0, n - 1}]] --> 137/12 (11.4167) Table of this Table[Sum[e, {k, 0, n - 1}], {n, 0, 10}] N@% -> {0, 1, 3, 11/2, 25/3, 137/12, 147/10, 363/20, 761/35, 7129/280, 7381/252} {0., 1., 3., 5.5, 8.33333, 11.4167, 14.7, 18.15, 21.7429, 25.4607, 29.2897} The sum: Sum[e, {k, 0, n - 1}] -> -n (-EulerGamma - PolyGamma[0, 1 + n]) """ Also see https://en.wikipedia.org/wiki/Coupon_collector%27s_problem This is a port of my Racket/Gamble model gamble_geometric_cereal_box.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. /* Here are results for n=10, 100, and 1000. n = 10 theoretical = 29.2897 approx = 29.298 mean = 29.1145 Histogram: 10.000: 30 # (0.003 / 0.003) 13.350: 436 ############### (0.044 / 0.047) 16.700: 893 ############################## (0.089 / 0.136) 20.050: 1329 ############################################ (0.133 / 0.269) 23.400: 1792 ############################################################ (0.179 / 0.448) 26.750: 1160 ####################################### (0.116 / 0.564) 30.100: 1020 ################################## (0.102 / 0.666) 33.450: 1093 ##################################### (0.109 / 0.775) 36.800: 555 ################### (0.056 / 0.831) 40.150: 434 ############### (0.043 / 0.874) 43.500: 427 ############## (0.043 / 0.917) 46.850: 216 ####### (0.022 / 0.938) 50.200: 171 ###### (0.017 / 0.956) 53.550: 146 ##### (0.015 / 0.970) 56.900: 81 ### (0.008 / 0.978) 60.250: 52 ## (0.005 / 0.983) 63.600: 62 ## (0.006 / 0.990) 66.950: 26 # (0.003 / 0.992) 70.300: 15 # (0.002 / 0.994) 73.650: 18 # (0.002 / 0.996) 77.000: 12 (0.001 / 0.997) 80.350: 14 (0.001 / 0.998) 83.700: 4 (0.000 / 0.999) 87.050: 4 (0.000 / 0.999) 90.400: 4 (0.000 / 0.999) 93.750: 0 (0.000 / 0.999) 97.100: 0 (0.000 / 0.999) 100.450: 2 (0.000 / 1.000) 103.800: 0 (0.000 / 1.000) 107.150: 1 (0.000 / 1.000) 110.500: 2 (0.000 / 1.000) 113.850: 0 (0.000 / 1.000) 117.200: 0 (0.000 / 1.000) 120.550: 0 (0.000 / 1.000) 123.900: 0 (0.000 / 1.000) 127.250: 0 (0.000 / 1.000) 130.600: 0 (0.000 / 1.000) 133.950: 0 (0.000 / 1.000) 137.300: 0 (0.000 / 1.000) 140.650: 1 (0.000 / 1.000) HPD intervals: HPD interval (0.84): 14.00000000000000..40.00000000000000 HPD interval (0.9): 13.00000000000000..44.00000000000000 HPD interval (0.95): 11.00000000000000..50.00000000000000 HPD interval (0.99): 11.00000000000000..66.00000000000000 HPD interval (0.999): 10.00000000000000..88.00000000000000 HPD interval (0.9999): 10.00000000000000..111.00000000000000 HPD interval (0.99999): 10.00000000000000..144.00000000000000 n = 100 theoretical = 518.738 approx = 518.739 mean = 517.89 Histogram: 265.000: 6 (0.001 / 0.001) 291.300: 44 ### (0.004 / 0.005) 317.600: 128 ######## (0.013 / 0.018) 343.900: 345 #################### (0.035 / 0.052) 370.200: 559 ################################# (0.056 / 0.108) 396.500: 719 ########################################## (0.072 / 0.180) 422.800: 924 ###################################################### (0.092 / 0.272) 449.100: 1020 ############################################################ (0.102 / 0.374) 475.400: 957 ######################################################## (0.096 / 0.470) 501.700: 877 #################################################### (0.088 / 0.558) 528.000: 867 ################################################### (0.087 / 0.645) 554.300: 676 ######################################## (0.068 / 0.712) 580.600: 584 ################################## (0.058 / 0.771) 606.900: 533 ############################### (0.053 / 0.824) 633.200: 390 ####################### (0.039 / 0.863) 659.500: 286 ################# (0.029 / 0.891) 685.800: 234 ############## (0.023 / 0.915) 712.100: 207 ############ (0.021 / 0.936) 738.400: 154 ######### (0.015 / 0.951) 764.700: 120 ####### (0.012 / 0.963) 791.000: 101 ###### (0.010 / 0.973) 817.300: 54 ### (0.005 / 0.978) 843.600: 41 ## (0.004 / 0.983) 869.900: 31 ## (0.003 / 0.986) 896.200: 30 ## (0.003 / 0.989) 922.500: 29 ## (0.003 / 0.992) 948.800: 22 # (0.002 / 0.994) 975.100: 15 # (0.002 / 0.995) 1001.400: 13 # (0.001 / 0.997) 1027.700: 6 (0.001 / 0.997) 1054.000: 4 (0.000 / 0.998) 1080.300: 7 (0.001 / 0.998) 1106.600: 4 (0.000 / 0.999) 1132.900: 6 (0.001 / 0.999) 1159.200: 0 (0.000 / 0.999) 1185.500: 1 (0.000 / 0.999) 1211.800: 0 (0.000 / 0.999) 1238.100: 1 (0.000 / 0.999) 1264.400: 1 (0.000 / 1.000) 1290.700: 4 (0.000 / 1.000) HPD intervals: HPD interval (0.84): 332.00000000000000..642.00000000000000 HPD interval (0.9): 327.00000000000000..698.00000000000000 HPD interval (0.95): 313.00000000000000..765.00000000000000 HPD interval (0.99): 288.00000000000000..926.00000000000000 HPD interval (0.999): 267.00000000000000..1126.00000000000000 HPD interval (0.9999): 265.00000000000000..1294.00000000000000 HPD interval (0.99999): 265.00000000000000..1317.00000000000000 n = 1000 theoretical = 7485.47 approx = 7485.47 mean = 7477.08 Histogram: 4564.000: 1 (0.000 / 0.000) 4868.525: 14 # (0.001 / 0.002) 5173.050: 57 ### (0.006 / 0.007) 5477.575: 156 ######## (0.016 / 0.023) 5782.100: 444 ####################### (0.044 / 0.067) 6086.625: 735 ###################################### (0.073 / 0.141) 6391.150: 975 ################################################### (0.098 / 0.238) 6695.675: 1070 ######################################################## (0.107 / 0.345) 7000.200: 1147 ############################################################ (0.115 / 0.460) 7304.725: 1050 ####################################################### (0.105 / 0.565) 7609.250: 901 ############################################### (0.090 / 0.655) 7913.775: 772 ######################################## (0.077 / 0.732) 8218.300: 640 ################################# (0.064 / 0.796) 8522.825: 482 ######################### (0.048 / 0.844) 8827.350: 402 ##################### (0.040 / 0.885) 9131.875: 314 ################ (0.031 / 0.916) 9436.400: 205 ########### (0.021 / 0.936) 9740.925: 167 ######### (0.017 / 0.953) 10045.450: 107 ###### (0.011 / 0.964) 10349.975: 104 ##### (0.010 / 0.974) 10654.500: 60 ### (0.006 / 0.980) 10959.025: 55 ### (0.005 / 0.986) 11263.550: 35 ## (0.004 / 0.989) 11568.075: 27 # (0.003 / 0.992) 11872.600: 22 # (0.002 / 0.994) 12177.125: 20 # (0.002 / 0.996) 12481.650: 10 # (0.001 / 0.997) 12786.175: 5 (0.001 / 0.998) 13090.700: 12 # (0.001 / 0.999) 13395.225: 4 (0.000 / 0.999) 13699.750: 0 (0.000 / 0.999) 14004.275: 2 (0.000 / 0.999) 14308.800: 3 (0.000 / 1.000) 14613.325: 1 (0.000 / 1.000) 14917.850: 0 (0.000 / 1.000) 15222.375: 0 (0.000 / 1.000) 15526.900: 0 (0.000 / 1.000) 15831.425: 0 (0.000 / 1.000) 16135.950: 0 (0.000 / 1.000) 16440.475: 1 (0.000 / 1.000) HPD intervals: HPD interval (0.84): 5776.00000000000000..8938.00000000000000 HPD interval (0.9): 5603.00000000000000..9332.00000000000000 HPD interval (0.95): 5404.00000000000000..10047.00000000000000 HPD interval (0.99): 5075.00000000000000..11683.00000000000000 HPD interval (0.999): 4855.00000000000000..13423.00000000000000 HPD interval (0.9999): 4564.00000000000000..14559.00000000000000 HPD interval (0.99999): 4564.00000000000000..16745.00000000000000 */ go ?=> Runs = 10000, member(N,[10,100,1000]), nl, println(n=N), println(theoretical=coupon_collectors_problem_theoretical(N)), println(approx=coupon_collectors_problem_approx(N)), reset_store(), run_model(Runs,$model(N),[mean,show_percentiles,show_hpd_intervals,show_histogram]), fail, nl. go => true. model(N) => % Note that we calculate for n-1 Vs = [ 1+geometric_dist(1 - (I / N)) : I in 0..N-1], V = Vs.sum, add("v",V).