# compute numerical second derivatives involving integrals
num.infoI <- function(coefficients, FUN, X, data,integrate)
{
FUN <- get(FUN, inherits = TRUE)
values <- FUN(coefficients, X, data,integrate)
p <- length(values)
Info <- matrix(0, p, p)
h <- rep(0, p)
delta <- cbind((abs(coefficients) + 1e-012) * 0.0001, rep(1e-012, p))
delta <- apply(delta, 1, max)
for(i in 1:p) {
h[i] <- delta[i]
new.values <- FUN(coefficients + h, X, data,integrate)
Info[, i] <- (new.values - values)/delta[i]
h[i] <- 0
}
Info
}