https://github.com/cran/nacopula
Raw File
Tip revision: e24ab1e1bade02a8136bd488815825e92fa6acd7 authored by Martin Maechler on 18 August 2010, 00:00:00 UTC
version 0.4-3
Tip revision: e24ab1e
AllClass.R
## 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 <http://www.gnu.org/licenses/>.

### TODO Should we try to be compatible to  'Copula' package ??

## This has advantage that arithmetic with scalars works "for free" already:
setClass("interval", contains =  "numeric", # of length 2
         representation(open  = "logical"),# of length 2
	 validity = function(object) {
	     if(length(rng <- object@.Data) != 2) "interval must be of length 2"
	     else if(length(object@open) != 2) "'open' must be of length 2"
	     else if(rng[2] < rng[1]) "'range[2]' must not be smaller than range[1]"
	     else TRUE
	 })

setClassUnion("maybeInterval", c("interval", "NULL"))

### Mother class of all (simple, *NON*-nested) Archimedean Copula Types
### for *any* dimension d
setClass("acopula",
	 representation(name = "character",
                        psi = "function",    # of (t, theta) -- the generator
                        psiInv = "function", # of (p, theta) -- psi_inverse: \psi^{-1}(p) = t
                        theta = "numeric", # value of theta or  'NA'  (for unspecified)
                        paraConstr = "function", # of (theta) ; constr(theta) |--> TRUE: "fulfilled"
                        ## when theta is one-dimensional, specifying the interval is more convenient:
                        paraInterval = "maybeInterval", # [.,.]  (.,.], etc ..
                        V0 = "function",	# of (n,theta) -- RNGenerator
                        tau= "function",	# of (theta)
                        tauInv = "function",    # of (tau)
                        lambdaL = "function",    # of (theta) lower bound  \lambda_l
                        lambdaLInv = "function", # of (lambda) - Inverse of \lambda_l
                        lambdaU = "function",    # of (theta)  - upper bound  \lambda_u
                        lambdaUInv = "function", # of (lambda) - Inverse of \lambda_u

                        ## Nesting properties if the child copulas are of the same family :
                        nestConstr = "function", # of (th0, th1) ; TRUE <==> "fulfilled"
                        V01= "function"	# of (V0,theta0,theta1)
                        ),
         prototype = prototype(theta = NA_real_),
	 validity = function(object) {
	     if (length(nm <- object@name) != 1 || nchar(nm) < 1)
		 return("'name' must be a string (of at least 1 character)")
             th <- tryCatch(validTheta(object), error = function(e) e)
             if(is(th, "error")) return(th $ message)

	     checkFun <- function(sName, nArgs, chkVectorizing=TRUE) {
		 f <- slot(object, sName)
		 if (length(formals(f)) < nArgs)
		     paste("slot '",sName,"' must be a function of at least ",nArgs,
			   " arguments", sep="")
		 else if(chkVectorizing && nArgs <= 2) {
		     ## test that the function can be called with NULL
		     r0 <- tryCatch(if(nArgs == 2) f(NULL, theta = th) else f(NULL),
                                    error = function(e) e)
		     if(is(r0, "error"))
			 sprintf("calling %s(NULL%s fails: %s", sName,
				 if(nArgs == 2) sprintf(", %g)", th) else ")",
				 r0$message)
		     else if(!identical(r0,
                                        {n0 <- numeric(0)
                                         if(nArgs == 2) f(n0, theta = th) else f(n0) }))
			 sprintf("%s(NULL ..) is not the same as %s(num(0) ..)",
				 sName, sName)
		     else TRUE
		 } else TRUE
	     } ## {checkFun}

             if(!isTRUE(tt <- checkFun("psi", 2)))	return(tt)
             if(!isTRUE(tt <- checkFun("psiInv", 2)))	return(tt)
             if(!isTRUE(tt <- checkFun("tau", 1)))	return(tt)

             if(!isTRUE(tt <- checkFun("paraConstr", 1, chkVect=FALSE))) return(tt)
             if(!isTRUE(tt <- checkFun("nestConstr", 2, chkVect=FALSE))) return(tt)

             ## ...
             ## ... (TODO)

             ## Check more :
	     if (object@psi(0, theta= 1/2) != 1)
		 return("psi(0, theta=1/2) != 1 -- seems an invalid generator")

             ## ....
             ## ....
             ## ....
             ## ....

	     ## 'else'	ok :
	     TRUE
	 })

## Construct 'paraConstr' slot automatically from 'paraInterval' :
setMethod("initialize", signature(.Object = "acopula"),
	  function(.Object, paraInterval, ...) {
	      if(!missing(paraInterval) && is(paraInterval, "interval")) {
		  .Object@paraConstr <- mkParaConstr(paraInterval)
		  .Object@paraInterval <- paraInterval
              }
              callNextMethod()
	  })


## A utility, not yet exported _ TODO ? _

setGeneric("validTheta", function(x, u) standardGeneric("validTheta"))
setMethod("validTheta", signature(x = "acopula"),
	  function(x) {
	      if(is((int <- x@paraInterval), "interval")) {
		  if(is.finite(d <- diff(as.numeric(int)))) # take mid-value in interval
		      int[1] + d/2
		  else if (all(iinf <- is.infinite(is.finite(int)))) ##	 (-Inf, Inf)
		      1/2
		  else { ## at least one end is finite
		      if(int[2] == Inf) int[1] + 1
		      else if(int[1] == -Inf) int[2] - 1
		      else stop("invalid paraInterval slot??")
		  }
	      } else  { ## need to use parameter-contraint function and "random search":
		  f.c <- x@paraConstr
		  for (th in c(1/2, (0:16)/16, round(exp(rnorm(20)),2)))
		      if(f.th <- f.c(th)) break
		  if(f.th) th else stop("could not find validTheta(<c>) for this copula")
	      }
	  })


### Nested Archimedean Copulas with *specified* dimension(s)
setClass("nacopula",
	 representation(copula = "acopula",
                        comp = "integer", # from 1:d -- of length in [0,d]
                        childCops = "list" #of nacopulas, possibly empty
                        ## TODO? nesting properties (V01,..) for specific mother-child relations
                        ),
         validity = function(object) {
             if(length(d <- dim(object)) != 1 || !is.numeric(d) || d <= 0)
                 return("invalid dim(.)")
             ic <- object@comp
             if((lc <- length(ic)) > d)
                 return("more components than the dimension")
	     if(!all("nacopula" == sapply(object@childCops, class)))
                 return("All 'childCops' elements must be 'nacopula' objects")
## FIXME: we need to recursively apply a "comp" function

             allC <- allComp(object)
             if(length(allC) != d)
                 return("must have d coordinates (from 'comps' and child copulas)")
             ##
             TRUE
         })


## FIXME: Maybe define a "strict_dim_nac" class which has the extra checks
## -----  for the "nacopula"s that appear as *SUB*-copulas, we do NOT want to
## require that their { components } are == { 1:d }
### Nested Archimedean Copulas with *specified* dimension(s)
setClass("outer_nacopula", contains = "nacopula",
         validity = function(object) {
             ## *Extra* checks in addition to those of "nacopula" :
             d <- dim(object)
             ic <- object@comp
             allC <- allComp(object)
             if(length(allC) != d)
                 return("must have d coordinates (from 'comps' and child copulas)")
             if(!all(sort(allC) == 1:d))
                 return(paste("The implicit coordinates are not identical to 1:d; instead\n  ",
                              paste(allC, collapse=", ")))
             ##
             TRUE
         })


## The dim() method is nicely defined  *recursive*ly :
setMethod("dim", signature(x = "nacopula"),
	  function(x) length(x@comp) + sum(unlist(lapply(x@childCops, dim))))

## also needed in validity above -- defined recursively, too :
allComp <- function(x) {
    stopifnot(is(x, "nacopula"))
    c(x@comp, unlist(lapply(x@childCops, allComp)))
}
back to top