This is the code block imported to perform various statistical functions. When run in standalone mode various tests results are returned. See the script header for accreditation for various functions.
'''
Distributions.py
Author
Dan Patterson
Dept of Geography and Environmental Studies
Carleton University, Ottawa, Canada
Dan_Patterson@carleton.ca
Date Feb 6, 2010
Modified June, 2013
Purpose:
Generate distribution related values for use with crosstabulation etc.
Sources:
cdf based on: http://www.johndcook.com/stand_alone_code.html
some of the gamma related functions derived from
Tom Loredo, translations of Numerical Recipes code
'''
#imports
import math, sys
#constants
sqrt_pi = math.sqrt(math.pi)
#gamma, lgamma and gamma incomplete methods.
def gamma_in(a, x, itmax=700, eps=3.e-7):
'''Series approx'n to the incomplete gamma function.'''
gln = math.lgamma(a)
if (x < 0.0):
raise bad_arg, x
if (x == 0.0):
return(0.0)
ap = a
sum = 1.0 / a
delta = sum
n = 1
while n <= itmax:
ap = ap + 1.0
delta = delta * x / ap
sum = sum + delta
if (abs(delta) < abs(sum)*eps):
return (sum * math.exp(-x + a*math.log(x) - gln), gln)
n = n + 1
raise max_iters, str((abs(delta), abs(sum)*eps))
def gamma_in_cf(a, x, itmax=200, eps=3.e-7):
'''Continued fraction approx'n of the incomplete gamma function.'''
gln = math.lgamma(a)
gold = 0.0
a0 = 1.0
a1 = x
b0 = 0.0
b1 = 1.0
fac = 1.0
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.0):
fac = 1.0 / a1
g = b1*fac
if (abs((g-gold)/g) < eps):
return (g*math.exp(-x + a*math.log(x) - gln), gln)
gold = g
n = n + 1
raise max_iters, str(abs((g-gold)/g))
def gamma_inc_p(a, x):
'''Incomplete gamma function, upper.'''
if (x < 0.0 or a <= 0.0):
raise ValueError, (a, x)
if (x < a+1.0):
return gamma_in(a,x)[0]
else:
return 1.0 - gamma_in_cf(a,x)[0]
def gamma_inc_q(a, x):
'''Incomplete gamma function, lower.'''
if (x < 0.0 or a <= 0.0):
raise ValueError, repr((a, x))
if (x < a+1.0):
return 1.0 - gamma_in(a,x)[0]
else:
return gamma_in_cf(a,x)[0]
# Chi squared distribution
def chisq_p(chisq, df):
'''Lower tail area of the chi**2 dist'n with df degrees of freedom.'''
half_df = 0.5 * df
half_chi = 0.5 * chisq
return gamma_inc_p(half_df, half_chi)
# ----------------------------------------------------------------
def chisq_q(chisq, df):
'''Upper tail area of the chi**2 dist'n with df degrees of freedom.'''
half_df = 0.5 * df
half_chi = 0.5 * chisq
return gamma_inc_q(half_df, half_chi)
def normal_cdf(z):
'''cumulative density function for gaussian (normal) distribution
input z-score
'''
# constants
a0, a1, a2, a3, a4 = [0.254829592, -0.284496736, 1.421413741,
-1.453152027, 1.061405429]
p = 0.3275911
# Save the sign of z
sign = 1
if z < 0:
sign = -1
z = abs(z) / math.sqrt(2.0)
# A&S formula 7.1.26
t = 1.0/(1.0 + p*z)
y = 1.0 - (((((a4*t + a3)*t) + a2)*t + a1)*t + a0)*t*math.exp(-z*z)
return 0.5 * (1.0 + sign*y)
def normal_inv_cdf(p):
'''Inverse of the normal cdf given a probability value
Source: http://www.johndcook.com/cpp_phi_inverse.html
'''
def log_one_plus_x(x):
if abs(x) > 1e-4:
return math.log(1.0 + x)
else:
return (-0.5*x + 1.0)*x
def ration_approx(p):
c = [2.515517, 0.802853, 0.010328, 1.432788, 0.189269, 0.001308]
return p - ((c[2]*p + c[1])*p + c[0]) / (((c[5]*p + c[4])*p +c[3])*p + 1.0)
#p = log_one_plus_x(p)
if p < 0.5:
return -ration_approx( math.sqrt(-2.0 * math.log(p)))
else:
return ration_approx( math.sqrt(-2.0 * math.log(1-p)))
def test_chisquare(chisq=None, df=None):
'''provide a chisquared value at df degrees of freedom'''
print "\nChi-squared test section"
if (chisq == None) or (df == None):
print "check the requirements"
sys.exit()
p_val_p = round(chisq_p(chisq, df), 4)
p_val_q = round(chisq_q(chisq, df), 4)
print "chi = %5s, df = % 3s, prob_p = %s, prob_q = %s" %( chisq, df, p_val_p, p_val_q)
def test_factorial_gamma(n_max=None):
'''factorial, gamma test'''
if n_max == None:
n_max = 10
print "\nFactorial test"
print "%12s %15s %15s" % ("n", "Factorial", "gamma")
n = 1
while n <= n_max:
print "%12s %15s %15s" % (n, math.factorial(n), math.gamma(n))
n += 1 #0.5
def test_gamma(lower=None, upper=None):
'''test gamma and lgamma'''
print "\nGamma, ln_gamma test"
print "%5s %15s %15s" % ("x", "gamma", "ln_gamma")
for x in range(lower,upper):
value1 = math.gamma(x)
value2 = math.lgamma(x)
print "%5s %15s %15s" %(x, value1, value2)
def test_z_cdf(zs=None):
'''test z-score and cdf functions'''
ps = [0.0000001, 0.00001, 0.001,
0.05, 0.15, 0.25, 0.35, 0.45,
0.5, 0.55, 0.65, 0.75, 0.85,
0.95, 0.999, 0.99999, 0.9999999]
if zs == None:
zs = [-4.0, -3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0]
print "\nTest z-score and cumulative density frequency"
print "\n%10s %12s" % ("z-score", "cdf")
for z in zs:
cdf = normal_cdf(z)
print "%10s %12.8f" % (z, cdf)
print "\n%10s %10s" %("prob.", "z")
for p in ps:
inv_cdf = normal_inv_cdf(p)
print "%10.8f %10.8f" % (p, inv_cdf)
#--------------------------------------------------------------------
if __name__ == "__main__":
test_z_cdf(zs=None)
test_factorial_gamma(10)
test_gamma(1,11)
test_chisquare(7.815, 3)
No comments:
Post a Comment