[spcopula-commits] r111 - pkg pkg/R pkg/man thoughts-Ben
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 11 18:05:37 CET 2013
Author: ben_graeler
Date: 2013-11-11 18:05:37 +0100 (Mon, 11 Nov 2013)
New Revision: 111
Modified:
pkg/DESCRIPTION
pkg/R/ClaytonGumbelCopula.R
pkg/R/joeBiCopula.R
pkg/R/spCopula.R
pkg/R/spVineCopula.R
pkg/R/spatio-temporalPreparation.R
pkg/R/stCopula.R
pkg/R/stVineCopula.R
pkg/man/condStVine.Rd
pkg/man/getStNeighbours.Rd
pkg/man/loglikByCopulasLags.Rd
pkg/man/spCopPredict.Rd
pkg/man/spVineCopula.Rd
pkg/man/stCopPredict.Rd
pkg/man/stCopula.Rd
pkg/man/stNeighbourhood.Rd
pkg/man/stVineCopula.Rd
thoughts-Ben/tawn3pCopula-stub.R
Log:
- moved VineCopula from Depends to Imports
- catched a couple of warnings in copula fitting procedures
- updated examples for spatio-temporal functions to be spatio-temporal
- chnaged the getStNeighbours function
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/DESCRIPTION 2013-11-11 17:05:37 UTC (rev 111)
@@ -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), R (>= 2.15.0)
-Imports: sp, spacetime (>= 1.0-2), methods
+Depends: copula (>= 0.999-7), R (>= 2.15.0)
+Imports: VineCopula (>= 1.1-2), sp, spacetime (>= 1.0-9), methods
URL: http://r-forge.r-project.org/projects/spcopula/
Collate:
Classes.R
Modified: pkg/R/ClaytonGumbelCopula.R
===================================================================
--- pkg/R/ClaytonGumbelCopula.R 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/R/ClaytonGumbelCopula.R 2013-11-11 17:05:37 UTC (rev 111)
@@ -71,7 +71,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("surClaytonCopula"),
function(copula, tau) {
- if(tau <= 0) warning("The survival Clayton copula can only represent positive dependence!")
+ if(tau <= 0)
+ return(NA)
linkVineCop.iTau(copula, max(1e-6,abs(tau)))
})
@@ -145,12 +146,12 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("r90ClaytonCopula"),
function(copula, tau) {
- if(tau >= 0) warning("The rotated Clayton copula can only represent negative dependence!")
+ if(tau >= 0)
+ return(NA)
linkVineCop.iTau(copula, min(-1e-6,-abs(tau)))
})
setMethod("tau",signature("r90ClaytonCopula"),linkVineCop.tau)
-
setMethod("tailIndex",signature("r90ClaytonCopula"),linkVineCop.tailIndex)
########################
@@ -205,7 +206,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("r270ClaytonCopula"),
function(copula, tau) {
- if(tau >= 0) warning("The rotated Clayton copula can only represent negative dependence!")
+ if(tau >= 0)
+ return(NA)
linkVineCop.iTau(copula, min(-1e-6,-abs(tau)))
})
@@ -286,7 +288,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("surGumbelCopula"),
function(copula, tau) {
- if(tau < 0) warning("The survival Gumbel copula can only represent non-negative dependence!")
+ if(tau < 0)
+ return(NA)
linkVineCop.iTau(copula, max(0,abs(tau)))
})
@@ -361,7 +364,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("r90GumbelCopula"),
function(copula, tau) {
- if(tau > 0) warning("The rotated Gumbel copula can only represent non-positive dependence!")
+ if(tau > 0)
+ return(NA)
linkVineCop.iTau(copula, min(0,-abs(tau)))
})
@@ -421,7 +425,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("r270GumbelCopula"),
function(copula, tau) {
- if(tau >= 0) warning("The rotated Gumbel copula can only represent negative dependence!")
+ if(tau >= 0)
+ return(NA)
linkVineCop.iTau(copula, min(-1e-6,-abs(tau)))
})
Modified: pkg/R/joeBiCopula.R
===================================================================
--- pkg/R/joeBiCopula.R 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/R/joeBiCopula.R 2013-11-11 17:05:37 UTC (rev 111)
@@ -68,7 +68,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("joeBiCopula"),
function(copula, tau) {
- if(tau <= 0) warning("The Joe copula can only represent positive dependence!")
+ if(tau <= 0)
+ return(NA)
linkVineCop.iTau(copula, max(1e-6,abs(tau)))
})
@@ -147,7 +148,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("surJoeBiCopula"),
function(copula, tau) {
- if(tau <= 0) warning("The survival Joe copula can only represent positive dependence!")
+ if(tau <= 0)
+ return(NA)
linkVineCop.iTau(copula, max(1e-6,abs(tau)))
})
@@ -224,7 +226,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("r90JoeBiCopula"),
function(copula, tau) {
- if(tau >= 0) warning("The rotated Joe copula can only represent negative dependence!")
+ if(tau >= 0)
+ return(NA)
linkVineCop.iTau(copula, min(-1e-6,-abs(tau)))
})
@@ -286,7 +289,8 @@
## Kendalls tau to parameter conversion
setMethod("iTau", signature("r270JoeBiCopula"),
function(copula, tau) {
- if(tau >= 0) warning("The rotated Joe copula can only represent negative dependence!")
+ if(tau >= 0)
+ return(NA)
linkVineCop.iTau(copula, min(-1e-6,-abs(tau)))
})
Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/R/spCopula.R 2013-11-11 17:05:37 UTC (rev 111)
@@ -487,7 +487,7 @@
loglik <- NULL
copulas <- list()
for (cop in families) {
- print(cop)
+ cat(cop at fullname,"\n")
tmploglik <- NULL
tmpCop <- list()
for(i in 1:length(bins$meanDists)) {
@@ -509,12 +509,16 @@
stop(paste(calcCor(NULL), "is not yet supported.")))
} else {
param <- moa(cop, bins$meanDists[i])
- cop at parameters[1:length(param)] <- param
+ if(!is.na(param))
+ cop at parameters[1:length(param)] <- param
}
}
}
- tmploglik <- c(tmploglik, sum(dCopula(bins$lagData[[i]], cop, log=T)))
+ if(is.na(param))
+ tmploglik <- c(tmploglik, NA)
+ else
+ tmploglik <- c(tmploglik, sum(dCopula(bins$lagData[[i]], cop, log=T)))
tmpCop <- append(tmpCop, cop)
}
loglik <- cbind(loglik, tmploglik)
@@ -531,7 +535,7 @@
fits <-lapply(families,
function(cop) {
- print(cop)
+ cat(cop at fullname,"\n")
lapply(bins$lagData,
function(x) {
tryCatch(fitCopula(cop, x, estimate.variance = FALSE),
Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/R/spVineCopula.R 2013-11-11 17:05:37 UTC (rev 111)
@@ -204,7 +204,7 @@
}
ePred <- integrate(condExp,0,1,subdivisions=10000L,stop.on.error=stop.on.error, ...)
- if(ePred$abs.error > 0.01)
+ if(ePred$abs.error > 0.05)
warning("Numerical integration in predExpectation performed at a level of absolute error of only ",
ePred$abs.error, " for location ",i,".")
predMean <- c(predMean, ePred$value)
Modified: pkg/R/spatio-temporalPreparation.R
===================================================================
--- pkg/R/spatio-temporalPreparation.R 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/R/spatio-temporalPreparation.R 2013-11-11 17:05:37 UTC (rev 111)
@@ -19,8 +19,8 @@
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,
+ new("stNeighbourhood", data=data, distances=distances, locations=ST,
+ dataLocs=STxDF, index=index, prediction=prediction, var=var,
sp=as(STxDF at sp, "Spatial"), time=STxDF at time[1],
endTime=STxDF at endTime[length(STxDF at endTime)])
}
@@ -48,36 +48,52 @@
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
+
+ if(!prediction) {
+ if(is.na(timeSteps)) {
+ timeSteps <- length(stData at time)+timeSpan
+ reSample <- function() (1-timeSpan):length(stData at time)
+ } else {
+ reSample <- function() sort(sample((1-timeSpan):length(stData at time), timeSteps))
+ }
+ nLocs <- length(ST at sp)*timeSteps
nghbrs <- getNeighbours(stData[,1], var=var, size=spSize, min.dist=min.dist)
+ } else {
+ nLocs <- length(ST)
+ nghbrs <- getNeighbours(stData[,1], ST at sp, var, spSize, prediction, min.dist)
+ timeNghbrs <- sapply(index(ST at time), function(x) which(x == index(stData at time)))
+ reSample <- function() timeNghbrs
+ timeSteps <- length(stData at time)+timeSpan
+ }
- stNeighData <- NULL
+ stNeighData <- matrix(NA, nLocs, (spSize-1)*length(t.lags)+1)
stDists <- array(NA,c(nLocs,(spSize-1)*length(t.lags),2))
stInd <- array(NA,c(nLocs,(spSize-1)*length(t.lags),2))
+
+ nTimeInst <- length(reSample())
+
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[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,][-1],
- tmpInst+t,var]@data[[1]],
- ncol=spSize-1, byrow=T))
- tmpInd <- cbind(tmpInd, matrix(rep(tmpInst+t,spSize-1),ncol=spSize-1))
+ timeInst <- reSample() # draw random time steps for each neighbourhood
+ stNeighData[(i-1)*timeSteps+(1:timeSteps),
+ 1:spSize] <- matrix(stData[nghbrs at index[i,], timeInst,
+ var, drop=F]@data[[1]],
+ ncol=spSize, byrow=T) # retrieve the top level data
+ tmpInd <- matrix(rep(timeInst, spSize-1), ncol=spSize-1)
+ for(j in 2:length(t.lags)) {
+ t <- t.lags[j]
+ stNeighData[(i-1)*timeSteps+(1:timeSteps),
+ (j-1)*(spSize-1)+2:(spSize)] <- matrix(stData[nghbrs at index[i,][-1],
+ timeInst+t,
+ var, drop=F]@data[[1]],
+ ncol=spSize-1, byrow=T)
+ tmpInd <- cbind(tmpInd, matrix(rep(timeInst+t,spSize-1),ncol=spSize-1))
}
- stNeighData <- rbind(stNeighData, tmpData) # bind data row-wise
+
stDists[(i-1)*timeSteps+1:timeSteps,,1] <- matrix(rep(nghbrs at distances[i,],
timeSteps*length(t.lags)),
byrow=T, ncol=length(t.lags)*(spSize-1)) # store sp distances
@@ -90,10 +106,12 @@
stInd[(i-1)*timeSteps+1:timeSteps,,2] <- tmpInd
}
- if (prediction)
+ if (prediction) {
dataLocs <- stData
- else
+ stNeighData <- stNeighData[,-1]
+ } else {
dataLocs <- NULL
+ }
return(stNeighbourhood(as.data.frame(stNeighData), stDists, stData, ST,
stInd, prediction, var))
}
Modified: pkg/R/stCopula.R
===================================================================
--- pkg/R/stCopula.R 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/R/stCopula.R 2013-11-11 17:05:37 UTC (rev 111)
@@ -5,18 +5,25 @@
## constructor ##
#################
-stCopula <- function(components, distances, t.lags, stDepFun, unit="m", t.res="day") {
- 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)){
- spCopList <- append(spCopList, getSpCop(components[[i]], distances[[i]], i))
- }
+stCopula <- function(components, t.lags, distances=NA, stDepFun, unit="m", t.res="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))
+ spCopList <- components
} else {
- for(i in 1:length(t.lags)){
- spCopList <- append(spCopList, spCopula(components[[i]], distances[[i]], unit=unit))
+ 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)){
+ spCopList <- append(spCopList, getSpCop(components[[i]], distances[[i]], i))
+ }
+ } else {
+ for(i in 1:length(t.lags)){
+ spCopList <- append(spCopList, spCopula(components[[i]], distances[[i]], unit=unit))
+ }
}
}
Modified: pkg/R/stVineCopula.R
===================================================================
--- pkg/R/stVineCopula.R 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/R/stVineCopula.R 2013-11-11 17:05:37 UTC (rev 111)
@@ -173,7 +173,7 @@
predQuantile <- NULL
for(i in 1:nrow(predNeigh at data)) { # i <-1
cat("[Predicting location ",i,".]\n", sep="")
- condSecVine <- condStVine(as.numeric(predNeigh at data[i,]), dists[i,], stVine)
+ condSecVine <- condStVine(as.numeric(predNeigh at data[i,]), dists[i,,,drop=F], stVine)
xVals <- attr(condSecVine,"xVals")
density <- condSecVine(xVals)
Modified: pkg/man/condStVine.Rd
===================================================================
--- pkg/man/condStVine.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/condStVine.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -35,27 +35,29 @@
\code{\linkS4class{stVineCopula}}, \code{\link{condSpVine}}
}
\examples{
-## the spatial version
-data(spCopDemo)
+# a spatio-temporal C-vine copula (with independent copulas in the upper vine)
-calcKTauPol <- fitCorFun(bins, degree=3)
+spCopT0 <- spCopula(components=list(claytonCopula(8), claytonCopula(4),
+ claytonCopula(2), claytonCopula(1),
+ claytonCopula(0.5), indepCopula()),
+ distances=c(100,200,300,400,500,600),
+ unit="km")
+spCopT1 <- spCopula(components=list(claytonCopula(4), claytonCopula(2),
+ claytonCopula(1), claytonCopula(0.5),
+ indepCopula()),
+ distances=c(100,200,300,400,500),
+ unit="km")
-spCop <- spCopula(components=list(normalCopula(0), tCopula(0, dispstr = "un"),
- frankCopula(1), normalCopula(0), claytonCopula(0),
- claytonCopula(0), claytonCopula(0), claytonCopula(0),
- claytonCopula(0), indepCopula()),
- distances=bins$meanDists,
- spDepFun=calcKTauPol, unit="m")
+stCop <- stCopula(components=list(spCopT0, spCopT1),
+ t.lags=-(0:1))
-spVineCop <- spVineCopula(spCop, vineCopula(4L))
+stVineCop <- stVineCopula(stCop, vineCopula(4L))
-dists <- list(c(473, 124, 116, 649))
+dists <- array(c(150, 250, 150, 250,0,0,-1,-1),dim=c(1,4,2))
condVar <- c(0.29, 0.55, 0.05, 0.41)
-condDensity <- condSpVine(condVar,dists,spVineCop)
+condDensity <- condStVine(condVar,dists,stVineCop)
curve(condDensity)
-mtext(paste("Dists:",paste(round(dists[[1]],0),collapse=", ")),line=0)
-mtext(paste("Cond.:",paste(round(condVar,2),collapse=", ")),line=1)
}
\keyword{ distribution }
\ No newline at end of file
Modified: pkg/man/getStNeighbours.Rd
===================================================================
--- pkg/man/getStNeighbours.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/getStNeighbours.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -42,11 +42,18 @@
See \code{\link{stNeighbourhood}} for the native constructor of a \code{\linkS4class{stNeighbourhood}} class. The pure spatial version can be found at \code{\link{getNeighbours}}.
}
\examples{
-## the spatial version:
library(sp)
-spdf <- data.frame(x=c(112,154,212,289,345),y=c(124,198,85,168,346),measure=rlnorm(5))
-coordinates(spdf) <- ~x+y
+library(spacetime)
-getNeighbours(spdf,size=4)
+sp <- SpatialPoints(matrix(c(181000,181100,333500,333600),2))
+time <- Sys.time()+60*60*24*c(0,1,2)
+data <- data.frame(var1=runif(6))
+
+stData <- STFDF(sp, time, data)
+stQuerry <- STF(SpatialPoints(matrix(c(181000,181200,333600,333600),2)),
+ time[2:3])
+
+getStNeighbours(stData=stData, ST=stQuerry, prediction=TRUE, spSize=3,
+ t.lags=-(0:1))
}
\keyword{ spatial }
\ No newline at end of file
Modified: pkg/man/loglikByCopulasLags.Rd
===================================================================
--- pkg/man/loglikByCopulasLags.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/loglikByCopulasLags.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -38,7 +38,7 @@
calcKTauPol <- fitCorFun(bins, degree=3)
-loglikTau <- loglikByCopulasLags(bins, calcKTauPol)
+loglikTau <- loglikByCopulasLags(bins, calcCor=calcKTauPol)
loglikTau$loglik
}
Modified: pkg/man/spCopPredict.Rd
===================================================================
--- pkg/man/spCopPredict.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/spCopPredict.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -58,7 +58,7 @@
meuse$rtZinc <- rank(meuse$zinc)/(length(meuse)+1)
-predMeuseNeigh <- getNeighbours(meuse[1:4,], meuse.grid[c(9:12,15:19,24:28,34:38),],
+predMeuseNeigh <- getNeighbours(meuse[1:4,], meuse.grid[c(9:12,16:19,25:28),],
"rtZinc", 5L, TRUE, -1)
qMar <- function(x) {
@@ -67,10 +67,12 @@
predMedian <- spCopPredict(predMeuseNeigh, spVineCop, list(q=qMar), "quantile", p=0.5)
+\dontrun{
spplot(predMedian,"quantile.0.5",
sp.layout=list("sp.points", meuse, pch = 19, col = "red"),
col.regions=bpy.colors())
}
+}
\keyword{ distribution }
\keyword{ prediction }
Modified: pkg/man/spVineCopula.Rd
===================================================================
--- pkg/man/spVineCopula.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/spVineCopula.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -35,6 +35,7 @@
distances=bins$meanDists,
spDepFun=calcKTauPol, unit="m")
+library(VineCopula)
RVM <- RVineMatrix(matrix(c(1,0,0,2,2,0,3,3,3),3,byrow=TRUE))
spVineCopula(spCop, vineCopula(RVM))
}
Modified: pkg/man/stCopPredict.Rd
===================================================================
--- pkg/man/stCopPredict.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/stCopPredict.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -32,44 +32,39 @@
\code{\link{condStVine}} and \code{\link{spCopPredict}} for the spatial version.
}
\examples{
-## the spatial version
-data(spCopDemo)
+library(sp)
+library(spacetime)
-calcKTauPol <- fitCorFun(bins, degree=3)
+spCopT0 <- spCopula(components=list(claytonCopula(8), claytonCopula(4),
+ claytonCopula(2), claytonCopula(1),
+ claytonCopula(0.5), indepCopula()),
+ distances=c(100,200,300,400,500,600),
+ unit="km")
+spCopT1 <- spCopula(components=list(claytonCopula(4), claytonCopula(2),
+ claytonCopula(1), claytonCopula(0.5),
+ indepCopula()),
+ distances=c(100,200,300,400,500),
+ unit="km")
-spCop <- spCopula(components=list(normalCopula(0), tCopula(0, dispstr = "un"),
- frankCopula(1), normalCopula(0), claytonCopula(0),
- claytonCopula(0), claytonCopula(0), claytonCopula(0),
- claytonCopula(0), indepCopula()),
- distances=bins$meanDists,
- spDepFun=calcKTauPol, unit="m")
+stCop <- stCopula(components=list(spCopT0, spCopT1),
+ t.lags=-(0:1))
-spVineCop <- spVineCopula(spCop, vineCopula(4L))
+stVineCop <- stVineCopula(stCop, vineCopula(4L))
-library(sp)
-data(meuse.grid)
-coordinates(meuse.grid) <- ~x+y
-gridded(meuse.grid) <- TRUE
+sp <- SpatialPoints(matrix(c(181000,181100,333500,333600),2))
+time <- Sys.time()+60*60*24*c(0,1,2)
+data <- data.frame(var1=runif(6))
-data(meuse)
-coordinates(meuse) <- ~x+y
+stData <- STFDF(sp, time, data)
+stQuerry <- STF(SpatialPoints(matrix(c(181000,181200,333600,333600),2)),
+ time[2:3])
-meuse$rtZinc <- rank(meuse$zinc)/(length(meuse)+1)
+stNeigh <- getStNeighbours(stData=stData, ST=stQuerry, prediction=TRUE, spSize=3,
+ t.lags=-(0:1))
-predMeuseNeigh <- getNeighbours(meuse[1:4,], meuse.grid[c(9:12,15:19,24:28,34:38),],
- "rtZinc", 5L, TRUE, -1)
-
-qMar <- function(x) {
- qlnorm(x,mean(log(meuse$zinc)),sd(log(meuse$zinc)))
+stCopPredict(stNeigh, stVineCop, list(q=qunif), "quantile", 0.5)
}
-predMedian <- spCopPredict(predMeuseNeigh, spVineCop, list(q=qMar), "quantile", p=0.5)
-
-spplot(predMedian,"quantile.0.5",
- sp.layout=list("sp.points", meuse, pch = 19, col = "red"),
- col.regions=bpy.colors())
-}
-
\keyword{ distribution }
\keyword{ prediction }
\keyword{ spatial }
Modified: pkg/man/stCopula.Rd
===================================================================
--- pkg/man/stCopula.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/stCopula.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -7,25 +7,17 @@
Constructor of a bivariate spatio-temporal copula \code{\linkS4class{stCopula}}.
}
\usage{
-stCopula(components, distances, t.lags, stDepFun, unit="m", t.res="day")
+stCopula(components, t.lags, distances=NA, stDepFun, unit="m", t.res="day")
}
\arguments{
- \item{components}{
-A list of bivariate spatial copulas (\code{\linkS4class{spCopula}}) to be used at each temporal lag.
+ \item{components}{A list of bivariate spatial copulas (\code{\linkS4class{spCopula}}) to be used at each temporal lag. Or a list of with lists of the spatial components per temporal lag together with the argument \code{distances}.}
+ \item{t.lags}{The temporal lags used in the spatio-temporal copula.}
+ \item{distances}{This and the follwoing 2 argumewnts are only neccessary when the provided \code{components} argument is not yet a list of \code{\linkS4class{spCopula}}s: A vector of the mean distances of the spatial lag classes.}
+ \item{stDepFun}{A list of spatial dependence functions; one per temporal lag. This argument is only needed when components is not yet a list of \code{\linkS4class{spCopula}}s.}
+ \item{unit}{The spatial unit, default: m (meters). This argument is only needed when components is not yet a list of \code{\linkS4class{spCopula}}s.}
+ \item{t.res}{The temporal resolution, default: day}
}
- \item{distances}{
-A vector of the mean distances of the spatial lag classes.
-}
-\item{t.lags}{The temporal lags used in the spatio-temporal copula.}
- \item{stDepFun}{
-A list of spatial dependence functions; one per temporal lag.
-}
- \item{unit}{
-The spatial unit, default: m (meters)
-}
-\item{t.res}{The temporal resolution, default: day}
-}
\value{
An instance of the spatio-temporal Copula class \code{\linkS4class{stCopula}}.
}
@@ -37,16 +29,23 @@
\code{\link{spCopula}}
}
\examples{
-data(spCopDemo)
+spCopT0 <- spCopula(components=list(claytonCopula(8), claytonCopula(4),
+ claytonCopula(2), claytonCopula(1),
+ claytonCopula(0.5), indepCopula()),
+ distances=c(100,200,300,400,500,600),
+ unit="km")
+spCopT1 <- spCopula(components=list(claytonCopula(4), claytonCopula(2),
+ claytonCopula(1), claytonCopula(0.5),
+ indepCopula()),
+ distances=c(100,200,300,400,500),
+ unit="km")
+spCopT2 <- spCopula(components=list(claytonCopula(2), claytonCopula(1),
+ claytonCopula(0.5), indepCopula()),
+ distances=c(100,200,300,400),
+ unit="km")
-calcKTauPol <- fitCorFun(bins, degree=3)
-
-spCop <- spCopula(components=list(normalCopula(0), tCopula(0, dispstr = "un"),
- frankCopula(1), normalCopula(0), claytonCopula(0),
- claytonCopula(0), claytonCopula(0), claytonCopula(0),
- claytonCopula(0), indepCopula()),
- distances=bins$meanDists,
- spDepFun=calcKTauPol, unit="m")
+stCop <- stCopula(components=list(spCopT0, spCopT1, spCopT2),
+ t.lags=-(0:2))
}
\keyword{spcopula}
\keyword{copula}
Modified: pkg/man/stNeighbourhood.Rd
===================================================================
--- pkg/man/stNeighbourhood.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/stNeighbourhood.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -25,16 +25,19 @@
\code{\linkS4class{stNeighbourhood}}, \code{\link{getStNeighbours}}
}
\examples{
-## the spatial version
library(sp)
-spdf <- data.frame(x=c(112,154,212,289),y=c(124,198,85,168),measure=rlnorm(4))
-coordinates(spdf) <- ~x+y
+library(spacetime)
-neigh <- getNeighbours(spdf,size=4)
-neigh
+sp <- SpatialPoints(matrix(c(181000,181100,333500,333600),2))
+time <- Sys.time()+60*60*24*c(0,1,2)
+data <- data.frame(var1=runif(6))
-# rebuilding neigh
-neighbourhood(neigh at data, neigh at distances, neigh at index, spdf, NULL, neigh at prediction, neigh at var)
+stData <- STFDF(sp, time, data)
+stQuerry <- STF(SpatialPoints(matrix(c(181000,181200,333600,333600),2)),
+ time[2:3])
+
+getStNeighbours(stData=stData, ST=stQuerry, prediction=TRUE, spSize=3,
+ t.lags=-(0:1))
}
\keyword{spatio-temporal}
\ No newline at end of file
Modified: pkg/man/stVineCopula.Rd
===================================================================
--- pkg/man/stVineCopula.Rd 2013-10-24 11:37:56 UTC (rev 110)
+++ pkg/man/stVineCopula.Rd 2013-11-11 17:05:37 UTC (rev 111)
@@ -20,21 +20,28 @@
Benedikt Graeler
}
\examples{
-## the spatial version
-# a spatial C-vine copula (with independent dummy copulas in the upper vine)
-data(spCopDemo)
+# a spatio-temporal C-vine copula (with independent copulas in the upper vine)
-calcKTauPol <- fitCorFun(bins, degree=3)
+spCopT0 <- spCopula(components=list(claytonCopula(8), claytonCopula(4),
+ claytonCopula(2), claytonCopula(1),
+ claytonCopula(0.5), indepCopula()),
+ distances=c(100,200,300,400,500,600),
+ unit="km")
+spCopT1 <- spCopula(components=list(claytonCopula(4), claytonCopula(2),
+ claytonCopula(1), claytonCopula(0.5),
+ indepCopula()),
+ distances=c(100,200,300,400,500),
+ unit="km")
+spCopT2 <- spCopula(components=list(claytonCopula(2), claytonCopula(1),
+ claytonCopula(0.5), indepCopula()),
+ distances=c(100,200,300,400),
+ unit="km")
-spCop <- spCopula(components=list(normalCopula(0), tCopula(0, dispstr = "un"),
- frankCopula(1), normalCopula(0), claytonCopula(0),
- claytonCopula(0), claytonCopula(0), claytonCopula(0),
- claytonCopula(0), indepCopula()),
- distances=bins$meanDists,
- spDepFun=calcKTauPol, unit="m")
+stCop <- stCopula(components=list(spCopT0, spCopT1, spCopT2),
+ t.lags=-(0:2))
-RVM <- RVineMatrix(matrix(c(1,0,0,2,2,0,3,3,3),3,byrow=TRUE))
-spVineCopula(spCop, vineCopula(RVM))
+library(VineCopula)
+stVineCopula(stCop, vineCopula(9L))
}
\keyword{ mulitvariate }
\keyword{ distribution }
\ No newline at end of file
Modified: thoughts-Ben/tawn3pCopula-stub.R
===================================================================
--- thoughts-Ben/tawn3pCopula-stub.R 2013-10-24 11:37:56 UTC (rev 110)
+++ thoughts-Ben/tawn3pCopula-stub.R 2013-11-11 17:05:37 UTC (rev 111)
@@ -9,11 +9,35 @@
beta <- param[2]
theta <- param[3]
(1-beta)*(t) + (1-alpha)*(1-t) + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
-# 1-beta + (beta-alpha)*t + ((alpha*t)^theta + (beta*(1-t))^theta)^(1/theta)
+
}
-curve(Atawn3p)
+ATawn <- function(copula, w) {
+ Atawn3p(w,copula at parameters)
+}
+setMethod("A",signature("tawn3pCopula"),ATawn)
+
+dAduTawn <- function(copula, w) {
+ alpha <- copula at parameters[1]
+ beta <- copula at parameters[2]
+ theta <- copula at parameters[3]
+
+ # 1st derivative
+ p1 <- (alpha*(alpha*(-(w-1)))^(theta-1)-beta*(beta*w)^(theta-1))
+ p2 <- ((alpha*(-(w-1)))^theta+(beta*w)^theta)^(1/theta-1)
+
+ # 2nd derivative
+ p3 <- (alpha*(-(w-1)))^(theta-2)
+ p4 <- (beta*w)^(theta-2)
+ p5 <- ((alpha*(-(w-1)))^theta+(beta*w)^theta)^(1/theta-2)
+
+ data.frame(der1=alpha-beta-p1*p2,
+ der2=alpha^2*beta^2*(theta-1)*p3*p4*p5)
+}
+
+setMethod("dAdu",signature("tawn3pCopula"),dAduTawn)
+
tawn3pCopula <- function (param = c(0.5, 0.5, 2)) {
# A(t) = (1-beta)*(1-t) + (1-alpha)*t + ((alpha*(1-t))^theta+(beta*t)^theta)^(1/theta)
# C(u1,u2) = exp(log(u1*u2) * A(log(u2)/log(u1*u2)))
@@ -70,54 +94,57 @@
persp(tawn3pCopula(c(0.25, 0.75, 2)), dCopula)
persp(tawn3pCopula(c(0.5, 1, 20)), pCopula)
-plot(1-rtTriples[,c(1,3)])
-
-dCopula(c(0.15,0.95), tawn3pCopula(c(0.5, 0.5, 2)))
-
-
-tawnFit <- fitCopula(tawn3pCopula(), 1-as.matrix(rtTriples[,c(1,3)]), hideWarnings=F,estimate.variance=F,
+tawnFit <- fitCopula(tawn3pCopula(c(0.25, 0.75, 2)), 1-as.matrix(rtTriples[,c(1,3)]), hideWarnings=F,estimate.variance=F,
start=c(0.9, 1, 8), method="mpl", lower=c(0,0,1), upper=c(1, 1, 10),
optim.method="L-BFGS-B",)
tawnFit at loglik # 742
-tawnFit at copula
+tawnCop <- tawnFit at copula
-fitCopula(gumbelCopula(5),1-as.matrix(rtTriples[,c(1,3)]))@loglik # 723
+par(mfrow=c(2,2))
+plot(rCopula(500,tawn3pCopula(c(tawnCop$par2, 1, tawnCop$par))),asp=1)
+plot(rCopula(500,tawnFit at copula),asp=1)
+plot(rCopula(500,cdfAFunCopula(aGevPar)),asp=1)
+plot(1-rtTriples[,c(3,1)], asp=1)
-dLeaf <- dCopula(as.matrix(rtTriples[,c(1,3)]), spcopula:::leafCopula())
-sum(log(dLeaf[dLeaf>0]))
+# fitCopula(gumbelCopula(5),1-as.matrix(rtTriples[,c(1,3)]))@loglik # 723
+#
+# dLeaf <- dCopula(as.matrix(rtTriples[,c(1,3)]), spcopula:::leafCopula())
+# sum(log(dLeaf[dLeaf>0]))
+#
+# persp(tawnFit at copula, dCopula)
+# contour(tawnFit at copula, dCopula, levels=c(0,0.5,1,2,4,8,100), asp=1)
+#
+# sum(dCopula(as.matrix(rtTriples[,c(1,3)]), cop13, log=T))
+# sum(dCopula(1-as.matrix(rtTriples[,c(1,3)]), tawnFit at copula, log=T))
-persp(tawnFit at copula, dCopula)
-contour(tawnFit at copula, dCopula, levels=c(0,0.5,1,2,4,8,100), asp=1)
-
-sum(dCopula(as.matrix(rtTriples[,c(1,3)]), cop13, log=T))
-sum(dCopula(1-as.matrix(rtTriples[,c(1,3)]), tawnFit at copula, log=T))
-
+par(mfrow=c(1,1))
plot(1-as.matrix(rtTriples[,c(1,3)]),asp=1,cex=0.5)
curve(x^(tawnFit at copula@parameters[1]),add=T, col="red")
abline(0,1,col="grey")
-copula:::fitCopula.ml
-
-tawn3pCopula()
-
-persp(tawn3pCopula(),dCopula)
-
-
###
-# h(t)
+# h(t), TUM thesis eq. (4.11)
+library(evd)
hist(log(1-rtTriples[,3])/log((1-rtTriples[,1])*(1-rtTriples[,3])),n=20,
- xlim=c(0,1), freq=F, add=T, col="blue")
+ xlim=c(0,1), freq=F, add=F, col="blue")
tSmpl <- log(1-rtTriples[,3])/log((1-rtTriples[,1])*(1-rtTriples[,3]))
- ((1-rtTriples[,1])-(1-rtTriples[,3]))*(0.5)+0.5
+# ((1-rtTriples[,1])-(1-rtTriples[,3]))*(0.5)+0.5
dlogNorm <- function(x) dlnorm(x, mean(log(tSmpl)), sd(log(tSmpl)))
+
aGevPar <- fgev(tSmpl)$estimate
-dGev <- function(x) dgev(x, 0.46938, 0.05057, 0.01720)
-dGamma <- function(x) dgamma(x, 60.40556, 120.92807)
-
+dGev <- function(x) dgev(x, aGevPar[1], aGevPar[2], aGevPar[3])
+
+optFun <- function(param) {
+ -sum(log(dgamma(tSmpl,param[1],param[2])))
+}
+aGammaPar <- optim(c(1,0.5),optFun)$par
+dGamma <- function(x) dgamma(x, aGammaPar[1], aGammaPar[2])
+
+par(mfrow=c(1,1))
hist(tSmpl, freq=F, xlim=c(0,1), n=20, add=F)
curve(dlogNorm, add=T)
curve(dGev, add=T, col="red")
@@ -131,10 +158,13 @@
Afit <- function(t) {
res <- t
res[res == 0] <- 1
- intFun <- function(z) (pgev(z, 0.46938, 0.05057, 0.01720)-z)/(z-z^2)
+
+ intFun <- function(z) (pgev(z, aGevPar[1], aGevPar[2], aGevPar[3])-z)/(z-z^2)
+
for(i in which(res != 1)) {
res[i] <- exp(integrate(intFun,0,t[i])$value)
}
+
return(res)
}
@@ -143,20 +173,143 @@
abline(1, -1, col="grey")
curve(Atawn3p,col="red",add=T)
+## understanding why some cdfs produce a convex A and some do not:
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/spcopula -r 111
More information about the spcopula-commits
mailing list