/* Birthday / Coinindcidence probability/quantile in Picat. This is port of R's birthday/coincidence functions (from stats library) """ Probability of coincidences Description: Computes answers to a generalised _birthday paradox_ problem. ‘pbirthday’ computes the probability of a coincidence and ‘qbirthday’ computes the smallest number of observations needed to have at least a specified probability of coincidence. Usage: qbirthday(prob = 0.5, classes = 365, coincident = 2) pbirthday(n, classes = 365, coincident = 2) Arguments: classes: How many distinct categories the people could fall into prob: The desired probability of coincidence n: The number of people coincident: The number of people to fall in the same category Details: The birthday paradox is that a very small number of people, 23, suffices to have a 50-50 chance that two or more of them have the same birthday. This function generalises the calculation to probabilities other than 0.5, numbers of coincident events other than 2, and numbers of classes other than 365. The formula used is approximate for ‘coincident > 2’. The approximation is very good for moderate values of ‘prob’ but less good for very small probabilities. Value: qbirthday: Minimum number of people needed for a probability of at least ‘prob’ that ‘k’ or more of them have the same one out of ‘classes’ equiprobable labels. pbirthday: Probability of the specified coincidence. """ Also, I added rbirthday for generating random variates. Note: The Picat PPL "standard" names are - however - recommended - birthday_dist(Classes,Coincident) - birthday_dist_pdf(Classes,Coincident,X) - birthday_dist_quantile(Classes,Coincident,X) See more in ppl_distributions.pi and ppl_distributions_test.pi Picat> pdf_all($birthday_dist(365,2),0.001,0.8).printf_list 2 0.002739726027397 3 0.008204165884781 4 0.01635591246655 5 0.027135573699793 6 0.040462483649111 7 0.056235703095975 8 0.074335292351669 9 0.094623833889167 10 0.116948177711078 11 0.141141378321733 12 0.167024788838064 13 0.194410275232429 14 0.223102512004973 15 0.252901319763686 16 0.28360400525285 17 0.315007665296561 18 0.34691141787179 19 0.379118526031537 20 0.41143838358058 21 0.443688335165206 22 0.47569530766255 23 0.507297234323986 24 0.538344257914529 25 0.568699703969464 26 0.598240820135939 27 0.626859282263242 28 0.654461472342399 29 0.680968537477777 30 0.706316242719269 31 0.730454633728644 32 0.753347527850321 33 0.774971854175772 34 0.795316864620154 35 0.814383238874715 Picat> quantile_all($birthday_dist(365,2)).printf_list 0.000001 2 0.00001 2 0.001 2 0.01 4 0.025 5 0.05 7 0.25 15 0.5 23 0.75 32 0.84 37 0.975 52 0.99 57 0.999 70 0.99999 89 0.999999 97 Cf my Gamble model gamble_birthday_probability.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. /* var : b Probabilities: 18: 0.0353000000000000 17: 0.0338000000000000 22: 0.0337000000000000 23: 0.0335000000000000 19: 0.0334000000000000 16: 0.0312000000000000 14: 0.0310000000000000 21: 0.0309000000000000 24: 0.0302000000000000 15: 0.0296000000000000 20: 0.0294000000000000 26: 0.0292000000000000 25: 0.0288000000000000 31: 0.0279000000000000 27: 0.0278000000000000 13: 0.0274000000000000 28: 0.0273000000000000 12: 0.0270000000000000 29: 0.0258000000000000 11: 0.0254000000000000 30: 0.0245000000000000 10: 0.0243000000000000 33: 0.0231000000000000 9: 0.0206000000000000 32: 0.0197000000000000 34: 0.0192000000000000 36: 0.0190000000000000 35: 0.0174000000000000 8: 0.0174000000000000 7: 0.0163000000000000 37: 0.0159000000000000 6: 0.0143000000000000 40: 0.0130000000000000 38: 0.0130000000000000 42: 0.0121000000000000 39: 0.0116000000000000 5: 0.0111000000000000 41: 0.0105000000000000 43: 0.0087000000000000 44: 0.0086000000000000 46: 0.0080000000000000 45: 0.0078000000000000 48: 0.0075000000000000 4: 0.0065000000000000 47: 0.0058000000000000 3: 0.0050000000000000 50: 0.0049000000000000 49: 0.0043000000000000 52: 0.0035000000000000 51: 0.0035000000000000 53: 0.0028000000000000 2: 0.0027000000000000 56: 0.0025000000000000 55: 0.0023000000000000 57: 0.0021000000000000 54: 0.0019000000000000 58: 0.0018000000000000 60: 0.0012000000000000 59: 0.0011000000000000 64: 0.0010000000000000 62: 0.0009000000000000 65: 0.0006000000000000 63: 0.0006000000000000 61: 0.0006000000000000 70: 0.0005000000000000 66: 0.0004000000000000 67: 0.0003000000000000 77: 0.0002000000000000 71: 0.0002000000000000 69: 0.0002000000000000 68: 0.0002000000000000 89: 0.0001000000000000 72: 0.0001000000000000 mean = 24.3852 Percentiles: (0.001 2) (0.01 4) (0.025 5) (0.05 7) (0.1 10) (0.25 15) (0.5 23) (0.75 32) (0.84 37) (0.9 41) (0.95 47) (0.975 52) (0.99 57.010000000000218) (0.999 70) (0.9999 77.001199999991513) (0.99999 87.800120000014431) HPD intervals: HPD interval (0.01): 5.00000000000000..5.00000000000000 HPD interval (0.84): 5.00000000000000..37.00000000000000 HPD interval (0.9): 5.00000000000000..42.00000000000000 HPD interval (0.95): 3.00000000000000..47.00000000000000 HPD interval (0.99): 2.00000000000000..57.00000000000000 HPD interval (0.99999): 2.00000000000000..89.00000000000000 Histogram (total 10000) 2.000: 77 ##### (0.008 / 0.008) 4.175: 176 ############ (0.018 / 0.025) 6.350: 306 #################### (0.031 / 0.056) 8.525: 380 ######################### (0.038 / 0.094) 10.700: 497 ################################ (0.050 / 0.144) 12.875: 544 #################################### (0.054 / 0.198) 15.050: 918 ############################################################ (0.092 / 0.290) 17.225: 691 ############################################# (0.069 / 0.359) 19.400: 628 ######################################### (0.063 / 0.422) 21.575: 646 ########################################## (0.065 / 0.486) 23.750: 637 ########################################## (0.064 / 0.550) 25.925: 858 ######################################################## (0.086 / 0.636) 28.100: 531 ################################### (0.053 / 0.689) 30.275: 524 ################################## (0.052 / 0.741) 32.450: 428 ############################ (0.043 / 0.784) 34.625: 366 ######################## (0.037 / 0.821) 36.800: 349 ####################### (0.035 / 0.856) 38.975: 376 ######################### (0.038 / 0.893) 41.150: 226 ############### (0.023 / 0.916) 43.325: 173 ########### (0.017 / 0.933) 45.500: 158 ########## (0.016 / 0.949) 47.675: 133 ######### (0.013 / 0.962) 49.850: 92 ###### (0.009 / 0.971) 52.025: 98 ###### (0.010 / 0.981) 54.200: 42 ### (0.004 / 0.985) 56.375: 46 ### (0.005 / 0.990) 58.550: 29 ## (0.003 / 0.993) 60.725: 18 # (0.002 / 0.995) 62.900: 15 # (0.002 / 0.996) 65.075: 20 # (0.002 / 0.998) 67.250: 5 (0.001 / 0.999) 69.425: 7 (0.001 / 0.999) 71.600: 3 (0.000 / 1.000) 73.775: 0 (0.000 / 1.000) 75.950: 2 (0.000 / 1.000) 78.125: 0 (0.000 / 1.000) 80.300: 0 (0.000 / 1.000) 82.475: 0 (0.000 / 1.000) 84.650: 0 (0.000 / 1.000) 86.825: 1 (0.000 / 1.000) var : p Probabilities: true: 0.5198000000000000 false: 0.4802000000000000 mean = 0.5198 Percentiles: (0.001 false) (0.01 false) (0.025 false) (0.05 false) (0.1 false) (0.25 false) (0.5 true) (0.75 true) (0.84 true) (0.9 true) (0.95 true) (0.975 true) (0.99 true) (0.999 true) (0.9999 true) (0.99999 true) HPD intervals: show_hpd_intervals: data is not numeric Histogram (total 10000) false : 4802 ####################################################### (0.480 / 0.480) true : 5198 ############################################################ (0.520 / 1.000) */ go ?=> reset_store, run_model(10_000,$model,[show_probs,mean, show_percentiles, show_hpd_intervals,hpd_intervals=[0.01,0.84,0.9,0.95,0.99,0.99999], show_histogram % min_accepted_samples=1000,show_accepted_samples=true ]), nl, % show_store_lengths,nl, % fail, nl. go => true. model() => Classes = 365, Coincident = 2, B = birthday_dist(Classes,Coincident), P = check(B <= 23), add("b",B), add("p",P).