Returns the minimal-width confidence interval [mu_low, mu_high] of
    confidence level alpha for a posterior distribution post on the parameter mu.
325343
def confidence_interval_min_width(mu, post, alpha=0.9):
    '''
    
    '''
    if not 0 < alpha < 1:
        raise ValueError("Confidence level must be in (0,1).")

    # choose a step size for the sliding confidence window
    alpha_step = 0.01

    # initialize the lower and upper limits
    mu_low = numpy.min(mu)
    mu_high = numpy.max(mu)

    # find the smallest window (by delta-mu) stepping by dalpha
    for ai in numpy.arange(0, 1 - alpha, alpha_step):
        ml = compute_lower_limit(mu, post, 1 - ai)
        mh = compute_upper_limit(mu, post, alpha + ai)
        if mh - ml < mu_high - mu_low:
            mu_low = ml
            mu_high = mh

    return mu_low, mu_high
1167254
def vp(self):
        
        return np.sqrt(self.mu * (2 / (self.rp) - 1 / self.kep.a))
566277
def get_pvalue(value, ci):
    
    from scipy.stats import norm
    se = (ci[1] - ci[0]) / (2.0 * 1.96)
    z = value / se
    pvalue = -2 * norm.cdf(-np.abs(z))
    return pvalue
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
518263
def _compute_stddevs(self, C, mag, rjb, imt, stddev_types):
        
        # aleatory uncertainty
        sigma_ale_m = np.interp(mag, [5.0, 5.5, 8.0],
                                [C['m50'], C['m55'], C['m80']])
        sigma_ale_rjb = np.interp(rjb, [5.0, 20.0], [C['r5'], C['r20']])
        sigma_ale = np.sqrt(sigma_ale_m ** 2 + sigma_ale_rjb ** 2)

        # epistemic uncertainty
        if imt.period < 1:
            sigma_epi = 0.36 + 0.07 * (mag - 6)
        else:
            sigma_epi = 0.34 + 0.06 * (mag - 6)

        sigma_total = np.sqrt(sigma_ale ** 2 + sigma_epi ** 2)

        stddevs = []
        for _ in stddev_types:
            stddevs.append(sigma_total)

        return stddevs
1023042
def calc_min_interval(x, alpha):
    

    n = len(x)
    cred_mass = 1.0 - alpha

    interval_idx_inc = int(np.floor(cred_mass * n))
    n_intervals = n - interval_idx_inc
    interval_width = x[interval_idx_inc:] - x[:n_intervals]

    if len(interval_width) == 0:
        print_('Too few elements for interval calculation')
        return [None, None]

    min_idx = np.argmin(interval_width)
    hdi_min = x[min_idx]
    hdi_max = x[min_idx + interval_idx_inc]
    return [hdi_min, hdi_max]
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
383341
def confidence_interval_hazard_(self):
        
        return self._compute_confidence_bounds_of_transform(self._hazard, self.alpha, self._ci_labels)
1045021
def rskew_normal(mu, tau, alpha, size=()):
    
    size_ = size or (1,)
    len_ = np.prod(size_)
    return flib.rskewnorm(
        len_, mu, tau, alpha, np.random.normal(size=2 * len_)).reshape(size)
265147
def mle_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"]])

        # calc posterior expectation
        upgrade = self.mle_covariance * self.jco.T * 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[islog,:] = 10.0**post_expt.loc[islog,:]
        return post_expt
667191
def LOO(self, kern, X, Y, likelihood, posterior, Y_metadata=None, K=None):
        
        g = posterior.woodbury_vector
        c = posterior.woodbury_inv
        c_diag = np.diag(c)[:, None]
        neg_log_marginal_LOO = 0.5*np.log(2*np.pi) - 0.5*np.log(c_diag) + 0.5*(g**2)/c_diag
        #believe from Predictive Approaches for Choosing Hyperparameters in Gaussian Processes
        #this is the negative marginal LOO
        return -neg_log_marginal_LOO
1313772
def rlaplace(mu, tau, size=None):
    

    u = np.random.uniform(-0.5, 0.5, size)
    return mu - np.sign(u) * np.log(1 - 2 * np.abs(u)) / tau
854337
def probability_lt(self, x):
        
        if self.mean is None:
            return
        return normdist(x=x, mu=self.mean, sigma=self.standard_deviation)
156937
def minimum(rule, artifact):
    
    m = 1.0
    for i in range(len(rule.R)):
        e = rule.R[i](artifact)
        if e is not None:
            if e < m:
                m = e
    return m
678885
def skew_normal_expval(mu, tau, alpha):
    
    delta = alpha / np.sqrt(1. + alpha ** 2)
    return mu + np.sqrt(2 / pi / tau) * delta
828845
def confidence_interval_dichotomous(
        point_estimate,
        sample_size,
        confidence=.95,
        bias=False,
        percentage=True,
        **kwargs
):
    """"""
    alpha = ppf((confidence + 1) / 2, sample_size - 1)
    p = point_estimate
    if percentage:
        p /= 100

    margin = sqrt(p * (1 - p) / sample_size)
    if bias:
        margin += .5 / sample_size
    if percentage:
        margin *= 100

    return (point_estimate - alpha * margin, point_estimate + alpha * margin)
242822
def _pval_from_bootci(boot, estimate):
    
    if estimate == 0:
        out = 1
    else:
        out = 2 * min(sum(boot > 0), sum(boot < 0)) / len(boot)
    return min(out, 1)