iLapCW: Improved Laplace approximation without nested optimisation
In iLaplace: Improved Laplace Approximation for Integrals of Unimodal Functions

Description Usage Arguments Details Value References Examples

View source: R/iLapCW.R

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.

1	iLapCW(fullOpt, ff, ff.gr, ff.hess, quad.data, ...)

`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 of`ff`, 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`

iLapCW 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 in iLap but with nested optimisations replaced by the approximations prposed by Cox & Wermuth (1990).

A double, the logarithm of the integral

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

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