Probability of reaching a threshold after several trials

I was wondering if there wasn't some better way to get to the bottom of this question I had regarding the chances of reaching a certain threshold of guessing numbers correctly over a series of trials. This isn't my specialty, so would appreciate answers that are like you were talking to a 5 year old or at least google-able :)

So say I have a game where I have to pick 10 out of 80 numbers and then 20 are drawn. For each one of my Picks that is chosen as one of the 20 I get a black marble.

I need to know the odds that if I can get AT LEAST 20 black marbles in 5 games, picking 10 numbers each game. Each game is discrete in that the picked numbers are replaced and the 20 numbers that are drawn are totally random

I can tell what the average amount of black marbles I will win per game (2.5) and therefore the average amount of black marbles I'll win in 5 games (12.5). I can do this because for a single game, hyper geometric distribution tells me the odds of each outcome of catches: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. I can then multiply each of these by their probabilities and sum them to get the average amount I will win in that game
e.g. $(HYPGEOMDIST(0, 10, 20, 80) * 0) + (HYPGEOMDIST(1, 10, 20, 80) * 1)$, etc

Here's how I would do this by hand if I had to (this is the part I'd like to simplify / correct my assumptions on!!):

Find each permutation of black marbles (0-10) I could win for each of the 5 games without regard to order.

Filter each permutation so that if the set of 5 games does not sum to 20 it is discarded
For the remaining sets find and multiply the HYPGEOMDIST for each of the 5 games in each set, then sum every set

Please excuse me if this has been asked somewhere else, I am unfamiliar with the terminology in order to search well enough on this site.

You have correctly identified that your problem is one of convolution of hypergeometrics.

Your two simplest choices appear to be either:

1. Direct numerical convolution of the probability functions. This is not difficult and given the probability function is on a grid, it can even be done in Excel fairly easily.

2. You could perhaps try an approximation to the hypergeometric, such as perhaps a normal approximation; if we're not too far in the tail for the sum it should be reasonable.

So let's try them:

1. Direct convolution:

(Calculation details omitted - I did it in R using the convolve function, though it could be done with fft and its inverse or filter as well, but they're all a bit fiddly. In this case it's probably easier to just do it directly as a calculation on the vector of probabilities.)

This is what the probability function of the convolution looks like:

The probability of getting $\geq 20$ on that p.f. is 0.009385. This is 'exact' to the numerical accuracy of the convolution (which should be to more than the number of figures I quoted, so the quoted figures should be accurate).

Here's the same graph with the p.f. of the hypergeometric and the convolution of 5 of them computed in Excel, in this case by repeated use of the SUMPRODUCT function. The answers are of course the same:

2. Normal approximation. Since it turns out (as we see from 1.) that we're quite far out in the tail, this won't be very accurate. But let's go ahead and try it:

From Wikipedia

For a single set of 20 draws ($N=80, K=20, n=10$, and $k$ is the number of 'hits' (black marbles you get).

Your mean for $k$ is: $n {K\over N} = 10 \times 20/80 = 2.5$ as you stated.

Your variance for $k$ is: $n{K\over N}{(N-K)\over N}{N-n\over N-1}$ $= 2.5 \cdot \frac{3}{4} \cdot \frac{70}{79} = 1.6614$ (s.d. 1.289)

Hence $S=Y_1+Y_2+\ldots +Y_5 \quad \dot{\sim}\,\, N(12.5,2.8822^2)$

With continuity correction:

$P(S\geq 19.5) \approx P(Z\geq (19.5-12.5)/2.8822) = P(Z\geq 2.4287) = 0.0076$

The z-value cutoff itself here is enough to hint to us that the approximation may not be very accurate.

(While not terrible accurate, the uncorrected version is significantly less accurate in this case, having about half the required probability.)

Given that the numerical convolution isn't all that much harder than this inaccurate approximation, in this case, it's worth the additional effort.

---

An alternative approximation that is often used is the binomial approximation to the hypergeometric; the usual one would take $n=10$ and $p=0.25$. In this case, that's quite convenient because the convolution of them will also be binomial (with $n=50$ and $p=0.25$); however, it tends to be inaccurate when you sample more than about 5-10% of the population, and in this situation it's a lot more than that, so we would anticipate that it won't be very good. Further, the usual form overestimates the variance (since it ignores the 'without replacement' aspect). In this case it yields a probability about 50% too high.

However, as Brunk et al.$^{(1)}$ point out, since several of the hypergeometric parameters can be interchanged without changing the probabilities, there are actually four quite natural binomial approximations to the hypergeometric, and they discuss another reference that says the one with the smallest "n" for the binomial is generally the best approximation. Brunk et al. give a detailed analysis which broadly seems to back this simple advice up.

In their notation, the total population is $N$, there are $M$ items in the population with the characteristic of interest and there are $n$ in the sample (i.e. my $K$ above is their $M$ and my $k$ is their $x$), then:

$$H(x;N,n,M) = \binom{M}{x}\binom{N-M}{n-x}/\binom{N}{n}\,,$$

where $0 \leq x \leq \min(M,n)$ and the usual binomial approximating probability is

$$B(x;n,M/N) = \binom{n}{x}(M/N)^r\, (1-M/N)^{n-r}\,,$$

for $0 \leq x \leq n$.

Since in the hypergeometric, $n$ and $M$ can be interchanged, as can the roles items-of-interest and items-not-of-interest,

$$H(x;N,n,M)=H(x;N,M,n)=H(n-x;N,N-M,n)=H(M-x;N,N-n,M)$$

leading to four binomial approximations to those probabilities:

$B(x;n,M/N)$
$B(x;M,n/N)$
$B(n-x;N-M,n/N)$
$B(M-x;N-n,M/N)$

In our situation, $N=80$, $M=20$ and $n=10$, so the four "sample sizes" are minimized by taking the first approximation, the one we already considered; this suggests we likely can't do better with the usual forms of binomial approximation.

--

An alternative possibility for a binomial approximation, since the probabilities are essentially zero near the maximum is to forget sticking to the given $n$ and instead choose $n$ and $p$ so as to approximately match the first two moments:

$np = 2.5$

$np(1-p)=1.6614$

Dividing the second equation by the first, $1-p = 0.66456$ or $p = 0.33544$, implying $n\approx 7.453$. If we seek to approximate the individual components of the convolution, the possibilities are to take $n = 7$ or $n=8$, with implied variances $1.607$ and $1.719$ (and of those we might anticipate $n=8$ to do slightly better since it gives more of a 'right tail'), but we could also consider rounding only to the nearest $0.2$ ($n=7.4$) so that we have integer $n$ once we have added the five components, yielding these approximations after the convolution:

(i) $n=35, p = 2.5/7$

(ii) $n=37, p = 2.5/7.4$

(iii) $n=40, p = 2.5/8$

> 1-pbinom(19,35,2.5/7)
[1] 0.007825608
> 1-pbinom(19,37,2.5/7.4)
[1] 0.008837474
> 1-pbinom(19,40,2.5/8)
[1] 0.01022883


With the smallest percentage error on the middle one, though the $n=8$ one is also "close". Nevertheless, the accuracy is not particularly good, though much better than our other approximations.

(1): Brunk, H. D., J. E. Holstein and F. Williams (1968),
A Comparison of Binomial Approximations to the Hypergeometric Distribution
The American Statistician, Vol. 22, No. 1 (Feb.), pp. 24-26