Python scipy.special.erfinv() Examples

def _function(self,
                 variable=None,
                 context=None,
                 params=None,
                 ):

        try:
            from scipy.special import erfinv
        except:
            raise FunctionError("The UniformToNormalDist function requires the SciPy package.")

        mean = self._get_current_function_param(DIST_MEAN, context)
        standard_deviation = self._get_current_function_param(STANDARD_DEVIATION, context)
        random_state = self._get_current_function_param('random_state', context)

        sample = random_state.rand(1)[0]
        result = ((np.sqrt(2) * erfinv(2 * sample - 1)) * standard_deviation) + mean

        return self.convert_output_type(result) 

def std_iqr(x):
    """Robust estimation of the standard deviation, based on the inter-quartile
    (IQR) distance of x.
    This computes the IQR 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 iqr
    from scipy.special import erfinv

    correction = 2 ** 0.5 * erfinv(0.5)
    return correction * iqr(x) 

def _rescale(self, samp, **kwargs):
        try:
            mode = kwargs['mode']
        except KeyError:
            mode = None

        if mode is None:
            if self.nmodes == 1:
                mode = 0
            else:
                mode = np.argwhere(self.cumweights - np.random.rand() > 0)[0][0]

        samp = erfinv(2. * samp - 1) * 2. ** 0.5

        # rotate and scale to the multivariate normal shape
        samp = self.mus[mode] + self.sigmas[mode] * np.einsum('ij,kj->ik',
                                                              samp * self.sqeigvalues[mode],
                                                              self.eigvectors[mode])
        return samp 

def test_erfinv() -> None:
    try:
        from scipy.special import erfinv as scipy_erfinv
    except:
        pytest.skip("scipy not installed skipping test for erf")

    x = np.linspace(-1.0 + 1.0e-4, 1 - 1.0e-4, 11)
    y_scipy = scipy_erfinv(x)

    # Test mx.nd
    y_mxnet = util.erfinv(mx.nd, mx.nd.array(x)).asnumpy()
    assert np.allclose(y_mxnet, y_scipy, rtol=1e-3)

    # Test mx.sym
    X = mx.symbol.Variable("x")
    func = util.erfinv(mx.sym, X)
    func_exec = func.bind(ctx=mx.cpu(), args={"x": mx.nd.array(x)})
    func_exec.forward()
    y_mxnet_sym = func_exec.outputs[0].asnumpy()
    assert np.allclose(y_mxnet_sym, y_scipy, rtol=1e-3)

    # Text np
    y_np = util.erfinv(np, x)
    assert np.allclose(y_np, y_scipy, rtol=1e-3) 

def date_uncertainty_due_to_rate(self, node, interval=(0.05, 0.095)):
        """use previously calculated variation of the rate to estimate
        the uncertainty in a particular numdate due to rate variation.

        Parameters
        ----------
        node : PhyloTree.Clade
            node for which the confidence interval is to be calculated
        interval : tuple, optional
            Array of length two, or tuple, defining the bounds of the confidence interval

        """
        if hasattr(node, "numdate_rate_variation"):
            from scipy.special import erfinv
            nsig = [np.sqrt(2.0)*erfinv(-1.0 + 2.0*x) if x*(1.0-x) else 0
                    for x in interval]
            l,c,u = [x[1] for x in node.numdate_rate_variation]
            return np.array([c + x*np.abs(y-c) for x,y in zip(nsig, (l,u))])

        else:
            return None 

def _using_cdf(self, x, exp):

        # Get the value of the cumulative probability function, instead of the survival function (1 - cdf),
        # because for extreme values sf(x) = 1 - cdf(x) = 1 due to numerical precision problems

        cdf = scipy.stats.poisson.cdf(x, exp)

        # print(cdf)

        out = np.zeros_like(x)

        idx = cdf >= 2 * self._epsilon

        # We can do a direct computation, because the numerical precision is sufficient
        # for this computation, as -sf = cdf - 1 is a representable number

        out[idx] = erfinv(2 * cdf[idx] - 1) * sqrt(2)

        # We use a lookup table with interpolation because the numerical precision would not
        # be sufficient to make the computation

        out[~idx] = -1 * self._interpolator(np.log10(cdf[~idx]))

        return out 

def _using_sf(self, x, exp):

        sf = scipy.stats.poisson.sf(x, exp)

        # print(sf)

        # return erfinv(2 * sf) * sqrt(2)

        return scipy.stats.norm.isf(sf) 

def animate(i):
    param = np.exp((float(i)-50)/10)
    #pp = ((1-(1-p)**param) + (p**param))/2

    pp = (erf( erfinv( p*2-1 ) * param )+1)/2

    fig.suptitle("power %f"%float(param))
    im1[0].set_ydata(pp)
    return im1 

def harmonic_1d_data():
    N = 201
    delta = 1. / (N+1)
    xgrid = 2 * np.arange(1, N+1) * delta - 1.
    x = np.sqrt(2) * erfinv(xgrid)
    return x.reshape(-1, 1) 

def _normalcdfinv(self, param, p):
        """The inverse CDF of the normal distribution, without bounds."""
        mu = self._mean[param]
        var = self._var[param]
        return mu + (2*var)**0.5 * erfinv(2*p - 1.) 

def rescale(self, val):
        """
        'Rescale' a sample from the unit line element to the appropriate LogNormal prior.

        This maps to the inverse CDF. This has been analytically solved for this case.
        """
        self.test_valid_for_rescaling(val)
        return np.exp(self.mu + np.sqrt(2 * self.sigma ** 2) * erfinv(2 * val - 1)) 

def rescale(self, val):
        """
        'Rescale' a sample from the unit line element to the appropriate truncated Gaussian prior.

        This maps to the inverse CDF. This has been analytically solved for this case.
        """
        self.test_valid_for_rescaling(val)
        return erfinv(2 * val * self.normalisation + erf(
            (self.minimum - self.mu) / 2 ** 0.5 / self.sigma)) * 2 ** 0.5 * self.sigma + self.mu 

def rescale(self, val):
        """
        'Rescale' a sample from the unit line element to the appropriate Gaussian prior.

        Parameters
        ----------
        val: Union[float, int, array_like]

        This maps to the inverse CDF. This has been analytically solved for this case.
        """
        self.test_valid_for_rescaling(val)
        return self.mu + erfinv(2 * val - 1) * 2 ** 0.5 * self.sigma 

def normal(sample, avg=0.0, std=1.0):
    return avg + std * np.sqrt(2) * erfinv(2 * sample - 1) 

def test_erfinv(self):
        i = special.erfinv(0)
        assert_equal(i,0) 

def normal(sample, avg=0.0, std=1.0):
    return avg + std * np.sqrt(2) * erfinv(2 * sample - 1) 

def unit2gaussian(x, mu, sigma):
    """
    mapping from uniform distribution on unit hypercube
    to truncated gaussian distribution on parameter space, 
    with mean 'mu' and std dev 'sigma'

    from Handley+15, eq. (A9)
    """
    return mu + SQRT2 * sigma * special.erfinv(2*x - 1) 

def spectral_rad_ppf(self, u):
        """Radial spectral ppf.

        Notes
        -----
        Not defined for 3D.
        """
        u = np.array(u, dtype=np.double)
        if self.dim == 1:
            return sps.erfinv(u) * np.sqrt(np.pi) / self.len_scale
        if self.dim == 2:
            return np.sqrt(np.pi) / self.len_scale * np.sqrt(-np.log(1.0 - u))
        return None 

def _zinnharvey(field, conn="high", mean=None, var=None):
    """
    Zinn and Harvey transformation to connect low or high values.

    Parameters
    ----------
    field : :class:`numpy.ndarray`
        Spatial Random Field with normal distributed values.
        As returned by SRF.
    conn : :class:`str`, optional
        Desired connectivity. Either "low" or "high".
        Default: "high"
    mean : :class:`float` or :any:`None`, optional
        Mean of the given field. If None is given, the mean will be calculated.
        Default: :any:`None`
    var : :class:`float` or :any:`None`, optional
        Variance of the given field.
        If None is given, the mean will be calculated.
        Default: :any:`None`

    Returns
    -------
        :class:`numpy.ndarray`
            Transformed field.
    """
    if mean is None:
        mean = np.mean(field)
    if var is None:
        var = np.var(field)
    field = np.abs((field - mean) / np.sqrt(var))
    field = 2 * erf(field / np.sqrt(2)) - 1
    field = np.sqrt(2) * erfinv(field)
    if conn == "high":
        field = -field
    return field * np.sqrt(var) + mean 

def test_erfinv(self):
        i = special.erfinv(0)
        assert_equal(i,0) 

def erfinv(x):
    """Elementwise inverse function of error function.

    .. note::
       Forward computation in CPU cannot be done if
       `SciPy <https://www.scipy.org/>`_ is not available.

    Args:
        x (:class:`~chainer.Variable` or :ref:`ndarray`): Input variable.

    Returns:
        ~chainer.Variable: Output variable.
    """
    return ErfInv().apply((x,))[0] 

def forward_gpu(self, x):
        self.retain_outputs((0,))
        return cuda.elementwise(
            'T x', 'T y',
            'y = erfinv(x)',
            'elementwise_erfinv',
        )(x[0]), 

def forward_cpu(self, x):
        if not available_cpu:
            raise ImportError('SciPy is not available. Forward computation'
                              ' of erfinv in CPU cannot be done. ' +
                              str(_import_error))
        self.retain_outputs((0,))
        return utils.force_array(special.erfinv(x[0]), dtype=x[0].dtype), 

def label(self):
        return 'erfinv' 

def _erfinv_cpu(x, dtype):
    from scipy import special
    return numpy.vectorize(special.erfinv, otypes=[dtype])(x) 

def xi_hat(self, xi_bar_hat):
		"""
		A priori SNR estimate.

		Argument/s:
			xi_bar_hat - mapped a priori SNR estimate.

		Returns:
			A priori SNR estimate.
		"""
		xi_db_hat = np.add(np.multiply(np.multiply(self.sigma, np.sqrt(2.0)),
			spsp.erfinv(np.subtract(np.multiply(2.0, xi_bar_hat), 1))), self.mu)
		return np.power(10.0, np.divide(xi_db_hat, 10.0))