[spcopula-commits] r127 - in pkg: . R data demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 19 20:33:05 CET 2014
Author: ben_graeler
Date: 2014-02-19 20:32:57 +0100 (Wed, 19 Feb 2014)
New Revision: 127
Modified:
pkg/NAMESPACE
pkg/R/Classes.R
pkg/R/spCopula.R
pkg/R/spatialPreparation.R
pkg/R/spatio-temporalPreparation.R
pkg/R/stCopula.R
pkg/data/spCopDemo.RData
pkg/demo/spCopula.R
pkg/man/calcBins.Rd
pkg/man/condCovariate.Rd
pkg/man/condStCoVarVine.Rd
pkg/man/condStVine.Rd
pkg/man/fitCorFun.Rd
pkg/man/fitSpCopula.Rd
pkg/man/getStNeighbours.Rd
pkg/man/loglikByCopulasLags.Rd
pkg/man/reduceNeighbours.Rd
pkg/man/spCopDemo.Rd
pkg/man/stCoVarVineCopula.Rd
pkg/man/stCopPredict.Rd
pkg/man/stCopula-class.Rd
pkg/man/stCopula.Rd
pkg/man/stNeighbourhood.Rd
pkg/man/stVineCopula.Rd
Log:
- redesign of spatial and spatio-temporal bins (i.e. dropping the lagData entry)
- renaming of arguments from t.lags to tlags and slots t.res to tres
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/NAMESPACE 2014-02-19 19:32:57 UTC (rev 127)
@@ -32,7 +32,7 @@
export(reduceNeighbours)
# fitting
-export(fitCorFun, loglikByCopulasLags, fitSpCopula, composeSpCopula)
+export(fitCorFun, loglikByCopulasLags, loglikByCopulasStLags, fitSpCopula, composeSpCopula)
export(tailDepFun, lowerTailDepFun, upperTailDepFun)
export(empTailDepFun, lowerEmpTailDepFun, upperEmpTailDepFun)
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/R/Classes.R 2014-02-19 19:32:57 UTC (rev 127)
@@ -163,14 +163,14 @@
############################
validStCopula <- function(object) {
- if(length(object at t.lags) != length(object at spCopList)) return("The length of the temporal distance vector must equal the number of spatial copulas.")
+ if(length(object at tlags) != length(object at spCopList)) return("The length of the temporal distance vector must equal the number of spatial copulas.")
return(TRUE) # validity of any spCopula in spCopList is tested by the constructor, I believe
}
setClass("stCopula", representation = representation("copula",
spCopList="list",
- t.lags="numeric",
- t.res="character"),
+ tlags="numeric",
+ tres="character"),
validity = validStCopula, contains = list("copula"))
###############################################
Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/R/spCopula.R 2014-02-19 19:32:57 UTC (rev 127)
@@ -448,7 +448,7 @@
}
}
-fitCorFun <- function(bins, degree=3, cutoff=NA, bounds=c(0,1),
+fitCorFun <- function(bins, degree=3, cutoff=NA, tlags, bounds=c(0,1),
cor.method=NULL, weighted=FALSE){
if(is.null(cor.method)) {
if(is.null(attr(bins,"cor.method")))
@@ -460,9 +460,10 @@
stop("The cor.method attribute of the bins argument and the argument cor.method do not match.")
}
- if(is.null(nrow(bins$lagCor)))
+ if(is.null(nrow(bins$lagCor))) # the spatial case
return(fitCorFunSng(bins, degree, cutoff, bounds, cor.method, weighted))
-
+
+ # the spatio-temporal case
degree <- rep(degree,length.out=nrow(bins$lagCor))
calcKTau <- list()
for(j in 1:nrow(bins$lagCor)) {
@@ -471,14 +472,21 @@
degree[j], cutoff, bounds,
cor.method, weighted)
}
- return(calcKTau)
+
+ corFun <- function(h, time, tlags=sort(tlags,decreasing=TRUE)) {
+ t <- which(tlags==time)
+ calcKTau[[time]](h)
+ }
+
+ attr(corFun, "tlags") <- sort(tlags, decreasing=TRUE)
+ return(corFun)
}
# towards b)
## loglikelihoods for a dynamic spatial copula
-loglikByCopulasLags.dyn <- function(bins, families, calcCor) {
+loglikByCopulasLags.dyn <- function(bins, lagData, families, calcCor) {
moa <- switch(calcCor(NULL),
kendall=function(copula, h) iTau(copula, calcCor(h)),
spearman=function(copula, h) iRho(copula, calcCor(h)),
@@ -496,18 +504,18 @@
if(class(cop)!="indepCopula") {
if(class(cop) == "asCopula") {
cop <- switch(calcCor(NULL),
- kendall=fitASC2.itau(cop, bins$lagData[[i]],
+ kendall=fitASC2.itau(cop, lagData[[i]],
tau=calcCor(bins$meanDists[i]))@copula,
- spearman=fitASC2.irho(cop, bins$lagData[[i]],
+ spearman=fitASC2.irho(cop, lagData[[i]],
rho=calcCor(bins$meanDists[i]))@copula,
stop(paste(calcCor(NULL), "is not yet supported.")))
param <- cop at parameters
} else {
if(class(cop) == "cqsCopula") {
cop <- switch(calcCor(NULL),
- kendall=fitCQSec.itau(cop, bins$lagData[[i]],
+ kendall=fitCQSec.itau(cop, lagData[[i]],
tau=calcCor(bins$meanDists[i]))@copula,
- spearman=fitCQSec.irho(cop, bins$lagData[[i]],
+ spearman=fitCQSec.irho(cop, lagData[[i]],
rho=calcCor(bins$meanDists[i]))@copula,
stop(paste(calcCor(NULL), "is not yet supported.")))
param <- cop at parameters
@@ -522,7 +530,7 @@
if(any(is.na(param)))
tmploglik <- c(tmploglik, NA)
else
- tmploglik <- c(tmploglik, sum(dCopula(bins$lagData[[i]], cop, log=T)))
+ tmploglik <- c(tmploglik, sum(dCopula(lagData[[i]], cop, log=T)))
tmpCop <- append(tmpCop, cop)
setTxtProgressBar(pb, i)
}
@@ -537,12 +545,12 @@
}
## loglikelihoods for a static spatial copula
-loglikByCopulasLags.static <- function(bins, families) {
+loglikByCopulasLags.static <- function(lagData, families) {
fits <-lapply(families,
function(cop) {
cat(cop at fullname,"\n")
- lapply(bins$lagData,
+ lapply(lagData,
function(x) {
tryCatch(fitCopula(cop, x, estimate.variance = FALSE),
error=function(e) return(NA))
@@ -571,28 +579,36 @@
return(list(loglik=loglik, copulas=copulas))
}
-# 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, families=c(normalCopula(0),
- tCopula(0,dispstr="un"),
- claytonCopula(0), frankCopula(1),
- gumbelCopula(1)),
- calcCor) {
- bins$lagData <- lapply(bins$lagData,
- function(pairs) {
- bool <- !is.na(pairs[,1]) & !is.na(pairs[,2])
- pairs[bool,]
- })
+##
+
+loglikByCopulasLags <- function(bins, data, families=c(normalCopula(),
+ tCopula(),
+ claytonCopula(), frankCopula(),
+ gumbelCopula()),
+ calcCor, lagSub=1:length(bins$meanDists)) {
+ var <- attr(bins, "variable")
+ lagData <- lapply(bins$lags[lagSub],
+ function(x) {
+ as.matrix((cbind(data[x[,1], var]@data,
+ data[x[,2], var]@data)))
+ })
+
+ lagData <- lapply(lagData,
+ function(pairs) {
+ bool <- !is.na(pairs[,1]) & !is.na(pairs[,2])
+ pairs[bool,]
+ })
+
if(missing(calcCor))
- return(loglikByCopulasLags.static(bins, families))
+ return(loglikByCopulasLags.static(lagData, families))
else
- return(loglikByCopulasLags.dyn(bins, families, calcCor))
+ return(loglikByCopulasLags.dyn(lapply(bins, function(x) x[lagSub]),
+ lagData, families, calcCor))
}
+
# towards d)
composeSpCopula <- function(bestFit, families, bins, calcCor, range=max(bins$meanDists)) {
nFits <- length(bestFit)
@@ -626,12 +642,12 @@
# 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,
+fitSpCopula <- function(bins, data, cutoff=NA,
families=c(normalCopula(), tCopula(),
claytonCopula(), frankCopula(),
gumbelCopula()), ...) {
calcCor <- fitCorFun(bins, cutoff=cutoff, ...)
- loglik <- loglikByCopulasLags(bins, families, calcCor)
+ loglik <- loglikByCopulasLags(bins, data, families, calcCor)
bestFit <- apply(apply(loglik$loglik, 1, rank),2,
function(x) which(x==length(families)))
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/R/spatialPreparation.R 2014-02-19 19:32:57 UTC (rev 127)
@@ -168,8 +168,10 @@
abline(h=c(-min(lagCor),0,min(lagCor)),col="grey")
}
- res <- list(np=np, meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=lags)
+# res <- list(np=np, meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=lags)
+ res <- list(np=np, meanDists = mDists, lagCor=lagCor, lags=lags)
attr(res,"cor.method") <- cor.method
+ attr(res,"variable") <- var
return(res)
}
@@ -177,7 +179,7 @@
# calc bins from a (conditional) neighbourhood
-calcNeighBins <- function(data, var=names(data), nbins=9, boundaries=NA,
+calcNeighBins <- function(data, var=data at var, nbins=9, boundaries=NA,
cutoff=NA, cor.method="kendall", plot=TRUE) {
dists <- data at distances
@@ -199,7 +201,6 @@
np <- numeric(nbins)
moa <- numeric(nbins)
meanDists <- numeric(nbins)
- lagData <- NULL
data <- as.matrix(data at data)
@@ -211,7 +212,7 @@
pairs <- rbind(pairs, data[bools[,col],c(1,1+col)])
}
- lagData <- append(lagData, list(pairs))
+# lagData <- append(lagData, list(pairs))
moa[i] <- corFun(pairs)
meanDists[i] <- mean(dists[bools])
np[i] <- sum(bools)
@@ -223,8 +224,11 @@
abline(h=c(-min(moa),0,min(moa)),col="grey")
}
- res <- list(np=np, meanDists = meanDists, lagCor=moa, lagData=lagData)
+# res <- list(np=np, meanDists = meanDists, lagCor=moa, lagData=lagData)
+ res <- list(np=np, meanDists = meanDists, lagCor=moa)
attr(res,"cor.method") <- switch(cor.method, fasttau="kendall", cor.method)
+ attr(res,"variable") <- var
+
return(res)
}
Modified: pkg/R/spatio-temporalPreparation.R
===================================================================
--- pkg/R/spatio-temporalPreparation.R 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/R/spatio-temporalPreparation.R 2014-02-19 19:32:57 UTC (rev 127)
@@ -63,13 +63,13 @@
# returns an neighbourhood object
##################################
-getStNeighbours <- function(stData, ST, spSize=4, t.lags=-(0:2),
+getStNeighbours <- function(stData, ST, spSize=4, tlags=-(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)
+ timeSpan <- min(tlags)
if(missing(ST) && !prediction)
ST=geometry(stData)
@@ -101,7 +101,7 @@
timeSteps <- length(stData at time)+timeSpan
}
- nStNeighs <- (spSize-1)*length(t.lags)
+ nStNeighs <- (spSize-1)*length(tlags)
stNeighData <- matrix(NA, nLocs, nStNeighs + 1 + length(coVar))
stDists <- array(NA,c(nLocs, nStNeighs, 2))
@@ -126,8 +126,8 @@
tmpInd <- matrix(rep(timeInst, spSize), ncol=spSize)
- for (j in 2:length(t.lags)) {
- t <- t.lags[j]
+ for (j in 2:length(tlags)) {
+ t <- tlags[j]
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)
@@ -135,15 +135,15 @@
}
# store spatial distances
- stDists[spInd,,1] <- matrix(rep(nghbrs at distances[i,], timeSteps*length(t.lags)),
+ stDists[spInd,,1] <- matrix(rep(nghbrs at distances[i,], timeSteps*length(tlags)),
byrow=T, ncol=nStNeighs)
# store temporal distances
- stDists[spInd,,2] <- matrix(rep(rep(t.lags,each=spSize-1), timeSteps),
+ stDists[spInd,,2] <- matrix(rep(rep(tlags,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)),
+ stInd[spInd,,1] <- matrix(rep(c(nghbrs at index[i, ], rep(nghbrs at index[i, -1], length(tlags)-1)),
timeSteps), ncol = nStNeighs + 1, byrow = T)
# store time indices
@@ -230,7 +230,7 @@
}
## to be redone
-# calcStNeighBins <- function(data, var="uniPM10", nbins=9, t.lags=-(0:2),
+# calcStNeighBins <- function(data, var="uniPM10", nbins=9, tlags=-(0:2),
# boundaries=NA, cutoff=NA, cor.method="fasttau") {
# dists <- data at distances[,,1]
#
@@ -248,8 +248,8 @@
# }
#
# lagData <- NULL
-# for(t.lag in t.lags) { # t.lag <- 0
-# tBool <- data at distances[,,2]==t.lag
+# for(tlag in tlags) { # tlag <- 0
+# tBool <- data at distances[,,2]==tlag
# tmpLagData <- NULL
# for(i in 1:nbins) { # i <- 1
# sBool <- (dists <= boundaries[i + 1] & dists > boundaries[i])
@@ -341,10 +341,10 @@
# 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
+# other -> temporal indexing as in spacetime/xts, the parameter tlags is set to 0 in this case.
+# tlags: numeric -> temporal shifts between obs
calcStBins <- function(data, var, nbins=15, boundaries=NA, cutoff=NA,
- instances=NA, t.lags=-(0:2), ...,
+ instances=NA, tlags=-(0:2), ...,
cor.method="fasttau", plot=FALSE) {
if(is.na(cutoff))
cutoff <- spDists(coordinates(t(data at sp@bbox)))[1,2]/3
@@ -363,11 +363,14 @@
}
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)
+ for (tlag in rev(tlags)) {
+ if(is.na(instances))
+ smplInd <- max(1,1-min(tlags)):min(lengthTime,lengthTime-min(tlags))
+ else
+ smplInd <- sort(sample(x=max(1,1-min(tlags)):min(lengthTime,lengthTime-min(tlags)),
+ size=min(instances,lengthTime-max(abs(tlags)))))
+
+ tempIndices <- cbind(smplInd+tlag, tempIndices)
tempIndices <- cbind(smplInd, tempIndices)
}
}
@@ -412,8 +415,10 @@
abline(h=c(-min(lagCor),0,min(lagCor)),col="grey")
}
- res <- list(meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=list(sp=spIndices, time=tempIndices))
+ # res <- list(meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=list(sp=spIndices, time=tempIndices))
+ res <- list(meanDists = mDists, lagCor=lagCor, lags=list(sp=spIndices, time=tempIndices))
attr(res,"cor.method") <- cor.method
+ attr(res, "variable") <- var
return(res)
}
Modified: pkg/R/stCopula.R
===================================================================
--- pkg/R/stCopula.R 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/R/stCopula.R 2014-02-19 19:32:57 UTC (rev 127)
@@ -5,23 +5,23 @@
## constructor ##
#################
-stCopula <- function(components, t.lags, distances=NA, stDepFun, unit="m", t.res="day") {
+stCopula <- function(components, tlags, distances=NA, stDepFun, unit="m", tres="day") {
if(all(sapply(components, function(x) class(x)=="spCopula"))) {
if(length(unique(sapply(components, function(x) x at unit))) >1 )
stop("All spatial copulas need to have the same distance unit.")
- stopifnot(length(t.lags) == length(components))
+ stopifnot(length(tlags) == length(components))
spCopList <- components
} else {
spCopList <- list()
if(!missing(stDepFun)) {
getSpCop <- function(comp,dist,time) spCopula(comp, dist,
- spDepFun=function(h) stDepFun(h,time), unit)
- for(i in 1:length(t.lags)){
+ spDepFun=function(h) stDepFun(h, time, 1:length(tlags)), unit)
+ for(i in 1:length(tlags)){
spCopList <- append(spCopList, getSpCop(components[[i]], distances[[i]], i))
}
} else {
- for(i in 1:length(t.lags)){
+ for(i in 1:length(tlags)){
spCopList <- append(spCopList, spCopula(components[[i]], distances[[i]], unit=unit))
}
}
@@ -35,7 +35,7 @@
new("stCopula", dimension=as.integer(2), parameters=param, param.names=param.names,
param.lowbnd=param.low, param.upbnd=param.up,
fullname="Spatio-Temporal Copula: distance and time dependent convex combination of bivariate copulas",
- spCopList=spCopList, t.lags=t.lags, t.res=t.res)
+ spCopList=spCopList, tlags=tlags, tres=tres)
}
## show method ##
@@ -47,8 +47,8 @@
cat("Copulas:\n")
for (i in 1:length(object at spCopList)) {
cmpCop <- object at spCopList[[i]]
- cat(" ", cmpCop at fullname, "at", object at t.lags[i],
- paste("[",object at t.res,"]",sep=""), "\n")
+ cat(" ", cmpCop at fullname, "at", object at tlags[i],
+ paste("[",object at tres,"]",sep=""), "\n")
show(cmpCop)
}
}
@@ -65,16 +65,16 @@
n <- nrow(u)
tDist <- unique(h[,2])
- if(any(is.na(match(tDist,copula at t.lags))))
+ if(any(is.na(match(tDist,copula at tlags))))
stop("Prediction time(s) do(es) not math the modelled time slices.")
if (length(tDist)==1) {
- res <- pSpCopula(u, copula at spCopList[[match(tDist, copula at t.lags)]], h[,1])
+ res <- pSpCopula(u, copula at spCopList[[match(tDist, copula at tlags)]], h[,1])
} else {
res <- numeric(n)
for(t in tDist) {
tmpInd <- h[,2]==t
- tmpCop <- copula at spCopList[[match(t, copula at t.lags)]]
+ tmpCop <- copula at spCopList[[match(t, copula at tlags)]]
res[tmpInd] <- pSpCopula(u[tmpInd,,drop=F], tmpCop, h[tmpInd,1])
}
}
@@ -95,16 +95,16 @@
n <- nrow(u)
tDist <- unique(h[,2])
- if(any(is.na(match(tDist,copula at t.lags))))
+ if(any(is.na(match(tDist,copula at tlags))))
stop("Prediction time(s) do(es) not math the modelled time slices.")
if (length(tDist)==1) {
- res <- dSpCopula(u, copula at spCopList[[match(tDist, copula at t.lags)]], log, h[,1])
+ res <- dSpCopula(u, copula at spCopList[[match(tDist, copula at tlags)]], log, h[,1])
} else {
res <- numeric(n)
for(t in tDist) {
tmpInd <- h[,2]==t
- tmpCop <- copula at spCopList[[match(t, copula at t.lags)]]
+ tmpCop <- copula at spCopList[[match(t, copula at tlags)]]
res[tmpInd] <- dSpCopula(u[tmpInd,,drop=F], tmpCop, log, h[tmpInd,1])
}
}
@@ -128,16 +128,16 @@
n <- nrow(u)
tDist <- unique(h[,2])
- if(any(is.na(match(tDist,copula at t.lags))))
+ if(any(is.na(match(tDist,copula at tlags))))
stop("Prediction time(s) do(es) not math the modelled time slices.")
if (length(tDist)==1) {
- res <- dduSpCopula(u, copula at spCopList[[match(tDist, copula at t.lags)]], h[,1])
+ res <- dduSpCopula(u, copula at spCopList[[match(tDist, copula at tlags)]], h[,1])
} else {
res <- numeric(n)
for(t in tDist) {
tmpInd <- h[,2]==t
- tmpCop <- copula at spCopList[[match(t, copula at t.lags)]]
+ tmpCop <- copula at spCopList[[match(t, copula at tlags)]]
res[tmpInd] <- dduSpCopula(u[tmpInd,,drop=F], tmpCop, h[tmpInd,1])
}
}
@@ -159,16 +159,16 @@
n <- nrow(u)
tDist <- unique(h[,2])
- if(any(is.na(match(tDist,copula at t.lags))))
+ if(any(is.na(match(tDist,copula at tlags))))
stop("Prediction time(s) do(es) not math the modelled time slices.")
if (length(tDist)==1) {
- res <- ddvSpCopula(u, copula at spCopList[[match(tDist,copula at t.lags)]], h[,1])
+ res <- ddvSpCopula(u, copula at spCopList[[match(tDist,copula at tlags)]], h[,1])
} else {
res <- numeric(n)
for(t in tDist) {
tmpInd <- h[,2]==t
- tmpCop <- copula at spCopList[[match(t, copula at t.lags)]]
+ tmpCop <- copula at spCopList[[match(t, copula at tlags)]]
res[tmpInd] <- ddvSpCopula(u[tmpInd,,drop=F], tmpCop, h[tmpInd,1])
}
}
@@ -179,7 +179,53 @@
function(u, copula, ...) ddvStCopula(matrix(u,ncol=2), copula, ...))
setMethod("ddvCopula", signature("matrix","stCopula"), ddvStCopula)
-# dropping a sptio-temporal tree
+# log-likelihood by copula for all spatio-temporal lags
+
+
+loglikByCopulasStLags <- function(stBins, data, families = c(normalCopula(),
+ tCopula(),
+ claytonCopula(),
+ frankCopula(),
+ gumbelCopula()),
+ calcCor, lagSub=1:length(stBins$meanDists)) {
+ nTimeLags <- dim(stBins$lagCor)[1]
+ var <- attr(stBins, "variable")
+
+ 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(stBins$lags[[1]][lagSub], retrieveData, tempIndices=stBins$lags[[2]])
+
+ tmpBins <- list(meanDists=stBins$meanDists[lagSub])
+ attr(tmpBins, "variable") <- var
+
+ loglikTau <- list()
+ for(j in 1:nTimeLags) {
+ tmpLagData <- lapply(lagData, function(x) x[,c(2*j-1,2*j)])
+ tmpLagData <- lapply(tmpLagData, function(pairs) {
+ bool <- !is.na(pairs[,1]) & !is.na(pairs[,2])
+ pairs[bool,]
+ })
+
+ if(missing(calcCor))
+ res <- loglikByCopulasLags.static(tmpLagData, families)
+ else
+ res <- loglikByCopulasLags.dyn(tmpBins, tmpLagData, families,
+ function(h) calcCor(h, j, 1:nTimeLags))
+ loglikTau[[paste("loglik",j,sep="")]] <- res
+ }
+
+ return(loglikTau)
+}
+
+# dropping a spatio-temporal tree
dropStTree <- function (stNeigh, dataLocs, stCop) {
stopifnot(class(stNeigh) == "stNeighbourhood")
Modified: pkg/data/spCopDemo.RData
===================================================================
(Binary files differ)
Modified: pkg/demo/spCopula.R
===================================================================
--- pkg/demo/spCopula.R 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/demo/spCopula.R 2014-02-19 19:32:57 UTC (rev 127)
@@ -26,7 +26,7 @@
meuse$marZinc <- pMar(meuse$zinc)
## lag classes ##
-bins <- calcBins(meuse,var="marZinc",nbins=10,cutoff=800)
+bins <- calcBins(meuse, var="marZinc", nbins=10, cutoff=800)
## calculate parameters for Kendall's tau function ##
# either linear
@@ -38,7 +38,7 @@
curve(calcKTauPol,0, 1000, col="purple",add=TRUE)
## find best fitting copula per lag class
-loglikTau <- loglikByCopulasLags(bins, calcKTauPol,
+loglikTau <- loglikByCopulasLags(bins, meuse, calcKTauPol,
families=c(normalCopula(0), tCopula(0),
claytonCopula(0), frankCopula(1),
gumbelCopula(1), joeBiCopula(1.5),
@@ -57,11 +57,16 @@
spDepFun=calcKTauPol, unit="m")
## compare spatial copula loglik by lag:
+lagData <- lapply(bins$lags, function(x) {
+ as.matrix((cbind(meuse[x[,1], "marZinc"]@data,
+ meuse[x[,2], "marZinc"]@data)))
+ })
+
spLoglik <- NULL
for(i in 1:length(bins$lags)) { # i <- 7
cat("Lag",i,"\n")
spLoglik <- c(spLoglik,
- sum((dCopula(u=bins$lagData[[i]], spCop,log=T,
+ sum((dCopula(u=lagData[[i]], spCop,log=T,
h=bins$lags[[i]][,3]))))
}
@@ -71,7 +76,7 @@
legend(6, 50,c("Spatial Copula", "best copula per lag", "Gaussian Copula",
"number of pairs"),
pch=c(1,16,5,50), col=c("black", "green", "red"))
-text(x=(1:10+0.5),y=spLoglik,lapply(bins$lagData,length))
+text(x=(1:10+0.5), y=spLoglik, lapply(lagData,length))
##
# spatial vine
Modified: pkg/man/calcBins.Rd
===================================================================
--- pkg/man/calcBins.Rd 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/man/calcBins.Rd 2014-02-19 19:32:57 UTC (rev 127)
@@ -34,7 +34,7 @@
\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{tlags}{a vector indicating the time lags to be investigated}
}
}
\item{cor.method}{
Modified: pkg/man/condCovariate.Rd
===================================================================
--- pkg/man/condCovariate.Rd 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/man/condCovariate.Rd 2014-02-19 19:32:57 UTC (rev 127)
@@ -38,7 +38,7 @@
time[2:3])
stNeigh <- getStNeighbours(stData=stData, ST=stQuerry,
- spSize=3, t.lags=-(0:1),
+ spSize=3, tlags=-(0:1),
var="var", coVar="coVar", prediction=TRUE)
condCovariate(stNeigh, function(x) gumbelCopula(7))
Modified: pkg/man/condStCoVarVine.Rd
===================================================================
--- pkg/man/condStCoVarVine.Rd 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/man/condStCoVarVine.Rd 2014-02-19 19:32:57 UTC (rev 127)
@@ -59,7 +59,7 @@
unit="km")
stCop <- stCopula(components=list(spCopT0, spCopT1, spCopT2),
- t.lags=-(0:2))
+ tlags=-(0:2))
# only a constant copula ius used for the covariate
stCVVC <- stCoVarVineCopula(function(x) gumbelCopula(7), stCop, vineCopula(5L))
Modified: pkg/man/condStVine.Rd
===================================================================
--- pkg/man/condStVine.Rd 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/man/condStVine.Rd 2014-02-19 19:32:57 UTC (rev 127)
@@ -49,7 +49,7 @@
unit="km")
stCop <- stCopula(components=list(spCopT0, spCopT1),
- t.lags=-(0:1))
+ tlags=-(0:1))
stVineCop <- stVineCopula(stCop, vineCopula(4L))
Modified: pkg/man/fitCorFun.Rd
===================================================================
--- pkg/man/fitCorFun.Rd 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/man/fitCorFun.Rd 2014-02-19 19:32:57 UTC (rev 127)
@@ -5,44 +5,30 @@
Automated fitting of a correlation function to the correlogram
}
\description{
-Polynomials of different degrees can be fitted to the correlogram calculated using \code{\link{calcBins}}. This function will be used to adjust the copula parameter in the spatial/spatio-temporal copula.
+Polynomials of different degrees can be fitted to the spatial/spatio-temporal correlogram calculated using \code{\link{calcBins}}. This function will be used to adjust the copula parameter in the spatial/spatio-temporal copula.
}
\usage{
-fitCorFun(bins, degree = 3, cutoff = NA, bounds = c(0, 1), cor.method = NULL,
- weighted = FALSE)
+fitCorFun(bins, degree = 3, cutoff = NA, tlags, bounds = c(0, 1),
+ cor.method = NULL, weighted = FALSE)
}
\arguments{
- \item{bins}{
-Typically the output of \code{\link{calcBins}}. Any \code{data.frame} with a columns \code{lagCor} and \code{meanDists} in the first two columns will do.
+ \item{bins}{Typically the output of \code{\link{calcBins}}. Any \code{data.frame} with a columns \code{lagCor} and \code{meanDists} in the first two columns will do.}
+ \item{degree}{The degree of polynomial to be fitted - recycled if needed.}
+ \item{cutoff}{Maximal distance to which lags should be included in the polynomial fit.}
+ \item{tlags}{The temporal lags used for the genration of \code{bins}.}
+ \item{bounds}{Bounds of the correlation values. The default is set [0,1] not allowing any negative relationship but perfect positive dependence.}
+ \item{cor.method}{The output of \code{\link{calcBins}} has an attribute \code{cor.method}, in case this is not present the parameter \code{cor.method} will be used. In case the parameter \code{cor.method} is not \code{NULL} and the attribute \code{cor.method} is present, they will be compared.}
+ \item{weighted}{shall the residuals be weighted by the number of points in the lag class?}
}
- \item{degree}{
-The degree of polynomial to be fitted.
-}
- \item{cutoff}{
-Maximal distance to which lags should be included in the polynomial fit.
-}
- \item{bounds}{
-Bounds of the correlation values. The default is set [0,1] not allowing any negative relationship but perfect positive dependence.
-}
- \item{cor.method}{
-The output of \code{\link{calcBins}} has an attribute \code{cor.method}, in case this is not present the parameter \code{cor.method} will be used. In case the parameter \code{cor.method} is not \code{NULL} and the attribute \code{cor.method} is present, they will be compared.
-}
- \item{weighted}{
- shall the residulas be weighted by the number of points in the lag class?
- }
-}
-\value{
-Returns a function that provides correlation estimate for every separating distance.
-}
-\author{
-Benedikt Graeler
-}
-\seealso{
-See also \code{\link{calcBins}} and \code{\link{spCopula}}.
-}
+\value{Returns a one/two-place function that provides correlation estimates for every separating spatial/spatio-temporal distance.}
+
+\author{Benedikt Graeler}
+
+\seealso{See also \code{\link{calcBins}} and \code{\link{spCopula}}.}
+
\examples{
-# a simplified bins object (from demo(spcopula_estimation))
+# a simplified bins object (from demo(spcopula))
bins <- list(meanDists=c(64, 128, 203, 281, 361, 442, 522, 602, 681, 760),
lagCor=c(0.57, 0.49, 0.32, 0.29, 0.15, 0.14, 0.10, -0.00, 0.03, -0.01))
attr(bins,"cor.method") <- "kendall"
Modified: pkg/man/fitSpCopula.Rd
===================================================================
--- pkg/man/fitSpCopula.Rd 2014-02-17 15:46:15 UTC (rev 126)
+++ pkg/man/fitSpCopula.Rd 2014-02-19 19:32:57 UTC (rev 127)
@@ -4,31 +4,26 @@
\title{
Spatial Copula Fitting
}
+
\description{
A bivariate spatial copula is composed out of a set of bivariate copulas. These are combined using a convex linear combination with weights based on distances where for copulas with a 1-1 correspondence of Kendall's tau or Spearman's rho a dependence function providing measures of association based on distances might be used. This function estimates a spatial dependence function, evaluates the log-likelihood per family and lag class, selects the best fits and composes a spatial bivariate copula.
}
\usage{
-fitSpCopula(bins, cutoff = NA, families = c(normalCopula(), tCopula(),
+fitSpCopula(bins, data, cutoff = NA, families = c(normalCopula(), tCopula(),
claytonCopula(), frankCopula(), gumbelCopula()), ...)
}
\arguments{
- \item{bins}{
-the bins to be used, typically output from \code{\link{calcBins}}.
+ \item{bins}{the bins to be used, typically output from \code{\link{calcBins}}.}
+ \item{data}{the spatial dataset that ahs been used to generate \code{bins}.}
+ \item{cutoff}{The maximal distance to be used in the fit.}
+ \item{families}{The set of families to be investigated.}
+ \item{\dots}{Passed on to the function \code{\link{fitCorFun}}.}
}
- \item{cutoff}{
-The maximal distance to be used in the fit.
-}
- \item{families}{
-The set of families to be investigated.
-}
- \item{\dots}{
-Passed on to the function \code{\link{fitCorFun}}.
-}
-}
\value{
A \code{\linkS4class{spCopula}} object.
}
+
\references{
Graeler, B. & E. Pebesma (2011): The pair-copula construction for spatial data: a new approach to model spatial dependency.
Poster at: Spatial Statistics 2011 - Mapping global change. Enschede, The Netherlands, 23-25 March 2011.
@@ -42,9 +37,20 @@
Take a look at \code{\link{fitCorFun}}, \code{\link{loglikByCopulasLags}}, \code{\link{composeSpCopula}} and \code{\linkS4class{spCopula}} to gain more control over the fitting procedure.
}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/spcopula -r 127
More information about the spcopula-commits
mailing list