iLap2d: Improved Laplace approximation for bivariate integrals of...
In iLaplace: Improved Laplace Approximation for Integrals of Unimodal Functions

Description Usage Arguments Value References Examples

View source: R/iLap2d.R

This function is similar to iLap except that it handles only bivariate integrals of user-written unimodal functions.

1	iLap2d(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 `iLap` for further 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`

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

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