def GaussianCdfInverse(p, mu=0, sigma=1):
"""Evaluates the inverse CDF of the gaussian distribution.
See http://en.wikipedia.org/wiki/Normal_distribution#Quantile_function
Args:
p: float
mu: mean parameter
sigma: standard deviation parameter
Returns:
float
"""
x = ROOT2 * erfinv(2 * p - 1)
return mu + x * sigma
python类erfinv()的实例源码
def GaussianCdfInverse(p, mu=0, sigma=1):
"""Evaluates the inverse CDF of the gaussian distribution.
See http://en.wikipedia.org/wiki/Normal_distribution#Quantile_function
Args:
p: float
mu: mean parameter
sigma: standard deviation parameter
Returns:
float
"""
x = ROOT2 * erfinv(2 * p - 1)
return mu + x * sigma
def normal_inverse_cdf(p: T.Tensor) -> T.Tensor:
"""
Elementwise inverse cumulative distribution function of the standard normal
distribution.
For the inverse CDF of a normal distributon with mean u and standard deviation sigma,
use u + sigma * normal_inverse_cdf(p).
Args:
p (tensor bounded in (0,1))
Returns:
tensor: Elementwise inverse cdf
"""
return SQRT2*erfinv(2*numpy.clip(p, a_min=EPSILON, a_max=1-EPSILON)-1)
def compute_cutoff_beta(T, frac=0.99):
"""Returns the velocity for which the fraction frac of a beam emitted from a thermionic
cathode with temperature T move more slowly in the longitudinal (z) direction.
Arguments:
T (float) : temperature of the cathode in K
frac (Optional[float]) : Fraction of beam with vz < cutoff. Defaults to 0.99.
Returns:
beta_cutoff (float) : cutoff velocity divided by c.
"""
sigma = np.sqrt(2*kb_J*T/m) # effective sigma for half-Gaussian distribution
multiplier = erfinv(frac) # the multiplier on sigma accounting for frac of the distribution
return multiplier*sigma/c
def prior_cdf(self, u):
"""Evaluate cumulative density function."""
value = (erfinv(2.0 * u - 1.0) * np.sqrt(2.)) * self._sigma + self._mu
value = (value - self._min_value) / (self._max_value - self._min_value)
return np.clip(value, 0.0, 1.0)
def normal(sample, avg=0.0, std=1.0):
return avg + std * np.sqrt(2) * erfinv(2 * sample - 1)
def fit(self, y):
x = y[~np.isnan(y)]
n = len(y)
self.y = erfinv(np.linspace(0., 1., n + 2)[1:-1])
self.s = np.sort(x)
def __init__(self, params):
super(GaussianRandomField, self).__init__(params) # don't forget this line in our classes!
# value of the correlation function at the origin
self.corr_func_at_origin = self.frac_volume * (1.0 - self.frac_volume)
# inverse slope of the normalized correlation function at the origin
beta = np.sqrt(2) * erfinv(2 * (1-self.frac_volume) - 1)
# second derivative of the field acf at the origin
acf_psi_doubleprime = -1.0 / 2 * ( (1.0/self.corr_length)**2 + 1.0 / 3 * (2 * np.pi / self.repeat_distance)**2 )
SSA_tilde = 2.0 / np.pi * np.exp( - beta**2 / 2) * np.sqrt( -acf_psi_doubleprime ) / self.frac_volume
self.inv_slope_at_origin = 4.0 * (1 - self.frac_volume) / SSA_tilde
def autocorrelation_function(self, r):
"""compute the real space autocorrelation function for the Gaussian random field model
"""
# compute the cut-level parameter beta
beta = np.sqrt(2) * erfinv(2 * (1-self.frac_volume) - 1)
# the covariance of the GRF
acf_psi = (np.exp(-r/self.corr_length) * (1 + r / self.corr_length)
* np.sinc(2*r/self.repeat_distance))
# integral discretization. henning says: the resolution 1e-2 is ad hoc, test required,
# the integrand has a (integrable) singularity for t=1 and acf_psi = 1, so an adaptive
# discretization seems preferable -> TODO
dt = 1e-2
t = np.arange(0, 1, dt)
# the gridded integrand, via change of integration variable
# compared to the wp-2 docu, to enable array-based computation
t_gridded, acf_psi_gridded = np.meshgrid(t, acf_psi)
integrand_gridded = (acf_psi_gridded / np.sqrt(1 - (t_gridded * acf_psi_gridded)**2)
* np.exp( - beta**2 / (1 + t_gridded * acf_psi_gridded)))
acf = 1.0 / (2 * np.pi) * np.trapz(integrand_gridded, x=t_gridded)
return acf
# ft not known analytically: deligate
# def ft_autocorrelation_function(self, k):
def ndpac(pha, amp, p, optimize):
"""Normalized direct Pac (Ozkurt, 2012).
Parameters
----------
pha : array_like
Array of phases of shapes (npha, ..., npts)
amp : array_like
Array of amplitudes of shapes (namp, ..., npts)
p : float
The p-value to use.
Returns
-------
pac : array_like
PAC of shape (npha, namp, ...)
"""
npts = amp.shape[-1]
# Normalize amplitude :
np.subtract(amp, np.mean(amp, axis=-1, keepdims=True), out=amp)
np.divide(amp, np.std(amp, axis=-1, keepdims=True), out=amp)
# Compute pac :
pac = np.abs(np.einsum('i...j, k...j->ik...', amp, np.exp(1j * pha),
optimize=optimize))
pac *= pac / npts
# Set to zero non-significant values:
xlim = erfinv(1 - p)**2
pac[pac <= 2 * xlim] = 0.
return pac
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 normal(sample, avg=0.0, std=1.0):
return avg + std * np.sqrt(2) * erfinv(2 * sample - 1)
def normal(sample, avg=0.0, std=1.0):
return avg + std * np.sqrt(2) * erfinv(2 * sample - 1)
def lineOfSightGrid( theta, phi, Ntheta, Nphi, sigmaTheta, pole=None ):
'''
assumes theta, phi are in Earth fixed coordinates.
generates a grid in line of sight coordinates around the corresponding theta_los, phi_los
returns that grid in Earth fixed coords
return thetaGRID, phiGRID
'''
theta_los, phi_los = ThetaPhi2LineOfSight(theta, phi, pole=pole)
# phi is easy, just wrap it around the azimuth
phiGRID = np.linspace(0, 2*np.pi, Nphi+1)[:-1]
# theta is a bit ad hoc, but should be a gaussian centered on the injected location
if Ntheta%2:
Ntheta += 1
thetaGRID = erfinv(np.linspace(0,1,Ntheta/2+1)[:-1])
thetaGRID = theta_los + sigmaTheta*np.concatenate((-thetaGRID[::-1], thetaGRID[1:]))
thetaGRID = thetaGRID[(thetaGRID>=0)*(thetaGRID<=np.pi)]
### rotate back to Earth-fixed
thetaGRID, phiGRID = np.meshgrid(thetaGRID, phiGRID)
thetaGRID, phiGRID = lineOfSight2ThetaPhi(thetaGRID.flatten(), phiGRID.flatten(), pole=pole)
phiGRID[phiGRID<0] += 2*np.pi ### we want these to be positive
return thetaGRID, phiGRID
def ka_erfinv_rank_transform(x):
'''
This is used on numeric variable, after doing this operation, one should do MM and SS on all dataset.
'''
mm = MinMaxScaler()
tmp = erfinv(np.clip(np.squeeze(mm.fit_transform(rankdata(x).reshape(-1,1))), 0, 0.999999999))
tmp = tmp - np.mean(tmp)
return tmp
def histClim( imData, cutoff = 0.01, bins_ = 512 ):
'''Compute display range based on a confidence interval-style, from a histogram
(i.e. ignore the 'cutoff' proportion lowest/highest value pixels)'''
if( cutoff <= 0.0 ):
return imData.min(), imData.max()
# compute image histogram
hh, bins_ = imHist(imData, bins_)
hh = hh.astype( 'float' )
# number of pixels
Npx = np.sum(hh)
hh_csum = np.cumsum( hh )
# Find indices where hh_csum is < and > Npx*cutoff
try:
i_forward = np.argwhere( hh_csum < Npx*(1.0 - cutoff) )[-1][0]
i_backward = np.argwhere( hh_csum > Npx*cutoff )[0][0]
except IndexError:
print( "histClim failed, returning confidence interval instead" )
from scipy.special import erfinv
sigma = np.sqrt(2) * erfinv( 1.0 - cutoff )
return ciClim( imData, sigma )
clim = np.array( [bins_[i_backward], bins_[i_forward]] )
if clim[0] > clim[1]:
clim = np.array( [clim[1], clim[0]] )
return clim
def normal(sample, avg=0.0, std=1.0):
return avg + std * np.sqrt(2) * erfinv(2 * sample - 1)
def erfinv(x: T.Tensor) -> T.Tensor:
"""
Elementwise error function of a tensor.
Args:
x: A tensor.
Returns:
tensor: Elementwise error function
"""
return special.erfinv(x)
