Robust implementations of log cdf of a standard normal.

     @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]]
     in C from Matthias Seeger.
556688
def logCdfNormal(z):
    
    if (abs(z) < ERF_CODY_LIMIT1):
        # Phi(z) approx  (1+y R_3(y^2))/2, y=z/sqrt(2)
        return np.log1p((z / M_SQRT2) * _erfRationalHelperR3(0.5 * z * z)) - M_LN2
    elif (z < 0.0):
        # Phi(z) approx N(z)Q(-z)/(-z), z<0
        return logPdfNormal(z) - np.log(-z) + np.log(_erfRationalHelper(-z))
    else:
        return np.log1p(-(np.exp(logPdfNormal(z))) * _erfRationalHelper(z) / z)
17899
def normal_log_q_initial(self,z):
                     
        means, scale = self.get_means_and_scales_from_q()
        return ss.norm.logpdf(z,loc=means,scale=scale)
1303097
@Deprecated
  public Counter<L> logProbabilityOf(RVFDatum<L, F> example) {
    // NB: this duplicate method is needed so it calls the scoresOf method
    // with an RVFDatum signature!!  Don't remove it!
    // JLS: type resolution of method parameters is static
    Counter<L> scores = scoresOf(example);
    Counters.logNormalizeInPlace(scores);
    return scores;
  }
686924
public static double cdf(double val, double mu, double sigma, double k) {
    final double x = (val - mu) / sigma;
    if(k > 0 || k < 0) {
      if(k * x > 1) {
        return k > 0 ? 1 : 0;
      }
      return FastMath.exp(-FastMath.exp(FastMath.log(1 - k * x) / k));
    }
    else { // Gumbel case:
      return FastMath.exp(-FastMath.exp(-x));
    }
  }
748837
def _initialize_likelihood_prior(self, positions, log_likelihoods, log_priors):
        
        func = SimpleCLFunction.from_string('''
            void compute(global mot_float_type* chain_position,
                         global mot_float_type* log_likelihood,
                         global mot_float_type* log_prior,
                         local mot_float_type* x_tmp,
                         void* data){
                
                bool is_first_work_item = get_local_id(0) == 0;
    
                if(is_first_work_item){
                    for(uint i = 0; i < ''' + str(self._nmr_params) + '''; i++){
                        x_tmp[i] = chain_position[i];
                    }
                    *log_prior = _computeLogPrior(x_tmp, data);
                }
                barrier(CLK_LOCAL_MEM_FENCE);
                    
                *log_likelihood = _computeLogLikelihood(x_tmp, data);
            }
        ''', dependencies=[self._get_log_prior_cl_func(), self._get_log_likelihood_cl_func()])

        kernel_data = {
            'chain_position': Array(positions, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'log_likelihood': Array(log_likelihoods, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'log_prior': Array(log_priors, 'mot_float_type', mode='rw', ensure_zero_copy=True),
            'x_tmp': LocalMemory('mot_float_type', self._nmr_params),
            'data': self._data
        }

        func.evaluate(kernel_data, self._nmr_problems,
                      use_local_reduction=all(env.is_gpu for env in self._cl_runtime_info.cl_environments),
                      cl_runtime_info=self._cl_runtime_info)
1005572
def _discover_bounds(cdf, tol=1e-7):
    
    class DistFromCDF(stats.distributions.rv_continuous):
        def cdf(self, x):
            return cdf(x)
    dist = DistFromCDF()
    # the ppf is the inverse cdf
    lower = dist.ppf(tol)
    upper = dist.ppf(1. - tol)
    return lower, upper
1234543
private static double logLikelihood(double a, double b, 
                                        double c, double d) {
        // Table set up as:
        //  a   b
        //  c   d
        double col1sum = a + c;
        double col2sum = b + d;
        double row1sum = a + b;
        double row2sum = c + d;
        double sum = row1sum + row2sum;
        
        // Calculate the expected values for a, b, c, d
        double aExp = (row1sum / sum) * col1sum;
        double bExp = (row1sum / sum) * col2sum;
        double cExp = (row2sum / sum) * col1sum;
        double dExp = (row2sum / sum) * col2sum;

        // log(0) = Infinity, which messes up the calcuation.  Therefore, check
        // whether the value is zero before calculating its contribution.
        double aVal = (a == 0) ? 0 : a * Math.log(a / aExp);
        double bVal = (b == 0) ? 0 : b * Math.log(b / bExp);
        double cVal = (c == 0) ? 0 : c * Math.log(c / cExp);
        double dVal = (d == 0) ? 0 : d * Math.log(d / dExp);

        return 2 * (aVal + bVal + cVal + dVal);            
    }
449346
public static double cdf(double x, double mu, double sigma, double skew) {
    x = (x - mu) / sigma;
    if(Math.abs(skew) > 0.) {
      double tmp = 1 - skew * x;
      if(tmp < 1e-15) {
        return (skew < 0.) ? 0. : 1.;
      }
      x = -FastMath.log(tmp) / skew;
    }
    return .5 + .5 * NormalDistribution.erf(x * MathUtil.SQRTHALF);
  }
3989
def MakeCdf(self, steps=101):
        """"""
        xs = [i / (steps - 1.0) for i in xrange(steps)]
        ps = [scipy.special.betainc(self.alpha, self.beta, x) for x in xs]
        cdf = Cdf(xs, ps)
        return cdf
116787
def rtruncated_normal(mu, tau, a=-np.inf, b=np.inf, size=None):
    

    sigma = 1. / np.sqrt(tau)
    na = utils.normcdf((a - mu) / sigma)
    nb = utils.normcdf((b - mu) / sigma)

    # Use the inverse CDF generation method.
    U = np.random.mtrand.uniform(size=size)
    q = U * nb + (1 - U) * na
    R = utils.invcdf(q)

    # Unnormalize
    return R * sigma + mu
1257515
public static double cdf(double val, double loc, double scale) {
    val = (val - loc) / scale;
    return 1. / (1. + FastMath.exp(-val));
  }
40035
public static double normalizeLogProps(double[] logProps) {
        double logPropSum = DoubleArrays.logSum(logProps);
        if (logPropSum == Double.NEGATIVE_INFINITY) {
            double uniform = FastMath.log(1.0 / (double)logProps.length);
            for (int d = 0; d < logProps.length; d++) {
                logProps[d] = uniform;
            }
        } else if (logPropSum == Double.POSITIVE_INFINITY) {
            int count = DoubleArrays.count(logProps, Double.POSITIVE_INFINITY);
            if (count == 0) {
                throw new RuntimeException("Unable to normalize since sum is infinite but contains no infinities: " + Arrays.toString(logProps));
            }
            double constant = FastMath.log(1.0 / (double) count);
            for (int d=0; d<logProps.length; d++) {
                if (logProps[d] == Double.POSITIVE_INFINITY) {
                    logProps[d] = constant;
                } else {
                    logProps[d] = Double.NEGATIVE_INFINITY;
                }
            }
        } else {
            for (int d = 0; d < logProps.length; d++) {
                logProps[d] -= logPropSum;
                assert(!Double.isNaN(logProps[d]));
            }            
        }
        return logPropSum;
    }
717701
def _cdf(self, xloc, dist, base, cache):
        """"""
        return evaluation.evaluate_forward(dist, base**xloc, cache=cache)
580199
public static double cdf(double x, double mu, double sigma, double lambda) {
    if(x == Double.NEGATIVE_INFINITY) {
      return 0.;
    }
    if(x == Double.POSITIVE_INFINITY) {
      return 1.;
    }
    final double u = lambda * (x - mu);
    final double v = lambda * sigma;
    final double v2 = v * v;
    final double logphi = FastMath.log(NormalDistribution.cdf(u, v2, v));
    return NormalDistribution.cdf(u, 0., v) - FastMath.exp(-u + v2 * .5 + logphi);
  }
1150228
def logspec2s2mfc(logspec, ncep=13):
    """"""
    nframes, nfilt = logspec.shape
    melcos = s2dctmat(nfilt, ncep, 1./nfilt)
    return numpy.dot(logspec, melcos.T) / nfilt
151831
public static double logcdf(double val, double k, double theta) {
    if(val < 0) {
      return Double.NEGATIVE_INFINITY;
    }
    double vt = val * theta;
    return (val == Double.POSITIVE_INFINITY) ? 0. : logregularizedGammaP(k, vt);
  }
1331573
function normcdf(x) {
  var a1 = 0.398942280444;
  var a2 = 0.399903438504;
  var a3 = 5.75885480458;
  var a4 = 29.8213557808;
  var a5 = 2.62433121679;
  var a6 = 48.6959930692;
  var a7 = 5.92885724438;
  var b0 = 0.398942280385;
  var b1 = 3.8052E-08;
  var b2 = 1.00000615302;
  var b3 = 3.98064794E-04;
  var b4 = 1.98615381364;
  var b5 = 0.151679116635;
  var b6 = 5.29330324926;
  var b7 = 4.8385912808;
  var b8 = 15.1508972451;
  var b9 = 0.742380924027;
  var b10 = 30.789933034;
  var b11 = 3.99019417011;
  var q;
  var y;
  var absx = Math.abs(x);

  if (absx <= 1.28) {
    y = 0.5 * x * x;
    q = 0.5 - absx * (a1 -
      a2 * y / (y + a3 -
      a4 / (y + a5 +
      a6 / (y + a7))));
  } else if (absx <= 12.7) {
    y = 0.5 * x * x;
    q = Math.exp(-y) * b0 / (absx - b1
      + b2 / (absx + b3
      + b4 / (absx - b5
      + b6 / (absx + b7
      - b8 / (absx + b9
      + b10 / (absx + b11))))));
  } else {
    q = 0.0;
  }

  if (x < 0.0)
  {
    return q;
  }
  return 1.0 - q;
}