The Binomial(10^{-6}, 10^{6}) distribution is shown below:

Though in theory this distribution extends from 0 to 10^{6}, in practice it is very much constrained to the lowest values. Calculating the probabilities for each possible value of this distribution is a problem, since one needs to evaluate 10^{6}!, which is beyond the capacity of most computers (Excel will calculate factorials up to 170, for example). The Stirling formula can be used to obtain a very good approximation of high factorials, but we still end up with the problem of having to deal with manipulating very high numbers. An easier approach is to use recursive formulas. These formulas relate the equation for the probability of the (*i*+1)^{th} value to the probability of the *i*^{th }value. Then one simply has to calculate the probability of any one value explicitly (a value chosen to give the simplest calculation) and then use the recursive formula to determine all other probabilities.

The binomial probability mass function gives:

{p(i)=\frac{n!}{i!(n-i)!}p^{i}(1-p)^{n-i}~~~\text{and}~~~p(i+1)=\frac{n!}{(i+1)!(n-i-1)!}p^{i+1}(1-p)^{n-i-1}} |

Thus,

{\frac{p\big(i+1\big)}{p\big(i\big)}=\frac{\frac{n!}{\big(i+1\big)!\big(n-i-1\big)!}p^{i+1}(1-p)^{n-i-1}}{\frac{n!}{i!\big(n-1\big)!}p^{i}\big(1-p\big)^{n-i}}=\frac{\big(n-i\big)p}{(i+1)(1-p)}} |

i.e.: {p(i+1)=\frac{(n-i)p}{(i+1)(1-p)}p(i)}

The binomial probability of zero successes is easily calculated:

{p(0)=(1-p)^{n}} |

so *p(1), p(2)*, etc. can be readily determined using the recursive formula:

{p(1)=(1-p)^{n}\frac{(n-0)p}{(0+1)(1-p)}=np(1-p)^{n-1}} |

{p(2)=np(1-p)^{n-1}\frac{(n-1)p}{(1+1)(1-p)}=\frac{n(n-1)p^{2}(1-p)^{n-2}}{2}} |

If this Binomial distribution is needed for a simulation, it is a simple enough task to use the *x* values and calculate probabilities as inputs into a Discrete distribution. The Discrete distribution then acts as if it was the required Binomial distribution. The spreadsheet of this model can be reached in the model Recursive Formulas.

The links to the Recursive Formulas software specific models are provided here:

If the binomial probability in this example had been extremely high, say 0.999 99, instead of 0.000 01, we would use the same technique but calculate backwards from p(1 000 000), i.e.

{p(n)=p^{n}} |

and work backwards using

{p(i)=\frac{(i+1)(1-p)}{(n-i)p}p(i+1)} |

Here are other useful recursive formulas for some of the most common discrete probability distributions:

{p(i+1)=\frac{\lambda}{i+1}p(i)} | {p(0)=e^{-\lambda}} | |

{p(i+1)=\frac{(s+i)(1-p)}{i+1}p(i)} | {p(s)=p^{s}} | |

{p(i+1)=(1-p)p(i)} | {p(1)=p} | |

{p(i+1)=\frac{(D-i)(n-i)}{(1+i)(1-D+M-n+i)}p(i)} | {p(0)=\frac{\binom{M-D}{n}}{\binom{M}{n}}} |

The formula for *p(*0*)* for the Hypergeometric distribution is unwieldy but can still be very accurately approximated without resorting to factorial calculations by using the Stirling formula to give:

{p(0)\approx \sqrt{\frac{(M-D)(M-n)}{M(M-D-n)}}\frac{(M-D)^{(M-D)}(M-n)^{(M-n)}}{(M-D-n)^{(M-D-n)}M^{M}}} |

This last formula will usually have to be calculated in log space first, then converted to real space at the end in order to avoid intermediary calculations of very large numbers. Another formula for *p(0)* which is only *very* slightly less accurate for larger *M* and generally more accurate for small *n* is provided by Cannon and Roe (1982):

{p(0)\approx =\Big(\frac{2M-2n-D+1}{2M-D+1}\Big)^{D}} |

(1)

This formula is produces by expanding the factorials for the equation for *p(0)* and cancelling common factors top and bottom to give:

{p(0)=\frac{\binom{M-D}{n}}{\binom{M}{n}}=\frac{(M-n)!(M-D)!}{M!(M-D-n)!}} |

(2)

{=\frac{(M-n-D+D)(M-n-D+(D+1))\dots(M-n-D+2)(M-n-D+1)}{(M-D+D)(M-D+(D+1))\dots(M-D+2)(M-D+1)}} |

We then take the first and last terms in each of the numerator and denominator, average the two terms and raise to the power of the number of terms (*D*) top and bottom:

{P(0)\approx \Big(\frac{(M-n-D+D)+(M-n-D+1)}{(M-D+D)+(M-D+1)}\Big)^{D}=\Big(\frac{2M-2n-D+1}{2M-D+1}\Big)^{D}} |

(3)

An alternative formulation is obtained by swapping around *n* and *D*, to give:

{p(0)\approx \Big(\frac{2M-2D-n+1}{2M-n+1}\Big)^{n}} |

This works since Equation (2) is symmetric in *n* and *D*, i.e. if *n* and *D* swap places the equation remains the same. Equation (1) is a better approximation when *n* > *D* and Equation (3) is better when *n* < *D*.