Computing Probability Functions

The binomial probability is
Direct evaluation of this expression is reasonable for small n. For larger n, multiplying out the factorials will lead to overflow problems.

Computer packages, including S-Plus and the statistics toolbox in Matlab, evaluate the probabilities using the representation

and use a log-gamma function to evaluate the log-factorials.

This is better than direct evaluation, but loss-of-precision can still occur, typically becoming noticeable for n about 10⁶. To overcome this limitation, I developed an alternative set of algorithms, based on saddle point methods and Stirling's formula, that significantly improve accuracy in large samples.

Similar problems occur for many common distributions. My DBinom library provides a stable implementation for the Binomial, Poisson, Hypergeometric, Gamma (incl. Chi-Square), Student's t, F, negative Binomial and Beta distributions.

• Source Code
• Documentation

The R programming environment incorporates a version of this code.

Example

The graph on the right shows an example, computing the binomial probability bin(3,n,2/n) for a range of values of n. As n increases to infinity, this probability should converge to the Poisson probability, p₀ = 4e^-2/3. The plot shows the number of digits agreement (measured by log₁₀|p/p₀-1|) for a range of n.

The saddle point algorithm (indicated by +) obtains the full accuracy of computer double precision (about 16 digits) for large n. The `classic' algorithm obtains similar results for n up to 10⁷, then loses precision after that.