from math import cos, e, factorial, log, pi, sin, sqrt
def integrate_sr(func_str, var, a, b, n = 200):
'''Composite Simpson's rule
Numerical Integration'''
h = float(b-a) / n
vars()[var] = a; ans = eval(func_str)
vars()[var] = b; ans += eval(func_str)
times, var_val = [(unused % 2) + 1 for unused in xrange(1, n)], [a + h * unused for unused in xrange(1, n)]
for unused in xrange(len(var_val)):
vars()[var] = var_val[unused]
ans += 2 * times[unused] * eval(func_str)
return ans * h / 3.
def integrate_g7(func_str, var, a, b):
'''Gauss Integration with 7 nodes
Numerical Integration'''
f = eval('lambda ' + var + ':' + func_str)
g = lambda unused:(b-a) / 2. * f((b-a) / 2. * unused + (a + b) / 2.)
nodes = [-0.949107912343, -0.741531185599, -0.405845151377, 0.0, 0.405845151377, 0.741531185599, 0.949107912343]
weight = [0.129484966169, 0.279705391489, 0.381830050505, 0.417959183673, 0.381830050505, 0.279705391489, 0.129484966169]
return sum([weight[x] * g(nodes[x]) for x in xrange(7)])
def integrate_k15(func_str, var, a, b):
'''Kronrod Integration with 15 nodes
Numerical Integration'''
f = eval('lambda ' + var + ':' + func_str)
g = lambda unused:(b-a) / 2. * f((b-a) / 2. * unused + (a + b) / 2.)
nodes = [-0.991455371121, -0.949107912343, -0.86486442336, -0.741531185599, -0.586087235468, -0.405845151377, -0.207784955008, 0.0, 0.207784955008, 0.405845151377, 0.586087235468, 0.741531185599, 0.86486442336, 0.949107912343, 0.991455371121]
weight = [0.0229353220105, 0.06309209263, 0.104790010322, 0.140653259716, 0.169004726639, 0.190350578065, 0.204432940075, 0.209482141085, 0.204432940075, 0.190350578065, 0.169004726639, 0.140653259716, 0.104790010322, 0.06309209263, 0.0229353220105]
return sum([weight[x] * g(nodes[x]) for x in xrange(15)])
def cumSum(l):
for x in xrange(1, len(l)):
l[x] += int(l[x-1])
return l
def delta_list(l):
m = []
for x in xrange(1, len(l)):
m += [l[x] -l[x-1]]
return m
def double_factorial(n):
if int(n) == n:
n = int(n)
if [0, 1].__contains__(n):
return 1
a = (n & 1) + 2
b = 1
while a <= n:
b *= a
a += 2
return float(b)
else:
return factorials(n / 2) * 2 ** (n / 2) * (pi / 2) ** (.25 * (-1 + cos(n * pi)))
def factorials(n):
'''Theoretically gets the factorial of any number, but since
recursions are limited, best use only integers or floats that have
decimals of .5'''
if n == int(n):
return float(factorial(n))
return pi ** (.5 * sin(n * pi) ** 2) * 2 ** (-n + .25 * (-1 + cos(2 * n * pi))) * double_factorial(2 * n)
def normalpdf(x, mu = 0, sd = 1):
return 1 / sqrt(2 * pi) * e ** (-.5 * ((x-mu) / sd) ** 2) / sd
def normalcdf(a, b, mu = 0, sd = 1):
return integrate_sr('normalpdf(x, ' + str(mu) + ', ' + str(sd) + ')', 'x', a, b)
def binompdf(n, p, k = -1):
if k != -1:
return nCr(n, k) * p ** k * (1-p) ** (n-k)
else:
return [nCr(n, k) * p ** k * (1-p) ** (n-k) for k in xrange(n + 1)]
def binomcdf(n, p, k = -1):
if k == -1:
return cumSum(binompdf(n, p))
else:
return sum(binompdf(n, p)[:k + 1])
def tpdf(t, df):
return (1 + float(t) ** 2 / df) ** (-.5 * (df + 1)) * factorials((df + 1) / 2.-1) / (sqrt(df * pi) * factorials(df / 2.-1))
def tcdf(a, b, df):
return integrate_sr('tpdf(x, ' + str(df) + ')', 'x', a, b)
def chi2pdf(x, k):
return .5 ** (k / 2.) * x ** (k / 2.-1) * e ** (-x / 2.) / factorials(k / 2.-1)
def chi2cdf(a, b, k):
return integrate_sr('chi2pdf(x, ' + str(k) + ')', 'x', a, b)
def poissonpdf(mean_occur, expected):
return e ** -mean_occur * mean_occur ** expected / factorials(expected)
def poissoncdf(mean_occur, max_expected):
return sum([poissonpdf(mean_occur, k) for k in xrange(max_expected + 1)])
def geometpdf(p, n):
if p > 1:
raise ValueError('P cannot be greater than 1')
return p * ((1 - p) ** (n - 1))
def geometcdf(p, n):
return sum([geometpdf(p, i) for i in xrange(1, n)])
def Fpdf(x, d1, d2):
x = float(x)
d1 = float(d1)
d2 = float(d2)
return ((d1 * x) ** d1 * d2 ** d2 / (d1 * x + d2) ** (d1 + d2)) ** .5 / (x * beta(d1 / 2, d2 / 2))
def Fcdf(a, b, d1, d2):
'''Doesnt work for for a = 0; use .01
Suggest use 100 for max b'''
return integrate_sr('Fpdf(x, ' + str(d1) + ', ' + str(d2) + ')', 'x', a, b)