Nothing
cov.sp <- function(coords, sp.type = "exponential",
sp.par = stop("specify sp.par argument"),
error.var = 0, smoothness = 0.5, finescale.var = 0, pcoords = NULL,
D = NULL, Dp = NULL, Dop = NULL)
{
#check validity of arguments
cov_sp_arg_check(coords, sp.type, sp.par, error.var, smoothness,
finescale.var, pcoords, D, Dp, Dop)
# calculate covariance matrix for observed responses
# calculate distance matrix if not supplied
if(is.null(D))
{
V <- simple.cov.sp(dist1(coords), sp.type, sp.par, error.var, smoothness, finescale.var)
}
else
{
V <- simple.cov.sp(D, sp.type, sp.par, error.var, smoothness, finescale.var)
}
#calculate covariance matrices for predicted responses and
#between observed responses and predicted responses
#calculate Dp and Dop if necessary
if(!is.null(pcoords))
{
if(is.null(Dp))
{
Vp <- simple.cov.sp(dist1(pcoords), sp.type, sp.par, error.var = 0, smoothness, finescale.var)
}
else
{
Vp <- simple.cov.sp(Dp, sp.type, sp.par, error.var = 0, smoothness, finescale.var)
}
if(is.null(Dop))
{
Vop <- simple.cov.sp(dist2(coords, pcoords), sp.type, sp.par, error.var = 0, smoothness, finescale.var)
}
else
{
Vop <- simple.cov.sp(Dop, sp.type, sp.par, error.var = 0, smoothness, finescale.var)
}
}
if(is.null(pcoords))
{
return(list(V = V))
}else
{
return(list(V = V, Vp = Vp, Vop = Vop))
}
}
simple.cov.sp <- function(D, sp.type, sp.par, error.var, smoothness, finescale.var)
{
if(sp.type == "exponential")
{
V <- sp.par[1]*exp(-D/sp.par[2])
}else if(sp.type == "gaussian")
{
V <- sp.par[1]*exp(-(D/sp.par[2])^2)
}else if(sp.type == "matern")
{
sD <- D/sp.par[2]
V <- (D > 0) * sp.par[1]*(2^(1-smoothness)/gamma(smoothness)*sD^smoothness*besselK(sD, nu = smoothness))
V[is.nan(V)] <- sp.par[1]
}else if(sp.type == "matern2")
{
sD <- 2 * sqrt(smoothness) * D/sp.par[2]
V <- (D > 0) * sp.par[1]*(2^(1-smoothness)/gamma(smoothness)*sD^smoothness*besselK(sD, nu = smoothness))
V[is.nan(V)] <- sp.par[1]
}else if(sp.type == "spherical")
{
sD <- D/sp.par[2]
V <- sp.par[1]*(1 - 1.5*sD + 0.5*(sD)^3)*(D < sp.par[2])
}
return(V + (D == 0)*(finescale.var + error.var))
}
cov.st <- function(coords, time, sp.type = "exponential", sp.par = stop("specify sp.par argument"),
error.var = 0, smoothness = 0.5, finescale.var = 0, t.type = "ar1", t.par = .5, pcoords = NULL,
ptime = NULL, D = NULL, Dp = NULL, Dop = NULL, T = NULL, Tp = NULL, Top = NULL)
{
if(!is.matrix(time)){ time <- matrix(time) }
if(!is.null(ptime))
{
if(!is.matrix(ptime)){ ptime <- matrix(ptime) }
}
cov_st_arg_check(coords, time, sp.type, sp.par,
error.var, smoothness, finescale.var, t.type, t.par, pcoords,
ptime, D, Dp, Dop, T, Tp, Top)
if(is.null(D)){ D <- dist1(coords) }
V.sp <- simple.cov.sp(D, sp.type, sp.par, error.var = 0, smoothness, finescale.var = 0)
if(is.null(T)){ T <- dist1(matrix(time)) }
V.time <- simple.cov.time(T, t.type, t.par)
#Only add finescale or measurement error when at the same location and time
V <- V.sp * V.time + (D==0) * (T==0) * (finescale.var + error.var)
if(!is.null(pcoords))
{
if(is.null(ptime)){ stop("ptime must be specified since pcoords specified") }
if(is.null(Dp)){ Dp <- dist1(pcoords) }
if(is.null(Dop)){ Dop <- dist2(coords, pcoords) }
if(is.null(Tp)){ Tp <- dist1(matrix(ptime)) }
if(is.null(Top)){ Top <- dist2(matrix(time), matrix(ptime)) }
Vp.sp <- simple.cov.sp(Dp, sp.type, sp.par, error.var = 0, smoothness, finescale.var = 0)
Vp.time <- simple.cov.time(Tp, t.type, t.par)
Vp <- Vp.sp*Vp.time + (Dp == 0)*finescale.var
Vop.sp <- simple.cov.sp(Dop, sp.type, sp.par, error.var = 0, smoothness, finescale.var = 0)
Vop.time <- simple.cov.time(Top, t.type, t.par)
Vop <- Vop.sp*Vop.time + (Dop == 0)*finescale.var
}
if(is.null(pcoords))
{
return(list(V = V))
}else
{
return(list(V = V, Vp = Vp, Vop = Vop))
}
}
simple.cov.time <- function(T, t.type, t.par)
{
if(t.type == "ar1")
{
V.time <- t.par^T
}
V.time
}
decomp.cov <- function(V, method = "eigen")
{
decomp_cov_check_arg(V, method)
if(method == "eigen")
{
eigenV <- eigen(V)
return(eigenV$vectors %*% diag(sqrt(pmax(eigenV$values,0))))
}else if(method == "chol")
{
return(t(chol(V)))
}else if(method == "svd")
{
svdV <- svd(V)
return(tcrossprod(svdV$u %*% diag(sqrt(svdV$d)), svdV$v))
}
}
#marginally faster than decomp.V for small matrices, but much slower for large
#decomp.cov2 <- function(V, method = "eigen")
#{
# decomp_cov_check_arg(V, method)
#
# if(method == "eigen")
# {
# .Call( "decomp_cov", Vs = V,
# method = 1, PACKAGE = "SpatialTools")
# }else if(method == "chol")
# {
# .Call( "decomp_cov", Vs = V,
# method = 2, PACKAGE = "SpatialTools")
# }else if(method == "svd")
# {
# .Call( "decomp_cov", Vs = V,
# method = 3, PACKAGE = "SpatialTools")
# }
#}
maxlik.cov.sp <- function(X, y, coords, sp.type = "exponential",
range.par = stop("specify range.par argument"),
error.ratio = stop("specify error.ratio argument"), smoothness = 0.5,
D = NULL, reml = TRUE, lower = NULL, upper = NULL,
control = list(trace = TRUE), optimizer="nlminb")
{
y <- as.vector(y)
maxlik_cov_sp_check_arg(X, y, coords, sp.type,
range.par, error.ratio, smoothness, D, reml, lower, upper)
if(is.null(D)){ D <- dist1(coords) }
#construct vector of initial parameters to pass to objective function
parms <- c(range.par, error.ratio)
if(sp.type == "matern"){ parms <- c(parms, smoothness) }
#define lower and upper bounds if not supplied by user
bounds <- bounds.cov.sp(sp.type)
if(is.null(lower)){ lower <- bounds$lower.bound }
if(is.null(upper)){ upper <- bounds$upper.bound }
if(optimizer=="nlminb")
{
min.out <- stats::nlminb(start=parms,
objective = logLik.cov.sp, X=X, y=y, D=D,
sp.type = sp.type, reml=reml, minus2 = TRUE, lower=lower, upper=upper,
control = control)
range.par <- min.out$par[1]
error.ratio <- min.out$par[2]
if(sp.type == "matern"){ smoothness <- min.out$par[3] }
}
V <- simple.cov.sp(D = D, sp.type = sp.type, sp.par = c(1, range.par),
error.var = error.ratio, smoothness = smoothness, finescale.var = 0)
ViX <- solve(V, X)
XtViX <- crossprod(ViX, X)
fitted <- X %*% solve(XtViX, crossprod(ViX, y))
r <- y - fitted
n <- length(y)
df.resid <- n - ncol(X)
rtVir <- crossprod(r, solve(V,r))
if(!reml)
{
sigmasq <- c(rtVir/n)
}else
{
sigmasq <- c(rtVir/df.resid)
}
error.var <- sigmasq * error.ratio
sp.par <- c(sigmasq, range.par)
return(list(sp.type = sp.type, sp.par = sp.par,
error.var = error.var, smoothness = smoothness, par = min.out$par,
objective = min.out$objective, convergence = min.out$convergence,
message = min.out$message, iterations = min.out$iter,
evaluations = min.out$eval))
}
logLik.cov.sp <- function(par, X, y, D, sp.type, reml = FALSE, minus2 = FALSE)
{
range.par <- par[1]
error.ratio <- par[2]
if(sp.type == "matern")
{
smoothness <- par[3]
}else
{
smoothness <- 0.5
}
V <- simple.cov.sp(D = D, sp.type = sp.type, sp.par = c(1, range.par),
error.var = error.ratio, smoothness = smoothness, finescale.var = 0)
ViX <- solve(V, X)
XtViX <- crossprod(ViX, X)
fitted <- X %*% solve(XtViX, crossprod(ViX, y))
r <- y - fitted
n <- length(y)
df.resid <- n - ncol(X)
rtVir <- crossprod(r, solve(V,r))
if(!reml)
{
sigmasq <- c(rtVir/n)
lik <- determinant(V, logarithm=TRUE)$modulus +
n + length(r) * log(2*pi) + length(r)*log(sigmasq)
}else
{
sigmasq <- c(rtVir/df.resid)
lik <- determinant(V, logarithm=TRUE)$modulus +
determinant(XtViX, logarithm=TRUE)$modulus +
(n - ncol(X)) + n * log(2*pi) + (n - ncol(X))*log(sigmasq)
}
if(!minus2)
{
lik <- -lik/2
}
return(lik)
}
maxlik.cov.st <- function(X, y, coords, time, sp.type = "exponential",
range.par = stop("specify range.par argument"),
error.ratio = stop("specify error.ratio argument"),
smoothness = 0.5, t.type = "ar1", t.par = .5, D = NULL, T = NULL,
reml = TRUE, lower = NULL, upper = NULL, control = list(trace = TRUE), optimizer="nlminb")
{
y <- as.vector(y)
maxlik_cov_st_check_arg(X = X, y = y, coords = coords, time = time,
sp.type = sp.type, range.par = range.par, error.ratio = error.ratio,
smoothness = smoothness, t.type = t.type,
D = D, T = T, reml = reml, lower = lower, upper = upper)
if(is.null(D)){ D <- dist1(coords) }
if(is.null(T)){ T <- dist1(matrix(time)) }
#construct vector of initial parameters to pass to objective function
parms <- c(range.par, error.ratio)
if(sp.type == "matern"){ parms <- c(parms, smoothness) }
parms <- c(parms, t.par)
#define lower and upper bounds if not supplied by user
bounds <- bounds.cov.sp(sp.type)
if(is.null(lower)){ lower <- bounds$lower.bound }
if(is.null(upper)){ upper <- bounds$upper.bound }
lower <- c(lower, 0)
upper <- c(upper, .999)
if(optimizer=="nlminb")
{
min.out <- stats::nlminb(start=parms,
objective = logLik.cov.st, X=X, y=y, D=D,
sp.type = sp.type, t.type = t.type, T=T,
reml=reml, minus2 = TRUE,
lower=lower, upper=upper, control = control)
range.par <- min.out$par[1]
error.ratio <- min.out$par[2]
if(sp.type == "matern")
{
smoothness <- min.out$par[3]
t.par <- min.out$par[4]
}else
{
t.par <- min.out$par[3]
}
}
V.sp <- simple.cov.sp(D, sp.type, sp.par = c(1, range.par),
error.var = 0, smoothness, finescale.var = 0)
V.time <- simple.cov.time(T, t.type, t.par)
V <- V.sp * V.time + diag(nrow(V.sp)) * error.ratio
n <- length(y)
ViX <- solve(V, X)
XtViX <- crossprod(ViX, X)
fitted <- X %*% solve(XtViX, crossprod(ViX, y))
r <- y - fitted
df.resid <- n - ncol(X)
rtVir <- crossprod(r, solve(V,r))
if(!reml)
{
sigmasq <- as.vector(rtVir/n)
}else
{
sigmasq <- as.vector(rtVir/df.resid)
}
error.var <- sigmasq * error.ratio
sp.par <- c(sigmasq, range.par)
return(list(sp.type = sp.type, sp.par = sp.par,
error.var = error.var, t.par = t.par,
smoothness = smoothness, par = min.out$par,
objective = min.out$objective, convergence = min.out$convergence,
message = min.out$message, iterations = min.out$iter,
evaluations = min.out$eval))
}
logLik.cov.st <- function(par, X, y, D, T, sp.type, t.type, reml = FALSE, minus2 = FALSE)
{
range.par <- par[1]
error.ratio <- par[2]
if(sp.type == "matern")
{
smoothness <- par[3]
t.par <- par[4]
}else
{
smoothness <- 0.5
t.par <- par[3]
}
V.sp <- simple.cov.sp(D, sp.type, sp.par = c(1, range.par),
error.var = 0, smoothness, finescale.var = 0)
V.time <- simple.cov.time(T, t.type, t.par)
V <- V.sp * V.time + diag(nrow(V.sp)) * error.ratio
ViX <- solve(V, X)
XtViX <- crossprod(ViX, X)
fitted <- X %*% solve(XtViX, crossprod(ViX, y))
r <- y - fitted
n <- length(y)
df.resid <- n - ncol(X)
rtVir <- crossprod(r, solve(V,r))
if(!reml)
{
sigmasq <- c(rtVir/n)
lik <- determinant(V, logarithm=TRUE)$modulus +
n + length(r) * log(2*pi) + length(r)*log(sigmasq)
}else
{
sigmasq <- c(rtVir/df.resid)
lik <- determinant(V, logarithm=TRUE)$modulus +
determinant(XtViX, logarithm=TRUE)$modulus +
(n - ncol(X)) + n * log(2*pi) + (n - ncol(X))*log(sigmasq)
}
if(!minus2)
{
lik <- -lik/2
}
return(lik)
}
bounds.cov.sp <- function(sp.type)
{
lower <- c(.001, 0)
upper <- c(Inf, 1)
if(sp.type == "matern")
{
lower <- c(lower, .001)
upper <- c(upper, 10)
}
return(list(lower.bound = lower, upper.bound = upper))
}
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.