Description Usage Arguments Details Value References Examples
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 |
fullOpt |
A list containing the minium (to be accesed via |
ff |
The minus logarithm of the integrand function (the |
ff.gr |
The gradient of |
ff.hess |
The Hessian matrix of |
quad.data |
Data for the Gaussian-Herimte quadratures; see "Details" |
... |
Additional arguments to be passed to |
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.
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 | # 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(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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.