https://github.com/cran/aster
Tip revision: aa47935123bfca8a22cbc8345d658d0c1713a289 authored by Charles J. Geyer on 14 December 2023, 15:20:02 UTC
version 1.1-3
version 1.1-3
Tip revision: aa47935
mlogl.R
mlogl <- function(parm, pred, fam, x, root, modmat, deriv = 0,
type = c("unconditional", "conditional"), famlist = fam.default(),
origin, origin.type = c("model.type", "unconditional", "conditional"))
{
type <- match.arg(type)
origin.type <- match.arg(origin.type)
stopifnot(is.element(fam, seq(along = famlist)))
setfam(famlist)
##### should probably be separate function not snarf-and-barf from aster
nnode <- length(fam)
stopifnot(is.numeric(x))
nind <- floor(length(x) / nnode)
stopifnot(length(x) == nind * nnode)
if (missing(origin)) {
origin <- .C(C_aster_default_origin,
nind = as.integer(nind),
nnode = as.integer(nnode),
fam = as.integer(fam),
theta = matrix(as.double(0), nind, nnode))$theta
if (type == "unconditional")
origin <- .C(C_aster_theta2phi,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
theta = origin,
phi = matrix(as.double(0), nind, nnode))$phi
} else {
stopifnot(is.numeric(origin))
storage.mode(origin) <- "double"
stopifnot(length(origin) == nind * nnode)
stopifnot(all(is.finite(origin)))
if (is.matrix(origin)) {
stopifnot(nrow(origin) == nind)
stopifnot(ncol(origin) == nnode)
} else {
origin <- matrix(origin, nind, nnode)
}
if (origin.type == "model.type")
origin.type <- type
if (type == "unconditional" && origin.type == "conditional")
origin <- .C(C_aster_theta2phi,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
theta = origin,
phi = matrix(as.double(0), nind, nnode))$phi
if (type == "conditional" && origin.type == "unconditional")
origin <- .C(C_aster_phi2theta,
nind = as.integer(nind),
nnode = as.integer(nnode),
pred = as.integer(pred),
fam = as.integer(fam),
phi = origin,
theta = matrix(as.double(0), nind, nnode))$theta
}
result <- mloglhelper(parm, pred, fam, x, root, modmat, origin,
deriv, type)
clearfam()
return(result)
}
mloglhelper <- function(parm, pred, fam, x, root, modmat, origin,
deriv = 0, type = c("unconditional", "conditional"))
{
type <- match.arg(type)
stopifnot(is.numeric(x))
stopifnot(is.numeric(root))
stopifnot(is.numeric(modmat))
stopifnot(is.numeric(parm))
stopifnot(is.numeric(pred))
stopifnot(is.numeric(fam))
stopifnot(is.numeric(deriv))
stopifnot(all(pred == as.integer(pred)))
stopifnot(all(fam == as.integer(fam)))
stopifnot(all(deriv == as.integer(deriv)))
stopifnot(all(pred < seq(along = pred)))
stopifnot(is.element(deriv, 0:2))
stopifnot(length(deriv) == 1)
ncoef <- length(parm)
nnode <- length(pred)
foo <- dim(modmat)
if (ncoef != foo[length(foo)])
stop("last dimension of modmat not length(parm)")
if (length(x) != prod(foo[- length(foo)]))
stop("product of all but last dimension of modmat not length(x)")
stopifnot(length(x) == length(root))
stopifnot(length(pred) == length(fam))
nind <- length(x) / nnode
if(nind != as.integer(nind))
stop("length(x) not multiple of length(pred)")
stopifnot(length(origin) == nind * nnode)
if (type == "unconditional") {
out <- .C(C_aster_mlogl_unco,
nind = as.integer(nind),
nnode = as.integer(nnode),
ncoef = as.integer(ncoef),
pred = as.integer(pred),
fam = as.integer(fam),
deriv = as.integer(deriv),
beta = as.double(parm),
root = as.double(root),
x = as.double(x),
origin = as.double(origin),
modmat = as.double(modmat),
value = double(1),
gradient = double(ncoef),
hessian = matrix(as.double(0), ncoef, ncoef))
} else {
out <- .C(C_aster_mlogl_cond,
nind = as.integer(nind),
nnode = as.integer(nnode),
ncoef = as.integer(ncoef),
pred = as.integer(pred),
fam = as.integer(fam),
deriv = as.integer(deriv),
beta = as.double(parm),
root = as.double(root),
x = as.double(x),
origin = as.double(origin),
modmat = as.double(modmat),
value = double(1),
gradient = double(ncoef),
hessian = matrix(as.double(0), ncoef, ncoef))
}
if (out$value == Inf) {
out$gradient <- NaN * out$gradient
out$hessian <- NaN * out$hessian
}
if (deriv == 0)
return(list(value = out$value))
if (deriv == 1)
return(list(value = out$value, gradient = out$gradient))
if (deriv == 2)
return(list(value = out$value, gradient = out$gradient,
hessian = out$hessian))
}