# Python math.gamma() Examples

The following are code examples for showing how to use math.gamma(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
def mvt_pdf(X, phi, K, nu):
"""
Multivariate student-t density:
output:
the density of the given element
input:
x = parameter (d dimensional numpy array or scalar)
mu = mean (d dimensional numpy array or scalar)
K = scale matrix (dxd numpy array)
nu = degrees of freedom
"""
d = X.shape[-1]
num = math.gamma((d + nu) / 2.) * pow(
1. + (1. / (nu - 2)) * ((X - phi).dot(np.linalg.inv(K)).dot(np.transpose(X - phi))), -(d + nu) / 2.)
denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi, d / 2.) * pow(np.linalg.det(K), 0.5)
return num / denom 
Example 2
def AR_1_Coeff(data):
"""
The autocovariance coefficient called as rho, for an AR(1) model can be calculated as shown here:

.. math::

\\rho(1) = \\frac{\\gamma(1)}{\\gamma(0)}

For further information look for example in "Zeitreihenanalyse", pages 17, by Matti Schneider, Sebastian Mentemeier,
SS 2010.

:param data: numerical list
:type data: list
:return: autocorrelation coefficient
:rtype: float
"""
return TimeSeries.acf(data, 1) / TimeSeries.acf(data, 0) 
Example 3
def test_cheb2_scheme(scheme):
def integrate_exact(k):
# \int_-1^1 x^k * sqrt(1 - x^2)
if k == 0:
return 0.5 * numpy.pi
if k % 2 == 1:
return 0.0
return (
numpy.sqrt(numpy.pi)
* ((-1) ** k + 1)
* math.gamma(0.5 * (k + 1))
/ math.gamma(0.5 * k + 2)
/ 4
)

degree = check_degree_1d(
lambda poly: scheme.integrate(poly, numpy.array([[-1.0], [1.0]])),
integrate_exact,
scheme.degree + 1,
)
assert degree >= scheme.degree
return 
Example 4
def test_integrate():
moments = quadpy.tools.integrate(lambda x: [x ** k for k in range(5)], -1, +1)
assert (moments == [2, 0, sympy.S(2) / 3, 0, sympy.S(2) / 5]).all()

lambda x: orthopy.line_segment.tree_legendre(x, 4, "monic", symbolic=True),
-1,
+1,
)
assert (moments == [2, 0, 0, 0, 0]).all()

# Example from Gautschi's "How to and how not to" article
lambda x: [x ** k * sympy.exp(-(x ** 3) / 3) for k in range(5)], 0, sympy.oo
)
S = numpy.vectorize(sympy.S)
gamma = numpy.vectorize(sympy.gamma)
n = numpy.arange(5)
reference = 3 ** (S(n - 2) / 3) * gamma(S(n + 1) / 3)
assert numpy.all([sympy.simplify(m - r) == 0 for m, r in zip(moments, reference)]) 
Example 5
def power_law_anisotropy(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
"""
equation (19) in Suyu+ 2010
:param r:
:return:
"""
# first term
theta_E = kwargs_profile['theta_E']
gamma = kwargs_profile['gamma']
r_ani = kwargs_anisotropy['r_ani']
a = 0.551 * kwargs_light['r_eff']
rho0_r0_gamma = self._rho0_r0_gamma(theta_E, gamma)
prefac1 = 4*np.pi * const.G * a**(-gamma) * rho0_r0_gamma / (3-gamma)
prefac2 = r * (r + a)**3/(r**2 + r_ani**2)
hyp1 = vel_util.hyp_2F1(a=2+gamma, b=gamma, c=3+gamma, z=1./(1+r/a))
hyp2 = vel_util.hyp_2F1(a=3, b=gamma, c=1+gamma, z=-a/r)
fac = r_ani**2/a**2 * hyp1 / ((2+gamma) * (r/a + 1)**(2+gamma)) + hyp2 / (gamma*(r/a)**gamma)
sigma2_dim_less = prefac1 * prefac2 * fac
return sigma2_dim_less * (self.cosmo.arcsec2phys_lens(1.) * const.Mpc / 1000)**2 
Example 6
def vel_disp(self, gamma, theta_E, r_eff, r_ani, rendering_number=1000):
"""
computes the averaged LOS velocity dispersion in the slit (convolved)

:param gamma: power-law slope of the mass profile (isothermal = 2)
:param theta_E: Einstein radius of the lens (in arcseconds)
:param r_eff: half light radius of the Hernquist profile (or as an approximation of any other profile to be described as a Hernquist profile
:param kwargs_aperture: keyword arguments describing the aperture of the collected spectral
:param rendering_number: number of spectral renderings drawn from the light distribution that go through the
slit of the observations

:return: LOS integrated velocity dispersion in units [km/s]
"""
sigma_s2_sum = 0
rho0_r0_gamma = self._rho0_r0_gamma(theta_E, gamma)
for i in range(0, rendering_number):
sigma_s2_draw = self.vel_disp_one(gamma, rho0_r0_gamma, r_eff, r_ani)
sigma_s2_sum += sigma_s2_draw
sigma_s2_average = sigma_s2_sum / rendering_number
return np.sqrt(sigma_s2_average) 
Example 7
def vel_disp_one(self, gamma, rho0_r0_gamma, r_eff, r_ani):
"""
computes one realisation of the velocity dispersion realized in the slit

:param gamma: power-law slope of the mass profile (isothermal = 2)
:param rho0_r0_gamma: combination of Einstein radius and power-law slope as equation (14) in Suyu+ 2010
:param r_eff: half light radius of the Hernquist profile (or as an approximation of any other profile to be described as a Hernquist profile
:param kwargs_aperture: keyword arguments describing the aperture of the collected spectral
:param FWHM: full width at half maximum of the seeing conditions, described as a Gaussian
:return: projected velocity dispersion of a single drawn position in the potential [km/s]
"""
a = 0.551 * r_eff
while True:
r = self.P_r(a)  # draw r
R, x, y = self.R_r(r)  # draw projected R
x_, y_ = self._psf.displace_psf(x, y)
bool = self.aperture.aperture_select(x_, y_)
if bool is True:
break
sigma_s2 = self.sigma_s2(r, R, r_ani, a, gamma, rho0_r0_gamma)
return sigma_s2 
Example 8
def __init__(self, params=None, core_model='moffat', aureole_model='power'):
"""
Parameters
----------
params : a dictionary containing keywords of PSF parameter
core_model : model of PSF core (moffat)
aureole_model : model of aureole ("power" or "multi-power")

"""
self.core_model = core_model
self.aureole_model = aureole_model

self.params = params

# Build attribute for parameters from dictionary keys
for key, val in params.items():
exec('self.' + key + ' = val')

if hasattr(self, 'fwhm'):
self.gamma = fwhm_to_gamma(self.fwhm, self.beta)
self.params['gamma'] = self.gamma

if hasattr(self, 'gamma'):
self.fwhm  = gamma_to_fwhm(self.gamma, self.beta)
self.params['fwhm'] = self.fwhm 
Example 9
def mvt_pdf(X, phi, K, nu):
"""
Multivariate student-t density:
output:
the density of the given element
input:
x = parameter (d dimensional numpy array or scalar)
mu = mean (d dimensional numpy array or scalar)
K = scale matrix (dxd numpy array)
nu = degrees of freedom
"""
d = X.shape[-1]
num = math.gamma((d + nu) / 2.) * pow(
1. + (1. / (nu - 2)) * ((X - phi).dot(np.linalg.inv(K)).dot(np.transpose(X - phi))), -(d + nu) / 2.)
denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi, d / 2.) * pow(np.linalg.det(K), 0.5)
return num / denom 
Example 10
 Project: tripp   Author: mjamesruggiero   File: hypothesis.py    BSD 3-Clause "New" or "Revised" License 5 votes
def B(alpha, beta):
"""a normalizing constant so that the total probability is 1"""
return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 
Example 11
def student_pdf_1d(X, phi, K, nu):
"""
Univariate student-t density:
output:
the density of the given element
input:
x = parameter scalar
mu = mean
var = scale matrix
nu = degrees of freedom
"""
num = math.gamma((1. + nu) / 2.) * pow(
1. + (1. / (nu - 2)) * (1. / K * (X - phi) ** 2), -(1. + nu) / 2.)
denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * K, 0.5)
return num / denom 
Example 12
def student_pdf_1d(X, mu, var, nu):
if nu > 50:
return gauss_pdf_1d(X, mu, var)
num = math.gamma((1. + nu) / 2.) * pow(
1. + (1. / (nu - 2)) * (1. / var * (X - mu) ** 2), -(1. + nu) / 2.)
denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * var, 0.5)
return num / denom 
Example 13
 Project: GPEXP   Author: goroda   File: utilities.py    GNU General Public License v2.0 5 votes
def distribFuncNDCircle(points, radius):
nDim = points.shape[1]
prob = 1.0/area
out = np.zeros((len(points)))
out[points[:,0]**2.0 + points[:,1]**2.0 < radius**2.0] = prob

return out 
Example 14
def ComputeScaleFunction(Mean, Shape):
return Mean/math.gamma(1 + 1.0/Shape) 
Example 15
def test_gamma(self):
self.assertAlmostEqual(math.gamma(0.5), math.sqrt(math.pi), places=15)
for i in xrange(1, 20):
self.assertEqual(math.factorial(i-1), math.gamma(i))
self.assertEqual(math.gamma(float('inf')), float('inf'))
self.assertRaises(ValueError, math.gamma, float('-inf'))
self.assertTrue(math.isnan(math.gamma(float('nan'))))
for i in xrange(0, -1001, -1):
self.assertRaises(ValueError, math.gamma, i) 
Example 16
def test_gamma_lgamma(self):
tolerance = 13
for x in itertools.count(0.001, 0.001):
if x > 5.0:
break
self.assertAlmostEqual(math.lgamma(x), math.log(math.gamma(x)), places=tolerance)
self.assertAlmostEqual(math.lgamma(x*x), math.log(math.gamma(x*x)), places=tolerance)
self.assertAlmostEqual(math.lgamma(2.0**x), math.log(math.gamma(2.0**x)), places=tolerance)

# Test negative values too, but not integers
if x % 1.0 != 0.0:
self.assertAlmostEqual(math.lgamma(-x), math.log(abs(math.gamma(-x))), places=tolerance)
self.assertAlmostEqual(math.lgamma(-x*x), math.log(abs(math.gamma(-x*x))), places=tolerance)
self.assertAlmostEqual(math.lgamma(-2.0**x), math.log(abs(math.gamma(-2.0**x))), places=tolerance) 
Example 17
def __init__(self, alpha):

self._alpha = np.array(alpha)
self._coef = gamma(np.sum(self._alpha)) / reduce(mul, [gamma(a) for a in self._alpha]) 
Example 18
def get_gamma(X, bandwidth):
x_norm = torch.sum(X ** 2, dim=1, keepdim=True)
x_t = torch.transpose(X, 0, 1)
x_norm_t = x_norm.view(1, -1)
t = x_norm + x_norm_t - 2.0 * torch.matmul(X, x_t)
dist2 = F.relu(Variable(t)).detach().data

d = dist2.cpu().numpy()
d = d[np.isfinite(d)]
d = d[d > 0]
median_dist2 = float(np.median(d))
gamma = 0.5 / median_dist2 / bandwidth
return gamma 
Example 19
def get_kernel_mat(x, landmarks, gamma):
d = pairwise_distances(x, landmarks)
k = torch.exp(d * -gamma)
k = k.view(x.shape[0], -1)
return k 
Example 20
def MMD(x, y, bandwidth=1.0):
y = y.detach()
gamma = get_gamma(y.detach(), bandwidth)
kxx = get_kernel_mat(x, x, gamma)
idx = torch.arange(0, x.shape[0], out=torch.LongTensor())
kxx = kxx * (1 - torch.eye(x.shape[0]).to(DEVICE))
kxx = torch.sum(kxx) / x.shape[0] / (x.shape[0] - 1)

kyy = get_kernel_mat(y, y, gamma)
idx = torch.arange(0, y.shape[0], out=torch.LongTensor())
kyy[idx, idx] = 0.0
kyy = torch.sum(kyy) / y.shape[0] / (y.shape[0] - 1)
kxy = torch.sum(get_kernel_mat(y, x, gamma)) / x.shape[0] / y.shape[0]
mmd = kxx + kyy - 2 * kxy
return mmd 
Example 21
def log_pdf(self, log_x):
return (self.alpha - 1) * log_x - self.beta * torch.exp(log_x) + self.alpha * np.log(self.beta) - np.log(math.gamma(self.alpha)) 
Example 22
def get_samples(self, num_samples):
rand_vars = np.random.gamma(self.alpha, 1.0 / self.beta, size=[num_samples, 1])  # numpy uses scale parameter
dtype=t_float).to(DEVICE) 
Example 23
 Project: pynam   Author: hbp-unibi   File: entropy.py    GNU General Public License v3.0 5 votes
def ncrr(x, y):
"""
Implementation of the binomial coefficient (ncr) function for real arguments
x and y (respectively corresponding to n and r).
"""
return math.gamma(x+1.0) / (math.gamma(y+1.0) * math.gamma(x-y+1.0)) 
Example 24
def B(alpha, beta):
"""a normalizing constant so that the total probability is 1"""
return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 
Example 25
 Project: mne-features   Author: mne-tools   File: univariate.py    BSD 3-Clause "New" or "Revised" License 5 votes
def _hurst_exp_helper(x, n_splits=20):
"""Helper function for :func:compute_hurst_exp.

Compute the Hurst exponent from a univariate time series. The Hurst
exponent is defined as the slope of the least-squares regression line
going through a cloud of n_splits points. Each point is obtained by
considering sub-series of x of n_splits different lenghts.

Parameters
----------
x : ndarray, shape (n_times,)

Returns
-------
output : ndarray, shape (n_splits,)
"""
n_times = x.shape[0]
_splits = np.floor(np.logspace(start=4, stop=np.log2(n_times / 2),
num=n_splits, base=2.))
splits = np.unique(_splits).astype(int)
reg = np.zeros((splits.size,))
for j, n in enumerate(splits):
a = x.copy()
d = int(floor(n_times / n))
a = np.lib.stride_tricks.as_strided(a, shape=(d, n),
strides=(n * a.strides[-1],
a.strides[-1]))
_rs = _hurst_exp_compute_rs(a)
_rs = _rs[~np.isnan(_rs)]
reg[j] = np.log(np.mean(_rs))
s = sum([sqrt((n - i) / i) for i in range(1, n)]) * ((n - 0.5) / n)
if n <= 340:
corr = (gamma((n - 1) / 2.) / (sqrt(np.pi) * gamma(n / 2.))) * s
else:
corr = ((n - 0.5) / n) * (1. / sqrt(np.pi * n / 2.)) * s
reg[j] -= (np.log(corr) - np.log(n) / 2)
return _slope_lstsq(np.log(splits), reg) 
Example 26
 Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    GNU General Public License v2.0 5 votes
def lower_incomplete_gamma2(a,x):
return gamma(a)-upper_incomplete_gamma2(a,x) 
Example 27
 Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    GNU General Public License v2.0 5 votes
def gammainc(a,x):
return lower_incomplete_gamma(a,x)/gamma(a) 
Example 28
 Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    GNU General Public License v2.0 5 votes
def gammaincc(a,x):
return upper_incomplete_gamma(a,x)/gamma(a) 
Example 29
 Project: python3_ios   Author: holzschu   File: bayes_update.py    BSD 3-Clause "New" or "Revised" License 5 votes
def beta_pdf(x, a, b):
return (x**(a-1) * (1-x)**(b-1) * math.gamma(a + b)
/ (math.gamma(a) * math.gamma(b))) 
Example 30
def test_update_exp_F():
for k in range(0,K):
BNMTF = bnmtf_vb_optimised(R,M,K,L,priors)
BNMTF.initialise(init_S,init_FG)
BNMTF.tauF = 4*numpy.ones((I,K))  # muF = [[0.5]], tauF = [[4.]]
BNMTF.update_exp_F(k) #-mu*sqrt(tau) = -0.5*2 = -1. lambda(1) = 0.241971 / (1-0.1587) = 0.2876155949126352. gamma = 0.37033832534958433
for i in range(0,I):
assert abs(BNMTF.expF[i,k] - (0.5 + 1./2. * 0.2876155949126352)) < 0.00001
assert abs(BNMTF.varF[i,k] - 1./4.*(1.-0.37033832534958433)) < 0.00001 
Example 31
def test_update_exp_S():
for k,l in itertools.product(xrange(0,K),xrange(0,L)):
BNMTF = bnmtf_vb_optimised(R,M,K,L,priors)
BNMTF.initialise(init_S,init_FG)
BNMTF.tauS = 4*numpy.ones((K,L)) # muS = [[1./3.]], tauS = [[4.]]
BNMTF.update_exp_S(k,l) #-mu*sqrt(tau) = -2./3., lambda(..) = 0.319448 / (1-0.2525) = 0.4273551839464883, gamma =
assert abs(BNMTF.expS[k,l] - (1./3. + 1./2. * 0.4273551839464883)) < 0.00001
assert abs(BNMTF.varS[k,l] - 1./4.*(1. - 0.4675359092102624)) < 0.00001 
Example 32
def test_update_exp_G():
for l in range(0,L):
BNMTF = bnmtf_vb_optimised(R,M,K,L,priors)
BNMTF.initialise(init_S,init_FG)
BNMTF.tauG = 4*numpy.ones((J,L)) # muG = [[1./4.]], tauG = [[4.]]
BNMTF.update_exp_G(l) #-mu*sqrt(tau) = -0.5., lambda(..) = 0.352065 / (1-0.3085) = 0.5091323210412148, gamma = 0.5137818808494219
for j in range(0,J):
assert abs(BNMTF.expG[j,l] - (1./4. + 1./2. * 0.5091323210412148)) < 0.0001
assert abs(BNMTF.varG[j,l] - 1./4.*(1. - 0.5137818808494219)) < 0.0001 
Example 33
def test_update_exp_U():
for k in range(0,K):
BNMF = bnmf_vb_optimised(R,M,K,priors)
BNMF.initialise()
BNMF.tauU = 4*numpy.ones((I,K)) # muU = [[0.5]], tauU = [[4.]]
BNMF.update_exp_U(k) #-mu*sqrt(tau) = -0.5*2 = -1. lambda(1) = 0.241971 / (1-0.1587) = 0.2876155949126352. gamma = 0.37033832534958433
for i in range(0,I):
assert abs(BNMF.expU[i,k] - (0.5 + 1./2. * 0.2876155949126352)) < 0.00001
assert abs(BNMF.varU[i,k] - 1./4.*(1.-0.37033832534958433)) < 0.00001 
Example 34
def factorial(n: Union[int, float]) -> Union[int, float]:
try:
import sympy.functions.combinatorial.factorials as fs
except ModuleNotFoundError:
if isinstance(n, int):
p = 1
for i in range(1, n + 1): p *= i
return p
else:
return math.gamma(n + 1)

return int(fs.factorial(n)) 
Example 35
def _alpha(n):
return math.pi**(n/2) / math.gamma(n/2 + 1) 
Example 36
def log_gamma(x):
"""
Returns the log-gamma value using Lanczos approximation formula

:param x: the point
:return: the log-gamma value for the point
"""
return math.lgamma(x) 
Example 37
def gamma(x):
"""
Returns the value of the gamma function: Gamma(x) = integral( t^(x-1) e^(-t),
t = 0 .. infinity)

:param x: the point
:return: the gamma value for the point
"""
return math.gamma(x) 
Example 38
def get_volume(radius, dim):
"""
Returns the volume of an N-ball of a certain radius.

:param dim: the number of dimensions to consider
:return: the resulting volume
"""
numerator = math.pow(math.pi, dim / 2.)
denominator = MathUtils.gamma((dim / 2.) + 1)
return radius2 * numerator / denominator 
Example 39
def plot(self, show_axes=True):
import matplotlib.pyplot as plt

ax = plt.gca()
plt.axis("equal")

if not show_axes:
ax.set_axis_off()

n = 2
I0 = 2 * math.factorial(n - 1) * math.pi ** (0.5 * n) / math.gamma(0.5 * n)

plot_disks(plt, self.points, self.weights, I0)
return 
Example 40
def plot(self, show_axes=True):
import matplotlib.pyplot as plt

ax = plt.gca()
plt.axis("equal")

if not show_axes:
ax.set_axis_off()

n = 2
I0 = 2 * math.factorial(n - 1) * math.pi ** (0.5 * n) / math.gamma(0.5 * n)

plot_disks(plt, self.points, self.weights, I0)
return 
Example 41
def integrate_monomial_over_enr(k):
if numpy.any(k % 2 == 1):
return 0
n = len(k)
return (
2
* math.factorial(sum(k) + n - 1)
* numpy.prod([math.gamma((kk + 1) / 2.0) for kk in k])
/ math.gamma((sum(k) + n) / 2)
)

# def simplex_monomials(degree):
#     exponents = numpy.concatenate(
#         [quadpy.helpers.partition(d, 2) for d in range(degree + 1)]
#     )
#
#     exact_vals = numpy.array(
#         [integrate_monomial_over_unit_simplex(k) for k in exponents]
#     )
#
#     def fun(x):
#         k = exponents.T
#         # <https://stackoverflow.com/a/46689653/353337>
#         s = x.shape[1:] + k.shape[1:]
#         return (
#             (x.reshape(x.shape[0], -1, 1) ** k.reshape(k.shape[0], 1, -1))
#             .prod(0)
#             .reshape(s)
#         )
#
#     return fun, exact_vals 
Example 42
def integrate_monomial_over_enr(k):
if numpy.any(k % 2 == 1):
return 0
n = len(k)
return (
2
* math.factorial(sum(k) + n - 1)
* numpy.prod([math.gamma((kk + 1) / 2.0) for kk in k])
/ math.gamma((sum(k) + n) / 2)
) 
Example 43
def integrate_monomial_over_enr2(k):
if numpy.any(k % 2 == 1):
return 0
return numpy.prod([math.gamma((kk + 1) / 2.0) for kk in k]) 
Example 44
def integrate_monomial_over_enr(k):
if numpy.any(k % 2 == 1):
return 0
n = len(k)
return (
2
* math.factorial(sum(k) + n - 1)
* numpy.prod([math.gamma((kk + 1) / 2.0) for kk in k])
/ math.gamma((sum(k) + n) / 2)
) 
Example 45
def gamma(x=('FloatPin', 0.0), Result=(REF, ('BoolPin', False))):
'''Return the Gamma function at x.'''
try:
Result(True)
return math.gamma(x)
except:
Result(False)
return -1 
Example 46
def B(alpha, beta):
"""a normalizing constant so that the total probability is 1"""
return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 
Example 47
def gamma(x):
if isinstance(x, LineValue): lx = x.get_value()
else: lx = x
if lx == NAN: return LineValue(NAN)
return LineValue(math.gamma(lx)) 
Example 48
def _rho0_r0_gamma(self, theta_E, gamma):
# equation (14) in Suyu+ 2010
return -1 * math.gamma(gamma/2.)/(np.sqrt(np.pi)*math.gamma((gamma-3)/2.)) * theta_E**gamma/self.cosmo.arcsec2phys_lens(theta_E) * self.cosmo.epsilon_crit * const.M_sun/const.Mpc**3  # units kg/m^3 
Example 49
def _rho0_r0_gamma(self, theta_E, gamma):
# equation (14) in Suyu+ 2010
return -1 * math.gamma(gamma/2) / (np.sqrt(np.pi)*math.gamma((gamma-3)/2.)) * theta_E ** gamma / \
self._cosmo.arcsec2phys_lens(theta_E) * self._cosmo.epsilon_crit * const.M_sun / const.Mpc ** 3 
Example 50
def sigma_r2(self, r, a, gamma, rho0_r0_gamma, r_ani):
"""
equation (19) in Suyu+ 2010
"""
# first term
prefac1 = 4*np.pi * const.G * a**(-gamma) * rho0_r0_gamma / (3-gamma)
prefac2 = r * (r + a)**3/(r**2 + r_ani**2)
hyp1 = vel_util.hyp_2F1(a=2+gamma, b=gamma, c=3+gamma, z=1./(1+r/a))
hyp2 = vel_util.hyp_2F1(a=3, b=gamma, c=1+gamma, z=-a/r)
fac = r_ani**2/a**2 * hyp1 / ((2+gamma) * (r/a + 1)**(2+gamma)) + hyp2 / (gamma*(r/a)**gamma)
return prefac1 * prefac2 * fac * (self._cosmo.arcsec2phys_lens(1.) * const.Mpc / 1000) ** 2 
Example 51
def gammainc_halfint(s, x) :
""" Lower incomplete gamma function =
integral from 0 to x of t ** (s-1) exp(-t) dt divided by gamma(s),
i.e., the fraction of gamma that you get if you integrate only until
x instead of all the way to infinity.
Implemented here only if s is a positive multiple of 0.5.
"""
# scipy equivalent: scipy.special.gammainc(s,x)

if s <= 0 :
raise ValueError('%s is not positive' % s)
if x < 0 :
raise ValueError('%s < 0' % x)
if s * 2 != int(s * 2) :
raise NotImplementedError('%s is not a multiple of 0.5' % s)

# Handle integers analytically
if s == int(s) :
term = 1
total = 1
for k in range(1, int(s)) :
term *= x /  k
total += term
return 1 - exp(-x) * total

# Otherwise s is integer + 0.5. Decrease to 0.5 using recursion formula:
result = 0.0
while s > 1 :
result -= x ** (s - 1) * exp(-x) / gamma(s)
s = s - 1
# Then use gammainc(0.5, x) = erf(sqrt(x))
result += erf(sqrt(x))
return result 
Example 52
def make_grid(self, image_size, pixel_scale=2.5):
""" Build grid for drawing """
self.image_size = image_size
self.yy, self.xx = np.mgrid[:image_size, :image_size]
self.cen = ((image_size-1)/2., (image_size-1)/2.)
self.pixel_scale = pixel_scale

for key, val in self.params.items():
if (key == 'gamma') | ('theta' in key):
val = val / pixel_scale
exec('self.' + key + '_pix' + ' = val') 
Example 53
def generate_core(self, folding_threshold=1.e-10):
""" Generate Galsim PSF of core. """
self.fwhm = self.gamma * 2. * np.sqrt(2**(1./self.beta)-1)
gsparams = galsim.GSParams(folding_threshold=folding_threshold)
psf_core = galsim.Moffat(beta=self.beta, fwhm=self.fwhm,
flux=1., gsparams=gsparams) # in arcsec
self.psf_core = psf_core
return psf_core 
Example 54
def draw_core2D_in_real(self, star_pos, Flux):
""" 2D drawing function of the core in real space given positions and flux (of core) of target stars """
from astropy.modeling.models import Moffat2D
gamma, alpha = self.gamma_pix, self.beta
Amps = np.array([moffat2d_Flux2Amp(gamma, alpha, Flux=flux)
for flux in Flux])
f_core_2d_s = np.array([Moffat2D(amplitude=amp, x_0=x0, y_0=y0,
gamma=gamma, alpha=alpha)
for ((x0,y0), amp) in zip(star_pos, Amps)])

return f_core_2d_s 
Example 55
def gamma_to_fwhm(gamma, beta):
""" in arcsec """
return gamma / fwhm_to_gamma(1, beta) 
Example 56
def moffat1d_normed(x, gamma, alpha):
""" Moffat for 1d array, flux normalized = 1 """
from astropy.modeling.models import Moffat1D
Mof_mod_1d = Moffat1D(amplitude=1, x_0=0, gamma=gamma, alpha=alpha)
y = Mof_mod_1d(x) / norm_mof
return y 
Example 57
def C_mof2Dto1D(r_core, beta):
""" gamma in pixel """
return 1. / (beta-1) * 2 * np.sqrt(np.pi) * r_core * math.gamma(beta) / math.gamma(beta - 1. / 2) 
Example 58
def C_mof1Dto2D(r_core, beta):
""" gamma in pixel """
return 1. / C_mof2Dto1D(r_core, beta) 
Example 59
 Project: recipes   Author: cwlseu   File: possion.py    GNU General Public License v3.0 5 votes
def bernonlli(p, n=1, k=1):
return C(n, k)*(p**k)*((1-p)**(n-k))

# gamma函数是阶乘在实数上的延拓，gamma(x+1) = x*gamma(x) 
Example 60
 Project: recipes   Author: cwlseu   File: possion.py    GNU General Public License v3.0 5 votes
def gamma(x):
if x <= 0:
raise ValueError("Gamma Parameter Error!")
if x == 0.5:
return math.pi**0.5
else:
return (x-1)*gamma(x-1) 
Example 61
def binomial(n, k):
"""Calculates the binomial coefficient."""
return math.gamma(n + 1) / (math.gamma(k + 1) * math.gamma(n - k + 1)) 
Example 62
def gamma(x):
return math.gamma(x) 
Example 63
def beta(a, b):
return (gamma(a) * gamma(b)) / gamma(a + b) 
Example 64
 Project: opytimizer   Author: gugarosa   File: distribution.py    GNU General Public License v3.0 5 votes
def generate_levy_distribution(beta=0.1, size=1):
"""Generates a n-dimensional array based on a Lévy distribution.

References:
X.-S. Yang and S. Deb. Computers & Operations Research.
Multiobjective Cuckoo Search for Design Optimization (2013).

Args:
beta (float): Skewness parameter.
size (int): Size of array.

Returns:
A Lévy distribution n-dimensional array.

"""

# Calculates the equation's numerator
num = gamma(1+beta) * sin(pi*beta/2)

# Calculates the equation's denominator
den = gamma((1+beta)/2) * beta * (2 ** ((beta-1)/2))

# Calculates the sigma for further distribution generation
sigma = (num/den) ** (1/beta)

# Calculates the 'u' distribution
u = r.generate_gaussian_random_number(size=size) * sigma

# Calculates the 'v' distribution
v = r.generate_gaussian_random_number(size=size)

# Finally, we can calculate the Lévy distribution
step = u / np.fabs(v) ** (1 / beta)

return step 
Example 65
def _factorial(n):
try:
return rmath.factorial(n)
except ValueError:
return rmath.gamma(n+1) 
Example 66
def fact(k: int) -> int:
"""Simple factorial.

>>> fact(1)
1
>>> fact(2)
2
>>> fact(3)
6
>>> fact(4)
24
"""
if k < 2:
return 1
return reduce(operator.mul, range(2,k+1))

# The implementation uses reduce( operator.mul, ... ) to compute
# the product of a sequence of integer values.
# We've included the @lru_cache because this is used often,
# and the small domain of possible values leads to some benefit
# from the cache.

# We could also use math.factorial(). In order to make use of
# the cache, we would need to do something like this.
#
# ..  parsed-literal::
#
#     fact = lru_cache(128)(math.factorial)
#
# This would create a similarly cached factorial function.

# Incomplete Gamma
# -----------------

# The incomplete (lower) gamma function is this:
#
# ..  math::
#
#     \gamma(s,z) = \sum_{k=0}^{\infty} \dfrac {(-1)^k} {k!} \; \dfrac {z^{s+k}} {s+k}
# 
Example 67
def Gamma_Half(k: float) -> float:
"""Gamma(k) with special case for k = n+1/2; k-1/2=n.

>>> import math
>>> round(Gamma_Half(2),1)
1.0
>>> round(Gamma_Half(3),1)
2.0
>>> round(Gamma_Half(4),1)
6.0
>>> round(Gamma_Half(5),1)
24.0

>>> round(Gamma_Half(.5), 7)
1.7724539
>>> round(math.sqrt(math.pi), 7)
1.7724539
>>> round(Gamma_Half(1.5), 7)
0.8862269
>>> round(math.sqrt(math.pi)/2, 7)
0.8862269
"""
ε = 1E-6
if abs(k-int(k)-.5) < ε:
n = int(k-.5)
return fact(2*n)/(4**n*fact(n))*math.sqrt(math.pi)
else:
return float(Gamma2(k))

# If the value is an :math:n+\dfrac{1}{2} \pm \epsilon, we'll use the special
# close-form expression. If the value is not close to :math:n+\dfrac{1}{2},
# we'll use a more general approximation.

# The math.gamma() Version
# ---------------------------

# Here's a test case for the math.gamma() function. 
Example 68
def student_pdf_1d(X, phi, K, nu):
"""
Univariate student-t density:
output:
the density of the given element
input:
x = parameter scalar
mu = mean
var = scale matrix
nu = degrees of freedom
"""
num = math.gamma((1. + nu) / 2.) * pow(
1. + (1. / (nu - 2)) * (1. / K * (X - phi) ** 2), -(1. + nu) / 2.)
denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * K, 0.5)
return num / denom 
Example 69
def student_pdf_1d(X, mu, var, nu):
if nu > 50:
return gauss_pdf_1d(X, mu, var)
num = math.gamma((1. + nu) / 2.) * pow(
1. + (1. / (nu - 2)) * (1. / var * (X - mu) ** 2), -(1. + nu) / 2.)
denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * var, 0.5)
return num / denom 
Example 70
def Slater_exec(Z_ind,orb):

'''

Define an executable Slater-type orbital wavefunction which takes only
the radius as an input argument. In this way, the usser specifies Z
and the orbital label string, and this generates a lambda function

:math:R(r) = (\\frac{2Z_{e}}{n_{e}})^{n_{e}} \\sqrt{\\frac{2Z_{e}}{n_{e}\\Gamma(2n_{e}+1)}} r^{n_{e}-1} e^{-\\frac{Z_{e} r}{n_{e}}}

*args*:

- **Z_ind**: int, atomic number

- **orb**: string, 'nlxx' as per standard orbital definitions used
throughout the library

*return*:

- executable function of position (in units of Angstrom)

***
'''
ne = ndic[int(orb[0])]
xi = Z_eff(Z_ind,orb)/ne
tmp = (2*xi)**ne*np.sqrt(2*xi/gamma(float(2*ne)+1))

def lambda_gen():

return lambda r: tmp*(r)**(ne-1)*np.exp(-xi*r)
return lambda_gen() 
Example 71
 Project: mmse-port   Author: rajivpoddar   File: mmse.py    BSD 2-Clause "Simplified" License 4 votes
def MMSESTSA(signal, fs, IS=0.25, W=1024, NoiseMargin=3, saved_params=None):
SP = 0.4
wnd = np.hamming(W)

y = segment(signal, W, SP, wnd)
Y = np.fft.fft(y, axis=0)
YPhase = np.angle(Y[0:int(np.fix(len(Y)/2))+1,:])
Y = np.abs(Y[0:int(np.fix(len(Y)/2))+1,:])
numberOfFrames = Y.shape[1]

NoiseLength = 9
NoiseCounter = 0
alpha = 0.99

NIS = int(np.fix(((IS * fs - W) / (SP * W) + 1)))
N = np.mean(Y[:,0:NIS].T).T

if saved_params != None:
NIS = 0
N = saved_params['N']
NoiseCounter = saved_params['NoiseCounter']

G = np.ones(N.shape)
Gamma = G

Gamma1p5 = math.gamma(1.5)
X = np.zeros(Y.shape)

for i in range(numberOfFrames):
Y_i = Y[:,i]

if i < NIS:
SpeechFlag = 0
NoiseCounter = 100
else:
SpeechFlag, NoiseCounter = vad(Y_i, N, NoiseCounter, NoiseMargin)

if SpeechFlag == 0:
N = (NoiseLength * N + Y_i) / (NoiseLength + 1)
LambdaD = (NoiseLength * LambdaD + (Y_i ** 2)) / (1 + NoiseLength)

gammaNew = (Y_i ** 2) / LambdaD
xi = alpha * (G ** 2) * Gamma + (1 - alpha) * np.maximum(gammaNew - 1, 0)

Gamma = gammaNew
nu = Gamma * xi / (1 + xi)

# log MMSE algo
#G = (xi/(1 + xi)) * np.exp(0.5 * expn(1, nu))

# MMSE STSA algo
G = (Gamma1p5 * np.sqrt(nu)) / Gamma * np.exp(-1 * nu / 2) * ((1 + nu) * bessel(0, nu / 2) + nu * bessel(1, nu / 2))
Indx = np.isnan(G) | np.isinf(G)
G[Indx] = xi[Indx] / (1 + xi[Indx])

X[:,i] = G * Y_i

output = OverlapAdd2(X, YPhase, W, SP * W)
return output, {'N': N, 'LambdaD': LambdaD, 'NoiseCounter': NoiseCounter} 
Example 72
 Project: data_algebra   Author: WinVector   File: SQLite.py    BSD 3-Clause "New" or "Revised" License 4 votes
def prepare_connection(self, conn):
# https://docs.python.org/3/library/sqlite3.html#sqlite3.Connection.create_function
# math fns
conn.create_function("acos", 1, math.acos)
conn.create_function("acosh", 1, math.acosh)
conn.create_function("asin", 1, math.asin)
conn.create_function("asinh", 1, math.asinh)
conn.create_function("atan", 1, math.atan)
conn.create_function("atanh", 1, math.atanh)
conn.create_function("ceil", 1, math.ceil)
conn.create_function("cos", 1, math.cos)
conn.create_function("cosh", 1, math.cosh)
conn.create_function("degrees", 1, math.degrees)
conn.create_function("erf", 1, math.erf)
conn.create_function("erfc", 1, math.erfc)
conn.create_function("exp", 1, math.exp)
conn.create_function("expm1", 1, math.expm1)
conn.create_function("fabs", 1, math.fabs)
conn.create_function("factorial", 1, math.factorial)
conn.create_function("floor", 1, math.floor)
conn.create_function("frexp", 1, math.frexp)
conn.create_function("gamma", 1, math.gamma)
conn.create_function("isfinite", 1, math.isfinite)
conn.create_function("isinf", 1, math.isinf)
conn.create_function("isnan", 1, math.isnan)
conn.create_function("lgamma", 1, math.lgamma)
conn.create_function("log", 1, math.log)
conn.create_function("log10", 1, math.log10)
conn.create_function("log1p", 1, math.log1p)
conn.create_function("log2", 1, math.log2)
conn.create_function("modf", 1, math.modf)
conn.create_function("sin", 1, math.sin)
conn.create_function("sinh", 1, math.sinh)
conn.create_function("sqrt", 1, math.sqrt)
conn.create_function("tan", 1, math.tan)
conn.create_function("tanh", 1, math.tanh)
conn.create_function("trunc", 1, math.trunc)
conn.create_function("atan2", 2, math.atan2)
conn.create_function("copysign", 2, math.copysign)
conn.create_function("fmod", 2, math.fmod)
conn.create_function("gcd", 2, math.gcd)
conn.create_function("hypot", 2, math.hypot)
conn.create_function("isclose", 2, math.isclose)
conn.create_function("ldexp", 2, math.ldexp)
conn.create_function("pow", 2, math.pow) 
Example 73
 Project: data_algebra   Author: WinVector   File: SQLite.py    BSD 3-Clause "New" or "Revised" License 4 votes
def prepare_connection(self, conn):
# https://docs.python.org/3/library/sqlite3.html#sqlite3.Connection.create_function
# math fns
conn.create_function("acos", 1, math.acos)
conn.create_function("acosh", 1, math.acosh)
conn.create_function("asin", 1, math.asin)
conn.create_function("asinh", 1, math.asinh)
conn.create_function("atan", 1, math.atan)
conn.create_function("atanh", 1, math.atanh)
conn.create_function("ceil", 1, math.ceil)
conn.create_function("cos", 1, math.cos)
conn.create_function("cosh", 1, math.cosh)
conn.create_function("degrees", 1, math.degrees)
conn.create_function("erf", 1, math.erf)
conn.create_function("erfc", 1, math.erfc)
conn.create_function("exp", 1, math.exp)
conn.create_function("expm1", 1, math.expm1)
conn.create_function("fabs", 1, math.fabs)
conn.create_function("factorial", 1, math.factorial)
conn.create_function("floor", 1, math.floor)
conn.create_function("frexp", 1, math.frexp)
conn.create_function("gamma", 1, math.gamma)
conn.create_function("isfinite", 1, math.isfinite)
conn.create_function("isinf", 1, math.isinf)
conn.create_function("isnan", 1, math.isnan)
conn.create_function("lgamma", 1, math.lgamma)
conn.create_function("log", 1, math.log)
conn.create_function("log10", 1, math.log10)
conn.create_function("log1p", 1, math.log1p)
conn.create_function("log2", 1, math.log2)
conn.create_function("modf", 1, math.modf)
conn.create_function("sin", 1, math.sin)
conn.create_function("sinh", 1, math.sinh)
conn.create_function("sqrt", 1, math.sqrt)
conn.create_function("tan", 1, math.tan)
conn.create_function("tanh", 1, math.tanh)
conn.create_function("trunc", 1, math.trunc)
conn.create_function("atan2", 2, math.atan2)
conn.create_function("copysign", 2, math.copysign)
conn.create_function("fmod", 2, math.fmod)
conn.create_function("gcd", 2, math.gcd)
conn.create_function("hypot", 2, math.hypot)
conn.create_function("isclose", 2, math.isclose)
conn.create_function("ldexp", 2, math.ldexp)
conn.create_function("pow", 2, math.pow) 
Example 74
def test_elbo():
I,J,K = 5,3,2
R = numpy.ones((I,J))
M = numpy.ones((I,J))
M[0,0], M[2,2], M[3,1] = 0, 0, 0 # size Omega = 12

lambdaU = 2*numpy.ones((I,K))
lambdaV = 3*numpy.ones((J,K))
alpha, beta = 3, 1
priors = { 'alpha':alpha, 'beta':beta, 'lambdaU':lambdaU, 'lambdaV':lambdaV }

expU = 5*numpy.ones((I,K))
expV = 6*numpy.ones((J,K))
varU = 11*numpy.ones((I,K))
varV = 12*numpy.ones((J,K))
exptau = 8.
explogtau = 9.

muU = 14*numpy.ones((I,K))
muV = 15*numpy.ones((J,K))
tauU = numpy.ones((I,K))/100.
tauV = numpy.ones((J,K))/101.
alpha_s = 20.
beta_s = 21.

# expU * expV = [[60]]
# (R - expU*expV)^2 = 12*59^2 = 41772
# Var[U*V] = 12*K*((11+5^2)*(12+6^2)-5^2*6^2) = 12*2*828 = 19872

# -muU*sqrt(tauU) = -14*math.sqrt(100) = -1.4
# -muV*sqrt(tauV) = -15*math.sqrt(101) = -1.4925557853149838
# cdf(-1.4) = 0.080756659233771066
# cdf(-1.4925557853149838) = 0.067776752211548219

ELBO = 12./2.*(explogtau - math.log(2*math.pi)) - 8./2.*(41772+19872) \
+ 5*2*(math.log(2.) - 2.*5.) + 3*2*(math.log(3.) - 3.*6.) \
+ 3.*numpy.log(1.) - numpy.log(math.gamma(3.)) + 2.*9. - 1.*8. \
- 20.*numpy.log(21.) + numpy.log(math.gamma(20.)) - 19.*9. + 21.*8. \
- 0.5*5*2*math.log(1./100.) + 0.5*5*2*math.log(2*math.pi) + 5*2*math.log(1.-0.080756659233771066) \
+ 0.5*5*2*1./100.*(11.+81.) \
- 0.5*3*2*math.log(1./101.) + 0.5*3*2*math.log(2*math.pi) + 3*2*math.log(1.-0.067776752211548219) \
+ 0.5*3*2*1./101.*(12.+81.)

BNMF = bnmf_vb_optimised(R,M,K,priors)
BNMF.expU = expU
BNMF.expV = expV
BNMF.varU = varU
BNMF.varV = varV
BNMF.exptau = exptau
BNMF.explogtau = explogtau
BNMF.muU = muU
BNMF.muV = muV
BNMF.tauU = tauU
BNMF.tauV = tauV
BNMF.alpha_s = alpha_s
BNMF.beta_s = beta_s
assert BNMF.elbo() == ELBO 
Example 75
def stroud_1967_5():

def f(x):
degree = 5
n = 7

lmbda, xi, mu, gamma = x
eta = 0
A = 1 / 9
B = 1 / 72
C = B

# data = [
#     (B, rd(n, [(+lmbda, 1), (+xi, n - 1)])),
#     (B, rd(n, [(-lmbda, 1), (-xi, n - 1)])),
#     (C, rd(n, [(+mu, 2), (+gamma, n - 2)])),
#     (C, rd(n, [(-mu, 2), (-gamma, n - 2)])),
#     (2 * A, numpy.full((1, n), eta)),
# ]
# points, weights = untangle(data)
# weights *= numpy.sqrt(numpy.pi) ** n

data = [
(B, rd(n, [(+lmbda, 1), (+xi, n - 1)])),
(B, rd(n, [(-lmbda, 1), (-xi, n - 1)])),
(C, rd(n, [(+mu, 2), (+gamma, n - 2)])),
(C, rd(n, [(-mu, 2), (-gamma, n - 2)])),
(2 * A, numpy.full((1, n), eta)),
]

points, weights = untangle(data)
weights *= numpy.sqrt(numpy.pi) ** n

A = numpy.concatenate(orthopy.enr2.tree(points.T, degree, symbolic=False))

out = numpy.dot(A, weights)
out[0] -= numpy.sqrt(numpy.sqrt(numpy.pi)) ** n

norm_v = numpy.sqrt(numpy.vdot(out, out))
return norm_v

x0 = [
2.009_505_637_083_749e00,
2.774_548_295_173_737e-01,
-1.062_215_595_206_724e00,
6.698_352_123_613_097e-01,
# 2.009_505_6,
# 0.277_454_83,
# -1.062_215_60,
# 0.669_835_21,
]

out = minimize(f, x0, method="Powell", tol=1.0e-20, options={"maxiter": 20000})
print(out.status, out.nfev, out.message, "Function value", out.fun)
assert out.success
print()
for x in out.x:
print(f"{x:.15e}")
return 
Example 76
def test_gautschi_how_to_and_how_not_to():
"""Test Gautschi's famous example from

W. Gautschi,
How and how not to check Gaussian quadrature formulae,
BIT Numerical Mathematics,
June 1983, Volume 23, Issue 2, pp 209–216,
<https://doi.org/10.1007/BF02218441>.
"""
points = numpy.array(
[
1.457697817613696e-02,
8.102669876765460e-02,
2.081434595902250e-01,
3.944841255669402e-01,
6.315647839882239e-01,
9.076033998613676e-01,
1.210676808760832,
1.530983977242980,
1.861844587312434,
2.199712165681546,
2.543839804028289,
2.896173043105410,
3.262066731177372,
3.653371887506584,
4.102376773975577,
]
)
weights = numpy.array(
[
3.805398607861561e-2,
9.622028412880550e-2,
1.572176160500219e-1,
2.091895332583340e-1,
2.377990401332924e-1,
2.271382574940649e-1,
1.732845807252921e-1,
9.869554247686019e-2,
3.893631493517167e-2,
9.812496327697071e-3,
1.439191418328875e-3,
1.088910025516801e-4,
3.546866719463253e-6,
3.590718819809800e-8,
5.112611678291437e-11,
]
)

# weight function exp(-t**3/3)
n = len(points)
moments = numpy.array(
[3.0 ** ((k - 2) / 3.0) * math.gamma((k + 1) / 3.0) for k in range(2 * n)]
)

errors_alpha, errors_beta = quadpy.tools.check_coefficients(moments, alpha, beta)

assert numpy.max(errors_alpha) > 1.0e-2
assert numpy.max(errors_beta) > 1.0e-2 
Example 77
def cdf(x: float, k: int) -> float:
"""χ² cumulative distribution function.

:param x: χ² value -- generally sum (obs[i]-exp[i])**2/exp[i]
for parallel sequences of observed and expected values.
:param k: degrees of freedom >= 1; generally len(data)-1

From http://en.wikipedia.org/wiki/Chi-squared_distribution

>>> round(cdf(0.004, 1), 2)
0.95
>>> round(cdf(10.83, 1), 3)
0.001
>>> round(cdf(3.94, 10), 2)
0.95
>>> round(cdf(29.59, 10), 3)
0.001
>>> expected=[0.95, 0.90, 0.80, 0.70, 0.50, 0.30, 0.20, 0.10, 0.05, 0.01, 0.001]
>>> chi2= [0.004, 0.02, 0.06, 0.15, 0.46, 1.07, 1.64, 2.71, 3.84, 6.64, 10.83]
>>> act= [round(x,3) for x in map(cdf, chi2, [1]*len(chi2))]
>>> act
[0.95, 0.888, 0.806, 0.699, 0.498, 0.301, 0.2, 0.1, 0.05, 0.01, 0.001]

From http://www.itl.nist.gov/div898/handbook/prc/section4/prc45.htm

>>> round(cdf(19.18, 6), 5)
0.00387
>>> round(cdf(12.5916, 6), 2)
0.05

From http://www.itl.nist.gov/div898/handbook/prc/section4/prc46.htm
>>> round(cdf(12.131, 4), 5)
0.01639
>>> round(cdf(9.488, 4), 2)
0.05

"""
return 1-gamma(k/2, x/2)/Gamma_Half(k/2)

# The calcuation is 1 minus the ratio the partial
# gamma to the full gamma.

# Unit Test Cases
# ===============

# We'll use the doctest comments in each function defined above.
# Additionally, we'll use some strings with doctest test cases.