Three-birthday paradox

[bonanova]

Most puzzle addicts know by now that in a room
containing just 23 random people there is a 50%
chance that at least two of them share a birthday.
But how many would there have to be to have an
even chance that at least three share a birthday?

[Logophobic]

Based on simulated trials, the answer is...

Spoiler

87

I recorded 100 runs of 100,000 with this number. The minimum was 50,046. The maximum was 50,316. The average was 50,177.64. Final result: 50.18%

Then I considered February 29, which was conveniently omitted in those trials. I repeated the process, allowing for leap day in its proper ratio. The minimum was 49,891. The maximum was 50,068. The average was 49,999.04. Final result: 50.00%

Very interesting. Given this number of people, the probability that at least three share a birthday might be ever-so-slightly less than 50%.

[bonanova]

What I got, too. I think doing the actual calculation is very much harder than for two people.

[plasmid]

This is a way to calculate probabilities exactly (if you're willing to ignore leap years), although it uses an iterative approach rather than providing a closed form solution.

Spoiler

Define S (singleton) as the number of birthdays that are held by at least one person in the room, D (duo) as the number of birthdays that are held by at least two people in the room, and T (triplet) as the number of birthdays that are held by at least three people in the room. Note that if a birthday is held by two people, that birthday will be counted in both the S and D terms. A room with N people and no more than three sharing a birthday can have its state described as a vector (S, D, T) where S + D + T = N.

If you start with an empty room and add one person, you go from state (0, 0, 0) to state (1, 0, 0) with 100% probability. If you have one person in the room already and you add a second person, you go from state (1, 0, 0) to state (2, 0, 0) with probability 364/365, or you go from state (1, 0, 0) to state (1, 1, 0) with probability 1/365. In general, if you add one person to any state (S, D, 0) then you have a D/365 chance of forming a triplet, an (S-D)/365 chance of adding a new duo from a singleton, and a (365-S)/365 chance of adding a new singleton.

Using this approach, you can make a spreadsheet where the rows represent S and columns represent D. In Excel or Open Office, use column A to label the rows (S) from 0 upward and row 1 to label the columns (D) from 0 upward. In cell B2 put 1 to represent 100% chance of having state (0, 0, 0) if no one is present. In each cell in column B below B2, put the formula representing getting an new birthday added from the state immediately above, which would be (365-D-S)/365 times the probability that you were in the preceding cell to begin with, so in spreadsheet format would be cell B3 with B2*(365-A2)/365 (remember that in column B the value of D is zero so it doesn't need to be included in the formula) and copy/paste that into all of the lower cells in column B. Since we know the answer to the famous two birthdays question is 23, the state (23, 0, 0) ought to have probability less than 0.5 if you set the spreadsheet up right. For the other cells in row 2, put zeros because you can never reach a state where D>S. In all other cells, the probability that they will be occupied is the sum of (the probability of having a singleton added from the state immediately above) + (the probability of having a duo added from the state to the left), which for cell C3 would be B3*(($A3-B$1)/365)+C2*((365-$A2)/365) which can be copy/pasted to all other cells. Note that the cells only represent states with probability (S, D, 0). What happens if you have a triplet so T is greater than 0? If that happens, the probability is no longer accounted for in the spreadsheet. And that means you can take the sum of probabilities along a diagonal from lower left to upper right to find the probability that out of N people there are no triplets that share a birthday.

Alternatively, if you're a real geek, you can instead use the following perl code so you don't have to sum along a diagonal in a spreadsheet. The probability of having no triplets drops to 48.89% when you have 88 people.

#!/usr/bin/perl
use strict;
use warnings;

# The @vector is an array of arrays such that
# $vector[S][D] will hold the probability of getting
# to state (S,D).
# This initializes @vector and makes $vector[0][0] be 1.
my @vector = ();
$vector[0] = [1];

# $N is the total number of people in the room
my $N = 0;

# The variable $sum is the sum of all probabilities
# where N people are in the room and the state
# can be desribed as (S, D, 0), which is the
# probability that there are no triplets. The loop
# will stop when the $sum for a given value of $N is
# less than half. I would have initialized it within
# the do-while loop but then it wouldn't work as an
# until criterion :(
my $sum = 0;

do {
  $N++;
  $sum = 0;
  
  # Because of the way the array-of-arrays is set up,
  # this will have to push zeros to the end of previous
  # arrays when you add new potential values for D.
  # I'm overkilling with this by including impossibly
  # high values of D because I'm lazy.
  # Just ignore this unless you're a hardcore perl geek.
  for (my $i=0; $i<$N+1; $i++) {
    $vector[$i][$N] = 0;
  }
  
  # Make $vector[N][0] be $vector[N-1][0] times
  # (365-(N-1))
  $vector[$N][0] = $vector[$N-1][0] * (365-($N-1))/365;
  $sum += $vector[$N][0];
  
  # Next consider the $vector[N-D][D] terms
  for (my $D=1; $D < 1+$N/2; $D++) {
    my $S = $N-$D;
    $vector[$S][$D] = $vector[$S][$D-1] * ($S-($D-1))/365 +
                      $vector[$S-1][$D] * (365-($S-1))/365;
    $sum += $vector[$S][$D];
  }
  
} while ($sum > 0.5);
print "Stopped at N=$N with probability of no triplets $sum\n";

 

 

Edited December 20, 2015 by plasmid
Interesting; just using [code]blah blah[/code] tags doesn't work well, but ctrl-rightclick in a code box lets you copy/paste into a popup that works fine

[bonanova]

Nice.
Closed form would probably be ugly.

[Logophobic]

I didn't fully trust the results of my first simulation because I was using a pRNG with has a relatively short cycle. On Saturday, I dug up code for a better pRNG and repeated the simulation. The result was far more accurate, giving a probability of 49.95% (leap day excluded). I then began an iterative sum approach similar to plasmid's above. I manually calculated the probabilities for up to 8 persons, then began to code a routine to sum the probabilities for up to 99 persons. I didn't have time to finish, so plasmid beat me to the post. The results, of course, were exactly the same:

Spoiler

 86 persons: 48.7826%
 87 persons: 49.9455%
 88 persons: 51.1065%
 89 persons: 52.2648%

The code (Excel VBA)

Spoiler

If you want a detailed explanation of this code, just let me know.

Sub IterativeSum()
  Dim p(9999) As Double, t As Double
  Dim s As Long, d As Long
  Dim i As Long, n As Long
  
  p(0) = 1
  For n = 0 To 98
    For s = n To 0 Step -2
      d = (n - s) \ 2
      i = 100 * d + s
      If d > 0 Then t = t + p(i) * d / 365
      If s > 0 Then p(i + 99) = p(i + 99) + p(i) * s / 365
      p(i + 1) = p(i + 1) + p(i) * (365 - s - d) / 365
    Next s
    Sheet1.Cells(n + 1, 1) = t
  Next n
End Sub

 

Edited December 21, 2015 by Logophobic
fixing the code per plasmid's tip

[Logophobic]

A bit more probabilities for at least three persons sharing a birthday.

Spoiler

 102 persons: 66.697%
 111 persons: 75.461%
 117 persons: 80.524%
 132 persons: 90.119% 

And some probabilities for at least four persons sharing a birthday. Because: Why not?

Spoiler

 138 persons: 20.113%
 148 persons: 25.324%
 162 persons: 33.491%
 189 persons: 50.523%

 

[bonanova]

Since we're baring our code, here is the APL code for the 3bd-case:
 

Spoiler

 

 Z←GROUPSIZE BIRTHDAY TRIALS ;I;P
⍝ LIKELIHOOD OF 3 (OR MORE) PEOPLE IN A GROUP OF GROUPSIZE SHARING A BIRTHDAY
Z←0×I←1 
⍝
LOOP:
Z←Z+2<⌈/+/P∘.=P←?GROUPSIZE⍴365
→(TRIALS>I←I+1)/LOOP
Z←(+/Z)÷I  

GROUPSUZE of 87 returned just greater than 0.5 for 10,000 trials. ? is rng.

 

 

[Logophobic]

I realized an error in my code for calculating the four-person probabilities. For the sake of correctness, here are the true (and extended) results:

Spoiler

138 persons: 20.114%
148 persons: 25.343%
162 persons: 33.635%
187 persons: 50.269%
212 persons: 66.974%
226 persons: 75.318%
235 persons: 80.062%
260 persons: 90.215%