[spcopula-commits] r97 - / pkg pkg/R pkg/man www
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 24 17:08:19 CEST 2013
Author: ben_graeler
Date: 2013-05-24 17:08:18 +0200 (Fri, 24 May 2013)
New Revision: 97
Added:
pkg/R/spatio-temporalPreparation.R
pkg/R/stVineCopula.R
pkg/man/condStVine.Rd
pkg/man/getStNeighbours.Rd
pkg/man/stCopPredict.Rd
pkg/man/stCopula.Rd
pkg/man/stNeighbourhood-class.Rd
pkg/man/stNeighbourhood.Rd
pkg/man/stVineCopula-class.Rd
pkg/man/stVineCopula.Rd
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/Classes.R
pkg/R/asCopula.R
pkg/R/cqsCopula.R
pkg/R/spCopula.R
pkg/R/spatialPreparation.R
pkg/R/stCopula.R
pkg/man/condSpVine.Rd
pkg/man/cqsCopula.Rd
pkg/man/neighbourhood-class.Rd
pkg/man/neighbourhood.Rd
pkg/man/spVineCopula.Rd
pkg/man/stCopula-class.Rd
spcopula_0.1-1.tar.gz
spcopula_0.1-1.zip
www/index.php
Log:
- added full spatio-temporal support (still needs tesating!)
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-05-21 08:26:21 UTC (rev 96)
+++ pkg/DESCRIPTION 2013-05-24 15:08:18 UTC (rev 97)
@@ -2,7 +2,7 @@
Type: Package
Title: copula driven spatial analysis
Version: 0.1-1
-Date: 2013-05-21
+Date: 2013-05-24
Author: Benedikt Graeler
Maintainer: Benedikt Graeler <ben.graeler at uni-muenster.de>
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.
@@ -20,6 +20,7 @@
spCopula.R
stCopula.R
spatialPreparation.R
+ spatio-temporalPreparation.R
wrappingCFunctions.R
linkingVineCopula.R
BB1copula.R
@@ -30,5 +31,6 @@
ClaytonGumbelCopula.R
vineCopulas.R
spVineCopula.R
+ stVineCopula.R
utilities.R
returnPeriods.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-05-21 08:26:21 UTC (rev 96)
+++ pkg/NAMESPACE 2013-05-24 15:08:18 UTC (rev 97)
@@ -9,8 +9,9 @@
export(joeBiCopula, surJoeBiCopula, r90JoeBiCopula, r270JoeBiCopula)
export(surClaytonCopula, r90ClaytonCopula, r270ClaytonCopula)
export(surGumbelCopula, r90GumbelCopula, r270GumbelCopula)
-export(vineCopula, spVineCopula)
-export(neighbourhood)
+export(spCopula, stCopula)
+export(vineCopula, spVineCopula, stVineCopula)
+export(neighbourhood, stNeighbourhood)
export(empiricalCopula, genEmpCop)
# general functions
@@ -20,18 +21,18 @@
export(invdduCopula, invddvCopula)
export(qCopula_u)
export(condSpVine,spCopPredict)
+export(condStVine,stCopPredict)
# tweaks
# export(setSizeLim)
# spatial
-export(getNeighbours)
+export(getNeighbours,getStNeighbours)
export(calcBins)
export(calcSpTreeDists, dropSpTree)
# fitting
export(fitCorFun, loglikByCopulasLags, fitSpCopula, composeSpCopula)
-export(spCopula)
# MRP functions
export(genEmpKenFun, genInvKenFun)
@@ -39,8 +40,8 @@
export(criticalPair, criticalTriple)
## classes
-exportClasses(asCopula, cqsCopula, neighbourhood, empiricalCopula)
-exportClasses(vineCopula, spCopula, stCopula, spVineCopula)
+exportClasses(asCopula, cqsCopula, neighbourhood, stNeighbourhood, empiricalCopula)
+exportClasses(spCopula, stCopula, vineCopula, spVineCopula, stVineCopula)
# wrappers to CDVine
exportClasses(BB1Copula, surBB1Copula, r90BB1Copula, r270BB1Copula)
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2013-05-21 08:26:21 UTC (rev 96)
+++ pkg/R/Classes.R 2013-05-24 15:08:18 UTC (rev 97)
@@ -50,7 +50,7 @@
return("Parameter and lower bound have non-equal length")
if (any(is.na(param) | param > upper | param < lower))
return("Parameter value out of bound")
- if (length(object at fixed >0)){
+ if (object at fixed != ""){
if(!("a" %in% object at fixed | "b" %in% object at fixed))
return("The slot fixed may only refer to \"a\" or \"b\".")
if ("a" %in% object at fixed & "b" %in% object at fixed)
@@ -208,6 +208,17 @@
setClassUnion("spVineCopula",c("mixedSpVineCopula","pureSpVineCopula"))
+#################################
+## Spatio-temporal Vine Copula ##
+#################################
+
+validStVineCopula <- function(object) {
+ return(validStCopula(object at stCop) & validObject(object at topCop))
+}
+
+setClass("stVineCopula", representation("copula", stCop="stCopula", topCop="copula"),
+ validity = validStVineCopula, contains=list("copula"))
+
########################################
## spatial classes providing the data ##
########################################
@@ -259,3 +270,32 @@
prediction="logical"),
validity = validNeighbourhood, contains = list("Spatial"))
+## 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)
+ 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])
+ return("Data and index have unequal number of rows.")
+ if (dimDists[2] != dimInd[2])
+ return("Data and index have unequal 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",
+ prediction="logical"),
+ validity = validStNeighbourhood, contains = list("ST"))
\ No newline at end of file
Modified: pkg/R/asCopula.R
===================================================================
--- pkg/R/asCopula.R 2013-05-21 08:26:21 UTC (rev 96)
+++ pkg/R/asCopula.R 2013-05-24 15:08:18 UTC (rev 97)
@@ -184,47 +184,36 @@
# method
# one of kendall or spearman according to the calculation of moa
-fitASC2.itau <- function(copula, data, estimate.variance) {
-tau <- cor(data,method="kendall")[1,2]
-esti <- fitASC2.moa(tau, data, method="itau")
-copula <- asCopula(esti)
-return(new("fitCopula",
- estimate = esti,
- var.est = matrix(NA),
- loglik = sum(log(dCopula(data, copula))),
- nsample = nrow(data),
- method = "Inversion of Kendall's tau and MLE",
- fitting.stats = list(convergence=as.integer(NA)),
- copula = copula))
+fitASC2.itau <- function(copula, data, estimate.variance, tau=NULL) {
+ if(is.null(tau))
+ tau <- VineCopula:::fasttau(data[,1],data[,2])
+ esti <- fitASC2.moa(tau, data, method="itau")
+ copula <- asCopula(esti)
+
+ new("fitCopula", estimate = esti, var.est = matrix(NA),
+ loglik = sum(log(dCopula(data, copula))), nsample = nrow(data),
+ method = "Inversion of Kendall's tau and MLE",
+ fitting.stats = list(convergence=as.integer(NA)), copula = copula)
}
-fitASC2.irho <- function(copula, data, estimate.variance){
-rho <- cor(data,method="spearman")[1,2]
-esti <- fitASC2.moa(rho, data, method="itau")
-copula <- asCopula(esti)
-return(new("fitCopula",
- estimate = esti,
- var.est = matrix(NA),
- loglik = sum(log(dCopula(data, copula))),
- nsample = nrow(data),
- method = "Inversion of Spearman's rho and MLE",
- fitting.stats = list(convergence=as.integer(NA)),
- copula = copula))
+fitASC2.irho <- function(copula, data, estimate.variance, rho=NULL){
+ if(is.null(rho))
+ rho <- cor(data,method="spearman")[1,2]
+ esti <- fitASC2.moa(rho, data, method="irho")
+ copula <- asCopula(esti)
+
+ new("fitCopula", estimate = esti, var.est = matrix(NA),
+ loglik = sum(log(dCopula(data, copula))), nsample = nrow(data),
+ method = "Inversion of Spearman's rho and MLE",
+ fitting.stats = list(convergence=as.integer(NA)), copula = copula)
}
fitASC2.moa <- function(moa, data, method="itau", tol=.Machine$double.eps^.5) {
smpl <- as.matrix(data)
+ iFun <- switch(method,
+ itau=function(p) iTauASC2(p,moa),
+ irho=function(p) iRhoASC2(p,moa))
- iTau <- function(p) {
- iTauASC2(p,moa)
- }
-
- iRho <- function(p) {
- iRhoASC2(p,moa)
- }
-
- iFun <- switch(method, itau=iTau, irho=iRho)
-
sec <- function (parameters) {
res <- NULL
for(param in parameters) {
@@ -235,9 +224,7 @@
b <- optimize(sec,c(-1,1), tol=tol)$minimum
- param <- c(iFun(b),b)
-
- return(param)
+ return(c(iFun(b),b))
}
# maximum log-likelihood estimation of a and b using optim
Modified: pkg/R/cqsCopula.R
===================================================================
--- pkg/R/cqsCopula.R 2013-05-21 08:26:21 UTC (rev 96)
+++ pkg/R/cqsCopula.R 2013-05-24 15:08:18 UTC (rev 97)
@@ -7,7 +7,7 @@
## ##
######################################################
-cqsCopula <- function (param=c(0,0), fixed=character(0)) {
+cqsCopula <- function (param=c(0,0), fixed="") {
new("cqsCopula", dimension = as.integer(2), parameters = param,
param.names = c("a", "b"), param.lowbnd = c(limA(param[2]),-1),
param.upbnd = c(1, 1),
@@ -272,8 +272,6 @@
copula=copula)
}
-
-
fitCQSec.moa <- function(moa, data, method="itau", tol=.Machine$double.eps^.5) {
smpl <- as.matrix(data)
@@ -298,9 +296,7 @@
b <- optimize(sec,c(-1,1), tol=tol)$minimum
- param <- c(iFun(b),b)
-
- return(param)
+ return(c(iFun(b),b))
}
# maximum log-likelihood estimation of a and b using optim
Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2013-05-21 08:26:21 UTC (rev 96)
+++ pkg/R/spCopula.R 2013-05-24 15:08:18 UTC (rev 97)
@@ -422,18 +422,7 @@
# cutoff -> maximal distance that should be considered for fitting
# bounds -> the bounds of the correlation function (typically c(0,1))
# method -> the measure of association, either "kendall" or "spearman"
-fitCorFun <- function(bins, degree=3, cutoff=NA, bounds=c(0,1),
- cor.method=NULL, weighted=FALSE) {
- if(is.null(cor.method)) {
- if(is.null(attr(bins,"cor.method")))
- stop("Neither the bins arguments has an attribute cor.method nor is the parameter cor.method provided.")
- else
- cor.method <- attr(bins,"cor.method")
- } else {
- if(!is.null(attr(bins,"cor.method")) && cor.method != attr(bins,"cor.method"))
- stop("The cor.method attribute of the bins argument and the argument cor.method do not match.")
- }
-
+fitCorFunSng <- function(bins, degree, cutoff, bounds, cor.method, weighted) {
if (weighted) {
bins <- as.data.frame(bins[c("np","meanDists","lagCor")])
if(!is.na(cutoff))
@@ -449,6 +438,9 @@
print(fitCor)
cat("Sum of squared residuals:",sum(fitCor$residuals^2),"\n")
+ if(cor.method=="fasttau")
+ cor.method <- "kendall"
+
function(x) {
if (is.null(x)) return(cor.method)
return(pmin(bounds[2], pmax(bounds[1],
@@ -456,7 +448,33 @@
}
}
+fitCorFun <- function(bins, degree=3, cutoff=NA, bounds=c(0,1),
+ cor.method=NULL, weighted=FALSE){
+ if(is.null(cor.method)) {
+ if(is.null(attr(bins,"cor.method")))
+ stop("Neither the bins arguments has an attribute cor.method nor is the parameter cor.method provided.")
+ else
+ cor.method <- attr(bins,"cor.method")
+ } else {
+ if(!is.null(attr(bins,"cor.method")) && cor.method != attr(bins,"cor.method"))
+ stop("The cor.method attribute of the bins argument and the argument cor.method do not match.")
+ }
+
+ if(is.null(nrow(bins$lagCor)))
+ return(fitCorFunSng(bins, degree, cutoff, bounds, cor.method, weighted))
+
+ degree <- rep(degree,length.out=nrow(bins$lagCor))
+ calcKTau <- list()
+ for(j in 1:nrow(bins$lagCor)) {
+ calcKTau[[paste("fun",j,sep="")]] <- fitCorFunSng(data.frame(meanDists=bins$meanDists,
+ lagCor=bins$lagCor[j,]),
+ degree[j], cutoff, bounds,
+ cor.method, weighted)
+ }
+ return(calcKTau)
+}
+
# towards b)
## loglikelihoods for a dynamic spatial copula
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2013-05-21 08:26:21 UTC (rev 96)
+++ pkg/R/spatialPreparation.R 2013-05-24 15:08:18 UTC (rev 97)
@@ -4,8 +4,8 @@
## ##
########################################################
-## neighbourhood constructor
-############################
+## spatial neighbourhood constructor
+####################################
neighbourhood <- function(data, distances, sp, dataLocs=NULL, index,
prediction, var) {
@@ -54,7 +54,7 @@
prediction=x at prediction)
}
-setMethod("[[", "neighbourhood", selectFromNeighbourhood)
+setMethod("[[", signature("neighbourhood","numeric","missing"), selectFromNeighbourhood)
## calculate neighbourhood from SpatialPointsDataFrame
@@ -84,7 +84,7 @@
allData <- NULL
for(i in 1:length(locations)) { # i <- 1
- tempDists <- spDistsN1(spData,locations[i,])
+ tempDists <- spDists(spData,locations[i,])
tempDists[tempDists < min.dist] <- Inf
spLocs <- order(tempDists)[1:(size-1)]
allLocs <- rbind(allLocs, spLocs)
Added: pkg/R/spatio-temporalPreparation.R
===================================================================
--- pkg/R/spatio-temporalPreparation.R (rev 0)
+++ pkg/R/spatio-temporalPreparation.R 2013-05-24 15:08:18 UTC (rev 97)
@@ -0,0 +1,98 @@
+###############################################################
+## ##
+## functions based on spacetime preparing the use of copulas ##
+## ##
+###############################################################
+
+## spatio-temporal neighbourhood constructor
+############################################
+
+stNeighbourhood <- function(data, distances, STxDF, ST=NULL,index,
+ prediction, var) {
+ data <- as.data.frame(data)
+ sizeN <- nrow(data)
+ dimDists <- dim(distances)
+
+ 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=".")
+ if (anyDuplicated(rownames(data))>0)
+ rownames <- 1:length(rownames)
+ new("stNeighbourhood", data=data, distances=distances, locations=STxDF,
+ dataLocs=ST, 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)])
+}
+
+## show
+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")
+}
+
+setMethod(show,signature("stNeighbourhood"),showStNeighbourhood)
+
+
+## calculate neighbourhood from SpatialPointsDataFrame
+
+# 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)))
+ stopifnot(min.dist>0 || prediction)
+
+ timeSpan <- min(t.lags)
+ if(missing(ST) && !prediction)
+ ST=stData
+ if(is.na(timeSteps))
+ timeSteps <- length(stData at time)+timeSpan
+
+ stopifnot(is(ST,"ST"))
+
+ nLocs <- length(ST at sp)*timeSteps
+
+ 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.")
+
+ if(prediction)
+ nghbrs <- getNeighbours(stData[,1], ST, var, spSize, prediction, min.dist)
+ else
+ nghbrs <- getNeighbours(stData[,1], var=var, size=spSize, min.dist=min.dist)
+
+ stNeighData <- NULL
+ stDists <- array(NA,c(nLocs,(spSize-1)*length(t.lags),2))
+ 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]],
+ 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,],
+ tmpInst+t,var]@data[[1]],
+ ncol=spSize-1, byrow=T))
+ tmpInd <- cbind(tmpInd, matrix(rep(tmpInst+t,spSize-1),ncol=spSize-1))
+ }
+ 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
+ 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)),
+ byrow=T, ncol=length(t.lags)*(spSize-1))
+ stInd[(i-1)*timeSteps+1:timeSteps,,2] <- tmpInd
+ }
+
+ if (prediction)
+ dataLocs <- stData
+ else
+ dataLocs <- NULL
+ return(stNeighbourhood(as.data.frame(stNeighData), stDists, stData, ST,
+ stInd, prediction, var))
+}
\ No newline at end of file
Modified: pkg/R/stCopula.R
===================================================================
--- pkg/R/stCopula.R 2013-05-21 08:26:21 UTC (rev 96)
+++ pkg/R/stCopula.R 2013-05-24 15:08:18 UTC (rev 97)
@@ -1,19 +1,10 @@
-# constructor
-# dimension = "integer" set to 2
-# parameters = "numeric" set of parameters
-# param.names = "character" appropriate names
-# param.lowbnd = "numeric" appropriate lower bounds
-# param.upbnd = "numeric" appropriate upper bounds
-# fullname = "character" messgae printed with "show"
-# components="list" list of copulas (will be automatically supplemented
-# by the independent copula)
-# distances="numeric" the linking distances + the range (will be assigned
-# to the independent copula)
-# unit="character" measurement unit of distance
-# depFun="function" a optional spatial dependence function providing
-# Kendalls tau or Spearman's rho to calib* or exact
-# parameters
+######################################
+## Spatio-Temporal Bivariate Copula ##
+######################################
+## constructor ##
+#################
+
stCopula <- function(components, distances, t.lags, stDepFun, unit="m", t.res="day") {
spCopList <- list()
@@ -40,7 +31,9 @@
spCopList=spCopList, t.lags=t.lags, t.res=t.res)
}
-## show method
+## show method ##
+#################
+
showStCopula <- function(object) {
cat(object at fullname, "\n")
cat("Dimension: ", object at dimension, "\n")
@@ -55,270 +48,126 @@
setMethod("show", signature("stCopula"), showStCopula)
-## spatial copula jcdf ##
+## spatial copula cdf ##
+########################
-# TODO: add again block support to the spatio-temporal copula
-# u
-# list containing two column matrix providing the transformed pairs, their respective
-# separation distances and time steps
-pStCopula <- function (u, copula) {
- if (!is.list(u) || !length(u)>=3) stop("Point pairs need to be provided with their separating spatial and temproal distances as a list.")
+pStCopula <- function (u, copula, h) {
+ stopifnot(ncol(h)==2)
+ stopifnot(nrow(h)==1 || nrow(h)==nrow(u))
- if(!is.matrix(u[[1]])) u[[1]] <- matrix(u[[1]],ncol=2)
- n <- nrow(u[[1]])
- h <- u[[2]]
- t.dist <- u[[3]]
+ n <- nrow(u)
+ tDist <- unique(h[,2])
- if(length(u)==4) {
- block <- u[[4]]
- if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
- } else block <- 1
-
- if(any(is.na(match(t.dist,copula at t.lags))))
+ if(any(is.na(match(tDist,copula at t.lags))))
stop("Prediction time(s) do(es) not math the modelled time slices.")
- if(length(h)>1 & length(h)!=n)
- stop("The spatial distance vector must either be of same length as rows in the data pairs or a single value.")
- if(length(t.dist)>1 & length(t.dist)!=n)
- stop("The temporal distances vector must either be of same length as rows in the data pairs or a single value.")
-
- if (length(t.dist)==1) {
- res <- pSpCopula(copula at spCopList[[match(t.dist,copula at t.lags)]],
- list(u[[1]], h))
+
+ if (length(tDist)==1) {
+ res <- pSpCopula(u, copula at spCopList[[match(tDist, copula at t.lags)]], h[,1])
} else {
- if(length(h)==1) h <- rep(h,n)
- res <- NULL
- for(i in 1:(n%/%block)) {
- res <- rbind(res, pSpCopula(copula at spCopList[[match(t.dist[i],copula at t.lags)]],
- list(u[[1]][((i-1)*block+1):(i*block),], h[i*block])))
+ res <- numeric(n)
+ 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])
}
}
-
- return(res)
+ res
}
-setMethod("pCopula", signature("numeric","stCopula"), pStCopula)
+setMethod(pCopula, signature("numeric","stCopula"),
+ function(u, copula, log, ...) pStCopula(matrix(u,ncol=2), copula, ...))
+setMethod(pCopula, signature("matrix","stCopula"), pStCopula)
-
## spatial Copula density ##
+############################
-# u
-# three column matrix providing the transformed pairs and their respective
-# separation distances
-dStCopula <- function (u, copula) {
- if (!is.list(u) || !length(u)>=3) stop("Point pairs need to be provided with their separating spatial and temproal distances as a list.")
+dStCopula <- function (u, copula, log, h) {
+ stopifnot(ncol(h)==2)
+ stopifnot(nrow(h)==1 || nrow(h)==nrow(u))
- if(!is.matrix(u[[1]])) u[[1]] <- matrix(u[[1]],ncol=2)
- n <- nrow(u[[1]])
- h <- u[[2]]
- t.dist <- u[[3]]
+ n <- nrow(u)
+ tDist <- unique(h[,2])
- if(length(u)==4) {
- block <- u[[4]]
- if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
- } else block <- 1
-
- if(any(is.na(match(t.dist,copula at t.lags))))
+ if(any(is.na(match(tDist,copula at t.lags))))
stop("Prediction time(s) do(es) not math the modelled time slices.")
- if(length(h)>1 & length(h)!=n)
- stop("The spatial distance vector must either be of same length as rows in the data pairs or a single value.")
- if(length(t.dist)>1 & length(t.dist)!=n)
- stop("The temporal distances vector must either be of same length as rows in the data pairs or a single value.")
- if (length(t.dist)==1) {
- res <- dSpCopula(copula at spCopList[[match(t.dist,copula at t.lags)]],
- list(u[[1]], h))
+ if (length(tDist)==1) {
+ res <- dSpCopula(u, copula at spCopList[[match(tDist, copula at t.lags)]], log, h[,1])
} else {
- if(length(h)==1) h <- rep(h,n)
- res <- NULL
- for(i in 1:(n%/%block)) {
- res <- rbind(res, dSpCopula(copula at spCopList[[match(t.dist[i],copula at t.lags)]],
- list(u[[1]][((i-1)*block+1):(i*block),], h[i*block])))
+ res <- numeric(n)
+ 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])
}
}
-
- return(res)
+ res
}
-setMethod("dCopula", signature("list","stCopula"), dStCopula)
+setMethod(dCopula, signature("numeric","stCopula"),
+ function(u, copula, log, ...) dStCopula(matrix(u,ncol=2), copula, log=log, ...))
+setMethod(dCopula, signature("matrix","stCopula"), dStCopula)
## partial derivatives ##
-## dduSpCopula
-###############
+## dduSpCopula ##
+#################
-dduStCopula <- function (u, copula) {
- if (!is.list(u) || !length(u)>=3) stop("Point pairs need to be provided with their separating spatial and temproal distances as a list.")
+dduStCopula <- function (u, copula, h) {
+ stopifnot(ncol(h)==2)
+ stopifnot(nrow(h)==1 || nrow(h)==nrow(u))
- if(!is.matrix(u[[1]])) u[[1]] <- matrix(u[[1]],ncol=2)
- n <- nrow(u[[1]])
- h <- u[[2]]
- t.dist <- u[[3]]
+ n <- nrow(u)
+ tDist <- unique(h[,2])
- if(length(u)==4) {
- t.block <- u[[4]]
- if (n%%t.block != 0) stop("The block size is not a multiple of the data length:",n)
- } else t.block <- 1
-
- if(any(is.na(match(t.dist,copula at t.lags))))
+ if(any(is.na(match(tDist,copula at t.lags))))
stop("Prediction time(s) do(es) not math the modelled time slices.")
- if(length(h)>1 & length(h)!=n)
- stop("The spatial distance vector must either be of same length as rows in the data pairs or a single value.")
- if(length(t.dist)>1 & length(t.dist)!=n)
- stop("The temporal distances vector must either be of same length as rows in the data pairs or a single value.")
- if (length(t.dist)==1) {
- res <- dduSpCopula(copula at spCopList[[match(t.dist,copula at t.lags)]],
- list(u[[1]], h, block=t.block))
+ if (length(tDist)==1) {
+ res <- dduSpCopula(u, copula at spCopList[[match(tDist, copula at t.lags)]], h[,1])
} else {
- if(length(h)==1) h <- rep(h,n)
- res <- NULL
- for(i in 1:(n%/%t.block)) {
- cop <- copula at spCopList[[match(t.dist[i*t.block],copula at t.lags)]]
- tmpPair <- u[[1]][((i-1)*t.block+1):(i*t.block),]
- res <- rbind(res, dduSpCopula(cop, list(tmpPair,h[i*t.block])))
+ res <- numeric(n)
+ 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])
}
}
-
- return(res)
+ res
}
-setMethod("dduCopula", signature("list","stCopula"), dduStCopula)
+setMethod("dduCopula", signature("numeric","stCopula"),
+ function(u, copula, ...) dduStCopula(matrix(u,ncol=2), copula, ...))
+setMethod("dduCopula", signature("matrix","stCopula"), dduStCopula)
-## ddvSpCopula
-###############
-ddvStCopula <- function (u, copula) {
- if (!is.list(u) || !length(u)>=3) stop("Point pairs need to be provided with their separating spatial and temproal distances as a list.")
+## ddvSpCopula ##
+#################
+
+ddvStCopula <- function (u, copula, h) {
+ stopifnot(ncol(h)==2)
+ stopifnot(nrow(h)==1 || nrow(h)==nrow(u))
- if(!is.matrix(u[[1]])) u[[1]] <- matrix(u[[1]],ncol=2)
- n <- nrow(u[[1]])
- h <- u[[2]]
- t.dist <- u[[3]]
+ n <- nrow(u)
+ tDist <- unique(h[,2])
- if(length(u)==4) {
- t.block <- u[[4]]
- if (n%%t.block != 0) stop("The block size is not a multiple of the data length:",n)
- } else t.block <- 1
-
- if(any(is.na(match(t.dist,copula at t.lags))))
+ if(any(is.na(match(tDist,copula at t.lags))))
stop("Prediction time(s) do(es) not math the modelled time slices.")
- if(length(h)>1 & length(h)!=n)
- stop("The spatial distance vector must either be of same length as rows in the data pairs or a single value.")
- if(length(t.dist)>1 & length(t.dist)!=n)
- stop("The temporal distances vector must either be of same length as rows in the data pairs or a single value.")
- if (length(t.dist)==1) {
- res <- ddvSpCopula(copula at spCopList[[match(t.dist,copula at t.lags)]],
- list(u[[1]], h, block=t.block))
+ if (length(tDist)==1) {
+ res <- ddvSpCopula(u, copula at spCopList[[match(tDist,copula at t.lags)]], h[,1])
} else {
- if(length(h)==1) h <- rep(h,n)
- res <- NULL
- for(i in 1:(n%/%t.block)) {
- cop <- copula at spCopList[[match(t.dist[i*t.block],copula at t.lags)]]
- tmpPair <- u[[1]][((i-1)*t.block+1):(i*t.block),]
- res <- rbind(res, ddvSpCopula(cop, list(tmpPair,h[i*t.block])))
+ res <- numeric(n)
+ 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])
}
}
-
- return(res)
+ res
}
-setMethod("ddvCopula", signature("list","stCopula"), ddvStCopula)
-
-# #############
-# ## ##
-# ## FITTING ##
-# ## ##
-# #############
-#
-# # two models:
-# # 1) Kendall's tau driven:
-# # fit curve through emp. Kendall's tau values, identify validity ranges for
-# # copula families deriving parameters from the fit, fade from one family to
-# # another at borders
-# # 2) convex-linear combination of copulas:
-# # fit one per lag, fade from one to another
-#
-# # towards the first model:
-#
-# # INPUT: the stBinning
-# # steps
-# # a) fit a curve
-# # b) estimate bivariate copulas per lag (limited to those with some 1-1-relation
-# # to Kendall's tau')
-# # INTERMEDIATE RESULT
-# # c) select best fits based on ... e.g. log-likelihood, visual inspection
-# # d) compose bivariate copulas to one spatial copula
-# # OUTPUT: a spatial copula parametrixued by distance through Kendall's tau
-#
-#
-# # towards b)
-# # bins -> typically output from calcBins
-# # calcTau -> a function on distance providing Kendall's tau estimates,
-# # families -> a vector of dummy copula objects of each family to be considered
-# # DEFAULT: c(normal, t_df=4, clayton, frank, gumbel
-# loglikByCopulasLags <- function(bins, calcTau, families=c(normalCopula(0), tCopula(0,dispstr="un"),
-# claytonCopula(0), frankCopula(1), gumbelCopula(1))) {
-# loglik <- NULL
-# for (cop in families) {
-# print(cop)
-# tmploglik <- NULL
-# for(i in 1:length(bins$meanDists)) {
-# cop at parameters[1] <- iTau(cop,tau=calcTau(bins$meanDists[i]))
-# tmploglik <- c(tmploglik, sum(log(dCopula(bins$lagData[[i]],cop))))
-# }
-# loglik <- cbind(loglik, tmploglik)
-# }
-#
-# colnames(loglik) <- sapply(families, function(x) class(x)[1])
-#
-# return(loglik)
-# }
-#
-# # towards d)
-# composeSpCop <- function(bestFit, families, bins, calcCor) {
-# nfits <- length(bestFit)
-# gaps <- which(diff(bestFit)!=0)
-#
-# if(missing(calcCor)) noCor <- nfits
-# else noCor <- min(which(calcCor(bins$meanDists)<=0), nfits)
-#
-# breaks <- sort(c(gaps, gaps+1, noCor))
-# breaks <- breaks[breaks<noCor]
-#
-# cops <- as.list(families[bestFit[breaks]])
-#
-# breaks <- unique(c(breaks, min(nfits,rev(breaks)[1]+1)))
-# distances <- bins$meanDists[breaks]
-#
-# if(length(breaks)==length(cops)) {
-# warning("Lags do not cover the entire range with positive correlation. ",
-# "An artifical one fading towards the indepCopula has been added.")
-# distances <- c(distances, rev(distances)[1]+diff(bins$meanDists[nfits+c(-1,0)]))
-# }
-#
-# if(missing(calcCor)) return(spCopula(components=cops, distances=distances, unit="m"))
-# else return(spCopula(components=cops, distances=distances, unit="m", spDepFun=calcCor))
-# }
-#
-# # in once
-#
-# # bins -> typically output from calcBins
-# # cutoff -> maximal distance that should be considered for fitting
-# # families -> a vector of dummy copula objects of each family to be considered
-# # DEFAULT: c(normal, t_df=4, clayton, frank, gumbel
-# # ...
-# # type -> the type of curve (by now only polynominals are supported)
-# # degree -> the degree of the polynominal
-# # bounds -> the bounds of the correlation function (typically c(0,1))
-# # method -> the measure of association, either "kendall" or "spearman"
-# fitSpCopula <- function(bins, cutoff=NA, families=c(normalCopula(0), tCopula(0,dispstr="un"),
-# claytonCopula(0), frankCopula(1), gumbelCopula(1)), ...) {
-# calcTau <- fitCorFun(bins, cutoff=cutoff, ...)
-# loglik <- loglikByCopulasLags(bins, calcTau=calcTau, families=families)
-#
-# bestFit <- apply(apply(loglik, 1, rank),2,function(x) which(x==length(families)))
-#
-# return(composeSpCop(bestFit, families, bins, calcTau))
-# }
+setMethod("ddvCopula", signature("numeric","stCopula"),
+ function(u, copula, ...) ddvStCopula(matrix(u,ncol=2), copula, ...))
+setMethod("ddvCopula", signature("matrix","stCopula"), ddvStCopula)
\ No newline at end of file
Added: pkg/R/stVineCopula.R
===================================================================
--- pkg/R/stVineCopula.R (rev 0)
+++ pkg/R/stVineCopula.R 2013-05-24 15:08:18 UTC (rev 97)
@@ -0,0 +1,202 @@
+#########################################
+## methods for the spatial vine copula ##
+#########################################
+
+## constructor ##
+#################
+
+stVineCopula <- function(stCop, topCop) {
+ stopifnot(is(stCop,"stCopula"))
+
+ new("stVineCopula", dimension = as.integer(topCop at dimension+1),
+ parameters=numeric(), param.names = character(), param.lowbnd = numeric(),
+ param.upbnd = numeric(),
+ fullname = paste("Spatio-temporal vine copula family with 1 spatio-temporal tree."),
+ stCop=stCop, topCop=topCop)
+}
+
+## show ##
+##########
+
+showStVineCopula <- function(object) {
+ dim <- object at dimension
+ cat(object at fullname, "\n")
+ cat("Dimension: ", dim, "\n")
+}
+
+setMethod("show", signature("stVineCopula"), showStVineCopula)
+
+## density ##
+#############
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/spcopula -r 97
More information about the spcopula-commits
mailing list