Birthday “Paradox” — with a different perspective

Background: Many people are familiar with the so-called Birthday "Paradox" that, in a room of 23 people, there is a better than 50/50 chance that two of them will share the same birthday. In its more general form for n people, the probability of no two people sharing the same birthday is $p(n) = \frac{365!}{ 365^n *(365-n)!}$. Similar calculations are used for understanding hash-space sizes, cryptographic attacks, etc.

Motivation: The reason for asking the following question is actually related to understanding a specific financial market behavior. However a variant on the "Birthday Paradox" problem fits exactly as an analogy and is likely to be of wider interest to people with different backgrounds. My question is therefore framed along the lines of the more familiar "Birthday Paradox", as follows.

Situation: There are a total of 60 people in a room. Of these, it turns out that there are 11 (eleven) PAIRS of people who share the same birthday, and one TRIPLE (i.e. group of 3 people) who have the same birthday. The remaining 60 - 2*11 - 3 = 35 people have different birthdays. Assuming a population in which any day is equally likely for a birthday (i.e. ignore Feb 29th & possible seasonal effects) and, given the specified distribution of birthdays, the questioner would like to understand how likely (or unlikely) it is that these 60 people were really chosen at random. This question was originally posed on another site where it was left unanswered, but the questioner was advised to re-state the question in the form that now follows below.

Question: "If 60 people are chosen at random from a population in which any day is equally likely to be a person's birthday, what is the probability that there are 11 days on which exactly 2 people share a birthday, one day on which exactly 3 of them share a birthday, and no days on which 4 or more people share a birthday?"

When in doubt, simulate. (I'm sure that you can actually put together a formula, but it will likely be painful to look at.) I'll work with the criterion being "at least one triple birthday" and "at least eleven double birthdays" -- changing the code below to checking for exactly so many birthdays is not hard.

n.sims <- 1e5 n.persons <- 60

counter <- 0
pb <- winProgressBar(max=n.sims)
for ( ii in 1:n.sims ) {
setWinProgressBar(pb,ii,paste(ii,"of",n.sims))
set.seed(ii)
birthdays <- sample(x=365,size=n.persons,replace=TRUE)
birthday.table <- table(birthdays)
if (    sum(birthday.table>=4) == 0  &
sum(birthday.table==3) >= 1  &
sum(birthday.table==2) >= 11 ) counter <- counter+1
}
close(pb)
counter/n.sims


Out of 100,000 simulations, I get six hits, for a $p$-value of $p=0.00006$. If you twiddle the set.seed() command, e.g., to set.seed(2*ii), you will get slightly different results (in this case $p=0.00004$), which serves as a sort of sensitivity analysis.

• (1) In what sense is this a "p-value"? There is no hypothesis test in evidence. (2) Your results are based only on 4 or 6 observations out of a large simulation. A good approximation is the Poisson; accordingly there's more than a 5% chance that the correct proportion is $7.5/10^5$ or greater and there is more than a 5% chance that the correct answer is $2.5/10^5$ or less. Therefore your simulations identify the answer only to within one half order of magnitude--not even one significant digit! You would need to make this 100 times bigger to identify even the first digit of $p$ reliably. – whuber Nov 7 '17 at 15:46
• @whuber: (1) It's a p value because the OP explicitly asked about the probability of getting an extreme result under the null hypothesis of independently picking people with uniformly distributed birthdays. (2) Point taken that you'd need a lot more simulations to get a "serious" result. Note that I did recommend looking at the sensitivity. – S. Kolassa - Reinstate Monica Nov 7 '17 at 15:52
• That's not what a p-value is, even though it sounds similar. A p-value is a number attached to a sample used to assess how consistent it is with a null hypothesis. You are just computing a probability. To establish a p-value in this circumstance you would have to propose a null hypothesis, construct a critical region, and then compute the chance of that critical region under the null--which you are definitely not doing. – whuber Nov 7 '17 at 16:15
• The sample is the hypothetical group of people with one "triple birthday" and at least eleven "double birthdays". I am assessing how consistent this is with the null hypothesis of having picked 60 people at random, with uniformly distributed birthdays. The critical region is the quadrant $x\geq 11$ and $y\geq 1$ in the "double/triple birthday space". I don't fully understand why the simulated probability is not a p-value. Where am I mistaken? – S. Kolassa - Reinstate Monica Nov 7 '17 at 16:20
• Yes, I believe that is mistaken. To compute a legitimate p-value you would also need to define what you mean by "extreme or more extreme" than this event. Such a definition of "more extreme" is not at all apparent and cannot be provided without specifying the alternative hypothesis, at which point it's likely the p-value you compute would be much greater than the number you are producing here. The error made in calling this a "p-value" is comparable to mistakenly computing the p-value of a Binomial outcome $X$ as the probability of $X$ itself. – whuber Nov 7 '17 at 16:30

Assuming you aim for the case that exactly (so not at least) 11 double birthdays and 1 triple birthdays occur:

\begin{align}p(60,n_2=11,n_3=1) &= \frac{\text{possibilities with 11 double birthdays and 1 triple birthdays}}{\text{all possibilities}}\\ &= \frac{\frac{60!}{35!22!3!} 21!! \,\cdot \,365 \cdot 364 \cdot \, ... \, \cdot (365-n+1+1\cdot11+2\cdot1)}{365^{60}}\end{align}

which my calculation in R

factorial(60)/(factorial(35)*factorial(22)*factorial(3))*
pracma::factorial2(21)*
cumprod(319:365)[47]/365^60


approximates as $3.64 * 10^{-5}$

Explanation of the terms in the calculation

• $365^{60}$ is the number of ways to select random 60 birthdays among equally probable 365 days.
• $365 \cdot 364 \cdot \, ... \, \cdot (365-n+1+1\cdot11+2\cdot1)$ is the number of ways to select random 60-11-2 unique days among equally probable 365 days.
• $\frac{60!}{35!22!3!}$ is the number of unique ways to partition 60 people in groups of 35, 22, and 3 (the numbers of people with single birthdays, double birthdays, and triple birthdays)
• $21!! = 21 \cdot 19 \cdot ... \cdot 3 \cdot 1$ is the number of ways that we can partition the 22 people in the double birthdays group into pairs.
• +1 this is the right answer. Just note that the OP specifically asked to ignore feb 29 (meaning there's 365 days within a year), and is $21!$ really the number of ways to partition 22 people into pairs or is it $22!$? – IWS Nov 7 '17 at 13:50
• I have changed it accordingly. I had ignored feb 29 in the sense to make it 'equally probable' instead of absent. The statement to ignore 29 feb is a bit ambiguous, but indeed in discussing the original birthday problem the OP mentions 365. – Sextus Empiricus Nov 7 '17 at 13:55
• It is 21!! (the double factorial). The first person of the pairs is a 'free' choice. So you select, without loss of generality, the first person in the first pair and have (22-1) options left to make a first pair, then (22-3) options to make the second pair, (22-5) options to make the third pair, etc. Alternatively you take 22! the number of ways to order the people and divide by $2^{11}$ and 11!, to correct for the number of ways to switch within and between the pairs (these have been counted multiple times by using 22!, the double factorial expression is more direct). – Sextus Empiricus Nov 7 '17 at 15:05
• Thanks for explaining! I did not really know/notice the double $!!$, and so I learned something useful today. – IWS Nov 7 '17 at 16:50
• Your question on the double !! made me think of a way to simplify the double factorial term which allows to create a simpler expression and opens the way for a, not so bad looking, general expression: $$p(n_1,n_2,n_3...n_i) = \frac{(n_t)! (365)!}{(365)^{n_t}} \frac{1}{(365-n_s)!} \prod \frac{1}{(i!)^{n_i}(n_i)!}$$ with $n_t = \sum i \cdot n_i$ the number of people and $n_s = \sum n_i$ the number of days that are a birthday. I imagine that, when we take the log and use Stirling's approximation, that we may get to an optimization formula that leads to the most probable state. – Sextus Empiricus Nov 7 '17 at 17:12