3
$\begingroup$

I have performed a calculation that a sample of records will contain at least 1 defective record assuming replacement. There are 20,000,000 records in the population, 7,000 are defective, and I will be taking 200,000 samples.

Assuming replacement, the probability that at least 1 of the 200,000 samples will contain a defective record is: 1 - (7,000/20,000,000)^200,000 = 1 - 3.9 x 10^(-31)

If I remove the assumption of replacing the records after each sample, I believe the calculation becomes: The product of [1- 7,000/(20,000,000 -n] with n=0 to n=(200,000 - 1)

My problem is, I don't now how to perform the product of the sequence calculation. In excel, I can create the 200,000 individual probabilitys and multiply them together, but I'm afraid rounding will be an issue (the probability was 1 - 2.7 x 10^(-31))

Can any suggest how I can perform this product of a sequence calculation more accurately?

$\endgroup$

1 Answer 1

6
$\begingroup$

The log Gamma function is the tool to use. Some caution is needed concerning the precision of the result, but even for large problems like this excellent precision can be attained with Excel. (This caution applies to all statistical packages that use double-precision floating point arithmetic.)


When sampling $k=200\ 000$ items uniformly out of $n=20\ 000\ 000$ items without replacement, where $d=7000$ of those items are "defective," the chance of not obtaining any defective item is the chance that the sample comes entirely from the non-defective items, of which there are $n-d$. Because all samples of size $k$ are equally probable, we only need to count the number of samples from $n-d$ items and divide by the total number of possible samples. Binomial coefficients (by definition) count samples, whence this chance is

$$\Pr(\text{no defective items in sample}) = \frac{\binom{n-d}{k}}{\binom{n}{k}}.$$

A (well-known) formula for binomial coefficients is

$$\binom{n}{k} = \frac{n!}{k!(n-k)!} = \frac{\Gamma(n+1)}{\Gamma(k+1)\Gamma(n-k+1)}.$$

When the numerator and denominator get large--which they quickly do--we must work with their logarithms to avoid numerical overflow on the computer. Using basic laws of logs, we obtain

$$\log\binom{n}{k} = \log \Gamma(n+1) - \log \Gamma(k+1) - \log \Gamma(n-k+1).$$

Therefore

$$\eqalign{ \log \frac{\binom{n-d}{k}}{\binom{n}{k}} &= \left[\log \Gamma(n-d+1) - \log \Gamma(k+1) - \log \Gamma(n-d-k+1)\right] \\ &- \left[\log \Gamma(n+1) - \log \Gamma(k+1) - \log \Gamma(n-k+1)\right]\\ &=\log\Gamma(n-d+1) -\log\Gamma(n-d-k+1) + \log\Gamma(n-k+1)-\log\Gamma(n+1). }$$

Excel offers a function to compute the natural logarithm of the Gamma function (which is why I rewrote the factorials in terms of $\Gamma$), GAMMALN. (We can trust this calculation because GAMMALN is easy to compute to high accuracy using its excellent asymptotic expansion and the recurrence relation $\Gamma(n+1)=n\Gamma(n)$; I strongly suspect that's how Excel does it.) Using this function, Excel finds the log to equal $-70.364\ 725\ldots$. Exponentiating this value gives $2.760\ 500\ 8\ldots\times 10^{-31}$ as the answer. Of course one will subtract the answer from $1$ at the end in order to express the complementary chance that at least one item in the sample is defective.


Because the final log of $-70. \ldots$ involved cancellation of values as large as $\log\Gamma(n+1)\approx 3\times 10^8$, some seven digits of precision were removed from the $15$ or so digits inherent in double-precision calculations. Therefore we cannot trust the answer to any more significant figures than I have written ($15-7=8$ of them). In this case the additional digits are $2.760\ 500\ \color{red}{804}\ldots$ whereas the correct answer (computed using exact arithmetic in Mathematica) is $2.760\ 500\ \color{red}{590}\ldots$, showing that this concern with loss of precision is justified.

$\endgroup$
1
  • 2
    $\begingroup$ Wow, this is a great explanation. Thank you so very much! $\endgroup$
    – Brandon
    Commented Sep 25, 2013 at 18:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.