R/kriging.R In geoFKF: Kriging Method for Spatial Functional Data

Documented in geo_fkflog_lik_rf

#' Maximum likelihood estimate for \eqn{\sigma^2}, \eqn{\phi} and \eqn{\rho}.
#'
#' This function maximum likelihood estimate for \eqn{\sigma^2}, \eqn{\phi}
#'  and \eqn{\rho} in the random field model for the covariance.
#'
#' @param m_coef Matrix where each column is an observed vector
#' @param m_coord Matrix where each observation contains the latitude and
#' longitude
#' @export
#' @return Return a list with
#' \describe{
#'    \item{par}{A vector with the estimates of \eqn{\sigma^2}, \eqn{\phi} and
#'    \eqn{\rho}.}
#'    \item{m_cov}{A matrix of covariances of the estimates.}
#' }
#' @examples
#'
#' m_data <- as.matrix(datasetCanada$m_data) #' m_coord <- as.matrix(datasetCanada$m_coord[, 1:2])
#'
#' p <- ceiling(1 + log2(nrow(m_data)))
#' m_coef <- sapply(seq_len(nrow(m_coord)), function(i) {
#'     coef_fourier(m_data[, i], p)
#' })
#' log_lik_rf(m_coef, m_coord)
log_lik_rf <- function(m_coef, m_coord) {

m_dist <- as.matrix(dist(m_coord))

m <- (nrow(m_coef) - 1) / 2

m_dif <- sapply(seq_len(2 * m + 1), function(j) abs(seq_len(2 * m + 1) - j))

fn <- function(arg) {
as.numeric(logLikMultiNorm(m_coef, m_dist, arg[1], arg[2], arg[3]))
}

inicial <- c(2, 0.5, 0.25) + rnorm(3, sd = 0.05)

fit <- nlminb(inicial, fn,
lower = c(0.001, 0.001, -0.999), upper = c(Inf, Inf, 0.999))

m_cov <- with(fit, solve(numDeriv::hessian(fn, par)))

if (min(diag(m_cov)) < 0) stop("Error in computing the hessian matrix.")

with(fit, list(par = par, m_cov = m_cov))
}

#' Kriging method for Spatial Functional Data.
#'
#' \code{geo_fkf} implements the kriging method for spatial functional datasets.
#'
#' @param m_data a tibble where each column or variable is data from a station
#' @param m_coord a tibble with two columns: latitude and longitude
#' @param new_loc a tible with one observation, where the columns or variables
#' are latitude and longitude
#' @param p order in the Fourier Polynomial
#' @param t a time series with values belonging to \eqn{[-\pi, \pi]} to
#' evaluate the estimate curve
#' @export
#'
#' @return a list with three entries: \code{estimates}, \code{Theta} and
#' \code{cov_params}
#' \describe{
#'   \item{estimates}{the estimate curve}
#'   \item{Theta}{weights (matrices) of the linear combination}
#'   \item{cov_params}{estimate \eqn{\sigma^2}, \eqn{\phi} and \eqn{\rho}}
#' }
#' @examples
#'
#'m_data <- as.matrix(datasetCanada$m_data) #'m_coord <- as.matrix(datasetCanada$m_coord[, 1:2])
#'pos <- sample.int(nrow(m_coord), 1)
#'log_pos <- !(seq_len(nrow(m_coord)) %in% pos)
#'new_loc <- m_coord[pos, ]
#'m_coord <- m_coord[log_pos, ]
#'m_data <- m_data[, log_pos]
#'
#'geo_fkf(m_data, m_coord, new_loc)
geo_fkf <- function(m_data, m_coord, new_loc, p,
t = seq(from = -pi, to = pi, by = 0.01)) {

n_loc <- nrow(m_coord)

if (missing(p)) p <- ceiling(1 + log2(nrow(m_data)))

m_coef <- sapply(seq_len(n_loc), function(i) {
x <- (as.matrix(m_data))[, i]
coef_fourier(x, p)
})

fit_cov <- log_lik_rf(as.matrix(m_coef), as.matrix(m_coord))
s2 <- fit_cov$par[1] phi <- fit_cov$par[2]
rho <- fit_cov\$par[3]

m_dist <- as.matrix(dist(rbind(as.matrix(m_coord),
as.vector(new_loc))))
m_dif <- sapply(seq_len(2 * p + 1),
function(j) abs(seq_len(2 * p + 1) - j))

m_a <- matrix(0, nrow = (n_loc + 1) * (2 * p + 1),
ncol = (n_loc + 1) * (2 * p + 1))

for (i in seq_len(n_loc + 1)) {
v_rows <- (n_loc * (2 * p + 1) + 1):((n_loc + 1) * (2 * p + 1))
v_cols <- ((i - 1) * (2 * p + 1) + 1):(i * (2 * p + 1))
m_a[v_rows, v_cols] <- m_a[v_cols, v_rows] <- diag(2 * p + 1)
}

for (i in 1:n_loc) {
for (j in i:n_loc) {
v_rows <- ((i - 1) * (2 * p + 1) + 1):(i * (2 * p + 1))
v_cols <- ((j - 1) * (2 * p + 1) + 1):(j * (2 * p + 1))
m_a[v_cols, v_rows] <-
m_a[v_rows, v_cols] <-
s2 * exp(-m_dist[i, j] / phi) * rho^m_dif * diag(2 * p + 1)
}
}

m_b <- matrix(1, nrow = (n_loc + 1) * (2 * p + 1), ncol = 1)
for (i in seq_len(n_loc)) {
v_rows <- ((i - 1) * (2 * p + 1) + 1):(i * (2 * p + 1))
m_b[v_rows] <- diag(s2 * exp(-m_dist[i, n_loc + 1] / phi) * rho^m_dif)
}

solution <- matrix(solve(m_a, m_b), nrow = 2 * p + 1)

c0 <- matrix(0, nrow = nrow(m_coef), ncol = 1)
for (i in seq_len(n_loc)) {
c0 <- c0 + diag(solution[, i]) %*% m_coef[, i]
}

list(estimates = fourier_b(c0, x = t),
Theta = solution[, seq_len(n_loc)],
cov_params = list(s2 = s2, phi = phi, rho = rho))
}


Try the geoFKF package in your browser

Any scripts or data that you put into this service are public.

geoFKF documentation built on Aug. 13, 2022, 1:05 a.m.