/* Birthday "paradox", two stage model in Picat. https://en.wikipedia.org/wiki/Birthday_problem """ In probability theory, the birthday problem asks for the probability that, in a set of n randomly chosen people, at least two will share a birthday. The birthday paradox refers to the counterintuitive fact that only 23 people are needed for that probability to exceed 50%. """ The "two stage" is the we use an "ounter" sampling function (birthday/2) to generate samples which is then used in the main model (model/1). 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. /* var : p Probabilities: 0.5: 0.3823529411764706 0.51: 0.3137254901960784 0.49: 0.3039215686274510 mean = 0.500098 Percentiles: (0.001 0.48999999999999999) (0.01 0.48999999999999999) (0.025 0.48999999999999999) (0.05 0.48999999999999999) (0.1 0.48999999999999999) (0.25 0.48999999999999999) (0.5 0.5) (0.75 0.51000000000000001) (0.84 0.51000000000000001) (0.9 0.51000000000000001) (0.95 0.51000000000000001) (0.975 0.51000000000000001) (0.99 0.51000000000000001) (0.999 0.51000000000000001) (0.9999 0.51000000000000001) (0.99999 0.51000000000000001) HPD intervals: HPD interval (0.5): 0.49000000000000..0.50000000000000 HPD interval (0.84): 0.49000000000000..0.51000000000000 HPD interval (0.9): 0.49000000000000..0.51000000000000 HPD interval (0.95): 0.49000000000000..0.51000000000000 HPD interval (0.99): 0.49000000000000..0.51000000000000 HPD interval (0.999): 0.49000000000000..0.51000000000000 HPD interval (0.9999): 0.49000000000000..0.51000000000000 HPD interval (0.99999): 0.49000000000000..0.51000000000000 Histogram (total 102) 0.490: 31 ################################################ (0.304 / 0.304) 0.490: 0 (0.000 / 0.304) 0.491: 0 (0.000 / 0.304) 0.491: 0 (0.000 / 0.304) 0.492: 0 (0.000 / 0.304) 0.492: 0 (0.000 / 0.304) 0.493: 0 (0.000 / 0.304) 0.493: 0 (0.000 / 0.304) 0.494: 0 (0.000 / 0.304) 0.494: 0 (0.000 / 0.304) 0.495: 0 (0.000 / 0.304) 0.495: 0 (0.000 / 0.304) 0.496: 0 (0.000 / 0.304) 0.496: 0 (0.000 / 0.304) 0.497: 0 (0.000 / 0.304) 0.497: 0 (0.000 / 0.304) 0.498: 0 (0.000 / 0.304) 0.498: 0 (0.000 / 0.304) 0.499: 0 (0.000 / 0.304) 0.499: 0 (0.000 / 0.304) 0.500: 39 ############################################################ (0.382 / 0.686) 0.500: 0 (0.000 / 0.686) 0.501: 0 (0.000 / 0.686) 0.502: 0 (0.000 / 0.686) 0.502: 0 (0.000 / 0.686) 0.502: 0 (0.000 / 0.686) 0.503: 0 (0.000 / 0.686) 0.504: 0 (0.000 / 0.686) 0.504: 0 (0.000 / 0.686) 0.504: 0 (0.000 / 0.686) 0.505: 0 (0.000 / 0.686) 0.506: 0 (0.000 / 0.686) 0.506: 0 (0.000 / 0.686) 0.506: 0 (0.000 / 0.686) 0.507: 0 (0.000 / 0.686) 0.508: 0 (0.000 / 0.686) 0.508: 0 (0.000 / 0.686) 0.508: 0 (0.000 / 0.686) 0.509: 0 (0.000 / 0.686) 0.510: 32 ################################################# (0.314 / 1.000) var : num people Probabilities: 24: 0.2352941176470588 23: 0.2254901960784314 22: 0.2254901960784314 21: 0.1176470588235294 25: 0.0980392156862745 20: 0.0686274509803922 26: 0.0196078431372549 28: 0.0098039215686275 mean = 22.8725 Percentiles: (0.001 20) (0.01 20) (0.025 20) (0.05 20) (0.1 21) (0.25 22) (0.5 23) (0.75 24) (0.84 24) (0.9 25) (0.95 25) (0.975 25.474999999999994) (0.99 26) (0.999 27.798000000000002) (0.9999 27.979800000000012) (0.99999 27.997980000000013) HPD intervals: HPD interval (0.5): 21.00000000000000..23.00000000000000 HPD interval (0.84): 20.00000000000000..24.00000000000000 HPD interval (0.9): 21.00000000000000..25.00000000000000 HPD interval (0.95): 20.00000000000000..25.00000000000000 HPD interval (0.99): 20.00000000000000..26.00000000000000 HPD interval (0.999): 20.00000000000000..28.00000000000000 HPD interval (0.9999): 20.00000000000000..28.00000000000000 HPD interval (0.99999): 20.00000000000000..28.00000000000000 Histogram (total 102) 20: 7 ################## (0.069 / 0.069) 21: 12 ############################## (0.118 / 0.186) 22: 23 ########################################################## (0.225 / 0.412) 23: 23 ########################################################## (0.225 / 0.637) 24: 24 ############################################################ (0.235 / 0.873) 25: 10 ######################### (0.098 / 0.971) 26: 2 ##### (0.020 / 0.990) 28: 1 ### (0.010 / 1.000) */ go ?=> MaxNumPeople = 100, % Maximum number of people in the room NumGen = 30, % Number of iterations per call in birthdays/2 println([max_num_people=MaxNumPeople,num_gen=NumGen]), reset_store, run_model(10_000,$model(MaxNumPeople,NumGen),[show_probs_trunc,mean,show_percentiles,show_hpd_intervals,show_histogram]), nl, show_store_lengths,nl, fail, nl. go => true. % Generate N birthdays and return the proportion (probability) that % there are at least one common birthday. birthdays(N,NumGen) = Res => Ps = [], foreach(_ in 1..NumGen) % Generate N birthdays Bs = random_integer1_n(365,N), % Check if there is at least common birthday Ps := Ps ++ [cond(N-Bs.remove_dups.len > 0,1,0)] end, Res = Ps.sum / NumGen. model(MaxNumPeople,NumGen) => NumPeople = random_integer(MaxNumPeople), P = birthdays(NumPeople,NumGen), % Ensure that there are (about) 50% change of at least one common birthday. observe( abs(P-0.5)<=0.02), if observed_ok then add("p",P), add("num people",NumPeople), end.