quantBisect {openEBGM}

Find quantiles of the posterior distribution

Description

quantBisect finds the desired quantile of the posterior distribution using the bisection method. Used to create credibility limits.

Usage

quantBisect(
  percent,
  theta_hat,
  N,
  E,
  qn,
  digits = 2,
  limits = c(-1e+05, 1e+05),
  max_iter = 2000
)

Arguments

percent

A numeric scalar between 1 and 99 for the desired percentile (e.g., 5 for 5th percentile).

theta_hat

A numeric vector of hyperparameter estimates (likely from autoHyper or from directly minimizing negLLsquash) ordered as: \alpha_1, \beta_1, \alpha_2, \beta_2, P.

N

A whole number vector of actual counts from processRaw.

E

A numeric vector of expected counts from processRaw.

qn

A numeric vector of posterior probabilities that \lambda came from the first component of the mixture, given N = n (i.e., the mixture fraction). See function Qn.

digits

A scalar whole number that determines the number of decimal places used when rounding the results.

limits

A whole number vector of length 2 for the upper and lower bounds of the search space.

max_iter

A whole number scalar for the maximum number of iterations. Used to prevent infinite loops.

Details

The hyperparameter estimates (theta_hat) are:

• \alpha_1, \beta_1: Parameter estimates of the first component of the prior distribution

• \alpha_2, \beta_2: Parameter estimates of the second component

• P: Mixture fraction estimate of the prior distribution

Although this function can find any quantile of the posterior distribution, it will often be used to calculate the 5th and 95th percentiles to create a 90% credibility interval.

The quantile is calculated by solving for x in the general equation F(x) = cutoff, or equivalently, F(x) - cutoff = 0, where F(x) is the cumulative distribution function of the posterior distribution and cutoff is the appropriate cutoff level (e.g., 0.05 for the 5th percentile).

Value

A numeric vector of quantile estimates.

Warning

The digits argument determines the tolerance for the bisection algorithm. The more decimal places you want returned, the longer the run time.

See Also

https://en.wikipedia.org/wiki/Bisection_method

autoHyper, exploreHypers, negLLsquash, negLL, negLLzero, and negLLzeroSquash for hyperparameter estimation.

processRaw for finding counts.

Qn for finding mixture fractions.

Other posterior distribution functions: Qn(), ebgm()

Examples

data.table::setDTthreads(2)  #only needed for CRAN checks
theta_init <- data.frame(
  alpha1 = c(0.5, 1),
  beta1  = c(0.5, 1),
  alpha2 = c(2,   3),
  beta2  = c(2,   3),
  p      = c(0.1, 0.2)
)
data(caers)
proc <- processRaw(caers)
squashed <- squashData(proc, bin_size = 300, keep_pts = 10)
squashed <- squashData(squashed, count = 2, bin_size = 13, keep_pts = 10)
theta_hat <- autoHyper(data = squashed, theta_init = theta_init)$estimates
qn <- Qn(theta_hat, N = proc$N, E = proc$E)
proc$QUANT_05 <- quantBisect(percent = 5, theta = theta_hat, N = proc$N,
                             E = proc$E, qn = qn)
## Not run: proc$QUANT_95 <- quantBisect(percent = 95, theta = theta_hat,
                                      N = proc$N, E = proc$E, qn = qn)
## End(Not run)
head(proc)