The Birthday Bet

Mark Senn

December 13, 2006

In a class of 23 students there is a 50.8% probability that two or more students have the same birthday.  This is a result from the classic statistics "birthday problem" (Dr. Math, 2006; Reese, 2005; Weisstein 2003).  (Year of birth, leap years, and non-uniform distribution of birthday issues are ignored.)

Birthday Combinations

First, for warm-up, how many possible combinations of birthdays are there for people in the class?  For one person there are 365 possibilities  For two people, 365 × 365 combinations.  For ten people

365 ^ 10

combinations.  For n people and d days, there are (d^n) possibilites.  So, for 365 days and 23 people there are

d = 365;
n = 23;
d^n

combinations.

Birthdays on Same Day

How likely is it that two or more people in a group of n people have their birthday on the same month and day?

The Hard Way

One can find the probability of two people, three people, and so on up to all n people having a birthday (Clough, 2006).  This is complicated, see Weisstein (2003)  for details.

The Easy Way

An easier way to solve the problem is to find the probability that everyone in the class has separate birthdays and subtract that from 1.  The first person has d/d ways to have a birthday, the second (d - 1)/ d ways to have a birthday on a different day than the first, the third has (d - 2)/d ways to have a birthday on a different day than the first and second, etc.  So, doing the math manually, the probability of everyone having a different birthday is

((d (d - 1) ⋯ (d - n + 
                    1))/(d^n))

=((d(d - 1) ⋯(
                    d - n + 1) (d - n)!)/((d^n)(
                    d - n)!))

=((d!)/((d^n)(
                  d - n)!)).

Or, Mathematica can be used to do the math (Agrawal, 2006).

Clear[d,n];
expr = Product[d-i, {i,0,n-1}] / d^n

fexp = FunctionExpand[
    expr,
    d > 0  &&  n > 0  &&  Element[{d,n}, Integers]
]

Mathematica apparently doesn't know that Gamma[x] can be simplified to (x-1)! if x is an integer so we define the rule for Mathematica to use.

Simplify[
    fexp,
    d > 0  &&  n > 0  &&  Element[{d,n}, Integers]
]  /.  {Gamma[x_] -> (x-1)!}

Does anyone know of a better way to do the above?

Define the Probability Function

Set d and define the probability function.

d = 365;
f[d_, n_] := 1 - d! / (d^n (d-n)!);

Does the function work as expected?

p2 = f[d,2]; p3 = f[d,3]; p365 = f[d,365]; p366 = f[d,366];

p2

Looks good for two people.  The probability of the second person having the same birthday as first should be 1/365.

p3

N[p3]

Looks reasonable for three people.

p365

Show numerical value to 100 digit precision.

N[p365,100]

The probability of two people having the same birthday in a group of 365 people being one doesn't make any sense.  Everyone could have a different birthday and it would require 366 people to have two people with the same birthday.  Maybe printing with 500 digit precision will give a clue what's going on.

N[p365,500]

Aha!  It was almost one, but not quite.

p366

This makes sense, with 366 people at least two people are guaranteed to have the same birthday.

What is the probability that exactly 366 people will be required?

p366365 = f[d,366] - f[d,365];
N[p366365]

So, 366 people are needed

N[1/p366365]

about 7 out of every (10^156) times.  Aside: there are around (10^80) atoms in the universe [Champion, 1998; Wikipedia, 2006] so the probability of choosing a specific atom from the universe at random (talk about a needle in a haystack!) is much better than choosing a group of 365 random people (remember this is a mathematical model and there is an inexhaustible supply of people) having different birthdays.

Make Some Plots

Plot[f[d,n], {n, 2, 366}]

The above plot wasn't very helpful.  Make a better plot using some options.

tplot = Plot[
    f[d,n],
    {n,2,364},
    AxesLabel    -> {"n","p"},
    GridLines    -> {{366}, {1}},
    PlotDivision -> 100,
    PlotPoints   -> 1000,
    PlotRange    -> {{0,366}, {0,1}},
    Ticks        -> {Range[0,366,50], Range[0,1,0.2]}
]

Let's get a better look at an area of interest in the graph.
Note that a new Plot[] command isn't needed---can just Show[] the old plot with some other options.

plot1 = Show[
        tplot,
        {
            GridLines -> {{80}, {1}},
            PlotRange -> {{0,80}, {0,1}},
            Ticks     -> {Range[0,80,10], Range[0,1,0.1]}
        }
    ]

For some specific probabilites, solve for n.  Normally "n -> ..." rules are produced.

points = Table[
    {v, FindRoot[f[365,n] == v, {n,5,50}]},
    {v, 0.1, 0.9, 0.1}
]

Show how to replace the "n -> ..." rules with the values for n.

points = Table[
    {v, n /. FindRoot[f[365,n] == v, {n,20,50}]},
    {v, 0.1, 0.9, 0.1}
]

Indicate the values of n for the specific probabilities just computed on the graph.

pp = Table[Range[9]];

For[i=1, i<=Length[points], i++,
    pp[[i]] = Graphics[
                 Line[{
                     {0,             points[[i,1]]},
                     {points[[i,2]], points[[i,1]]},
                     {points[[i,2]], 0}
                 }]
             ]
]

Show[plot1, Table[pp[[i]], {i,1,9}], GridLines->{{80},{1}}]

Given Probabilities, How Many People are Needed?

Find the number of people required for some specific probabilities.

n = .;  (* n is currently 23, make it an unknown *)
FindRoot[f[365,n] == 0.5,      {n,10, 50}]
FindRoot[f[365,n] == 0.99,     {n,40,364}]
FindRoot[f[365,n] == 0.999,    {n,40,364}]
FindRoot[f[365,n] == 0.9999,   {n,40,364}]
FindRoot[f[365,n] == 0.99999,  {n,40,364}]
FindRoot[f[365,n] == 0.999999, {n,40,364}]

Given Number of People, What are Probabilities?

f[d,23] //N

(* check around 365/2 = 182.5 people *)
f[d,183] //N

N[f[d,183], 100]

Run a Simulation

Let's run a computer simulation and see how it compares with the theoretical results.

(* Make the number of days in a year a variable. *)
(* Setting d to 7 will give probabilites for     *)
(* birthdays on same day of week.                *)
d = 365;

(* total birthday count for all runs *)
totals = Table[0, {d+1}];

(* find this number of matching birthdays *)
runs = 100;

For[trial = 1,  trial <= runs,  trial++,

    (* all birthdays have occurred zero times during this trial *)
    count = Table[0, {d}];
    
    (* we haven't found a birthday match yet *)
    match = 0;
    
    (* might have to generate up to d + 1 birthdays for match *)
    For[i = 1,  i <= d + 1,  i++,
    
        birthday = Random[Integer, {1,d}];
        
        (* if already have a birthday on that day we're done *)
        If [count[[birthday]] == 1,
            match = i;
            Break[];
        ];
        
        (* add one to the number of birthdays on this day *)
        count[[birthday]] = count[[birthday]] + 1;
        (* or we could have used the shorter equivalent count[[birthay]]++ *)
        
    ];
    
    (* something's wrong if no matching birthday here *)
    If [match == 0,
        Print["Couldn't find a birthday match."];
        Print["Something is wrong with program."];
        Abort[];
    ];
        
    (* Increment number of days needed to find a matching birthday. *)
    totals[[match]]++;
                
];

totals

ListPlot[totals];

Convert the distribution to a cumulative distribution.

For [i = 2,  i <= d + 1,  i++,
    totals[[i]] = totals[[i-1]] + totals[[i]];
];

totals = totals / runs;

plot2 = ListPlot[totals];

Compare theory with simulation

p = Show[
        {plot1, plot2},
        {
            GridLines -> {{80}, {1}},
            PlotRange -> {{0,80}, {0,1}},
            Ticks     -> {
                Range[0,80,10],
                Range[0,1,0.1]
            }
        }
    ]

The theory and simulation results appear to match reasonably closely.  Try changing the number of runs from 100 to 1000 and do Kernel => Evaluation => Evaluate Notebook to see if the results match better.  Do the same thing for 10 000 runs (this may take a little while to run).

References

Agrawal, Jyota.  Wolfram Research Technical Support tracking number 33367.  2006-12-13.

Champion, Matthew.  "How many atoms make up the universe?".  Last modified 1998-09-11.  Retrieved 2006-11-21.

Clough, Tim.  Personal communication.  2006-11-12.

Dr. Math.  "The Birthday Problem".  Retrieved 2006-11-21.

Reese, George.  "The Birthday Problem :  
              A short lesson in probability."  Last modified 2005-08-24.  Retrieved 2006-11-20.

Weisstein, Eric W.  "Birthday Problem".  From MathWorld--A Wolfram Web Resource.  Last modified: 2003-08-17.  Retrieved 2006-11-20.

Wikipedia, the free encyclopedia.  "Observable universe".  Last Modified: 2006-11-19.  Retrieved 2006-11-21.