[spcopula-commits] r99 - / pkg pkg/R pkg/demo pkg/man www
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 10 09:07:34 CEST 2013
Author: ben_graeler
Date: 2013-07-10 09:07:34 +0200 (Wed, 10 Jul 2013)
New Revision: 99
Added:
pkg/demo/pureSpVineCopula.R
Modified:
pkg/NAMESPACE
pkg/R/Classes.R
pkg/R/spCopula.R
pkg/R/spVineCopula.R
pkg/R/spatialPreparation.R
pkg/R/spatio-temporalPreparation.R
pkg/R/stCopula.R
pkg/R/stVineCopula.R
pkg/demo/00Index
pkg/demo/spCopula.R
pkg/man/calcBins.Rd
pkg/man/getNeighbours.Rd
pkg/man/getStNeighbours.Rd
pkg/man/neighbourhood-class.Rd
pkg/man/neighbourhood.Rd
pkg/man/spCopPredict.Rd
pkg/man/stCopPredict.Rd
pkg/man/stNeighbourhood.Rd
spcopula_0.1-1.tar.gz
spcopula_0.1-1.zip
www/index.php
Log:
- changed neighbourhood structure
- new, additional demo
- updated spatio-temporal structures (more to come)
- more on multi-level spatial/spatio-temporal copulas (still more to come)
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/NAMESPACE 2013-07-10 07:07:34 UTC (rev 99)
@@ -50,6 +50,4 @@
exportClasses(BB8Copula, surBB8Copula, r90BB8Copula, r270BB8Copula)
exportClasses(joeBiCopula, surJoeBiCopula, r90JoeBiCopula, r270JoeBiCopula)
exportClasses(surClaytonCopula, r90ClaytonCopula, r270ClaytonCopula)
-exportClasses(surGumbelCopula, r90GumbelCopula, r270GumbelCopula)
-
-# useDynLib("spcopula")
\ No newline at end of file
+exportClasses(surGumbelCopula, r90GumbelCopula, r270GumbelCopula)
\ No newline at end of file
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/R/Classes.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -240,34 +240,34 @@
# 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) {
- sizeN <- ncol(object at distances)+1
- nVars <- length(object at var)
- if (object at prediction & is.null(object at dataLocs))
- return("The locations of the data have to provided for the estimation procedure.")
+ 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 (ncol(object at data) %% (sizeN-object at prediction) != 0)
- return("Data and distances have non matching number of columns.")
if (nrow(object at data) != nrow(object at index))
return("Data and index have unequal number of rows.")
- if (ncol(object at distances) != ncol(object at index))
+ # check for columns
+ if (ncol(object at data) != ncol(object at distances)+1)
+ 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) != (sizeN-object at prediction) * nVars)
- return(paste("Number of columns in data does not equal the product of the neighbourhood's size (",sizeN,") with number of variables (",nVars,").",sep=""))
else
return(TRUE)
}
-setClassUnion("optionalDataLocs",c("NULL","Spatial"))
+setClassUnion("optionalLocs",c("NULL","Spatial"))
setClass("neighbourhood",
representation = representation(data = "data.frame",
distances="matrix",
index="matrix",
- locations="Spatial",
- dataLocs="optionalDataLocs",
- var="character",
- prediction="logical"),
+ dataLocs="Spatial",
+ predLocs="optionalLocs",
+ prediction="logical",
+ var="character"),
validity = validNeighbourhood, contains = list("Spatial"))
## ST neighbourhood
Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/R/spCopula.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -213,7 +213,7 @@
} else {
if(class(lowerCop) != "indepCopula") {
lowerParam <- calibPar(lowerCop, h)
- lowerCop at parameters <- lowerParam
+ lowerCop at parameters[length(lowerParam)] <- lowerParam
}
return(fun(pairs, lowerCop, ...))
}
@@ -608,23 +608,31 @@
## dropping a spatial tree, returning a conditional neighbourhood
dropSpTree <- function(neigh, spCop) {
- u1 <- NULL
- h1 <- NULL
+ u1 <- matrix(NA,nrow(neigh at data),ncol(neigh at data)-1)
+ h1 <- matrix(NA,nrow(neigh at distances),ncol(neigh at distances)-1)
+
for(i in 1:ncol(neigh at distances)) {
- u1 <- cbind(u1, dduCopula(as.matrix(neigh at data[,c(1,1+i)]), spCop,
- neigh at distances[,i]))
+ u1[,i] <- dduCopula(as.matrix(neigh at data[,c(1,1+i)]), spCop,
+ neigh at distances[,i])
if (i < ncol(neigh at distances)) {
- h1 <- cbind(h1, apply(neigh at index[,c(1,i+1)],1,
- function(x) spDists(neigh at locations[x,])[1,2]))
+ h1[,i] <- apply(neigh at index[,c(2,2+i)],1,
+ function(x) spDists(neigh at dataLocs[x,])[1,2])
}
}
varSplit <- strsplit(neigh at var,"|",fixed=TRUE)[[1]]
cond <- suppressWarnings(as.numeric(varSplit[length(varSplit)]))
- if(is.na(cond))
- cond <- paste(neigh at var,"|0",sep="")
- else
- cond <- paste(neigh at var,cond+1,sep="")
- return(neighbourhood(u1, h1, neigh at locations, neigh at dataLocs, neigh at index[,-1], neigh at prediction,
- cond))
+ if(is.na(cond)) {
+ var <- paste(neigh at var,"|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="")
+ 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))
}
\ No newline at end of file
Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/R/spVineCopula.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -103,7 +103,7 @@
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,spTree+i)],1,
- function(x) spDists(data at locations[x,])[1,2]))
+ function(x) spDists(data at dataLocs[x,])[1,2]))
}
}
u0 <- u1
@@ -142,11 +142,6 @@
# deriving all spatial tree distances
calcSpTreeDists <- function(neigh, n.trees) {
- if(!neigh at prediction)
- data <- neigh at locations
- else
- data <- neigh at dataLocs
-
condDists <- list(n.trees)
condDists[[1]] <- neigh at distances
if(n.trees==1)
@@ -154,8 +149,8 @@
for (spTree in 1:(n.trees-1)) {
h1 <- NULL
for(i in 1:(ncol(neigh at distances)-spTree)) {
- h1 <- cbind(h1,apply(neigh at index[,c(spTree,i+spTree),drop=F],1,
- function(x) spDists(data[x,])[1,2]))
+ 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]))
dimnames(h1) <- NULL
}
condDists[[spTree+1]] <- h1
@@ -201,7 +196,7 @@
predMean <- NULL
for(i in 1:nrow(predNeigh at data)) { # i <-1
cat("[Predicting location ",i,".]\n", sep="")
- condSecVine <- condSpVine(as.numeric(predNeigh at data[i,]),
+ condSecVine <- condSpVine(as.numeric(predNeigh at data[i,-1]),
lapply(dists,function(x) x[i,]), spVine)
condExp <- function(x) {
@@ -214,14 +209,14 @@
ePred$abs.error, " for location ",i,".")
predMean <- c(predMean, ePred$value)
}
- if ("data" %in% slotNames(predNeigh at locations)) {
- res <- predNeigh at locations
+ if ("data" %in% slotNames(predNeigh at predLocs)) {
+ res <- predNeigh at predLocs
res at data[["expect"]] <- predMean
return(res)
} else {
predMean <- data.frame(predMean)
colnames(predMean) <- "expect"
- return(addAttrToGeom(predNeigh at locations, predMean, match.ID=FALSE))
+ return(addAttrToGeom(predNeigh at predLocs, predMean, match.ID=FALSE))
}
}
@@ -232,7 +227,7 @@
predQuantile <- NULL
for(i in 1:nrow(predNeigh at data)) { # i <-1
cat("[Predicting location ",i,".]\n", sep="")
- condSecVine <- condSpVine(as.numeric(predNeigh at data[i,]),
+ condSecVine <- condSpVine(as.numeric(predNeigh at data[i,-1]),
lapply(dists,function(x) x[i,]), spVine)
xVals <- attr(condSecVine,"xVals")
@@ -253,14 +248,14 @@
predQuantile <- c(predQuantile, margin$q(xVals[lower]+xRes))
}
- if ("data" %in% slotNames(predNeigh at locations)) {
- res <- predNeigh at locations
+ if ("data" %in% slotNames(predNeigh at predLocs)) {
+ res <- predNeigh at 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 locations, predQuantile, match.ID=FALSE))
+ return(addAttrToGeom(predNeigh at predLocs, predQuantile, match.ID=FALSE))
}
}
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/R/spatialPreparation.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -7,18 +7,15 @@
## spatial neighbourhood constructor
####################################
-neighbourhood <- function(data, distances, sp, dataLocs=NULL, index,
+neighbourhood <- function(data, distances, index, dataLocs, predLocs=NULL,
prediction, var) {
sizeN <- ncol(distances)+1
data <- as.data.frame(data)
- colnames(data) <- paste(paste("N", rep((0+as.numeric(prediction)):(sizeN-1),
- each=length(var)), sep=""),
- rep(var,(sizeN-prediction)),sep=".")
if (anyDuplicated(rownames(data))>0)
rownames <- 1:length(rownames)
- new("neighbourhood", data=data, distances=distances, locations=sp,
- dataLocs=dataLocs, bbox=sp at bbox, proj4string=sp at proj4string, index=index,
- var=var, prediction=prediction)
+ new("neighbourhood", data=data, distances=distances, index=index,
+ dataLocs=dataLocs, predLocs=predLocs, prediction=prediction, var=var,
+ bbox=dataLocs at bbox, proj4string=dataLocs at proj4string)
}
## show
@@ -31,13 +28,13 @@
setMethod(show,signature("neighbourhood"),showNeighbourhood)
## names (from sp)
-setMethod(names, signature("neighbourhood"), namesNeighbourhood <- function(x) x at var)
+setMethod(names, signature("neighbourhood"), function(x) x at var)
## 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 locations,
+ spdf <- SpatialPointsDataFrame(coords=obj at dataLocs,
data=obj at data[,pattern,drop=FALSE],
proj4string=obj at proj4string, bbox=obj at bbox)
spplot(spdf, ...)
@@ -49,7 +46,7 @@
newSp <- x at locations[i,]
new("neighbourhood", data=x at data[i,,drop=F],
distances=x at distances[i,,drop=F],
- locations=newSp, dataLocs=x at dataLocs, bbox=newSp at bbox,
+ 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)
}
@@ -61,45 +58,44 @@
# returns an neighbourhood object
##################################
-getNeighbours <- function(spData, locations, var=names(spData)[1], size=5,
+getNeighbours <- function(dataLocs, predLocs, var=names(dataLocs)[1], size=5,
prediction=FALSE, min.dist=0.01) {
- stopifnot((!prediction && missing(locations)) || (prediction && !missing(locations)))
+ stopifnot((!prediction && missing(predLocs)) || (prediction && !missing(predLocs)))
stopifnot(min.dist>0 || prediction)
- if(missing(locations) && !prediction)
- locations=spData
+ if(missing(predLocs) && !prediction)
+ predLocs=dataLocs
- stopifnot(is(locations,"Spatial"))
+ stopifnot(is(predLocs,"Spatial"))
- nLocs <- length(locations)
-
- if(any(is.na(match(var,names(spData)))))
+ if(any(is.na(match(var,names(dataLocs)))))
stop("At least one of the variables is unkown or is not part of the data.")
- size <- min(size,length(spData)+prediction)
-
- allDists <- NULL
- allLocs <- NULL
- allData <- NULL
+ nLocs <- length(predLocs)
+ size <- min(size, length(dataLocs)+prediction)
- for(i in 1:length(locations)) { # i <- 1
- tempDists <- spDists(spData,locations[i,])
+ 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 <- rbind(allLocs, spLocs)
- allDists <- rbind(allDists, tempDists[spLocs])
+ spLocs <- order(tempDists)[1:(size - 1)]
- if(!prediction)
- spLocs <- c(i,spLocs)
- allData <- rbind(allData, as.vector(spData[spLocs, var, drop=F]@data[[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)
- dataLocs <- spData
- else
- dataLocs <- NULL
- return(neighbourhood(allData, allDists, locations, dataLocs,
- allLocs, prediction, var))
+
+ 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))
}
#############
@@ -148,6 +144,8 @@
if(any(is.na(boundaries))) {
boundaries <- ((1:nbins) * cutoff/nbins)
}
+
+ nbins <- length(boundaries)-1
lags <- calcSpLagInd(data, boundaries)
@@ -196,14 +194,16 @@
boundaries <- unique(c(0,boundaries))
}
- np <- numeric(0)
- moa <- numeric(0)
+ nbins <- length(boundaries)-1
+
+ np <- numeric(nbins)
+ moa <- numeric(nbins)
+ meanDists <- numeric(nbins)
lagData <- NULL
- meanDists <- numeric(0)
data <- as.matrix(data at data)
- for ( i in 1:nbins) {
+ for (i in 1:nbins) {
bools <- (dists <= boundaries[i+1] & dists > boundaries[i])
pairs <- NULL
@@ -212,14 +212,14 @@
}
lagData <- append(lagData, list(pairs))
- moa <- c(moa, corFun(pairs))
- meanDists <- c(meanDists, mean(dists[bools]))
- np <- c(np, sum(bools))
+ moa[i] <- corFun(pairs)
+ meanDists[i] <- mean(dists[bools])
+ np[i] <- sum(bools)
}
if(plot) {
plot(meanDists, moa, xlab="distance", ylab=paste("correlation [",cor.method,"]",sep=""),
- ylim=1.05*c(-abs(min(moa)),max(moa)), xlim=c(0,max(meanDists)))
+ ylim=1.05*c(-abs(min(moa, na.rm=T)),max(moa, na.rm=T)), xlim=c(0,max(meanDists,na.rm=T)))
abline(h=c(-min(moa),0,min(moa)),col="grey")
}
Modified: pkg/R/spatio-temporalPreparation.R
===================================================================
--- pkg/R/spatio-temporalPreparation.R 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/R/spatio-temporalPreparation.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -35,7 +35,7 @@
setMethod(show,signature("stNeighbourhood"),showStNeighbourhood)
-## calculate neighbourhood from SpatialPointsDataFrame
+## calculate neighbourhood from ST
# returns an neighbourhood object
##################################
@@ -78,13 +78,13 @@
}
stNeighData <- rbind(stNeighData, tmpData) # bind data row-wise
stDists[(i-1)*timeSteps+1:timeSteps,,1] <- matrix(rep(nghbrs at distances[i,],
- timeSteps*(spSize-1)),
- byrow=T, ncol=length(t.lags)*(spSize-1)) # store sp distances
+ 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,],
- timeSteps*(spSize-1)),
+ timeSteps*length(t.lags)),
byrow=T, ncol=length(t.lags)*(spSize-1))
stInd[(i-1)*timeSteps+1:timeSteps,,2] <- tmpInd
}
@@ -95,4 +95,112 @@
dataLocs <- NULL
return(stNeighbourhood(as.data.frame(stNeighData), stDists, stData, ST,
stInd, prediction, var))
-}
\ No newline at end of file
+}
+
+calcStNeighBins <- function(data, var="uniPM10", nbins=9, t.lags=-(0:2),
+ boundaries=NA, cutoff=NA, cor.method="fasttau") {
+# dists <- data at distances[,,1]
+#
+# corFun <- switch(cor.method,
+# fasttau=function(x) VineCopula:::fasttau(x[,1],x[,2]),
+# function(x) cor(x,method=cor.method)[1,2])
+#
+# if (any(is.na(boundaries)))
+# boundaries <- quantile(as.vector(dists), probs=c(1:nbins/nbins))
+# if(!is.na(cutoff)) {
+# boundaries <- boundaries[boundaries < cutoff]
+# boundaries <- unique(c(0,boundaries,cutoff))
+# } else {
+# boundaries <- unique(c(0,boundaries))
+# }
+#
+# lagData <- NULL
+# for(t.lag in t.lags) { # t.lag <- 0
+# tBool <- data at distances[,,2]==t.lag
+# tmpLagData <- NULL
+# for(i in 1:nbins) { # i <- 1
+# sBool <- (dists <= boundaries[i + 1] & dists > boundaries[i])
+# bool <- tBool & sBool
+# pairs <- NULL
+# for (col in 1:(dim(tBool)[2])) { # col <- 1
+# if(!any(bool[, col]))
+# next
+# sInd <- data at index[bool[, col], c(1, 1 + col),1]
+# tInd <- data at index[bool[, col], c(1, 1 + col),2]
+# p1 <- apply(cbind(sInd[,1], tInd[,1]),1,
+# function(x) data at locations[x[1], x[2],var])
+# p2 <- apply(cbind(sInd[,2], tInd[,2]),1,
+# function(x) data at locations[x[1], x[2],var])
+# pairs <- rbind(pairs, cbind(p1,p2))
+# }
+# tmpLagData <- append(tmpLagData,list(pairs))
+# }
+# lagData <- append(lagData,list(tmpLagData))
+#
+# }
+#
+# 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)) {
+# cors <- c(cors, VineCopula:::fasttau(binnedData[, 2 * i - 1], binnedData[, 2 * i]))
+# }
+# return(cors)
+# }
+# calcCor <- switch(cor.method, fasttau = calcTau, calcStats)
+# lagCor <- sapply(lagData, calcCor)
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+# np <- numeric(0)
+# moa <- numeric(0)
+# lagData <- NULL
+# meanDists <- numeric(0)
+#
+# data <- as.matrix(data at data)
+#
+# for ( i in 1:nbins) {
+# bools <- (dists <= boundaries[i+1] & dists > boundaries[i])
+#
+# pairs <- NULL
+# for(col in 1:(dim(bools)[2])) {
+# pairs <- rbind(pairs, data[bools[,col],c(1,1+col)])
+# }
+#
+# lagData <- append(lagData, list(pairs))
+# moa <- c(moa, corFun(pairs))
+# meanDists <- c(meanDists, mean(dists[bools]))
+# np <- c(np, sum(bools))
+# }
+#
+# if(plot) {
+# plot(meanDists, moa, xlab="distance", ylab=paste("correlation [",cor.method,"]",sep=""),
+# ylim=1.05*c(-abs(min(moa)),max(moa)), xlim=c(0,max(meanDists)))
+# abline(h=c(-min(moa),0,min(moa)),col="grey")
+# }
+#
+# res <- list(np=np, meanDists = meanDists, lagCor=moa, lagData=lagData)
+# attr(res,"cor.method") <- switch(cor.method, fasttau="kendall", cor.method)
+# return(res)
+}
+
+setMethod(calcBins, signature="stNeighbourhood", calcStNeighBins)
\ No newline at end of file
Modified: pkg/R/stCopula.R
===================================================================
--- pkg/R/stCopula.R 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/R/stCopula.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -68,7 +68,7 @@
for(t in tDist) {
tmpInd <- h[,2]==t
tmpCop <- copula at spCopList[[match(t, copula at t.lags)]]
- res[tmpInd] <- pSpCopula(u[tmpInd,], tmpCop, h[tmpInd,1])
+ res[tmpInd] <- pSpCopula(u[tmpInd,,drop=F], tmpCop, h[tmpInd,1])
}
}
res
@@ -98,7 +98,7 @@
for(t in tDist) {
tmpInd <- h[,2]==t
tmpCop <- copula at spCopList[[match(t, copula at t.lags)]]
- res[tmpInd] <- dSpCopula(u[tmpInd,], tmpCop, log, h[tmpInd,1])
+ res[tmpInd] <- dSpCopula(u[tmpInd,,drop=F], tmpCop, log, h[tmpInd,1])
}
}
res
@@ -131,7 +131,7 @@
for(t in tDist) {
tmpInd <- h[,2]==t
tmpCop <- copula at spCopList[[match(t, copula at t.lags)]]
- res[tmpInd] <- dduSpCopula(u[tmpInd,], tmpCop, h[tmpInd,1])
+ res[tmpInd] <- dduSpCopula(u[tmpInd,,drop=F], tmpCop, h[tmpInd,1])
}
}
res
@@ -162,7 +162,7 @@
for(t in tDist) {
tmpInd <- h[,2]==t
tmpCop <- copula at spCopList[[match(t, copula at t.lags)]]
- res[tmpInd] <- ddvSpCopula(u[tmpInd,], tmpCop, h[tmpInd,1])
+ res[tmpInd] <- ddvSpCopula(u[tmpInd,,drop=F], tmpCop, h[tmpInd,1])
}
}
res
Modified: pkg/R/stVineCopula.R
===================================================================
--- pkg/R/stVineCopula.R 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/R/stVineCopula.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -81,15 +81,18 @@
stopifnot(class(data)=="stNeighbourhood")
stopifnot(copula at dimension == ncol(data at data))
- u0 <- as.matrix(data at data) # previous level's (contitional) data
+ u0 <- as.matrix(data at data) # previous level's (conditional) data
h0 <- data at distances # previous level's distances
l0 <- rep(0,nrow(u0)) # spatial density
u1 <- NULL # current level of conditional data
+ cat("[Margin ")
for(i in 1:dim(h0)[2]) { # i <- 1
l0 <- l0 + dCopula(u0[,c(1,i+1)], copula at stCop, h=h0[,i,], log=T)
+ cat(i,", ")
u1 <- cbind(u1, dduCopula(u0[,c(1,i+1)], copula at stCop, h=h0[,i,]))
}
u0 <- u1
+ cat(".]\n")
cat("[Estimating a",ncol(u0),"dimensional copula at the top.]\n")
vineCopFit <- fitCopula(copula at topCop, u0, method, estimate.variance)
@@ -199,4 +202,37 @@
switch(method,
quantile=stCopPredict.quantile(predNeigh, stVine, margin, p),
expectation=stCopPredict.expectation(predNeigh, stVine, margin, ...))
+}
+
+dropStTree <- function(neigh, stCop) {
+ stopifnot(class(neigh)=="stNeighbourhood")
+
+ u0 <- as.matrix(neigh at data) # previous level's (conditional) data
+ h0 <- neigh at distances # previous level's distances
+ u1 <- NULL # current level of conditional data
+ h1s <- NULL # upcoming distances
+ h1t <- NULL # upcoming distances
+ cat("[Margin ")
+ for(i in 1:dim(h0)[2]) { # i <- 1
+ cat(i,", ")
+ u1 <- cbind(u1, dduCopula(u0[,c(1,i+1)], stCop, h=h0[,i,]))
+ if (i < ncol(neigh at distances)) {
+ h1s <- cbind(h1s, apply(neigh at index[, c(1, i + 1),1], 1,
+ function(x) spDists(neigh at locations@sp[x, ])[1, 2]))
+ h1t <- cbind(h1t, apply(neigh at index[, c(1, i + 1),2], 1,
+ function(x) diff(x)))
+ }
+ }
+ h1 <- array(dim=c(dim(h1s),2))
+ h1[,,1] <- h1s
+ h1[,,2] <- h1t
+
+ varSplit <- strsplit(neigh at var, "|", fixed = TRUE)[[1]]
+ cond <- suppressWarnings(as.numeric(varSplit[length(varSplit)]))
+ if (is.na(cond))
+ cond <- paste(neigh at var, "|0", sep = "")
+ else cond <- paste(neigh at var, cond + 1, sep = "")
+ return(stNeighbourhood(data=u1, distances=h1, STxDF=neigh at locations,
+ ST=neigh at dataLocs, index=neigh at index[, -1,],
+ prediction=neigh at prediction, var=cond))
}
\ No newline at end of file
Modified: pkg/demo/00Index
===================================================================
--- pkg/demo/00Index 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/demo/00Index 2013-07-10 07:07:34 UTC (rev 99)
@@ -1,2 +1,3 @@
MRP The MRP demo gives insight in the code used in the paper: Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation, by Gräler et al. (2013), HESS-17-1281-2013.
-spCopula A demo illustrating the estiamtion of a spatial copula for a SpatialPointsDataFrame.
+spCopula A demo illustrating the estiamtion of a single spatial tree vine copula for a SpatialPointsDataFrame.
+pureSpVineCopula A demo illustrating the estiamtion of a pure spatial vine copula for a SpatialPointsDataFrame.
Added: pkg/demo/pureSpVineCopula.R
===================================================================
--- pkg/demo/pureSpVineCopula.R (rev 0)
+++ pkg/demo/pureSpVineCopula.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -0,0 +1,198 @@
+## librarys ##
+library(spcopula)
+# library(evd)
+
+## meuse - spatial poionts data.frame ##
+data(meuse)
+coordinates(meuse) = ~x+y
+
+spplot(meuse,"zinc", col.regions=bpy.colors(5))
+
+## margins ##
+hist(meuse[["zinc"]],freq=F,n=30,ylim=c(0,0.0035),
+ main="Histogram of zinc", xlab="zinc concentration")
+
+ #gevEsti <- fgev(meuse[["zinc"]])$estimate
+meanLog <- mean(log(meuse[["zinc"]]))
+sdLog <- sd(log(meuse[["zinc"]]))
+# curve(dgev(x,gevEsti[1], gevEsti[2], gevEsti[3]),add=T,col="red")
+curve(dlnorm(x,meanLog,sdLog),add=T,col="green")
+
+pMar <- function(q) plnorm(q, meanLog, sdLog)
+qMar <- function(p) qlnorm(p, meanLog, sdLog)
+dMar <- function(x) dlnorm(x, meanLog, sdLog)
+
+# pMar <- function(q) pgev(q, gevEsti[1], gevEsti[2], gevEsti[3])
+# qMar <- function(p) qgev(p, gevEsti[1], gevEsti[2], gevEsti[3])
+# dMar <- function(x) dgev(x, gevEsti[1], gevEsti[2], gevEsti[3])
+
+meuse$rtZinc <- rank(meuse$zinc)/(length(meuse)+1)
+
+## lag classes ##
+bins <- calcBins(meuse,var="rtZinc", nbins=10, cutoff=800)
+
+## calculate parameters for Kendall's tau function ##
+calcKTau <- fitCorFun(bins, degree=2)
+curve(calcKTau,0, 1000, col="purple",add=TRUE)
+
+families <- list(normalCopula(0), tCopula(0,df=2.15), claytonCopula(0),
+ gumbelCopula(1), frankCopula(1), joeBiCopula(1.5),
+ surClaytonCopula(1), surGumbelCopula(1), surJoeBiCopula(1.5))
+
+## find best fitting copula per lag class
+loglikTau <- loglikByCopulasLags(bins, families, calcKTau)
+bestFitTau <- apply(apply(loglikTau$loglik, 1, rank, na.last=T), 2,
+ function(x) which(x==9))
+colnames(loglikTau$loglik)[bestFitTau]
+
+## set up a first bivariate spatial Copula
+###########################################
+spCop <- spCopula(c(families[bestFitTau[1:7]],indepCopula()),
+ distances=bins$meanDists[1:8],
+ spDepFun=calcKTau, unit="m")
+
+## estimation neighbourhood for the pure spatial vine copula
+#############################################################
+vineDim <- 20L
+meuseNeigh1 <- getNeighbours(dataLocs=meuse,var="rtZinc",size=vineDim)
+
+## second spatial tree
+#######################
+meuseNeigh2 <- dropSpTree(meuseNeigh1, spCop)
+bins2 <- calcBins(meuseNeigh2, boundaries=c(0,2:15)*50, plot=F)
+points(bins2$meanDists, bins2$lagCor, pch=2)
+bins2$np
+calcKTau2 <- fitCorFun(bins2, degree=3,cutoff=500)
+curve(calcKTau2,0, 800, col="green",add=TRUE)
+curve(calcKTau, 0, 800, col="purple",add=TRUE)
+
+loglikTau2 <- loglikByCopulasLags(bins2, families) #, calcKTau2)
+bestFitTau2 <- apply(apply(loglikTau2$loglik, 1, rank, na.last=T), 2,
+ function(x) which(x==9))
+
+colnames(loglikTau2$loglik)[bestFitTau2]
+
+## set up the second bivariate spatial Copula
+##############################################
+spCop2 <- spCopula(c(families[bestFitTau2[1:5]],indepCopula()),
+ distances=bins2$meanDists[1:6],
+ spDepFun=calcKTau2, unit="m")
+
+## third spatial tree
+######################
+meuseNeigh3 <- dropSpTree(meuseNeigh2, spCop2)
+bins3 <- calcBins(meuseNeigh3)
+calcKTau3 <- fitCorFun(bins3, degree=1,cutoff=400)
+curve(calcKTau3,0, 500, col="purple",add=TRUE)
+
+loglikTau3 <- loglikByCopulasLags(bins3, families, calcKTau3)
+bestFitTau3 <- apply(apply(loglikTau3$loglik, 1, rank, na.last=T), 2,
+ function(x) which(x==9))
+colnames(loglikTau3$loglik)[bestFitTau3]
+
+## set up the third bivariate spatial Copula
+#############################################
+spCop3 <- spCopula(c(families[bestFitTau3[1:6]],indepCopula()),
+ distances=bins3$meanDists[1:7],
+ spDepFun=calcKTau3, unit="m")
+
+## fourth spatial tree
+#######################
+meuseNeigh4 <- dropSpTree(meuseNeigh3, spCop3)
+bins4 <- calcBins(meuseNeigh4)
+calcKTau4 <- fitCorFun(bins4, degree=1,cutoff=400)
+curve(calcKTau4,0, 500, col="purple",add=TRUE)
+
+loglikTau4 <- loglikByCopulasLags(bins4, families, calcKTau4)
+bestFitTau4 <- apply(apply(loglikTau4$loglik, 1, rank, na.last=T), 2,
+ function(x) which(x==9))
+colnames(loglikTau4$loglik)[bestFitTau4]
+
+## set up the third bivariate spatial Copula
+#############################################
+spCop4 <- spCopula(c(families[bestFitTau4[1:4]],indepCopula()),
+ distances=bins4$meanDists[1:5],
+ spDepFun=calcKTau4, unit="m")
+
+## pure spatial vine
+#####################
+meuseSpVine <- spVineCopula(list(spCop, spCop2, spCop3, spCop4))
+
+## neighbourhood for cross-validation using a 5 dim pure spatial vine copula
+#############################################################################
+vineDim <- 5L
+
+meuse$lnZinc <- pMar(meuse$zinc)
+meuseNeigh <- getNeighbours(dataLocs=meuse, predLocs=meuse, prediction=T,
+ min.dist=10, var="lnZinc", size=vineDim)
+
+# meuse$evZinc <- pMar(meuse$zinc)
+# meuseNeigh <- getNeighbours(dataLocs=meuse, predLocs=meuse, prediction=T,
+# min.dist=10, var="evZinc", size=vineDim)
+
+meuseNeigh <- getNeighbours(dataLocs=meuse, predLocs=meuse, prediction=T,
+ min.dist=10, var="rtZinc", size=vineDim)
+
+## leave-one-out x-validation
+##############################
+
+time <- proc.time() # ~240 s
+predMedian <- spCopPredict(meuseNeigh, meuseSpVine, list(q=qMar), "quantile")
+predMean <- spCopPredict(meuseNeigh, meuseSpVine, list(q=qMar), "expectation")
+proc.time() - time
+
+c(mean(abs(predMean$expect-meuse$zinc)),
+ mean(predMean$expect-meuse$zinc),
+ sqrt(mean((predMean$expect-meuse$zinc)^2)))
+
+c(mean(abs(predMedian$quantile.0.5-meuse$zinc)),
+ mean(predMedian$quantile.0.5-meuse$zinc),
+ sqrt(mean((predMedian$quantile.0.5-meuse$zinc)^2)))
+
+plot(predMean$expect, meuse$zinc,asp=1)
+abline(0,1)
+
+plot(predMedian$quantile.0.5, meuse$zinc,asp=1)
+abline(0,1)
+
+## copula, evd:
+# Median:
+# MAE: 164
+# ME: -82
+# RMSE: 265
+# Mean:
+# MAE: 278
+# ME: 196
+# RMSE: 507
+
+## copula, lnorm, rt:
+# Median:
+# MAE: 167
+# ME: -61
+# RMSE: 263
+# Mean:
+# MAE: 171
+# ME: 0
+# RMSE: 247
+
+## copula, lnorm, pMar:
+# Median:
+# MAE: 163
+# ME: -60
+# RMSE: 263
+# Mean:
+# MAE: 167
+# ME: 2
+# RMSE: 248
+
+
+## kriging results:
+# same neighbourhood size 5L:
+# MAE: 158.61
+# BIAS: -4.24
+# RMSE: 239.85
+#
+# global kriging:
+# MAE: 148.85
+# BIAS: -3.05
+# RMSE: 226.15
\ No newline at end of file
Modified: pkg/demo/spCopula.R
===================================================================
--- pkg/demo/spCopula.R 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/demo/spCopula.R 2013-07-10 07:07:34 UTC (rev 99)
@@ -89,7 +89,7 @@
##
# leave-one-out x-validation
-time <- proc.time() # ~60 s
+time <- proc.time() # ~100 s
predMedian <- NULL
predMean <- NULL
for(loc in 1:nrow(meuseNeigh at data)) { # loc <- 145
Modified: pkg/man/calcBins.Rd
===================================================================
--- pkg/man/calcBins.Rd 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/man/calcBins.Rd 2013-07-10 07:07:34 UTC (rev 99)
@@ -5,6 +5,7 @@
\alias{calcBins,Spatial-method}
\alias{calcBins,STFDF-method}
\alias{calcBins,neighbourhood-method}
+\alias{calcBins,stNeighbourhood-method}
\title{
A function calculating the spatial/spatio-temporal bins
Modified: pkg/man/getNeighbours.Rd
===================================================================
--- pkg/man/getNeighbours.Rd 2013-06-19 08:08:51 UTC (rev 98)
+++ pkg/man/getNeighbours.Rd 2013-07-10 07:07:34 UTC (rev 99)
@@ -8,13 +8,13 @@
This function calculates a local neighbourhood to be used for fitting of spatial/spatio-temporal vine copulas and for prediction using spatial/spatio-temporal vine copulas.
}
\usage{
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/spcopula -r 99
More information about the spcopula-commits
mailing list