probability distribution function in python
Tom Loredo
loredo at spacenet.tn.cornell.edu
Wed Apr 11 14:14:37 EDT 2001
Howdy-
Here is some "pure Python" (rather than an extension) implementing
some functions from the *Numerical Recipes* books related to
factorials. Use at your own risk! (i.e., not yet thoroughly tested).
-Tom Loredo
from math import *
import Numeric
N = Numeric
#============= Exceptions ===============
max_iters = 'Too many iterations: '
#============= Global constants ===============
rt2 = sqrt(2.)
gammln_cof = N.array([76.18009173, -86.50532033, 24.01409822,
-1.231739516e0, 0.120858003e-2, -0.536382e-5])
gammln_stp = 2.50662827465
#============= Gamma, Incomplete Gamma ===========
def gammln(xx):
"""Logarithm of the gamma function."""
global gammln_cof, gammln_stp
x = xx - 1.
tmp = x + 5.5
tmp = (x + 0.5)*log(tmp) - tmp
ser = 1.
for j in range(6):
x = x + 1.
ser = ser + gammln_cof[j]/x
return tmp + log(gammln_stp*ser)
def gser(a, x, itmax=700, eps=3.e-7):
"""Series approx'n to the incomplete gamma function."""
gln = gammln(a)
if (x < 0.):
raise bad_arg, x
if (x == 0.):
return(0.)
ap = a
sum = 1. / a
delta = sum
n = 1
while n <= itmax:
ap = ap + 1.
delta = delta * x / ap
sum = sum + delta
if (abs(delta) < abs(sum)*eps):
return (sum * exp(-x + a*log(x) - gln), gln)
n = n + 1
raise max_iters, str((abs(delta), abs(sum)*eps))
def gcf(a, x, itmax=200, eps=3.e-7):
"""Continued fraction approx'n of the incomplete gamma function."""
gln = gammln(a)
gold = 0.
a0 = 1.
a1 = x
b0 = 0.
b1 = 1.
fac = 1.
n = 1
while n <= itmax:
an = n
ana = an - a
a0 = (a1 + a0*ana)*fac
b0 = (b1 + b0*ana)*fac
anf = an*fac
a1 = x*a0 + anf*a1
b1 = x*b0 + anf*b1
if (a1 != 0.):
fac = 1. / a1
g = b1*fac
if (abs((g-gold)/g) < eps):
return (g*exp(-x+a*log(x)-gln), gln)
gold = g
n = n + 1
raise max_iters, str(abs((g-gold)/g))
def gammp(a, x):
"""Incomplete gamma function."""
if (x < 0. or a <= 0.):
raise ValueError, (a, x)
if (x < a+1.):
return gser(a,x)[0]
else:
return 1.-gcf(a,x)[0]
def gammq(a, x):
"""Incomplete gamma function."""
if (x < 0. or a <= 0.):
raise ValueError, repr((a, x))
if (x < a+1.):
return 1.-gser(a,x)[0]
else:
return gcf(a,x)[0]
def factrl(n, ntop=0, prev=N.ones((33),N.Float)):
"""Factorial of n.
The first 33 values are stored as they are calculated to
speed up subsequent calculations."""
if n < 0:
raise ValueError, 'Negative argument!'
elif n <= ntop:
return prev[n]
elif n <= 32:
for j in range(ntop+1, n+1):
prev[j] = j * prev[j-1]
ntop = n
return prev[n]
else:
return exp(gammln(n+1.))
def factln(n, prev=N.array(101*(-1.,))):
"""Log factorial of n.
Values for n=0 to 100 are stored as they are calculated to
speed up subsequent calls."""
if n < 0:
raise ValueError, 'Negative argument!'
elif n <= 100:
if prev[n] < 0:
prev[n] = gammln(n+1.)
return prev[n]
else:
return gammln(n+1.)
def combln(Ntot, n):
"""Log of number of combinations of n items from Ntot total."""
return factln(Ntot) - factln(n) - factln(Ntot-n)
More information about the Python-list
mailing list