# Python scipy.special.gamma() Examples

The following are 30 code examples of scipy.special.gamma(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module , or try the search function .
Example #1
 Source File: reciprocal_kprime.py    From burnman with GNU General Public License v2.0 6 votes
```def _upper_incomplete_gamma(z, a):
"""
An implementation of the non-regularised upper incomplete gamma
function. Computed using the relationship with the regularised
lower incomplete gamma function (scipy.special.gammainc).
Uses the recurrence relation wherever z<0.
"""
n = int(-np.floor(z))
if n > 0:
z = z + n
u_gamma = (1. - gammainc(z, a))*gamma(z)

for i in range(n):
z = z - 1.
u_gamma = (u_gamma - np.power(a, z)*np.exp(-a))/z
return u_gamma
else:
return (1. - gammainc(z, a))*gamma(z) ```
Example #2
 Source File: test_Calculus.py    From PySpice with GNU General Public License v3.0 6 votes
```def compute_finite_difference_coefficients(derivative_order, grid_size):

# from http://www.scientificpython.net/pyblog/uniform-finite-differences-all-orders

n = 2*grid_size -1
A = np.tile(np.arange(grid_size), (n,1)).T
B = np.tile(np.arange(1-grid_size,grid_size), (grid_size,1))
M = (B**A)/gamma(A+1)

r = np.zeros(grid_size)
r[derivative_order] = 1

D = np.zeros((grid_size, grid_size))
for k in range(grid_size):
indexes = k + np.arange(grid_size)
D[:,k] = solve(M[:,indexes], r)

return D

#################################################################################################### ```
Example #3
```def _pdf_skip(self, x, dfn, dfd, nc):
# ncf.pdf(x, df1, df2, nc) = exp(nc/2 + nc*df1*x/(2*(df1*x+df2))) *
#             df1**(df1/2) * df2**(df2/2) * x**(df1/2-1) *
#             (df2+df1*x)**(-(df1+df2)/2) *
#             gamma(df1/2)*gamma(1+df2/2) *
#             L^{v1/2-1}^{v2/2}(-nc*v1*x/(2*(v1*x+v2))) /
#             (B(v1/2, v2/2) * gamma((v1+v2)/2))
n1, n2 = dfn, dfd
term = -nc/2+nc*n1*x/(2*(n2+n1*x)) + sc.gammaln(n1/2.)+sc.gammaln(1+n2/2.)
term -= sc.gammaln((n1+n2)/2.0)
Px = np.exp(term)
Px *= n1**(n1/2) * n2**(n2/2) * x**(n1/2-1)
Px *= (n2+n1*x)**(-(n1+n2)/2)
Px *= sc.assoc_laguerre(-nc*n1*x/(2.0*(n2+n1*x)), n2/2, n1/2-1)
Px /= sc.beta(n1/2, n2/2)
# This function does not have a return.  Drop it for now, the generic
# function seems to work OK. ```
Example #4
```def _pdf(self, x, df, nc):
# nct.pdf(x, df, nc) =
#                               df**(df/2) * gamma(df+1)
#                ----------------------------------------------------
#                2**df*exp(nc**2/2) * (df+x**2)**(df/2) * gamma(df/2)
n = df*1.0
nc = nc*1.0
x2 = x*x
ncx2 = nc*nc*x2
fac1 = n + x2
trm1 = n/2.*np.log(n) + sc.gammaln(n+1)
trm1 -= n*np.log(2)+nc*nc/2.+(n/2.)*np.log(fac1)+sc.gammaln(n/2.)
Px = np.exp(trm1)
valF = ncx2 / (2*fac1)
trm1 = np.sqrt(2)*nc*x*sc.hyp1f1(n/2+1, 1.5, valF)
trm1 /= np.asarray(fac1*sc.gamma((n+1)/2))
trm2 = sc.hyp1f1((n+1)/2, 0.5, valF)
trm2 /= np.asarray(np.sqrt(fac1)*sc.gamma(n/2+1))
Px *= trm1+trm2
return Px ```
Example #5
```def _munp(self, n, beta, m):
"""
Returns the n-th non-central moment of the crystalball function.
"""
N = 1.0 / (m/beta / (m-1) * np.exp(-beta**2 / 2.0) + _norm_pdf_C * _norm_cdf(beta))

def n_th_moment(n, beta, m):
"""
Returns n-th moment. Defined only if n+1 < m
Function cannot broadcast due to the loop over n
"""
A = (m/beta)**m * np.exp(-beta**2 / 2.0)
B = m/beta - beta
rhs = 2**((n-1)/2.0) * sc.gamma((n+1)/2) * (1.0 + (-1)**n * sc.gammainc((n+1)/2, beta**2 / 2))
lhs = np.zeros(rhs.shape)
for k in range(n + 1):
lhs += sc.binom(n, k) * B**(n-k) * (-1)**k / (m - k - 1) * (m/beta)**(-m + k + 1)
return A * lhs + rhs

return N * _lazywhere(np.atleast_1d(n + 1 < m),
(n, beta, m),
np.vectorize(n_th_moment, otypes=[np.float]),
np.inf) ```
Example #6
```def mean(t, a, b):
# TODO this is not tested yet.
# tests:
#    cemean(0., a, b)==mean(a, b, p)
#    mean(t, a, 1., p)==mean(0., a, 1., p) == a
# conditional excess mean
# E[Y|y>t]
# (conditional mean age at failure)
# http://reliabilityanalyticstoolkit.appspot.com/conditional_distribution
from scipy.special import gamma
from scipy.special import gammainc
# Regularized lower gamma
print('not tested')

v = 1. + 1. / b
gv = gamma(v)
L = np.power((t + .0) / a, b)
cemean = a * gv * np.exp(L) * (1 - gammainc(v, t / a) / gv)

return cemean ```
Example #7
```def volume_unit_ball(n_dims):
return (pi ** (.5 * n_dims)) / gamma(.5 * n_dims + 1) ```
Example #8
 Source File: test_pre.py    From skan with BSD 3-Clause "New" or "Revised" License 5 votes
```def ball_volume(ndim, radius):
"""Return the volume of a ball of dimension `ndim` and radius `radius`."""
n = ndim
return np.pi ** (n / 2) / special.gamma(n / 2 + 1) * r ** n ```
Example #9
```def levy(self, D):
r"""Levy function.

Returns:
float: Next Levy number.
"""
sigma = (Gamma(1 + self.beta) * sin(pi * self.beta / 2) / (Gamma((1 + self.beta) / 2) * self.beta * 2 ** ((self.beta - 1) / 2))) ** (1 / self.beta)
return 0.01 * (self.normal(0, 1, D) * sigma / fabs(self.normal(0, 1, D)) ** (1 / self.beta)) ```
Example #10
```def factorial(n):
return gamma(n + 1) ```
Example #11
```def _hs(k, cs, rho, omega):
c0 = (cs * cs * (1 + rho * rho) / (1 - rho * rho) /
(1 - 2 * rho * rho * cos(2 * omega) + rho ** 4))
gamma = (1 - rho * rho) / (1 + rho * rho) / tan(omega)
ak = abs(k)
return c0 * rho ** ak * (cos(omega * ak) + gamma * sin(omega * ak)) ```
Example #12
```def from_spline(cls, tck, extrapolate=None):
"""
Construct a piecewise polynomial from a spline

Parameters
----------
tck
A spline, as returned by `splrep` or a BSpline object.
extrapolate : bool or 'periodic', optional
If bool, determines whether to extrapolate to out-of-bounds points
based on first and last intervals, or to return NaNs.
If 'periodic', periodic extrapolation is used. Default is True.
"""
if isinstance(tck, BSpline):
t, c, k = tck.tck
if extrapolate is None:
extrapolate = tck.extrapolate
else:
t, c, k = tck

cvals = np.empty((k + 1, len(t)-1), dtype=c.dtype)
for m in xrange(k, -1, -1):
y = fitpack.splev(t[:-1], tck, der=m)
cvals[k - m, :] = y/spec.gamma(m+1)

return cls.construct_fast(cvals, t, extrapolate) ```
Example #13
```def fromspline(cls, xk, cvals, order, fill=0.0):
# Note: this spline representation is incompatible with FITPACK
N = len(xk)-1
sivals = np.empty((order+1, N), dtype=float)
for m in xrange(order, -1, -1):
fact = spec.gamma(m+1)
res = _fitpack._bspleval(xk[:-1], xk, cvals, order, m)
res /= fact
sivals[order-m, :] = res
return cls(sivals, xk, fill=fill)

# The 3 private functions below can be called by splmake(). ```
Example #14
```def _pdf(self, x, a, b):
#                     gamma(a+b) * x**(a-1) * (1-x)**(b-1)
# beta.pdf(x, a, b) = ------------------------------------
#                              gamma(a)*gamma(b)
return np.exp(self._logpdf(x, a, b)) ```
Example #15
```def _pdf(self, x, df):
#                   x**(df-1) * exp(-x**2/2)
# chi.pdf(x, df) =  -------------------------
#                   2**(df/2-1) * gamma(df/2)
return np.exp(self._logpdf(x, df)) ```
Example #16
```def _stats(self, df):
mu = np.sqrt(2)*sc.gamma(df/2.0+0.5)/sc.gamma(df/2.0)
mu2 = df - mu*mu
g1 = (2*mu**3.0 + mu*(1-2*df))/np.asarray(np.power(mu2, 1.5))
g2 = 2*df*(1.0-df)-6*mu**4 + 4*mu**2 * (2*df-1)
g2 /= np.asarray(mu2**2.0)
return mu, mu2, g1, g2 ```
Example #17
```def _pdf(self, x, df):
# chi2.pdf(x, df) = 1 / (2*gamma(df/2)) * (x/2)**(df/2-1) * exp(-x/2)
return np.exp(self._logpdf(x, df)) ```
Example #18
```def _rvs(self, a):
sz, rndm = self._size, self._random_state
u = rndm.random_sample(size=sz)
gm = gamma.rvs(a, size=sz, random_state=rndm)
return gm * np.where(u >= 0.5, 1, -1) ```
Example #19
```def _pdf(self, x, a):
# dgamma.pdf(x, a) = 1 / (2*gamma(a)) * abs(x)**(a-1) * exp(-abs(x))
ax = abs(x)
return 1.0/(2*sc.gamma(a))*ax**(a-1.0) * np.exp(-ax) ```
Example #20
```def _munp(self, n, c):
return sc.gamma(1.0+n*1.0/c) ```
Example #21
```def _munp(self, n, c):
val = sc.gamma(1.0+n*1.0/c)
if int(n) % 2:
sgn = -1
else:
sgn = 1
return sgn * val ```
Example #22
```def _munp(self, n, c):
def __munp(n, c):
val = 0.0
k = np.arange(0, n + 1)
for ki, cnk in zip(k, sc.comb(n, k)):
val = val + cnk * (-1) ** ki / (1.0 - c * ki)
return np.where(c * n < 1, val * (-1.0 / c) ** n, np.inf)
return _lazywhere(c != 0, (c,),
lambda c: __munp(n, c),
sc.gamma(n + 1)) ```
Example #23
```def _stats(self, c):
g = lambda n: sc.gamma(n*c + 1)
g1 = g(1)
g2 = g(2)
g3 = g(3)
g4 = g(4)
g2mg12 = np.where(abs(c) < 1e-7, (c*np.pi)**2.0/6.0, g2-g1**2.0)
gam2k = np.where(abs(c) < 1e-7, np.pi**2.0/6.0,
sc.expm1(sc.gammaln(2.0*c+1.0)-2*sc.gammaln(c + 1.0))/c**2.0)
eps = 1e-14
gamk = np.where(abs(c) < eps, -_EULER, sc.expm1(sc.gammaln(c + 1))/c)

m = np.where(c < -1.0, np.nan, -gamk)
v = np.where(c < -0.5, np.nan, g1**2.0*gam2k)

# skewness
sk1 = _lazywhere(c >= -1./3,
(c, g1, g2, g3, g2mg12),
lambda c, g1, g2, g3, g2gm12:
np.sign(c)*(-g3 + (g2 + 2*g2mg12)*g1)/g2mg12**1.5,
fillvalue=np.nan)
sk = np.where(abs(c) <= eps**0.29, 12*np.sqrt(6)*_ZETA3/np.pi**3, sk1)

# kurtosis
ku1 = _lazywhere(c >= -1./4,
(g1, g2, g3, g4, g2mg12),
lambda g1, g2, g3, g4, g2mg12:
(g4 + (-4*g3 + 3*(g2 + g2mg12)*g1)*g1)/g2mg12**2,
fillvalue=np.nan)
ku = np.where(abs(c) <= (eps)**0.23, 12.0/5.0, ku1-3.0)
return m, v, sk, ku ```
Example #24
```def _munp(self, n, c):
k = np.arange(0, n+1)
vals = 1.0/c**n * np.sum(
sc.comb(n, k) * (-1)**k * sc.gamma(c*k + 1),
axis=0)
return np.where(c*n > -1, vals, np.inf) ```
Example #25
```def _pdf(self, x, a):
# gamma.pdf(x, a) = x**(a-1) * exp(-x) / gamma(a)
return np.exp(self._logpdf(x, a)) ```
Example #26
```def _fitstart(self, data):
# The skewness of the gamma distribution is `4 / np.sqrt(a)`.
# We invert that to estimate the shape `a` using the skewness
# of the data.  The formula is regularized with 1e-8 in the
# denominator to allow for degenerate data where the skewness
# is close to 0.
a = 4 / (1e-8 + _skew(data)**2)
return super(gamma_gen, self)._fitstart(data, args=(a,)) ```
Example #27
```def _pdf(self, x, a, c):
# gengamma.pdf(x, a, c) = abs(c) * x**(c*a-1) * exp(-x**c) / gamma(a)
return np.exp(self._logpdf(x, a, c)) ```
Example #28
```def _munp(self, n, a, c):
# Pochhammer symbol: sc.pocha,n) = gamma(a+n)/gamma(a)
return sc.poch(a, n*1.0/c) ```
Example #29
```def _munp(self, n):
if n == 1:
return 2*np.log(2)
if n == 2:
return np.pi*np.pi/3.0
if n == 3:
return 9*_ZETA3
if n == 4:
return 7*np.pi**4 / 15.0
return 2*(1-pow(2.0, 1-n))*sc.gamma(n+1)*sc.zeta(n, 1) ```
Example #30
```def _pdf(self, x, a):