[spcopula-commits] r90 - / pkg pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 4 19:50:25 CEST 2013
Author: ben_graeler
Date: 2013-04-04 19:50:25 +0200 (Thu, 04 Apr 2013)
New Revision: 90
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/Classes.R
pkg/R/spVineCopula.R
pkg/R/spatialPreparation.R
pkg/R/vineCopulas.R
spcopula_0.1-1.tar.gz
Log:
- added spCopPredict functions
- cleaning
- getNeighbours can be used for prediction as well, twice as fast as before
- all fitCopula retun now a fitCopula object
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-03-27 19:13:09 UTC (rev 89)
+++ pkg/DESCRIPTION 2013-04-04 17:50:25 UTC (rev 90)
@@ -2,7 +2,7 @@
Type: Package
Title: copula driven spatial analysis
Version: 0.1-1
-Date: 2013-02-20
+Date: 2013-04-04
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.
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-03-27 19:13:09 UTC (rev 89)
+++ pkg/NAMESPACE 2013-04-04 17:50:25 UTC (rev 90)
@@ -19,6 +19,7 @@
export(dduCopula,ddvCopula)
export(invdduCopula, invddvCopula)
export(qCopula_u)
+export(condSpVine,spCopPredict)
# tweaks
# export(setSizeLim)
Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R 2013-03-27 19:13:09 UTC (rev 89)
+++ pkg/R/Classes.R 2013-04-04 17:50:25 UTC (rev 90)
@@ -147,9 +147,12 @@
else return (TRUE)
}
+setOldClass("RVineMatrix")
+
setClass("vineCopula",
representation = representation(copulas="list", dimension="integer",
- RVM="list"),
+ RVM="RVineMatrix"),
+ prototype = prototype(RVM=structure(list(),class="RVineMatrix")),
validity = validVineCopula,
contains = list("copula")
)
@@ -188,18 +191,20 @@
validNeighbourhood <- function(object) {
sizeN <- ncol(object at distances)+1
nVars <- length(object at varNames)
- if (sizeN > sizeLim) return("The limting size of the neighbourhood is exceeded. Increase the constant sizeLim if needed.")
if (nrow(object at data) != nrow(object at distances)) return("Data and distances have unequal number of rows.")
- if (ncol(object at data) %% sizeN != 0) return("Data and distances have non matching number of columns.")
-# if (nrow(object at data) != nrow(object at coords) ) return("Data and sp at coordinates have unequal number of rows.")
+ if (ncol(object at data) %% (sizeN-object at prediction) != 0) return("Data and distances have non matching number of columns.")
if (nrow(object at data) != nrow(object at index)) return("Data and index have unequal number of rows.")
- if (sizeN != ncol(object at index)) return("Data and index have unequal number of columns.")
- if (ncol(object at data) != sizeN * nVars) return(paste("Number of columns in data does not equal the product of the neighbourhood's size (",sizeN,") with number of variables (",nVars,").",sep=""))
+ if (ncol(object at distances) != ncol(object at index)) return("Data and index have unequal number of columns.")
+ if (ncol(object at data) != (sizeN-object at prediction) * nVars) return(paste("Number of columns in data does not equal the product of the neighbourhood's size (",sizeN,") with number of variables (",nVars,").",sep=""))
else return(TRUE)
}
setClass("neighbourhood",
- representation = representation(data = "data.frame", distances="matrix", "SpatialPoints", index="matrix", varNames="character"),
- validity = validNeighbourhood,
- contains = list("SpatialPoints"))
+ representation = representation(data = "data.frame",
+ distances="matrix",
+ index="matrix",
+ locations="Spatial",
+ varNames="character",
+ prediction="logical"),
+ validity = validNeighbourhood, contains = list("Spatial"))
Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R 2013-03-27 19:13:09 UTC (rev 89)
+++ pkg/R/spVineCopula.R 2013-04-04 17:50:25 UTC (rev 90)
@@ -47,9 +47,14 @@
dspVine(matrix(u,ncol=copula at dimension), copula at spCop, copula at vineCop, log=log, ...)
})
+setMethod("dCopula",signature=signature("data.frame","spVineCopula"),
+ function(u, copula, log, ...) {
+ dspVine(as.matrix(u), copula at spCop, copula at vineCop, log=log, ...)
+ })
+
# fiiting the spatial vine for a given spatial copula
-fitSpVine <- function(copula, data) {
+fitSpVine <- function(copula, data, method) {
stopifnot(class(data)=="neighbourhood")
stopifnot(copula at dimension == ncol(data at data))
@@ -60,9 +65,83 @@
copula=copula at spCop, h=data at distances[,i]))
}
- vineCop <- fitCopula(copula at vineCop, secLevel)
+ vineCopFit <- fitCopula(copula at vineCop, secLevel, method)
- return(spVineCopula(copula at spCop, vineCop))
+ spVineCop <- spVineCopula(copula at spCop, vineCopFit at copula)
+
+ return(new("fitCopula", estimate = spVineCop at parameters, var.est = matrix(NA),
+ method = method,
+ loglik = sum(dCopula(data at data, spVineCop, h=data at distances, log=T)),
+ fitting.stats=list(convergence = as.integer(NA)),
+ nsample = nrow(data at data), copula=spVineCop))
}
-setMethod("fitCopula",signature=signature("spVineCopula"),fitSpVine)
\ No newline at end of file
+setMethod("fitCopula",signature=signature("spVineCopula"),fitSpVine)
+
+# conditional spatial vine
+condSpVine <- function(condVar, dists, spVine, n=100) {
+ rat <- 0.2/(1:(n/2))-(0.1/((n+1)/2))
+ xVals <- unique(sort(c(rat,1-rat,1:(n-1)/(n))))
+ xLength <- length(xVals)
+ repCondVar <- matrix(condVar, ncol=length(condVar), nrow=xLength, byrow=T)
+ density <- dCopula(cbind(xVals, repCondVar), spVine, h=dists)
+
+ linAppr <- approxfun(c(0,xVals,1), density[c(1,1:xLength,xLength)] ,yleft=0, yright=0)
+ int <- integrate(linAppr,lower=0, upper=1)$value
+
+ return(function(u) linAppr(u)/int)
+}
+
+# interpolation
+
+spCopPredict.expectation <- function(predNeigh, spVine, margin) {
+ stopifnot(is.function(margin$d))
+ stopifnot(is.function(margin$p))
+
+ predMean <- NULL
+ for(i in 1:nrow(predNeigh at data)) { # i <-1
+ condSecVine <- condSpVine(as.numeric(predNeigh at data[i,]), predNeigh at distances[i,], spVine)
+
+ condExp <- function(x) {
+ condSecVine(margin$p(x))*margin$d(x)*x
+ }
+
+ predMean <- c(predMean, integrate(condExp,0,3000,subdivisions=1e6)$value)
+ }
+ if ("data" %in% slotNames(predNeigh at locations)) {
+ res <- predNeigh at locations
+ res at data[["expect"]] <- predMean
+ return(res)
+ } else {
+ predMean <- data.frame(predMean)
+ colnames(predMean) <- "expect"
+ return(addAttrToGeom(predNeigh at locations, predMean, match.ID=FALSE))
+ }
+}
+
+spCopPredict.quantile <- function(predNeigh, spVine, margin, p=0.5) {
+ stopifnot(is.function(margin$q))
+
+ predQuantile <- NULL
+ for(i in 1:nrow(predNeigh at data)) { # i <-1
+ condSecVine <- condSpVine(as.numeric(predNeigh at data[i,]), predNeigh at distances[i,], spVine)
+ pPred <- optimise(function(x) abs(integrate(condSecVine,0,x)$value-p),
+ c(0,1))$minimum
+ predQuantile <- c(predQuantile, margin$q(pPred))
+ }
+ if ("data" %in% slotNames(predNeigh at locations)) {
+ res <- predNeigh at locations
+ res at data[[paste("quantile.",p,sep="")]] <- predQuantile
+ return(res)
+ } else {
+ predQuantile <- data.frame(predQuantile)
+ colnames(predQuantile) <- paste("quantile.",p,sep="")
+ return(addAttrToGeom(predNeigh at locations, predQuantile, match.ID=FALSE))
+ }
+}
+
+spCopPredict <- function(predNeigh, spVine, margin, method="quantile", p=0.5) {
+ switch(method,
+ quantile=spCopPredict.quantile(predNeigh, spVine, margin, p),
+ expectation=spCopPredict.expectation(predNeigh, spVine, margin))
+}
\ No newline at end of file
Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R 2013-03-27 19:13:09 UTC (rev 89)
+++ pkg/R/spatialPreparation.R 2013-04-04 17:50:25 UTC (rev 90)
@@ -11,12 +11,14 @@
# sp="SpatialPoints" SpatialPoints object providing the coordinates
# index="matrix" linking the obs. in data to the coordinates
-neighbourhood <- function(data, distances, sp, index){
- varNames <- names(data)[[1]]
+neighbourhood <- function(data, distances, sp, index, prediction, varNames){
sizeN <- ncol(distances)+1
data <- as.data.frame(data)
- colnames(data) <- paste(paste("N",rep(0:(sizeN-1),each=length(varNames)),sep=""),rep(varNames,sizeN),sep=".")
- new("neighbourhood", data=data, distances=distances, coords=sp at coords, bbox=sp at bbox, proj4string=sp at proj4string, index=index, varNames=varNames)
+ colnames(data) <- paste(paste("N", rep((0+prediction):(sizeN-1), each=length(varNames)), sep=""),
+ rep(varNames,(sizeN-prediction)),sep=".")
+ new("neighbourhood", data=data, distances=distances, locations=sp,
+ bbox=sp at bbox, proj4string=sp at proj4string, index=index, varNames=varNames,
+ prediction=prediction)
}
## show
@@ -45,60 +47,58 @@
## calculate neighbourhood from SpatialPointsDataFrame
# returns an neighbourhood object
-# spData spatialPointsDataFrame
-# var one or multiple variable names, all is the default
-# size the size of the neighbourhood, default of 5
-# dep denoting a subset of dependent locations (default NULL: all locations will be used)
-# indep denoting a subset of independent locations (default NULL: all locations will be used)
-# no location will be paired with itself
-getNeighbours <- function(spData,var=names(spData),size=4,dep=NULL,indep=NULL,min.dist=10){
-nLocs <- length(spData)
-distMat <- spDists(spData)
-if(min.dist>0) distMat[distMat<min.dist] <- Inf
-if ( any(is.na( match(var,names(spData)) )) )
- stop("At least one of the variables is unkown or is not part of the data.")
-if(is.null(dep) & !is.null(indep)) dep <- 1:nLocs[-indep]
-if(!is.null(dep) & is.null(indep)) indep <- 1:nLocs[-dep]
-if(!is.null(dep) & !is.null(indep)) {
- cat("Reduced distance matrix is used: (",paste(dep,collapse=", "),") x (",paste(indep,collapse=", "),")",sep="")
-} else {
- dep <- 1:nLocs
- indep <- 1:nLocs
-}
+# spData spatialPointsDataFrame
+# locations the prediction locations, for fitting, locations=spData
+# var one or multiple variable names, all is the default
+# size the size of the neighbourhood, note that for prediction the size
+# is one less than for the copula estimation (default of 5)
+# min.dist the minimum distance between neighbours, must be positive
-size <- min(size,length(indep)-1)
-if (size > sizeLim) {
- stop(paste("Evaluation of copulas might take a long time for more than",
- sizeLim," neighbours. Increase sizeLim if you want to evaluate neighbourhoods with",
- size,"locations."))
-}
+getNeighbours <- function(spData, locations, var=names(spData), size=5,
+ prediction=FALSE, min.dist=0.01) {
+
+ stopifnot((!prediction && missing(locations)) || (prediction && !missing(locations)))
+ stopifnot(min.dist>0 || prediction)
+
+ if(missing(locations) && !prediction)
+ locations=spData
+
+ stopifnot(is(locations,"Spatial"))
+
+ nLocs <- length(locations)
+
+ if(any(is.na(match(var,names(spData)))))
+ stop("At least one of the variables is unkown or is not part of the data.")
-lData <- vector("list",size)
-index <- NULL
-dists <- NULL
+ size <- min(size,length(spData))
-for (i in dep) {
- nbrs <- matrix(Inf,ncol=2,nrow=size-1)
- ind <- logical(nLocs)
- ind[indep] <- TRUE
- ind[i] <- FALSE
- for (j in (1:nLocs)[ind]) {
- tmpDist <- distMat[i,j]
- if (any(tmpDist < nbrs[,1])) {
- nbrs[size-1,] <- c(tmpDist,j)
- nbrs <- nbrs[order(nbrs[,1]),]
- }
+ allDists <- NULL
+ allLocs <- NULL
+ allData <- NULL
+
+ for(i in 1:length(locations)) { # i <- 1
+ tempDists <- spDistsN1(spData,locations[i,])
+ tempDists[tempDists < min.dist] <- Inf
+ spLocs <- order(tempDists)[1:(size-1)]
+ allLocs <- rbind(allLocs, spLocs)
+ allDists <- rbind(allDists, tempDists[spLocs])
+
+ if(!prediction)
+ spLocs <- c(i,spLocs)
+ allData <- rbind(allData, as.vector(spData[spLocs, var, drop=F]@data[[1]]))
}
- lData[[1]] <- rbind(lData[[1]],spData at data[i,var,drop=FALSE])
- for (nbr in 2:size) {
- lData[[nbr]] <- rbind(lData[[nbr]],spData at data[nbrs[nbr-1,2],var,drop=FALSE])
- }
- index <- rbind(index, c(i,nbrs[,2]))
- dists <- rbind(dists, c(nbrs[,1]))
+
+ return(neighbourhood(allData, allDists, locations,
+ allLocs, prediction, var))
}
-return(neighbourhood(as.data.frame(lData), dists, SpatialPoints(spData), index))
-}
+# data(meuse)
+# coordinates(meuse) <- ~x+y
+# meuseNeigh <- getNeighbours(meuse,var="zinc",size=5)
+# str(meuseNeigh)
+#
+# meuseNeigh <- addAttrToGeom(meuseNeigh at locations,data.frame(rnd=runif(155)),match.ID=F)
+# str(meuseNeigh)
#############
## BINNING ##
Modified: pkg/R/vineCopulas.R
===================================================================
--- pkg/R/vineCopulas.R 2013-03-27 19:13:09 UTC (rev 89)
+++ pkg/R/vineCopulas.R 2013-04-04 17:50:25 UTC (rev 90)
@@ -14,13 +14,13 @@
RVM <- D2RVine(1:RVM,rep(0,RVM*(RVM-1)/2),rep(0,RVM*(RVM-1)/2))
}
- # handling non S4-class as sub-element in a S4-class
stopifnot(class(RVM)=="RVineMatrix")
- class(RVM) <- "list"
ltr <- lower.tri(RVM$Matrix)
copDef <- cbind(RVM$family[ltr], RVM$par[ltr], RVM$par2[ltr])
- copulas <- rev(apply(copDef,1, function(x) copulaFromFamilyIndex(x[1],x[2],x[3])))
+ copulas <- rev(apply(copDef,1, function(x) {
+ copulaFromFamilyIndex(x[1],x[2],x[3])
+ }))
new("vineCopula", copulas=copulas, dimension = as.integer(nrow(RVM$Matrix)),
RVM=RVM, parameters = numeric(),
@@ -45,7 +45,7 @@
dRVine <- function(u, copula, log=F) {
RVM <- copula at RVM
- class(RVM) <- "RVineMatrix"
+# class(RVM) <- "RVineMatrix"
vineLoglik <- RVineLogLik(u, RVM, separate=T)$loglik
if(log)
return(vineLoglik)
@@ -54,111 +54,15 @@
}
setMethod("dCopula", signature("numeric","vineCopula"),
- function(u, copula, log, ...) dRVine(matrix(u, ncol=copula at dimension), copula, log, ...))
+ function(u, copula, log, ...) {
+ dRVine(matrix(u, ncol=copula at dimension), copula, log, ...)
+ })
setMethod("dCopula", signature("matrix","vineCopula"), dRVine)
+setMethod("dCopula", signature("data.frame","vineCopula"),
+ function(u, copula, log, ...) {
+ dRVine(as.matrix(u), copula, log, ...)
+ })
-# ## d-vine structure
-#
-# # copula <- vineFit
-# # u <- empVine
-# # empCopVine
-#
-# # dDvine(vineFit, empVine,log=T)
-#
-# dDvine <- function(copula, u, log=FALSE){
-# dim <- copula at dimension
-# tmp <- u
-# u <- NULL
-# u[[1]] <- matrix(tmp,ncol=dim)
-#
-# den <- rep(1,nrow(u[[1]]))
-#
-# newU <- NULL
-# for (i in 1:(dim-1)) {
-# tmpCop <- copula at copulas[[i]]
-# tmpU <- u[[1]][,i:(i+1)]
-# if(log)
-# den <- den + dCopula(tmpU, tmpCop,log=T)
-# else
-# den <- den*dCopula(tmpU,tmpCop,log=F)
-# if (i == 1) {
-# newU <- cbind(newU, ddvCopula(tmpU, tmpCop))
-# } else {
-# newU <- cbind(newU, dduCopula(tmpU, tmpCop))
-# }
-# if (1<i & i<(dim-1)) {
-# newU <- cbind(newU, ddvCopula(tmpU, tmpCop))
-# }
-# }
-# u[[2]] <- newU
-#
-# used <- dim-1
-# for (l in 2:(dim-1)) {
-# newU <- NULL
-# for (i in 1:(dim-l)) {
-# # cat(used+i,"\n")
-# tmpCop <- copula at copulas[[used+i]]
-# tmpU <- u[[l]][,(i*2-1):(i*2)]
-# if(log)
-# den <- den + dCopula(tmpU, tmpCop,log=T)
-# else
-# den <- den*dCopula(tmpU, tmpCop, log=F)
-# if (l < dim-1) {
-# if (i == 1) {
-# newU <- cbind(newU,ddvCopula(tmpU, tmpCop))
-# } else {
-# newU <- cbind(newU,dduCopula(tmpU, tmpCop))
-# }
-# if (1<i & i<(dim-1)) {
-# newU <- cbind(newU,ddvCopula(tmpU, tmpCop))
-# }
-# }
-# }
-# u[[l+1]] <- newU
-# used <- used + dim - l
-# }
-#
-# return(den)
-# }
-#
-# ## c-vine structure
-#
-# dCvine <- function(copula, u) {
-# # cat("c-vine \n")
-# dim <- copula at dimension
-# tmp <- u
-# u <- NULL
-# u[[1]] <- matrix(tmp,ncol=dim)
-#
-# den <- rep(1,nrow(u[[1]]))
-#
-# used <- 0 # already used copulas
-#
-# for (l in 1:(dim-1)) {
-# newU <- NULL
-# for (i in 1:(dim-l)) {
-# # cat(used+i,"\n")
-# tmpCop <- copula at copulas[[used+i]]
-# tmpU <- u[[l]][,c(1,(i+1))]
-# den <- den*dCopula(tmpU, tmpCop)
-# if(l < (dim-1)) newU <- cbind(newU,dduCopula(tmpU, tmpCop))
-# }
-# if(l < (dim-1)) {
-# u[[l+1]] <- newU
-# used <- used + dim - l
-# }
-# }
-#
-# return(den)
-# }
-#
-# ##
-#
-# dvineCopula <- function(u, copula, log=F) {
-# den <- switch(getNumType(copula),dCvine ,dDvine)
-# return(den(copula, u, log))
-# }
-
## jcdf ##
pvineCopula <- function(u, copula) {
empCop <- genEmpCop(copula, 1e5)
@@ -166,35 +70,17 @@
return(pCopula(u, empCop))
}
-setMethod("pCopula", signature("numeric","vineCopula"), pvineCopula)
+setMethod("pCopula", signature("numeric","vineCopula"),
+ function(u,copula) {
+ pvineCopula(matrix(u, ncol=copula at dimension),copula)
+ })
+setMethod("pCopula", signature("data.frame","vineCopula"),
+ function(u,copula) pvineCopula(as.matrix(u),copula))
setMethod("pCopula", signature("matrix","vineCopula"), pvineCopula)
-
-## random numbers
-# linkVineCopSim <- function(n, copula) {
-# numType <- getNumType(copula)
-#
-# getFamily <- function(copula) {
-# if("family" %in% slotNames(copula)) numFam <- copula at family
-# else {
-# numFam <- switch(class(copula)[1], normalCopula=1, tCopula=2, claytonCopula=3, gumbelCopula=4, frankCopula=5)
-# }
-# }
-#
-# par1 <- unlist(lapply(copula at copulas,function(x) x at parameters[1]))
-# par2 <- unlist(lapply(copula at copulas,function(x) x at parameters[2]))
-# par2[is.na(par2)] <- 0
-# numFam <- unlist(lapply(copula at copulas,getFamily))
-# tcops <- which(numFam==2) #? length(which(5==3))
-# if(length(tcops)>0)
-# par2[tcops] <- unlist(lapply(copula at copulas[tcops], function(x) x at df))
-#
-# return(RVineSim(n, C2RVine(1:copula at dimension, numFam, par1, par2)))
-# }
-
rRVine <- function(n, copula) {
RVM <- copula at RVM
- class(RVM) <- "RVineMatrix"
+# class(RVM) <- "RVineMatrix"
RVineSim(n, RVM)
}
@@ -204,9 +90,16 @@
fitVineCop <- function(copula, data, method) {
stopifnot(copula at dimension==ncol(data))
if("StructureSelect" %in% method)
- vineCopula(RVineStructureSelect(data, indeptest="indeptest" %in% method))
+ vineCop <- vineCopula(RVineStructureSelect(data, indeptest="indeptest" %in% method))
else
- vineCopula(RVineCopSelect(data, Matrix=copula at RVM$Matrix, indeptest="indeptest" %in% method))
+ vineCop <- vineCopula(RVineCopSelect(data, Matrix=copula at RVM$Matrix,
+ indeptest="indeptest" %in% method))
+
+ return(new("fitCopula", estimate = vineCop at parameters, var.est = matrix(NA),
+ method = method,
+ loglik = RVineLogLik(data, vineCop at RVM)$loglik,
+ fitting.stats=list(convergence = as.integer(NA)),
+ nsample = nrow(data), copula=vineCop))
}
setMethod("fitCopula", signature=signature("vineCopula"), fitVineCop)
\ No newline at end of file
Modified: spcopula_0.1-1.tar.gz
===================================================================
(Binary files differ)
More information about the spcopula-commits
mailing list