Nothing
### penmlogl is
###
### - l(a + M alpha + Z D^(1 / 2) c) + (1 / 2) c^T c
###
### in the notation of equation (26) of inst/doc/re.tex
penmlogl <- function(parm, sigma, fixed, random, obj, y, origin, deriv = 2) {
stopifnot(inherits(obj, "aster"))
stopifnot(is.numeric(fixed))
stopifnot(is.finite(fixed))
stopifnot(is.matrix(fixed))
stopifnot(nrow(fixed) == length(obj$x))
stopifnot(is.numeric(sigma))
stopifnot(is.finite(sigma))
stopifnot(is.numeric(parm))
stopifnot(is.finite(parm))
if (! is.list(random))
random <- list(random)
for (i in seq(along = random)) {
if (! is.numeric(random[[i]]))
stop("random or component of random (if list) not numeric")
if (! all(is.finite(random[[i]])))
stop("random or component of random (if list) not all finite")
if (! is.matrix(random[[i]]))
stop("random or component of random (if list) not matrix")
if (nrow(random[[i]]) != length(obj$x))
stop("nrow of random or component of random (if list) not length of response vector")
}
if (length(sigma) != length(random))
stop("length(sigma) != length(random) if list or != 1 otherwise")
if (length(parm) != ncol(fixed) + sum(sapply(random, ncol)))
stop("length(parm) != number of columns of all model matrices, fixed and random")
if (missing(y)) {
y <- obj$x
} else {
stopifnot(is.numeric(y))
stopifnot(is.matrix(y))
stopifnot(identical(dim(y), dim(obj$root)))
stopifnot(all(is.finite(as.vector(y))))
}
if (! missing(origin)) {
stopifnot(is.numeric(origin))
stopifnot(is.matrix(origin))
stopifnot(identical(dim(origin), dim(obj$root)))
stopifnot(all(is.finite(as.vector(origin))))
}
stopifnot(deriv %in% 0:2)
scalevec <- rep(1, ncol(fixed))
penaltyvec <- rep(0, ncol(fixed))
modmat <- fixed
for (i in seq(along = random)) {
foo <- sigma[i] * random[[i]]
modmat <- cbind(modmat, foo)
scalevec <- c(scalevec, rep(sigma[i], ncol(foo)))
penaltyvec <- c(penaltyvec, rep(1, ncol(foo)))
}
### note: despite documentation of the mlogl function, it actually
### works to have modmat a matrix rather than a 3-way array
bar <- mlogl(parm, obj$pred, obj$fam, y, obj$root, modmat, deriv = deriv,
famlist = obj$famlist, origin = origin)
val <- bar$value + sum(penaltyvec * parm^2) / 2
if (deriv == 0)
return(list(value = val, argument = parm, scale = scalevec))
if (deriv >= 1)
grad <- bar$gradient + penaltyvec * parm
if (deriv == 1)
return(list(value = val, gradient = grad, argument = parm,
scale = scalevec, mlogl.gradient = bar$gradient))
## if we get here deriv == 2
hess <- bar$hessian + diag(penaltyvec)
return(list(value = val, gradient = grad, hessian = hess,
argument = parm, scale = scalevec, mlogl.hessian = bar$hessian,
mlogl.gradient = bar$gradient))
}
penmlogl2 <- function(parm, alpha, sigma, fixed, random, obj, y, origin) {
pout <- penmlogl(c(alpha, parm), sigma, fixed, random, obj, y, origin)
nfix <- ncol(fixed)
idx <- seq(along = pout$gradient) > nfix
pout$gradient <- pout$gradient[idx]
tmp <- pout$hessian[idx, , drop = FALSE]
pout$hessian <- tmp[ , idx, drop = FALSE]
pout$scale <- pout$scale[idx]
return(pout)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.