[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