iLaplace (version 1.1.0)

iLap2d: Improved Laplace approximation for bivariate integrals of unimodal functions

Description

This function is similar to iLap except that it handles only bivariate integrals of user-written unimodal functions.

Usage

iLap2d(fullOpt, ff, ff.gr, ff.hess, quad.data, ...)

Arguments

fullOpt

A list containing the minium (to be accesed via fullOpt$par), the value of the function at the minimum (to be accessed via fullOpt$objective) and the Hessian matrix at the minimum (to be accessed via fullOpt$hessian

ff

The minus logarithm of the integrand function (the h function, see iLap for further details).

ff.gr

The gradient of ff, having the exact same arguments as ff

ff.hess

The Hessian matrix offf, having the exact same arguments as ff

quad.data

Data for the Gaussian-Herimte quadratures; see "Details"

...

Additional arguments to be passed to ff, ff.gr and ff.hess

Value

a double, the logarithm of the integral

References

Ruli E., Sartori N. & Ventura L. (2015) Improved Laplace approximation for marignal likelihoods. http://arxiv.org/abs/1502.06440

Examples

# The negative integrand function in log
# is the negative log-density of the multivariate
# Student-t density centred at 0 with unit scale matrix
ff <- function(x, df) {
       d <- length(x)
       S <- diag(1, d, d)
       S.inv <- solve(S)
       Q <- colSums((S.inv %*% x) * x)
       logDet <- determinant(S)$modulus
       logPDF <- (lgamma((df + d)/2) - 0.5 * (d * logb(pi * df) +
       logDet) - lgamma(df/2) - 0.5 * (df + d) * logb(1 + Q/df))
       return(-logPDF)
       }

# the gradient of ff
ff.gr <- function(x, df){
            m <- length(x)
            kr = 1 + crossprod(x,x)/df
            return((m+df)*x/(df*kr))
            }

# the Hessian matrix of ff
ff.hess <- function(x, df) {
m <- length(x)
kr <- as.double(1 + crossprod(x,x)/df)
ll <- -(df+m)*2*tcrossprod(x,x)/(df*kr)^2.0
dd = (df+m)*(kr - 2*x^2/df)/(df*kr^2.0)
diag(ll) = dd;
return(ll)
}

dgf = 5
opt <- nlminb(rep(1,2), ff, gradient = ff.gr, hessian = ff.hess, df = dgf)
opt$hessian <- ff.hess(opt$par, df = dgf);
quad.data = gaussHermiteData(50)

# The improved Laplace approximation (the truth equals 0.0)
iLap <- iLap2d(fullOpt = opt, ff = ff, ff.gr = ff.gr,
               ff.hess = ff.hess, quad.data = quad.data,
               df = dgf)
# The standard Laplace approximation (the truth equals 0.0)
Lap <- log(2*pi) - opt$objective - 0.5*determinant(opt$hessian)$mod;