# Non-uniform birthday problem in Mathematica

Suppose my year has $$d$$ days. I have a vector of probabilities of being born on those days as $$h=p_1,\ldots,p_d$$, and want to estimate the smallest $$k$$ such that probability of of birthday collision (2 or more people sharing the same birthday), is at least 50%

What's the fastest way of doing this in Mathematica? Function below birthdayThreshold[h] below estimates this by sampling from MultinomialDistribution, but this approach is very slow, what's the best way to speed this up? In particular for large $$d \gtrsim 10000$$

Related discussion

ClearAll["Global*"];
spectrum2R[h_] := Total[h]^2/Total[h^2];
birthdayProbabilities[h2_] := (
h = h2/Total[h2];
maxLength = Ceiling[2*Sqrt[spectrum2R[h]]];
dist = MultinomialDistribution[1, h];
data = RandomVariate[dist, {numSamples, maxLength}];
collisionSeries[drawSeries_] :=
Total[UnitStep[# - 2]] & /@ Accumulate@drawSeries;
collisionCounts = collisionSeries /@ data; (*
numSamples x maxLength *)
collisionBoole = N@Map[UnitStep[# - 1] &, collisionCounts, {2}];
Mean[collisionBoole]
);
birthdayThreshold[h_] := (
probs = birthdayProbabilities[h];
First@Flatten@Position[# > .5 & /@ probs, True, 1, 1]
);

normalize[h_] := h/Total[h];
numSamples = 100;
d = 100;
plist = Table[1 + 1./r, {r, 1, d, d/100}];
dists = Table[normalize[Table[i^-p, {i, 1, d}]], {p, plist}];
ListPlot[birthdayThreshold /@ dists] // Timing


• – JimB
Commented May 4, 2023 at 16:01
• Hint: Surely you order the days by descending probability of carrying birthdays. Then, compute a hypergeometric distribution with successively larger sets of birthdays until the criterion is reached. Commented May 4, 2023 at 16:15
• I don't see why there needs to be a solution, in general. Suppose there are $d = 3$ days and $P_i = 1/3$. There is no $k>1$ for which "the smallest $k$ such that probability of $k$ people sharing the same birthday is at least 50%." (I presume the $k=1$ case is excluded as it makes little sense to say a person 'shares' his own birthday.") Commented May 4, 2023 at 16:55
• I'm probably being thick: "the smallest $k$ such that probability of $k$ people sharing the same birthday is at least 50%". Should that be "the smallest number of randomly selected people such that the probability of exactly $k$ people sharing the same birthday is at least 50%". In other words, I'm confused by the use of $k$ twice in the same sentence for what seems to be two different meanings.
– JimB
Commented May 4, 2023 at 22:17
• @JimB doh, that's bad wording on my part. Classical birthday problem asks for size of group such that 2 or more people in the group are likely to share the same birthday. I wanted to know the extension of this to non-uniform day probabilities Commented May 5, 2023 at 5:50

No simulations are needed as it's faster to obtain the (usually) exact result. By "usually" I mean using machine precision numbers for the multinomial probability vectors as opposed to rationalizing those probabilities. (My previous answer used Rationalize for the exact solution and sped up the simulation solution but it turns out neither of those are needed.)

Given a multinomial probability vector the following function returns the smallest number of samples ($$k$$) needed to have at least a 50% chance of observing two samples from the same category.

minK50[p_] := Module[{prob, k},
prob = AugmentedSymmetricPolynomial[{2}, p];
k = 2;
While[prob < 1/2,
k = k + 1;
prob = prob + (k - 1) AugmentedSymmetricPolynomial[Join[{2}, ConstantArray[1, k - 2]], p]];
{k, prob}
]


Using your set-up of 100 different multinomial probability vectors but with d = 10000:

normalize[h_] := h/Total[h];
d = 10000;
plist = Table[1 + 1./r, {r, 1, d, d/100}];
dists = Table[normalize[Table[i^-p, {i, 1, d}]], {p, plist}] // N;

(t = Table[{i, minK50[dists[[i]]]}, {i, Length[dists]}];) // Timing
(* {8.625, Null} *)

ListPlot[t[[All, 2, 1]]]


Here's an approach to estimate the minimum number of people needed to have at least $$k$$ people with the same birthday with a probability of at least 0.5 (which does not exactly coincide with how you've currently described the problem).

This ends up being about 20 times faster than your current code. Specifically it gets 10 times the simulations in about half of the time of your code.

normalize[h_] := h/Total[h];
numSamples = 1000;
d = 100;
plist = Table[1 + 1./r, {r, 1, d, d/100}];
dists = Table[normalize[Table[i^-p, {i, 1, d}]], {p, plist}];

k = 4; (* k is number of selected people with the same birthday *);

SeedRandom[12345];
Timing[
results = ConstantArray[{0, 0}, Length[dists]];
Do[
counts = ConstantArray[0, Length[dists[[i]]]];
Do[x = RandomChoice[dists[[i]] -> Range[Length[dists[[i]]]],  Length[dists[[i]]]];
(* Find the number of people selected where k people have the same birthday *)
m = 1;
While[Max[Tally[x[[1 ;; m]]][[All, 2]]] < k, m = m + 1];
Do[counts[[j]] = counts[[j]] + 1, {j, m, Length[dists[[i]]]}], {n, 1, numSamples}];
counts = counts/numSamples;
results[[i]] = Select[Transpose[{Range[Length[dists[[i]]]], counts}], #[[2]] >= 1/2 &][[1]],
{i, 1, Length[dists]}]
]
(* {12.3594, Null} *)

ListPlot[results[[All, 1]], PlotRange -> All]


For k = 2 (and other values of k) one can find the minimum values and the associated probabilities exactly but in about the same timing as your code.

results = ConstantArray[{0, 0}, Length[dists]];
Do[
dd = Rationalize[dists[[i]], 0];
prob = AugmentedSymmetricPolynomial[{2}, dd];
m = 2;
While[prob < 1/2,
m = m + 1;
prob =
prob + (m - 1) AugmentedSymmetricPolynomial[Join[{2}, ConstantArray[1, m - 2]], dd]
];
results[[i]] = {m, prob // N},
{i, 100}] // Timing
(* {21.8281, Null} *)
ListPlot[results[[All, 1]], PlotRange -> All]


• I'm going to modify my answer by removing the Rationalize part because it doesn't appear to be needed when using AugmentedSymmetricPolynomial (but one does need it if those augmented symmetric polynomials are converted with PowerSymmetricPolynomial. Also, without the Rationalize, it works relatively quickly for d = 10000`.
– JimB
Commented May 6, 2023 at 2:22