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