Nothing
#======================================================================================
# Local likelihood estimation for covariance functions with spatially-varying
# parameters: the convoSPAT() package for R
# Mark D. Risser / The Ohio State University / 2014-2015
#======================================================================================
#======================================================================================
# Local/Global Parameter Estimation
#======================================================================================
#======================================================================================
# Global parameter estimation
#======================================================================================
# For a given correlation matrix, the following functions calculate the global
# parameters based on what is NOT specified to be spatially-varying.
#======================================================================================
#===================================
# Calculate spatial covariance
#ROxygen comments ----
#' Calculate spatial covariance.
#'
#' This function replaces the geoR function \code{cov.spatial}, which
#' is now defunct.
#' Options available in this package are: "\code{exponential}",
#' \code{"matern"}, and \code{"gaussian"}.
#'
#' @param Dist.mat A matrix of scaled distances.
#' @param cov.model A string specifying the model for the correlation
#' function; defaults to \code{"exponential"}.
#' Options available in this package are: "\code{exponential}",
#' \code{"matern"}, and \code{"gaussian"}.
#' @param cov.pars Fixed values; not used in the function.
#' @param kappa Scalar; value of the smoothness parameter.
#'
#' @return This function returns a correlation matrix.
#'
#' @examples
#' Distmat <- as.matrix(dist(matrix(runif(20), ncol = 2), diag = TRUE, upper = TRUE))
#' C <- cov_spatial( Dist.mat = Distmat )
#'
#' @export
cov_spatial <- function( Dist.mat, cov.model = "exponential", cov.pars = c(1,1), kappa = 0.5 ){
# Make sure cov.model is one of the permissible options
if( cov.model != "matern" & cov.model != "gaussian" & cov.model != "exponential" ){
stop("Please specify a valid covariance model (matern, gaussian, or exponential).")
}
if(cov.model == "exponential"){
C <- exp(-Dist.mat)
}
if( cov.model == "matern"){
C <- (exp(lgamma(kappa)) * 2^(kappa - 1))^(-1) * Dist.mat^kappa * besselK(Dist.mat, kappa)
C[is.na(C)] <- 1
}
if( cov.model == "gaussian"){
C <- exp(-(Dist.mat^2))
}
return(C)
}
#===================================
# First, models without kappa
# Estimates: nugget, variance
#ROxygen comments ----
#' Constructor functions for global parameter estimation.
#'
#' This function generates another function to be used within \code{optim} to
#' obtain maximum likelihood estimates of
#' global variance parameters tausq, sigmasq with a fixed correlation
#' matrix (smoothness is fixed).
#'
#' @param data A vector or matrix of data to use in the likelihood
#' calculation.
#' @param Xmat The design matrix for the mean model.
#' @param Corr The correlation matrix.
#' @param nugg2.var Fixed values for the covariance of the second nugget term.
#'
#' @return This function returns another function for use in \code{optim}.
#'
#' @examples
#' \dontrun{
#' make_global_loglik1( data, Xmat, Corr, nugg2.var )
#' }
#'
#' @export
make_global_loglik1 <- function( data, Xmat, Corr, nugg2.var ){
fixed <- c(FALSE, FALSE)
params <- fixed
function(p){
params[!fixed] <- p
# Parameters
tausq <- params[1]
sigmasq <- params[2]
N <- dim(data)[1]
p <- dim(data)[2]
Edecomp <- eigen(sigmasq*Corr + diag(rep(tausq,N)) + nugg2.var)
Dmat.temp <- diag(Edecomp$values)
Vmat.temp <- Edecomp$vectors
Cov.inv <- Vmat.temp %*% diag( 1/diag(Dmat.temp) ) %*% t(Vmat.temp)
XCX <- t(Xmat) %*% Cov.inv %*% Xmat
XCX.inv <- chol2inv( chol(XCX) )
Ptemp <- Cov.inv - Cov.inv %*% Xmat %*% XCX.inv %*% t(Xmat) %*% Cov.inv
# Default is to minimize -- want the maximum. Calculate the neg loglikelihood.
loglikelihood <- 0.5*( p*sum(log(diag(Dmat.temp))) + p*log(det(XCX)) + sum(diag(t(data)%*%Ptemp%*%data)) )
if(abs(loglikelihood) == Inf){ loglikelihood <- 100000 }
return(loglikelihood)
}
}
# Estimates: variance
#ROxygen comments ----
#' Constructor functions for global parameter estimation.
#'
#' This function generates another function to be used within \code{optim} to
#' obtain maximum likelihood estimates of
#' global variance parameter sigmasq with a fixed correlation
#' matrix (smoothness is fixed). The nugget variance is taken
#' to be spatially-varing.
#'
#' @param data A vector or matrix of data to use in the likelihood
#' calculation.
#' @param Xmat The design matrix for the mean model.
#' @param Corr The correlation matrix.
#' @param obs.nuggets A vector containing the spatially-varying nuggets
#' corresponding to each data location.
#' @param nugg2.var Fixed values for the covariance of the second nugget term.
#'
#' @return This function returns another function for use in \code{optim}.
#'
#' @examples
#' \dontrun{
#' make_global_loglik2( data, Xmat, Corr, obs.nuggets, nugg2.var )
#' }
#'
#' @export
make_global_loglik2 <- function( data, Xmat, Corr, obs.nuggets, nugg2.var ){
fixed <- c(FALSE)
params <- fixed
function(p){
params[!fixed] <- p
# Parameters
sigmasq <- params[1]
N <- dim(data)[1]
p <- dim(data)[2]
Edecomp <- eigen(sigmasq*Corr + diag(obs.nuggets) + nugg2.var)
Dmat.temp <- diag(Edecomp$values)
Vmat.temp <- Edecomp$vectors
Cov.inv <- Vmat.temp %*% diag( 1/diag(Dmat.temp) ) %*% t(Vmat.temp)
XCX <- t(Xmat) %*% Cov.inv %*% Xmat
XCX.inv <- chol2inv( chol(XCX) )
Ptemp <- Cov.inv - Cov.inv %*% Xmat %*% XCX.inv %*% t(Xmat) %*% Cov.inv
# Default is to minimize -- want the maximum. Calculate the neg loglikelihood.
loglikelihood <- 0.5*( p*sum(log(diag(Dmat.temp))) + p*log(det(XCX)) + sum(diag(t(data)%*%Ptemp%*%data)) )
if(abs(loglikelihood) == Inf){ loglikelihood <- 100000 }
return(loglikelihood)
}
}
# Estimates: nugget
#ROxygen comments ----
#' Constructor functions for global parameter estimation.
#'
#' This function generates another function to be used within \code{optim} to
#' obtain maximum likelihood estimates of
#' global variance parameter tausq with a fixed correlation
#' matrix (smoothness is fixed). The process variance is taken
#' to be spatially-varing.
#'
#' @param data A vector or matrix of data to use in the likelihood
#' calculation.
#' @param Xmat The design matrix for the mean model.
#' @param Corr The correlation matrix matrix.
#' @param obs.variance A vector containing the spatially-varying variance
#' corresponding to each data location.
#' @param nugg2.var Fixed values for the covariance of the second nugget term.
#'
#' @return This function returns another function for use in \code{optim}.
#'
#' @examples
#' \dontrun{
#' make_global_loglik3( data, Xmat, Corr, obs.variance, nugg2.var )
#' }
#'
#' @export
make_global_loglik3 <- function( data, Xmat, Corr, obs.variance, nugg2.var ){
fixed <- c(FALSE)
params <- fixed
function(p){
params[!fixed] <- p
tausq <- params[1]
N <- dim(data)[1]
p <- dim(data)[2]
Edecomp <- eigen(diag(sqrt(obs.variance)) %*% Corr %*% diag(sqrt(obs.variance))
+ diag(rep(tausq,N)) + nugg2.var )
Dmat.temp <- diag(Edecomp$values)
Vmat.temp <- Edecomp$vectors
Cov.inv <- Vmat.temp %*% diag( 1/diag(Dmat.temp) ) %*% t(Vmat.temp)
XCX <- t(Xmat) %*% Cov.inv %*% Xmat
XCX.inv <- chol2inv( chol(XCX) )
Ptemp <- Cov.inv - Cov.inv %*% Xmat %*% XCX.inv %*% t(Xmat) %*% Cov.inv
# Default is to minimize -- want the maximum. Calculate the neg loglikelihood.
loglikelihood <- 0.5*( p*sum(log(diag(Dmat.temp))) + p*log(det(XCX)) + sum(diag(t(data)%*%Ptemp%*%data)) )
if(abs(loglikelihood) == Inf){ loglikelihood <- 100000 }
return(loglikelihood)
}
}
#===================================
#===================================
# Next, models with kappa
# Estimates: nugget, variance, kappa
#ROxygen comments ----
#' Constructor functions for global parameter estimation.
#'
#' This function generates another function to be used within \code{optim} to
#' obtain maximum likelihood estimates of
#' global variance parameters tausq, sigmasq, and nu.
#'
#' @param data A vector or matrix of data to use in the likelihood
#' calculation.
#' @param Xmat The design matrix for the mean model.
#' @param cov.model String; the covariance model.
#' @param Scalemat Matrix; contains the scaling quantities from the
#' covariance function.
#' @param Distmat Matrix; contains the scaled distances.
#' @param nugg2.var Fixed values for the covariance of the second nugget term.
#'
#' @return This function returns another function for use in \code{optim}.
#'
#' @examples
#' \dontrun{
#' make_global_loglik1_kappa( data, Xmat, cov.model, Scalemat, Distmat, nugg2.var )
#' }
#'
#' @export
make_global_loglik1_kappa <- function( data, Xmat, cov.model, Scalemat, Distmat, nugg2.var ){
fixed <- c(FALSE, FALSE, FALSE)
params <- fixed
function(p){
params[!fixed] <- p
# Parameters
tausq <- params[1]
sigmasq <- params[2]
kapp <- params[3]
Unscl.corr <- cov_spatial( Distmat, cov.model = cov.model,
cov.pars = c(1,1), kappa = kapp )
Corr <- Scalemat*Unscl.corr
N <- dim(data)[1]
p <- dim(data)[2]
Edecomp <- eigen(sigmasq*Corr + diag(rep(tausq,N)) + nugg2.var)
Dmat.temp <- diag(Edecomp$values)
Vmat.temp <- Edecomp$vectors
Cov.inv <- Vmat.temp %*% diag( 1/diag(Dmat.temp) ) %*% t(Vmat.temp)
XCX <- t(Xmat) %*% Cov.inv %*% Xmat
XCX.inv <- chol2inv( chol(XCX) )
Ptemp <- Cov.inv - Cov.inv %*% Xmat %*% XCX.inv %*% t(Xmat) %*% Cov.inv
# Default is to minimize -- want the maximum. Calculate the neg loglikelihood.
loglikelihood <- 0.5*( p*sum(log(diag(Dmat.temp))) + p*log(det(XCX)) + sum(diag(t(data)%*%Ptemp%*%data)) )
if(abs(loglikelihood) == Inf){ loglikelihood <- 100000 }
return(loglikelihood)
}
}
# Estimates: variance, kappa
#ROxygen comments ----
#' Constructor functions for global parameter estimation.
#'
#' This function generates another function to be used within \code{optim} to
#' obtain maximum likelihood estimates of
#' global variance parameters sigmasq and nu. The nugget variance is
#' taken to be spatially-varying.
#'
#' @param data A vector or matrix of data to use in the likelihood
#' calculation.
#' @param Xmat The design matrix for the mean model.
#' @param cov.model String; the covariance model.
#' @param Scalemat Matrix; contains the scaling quantities from the
#' covariance function.
#' @param Distmat Matrix; contains the scaled distances.
#' @param obs.nuggets A vector containing the spatially-varying nuggets
#' corresponding to each data location.
#' @param nugg2.var Fixed values for the covariance of the second nugget term.
#'
#' @return This function returns another function for use in \code{optim}.
#'
#' @examples
#' \dontrun{
#' make_global_loglik2_kappa( data, Xmat, cov.model, Scalemat, Distmat, obs.nuggets, nugg2.var )
#' }
#'
#' @export
make_global_loglik2_kappa <- function( data, Xmat, cov.model, Scalemat, Distmat, obs.nuggets, nugg2.var ){
fixed <- c(FALSE, FALSE)
params <- fixed
function(p){
params[!fixed] <- p
# Parameters
sigmasq <- params[1]
kapp <- params[2]
Unscl.corr <- cov_spatial( Distmat, cov.model = cov.model,
cov.pars = c(1,1), kappa = kapp )
Corr <- Scalemat*Unscl.corr
N <- dim(data)[1]
p <- dim(data)[2]
Edecomp <- eigen(sigmasq * Corr + diag(obs.nuggets) + nugg2.var)
Dmat.temp <- diag(Edecomp$values)
Vmat.temp <- Edecomp$vectors
Cov.inv <- Vmat.temp %*% diag( 1/diag(Dmat.temp) ) %*% t(Vmat.temp)
XCX <- t(Xmat) %*% Cov.inv %*% Xmat
XCX.inv <- chol2inv( chol(XCX) )
Ptemp <- Cov.inv - Cov.inv %*% Xmat %*% XCX.inv %*% t(Xmat) %*% Cov.inv
# Default is to minimize -- want the maximum. Calculate the neg loglikelihood.
loglikelihood <- 0.5*( p*sum(log(diag(Dmat.temp))) + p*log(det(XCX)) + sum(diag(t(data)%*%Ptemp%*%data)) )
if(abs(loglikelihood) == Inf){ loglikelihood <- 100000 }
return(loglikelihood)
}
}
# Estimates: nugget, kappa
#ROxygen comments ----
#' Constructor functions for global parameter estimation.
#'
#' This function generates another function to be used within \code{optim} to
#' obtain maximum likelihood estimates of
#' global variance parameters tausq and nu. The process variance is
#' taken to be spatially-varying.
#'
#' @param data A vector or matrix of data to use in the likelihood
#' calculation.
#' @param Xmat The design matrix for the mean model.
#' @param cov.model String; the covariance model.
#' @param Scalemat Matrix; contains the scaling quantities from the
#' covariance function.
#' @param Distmat Matrix; contains the scaled distances.
#' @param obs.variance A vector containing the spatially-varying variance
#' corresponding to each data location.
#' @param nugg2.var Fixed values for the covariance of the second nugget term.
#'
#' @return This function returns another function for use in \code{optim}.
#'
#' @examples
#' \dontrun{
#' make_global_loglik3_kappa( data, Xmat, cov.model, Scalemat, Distmat, obs.variance, nugg2.var )
#' }
#'
#' @export
make_global_loglik3_kappa <- function( data, Xmat, cov.model, Scalemat, Distmat, obs.variance, nugg2.var ){
fixed <- c(FALSE, FALSE)
params <- fixed
function(p){
params[!fixed] <- p
# Parameters
tausq <- params[1]
kapp <- params[2]
Unscl.corr <- cov_spatial( Distmat, cov.model = cov.model,
cov.pars = c(1,1), kappa = kapp )
Corr <- Scalemat*Unscl.corr
N <- dim(data)[1]
p <- dim(data)[2]
Edecomp <- eigen(diag(sqrt(obs.variance)) %*% Corr %*% diag(sqrt(obs.variance))
+ diag(rep(tausq,N)) + nugg2.var )
Dmat.temp <- diag(Edecomp$values)
Vmat.temp <- Edecomp$vectors
Cov.inv <- Vmat.temp %*% diag( 1/diag(Dmat.temp) ) %*% t(Vmat.temp)
XCX <- t(Xmat) %*% Cov.inv %*% Xmat
XCX.inv <- chol2inv( chol(XCX) )
Ptemp <- Cov.inv - Cov.inv %*% Xmat %*% XCX.inv %*% t(Xmat) %*% Cov.inv
# Default is to minimize -- want the maximum. Calculate the neg loglikelihood.
loglikelihood <- 0.5*( p*sum(log(diag(Dmat.temp))) + p*log(det(XCX)) + sum(diag(t(data)%*%Ptemp%*%data)) )
if(abs(loglikelihood) == Inf){ loglikelihood <- 100000 }
return(loglikelihood)
}
}
# Estimates: kappa
#ROxygen comments ----
#' Constructor functions for global parameter estimation.
#'
#' This function generates another function to be used within \code{optim} to
#' obtain maximum likelihood estimates of
#' global variance parameters nu. The process variance
#' and nugget variance are taken to be spatially-varying.
#'
#' @param data A vector or matrix of data to use in the likelihood
#' calculation.
#' @param Xmat The design matrix for the mean model.
#' @param cov.model String; the covariance model.
#' @param Scalemat Matrix; contains the scaling quantities from the
#' covariance function.
#' @param Distmat Matrix; contains the scaled distances.
#' @param obs.variance A vector containing the spatially-varying variance
#' corresponding to each data location.
#' @param obs.nuggets A vector containing the spatially-varying nuggets
#' corresponding to each data location.
#' @param nugg2.var Fixed values for the covariance of the second nugget term.
#'
#' @return This function returns another function for use in \code{optim}.
#'
#' @examples
#' \dontrun{
#' make_global_loglik4_kappa( data, Xmat, cov.model, Scalemat, Distmat,
#' obs.variance, obs.nuggets, nugg2.var )
#' }
#'
#' @export
make_global_loglik4_kappa <- function( data, Xmat, cov.model, Scalemat, Distmat, obs.variance,
obs.nuggets, nugg2.var ){
fixed <- c(FALSE)
params <- fixed
function(p){
params[!fixed] <- p
# Parameters
kapp <- params[1]
Unscl.corr <- cov_spatial( Distmat, cov.model = cov.model,
cov.pars = c(1,1), kappa = kapp )
Corr <- Scalemat*Unscl.corr
N <- dim(data)[1]
p <- dim(data)[2]
Edecomp <- eigen(diag( sqrt(obs.variance) ) %*% Corr %*% diag( sqrt(obs.variance) ) + diag(obs.nuggets) + nugg2.var)
Dmat.temp <- diag(Edecomp$values)
Vmat.temp <- Edecomp$vectors
Cov.inv <- Vmat.temp %*% diag( 1/diag(Dmat.temp) ) %*% t(Vmat.temp)
XCX <- t(Xmat) %*% Cov.inv %*% Xmat
XCX.inv <- chol2inv( chol(XCX) )
Ptemp <- Cov.inv - Cov.inv %*% Xmat %*% XCX.inv %*% t(Xmat) %*% Cov.inv
# Default is to minimize -- want the maximum. Calculate the neg loglikelihood.
loglikelihood <- 0.5*( p*sum(log(diag(Dmat.temp))) + p*log(det(XCX)) + sum(diag(t(data)%*%Ptemp%*%data)) )
if(abs(loglikelihood) == Inf){ loglikelihood <- 100000 }
return(loglikelihood)
}
}
#======================================================================================
#======================================================================================
# Function to estimate the local models
#======================================================================================
# Using a subset of the data, calculate MLEs of covariance and mean
# parameters. Options for
# (1) method: ml vs. reml
# (2) smoothness parameter: no kappa vs. estimate kappa vs. fixed kappa
# (3) locally isotropic vs. locally anisotropic
# (4) fixed tausq vs. estimate tausq
#======================================================================================
#ROxygen comments ----
#' Constructor functions for local parameter estimation.
#'
#' This function generates another function to be used within \code{optim} to
#' obtain maximum likelihood estimates of covariance (and possibly mean) parameters.
#' The function includes options for
#' (1) maximum likelihood (\code{"ml"}) vs. restricted maximum likelihood
#' (\code{"reml"}),
#' (2) smoothness (\code{kappa}): models without smoothness vs. estimating the
#' smoothness vs. using fixed smoothness,
#' (3) locally isotropic vs. locally anisotropic, and
#' (4) fixed nugget variance (\code{tausq}): fixed vs. estimated.
#'
#' @param locations A matrix of locations.
#' @param cov.model String; the covariance model.
#' @param data A vector or matrix of data to use in the likelihood
#' calculation.
#' @param Xmat The design matrix for the mean model.
#' @param nugg2.var Fixed values for the variance/covariance of the second nugget term; defaults
#' to a matrix of zeros.
#' @param tausq Scalar; fixed value for the nugget variance (when
#' \code{fix.tausq = TRUE}).
#' @param kappa Scalar; fixed value for the smoothness (when \code{fix.kappa = TRUE}).
#' @param fixed Logical vector of \code{FALSE} values; length corresponds to the number
#' of parameters to be estimated.
#' @param method Indicates the estimation method, either maximum likelihood (\code{"ml"})
#' or restricted maximum likelihood (\code{"reml"}).
#' @param local.aniso Logical; indicates if the local covariance should be
#' anisotropic (\code{TRUE}) or isotropic (\code{FALSE}). Defaults to \code{TRUE}.
#' @param fix.kappa Logical; indicates if the kappa parameter should be
#' fixed (\code{TRUE}) or estimated (\code{FALSE}). Defaults to \code{FALSE}
#' (only valid for \code{cov.model = "matern"} and \code{cov.model = "cauchy"}).
#' @param fix.tausq Logical; indicates whether the default nugget term
#' (tau^2) should be fixed (\code{TRUE}) or estimated (\code{FALSE}). Defaults to
#' \code{FALSE}.
#'
#' @return This function returns another function for use in \code{optim}.
#'
#' @examples
#' \dontrun{
#' make_local_lik( locations, cov.model, data, Xmat )
#' }
#'
#' @export
#' @importFrom StatMatch mahalanobis.dist
make_local_lik <- function( locations, cov.model, data, Xmat,
nugg2.var = matrix(0, nrow(locations), nrow(locations)),
tausq = 0,
kappa = 0.5, fixed = rep(FALSE, 6),
method = "reml", local.aniso = TRUE,
fix.tausq = FALSE, fix.kappa = FALSE ){
params <- fixed
function(p) {
params[!fixed] <- p
if( local.aniso ){ # Locally anisotropic
lam1 <- params[1]
lam2 <- params[2]
eta <- params[3]
if( cov.model == "matern" || cov.model == "cauchy" ){ # For covariance models with kappa
if( !fix.kappa ){ # Estimate kappa
if( !fix.tausq ){ # Estimate tausq
tausq <- params[4]
sigmasq <- params[5]
kappa <- params[6]
if( method == "ml" ) beta <- params[-c(1:6)]
}
if( fix.tausq ){ # Fix tausq
tausq <- tausq
sigmasq <- params[4]
kappa <- params[5]
if( method == "ml" ) beta <- params[-c(1:5)]
}
}
if( fix.kappa ){ # Fix kappa
if( !fix.tausq ){ # Estimate tausq
tausq <- params[4]
sigmasq <- params[5]
kappa <- kappa
if( method == "ml" ) beta <- params[-c(1:5)]
}
if( fix.tausq ){ # Fix tausq
tausq <- tausq
sigmasq <- params[4]
kappa <- kappa
if( method == "ml" ) beta <- params[-c(1:4)]
}
}
}
if( cov.model != "matern" & cov.model != "cauchy" ){ # For covariance models without kappa
if( !fix.tausq ){ # Estimate tausq
tausq <- params[4]
sigmasq <- params[5]
if( method == "ml" ) beta <- params[-c(1:5)]
}
if( fix.tausq ){ # Fix tausq
tausq <- tausq
sigmasq <- params[4]
if( method == "ml" ) beta <- params[-c(1:4)]
}
}
}
if( !local.aniso ){ # Locally isotropic
lam1 <- params[1]
lam2 <- params[1]
eta <- 0
if( cov.model == "matern" || cov.model == "cauchy" ){ # For covariance models with kappa
if( !fix.kappa ){ # Estimate kappa
if( !fix.tausq ){ # Estimate tausq
tausq <- params[2]
sigmasq <- params[3]
kappa <- params[4]
if( method == "ml" ) beta <- params[-c(1:4)]
}
if( fix.tausq ){ # Fix tausq
tausq <- tausq
sigmasq <- params[2]
kappa <- params[3]
if( method == "ml" ) beta <- params[-c(1:3)]
}
}
if( fix.kappa ){ # Fix kappa
if( !fix.tausq ){ # Estimate tausq
tausq <- params[2]
sigmasq <- params[3]
kappa <- kappa
if( method == "ml" ) beta <- params[-c(1:3)]
}
if( fix.tausq ){ # Fix tausq
tausq <- tausq
sigmasq <- params[2]
kappa <- kappa
if( method == "ml" ) beta <- params[-c(1:2)]
}
}
}
if( cov.model != "matern" & cov.model != "cauchy" ){ # For covariance models without kappa
if( !fix.tausq ){ # Estimate tausq
tausq <- params[2]
sigmasq <- params[3]
if( method == "ml" ) beta <- params[-c(1:3)]
}
if( fix.tausq ){ # Fix tausq
tausq <- tausq
sigmasq <- params[2]
if( method == "ml" ) beta <- params[-c(1:2)]
}
}
}
# Setup and covariance calculation
N <- dim(locations)[1]
m <- dim(data)[2]
Pmat <- matrix(c(cos(eta), -sin(eta), sin(eta), cos(eta)), nrow = 2, byrow = T)
Dmat <- diag(c(lam1, lam2))
Sigma <- Pmat %*% Dmat %*% t(Pmat)
distances <- mahalanobis.dist(data.x = locations, vc = Sigma)
NS.cov <- nugg2.var + tausq * diag(N) + sigmasq * cov_spatial(distances, cov.model = cov.model,
cov.pars = c(1, 1), kappa = kappa)
possibleError <- tryCatch(
cov.chol <- chol(NS.cov),
error=function(e) e
)
# Check for error before calculating the likelihood
if(inherits(possibleError, "error")){
loglikelihood <- 1e+06
} else{
# Likelihood calculation
if( method == "ml" ){
tmp1 <- backsolve(cov.chol, data - Xmat %*% beta, transpose = TRUE)
ResCinvRes <- t(tmp1) %*% tmp1
loglikelihood <- m * sum(log(diag(cov.chol))) + 0.5 * sum(diag(ResCinvRes))
}
if( method == "reml" ){
tmp1 <- backsolve(cov.chol, Xmat, transpose = TRUE)
XCinvX <- t(tmp1) %*% tmp1
xcx.chol <- chol(XCinvX)
tmp2 <- backsolve(cov.chol, data, transpose = TRUE)
ZCinvZ <- t(tmp2) %*% tmp2
XCinvZ <- t(Xmat) %*% backsolve(cov.chol, tmp2)
tmp3 <- backsolve(xcx.chol, XCinvZ, transpose = TRUE)
qf <- t(tmp3) %*% tmp3
loglikelihood <- m * sum(log(diag(cov.chol))) + m * sum(log(diag(xcx.chol))) + 0.5 * sum(diag(ZCinvZ)) - 0.5 * sum(diag(qf))
}
}
if (abs(loglikelihood) == Inf) { loglikelihood <- 1e+06 } # Make sure not Inf
return(loglikelihood)
}
}
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.