# Can this multinomial probability query performance be improved?

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?

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


0.00363067

181.405

49964.7

True

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


0.0102418

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


{0.918627,0.727751}

{15.452,0.727751}

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