How to get probabilities for multinomial & hypergeometric distribution ranges more quickly?

Question

I often need to calculate the probability of a result from a multinomial distribution falling within some range, e.g., the maximum of the categories having an upper bound, or a lower bound for the minimum, or both.

For example, posed in a classic "balls and bins" problem style, one might ask "If 50 balls are thrown at 4 bins, with the probabilities of landing in the respective bins of {1/2, 3/8, 1/16, 1/16}, what's the probability that no bin has less than 5 balls or more than 20 balls?"

In Mathematica, one way of evaluating this would be :

Probability[Min[a, b, c, d] >= 5 && Max[a, b, c, d] <= 20, 
         {a, b, c, d} \[Distributed] MultinomialDistribution[50, {1/2, 3/8, 1/16, 1/16}]]

But, as anyone who's done such things in Mathematica knows, this will get untenable pretty quickly - the above takes 5 minutes on my loungebook. Adding categories and/or more trials explodes the time exponentially.

Sometimes the problem is even uglier:

Probability[Min[a, b, c, d] >= 11 && Max[a, b, c, d] <= 20, 
         {a, b, c, d} \[Distributed] MultinomialDistribution[40, {1/2, 3/8, 1/16, 1/16}]]

Mathematica will merrily churn along, finally returning 0, something that can be seen readily by simple inspection: you can't have at least 11 in each category with only 40 trials.

Is there a more efficient / faster way of getting such results? And the same for the multivariate hypergeometric?

Answer 1

Yes, there is!

One way to do such things more quickly:

multinomialRangedIP[n_, ps_, min_, max_] := 
 Total[PDF[MultinomialDistribution[n, ps], Permutations@#] & /@ 
   IntegerPartitions[n, {Length@ps}, Range[min, max]], 2]

Using the OP example,

multinomialRangedIP[50, N@{1/2, 3/8, 1/16, 1/16}, 5, 20]

is over 200X faster. That advantage grows as the number of trials/categories increases.

However, with large enough cases, even this becomes impractical (although it will always return a 0 result quickly for impossible configurations).

Another technique I use for larger cases:

multinomialRanged[n_, ps_, min_, max_] :=
  Module[{ups = MapIndexed[(#1/Total[ps[[First@#2 ;;]]]) &, ps], sa, lp = Length@ps, pr},

   pr[nx_, sk_, sk1_, pk_, minx_: 0, maxx_: n] := 
    Piecewise[{{0, sk < sk || sk - sk1 > maxx || sk - sk1 < minx},
      {Binomial[nx - sk1, sk - sk1]*pk^(sk - sk1)*(1 - pk)^(nx - sk),True}}];

   sa = Append[
     Table[Array[pr[n, #2 - 1, #1 - 1, ups[[k]], min, max] &, {n + 1, n + 1}], {k, 1, lp - 1}], 
     Reverse@Array[{Boole[# - 1 >= min && # - 1 <= max]} &, n + 1]];
   sa[[1]] = sa[[1, 1]];

   First@(Dot @@ sa)];

Comparing

multinomialRangedIP[60, N@{1/4, 1/4, 3/8, 1/16, 1/32, 1/32}, 5, 20] 
multinomialRanged[60, N@{1/4, 1/4, 3/8, 1/16, 1/32, 1/32}, 5, 20] 

Shows an over 100X speed advantage for multinomailRanged over multinomailRangedIP, again with the advantage increasing with increasing problem size. I have no idea how long Mathematica would take using native functionality - I ran out of patience...

There's certainly some optimization left in the latter code, e.g. the structure of the array has high density of zeroes, so time could probably be shaved with a custom Dot, but once it was fast enough for my needs, I left it as is.

The same technique in the latter can be applied to ranged probabilities for the multivariate hypergeometric:

hypergeometricRanged[n_, ps_, min_, max_] := 
  Module[{vkseq, t, prh, sa},

   vkseq = Accumulate@Prepend[ps, 0];
   t = Total@ps;

   sa = PadLeft[
     Table[If[min <= Subtract[y, x] <= max, 
       With[{sk = Subtract[y, 1], sk1 = Subtract[x, 1], vk = vkseq[[z + 1]], vk1 = vkseq[[z]]}, 
        If[(prh = Binomial[t - vk1, n - sk1]) == 0, 0, 
         Binomial[vk - vk1, sk - sk1] Binomial[t - vk, n - sk]/prh]], 
       0], {z, Length@ps - 1}, {x, n + 1}, {y, x, n + 1}]];

   sa[[1]] = sa[[1, 1]];

   sa = Append[sa, 
     Transpose[{Array[Boole[min <= n + 1 - # <= max] &, n + 1]}]];

   First@(Dot @@ sa)];

On a fairly trivial example:

hypergeometricRanged[30, {20, 20, 20, 10, 5}, 2, 12]

Probability[Min[a, b, c, d, e] >= 2 && Max[a, b, c, d, e] <= 12,
            {a, b, c, d, e} \[Distributed] 
                MultivariateHypergeometricDistribution[30, {20, 20, 20, 10, 5}]]

this was ~8500X faster, and as before, advantage grows dramatically with problem size.

N.B.: I use these in routines that call with well-formed cases, so no error checking is done - call them with really whacky things (like negative values) at your own peril (or add E/C code).