def rvs(df, gamma, n):
""" Generates random variables from a Skew t distribution
Parameters
----------
df : float
degrees of freedom parameter
gamma : float
skewness parameter
n : int or list
Number of simulations to perform; if list input, produces array
"""
if type(n) == list:
u = np.random.uniform(size=n[0]*n[1])
result = Skewt.ppf(q=u, df=df, gamma=gamma)
result = np.split(result,n[0])
return np.array(result)
else:
u = np.random.uniform(size=n)
if isinstance(df, np.ndarray) or isinstance(gamma, np.ndarray):
return np.array([Skewt.ppf(q=np.array([u[i]]), df=df[i], gamma=gamma[i])[0] for i in range(n)])
else:
return Skewt.ppf(q=u, df=df, gamma=gamma)
python类gamma()的实例源码
def logpdf_internal_prior(x, df, loc=0.0, scale=1.0, gamma = 1.0):
if x-loc < 0.0:
return np.log(2.0) - np.log(gamma + 1.0/gamma) + ss.t.logpdf(x=gamma*x, loc=loc*gamma,df=df, scale=scale)
else:
return np.log(2.0) - np.log(gamma + 1.0/gamma) + ss.t.logpdf(x=x/gamma, loc=loc/gamma,df=df, scale=scale)
def markov_blanket(y, mean, scale, shape, skewness):
""" Markov blanket for each likelihood term
Parameters
----------
y : np.ndarray
univariate time series
mean : np.ndarray
array of location parameters 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
----------
- Markov blanket 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
return Skewt.logpdf_internal(x=y, df=shape, loc=mean, gamma=skewness, scale=scale)
def neg_loglikelihood(y, mean, scale, shape, skewness):
""" Negative loglikelihood function
Parameters
----------
y : np.ndarray
univariate time series
mean : np.ndarray
array of location parameters 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
----------
- Negative loglikelihood 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
return -np.sum(Skewt.logpdf_internal(x=y, df=shape, loc=mean, gamma=skewness, scale=scale))
def reg_score_function(X, y, mean, scale, shape, skewness):
""" GAS Skew t Regression Update term using gradient only - native Python function
Parameters
----------
X : float
datapoint for the right hand side variable
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)*X)/(np.power(skewness*scale,2) + (np.power(y-mean,2)/shape))
else:
return ((shape+1)/shape)*((y-mean)*X)/(np.power(scale,2) + (np.power(skewness*(y-mean),2)/shape))
def second_order_score(y, mean, scale, shape, skewness):
""" GAS Skew t Update term potentially using second-order information - 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
----------
- Adjusted 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))
def cdf(x, df, loc=0.0, scale=1.0, gamma = 1.0):
"""
CDF function for Skew t distribution
"""
result = np.zeros(x.shape[0])
result[x<0] = 2.0/(np.power(gamma,2) + 1.0)*ss.t.cdf(gamma*(x[x-loc < 0]-loc[x-loc < 0])/scale, df=df)
result[x>=0] = 1.0/(np.power(gamma,2) + 1.0) + 2.0/((1.0/np.power(gamma,2)) + 1.0)*(ss.t.cdf((x[x-loc >= 0]-loc[x-loc >= 0])/(gamma*scale), df=df)-0.5)
return result
def ppf(q, df, loc=0.0, scale=1.0, gamma = 1.0):
"""
PPF function for Skew t distribution
"""
result = np.zeros(q.shape[0])
probzero = Skewt.cdf(x=np.zeros(1),loc=np.zeros(1),df=df,gamma=gamma)
result[q<probzero] = 1.0/gamma*ss.t.ppf(((np.power(gamma,2) + 1.0) * q[q<probzero])/2.0,df)
result[q>=probzero] = gamma*ss.t.ppf((1.0 + 1.0/np.power(gamma,2))/2.0*(q[q >= probzero] - probzero) + 0.5, df)
return result
# Optional Cythonized recursions below for GAS Skew t models
def plot_predict_is(self, h=5, fit_once=True, fit_method='MLE', **kwargs):
""" Plots forecasts with the estimated model against data (Simulated prediction with data)
Parameters
----------
h : int (default : 5)
How many steps to forecast
fit_once : boolean
(default: True) Fits only once before the in-sample prediction; if False, fits after every new datapoint
fit_method : string
Which method to fit the model with
Returns
----------
- Plot of the forecast against data
"""
import matplotlib.pyplot as plt
import seaborn as sns
figsize = kwargs.get('figsize',(10,7))
plt.figure(figsize=figsize)
date_index = self.index[-h:]
predictions = self.predict_is(h, fit_method=fit_method, fit_once=fit_once)
data = self.data[-h:]
t_params = self.transform_z()
loc = t_params[-2] + t_params[-1]*predictions.values.T[0] + (t_params[-4] - (1.0/t_params[-4]))*predictions.values.T[0]*(np.sqrt(t_params[-3])*sp.gamma((t_params[-3]-1.0)/2.0))/(np.sqrt(np.pi)*sp.gamma(t_params[-3]/2.0))
plt.plot(date_index, np.abs(data-loc), label='Data')
plt.plot(date_index, predictions, label='Predictions', c='black')
plt.title(self.data_name)
plt.legend(loc=2)
plt.show()
def logpdf(x, shape, loc=0.0, scale=1.0, skewness = 1.0):
m1 = (np.sqrt(shape)*sp.gamma((shape-1.0)/2.0))/(np.sqrt(np.pi)*sp.gamma(shape/2.0))
loc = loc + (skewness - (1.0/skewness))*scale*m1
result = np.zeros(x.shape[0])
result[x-loc<0] = np.log(2.0) - np.log(skewness + 1.0/skewness) + ss.t.logpdf(x=skewness*x[(x-loc) < 0], loc=loc[(x-loc) < 0]*skewness,df=shape, scale=scale[(x-loc) < 0])
result[x-loc>=0] = np.log(2.0) - np.log(skewness + 1.0/skewness) + ss.t.logpdf(x=x[(x-loc) >= 0]/skewness, loc=loc[(x-loc) >= 0]/skewness,df=shape, scale=scale[(x-loc) >= 0])
return result
def tv_variate_exp(df):
return (sqrt(df)*gamma((df-1.0)/2.0))/(sqrt(pi)*gamma(df/2.0))
def ball_volume(ndim, radius):
"""Return the volume of a ball of dimension `ndim` and radius `radius`."""
n = ndim
r = radius
return np.pi ** (n / 2) / special.gamma(n / 2 + 1) * r ** n
def mean(a, b):
""" Continuous mean. at most 1 step below discretized mean
`E[T ] <= E[Td] + 1` true for positive distributions.
"""
from scipy.special import gamma
return a * gamma(1.0 + 1.0 / b)
def mean(a, b):
"""Continuous mean. Theoretically at most 1 step below discretized mean
`E[T ] <= E[Td] + 1` true for positive distributions.
:param a: Alpha
:param b: Beta
:return: `a * gamma(1.0 + 1.0 / b)`
"""
from scipy.special import gamma
return a * gamma(1.0 + 1.0 / b)
def gammaDist(x,params):
'''Gamma distribution function
M,K = params, where K is average photon counts <x>,
M is the number of coherent modes,
In case of high intensity, the beam behavors like wave and
the probability density of photon, P(x), satify this gamma function.
'''
K,M = params
K = float(K)
M = float(M)
coeff = np.exp(M*np.log(M) + (M-1)*np.log(x) - gammaln(M) - M*np.log(K))
Gd = coeff*np.exp(-M*x/K)
return Gd
def gamma_dist(bin_values, K, M):
"""
Gamma distribution function
Parameters
----------
bin_values : array
scattering intensities
K : int
average number of photons
M : int
number of coherent modes
Returns
-------
gamma_dist : array
Gamma distribution
Notes
-----
These implementations are based on the references under
nbinom_distribution() function Notes
: math ::
P(K) =(\frac{M}{<K>})^M \frac{K^(M-1)}{\Gamma(M)}\exp(-M\frac{K}{<K>})
"""
#gamma_dist = (stats.gamma(M, 0., K/M)).pdf(bin_values)
x= bin_values
coeff = np.exp(M*np.log(M) + (M-1)*np.log(x) - gammaln(M) - M*np.log(K))
gamma_dist = coeff*np.exp(-M*x/K)
return gamma_dist
def poisson(x,K):
'''Poisson distribution function.
K is average photon counts
In case of low intensity, the beam behavors like particle and
the probability density of photon, P(x), satify this poisson function.
'''
K = float(K)
Pk = np.exp(-K)*power(K,x)/gamma(x+1)
return Pk
def Schultz_Zimm(x,u,sigma):
'''http://sasfit.ingobressler.net/manual/Schultz-Zimm
See also The size distribution of ‘gold standard’ nanoparticles
Anal Bioanal Chem (2009) 395:1651–1660
DOI 10.1007/s00216-009-3049-5
'''
k = 1.0/ (sigma)**2
return 1.0/u * (x/u)**(k-1) * k**k*np.exp( -k*x/u)/gamma(k)
def gamma_dist(bin_values, K, M):
"""
Gamma distribution function
Parameters
----------
bin_values : array
scattering intensities
K : int
average number of photons
M : int
number of coherent modes
Returns
-------
gamma_dist : array
Gamma distribution
Notes
-----
These implementations are based on the references under
nbinom_distribution() function Notes
: math ::
P(K) =(\frac{M}{<K>})^M \frac{K^(M-1)}{\Gamma(M)}\exp(-M\frac{K}{<K>})
"""
#gamma_dist = (stats.gamma(M, 0., K/M)).pdf(bin_values)
x= bin_values
coeff = np.exp(M*np.log(M) + (M-1)*np.log(x) - gammaln(M) - M*np.log(K))
gamma_dist = coeff*np.exp(-M*x/K)
return gamma_dist
def poisson(x,K):
'''Poisson distribution function.
K is average photon counts
In case of low intensity, the beam behavors like particle and
the probability density of photon, P(x), satify this poisson function.
'''
K = float(K)
Pk = np.exp(-K)*power(K,x)/gamma(x+1)
return Pk
