[spcopula-commits] r110 - in pkg: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 24 13:37:57 CEST 2013
Author: ben_graeler
Date: 2013-10-24 13:37:56 +0200 (Thu, 24 Oct 2013)
New Revision: 110
Added:
pkg/demo/stVineCopFit.R
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/asCopula.R
pkg/R/cqsCopula.R
pkg/R/spCopula.R
pkg/R/spatialPreparation.R
pkg/R/spatio-temporalPreparation.R
pkg/R/stVineCopula.R
pkg/demo/00Index
pkg/man/calcBins.Rd
pkg/man/getNeighbours.experimental.Rd
Log:
- added spatio-temporal demo stVineCopFit.R
- applied a few chages from the spatial case to the spatio-temporal one
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/DESCRIPTION 2013-10-24 11:37:56 UTC (rev 110)
@@ -2,7 +2,7 @@
Type: Package
Title: copula driven spatial analysis
Version: 0.1-1
-Date: 2013-09-18
+Date: 2013-10-24
Authors at R: c(person("Benedikt", "Graeler", role = c("aut", "cre"), email =
"ben.graeler at uni-muenster.de"), person("Marius", "Appel",
role = "ctb"))
@@ -10,8 +10,8 @@
Description: This package provides a framework to analyse via copulas spatial and spatio-temporal data provided in the format of the spacetime package. Additionally, support for calculating different multivariate return periods is implemented.
License: GPL-2
LazyLoad: yes
-Depends: copula (>= 0.999-7), VineCopula (>= 1.1-2), methods, R (>= 2.15.0)
-Imports: sp, spacetime (>= 1.0-2)
+Depends: copula (>= 0.999-7), VineCopula (>= 1.1-2), R (>= 2.15.0)
+Imports: sp, spacetime (>= 1.0-2), methods
URL: http://r-forge.r-project.org/projects/spcopula/
Collate:
Classes.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/NAMESPACE 2013-10-24 11:37:56 UTC (rev 110)
@@ -1,5 +1,6 @@
import(copula, VineCopula)
import(sp, spacetime)
+import(methods)
# constructor
export(asCopula, cqsCopula)
Modified: pkg/R/asCopula.R
===================================================================
--- pkg/R/asCopula.R 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/R/asCopula.R 2013-10-24 11:37:56 UTC (rev 110)
@@ -15,7 +15,7 @@
## density ##
-dASC2 <- function (u, copula) {
+dASC2 <- function (u, copula, log=FALSE) {
a <- copula at parameters[1]
b <- copula at parameters[2]
@@ -24,12 +24,15 @@
u1 <- u[, 1]
u2 <- u[, 2]
- return(pmax(a * u2 * (((12 - 9 * u1) * u1 - 3) * u2 + u1 * (6 * u1 - 8) + 2) + b * (u2 * ((u1 * (9 * u1 - 12) + 3) * u2 + (12 - 6 * u1) * u1 - 4) - 2 * u1 + 1) + 1,0))
+ if(log)
+ return(log(pmax(a * u2 * (((12 - 9 * u1) * u1 - 3) * u2 + u1 * (6 * u1 - 8) + 2) + b * (u2 * ((u1 * (9 * u1 - 12) + 3) * u2 + (12 - 6 * u1) * u1 - 4) - 2 * u1 + 1) + 1,0)))
+ else
+ return(pmax(a * u2 * (((12 - 9 * u1) * u1 - 3) * u2 + u1 * (6 * u1 - 8) + 2) + b * (u2 * ((u1 * (9 * u1 - 12) + 3) * u2 + (12 - 6 * u1) * u1 - 4) - 2 * u1 + 1) + 1,0))
}
setMethod("dCopula", signature("numeric","asCopula"),
- function(u, copula, ...) {
- dASC2(matrix(u,ncol=copula at dimension),copula)
+ function(u, copula, log) {
+ dASC2(matrix(u,ncol=copula at dimension),copula, log)
})
setMethod("dCopula", signature("matrix","asCopula"), dASC2)
Modified: pkg/R/cqsCopula.R
===================================================================
--- pkg/R/cqsCopula.R 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/R/cqsCopula.R 2013-10-24 11:37:56 UTC (rev 110)
@@ -1,5 +1,5 @@
## make fitCopula generic
-setGeneric("fitCopula")
+setGeneric("fitCopula",fitCopula)
######################################################
## ##
@@ -233,7 +233,7 @@
# method
# one of kendall or spearman according to the calculation of moa
-fitCQSec.itau <- function(copula, data, estimate.variance, tau=NULL) {
+fitCQSec.itau <- function(copula, data, estimate.variance=FALSE, tau=NULL) {
if(is.null(tau))
tau <- VineCopula:::fasttau(data[,1],data[,2])
if(copula at fixed=="a")
@@ -253,7 +253,7 @@
}
-fitCQSec.irho <- function(copula, data, estimate.variance, rho=NULL){
+fitCQSec.irho <- function(copula, data, estimate.variance=FALSE, rho=NULL){
if(is.null(rho))
rho <- cor(data,method="spearman")[1,2]
if(copula at fixed=="a")
Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/R/spCopula.R 2013-10-24 11:37:56 UTC (rev 110)
@@ -492,19 +492,36 @@
tmpCop <- list()
for(i in 1:length(bins$meanDists)) {
if(class(cop)!="indepCopula") {
- param <- moa(cop, bins$meanDists[i])
- cop at parameters[1:length(param)] <- param
+ if(class(cop) == "asCopula") {
+ cop <- switch(calcCor(NULL),
+ kendall=fitASC2.itau(cop, bins$lagData[[i]],
+ tau=calcCor(bins$meanDists[i]))@copula,
+ spearman=fitASC2.irho(cop, bins$lagData[[i]],
+ rho=calcCor(bins$meanDists[i]))@copula,
+ stop(paste(calcCor(NULL), "is not yet supported.")))
+ } else {
+ if(class(cop) == "cqsCopula") {
+ cop <- switch(calcCor(NULL),
+ kendall=fitCQSec.itau(cop, bins$lagData[[i]],
+ tau=calcCor(bins$meanDists[i]))@copula,
+ spearman=fitCQSec.irho(cop, bins$lagData[[i]],
+ rho=calcCor(bins$meanDists[i]))@copula,
+ stop(paste(calcCor(NULL), "is not yet supported.")))
+ } else {
+ param <- moa(cop, bins$meanDists[i])
+ cop at parameters[1:length(param)] <- param
+ }
+ }
}
- tmploglik <- c(tmploglik, sum(dCopula(bins$lagData[[i]],cop, log=T)))
+ tmploglik <- c(tmploglik, sum(dCopula(bins$lagData[[i]], cop, log=T)))
tmpCop <- append(tmpCop, cop)
}
loglik <- cbind(loglik, tmploglik)
- copulas <- append(copulas,tmpCop)
+ copulas[[class(cop)]] <- tmpCop
}
colnames(loglik) <- sapply(families, function(x) class(x)[1])
- names(copulas) <- colnames(loglik)
return(list(loglik=loglik, copulas=copulas))
}
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/R/spatialPreparation.R 2013-10-24 11:37:56 UTC (rev 110)
@@ -182,15 +182,15 @@
# the generic calcBins, calculates bins for spatial and spatio-temporal data
setGeneric("calcBins", function(data, var, nbins=15, boundaries=NA, cutoff=NA,
- cor.method="kendall", plot=TRUE, ...) {
+ ..., cor.method="fasttau", plot=TRUE) {
standardGeneric("calcBins")
})
## calculating the spatial bins
################################
-calcSpBins <- function(data, var=names(data), nbins=15, boundaries=NA,
- cutoff=NA, cor.method="kendall", plot=TRUE) {
+calcSpBins <- function(data, var, nbins=15, boundaries=NA, cutoff=NA, ...,
+ cor.method="fasttau", plot=TRUE) {
if(is.na(cutoff)) {
cutoff <- spDists(coordinates(t(data at bbox)))[1,2]/3
@@ -289,9 +289,8 @@
# 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=10, t.lags=-(0:2), cor.method="fasttau",
- plot=FALSE) {
-
+ 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))
@@ -360,4 +359,4 @@
return(res)
}
-setMethod(calcBins, signature("STFDF"), calcStBins)
+setMethod(calcBins, signature(data="STFDF"), calcStBins)
Modified: pkg/R/spatio-temporalPreparation.R
===================================================================
--- pkg/R/spatio-temporalPreparation.R 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/R/spatio-temporalPreparation.R 2013-10-24 11:37:56 UTC (rev 110)
@@ -39,6 +39,7 @@
# 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) {
stopifnot((!prediction && missing(ST)) || (prediction && !missing(ST)))
@@ -67,11 +68,11 @@
stInd <- array(NA,c(nLocs,(spSize-1)*length(t.lags),2))
for(i in 1:nrow(nghbrs at index)){ # i <- 1
tmpInst <- sample((1-timeSpan):length(stData at time), timeSteps) # draw random time steps for each neighbourhood
- tmpData <- matrix(stData[c(i, nghbrs at index[i,]), tmpInst, var]@data[[1]],
+ tmpData <- matrix(stData[nghbrs at index[i,], tmpInst, var]@data[[1]],
ncol=spSize, byrow=T) # retrieve the top level data
tmpInd <- matrix(rep(tmpInst,spSize-1),ncol=spSize-1)
for(t in t.lags[-1]) {
- tmpData <- cbind(tmpData, matrix(stData[nghbrs at index[i,],
+ tmpData <- cbind(tmpData, matrix(stData[nghbrs at index[i,][-1],
tmpInst+t,var]@data[[1]],
ncol=spSize-1, byrow=T))
tmpInd <- cbind(tmpInd, matrix(rep(tmpInst+t,spSize-1),ncol=spSize-1))
@@ -83,7 +84,7 @@
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,],
+ 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
Modified: pkg/R/stVineCopula.R
===================================================================
--- pkg/R/stVineCopula.R 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/R/stVineCopula.R 2013-10-24 11:37:56 UTC (rev 110)
@@ -88,11 +88,11 @@
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,", ")
+ cat(i,", ", sep="")
u1 <- cbind(u1, dduCopula(u0[,c(1,i+1)], copula at stCop, h=h0[,i,]))
}
u0 <- u1
- cat(".]\n")
+ cat("]\n")
cat("[Estimating a",ncol(u0),"dimensional copula at the top.]\n")
vineCopFit <- fitCopula(copula at topCop, u0, method, estimate.variance)
Modified: pkg/demo/00Index
===================================================================
--- pkg/demo/00Index 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/demo/00Index 2013-10-24 11:37:56 UTC (rev 110)
@@ -1,3 +1,4 @@
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 single spatial tree vine copula for a SpatialPointsDataFrame.
pureSpVineCopula A demo illustrating the estiamtion of a pure spatial vine copula for a SpatialPointsDataFrame.
+stVineCopFit A demo corresponding to the vignette estimating a spatio-temporal vine copula.
Added: pkg/demo/stVineCopFit.R
===================================================================
--- pkg/demo/stVineCopFit.R (rev 0)
+++ pkg/demo/stVineCopFit.R 2013-10-24 11:37:56 UTC (rev 110)
@@ -0,0 +1,67 @@
+## stVineCopula estimation following the vignette, but using a randomly
+## selected smaller subset of the original data to reduce calculation demands.
+## Thus, results are likely to differ (a little) from the original study.
+
+library(spcopula)
+load(url("http://ifgi.uni-muenster.de/~b_grae02/publications/EU_RB_2005.RData"))
+
+## spatio-temporal copula
+# binning, using only 90 out of all temporal instances
+stBins <- calcBins(EU_RB_2005, "rtPM10", nbins=40, t.lags=-(0:2),
+ instances=10)
+
+calcKTau <- fitCorFun(stBins,c(3,3,3))
+
+families <- c(normalCopula(0), tCopula(0), claytonCopula(0),
+ frankCopula(1), gumbelCopula(1), joeBiCopula(),
+ cqsCopula(), asCopula())
+
+loglikTau <- list()
+for(j in 1:3) { # j <-1
+ tmpBins <- list(meanDists=stBins$meanDists,
+ lagData=lapply(stBins$lagData, function(x) x[,c(2*j-1,2*j)]))
+ loglikTau[[j]] <- loglikByCopulasLags(tmpBins, families, calcKTau[[j]])
+}
+
+bestFitTau <- list()
+bestFitTau[[1]] <- apply(apply(loglikTau[[1]]$loglik[,-8], 1, rank),
+ 2, which.max)
+bestFitTau[[2]] <- apply(apply(loglikTau[[2]]$loglik, 1, rank),
+ 2, which.max)
+bestFitTau[[3]] <- apply(apply(loglikTau[[3]]$loglik, 1, rank),
+ 2, which.max)
+
+# gather the fitted copulas and representative distances
+listCops <- NULL
+
+for(t.level in 1:3) {
+ listCops[[t.level]] <- sapply(1:35,
+ function(i) {
+ tmpTau <- calcKTau[[t.level]](stBins$meanDists[i])
+ if (tmpTau < 0.05) {
+ return(indepCopula(dim=2))
+ } else {
+ return(loglikTau[[t.level]]$copulas[[bestFitTau[[t.level]][i]]][[i]])
+ }})
+}
+
+listDists <- NULL
+listDists[[1]] <- stBins$meanDists[1:35]
+listDists[[2]] <- stBins$meanDists[1:35]
+listDists[[3]] <- stBins$meanDists[1:35]
+
+# build the spatio-temporal copula
+stConvCop <- stCopula(components=listCops, distances=listDists,
+ t.lags=c(0,-1,-2))
+
+# get the neighbours
+stNeigh <- getStNeighbours(EU_RB_2005, var="rtPM10", spSize=4,
+ t.lags=-(0:2), timeSteps=10, min.dist=10)
+
+
+stVineFit <- fitCopula(stVineCopula(stConvCop,vineCopula(9L)), stNeigh,
+ method="indeptest")
+stVine <- stVineFit at copula
+
+# retrieve the log-likelihood
+stVineFit at loglik
\ No newline at end of file
Modified: pkg/man/calcBins.Rd
===================================================================
--- pkg/man/calcBins.Rd 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/man/calcBins.Rd 2013-10-24 11:37:56 UTC (rev 110)
@@ -14,8 +14,8 @@
The (spatio-temporal) space is subdivided into pairs of observations that belong to certain spatial/spatio-temporal distance classes. For each distance class, the mean separating distance of all pairs involved is calculated alongside a correlation measure. The spatial/spatio-temporal correlogram is plotted by default.
}
\usage{
-calcBins(data, var, nbins = 15, boundaries = NA, cutoff = NA, cor.method="kendall",
- plot=TRUE, ...)
+calcBins(data, var, nbins = 15, boundaries = NA, cutoff = NA,
+ ..., cor.method="fasttau", plot=TRUE)
}
\arguments{
@@ -34,16 +34,16 @@
\item{cutoff}{
The largest spatial distance to be investigated, if set to \code{NA}, one third of the bounding box's diagonal will be used.
}
- \item{cor.method}{
- defining the correlation method that should be used. Possible choices are \code{kendall} (default), \code{pearson}, \code{spearman} and \code{fasttau} that re-uses a very fast implementation of Kendall's tau from \code{\link{VineCopula-package}}.
- }
- \item{plot}{Whether the correlogram should be plotted.}
- \item{\dots}{Additional arguments for the spatio-temporal case:
- \describe{
+\item{\dots}{Additional arguments for the spatio-temporal case:
+\describe{
\item{instances}{To reduce the data size or circumvent unwanted autocorrelation effects, one might provide a number of randomly selected time instances from the spatio-temporal \code{data.frame}. If this parameter is set to \code{NA}, the complete time series will be used, if different from a single number, \code{instances} will be passed on as to index time.}
\item{t.lags}{a vector indicating the time lags to be investigated}
}
}
+ \item{cor.method}{
+ defining the correlation method that should be used. Possible choices are \code{kendall}, \code{pearson}, \code{spearman} and \code{fasttau} (default) that re-uses a very fast implementation of Kendall's tau from \code{\link{VineCopula-package}}.
+ }
+ \item{plot}{Whether the correlogram should be plotted.}
}
\value{A list containing the following entries
\item{meanDists}{A vector holding the mean separating distance for each lag.}
Modified: pkg/man/getNeighbours.experimental.Rd
===================================================================
--- pkg/man/getNeighbours.experimental.Rd 2013-09-18 13:13:11 UTC (rev 109)
+++ pkg/man/getNeighbours.experimental.Rd 2013-10-24 11:37:56 UTC (rev 110)
@@ -8,7 +8,8 @@
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{
-getNeighbours.experimental(dataLocs, predLocs, var = names(dataLocs)[1], size = 5, prediction=FALSE, min.dist = 0.01)
+getNeighbours.experimental(dataLocs, predLocs, var = names(dataLocs)[1], size = 5,
+prediction=FALSE, min.dist = 0.01)
}
\arguments{
\item{dataLocs}{
More information about the spcopula-commits
mailing list