library(aster) library(numDeriv) # needed because of the change in R function "sample" in R-devel suppressWarnings(RNGversion("3.5.2")) set.seed(42) nind <- 25 nnode <- 5 ncoef <- nnode + 1 famlist <- fam.default() fam <- c(1, 1, 2, 3, 3) pred <- c(0, 1, 1, 2, 3) modmat <- array(0, c(nind, nnode, ncoef)) modmat[ , , 1] <- 1 for (i in 2:nnode) modmat[ , i, i] <- 1 modmat[ , , ncoef] <- rnorm(nind * nnode) beta <- rnorm(ncoef) / 10 theta <- matrix(modmat, ncol = ncoef) %*% beta theta <- matrix(theta, ncol = nnode) root <- sample(1:3, nind * nnode, replace = TRUE) root <- matrix(root, nind, nnode) x <- raster(theta, pred, fam, root) out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 2, type = "conditional") my.value <- 0 for (j in 1:nnode) { ifam <- fam[j] k <- pred[j] if (k > 0) xpred <- x[ , k] else xpred <- root[ , j] for (i in 1:nind) my.value <- my.value - sum(x[i, j] * theta[i, j] - xpred[i] * famfun(famlist[[ifam]], 0, theta[i, j])) } all.equal(out$value, my.value) foo <- function(beta) { stopifnot(is.numeric(beta)) stopifnot(is.finite(beta)) mlogl(beta, pred, fam, x, root, modmat, type = "conditional")$value } g <- grad(foo, beta) all.equal(g, out$gradient) h <- hessian(foo, beta) all.equal(h, out$hessian) ########## objfun <- function(beta) { out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 1, type = "conditional") result <- out$value attr(result, "gradient") <- out$gradient return(result) } out1 <- nlm(objfun, beta, fscale = nind) ########## objfun <- function(beta) { out <- mlogl(beta, pred, fam, x, root, modmat, deriv = 2, type = "conditional") result <- out$value attr(result, "gradient") <- out$gradient attr(result, "hessian") <- out$hessian return(result) } out <- nlm(objfun, beta, fscale = nind) all.equal(out1$minimum, out$minimum) all.equal(out1$estimate, out$estimate, tolerance = 1e-4) ########## expected Fisher information ########## aster:::setfam(fam.default()) fout <- .C(aster:::C_aster_fish_cond, nind = as.integer(nind), nnode = as.integer(nnode), ncoef = as.integer(ncoef), pred = as.integer(pred), fam = as.integer(fam), beta = as.double(out$estimate), root = as.double(root), x = as.double(x), modmat = as.double(modmat), fish = matrix(as.double(0), ncoef, ncoef)) mout <- mlogl(out$estimate, pred, fam, x, root, modmat, deriv = 2, type = "conditional") aster:::setfam(fam.default()) theta <- .C(aster:::C_aster_mat_vec_mult, nrow = as.integer(nind * nnode), ncol = as.integer(ncoef), a = as.double(modmat), b = as.double(out$estimate), c = matrix(as.double(0), nind, nnode))$c ctau <- .C(aster:::C_aster_theta2ctau, nind = as.integer(nind), nnode = as.integer(nnode), pred = as.integer(pred), fam = as.integer(fam), theta = as.double(theta), ctau = matrix(as.double(0), nind, nnode))$ctau tau <- .C(aster:::C_aster_ctau2tau, nind = as.integer(nind), nnode = as.integer(nnode), pred = as.integer(pred), fam = as.integer(fam), root = as.double(root), ctau = as.double(ctau), tau = matrix(as.double(0), nind, nnode))$tau psiDoublePrime <- .C(aster:::C_aster_theta2whatsis, nind = as.integer(nind), nnode = as.integer(nnode), pred = as.integer(pred), fam = as.integer(fam), deriv = as.integer(2), theta = as.double(theta), result = matrix(as.double(0), nind, nnode))$result my.hessian <- theta * NaN my.fish <- theta * NaN for (i in 1:nind) { for (j in 1:nnode) { k <- pred[j] if (k > 0) { my.hessian[i, j] <- x[i, k] * psiDoublePrime[i, j] my.fish[i, j] <- tau[i, k] * psiDoublePrime[i, j] } else { my.hessian[i, j] <- root[i, j] * psiDoublePrime[i, j] my.fish[i, j] <- root[i, j] * psiDoublePrime[i, j] } } } modmat <- matrix(as.double(modmat), ncol = ncoef) my.hessian <- as.numeric(my.hessian) my.fish <- as.numeric(my.fish) my.hessian <- t(modmat) %*% diag(my.hessian) %*% modmat my.fish <- t(modmat) %*% diag(my.fish) %*% modmat all.equal(my.hessian, mout$hessian) all.equal(my.fish, fout$fish)