Python源码示例：scipy.stats.norm.ppf()
示例1
def _t_value(self):
r"""
Returns the critical t-statistic given the input alpha-level (defaults to 0.05).
Returns
-------
tval : float
The critical t-value for using in computing the Least Significant Difference.
Notes
-----
Scipy's :code:`t.ppf` method is used to compute the critical t-value.
"""
tval = t.ppf(1 - self.alpha / 2, self.n - self.k)
return tval
示例2
def test_calc_bias_correction_bca(self):
# There are 100 bootstrap replicates, already in ascending order for
# each column. If we take row 51 to be the mle, then 50% of the
# replicates are less than the mle, and we should have bias = 0.
expected_result = np.zeros(self.mle_params.size)
# Alias the function to be tested.
func = bc.calc_bias_correction_bca
# Perform the desired test
func_result = func(self.bootstrap_replicates, self.mle_params)
self.assertIsInstance(func_result, np.ndarray)
self.assertEqual(func_result.shape, expected_result.shape)
npt.assert_allclose(func_result, expected_result)
# Create a fake mle that should be higher than 95% of the results
fake_mle = self.bootstrap_replicates[95]
expected_result_2 = norm.ppf(0.95) * np.ones(self.mle_params.size)
func_result_2 = func(self.bootstrap_replicates, fake_mle)
self.assertIsInstance(func_result_2, np.ndarray)
self.assertEqual(func_result_2.shape, expected_result_2.shape)
npt.assert_allclose(func_result_2, expected_result_2)
return None
示例3
def update(self, model):
self.model = model
self.sn2 = self.model.get_noise()
# Sample representer points
self.sampling_acquisition.update(model)
self.sample_representer_points()
# Omega values which are needed for the innovations
# by sampling from a uniform grid
self.W = norm.ppf(np.linspace(1. / (self.Np + 1),
1 - 1. / (self.Np + 1),
self.Np))[np.newaxis, :]
# Compute current posterior belief at the representer points
self.Mb, self.Vb = self.model.predict(self.zb, full_cov=True)
self.pmin = mc_part.joint_pmin(self.Mb, self.Vb, self.Nf)
self.logP = np.log(self.pmin)
示例4
def testSecurityNormInvValueHolder(self):
mm1 = SecurityNormInvValueHolder('open')
mm2 = SecurityNormInvValueHolder('open', fullAcc=True)
for i in range(len(self.aapl['close'])):
data = dict(aapl=dict(open=norm.cdf(self.aapl['open'][i])),
ibm=dict(open=norm.cdf(self.ibm['open'][i])))
mm1.push(data)
mm2.push(data)
value1 = mm1.value
value2 = mm2.value
for name in value1.index():
expected = norm.ppf(data[name]['open'])
calculated = value1[name]
self.assertAlmostEqual(expected, calculated, 6, 'at index {0}\n'
'expected: {1:.12f}\n'
'calculat: {2:.12f}'
.format(i, expected, calculated))
calculated = value2[name]
self.assertAlmostEqual(expected, calculated, 12, 'at index {0}\n'
'expected: {1:.12f}\n'
'calculat: {2:.12f}'
.format(i, expected, calculated))
示例5
def liptak(pvalues):
r"""
Apply Liptak's combining function
.. math:: \sum_i \Phi^{-1}(1-p_i)
where $\Phi^{-1}$ is the inverse CDF of the standard normal distribution.
Parameters
----------
pvalues : array_like
Array of p-values to combine
Returns
-------
float
Liptak's combined test statistic
"""
return np.sum(norm.ppf(1 - pvalues))
示例6
def predictions(self, image, forward_batch_size=32):
from scipy.stats import norm
image, _ = self._process_input(image)
image_batch = np.vstack([[image]] * self._iterations)
noise = np.random.normal(scale=self._std, size=image_batch.shape).astype(np.float32)
image_batch += noise
predictions = self._model.batch_predictions(image_batch)
logits = np.argmax(predictions, axis=1)
one_hot = np.zeros([self._iterations, self._num_classes])
logits_one_hot = np.eye(self._num_classes)[logits]
one_hot += logits_one_hot
one_hot = np.sum(one_hot, axis=0)
ranks = sorted(one_hot / np.sum(one_hot))[::-1]
qi = ranks[0] - 1e-9
qj = ranks[1] + 1e-9
bound = self._std / 2. * (norm.ppf(qi) - norm.ppf(qj))
return np.argmax(one_hot), bound
示例7
def std_mad(x):
"""Robust estimation of the standard deviation, based on the Corrected Median
Absolute Deviation (MAD) of x.
This computes the MAD of x, and applies the Gaussian distribution
correction, making it a consistent estimator of the standard-deviation
(when the sample looks Gaussian with outliers).
Parameters
----------
x : `np.ndarray`
Input vector
Returns
-------
output : `float`
A robust estimation of the standard deviation
"""
from scipy.stats import norm
correction = 1 / norm.ppf(3 / 4)
return correction * np.median(np.abs(x - np.median(x)))
示例8
def bias_corrected_ci(estimate, samples, conf=95):
"""
Return the bias-corrected bootstrap confidence interval for an estimate
:param estimate: Numerical estimate in the original sample
:param samples: Nx1 array of bootstrapped estimates
:param conf: Level of the desired confidence interval
:return: Bias-corrected bootstrapped LLCI and ULCI for the estimate.
"""
# noinspection PyUnresolvedReferences
ptilde = ((samples < estimate) * 1).mean()
Z = norm.ppf(ptilde)
Zci = z_score(conf)
Zlow, Zhigh = -Zci + 2 * Z, Zci + 2 * Z
plow, phigh = norm._cdf(Zlow), norm._cdf(Zhigh)
llci = np.percentile(samples, plow * 100, interpolation="lower")
ulci = np.percentile(samples, phigh * 100, interpolation="higher")
return llci, ulci
示例9
def _fisher_confint(self, alpha: float, observed: bool = False) -> List[float]:
"""Compute the Fisher information confidence interval for the MLE of the previous run.
Args:
alpha: Specifies the (1 - alpha) confidence level (0 < alpha < 1).
observed: If True, the observed Fisher information is used to construct the
confidence interval, otherwise the expected Fisher information.
Returns:
The Fisher information confidence interval.
"""
shots = self._ret['shots']
mle = self._ret['ml_value']
# approximate the standard deviation of the MLE and construct the confidence interval
std = np.sqrt(shots * self._compute_fisher_information(observed))
ci = mle + norm.ppf(1 - alpha / 2) / std * np.array([-1, 1])
# transform the confidence interval from [0, 1] to the target interval
return [self.a_factory.value_to_estimation(bound) for bound in ci]
示例10
def _init_model(self, C, eta):
"""
Initialize model.
"""
logger.info("init model starts")
self.model["mu"] = defaultdict() # model parameter mean
self.model["S"] = defaultdict() # model parameter covariance
self.model["C"] = C # PA parameter
self.model["eta"] = eta # confidence parameter
self.model["phi"] = norm.ppf(norm.cdf(eta)) # inverse of cdf(eta)
self.model["phi_2"] = np.power(self.model["phi"], 2)
self.model["psi"] = 1 + self.model["phi_2"] / 2
self.model["zeta"] = 1 + self.model["phi_2"]
logger.info("init model finished")
pass
示例11
def prewarp(self, xx):
"""Extra work needed to get variables into the Gaussian space
representation."""
xxw = {}
for arg_name, vv in xx.items():
assert np.isscalar(vv)
space = self.space[arg_name]
if space is not None:
# Warp so we think it is apriori uniform in [a, b]
vv = space.warp(vv)
assert vv.size == 1
# Now make uniform on [0, 1], also unpack warped to scalar
(lb, ub), = space.get_bounds()
vv = linear_rescale(vv.item(), lb, ub, 0, 1)
# Now make std Gaussian apriori
vv = norm.ppf(vv)
assert np.isscalar(vv)
xxw[arg_name] = vv
return xxw
示例12
def prewarp(self, xx):
"""Extra work needed to get variables into the Gaussian space
representation."""
xxw = {}
for arg_name, vv in xx.items():
assert np.isscalar(vv)
space = self.space[arg_name]
if space is not None:
# Warp so we think it is apriori uniform in [a, b]
vv = space.warp(vv)
assert vv.size == 1
# Now make uniform on [0, 1], also unpack warped to scalar
(lb, ub), = space.get_bounds()
vv = linear_rescale(vv.item(), lb, ub, 0, 1)
# Now make std Gaussian apriori
vv = norm.ppf(vv)
assert np.isscalar(vv)
xxw[arg_name] = vv
return xxw
示例13
def _confidence_interval_by_alpha(cls, p_hat, n, alpha, method='wald'):
"""Compute confidence interval for estimate of Bernoulli parameter p.
Args:
p_hat: maximum likelihood estimate of p
n: samples observed
alpha: the probability that the true p falls outside the CI
Returns:
left, right
"""
prob = 1 - 0.5 * alpha
z = norm.ppf(prob)
compute_ci = cls._confidence_interval_by_z_wald if method == 'wald' else cls._confidence_interval_by_z_wilson
return compute_ci(p_hat, n, z)
示例14
def test_conditional_value_at_risk_mc(self):
for mu, sigma, alpha in [(1, 1, 0.05), (0.4, 0.1, 0.02), (0.1, 2, 0.01)]:
# prepare estimator dummy
mu1 = np.array([mu])
sigma1 = np.identity(n=1) * sigma
est = GaussianDummy(mean=mu1, cov=sigma1**2, ndim_x=1, ndim_y=1, has_pdf=True)
est.fit(None, None)
CVaR_true = mu - sigma/alpha * norm.pdf(norm.ppf(alpha))
CVaR_est = est.conditional_value_at_risk(x_cond=np.array([[0],[1]]), alpha=alpha)
print("CVaR True (%.2f, %.2f):"%(mu, sigma), CVaR_true)
print("CVaR_est (%.2f, %.2f):"%(mu, sigma), CVaR_est)
print("VaR (%.2f, %.2f):"%(mu, sigma), est.value_at_risk(x_cond=np.array([[0],[1]]), alpha=alpha))
self.assertAlmostEqual(CVaR_est[0], CVaR_true, places=2)
self.assertAlmostEqual(CVaR_est[1], CVaR_true, places=2)
示例15
def value_at_risk(self, x_cond, alpha=0.01, **kwargs):
""" Computes the Value-at-Risk (VaR) of the fitted distribution. Only if ndim_y = 1
Args:
x_cond: different x values to condition on - numpy array of shape (n_values, ndim_x)
alpha: quantile percentage of the distribution
Returns:
VaR values for each x to condition on - numpy array of shape (n_values)
"""
assert self.ndim_y == 1, "Value at Risk can only be computed when ndim_y = 1"
assert x_cond.ndim == 2
VaR = norm.ppf(alpha, loc=x_cond, scale=self._std(x_cond))[:,0]
assert VaR.shape == (x_cond.shape[0],)
return VaR
示例16
def value_at_risk(self, x_cond, alpha=0.01, **kwargs):
""" Computes the Value-at-Risk (VaR) of the fitted distribution. Only if ndim_y = 1
Args:
x_cond: different x values to condition on - numpy array of shape (n_values, ndim_x)
alpha: quantile percentage of the distribution
Returns:
VaR values for each x to condition on - numpy array of shape (n_values)
"""
assert self.ndim_y == 1, "Value at Risk can only be computed when ndim_y = 1"
assert x_cond.ndim == 2
VaR = norm.ppf(alpha, loc=self._mean(x_cond), scale=self._std(x_cond))[:,0]
assert VaR.shape == (x_cond.shape[0],)
return VaR
示例17
def conditional_value_at_risk(self, x_cond, alpha=0.01, **kwargs):
""" Computes the Conditional Value-at-Risk (CVaR) / Expected Shortfall of the fitted distribution. Only if ndim_y = 1
Args:
x_cond: different x values to condition on - numpy array of shape (n_values, ndim_x)
alpha: quantile percentage of the distribution
n_samples: number of samples for monte carlo model_fitting
Returns:
CVaR values for each x to condition on - numpy array of shape (n_values)
"""
assert self.ndim_y == 1, "Value at Risk can only be computed when ndim_y = 1"
x_cond = self._handle_input_dimensionality(x_cond)
assert x_cond.ndim == 2
mean = self._mean(x_cond)
sigma = self._std(x_cond)
CVaR = (mean - sigma * (1/alpha) * norm.pdf(norm.ppf(alpha)))[:,0]
assert CVaR.shape == (x_cond.shape[0],)
return CVaR
示例18
def _transform_col(self, x, col):
"""Normalize one numerical column.
Args:
x (pandas.Series): a numerical column to normalize
col (int): column index
Returns:
A normalized feature vector.
"""
return norm.ppf(self.ecdfs[col](x.values) * .998 + .001)
示例19
def __setattr__(self, name, value):
if (name == "conf_coef"):
value = norm_dist.ppf(value / 100.)
super().__setattr__(name, value)
示例20
def __setattr__(self, name, value):
if (name == "conf_coef"):
value = norm_dist.ppf(value / 100.)
super().__setattr__(name, value)
示例21
def piece_linear(self, hyper, M, prob_R):
'''
model: straight line
'''
c, slope, sigma, trans = self.split_hyper_linear(hyper)
R = np.zeros_like(M)
for i in range(4):
ind = self.indicate(M, trans, i)
mu = c[i] + M[ind]*slope[i]
R[ind] = norm.ppf(prob_R[ind], mu, sigma[i])
return R
# # Unused functions
# def classification(self, logm, trans ):
# '''
# classify as four worlds
# '''
# count = np.zeros(4)
# sample_size = len(logm)
# for iclass in range(4):
# for isample in range(sample_size):
# ind = self.indicate( logm[isample], trans[isample], iclass)
# count[iclass] = count[iclass] + ind
# prob = count / np.sum(count) * 100.
# print 'Terran %(T).1f %%, Neptunian %(N).1f %%, Jovian %(J).1f %%, Star %(S).1f %%' \
# % {'T': prob[0], 'N': prob[1], 'J': prob[2], 'S': prob[3]}
# return None
#
# def ProbRGivenM(self, radii, M, hyper):
# '''
# p(radii|M)
# '''
# c, slope, sigma, trans = self.split_hyper_linear(hyper)
# prob = np.zeros_like(M)
# for i in range(4):
# ind = self.indicate(M, trans, i)
# mu = c[i] + M[ind]*slope[i]
# sig = sigma[i]
# prob[ind] = norm.pdf(radii, mu, sig)
# prob = prob/np.sum(prob)
# return prob
示例22
def calc_bias_correction_abc(acceleration, total_curvature):
"""
Calculate the bias correction constant for the approximate bootstrap
confidence (ABC) intervals.
Parameters
----------
acceleration : 1D ndarray of scalars.
Should contain the ABC intervals' estimated acceleration constants.
total_curvature : 1D ndarray of scalars.
Should denote the ABC intervals' computred total curvature values.
Returns
-------
bias_correction : 1D ndarray of scalars.
Contains the computed bias correction for the MLE estimates that the
ABC interval is being computed for.
References
----------
Efron, Bradley, and Robert J. Tibshirani. An Introduction to the Bootstrap.
CRC press, 1994. Section 22.6, Equation 22.40, line 1.
"""
inner_arg = 2 * norm.cdf(acceleration) * norm.cdf(-1 * total_curvature)
bias_correction = norm.ppf(inner_arg)
return bias_correction
示例23
def calc_bias_correction_bca(bootstrap_replicates, mle_estimate):
"""
Calculate the bias correction for the Bias Corrected and Accelerated (BCa)
bootstrap confidence intervals.
Parameters
----------
bootstrap_replicates : 2D ndarray.
Each row should correspond to a different bootstrap parameter sample.
Each column should correspond to an element of the parameter vector
being estimated.
mle_estimate : 1D ndarray.
The original dataset's maximum likelihood point estimate. Should have
one elements for each component of the estimated parameter vector.
Returns
-------
bias_correction : 1D ndarray.
There will be one element for each element in `mle_estimate`. Elements
denote the bias correction factors for each component of the parameter
vector.
References
----------
Efron, Bradley, and Robert J. Tibshirani. An Introduction to the Bootstrap.
CRC press, 1994. Section 14.3, Equation 14.14.
"""
numerator = (bootstrap_replicates < mle_estimate[None, :]).sum(axis=0)
denominator = float(bootstrap_replicates.shape[0])
bias_correction = norm.ppf(numerator / denominator)
return bias_correction
示例24
def get_score_df(self):
'''
:return: pd.DataFrame
'''
term_freq_df = self.term_ranker_.get_ranks('')
cat_freq_df = pd.DataFrame({
'cat': term_freq_df[self.category_name],
'ncat': term_freq_df[self.not_category_names].sum(axis=1),
})
if self.neutral_category_names:
cat_freq_df['neut'] = term_freq_df[self.neutral_category_names].sum(axis=1)
cat_freq_df['all'] = cat_freq_df.sum(axis=1)
N = cat_freq_df['all'].sum()
catN = cat_freq_df['cat'].sum()
ncatN = cat_freq_df['ncat'].sum()
cat_freq_df['cat_pct'] = cat_freq_df['cat'] * 1. / catN
cat_freq_df['ncat_pct'] = cat_freq_df['ncat'] * 1. / ncatN
def row_beta_posterior(row):
return pd.Series({
'cat_p': beta(row['all'], N - row['all']).sf(row['cat'] * 1. / catN),
'ncat_p': beta(row['all'], N - row['all']).sf(row['ncat'] * 1. / ncatN),
})
p_val_df = cat_freq_df.apply(row_beta_posterior, axis=1)
cat_freq_df['cat_p'] = p_val_df['cat_p']
cat_freq_df['ncat_p'] = p_val_df['ncat_p']
cat_freq_df['cat_z'] = norm.ppf(p_val_df['cat_p'])
cat_freq_df['ncat_z'] = norm.ppf(p_val_df['ncat_p'])
cat_freq_df['score'] = None
cat_freq_df['score'][cat_freq_df['cat_pct'] == cat_freq_df['ncat_pct']] = 0
cat_freq_df['score'][cat_freq_df['cat_pct'] < cat_freq_df['ncat_pct']] = cat_freq_df['ncat_z']
cat_freq_df['score'][cat_freq_df['cat_pct'] > cat_freq_df['ncat_pct']] = -cat_freq_df['cat_z']
return cat_freq_df
示例25
def _clopper_pearson_interval(self):
r"""
Computes the Clopper-Pearson 'exact' confidence interval.
References
----------
Wikipedia contributors. (2018, July 14). Binomial proportion confidence interval.
In Wikipedia, The Free Encyclopedia. Retrieved 00:40, August 15, 2018,
from https://en.wikipedia.org/w/index.php?title=Binomial_proportion_confidence_interval&oldid=850256725
"""
p = self.x / self.n
if self.alternative == 'less':
lower_bound = 0.0
upper_bound = beta.ppf(1 - self.alpha, self.x + 1, self.n - self.x)
elif self.alternative == 'greater':
upper_bound = 1.0
lower_bound = beta.ppf(self.alpha, self.x, self.n - self.x + 1)
else:
lower_bound = beta.ppf(self.alpha / 2, self.x, self.n - self.x + 1)
upper_bound = beta.ppf(1 - self.alpha / 2, self.x + 1, self.n - self.x)
clopper_pearson_interval = {
'probability of success': p,
'conf level': 1 - self.alpha,
'interval': (lower_bound, upper_bound)
}
return clopper_pearson_interval
示例26
def _normal_scores(self):
r"""
Calculates the normal scores used in the Van der Waerden test.
Returns
-------
score_matrix : array-like
Numpy ndarray representing the data matrix with ranked observations and computed normal test scores.
Notes
-----
Let :math:`n_j`, be the number of samples for each of the :math:`k` groups where :math:`j` is the j-th group.
:math:`N` is the number of total samples in all groups, while :math:`X_{ij}` is the i-th value of the j-th
group. The normal scores used in the Van der Waerden test are calculated as:
.. math::
A_{ij} = \phi^{-1} \left( \frac{R \left( X_{ij} \right)}{N + 1} \right)
References
----------
Conover, W. J. (1999). Practical Nonparameteric Statistics (Third ed.). Wiley.
Wikipedia contributors. "Van der Waerden test." Wikipedia, The Free Encyclopedia.
Wikipedia, The Free Encyclopedia, 8 Feb. 2017. Web. 8 Mar. 2020.
"""
aij = norm.ppf(list(self.ranked_matrix[:, 2] / (self.n + 1)))
score_matrix = np.column_stack([self.ranked_matrix, aij])
return score_matrix
示例27
def _forecast_conf_int(self, forecast, fcasterr, alpha):
const = norm.ppf(1 - alpha / 2.)
conf_int = np.c_[forecast - const * fcasterr,
forecast + const * fcasterr]
return conf_int
示例28
def _forecast_conf_int(self, forecast, fcerr, alpha):
const = norm.ppf(1 - alpha/2.)
conf_int = np.c_[forecast - const*fcerr, forecast + const*fcerr]
return conf_int
示例29
def mad(a, c=Gaussian.ppf(3/4.), axis=0, center=np.median):
# c \approx .6745
"""
The Median Absolute Deviation along given axis of an array
Parameters
----------
a : array-like
Input array.
c : float, optional
The normalization constant. Defined as scipy.stats.norm.ppf(3/4.),
which is approximately .6745.
axis : int, optional
The defaul is 0. Can also be None.
center : callable or float
If a callable is provided, such as the default `np.median` then it
is expected to be called center(a). The axis argument will be applied
via np.apply_over_axes. Otherwise, provide a float.
Returns
-------
mad : float
`mad` = median(abs(`a` - center))/`c`
"""
a = np.asarray(a)
if callable(center):
center = np.apply_over_axes(center, a, axis)
return np.median((np.fabs(a-center))/c, axis=axis)
示例30
def hall_sheather(n, q, alpha=.05):
z = norm.ppf(q)
num = 1.5 * norm.pdf(z)**2.
den = 2. * z**2. + 1.
h = n**(-1. / 3) * norm.ppf(1. - alpha / 2.)**(2./3) * (num / den)**(1./3)
return h
