I'd like to get the probability that some N samples of M equiprobable events does not have any case of the M possible individual outcomes having more than some maximum value.

The code I'm using on Mathematica 11 is:

Probability[Max[Array[x, bins]] <= max, Array[x, bins] \[Distributed] 
             MultinomialDistribution[balls, ConstantArray[1/bins, bins]]]

But, this gets very slow with cases I'm interested in. For example:

{balls, bins, max} = {100, 10, 20}

has been running for several minutes on my reasonably fast 4-core Intel laptop.

Can the performance be improved in any way to get cases like this or similar magnitude of values to finish in less time?


1 Answer 1


Here's an old function I cooked up for this very scenario: PMF of the maximum of categories in a multinomial with equal category probabilities:

pMax[balls_, bins_, max_, p_: Infinity] := Module[{h},
   h[0, n_, m_] = N[1, p];
   h[s_, n_, m_] := h[s, n, m] = Sum[(n*x + x - s) (h[s - x, n, m]/x!), {x, Min[s, m]}]/s;
   balls! h[balls, bins, max]/bins^balls];

{balls, bins, max} = {50, 7, 10};

r1 = pMax[balls, bins, max]; // AbsoluteTiming // First

r2 = Probability[Max[Array[r, bins]] <= max, 
     Array[r, bins] \[Distributed] 
      MultinomialDistribution[balls, ConstantArray[1/bins, bins]]]; //
   AbsoluteTiming // First


r1 == r2





So, about 50,000X faster than built-in on a smaller example than the OP (I got tired of waiting for the built-in to finish on that example, but it was over half an hour running before I killed it, on a fairly beefy laptop.)

The OP case:

pMax[100, 10, 20]; // AbsoluteTiming//First


This can handle pretty big cases.

You may need to increase $RecursionLimit for larger cases.

The last optional argument allows use of inexact calculations. In use, its value is generally of the order of the # of balls, but experiment - it varies with case characteristics. You'll get the indication of "No significant digits are available to display" warning if too small. In any case, a nice performance boost can be had then if exact results are not required:

Block[{$RecursionLimit = 50000}, pMax[5000, 300, 30, 5000] // N[#] &] // AbsoluteTiming
Block[{$RecursionLimit = 50000}, pMax[5000, 300, 30] // N[#] &] // AbsoluteTiming



Caveat: This is an old quick-and-dirty, it has no real error checking, so feed it nonsense values at your own peril.


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.