maggiexyz/PDL-Stats — pull request #1 by run4flat (2012-10-31 16:34:48+00:00)
Title: Use Stirling's approximation to expand range of Poisson pmf

The Poisson distribution is well behaved for large input values,
but the factorial form is not numerically stable and has a hard-coded
limit for its input values. Using Stirling's Approximation for the
logarithm of the factorial, we can rewrite the probability as a
sum of the logarithm of its terms, which is numerically stable.

This commit includes two new Poisson distribution functions as well
as modifications to the current distribution function. The factorial
form is called pmf_poisson_factorial and is a slight variant on the
original code that sets all probabilities for large input x to zero.
Perhaps a value of -1 would be better. The Stirling form is called
pmf_poisson_stirling and uses Stirling's approximation for all
values of x. This leads to numerically inaccurate results for very
small values of the input argument (though the inaccuracy is minor).
The Combined form is called pmf_poisson and uses the factorial form
for small values of x (roughly < 85) and the Stirling form for large
values of x.

Results from the Stirling and factorial forms differ, but that
difference falls within numerical bit noise for x roughly at
50.

This commit does not include tests. It is intended to demonstrate
the different forms. If accepted, tests should be fairly easy and
should follow shortly.