def cov_ellipse(cov, q=None, nsig=None, **kwargs):
""" Code is slightly modified, but essentially borrowed from:
https://stackoverflow.com/questions/18764814/make-contour-of-scatter
"""
if q is not None:
q = np.asarray(q)
elif nsig is not None:
q = 2 * norm.cdf(nsig) - 1
else:
raise ValueError('Either `q` or `nsig` should be specified')
r2 = chi2.ppf(q, 2)
val, vec = np.linalg.eigh(cov)
width, height = 2 * np.sqrt(val[:, None] * r2)
rotation = np.degrees(np.arctan2(*vec[::-1, 0]))
return width, height, rotation
python类cdf()的实例源码
def _erpac(xp, xa, n_perm, n_jobs):
"""Sub erpac function
[xp] = [xa] = (npts, ntrials)
"""
npts, ntrials = xp.shape
# Compute ERPAC
xerpac = np.zeros((npts,))
for t in range(npts):
xerpac[t] = circ_corrcc(xp[t, :], xa[t, :])[0]
# Compute surrogates:
data = Parallel(n_jobs=n_jobs)(delayed(_erpacSuro)(
xp, xa, npts, ntrials) for pe in range(n_perm))
suro = np.array(data)
# Normalize erpac:
xerpac = (xerpac - suro.mean(0))/suro.std(0)
# Get p-value:
pvalue = norm.cdf(-np.abs(xerpac))*2
return xerpac, pvalue
def _mu(distr, z, eta):
"""The non-linearity (inverse link)."""
if distr in ['softplus', 'gamma']:
mu = np.log1p(np.exp(z))
elif distr == 'poisson':
mu = z.copy()
intercept = (1 - eta) * np.exp(eta)
mu[z > eta] = z[z > eta] * np.exp(eta) + intercept
mu[z <= eta] = np.exp(z[z <= eta])
elif distr == 'gaussian':
mu = z
elif distr == 'binomial':
mu = expit(z)
elif distr == 'probit':
mu = norm.cdf(z)
return mu
def ProbabilityImprovement(self, tau, mean, std):
"""
Probability of Improvement acquisition function.
Parameters
----------
tau: float
Best observed function evaluation.
mean: float
Point mean of the posterior process.
std: float
Point std of the posterior process.
Returns
-------
float
Probability of improvement.
"""
z = (mean - tau - self.eps) / (std + self.eps)
return norm.cdf(z)
def ExpectedImprovement(self, tau, mean, std):
"""
Expected Improvement acquisition function.
Parameters
----------
tau: float
Best observed function evaluation.
mean: float
Point mean of the posterior process.
std: float
Point std of the posterior process.
Returns
-------
float
Expected improvement.
"""
z = (mean - tau - self.eps) / (std + self.eps)
return (mean - tau) * norm.cdf(z) + std * norm.pdf(z)[0]
def tExpectedImprovement(self, tau, mean, std, nu=3.0):
"""
Expected Improvement acquisition function. Only to be used with `tStudentProcess` surrogate.
Parameters
----------
tau: float
Best observed function evaluation.
mean: float
Point mean of the posterior process.
std: float
Point std of the posterior process.
Returns
-------
float
Expected improvement.
"""
gamma = (mean - tau - self.eps) / (std + self.eps)
return gamma * std * t.cdf(gamma, df=nu) + std * (1 + (gamma ** 2 - 1)/(nu - 1)) * t.pdf(gamma, df=nu)
def EI(params, means, stand_devi, parms_done, y, n_iterations, k):
"""
Expected Improvement acquisition function
:param params: test data
:param means: GP posterior mean
:param stand_devi: standard deviation
:param parms_done: training data
:param y: training targets
:return: next point that need to pick up
"""
s = 0.0005 # small value
max_mean = np.max(y)
f_max = max_mean + s
z = (means - f_max) / stand_devi
EI_vector = (means - f_max) * norm.cdf(z) + stand_devi * norm.pdf(z)
max_index = np.where(EI_vector == np.max(EI_vector))
next_point = params[max_index]
if k == n_iterations-1:
plt.subplot(2, 1, 2)
plt.plot(params, EI_vector, label='EI')
plt.legend(loc=3)
return next_point
def task(seasons, parameters):
##seasonIndexList = getSeasonIndexList(leagueId)
seasons = filterSeasonsByParams(seasons, parameters)
printTitleBox(parameters.getLevelId())
predictedTotal = 0
actualTotal = 0
for season in seasons:
print "%s %i teams," % (season.getSeasonId(), season.getTotalNumberOfTeams())
print " AWQI, APQI, AGCI, AOQI, ADQI, CPQI, Points, Points Pct \nAverages: %2.2f, %2.2f, %2.2f, %2.2f, %2.2f, %2.2f, %2.2f, %2.2f\nSigmas: %2.2f, %2.2f, %2.2f, %2.2f, %2.2f, %2.2f, %2.2f, %2.2f" % (season.getTeamStatAverage(team.Team.getAWQI), season.getTeamStatAverage(team.Team.getAPQI), season.getTeamStatAverage(team.Team.getAGCI), season.getTeamStatAverage(team.Team.getAOQI), season.getTeamStatAverage(team.Team.getADQI), season.getTeamStatAverage(team.Team.getCPQI), season.getTeamStatAverage(team.Team.getSeasonPointsTotal), season.getTeamStatAverage(team.Team.getPointsPercentage), standardDeviation(season.getTeamStatList(team.Team.getAWQI)), standardDeviation(season.getTeamStatList(team.Team.getAPQI)), standardDeviation(season.getTeamStatList(team.Team.getAGCI)), standardDeviation(season.getTeamStatList(team.Team.getAOQI)), standardDeviation(season.getTeamStatList(team.Team.getADQI)), standardDeviation(season.getTeamStatList(team.Team.getCPQI)), standardDeviation(season.getTeamStatList(team.Team.getSeasonPointsTotal)), standardDeviation(season.getTeamStatList(team.Team.getPointsPercentage)))
oneSigma = standardDeviation(season.getTeamStatList(team.Team.getCPQI))
oneGoalSigma = 0.987/oneSigma
prediction = (1.000-norm.cdf(oneGoalSigma))*season.getTotalNumberOfTeams()
actual = season.getTeamsAboveOneGoalCutoff()
print "One goal cutoff: %f sigma, above teams predicted, %.3f, actual %i\n" % (norm.cdf(oneGoalSigma), prediction, actual)
predictedTotal += prediction
actualTotal += actual
print "%s oneGoal Team totals: predicted %i, actual %i" % (parameters.getLevelId(), predictedTotal, actualTotal)
def get_p_vals(self, X):
'''
Imputes p-values from the Z-scores of `ScaledFScore` scores. Assuming incorrectly
that the scaled f-scores are normally distributed.
Parameters
----------
X : np.array
Array of word counts, shape (N, 2) where N is the vocab size. X[:,0] is the
positive class, while X[:,1] is the negative class.
Returns
-------
np.array of p-values
'''
f_scores = ScaledFScore.get_scores(X[:,0], X[:,1], self.scaler_algo, self.beta)
z_scores = (f_scores - np.mean(f_scores))/(np.std(f_scores)/np.sqrt(len(f_scores)))
return norm.cdf(z_scores)
def forward_transform(self, X):
tX = X[:, self.data["cols"]]
if(self.algo == "min-max"):
return (tX-self.data["min"])/(self.data["max"]-self.data["min"])
elif(self.algo == "normal"):
return (tX-self.data["mu"])/self.data["std"]
elif(self.algo == "inv-normal"):
return norm.cdf((tX-self.data["mu"])/self.data["std"])
elif(self.algo == "auto-normal"):
tX = (tX-self.data["min"])/(self.data["max"]-self.data["min"])
lm = self.data['boxcox'][None, :]
boxcox = lambda x: (np.sign(x)*np.abs(x)**lm-1)/lm
return (boxcox(tX)-self.data["mu"])/self.data["std"]
elif(self.algo == "auto-inv-normal"):
tX = (tX-self.data["min"])/(self.data["max"]-self.data["min"])
lm = self.data['boxcox'][None, :]
boxcox = lambda x: (np.sign(x)*np.abs(x)**lm-1)/lm
return norm.cdf(boxcox(tX), self.data["mu"], self.data["std"])
def updateInterface1D(self,solver,bbox,history,bounds):
x = np.linspace(0, 1, 80)
z=map(lambda y:bbox.queryAt([y]),x)
xx=np.atleast_2d(x).T
z_pred, sigma2_pred = solver.gp.predict(xx, eval_MSE=True)
self.ax_scat.clear()
self.ax_approx.clear()
self.ax_eimap.clear()
self.ax_scat.plot(x,np.array(z))
self.ax_scat.scatter(np.array(history)[:,0],np.array(history)[:,1])
self.ax_approx.plot(x,np.array(z_pred))
self.ax_approx.fill_between(x,np.array(z_pred)+np.array(np.sqrt(sigma2_pred)),np.array(z_pred)-np.array(np.sqrt(sigma2_pred)),alpha=0.2)
self.ax_approx.plot(x)
target=min(np.array(history)[:,1])
mean,variance=solver.gp.predict(xx,eval_MSE=True)
z=(target-mean)/np.sqrt(variance)
self.ax_approx.plot(x,np.sqrt(variance)*(z*norm.cdf(z)+norm.pdf(z)))
def update_interface_1D(ax,ax2,solver,bbox,history):
x = np.linspace(0, 1, 80)
z=map(lambda y:bbox.queryAt([y]),x)
xx=np.atleast_2d(x).T
z_pred, sigma2_pred = solver.gp.predict(xx, eval_MSE=True)
ax.clear()
ax2.clear()
ax.plot(x,np.array(z))
ax.scatter(np.array(history)[:,0],np.array(history)[:,1])
ax2.plot(x,np.array(z_pred))
ax2.fill_between(x,np.array(z_pred)+np.array(np.sqrt(sigma2_pred)),np.array(z_pred)-np.array(np.sqrt(sigma2_pred)),alpha=0.2)
ax.set_xlim([0,1])
ax2.set_xlim([0,1])
ax2.plot(x)
target=min(np.array(history)[:,1])
mean,variance=solver.gp.predict(xx,eval_MSE=True)
z=(target-mean)/np.sqrt(variance)
ax2.plot(x,np.sqrt(variance)*(z*norm.cdf(z)+norm.pdf(z)))
def update_interface_2D(ax,ax2,ax3,solver,bbox,history):
x = np.linspace(0, 1, 200)
y = np.linspace(0, 1, 200)
x,y=np.meshgrid(x, y)
xx=x.ravel()
yy=y.ravel()
#z=map(lambda x:bbox.queryAt(x),[np.array(i) for i in zip(xx,yy)])
points_pred= np.array(map(lambda s:np.array(s),zip(xx,yy)))
z_pred, sigma2_pred = solver.gp.predict(points_pred, eval_MSE=True)
#target=min(map(lambda x:bbox.queryAt(x),[np.array(i) for i in history]))
target=min(list(np.array(history)[:,1]))
u=(target-z_pred)/np.sqrt(sigma2_pred)
ax.clear()
ax2.clear()
#c1=ax.contourf(x,y,np.array(z).reshape(-1,len(x[0])))
tt=np.array(map(np.asarray,np.array(history).reshape(-1,2)[:,0]))
ax.scatter(tt[:,0],tt[:,1])
c2=ax2.contourf(x,y,np.array(z_pred).reshape(-1,len(x[0])))
c3=ax3.contourf(x,y,np.array(np.sqrt(sigma2_pred)*(u*norm.cdf(u)+norm.pdf(u))).reshape(-1,len(x[0])))
ax2.scatter(tt[:,0],tt[:,1])
ax3.scatter(tt[:,0],tt[:,1])
#c1.set_clim(min(z),max(z))
#c2.set_clim(min(z),max(z))
def cumulative_neg_binom(x, n, p):
"""
Get the cumulative probability distribution for a negative
binomial, over an array of points x that are log-scaled.
:param x: Points at which to calculate probability density
:type x: :class:`~numpy.ndarray`
:param int n: Number of trials (see :data:`~scipy.stats.nbinom`)
:param int p: Probability of success (see :data:`~scipy.stats.nbinom`)
:returns: Cumulative probability distribution over x
"""
# x is in log, so transform it back
x = list(10 ** x[1:])
# Add the point 0.0
x = [0.0] + x
return nbinom.cdf(x, n, p)
def normal(x, loc, scale):
"""
Get the probability density for a normal distribution.
:param x: Edges of bins within which to calculate probability density
:type x: :class:`~numpy.ndarray`
:param int loc: Mean of the distribution (see :data:`~scipy.stats.norm`)
:param int scale: Standard deviation of the distribution \
(see :data:`~scipy.stats.norm`)
:returns: Probability density over x. Return \
array has one less value than x, as the first value of \
the return array is the probability density between x[0] \
and x[1].
"""
norm_y = norm.cdf(x, loc, scale)
norm_y = un_cumulative(norm_y)
return sum_to_1(norm_y)
def _BlackScholesCall(S, K, T, sigma, r, q):
d1 = (log(S/K) + (r - q + (sigma**2)/2)*T)/(sigma*sqrt(T))
d2 = d1 - sigma*sqrt(T)
return S*exp(-q*T)*norm.cdf(d1) - K*exp(-r*T)*norm.cdf(d2)
def _BlackScholesPut(S, K, T, sigma, r, q):
d1 = (log(S/K) + (r - q + (sigma**2)/2)*T)/(sigma*sqrt(T))
d2 = d1 - sigma*sqrt(T)
return K*exp(-r*T)*norm.cdf(-d2) - S*exp(-q*T)*norm.cdf(-d1)
def expected_improvement(self, optimiser, x):
mu,std = optimiser.predict([x], return_std=True)
current_best = max([score for score, params in self.hyperparam_history])
gamma = (current_best - mu[0])/std[0]
exp_improv = std[0] * (gamma * norm.cdf(gamma) + norm.pdf(gamma))
return -1 * exp_improv
def probability_of_improvement(self, optimiser, x):
mu,std = optimiser.predict([x], return_std=True)
current_best = max([score for score, params in self.hyperparam_history])
gamma = (current_best - mu[0])/std[0]
return -1 * norm.cdf(gamma)
def tick_values(self, vmin, vmax) :
eps = 0.
return norm.cdf(AutoLocator.tick_values( self , norm.ppf( max(eps,vmin) ) , norm.ppf(min(1-eps,vmax)) ))
def transform_non_affine(self, a):
return norm.cdf(a)
def cdf(self, samples):
'''
Calculates the cumulative distribution function.
Parameters
----------
samples : array_like
n-by-2 matrix of samples where n is the number of samples.
Returns
-------
ndarray
Cumulative distribution function evaluated at `samples`.
'''
return np.exp(self.logcdf(samples))
def _ccdf(self, samples):
vals = np.zeros(samples.shape[0])
# Avoid subtraction of infinities
neqz = np.bitwise_and(np.any(samples > 0.0, axis=1),
np.any(samples < 1.0, axis=1))
nrvs = norm.ppf(samples[neqz, :])
vals[neqz] = norm.cdf((nrvs[:, 0] - self.theta * nrvs[:, 1])
/ np.sqrt(1 - self.theta**2))
vals[np.invert(neqz)] = norm.cdf(0.0)
return vals
def _ppcf(self, samples):
nrvs = norm.ppf(samples)
vals = norm.cdf(nrvs[:, 0] * np.sqrt(1 - self.theta**2)
+ self.theta * nrvs[:, 1])
return vals
def itransform(self, y_transformed):
# FIXME: Interpolate rather than solve for speed?
ycdf = norm.cdf(y_transformed)
y = [brentq(self._obj, a=self.lb, b=self.ub, args=(yi,))
for yi in ycdf]
return np.array(y)
def _get_levels(self):
sigma2d = self.parent.config["sigma2d"]
if sigma2d:
levels = 1.0 - np.exp(-0.5 * self.parent.config["sigmas"] ** 2)
else:
levels = 2 * norm.cdf(self.parent.config["sigmas"]) - 1.0
return levels
def obrien_fleming(information_fraction, alpha=0.05):
"""
Calculate an approximation of the O'Brien-Fleming alpha spending function.
Args:
information_fraction (scalar or array_like): share of the information
amount at the point of evaluation, e.g. the share of the maximum
sample size
alpha: type-I error rate
Returns:
float: redistributed alpha value at the time point with the given
information fraction
"""
return (1 - norm.cdf(norm.ppf(1 - alpha / 2) / np.sqrt(information_fraction))) * 2
def calculate_prob(self, multiplicity: int, avg_coverage: float):
mu = self.mu * multiplicity
sigma = self.sigma * multiplicity
difference = math.fabs(mu - avg_coverage)
return 2*norm.cdf(mu - difference, loc=mu, scale=sigma)
def testCumulative(self):
from scipy.stats import norm
x = numpy.linspace(-5., 5., 200)
samples = numpy.random.normal(size=100000)
cumula = Cumulative(samples)
c = cumula(x)
cx = norm.cdf(x)
for i in range(199):
self.assertAlmostEqual(cx[i], c[i], delta=1e-2)
def CalcSupervDelta(self, Superv_Vol=None):
if not Superv_Vol:
return 1 if self.BuySell=="Buy" else -1
if (self.UnderlyingPrice * self.StrikePrice < 0):
num = self.UnderlyingPrice - self.StrikePrice
else:
num = (log(self.UnderlyingPrice/self.StrikePrice)+0.5*Superv_Vol**2*self.Si)
den = Superv_Vol*self.Si**0.5
temp = num/den
flip = 1 if self.BuySell=="Buy" else -1
return flip*norm.cdf(temp) if self.OptionType == 'Call' else flip*-norm.cdf(-temp)
