compute_quantiles {aghq}

Quantiles

Description

Compute marginal quantiles using AGHQ. This function works by first approximating the CDF using aghq::compute_pdf_and_cdf and then inverting the approximation numerically.

Usage

compute_quantiles(
  obj,
  q = c(0.025, 0.975),
  transformation = default_transformation(),
  ...
)

## Default S3 method:
compute_quantiles(
  obj,
  q = c(0.025, 0.975),
  transformation = default_transformation(),
  interpolation = "auto",
  ...
)

## S3 method for class 'list'
compute_quantiles(
  obj,
  q = c(0.025, 0.975),
  transformation = default_transformation(),
  ...
)

## S3 method for class 'aghq'
compute_quantiles(
  obj,
  q = c(0.025, 0.975),
  transformation = obj$transformation,
  ...
)

Arguments

obj

Either the output of aghq::aghq(), its list of marginal distributions (element marginals), or an individual data.frame containing one of these marginal distributions as output by aghq::marginal_posterior().

q

Numeric vector of values in (0,1). The quantiles to compute.

transformation

Optional. Calculate marginal quantiles for a transformation of the parameter whose posterior was normalized using adaptive quadrature. transformation is either: a) an aghqtrans object returned by aghq::make_transformation, or b) a list that will be passed to that function internally. See ?aghq::make_transformation for details. Note that since g has to be monotone anyways, this just returns sort(g(q)) instead of q.

...

Used to pass additional arguments.

interpolation

Which method to use for interpolating the marginal posterior, 'polynomial' (default) or 'spline'? If k > 3 then the polynomial may be unstable and you should use the spline, but the spline doesn't work unless k > 3 so it's not the default. See interpolate_marginal_posterior().

Value

A named numeric vector containing the quantiles you asked for, for the variable whose marginal posterior you provided.

See Also

Other summaries: compute_pdf_and_cdf(), interpolate_marginal_posterior(), marginal_posterior()

Examples

logfteta2d <- function(eta,y) {
  # eta is now (eta1,eta2)
  # y is now (y1,y2)
  n <- length(y)
  n1 <- ceiling(n/2)
  n2 <- floor(n/2)
  y1 <- y[1:n1]
  y2 <- y[(n1+1):(n1+n2)]
  eta1 <- eta[1]
  eta2 <- eta[2]
  sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
    sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
  fn = objfunc2d,
  gr = function(x) numDeriv::grad(objfunc2d,x),
  he = function(x) numDeriv::hessian(objfunc2d,x)
)
opt_sparsetrust_2d <- optimize_theta(funlist2d,c(1.5,1.5))
margpost <- marginal_posterior(opt_sparsetrust_2d,3,1) # margpost for theta1
etaquant <- compute_quantiles(margpost)
etaquant
# lambda = exp(eta)
exp(etaquant)
# Compare to truth
qgamma(.025,1+sum(y1),1+n1)
qgamma(.975,1+sum(y1),1+n1)