iLaplace (version 1.1.0)

iLap_par: Improved Laplace approximation in parallel

Description

Does the same as iLap but in parallel

Usage

iLap_par(fullOpt, ff, ff.gr, ff.hess, quad.data, control = list(n.cores = 1), ...)

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"

control

A named list of control parameters with elements n.cores which sets the number of cores to be used for the parallel computiations. See "Details" for more information.

...

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

Value

A double, the logarithm of the integral

Details

See iLap.

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 <- iLap_par(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)