[spcopula-commits] r162 - / pkg pkg/R pkg/man pkg/tests/Examples

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Feb 7 12:35:01 CET 2017


Author: ben_graeler
Date: 2017-02-07 12:35:01 +0100 (Tue, 07 Feb 2017)
New Revision: 162

Added:
   pkg/R/hkCopula.R
   pkg/man/hkCopula-class.Rd
   pkg/man/hkCopula.Rd
   pkg/man/rCopula_y.Rd
Modified:
   /
   pkg/DESCRIPTION
   pkg/NAMESPACE
   pkg/R/Classes.R
   pkg/R/asCopula.R
   pkg/R/cqsCopula.R
   pkg/R/empiricalCopula.R
   pkg/R/mixtureCopula.R
   pkg/R/partialDerivatives.R
   pkg/R/spCopula.R
   pkg/R/spVineCopula.R
   pkg/R/stCoVarVineCopula.R
   pkg/R/stCopula.R
   pkg/R/stVineCopula.R
   pkg/man/asCopula.Rd
   pkg/man/cqsCopula.Rd
   pkg/man/mixtureCopula.Rd
   pkg/tests/Examples/spcopula-Ex.Rout.save
Log:
keeping track with copula, adds some functionality for Hierarchical Kendall Copulas


Property changes on: 
___________________________________________________________________
Modified: svn:ignore
   - .Rd2pdf6724
.Rhistory
.Rproj.user
Meuse_spcopula_estimation.R
pkg.Rcheck
pkg.pdf
spcopula.Rcheck
spcopula.Rproj
spCopDemo_old_Mar_2013.RData
.RData
inst
spcopula_0.2-1.tar.gz

   + .Rd2pdf6724
.Rhistory
.Rproj.user
Meuse_spcopula_estimation.R
pkg.Rcheck
pkg.pdf
spcopula.Rcheck
spcopula.Rproj
spCopDemo_old_Mar_2013.RData
.RData
inst
spcopula_0.2-1.tar.gz
copula_example_template.R
spcopula_0.2-2.tar.gz


Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/DESCRIPTION	2017-02-07 11:35:01 UTC (rev 162)
@@ -1,13 +1,13 @@
 Package: spcopula
 Type: Package
-Title: Copula Driven Spatio-Temporal Analysis
-Version: 0.2-1
-Date: 2016-09-12
+Title: Copula Driven Analysis - Multivariate, Spatial, Spatio-Temporal
+Version: 0.2-2
+Date: 2017-02-07
 Authors at R: c(person("Benedikt", "Graeler", role = c("aut", "cre"),
-                    email = "ben.graeler at uni-muenster.de"),
+                    email = "b.graeler at 52north.org"),
              person("Marius", "Appel",role = "ctb"))
-Maintainer: Benedikt Graeler <ben.graeler at uni-muenster.de>
-Description: We provide a framework to analyse spatial and spatio-temporal data via copulas and vine copulas. The data needs to be provided in the form of the sp and spacetime package respectively. Additionally, support for calculating different multivariate return periods based on copulas and vine copulas is implemented.
+Maintainer: Benedikt Graeler <b.graeler at 52north.org>
+Description: A framework to analyse multivariate, spatial and spatio-temporal data via copulas, vine copulas and hierarchical Kendall copulas is provided. The spatial/spatio-temporal data needs to be provided in the form of the sp and spacetime package respectively. Additionally, support for calculating different multivariate return periods based on copulas and vine copulas is implemented. As of version 0.2-2, Hierarchical Kendall Copulas are also provided.
 License: GPL-3
 LazyLoad: yes
 Depends: copula (>= 0.999-15), R (>= 3.1.0), VineCopula (>= 2.0.4)
@@ -35,4 +35,5 @@
   tawn3pCopula.R
   tailDependenceFunctions.R
   KendallDistribution.R
+  hkCopula.R
   zzz.R

Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/NAMESPACE	2017-02-07 11:35:01 UTC (rev 162)
@@ -4,10 +4,10 @@
 import(methods)
 
 importFrom("graphics", "abline", "smoothScatter")
-  importFrom("stats", "D", "approxfun", "cor", "ecdf", "integrate", "lm",
-             "optim", "optimise", "optimize", "pnorm", "predict", "pt",
-             "qnorm", "qt", "quantile", "runif", "uniroot", "var")
-  importFrom("utils", "setTxtProgressBar", "txtProgressBar", "tail")
+importFrom("stats", "D", "approxfun", "cor", "ecdf", "integrate", "lm",
+           "optim", "optimise", "optimize", "pnorm", "predict", "pt",
+           "qnorm", "qt", "quantile", "runif", "uniroot", "var")
+importFrom("utils", "setTxtProgressBar", "txtProgressBar", "tail", "capture.output", "str")
 
 importMethodsFrom(VineCopula, fitCopula)
 importMethodsFrom(VineCopula, dduCopula,ddvCopula)
@@ -46,13 +46,14 @@
 export(neighbourhood, stNeighbourhood)
 export(empiricalCopula, genEmpCop, empSurCopula, genEmpSurCop)
 export(mixtureCopula)
+export(hkCopula)
 
 # general functions
 export(rankTransform, dependencePlot, unitScatter, univScatter)
 export(fitCopula)
 export(dduCopula,ddvCopula)
 export(invdduCopula, invddvCopula)
-export(qCopula_u, qCopula_v)
+export(qCopula_u, qCopula_v, rCopula_y)
 export(condSpVine,spCopPredict)
 export(condStVine,stCopPredict)
 export(condStCoVarVine, condCovariate)
@@ -82,4 +83,5 @@
 exportClasses(asCopula, cqsCopula, tawn3pCopula, neighbourhood, stNeighbourhood, empiricalCopula, empSurCopula)
 exportClasses(spCopula, stCopula, spVineCopula, stVineCopula)
 exportClasses(stCoVarVineCopula)
-exportClasses(mixtureCopula)
\ No newline at end of file
+exportClasses(mixtureCopula)
+exportClasses(hkCopula)
\ No newline at end of file

Modified: pkg/R/Classes.R
===================================================================
--- pkg/R/Classes.R	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/R/Classes.R	2017-02-07 11:35:01 UTC (rev 162)
@@ -1,87 +1,4 @@
-## an asymmetric copula with cubic and quadratic sections
-
-validAsCopula = function(object) {
-  if (object at dimension != 2)
-    return("Only copulas with cubic quadratic sections of dimension 2 are supported.")
-  param <- object at parameters
-  upper <- object at param.upbnd
-  lower <- object at param.lowbnd
-  if (length(param) != length(upper))
-    return("Parameter and upper bound have non-equal length")
-  if (length(param) != length(lower))
-    return("Parameter and lower bound have non-equal length")
-  if (any(is.na(param) | param > upper | param < lower))
-    return("Parameter value out of bound")
-  else return (TRUE)
-}
-
-# the lower bound of the parameter a dependening on the parameter b
-limA <- function (b) {
-  stopifnot(abs(b) <= 1)
-  0.5*(-sqrt(-3*b^2+6*b+9)+b-3)
-}
-
-# the lower and upper bound of the parameter b dependening on the parameter a
-limB <- function (a) {
-  stopifnot(a <=1 & a >= -3)
-  if(a>-2)
-    return(c(-1,1))
-  pmax(pmin(0.5*(c(-1,1)*(sqrt(3)*sqrt(-a^2-2*a+3))+a+3),1),-1)
-}
-
-setClass("asCopula",
-  representation = representation("copula"),
-  validity = validAsCopula,
-  contains = list("copula")
-)
-
 ####
-## a symmetric copula with cubic and quadratic sections
-
-validCqsCopula <- function(object) {
-  if (object at dimension != 2)
-    return("Only copulas with cubic quadratic sections of dimension 2 are supported.")
-  param <- object at parameters
-  upper <- object at param.upbnd
-  lower <- object at param.lowbnd
-  if (length(param) != length(upper))
-    return("Parameter and upper bound have non-equal length")
-  if (length(param) != length(lower))
-    return("Parameter and lower bound have non-equal length")
-  if (any(is.na(param) | param > upper | param < lower))
-    return("Parameter value out of bound")
-  if (object at fixed != ""){
-    if(!("a" %in% object at fixed | "b" %in% object at fixed))
-      return("The slot fixed may only refer to \"a\" or \"b\".")
-    if ("a" %in% object at fixed & "b" %in% object at fixed)
-      return("Only one of the parameters may be kept fixed.")
-  }
-  else return (TRUE)
-}
-
-setClass("cqsCopula",
-  representation = representation("copula",fixed="character"),
-  validity = validCqsCopula,
-  contains = list("copula")
-)
-
-####
-## an empirical copula representation
-
-validEmpCopula <- function(object) {
-  if(ncol(object at sample) != object at dimension)
-    return("Dimension of the copula and the sample do not match.")
-  else
-    return(TRUE)
-}
-
-setClass("empiricalCopula",
-         representation = representation("copula", sample="matrix"),
-         validity = validEmpCopula,
-         contains = list("copula")
-)
-
-####
 ## an empirical survival copula representation
 
 validEmpSurCopula <- function(object) {
@@ -115,6 +32,32 @@
          contains = list("copula")
 )
 
+#####
+## Hierarchical Kendall Copulas
+
+validHkCopula <- function(object) {
+  stopifnot(all(sapply(object at clusterCops, function(x) inherits(x[[1]], "copula"))))
+  if (object at dimension != sum(sapply(object at clusterCops, function(x) x[[1]]@dimension))+object at nestingCop@dimension-length(object at clusterCops))
+    return("The dimensions of the hierarchical Kendall copula do not match.")
+  
+  if(length(object at clusterCops) != length(object at kenFuns))
+    return("Each cluster copula needs to have its Kendall function in 'kenFuns'.")
+  
+  else return (TRUE)
+}
+
+setClass("hkCopula",
+         representation = representation("copula",
+                                         nestingCop = "copula",
+                                         clusterCops = "list",
+                                         kenFuns = "list"),
+         validity = validHkCopula,
+         contains = list("copula")
+)
+
+
+
+
 ## 
 ## the spatial copula
 ##
@@ -158,10 +101,10 @@
       check.lower <- c(check.lower, is.na(object at calibMoa(object at components[[i]], c(0,object at distances)[i])))
     }
     if(sum(check.upper>0)) return(paste("Reconsider the upper boundary conditions of the following copula(s): \n",
-                                        paste(sapply(object at components[check.upper], function(x) x at fullname), 
+                                        paste(sapply(object at components[check.upper], function(x) describeCop(x, "very short")), 
                                               "at", object at distances[check.upper],collapse="\n")))
     if(sum(check.lower>0)) return(paste("Reconsider the lower boundary conditions of the following copula(s): \n",
-                                        paste(sapply(object at components[check.lower], function(x) x at fullname), 
+                                        paste(sapply(object at components[check.lower], function(x) describeCop(x, "very short")), 
                                               "at", object at distances[check.lower],collapse="\n")))
   }
   return(TRUE)
@@ -280,4 +223,5 @@
                                          var="character", 
                                          coVar="character",
                                          prediction="logical"),
-         validity = validStNeighbourhood)
\ No newline at end of file
+         validity = validStNeighbourhood)
+

Modified: pkg/R/asCopula.R
===================================================================
--- pkg/R/asCopula.R	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/R/asCopula.R	2017-02-07 11:35:01 UTC (rev 162)
@@ -5,6 +5,41 @@
 ##########################
 # (see Example 3.16 in: Nelsen, Roger B. (2006): An Introduction to Copulas, second edition, Springer)
 
+validAsCopula = function(object) {
+  if (object at dimension != 2)
+    return("Only copulas with cubic quadratic sections of dimension 2 are supported.")
+  param <- object at parameters
+  upper <- object at param.upbnd
+  lower <- object at param.lowbnd
+  if (length(param) != length(upper))
+    return("Parameter and upper bound have non-equal length")
+  if (length(param) != length(lower))
+    return("Parameter and lower bound have non-equal length")
+  if (any(is.na(param) | param > upper | param < lower))
+    return("Parameter value out of bound")
+  else return (TRUE)
+}
+
+# the lower bound of the parameter a dependening on the parameter b
+limA <- function (b) {
+  stopifnot(abs(b) <= 1)
+  0.5*(-sqrt(-3*b^2+6*b+9)+b-3)
+}
+
+# the lower and upper bound of the parameter b dependening on the parameter a
+limB <- function (a) {
+  stopifnot(a <=1 & a >= -3)
+  if(a>-2)
+    return(c(-1,1))
+  pmax(pmin(0.5*(c(-1,1)*(sqrt(3)*sqrt(-a^2-2*a+3))+a+3),1),-1)
+}
+
+setClass("asCopula",
+         representation = representation("copula"),
+         validity = validAsCopula,
+         contains = list("copula")
+)
+
 # constructor
 asCopula <- function (param=c(0,0)) {
   val <- new("asCopula", dimension = as.integer(2), parameters = param, 
@@ -13,6 +48,24 @@
   return(val)
 }
 
+## printing
+setMethod("describeCop", c("asCopula", "character"),
+          function(x, kind = c("short", "very short", "long"), prefix = "", ...) {
+            kind <- match.arg(kind)
+            if(kind == "very short") # e.g. for show() which has more parts
+              return(paste0(prefix, "AS-CQS copula"))
+            
+            name <- "asymmetric"
+            d <- dim(x)
+            ch <- paste0(prefix, name, " copula, dim. d = ", d)
+            switch(kind <- match.arg(kind),
+                   short = ch,
+                   long = paste0(ch, "\n", prefix, " param.: ",
+                                 capture.output(str(x at parameters,
+                                                    give.head=FALSE))),
+                   stop("invalid 'kind': ", kind))
+          })
+
 ## density ##
 
 dASC2 <- function (u, copula, log=FALSE) {
@@ -184,11 +237,14 @@
 fitCopulaASC2 <- function (copula, data, method = "ml", start=c(0,0),
                            lower=c(-3,-1), upper=c(1,1), 
                            optim.method="L-BFGS-B", optim.control=list(),
-                           estimate.variance = FALSE) {
+                           estimate.variance = FALSE, call) {
+  if(missing(call)) 
+    call <- match.call()
+
   fit <- switch(method, 
-                ml=fitASC2.ml(copula, data, start, lower, upper, optim.control, optim.method),
-                itau=fitASC2.itau(copula, data, estimate.variance),
-                irho=fitASC2.irho(copula, data, estimate.variance),
+                ml=fitASC2.ml(copula, data, start, lower, upper, optim.control, optim.method, call),
+                itau=fitASC2.itau(copula, data, estimate.variance, call),
+                irho=fitASC2.irho(copula, data, estimate.variance, call),
                 stop("Implemented methods for copulas in the spCopula package are: ml, itau, and irho."))
   return(fit)
 }
@@ -207,27 +263,33 @@
 # method
 #  one of kendall or spearman according to the calculation of moa
 
-fitASC2.itau <- function(copula, data, estimate.variance, tau=NULL) {
+fitASC2.itau <- function(copula, data, estimate.variance, tau=NULL, call) {
+  if(missing(call)) 
+    call <- match.call()
   if(is.null(tau))
     tau <- TauMatrix(data)[1,2]
+  
   esti <- fitASC2.moa(tau, data, method="itau")
   copula <- asCopula(esti)
   
   new("fitCopula", estimate = esti, var.est = matrix(NA), 
       loglik = sum(log(dCopula(data, copula))), nsample = nrow(data),
-      method = "Inversion of Kendall's tau and MLE", 
+      method = "Inversion of Kendall's tau and MLE", call = call,
       fitting.stats = list(convergence=as.integer(NA)), copula = copula)
 }
 
 fitASC2.irho <- function(copula, data, estimate.variance, rho=NULL){
+  if(missing(call)) 
+    call <- match.call()
   if(is.null(rho))
     rho <- cor(data,method="spearman")[1,2]
+  
   esti <- fitASC2.moa(rho, data, method="irho")
   copula <- asCopula(esti)
   
   new("fitCopula", estimate = esti, var.est = matrix(NA), 
       loglik = sum(log(dCopula(data, copula))), nsample = nrow(data),
-      method = "Inversion of Spearman's rho and MLE",
+      method = "Inversion of Spearman's rho and MLE", call=call,
       fitting.stats = list(convergence=as.integer(NA)), copula = copula)
 }
 
@@ -252,8 +314,11 @@
 
 # maximum log-likelihood estimation of a and b using optim
 
-fitASC2.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:")
+fitASC2.ml <- function(copula, data, start, lower, upper, optim.control, optim.method, call) { 
+  if(missing(call)) 
+    call <- match.call()
+  if(length(start)!=2) 
+    stop("Start values need to have same length as parameters:")
   
   optFun <- function(param=c(0,0)) {
     if(any(param > 1) | param[2] < -1 | param[1] < limA(param[2])) return(1)
@@ -269,6 +334,7 @@
              loglik = -optimized$value,
              nsample = nrow(data),
              method = "Numerical MLE over the full range.", 
+             call = call,
              fitting.stats = optimized,
              copula = asCopula(optimized$par)))
 }

Modified: pkg/R/cqsCopula.R
===================================================================
--- pkg/R/cqsCopula.R	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/R/cqsCopula.R	2017-02-07 11:35:01 UTC (rev 162)
@@ -3,7 +3,36 @@
 ## a symmetric copula with cubic quadratic sections ##
 ##                                                  ##
 ######################################################
+# (see Example 3.16 in: Nelsen, Roger B. (2006): An Introduction to Copulas, second edition, Springer)
 
+validCqsCopula <- function(object) {
+  if (object at dimension != 2)
+    return("Only copulas with cubic quadratic sections of dimension 2 are supported.")
+  param <- object at parameters
+  upper <- object at param.upbnd
+  lower <- object at param.lowbnd
+  if (length(param) != length(upper))
+    return("Parameter and upper bound have non-equal length")
+  if (length(param) != length(lower))
+    return("Parameter and lower bound have non-equal length")
+  if (any(is.na(param) | param > upper | param < lower))
+    return("Parameter value out of bound")
+  if (object at fixed != ""){
+    if(!("a" %in% object at fixed | "b" %in% object at fixed))
+      return("The slot fixed may only refer to \"a\" or \"b\".")
+    if ("a" %in% object at fixed & "b" %in% object at fixed)
+      return("Only one of the parameters may be kept fixed.")
+  }
+  else return (TRUE)
+}
+
+setClass("cqsCopula",
+         representation = representation("copula",fixed="character"),
+         validity = validCqsCopula,
+         contains = list("copula")
+)
+
+# constructor 
 cqsCopula <- function (param=c(0,0), fixed="") {
   new("cqsCopula", dimension = as.integer(2), parameters = param, 
       param.names = c("a", "b"), param.lowbnd = c(limA(param[2]),-1),
@@ -11,6 +40,24 @@
       fullname = "copula family with cubic quadratic sections", fixed=fixed)
 }
 
+## printing
+setMethod("describeCop", c("cqsCopula", "character"),
+          function(x, kind = c("short", "very short", "long"), prefix = "", ...) {
+            kind <- match.arg(kind)
+            if(kind == "very short") # e.g. for show() which has more parts
+              return(paste0(prefix, "CQS copula"))
+            
+            name <- "cubic-quadratic sections"
+            d <- dim(x)
+            ch <- paste0(prefix, name, " copula, dim. d = ", d)
+            switch(kind <- match.arg(kind),
+                   short = ch,
+                   long = paste0(ch, "\n", prefix, " param.: ",
+                                 capture.output(str(x at parameters,
+                                                    give.head=FALSE))),
+                   stop("invalid 'kind': ", kind))
+          })
+
 ## density ##
 dCQSec <- function (u, copula, log=F) {
   a <- copula at parameters[1]

Modified: pkg/R/empiricalCopula.R
===================================================================
--- pkg/R/empiricalCopula.R	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/R/empiricalCopula.R	2017-02-07 11:35:01 UTC (rev 162)
@@ -4,6 +4,21 @@
 ##                                    ##
 ########################################
 
+# validity
+validEmpCopula <- function(object) {
+  if(ncol(object at sample) != object at dimension)
+    return("Dimension of the copula and the sample do not match.")
+  else
+    return(TRUE)
+}
+
+# class definition
+setClass("empiricalCopula",
+         representation = representation("copula", sample="matrix"),
+         validity = validEmpCopula,
+         contains = list("copula")
+)
+
 # constructor
 empiricalCopula <- function (sample=NULL, copula) {
   if(is.null(sample) && missing(copula))
@@ -22,7 +37,7 @@
   new("empiricalCopula", dimension = copula at dimension, 
       parameters = copula at parameters, param.names = copula at param.names, 
       param.lowbnd = copula at param.lowbnd, param.upbnd = copula at param.upbnd, 
-      fullname = paste("Empirical copula derived from",copula at fullname),
+      fullname = paste("Empirical copula derived from", describeCop(copula, "very short")),
       sample=sample)
 }
 
@@ -33,6 +48,24 @@
   empiricalCopula(rCopula(sample.size, copula), copula)
 }
 
+# printing
+setMethod("describeCop", c("empiricalCopula", "character"),
+          function(x, kind = c("short", "very short", "long"), prefix = "", ...) {
+            kind <- match.arg(kind)
+            name <- "empirical"
+            if(kind == "very short") # e.g. for show() which has more parts
+              return(paste0(prefix, name, " copula"))
+            ## else
+            d <- dim(x)
+            ch <- paste0(prefix, name, " copula, dim. d = ", d)
+            switch(kind <- match.arg(kind),
+                   short = ch,
+                   long = paste0(ch, "\n", prefix, " param.: ",
+                                 capture.output(str(x at parameters,
+                                                    give.head=FALSE))),
+                   stop("invalid 'kind': ", kind))
+          })
+
 ## density, not yet needed and hence not implemented ##
 
 ## jcdf ##
@@ -88,7 +121,7 @@
   new("empSurCopula", dimension = copula at dimension, 
       parameters = copula at parameters, param.names = copula at param.names, 
       param.lowbnd = copula at param.lowbnd, param.upbnd = copula at param.upbnd, 
-      fullname = paste("Empirical survival copula derived from",copula at fullname),
+      fullname = paste("Empirical survival copula derived from", describeCop(copula, "very short")),
       sample=sample)
 }
 

Added: pkg/R/hkCopula.R
===================================================================
--- pkg/R/hkCopula.R	                        (rev 0)
+++ pkg/R/hkCopula.R	2017-02-07 11:35:01 UTC (rev 162)
@@ -0,0 +1,232 @@
+## Hierarchical Kendall Copulas as defined in Brechmann, Eike Christian. 
+## "Hierarchical Kendall copulas: Properties and inference." Canadian Journal 
+## of Statistics 42.1 (2014): 78-108.
+
+# Slots:
+#   
+# Name:     nestingCop   clusterCops   kenFuns   dimension   parameters  param.names param.lowbnd  param.upbnd     fullname
+# Class:        copula          list      list     integer      numeric    character      numeric      numeric    character
+
+## nestingCop copula
+## clusterCop list of list (copula, ind)
+
+# easy constructor
+
+hkCopula <- function(nestingCop, clusterCops, kenFuns=NULL) {
+  if (is.null(kenFuns)) {
+    kenFuns <- lapply(clusterCops, function(copInd) getKendallDistr(copInd[[1]]))
+  }
+  
+  new("hkCopula", 
+      nestingCop = nestingCop, 
+      clusterCops = clusterCops,
+      kenFuns = kenFuns,
+      dimension = as.integer(sum(sapply(clusterCops, function(x) x[[1]]@dimension))+nestingCop at dimension-length(clusterCops)),
+      parameters = NA_real_,
+      param.names =NA_character_,
+      param.lowbnd = NA_real_,
+      param.upbnd = NA_real_, 
+      fullname = "Hierarchical Kendall Copula")
+}
+
+showHkCopula <- function(object) {
+  cat(object at fullname, "\n")
+  cat("Dimension: ", object at dimension, "\n")
+  cat("Nesting copula:\n")
+  show(object at nestingCop)
+  cat("Cluster copulas:\n")
+  for (i in 1:length(object at clusterCops)) {
+    cmpCop <- object at clusterCops[[i]][[1]]
+    cat("  ", describeCop(cmpCop, "very short"), 
+        "of dimension", cmpCop at dimension,
+        "for indices", object at clusterCops[[i]][[2]], "\n")
+  }
+}
+
+setMethod("show", signature("hkCopula"), showHkCopula)
+
+## density
+
+dHkCop  <- function(u, copula, log=F, ...) {
+  stopifnot(ncol(u) == copula at dimension)
+
+  lik <- NULL
+  kenVal <- NULL
+  
+  for (i in 1:length(copula at clusterCops)) {
+    cop <- copula at clusterCops[[i]][[1]]
+    ind <- copula at clusterCops[[i]][[2]]
+    ken <- copula at kenFuns[[i]]
+    
+    lik <- cbind(lik, dCopula(u[, ind], cop, log=log))
+    kenVal <- cbind(kenVal, ken(pCopula(u[, ind], cop)))
+  }
+  
+  if (ncol(kenVal) < copula at nestingCop@dimension) {
+    kenVal <- cbind(kenVal, u[, -sapply(copula at clusterCops, function(x) x[[2]])])
+  }
+  
+  lik <- cbind(lik, dCopula(kenVal, copula at nestingCop, log=log))
+  
+  if (log)
+    return(apply(lik, 1, sum))
+  
+  return(apply(lik, 1, prod))
+}
+
+setMethod("dCopula", signature("matrix", "hkCopula"), dHkCop)
+setMethod("dCopula", signature("numeric", "hkCopula"), 
+          function(u, copula, log, ...) dHkCop(matrix(u, ncol = copula at dimension), copula, log, ...))
+
+rHkCop <- function(n, copula, ...) {
+  smpl <- matrix(NA, n, copula at dimension)
+
+  nestSmpl <- rCopula(n, copula at nestingCop)
+  
+  for (i in 1:length(copula at clusterCops)) {
+    cop <- copula at clusterCops[[i]][[1]]
+    ind <- copula at clusterCops[[i]][[2]]
+    ken <- copula at kenFuns[[i]]
+    invKen <- genInvKenFun(ken)
+    
+    smpl[,ind] <- rCopula_y(invKen(nestSmpl[,i]), cop)
+  }
+  
+  if (ncol(nestSmpl) > length(copula at clusterCops)) {
+    smpl[,-sapply(copula at clusterCops, function(x) x[[2]])] <- nestSmpl[, -c(1:length(copula at clusterCops))]
+  }
+  
+  return(smpl)
+}
+
+setMethod(rCopula, signature = c("numeric","hkCopula"), rHkCop)
+
+setMethod(pCopula, signature = c("numeric","hkCopula"), 
+          function(u, copula, ...) stop("Please use an empirical representation (i.e. \"genEmpCop\" applied to a sample of this copula)."))
+
+setMethod(pCopula, signature = c("matrix","hkCopula"), 
+          function(u, copula, ...) stop("Please use an empirical representation (i.e. \"genEmpCop\" applied to a sample of this copula)."))
+
+
+rCop_y <- function(y, copula, n=1, n.disc = 1e2) {
+  stopifnot(copula at dimension == 2)
+  n.y <- length(y)
+  stopifnot(n.y == 1 | n == 1)
+
+  smpl <- matrix(NA, n.y*n, 2)
+  
+  for (i in 1:n.y) { # i <- 1
+    condVals <- seq(y[i], 1-(1-y[i])/n.disc, length.out = n.disc)
+    uv <- qCopula_v(copula, rep(y[i], n.disc), condVals)
+    uv <- rbind(uv, qCopula_u(copula, rep(y[i], n.disc), condVals))
+    
+    uv <- uv[order(uv[,1]),]
+    
+    dSeq <- cumsum(c(0, apply((uv[-nrow(uv),]-uv[-1,])^2, 1, function (x) sqrt(sum(x)))))
+    probs <- dCopula(uv, copula)
+    
+    apFun <- approxfun(dSeq, probs, rule = 2)
+    probCor <- integrate(apFun, 0, max(dSeq))$value
+    
+    rContour <- runif(n, 0, probCor)
+    
+    funAppConPoint <- function(rCont) {
+      invCDFContour <- function(x) {
+        abs(integrate(apFun, 0, x)$value - rCont)
+      } 
+      
+      lContour <- optimise(invCDFContour, c(0, max(dSeq)))$minimum
+      
+      dSeqInt <- findInterval(lContour, dSeq)
+      
+      lSeq <- sqrt(sum((uv[dSeqInt,]-uv[dSeqInt+1,])^2))
+      
+      uv[dSeqInt,] + (lContour - dSeq[dSeqInt])/lSeq * (uv[dSeqInt+1,]-uv[dSeqInt,])
+    }
+    
+    if (n == 1) {
+      appConPoint <- funAppConPoint(rContour)
+      
+      if (appConPoint[1] > appConPoint[2]) {
+        smpl[i,] <- qCopula_u(copula, y[i], appConPoint[1])
+      } else {
+        smpl[i,] <- qCopula_v(copula, y[i], appConPoint[2])
+      }
+    } else {
+      appConPoint <- t(sapply(rContour, funAppConPoint))
+      
+      boolLower <- appConPoint[,1] > appConPoint[,2]
+      smpl[boolLower,] <- qCopula_u(copula, rep(y, sum(boolLower)), appConPoint[boolLower, 1])
+      smpl[!boolLower,] <- qCopula_v(copula, rep(y, sum(!boolLower)), appConPoint[!boolLower, 2])
+    }
+    
+    # plot(uv, type="l", xlim=c(uv[dSeqInt+c(0,1)]+c(-1,1)/1000), asp=1)
+    # points(uv[dSeqInt+c(0,1),], col=c("red", "purple"))
+    # points(matrix(appConPoint, nrow = 1), col="green")
+    # points(matrix(smpl, nrow = 1), col="green", pch=2)
+  }
+
+  return(smpl)  
+}
+
+setGeneric("rCopula_y", function(y, copula, n=1, n.disc=1e2) NULL)
+setMethod("rCopula_y", signature("numeric", "copula"), rCop_y)
+
+# ## attic
+# hkCop <- new("hkCopula", nestingCop=normalCopula(0.6),
+#              clusterCops=list(list(cop=frankCopula(3), ind=c(1,2))),
+#              kenFuns = list(getKendallDistr(frankCopula(3))),
+#              dimension = 3L,
+#              parameters = NA_real_,
+#              param.names ="",
+#              param.lowbnd = NA_real_,
+#              param.upbnd = NA_real_,
+#              fullname = "Hierarchical Kendall Copula")
+# 
+# hkCop4D <- new("hkCopula", nestingCop=normalCopula(0.6),
+#                clusterCops=list(list(cop=frankCopula(3), ind=c(1,2)),
+#                                 list(cop=gumbelCopula(5), ind=c(3,4))),
+#                kenFuns = list(getKendallDistr(frankCopula(3)), 
+#                               getKendallDistr(gumbelCopula(5))),
+#                dimension = 4L,
+#                parameters = c(0),
+#                param.names ="",
+#                param.lowbnd = c(0),
+#                param.upbnd = 0, 
+#                fullname = "Hierarchical Kendall Copula")
+# 
+# rHkCop(10, hkCop)
+# 
+# smplRHkCop3D <- rCopula(100, hkCop)
+# 
+# library(rgl)
+# plot3d(smplRHkCop3D)
+# kenHKcop <- genEmpKenFun(hkCop, sample = smplRHkCop3D)
+# 
+# curve(kenHKcop)
+# 
+# showMethods("pCopula")
+# 
+# plot(rCopula_y(0.9, gumbelCopula(5), 100))
+# plot(rCopula_y(0.9, tawn3pCopula(c(0.75,.25,5)), 100))
+# plot(rCopula_y(0.9, normalCopula(0.4), 100), asp=1)
+# points(rCopula_y(0.9, normalCopula(0.8), 100), asp=1, pch=2)
+# points(rCopula_y(0.9, normalCopula(-0.8), 100), asp=1, pch=3)
+# 
+# points(rCopula_y(0.4, normalCopula(-0.3), 100), asp=1, pch=4)
+# abline(1.9,-1, col="red")
+# 
+# contour(normalCopula(-0.3), pCopula, asp=1)
+# 
+# sum(dHkCop(rHkCop(10, hkCop), hkCop))/10
+# sum(dHkCop(rHkCop(10, hkCop4D), hkCop4D))/10
+# 
+# sum(dHkCop(matrix(runif(3000), 1000), hkCop))/1000
+# sum(dHkCop(matrix(runif(4000), 1000), hkCop4D))/1000
+# 
+# par(mfrow=c(2,1))
+# hist(dHkCop(matrix(runif(4*1e5),ncol = 4), hkCop4D), n=4000, xlim=c(0,10))
+# hist(dHkCop(matrix(runif(3*1e5),ncol = 3), hkCop), n=400, xlim=c(0,10))
+# 
+# sum(dHkCop(matrix(runif(4*1e5),ncol = 4), hkCop4D))/1e5
+# sum(dHkCop(matrix(runif(3*1e5),ncol = 3), hkCop))/1e5
\ No newline at end of file

Modified: pkg/R/mixtureCopula.R
===================================================================
--- pkg/R/mixtureCopula.R	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/R/mixtureCopula.R	2017-02-07 11:35:01 UTC (rev 162)
@@ -30,7 +30,8 @@
       param.names = c(memberCops[[1]]@param.names, memberCops[[2]]@param.names, "mixLambda"),
       param.lowbnd = c(memberCops[[1]]@param.lowbnd, memberCops[[2]]@param.lowbnd, 0),
       param.upbnd = c(memberCops[[1]]@param.upbnd, memberCops[[2]]@param.upbnd, 1), 
-      fullname = paste("mixture of a", memberCops[[1]]@fullname, "and a", memberCops[[2]]@fullname))
+      fullname = paste("mixture of a", describeCop(memberCops[[1]], "very short"),
+                       "and a", describeCop(memberCops[[2]], "very short")))
 }
 
 ## density ##
@@ -71,14 +72,14 @@
 # invddu
 invdduMixCop <- function (u, copula, y) {
   stopifnot(length(u) == length(y)) 
- 
+  
   opti <- function(ind) {
     optFun <- function(v) {
       (dduCopula(cbind(u[ind], v), copula) - y[ind])^2
     }
     optimise(optFun, c(0,1))$minimum
   }
-
+  
   sapply(1:length(y), opti)  
 }
 
@@ -131,17 +132,16 @@
     lower <- copula at param.lowbnd
   if(is.null(upper))
     upper <- copula at param.upbnd
-    
+  
   optFun <- function(parSet) {
     cop <- mixtureCopula(parSet, copula at memberCops)
-    cat(cop at parameters, "\n")
     -sum(log(dCopula(data, cop)))
   }
   
   optOut <- optim(start, optFun, method = optim.method, 
                   lower = lower, upper = upper, 
                   control = optim.control)
-
+  
   new("fitCopula",
       copula = mixtureCopula(optOut$par, copula at memberCops),
       estimate = optOut$par,
@@ -165,6 +165,11 @@
             (1-mixLambda) * tau(copula at memberCops[[1]], ...) + mixLambda * tau(copula at memberCops[[2]], ...)
           })
 
+setMethod("rho", signature = c(copula = "mixtureCopula"),
+          function(copula, ...) {
+            mixLambda <- tail(copula at parameters, 1)
+            (1-mixLambda) * rho(copula at memberCops[[1]], ...) + mixLambda * rho(copula at memberCops[[2]], ...)
+          })
 
 setMethod("lambda", signature = c(copula = "mixtureCopula"),
           function(copula, ...) {

Modified: pkg/R/partialDerivatives.R
===================================================================
--- pkg/R/partialDerivatives.R	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/R/partialDerivatives.R	2017-02-07 11:35:01 UTC (rev 162)
@@ -266,7 +266,7 @@
 ##########################
 
 ## Wolfram alpha:
-# -e^α/(e^(α u) - e^α) - ((e^α - 1) e^(α + α u))/((e^(α u) - e^α) (e^α - e^(α + α u) + e^(α u + α v) - e^(α + α v)))
+# -e^a/(e^(a u) - e^a) - ((e^a - 1) e^(a + a u))/((e^(a u) - e^a) (e^a - e^(a + a u) + e^(a u + a v) - e^(a + a v)))
 
 dduFrank <- function(u, copula){
   rho <- copula at parameters

Modified: pkg/R/spCopula.R
===================================================================
--- pkg/R/spCopula.R	2017-02-01 15:31:26 UTC (rev 161)
+++ pkg/R/spCopula.R	2017-02-07 11:35:01 UTC (rev 162)
@@ -56,7 +56,7 @@
                                  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,
       fullname="Spatial Copula: distance dependent convex combination of bivariate copulas",
@@ -70,14 +70,11 @@
   cat("Copulas:\n")
   for (i in 1:length(object at components)) {
     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")
-#     }
+    cat("  ", describeCop(cmpCop, "very short"), "at", object at distances[i], 
+        paste("[",object at unit,"]",sep=""), "\n")
   }
-  if(!is.null(object at calibMoa(normalCopula(0),0))) cat("A spatial dependence function is used. \n")
+  if(!is.null(object at calibMoa(normalCopula(0),0))) 
+    cat("A spatial dependence function is used. \n")
 }
 
 setMethod("show", signature("spCopula"), showCopula)
@@ -555,7 +552,7 @@
   loglik <- NULL
   copulas <- list()
   for (cop in families) {
-    cat(cop at fullname,"\n")
+    cat(describeCop(cop, "very short"),"\n")
     tmploglik <- NULL
     tmpCop <- list()
     
@@ -609,7 +606,7 @@
   
   fits <-lapply(families, 
                 function(cop) {
-                  cat(cop at fullname,"\n")
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/spcopula -r 162


More information about the spcopula-commits mailing list