iLaplace (version 1.1.0)

iLapCW: Improved Laplace approximation without nested optimisation

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

iLapCW(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

Value

A double, the logarithm of the integral

Details

iLapCW approximates integrals of the type $$I = \int_{x\in\mathcal{R}^d}\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 in iLap but with nested optimisations replaced by the approximations prposed by Cox & Wermuth (1990).

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.

Cox, D.R and Wermuth, W. (1990). An approximation to maximum likelihood estimates in reduced models. Biometrika 77, 747-761

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:8
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 <- iLapCW(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)