Returns the upper limit mu_high of confidence level alpha for a
    posterior distribution post on the given parameter mu.
    The posterior need not be normalized.
965185
def compute_upper_limit(mu_in, post, alpha=0.9):
    
    if 0 < alpha < 1:
        dp = integral_element(mu_in, post)
        high_idx = bisect.bisect_left(dp.cumsum() / dp.sum(), alpha)
        # if alpha is in (0,1] and post is non-negative, bisect_left
        # will always return an index in the range of mu since
        # post.cumsum()/post.sum() will always begin at 0 and end at 1
        mu_high = mu_in[high_idx]
    elif alpha == 1:
        mu_high = numpy.max(mu_in[post > 0])
    else:
        raise ValueError("Confidence level must be in (0,1].")

    return mu_high
648627
def normal_prior(value, mean, sigma):
    
    return -0.5 * (2 * np.pi * sigma) - (value - mean) ** 2 / (2.0 * sigma)
1153674
def _maximization_step(X, posteriors):
    

    ### Prior probabilities or class weights
    sum_post_proba = np.sum(posteriors, axis=0)
    prior_proba = sum_post_proba / (sum_post_proba.sum() + Epsilon)
    
    ### means
    means = np.dot(posteriors.T, X) / (sum_post_proba[:, np.newaxis] + Epsilon)

    ### covariance matrices
    n_components = posteriors.shape[1]
    n_features = X.shape[1]
    covars = np.empty(shape=(n_components, n_features, n_features), dtype=float)
  
    for i in range(n_components):
        post_i = posteriors[:, i]
        mean_i = means[i]
        diff_i = X - mean_i

        with np.errstate(under='ignore'):
            covar_i = np.dot(post_i * diff_i.T, diff_i) / (post_i.sum() + Epsilon)
        covars[i] = covar_i + Lambda * np.eye(n_features)


    _validate_params(prior_proba, means, covars)
    return(prior_proba, means, covars)
1190905
def lnprior(pars):
    
    # Limit norm and B to be positive
    logprob = (
        naima.uniform_prior(pars[0], 0.0, np.inf)
        + naima.uniform_prior(pars[1], -1, 5)
        + naima.uniform_prior(pars[3], 0, np.inf)
    )

    return logprob
1355030
def alpha_abs(self, R, amp, sigma):
        
        sigma_x, sigma_y = sigma, sigma
        amp_density = self._amp2d_to_3d(amp, sigma_x, sigma_y)
        alpha = self.mass_2d(R, amp_density, sigma) / np.pi / R
        return alpha
335624
def bayesianUpperLimit2(self, alpha, steps=1.e5, plot=False):
        
        cut = ((self.y / 2.) > -30.) # Numeric stability
        try:
            f = scipy.interpolate.interp1d(self.x[cut], self.y[cut], kind='cubic')
        except:
            f = scipy.interpolate.interp1d(self.x[cut], self.y[cut], kind='linear')
        x = np.linspace(0., np.max(self.x[cut]), steps)
        y = np.exp(f(x) / 2.)
        #forbidden = np.nonzero((y / np.exp(self.vertex_y / 2.)) < 1.e-10)[0]
        forbidden = np.nonzero((y / self.vertex_y) < 1.e-10)[0]
        if len(forbidden) > 0:
            index = forbidden[0] # Numeric stability
            x = x[0: index]
            y = y[0: index]
        cdf = np.cumsum(y)
        cdf /= cdf[-1]
        cdf_reflect = scipy.interpolate.interp1d(cdf, x)

        return cdf_reflect(alpha)
604810
def fit_mle(self, data, init_vals=(80, 80)):
        

        if len(data) > 1:
            mu = np.mean(data)
            var = np.var(data)
            theta0 = var / mu
            alpha0 = mu / theta0
        else:
            alpha0 = init_vals[0]
            theta0 = init_vals[1]

        def mle(params):
            return -np.sum(np.log(self.pmf(data, params[0], params[1])))

        # Bounded fmin?
        alpha, theta = optim.fmin(mle, x0=[alpha0, theta0], disp=0)

        return alpha, theta
275112
def EvalGaussianPdf(x, mu, sigma):
    
    return scipy.stats.norm.pdf(x, mu, sigma)
1085589
private double computeConfidence(int support, int samples) {
    final double z = NormalDistribution.standardNormalQuantile(alpha);
    final double eprob = support / (double) samples;
    return Math.max(0., eprob - z * FastMath.sqrt((eprob * (1 - eprob)) / samples));
  }
743123
def confidenceInterval(self, alpha=0.6827, steps=1.e5, plot=False):
        
        x_dense, y_dense = self.densify()
        y_dense -= np.max(y_dense) # Numeric stability
        f = scipy.interpolate.interp1d(x_dense, y_dense, kind='linear')
        x = np.linspace(0., np.max(x_dense), steps)
        # ADW: Why does this start at 0, which often outside the input range?
        # Wouldn't starting at xmin be better:
        #x = np.linspace(np.min(x_dense), np.max(x_dense), steps)
        pdf = np.exp(f(x) / 2.)
        cut = (pdf / np.max(pdf)) > 1.e-10
        x = x[cut]
        pdf = pdf[cut]

        sorted_pdf_indices = np.argsort(pdf)[::-1] # Indices of PDF in descending value
        cdf = np.cumsum(pdf[sorted_pdf_indices])
        cdf /= cdf[-1]
        sorted_pdf_index_max = np.argmin((cdf - alpha)**2)
        x_select = x[sorted_pdf_indices[0: sorted_pdf_index_max]]
            
        return np.min(x_select), np.max(x_select)
1363963
def norm(x, mu, sigma=1.0):
    
    return stats.norm(loc=mu, scale=sigma).pdf(x)
260418
def hpd_threshold(mu_in, post, alpha, tol):
    '''
    
    '''
    norm_post = normalize_pdf(mu_in, post)
    # initialize bisection search
    p_minus = 0.0
    p_plus = max(post)
    while abs(hpd_coverage(mu_in, norm_post, p_minus) -
              hpd_coverage(mu_in, norm_post, p_plus)) >= tol:
        p_test = (p_minus + p_plus) / 2.
        if hpd_coverage(mu_in, post, p_test) >= alpha:
            # test value was too low or just right
            p_minus = p_test
        else:
            # test value was too high
            p_plus = p_test
    # p_minus never goes above the required threshold and p_plus never goes below
    # thus on exiting p_minus is at or below the required threshold and the
    # difference in coverage is within tolerance
    return p_minus
397403
private static long specialNPrimeB(final long numSamplesI, final double p, final double delta) {
    final double q, numSamplesF, oneMinusDelta;
    double tot, curTerm;
    long m;
    assertTrue(numSamplesI >= 1);
    assertTrue((0.0 < p) && (p < 1.0));
    assertTrue((0.0 < delta) && (delta < 1.0));
    q = 1.0 - p;
    oneMinusDelta = 1.0 - delta;
    numSamplesF = numSamplesI;
    curTerm = Math.pow(p, numSamplesF);  // curTerm = posteriorProbability (k, k, p)
    assertTrue(curTerm > 1e-100); // sanity check for non-use of logarithms
    tot = curTerm;
    m = numSamplesI;
    while (tot < oneMinusDelta) {
      curTerm = (curTerm * q * (m)) / ((m + 1) - numSamplesI);
      tot += curTerm;
      m += 1;
    }
    return (m); // don't need to back up
  }
941429
def q_prior(q, m=1, gamma=0.3, qmin=0.1):
    
    if q < qmin or q > 1:
        return 0
    C = 1/(1/(gamma+1)*(1 - qmin**(gamma+1)))
    return C*q**gamma
1167254
def vp(self):
        
        return np.sqrt(self.mu * (2 / (self.rp) - 1 / self.kep.a))
1199316
def _get_SRF_phi(self, imt_per):
        
        if imt_per < 0.6:
            srf = 0.8
        elif 0.6 <= imt_per < 1:
            srf = self._interp_function(0.7, 0.8, 1, 0.6, imt_per)
        elif 1 <= imt_per <= 10:
            srf = self._interp_function(0.6, 0.7, 10, 1, imt_per)
        else:
            srf = 1

        return srf
955951
def map_parameter_estimate(self):
        
        res = self.pst.res
        assert res is not None
        # build the prior expectation parameter vector
        prior_expt = self.pst.parameter_data.loc[:,["parval1"]].copy()
        islog = self.pst.parameter_data.partrans == "log"
        prior_expt.loc[islog] = prior_expt.loc[islog].apply(np.log10)
        prior_expt = Matrix.from_dataframe(prior_expt)
        prior_expt.col_names = ["prior_expt"]
        # build the residual vector
        res_vec = Matrix.from_dataframe(res.loc[:,["residual"]])

        # form the terms of Schur's complement
        b = self.parcov * self.jco.T
        c = ((self.jco * self.parcov * self.jco.T) + self.obscov).inv
        bc = Matrix((b * c).x, row_names=b.row_names, col_names=c.col_names)

        # calc posterior expectation
        upgrade = bc * res_vec
        upgrade.col_names = ["prior_expt"]
        post_expt = prior_expt + upgrade

        # post processing - back log transform
        post_expt = pd.DataFrame(data=post_expt.x,index=post_expt.row_names,
                                 columns=["post_expt"])
        post_expt.loc[:,"prior_expt"] = prior_expt.x.flatten()
        post_expt.loc[islog,:] = 10.0**post_expt.loc[islog,:]
        # unc_sum = self.get_parameter_summary()
        # post_expt.loc[:,"standard_deviation"] = unc_sum.post_var.apply(np.sqrt)
        post_expt.sort_index(inplace=True)
        return post_expt