[spcopula-commits] r124 - in pkg: R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 13 21:05:43 CET 2014
Author: ben_graeler
Date: 2014-02-13 21:05:42 +0100 (Thu, 13 Feb 2014)
New Revision: 124
Modified:
pkg/R/Classes.R
pkg/R/spCopula.R
pkg/R/spVineCopula.R
pkg/R/spatialGaussianCopula.R
pkg/R/spatialPreparation.R
pkg/R/spatio-temporalPreparation.R
pkg/R/stVineCopula.R
pkg/demo/spCopula.R
pkg/man/calcSpTreeDists.Rd
pkg/man/getNeighbours.Rd
pkg/man/getStNeighbours.Rd
pkg/man/neighbourhood-class.Rd
pkg/man/neighbourhood.Rd
pkg/man/spCopPredict.Rd
pkg/man/spGaussCopPredict.Rd
pkg/man/spGaussLogLik.Rd
pkg/man/spVineCopula-class.Rd
pkg/man/spVineCopula.Rd
pkg/man/stCopPredict.Rd
pkg/man/stNeighbourhood-class.Rd
pkg/man/stNeighbourhood.Rd
Log:
- redesign of spatial and spatio-temporal neighbourhoods (make them lighter: skiped several slots)
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2014-02-13 09:38:55 UTC (rev 123)
+++ pkg/R/Classes.R 2014-02-13 20:05:42 UTC (rev 124)
@@ -173,9 +173,9 @@
t.res="character"),
validity = validStCopula, contains = list("copula"))
-####################
-## vine copulas ##
-####################
+###############################################
+## vine copulas, happens now in VineCopula ##
+###############################################
# validVineCopula = function(object) {
# dim <- object at dimension
@@ -235,77 +235,58 @@
## neighbourhood:
-sizeLim <- 25 # a constant
-# setSizeLim <- function(x) {
-# env <- parent.env(environment())
-# unlockBinding("neighbourLim",env)
-# assign("neighbourLim", x,envir=env)
-# lockBinding("neighbourLim",env)
-# }
-
-# a class combining two matrices holding the data and the corresponding
-# distances as well a slot for the coordinates refernce system and an attribute
-# if the data is already transformed to uniform on [0,1] distributed variables
-# data: a list of data.frames holding the data per neighbour. each neighbour needs to have the same number of variables in the same order
-# sp: an optional slot providing the coordinates of locations
-# index: a matrix linking the data entries with the coordinates of the locations
validNeighbourhood <- function(object) {
if(length(var)>1)
return("Only a single variable name is supported.")
- if (object at prediction & is.null(object at predLocs))
- return("The prediction locations have to be provided for the prediction procedure.")
# check for number of rows
if (nrow(object at data) != nrow(object at distances))
return("Data and distances have unequal number of rows.")
if (nrow(object at data) != nrow(object at index))
return("Data and index have unequal number of rows.")
# check for columns
- if (ncol(object at data) != ncol(object at distances)+1)
+ if (ncol(object at data) != ncol(object at distances) + 1 + length(object at coVar))
return("Data and distances have non matching number of columns.")
- if (ncol(object at data) != ncol(object at index))
- return("Data and index have unequal number of columns.")
+ if (ncol(object at data) != ncol(object at index) + length(object at coVar))
+ return("Data and index have non matching number of columns.")
else
return(TRUE)
}
-setClassUnion("optionalLocs",c("NULL","Spatial"))
-
setClass("neighbourhood",
representation = representation(data = "data.frame",
distances="matrix",
index="matrix",
- dataLocs="Spatial",
- predLocs="optionalLocs",
- prediction="logical",
- var="character"),
- validity = validNeighbourhood, contains = list("Spatial"))
+ var="character",
+ coVar="character",
+ prediction="logical"),
+ validity = validNeighbourhood)
## ST neighbourhood
-
validStNeighbourhood <- function(object) {
- sizeN <- nrow(object at data)
- if (object at prediction & is.null(object at dataLocs))
- return("The spatio-temporal locations of the data have to be provided for the estimation procedure.")
dimDists <- dim(object at distances)
+ dimInd <- dim(object at index)
+
+ stopifnot(length(dimDists)==3)
+ stopifnot(length(dimInd)==3)
+ stopifnot(dimDists[3] == dimInd[3])
+
if (nrow(object at data) != dimDists[1])
return("Data and distances have unequal number of rows.")
- dimInd <- dim(object at index)
- if (nrow(object at data) != dimInd[1])
+ if (nrow(object at data) != dimInd[1])
return("Data and index have unequal number of rows.")
- if (dimDists[2] != dimInd[2])
- return("Data and index have unequal number of columns.")
+ if (ncol(object at data) + object at prediction != dimInd[2] + length(object at coVar))
+ return("Data and index have non matching number of columns.")
+ if (dimDists[2]+1 != dimInd[2])
+ return("Data and index have non matching number of columns.")
else
return(TRUE)
}
-setClassUnion("optionalST",c("NULL","ST"))
-
setClass("stNeighbourhood",
representation = representation(data = "data.frame",
distances="array",
index="array",
- locations="ST",
- dataLocs="optionalST",
var="character",
+ coVar="character",
prediction="logical"),
- validity = validStNeighbourhood, contains = list("ST"))
\ No newline at end of file
+ validity = validStNeighbourhood)
\ No newline at end of file
Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2014-02-13 09:38:55 UTC (rev 123)
+++ pkg/R/spCopula.R 2014-02-13 20:05:42 UTC (rev 124)
@@ -657,15 +657,16 @@
cond <- suppressWarnings(as.numeric(varSplit[length(varSplit)]))
if(is.na(cond)) {
var <- paste(neigh at var,"|0",sep="")
+ coVar <- paste(neigh at coVar,"|0",sep="")
colnames(u1) <- paste(paste("N", rep(1:(ncol(u1)), each=length(var)), sep=""),
rep(var,ncol(u1)),sep=".")
} else {
var <- paste(neigh at var,cond+1,sep="")
+ coVar <- neigh at coVar
colnames(u1) <- paste(paste("N", rep(cond:(ncol(u1)+cond-1)+2,
each=length(var)), sep=""),
rep(var,ncol(u1)),sep=".")
}
return(neighbourhood(data=u1, distances=h1, index=neigh at index[,-1],
- dataLocs=neigh at dataLocs, predLocs=neigh at predLocs,
- prediction=neigh at prediction, var=var))
+ var=var, coVar=coVar, prediction=neigh at prediction))
}
\ No newline at end of file
Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R 2014-02-13 09:38:55 UTC (rev 123)
+++ pkg/R/spVineCopula.R 2014-02-13 20:05:42 UTC (rev 124)
@@ -86,13 +86,18 @@
dspVine(as.matrix(u), copula at spCop, NULL, log=log, ...)
})
-# fiiting the spatial vine for a given list of spatial copulas
-fitSpVine <- function(copula, data, method, estimate.variance=F) {
- stopifnot(class(data)=="neighbourhood")
- stopifnot(copula at dimension == ncol(data at data))
+# fitting the spatial vine for a given list of spatial copulas
+fitSpVine <- function(copula, data, method, estimate.variance=FALSE) {
+ stopifnot(is.list(data))
+ stopifnot(length(data)==2)
+ neigh <- data[[1]]
+ dataLocs <- data[[2]]
- u0 <- as.matrix(data at data) # previous level's (contitional) data
- h0 <- data at distances # previous level's distances
+ stopifnot(class(neigh)=="neighbourhood")
+ stopifnot(copula at dimension == ncol(neigh at data))
+
+ u0 <- as.matrix(neigh at data) # previous level's (contitional) data
+ h0 <- neigh at distances # previous level's distances
l0 <- rep(0,nrow(u0)) # spatial density
for(spTree in 1:length(copula at spCop)) {
cat("[Dropping ", spTree, ". spatial tree.]\n",sep="")
@@ -101,9 +106,9 @@
for(i in 1:ncol(h0)) { # i <- 1
l0 <- l0 + dCopula(u0[,c(1,i+1)], copula at spCop[[spTree]], h=h0[,i], log=T)
u1 <- cbind(u1, dduCopula(u0[,c(1,i+1)], copula at spCop[[spTree]], h=h0[,i]))
- if (i < ncol(h0)) {
- h1 <- cbind(h1,apply(data at index[,c(spTree+1,spTree+i+1)],1,
- function(x) spDists(data at dataLocs[x,])[1,2]))
+ if (i < ncol(h0) & spTree < length(copula at spCop)) {
+ h1 <- cbind(h1,apply(neigh at index[,c(spTree+1,spTree+i+1)],1,
+ function(x) spDists(dataLocs[x,])[1,2]))
}
}
u0 <- u1
@@ -135,13 +140,13 @@
method = sapply(method,paste,collapse=", "),
loglik = sum(l0)+loglik,
fitting.stats=list(convergence = as.integer(NA)),
- nsample = nrow(data at data), copula=spVineCop))
+ nsample = nrow(neigh at data), copula=spVineCop))
}
-setMethod("fitCopula",signature=signature("spVineCopula"),fitSpVine)
+setMethod("fitCopula", signature=signature("spVineCopula"), fitSpVine)
# deriving all spatial tree distances
-calcSpTreeDists <- function(neigh, n.trees) {
+calcSpTreeDists <- function(neigh, dataLocs, n.trees) {
condDists <- list(n.trees)
condDists[[1]] <- neigh at distances
if(n.trees==1)
@@ -150,7 +155,7 @@
h1 <- NULL
for(i in 1:(ncol(neigh at distances)-spTree)) {
h1 <- cbind(h1,apply(neigh at index[,c(spTree+1,spTree+i+1),drop=F],1,
- function(x) spDists(neigh at dataLocs[x,])[1,2]))
+ function(x) spDists(dataLocs[x,])[1,2]))
dimnames(h1) <- NULL
}
condDists[[spTree+1]] <- h1
@@ -185,11 +190,13 @@
# interpolation
-spCopPredict.expectation <- function(predNeigh, spVine, margin, ..., stop.on.error=F) {
+spCopPredict.expectation <- function(data, spVine, margin, ..., stop.on.error=F) {
stopifnot(is.function(margin$q))
- dists <- calcSpTreeDists(predNeigh,length(spVine at spCop))
+ predNeigh <- data[[1]]
+ dists <- calcSpTreeDists(predNeigh, data[[2]], length(spVine at spCop))
+
predMean <- NULL
pb <- txtProgressBar(0,nrow(predNeigh at data), 0, width=getOption("width")-10, style=3)
@@ -210,21 +217,24 @@
}
close(pb)
- if ("data" %in% slotNames(predNeigh at predLocs)) {
+ predLocs <- data[[3]]
+ if ("data" %in% slotNames(predLocs)) {
res <- predNeigh at predLocs
res at data[["expect"]] <- predMean
return(res)
} else {
predMean <- data.frame(predMean)
colnames(predMean) <- "expect"
- return(addAttrToGeom(predNeigh at predLocs, predMean, match.ID=FALSE))
+ return(addAttrToGeom(predLocs, predMean, match.ID=FALSE))
}
}
-spCopPredict.quantile <- function(predNeigh, spVine, margin, p=0.5) {
+spCopPredict.quantile <- function(data, spVine, margin, p=0.5) {
stopifnot(is.function(margin$q))
- dists <- calcSpTreeDists(predNeigh,length(spVine at spCop))
+ predNeigh <- data[[1]]
+ dists <- calcSpTreeDists(predNeigh, data[[2]], length(spVine at spCop))
+
predQuantile <- NULL
pb <- txtProgressBar(0, nrow(predNeigh at data), 0, width=getOption("width")-10, style=3)
for(i in 1:nrow(predNeigh at data)) { # i <-1
@@ -241,31 +251,29 @@
b <- density[lower]
xRes <- -b/m+sign(m)*sqrt(b^2/m^2+2*(p-int[lower])/m)
-# pPred <- optimise(function(x) abs(integrate(condSecVine, 0, x,
-# subdivisions=10000L,
-# abs.tol=1e-6)$value-p), c(0,1))
-# if(pPred$objective > 1e-4)
-# warning("Numerical evaluation in predQuantile achieved an objective of only ",
-# pPred$objective, " where 0 has been sought for location ",i,".")
predQuantile <- c(predQuantile, margin$q(xVals[lower]+xRes))
}
close(pb)
- if ("data" %in% slotNames(predNeigh at predLocs)) {
- res <- predNeigh at predLocs
+ predLocs <- data[[3]]
+ if ("data" %in% slotNames(predLocs)) {
+ res <- predLocs
res at data[[paste("quantile.",p,sep="")]] <- predQuantile
return(res)
} else {
predQuantile <- data.frame(predQuantile)
colnames(predQuantile) <- paste("quantile.",p,sep="")
- return(addAttrToGeom(predNeigh at predLocs, predQuantile, match.ID=FALSE))
+ return(addAttrToGeom(predLocs, predQuantile, match.ID=FALSE))
}
}
-spCopPredict <- function(predNeigh, spVine, margin, method="quantile", p=0.5, ...) {
+spCopPredict <- function(data, spVine, margin, method="quantile", p=0.5, ...) {
+ stopifnot(is.list(data))
+ stopifnot(length(data)==3)
+
switch(method,
- quantile=spCopPredict.quantile(predNeigh, spVine, margin, p),
- expectation=spCopPredict.expectation(predNeigh, spVine, margin, ...))
+ quantile=spCopPredict.quantile(data, spVine, margin, p),
+ expectation=spCopPredict.expectation(data, spVine, margin, ...))
}
# draw from a spatial vine
Modified: pkg/R/spatialGaussianCopula.R
===================================================================
--- pkg/R/spatialGaussianCopula.R 2014-02-13 09:38:55 UTC (rev 123)
+++ pkg/R/spatialGaussianCopula.R 2014-02-13 20:05:42 UTC (rev 124)
@@ -1,10 +1,10 @@
## spatial Gaussian Copula
# "density" evaluation
-spGaussLogLik <- function(corFun, neigh, log=T) {
+spGaussLogLik <- function(corFun, neigh, dataLocs, log=T) {
neighDim <- ncol(neigh at data)
- allDataDists <- spDists(neigh at dataLocs)
+ allDataDists <- spDists(dataLocs)
pb <- txtProgressBar(0, nrow(neigh at data), 0, width = getOption("width") - 10, style = 3)
@@ -29,12 +29,12 @@
}
# interpolation based on a valid corelogram function
-spGaussCopPredict <- function(corFun, predNeigh, margin, p=0.5, ..., n=1000) {
+spGaussCopPredict <- function(corFun, predNeigh, dataLocs, predLocs, margin, p=0.5, ..., n=1000) {
stopifnot(is.list(margin))
stopifnot(is(margin$q, "function"))
neighDim <- ncol(predNeigh at data)
- allDataDists <- spDists(predNeigh at dataLocs)
+ allDataDists <- spDists(dataLocs)
pb <- txtProgressBar(0, nrow(predNeigh at data), 0, width = getOption("width") - 10, style = 3)
@@ -73,15 +73,15 @@
}
close(pb)
- if ("data" %in% slotNames(predNeigh at predLocs)) {
- res <- predNeigh at predLocs
+ if ("data" %in% slotNames(predLocs)) {
+ res <- predLocs
res at data[[paste("quantile.", p, sep = "")]] <- predQuantile
return(res)
}
else {
predQuantile <- data.frame(predQuantile)
colnames(predQuantile) <- paste("quantile.", p, sep = "")
- return(addAttrToGeom(predNeigh at predLocs, predQuantile,
+ return(addAttrToGeom(predLocs, predQuantile,
match.ID = FALSE))
}
}
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2014-02-13 09:38:55 UTC (rev 123)
+++ pkg/R/spatialPreparation.R 2014-02-13 20:05:42 UTC (rev 124)
@@ -7,110 +7,59 @@
## spatial neighbourhood constructor
####################################
-neighbourhood <- function(data, distances, index, dataLocs, predLocs=NULL,
- prediction, var) {
+neighbourhood <- function(data, distances, index, var, coVar=character(), prediction=FALSE) {
sizeN <- ncol(distances)+1
data <- as.data.frame(data)
if (anyDuplicated(rownames(data))>0)
rownames <- 1:length(rownames)
+
new("neighbourhood", data=data, distances=distances, index=index,
- dataLocs=dataLocs, predLocs=predLocs, prediction=prediction, var=var,
- bbox=dataLocs at bbox, proj4string=dataLocs at proj4string)
+ var=var, coVar=coVar, prediction=prediction)
}
## show
showNeighbourhood <- function(object){
cat("A set of neighbourhoods consisting of", ncol(object at distances)+1, "locations each \n")
- cat("with",nrow(object at data),"rows of observations for:\n")
- cat(object at var,"\n")
+ if (length(object at var)>0) {
+ cat("with", nrow(object at data), "rows of observations for:\n")
+ cat(object at var, "\n")
+ } else {
+ cat("without data \n")
+ }
+ if(length(object at coVar)>0)
+ cat("with covariate", object at coVar, "\n")
}
-setMethod(show,signature("neighbourhood"),showNeighbourhood)
+setMethod(show, signature("neighbourhood"), showNeighbourhood)
## names (from sp)
-setMethod(names, signature("neighbourhood"), function(x) x at var)
+setMethod(names, signature("neighbourhood"), function(x) c(x at var,x at coVar))
-## spplot ##
-spplotNeighbourhood <- function(obj, zcol=names(obj), ..., column=0) {
- stopifnot(all(column<ncol(obj at data)))
- pattern <- paste(paste("N", column, ".", sep=""), zcol, sep="")
- spdf <- SpatialPointsDataFrame(coords=obj at dataLocs,
- data=obj at data[,pattern,drop=FALSE],
- proj4string=obj at proj4string, bbox=obj at bbox)
- spplot(spdf, ...)
-}
-
-setMethod(spplot, signature("neighbourhood"), spplotNeighbourhood)
-
selectFromNeighbourhood <- function(x, i) {
- newSp <- x at dataLocs[i,]
new("neighbourhood", data=x at data[i,,drop=F],
- distances=x at distances[i,,drop=F],
- dataLocs=newSp, predLocs=x at predLocs, bbox=newSp at bbox,
- proj4string=x at proj4string, index=x at index[i,,drop=F], var=x at var,
- prediction=x at prediction)
+ distances=x at distances[i,,drop=F], index=x at index[i,,drop=F],
+ var=x at var, coVar=x at coVar, prediction=x at prediction)
}
-setMethod("[[", signature("neighbourhood","numeric","missing"), selectFromNeighbourhood)
+setMethod("[", signature("neighbourhood","numeric"), selectFromNeighbourhood)
## calculate neighbourhood from SpatialPointsDataFrame
-
-# returns an neighbourhood object
-##################################
-
-# getNeighbours <- function(dataLocs, predLocs, var=names(dataLocs)[1], size=5,
-# prediction=FALSE, min.dist=0.01) {
-#
-# stopifnot((!prediction && missing(predLocs)) || (prediction && !missing(predLocs)))
-# stopifnot(min.dist>0 || prediction)
-#
-# if(missing(predLocs) && !prediction)
-# predLocs=dataLocs
-#
-# stopifnot(is(predLocs,"Spatial"))
-#
-# if(any(is.na(match(var,names(dataLocs)))))
-# stop("At least one of the variables is unkown or is not part of the data.")
-#
-# nLocs <- length(predLocs)
-# size <- min(size, length(dataLocs)+prediction)
-#
-# allLocs <- matrix(NA,nLocs,size)
-# allDists <- matrix(NA,nLocs,size-1)
-# allData <- matrix(NA,nLocs,size)
-# for (i in 1:nLocs) {
-# tempDists <- spDists(dataLocs, predLocs[i, ])
-# tempDists[tempDists < min.dist] <- Inf
-# spLocs <- order(tempDists)[1:(size - 1)]
-#
-# allLocs[i,] <- c(i, spLocs)
-# allDists[i,] <- tempDists[spLocs]
-# allData[i,(prediction+1):size] <- dataLocs[c(i[!prediction],spLocs),
-# var, drop = F]@data[[1]]
-# }
-#
-# if (!prediction)
-# predLocs <- NULL
-# colnames(allData) <- paste(paste("N", rep(0:(size-1),
-# each=length(var)), sep=""),
-# rep(var,size),sep=".")
-# return(neighbourhood(allData, allDists, allLocs, dataLocs,
-# predLocs, prediction, var))
-# }
-#
-
-
-getNeighbours <- function (dataLocs, predLocs, var = names(dataLocs)[1], size = 5,
- prediction = FALSE, min.dist = 0.01)
-{
- stopifnot((!prediction && missing(predLocs)) || (prediction &&
- !missing(predLocs)))
+getNeighbours <- function (dataLocs, predLocs, size = 5,
+ var = names(dataLocs)[1], coVar=character(),
+ prediction = FALSE, min.dist = 0.01) {
+ stopifnot((!prediction && missing(predLocs)) || (prediction && !missing(predLocs)))
stopifnot(min.dist > 0 || prediction)
+
if (missing(predLocs) && !prediction)
predLocs = dataLocs
+
stopifnot(is(predLocs, "Spatial"))
- if (any(is.na(match(var, names(dataLocs)))))
- stop("At least one of the variables is unkown or is not part of the data.")
+
+ if ("data" %in% slotNames(dataLocs)) {
+ if (any(is.na(match(var, names(dataLocs)))))
+ stop("The variables is not part of the data.")
+ }
+
nLocs <- length(predLocs)
size <- min(size, length(dataLocs) + prediction)
@@ -128,18 +77,25 @@
as.integer(prediction),
PACKAGE="spcopula")
- if (!prediction) allData <- matrix(dataLocs[result$allLocs, var, drop = F]@data[[1]], nLocs, size)
- else allData <- matrix(c(rep(NA,nLocs),dataLocs[result$allLocs[(nLocs+1):(nLocs*size)], var, drop = F]@data[[1]]), nLocs, size)
+ if ("data" %in% slotNames(dataLocs)) {
+ if (!prediction) {
+ allData <- matrix(dataLocs[result$allLocs, var, drop = F]@data[[1]],
+ nLocs, size)
+ } else {
+ allData <- matrix(c(rep(NA,nLocs),
+ dataLocs[result$allLocs[(nLocs+1):(nLocs*size)], var, drop = F]@data[[1]]),
+ nLocs, size)
+ }
+ colnames(allData) <- paste(paste("N", rep(0:(size - 1), each = length(var)), sep = ""),
+ rep(var, size), sep = ".")
+ } else {
+ allData <- as.data.frame(matrix(NA, nLocs, size + length(coVar)))
+ var <- character()
+ }
- if (!prediction)
- predLocs <- NULL
-
- colnames(allData) <- paste(paste("N", rep(0:(size - 1), each = length(var)),
- sep = ""), rep(var, size), sep = ".")
-
-
- return(neighbourhood(allData, matrix(result$allDists,nLocs, size - 1), matrix(result$allLocs, nLocs, size), dataLocs,
- predLocs, prediction, var))
+ return(neighbourhood(data=allData, distances=matrix(result$allDists, nLocs, size - 1),
+ index=matrix(result$allLocs, nLocs, size), var=var, coVar=coVar,
+ prediction=prediction))
}
#############
@@ -272,84 +228,4 @@
return(res)
}
-setMethod(calcBins, signature="neighbourhood", calcNeighBins)
-
-# instances: number -> number of randomly choosen temporal intances
-# NA -> all observations
-# other -> temporal indexing as in spacetime/xts, the parameter t.lags is set to 0 in this case.
-# t.lags: numeric -> temporal shifts between obs
-calcStBins <- function(data, var, nbins=15, boundaries=NA, cutoff=NA,
- instances=NA, t.lags=-(0:2), ...,
- cor.method="fasttau", plot=FALSE) {
- if(is.na(cutoff))
- cutoff <- spDists(coordinates(t(data at sp@bbox)))[1,2]/3
- if(is.na(boundaries))
- boundaries <- ((1:nbins) * cutoff / nbins)
- if(is.na(instances))
- instances=length(data at time)
-
- spIndices <- calcSpLagInd(data at sp, boundaries)
-
- mDists <- sapply(spIndices,function(x) mean(x[,3]))
-
- lengthTime <- length(data at time)
- if (!is.numeric(instances) | !length(instances)==1) {
- tempIndices <- cbind(instances, instances)
- }
- else {
- tempIndices <- NULL
- for (t.lag in rev(t.lags)) {
-# smplInd <- max(1,1-min(t.lags)):min(lengthTime,lengthTime-min(t.lags))
- smplInd <- sample(x=max(1,1-min(t.lags)):min(lengthTime,lengthTime-min(t.lags)),
- size=min(instances,lengthTime-max(abs(t.lags))))
- tempIndices <- cbind(smplInd+t.lag, tempIndices)
- tempIndices <- cbind(smplInd, tempIndices)
- }
- }
-
- retrieveData <- function(spIndex, tempIndices) {
- binnedData <- NULL
- for (i in 1:(ncol(tempIndices)/2)) {
- binnedData <- cbind(binnedData,
- as.matrix((cbind(data[spIndex[,1], tempIndices[,2*i-1], var]@data,
- data[spIndex[,2], tempIndices[,2*i], var]@data))))
- }
- return(binnedData)
- }
-
- lagData <- lapply(spIndices, retrieveData, tempIndices=tempIndices)
-
- calcStats <- function(binnedData) {
- cors <- NULL
- for(i in 1:(ncol(binnedData)/2)) {
- cors <- c(cors, cor(binnedData[,2*i-1], binnedData[,2*i], method=cor.method, use="pairwise.complete.obs"))
- }
- return(cors)
- }
-
- calcTau <- function(binnedData) {
- cors <- NULL
- for(i in 1:(ncol(binnedData)/2)) {
- tmpData <- binnedData[,2*i+c(-1,0)]
- tmpData <- tmpData[!apply(tmpData, 1, function(x) any(is.na(x))),]
- cors <- c(cors, TauMatrix(tmpData)[1,2])
- }
- return(cors)
- }
-
- calcCor <- switch(cor.method, fasttau=calcTau, calcStats)
-
- lagCor <- sapply(lagData, calcCor)
-
- if(plot) {
- plot(mDists, as.matrix(lagCor)[1,], xlab="distance",ylab=paste("correlation [",cor.method,"]",sep=""),
- ylim=1.05*c(-abs(min(lagCor)),max(lagCor)), xlim=c(0,max(mDists)))
- abline(h=c(-min(lagCor),0,min(lagCor)),col="grey")
- }
-
- res <- list(meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=list(sp=spIndices, time=tempIndices))
- attr(res,"cor.method") <- cor.method
- return(res)
-}
-
-setMethod(calcBins, signature(data="STFDF"), calcStBins)
+setMethod(calcBins, signature="neighbourhood", calcNeighBins)
\ No newline at end of file
Modified: pkg/R/spatio-temporalPreparation.R
===================================================================
--- pkg/R/spatio-temporalPreparation.R 2014-02-13 09:38:55 UTC (rev 123)
+++ pkg/R/spatio-temporalPreparation.R 2014-02-13 20:05:42 UTC (rev 124)
@@ -7,52 +7,79 @@
## spatio-temporal neighbourhood constructor
############################################
-stNeighbourhood <- function(data, distances, STxDF, ST=NULL,index,
- prediction, var) {
+stNeighbourhood <- function(data, distances, index, var, coVar=character(), prediction=FALSE) {
data <- as.data.frame(data)
sizeN <- nrow(data)
+
dimDists <- dim(distances)
+ dimInd <- dim(index)
- stopifnot(dimDists[1]==sizeN)
- stopifnot(dimDists[2]==ncol(data)-(!prediction))
- stopifnot(dimDists[3]==2)
- colnames(data) <- paste(paste("N", (0+prediction):dimDists[2], sep=""),var,sep=".")
+ stopifnot(length(dimDists) == 3)
+ stopifnot(length(dimInd) == 3)
+
+ stopifnot(dimDists[1] == sizeN)
+ stopifnot(dimInd[1] == dimDists[1])
+
+ stopifnot(((dimDists[2] + !prediction) + length(coVar)) == ncol(data))
+ stopifnot(dimInd[2] == dimDists[2]+1)
+
+ stopifnot(dimDists[3] == 2)
+ stopifnot(dimInd[3] == dimDists[3])
+
+ colnames(data) <- paste(paste("N", (0+prediction):dimDists[2], sep=""), var, sep=".")
+
+ if(length(coVar)>0)
+ colnames(data)[dimDists[2] + 1:length(coVar)] <- paste("N0", coVar)
+
if (anyDuplicated(rownames(data))>0)
rownames <- 1:length(rownames)
- new("stNeighbourhood", data=data, distances=distances, locations=ST,
- dataLocs=STxDF, index=index, prediction=prediction, var=var,
- sp=as(STxDF at sp, "Spatial"), time=STxDF at time[1],
- endTime=STxDF at endTime[length(STxDF at endTime)])
+
+ new("stNeighbourhood", data=data, distances=distances, index=index,
+ var=var, coVar=coVar, prediction=prediction)
}
## show
-showStNeighbourhood <- function(object){
+showStNeighbourhood <- function(object) {
cat("A set of spatio-temporal neighbourhoods consisting of", dim(object at distances)[2]+1, "locations each \n")
cat("with",nrow(object at data),"rows of observations for:\n")
cat(object at var,"\n")
+ if(length(object at coVar > 0))
+ cat("with covariate", object at coVar, "\n")
}
-setMethod(show,signature("stNeighbourhood"),showStNeighbourhood)
+setMethod(show, signature("stNeighbourhood"), showStNeighbourhood)
+# select
+selectFromStNeighbourhood <- function(x, i) {
+ new("stNeighbourhood", data=x at data[i,,drop=F],
+ distances=x at distances[i,,,drop=F], index=x at index[i,,,drop=F],
+ var=x at var, coVar=x at coVar, prediction=x at prediction)
+}
+setMethod("[", signature("stNeighbourhood","numeric"), selectFromStNeighbourhood)
+
## calculate neighbourhood from ST
# returns an neighbourhood object
##################################
-getStNeighbours <- function(stData, ST, var=names(stData at data)[1], spSize=4,
- t.lags=-(0:2), timeSteps=NA, prediction=FALSE, min.dist=0.01) {
+getStNeighbours <- function(stData, ST, spSize=4, t.lags=-(0:2),
+ var=names(stData at data)[1], coVar=character(),
+ timeSteps=NA, prediction=FALSE, min.dist=0.01) {
stopifnot((!prediction && missing(ST)) || (prediction && !missing(ST)))
stopifnot(min.dist>0 || prediction)
timeSpan <- min(t.lags)
if(missing(ST) && !prediction)
- ST=stData
+ ST=geometry(stData)
stopifnot(is(ST,"ST"))
if(any(is.na(match(var,names(stData at data)))))
- stop("At least one of the variables is unkown or is not part of the data.")
+ stop("The variables is not part of stData.")
+ if(length(coVar)>0)
+ if(any(is.na(match(coVar,names(stData at data)))))
+ stop("The covariate is not part of stData.")
if(!prediction) {
if(is.na(timeSteps)) {
@@ -62,58 +89,73 @@
reSample <- function() sort(sample((1-timeSpan):length(stData at time), timeSteps))
}
nLocs <- length(ST at sp)*timeSteps
- nghbrs <- getNeighbours(stData[,1], var=var, size=spSize, min.dist=min.dist)
+ nghbrs <- getNeighbours(dataLocs=geometry(stData at sp), var=character(), size=spSize,
+ min.dist=min.dist)
} else {
nLocs <- length(ST)
- nghbrs <- getNeighbours(stData[,1], ST at sp, var, spSize, prediction, min.dist)
+ nghbrs <- getNeighbours(dataLocs=geometry(stData at sp), predLocs=geometry(ST at sp),
+ size=spSize, var=character(), prediction=prediction,
+ min.dist=min.dist)
timeNghbrs <- sapply(index(ST at time), function(x) which(x == index(stData at time)))
reSample <- function() timeNghbrs
timeSteps <- length(stData at time)+timeSpan
}
- stNeighData <- matrix(NA, nLocs, (spSize-1)*length(t.lags)+1)
- stDists <- array(NA,c(nLocs,(spSize-1)*length(t.lags),2))
- stInd <- array(NA,c(nLocs,(spSize-1)*length(t.lags),2))
+ nStNeighs <- (spSize-1)*length(t.lags)
+ stNeighData <- matrix(NA, nLocs, nStNeighs + 1 + length(coVar))
+ stDists <- array(NA,c(nLocs, nStNeighs, 2))
+ stInd <- array(NA,c(nLocs, nStNeighs + 1, 2))
+
nTimeInst <- length(reSample())
- for(i in 1:nrow(nghbrs at index)){ # i <- 1
+ for (i in 1:nrow(nghbrs at index)) {
timeInst <- reSample() # draw random time steps for each neighbourhood
- stNeighData[(i-1)*timeSteps+(1:timeSteps),
- 1:spSize] <- matrix(stData[nghbrs at index[i,], timeInst,
- var, drop=F]@data[[1]],
- ncol=spSize, byrow=T) # retrieve the top level data
- tmpInd <- matrix(rep(timeInst, spSize-1), ncol=spSize-1)
- for(j in 2:length(t.lags)) {
+ spInd <- (i-1)*timeSteps+(1:timeSteps)
+
+ stNeighData[spInd, 1:spSize] <- matrix(stData[nghbrs at index[i,], timeInst,
+ var, drop=F]@data[[1]],
+ ncol=spSize, byrow=T)
+ # add covariate(s) to the last column(s)
+ if (length(coVar) > 0) {
+ coVarCols <- nStNeighs + 1 + (1:length(coVar))
+ stNeighData[spInd, coVarCols] <- matrix(stData[nghbrs at index[i,1], timeInst,
+ coVar, drop=F]@data[[1]],
+ ncol=length(coVar), byrow=T)
+ }
+
+ tmpInd <- matrix(rep(timeInst, spSize), ncol=spSize)
+
+ for (j in 2:length(t.lags)) {
t <- t.lags[j]
- stNeighData[(i-1)*timeSteps+(1:timeSteps),
- (j-1)*(spSize-1)+2:(spSize)] <- matrix(stData[nghbrs at index[i,][-1],
- timeInst+t,
- var, drop=F]@data[[1]],
- ncol=spSize-1, byrow=T)
- tmpInd <- cbind(tmpInd, matrix(rep(timeInst+t,spSize-1),ncol=spSize-1))
+ stNeighData[spInd, (j-1)*(spSize-1)+2:(spSize)] <- matrix(stData[nghbrs at index[i,][-1],
+ timeInst+t, var, drop=F]@data[[1]],
+ ncol=spSize-1, byrow=T)
+ tmpInd <- cbind(tmpInd, matrix(rep(timeInst+t,spSize-1), ncol=spSize-1))
}
-
- stDists[(i-1)*timeSteps+1:timeSteps,,1] <- matrix(rep(nghbrs at distances[i,],
- timeSteps*length(t.lags)),
- byrow=T, ncol=length(t.lags)*(spSize-1)) # store sp distances
- stDists[(i-1)*timeSteps+1:timeSteps,,2] <- matrix(rep(rep(t.lags,each=spSize-1),
- timeSteps),
- byrow=T, ncol=length(t.lags)*(spSize-1)) # store tmp distances
- stInd[(i-1)*timeSteps+1:timeSteps,,1] <- matrix(rep(nghbrs at index[i,][-1],
- timeSteps*length(t.lags)),
- byrow=T, ncol=length(t.lags)*(spSize-1))
- stInd[(i-1)*timeSteps+1:timeSteps,,2] <- tmpInd
+
+ # store spatial distances
+ stDists[spInd,,1] <- matrix(rep(nghbrs at distances[i,], timeSteps*length(t.lags)),
+ byrow=T, ncol=nStNeighs)
+
+ # store temporal distances
+ stDists[spInd,,2] <- matrix(rep(rep(t.lags,each=spSize-1), timeSteps),
+ byrow=T, ncol=nStNeighs)
+
+ # store space indices
+ stInd[spInd,,1] <- matrix(rep(c(nghbrs at index[i, ], rep(nghbrs at index[i, -1], length(t.lags)-1)),
+ timeSteps), ncol = nStNeighs + 1, byrow = T)
+
+ # store time indices
+ stInd[spInd,,2] <- tmpInd
}
if (prediction) {
- dataLocs <- stData
stNeighData <- stNeighData[,-1]
} else {
dataLocs <- NULL
}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/spcopula -r 124
More information about the spcopula-commits
mailing list