def skewgauss(x, sigma, mu, alpha, a):
"""Skewed Gaussian.
Args:
x (np.ndarray): Dataset.
sigma (float): Standard deviation.
mu (float): Mean.
alpha (float): Skewness.
a (float): Normalization factor.
Returns:
skewed gaussian (np.ndarray): Skewed Gaussian.
"""
normpdf = (1 / (sigma * np.sqrt(2 * np.pi))) * np.exp(-(x - mu)**2 / (2 * sigma**2))
normcdf = 0.5 * (1 + erf((alpha * ((x - mu) / sigma)) / (np.sqrt(2))))
return 2 * a * normpdf * normcdf
python类erf()的实例源码
def test_time(res, off, t, signal):
"""Time domain analytical half-space solution.
- Source at x = y = z = 0 m
- Receiver at y = z = 0 m; x = off
- Resistivity of halfspace res
- Times t, t > 0 s
- Impulse response if signal = 0
- Switch-on response if signal = 1
"""
tau = np.sqrt(mu_0*off**2/(res*t))
fact1 = res/(2*np.pi*off**3)
fact2 = tau/np.sqrt(np.pi)
if signal == 0:
return fact1*tau**3/(4*t*np.sqrt(np.pi))*np.exp(-tau**2/4)
else:
resp = fact1*(2 - special.erf(tau/2) + fact2*np.exp(-tau**2/4))
if signal < 0:
DC = test_time(res, off, 1000000, 1)
resp = DC-resp
return resp
# Time-domain solution
def particle(x0, y0, d):
sigma = 0.42*d
C = np.pi/8.*sigma**2
out = (erf((x - x0 + 0.5)/(sigma*np.sqrt(2))) - erf((x - x0 - 0.5)/(sigma*np.sqrt(2))))
out *= (erf((y - y0 + 0.5)/(sigma*np.sqrt(2))) - erf((y - y0 - 0.5)/(sigma*np.sqrt(2))))
return C*out
def sequential_colormap(x):
from numpy import array
from scipy.special import erf
x = array(x)
if any(x < 0) or any(x > 1):
raise ValueError('x must be between 0 and 1 inclusive.')
red = 1.000 - 0.392 * (1 + erf((x - 0.869) / 0.255))
grn = 1.021 - 0.456 * (1 + erf((x - 0.527) / 0.376))
blu = 1.000 - 0.493 * (1 + erf((x - 0.272) / 0.309))
return array([red, grn, blu]).T
def __init__(self, x, min_limit=-np.inf, max_limit=np.inf, weights=1.0):
self.points = x
self.N = x.size
self.min_limit = min_limit
self.max_limit = max_limit
self.sigma = get_sigma(x, min_limit=min_limit, max_limit=max_limit)
self.weights = (2
/ (erf((max_limit - x) / (np.sqrt(2.) * self.sigma))
- erf((min_limit - x) / (np.sqrt(2.) * self.sigma)))
* weights)
self.W_sum = np.sum(self.weights)
def Erf2D(x, y, xc, yc, a, b, sigma):
'''
A simple 2D error function PSF model.
'''
c = np.sqrt(2) * sigma
ex = (erf((x + 1 - xc) / c) - erf((x - xc) / c))
ey = (erf((y + 1 - yc) / c) - erf((y - yc) / c))
return a * ex * ey + b
def _cdf_y_gt_f(self, y, f):
return (erf((y - f) / (self.l * np.sqrt(2))) * self.l * SQRTPI2 + f) \
/ (f + self.l * SQRTPI2)
def _cdf_f_lt_0(self, y, f):
return erf((y - f) / (self.l * np.sqrt(2)))
def gaussian_sum(self, centers, weights, settings):
"""Calculates a discrete version of a sum of Gaussian distributions.
The calculation is done through the cumulative distribution function
that is better at keeping the integral of the probability function
constant with coarser grids.
The values are normalized by dividing with the maximum value of a
gaussian with the given standard deviation.
Args:
centers (1D np.ndarray): The means of the gaussians.
weights (1D np.ndarray): The weights for the gaussians.
settings (dict): The grid settings
Returns:
Value of the gaussian sums on the given grid.
"""
start = settings["min"]
stop = settings["max"]
sigma = settings["sigma"]
n = settings["n"]
max_val = 1/(sigma*math.sqrt(2*math.pi))
dx = (stop - start)/(n-1)
x = np.linspace(start-dx/2, stop+dx/2, n+1)
pos = x[np.newaxis, :] - centers[:, np.newaxis]
y = weights[:, np.newaxis]*1/2*(1 + erf(pos/(sigma*np.sqrt(2))))
f = np.sum(y, axis=0)
f /= max_val
f_rolled = np.roll(f, -1)
pdf = (f_rolled - f)[0:-1]/dx # PDF is the derivative of CDF
return pdf
def _cdf(z):
"""Cumulative density function (CDF) of the standard normal distribution."""
return .5 * (1. + erf(z/np.sqrt(2.)))
def truncated_normal(shape=None, mu=0., sigma=1., x_min=None, x_max=None):
"""
Generates random variates from a lower-and upper-bounded normal distribution
@param shape: shape of the random sample
@param mu: location parameter
@param sigma: width of the distribution (sigma >= 0.)
@param x_min: lower bound of variate
@param x_max: upper bound of variate
@return: random variates of lower-bounded normal distribution
"""
from scipy.special import erf, erfinv
from numpy.random import standard_normal
from numpy import inf, sqrt
if x_min is None and x_max is None:
return standard_normal(shape) * sigma + mu
elif x_min is None:
x_min = -inf
elif x_max is None:
x_max = inf
x_min = max(-1e300, x_min)
x_max = min(+1e300, x_max)
var = sigma ** 2 + 1e-300
sigma = sqrt(2 * var)
a = erf((x_min - mu) / sigma)
b = erf((x_max - mu) / sigma)
return probability_transform(shape, erfinv, a, b) * sigma + mu
def anomalies(log_z, row_label=None, prefix=''):
from scipy.special import erf
ns = log_z.shape[1]
if row_label is None:
row_label = list(map(str, range(ns)))
a_score = np.sum(log_z[:, :, 0], axis=0)
mean, std = np.mean(a_score), np.std(a_score)
a_score = (a_score - mean) / std
percentile = 1. / ns
anomalies = np.where(0.5 * (1 - erf(a_score / np.sqrt(2)) ) < percentile)[0]
f = safe_open(prefix + '/text_files/anomalies.txt', 'w+')
for i in anomalies:
f.write(row_label[i] + ', %0.1f\n' % a_score[i])
f.close()
def flow_para_function_with_vibration( x, beta, relaxation_rate, flow_velocity, freq, amp, baseline=1):
vibration_part = (1 + amp*np.cos( 2*np.pi*freq* x) )
Diff_part= np.exp(-2 * relaxation_rate * x)
Flow_part = np.pi**2/(16*x*flow_velocity) * abs( erf( np.sqrt( 4/np.pi * 1j* x * flow_velocity ) ) )**2
return beta* vibration_part* Diff_part * Flow_part + baseline
def flow_para_function( x, beta, relaxation_rate, flow_velocity, baseline=1):
'''flow_velocity: q.v (q vector dot v vector = q*v*cos(angle) )'''
Diff_part= np.exp(-2 * relaxation_rate * x)
Flow_part = np.pi**2/(16*x*flow_velocity) * abs( erf( np.sqrt( 4/np.pi * 1j* x * flow_velocity ) ) )**2
return beta*Diff_part * Flow_part + baseline
def flow_para_function_explicitq( x, beta, relaxation_rate, flow_velocity, baseline=1, qr=1, q_ang=0 ):
'''Nov 9, 2017 Basically, make q vector to (qr, angle), relaxation_rate is actually a diffusion rate
flow_velocity: q.v (q vector dot v vector = q*v*cos(angle) )'''
Diff_part= np.exp(-2 * relaxation_rate* qr**2 * x)
Flow_part = np.pi**2/(16*x*flow_velocity*qr* np.cos(q_ang) ) * abs( erf( np.sqrt( 4/np.pi * 1j* x * flow_velocity * qr* np.cos(q_ang) ) ) )**2
return beta*Diff_part * Flow_part + baseline
def stretched_flow_para_function( x, beta, relaxation_rate, alpha, flow_velocity, baseline=1):
'''
flow_velocity: q.v (q vector dot v vector = q*v*cos(angle) )
'''
Diff_part= np.exp(-2 * (relaxation_rate * x)**alpha )
Flow_part = np.pi**2/(16*x*flow_velocity) * abs( erf( np.sqrt( 4/np.pi * 1j* x * flow_velocity ) ) )**2
return beta*Diff_part * Flow_part + baseline
def shape(xs, center, amplitude, sigma, gamma):
norm = (amplitude) / (sigma * sqrt(2 * pi)) * \
exp(-((xs - center) ** 2) / (2 * sigma ** 2))
skew = (1 + erf((gamma * (xs - center)) / (sigma * sqrt(2))))
return norm * skew
def evaluate(self, F, xo, yo, s):
return (F / 4 *
((erf((self.x - xo + 0.5) / (np.sqrt(2) * s)) -
erf((self.x - xo - 0.5) / (np.sqrt(2) * s))) *
(erf((self.y - yo + 0.5) / (np.sqrt(2) * s)) -
erf((self.y - yo - 0.5) / (np.sqrt(2) * s)))))
def errfunc(c):
arg = 1.00 * c / sqrt(2.00)
return special.erf(arg);
def plot(begin, end, step):
c = cvector(begin, end, step)
erf = errfunc(c)
print(erf);
plt.plot(c, erf);
plt.xlabel('0 <= c <= 5');
plt.ylabel('I(c)')
plt.show()
def calc_theta_i(mag, mag_err, maxmag, limmag):
"""
Calculate theta_i. This is reproduced from calclambda_chisq_theta_i.pr
parameters
----------
mag:
mag_err:
maxmag:
limmag:
returns
-------
theta_i:
"""
theta_i = np.ones((len(mag)))
eff_lim = np.clip(maxmag,0,limmag)
dmag = eff_lim - mag
calc = dmag < 5.0
N_calc = np.count_nonzero(calc==True)
if N_calc > 0: theta_i[calc] = 0.5 + 0.5*erf(dmag[calc]/(np.sqrt(2)*mag_err[calc]))
hi = mag > limmag
N_hi = np.count_nonzero(hi==True)
if N_hi > 0:
theta_i[hi] = 0.0
return theta_i
def StandardGaussianCdf(x):
"""Evaluates the CDF of the standard Gaussian distribution.
See http://en.wikipedia.org/wiki/Normal_distribution
#Cumulative_distribution_function
Args:
x: float
Returns:
float
"""
return (erf(x / ROOT2) + 1) / 2
def StandardGaussianCdf(x):
"""Evaluates the CDF of the standard Gaussian distribution.
See http://en.wikipedia.org/wiki/Normal_distribution
#Cumulative_distribution_function
Args:
x: float
Returns:
float
"""
return (erf(x / ROOT2) + 1) / 2
def erf(x: T.Tensor) -> T.Tensor:
"""
Elementwise error function of a tensor.
Args:
x: A tensor.
Returns:
tensor: Elementwise error function
"""
return special.erf(x)
def normal_cdf(x: T.Tensor) -> T.Tensor:
"""
Elementwise cumulative distribution function of the standard normal distribution.
For the CDF of a normal distributon with mean u and standard deviation sigma, use
normal_cdf((x-u)/sigma).
Args:
x (tensor)
Returns:
tensor: Elementwise cdf
"""
return 0.5 * (1 + erf(x/SQRT2))
def Edeta(deta,x):
if deta != 0:
g = 0.5*(1+erf(x/deta))
return g
else:
g = np.zeros(x.shape)
g[x >= 0 ] =1
return g
def predict_proba(self, X):
"""Return probability estimates for the test vector X.
Parameters
----------
X : array-like, shape = (n_samples, n_features)
Returns
-------
C : array-like, shape = (n_samples, n_classes)
Returns the probability of the samples for each class in
the model. The columns correspond to the classes in sorted
order, as they appear in the attribute ``classes_``.
"""
check_is_fitted(self, ["X_train_", "y_train_", "pi_", "W_sr_", "L_"])
# Based on Algorithm 3.2 of GPML
K_star = self.kernel_(self.X_train_, X) # K_star =k(x_star)
f_star = K_star.T.dot(self.y_train_ - self.pi_) # Line 4
v = solve(self.L_, self.W_sr_[:, np.newaxis] * K_star) # Line 5
# Line 6 (compute np.diag(v.T.dot(v)) via einsum)
var_f_star = self.kernel_.diag(X) - np.einsum("ij,ij->j", v, v)
# Line 7:
# Approximate \int log(z) * N(z | f_star, var_f_star)
# Approximation is due to Williams & Barber, "Bayesian Classification
# with Gaussian Processes", Appendix A: Approximate the logistic
# sigmoid by a linear combination of 5 error functions.
# For information on how this integral can be computed see
# blitiri.blogspot.de/2012/11/gaussian-integral-of-error-function.html
alpha = 1 / (2 * var_f_star)
gamma = LAMBDAS * f_star
integrals = np.sqrt(np.pi / alpha) \
* erf(gamma * np.sqrt(alpha / (alpha + LAMBDAS**2))) \
/ (2 * np.sqrt(var_f_star * 2 * np.pi))
pi_star = (COEFS * integrals).sum(axis=0) + .5 * COEFS.sum()
return np.vstack((1 - pi_star, pi_star)).T
def compute_crossing_fraction(cathode_temp, phi, zmesh):
"""
Returns the % of particles expected to cross the gap
Arguments:
cathode_temp (float) : cathode temperature in K
phi (scipy) : Interpolation (scipy.interpolate.interpolate.interp1d) of Phi from initial solve
zmesh (ndArray) : 1D array with positions of the z-coordinates of the grid
Returns:
e_frac (float) : fraction of beam that overcomes the peak potential barrier phi
"""
# Conditions at cathode edge
var_xy = kb_J*cathode_temp/m # Define the variance of the distribution in the x,y planes
var_z = 2*kb_J*cathode_temp/m # Variance in z plane is twice that in the horizontal
#Phi is in eV
vz_phi = np.sqrt(2.*e*np.max(phi(zmesh))/m_e)
#cumulative distribution function for a Maxwellian
CDF_Maxwell = lambda v, a: erf(v/(np.sqrt(2)*a)) - np.sqrt(2./np.pi)*v*np.exp(-1.*v**2/(2.*a**2))/a
e_frac = CDF_Maxwell(vz_phi,np.sqrt(var_z)) #fraction with velocity below vz_phi
return (1.-e_frac)*100. #Multiply by 100 to get % value
def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, t, current=1., length=1., orientation='X', kappa=0., epsr=1.):
"""
Computing the analytic electric fields (E) from an electrical dipole in a wholespace
- You have the option of computing E for multiple times at a single reciever location
or a single time at multiple locations
:param numpy.array XYZ: reciever locations at which to evaluate E
:param numpy.array srcLoc: [x,y,z] triplet defining the location of the electric dipole source
:param float sig: value specifying the conductivity (S/m) of the wholespace
:param numpy.array t: array of times at which to measure
:param float current: size of the injected current (A), default is 1.0 A
:param float length: length of the dipole (m), default is 1.0 m
:param str orientation: orientation of dipole: 'X', 'Y', or 'Z'
:param float kappa: magnetic susceptiblity value (unitless), default is 0.
:param float epsr: relative permitivitty value (unitless), default is 1.0
:rtype: numpy.array
:return: Ex, Ey, Ez: arrays containing all 3 components of E evaluated at the specified locations and times.
"""
mu = mu_0*(1+kappa)
epsilon = epsilon_0*epsr
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
# Check
if XYZ.shape[0] > 1 & t.shape[0] > 1:
raise Exception("I/O type error: For multiple field locations only a single time can be specified.")
dx = XYZ[:,0]-srcLoc[0]
dy = XYZ[:,1]-srcLoc[1]
dz = XYZ[:,2]-srcLoc[2]
r = np.sqrt( dx**2. + dy**2. + dz**2.)
theta = np.sqrt((mu*sig)/(4*t))
front = current * length / (4.* pi * sig * r**3)
mid = 3 * erf(theta*r) - (4/np.sqrt(pi) * (theta)**3 * r**3 + 6/np.sqrt(pi) * theta * r) * np.exp(-(theta)**2 * (r)**2)
extra = (erf(theta*r) - (4/np.sqrt(pi) * (theta)**3 * r**3 + 2/np.sqrt(pi) * theta * r) * np.exp(-(theta)**2 * (r)**2))
if orientation.upper() == 'X':
Ex = front*(dx**2 / r**2)*mid - front*extra
Ey = front*(dx*dy / r**2)*mid
Ez = front*(dx*dz / r**2)*mid
return Ex, Ey, Ez
elif orientation.upper() == 'Y':
# x--> y, y--> z, z-->x
Ey = front*(dy**2 / r**2)*mid - front*extra
Ez = front*(dy*dz / r**2)*mid
Ex = front*(dy*dx / r**2)*mid
return Ex, Ey, Ez
elif orientation.upper() == 'Z':
# x --> z, y --> x, z --> y
Ez = front*(dz**2 / r**2)*mid - front*extra
Ex = front*(dz*dx / r**2)*mid
Ey = front*(dz*dy / r**2)*mid
return Ex, Ey, Ez
def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, t, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing the analytic magnetic fields (H) from an electrical dipole in a wholespace
- You have the option of computing H for multiple times at a single reciever location
or a single time at multiple locations
:param numpy.array XYZ: reciever locations at which to evaluate H
:param numpy.array srcLoc: [x,y,z] triplet defining the location of the electric dipole source
:param float sig: value specifying the conductivity (S/m) of the wholespace
:param numpy.array t: array of times at which to measure
:param float current: size of the injected current (A), default is 1.0 A
:param float length: length of the dipole (m), default is 1.0 m
:param str orientation: orientation of dipole: 'X', 'Y', or 'Z'
:param float kappa: magnetic susceptiblity value (unitless), default is 0.
:param float epsr: relative permitivitty value (unitless), default is 1.0
:rtype: numpy.array
:return: Hx, Hy, Hz: arrays containing all 3 components of H evaluated at the specified locations and times.
"""
mu = mu_0*(1+kappa)
epsilon = epsilon_0*epsr
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
# Check
if XYZ.shape[0] > 1 & t.shape[0] > 1:
raise Exception("I/O type error: For multiple field locations only a single time can be specified.")
dx = XYZ[:, 0]-srcLoc[0]
dy = XYZ[:, 1]-srcLoc[1]
dz = XYZ[:, 2]-srcLoc[2]
r = np.sqrt(dx**2. + dy**2. + dz**2.)
theta = np.sqrt((mu*sig)/(4*t))
front = (current * length) / (4.*pi*(r)**3)
mid = erf(theta*r) - (2/np.sqrt(pi)) * theta * r * np.exp(-(theta)**2 * (r)**2)
if orientation.upper() == 'X':
Hy = front * mid * -dz
Hz = front * mid * dy
Hx = np.zeros_like(Hy)
return Hx, Hy, Hz
elif orientation.upper() == 'Y':
Hx = front * mid * dz
Hz = front * mid * -dx
Hy = np.zeros_like(Hx)
return Hx, Hy, Hz
elif orientation.upper() == 'Z':
Hx = front * mid * -dy
Hy = front * mid * dx
Hz = np.zeros_like(Hx)
return Hx, Hy, Hz
