iLap: Improved Laplace approximation for integrals of unimodal...

In erlisR/iLaplace: Improved Laplace Approximation for Integrals of Unimodal Functions

Description

This function implements the improved Laplace approximation of Ruli et al. (2015) for multivariate integrals of user-written unimodal functions. See "Details" below for more information.

Usage

iLap(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 "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

Details

iLap approximates integrals of the type

I = \int\exp(-h(x)) dx

where -h() is a concave and unimodal function, with x being d dimensional real vector (d>1). The approximation of I is obtained as the ratio between the unormalised kernel -h(x) and an approximate density function f(x), both evaluated at the modal value x = \hat{x}. The approximate density function f(x) is obtained by resorting to the Laplace approximation for marginal densities. The minimisations are performed with nlminb by suppling the gradient ff.gr and Hessian matrix ff.hess of f(x). The normalisation of the univariate components is perforemd via Gaussian-Hermite quadratures as implemented in the function aghQuad. The Gaussian-Quadrature data, to be provided via the argument quad.data, can be computed with the function gaussHermiteData for a desired number of quadrature points. See "Examples" below.

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

Liu, Q. and Pierce, D. A. (1994). A Note on Gauss-Hermite Quadrature. Biometrika 81, 624-629.

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)
}

df = 5
dims <- 5:15
normConts <- sapply(dims, function(mydim) {
opt <- nlminb(rep(1,mydim), ff, gradient = ff.gr, hessian = ff.hess, df = df)
opt$hessian <- ff.hess(opt$par, df = df);
quad.data = gaussHermiteData(50)
iLap <- iLap(opt, ff, ff.gr, ff.hess, quad.data = quad.data, df = df);
Lap <- mydim*log(2*pi)/2 - opt$objective - 0.5*determinant(opt$hessian)$mod;
return(c(iLap = iLap, Lap = Lap))
})
# plot the results
## Not run: 
plot(dims, normConts[1,], pch="*", cex = 1.6,
 ylim = c(-5, 0)) #improved Laplace
lines(dims, normConts[2,], type = "p", pch = "+") #standard Laplace
abline(h = 0) # the true value

## End(Not run)

## Not run: 
## See also the examples provided in the pacakge iLaplaceExamples, which is
## an auxiliary R pacakge for iLaplace. To download it (be sure you have
## the devtools package) run from R
## devtools::install_github(erlisR/iLaplaceExamples)
## or download the source at \url{https://github.com/erlisR/iLaplaceExamples}.


## End(Not run)

erlisR/iLaplace documentation built on May 16, 2019, 8:48 a.m.