Description Usage Arguments Details Value References Examples
Does the same as iLap
but in parallel
1 2 |
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" |
control |
A named list of control parameters with elements |
... |
Additional arguments to be passed to |
See iLap
.
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
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: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_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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.