Python math.erf() Examples

def sum_of_erf(mu, sigma, N=1000):
    """
    Compute the sum of erf term

    :param mu: The Gaussian mean
    :param sigma: The Gaussian sigma
    :param N: The number of iterations in the sum
    """
    from math import sqrt, erf

    sum1 = 0
    sum2 = N * erf((N + 1 - mu) / (sqrt(2) * sigma))
    sum3 = N * erf((-N - mu) / (sqrt(2) * sigma))
    for i in range(-N, N):
        sum1 += erf((i + 1 - mu) / (sqrt(2) * sigma))
    return -0.5 * (sum1 + sum2 + sum3) 

def expimp(self, x):
        '''
        Returns the expected improvement at the design vector X in the model
        :param x: A real world coordinates design vector
        :return EI: The expected improvement value at the point x in the model
        '''
        S = self.predicterr_normalized(x)
        y_min = np.min(self.y)
        if S <= 0.:
            EI = 0.
        elif S > 0.:
            EI_one = ((y_min - self.predict_normalized(x)) * (0.5 + 0.5*m.erf((
                      1./np.sqrt(2.))*((y_min - self.predict_normalized(x)) /
                                       S))))
            EI_two = ((S * (1. / np.sqrt(2. * np.pi))) * (np.exp(-(1./2.) *
                      ((y_min - self.predict_normalized(x))**2. / S**2.))))
            EI = EI_one + EI_two
        return EI 

def integral_gaussian(a, b, mu, sigma):
    """
    Computes the integral of a gaussian centered
    in mu with a given sigma
    :param a:
    :param b:
    :param mu:
    :param sigma:
    :return:
    """

    # Integral from -\infty to a
    val_floor = 0.5 * (1 + math.erf((a - mu) / (sigma * math.sqrt(2.0))))

    # Integral from -\infty to b
    val_ceil = 0.5 * (1 + math.erf((b - mu) / (sigma * sqrt(2.0))))

    return val_ceil - val_floor 

def expimp(self, x):
        '''
        Returns the expected improvement at the design vector X in the model
        :param x: A real world coordinates design vector
        :return EI: The expected improvement value at the point x in the model
        '''
        S = self.predicterr_normalized(x)
        y_min = np.min(self.y)
        if S <= 0.:
            EI = 0.
        elif S > 0.:
            EI_one = ((y_min - self.predict_normalized(x)) * (0.5 + 0.5*m.erf((
                      1./np.sqrt(2.))*((y_min - self.predict_normalized(x)) /
                                       S))))
            EI_two = ((S * (1. / np.sqrt(2. * np.pi))) * (np.exp(-(1./2.) *
                      ((y_min - self.predict_normalized(x))**2. / S**2.))))
            EI = EI_one + EI_two
        return EI 

def ndtr(a):
    """
    Returns the area under the Gaussian probability density function,
    integrated from minus infinity to x.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ndtr.html#scipy.special.ndtr

    """
    sqrth = math.sqrt(2) / 2
    x = float(a) * sqrth
    z = abs(x)
    if z < sqrth:
        y = 0.5 + 0.5 * math.erf(x)
    else:
        y = 0.5 * math.erfc(z)
        if x > 0:
            y = 1 - y
    return y 

def erfinv(y):
    # courtesy of
    # https://github.com/antelopeusersgroup/antelope_contrib/blob/master/lib/location/libgenloc/erfinv.c#L15
    a = [0.886226899, -1.645349621, 0.914624893, -0.140543331]
    b = [-2.118377725, 1.442710462, -0.329097515, 0.012229801]
    c = [-1.970840454, -1.624906493, 3.429567803, 1.641345311]
    d = [3.543889200, 1.637067800]
    if y > 1:
        return None
    if y == 1:
        return math.copysign(1, y) * float("inf")
    if y <= 0.7:
        z = y * y
        num = (((a[3] * z + a[2]) * z + a[1]) * z + a[0])
        dem = ((((b[3] * z + b[2]) * z + b[1]) * z + b[0]) * z + 1.0)
        x = y * num / dem
    else:
        z = math.sqrt(-math.log((1.0 - abs(y)) / 2.0))
        num = ((c[3] * z + c[2]) * z + c[1]) * z + c[0]
        dem = (d[1] * z + d[0]) * z + 1.0
        x = (math.copysign(1.0, y)) * num / dem
    for _ in [1, 2]:
        x = x - (erf(x) - y) / ((2 / math.sqrt(trig.c_pi)) * math.exp(-x * x))
    return x 

def integral_default_gaussian(y, x):
    """Calculate the value of the integral from x to y of the erf function

    Parameters
    ----------
    y : float
        upper limit
    x : float
        lower limit

    Returns
    -------
    float
        Calculated value of integral

    """
    return 0.5 * (math.erf(x) - math.erf(y)) 

def erf_local(x):
    # save the sign of x
    sign = 1 if x >= 0 else -1
    x = abs(x)

    # constants
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    # A&S formula 7.1.26
    t = 1.0/(1.0 + p*x)
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)
    return sign*y # erf(-x) = -erf(x) 

def phi(x, m, s):
        """ Cumulative Probability """
        return 1. / 2 * (1 + math.erf((x - m) / s / math.sqrt(2))) 

def erf(x=('FloatPin', 0.0)):
        '''Return the error function at `x`.'''
        return math.erf(x) 

def Φ(x, µ, σ):
    """ Cumulative Probability """
    return 1 / 2 * (1 + math.erf((x - µ) / σ / math.sqrt(2))) 

def normal_cdf(x, mu=0, sigma=1.0):
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2 

def normal_cdf(x, mean = 0, std = 1):
    return 1/2*(1+ math.erf((x-mean)/(std*np.sqrt(2)))) 

def StandardNormalCdf(x):
    """Evaluates the CDF of the standard Normal distribution.
    
    See http://en.wikipedia.org/wiki/Normal_distribution
    #Cumulative_distribution_function

    Args:
        x: float
                
    Returns:
        float
    """
    return (math.erf(x / ROOT2) + 1) / 2 

def phi(n):
    """
    Computes the cumulative distribution function of the standard,
    normal distribution for values <= n.

    >>> round(phi(0), 3)
    0.5
    >>> round(phi(-1), 3)
    0.159
    >>> round(phi(+1), 3)
    0.841

    """
    return (1+erf(n/sqrt(2)))/2 

def normal_cdf(x: float, mu: float = 0, sigma: float = 1) -> float:
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2 

def normal_cdf(x, mu=0,sigma=1):
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2 

def normal_cdf(x, mu=0,sigma=1):
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2 

def qv_success(self):
        """Return whether each depth was successful (>2/3 with confidence
        greater than 97.5) and the confidence

        Returns:
            list: List of lenth depth with eact element a 3 list with
            - success True/False
            - confidence
        """

        success_list = []

        for depth_ind, _ in enumerate(self._depths):
            success_list.append([False, 0.0])
            hmean = self._ydata[0][depth_ind]
            if hmean > 2/3:
                cfd = 0.5 * (1 +
                             math.erf((hmean - 2/3)
                                      / (1e-10 +
                                         self._ydata[1][depth_ind])/2**0.5))
                success_list[-1][1] = cfd

                if cfd > 0.975:
                    success_list[-1][0] = True

        return success_list 

def Φ(x, µ, σ):
    """ Cumulative Probability """
    return 1 / 2 * (1 + math.erf((x - µ) / σ / math.sqrt(2))) 

def erf(*args, **kwargs):
            import mpmath
            return mpmath.erf(*args, **kwargs) 

def phi(x, m, s):
    """ Cumulative Probability """
    return 1. / 2 * (1 + math.erf((x - m) / s / math.sqrt(2))) 

def phi(x, m, s):
    """ Cumulative Probability """
    return 1. / 2 * (1 + math.erf((x - m) / s / math.sqrt(2))) 

def testHalfBasic(self):
    x = np.arange(-3, 3).reshape(1, 3, 2).astype(np.float16)
    y = (x + .5).astype(np.float16)    # no zero
    z = (x + 15.5).astype(np.float16)  # all positive
    self._compareBoth(x, np.abs, tf.abs)
    self._compareBoth(x, np.abs, _ABS)
    self._compareBoth(x, np.negative, tf.neg)
    self._compareBoth(x, np.negative, _NEG)
    self._compareBoth(y, self._inv, tf.inv)
    self._compareBoth(x, np.square, tf.square)
    self._compareBoth(z, np.sqrt, tf.sqrt)
    self._compareBoth(z, self._rsqrt, tf.rsqrt)
    self._compareBoth(x, np.exp, tf.exp)
    self._compareBoth(z, np.log, tf.log)
    self._compareBoth(z, np.log1p, tf.log1p)
    self._compareBoth(x, np.tanh, tf.tanh)
    self._compareBoth(x, self._sigmoid, tf.sigmoid)
    self._compareBoth(y, np.sign, tf.sign)
    self._compareBoth(x, np.sin, tf.sin)
    self._compareBoth(x, np.cos, tf.cos)
    self._compareBoth(
        y,
        np.vectorize(self._replace_domain_error_with_inf(math.lgamma)),
        tf.lgamma)
    self._compareBoth(x, np.vectorize(math.erf), tf.erf)
    self._compareBoth(x, np.vectorize(math.erfc), tf.erfc)

    self._compareBothSparse(x, np.abs, tf.abs)
    self._compareBothSparse(x, np.negative, tf.neg)
    self._compareBothSparse(x, np.square, tf.square)
    self._compareBothSparse(z, np.sqrt, tf.sqrt, tol=1e-3)
    self._compareBothSparse(x, np.tanh, tf.tanh)
    self._compareBothSparse(y, np.sign, tf.sign)
    self._compareBothSparse(x, np.vectorize(math.erf), tf.erf, tol=1e-3) 

def testFloatEmpty(self):
    x = np.empty((2, 0, 5), dtype=np.float32)
    self._compareBoth(x, np.abs, tf.abs)
    self._compareBoth(x, np.abs, _ABS)
    self._compareBoth(x, np.negative, tf.neg)
    self._compareBoth(x, np.negative, _NEG)
    self._compareBoth(x, self._inv, tf.inv)
    self._compareBoth(x, np.square, tf.square)
    self._compareBoth(x, np.sqrt, tf.sqrt)
    self._compareBoth(x, self._rsqrt, tf.rsqrt)
    self._compareBoth(x, np.exp, tf.exp)
    self._compareBoth(x, np.log, tf.log)
    self._compareBoth(x, np.log1p, tf.log1p)
    self._compareBoth(x, np.tanh, tf.tanh)
    self._compareBoth(x, self._sigmoid, tf.sigmoid)
    self._compareBoth(x, np.sign, tf.sign)
    self._compareBoth(x, np.sin, tf.sin)
    self._compareBoth(x, np.cos, tf.cos)
    # Can't use vectorize below, so just use some arbitrary function
    self._compareBoth(x, np.sign, tf.lgamma)
    self._compareBoth(x, np.sign, tf.erf)
    self._compareBoth(x, np.sign, tf.erfc)
    self._compareBoth(x, np.tan, tf.tan)
    self._compareBoth(x, np.arcsin, tf.asin)
    self._compareBoth(x, np.arccos, tf.acos)
    self._compareBoth(x, np.arctan, tf.atan)

    self._compareBothSparse(x, np.abs, tf.abs)
    self._compareBothSparse(x, np.negative, tf.neg)
    self._compareBothSparse(x, np.square, tf.square)
    self._compareBothSparse(x, np.sqrt, tf.sqrt, tol=1e-3)
    self._compareBothSparse(x, np.tanh, tf.tanh)
    self._compareBothSparse(x, np.sign, tf.sign)
    self._compareBothSparse(x, np.sign, tf.erf) 

def WP(RA, RB, VA, VB):
    return (math.erf((RB - RA) / math.sqrt(2 * (VA * VA + VB * VB))) + 1) / 2.0 

def _compute_harmonic_volume(radius, spring_constant, beta):
    """Compute the volume of an harmonic potential from 0 to radius.

    Parameters
    ----------
    radius : simtk.unit.Quantity
        The upper limit on the distance (units of length).
    spring_constant : simtk.unit.Quantity
        The spring constant of the harmonic potential (units of
        energy/mole/length^2).
    beta : simtk.unit.Quantity
        Thermodynamic beta (units of mole/energy).

    Returns
    -------
    volume : simtk.unit.Quantity
        The volume of the harmonic potential (units of length^3).

    """
    # Turn everything to consistent dimension-less units.
    length_unit = unit.nanometers
    energy_unit = unit.kilojoules_per_mole
    radius /= length_unit
    beta *= energy_unit
    spring_constant /= energy_unit/length_unit**2

    bk = beta * spring_constant
    bk_2 = bk / 2
    bkr2_2 = bk_2 * radius**2
    volume = math.sqrt(math.pi/2) * math.erf(math.sqrt(bkr2_2)) / bk**(3.0/2)
    volume -= math.exp(-bkr2_2) * radius / bk
    return 4 * math.pi * volume * length_unit**3 

def test_erf_erfc(self):
        tolerance = 15
        for x in itertools.count(0.0, 0.001):
            if x > 5.0:
                break
            self.assertAlmostEqual(math.erf(x), -math.erf(-x), places=tolerance)
            self.assertAlmostEqual(math.erfc(x), 2.0 - math.erfc(-x), places=tolerance)

            self.assertAlmostEqual(1.0 - math.erf(x), math.erfc(x), places=tolerance)
            self.assertAlmostEqual(1.0 - math.erf(-x), math.erfc(-x), places=tolerance)
            self.assertAlmostEqual(1.0 - math.erfc(x), math.erf(x), places=tolerance)
            self.assertAlmostEqual(1.0 - math.erfc(-x), math.erf(-x), places=tolerance) 

def erf(*args, **kwargs):
            from scipy.special import erf
            return erf(*args, **kwargs)


#    from scipy.special import lambertw, ellipe, gammaincc, gamma # fluids
#    from scipy.special import i1, i0, k1, k0, iv # ht
#    from scipy.special import hyp2f1    
#    if erf is None:
#        from scipy.special import erf 

def _no_grad_trunc_normal_(tensor, mean, std, a, b):
    # Cut & paste from PyTorch official master until it's in a few official releases - RW
    # Method based on https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf
    def norm_cdf(x):
        # Computes standard normal cumulative distribution function
        return (1. + math.erf(x / math.sqrt(2.))) / 2.

    if (mean < a - 2 * std) or (mean > b + 2 * std):
        warnings.warn("mean is more than 2 std from [a, b] in nn.init.trunc_normal_. "
                      "The distribution of values may be incorrect.",
                      stacklevel=2)

    with torch.no_grad():
        # Values are generated by using a truncated uniform distribution and
        # then using the inverse CDF for the normal distribution.
        # Get upper and lower cdf values
        l = norm_cdf((a - mean) / std)
        u = norm_cdf((b - mean) / std)

        # Uniformly fill tensor with values from [l, u], then translate to
        # [2l-1, 2u-1].
        tensor.uniform_(2 * l - 1, 2 * u - 1)

        # Use inverse cdf transform for normal distribution to get truncated
        # standard normal
        tensor.erfinv_()

        # Transform to proper mean, std
        tensor.mul_(std * math.sqrt(2.))
        tensor.add_(mean)

        # Clamp to ensure it's in the proper range
        tensor.clamp_(min=a, max=b)
        return tensor