[spcopula-commits] r94 - / pkg pkg/R pkg/man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Apr 23 16:17:12 CEST 2013


Author: ben_graeler
Date: 2013-04-23 16:17:11 +0200 (Tue, 23 Apr 2013)
New Revision: 94

Modified:
   pkg/DESCRIPTION
   pkg/R/Classes.R
   pkg/R/partialDerivatives.R
   pkg/R/returnPeriods.R
   pkg/R/spCopula.R
   pkg/R/spVineCopula.R
   pkg/R/spatialPreparation.R
   pkg/man/condSpVine.Rd
   pkg/man/fitCorFun.Rd
   pkg/man/spCopPredict.Rd
   spcopula_0.1-1.tar.gz
   spcopula_0.1-1.zip
Log:
- fine tuning of condSpVine
- weighted estimate in fitCorFun

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/DESCRIPTION	2013-04-23 14:17:11 UTC (rev 94)
@@ -2,7 +2,7 @@
 Type: Package
 Title: copula driven spatial analysis
 Version: 0.1-1
-Date: 2013-04-11
+Date: 2013-04-23
 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/R/Classes.R
===================================================================
--- pkg/R/Classes.R	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/Classes.R	2013-04-23 14:17:11 UTC (rev 94)
@@ -95,8 +95,10 @@
   check.upper <- NULL
   check.lower <- NULL
   
+  nComp <- length(object at components)
   if(!is.null(object at calibMoa(normalCopula(0),0))) {
-    for (i in 1:(length(object at components)-1)) {
+    nonIndep <- sapply(object at components[-nComp], function(x) class(x) != "indepCopula")
+    for (i in (1:(nComp-1))[nonIndep]) {
       check.upper <- c(check.upper, is.na(object at calibMoa(object at components[[i]], object at distances[i+1])))
       check.lower <- c(check.lower, is.na(object at calibMoa(object at components[[i]], c(0,object at distances)[i])))
     }

Modified: pkg/R/partialDerivatives.R
===================================================================
--- pkg/R/partialDerivatives.R	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/partialDerivatives.R	2013-04-23 14:17:11 UTC (rev 94)
@@ -218,7 +218,7 @@
   u1 <- u[,1]
   u2 <- u[,2]
 
-  pcopula(gumbelCopula(rho),u) * ((-log(u1))^rho+(-log(u2))^rho)^(1/rho-1) * (-log(u1))^(rho-1)/u1
+  pCopula(u,gumbelCopula(rho)) * ((-log(u1))^rho+(-log(u2))^rho)^(1/rho-1) * (-log(u1))^(rho-1)/u1
 }
 
 setMethod("dduCopula", signature("numeric","gumbelCopula"),
@@ -237,7 +237,7 @@
   u1 <- u[,1]
   u2 <- u[,2]
 
-  pcopula(gumbelCopula(rho),u) * ((-log(u2))^rho+(-log(u1))^rho)^(1/rho-1) * (-log(u2))^(rho-1)/u2
+  pCopula(u,gumbelCopula(rho)) * ((-log(u2))^rho+(-log(u1))^rho)^(1/rho-1) * (-log(u2))^(rho-1)/u2
 }
 
 setMethod("ddvCopula", signature("numeric","gumbelCopula"),

Modified: pkg/R/returnPeriods.R
===================================================================
--- pkg/R/returnPeriods.R	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/returnPeriods.R	2013-04-23 14:17:11 UTC (rev 94)
@@ -123,4 +123,4 @@
   return(cbind(u,params))
 }
 
-setMethod("qCopula_u",signature("copula"),qCopula_u.def)
\ No newline at end of file
+setMethod("qCopula_u", signature("copula"), qCopula_u.def)

Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/spCopula.R	2013-04-23 14:17:11 UTC (rev 94)
@@ -28,7 +28,8 @@
                        id=function(copula, h) return(h))
     
     for (i in 1:length(components)) {
-      param <- try(calibMoa(components[[i]], distances[i]),T)
+      if(class(components[[i]]) != "indepCopula")
+        param <- try(calibMoa(components[[i]], distances[i]),T)
       if (class(param) == "numeric")
         components[[i]]@parameters[1] <- param # take care of non single parameter copulas
     }
@@ -36,10 +37,26 @@
 
 #   components <- append(components,indepCopula())
   
-  param       <- unlist(lapply(components, function(x) x at parameters))
-  param.names <- unlist(lapply(components, function(x) x at param.names))
-  param.low   <- unlist(lapply(components, function(x) x at param.lowbnd))
-  param.up    <- unlist(lapply(components, function(x) x at param.upbnd))
+  param       <- unlist(lapply(components, 
+                               function(x) {
+                                 if(class(x)=="indepCopula") 
+                                   return(NA)
+                                 x at parameters}))
+  param.names <- unlist(lapply(components, 
+                               function(x) {
+                                 if(class(x)=="indepCopula") 
+                                   return(NA)
+                                 x at param.names}))
+  param.low   <- unlist(lapply(components, 
+                               function(x) {
+                                 if(class(x)=="indepCopula") 
+                                   return(NA)
+                                 x at param.lowbnd}))
+  param.up    <- unlist(lapply(components, 
+                               function(x) {
+                                 if(class(x)=="indepCopula") 
+                                   return(NA)
+                                 x at param.upbnd}))
      
   new("spCopula", dimension=as.integer(2), parameters=param, param.names=param.names,
       param.lowbnd=param.low, param.upbnd=param.up,
@@ -75,6 +92,7 @@
   n.dists <- length(dists)
   calibPar <- copula at calibMoa
   
+  # data is sorted to be ascending in distance
   res <- numeric(0)
   sel <- which(h < dists[1])
   if(sum(sel)>0) {
@@ -136,7 +154,7 @@
   sel <- which(h >= dists[n.dists])
   if(sum(sel)>0) {
     res <- c(res, fun(pairs[which(h >= dists[n.dists]),],
-                      copula at components[[n.dists]],, ...))
+                      copula at components[[n.dists]], ...))
   }
 
   return(res)
@@ -152,42 +170,39 @@
     tmpCop <- copula at components[[1]]
     if(class(tmpCop) != "indepCopula")
       tmpCop at parameters[1] <- calibPar(tmpCop, h)
-    res <- fun(pairs, tmpCop, ...)
+    return(fun(pairs, tmpCop, ...))
   }
   
-  if (n.dists >= 2) {
-    for ( i in 2:n.dists ) {
-      low  <- dists[i-1]
-      high <- dists[i]
-      if (h >= low & h < high) {
-        lowerCop <- copula at components[[i-1]]
-        upperCop <- copula at components[[i]]
-        if (class(lowerCop) != class(upperCop)) {
-          if(class(lowerCop) != "indepCopula")
-            lowerCop at parameters[1] <- calibPar(lowerCop, h)
-          if(class(upperCop) != "indepCopula")
-            upperCop at parameters[1] <- calibPar(upperCop, h)
+  if(h >= dists[n.dists]) {
+    return(fun(pairs, copula at components[[n.dists]], ...))
+  }
+  
+  for (i in 2:n.dists) {
+    low  <- dists[i-1]
+    high <- dists[i]
+    if(low <= h & h < high) {
+      lowerCop <- copula at components[[i-1]]
+      upperCop <- copula at components[[i]]
+      if (class(lowerCop) != class(upperCop)) {
+        if(class(lowerCop) != "indepCopula")
+          lowerCop at parameters[1] <- calibPar(lowerCop, h)
+        if(class(upperCop) != "indepCopula")
+          upperCop at parameters[1] <- calibPar(upperCop, h)
 
-          lowerVals <- fun(pairs, lowerCop)
-          upperVals <- fun(pairs, upperCop)
+        lowerVals <- fun(pairs, lowerCop)
+        upperVals <- fun(pairs, upperCop)
 
-          res <- (high-h)/(high-low)*lowerVals + (h-low)/(high-low)*upperVals
-          if(do.logs)
-            res <- log(res)
-        } else {
-          if(class(lowerCop) != "indepCopula")
-            lowerCop at parameters <- calibPar(lowerCop, h)
-          res <- fun(pairs, lowerCop, ...)
-        }
+        res <- (high-h)/(high-low)*lowerVals + (h-low)/(high-low)*upperVals
+        if(do.logs)
+          return(log(res))
+        return(res)
+      } else {
+        if(class(lowerCop) != "indepCopula")
+          lowerCop at parameters <- calibPar(lowerCop, h)
+        return(fun(pairs, lowerCop, ...))
       }
     }
   }
-  
-  if(h >= dists[n.dists]) {
-    res <- fun(pairs, copula at components[[n.dists]], ...)
-  }
-  
-  return(res)
 }
 
 
@@ -232,106 +247,63 @@
 }
 
 
-# u 
-#   two column matrix providing the transformed pairs and their respective 
-#   separation distances
-# h providing the separating distance(s)
-# block
-#   block distances, pairs are assumed to be ordered block wise
+## spatial copula CDF
+######################
+
 pSpCopula <- function (u, copula, h, block=1) {
-  if (missing(h)) stop("Point pairs need to be provided with their separating distance \"h\".")
-  
-  n <- nrow(u)
-  
-  if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
-  
-  if(length(h)>1 && length(h)!=n) {
+  if (missing(h)) 
+    stop("Point pairs need to be provided with their separating distance \"h\".")
+  if(length(h)>1 && length(h)!=nrow(u))
     stop("The distance vector must either be of the same length as rows in the data pairs or a single value.")
-  }
   
+  if(is.null(copula at calibMoa(normalCopula(0),0)))
+    return(spConCop(pCopula, copula, u, rep(h,length.out=nrow(u))))
   
-  if(is.null(copula at calibMoa(normalCopula(0),0))) {
-    res <- spConCop(pCopula, copula, u, rep(h,length.out=nrow(pairs)))
-  } else {
-    if(length(h)>1) {
-      if (block == 1){
-        ordering <- order(h)
+  if(length(h)>1) {
+    ordering <- order(h)
         
-        # ascending sorted pairs allow for easy evaluation
-        u <- u[ordering,,drop=FALSE] 
-        h <- h[ordering]
+    # ascending sorted pairs allow for easy evaluation
+    u <- u[ordering,,drop=FALSE] 
+    h <- h[ordering]
         
-        res <- spDepFunCop(pCopula, copula, u, h)
+    res <- spDepFunCop(pCopula, copula, u, h)
         
-        # reordering the values
-        res <- res[order(ordering)]
-      } else {
-        res <- NULL
-        for(i in 1:(n%/%block)) {
-          res <- c(res, spDepFunCopSnglDist(pCopula, copula, 
-                                            u[((i-1)*block+1):(i*block),], 
-                                            h[i*block]))
-        }
-      }
-    } else {
-      res <- spDepFunCopSnglDist(pCopula, copula, u, h)
-    }
-  }
-  
-  return(res)
+    # reordering the values
+    return(res[order(ordering)])
+  } else
+      return(spDepFunCopSnglDist(pCopula, copula, u, h))
 }
 
 setMethod(pCopula, signature("numeric","spCopula"), 
           function(u, copula, ...) pSpCopula(matrix(u,ncol=2),copula, ...))
 setMethod(pCopula, signature("matrix","spCopula"), pSpCopula)
 
-## spatial Copula density ##
+## spatial Copula density 
+##########################
 
-# u 
-#   three column matrix providing the transformed pairs and their respective 
-#   separation distances
-dSpCopula <- function (u, copula, log, h, block=1) {
-  if (missing(h)) stop("Point pairs need to be provided with their separating distance \"h\".")
-  
-  n <- nrow(u)
-  
-  if (n%%block != 0) stop("The block size is not a multiple of the data length:",n)
-  
-  if(length(h)>1 && length(h)!=n) {
+dSpCopula <- function (u, copula, log, h) {
+  if (missing(h)) 
+    stop("Point pairs need to be provided with their separating distance \"h\".")
+  if(length(h)>1 && length(h)!=nrow(u))
     stop("The distance vector must either be of the same length as rows in the data pairs or a single value.")
-  }
   
-  if(is.null(copula at calibMoa(normalCopula(0),0))){
-    res <- spConCop(dCopula, copula, u, rep(h, length.out=n), 
-                    do.logs=log, log=log)
-  }
-  else {
-    if(length(h)>1) {
-      if (block == 1){
-        ordering <- order(h)
+  if(is.null(copula at calibMoa(normalCopula(0),0)))
+    return(spConCop(dCopula, copula, u, rep(h, length.out=nrow(u)), do.logs=log,
+                    log=log))
+
+  if(length(h)>1) {
+    ordering <- order(h)
+    
+    # ascending sorted pairs allow for easy evaluation
+    u <- u[ordering,,drop=FALSE] 
+    h <- h[ordering]
         
-        # ascending sorted pairs allow for easy evaluation
-        u <- u[ordering,,drop=FALSE] 
-        h <- h[ordering]
+    res <- spDepFunCop(dCopula, copula, u, h, do.logs=log, log=log)
         
-        res <- spDepFunCop(dCopula, copula, u, h, do.logs=log, log=log)
-        
-        # reordering the values
-        res <- res[order(ordering)]
-      } else {
-        res <- NULL
-        for(i in 1:(n%/%block)) {
-          res <- c(res, spDepFunCopSnglDist(dCopula, copula, 
-                                            u[((i-1)*block+1):(i*block),], 
-                                            h[i*block], do.logs=log, log=log))
-        }
-      }
-    } else {
-      res <- spDepFunCopSnglDist(dCopula, copula, u, h, do.logs=log, log=log)
-    }
-  }
-  
-  return(res)
+    # reordering the values
+    return(res[order(ordering)])
+  } else
+      return(spDepFunCopSnglDist(dCopula, copula, u, h, do.logs=log, log=log))
 }
 
 setMethod(dCopula, signature("numeric","spCopula"), 
@@ -339,56 +311,31 @@
 setMethod(dCopula, signature("matrix","spCopula"), dSpCopula)
 
 ## partial derivatives ##
-
 ## dduSpCopula
 ###############
 
-dduSpCopula <- function (u, copula, h, block=1) {
-  if (missing(h)) stop("Point pairs need to be provided with their separating distance h.")
-  
-  n <- nrow(u)
-  
-  if(length(h)>1 && length(h)!=n) {
+dduSpCopula <- function (u, copula, h) {
+  if (missing(h)) 
+    stop("Point pairs need to be provided with their separating distance h.")
+  if(length(h)>1 && length(h)!=nrow(u))
     stop("The distance vector must either be of the same length as rows in the data pairs or a single value.")
-  }
 
   if(is.null(copula at calibMoa(normalCopula(0),0)))
-    res <- spConCop(dduCopula, copula, u, rep(h, length.out=n))
+    return(spConCop(dduCopula, copula, u, rep(h, length.out=nrow(u))))
   
-  else {
-    if(length(h)>1) {
-      if (block == 1){
-        ordering <- order(h)
+  if(length(h)>1) {
+    ordering <- order(h)
         
-        # ascending sorted pairs allow for easy evaluation
-        u <- u[ordering,,drop=FALSE] 
-        h <- h[ordering]
+    # ascending sorted pairs allow for easy evaluation
+    u <- u[ordering, , drop=FALSE] 
+    h <- h[ordering]
         
-        res <- spDepFunCop(dduCopula, copula, u, h)      
+    res <- spDepFunCop(dduCopula, copula, u, h)      
         
-        # reordering the values
-        res <- res[order(ordering)]
-      } else {
-        res <- NULL
-        for(i in 1:(n%/%block)) {
-          res <- c(res, spDepFunCopSnglDist(dduCopula, copula, 
-                                            u[((i-1)*block+1):(i*block),], 
-                                            h[i*block]))
-        }
-      }
-    } else {
-      res <- spDepFunCopSnglDist(dduCopula, copula, u, h)
-    }
-  }
-  
-  if(any(res < 0) || any(res > 1)) {
-    warning("Partial derivative produced values outside of [0,1], corrections will be applied to:\n",
-            paste(res[res<0],collapse=" "),paste(res[res>1],collapse=" "))
-    res[res<0] <- 0
-    res[res>1] <- 1
-  }
-  
-  return(res)
+    # reordering the values
+    return(res[order(ordering)])
+  } else
+    return(spDepFunCopSnglDist(dduCopula, copula, u, h))
 }
 
 setMethod("dduCopula", signature("matrix","spCopula"), dduSpCopula)
@@ -398,53 +345,28 @@
 ## ddvSpCopula
 ###############
 
-
 ddvSpCopula <- function (u, copula, h, block=1) {
-  if (missing(h)) stop("Point pairs need to be provided with their separating distance h.")
-  
-  n <- nrow(u)
-  
-  if(length(h)>1 && length(h)!=n) {
+  if (missing(h)) 
+    stop("Point pairs need to be provided with their separating distance h.")
+  if(length(h)>1 && length(h)!=nrow(u))
     stop("The distance vector must either be of the same length as rows in the data pairs or a single value.")
-  }
-  
+    
   if(is.null(copula at calibMoa(normalCopula(0),0)))
-    res <- spConCop(dduCopula, copula, u, rep(h, length.out=n))
+    return(spConCop(ddvCopula, copula, u, rep(h, length.out=nrow(u))))
   
-  else {
-    if(length(h)>1) {
-      if (block == 1){
-        ordering <- order(h)
+  if(length(h)>1) {
+    ordering <- order(h)
         
-        # ascending sorted pairs allow for easy evaluation
-        u <- u[ordering,,drop=FALSE] 
-        h <- h[ordering]
+    # ascending sorted pairs allow for easy evaluation
+    u <- u[ordering,,drop=FALSE] 
+    h <- h[ordering]
         
-        res <- spDepFunCop(ddvCopula, copula, u, h)      
+    res <- spDepFunCop(ddvCopula, copula, u, h)      
         
-        # reordering the values
-        res <- res[order(ordering)]
-      } else {
-        res <- NULL
-        for(i in 1:(n%/%block)) {
-          res <- c(res, spDepFunCopSnglDist(ddvCopula, copula, 
-                                            u[((i-1)*block+1):(i*block),], 
-                                            h[i*block]))
-        }
-      }
-    } else {
-      res <- spDepFunCopSnglDist(ddvCopula, copula, u, h)
-    }
-  }
-  
-  if(any(res < 0) || any(res > 1)) {
-    warning("Partial derivative produced values outside of [0,1], corrections will be applied to:\n",
-            paste(res[res<0],collapse=" "),paste(res[res>1],collapse=" "))
-    res[res<0] <- 0
-    res[res>1] <- 1
-  }
-  
-  return(res)
+    # reordering the values
+    return(res[order(ordering)])
+  } else
+      return(spDepFunCopSnglDist(ddvCopula, copula, u, h))
 }
 
 setMethod("ddvCopula", signature("matrix","spCopula"), ddvSpCopula)
@@ -485,21 +407,29 @@
 # bounds -> the bounds of the correlation function (typically c(0,1))
 # method -> the measure of association, either "kendall" or "spearman"
 fitCorFun <- function(bins, degree=3, cutoff=NA, bounds=c(0,1), 
-                      cor.method=NULL) {
+                      cor.method=NULL, weighted=FALSE) {
   if(is.null(cor.method)) {
     if(is.null(attr(bins,"cor.method")))
       stop("Neither the bins arguments has an attribute cor.method nor is the parameter cor.method provided.") 
     else 
       cor.method <- attr(bins,"cor.method")
-  } else 
+  } else {
     if(!is.null(attr(bins,"cor.method")) && cor.method != attr(bins,"cor.method"))
       stop("The cor.method attribute of the bins argument and the argument cor.method do not match.")
-    
-  bins <- as.data.frame(bins[1:2])
-  if(!is.na(cutoff)) bins <- bins[which(bins[[1]] <= cutoff),]
+  }
+
+  if (weighted) {
+    bins <- as.data.frame(bins[c("np","meanDists","lagCor")])
+    if(!is.na(cutoff)) 
+      bins <- bins[bins$meanDists <= cutoff,]
+    fitCor <- lm(lagCor ~ poly(meanDists, degree), data = bins, weights=bins$np)
+  } else {
+    bins <- as.data.frame(bins[c("meanDists","lagCor")])
+    if(!is.na(cutoff)) 
+      bins <- bins[bins$meanDists <= cutoff,]
+    fitCor <- lm(lagCor ~ poly(meanDists, degree), data = bins)
+  }
   
-  fitCor <- lm(lagCor ~ poly(meanDists, degree), data = bins)
-  
   print(fitCor)
   cat("Sum of squared residuals:",sum(fitCor$residuals^2),"\n")
   
@@ -564,35 +494,6 @@
   }
 }
 
-# older implementation:
-# composeSpCop <- function(bestFit, families, bins, calcCor) {
-#   nfits <- length(bestFit)
-#   gaps <- which(diff(bestFit)!=0)
-# 
-#   if(missing(calcCor)) noCor <- nfits
-#   else noCor <- min(which(calcCor(bins$meanDists)<=0), nfits)
-#   
-#   breaks <- sort(c(gaps, gaps+1, noCor))
-#   breaks <- breaks[breaks<noCor]
-#   
-#   cops <- as.list(families[bestFit[breaks]])
-#   
-#   breaks <- unique(c(breaks, min(nfits,rev(breaks)[1]+1)))
-#   distances <- bins$meanDists[breaks]
-#   
-#   if(length(breaks)==length(cops)) {
-#     warning("Lags do not cover the entire range with positive correlation. ", 
-#              "An artifical one fading towards the indepCopula has been added.")
-#     distances <- c(distances, 
-#                    rev(distances)[1]+diff(bins$meanDists[nfits+c(-1,0)]))
-#   }
-# 
-#   if(missing(calcCor)) return(spCopula(components=cops, distances=distances, 
-#                                        unit="m"))
-#   else return(spCopula(components=cops, distances=distances, unit="m", 
-#                        spDepFun=calcCor))
-# }
-
 # in once
 
 # bins   -> typically output from calcBins

Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/spVineCopula.R	2013-04-23 14:17:11 UTC (rev 94)
@@ -24,13 +24,14 @@
   l0 <- rep(0,nrow(u)) # level 0 (spatial) density
   u0 <- NULL # level 0 conditional data
   
-  if(!is.matrix(h)) h <- matrix(h, ncol=length(h))
+  if(!is.matrix(h)) 
+    h <- matrix(h, ncol=length(h))
   
   for(i in 1:(ncol(u)-1)) { # i <- 1
-    l0 <- l0+dCopula(as.matrix(u[,c(1,i+1)]), spCop, h=h[,i], log=T)
-    u0 <- cbind(u0, dduCopula(as.matrix(u[,c(1,i+1)]), spCop, h=h[,i]))
+    l0 <- l0 + dCopula(u[,c(1,i+1)], spCop, h=h[,i], log=T)
+    u0 <- cbind(u0, dduCopula(u[,c(1,i+1)], spCop, h=h[,i]))
   }
-  
+
   l1 <- dCopula(u0, vine, log=T)
   if(log)
     return(l0+l1)
@@ -65,6 +66,8 @@
                                 copula=copula at spCop, h=data at distances[,i]))
   }
   
+  cat(summary(as.data.frame(secLevel)))
+    
   vineCopFit <- fitCopula(copula at vineCop, secLevel, method) 
   
   spVineCop <- spVineCopula(copula at spCop, vineCopFit at copula)
@@ -79,22 +82,29 @@
 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)
+condSpVine <- function (condVar, dists, spVine, n = 1000) {
+  # add some points in the tails
+  rat <- 29:1%x%c(1e-6,1e-5,1e-4,1e-3)
+  xVals <- unique(sort(c(rat, 1 - rat, 1:(n - 1)/(n))))
+  nx <- length(xVals)
   
-  linAppr <- approxfun(c(0,xVals,1), density[c(1,1:xLength,xLength)] ,yleft=0, yright=0)
-  int <- integrate(linAppr,lower=0, upper=1)$value
+  repCondVar <- matrix(condVar, ncol = length(condVar), nrow = nx, byrow = T)
+  density <- dCopula(cbind(xVals, repCondVar), spVine, h = dists)
   
+  # the 1-e7 corners linearily to [0,1], but keep non-negative
+  density <- c(max(0,2*density[1]-density[2]),
+               density, max(0,2*density[nx]-density[nx-1]))
+  linAppr <- approxfun(c(0, xVals, 1), density)
+  
+  # sum up the denstiy to rescale
+  int <- sum(diff(c(0,xVals,1))*(0.5*diff(density)+density[-(nx+2)]))
   return(function(u) linAppr(u)/int)
 }
 
 # interpolation
 
-spCopPredict.expectation <- function(predNeigh, spVine, margin) {
+spCopPredict.expectation <- function(predNeigh, spVine, margin, range) {
+  stopifnot(!is.null(range))
   stopifnot(is.function(margin$d))
   stopifnot(is.function(margin$p))
   
@@ -106,7 +116,7 @@
       condSecVine(margin$p(x))*margin$d(x)*x
     }
     
-    predMean <- c(predMean, integrate(condExp,0,3000,subdivisions=1e6)$value)
+    predMean <- c(predMean, integrate(condExp,range[1],range[2],subdivisions=1e6)$value)
   }
   if ("data" %in% slotNames(predNeigh at locations)) {
     res <- predNeigh at locations
@@ -125,10 +135,16 @@
   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))
+    pPred <- optimise(function(x) abs(integrate(condSecVine, 0, x, 
+                                                subdivisions=10000L, 
+                                                abs.tol=1e-6)$value-p),
+                      c(0,1))
+    if(pPred$objective > 1e-6)
+      warning("Numerical evaluation in predQuantile achieved an obkective of only ",
+              pPred$objective, " where 0 has been sought.")
+    predQuantile <- c(predQuantile, margin$q(pPred$minimum))
   }
+  
   if ("data" %in% slotNames(predNeigh at locations)) {
     res <- predNeigh at locations
     res at data[[paste("quantile.",p,sep="")]] <- predQuantile
@@ -140,8 +156,8 @@
   }
 }
 
-spCopPredict <- function(predNeigh, spVine, margin, method="quantile", p=0.5) {
+spCopPredict <- function(predNeigh, spVine, margin, method="quantile", p=0.5, range=NULL) {
   switch(method,
          quantile=spCopPredict.quantile(predNeigh, spVine, margin, p),
-         expectation=spCopPredict.expectation(predNeigh, spVine, margin))
+         expectation=spCopPredict.expectation(predNeigh, spVine, margin, range))
 }
\ No newline at end of file

Modified: pkg/R/spatialPreparation.R
===================================================================
--- pkg/R/spatialPreparation.R	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/R/spatialPreparation.R	2013-04-23 14:17:11 UTC (rev 94)
@@ -1,15 +1,11 @@
-##################################################################
-##                                                              ##
-## dedicated functions based on sp preparing the use of copulas ##
-##                                                              ##
-##################################################################
+########################################################
+##                                                    ##
+## functions based on sp preparing the use of copulas ##
+##                                                    ##
+########################################################
 
-## neighbourhood
-# constructor
-# data = "matrix"	a matrix or array providing the data
-# distances="matrix"	a matrix providing the distances
-# sp="SpatialPoints"	SpatialPoints object providing the coordinates
-# index="matrix"	linking the obs. in data to the coordinates
+## neighbourhood constructor
+############################
 
 neighbourhood <- function(data, distances, sp, index, prediction, var){
   sizeN <- ncol(distances)+1
@@ -31,7 +27,6 @@
 setMethod(show,signature("neighbourhood"),showNeighbourhood)
 
 ## names (from sp)
-
 setMethod(names, signature("neighbourhood"), namesNeighbourhood <- function(x) x at var)
 
 ## spplot ##
@@ -49,12 +44,7 @@
 ## calculate neighbourhood from SpatialPointsDataFrame
 
 # returns an neighbourhood object
-# 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
+##################################
 
 getNeighbours <- function(spData, locations, var=names(spData)[1], size=5, 
                           prediction=FALSE, min.dist=0.01) {
@@ -94,14 +84,6 @@
                        allLocs, prediction, var))
 }
 
-# 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 ##
 #############
@@ -131,40 +113,46 @@
 }
 
 # the generic calcBins, calculates bins for spatial and spatio-temporal data
+setGeneric("calcBins", function(data, var, nbins=15, boundaries=NA, cutoff=NA,
+                                cor.method="kendall", plot=T, ...) {
+                         standardGeneric("calcBins") 
+                         })
 
-setGeneric("calcBins", function(data, var, nbins=15, boundaries=NA, cutoff=NA, cor.method="kendall", plot=T, ...) standardGeneric("calcBins") )
+## calculating the spatial bins
+################################
 
-# calculating the spatial bins
-# 
-# data denotes the spatial data object
-# var denotes the only variable name used
-# cor.method is passed on to cor() (default="kendall")
-# if plot=TRUE (default), the correlation measures are plotted agaisnt the mean lag separation distance
-# 
-calcSpBins <- function(data, var=names(data), nbins=15, boundaries=NA, cutoff=NA, cor.method="kendall", plot=TRUE) {
+calcSpBins <- function(data, var=names(data), nbins=15, boundaries=NA, 
+                       cutoff=NA, cor.method="kendall", plot=TRUE) {
 
+  if(is.na(cutoff)) {
+    cutoff <- spDists(coordinates(t(data at bbox)))[1,2]/3
+  }
   if(is.na(boundaries)) {
-    diagonal <- spDists(coordinates(t(data at bbox)))[1,2]
-    boundaries <- ((1:nbins) * min(cutoff,diagonal/3,na.rm=T) / nbins)
+    boundaries <- ((1:nbins) * cutoff/nbins)
   }
   
   lags <- calcSpLagInd(data, boundaries)
     
-  mDists <- sapply(lags,function(x) mean(x[,3]))
+  mDists <- sapply(lags, function(x) mean(x[,3]))
+  np <- sapply(lags, function(x) length(x[,3]))
   lagData <- lapply(lags, function(x) as.matrix((cbind(data[x[,1],var]@data, data[x[,2],var]@data))))
   
   if(cor.method == "fasttau")
     lagCor <- sapply(lagData, function(x) VineCopula:::fasttau(x[,1], x[,2]))
-  else 
+  if(cor.method %in% c("kendall","spearman","perarson"))
     lagCor <- sapply(lagData, function(x) cor(x,method=cor.method)[1,2])
-  
+  if(cor.method == "normVariogram")  
+    lagCor <- sapply(lagData, function(x) 1-cor(x,method="pearson")[1,2])
+  if(cor.method == "variogram")  
+    lagCor <- sapply(lagData, function(x) 0.5*mean((x[,1]-x[,2])^2,na.rm=T))
+    
   if(plot) { 
     plot(mDists, lagCor, xlab="distance",ylab=paste("correlation [",cor.method,"]",sep=""), 
          ylim=1.05*c(-abs(min(lagCor)),max(lagCor)), xlim=c(0,max(mDists)))
     abline(h=c(-min(lagCor),0,min(lagCor)),col="grey")
   }
   
-  res <- list(meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=lags)
+  res <- list(np=np, meanDists = mDists, lagCor=lagCor, lagData=lagData, lags=lags)
   attr(res,"cor.method") <- cor.method
   return(res)
 }

Modified: pkg/man/condSpVine.Rd
===================================================================
--- pkg/man/condSpVine.Rd	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/man/condSpVine.Rd	2013-04-23 14:17:11 UTC (rev 94)
@@ -8,7 +8,7 @@
 A spatial vine copula is conditioned under the observations of all but one neighbour generating a conditional uivariate distribution used ofr prediction.
 }
 \usage{
-condSpVine(condVar, dists, spVine, n = 100)
+condSpVine(condVar, dists, spVine, n = 1000)
 }
 \arguments{
   \item{condVar}{

Modified: pkg/man/fitCorFun.Rd
===================================================================
--- pkg/man/fitCorFun.Rd	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/man/fitCorFun.Rd	2013-04-23 14:17:11 UTC (rev 94)
@@ -8,7 +8,7 @@
 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.
 }
 \usage{
-fitCorFun(bins, degree = 3, cutoff = NA, bounds = c(0, 1), cor.method = NULL)
+fitCorFun(bins, degree = 3, cutoff = NA, bounds = c(0, 1), cor.method = NULL, weighted = FALSE)
 }
 \arguments{
   \item{bins}{
@@ -26,6 +26,9 @@
   \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.

Modified: pkg/man/spCopPredict.Rd
===================================================================
--- pkg/man/spCopPredict.Rd	2013-04-15 11:39:18 UTC (rev 93)
+++ pkg/man/spCopPredict.Rd	2013-04-23 14:17:11 UTC (rev 94)
@@ -8,7 +8,7 @@
 A spatial vine copula is used to predict values at unobserved locations conditioned on observations of a local neighbourhood.
 }
 \usage{
-spCopPredict(predNeigh, spVine, margin, method = "quantile", p = 0.5)
+spCopPredict(predNeigh, spVine, margin, method = "quantile", p = 0.5, range = NULL)
 }
 \arguments{
   \item{predNeigh}{
@@ -26,6 +26,9 @@
   \item{p}{
 only used for the quantile predictor indicating the desired fraction the quantile should correspond to.
 }
+  \item{range}{
+  the range of integartion to be used in the numerical integration in \code{"expectation"} by \code{\link{integrate}}.
+  }
 }
 \details{
 Predictions are done based on \code{\link{condSpVine}} through numerical integration/optimisation.

Modified: spcopula_0.1-1.tar.gz
===================================================================
(Binary files differ)

Modified: spcopula_0.1-1.zip
===================================================================
(Binary files differ)



More information about the spcopula-commits mailing list