def chin(self):
"""approximation of the expected length.
The exact value could be computed by::
from scipy.special import gamma
return 2**0.5 * gamma((self.dimension+1) / 2) / gamma(self.dimension / 2)
The approximation obeys ``chin < chin_hat < (1 + 5e-5) * chin``.
"""
values = {1: 0.7978845608028656, 2: 1.2533141373156,
3: 1.59576912160574, 4: 1.87997120597326}
try:
val = values[self.dimension]
except KeyError:
# for dim > 4 we have chin < chin_hat < (1 + 5e-5) * chin
N = self.dimension
val = N**0.5 * (1 - 1. / (4 * N) + 1. / (26 * N**2)) # was: 21
return val
python类gamma()的实例源码
def chin(self):
"""approximation of the expected length.
The exact value could be computed by::
from scipy.special import gamma
return 2**0.5 * gamma((self.dimension+1) / 2) / gamma(self.dimension / 2)
The approximation obeys ``chin < chin_hat < (1 + 5e-5) * chin``.
"""
values = {1: 0.7978845608028656, 2: 1.2533141373156,
3: 1.59576912160574, 4: 1.87997120597326}
try:
val = values[self.dimension]
except KeyError:
# for dim > 4 we have chin < chin_hat < (1 + 5e-5) * chin
N = self.dimension
val = N**0.5 * (1 - 1. / (4 * N) + 1. / (26 * N**2)) # was: 21
return val
def chin(self):
"""approximation of the expected length.
The exact value could be computed by::
from scipy.special import gamma
return 2**0.5 * gamma((self.dimension+1) / 2) / gamma(self.dimension / 2)
The approximation obeys ``chin < chin_hat < (1 + 5e-5) * chin``.
"""
values = {1: 0.7978845608028656, 2: 1.2533141373156,
3: 1.59576912160574, 4: 1.87997120597326}
try:
val = values[self.dimension]
except KeyError:
# for dim > 4 we have chin < chin_hat < (1 + 5e-5) * chin
N = self.dimension
val = N**0.5 * (1 - 1. / (4 * N) + 1. / (26 * N**2)) # was: 21
return val
def kolmogorov(self, r0):
self.covariance = np.zeros((self.nZernike,self.nZernike))
for i in range(self.nZernike):
ni, mi = wf.nollIndices(i+self.zero)
for j in range(self.nZernike):
nj, mj = wf.nollIndices(j+self.zero)
if (even(i - j)):
if (mi == mj):
phase = (-1.0)**(0.5*(ni+nj-2*mi))
t1 = np.sqrt((ni+1)*(nj+1)) * np.pi**(8.0/3.0) * 0.0072 * (self.telescopeD / r0)**(5.0/3.0)
t2 = sp.gamma(14./3.0) * sp.gamma(0.5*(ni+nj-5.0/3.0))
t3 = sp.gamma(0.5*(ni-nj+17.0/3.0)) * sp.gamma(0.5*(nj-ni+17.0/3.0)) * sp.gamma(0.5*(ni+nj+23.0/3.0))
self.covariance[i,j] = phase * t1 * t2 / t3
self.coeff = np.random.multivariate_normal(np.zeros(self.nZernike), self.covariance, size=(self.nSteps, self.nImages))
def logpdf(self, mu):
"""
Log PDF for Skew t prior
Parameters
----------
mu : float
Latent variable for which the prior is being formed over
Returns
----------
- log(p(mu))
"""
if self.transform is not None:
mu = self.transform(mu)
return self.logpdf_internal_prior(mu, df=self.df0, loc=self.loc0, scale=self.scale0, gamma=self.gamma0)
def setup():
""" Returns the attributes of this family
Notes
----------
- scale notes whether family has a variance parameter (sigma)
- shape notes whether family has a tail thickness parameter (nu)
- skewness notes whether family has a skewness parameter (gamma)
- mean_transform is a function which transforms the location parameter
- cythonized notes whether the family has cythonized routines
Returns
----------
- model name, link function, scale, shape, skewness, mean_transform, cythonized
"""
name = "Skewt"
link = np.array
scale = True
shape = True
skewness = True
mean_transform = np.array
cythonized = True
return name, link, scale, shape, skewness, mean_transform, cythonized
def pdf(self, mu):
"""
PDF for Skew t prior
Parameters
----------
mu : float
Latent variable for which the prior is being formed over
Returns
----------
- p(mu)
"""
if self.transform is not None:
mu = self.transform(mu)
return self.pdf_internal(mu, df=self.df0, loc=self.loc0, scale=self.scale0, gamma=self.gamma0)
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
def poisson_dist(bin_values, K):
"""
Poisson Distribution
Parameters
---------
K : int
average counts of photons
bin_values : array
scattering bin values
Returns
-------
poisson_dist : array
Poisson Distribution
Notes
-----
These implementations are based on the references under
nbinom_distribution() function Notes
:math ::
P(K) = \frac{<K>^K}{K!}\exp(-<K>)
"""
#poisson_dist = stats.poisson.pmf(K, bin_values)
K = float(K)
poisson_dist = np.exp(-K) * np.power(K, bin_values)/gamma(bin_values + 1)
return poisson_dist
def poisson_dist(bin_values, K):
"""
Poisson Distribution
Parameters
---------
K : int
average counts of photons
bin_values : array
scattering bin values
Returns
-------
poisson_dist : array
Poisson Distribution
Notes
-----
These implementations are based on the references under
nbinom_distribution() function Notes
:math ::
P(K) = \frac{<K>^K}{K!}\exp(-<K>)
"""
#poisson_dist = stats.poisson.pmf(K, bin_values)
K = float(K)
poisson_dist = np.exp(-K) * np.power(K, bin_values)/gamma(bin_values + 1)
return poisson_dist
def estimate_bayes_factor(traces, logp, r=0.05):
"""
Esitmate Odds ratios in a random subsample of the chains in MCMC
AstroML (see Eqn 5.127, pg 237)
"""
ndim, nsteps = traces.shape # [ndim,number of steps in chain]
# compute volume of a n-dimensional (ndim) sphere of radius r
Vr = np.pi ** (0.5 * ndim) / gamma(0.5 * ndim + 1) * (r ** ndim)
# use neighbor count within r as a density estimator
bt = BallTree(traces.T)
count = bt.query_radius(traces.T, r=r, count_only=True)
# BF = N*p/rho
bf = logp + np.log(nsteps) + np.log(Vr) - np.log(count) #log10(bf)
p25, p50, p75 = np.percentile(bf, [25, 50, 75])
return p50, 0.7413 * (p75 - p25)
########################################################################
def g(self, x):
"""Fallout Deposition Distribution Function.
Throughout the growth and transport of the radioactive cloud there is a
continual fall of particles back to the ground. WSEG states that there must be some
function "g(t)" which describes the fractional rate of activity arrival on the ground
everywhere at some time t. The integral of this function, G(t), represents the
fraction of activity down at time t. This g(t) function will be independent of the
horizontal activity distribution and therefore independent of the growth of a with
time. On the other hand g(t) will be dependent on the initial vertical distribution
and the activity/size distribution which determines particle fall rate. This
arbitrary choice of g(t) is based on Rand calculations which assume an activity/size
distribution given by activity_size_distribution(). These calculations are neither
shown nor referenced in the original 1959 WSEG model. If the activity/size
distribution for a given set of initial conditions is different
than that given by activity_size_distribution(), the form of g(t) should change.
This is not possible under the WSEG model where the function g(t) is fixed. The only
possible compensation for various activity/size distributions results because T_c
varies with yield."""
return np.exp(-(np.abs(x) / self.L)**self.n) / (self.L * gamma(1 + 1 / self.n))
def K2(bn, data, ed=None):
"""
K2 is bayesian posterior probability of structure given the data,
where N'ijk = 1.
"""
counts_dict = mle_fast(bn, data, counts=True, np=True)
a_ijk = []
k2 = 1
for rv, value in counts_dict.items():
nijk = value['cpt']
nijk_prime = 1
k2 *= gamma(nijk+nijk_prime)/gamma(nijk_prime)
nij_prime = nijk_prime*(len(cpt)/bn.card(rv))
nij = np.mean(nijk.reshape(-1, bn.card(rv)), axis=1) # sum along parents
k2 *= gamma(nij_prime) / gamma(nij+nij_prime)
return k2
def compute_by_noise_pow(self, signal, n_pow):
s_spec = np.fft.fftpack.fft(signal * self._window)
s_amp = np.absolute(s_spec)
s_phase = np.angle(s_spec)
gamma = self._calc_aposteriori_snr(s_amp, n_pow)
xi = self._calc_apriori_snr(gamma)
self._prevGamma = gamma
nu = gamma * xi / (1.0 + xi)
self._G = (self._gamma15 * np.sqrt(nu) / gamma) * np.exp(-nu / 2.0) *\
((1.0 + nu) * spc.i0(nu / 2.0) + nu * spc.i1(nu / 2.0))
idx = np.less(s_amp ** 2.0, n_pow)
self._G[idx] = self._constant
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = xi[idx] / (xi[idx] + 1.0)
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = self._constant
self._G = np.maximum(self._G, 0.0)
amp = self._G * s_amp
amp = np.maximum(amp, 0.0)
amp2 = self._ratio * amp + (1.0 - self._ratio) * s_amp
self._prevAmp = amp
spec = amp2 * np.exp(s_phase * 1j)
return np.real(np.fft.fftpack.ifft(spec))
def compute_by_noise_pow(self, signal, n_pow):
s_spec = np.fft.fftpack.fft(signal * self._window)
s_amp = np.absolute(s_spec)
s_phase = np.angle(s_spec)
gamma = self._calc_aposteriori_snr(s_amp, n_pow)
# xi = self._calc_apriori_snr2(gamma,n_pow)
xi = self._calc_apriori_snr(gamma)
self._prevGamma = gamma
u = 0.5 - self._mu / (4.0 * np.sqrt(gamma * xi))
self._G = u + np.sqrt(u ** 2.0 + self._tau / (gamma * 2.0))
idx = np.less(s_amp ** 2.0, n_pow)
self._G[idx] = self._constant
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = xi[idx] / (xi[idx] + 1.0)
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = self._constant
self._G = np.maximum(self._G, 0.0)
amp = self._G * s_amp
amp = np.maximum(amp, 0.0)
amp2 = self._ratio * amp + (1.0 - self._ratio) * s_amp
self._prevAmp = amp
spec = amp2 * np.exp(s_phase * 1j)
return np.real(np.fft.fftpack.ifft(spec))
def volume_of_the_unit_ball(d):
""" Volume of the d-dimensional unit ball.
Parameters
----------
d : int
dimension.
Returns
-------
vol : float
volume.
"""
vol = pi**(d/2) / gamma(d/2+1) # = 2 * pi^(d/2) / ( d*gamma(d/2) )
return vol
def K(self, X, Xstar):
"""
Computes covariance function values over `X` and `Xstar`.
Parameters
----------
X: np.ndarray, shape=((n, nfeatures))
Instances
Xstar: np.ndarray, shape=((n, nfeatures))
Instances
Returns
-------
np.ndarray
Computed covariance matrix.
"""
r = l2norm_(X, Xstar)
bessel = kv(self.v, np.sqrt(2 * self.v) * r / self.l)
f = 2 ** (1 - self.v) / gamma(self.v) * (np.sqrt(2 * self.v) * r / self.l) ** self.v
res = f * bessel
res[np.isnan(res)] = 1
res = self.sigmaf * res + self.sigman * kronDelta(X, Xstar)
return (res)
def K(self, X, Xstar):
"""
Computes covariance function values over `X` and `Xstar`.
Parameters
----------
X: np.ndarray, shape=((n, nfeatures))
Instances
Xstar: np.ndarray, shape=((n, nfeatures))
Instances
Returns
-------
np.ndarray
Computed covariance matrix.
"""
r = l2norm_(X, Xstar)
return self.sigmaf * (np.exp(-(r / self.l) ** self.gamma)) + \
self.sigman * kronDelta(X, Xstar)
def estimate_MSE_vs_alpha(transform, Ms, alphas, K):
# Without loss of generality
Z = 1
tZ = transform(Z)
# Estimate MSEs by constructing estimators K times
MSEs = np.empty((len(Ms), len(alphas)))
MSEs_stdev = np.empty((len(Ms), len(alphas)))
for Mi, M in enumerate(Ms):
# Compute means (K x alphas) in a loop, as otherwise
# this runs out of memory with K = 100,000.
means = np.empty((K, len(alphas)))
for ai, alpha in enumerate(alphas):
Ws = np.power(np.random.exponential(1.0, size=(K, M)), alpha) # (K, M)
means[:, ai] = np.mean(Ws, axis=1)
print_progress('M = %d: done %.0f%%' % (M, 100.0 * (ai+1) / len(alphas)))
print('')
g = np.power(gamma(1.0 + alphas), 1.0 / alphas) # (alphas)
tZ_hats = transform(g * np.power(means, -1.0/alphas)) # (K, alphas)
SEs = (tZ_hats - tZ) ** 2 # (K, alphas)
MSEs[Mi] = np.mean(SEs, axis=0) # (alphas)
MSEs_stdev[Mi] = np.std(SEs, axis=0) / np.sqrt(K) # (alphas)
return MSEs, MSEs_stdev
def simple_multivariate_t_distribution(self, x, mu, sigma, df, d):
'''
Multivariate t-student 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)
Sigma = scale matrix (dxd numpy array)
df = degrees of freedom
d: dimension
'''
num = special.gamma((df + d)/2)
xSigma = np.dot((x - mu), np.linalg.inv(sigma))
xSigma = np.array([xSigma[i, i] for i in range(self.K)])
denom = special.gamma(df / 2) * np.power(df * np.pi, d / 2.0)\
* np.power(np.linalg.det(sigma), 1 / 2.0)\
* np.power(1 + (1. / df)\
* np.sum(xSigma * (x - mu), axis=1), (d + df) / 2)
result = num / denom
return result
def matern(h, a, C0, s, b=0):
"""
The Matérn model.
For Matérn function see:
Minasny, B., & McBratney, A. B. (2005). The Matérn function as a general model for soil variograms.
Geoderma, 128(3–4 SPEC. ISS.), 192–207. http://doi.org/10.1016/j.geoderma.2005.04.003.
:param h: lag
:param a: range
:param C0: sill
:param s: smoothness parameter
:param b: nugget
:return:
"""
# prepare parameters
r = a
C0 -= b
return b + C0 * (1 - ( (1 / (np.power(2, s - 1) * special.gamma(s))) * np.power(h / r, s) * special.kv(s, h / r) ))
# --- Adaptions using no nugget effect --- #
def spherical_beta_logp(x, lon_lat, alpha):
if x[1] < -90. or x[1] > 90.:
raise ZeroProbability
return -np.inf
if alpha == 1.0:
return np.log(1. / 4. / np.pi)
mu = np.array([np.cos(lon_lat[1] * d2r) * np.cos(lon_lat[0] * d2r),
np.cos(lon_lat[1] * d2r) * np.sin(lon_lat[0] * d2r),
np.sin(lon_lat[1] * d2r)])
test_point = np.transpose(np.array([np.cos(x[1] * d2r) * np.cos(x[0] * d2r),
np.cos(x[1] * d2r) *
np.sin(x[0] * d2r),
np.sin(x[1] * d2r)]))
thetas = np.arccos(np.dot(test_point, mu))
normalization = sp.gamma(alpha + 0.5) / \
sp.gamma(alpha) / np.sqrt(np.pi) / np.pi / 2.
logp_elem = np.log(np.sin(thetas)) * (2. * alpha - 2) + \
np.log(normalization)
logp = logp_elem.sum()
return logp
def multivariatet(mu, Sigma, N, M, rng):
"""Return a sample (or samples) from the multivariate t distribution.
This function is adopted from
http://kennychowdhary.me/2013/03/python-code-to-generate-samples-from-multivariate-t/.
:param mu: Mean.
:type mu: ndarray, shape=(n_dim,), dtype=float
:param Sigma: Scaling matrix.
:type Sigma: ndarray, shape=(n_dim, n_dim), dtype=float
:param float N: Degrees of freedom.
:param int M: Number of samples to produce.
:param np.random.RandomState rng: Random number generator.
:return: M samples of (n_dim)-dimensional multivariate t distribution.
:rtype: ndarray, shape=(n_samples, n_dim), dtype=float
"""
d = len(Sigma)
g = np.tile(rng.gamma(N/2., 2./N, M), (d, 1)).T
Z = rng.multivariate_normal(np.zeros(d), Sigma, M)
return mu + Z / np.sqrt(g)
def nu2phr(nu):
phr = np.sqrt(2.0 / nu) * spesh.gamma((nu + 1.0) / 2.0) / \
spesh.gamma(nu / 2.0)
phr = sps.t.pdf(0.0, nu) / sps.norm.pdf(0.0)
return phr
# read in data
def _sample_alpha(self, n=1):
eps = 1e-5
loglikelihood = lambda alpha: (alpha - 1) * np.log(self.a_0+eps) + alpha * self.c_0 * np.log(self.beta+eps) \
- self.b_0 * np.log(special.gamma(alpha))
likelihood = lambda alpha: np.exp(loglikelihood(alpha))
stop = 15
alpha_space = np.linspace(0, stop, num=old_div(stop, 0.001))
alpha_dist = likelihood(alpha_space)
alpha_dist = old_div(alpha_dist, alpha_dist.sum())
return np.random.choice(a=alpha_space, p=alpha_dist, size=n)
def factorial(n):
"""
Compute the factorial of a number
Args:
n (real): input number
Returns:
real: n!
"""
return gamma(n+1)
def plot_fit(self, **kwargs):
"""
Plots the fit of the model against the data
"""
import matplotlib.pyplot as plt
import seaborn as sns
figsize = kwargs.get('figsize',(10,7))
plt.figure(figsize=figsize)
date_index = self.index[max(self.ar, self.ma):self.data.shape[0]]
mu, Y = self._model(self.latent_variables.get_z_values())
# Catch specific family properties (imply different link functions/moments)
if self.model_name2 == "Exponential":
values_to_plot = 1.0/self.link(mu)
elif self.model_name2 == "Skewt":
t_params = self.transform_z()
model_scale, model_shape, model_skewness = self._get_scale_and_shape(t_params)
m1 = (np.sqrt(model_shape)*sp.gamma((model_shape-1.0)/2.0))/(np.sqrt(np.pi)*sp.gamma(model_shape/2.0))
additional_loc = (model_skewness - (1.0/model_skewness))*model_scale*m1
values_to_plot = mu + additional_loc
else:
values_to_plot = self.link(mu)
plt.plot(date_index, Y, label='Data')
plt.plot(date_index, values_to_plot, label='ARIMA model', c='black')
plt.title(self.data_name)
plt.legend(loc=2)
plt.show()
def plot_fit(self, **kwargs):
"""
Plots the fit of the model against the data
"""
import matplotlib.pyplot as plt
import seaborn as sns
figsize = kwargs.get('figsize',(10,7))
plt.figure(figsize=figsize)
date_index = self.index[max(self.ar, self.ma):self.data_length]
mu, Y = self._model(self.latent_variables.get_z_values())
# Catch specific family properties (imply different link functions/moments)
if self.model_name2 == "Exponential":
values_to_plot = 1.0/self.link(mu)
elif self.model_name2 == "Skewt":
t_params = self.transform_z()
model_scale, model_shape, model_skewness = self._get_scale_and_shape(t_params)
m1 = (np.sqrt(model_shape)*sp.gamma((model_shape-1.0)/2.0))/(np.sqrt(np.pi)*sp.gamma(model_shape/2.0))
additional_loc = (model_skewness - (1.0/model_skewness))*model_scale*m1
values_to_plot = mu + additional_loc
else:
values_to_plot = self.link(mu)
plt.plot(date_index, Y, label='Data')
plt.plot(date_index, values_to_plot, label='ARIMA model', c='black')
plt.title(self.data_name)
plt.legend(loc=2)
plt.show()
def __init__(self, loc=0.0, scale=1.0, df=8.0, gamma=1.0, transform=None, **kwargs):
"""
Parameters
----------
loc : float
Location parameter for the Skew t distribution
scale : float
Scale parameter for the Skew t distribution
df : float
Degrees of freedom parameter for the Skew t distribution
gamma : float
Skewness parameter (1.0 is skewed; under 1.0, -ve skewed; over 1.0, +ve skewed)
transform : str
Whether to apply a transformation to the location variable - e.g. 'exp' or 'logit'
"""
super(Skewt, self).__init__(transform)
self.loc0 = loc
self.scale0 = scale
self.df0 = df
self.gamma0 = gamma
self.covariance_prior = False
self.gradient_only = kwargs.get('gradient_only', False) # used for GAS t models
if self.gradient_only is True:
self.score_function = self.first_order_score
else:
self.score_function = self.second_order_score
def first_order_score(y, mean, scale, shape, skewness):
""" GAS Skew t Update term using gradient only - native Python function
Parameters
----------
y : float
datapoint for the time series
mean : float
location parameter for the Skew t distribution
scale : float
scale parameter for the Skew t distribution
shape : float
tail thickness parameter for the Skew t distribution
skewness : float
skewness parameter for the Skew t distribution
Returns
----------
- Score of the Skew t family
"""
m1 = (np.sqrt(shape)*sp.gamma((shape-1.0)/2.0))/(np.sqrt(np.pi)*sp.gamma(shape/2.0))
mean = mean + (skewness - (1.0/skewness))*scale*m1
if (y-mean)>=0:
return ((shape+1)/shape)*(y-mean)/(np.power(skewness*scale,2) + (np.power(y-mean,2)/shape))
else:
return ((shape+1)/shape)*(y-mean)/(np.power(scale,2) + (np.power(skewness*(y-mean),2)/shape))
