4
$\begingroup$

Bug introduced in 10.4 or earlier and persisting through 12.0.0 or later


Two simple kinds of simulating the one-sided exact Fisher-test of obtaining at least 86 successes in 129 draws in a population of size 299 containing 167 successes

1st Version is rather slow and yields too large a probability

t=Table[RandomVariate[HypergeometricDistribution[129,167,299]],
1000000];//AbsoluteTiming
p1=Count[t,u_/;u>=86]/1000000.
{136.942,Null}
0.001054

2ndVersion is fast and its result is pretty similar to the theoretical probability, see below

sf=Join[s=Table[1,167],f=Table[0,132]];
t=Table[RandomSample[sf,129]//Total,10000000];//AbsoluteTiming
p2=Count[t,u_/;u>=86]/10000000.
{18.5029,Null}
0.000744

Theoretical result

p3=1-CDF[HypergeometricDistribution[129,167,299],86-1]//N
0.00074213

Comparisons

p1/p3
p2/p3
1.42024
1.00252

What is the defect in the 1st version ?

![Behavior of 10^6 simulations of HypergeometricDistribution[129,167,299]. Thick gray line:exact HypergeometricDistribution[129,167,299]. Dashed lines:two bulk-generated realizations.
Red line:singly-generated realization (using RandomVariate[...] but not RandomSample[...] !!) . ]
1

$\endgroup$
  • $\begingroup$ RandomVariate[HypergeometricDistribution[129, 167, 299], 1000000] is considerably faster than both your methods. Seems more accurate, too, but I'm not sure why p1 is so far off. -- Also p1 = Mean@N@UnitStep[t - 86] is much faster, but neither way takes much time. $\endgroup$ – Michael E2 May 21 at 12:37
  • 1
    $\begingroup$ TY. Your code is very fast and yields a reasonably precise solution, which now, however, is systematically somewhat too small. I submitted my question to WM-support and will let you know their answer - if any. Seems there is some kind of error propagation or systematic error increase/decrease depending on the kind we code. HS. t=RandomVariate[HypergeometricDistribution[129,167,299],1000 000 000]; (p1=Mean@N@UnitStep[t-86]//N) p1/p3 0.000742018 $\endgroup$ – Hagen Scherb May 21 at 13:51
  • $\begingroup$ I added the bugs tag. There hasn't been much community interest/discussion, but it can be removed if the issue turns out to be something else. $\endgroup$ – Michael E2 May 22 at 10:50
4
$\begingroup$

It seems a bug, perhaps. It makes a difference how the random variates are generated, a disturbing difference. The variates can be generated singly or in bulk:

Table[RandomVariate[HypergeometricDistribution[129, 167, 299]], k] (* singly *)
RandomVariate[HypergeometricDistribution[129, 167, 299], k]        (* bulk *)

For k up to the number of draws 129, these will generate the same variates. But for k > 129 they are different and have different properties. It seems that generating variates singly is faulty.

I'll use SeedRandom[] and BlockRandom[] to make certain results reproducible. First, some basic parameters of the distribution:

mean = N@Mean@HypergeometricDistribution[129, 167, 299]
sd = N@StandardDeviation@HypergeometricDistribution[129, 167, 299]
(*
  72.0502
  4.25976
*)

Here is a comparator of the two methods of generating variates. It also sows the the generated samples, which we will sometimes use for further analysis. We can see that the two methods generate the same samples up to k == 129 and different samples for k >= 130:

cmp[k_] := 
 BlockRandom[
   Sow@Table[RandomVariate[HypergeometricDistribution[129, 167, 299]],k]] ==
  BlockRandom[
   Sow@RandomVariate[HypergeometricDistribution[129, 167, 299], k]]

And @@ Table[SeedRandom[k]; cmp[k], {k, 129}]
(*  True  *)

Or @@ Table[SeedRandom[k]; cmp[k], {k, 130, 200}]
(*  False  *)

Here is a utility for checking a sample for a given size k. It compares the mean of the same with the expected mean:

check[k_] := Module[{sameq, t1, t2},
   {sameq, {{t1, t2}}} = Reap@cmp[k];
   <|"same" -> sameq,
    "Z[singly]" -> (Mean@N@t1 - mean) Sqrt[k]/sd,
    "Z[bulk]" -> (Mean@N@t2 - mean) Sqrt[k]/sd|>
   ];

For small samples, it's hard to say there's a problem:

SeedRandom[1];
check[129]
(*  <|"same" -> True, "Z[singly]" -> 1.74712, "Z[bulk]" -> 1.74712|>  *)

SeedRandom[1];
check[130]
(*  <|"same" -> False, "Z[singly]" -> 1.657, "Z[bulk]" -> -0.607834|>  *)

For larger k, we begin to get unlikely deviations from the mean:

SeedRandom[1];
check[2000]
(*  <|"same" -> False, "Z[singly]" -> 2.5284, "Z[bulk]" -> 0.706897|>  *)

SeedRandom[1];
check[10^6]  (* OP's sample size. Beware: it takes a couple of minutes *)
(*  <|"same" -> False, "Z[singly]" -> 53.5724, "Z[bulk]" -> -1.31421|>  *)

I think that suggests a bug in the way iterates are generated when generated in sample sizes under 130. The OP notes in a comment that for larger sizes, the mean is "systematically somewhat too small." That's not the case for k = 2000, but for 10^6 it is somewhat unusually small. It's also rather small for k = 10^5 and the OP checked it for k = 10^9, but I don't propose to check thoroughly.

Update: It seems that the switch between the two different algorithms for RandomVariate[HypergeometricDistribution[a, b, c], k] occurs when k is greater than any of a, b or c - b.

$\endgroup$
  • 1
    $\begingroup$ I received the following message from WTS: Hello, Thank you for contacting Wolfram Technical Support. I understand that the two calculations below are giving inconsistent results with the first calculation giving the wrong result …. . I have reported this issue to our developers. Thank you for bringing this to our attention. Regards, Wolfram Technical Support, Wolfram Research Inc. $\endgroup$ – Hagen Scherb May 22 at 18:57
  • $\begingroup$ BTW I regularly use 11.3, but same issue in 10.4. $\endgroup$ – Hagen Scherb May 22 at 19:08
  • $\begingroup$ @HagenScherb I updated the bugs header in the question to reflect the version info. $\endgroup$ – Michael E2 May 22 at 19:28

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.