10
$\begingroup$

2016-04-09: This bug is indeed fixed (thank you, Stefan R.!) in 10.4. It may have been fixed in 10.3 or 10.3.1, but I haven't tried these versions.

dist=MultivariateHypergeometricDistribution[5,ConstantArray[2,10]];
prob=Probability[a==1,Distributed[{a,b,c,d,e,f,g,h,i,j},dist]]

We stuff an urn with coloured balls, ten colours, two balls of each colour, for a total of twenty. We draw five balls without replacing them; the list {a,…,j} then (loosely speaking) represents what we just have drawn. The probability of picking exactly one ball of colour "a" P(a==1) then comes out to 15/38, easily verified and computed in a snap by Mathematica.

What about P(a==1&&b==1), the probability of getting exactly one "a" and one "b" ball each? You can probably model this in your head, do the calculation with pencil and paper – and finish long before Mathematica 10.2, which hasn't finished yet in the time it took me to type all this.

Or just checking for P(a==1&&b<0)? Neat, the answer is instantaneous and correct (p==0). And since the number of balls of colour "b" must obviously be at least 0, let's try this, too:

dist=MultivariateHypergeometricDistribution[5,{2,2,2,2,2,2,2,2,2,2}];
prob=Probability[a==1&&b>=0,Distributed[{a,b,c,d,e,f,g,h,i,j},dist]]

P(a==1||b>=0)? No luck. P(a==b)? Perhaps next week. P(a!=b)? Nope.

Apparently, the cutoff is at five dimensions, that is, when the length of the list that is passed to MultivariateHypergeometricDistribution[n,list] becomes greater than five. Up to five, everything works as expected.

My question is, I fear, obvious and not very specific: What am I doing wrong? Or must I keep recalling my high-school maths and calculate my probabilities on an abacus?

$\endgroup$
  • 2
    $\begingroup$ Not doing anything "wrong", this is an area of weakness in an otherwise quite nice set of facilities: sometimes you just have to use manual intervention to transform the problem... $\endgroup$ – ciao Sep 14 '15 at 23:08
  • $\begingroup$ @ciao Thanks, that is a relief. Maybe in MMA 11 … $\endgroup$ – Felix Kasza Sep 14 '15 at 23:13
  • $\begingroup$ See my CW answer.... might be of use. $\endgroup$ – ciao Sep 14 '15 at 23:16
  • 2
    $\begingroup$ Thanks for this question! We've looked into it and found that we can substantially speed up this computation. The next version of Mathematica will be able to find this probability in less than a second. $\endgroup$ – Stefan R Sep 16 '15 at 15:41
  • $\begingroup$ @StefanR Thank you very, very much! $\endgroup$ – Felix Kasza Sep 17 '15 at 12:24
9
$\begingroup$

This is one of my few gripes re: an otherwise quite nice probability functionality. Sometimes, performance is inexplicably poor.

Usually, trivial manual transformation/intervention can get the results desired speedily, e.g., your example cases:

ClearAll["Global`*"]

dist = MultivariateHypergeometricDistribution[5, ConstantArray[2, 10]];

PDF[dist, 
  Cases[Join @@ 
    Permutations /@ IntegerPartitions[5, {10}, Range[0, 5]], 
        {a_, b_, c_, d_, e_, f_, g_, h_, i_, j_} /; a == 1]] // Tr

PDF[dist, 
  Cases[Join @@ 
    Permutations /@ IntegerPartitions[5, {10}, Range[0, 5]], 
    {a_, b_, c_, d_, e_, f_, g_, h_, i_, j_} /; a == 1 && b >= 0]] // Tr

PDF[dist, 
  Cases[Join @@ 
    Permutations /@ IntegerPartitions[5, {10}, Range[0, 5]], 
       {a_, b_, c_, d_, e_, f_, g_, h_, i_, j_} /; a == b]] // Tr

PDF[dist, 
  Cases[Join @@ 
    Permutations /@ IntegerPartitions[5, {10}, Range[0, 5]], 
    {a_, b_, c_, d_, e_, f_, g_, h_, i_, j_} /; a != b]] // Tr

All return quickly.

When only certain categories are involved, further transformation will speed things hugely:

dist = MultivariateHypergeometricDistribution[5, Append[ConstantArray[2, 2], 16]];

PDF[dist, 
  Cases[Join @@ 
    Permutations /@ IntegerPartitions[5, {3}, Range[0, 5]], {a_, b_, c_} /; a != b]] // Tr

Obviously, for cases with equal categories, further speed can be had by dispensing with permutations and applying appropriate factor to each result for number of permutations.

Lastly, for ranged queries, see How to get probabilities for multinomial & hypergeometric distribution ranges more quickly?, a little sorcery I cooked up...

$\endgroup$
  • $\begingroup$ Thank you! I had already found your article on ranges with multivariate hypergeometric distributions and was prepared to solve my problem "by hand", a task which you were so kind as to anticipate, but ranges were not quite what I needed. Also, I wanted to be sure that I was dealing with an idiosyncrasy in MMA and not just with my own daftness … $\endgroup$ – Felix Kasza Sep 15 '15 at 6:04
7
$\begingroup$

Although @ciao's solution gets you to the answer, I would like to offer perhaps another angle at it.

Given a tuple $\{X_1, X_2, \ldots, X_n\}$ that follows a multivariate hypergeometric distribution with parameters $N$, $\{M_1, \ldots, M_n\}$, the tuple $\{X_1, X_2, \sum_{k=3}^n X_k \}$ also follows a multivariate hypergeometric with parameters $N$ and $\{M_1, M_2, \sum_{k=3}^n M_k\}$.

Hence the problem can be reformulated with fewer variables, which leads to a faster solution:

AbsoluteTiming[
 Probability[a == 1 && b >= 0, 
  Distributed[{a, b, c}, 
   MultivariateHypergeometricDistribution[
    5, {2, 2, Plus[2, 2, 2, 2, 2, 2, 2, 2]}]]]]

(* Out[29]= {0.0504801, 15/38} *)

AbsoluteTiming[
 Probability[a == b, 
  Distributed[{a, b, c}, 
   MultivariateHypergeometricDistribution[
    5, {2, 2, Plus[2, 2, 2, 2, 2, 2, 2, 2]}]]]]

(* Out[30]= {0.0075801, 138/323} *)

AbsoluteTiming[
 Probability[a > b, 
  Distributed[{a, b, c}, 
   MultivariateHypergeometricDistribution[
    5, {2, 2, Plus[2, 2, 2, 2, 2, 2, 2, 2]}]]]]

(* Out[32]= {0.00749303, 185/646} *)
$\endgroup$
  • $\begingroup$ This is the second example in my answer.... $\endgroup$ – ciao Sep 16 '15 at 0:54
  • $\begingroup$ The transformation, which @ciao also hinted at, is elegant and perfectly correct. Regrettably, it does not cover my original case; only looking at two out of ten categories was just a simple case that still demonstrated the sub-optimal behaviour of Mathematica. For my, admittedly rather odd, case, I think ciao's original answer is currently the best option I have. But I am grateful for making explicit and reminding me of the possibility of collapsing unneeded categories! $\endgroup$ – Felix Kasza Sep 16 '15 at 9:59

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.