Description Usage Arguments Value References Examples
This function is similar to iLap
except that it handles only bivariate integrals of user-written unimodal functions.
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 |
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
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 | # 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)
}
dgf = 5
opt <- nlminb(rep(1,2), ff, gradient = ff.gr, hessian = ff.hess, df = dgf)
opt$hessian <- ff.hess(opt$par, df = dgf);
quad.data = gaussHermiteData(50)
# The improved Laplace approximation (the truth equals 0.0)
iLap <- iLap2d(fullOpt = opt, ff = ff, ff.gr = ff.gr,
ff.hess = ff.hess, quad.data = quad.data,
df = dgf)
# The standard Laplace approximation (the truth equals 0.0)
Lap <- log(2*pi) - opt$objective - 0.5*determinant(opt$hessian)$mod;
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.