我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用scipy.special.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 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 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))
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)
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)
