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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Feb 9 17:07:24 CET 2015


Author: ben_graeler
Date: 2015-02-09 17:07:24 +0100 (Mon, 09 Feb 2015)
New Revision: 139

Added:
   pkg/R/KendallDistribution.R
   pkg/man/kendall-methods.Rd
Modified:
   /
   pkg/DESCRIPTION
   pkg/NAMESPACE
   pkg/R/returnPeriods.R
Log:
- added a new kendall function up to dimension 4 for Archimedean 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

   + .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


Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2014-10-01 09:31:41 UTC (rev 138)
+++ pkg/DESCRIPTION	2015-02-09 16:07:24 UTC (rev 139)
@@ -2,7 +2,7 @@
 Type: Package
 Title: copula driven spatial analysis
 Version: 0.2-1
-Date: 2014-10-01
+Date: 2015-02-06
 Authors at R: c(person("Benedikt", "Graeler", role = c("aut", "cre"),
                     email = "ben.graeler at uni-muenster.de"),
              person("Marius", "Appel",role = "ctb"))
@@ -10,7 +10,7 @@
 Description: This package provides a framework to analyse via copulas spatial and spatio-temporal data provided in the format of sp and spacetime package respectively. Additionally, support for calculating different multivariate return periods is implemented.
 License: GPL-3
 LazyLoad: yes
-Depends: copula (>= 0.999-7), R (>= 2.15.0), VineCopula (>= 1.2-1)
+Depends: copula (>= 0.999-12), VineCopula (>= 1.4), R (>= 3.1.0)
 Imports: methods, sp, spacetime (>= 1.0-9)
 Suggests: evd
 URL: http://r-forge.r-project.org/projects/spcopula/
@@ -33,4 +33,5 @@
   spatialGaussianCopula.R
   tawn3pCopula.R
   tailDependenceFunctions.R
+  KendallDistribution.R
   zzz.R

Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE	2014-10-01 09:31:41 UTC (rev 138)
+++ pkg/NAMESPACE	2015-02-09 16:07:24 UTC (rev 139)
@@ -40,6 +40,7 @@
 export(genEmpKenFun, genInvKenFun)
 export(kendallRP, criticalLevel, kendallDistribution, getKendallDistr)
 export(criticalPair, criticalTriple)
+export(kendall)
 
 ## classes
 exportClasses(asCopula, cqsCopula, tawn3pCopula, neighbourhood, stNeighbourhood, empiricalCopula)

Added: pkg/R/KendallDistribution.R
===================================================================
--- pkg/R/KendallDistribution.R	                        (rev 0)
+++ pkg/R/KendallDistribution.R	2015-02-09 16:07:24 UTC (rev 139)
@@ -0,0 +1,238 @@
+# derivation of kendall distributions in higher dimensions for Archimedean copulas
+
+###########
+## Frank ##
+###########
+
+# generator
+genFrank <- function(t, theta) copFrank at iPsi(t, theta)
+# -log( (exp(-theta*u)-1) / (exp(-theta)-1) )
+
+# use series expansion for small u?
+d1genFrank <- function(t, theta) {
+  theta / (1 - exp(theta * t))
+}
+
+d2genFrank <- function(t, theta) {
+  (theta^2 * exp(theta * t)) / (1 - exp(theta * t))^2
+}
+
+d3genFrank <- function(t, theta) {
+  -(theta^3 * exp(theta*t) * (exp(theta * t) + 1))/(exp(theta * t)-1)^3
+}
+
+## inverse generator
+igenFrank <- function(s, theta) copFrank at psi(s, theta)
+# -log(1-(1-exp(-theta))*exp(-t))/theta
+
+d1igenFrank <- function(s, theta) {
+  eth <- exp(theta)
+  (1 - eth) / (theta * (-eth + 1 + exp(theta + s)))
+}
+
+d2igenFrank <- function(s, theta) {
+  eth <- exp(theta)
+  eths <- exp(theta + s)
+  ((eth-1) * eths) / (theta * (-eth + eths + 1)^2)
+}
+
+d3igenFrank <- function(s, theta) {
+  eth <- exp(theta)
+  eths <- exp(theta + s)
+  ((eth - 1) * eths) / (theta * (-eth + eths + 1)^2)-(2*(eth-1) * exp(2 * theta + 2 * s))/(theta * (-eth + eths + 1)^3)  
+}
+
+kdFrank <- function(t, cop) {
+  # K^d (t) = t + sum_{i=1}^{d-1} (−1)^i/i! \varphi(t)^i (\varphi^{-1})^(i)(\varphi(t))
+  stopifnot(cop at dimension <=4)
+  .theta <- cop at parameters
+  .val <- 0 < t & t < 1
+  
+  gt <- genFrank(t[.val], .theta)
+  sum1 <- 0
+  for (i in 1:(cop at dimension-1)) {
+    digen <- switch(i,
+                    d1igenFrank(gt, .theta),
+                    d2igenFrank(gt, .theta),
+                    d3igenFrank(gt, .theta))
+    sum1 <- sum1 + (-1)^i/factorial(i) * gt^i * digen
+  }
+  
+  res <- t
+  res[.val] <- res[.val] + sum1
+  res
+}
+
+setMethod("kendall", signature("numeric", "frankCopula"), function(t, copula) kdFrank(t, copula))
+
+###################
+## Gumbel Copula ##
+###################
+# generator
+genGumbel <- function(t, theta) copGumbel at iPsi(t, theta)
+# -log(t)^theta
+
+# d1genGumbel <- function(t, theta) {}
+# d2genGumbel <- function(t, theta) {}
+# d3genGumbel <- function(t, theta) {}
+
+## inverse generator
+igenGumbel <- function(s, theta) copGumbel at psi(s, theta)
+# exp(-s^(1/theta))
+
+d1igenGumbel <- function(s, theta) {
+  -(exp(-s^(1/theta)) * s^(1/theta-1))/theta
+}
+
+d2igenGumbel <- function(s, theta) {
+  s1th <- s^(1/theta)
+  (exp(-s1th) * s^(1/theta-2) * (theta+s1th-1))/theta^2
+}
+
+d3igenGumbel <- function(s, theta) {
+  s1th <- s^(1/theta)
+  ems1th <- exp(-s1th)
+  s2th3 <- s^(2 / theta - 3)
+  -(ems1th * (theta + s1th - 1) * s2th3) / (theta^3 + ems1th * s2th3) / theta^3 + ( (1 / theta - 2) * ems1th * (theta + s1th-1) * s^(1 / theta - 3)) / theta^2
+}
+
+kdGumbel <- function(t, cop) {
+  # K^d (t) = t + sum_{i=1}^{d-1} (−1)^i/i! \varphi(t)^i (\varphi^{-1})^(i)(\varphi(t))
+  stopifnot(cop at dimension <=4)
+  .theta <- cop at parameters
+  .val <- 0 < t & t < 1
+  
+  gt <- genGumbel(t[.val], .theta)
+  sum1 <- 0
+  for (i in 1:(cop at dimension-1)) {
+    digen <- switch(i,
+                    d1igenGumbel(gt, .theta),
+                    d2igenGumbel(gt, .theta),
+                    d3igenGumbel(gt, .theta))
+    sum1 <- sum1 + (-1)^i/factorial(i) * gt^i * digen
+  }
+  
+  res <- t
+  res[.val] <- res[.val] + sum1
+  res
+}
+
+setMethod("kendall", signature("numeric", "gumbelCopula"), function(t, copula) kdGumbel(t, copula))
+
+#############
+## Clayton ##
+#############
+
+# generator
+genClayton <- function(t, theta) copClayton at iPsi(t, theta)
+# u^(-theta) - 1
+
+# d1genClayton <- function(t, theta) {}
+# d2genClayton <- function(t, theta) {}
+# d3genClayton <- function(t, theta) {}
+
+## inverse generator
+igenClayton <- function(s, theta) copClayton at psi(s, theta)
+# (1 + t)^(-1/theta)
+
+d1igenClayton <- function(s, theta) {
+  -(s+1)^(-(theta+1)/theta)/theta
+}
+
+d2igenClayton <- function(s, theta) {
+  ((theta + 1) * (s + 1)^(-1 / theta - 2)) / theta^2
+}
+
+d3igenClayton <- function(s, theta) {
+  -((theta + 1) * (2 * theta + 1) * (s + 1)^(-1 / theta - 3)) / theta^3
+}
+
+kdClayton <- function(t, cop) {
+  # K^d (t) = t + sum_{i=1}^{d-1} (−1)^i/i! \varphi(t)^i (\varphi^{-1})^(i)(\varphi(t))
+  stopifnot(cop at dimension <=4)
+  .theta <- cop at parameters
+  .val <- 0 < t & t < 1
+  
+  gt <- genClayton(t[.val], .theta)
+  sum1 <- 0
+  for (i in 1:(cop at dimension-1)) {
+    digen <- switch(i,
+                    d1igenClayton(gt, .theta),
+                    d2igenClayton(gt, .theta),
+                    d3igenClayton(gt, .theta))
+    sum1 <- sum1 + (-1)^i/factorial(i) * gt^i * digen
+  }
+  
+  res <- t
+  res[.val] <- res[.val] + sum1
+  res
+}
+
+setMethod("kendall", signature("numeric", "claytonCopula"), function(t, copula) kdClayton(t, copula))
+
+#########
+## Joe ##
+#########
+
+
+# generator
+genJoe <- function(t, theta) copJoe at iPsi(t, theta)
+# -log(1-(1 - t)^theta)
+
+# d1genClayton <- function(t, theta) {}
+# d2genClayton <- function(t, theta) {}
+# d3genClayton <- function(t, theta) {}
+
+## inverse generator
+igenJoe <- function(s, theta) copJoe at psi(s, theta)
+# 1 - (1 - exp(-s))^(1 / theta)
+
+d1igenJoe <- function(s, theta) {
+  ( -expm1(-s))^(1/theta)/(theta-theta * exp(s))
+}
+
+d2igenJoe <- function(s, theta) {
+  exps <- exp(s)
+  expm1ms <- expm1(-s)
+  expm1s <- expm1(s)
+  
+  ((-expm1ms)^(1/theta) * (theta * exps - 1))/(theta * expm1s)^2
+}
+
+# (e^(-s) (1-e^(-s))^(1/theta-1))/(theta (theta-theta e^s))
+# (theta e^s (1-e^(-s))^(1/theta))/(theta-theta e^s)^2
+
+d3igenJoe <- function(s, theta) {
+  exps <- exp(s)
+  expm1ms <- expm1(-s)
+  expm1s <- expm1(s)
+  
+  ds1 <- (-expm1ms)^(1 / theta) * (2 * theta * exps - 1)
+  ds2 <- -theta * (exps * (-expm1ms)^(1/theta) * (theta + theta * exps - 1))
+  
+  (ds1 + ds2) / (theta * expm1s)^3 
+}
+
+kdJoe <- function(t, cop) {
+  # K^d (t) = t + sum_{i=1}^{d-1} (−1)^i/i! \varphi(t)^i (\varphi^{-1})^(i)(\varphi(t))
+  stopifnot(cop at dimension <= 4)
+  .theta <- cop at parameters
+  .val <- 0 < t & t < 1
+  
+  gt <- genJoe(t[.val], .theta)
+  sum1 <- 0
+  for (i in 1:(cop at dimension-1)) {
+    digen <- switch(i,
+                    d1igenJoe(gt, .theta),
+                    d2igenJoe(gt, .theta),
+                    d3igenJoe(gt, .theta))
+    sum1 <- sum1 + (-1)^i/factorial(i) * gt^i * digen
+  }
+  
+  res <- t
+  res[.val] <- res[.val] + sum1
+  res
+}
+
+setMethod("kendall", signature("numeric", "joeCopula"), function(t, copula) kdJoe(t, copula))
+setMethod("kendall", signature("numeric", "joeBiCopula"), function(t, copula) kdJoe(t, copula))
\ No newline at end of file

Modified: pkg/R/returnPeriods.R
===================================================================
--- pkg/R/returnPeriods.R	2014-10-01 09:31:41 UTC (rev 138)
+++ pkg/R/returnPeriods.R	2015-02-09 16:07:24 UTC (rev 139)
@@ -2,7 +2,11 @@
   if(is.null(sample)) sample <- rCopula(1e6,copula)
   # as empirical copula:
   # copula <- genEmpCop(copula, sample)
-  ken <- pCopula(sample, copula)
+  if(missing(copula)) {
+    ken <- mapply(function(x,y) sum(x > sample[,1] & y > sample[,2])/length(x), sample[,1], sample[,2])
+  } else {
+    ken <- pCopula(sample, copula)
+  }
   
   empKenFun <- function(tlevel) {
     sapply(tlevel,function(t) sum(ken<=t))/nrow(sample)
@@ -165,6 +169,11 @@
 
 ## kendall distribution
 
+# generic kendall function
+kendall <- function(t, copula) {
+  standardGeneric(t, copula)
+}
+
 # empirical default
 getKendallDistr <- function(copula, sample=NULL) {
   standardGeneric("getKendallDistr")
@@ -187,7 +196,7 @@
 ## 
 
 kendallDistribution <- function(copula, t) {
-  stop("There is no analytical expression implemented for this copula family. See 'getKendallDstr' for a numerical solution instead.")
+  stop("There is no analytical expression implemented for this copula family. See 'getKendallDistr' for a numerical solution instead.")
 }
 
 setGeneric("kendallDistribution")
@@ -204,7 +213,8 @@
   return(kt)  
 }
 
-setMethod("kendallDistribution", signature("claytonCopula"), kendall.Clayton)
+setMethod("kendallDistribution", signature("claytonCopula"), kendall.Clayton) # for easy backwards compatibility
+setMethod("kendall", signature("numeric", "claytonCopula"), function(t, copula) kendall.Clayton(copula, t))
 
 setMethod("getKendallDistr", signature("claytonCopula"), 
           function(copula) return(function(t) kendall.Clayton(copula, t)))
@@ -221,7 +231,8 @@
   return(kt)  
 }
 
-setMethod("kendallDistribution", signature("gumbelCopula"), kendall.Gumbel)
+setMethod("kendallDistribution", signature("gumbelCopula"), kendall.Gumbel) # for easy backwards compatibility
+setMethod("kendall", signature("numeric","gumbelCopula"), function(t, copula) kendall.Gumbel(copula, t)) 
 
 setMethod("getKendallDistr", signature("gumbelCopula"), 
           function(copula) return(function(t) kendall.Gumbel(copula, t)))
@@ -238,7 +249,8 @@
   return(kt)  
 }
 
-setMethod("kendallDistribution", signature("frankCopula"), kendall.Frank)
+setMethod("kendallDistribution", signature("frankCopula"), kendall.Frank) # for easy backwards compatibility
+setMethod("kendall", signature("numeric", "frankCopula"),  function(t, copula) kendall.Frank)
 
 setMethod("getKendallDistr", signature("frankCopula"), 
           function(copula) return(function(t) kendall.Frank(copula, t)))
@@ -259,6 +271,11 @@
 }
 
 setMethod("kendallDistribution", signature("BB1Copula"), kendall.BB1)
+setMethod("kendall", signature("numeric", "BB1Copula"), 
+          function(t, copula) {
+            stopifnot(copula at dimension <= 2) 
+            kendall.BB1(copula, t)
+          })
 
 setMethod("getKendallDistr", signature("BB1Copula"), function(copula) return(function(t) kendall.BB1(copula, t)) )
 
@@ -276,7 +293,12 @@
   return(kt)  
 }
 
-setMethod("kendallDistribution", signature("BB6Copula"), kendall.BB6)
+setMethod("kendallDistribution", signature("BB6Copula"), kendall.BB6) # for easy backwards compatibility
+setMethod("kendall", signature("numeric", "BB6Copula"), 
+          function(t, copula) {
+            stopifnot(copula at dimension <= 2) 
+            kendall.BB6(copula, t)
+          })
 
 setMethod("getKendallDistr", signature("BB6Copula"), 
           function(copula) return(function(t) kendall.BB6(copula, t)))
@@ -296,6 +318,11 @@
 }
 
 setMethod("kendallDistribution", signature("BB7Copula"), kendall.BB7)
+setMethod("kendall", signature("numeric", "BB7Copula"), 
+          function(t, copula) {
+            stopifnot(copula at dimension <= 2) 
+            kendall.BB7(copula, t)
+          })
 
 setMethod("getKendallDistr", signature("BB7Copula"), 
           function(copula) return(function(t) kendall.BB7(copula, t)))
@@ -313,23 +340,32 @@
   return(kt)  
 }
 
-setMethod("kendallDistribution", signature("BB8Copula"), kendall.BB8)
+setMethod("kendallDistribution", signature("BB8Copula"), kendall.BB8) # for easy backwards compatibility
+setMethod("kendall", signature("numeric", "BB8Copula"), 
+          function(t, copula) {
+            stopifnot(copula at dimension <= 2) 
+            kendall.BB8(copula, t)
+          })
+
 setMethod("getKendallDistr", signature("BB8Copula"), 
           function(copula) return(function(t) kendall.BB8(copula, t)))
 
 # BiJoe
-## kendall distribution/measure, taken from VineCopula:::obs.stat
-kendall.Joe <- function(copula, t){
-  par = copula at parameters[1]
+kendall.Joe <- function(copula, t) kdJoe(t, copula)
   
-  kt <- rep(NA,length(t))
-  kt <- t - (log(1 - (1 - t)^par) * (1 - (1 - t))^par)/(par * (1 - t)^(par - 1))
-  kt[t==1] <- 1
-  kt[t==0] <- 0
-  return(kt)  
-}
+  ## kendall distribution/measure, taken from VineCopula:::obs.stat
+#   {
+#   par = copula at parameters[1]
+#   
+#   kt <- rep(NA,length(t))
+#   kt <- t - (log(1 - (1 - t)^par) * (1 - (1 - t))^par)/(par * (1 - t)^(par - 1))
+#   kt[t==1] <- 1
+#   kt[t==0] <- 0
+#   return(kt)  
+# }
 
 setMethod("kendallDistribution", signature("joeBiCopula"), kendall.Joe)
+
 setMethod("getKendallDistr", signature("joeBiCopula"), 
           function(copula) return(function(t) kendall.Joe(copula, t)))
 

Added: pkg/man/kendall-methods.Rd
===================================================================
--- pkg/man/kendall-methods.Rd	                        (rev 0)
+++ pkg/man/kendall-methods.Rd	2015-02-09 16:07:24 UTC (rev 139)
@@ -0,0 +1,31 @@
+\name{kendall-methods}
+\docType{methods}
+\alias{kendall}
+\alias{kendall-methods}
+\alias{kendall,ANY,ANY-method}
+\alias{kendall,numeric,claytonCopula-method}
+\alias{kendall,numeric,frankCopula-method}
+\alias{kendall,numeric,gumbelCopula-method}
+\alias{kendall,numeric,BB1Copula-method}
+\alias{kendall,numeric,BB6Copula-method}
+\alias{kendall,numeric,BB7Copula-method}
+\alias{kendall,numeric,BB8Copula-method}
+\alias{kendall,numeric,joeBiCopula-method}
+\alias{kendall,numeric,joeCopula-method}
+
+\title{ Methods for Function \code{kendall} }
+\description{
+  Methods for function \code{kendall} in case of Archimedian copulas up to 4 dimensions
+}
+\section{Methods}{
+\describe{
+
+\item{\code{signature(t = "ANY", copula = "ANY")}}{nothing}
+
+\item{\code{signature(t = "numeric", copula = "claytonCopula")}}{Kendall distribtuion for Clayton family.}
+
+\item{\code{signature(t = "numeric", copula = "frankCopula")}}{Kendall distribtuion for Frank family.}
+
+\item{\code{signature(t = "numeric", copula = "gumbelCopula")}}{Kendall distribtuion for Gumbel family.}
+}}
+\keyword{methods}



More information about the spcopula-commits mailing list