3
$\begingroup$

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

enter image description here

$\endgroup$
5
  • 1
    $\begingroup$ Related: math.stackexchange.com/questions/177692/…. $\endgroup$
    – JimB
    Commented May 4, 2023 at 16:01
  • 1
    $\begingroup$ 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. $\endgroup$ Commented May 4, 2023 at 16:15
  • $\begingroup$ 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.") $\endgroup$ Commented May 4, 2023 at 16:55
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Commented May 4, 2023 at 22:17
  • $\begingroup$ @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 $\endgroup$ Commented May 5, 2023 at 5:50

2 Answers 2

0
$\begingroup$

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]]]

Minimum values of samples to achieve at least a 50% chance for at least one match

$\endgroup$
3
$\begingroup$

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]

Results from simulations

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]

True minimum numbers of people to be selected for each of the 100 scenarios

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Commented May 6, 2023 at 2:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.