Normalize the log-posterior distribution using Adaptive Gauss-Hermite Quadrature. This function takes in a function and its gradient and Hessian, and returns a list of information about the normalized posterior, with methods for summarizing and plotting.
aghq(
ff,
k,
startingvalue,
transformation = default_transformation(),
optresults = NULL,
basegrid = NULL,
control = default_control(),
...
)An object of class aghq which is a list containing elements:
normalized_posterior: The output of the normalize_logpost function, which
itself is a list with elements:
nodesandweights: a dataframe containing the nodes and weights for the adaptive quadrature rule, with the un-normalized and normalized log posterior evaluated at the nodes.
thegrid: a NIGrid object from the mvQuad package, see ?mvQuad::createNIGrid.
lognormconst: the actual result of the quadrature: the log of the normalizing constant of the posterior.
marginals: a list of the same length as startingvalue of which element j
is the result of calling aghq::marginal_posterior with that j. This is
a tbl_df/tbl/data.frame containing the normalized log marginal posterior
for theta_j evaluated at the original quadrature points. Has columns
"thetaj","logpost_normalized","weights", where j is the j you specified.
optresults: information and results from the optimization of the log posterior, the
result of calling aghq::optimize_theta. This a list with elements:
ff: the function list that was provided
mode: the mode of the log posterior
hessian: the hessian of the log posterior at the mode
convergence: specific to the optimizer used, a message indicating whether it converged
control: the control parameters passed
A list with three elements:
fn: function taking argument theta and returning a numeric
value representing the log-posterior at theta
gr: function taking argument theta and returning a numeric
vector representing the gradient of the log-posterior at theta
he: function taking argument theta and returning a numeric
matrix representing the hessian of the log-posterior at theta
The user may wish to use numDeriv::grad and/or numDeriv::hessian to
obtain these. Alternatively, the user may consider the TMB package. This
list is deliberately formatted to match the output of TMB::MakeADFun.
Integer, the number of quadrature points to use. I suggest at least 3. k = 1 corresponds to a Laplace approximation.
Value to start the optimization. ff$fn(startingvalue),
ff$gr(startingvalue), and ff$he(startingvalue) must all return
appropriate values without error.
Optional. Do the quadrature for parameter theta, but
return summaries and plots for parameter g(theta).
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.
Optional. A list of the results of the optimization of the log
posterior, formatted according to the output of aghq::optimize_theta. The
aghq::aghq function handles the optimization for you; passing this list
overrides this, and is useful for when you know your optimization is too difficult to be
handled by general-purpose software. See the software paper for several examples of this.
If you're unsure whether this option is needed for your problem then it probably is not.
Optional. Provide an object of class NIGrid from the mvQuad
package, representing the base quadrature rule that will be adapted. This is only
for users who want more complete control over the quadrature, and is not necessary
if you are fine with the default option which basically corresponds to
mvQuad::createNIGrid(length(theta),'GHe',k,'product'). Note: the mvQuad
functions used within aghq operate on grids in memory, so your basegrid
object will be changed after you run aghq.
A list with elements
method: optimization method to use:
'sparse_trust' (default): trustOptim::trust.optim with method = 'sparse'
'SR1' (default): trustOptim::trust.optim with method = 'SR1'
'trust': trust::trust
'BFGS': optim(...,method = "BFGS")
optcontrol: optional: a list of control parameters to pass to the
internal optimizer you chose. The aghq package uses sensible defaults.
Additional arguments to be passed to ff$fn, ff$gr, and ff$he.
When k = 1 the AGHQ method is a Laplace approximation, and you should use
the aghq::laplace_approximation function, since some of the methods for
aghq objects won't work with only one quadrature point. Objects of
class laplace have different methods suited to this case. See ?aghq::laplace_approximation.
Other quadrature:
get_hessian(),
get_log_normconst(),
get_mode(),
get_nodesandweights(),
get_numquadpoints(),
get_opt_results(),
get_param_dim(),
laplace_approximation(),
marginal_laplace_tmb(),
marginal_laplace(),
nested_quadrature(),
normalize_logpost(),
optimize_theta(),
plot.aghq(),
print.aghqsummary(),
print.aghq(),
print.laplacesummary(),
print.laplace(),
print.marginallaplacesummary(),
summary.aghq(),
summary.laplace(),
summary.marginallaplace()
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)
)
thequadrature <- aghq(funlist2d,3,c(0,0))
Run the code above in your browser using DataLab