# Can this probability function speed be improved?

I have a process that as part of it requires a certain probability query, specifically, given a multinomial distribution with f equiprobable categories and r trials, some specific category has a given value v and that value is the (not necessarily unique) maximum across all categories.

This can obviously be done directly in Mathematica with

Probability[a[1] == v && v == Max[Array[a, f]],
Array[a, f] \[Distributed] MultinomialDistribution[r, ConstantArray[1/f, f]]]


but it is quite sluggish, for example with {r,f,v}={25,12,6} it takes some... well, I don't know, I aborted it after several minutes running.

I've come up with

f[r_, f_, v_] :=
Block[{t},
Coefficient[
Expand[Sum[t^j/j!, {j, 0, v}]^(f - 1)] (f - 1)^(-(r - v))*
(r - v)!*Binomial[r, v]/f^v (1 - 1/f)^(r - v), t^(r - v)]];


which is fairly snappy.

Can this query be done faster? I'd like to handle cases of r up to ~200, f up to ~25 and v generally ~r/3.

Edit/Update:

Using a technique from an old answer of mine, the new function

f2[r_, f_, v_, pf_ : Infinity] := Module[{pMax},
pMax[q_, b_, max_, p_] := 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;
q! h[q, b, max]/b^q];
pMax[r - v, f - 1, v, pf] Binomial[r, v]/f^v (1 - 1/f)^(r - v)];


is notably faster on larger cases and comparable on small ones. For example, with {r,f,v}={25,12,15} this takes ~0.36ms, vs ~2.5ms for my original function. By way of comparison, the direct calculation in Mathematica took ~2800 seconds.

If no better result arrives, I'll self-answer with this function.

• Here's is a (extremely?) minor speed improvement: Change the Coefficient[] term to If[r > v f, 0, Coefficient[Expand[Sum[t^j/j!, {j, 0, v}]^(f - 1)] as the probability is zero if $r>v f$. That way the Expand doesn't have to happen every time. The only other suggestion I have is to look at the coefficients in the expansion to see if something matches in oeis.org. – JimB Aug 30 '20 at 16:16
• @JimB - already done in auxiliary definitions I left out of OP for brevity. I'd already poked at OEIS with some tantalizing but eventually faulty hits, but will revisit. Thanks for the suggestions! – ciao Aug 30 '20 at 17:43
• Here's one more desperation level suggestion: Some of the extreme results might be easier to simplify. For example, if $2v>n$, then there can be only one category with a frequency of $v$. The function that matches your function under that condition is f2[r_, f_, v_] := Binomial[r, v] (f - 1)^(r - v)/f^v. A little more complicated would be if $3v>n$ which would mean (I think) that there could only be at most 2 categories with the same value of $v$. – JimB Aug 30 '20 at 18:04
• You're probably way ahead of me on this but maybe the following which uses a recurrence relation (rather than getting the associated coefficient) on a slightly different problem (not conditioning on a specific category having the maximum) might help: pnas.org/content/104/27/11189. – JimB Aug 30 '20 at 22:43
• @JimB ;-) - I am testing right now with that very technique from an old answer of mine that I remembered about an hour ago: mathematica.stackexchange.com/a/186535/11467 It's looking promising for large cases, silly I didn't recognize I'd seen this kind of problem before. – ciao Aug 30 '20 at 22:53

This is just an extended comment: Because the question is about speeding things up, I think you'll need to consider merging different functions depending on what values of $$v$$ are of most interest. Here is an example of the timing:

f[r_, f_, v_] := Block[{t}, Coefficient[Expand[Sum[t^j/j!, {j, 0, v}]^(f - 1)] (f - 1)^(-(r - v))*
(r - v)!*Binomial[r, v]/f^v (1 - 1/f)^(r - v), t^(r - v)]];

f2[r_, f_, v_, pf_ : Infinity] := Module[{pMax},
pMax[q_, b_, max_, p_] := 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;
q! h[q, b, max]/b^q];
pMax[r - v, f - 1, v, pf] Binomial[r, v]/f^v (1 - 1/f)^(r - v)];

(* This function only useful when 2v > r *)
f3[r_, f_, v_] := If[2 v > r, Binomial[r, v] (f - 1)^(r - v)/f^r]

r = 200;
g = 25;  (* Really this is f but I get confused between the constant f and the function f *)
data = ConstantArray[{0, 0, 0, 0}, r];
Do[data[[v]] = {v, If[2 v > r, AbsoluteTiming[Binomial[r, v] (g - 1)^(r - v)/g^v][[1]], 10],
AbsoluteTiming[f[r, g, v]][[1]],
AbsoluteTiming[f2[r, g, v]][[1]]}, {v, 1, r - 1}]

ListLogPlot[{data[[All, {1, 3}]], data[[All, {1, 4}]], data[[All, {1, 2}]]},
PlotRange -> {Automatic, {0, 9}}, Frame -> True,
FrameLabel -> (Style[#, Bold, 18] &) /@ {"v", "Seconds"},
PlotLegends -> {"f", "f2", "f3"}]


• Don't know how I missed this (excellent) response. Mea culpa, +1 and accept. – ciao Apr 9 at 22:35