# Inconsistent simple simulation outcomes?

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 ?

• 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. Commented May 21, 2019 at 12:37
• 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 Commented May 21, 2019 at 13:51
• 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. Commented May 22, 2019 at 10:50

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.

• 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. Commented May 22, 2019 at 18:57
• BTW I regularly use 11.3, but same issue in 10.4. Commented May 22, 2019 at 19:08
• @HagenScherb I updated the bugs header in the question to reflect the version info. Commented May 22, 2019 at 19:28