MLEfast.R
MLE.objective.fn <- function(ltheta, info, value = TRUE) {
# marginal process covariance matrix
y <- as.matrix(info$y)
x <- info$x
smoothness <- info$smoothness
ngrid <- info$ngrid
# number of reps
M <- ncol(y)
Tmatrix <- fields.mkpoly(x, 2)
qr.T <- qr(Tmatrix)
N <- nrow(y)
Q2 <- qr.yq2(qr.T, diag(1, N))
ys <- t(Q2) %*% y
N2 <- length(ys)
theta <- exp(ltheta)
K <- Matern(rdist(x, x)/theta, smoothness = smoothness)
Ke <- eigen(t(Q2) %*% K %*% Q2, symmetric = TRUE)
u2 <- t(Ke$vectors) %*% ys
# mean over replicates -- mean square for coefficients for
# a particular eigenfunction.
u2.MS <- c(rowMeans(u2^2))
D2 <- Ke$values
N2 <- length(D2)
# grid of lambda based on spacing of eigenvalues
ngrid <- min(ngrid, N2)
lambda.grid <- exp(seq(log(D2[1]), log(D2[N2]), , ngrid))
trA <- minus.pflike <- rep(NA, ngrid)
#grid search followed by golden section
# -log likelihood
temp.fn <- function(llam, info) {
lam.temp <- exp(llam)
u2 <- info$u2.MS
D2 <- info$D2
N2 <- length(u2.MS)
# MLE of rho
rho.MLE <- (sum((u2.MS)/(lam.temp + D2)))/N2
# ln determinant
lnDetCov <- sum(log(lam.temp + D2))
-1 * M * (-N2/2 - log(2 * pi) * (N2/2) - (N2/2) * log(rho.MLE) -
(1/2) * lnDetCov)
}
# information list for calling golden section search.
info <- list(D2 = D2, u2 = u2.MS, M = M)
out <- golden.section.search(f = temp.fn, f.extra = info,
gridx = log(lambda.grid), tol = 1e-07)
minus.LogProfileLike <- out$fmin
lambda.MLE <- exp(out$x)
rho.MLE <- (sum((u2.MS)/(lambda.MLE + D2)))/N2
sigma.MLE <- sqrt(lambda.MLE * rho.MLE)
trA <- sum(D2/(lambda.MLE + D2))
pars <- c(rho.MLE, theta, sigma.MLE, trA)
names(pars) <- c("rho", "theta", "sigma", "trA")
if (value) {
return(minus.LogProfileLike)
}
else {
return(list(minus.lPlike = minus.LogProfileLike, lambda.MLE = lambda.MLE,
pars = pars, mle.grid = out$coarse.search))
}
}
MLE.Matern.fast <- function(x, y, smoothness, theta.grid = NULL,
ngrid = 20, verbose = FALSE, m = 2, ...) {
# remove missing values and print out a warning
bad <- is.na(y)
if (sum(bad) > 0) {
cat("removed ", sum(bad), " NAs", fill = TRUE)
x <- x[!bad, ]
y <- y[!bad]
}
# list to pass to the objective function
# NOTE: large ngrid here is very cheap after the eigen decomposition
# has been done.
info <- list(x = x, y = y, smoothness = smoothness, ngrid = 80)
# if grid for ranges is missing use some quantiles of
# pairwise distances among data.
if (is.null(theta.grid)) {
theta.range <- quantile(rdist(x, x), c(0.03, 0.97))
theta.grid <- seq(theta.range[1], theta.range[2], , ngrid)
}
if (length(theta.grid) == 2) {
theta.grid <- seq(theta.grid[1], theta.grid[2], , ngrid)
}
else {
ngrid <- length(theta.grid)
}
# grid search/golden section search
# note that search is in log scale.
out <- golden.section.search(f = MLE.objective.fn, f.extra = info,
gridx = log(theta.grid))
theta.MLE <- exp(out$x)
REML <- -out$fmin
# one final call with the theta.MLE value to recover MLEs for rho and sigma
out2 <- MLE.objective.fn(log(theta.MLE), info, value = FALSE)
return(list(smoothness = smoothness, pars = out2$pars[1:3],
REML = REML, trA = out2$pars[4], REML.grid = cbind(theta.grid,
-1 * out$coarse.search[, 2])))
}