[spcopula-commits] r112 - pkg/R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Nov 19 09:57:27 CET 2013


Author: ben_graeler
Date: 2013-11-19 09:57:27 +0100 (Tue, 19 Nov 2013)
New Revision: 112

Modified:
   pkg/R/Classes.R
   pkg/R/cqsCopula.R
   pkg/R/linkingVineCopula.R
   pkg/R/partialDerivatives.R
   pkg/R/spCopula.R
   pkg/R/spVineCopula.R
   pkg/R/utilities.R
Log:
- correction in partial derivative of cqsCopula 
- typo in spCopula for single distances

Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R	2013-11-11 17:05:37 UTC (rev 111)
+++ pkg/R/Classes.R	2013-11-19 08:57:27 UTC (rev 112)
@@ -129,7 +129,7 @@
     nonIndep <- sapply(object at components[-nComp], function(x) class(x) != "indepCopula")
     for (i in (1:(nComp-1))[nonIndep]) {
       upParam <- object at calibMoa(object at components[[i]], object at distances[i+1])
-      if(is.na(upParam)) {
+      if(any(is.na(upParam))) {
         check.upper <- c(check.upper, TRUE)
       } else {
         if (class(object at components[[i]]) == "frankCopula" && upParam == 0) {

Modified: pkg/R/cqsCopula.R
===================================================================
--- pkg/R/cqsCopula.R	2013-11-11 17:05:37 UTC (rev 111)
+++ pkg/R/cqsCopula.R	2013-11-19 08:57:27 UTC (rev 112)
@@ -131,15 +131,13 @@
   c0 <- -y
 
   v <- solveCubicEq(c3,c2,c1,c0)
-
-  return(v)
   
-#   filter <- function(vec){
-#     vec <- vec[!is.na(vec)]
-#     return(vec[vec >= 0 & vec <= 1])
-#   }
-# 
-#   return(apply(v,1,filter))
+  filter <- function(vec){
+    vec <- vec[!is.na(vec)]
+    return(vec[vec >= 0 & vec <= 1])
+  }
+
+  return(apply(v,1,filter))
 }
 
 setMethod("invdduCopula", signature("numeric","cqsCopula","numeric"), invdduCQSec)
@@ -180,14 +178,12 @@
 
   u <- solveCubicEq(c3,c2,c1,c0)
   
-  return(u)
+  filter <- function(vec){
+    vec <- vec[!is.na(vec)]
+    return(vec[vec >= 0 & vec <= 1])
+  }
 
-#   filter <- function(vec){
-#     vec <- vec[!is.na(vec)]
-#     return(vec[vec >= 0 & vec <= 1])
-#   }
-# 
-#   return(apply(u,1,filter))
+  return(apply(u,1,filter))
 }
 
 setMethod("invddvCopula", signature("numeric","cqsCopula","numeric"), invddvCQSec)
@@ -198,7 +194,7 @@
   u <- runif(n, min = 0, max = 1)
   y <- runif(n, min = 0, max = 1)
     
-  res <- cbind(u, invdduCQSec(u, copula, y))
+  res <- cbind(u, spcopula:::invdduCQSec(u, copula, y))
   colnames(res) <- c("u","v")
     
   return(res)
@@ -304,7 +300,7 @@
 fitCQSec.ml <- function(copula, data, start, lower, upper, optim.control, optim.method) { 
   if(length(start)!=2) stop("Start values need to have same length as parameters.")
   
-  if (length(copula at fixed)==0) {
+  if (copula at fixed=="") {
     optFun <- function(param=c(0,0)) {
       if(any(param > 1) | param[2] < -1 | param[1] < limA(param[2])) 
         return(100)

Modified: pkg/R/linkingVineCopula.R
===================================================================
--- pkg/R/linkingVineCopula.R	2013-11-11 17:05:37 UTC (rev 111)
+++ pkg/R/linkingVineCopula.R	2013-11-19 08:57:27 UTC (rev 112)
@@ -19,14 +19,16 @@
   constr(c(par,par2))
 }
 
-#####################################################
+#########################################################
 ## generic wrapper functions to the VineCopula package ##
-#####################################################
+#########################################################
 
 # density from BiCopPDF
 linkVineCop.PDF <- function (u, copula, log=FALSE) {
   param <- copula at parameters
-  if(length(param)==1) param <- c(param,0)
+
+  if(length(param)==1) 
+    param <- c(param,0)
   n <- nrow(u)
   fam <- copula at family
 
@@ -34,8 +36,10 @@
 #   coplik = .C("LL_mod_seperate", as.integer(fam), as.integer(n), as.double(u[,1]), 
 #               as.double(u[,2]), as.double(param[1]), as.double(param[2]), 
 #               as.double(rep(0, n)), PACKAGE = "VineCopula")[[7]]
-  if(log) return(coplik)
-  else return(exp(coplik))
+  if(log) 
+    return(coplik)
+  else 
+    return(exp(coplik))
 }
 
 # cdf from BiCopCDF

Modified: pkg/R/partialDerivatives.R
===================================================================
--- pkg/R/partialDerivatives.R	2013-11-11 17:05:37 UTC (rev 111)
+++ pkg/R/partialDerivatives.R	2013-11-19 08:57:27 UTC (rev 112)
@@ -101,9 +101,23 @@
 ## independent copula ##
 ########################
 
+## Kendall's tau
+################
+setMethod("tau", signature("indepCopula"), function(copula, ...) return(0))
+
+## Spearman's rho
+#################
+setMethod("rho", signature("indepCopula"), function(copula, ...) return(0))
+
+## indepCopula as evCopula derivatives of A
+###########################################
+setMethod("dAdu", signature("indepCopula"), 
+          function(copula, w) {
+            data.frame(der1=rep(0, length(w)), der2=rep(0, length(w)))
+          })
+
 ## partial derivative d/du
 ##########################
-
 setMethod("dduCopula", signature("numeric","indepCopula"),
           function(u, copula, ...) {
             matrix(u,ncol=copula at dimension)[,2]
@@ -112,7 +126,6 @@
 
 ## inverse of the partial derivative d/du
 #########################################
-
 invdduIndep <- function(u, copula, y){
   return(y)
 }
@@ -121,7 +134,6 @@
 
 ## partial derivative d/dv
 ##########################
-
 setMethod("ddvCopula", signature("numeric","indepCopula"),
           function(u, copula, ...) {
             matrix(u,ncol=copula at dimension)[,1]
@@ -130,7 +142,6 @@
 
 ## inverse of the partial derivative d/dv
 #########################################
-
 invddvIndep <- function(v, copula, y){
   return(y)
 }
@@ -144,7 +155,6 @@
 
 ## partial derivative d/du
 ##########################
-
 dduClayton <- function(u, copula){
   rho <- copula at parameters
 

Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R	2013-11-11 17:05:37 UTC (rev 111)
+++ pkg/R/spCopula.R	2013-11-19 08:57:27 UTC (rev 112)
@@ -73,10 +73,10 @@
     cmpCop <- object at components[[i]]
     cat("  ", cmpCop at fullname, "at", object at distances[i], 
       paste("[",object at unit,"]",sep=""), "\n")
-    if (length(cmpCop at parameters) > 0) {
-      for (i in (1:length(cmpCop at parameters))) 
-        cat("    ", cmpCop at param.names[i], " = ", cmpCop at parameters[i], "\n")
-    }
+#     if (length(cmpCop at parameters) > 0) {
+#       for (i in (1:length(cmpCop at parameters))) 
+#         cat("    ", cmpCop at param.names[i], " = ", cmpCop at parameters[i], "\n")
+#     }
   }
   if(!is.null(object at calibMoa(normalCopula(0),0))) cat("A spatial dependence function is used. \n")
 }
@@ -213,7 +213,7 @@
       } else {
         if(class(lowerCop) != "indepCopula") {
           lowerParam <- calibPar(lowerCop, h)
-          lowerCop at parameters[length(lowerParam)] <- lowerParam
+          lowerCop at parameters[1:length(lowerParam)] <- lowerParam
         }
         return(fun(pairs, lowerCop, ...))
       }
@@ -499,6 +499,7 @@
                         spearman=fitASC2.irho(cop, bins$lagData[[i]],
                                               rho=calcCor(bins$meanDists[i]))@copula,
                         stop(paste(calcCor(NULL), "is not yet supported.")))
+          param <- cop at parameters
         } else {
           if(class(cop) == "cqsCopula") {
             cop <- switch(calcCor(NULL),
@@ -507,6 +508,7 @@
                           spearman=fitCQSec.irho(cop, bins$lagData[[i]],
                                                 rho=calcCor(bins$meanDists[i]))@copula,
                           stop(paste(calcCor(NULL), "is not yet supported.")))
+            param <- cop at parameters
           } else {
             param <- moa(cop, bins$meanDists[i])
             if(!is.na(param))
@@ -515,7 +517,7 @@
         }
       }
       
-      if(is.na(param))
+      if(any(is.na(param)))
         tmploglik <- c(tmploglik, NA)
       else 
         tmploglik <- c(tmploglik, sum(dCopula(bins$lagData[[i]], cop, log=T)))

Modified: pkg/R/spVineCopula.R
===================================================================
--- pkg/R/spVineCopula.R	2013-11-11 17:05:37 UTC (rev 111)
+++ pkg/R/spVineCopula.R	2013-11-19 08:57:27 UTC (rev 112)
@@ -186,16 +186,15 @@
 # interpolation
 
 spCopPredict.expectation <- function(predNeigh, spVine, margin, ..., stop.on.error=F) {
-#   stopifnot(!is.null(range))
-#   stopifnot(is.function(margin$d))
-#   stopifnot(is.function(margin$p))
   stopifnot(is.function(margin$q))
   
   dists <- calcSpTreeDists(predNeigh,length(spVine at spCop))
   
   predMean <- NULL
+  
+  pb <- txtProgressBar(0,nrow(predNeigh at data), 0, width=getOption("width")-10, style=3)
   for(i in 1:nrow(predNeigh at data)) { # i <-1
-    cat("[Predicting location ",i,".]\n", sep="")
+    setTxtProgressBar(pb, i)
     condSecVine <- condSpVine(as.numeric(predNeigh at data[i,-1]), 
                               lapply(dists,function(x) x[i,]), spVine)
     
@@ -203,12 +202,14 @@
       margin$q(x)*condSecVine(x)
     }
     
-    ePred <- integrate(condExp,0,1,subdivisions=10000L,stop.on.error=stop.on.error, ...)
+    ePred <- integrate(condExp,0+.Machine$double.eps,1-.Machine$double.eps,subdivisions=10000L,stop.on.error=stop.on.error, ...)
     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)
   }
+  close(pb)
+  
   if ("data" %in% slotNames(predNeigh at predLocs)) {
     res <- predNeigh at predLocs
     res at data[["expect"]] <- predMean
@@ -225,8 +226,9 @@
   dists <- calcSpTreeDists(predNeigh,length(spVine at spCop))
   
   predQuantile <- NULL
+  pb <- txtProgressBar(0, nrow(predNeigh at data), 0, width=getOption("width")-10, style=3)
   for(i in 1:nrow(predNeigh at data)) { # i <-1
-    cat("[Predicting location ",i,".]\n", sep="")
+    setTxtProgressBar(pb, i)
     condSecVine <- condSpVine(as.numeric(predNeigh at data[i,-1]), 
                               lapply(dists,function(x) x[i,]), spVine)
     
@@ -247,6 +249,7 @@
 #               pPred$objective, " where 0 has been sought for location ",i,".")
     predQuantile <- c(predQuantile, margin$q(xVals[lower]+xRes))
   }
+  close(pb)
   
   if ("data" %in% slotNames(predNeigh at predLocs)) {
     res <- predNeigh at predLocs

Modified: pkg/R/utilities.R
===================================================================
--- pkg/R/utilities.R	2013-11-11 17:05:37 UTC (rev 111)
+++ pkg/R/utilities.R	2013-11-19 08:57:27 UTC (rev 112)
@@ -2,7 +2,7 @@
 
 
 # ranks are automatically removed and NAs are by default randomly distributed
-rankTransform <- function(u,v=NULL, ties.method="average") {
+rankTransform <- function(u,v=NULL, na.last=TRUE, ties.method="average") {
   if(!(is.matrix(u) | is.data.frame(u))) {
     if (is.null(v))
       stop("u must either be a matrix with at least 2 columns or u and v must be given.")
@@ -10,8 +10,7 @@
       u <- cbind(u,v)
   }
 
-  bool <- apply(u,1,function(row) !any(is.na(row)))
-  res <- apply(u[bool,],2,rank,ties.method=ties.method)/(sum(bool)+1)
+  res <- apply(u,2,function(x) rank(x,na.last,ties.method)/(sum(!is.na(x))+1))
   if(is.data.frame(u))
     return(as.data.frame(res))
   return(res)



More information about the spcopula-commits mailing list