Nothing
corMatern <- function(value = c(1, 0.5), form = ~1, nugget=FALSE, nuScaled=FALSE,
metric = c("euclidean", "maximum", "manhattan"), fixed = FALSE) {
attr(value, "formula") <- form
attr(value, "nugget") <- nugget
attr(value, "nuScaled") <- nuScaled
attr(value, "metric") <- match.arg(metric)
attr(value, "fixed") <- fixed
class(value) <- c("corMatern", "corStruct") ## corSpatial removed
value
}
getCovariate.corMatern <- function(object, form = formula(object), data) {
if (is.null(covar <- attr(object, "covariate"))) { # need to calculate it
if (missing(data)) {
stop("need data to calculate covariate")
}
covForm <- nlme::getCovariateFormula(form)
if (length(all.vars(covForm)) > 0) { # covariate present
if (attr(terms(covForm), "intercept") == 1) {
covForm <- eval(parse(text = paste0("~", .DEPARSE(covForm[[2]]),"-1")))
}
covar <- as.data.frame(unclass(model.matrix(covForm,
model.frame(covForm, data, drop.unused.levels = TRUE))))
} else {
covar <- NULL
}
if (!is.null(nlme::getGroupsFormula(form))) { # by groups
grps <- getGroups(object, data = data)
if (is.null(covar)) {
distmats <- lapply(split(grps, grps), function(x) (.dist_fn(1:length(x))))
} else {
locfn <- function(el, metric) {
el <- as.matrix(el)
if (nrow(el) > 1) {
.dist_fn(el, metric)
} else {
numeric(0)
}
}
distmats <- lapply(split(covar, grps), locfn, metric = attr(object, "metric"))
}
covar <- lapply(distmats,as.vector)
covar <- covar[sapply(covar, length) > 0] # no 1-obs groups ## takes only list elements of nonzero length
} else { # no groups
if (is.null(covar)) {
distmats <- .dist_fn(1:nrow(data))
} else {
distmats <- .dist_fn(as.matrix(covar),
method = attr(object, "metric"))
}
covar <- lapply(distmats,as.vector)
}
if (any(unlist(covar) == 0)) {
stop("cannot have zero distances in \"corMatern\"")
}
attr(covar,"distmats") <- distmats
}
covar
}
Initialize.corMatern <- function(object, data, ...) {
if (!is.null(attr(object, "minD"))) { #already initialized
return(object)
}
form <- formula(object)
## obtaining the groups information, if any
if (!is.null(nlme::getGroupsFormula(form))) {
attr(object, "groups") <- getGroups(object, form, data = data)
attr(object, "Dim") <- nlme::Dim(object, attr(object, "groups"))
} else {
# no groups
attr(object, "Dim") <- nlme::Dim(object, as.factor(rep(1, nrow(data))))
}
## obtaining the covariate(s)
## is this where the distance matrix is actually computed ?
attr(object, "covariate") <- getCovariate.corMatern(object, data = data)
nug <- attr(object, "nugget")
#nusc <- attr(object, "nuScaled")
val <- as.vector(object) # how many parameters?
if (length(val) > 0) { # is initialized
if (val[1] <= 0) { # test for values of range
stop("'range' must be > 0 in \"corMatern\" initial value")
}
if (val[2] <= 0) { # test for values of nu
stop("'nu' must be > 0 in \"corMatern\" initial value")
}
if (nug) { # with nugget effect
if (length(val) == 2) { # assuming nugget effect not given
val <- c(val, 0.1) # setting it to 0.1
} else {
if (length(val) != 3) {
stop("initial value for \"corMatern\" parameters of wrong dimension")
}
}
if ((val[3] <= 0) || (val[3] >= 1)) {
stop("initial value of nugget ratio must be between 0 and 1")
}
} else { # only range and nu
if (length(val) != 2) {
stop("initial value for \"corMatern\" parameters of wrong dimension")
}
}
} else { ## val of length 0
val <- min(unlist(attr(object, "covariate"))) / 0.9
val <- c(val,0.5) ## FR 18/03/13
if (nug) val <- c(val, 0.1)
}
val[1] <- log(val[1]) ## log(range=1/rho)
val[2] <- log(val[2]) ## log(nu)
if (nug) val[3] <- log(val[3]/(1 - val[3]))
oldAttr <- attributes(object)
object <- val
attributes(object) <- oldAttr
attr(object, "minD") <- min(unlist(attr(object, "covariate")))
attr(object, "factor") <- corFactor.corMatern(object)
attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
object
} ## Initialize.corMatern
logDet.corMatern <- function(object, covariate = getCovariate(object), ...) {
if (!is.null(aux <- attr(object, "logDet"))) {
return(aux)
}
if (is.null(aux <- attr(object, "factor"))) {
## getting the transpose sqrt factor
aux <- corMatrix(object, covariate = covariate, corr = FALSE)
}
if (is.null(aux1 <- attr(aux, "logDet"))) {
## checking for logDet attribute; if not present, get corr matrix
aux <- corMatrix(object, covariate)
if (data.class(aux) == "list") { # by group
sum(log(abs(unlist(lapply(aux, function(el) svd(el)$d)))))/2
} else {
sum(log(abs(svd(aux)$d)))/2
}
} else {
-aux1
}
}
corMatrix.corMatern <- function(object, covariate = getCovariate(object),
corr = TRUE, ## returns matrix, else returna cholesky factor + the logDet
...) {
if (data.class(covariate) == "list") { ## groups are defined
if (is.null(names(covariate))) {
names(covariate) <- 1:length(covariate)
}
corD <- nlme::Dim(object, rep(names(covariate),
unlist(lapply(covariate,
function(el) round((1 + sqrt(1 + 8 * length(el)))/2)))))
} else { ## no groups
corD <- nlme::Dim(object, rep(1, round((1 + sqrt(1 + 8* length(covariate)))/2)))
}
par <- as.vector(object)
rho <- exp(-par[1]);
nu <- exp(par[2]);
if (attr(object, "nugget")) {
aux <- exp(par[3]);
Nugget <- 1 - 1 / (1.0 + aux)
} else Nugget <- 0.
lD <- NULL
val <- vector("list",length(covariate))
distmats <- attr(covariate,"distmats")
val <- vector("list",length(distmats))
if ( ! corr) lD <- numeric(length(covariate))
for (it in seq_along(distmats)) {
tmp <- MaternCorr(distmats[[it]], nu=nu, rho=rho,Nugget=Nugget)
tmp <- proxy::as.matrix(tmp, diag=1)
if (corr) { ## returns a list of correlation matrices
val[[it]] <- tmp
} else {
tmp <- mat_sqrt(m=tmp)
val[[it]] <- solve(tmp) # t(solve(chol(tmp)))
lD[it] <- determinant(val[[it]])$modulus
lD <- sum(lD)
}
}
if (corD[["M"]] > 1) {
val <- lapply(val, function(el) {
nel <- round(sqrt(length(el)))
array(el, c(nel, nel))
})
names(val) <- names(corD[["len"]])
val <- as.list(val)
} else {
val <- array(val, c(corD[["N"]], corD[["N"]]))
}
attr(val, "logDet") <- lD
val
} ## corMatrix.corMatern
coef.corMatern <- function(object, unconstrained = TRUE, ...) {
##cat("tata",object[1],object[2],"\n")
if (attr(object, "fixed") && unconstrained) {
return(numeric(0))
}
val <- as.vector(object)
if (length(val) == 0) { # uninitialized
return(val)
}
if (!unconstrained) {
val <- exp(val)
if (attr(object, "nugget")) val[3] <- val[3]/(1+val[3])
}
if (attr(object, "nugget")) names(val) <- c("range", "nu", "nugget")
else names(val) <- c("range","nu")
val
}
"coef<-.corMatern" <- function(object, ..., value) {
if (length(value) != length(object)) {
stop("cannot change the length of the parameter after initialization")
}
# cat("tutu",object[1],"->",value[1],object[2],"->",value[2],"\n")
object[] <- value
aux <- corFactor.corMatern(object, ...)
attr(object, "factor") <- aux
attr(object, "logDet") <- - attr(aux,"logDet")
object
} ## "coef<-.corMatern"
corFactor.corMatern <- function(object, ...) {
corD <- nlme::Dim(object)
## parameter assumed in unconstrained form = log transf for unconstrained range
par <- as.vector(object)
rho <- exp(-par[1]);
nu <- exp(par[2]);
if (attr(object, "nugget")) {
aux <- exp(par[3]);
Nugget <- 1 - 1 / (1.0 + aux)
} else Nugget <- 0.
distmats <- attr(attr(object,"covariate"),"distmats")
lD <- 0
val <- vector("list",length(distmats))
for (it in seq_along(distmats)) {
tmp <- MaternCorr(distmats[[it]], nu=nu, rho=rho,Nugget=Nugget)
tmp <- proxy::as.matrix(tmp, diag=1)
tmp <- mat_sqrt(m=tmp)
tmp <- solve(tmp) # t(solve(chol(tmp)))
lD <- lD + determinant(tmp)$modulus[[1L]]
val[[it]] <- as.vector(tmp)
}
val <- unlist(val)
attr(val, "logDet") <- lD
return(val)
} ## corFactor.corMatern
.Dim.corMatern <- function(object, groups, ...) {
if (missing(groups)) return(attr(object, "Dim"))
ugrp <- unique(groups)
groups <- factor(groups, levels = ugrp)
len <- table(groups)
val <- list(N = length(groups),
M = length(len),
maxLen = max(len),
sumLenSq = sum(len^2),
len = len,
start = NA,
spClass=0)
val[["start"]] <- c(0, 4, cumsum(val[["len"]] * (val[["len"]] - 1)/2)[-val[["M"]]])
val
}
recalc.corMatern <- function(object, conLin, ...) {
## parameter assumed in unconstrained form = log transf for unconstrained range
rho <- exp(-object[1]);
nu <- exp(object[2]);
if (attr(object, "nugget")) {
aux <- exp(object[3]);
Nugget <- 1 - 1 / (1.0 + aux) # /* 1 - nugget */
} else Nugget <- 0
covariate <- getCovariate.corMatern(object)
distmats <- attr(covariate,"distmats")
val <- vector("list",length(distmats))
cumranges <- cumsum(c(0,unlist(lapply(distmats,nrow))))
for (it in seq_along(distmats)) {
tmp <- MaternCorr(distmats[[it]], nu=nu, rho=rho,Nugget=Nugget)
tmp <- proxy::as.matrix(tmp, diag=1)
tmp <- mat_sqrt(m=tmp)
tmp <- solve(tmp) # t(solve(chol(tmp)))
yrange <- (cumranges[it]+1L):(cumranges[it+1L])
conLin[["Xy"]][yrange,] <- tmp %*% conLin[["Xy"]][yrange,,drop=FALSE]
}
conLin[["logLik"]] <- conLin[["logLik"]] +attr(attr(object,"factor"),"logDet")
conLin
} ## recalc.corMatern
Variogram.corMatern <- function(object, distance = NULL, sig2 = 1, length.out = 50, FUN, ...) {
if (is.null(distance)) {
rangeDist <- range(unlist(getCovariate(object)))
distance <- seq(rangeDist[1], rangeDist[2], length = length.out)
}
params <- coef(object, unconstrained = FALSE)
if (length(params) == 2) { # no nugget effect
range <- params[1]
nu <- params[2]
nugg <- 0
} else { # nugget effect
range <- params[1]
nu <- params[2]
nugg <- params[3]
}
par <- as.vector(object)
rho <- exp(-par[1]);
nu <- exp(par[2]);
zz <- 1-MaternCorr(distance, nu=nu, rho=rho,Nugget=0)
val <- data.frame(variog = sig2 * (nugg + (1 - nugg) * zz),
dist = distance)
class(val) <- c("Variogram", "data.frame")
val
} ## Variogram.corMatern
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.