7
$\begingroup$

Edit high-level summary of problem/solutions in Mathematica community post


Suppose I draw (with replacement) $k$ colored balls from bins with probabilities $h_1,\ldots,h_d$ until some color is drawn more than once.

What's the fastest way to get expected value of $k$?

Adapting JimB's code from this answer is below, wondering if it can be made much faster. Exact solution using AugmentedSymmetricPolynomials is too slow to be practical for this problem size.

ClearAll["Global`*"];
d = 10000;
p = 1.;
h = Table[i^(-p), {i, d}];
h = h/Total[h];

meanTimeToCollision = Compile[{{h2, _Real, 1}},
   Module[{s, x, mean, k, nsim},
    h = h2/Total[h2];
    nsim = 1000;
    x = RandomChoice[h -> Range[Length[h]], {Floor[nsim], d + 1}];
    mean = 0.;
    Do[k = 2;
     While[
      Select[x[[i, 1 ;; k - 1]], # == x[[i, k]] &] == {} && 
       k < (d + 1), k = k + 1];
     mean = mean + k,
     {i, nsim}
     ];
    mean = mean/nsim
    ]
   ];
meanTimeToCollision[h] // Timing

This code lets me estimate these times for various distributions up to $d=1000$ reasonably fast, although I'm wishing for a 10x faster procedure to plot it for $d=10000$ where asymptotics are more obvious

enter image description here

Notebook

$\endgroup$
11
  • $\begingroup$ Naive question: Is this problem in some way different from your recent weighted birthday question? (I have no objection if they’re the same, just checking whether my sanity left town for the weekend.) $\endgroup$ Commented May 6, 2023 at 15:43
  • $\begingroup$ @Daniel Lichtblau: Hi Danny, Doesn’t the birthday problem have (approximately) uniform probabilities of \sigmeq 1/365? $\endgroup$ Commented May 6, 2023 at 15:53
  • $\begingroup$ Hey Craig, yes it does. I should add a link to the prior question, which posited nonuniform probabilities. Except I’m doing this on a phone, and adding links is…real work (not for young people, of course). $\endgroup$ Commented May 6, 2023 at 16:00
  • $\begingroup$ OK, here’s that link (maybe, if this works…). $\endgroup$ Commented May 6, 2023 at 16:02
  • 1
    $\begingroup$ @wolfies thanks for the comment (I have your book!). $R$ is the "effective number of bins", and is equal to $d$ when $p_1=p_2=\ldots=p_d$, the underlying motivation is is given here for a related problem. I used term "collision" because of similarity of this problem to "hash collision" problem $\endgroup$ Commented May 7, 2023 at 9:07

5 Answers 5

3
$\begingroup$

Edit: Added an estimate of the standard error of the sample mean because an estimate without a measure of precision is at best of unknown value.

Because a "collision" will typically happen long before $d=10,000$, initially obtaining a sample of smaller size than 10,000 followed by a supplementary sample if necessary can considerably reduce the computation time.

ClearAll["Global`*"];
d = 10000;
p = 1.;
h = Table[i^(-p), {i, d}];
h = h/Total[h];
meanTimeToCollision = Compile[{{h2, _Real, 1}}, 
   Module[{s, x, mean, k, nsim, threshold, y, se, k2}, 
    h = h2/Total[h2];
    threshold = 50;
    nsim = 1000;
    x = RandomChoice[h -> Range[Length[h]], {Floor[nsim], Min[threshold, d + 1]}];
    mean = 0.;
    k2 = 0;
    Do[k = 2;
     While[
      Select[x[[i, 1 ;; k - 1]], # == x[[i, k]] &] == {} && 
       k < threshold, k = k + 1];
     If[k == threshold && Select[x[[i, 1 ;; k - 1]], # == x[[i, k]] &] == {}, 
      y = RandomChoice[h -> Range[Length[h]], Floor[d + 1 - threshold]];
      x[[i]] = Join[x[[i]], y];
      While[Select[x[[i, 1 ;; k - 1]], # == x[[i, k]] &] == {} && 
        k < (d + 1), k = k + 1]];
     mean = mean + k;
     k2 = k2 + k^2, {i, nsim}];
    mean = mean/nsim;
    se = Sqrt[(k2 - nsim*mean^2)/(nsim (nsim - 1))];
    {mean, se}]];
meanTimeToCollision[h] // Timing
(* {0.015625, {11.857, 0.201656}} *)

I've arbitrarily chosen 50 as the initial sample size threshold. Something else might be better. And maybe sequential thresholds might make it more efficient.

$\endgroup$
7
  • $\begingroup$ can you reality check my mn function. Essentially it just draws from distribution till encounters something already drawn and calculates length of draws. Any way +1 $\endgroup$
    – ubpdqn
    Commented May 8, 2023 at 18:36
  • $\begingroup$ @ubpdqn Gladly. But first I need to figure out NestWhile as I've not used that function before. I'll also check it against the exact values. $\endgroup$
    – JimB
    Commented May 9, 2023 at 3:27
  • $\begingroup$ Thanks. NeetWhile is continued f(f(…f(x)) while condition met. Stop when not met. Outputs ‘nesting’ were condition not met. A better condition than my messy Tally is: Length[#]==Length[Union@#]&. This iterates while no duplicates in output…@kglr neat trick $\endgroup$
    – ubpdqn
    Commented May 9, 2023 at 3:42
  • $\begingroup$ Thanks for the tips, this was the fasted method for $d<100$ (also using MKL as random number generator), kglr's solution was faster for d>1000, community post $\endgroup$ Commented May 9, 2023 at 12:47
  • $\begingroup$ Because @kglr gave the faster approach for your objective of large $n$, that answer should be the one accepted. $\endgroup$
    – JimB
    Commented May 9, 2023 at 16:19
4
$\begingroup$
SeedRandom[Method -> "MersenneTwister"];
d = 10000;
p = 1.;
h = Table[i^(-p), {i, d}];
h = h/Total[h];

Note that, with nsim = d = 10000, generating the random 10000$\times$10000 matrix already takes more than half of the total processing time in OP's meanTimeToCollision.

SeedRandom[1];
RandomChoice[h -> Range[10000], {10000, 1 + 10000}]; // Timing // First
6.41687

In the following, singledraw generates and processes a single sample RandomChoice[h -> Range[10000], 1+10000]. We use singledraw with ParallelTable to get the desired speed-up:

ClearAll[singledraw]
singledraw = Module[{rc = RandomChoice[# -> Range[Length[#]], Length[#] + 1], 
     ctime = 0., k = 1}, 
  While[++k == Length[Union@rc[[;; k]]]]; ctime += k] &;

SeedRandom[1];
singledraw[h] // Timing
{0.004907, 9.}

SeedRandom[1];
sample = ParallelTable[singledraw[h], 10000, ProgressReporting -> False]; //
   Timing // First
 0.740078
Mean @ sample
 12.1707
dist = FindDistribution[sample - 2]
MixtureDistribution[{0.705792, 0.294208}, 
 {NegativeBinomialDistribution[3, 0.296478], 
  NegativeBinomialDistribution[8, 0.327494]}]

Show[Histogram[sample, {1}, "PDF"], 
 Plot[PDF[SmoothKernelDistribution[sample, 
    "Oversmooth", {"Bounded", {2, \[Infinity]}, "Gaussian"}], x], 
   {x, 2, 50}], 
 DiscretePlot[PDF[dist, x - 2], {x, 2, 50}, 
  PlotStyle -> Directive[Red, PointSize[Medium]]], 
 ImageSize -> Large]

enter image description here

$\endgroup$
8
  • 1
    $\begingroup$ clever use of Union to detect prior colour drawn. +1 Much better than my use of Tally. $\endgroup$
    – ubpdqn
    Commented May 8, 2023 at 19:00
  • 1
    $\begingroup$ I find similarity of results reassuring: see here $\endgroup$
    – ubpdqn
    Commented May 8, 2023 at 19:22
  • $\begingroup$ @ubpdqn, i think we are using the same idea with different implementations. $\endgroup$
    – kglr
    Commented May 8, 2023 at 19:25
  • 2
    $\begingroup$ this is the fastest solution for $d=10000$, about 3x faster on my laptop. For $d=100$ the other solution is faster. Summarized in community post $\endgroup$ Commented May 9, 2023 at 12:48
  • 1
    $\begingroup$ +1 I do have a minor suggestion about the figure. As the positive probabilities start at 2, the following gives a slightly better match to the data: dist = FindDistribution[sample - 2]; Show[Histogram[sample, {1}, "PDF"], Plot[PDF[SmoothKernelDistribution[sample, "Oversmooth", {"Bounded", {2, \[Infinity]}, "Gaussian"}], x], {x, 2, 50}], DiscretePlot[PDF[dist, x - 2], {x, 2, 50}, PlotStyle -> Directive[Red, PointSize[Medium]]], ImageSize -> Large]. $\endgroup$
    – JimB
    Commented May 10, 2023 at 17:04
3
$\begingroup$

Updated answer Refer to functions in original answer.

     pl = ListPlot[
Table[{1/(prob[j].prob[j]), mn[1000,
prob[j]],  {j,
Join[Range [10, 500, 10], Range [1000,
10000, 1000]]}],
PlotRange -> Full]

Using R as defined in hyperlink:

 Show[pl, ListPlot[Table[{1 / (prob[j].
prob[j]), Sqrt[Pi / (2 prob[j].
prob[j])}, (j, 10, 10000, 10}1,
Joined -› True, PlotRange -> All]]

enter image description here

Interpreting question as number of draws till two sequential balls of same colour are drawn:

If my interpretation is still incorrect please let me know.

Setting up

prob[n_] := Module[{h = Table[1./i, {i, n}]},
  h = h/Total[h]]
f[x_, p_] := Append[x, RandomChoice[p -> Range[Length[p]]]]

Example draws:

 example[p_] := 
 NestWhile[f[#, p] &, Nest[f[#, p] &, {}, 2], #[[-1]] != #[[-2]] &]

For example: example[prob[10000]] yields: {235, 3, 67, 8, 204, 2486, 1785, 1668, 2406, 89, 6, 309, 60, 120, 1,
386, 1, 2, 92, 1, 21, 4329, 257, 19, 2, 26, 3, 4326, 1317, 540, 2648,
14, 2224, 2, 1, 1}

Simulating:

sim[n_, p_] := 
 Mean[Table[
    Length[NestWhile[f[#, p] &, 
      Nest[f[#, p] &, {}, 2], #[[-1]] != #[[-2]] &]], n]] // N

Sample sizes of 1000:

ListPlot[
 Table[{j, sim[1000, prob[j]]}, {j, 
   Join[Range[10, 500, 10], Range[1000, 10000, 1000]]}], 
 PlotRange -> Full]

Note this does not follow the function in OP except for small values of d.

enter image description here Original answer

This original answer was number of draws till a ball of same colour than previous draws.

Generalizing probability distribution:

prob[n_] := Module[{h = Table[1./i, {i, n}]},
  h = h/Total[h]]

Calculating means from simulation using NestWhile

f[x_, p_] := Append[x, RandomChoice[p -> Range[Length[p]]]]
mn[n_, p_] := 
 Mean[Table[
    Length[NestWhile[f[#, p] &, {}, Max[Tally[#][[All, 2]]] <= 1 &]], 
    n]] // N

Simulations using 1000 samples:

ListPlot[
 Table[{j, mn[1000, prob[j]]}, {j, 
   Join[Range[10, 500, 10], Range[1000, 10000, 1000]]}], 
 PlotRange -> Full]

enter image description here

$\endgroup$
6
  • $\begingroup$ Hm, didn't mean consecutive, consecutive same color seems like a simpler problem, rephrased the original question for clarity $\endgroup$ Commented May 8, 2023 at 8:34
  • $\begingroup$ @YaroslavBulatov see update. My original answer was the case of number of draws till drawing a colour already drawn. Using your R definition I have updated. You can of course use larger sample sizes. $\endgroup$
    – ubpdqn
    Commented May 8, 2023 at 9:29
  • $\begingroup$ hm, trying your code for the example in original post gives me something that's off by factor of 4 from true answer $\endgroup$ Commented May 8, 2023 at 12:11
  • $\begingroup$ @YaroslavBulatov use function mn NOT sim. If you read my answer. I first answered question based on number of draws till previous colour: mn. If you look at code sim it is clearly number of draws till consecutive same colour…which of course will be longer. As you ignored my original answer I assumed I interpreted incorrectly. If you look at ORIGINAL answer, I suggest re-run: mn. $\endgroup$
    – ubpdqn
    Commented May 8, 2023 at 17:06
  • $\begingroup$ Running mn[100000,prob[10000]] yields 12.2641. Closer to your ‘true’ answer. Note this a sample size of 100000. $\endgroup$
    – ubpdqn
    Commented May 8, 2023 at 17:18
2
$\begingroup$

Here is a different strategy. I am not sure if it is faster or not, but put it out there to consider.

Create a random sequence with weights:

roll[number_, length_] := 
 RandomChoice[1/Range[1, number] -> Range[1, number], length]

Look for sequential pairs by rotating a list of a given size:

RepeatedTiming[
With[{try = roll[100000, 1000] },
 FirstPosition[Most[Transpose[{try, RotateLeft[try]}]], {a_, a_}]
 ]
]

About 0.1 seconds.

If the number of colors is very large then use sequential rolls

findFirstPair[number_, trialLength_] := 
 Block[{count = -1, position = Missing["NotFound"]}, 
  While[position == Missing["NotFound"], count++;
   position = 
    With[{try = roll[number, trialLength]}, 
     FirstPosition[
      Most[Transpose[{try, RotateLeft[try]}]], {a_, a_}]]];
  count trialLength + position[[1]]
  ]

For example:

findFirstPair[1000000, 100]

Find a mean of 100 attempts

AbsoluteTiming[
 Mean[Table[findFirstPair[1000000, 1000], {100}]]
 ]

About 10 seconds. I didn't try to find the optimal trialLength for a given number of colors.

$\endgroup$
2
$\begingroup$
SeedRandom[Method -> "MersenneTwister"];
d = 10000; p = 1.; h = Normalize[Table[i^(-p), {i, d}], Total];

We process RandomChoice[h -> Range[Length[h]]] output to transform it to a draw from MultinomialDistribution[1, h].

(Somehow, more direct approach using RandomVariate[MultinomialDistribution[1, h]] is much slower).

To find the time to first collision we accumulate draws until the maximum of accumulated draws becomes 2:

ClearAll[oneHot, drawMND, ttcMND]

oneHot = ReplaceAll[x_Integer :> UnitVector[Length @ h, x]];

drawMND := oneHot @ RandomChoice[h -> Range[Length[h]]];

ttcMND := Module[{ttc = 0., k = 1, r0 = drawMND}, 
  While[Max[r0 += drawMND] < 2, k++]; ttc += k; ttc]

SeedRandom[1];
ttcMND // Timing
{0.002681, 8.}

Use ttcMND in ParallelTable to generate a sample of size 10000:

SeedRandom[1];
sampleMND = ParallelTable[ttcMND, 10000, ProgressReporting -> False]; // Timing
{0.746258, Null}
Through[{Min, Mean, Median, Max} @ sampleMND]
{1., 11.2721, 10., 49.}
dst = FindDistribution[sampleMND]
NegativeBinomialDistribution[3, 0.213528]
Show[Histogram[sampleMND, {1}, "PDF"], 
 SmoothHistogram[sampleMND], 
 DiscretePlot[PDF[dst, x], {x, 0, 50}, 
  PlotStyle -> Directive[Red, PointSize[Medium]]], 
 ImageSize -> Large]

enter image description here

$\endgroup$

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.