Basically, my question is the same as this one, except I need more than the $k = 0$ special case:
Given a sum of independent random variables each following a hypergeometric distribution, is there any efficient way to compute the PMF for that mixture?
At the very least, are there any tricks that might make a numerical evaluation less painful than a straightforward convolution (for cases where the number of variables and/or population size is very high)?
Sampling in the ballpark of half the population doesn't usually lead to problems with the normal approximation. Here's an example for the distribution of the number of white balls drawn from a population of 300 white balls and 700 black balls, sampling 500 balls without replacement, along with a normal distribution with the same mean and variance as the hypergeometric.
x=dhyper(120:180,300,700,500)
plot(120:180,x,type="h")
This suggests that as long as the number of each kind of ball are not too large or small and the total population size is reasonably large, just using normal approximations (possibly with continuity correction, depending on circumstances) may be quite feasible.
In cases where the normal approximation on the individual hypergeometric components isn't reasonable it may still be that a normal approximation to the sum may be adequate if there are enough terms in the sum.
If the proportion of white balls (/"successes") is very large or vary small compared to the population, a binomial or even a Poisson approximation (to the black balls if the proportion of white is large) may be adequate, which suggests either a Poisson-binomial or Poisson approximation to the sum may be reasonable in some situations. In other cases, a moment-matched (possibly shifted-) binomial may be adequate.
If none of those are adequate you may have to fall back on convolution. Assuming all terms have different parameters, so that it doesn't admit some shortcuts, $f=f_1*f_2*...*f_n$ can be approached in a number of different ways. Let $\hat{f}_i=\mathcal{F}(f_i)$ be the Fourier transform of the $i\,$th term; then the Fourier transform of the convolution is $\hat{f}=\prod_i \hat{f}_i$. These products might be generated one at a time (that is, looping through and taking the product of each with the next), or - especially if you have multiple cores that can run in parallel, they might be generated pairwise, then pairs combined and so on in turn, before inverting back to the resulting density, $f=\mathcal{F}^{-1}(\hat{f})$.
I don't think the use of the Fourier transform should be slow unless there are so many components that suggestion 2. should probably work.
I commented above asking if someone knows how to do a convolution in closed form. Based on @Glen_b's comment about doing it computationally, I figured I'd add an example in Python in case it is helpful.
from scipy.stats import hypergeom
import numpy as np
import matplotlib.pyplot as plt
def distr_of_indep_xsum(x1, x2, p1, p2, trim=False):
# calculate probablities of sums of x1, x2 happening independently
mx = x1.max() + x2.max()
# consider values from 0 to mx so can just use index
xvals = np.arange(mx + 1)
pvals = np.zeros(mx + 1)
for xx1, pp1 in zip(x1, p1):
for xx2, pp2, in zip(x2, p2):
xsum = xx1 + xx2
# independent probabilities are multiplied
pvals[xsum] += (pp1*pp2)
if trim:
max_nonzero = np.where(pvals > 0)[0].max()
xvals = xvals[:(max_nonzero+1)]
pvals = pvals[:(max_nonzero+1)]
return xvals, pvals / pvals.sum()
def hgprobs(ndraws, nitems, nitems_same):
# note: the notation in the function differs from that in wikipedia, so I am making [![enter image description here][1]][1]it clearer
x = np.array(list(range(0, ndraws + 1))).astype(int)
p = np.array([hypergeom.pmf(xx, nitems, nitems_same, ndraws) for xx in x])
return x, p
def hypergeometric_convolution():
# two hypergeometric variables
x1, p1 = hgprobs(ndraws=3, nitems=10, nitems_same=5)
x2, p2 = hgprobs(ndraws=12, nitems=20, nitems_same=10)
xsum, psum = distr_of_indep_xsum(x1=x1, x2=x2, p1=p1, p2=p2)
plt.bar(xsum, psum)
plt.show()
hypergeometric_convolution()
