# Extending the birthday paradox to more than 2 people

In the traditional Birthday Paradox the question is "what are the chances that two or more people in a group of $n$ people share a birthday". I'm stuck on a problem which is an extension of this.

Instead of knowing the probability that two people share a birthday, I need to extend the question to know what is the probability that $x$ or more people share a birthday. With $x=2$ you can do this by calculating the probability that no two people share a birthday and subtract that from $1$, but I don't think I can extend this logic to larger numbers of $x$.

To further complicate this I also need a solution which will work for very large numbers for $n$ (millions) and $x$ (thousands).

• I presume that it's bioinformatics problem Aug 5, 2010 at 14:41
• It is actually a bioinformatics problem, but since it boils down to the same concept as the birthday paradox I thought I'd save the irrelevant specifics! Aug 5, 2010 at 14:50
• Normally I would agree with you, but in this case the specifics might matter since there could already be a bioconductor package that does what you ask. Aug 5, 2010 at 14:53
• If you really want to know, it's a pattern finding problem where I'm trying to accurately estimate the probability of a given level of enrichment of a subsequence within a set of larger sequences. I therefore have a set of subsequences with associated counts and I know how many subsequences I observed and how many theoretically observable sequences are available. If I saw a particular sequence 10 times out of 10,000 observations I need to know how likely that was to have occurred by chance. Aug 5, 2010 at 15:21
• Almost eight years later, I posted an answer to this problem at stats.stackexchange.com/questions/333471. The code there does not work for large $n,$ though, because it takes quadratic time in $n$.
– whuber
Mar 17, 2018 at 21:58

This is a counting problem: there are $b^n$ possible assignments of $b$ birthdays to $n$ people. Of those, let $q(k; n, b)$ be the number of assignments for which no birthday is shared by more than $k$ people but at least one birthday actually is shared by $k$ people. The probability we seek can be found by summing the $q(k;n,b)$ for appropriate values of $k$ and multiplying the result by $b^{-n}$.

These counts can be found exactly for values of $n$ less than several hundred. However, they will not follow any straightforward formula: we have to consider the patterns of ways in which birthdays can be assigned. I will illustrate this in lieu of providing a general demonstration. Let $n = 4$ (this is the smallest interesting situation). The possibilities are:

• Each person has a unique birthday; the code is {4}.
• Exactly two people share a birthday; the code is {2,1}.
• Two people have one birthday and the other two have another; the code is {0,2}.
• Three people share a birthday; the code is {1,0,1}.
• Four people share a birthday; the code is {0,0,0,1}.

Generally, the code $\{a, a, \ldots\}$ is a tuple of counts whose $k^\text{th}$ element stipulates how many distinct birthdates are shared by exactly $k$ people. Thus, in particular,

$$1 a + 2a + ... + k a[k] + \ldots = n.$$

Note, even in this simple case, that there are two ways in which the maximum of two people per birthday is attained: one with the code $\{0,2\}$ and another with the code $\{2,1\}$.

We can directly count the number of possible birthday assignments corresponding to any given code. This number is the product of three terms. One is a multinomial coefficient; it counts the number of ways of partitioning $n$ people into $a$ groups of $1$, $a$ groups of $2$, and so on. Because the sequence of groups does not matter, we have to divide this multinomial coefficient by $a!a!\cdots$; its reciprocal is the second term. Finally, line up the groups and assign them each a birthday: there are $b$ candidates for the first group, $b-1$ for the second, and so on. These values have to be multiplied together, forming the third term. It is equal to the "factorial product" $b^{(a+a+\cdots)}$ where $b^{(m)}$ means $b(b-1)\cdots(b-m+1)$.

There is an obvious and fairly simple recursion relating the count for a pattern $\{a, \ldots, a[k]\}$ to the count for the pattern $\{a, \ldots, a[k-1]\}$. This enables rapid calculation of the counts for modest values of $n$. Specifically, $a[k]$ represents $a[k]$ birthdates shared by exactly $k$ people each. After these $a[k]$ groups of $k$ people have been drawn from the $n$ people, which can be done in $x$ distinct ways (say), it remains to count the number of ways of achieving the pattern $\{a, \ldots, a[k-1]\}$ among the remaining people. Multiplying this by $x$ gives the recursion.

I doubt there is a closed form formula for $q(k; n, b)$, which is obtained by summing the counts for all partitions of $n$ whose maximum term equals $k$. Let me offer some examples:

With $b=5$ (five possible birthdays) and $n=4$ (four people), we obtain

\eqalign{ q(1) &= q(1;4,5) &= 120 \\ q(2) &= 360 + 60 &= 420 \\ q(3) &&= 80 \\ q(4) &&= 5.\\ }

Whence, for example, the chance that three or more people out of four share the same "birthday" (out of $5$ possible dates) equals $(80 + 5)/625 = 0.136$.

As another example, take $b = 365$ and $n = 23$. Here are the values of $q( k;23,365)$ for the smallest $k$ (to six sig figs only):

\eqalign{ k=1: &0.49270 \\ k=2: &0.494592 \\ k=3: &0.0125308 \\ k=4: &0.000172844 \\ k=5: &1.80449E-6 \\ k=6: &1.48722E-8 \\ k=7: &9.92255E-11 \\ k=8: &5.45195E-13. }

Using this technique, we can readily compute that there is about a 50% chance of (at least) a three-way birthday collision among 87 people, a 50% chance of a four-way collision among 187, and a 50% chance of a five-way collision among 310 people. That last calculation starts taking a few seconds (in Mathematica, anyway) because the number of partitions to consider starts getting large. For substantially larger $n$ we need an approximation.

One approximation is obtained by means of the Poisson distribution with expectation $n/b$, because we can view a birthday assignment as arising from $b$ almost (but not quite) independent Poisson variables each with expectation $n/b$: the variable for any given possible birthday describes how many of the $n$ people have that birthday. The distribution of the maximum is therefore approximately $F(k)^b$ where $F$ is the Poisson CDF. This is not a rigorous argument, so let's do a little testing. The approximation for $n = 23$, $b = 365$ gives

\eqalign{ k=1: &0.498783 \\ k=2: &0.496803\\ k=3: &0.014187\\ k=4: &0.000225115. }

By comparing with the preceding you can see that the relative probabilities can be poor when they are small, but the absolute probabilities are reasonably well approximated to about 0.5%. Testing with a wide range of $n$ and $b$ suggests the approximation is usually about this good.

To wrap up, let's consider the original question: take $n = 10,000$ (the number of observations) and $b = 1\,000\,000$ (the number of possible "structures," approximately). The approximate distribution for the maximum number of "shared birthdays" is

\eqalign{ k=1: &0 \\ k=2: &0.8475+\\ k=3: &0.1520+\\ k=4: &0.0004+\\ k\gt 4: &\lt 1E-6. }

(This is a fast calculation.) Clearly, observing one structure 10 times out of 10,000 would be highly significant. Because $n$ and $b$ are both large, I expect the approximation to work quite well here.

Incidentally, as Shane intimated, simulations can provide useful checks. A Mathematica simulation is created with a function like

simulate[n_, b_] := Max[Last[Transpose[Tally[RandomInteger[{0, b - 1}, n]]]]];

which is then iterated and summarized, as in this example which runs 10,000 iterations of the $n = 10000$, $b = 1\,000\,000$ case:

Tally[Table[simulate[10000, 1000000], {n, 1, 10000}]] // TableForm

Its output is

2 8503

3 1493

4 4

These frequencies closely agree with those predicted by the Poisson approximation.

• What a fantastic answer, thank you very much @whuber. Aug 21, 2014 at 20:21
• "There is an obvious and fairly simple recursion" — Namely? Oct 14, 2016 at 20:28
• @Kodiologist I inserted a brief description of the idea.
– whuber
Oct 14, 2016 at 20:56
• +1 but where in the original question did you see that n=10000 and b=1mln? The OP looks like it is asking about n=1mln and k=10000, with b unspecified (presumably b=365). Not that it matters at this point :) Oct 14, 2016 at 21:20
• @amoeba After all this time (six years, 1600 answers, and closely reading tens of thousands of posts) I cannot recall, but most likely I misinterpreted the last line. In my defense, note that if we read it literally the answer is immediate (upon applying a version of the Pigeonhole Principle): it is certain that among $n$=millions of people there will be at least one birthday that is shared among at least $x$=thousands of them!
– whuber
Oct 14, 2016 at 21:25

It is always possible to solve this problem with a monte-carlo solution, although that's far from the most efficient. Here's a simple example of the 2 person problem in R (from a presentation I gave last year; I used this as an example of inefficient code), which could be easily adjusted to account for more than 2:

birthday.paradox <- function(n.people, n.trials) {
matches <- 0
for (trial in 1:n.trials) {
birthdays <- cbind(as.matrix(1:365), rep(0, 365))
for (person in 1:n.people) {
day <- sample(1:365, 1, replace = TRUE)
if (birthdays[birthdays[, 1] == day, 2] == 1) {
matches <- matches + 1
break
}
birthdays[birthdays[, 1] == day, 2] <- 1
}
birthdays <- NULL
}
print(paste("Probability of birthday matches = ", matches/n.trials))
}

• I am not sure if the multiple types solution will work here.
– user28
Aug 5, 2010 at 14:32
• I think that generalisation still only works for 2 or more people sharing a birthday - just that you can have different sub-classes of people. Aug 5, 2010 at 14:45

This is an attempt at a general solution. There may be some mistakes so use with caution!

First some notation:

$P(x,n)$ be the probability that $x$ or more people share a birthday among $n$ people,

$P(y|n)$ be the probability that exactly $y$ people share a birthday among $n$ people.

Notes:

1. Abuse of notation as $P(.)$ is being used in two different ways.

2. By definition $y$ cannot take the value of 1 as it does not make any sense and $y$ = 0 can be interpreted to mean that no one shares a common birthday.

Then the required probability is given by:

$P(x,n) = 1 - P(0|n) - P(2|n) - P(3|n) .... - P(x-1|n)$

Now,

$P(y|n) = {n \choose y} (\frac{365}{365})^y \ \prod_{k=1}^{k=n-y}(1 -\frac{k}{365})$

Here is the logic: You need the probability that exactly $y$ people share a birthday.

Step 1: You can pick $y$ people in ${n \choose y}$ ways.

Step 2: Since they share a birthday it can be any of the 365 days in a year. So, we basically have 365 choices which gives us $(\frac{365}{365})^y$.

Step 3: The remaining $n-y$ people should not share a birthday with the first $y$ people or with each other. This reasoning gives us $\prod_{k=1}^{k=n-y}(1 -\frac{k}{365})$.

You can check that for $x$ = 2 the above collapses to the standard birthday paradox solution.

• Will this solution suffer from the curse of dimensionality? If instead of n=365, n=10^6 is this solution still feasible? Aug 5, 2010 at 15:38
• Some approximations may have to be used to deal with high dimensions. Perhaps, use Stirling's approximation for factorials in the binomial coefficient. To deal with the product terms you could take logs and compute the sums instead of the products and then take the anti-log of the sum.
– user28
Aug 5, 2010 at 15:45
• There are also several other forms of approximations possible using for example the Taylor series expansion for the exponential function. See the wiki page for these approximations: en.wikipedia.org/wiki/Birthday_problem#Approximations
– user28
Aug 5, 2010 at 16:15
• Suppose y=2, n=4, and there are just two birthdays. Your formula, adapted by replacing 365 by 2, seems to say the probability that exactly 2 people share a birthday is Comb(4,2)*(2/2)^2*(1-1/2)*(1-2/2) = 0. (In fact, it's easy to see--by brute force enumeration if you like--that the probabilities that 2, 3, or 4 people share a "birthday" are 6/16, 8/16, and 2/16, respectively.) Indeed, whenever n-y >= 365, your formula yields 0, whereas as n gets large and y is fixed the probability should increase to a non-zero maximum before n reaches 365*y and then decrease, but never down to 0.
– whuber
Aug 18, 2010 at 15:05
• @svadali I did not replace 365 by $n$; I replaced it by 2, which gives a simple problem we can check by hand. At any rate, Step 2 has a typo and the reasoning in Step 3 is incorrect. In fact, the approach is flawed because it does not consider various complications that arise with multiple "collisions," such as two people sharing one birthday and another two sharing another birthday. I get the value of 6/16 (four people, two possible birthdays, at most two per birthday) because of the 16 possible configurations, 6 have two people share one birthday and the others sharing the other birthday.
– whuber
Feb 22, 2011 at 17:06