5
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Aug 30 '20 at 16:16
  • $\begingroup$ @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! $\endgroup$
    – ciao
    Aug 30 '20 at 17:43
  • 1
    $\begingroup$ 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$. $\endgroup$
    – JimB
    Aug 30 '20 at 18:04
  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Aug 30 '20 at 22:43
  • $\begingroup$ @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. $\endgroup$
    – ciao
    Aug 30 '20 at 22:53
5
$\begingroup$

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

Plot of timing for each of the 3 functions

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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