[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