Movatterモバイル変換
[0]ホーム
probability distribution function in python
Tom Loredoloredo at spacenet.tn.cornell.edu
Wed Apr 11 14:14:37 EDT 2001
Howdy-Here is some "pure Python" (rather than an extension) implementingsome functions from the *Numerical Recipes* books related tofactorials. Use at your own risk! (i.e., not yet thoroughly tested).-Tom Loredofrom math import *import NumericN = 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_stpx = xx - 1.tmp = x + 5.5tmp = (x + 0.5)*log(tmp) - tmpser = 1.for j in range(6):x = x + 1.ser = ser + gammln_cof[j]/xreturn 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, xif (x == 0.):return(0.)ap = asum = 1. / adelta = sumn = 1while n <= itmax:ap = ap + 1.delta = delta * x / apsum = sum + deltaif (abs(delta) < abs(sum)*eps):return (sum * exp(-x + a*log(x) - gln), gln)n = n + 1raise 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 = xb0 = 0.b1 = 1.fac = 1.n = 1while n <= itmax:an = nana = an - aa0 = (a1 + a0*ana)*facb0 = (b1 + b0*ana)*facanf = an*faca1 = x*a0 + anf*a1b1 = x*b0 + anf*b1if (a1 != 0.):fac = 1. / a1g = b1*facif (abs((g-gold)/g) < eps):return (g*exp(-x+a*log(x)-gln), gln)gold = gn = n + 1raise 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-listmailing list
[8]ページ先頭