Take the following:
p= N@{459/10703, 1/973, 95/10703, 635/10703, 179/10703, 565/10703,
54/973, 474/10703, 794/10703, 548/10703, 52/1529, 1/10703,
61/973, 86/1529, 775/10703, 162/10703, 160/10703, 870/10703,
157/10703, 816/10703, 691/10703, 471/10703, 192/10703,
10/973, 307/10703};
q={5, 6, 8, 20, 5, 14, 15, 12, 14, 5, 16, 14, 9, 8, 14, 9, 17, 18,
6, 6, 7, 11, 14, 17, 15};
n=200;
CDF[MultinomialDistribution[n, p], q] // AbsoluteTiming
I gave up waiting and aborted.
Any ideas for a faster method in native Mathematica for multinomial distribution CDF?
Edit: I'll be adding a bounty as soon as it's available to stimulate ideas/answers.
Here's a more trivial test case to use for those interested with results/timings from loungebook. Even with this much simpler case, better than five orders of magnitude performance gain.
probs = {.15, .1, .05, .2, .35, .12, .03};
counts = {18, 16, 13, 20, 21, 12, 10};
n = 49;
ClearSystemCache[];
fast = multinomFastCDF2[n, probs, counts] // AbsoluteTiming
ClearSystemCache[];
mma = CDF[MultinomialDistribution[n, probs], counts] // AbsoluteTiming
Last@fast == Last@mma
First@mma/First@fast
{0.00393243, 0.897792}
{2526.37, 0.897792}
True
642446.
Answer
Came up with this, using N[{...},14]
on the setting of p
in the OP to get sufficient result precision, it finishes in a few tenths on the loungebook:
multinomFastCDF2[n_Integer, p_, q_] := If[n < 0 || n > Tr[q], 0,
FullSimplify[E n! Last@Fold[Take[ListConvolve[##, {1, -1}, 0], UpTo[n + 1]] &,
MapThread[Divide[Gamma[#2 + 1, #1], Gamma[#2 + 1]]
Normalize[Divide[E^-#1 #1^Range[0, #2], Range[0, #2]!], Total] &,
{p,q}]]]];
Arguments are the number of trials, category probabilities, and CDF vector to evaluate.
If anyone can come up with better, I'm all eyes...
Comments
Post a Comment