Simulating the Birthday Problem with data derived probabilities

You've probably heard of the Birthday Paradox: it only takes a small gathering of people before it's quite likely that two of them share the same birthday. You can solve the problem analytically or with simulation, but usually in either case simplifying assumptions are made (no-one born on February 29, for example). Joe Rickert uses Revolution R Enterprise 6 to investigate whether these assumptions make a difference, and finds evidence for a "Summer Baby Boom" along the way. - Ed.

William Feller’s “An Introduction to Probability Theory and Its Applications” (Wiley 1950) is a beautiful book. Written with remarkable clarity, it is still very much worth reading. On page 33 (3^rd ed. 1968) Feller gives a formula for calculating the probability that at least two people in a random sample of size r share a common birthday. Assuming 365 days in a year and a uniform probability of being born on any given day, Feller uses a simple combinatorial argument to derive his “first approximation” solution and notes that when r = 23 there is a better than 50% chance of a match.  He also provides a formula for a numerical approximation suggesting that for “large” r the calculation should be done with logarithms. Apparently, in Feller’s time calculating the probability of at least one matching birthday for a given sample size was a bit of a pain.

Nevertheless, it is astounding that of all the problems and examples in Feller’s two-volume work the “birthday problem”, as we now call it, should have reached urban legend fame as an example of how non-intuitive probability can be. It’s on YouTube, NPR’s “Math Guy” has discussed it, and Trumbo et al. published a very elegant R based presentation in the 2004 Proceedings of the American Statistical Association, Session 345. In a recent blog post, Laura McLay poses the question: Have you simulated the Birthday Problem with and unequal birthday distribution?”. And so it is, having been given four years (1985 to 1988) of day of the week birth data (which my colleague, Sue Ranney, imported from the CDC for a talk she gave at Use-R last year) I find myself compelled to do the simulation. The small portion of the CDC file that I have looks like this:

Number of observations: 15252733
Number of variables: 5
Number of blocks: 32
Variable information:

Var 1: DOB_DAY, Day of Month
       Type: integer, Low/High: (1, 31)
Var 2: DOB_MM, Birth Month
       Type: integer, Low/High: (1, 12)
Var 3: DOB_YY, Birth Year
       Type: integer, Low/High: (1985, 1988)
Var 4: WEEKDAY, Day of Week Child Born
       7 factor levels: Sunday Monday Tuesday Wednesday Thursday Friday Saturday
Var 5: SEX, Sex of Infant
       2 factor levels: Male Female

15M records makes this an ideal job for Revolution’s RevoScaleR package. (The code for all of the data preparation and probability calculations are in the script Heat Map.R.) The first step was to transform this information into convenient form for tabulating births by day of year and estimation probabilities. The data transformations including the calculation to assign births to a 366 day year were accomplished using the transforms option of the rxDataStep function. The logic computing day of year (DOY) is

DOY = DateVar - startOfYear + 1,
DOY = ifelse(DOB_YY != 1988 & DOB_MM > 2, DOY + 1, DOY)

where the variables are as given above and startOfYear contains the date of the first of the year for each observation

Tabulating number of births by day of year was done with the function rxCube that produces a cross tabulation for factor variables in long form. On my 4-core Dell laptop with 8GB of RAM it only takes about .3 seconds to tabulate 15 million records. Dividing by the total number of births gives an estimate of the probability of birth for each day. The mean probability (shown below) is very close to 1/366 as would be expected:

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
0.0006424 0.0026360 0.0027260 0.0027320 0.0028410 0.0031970

Figure 1 plots a regression line on the data that shows that the probabilities trend slightly upward as the year progresses. The heat map in Figure 2 shows the same features of the NPR heat map mentioned in Laura’s post: low probabilities in the beginning of January, high probabilities at the end of December and a booming baby season in the summer. Figure 3, a plot of the differences between male and female birth probabilities shows some interesting structure. The differences are small numbers, but it does appear as if the year starts out favoring girls, gives boys the edge in the summer and then goes girls again.

Figure 1

Figure 2

Figure 3

A brute force simulation of the birthday problem now becomes a simple exercise of using R’s sample function to simulate from the estimated probability distribution. The script Birthday Problem contains the code to do this using the RevoScaleR function rxExec to run the simulation in parallel.

rxExec(pbirthday2, rxElemArg(myArgs), taskChunkSize=20)

(bpirthday2 is a function that performs the simulation and rxElemArg passes the required arguments to pbirthday2.) rxExec allows you to run arbitrary R functions in a distributed fashion using all of available computational resources. In my case, the following line of code established a local compute context instructing rxExec to distribute the computations across the 4 available cores of my laptop.

rxSetComputeContext(RxLocalParallel())

For someone who happened to have a cluster available, running the same code on the cluster would only require changing the compute context. The following code snippets respectively show how to set the compute context to be a Microsoft HPC cluster or a Linux cluster.

Code snippet for Microsoft HPC cluster

myCluster <- RxHpcServer(headNode="cluster-head2",
                revoPath="C:\\Revolution\\ . . .",
                shareDir="\\AllShare\\myName",
                dataPath="C:\\data")
rxSetComputeContext(myCluster)

Code snippet for Linux cluster

myCluster <- RxLsfCluster(shareDir="/usr/share/lsf/shares/myName",
                dataPath="/usr/local/data")
rxSetComputeContext(myCluster)

Figure 4 shows the cumulative distribution of the probability of a match built by simulating sample sizes of 2 to 100 people. Note that the color coding indicates in which node (core) that part of the calculation would be done.

Figure 4

Using 10,000 random draws for each point it took about 30 seconds to run this simulation on my laptop. The probability of at least one match for a sample size of 23 people comes out to be .5087, not far from the value: 0. 507297234 that Trumbo got using R and Feller’s approximate analytic solution. This is hardly the last word on the Birthday Problem. (It would be nice to look at the number of births as a time series over more years of data.) Nevertheless, Feller’s first approximation solution is still looking pretty good.

Comments

[Carl Witthoft]

Is the Feb 29 valley adjusted by 4x to account for the rate of occurrence of Feb29s?
Meanwhile, I do find it amusing that the birthrate is slightly lower (heatmap) on the 4th of July and Christmas and New Year's day. I bet that first dip applies only in the USA :-) .
And, finally, a pedantic note: given that the calendar year is cyclic and the end of the year is selected arbitrarily, it's a bad idea to attempt a linear fit, especially after you've overlaid several years' data. What does the fit look like if you plot "true time," i.e. from day 1 of first year to day 365 of the last year?

[Eric]

Link to punkrockor blog post is broken.

[Rick]

Where is the data? The CDC files linked only include Birth Month and Year, and not the day.

[Chris Mulligan]

The CDC data includes birth day for 1969 to 1988.

I actually did a very similar simulation a few days ago using that full data and found that they were nearly identical. Slightly more likely in reality than the simulation, but only 0.14% more likely at n=23 (and n=23 is still the minimum group size necessary for >= 50%.

Post: http://chmullig.com/2012/06/births-by-day-of-year/

[David Smith]

@Chris -- thanks! I really liked your time series chart, which I reproduced in this post.

[Derek O'Connor]

Does Revolution R have a Birthday Problem?

Here is the Matlab answer:

http://www.mathworks.co.uk/matlabcentral/answers/879-does-matlab-have-a-birthday-problem

http://www.scribd.com/doc/48243415/O-Connor-The-Birthday-Paradox-and-Random-Number-Generation