tensor_predictors/tensorPredictors/R/GMLM.R

65 lines
2.0 KiB
R

#' Fitting Generalized Multi-Linear Models
#'
#' @export
GMLM <- function(...) {
stop("Not Implemented")
}
#' @export
GMLM.default <- function(X, Fy, sample.axis = 1L,
family = "normal",
...,
eps = sqrt(.Machine$double.eps),
logger = NULL
) {
stopifnot(exprs = {
length(sample.axis) == 1L
(1L <= sample.axis) && (sample.axis <= length(dim(X)))
(dim(X) == dim(Fy))[sample.axis]
})
# rearrange `X`, `Fy` such that the last axis enumerates observations
axis.perm <- c(seq_along(dim(X))[-sample.axis], sample.axis)
X <- aperm(X, axis.perm)
Fy <- aperm(Fy, axis.perm)
# setup family specific GLM (pseudo) "inverse" link
family <- make.gmlm.family(family)
# wrap logger in callback for NAGD
callback <- if (is.function(logger)) {
function(iter, params) {
do.call(logger, c(list(iter), params))
}
}
# optimize likelihood using Nesterov Excelerated Gradient Descent
params.fit <- NAGD(
fun.loss = function(params) {
# scaled negative log-likelihood
# eta1 alphas Omegas
family$log.likelihood(X, Fy, params[[1]], params[[2]], params[[3]])
},
fun.grad = function(params) {
# gradient of the scaled negative log-likelihood
# eta1 alphas Omegas
family$grad(X, Fy, params[[1]], params[[2]], params[[3]])
},
params = family$initialize(X, Fy), # initial parameter estimates
fun.lincomb = function(a, lhs, b, rhs) {
list(
a * lhs[[1]] + b * rhs[[1]],
Map(function(l, r) a * l + b * r, lhs[[2]], rhs[[2]]),
Map(function(l, r) a * l + b * r, lhs[[3]], rhs[[3]])
)
},
fun.norm2 = function(params) {
sum(unlist(params)^2)
},
callback = callback,
...
)
structure(params.fit, names = c("eta1", "alphas", "Omegas"))
}