[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