## Copyright (C) 2010 Marius Hofert and Martin Maechler ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software ## Foundation; either version 3 of the License, or (at your option) any later ## version. ## ## This program is distributed in the hope that it will be useful, but WITHOUT ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS ## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more ## details. ## ## You should have received a copy of the GNU General Public License along with ## this program; if not, see . require(nacopula) set.seed(1) if(!dev.interactive()) pdf("copula-play.pdf") ##' ==== testing psi ==== myCop <- setTheta(copAMH, value = 0.5) # is maybe more natural setGeneric("psi", function(cop) standardGeneric("psi")) setMethod(psi, "acopula", function(cop) { function(t) cop@psi(t, theta = cop@theta) }) psi(myCop) # is a function psi(myCop)(0:4) curve(psi(myCop)(x), 0, 4) ##' but this can also be done directly [ => same curve "on top" :] curve(myCop@psi(x, theta = myCop@theta), 0, 4, col = 2, add = TRUE) ##' ==== testing Kendall's tau ==== p.Tau <- function(cop, n = 201, xlim = pmin(paraI, 50), ...) { stopifnot(is(cop, "acopula")) paraI <- cop@paraInterval theta <- seq(xlim[1], xlim[2], length.out = n) tit <- substitute(tau[NAME](theta), list(NAME = cop@name)) plot(theta, cop@tau(theta), type = "l", main = tit, ...) abline(h = c(0,1), lty = 3, col = "gray20") } p.Tau(copAMH) p.Tau(copClayton) p.Tau(copFrank, xlim = c(0, 80), ylim= 0:1) # fast via debye_1() p.Tau(copGumbel) p.Tau(copJoe, ylim = 0:1, yaxs="i") ##' ==== test function ==== ##' procedure to do several measurements ##' @param cop copula ##' @param theta1 parameter theta1 ##' @param thetavec vector of parameters ##' @param i10 values where psi is evaluated ##' @param nRnd number of generated V0's and V01's ##' @param t01 values where psiinv is evaluated ##' @param lambdaLvec vector of lower tail-dependence coefficients ##' @param lambdaUvec vector of upper tail-dependence coefficients ##' @return list of measurements ##' @author Marius Hofert, Martin Maechler tstCop <- function(cop, theta1 = cop@theta, thetavec = cop@theta, i10 = 1:10, nRnd = 50, t01 = (1:63)/64, # exact binary fractions lambdaLvec = NA_real_, lambdaUvec = NA_real_) { stopifnot(is(cop, "acopula")) cat0 <- function(...) cat(..., "\n", sep = "") theta0 <- cop@theta CT <- list() cat0(sprintf("(1) copula family: %10s, theta0 = %g", cop@name, theta0)) cat("\n(2) values of psi at i10:\n") CT <- c(CT, list(psi = system.time(p.i <- cop@psi(i10,theta = theta0)))) print(p.i) cat("check if psi(Inf)=0: ") stopifnot(cop@psi(Inf, theta = theta0)==0) cat0("TRUE") cat("check if psiInv(numeric(0)) is numeric(0): ") n0 <- numeric(0) stopifnot(identical(n0, cop@psiInv(n0, theta = theta0))) cat0("TRUE") cat("check if psiInv(0)=Inf: ") stopifnot(cop@psiInv(0, theta = theta0)==Inf) cat0("TRUE") cat0("values of psiInv at t01:\n") CT <- c(CT, list(psiI = system.time(pi.t <- cop@psiInv(t01, theta = theta0)))) print(pi.t) CT[["psiI"]] <- CT[["psiI"]] + system.time(pi.pi <- cop@psiInv(p.i,theta = theta0)) CT[["psi" ]] <- CT[["psi" ]] + system.time(p.pit <- cop@psi(pi.t, theta = theta0)) cat0("check if psiInv(psi(i10))==i10: ", all.equal(pi.pi, i10)) cat0("check if psi(psiInv(t01))==t01: ", all.equal(p.pit, t01)) cat("\n(3) parameter interval:\n") print(cop@paraInterval) cat0("theta1=",theta1) cat0("nesting condition for theta0 and theta1 fulfilled: ", cop@nestConstr(theta0,theta1)) CT <- c(CT, list(V0 = system.time(V0 <- cop@V0(nRnd,theta0)))) cat0("\n(4) ",nRnd," generated V0's:") print(summary(V0)) CT <- c(CT, list(V01 = system.time(V01 <- cop@V01(V0,theta0,theta1)))) cat0(nRnd," generated V01's:") print(summary(V01)) nt <- length(thetavec) cat("\n(5) tau at thetavec:\n") CT <- c(CT, list(tau = system.time(ta <- cop@tau(thetavec)))) print(ta) CT <- c(CT, list(tauI = system.time(ta.I <- cop@tauInv(ta)))) cat0("check if tauInv(tau(thetavec))==thetavec: ", all.equal(ta.I, thetavec)) lambdaLvec <- rep(as.double(lambdaLvec), length.out= nt) lambdaUvec <- rep(as.double(lambdaUvec), length.out= nt) cat("\n(6) lambdaL at thetavec:\n") CT <- c(CT, list(lambdaL = system.time(lT <- cop@lambdaL(thetavec)))) CT <- c(CT, list(lT.I = system.time(lT.I <- cop@lambdaLInv(lT)))) print(lT) cat0("check if lambdaLInv(lambdaL(thetavec))==lambdaLvec: ", all.equal(lT.I, lambdaLvec)) cat("\n(7) lambdaU at thetavec:\n") CT <- c(CT, list(lambdaU = system.time(uT <- cop@lambdaU(thetavec)))) CT <- c(CT, list(uT.I = system.time(uT.I <- cop@lambdaUInv(uT)))) print(uT) cat0("check if lambdaUInv(lambdaU(thetavec))==lambdaUvec: ", all.equal(uT.I, lambdaUvec)) class(CT) <- "proc_time_list" CT } ##' print() method for the tstCop() results print.proc_time_list <- function (x, ...) { stopifnot(is.list(x), !is.null(nx <- names(x))) for(nm in nx) if(!all(x[[nm]] == 0, na.rm=TRUE)) { cat(nm,":\n"); print(x[[nm]], ...) } invisible(x) } ##' ==== copAMH ==== myAMH <- setTheta(copAMH, 0.7135001) thetavec <- c(0.1,0.3,0.5,0.7,0.9) tstCop(myAMH, 0.9429679, thetavec = thetavec) ##' ==== copClayton ==== myClayton <- setTheta(copClayton, 0.5) thetavec <- c(0.5,1,2,5,10) tstCop(myClayton, 2, thetavec, lambdaL = thetavec, lambdaU = NA) ##' ==== copFrank === myFrank <- setTheta(copFrank, 1.860884) thetavec <- c(0.5,1,2,5,10) tstCop(myFrank, 5.736283, thetavec) ##' with a slightly more extensive test: tau.th <- c(0.055417, 0.11002, 0.21389, 0.4567, 0.66578) tau.F <- myFrank@tau(thetavec) stopifnot(all.equal(tau.th, tau.F, tol = 0.0001), all.equal(.9999, copFrank@tau(copFrank@tauInv(0.9999))), all.equal(myFrank@tauInv(tau.F, tol = 1e-14), thetavec, tol=1e-11)) ##' ==== copGumbel === myGumbel <- setTheta(copGumbel, 1.25) thetavec <- c(1,2,4,6,10) tstCop(myGumbel,2, thetavec, lambdaL = NA, lambdaU = thetavec) ##' ==== copJoe === myJoe <- setTheta(copJoe, 1.25) thetavec <- c(1.1,2,4,6,10) tstCop(myJoe, 2, thetavec, lambdaL = NA, lambdaU = thetavec)