https://github.com/cran/aster
Tip revision: cd7e4fc006dc5296865fa6523ce7d087d86d3ca8 authored by Charles J. Geyer on 20 October 2012, 00:00:00 UTC
version 0.8-20
version 0.8-20
Tip revision: cd7e4fc
quickle.R
### implements (7) of the design doc
quickle <- function(alphanu, bee, fixed, random, obj, y, origin, zwz,
deriv = 0)
{
stopifnot(inherits(obj, "aster"))
if (missing(y)) {
y <- obj$x
} else {
stopifnot(is.matrix(y))
stopifnot(is.numeric(y))
stopifnot(is.finite(y))
stopifnot(dim(y) == dim(obj$x))
}
if (missing(origin)) {
origin <- obj$origin
} else {
stopifnot(is.matrix(origin))
stopifnot(is.numeric(origin))
stopifnot(is.finite(origin))
stopifnot(dim(origin) == dim(obj$origin))
}
stopifnot(is.matrix(fixed))
stopifnot(is.numeric(fixed))
stopifnot(is.finite(fixed))
nfix <- ncol(fixed)
stopifnot(is.matrix(random) | is.list(random))
if (! is.list(random))
random <- list(random)
for (i in seq(along = random)) {
foo <- random[[i]]
if (! is.matrix(foo))
stop("random not matrix or list of matrices")
if (! is.numeric(foo))
stop("random not numeric matrix or list of such")
if (! all(is.finite(foo)))
stop("some random effects model matrix not all finite")
if (nrow(foo) != nrow(fixed))
stop("fixed and random effects model matrices with different nrow")
}
nrand <- sapply(random, ncol)
stopifnot(is.matrix(zwz))
stopifnot(is.numeric(zwz))
stopifnot(is.finite(zwz))
if (any(dim(zwz) != sum(nrand)))
stop("zwz not square matrix with dimension = number of random effects")
stopifnot(length(deriv) == 1)
stopifnot(deriv %in% c(0, 1, 2))
stopifnot(is.vector(alphanu))
stopifnot(is.numeric(alphanu))
stopifnot(is.finite(alphanu))
if (length(alphanu) != nfix + length(nrand))
stop("alphanu wrong length")
stopifnot(is.vector(bee))
stopifnot(is.numeric(bee))
stopifnot(is.finite(bee))
if (length(bee) != sum(nrand))
stop("bee wrong length")
idx <- seq(along = alphanu)
is.alpha <- idx <= nfix
is.nu <- nfix < idx
alpha <- alphanu[is.alpha]
nu <- alphanu[is.nu]
dee <- rep(nu, times = nrand)
if (any(nu < 0))
return(list(value = Inf, alpha = alpha, bee = bee, nu = nu))
modmat <- cbind(fixed, Reduce(cbind, random))
### note: despite documentation of the mlogl function, it actually
### works to have modmat a matrix rather than a 3-way array
mout <- mlogl(c(alpha, bee), obj$pred, obj$fam, y, obj$root, modmat,
deriv = 2, famlist = obj$famlist, origin = origin)
idx <- seq(along = mout$gradient)
is.alpha <- idx <= nfix
is.bee <- nfix < idx
if (missing(origin))
mymlogl <- function(bee)
mlogl(c(alpha, bee), obj$pred, obj$fam, y, obj$root,
modmat, deriv = 2, famlist = obj$famlist)
else
mymlogl <- function(bee)
mlogl(c(alpha, bee), obj$pred, obj$fam, y, obj$root,
modmat, deriv = 2, famlist = obj$famlist, origin = origin)
objfun <- function(bee) {
### note: despite documentation of the mlogl function, it actually
### works to have modmat a matrix rather than a 3-way array
mout <- mymlogl(bee)
val <- mout$value + sum(bee^2 / dee) / 2
grad <- mout$gradient[is.bee] + bee / dee
hess <- mout$hessian
hess <- hess[is.bee, , drop = FALSE]
hess <- hess[ , is.bee, drop = FALSE]
### see note on help page for diag !!!
hess <- hess + diag(1 / dee, nrow = length(dee))
return(list(value = val, gradient = grad, hessian = hess))
}
tout <- trust(objfun, bee, rinit = 1, rmax = 10, iterlim = 1000)
stopifnot(tout$converged)
bee <- tout$argument
mout <- mlogl(c(alpha, bee), obj$pred, obj$fam, y, obj$root, modmat,
deriv = 2, famlist = obj$famlist, origin = origin)
a <- sqrt(dee)
bigh <- zwz * outer(a, a) + diag(length(a))
bigh.chol <- chol(bigh)
val <- mout$value + sum(bee^2 / dee) / 2 + sum(log(diag(bigh.chol)))
if (deriv == 0)
return(list(value = val, alpha = alpha, bee = bee, nu = nu))
pa <- mout$gradient[is.alpha]
pb <- mout$gradient[is.bee] + bee / dee
bigh.inv <- chol2inv(bigh.chol)
idx <- rep(seq(along = nu), times = nrand)
pn <- rep(NaN, length(nu))
for (k in seq(along = nu)) {
eek <- as.numeric(idx == k)
pn[k] <- sum(bigh.inv * zwz * outer(a, eek / a)) / 2 -
sum(bee^2 / dee^2 * eek) / 2
}
if (deriv == 1)
return(list(value = val, gradient = c(pa, pn),
alpha = alpha, bee = bee, nu = nu))
foo <- mout$hessian[is.alpha, , drop = FALSE]
paa <- foo[ , is.alpha, drop = FALSE]
pab <- foo[ , is.bee, drop = FALSE]
foo <- mout$hessian[is.bee, , drop = FALSE]
pbb <- foo[ , is.bee, drop = FALSE] + diag(1 / dee, nrow = length(dee))
pan <- matrix(0, length(alpha), length(nu))
pbn <- matrix(NaN, length(bee), length(nu))
pnn <- matrix(NaN, length(nu), length(nu))
bigh.inverse <- chol2inv(bigh.chol)
for (k in seq(along = nu)) {
eek <- as.numeric(idx == k)
pbn[ , k] <- (- bee * eek / dee^2)
fook <- zwz * outer(a, eek / a)
fook <- bigh.inverse %*% fook
for (m in seq(along = nu)) {
eem <- as.numeric(idx == m)
foom <- zwz * outer(a, eem / a)
foom <- bigh.inverse %*% foom
pnn[k, m] <- sum(bee^2 * eek * eem / dee^3) -
sum(fook * t(foom)) / 2
}
}
pbb.inv <- chol2inv(chol(pbb))
poo <- rbind(cbind(paa, pan), cbind(t(pan), pnn))
pob <- rbind(pab, t(pbn))
hess <- poo - pob %*% pbb.inv %*% t(pob)
return(list(value = val, gradient = c(pa, pn), hessian = hess,
alpha = alpha, bee = bee, nu = nu, pbb.inv = pbb.inv,
pba = t(pab), pbn = pbn))
}