## Copyright (C) 2010--2011 Marius Hofert and Martin Maechler
##
## This program is free software; you can redistribute it and/or modify it under
## 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/>.

paste0 <- function(...) paste(..., sep="")

#### Functions and Methods for "acopula" objects
#### class definition in ./AllClass.R

### Ali-Mikhail-Haq == "AMH" ###################################################

##' @title Ali-Mikhail-Haq ("AMH")'s  tau(theta)
##' @param th
##' @return 1 - 2*((1-th)*(1-th)*log(1-th)+th)/(3*th*th)
##' numerically accurately, for both limits  th -> 0  &  th -> 1
##' @author Martin Maechler
tauAMH <- function(th) {
if(length(th) == 0) return(numeric(0)) # to work with NULL
r <- th
na <- is.na(th)
lrg <- (th > 0.01) & !na ## for "large" theta (<=> theta >= 0.01), use standard formula
f. <- function(t) {
## 1 - 2*((1-t)*(1-t)*log1p(-t) + t) / (3*t*t)
r <- t
r[i1 <- (1-t) == 0] <- 1/3
t <- t[s <- !i1]
r[s] <- 1 - 2*((1-t)^2 *log1p(-t) + t) / (3*t*t)
r
}
r[lrg] <- f.(th[lrg])
## Otherwise use "smart" Taylor polynomials: for given order, *drop* the last two..
## with the following coefficients:
##_ x <- polynomial()
##_ MASS::fractions(as.vector(mkAMHtau.0series(12,TRUE) * 9/(2*x)))
##  1  1/4  1/10  1/20  1/35  1/56  1/84  1/120  1/165  1/220

if(any(!lrg & !na)) {
l1 <- !lrg & !na & (ll1 <- th > 2e-4) ## --> k = 7
r[l1] <- (function(x) 2*x/9*(1+ x*(1/4 + x/10*(1 + x*(1/2 + x/3.5)))))(th[l1])
l2 <- !ll1 & !na & (ll2 <- th > 1e-5)	 ## --> k = 6
r[l2] <- (function(x) 2*x/9*(1+ x*(1/4 + x/10*(1 + x/2))))(th[l2])
l3 <- !ll2 & !na & (ll <- th > 2e-8)	## --> k = 5
r[l3] <- (function(x) 2*x/9*(1+ x*(1/4 + x/10)))(th[l3])
irem <- which(!ll & !na)## rest: th <= 2e-8 : k == 4
r[irem] <- (function(x) 2*x/9*(1+ x/4))(th[irem])
}
r
}

### Clayton ####################################################################

##' Note: this is mainly to show that this function can be very well
##' approximated much more simply by just using m <- round(V0).
##'
##' @title Optimal constant for fast rejection
##' @param V0 numeric vector >= 0
##' @return optimal constant m for the fast rejection algorithm
##' @author Martin Maechler (based on Marius Hofert's code)
m.opt.retst <- function(V0) {
n <- length(V0)
fV <- floor(V0)
cV <- ceiling(V0)
v1 <- fV*exp(V0/fV)
v2 <- cV*exp(V0/cV)

m <- integer(n)
l1 <- (V0 <= 1)
m[which(l1)] <- 1L

i2 <- which(!l1) ## those with V0 > 1
l3 <- (v1[i2] <= v2[i2])
i3 <- i2[l3]
m[i3] <- fV[i3]
i4 <- i2[!l3]
m[i4] <- cV[i4]
m
}

### fast rejection for fixed m, R version

##' Sample a random variate St ~ \tilde{S}(alpha, 1, (cos(alpha*pi/2)
##' *V_0)^{1/alpha}, V_0*I_{alpha = 1}, I_{alpha != 1}; 1) with
##' Laplace-Stieltjes transform exp(-V_0((1+t)^alpha-1)), see Nolan's book for
##' the parametrization, via an m-fold sum of random variates from
##' \tilde{S}(alpha, 1, (cos(alpha*pi/2)*V_0/m)^{1/alpha}, (V_0/m)
##' *I_{alpha = 1}, I_{alpha != 1}; 1) with Laplace-Stieltjes transform
##' exp(-(V_0/m)*((1+t)^alpha-1)). This is a building block for the fast rejection.
##'
##' @title Sample an exponentially tilted stable distribution as an m-fold sum
##' @param m number of summands, any positive integer
##' @param V0 random variate
##' @param alpha parameter in (0,1]
##' @return St
##' @author Marius Hofert, Martin Maechler
retstablerej <- function(m,V0,alpha) {
gamm. <- (cos(alpha*pi/2)*V0/m)^(1/alpha)
sum(unlist(lapply(integer(m),
function(.) {
## apply standard rejection for sampling
## \tilde{S}(alpha, 1, (cos(alpha*pi/2)
##	*V_0/m)^{1/alpha}, (V_0/m)*I_{alpha = 1},
## h*I_{alpha != 1}; 1) with Laplace-Stieltjes
## transform exp(-(V_0/m)*((h+t)^alpha-h^alpha))
repeat {
V__ <- rstable1(1, alpha, beta=1, gamma = gamm.)
if(runif(1) <= exp(-V__))
return(V__)
}})
## on acceptance, St_k ~ \tilde{S}(alpha, 1, (cos(alpha*pi/2)
## *V_0/m)^{1/alpha}, (V_0/m)*I_{alpha = 1}, h*I_{alpha != 1};
## 1) with Laplace-Stieltjes transform
## exp(-(V_0/m)*((h+t)^alpha-h^alpha))
))
}

### fast rejection, R version

##' Sample a vector of random variates St ~ \tilde{S}(alpha, 1, (cos(alpha*pi/2)
##' *V_0)^{1/alpha}, V_0*I_{alpha = 1}, h*I_{alpha != 1}; 1) with
##' Laplace-Stieltjes transform exp(-V_0((h+t)^alpha-h^alpha)), see Nolan's book for
##' the parametrization, with the fast rejection. This procedure calls retstablerej.
##'
##' @title Sampling an exponentially tilted stable distribution
##' @param alpha parameter in (0,1]
##' @param V0 vector of random variates
##' @param h non-negative real number
##' @return vector of variates St
##' @author Marius Hofert, Martin Maechler
retstableR <- function(alpha, V0, h=1) {
n <- length(V0)
stopifnot(is.numeric(alpha), length(alpha) == 1,
0 <= alpha, alpha <= 1) ## <- alpha > 1 ==> cos(pi/2 *alpha) < 0
## case alpha==1
if(alpha == 1 || n == 0) { # alpha == 1 => St corresponds to a point mass at V0 with
return(V0) # Laplace-Stieltjes transform exp(-V0*t)
}
## else alpha != 1 : call fast rejection algorithm with optimal m
m <- m.opt.retst(V0)
mapply(retstablerej, m=m, V0=V0, alpha=alpha)
}

### state-of-the-art: fast rejection + Luc's algorithm, C version

##' Sample random variates St ~ \tilde{S}(alpha, 1, (cos(alpha*pi/2)
##' *V_0)^{1/alpha}, V_0*I_{alpha = 1}, h*I_{alpha != 1}; 1) with
##' Laplace-Stieltjes transform exp(-V_0((h+t)^alpha-h^alpha)), see Nolan's book for
##' the parametrization, with the fast rejection.
##' This procedure is more efficient than retstableR since it calls the C
##' function retstable_c and uses both the fast rejection and Luc Devroye's algorithm.
##'
##' @title Efficiently sampling an exponentially tilted stable distribution
##' @param alpha parameter in (0,1]
##' @param V0 vector of random variates
##' @param h non-negative real number
##' @param method which method to call ("Marius Hofert", "Luc Devroye")
##' @return vector of variates St
##' @author Martin Maechler
retstableC <- function(alpha, V0, h = 1, method = NULL) {
n <- length(V0)
stopifnot(is.numeric(alpha), length(alpha) == 1,
0 < alpha, alpha <= 1,
is.numeric(h), length(h) == 1, h > 0)
if(alpha == 1 || n == 0) {
## alpha == 1 => St corresponds to a point mass at V0 with
V0           # Laplace-Stieltjes transform exp(-V0*t)
}
else {
if(is.null(method)) {
if(any(diff(V.is.sml <- V0 * h^alpha < 4))) { ## use *both* methods
r <- numeric(n)
r[ V.is.sml] <- .Call(retstable_c, V0[ V.is.sml], h = h, alpha, "MH")
r[!V.is.sml] <- .Call(retstable_c, V0[!V.is.sml], h = h, alpha, "LD")
return(r)
}
else
method <- if(V.is.sml[1]) "MH" else "LD"
}
else
method <- match.arg(method, c("MH","LD"))
.Call(retstable_c, V0, h = h, alpha, method)
}
}

## switch to make fast C code the default
if(FALSE)
retstable <- retstableR
retstable <- retstableC # retstable is by default retstableC

### Frank ######################################################################

### sampling a logarithmic distribution, R version

##' Random number generator for a Log(p) distribution with the algorithm "LK" of
##' Kemp (1981), R version.
##'
##' @title Sample a Log(p) distribution
##' @param n number of random variates to be generated
##' @param p parameter in (0,1)
##' @param Ip = 1 - p_ (possibly more accurate) -- use, if p ~= 1
##' @return vector of random variates from Log(p)
##' @author Marius Hofert, Martin Maechler
rlogR <- function(n, p, Ip = 1-p) {
if(missing(p)) p <- 1 - Ip
stopifnot((n <- as.integer(n)) >= 0,
0 < p, p <= 1, 0 < Ip)
vec <- numeric(n)
if(n >= 1) {
u <- runif(n)
l1 <- u > p
vec[l1] <- 1
i2 <- which( !l1 ) # of shorter length, say n2
q2 <- 1-(Iq2 <- Ip^runif(length(i2))) # length n2
l3 <- u[i2] < q2*q2
i3 <- i2[l3]
vec[i3] <- floor(1+abs(log(u[i3])/log1p(-Iq2[l3])))
l4 <- u[i2] > q2
vec[i2[l4]] <- 1
l5 <- ! (l3 | l4) # q2^2 <= u[i2] <= q2
vec[i2[l5]] <- 2
}
vec
}

### state-of-the art: sampling a logarithmic distribution, C version

##' Random number generator for a Log(p) distribution with the algorithm "LK" of
##' Kemp (1981), C version.
##'
##' @title Efficiently sampling a Log(p) distribution
##' @param n number of random variates to be generated
##' @param p parameter in (0,1)
##' @param Ip = 1 - p_ (possibly more accurate)
##' @return vector of random variates from Log(p)
##' @author Martin Maechler
rlog <- function(n, p, Ip = 1-p) {
if(missing(p)) p <- 1 - Ip
stopifnot(n >= 0, 0 < p, p <= 1, 0 < Ip)
.Call(rLog_vec_c, n, p, Ip)
}

### state-of-the art: sampling F01Frank, C version

##' Generate a vector of variates V01 ~ F01 with Laplace-Stieltjes transform
##' ((1-(1-exp(-t)*(1-e^(-theta1)))^alpha)/(1-e^(-theta0)))^V0.
##'
##' @title Efficiently sampling V01 for Frank
##' @param V0 vector of random variates from F0
##' @param theta0 parameter theta0 in (0,infinity)
##' @param theta1 parameter theta1 in [theta0, infinity)
##' @param rej method switch: if V0*theta_0*p0^(V0-1) > rej a rejection
##'        from F01 of Joe is applied (otherwise, the sum is
##'        sampled via a logarithmic envelope for the summands)
##' @param approx largest number of summands before asymptotics is used
##' @return vector of random variates V01
##' @author Marius Hofert
rF01Frank <- function(V0, theta0, theta1, rej, approx) {
.Call(rF01Frank_vec_c, V0, theta0, theta1, rej, approx)
}

### wrapper for inner distribution F for Frank

##' Generate a vector of variates V ~ F with Laplace-Stieltjes transform
##' (1-(1-exp(-t)*(1-e^(-theta1)))^alpha)/(1-e^(-theta0)).
##'
##' @title Sampling F for Frank
##' @param n number of variates from F
##' @param theta0 parameter theta0 in (0,infinity)
##' @param theta1 parameter theta1 in [theta0, infinity)
##' @param rej method switch: if theta_0 > rej a rejection
##'        from Joe (Sibuya distribution) is applied (otherwise, a logarithmic
##' 	   envelope is used)
##' @return vector of random variates V
##' @author Marius Hofert
rFFrank <- function(n, theta0, theta1, rej)
rF01Frank(rep(1,n),theta0,theta1,rej,1) # approx = 1 (dummy)

dDiagFrank <- function(u, theta, d, log=FALSE,
method = c("auto", paste("poly",1:4,sep=""),
"polyFull",  "m1", "MH", "MMH"))
{
stopifnot(length(d <- as.integer(d)) == 1, d >= 1)
r <- ut <- u*theta
Ie <- -expm1(-ut) # 1 - exp(-u * theta)
## L <- ut > log(2)*.Machine$double.digits # <==> exp(-ut) < Mach..eps method <- match.arg(method) if(method == "auto") method <- { if(d <= 3) "poly1" else # d in {1,2,3} if(d <= 6) paste("poly", d-2, sep="") # d in {4,5,6} else ## this is not really good: but (u*th) is *vector* "polyFull" } if(substr(method, 1,4) == "poly") { e.ut <- exp(-ut) ep <- (e.ut - exp(ut-theta))/Ie # "FIXME": maybe improve, testing u > 1/2 d1 <- d-1 D <- d + d1*ep # D = d + (d-1) * (exp(-theta*u)-exp(theta*u-theta)) if(d > 2) { # <==> d1 = d-1 >= 2 f <- 1 + ep # 1+exp(-theta*u)-exp(theta*u-theta) delt <- e.ut * f # exp(-theta*u) * (1+exp(-theta*u)-exp(theta*u-theta)) D <- D + f* delt * switch(method, "poly1" = d1*(d-2L)/2, "poly2" = d1*(d-2L)/2 *(1 + (d-3L)/3 * delt), "poly3" = d1*(d-2L)/2 *(1 + (d-3L)/3 * delt * (1 + (d-4L)/4 * delt)), "poly4" = d1*(d-2L)/2 *(1 + (d-3L)/3 * delt * (1 + (d-4L)/4 * delt*(1 + (d-5L)/5 * delt))), "polyFull" = { ## full polynomial formula if(is(ut, "mpfr")) sfsmisc::polyn.eval(chooseMpfr.all(d1)[-1], x = delt) else polynEval(choose (d1, 2:d1), x = delt) }, stop("invalid poly* method: ", method, ". Should never happen") ) } if(log) log(d) - log(D) else d / D } else switch(method, "m1" = { h <- -expm1(-theta) # = 1 - exp(-theta) D <- (h/Ie)^(d-1) - Ie # D := [d]enominator if(log) log(d) -ut -log(D) else d*exp(-ut) / D }, "MH" = { # Marius Hofert h <- -expm1(-theta) # = 1-exp(-theta) x <- Ie/h # = (1-exp(-theta*u)) / (1-exp(-theta)) if(log) log(d)+(d-1)*log(x)+log((1-h*x)/(1-h*x^d)) else d*x^(d-1)*(1-h*x)/(1-h*x^d) }, "MMH" = { # Martin's version of "MH" h <- -expm1(-theta) # = 1 - exp(-theta) x <- Ie/h #-> h*x = Ie r <- ## log( (1-h*x)/(1-h*x^d) ); 1-h*x = 1-Ie = exp(-u * theta) = exp(-ut) -ut - log1p(-h*x^d) if(log) log(d)+(d-1)*log(x) + r else d*x^(d-1)*exp(r) }, stop("impossible method: ", method,". Should never happen")) } ### Gumbel ##################################################################### ### compute the coefficients for polyG ##' Compute the coefficients a_{dk}(theta) involved in the generator derivatives ##' and the copula density of a Gumbel copula ##' ##' @title Coefficients of the polynomial involved in the generator derivatives ##' and density for Gumbel ##' @param d number of coefficients, d >= 1 ##' @param alpha parameter (1/theta) in (0,1]; ##' you maye use mpfr(alph, precBits = <n_prec>) for higher precision "Rmpfr*" methods ##' @param method a character string, one of ##' "sort": compute coefficients via exp(log()) pulling out the maximum, and sort ##' "horner": uses polynEval() ##' "direct": brute force approach ##' "dsSib.<FOO>": uses dsumSibuya(..., method= "<FOO>") ##' @param log boolean which determines if the logarithm is returned ##' @param verbose logical for method == sort ##' @return a_k(theta, d) = (-1)^{d-k}\sum_{j=k}^d alpha^j * s(d,j) * S(j,k), k in ##' {1,..,d} ##' @author Marius Hofert and Martin Maechler ##' Note: - This function is known to cause numerical problems, e.g., for d=100, ##' alpha=0.8 ##' - sign(s(n,k)) = (-1)^{n-k} coeffG <- function(d, alpha, method = c("sort", "horner", "direct", "dsumSibuya", paste("dsSib", eval(formals(dsumSibuya)$method), sep=".")),
log = FALSE, verbose = FALSE)
{
stopifnot(is.numeric(d), length(d) == 1, d >= 1, length(alpha) == 1,
0 < alpha, alpha <= 1)
a <- numeric(d) # for the a_k(theta, d)'s
method <- match.arg(method)
method <- "dsSib.log"
}
Meth <-
if(grepl("^dsSib", method)) {
meth.dsSib <- sub("^dsSib\\.(.*)", "\\1", method)
"dsSib"
} else method
switch(Meth,
"sort" = {
ls <- log(abs(Stirling1.all(d))) # log(|s(d, i)|), i=1:d
lS <- lapply(1:d, function(n) log(Stirling2.all(n)))
##-> lS[[n]][i] contains log(S(n,i)), i = 1,..,n
wrong.sign <- integer()
for(k in 1:d) { # deal with a_k(theta, d)
j <- k:d
## compute b_j = log(alpha^j*|s(d,j)|*S(j,k)), j = k,..,d
b <- j * log(alpha) + ls[j] +
unlist(lapply(j, function(i) lS[[i]][k]))
b.max <- max(b) # max_{j=k,..,d} b_j
exponents <- b - b.max # exponents
## compute critical sum (known to be positive)
exps <- exp(exponents) # (d-k+1)-many
even <- if(k == d) NULL else seq(2, d-k+1, by=2)
odd <- seq(1, d-k+1, by=2)
sum.neg <- sum(sort(exps[even]))
sum.pos <- sum(sort(exps[odd]))
sum. <- sum.pos - sum.neg
a[k] <- if(log) b.max + log(sum.) else exp(b.max)*sum.
if(sum.neg > sum.pos) {
if(verbose) message("sum.neg > sum.pos for k = ", k)
wrong.sign <- c(wrong.sign, k)
}
}
if(length(wrong.sign))
attr(a, "wrong.signs") <- wrong.sign
a
},
"dsSib" = {
## a_k(theta,d) = d!/k! * dsumSibuya(d, k, alpha)
k <- 1:d
ck <- ## c_k := d!/k!
if(log) c(0,cumsum(log(d:2)))[d:1]
else c(1,cumprod(d:2))[d:1]
p <- dsumSibuya(d, k, alpha, method = meth.dsSib, log=log)
##   ----------              -------------------
if(log) p + ck else p * ck
},
"horner" = {
s.abs <- abs(Stirling1.all(d))
k <- 1:d
S <- lapply(k, Stirling2.all)## S[[n]][i] contains S(n,i), i = 1,...,n
pol <- vapply(k, function(k.) {
j <- 0:(d-k.)
## compute coefficients (c_k = |s(d,j+k)|*S(j+k,k))
c.j <- s.abs[j+k.] * vapply(j, function(i) S[[i+k.]][k.], 1.)
polynEval(c.j, -alpha)
}, NA_real_)

if(log) k*log(alpha) + log(pol) else alpha^k * pol
},
"direct" = {
s <- Stirling1.all(d) # s(d,1), ..., s(d,d)
k <- 1:d
S <- lapply(k, Stirling2.all)## S[[n]][i] contains S(n,i), i = 1,...,n
vapply(k, function(k.) {
j <- k.:d
## extract a column of Stirling2 numbers:
S. <- vapply(j, function(i) S[[i]][k.], 1.)
sm <- sum(alpha^j * s[j]*S.)
if(log) log(abs(sm)) else (-1)^(d-k.)*sm
}, NA_real_)
},
stop(sprintf("unsupported method '%s' in coeffG", method))
) ## switch()
} ## coeffG()

### compute the polynomial for Gumbel

##' Compute the polynomial involved in the generator derivatives and the
##' copula density of a Gumbel copula
##'
##' @title Polynomial involved in the generator derivatives and density for Gumbel
##' @param lx = log(x); where x: evaluation point (vector);
##'        e.g., for copGumbel@dacopula, lx = alpha*log(rowSums(psiInv(u)))
##'        where u = (u_1,..,u_d) is the evaluation point of the density of Joe's copula)
##' @param alpha parameter in (0,1]   alpha := 1/theta = 1 - tau
##' @param d number of summands, >= 1
##' @param method a string, one of
##'   "default"         uses a combination of the other methods
##'   "pois.direct"     uses ppois directly
##'   "pois"            uses ppois with pulling out max
##'   "stirling"        uses the representation via Stirling numbers and once horner
##'   "stirling.horner" uses the representation via Stirling numbers and twice horner
##'    .....            all the method from coeffG(), see there
##' @param log boolean which determines if the logarithm is returned
##' @return \sum_{k=1}^d  a_{dk}(\theta)  x ^ k
##'       = \sum_{k=1}^d  a_{dk} *     exp(lx*k)
##'  where a_{dk}(theta)
##'       = (-1)^{d-k}\sum_{j=k}^d \theta^{-j} s(d,j) S(j,k)
##'       = (d!/k!)\sum_{l=1}^k (-1)^{d-l} \binom{k}{l}\binom{\alpha l}{d}
##' @author Marius Hofert
polyG <- function(lx, alpha, d, method=
c("default", "pois", "pois.direct", "stirling", "stirling.horner",
coeffG.methods),
log=FALSE)
{
stopifnot(length(alpha)==1, 0 < alpha, alpha <= 1,
d == as.integer(d), d >= 1)
k <- 1:d
## with a fixed match.arg(): method <- match.arg(method) ## -- R 2.15. ??
method <- match.arg(method, choices = eval(formals()[["method"]]))
Meth <- if(method %in% coeffG.methods) "coeffG" else method
switch(Meth,
"default" =
{
## default method: the following is clearly too simple
## ------- should even use different methods for different  x := exp(lx)

Recall(lx, alpha=alpha, d=d,
method = if(d <= 100) {
if(alpha <= 0.54) "stirling"
else if(alpha <= 0.77) "pois.direct"
else "dsSib.log"
} else "pois", # slower but more stable, e.g., for d=150
log=log)
},
"pois" =
{
## build list of b's
n <- length(lx)
x <- exp(lx)                                   # e^lx = x
lppois <- outer(d-k, x, FUN=ppois, log.p=TRUE) # a (d x n)-matrix; log(ppois(d-k, x))
llx <- k %*% t(lx)           # also a (d x n)-matrix; k*lx
labsPoch <- vapply(k, function(j) sum(log(abs(alpha*j-(k-1L)))), NA_real_) # log|(alpha*j)_d|, j=1,..,d
lfac <- lfactorial(k)        # log(j!), j=1,..,d
## build matrix of exponents
lxabs <- llx + lppois + rep(labsPoch - lfac, n) + rep(x, each = d)
res <- lssum(lxabs, signFF(alpha, k, d), strict=FALSE)
if(log) res else exp(res)
},
"pois.direct" =
{
## build coefficients
xfree <- lchoose(alpha*k,d) + lfactorial(d) - lfactorial(k)
x <- exp(lx)
lppois <- outer(d-k, x, FUN=ppois, log.p=TRUE) # (length(x),d)-matrix
klx <- lx %*% t(k)
exponents <- exp(t(x+klx)+lppois+xfree) # (d,length(x))-matrix
res <- as.vector(signFF(alpha, k, d) %*% exponents)
if(log) log(res) else res
},
"stirling" =
{
## implementation of \sum_{k=1}^d a_{dk}(\theta) x^k
## = (-1)^{d-1} * x * \sum_{k=1}^d alpha^k * s(d,k) * \sum_{j=1}^k S(k,j) * (-x)^{j-1}
## = (-1)^{d-1} * x * \sum_{k=1}^d alpha^k * s(d,k) * polynEval(...)
## inner function is evaluated via polynEval
x <- exp(lx)
s <- Stirling1.all(d) # s(d,1), ..., s(d,d)
S <- lapply(k, Stirling2.all) # S[[l]][n] contains S(l,n), n = 1,...,l
lst <- lapply(k, function(k.) (-1)^(d-1)*x*alpha^k.*s[k.]*polynEval(S[[k.]],-x))
res <- rowSums(matrix(unlist(lst), nrow=length(x)))
if(log) log(res) else res
},
"stirling.horner" =
{
## implementation of \sum_{k=1}^d a_{dk}(\theta) x^k
## = (-1)^{d-1} * x * \sum_{k=1}^d alpha^k * s(d,k) * \sum_{j=1}^k S(k,j) * (-x)^{j-1}
## = (-1)^{d-1} * x * \sum_{k=1}^d alpha^k * s(d,k) * polynEval(...)
## polynEval is used twice
x <- exp(lx)
s <- Stirling1.all(d) # s(d,1), ..., s(d,d)
S <- lapply(k, Stirling2.all) # S[[l]][n] contains S(l,n), n = 1,...,l
len <- length(x)
poly <- matrix(unlist(lapply(k, function(k.) polynEval(S[[k.]],-x))), nrow=len) # (len,d)-matrix
res <- (-1)^(d-1)*alpha*x* vapply(1:len, function(i) polynEval(s*poly[i,], alpha), 1.)
if(log) log(res) else res
## the following code was *not* faster
## poly <- t(sapply(k, function(k.) polynEval(S[[k.]],-x))) # (d,len(x))-matrix
## coeff <- if(length(x)==1) t(s*poly) else s*poly
## res <- (-1)^(d-1)*alpha*x*apply(coeff, 2, polynEval, x=alpha)
},
"coeffG" = ## <<< all the different 'coeffG' methods --------------------
{
## note: these methods are all known to show numerical deficiencies
if(d > 220) stop("d > 220 not yet supported") # would need Stirling2.all(d, log=TRUE)
## compute the log of the coefficients:
l.a.dk <- coeffG(d, alpha, method=method, log = TRUE)
##	     ~~~~~~			     ~~~~~~~~~~
## evaluate the sum
## for this, create a matrix B with (k,i)-th entry
## B[k,i] = log(a_{dk}(theta)) + k * lx[i],
##          where k in {1,..,d}, i in {1,..,n} [n = length(lx)]
logx <- l.a.dk + k %*% t(lx)
if(log) {
## compute log(colSums(exp(B))) stably (no overflow) with the idea of
## pulling out the maxima
lsum(logx)
}else colSums(exp(logx))
},
stop(sprintf("unsupported method '%s' in polyG", method))
) # end{switch}
}

### Joe ########################################################################

### sampling a Sibuya(alpha) distribution, R version

##' Sample V from a Sibuya(alpha) distribution with cdf F(n) = 1-1/(n*B(n,1-alpha)),
##' n in IN, with Laplace-Stieltjes transform 1-(1-exp(-t))^alpha via the
##' algorithm of Hofert (2011), Proposition 3.2. R version.
##'
##' @param n  sample size
##' @param alpha parameter in (0,1]
##' @return vector of random variates V
##' @author Marius Hofert, Martin Maechler
stopifnot((n <- as.integer(n)) >= 0, length(alpha)==1, 0 < alpha, alpha <= 1)
V <- numeric(n)
if(n >= 1) {
if(alpha == 1) {
V[] <- 1
} else {
u <- runif(n)
## FIXME(MM): (for alpha not too close to 1): re-express using 1-u
l1 <- u <= alpha
V[l1] <- 1
i2 <- which(!l1)
Ginv <- ((1-u[i2])*gamma(1-alpha))^(-1/alpha)
floorGinv <- floor(Ginv)
l3 <- (1-1/(floorGinv*beta(floorGinv,1-alpha)) < u[i2])
V[i2[l3]] <- ceiling(Ginv[l3])
i4 <- which(!l3)
V[i2[i4]] <- floorGinv[i4]
}
}
V
}

### state-of-the-art: sampling a Sibuya(alpha) distribution, C version

##' Sample V from a Sibuya(alpha) distribution with cdf F(n) = 1-1/(n*B(n,1-alpha)),
##' n in IN, with Laplace-Stieltjes transform 1-(1-exp(-t))^alpha via the
##' algorithm of Hofert (2011), Proposition 3.2. C version.
##'
##' @title Efficiently sampling Sibuya(alpha) distributions
##' @param n sample size (has to be numeric, >= 0)
##' @param alpha parameter in (0,1]
##' @return vector of random variates V
##' @author Martin Maechler
}

##' Probability mass function of a Sibuya(alpha) distribution
##'
##' @title Probability mass function of a Sibuya(alpha) distribution
##' @param x evaluation point [integer]
##' @param alpha parameter alpha
##' @param log boolean which determines if the logarithm is returned
##' @return p_x = choose(alpha, x) * (-1)^(x-1)
##' @author Marius Hofert and Martin Maechler
if(log) lchoose(alpha, x) else abs(choose(alpha, x))

##' Distribution function of a Sibuya(alpha) distribution
##'
##' @title Distribution function of a Sibuya(alpha) distribution
##' @param x evaluation point [integer]
##' @param alpha parameter alpha
##' @param lower.tail if TRUE, probabilities are P[X ≤ x], otherwise, P[X > x]
##' @param log.p boolean which determines if the logarithm is returned
##' @return F(x) = 1 - (-1)^x * choose(alpha-1, x)
##' @author Marius Hofert and Martin Maechler
pSibuya <- function(x, alpha, lower.tail=TRUE, log.p=FALSE)
{
## F(x) = 1 - 1/(x*Beta(x,1-alpha)) = 1 - (x*beta(x, 1-alpha))^(-1)
if(log.p) {
if(lower.tail) # log(1 - 1/(x*beta(x, 1-alpha)))
log1p(-1/(x*beta(x, 1-alpha)))
else ## log(1/(x*beta(x, 1-alpha))) = - log(x * beta(..)) =
-log(x) - lbeta(x, 1-alpha)
} else { ## no log
xb <- 1/(x*beta(x, 1-alpha))
if(lower.tail) 1 - xb else xb
}
}

### state-of-the art: sampling F01Joe, C version

##' Generate a vector of variates V01 ~ F01 with Laplace-Stieltjes transform
##' ((1-(1-exp(-t))^alpha))^V0. Bridge to R. Used, e.g., to draw several variates
##' from rF01Joe.
##'
##' @title Sampling F01 for Joe's family
##' @param V0 vector of random variates from F0
##' @param parameter alpha = theta0/theta1 in (0,1]
##' @param approx largest number of summands before asymptotics is used
##' @return vector of random variates V01
##' @author Marius Hofert
rF01Joe <- function(V0, alpha, approx) {
.Call(rF01Joe_vec_c, V0, alpha, approx)
}

### wrapper for inner distribution F for Joe

##' Generate a vector of variates V ~ F with Laplace-Stieltjes transform
##' 1-(1-exp(-t))^alpha.
##'
##' @title Sampling F for Joe
##' @param n number of variates from F
##' @param parameter alpha = theta0/theta1 in (0,1]
##' @return vector of random variates V
##' @author Marius Hofert
rFJoe <- function(n, alpha) rSibuya(n, alpha)

### polynomial evaluation for Joe

##' Inner probability mass function for a nested Joe copula, i.e. a Sibuya sum
##' also used in coeffG() -> polyG() for Gumbel's copula
##'
##' @title Inner probability mass function for a nested Joe copula
##' @param x vector (or number) of natural numbers
##' @param n vector (or number) of natural numbers
##' @param alpha parameter in (0,1]
##' @param method method applied
##'        log:      proper log computation based on lssum
##'        direct:   brute-force evaluation of the sum and its log
##'        Rmpfr, RmpfrM:    multi-precision; the latter *return* multi-prec.
##'        diff:     via forward differences
##'        exp.log:  similar to method = "log", but without *proper/intelligent* log
##' @param log boolean which determines if the logarithm is returned
##' @return p_{x,n} = \sum_{j=1}^n choose(n,j)*choose(alpha*j,x)*(-1)^(x-j)
##'         which is a probability mass function in x on IN with generating function
##'         g(z) = (1-(1-z)^alpha)^n
##' @author Marius Hofert and Martin Maechler
##' note: - p_{x,n} = 0 for x < n;  p_{n,n} = alpha^n
##'       - numerically challenging, e.g., dsumSibuya(100, 96, 0.01) < 0 for all methods
method= c("log", "direct", "Rmpfr", "Rmpfr0", "RmpfrM", "diff", "exp.log"),
log=FALSE)
{
stopifnot(x == round(x), n == round(n), n >= 0, length(alpha) == 1,
0 < alpha, alpha <= 1)
if((l.x <- length(x)) == 0 || (l.n <- length(n)) == 0)
return(numeric())
## "FIXME": from coeffG(), always have  {x = d, n = 1:d} -- can be faster *not* recycling
if((len <- l.x) != l.n) { ## do recycle to common length
len <- max(l.x, l.n)
if(l.x < len) {
x. <- x
x <- rep(x, length.out = len)
}
else ## if(l.n < len)
n <- rep(n, length.out = len)
}
if(alpha == 1)
return(x == n)
method <- if(missing(method) && is(alpha, "mpfr"))
"Rmpfr" else match.arg(method)
switch(method,
"log" =
{
## computes *proper* log based on lssum

## determine the matrix of signs of choose(alpha*j,x)*(-1)^(x-j),
## j in {1,..,m} -- which notably do *not* depend on x !
signs <- signFF(alpha, seq_len(max(n)), x)
## for one pair of x and n:
f.one <- function(x,n) {
if(x < n) return(-Inf)	# = log(0)
j <- seq_len(n)
lxabs <- lchoose(n, j) + lchoose(alpha*j, x)
## *NON*-strict -- otherwise need try() :
lssum(as.matrix(lxabs), signs[j], strict=FALSE)
}
S <- mapply(f.one, x, n, USE.NAMES = FALSE)
if(log) S else exp(S)
},
"direct" =
{
## brute force evaluation of the sum and its log
f.one <- function(x,n) {
if(x < n) return(0)
j <- seq_len(n)
sum(choose(n,j)*choose(alpha*j,x)*(-1)^(x-j))
}
S <- mapply(f.one, x, n, USE.NAMES = FALSE)
if(log) log(S) else S
},
"Rmpfr" =,
"Rmpfr0" =,
"RmpfrM" =
{
## as "direct" but using high-precision arithmetic, where
## the precision should be set via alpha = mpfr(*, precBits= .)
stopifnot(require(Rmpfr))
if(!is(alpha, "mpfr"))
alpha <- mpfr(alpha, precBits = max(100, min(x, 10000)))
pr. <- getPrec(alpha)
mpfr.0 <- mpfr(0, precBits = pr.)
mayRecall <- (method != "Rmpfr0")
## FIXME
if(FALSE && l.x == 1 && l.n == x. && all(n == seq_len(x.))) { ## Special case -- from coeffG()
message("fast special case ..") ## <- just for now
## change notation (for now)
j <- n  # == 1:d
n <- x. # == d
c.n <- chooseMpfr.all(n)
ca.j <- chooseMpfr(alpha*j,n)*(-1)^(n-j)
f.sp <- function(j) {
j. <- seq_len(j)
sum(c.n[j.] * ca.j[j.])
}
stop("fast special case -- is still wrong !")
S <- new("mpfr", unlist(lapply(j, f.sp)))
} else { ## usual case
f.one <- function(x,n) {
if(x < n) return(mpfr.0)
j <- seq_len(n)
## now do this in parts {to analyze}:
## sum(chooseMpfr.all(n)*chooseMpfr(alpha*j,x)*(-1)^(x-j))
c.n <- chooseMpfr.all(n)
ca.j <- chooseMpfr(alpha*j,x) * (-1)^(x-j)
trms <- c.n * ca.j
S <- sum(trms) # <-- sum of *huge* (partly alternating) terms getting almost 0
if(mayRecall) {
recall <- (S <= 0)
if(recall) { ## complete loss of precision -- recall with higher precision
newPrec <- round(1.5 * pr.)
MSG <- " |--> S < 0"
} else {		# S > 0 :
bitLoss <- round(as.numeric(log2(max(abs(trms)) / S )), 1)
recall <- (pr. - bitLoss < 20) ## less than 20 bits precision left
if(recall) {
newPrec <- pr. + 30
MSG <- paste("-> bit loss =", bitLoss)
}
}
if(recall) {
message(sprintf("dsumSibuya(%d, %d, alpha= %g, method='%s'): %s ==> recalled with new prec. %d",
x, n, as.numeric(alpha), method, MSG, newPrec))
method="RmpfrM", log=FALSE)
} else S
} else S
}
S <- new("mpfr", mapply(f.one, x, n, USE.NAMES=FALSE))
}
if(method == "RmpfrM")
if(log) log(S) else S
else
as.numeric(if(log) log(S) else S)
},
"diff" =
{
S. <- mapply(function(x,n) diff(choose(n:0*alpha, x), differences=n) * (-1)^x,
x, n, USE.NAMES = FALSE)
if(log) log(S.) else S.
},
"exp.log" =
{
## similar to method = "log", but without *proper/intelligent* log
## and inefficient due to the signs (old version)
f.one <- function(x,n) {
if(x < n) return(0)
j <- seq_len(n) ## indices of the summands
signs <- (-1)^(j+x)
## determine the signs of choose(j*alpha,x) for each component of j
to.subtract <- 0:(x-1)
sig.choose <-
vapply(j, function(l) prod(sign(l*alpha-to.subtract)),
1., USE.NAMES=FALSE)
signs <- signs*sig.choose
binom.coeffs <- exp(lchoose(n,j) + lchoose(j*alpha,x))
sum(signs*binom.coeffs)
}
S <- mapply(f.one, x, n, USE.NAMES = FALSE)
if(log) log(S) else S
},
## otherwise
stop(sprintf("unsupported method '%s' in dsumSibuya", method)))
}

### polynomial evaluation for Joe

##' Compute the polynomial involved in the generator derivatives and the
##' copula density of a Joe copula
##'
##' @title Polynomial involved in the generator derivatives and density for Joe
##' @param lx (log) evaluation point (lx is meant to be log(x) for some x which
##'        was used earlier; e.g., for copJoe@dacopula, lx = log(h(u)/(1-h(u))) for
##'        h(u) = \prod_{j=1}^d(1-(1-u_j)^theta), where u = (u_1,..,u_d) is the
##'        evaluation point of the density of Joe's copula)
##' @param alpha parameter alpha ( := 1/theta ) in (0,1]
##' @param d number of summands
##' @param method different methods, can be
##'        "log.poly" intelligent log version
##'        "poly"     brute force log version
##' @param log boolean which determines if the logarithm is returned
##' @return \sum_{k=1}^d a_{d,k}(theta) exp((k-1)*lx) = \sum_{k=0}^{d-1} a_{d,k+1}(theta) x^k
##'         where a_{d,k}(theta) = S(d,k)*(k-1-alpha)_{k-1} = S(d,k)*Gamma((1:d)-alpha)/Gamma(1-alpha)
##' @author Marius Hofert and Martin Maechler
polyJ <- function(lx, alpha, d, method=c("log.poly","log1p","poly"), log=FALSE) {
stopifnot(length(alpha)==1, 0 < alpha, alpha <= 1)
## compute the log of the coefficients a_{dk}(theta)
if(d > 220) stop("d > 220 not yet supported")# would need Stirling2.all(d, log=TRUE)
k <- 1:d
l.a.k <- log(Stirling2.all(d)) + lgamma(k-alpha) - lgamma(1-alpha) # log(a_{dk}(theta)), k = 1,..,d
## evaluate polynomial via exp( log(<poly>) )
## for this, create a matrix B with (k,i)-th entry B[k,i] = log(a_{dk}(theta)) + (k-1) * lx[i],
## where k in {1,..,d}, i in {1,..,n} [n = length(lx)]
B <- l.a.k + (k-1) %*% t(lx)
method <- match.arg(method)
switch(method,
"log.poly" = {
## stably compute log(colSums(exp(B))) (no overflow)
## Idea:
## (1) let b_k := log(a_{dk}(theta)) + (k-1)*lx and b_{max} := argmax{b_k}.
## (2) \sum_{k=1}^d a_{dk}(theta)\exp((k-1)*lx) = \sum_{k=1}^d \exp(log(a_{dk}(theta))
##     + (k-1)*lx) = \sum_{k=1}^d \exp(b_k) = \exp(b_{max})*\sum_{k=1}^d
##     \exp(b_k-b_{max})
## (3) => log(\sum...) = b_{max} + log(\sum_{k=1}^d \exp(b_k-b_{max}))
if(log) lsum(B) else exp(lsum(B))
},
"log1p" = {
## use log(1 + sum(<smaller>)) = log1p(sum(<smaller>)),
## but we don't expect it to make a difference
im <- apply(B, 2, which.max) # indices (vector) of maxima
n <- length(lx) ; d1 <- d-1L
max.B <- B[cbind(im, seq_len(n))] # get max(B[,i])_{i=1,..,n} == apply(B, 2, max)
B.wo.max <- matrix(B[unlist(lapply(im, function(j) k[-j])) +
d*rep(0:(n-1), each = d1)], d1, n) # matrix B without maxima
res <- max.B + log1p(colSums(exp(B.wo.max - rep(max.B, each = d1))))
if(log) res else exp(res)
},
"poly" = {
## brute force ansatz
res <- colSums(exp(B))
if(log) log(res) else res
},
stop(sprintf("unsupported method '%s' in polyJ", method)))
}

##' Circular/Rational function	(1 - x^d)/(1 - x) for x ~~ 1, i.e.,
##' compute (1 - x^d)/(1 - x) = (1 - (1-e)^d) / e   for	 e = 1-x (<< 1) and integer d
##'
##' @title Circular/Rational function  (1 - (1-e)^d) / e  {incl. limit e -> 0}
##' @param e numeric vector in [0, 1]
##' @param d integer (scalar), >= 1
##' @return (1 - (1-e)^d) / e
##' @author Martin Maechler, Date: 25 Jul 2011
circRat <- function(e, d)
{
### TODO (?):  improve "log=TRUE", for e ~= 1:	log(circRat(e, d)) = log1p(-x^d) - log(e)
stopifnot(length(d) == 1, d == as.integer(d), d >= 1)
if(d <= 6)
switch(d,
1-0*e, ## <<- d = 1
2-e,   ## <<- d = 2: 1 2 1
3-e*(3-e),#   d = 3: 1 3 3 1
4-e*(6-e*(4-e)),		       ## d = 4: 1 4 6 4 1
5-e*(10-e*(10-e*(5-e))),	       ## d = 5: 1 5 10 10 5 1
6-e*(15-e*(20-e*(15-e*(6-e))))  ## d = 6: 1 6 15 20 15 6 1
)
else { ## d >= 7 ---------------
r <- e
eps <- .Machine$double.eps if(any(l1 <- ((d1 <- (d-1)/2)*e < eps))) r[l1] <- d if(any(l2 <- !l1 & ((d2 <- (d-2)/3)*(e2 <- e*e) < eps))) r[l2] <- d*(1 - d1*e[l2]) if(any(l3 <- !l1 & !l2 & ((d3 <- (d-3)/4)*e*e2 < eps))) r[l3] <- d*(1 - d1*e[l3]*(1 - d2*e[l3])) ## and for the remaining ones, we afford a little precision loss: if(any(lrg <- !l1 & !l2 & !l3)) { e <- e[lrg] r[lrg] <- (1 - (1-e)^d)/e } r } } ### Misc ####################################################################### ##' Conditional copula function C(u[,d]|u[,1],...,u[,d-1]) ##' ##' @title Conditional copula function ##' @param u (n x d)-matrix of evaluation points (first d-1 colums are conditioned on) ##' @param cop an outer_nacopula ##' @param n.MC Monte Carlo sample size ##' @param log if TRUE the logarithm of the conditional copula is returned ##' @author Marius Hofert cacopula <- function(u, cop, n.MC=0, log=FALSE) { stopifnot(is(cop, "outer_nacopula")) if(length(cop@childCops)) stop("currently, only Archimedean copulas are provided") if(!is.matrix(u)) u <- rbind(u) stopifnot(0 <= u, u <= 1) ## note: for some (but not all) families, cacopula also makes sense on ## the boundaries since the corresponding limits exist th <- cop@copula@theta stopifnot(cop@copula@paraConstr(th)) dim. <- dim(u) n <- dim.[1] d <- dim.[2] psiI <- cop@copula@psiInv(u, theta=th) arg.denom <- rowSums(psiI[,1:(d-1), drop=FALSE]) arg.num <- arg.denom + psiI[,d] logD <- cop@copula@psiDabs(c(arg.num, arg.denom), theta=th, degree=d-1, n.MC=n.MC, log=TRUE) res <- logD[1:n]-logD[(n+1):(2*n)] if(log) res else exp(res) } ##' Function which computes psiDabs via Monte Carlo ##' ##' @title Computing the absolute value of the generator derivatives via Monte Carlo ##' @param t evaluation points ##' @param family Archimedean family (name or object) ##' @param theta parameter value ##' @param degree order of derivative ##' @param n.MC Monte Carlo sample size ##' @param method different methods ##' log: proper log using lsum ##' direct: direct evaluation of the sum ##' pois.direct: directly uses the Poisson density ##' pois: intelligently uses the Poisson density with lsum ##' @param log if TRUE the log of psiDabs is returned ##' @author Marius Hofert ##' Note: psiDabsMC(0) is always finite, although, theoretically, psiDabs(0) may ##' be Inf (e.g., for Gumbel and Joe) psiDabsMC <- function(t, family, theta, degree=1, n.MC, method=c("log", "direct", "pois.direct", "pois"), log = FALSE)# not yet: is.log.t = FALSE) { res <- numeric(length(t)) V <- getAcop(family)@V0(n.MC, theta) method <- match.arg(method) switch(method, ## the following is not faster than "log": ## "default" = { # basically, use "direct" if numerically not critical and "log" otherwise ## lx <- -V %*% t(t) + degree*log(V) ## explx <- exp(lx) # (n.MC, n)-matrix containing the summands ## explx0 <- explx==0 # can be TRUE due to t == Inf or t finite but too large ## t.too.large <- unlist(lapply(1:n, function(x) any(explx0))) # boolean vector of length n indicating which column of explx contains zeros ## r1 <- colMeans(explx[,!t.too.large, drop=FALSE]) ## res[!t.too.large] <- if(log) log(r1) else r1 ## r2 <- lsum(lx[,t.too.large, drop=FALSE] - log(n.MC)) ## res[t.too.large] <- if(log) r2 else exp(r2) ## res[is.infinite(t)] <- if(log) -Inf else 0 ## res ## }, "log" = { # intelligent log iInf <- is.infinite(t) res[iInf] <- -Inf # log(0) if(any(!iInf)) res[!iInf] <- lsum(-V %*% t(t[!iInf]) + degree*log(V) - log(n.MC)) if(log) res else exp(res) }, "direct" = { # direct method lx <- -V %*% t(t) + degree*log(V) res <- colMeans(exp(lx)) # can be all zeros if lx is too small [e.g., if t is too large] if(log) log(res) else res }, "pois.direct" = { poi <- dpois(degree, lambda=V %*% t(t)) res <- factorial(degree)/t^degree * colMeans(poi) if(log) log(res) else res }, "pois" = { iInf <- is.infinite(t) res[iInf] <- -Inf # log(0) if(any(!iInf)) { t. <- t[!iInf] lpoi <- dpois(degree, lambda=V %*% t(t.), log=TRUE) # (n.MC, length(t.))-matrix b <- -log(n.MC) + lfactorial(degree) - degree*rep(log(t.), each=n.MC) + lpoi # (n.MC, length(t.))-matrix res[!iInf] <- lsum(b) } if(log) res else exp(res) }, stop(sprintf("unsupported method '%s' in psiDabsMC", method))) } ### Non-numerics ############################################################### ##' Function for setting the parameter in an acopula ##' ##' @title Settting the parameter in an acopula ##' @param x acopula ##' @param value parameter value ##' @param na.ok logical indicating if NA values are ok for theta ##' @return acopula with theta set to value ##' @author Martin Maechler setTheta <- function(x, value, na.ok = TRUE) { stopifnot(is(x, "acopula"), is.numeric(value) | (ina <- is.na(value))) if(ina) { if(!na.ok) stop("NA value, but 'na.ok' is not TRUE") value <- NA_real_ } if(ina || x@paraConstr(value)) ## parameter constraints are fulfilled x@theta <- value else stop("theta (=", format(value), ") does not fulfill paraConstr()") x } ##' Construct "paraConstr" function from an "interval" ##' ##' @title Construct "paraConstr" function from an "interval" ##' @param int interval ##' @return parameter constraint function ##' @author Martin Maechler mkParaConstr <- function(int) { stopifnot(is(int, "interval")) # for now is.o <- int@open eL <- substitute(LL <= theta, list(LL = int[1])); if(is.o[1]) eL[[1]] <- as.symbol("<") eR <- substitute(theta <= RR, list(RR = int[2])); if(is.o[2]) eR[[1]] <- as.symbol("<") bod <- substitute(length(theta) == 1 && LEFT && RIGHT, list(LEFT = eL, RIGHT= eR)) as.function(c(alist(theta=), bod), parent.env(environment())) ## which is a fast version of ## r <- function(theta) {} ## environment(r) <- parent.env(environment()) ## body(r) <- bod ## r } printAcopula <- function(x, slots = TRUE, indent = 0, digits = getOption("digits"), width = getOption("width"), ...) { cl <- class(x) cld <- getClassDef(cl) stopifnot(indent >= 0, extends(cld, "acopula")) ch.thet <- { if(!all(is.na(x@theta)))## show theta paste(", theta= (", paste(sapply(x@theta, format, digits=digits), collapse=", "), ")", sep="") else "" } bl <- paste(rep.int(" ",indent), collapse="") cat(sprintf('%sArchimedean copula ("%s"), family "%s"%s\n', bl, cl, x@name, ch.thet)) if(slots) { nms <- slotNames(cld) nms <- nms[!(nms %in% c("name", "theta"))] i2 <- indent+2 cat(bl, " It contains further slots, named\n", paste(strwrap(paste(dQuote(nms),collapse=", "), width = 0.95 * (width-2), indent=i2, exdent=i2), collapse="\n"), "\n", sep="") } invisible(x) } setMethod(show, "acopula", function(object) printAcopula(object)) ## This is now exported & has help file --> ../man/printNacopula.Rd : printNacopula <- function(x, labelKids = NA, deltaInd = if(identical(labelKids,FALSE)) 5 else 3, indent.str="", digits = getOption("digits"), width = getOption("width"), ...) { cl <- class(x) stopifnot(deltaInd >= 0, is.character(indent.str), length(indent.str) == 1, extends(cl, "nacopula")) mkBlanks <- function(n) paste(rep.int(" ", n), collapse="") bl <- mkBlanks(nIS <- nchar(indent.str)) ## cat(sprintf(" __deltaInd = %d, nIS = %d__ ", deltaInd, nIS)) ch1 <- sprintf("%sNested Archimedean copula (\"%s\"), with ", indent.str, cl) ch2 <- if(length(c.j <- x@comp)) { sprintf("slot \n%s'comp' = %s", bl, paste("(",paste(c.j, collapse=", "),")", sep="")) } else "empty slot 'comp'" cat(ch1, ch2, sprintf(" and root\n%s'copula' = ", bl), sep="") printAcopula(x@copula, slots=FALSE, digits=digits, width=width, ...) nk <- length(kids <- x@childCops) if(nk) { cat(sprintf("%sand %d child copula%s\n", bl, nk, if(nk > 1)"s" else "")) doLab <- if(is.na(labelKids)) nk > 1 else as.logical(labelKids) if(doLab) { hasNms <- !is.null(nms <- names(kids)) lab <- if(hasNms) paste0(nms,": ") else paste0(seq_len(nk),") ") } bl <- mkBlanks(nIS + deltaInd) for(ic in seq_along(kids)) printNacopula(kids[[ic]], deltaInd=deltaInd, indent.str = paste0(bl, if(doLab) lab[ic]), labelKids=labelKids, digits=digits, width=width, ...) } else cat(sprintf("%sand *no* child copulas\n", bl)) invisible(x) } setMethod(show, "nacopula", function(object) printNacopula(object)) ##' Get one of our "acopula" family objects by name ##' ##' @title Get one of our "acopula" family objects by name ##' @param family either character string (short or longer form of ##' copula family name) or an "acopula" family object ##' @param check logical indicating if the class of the return value should ##' be checked. ##' @return one of our "acopula" objects ##' @author Martin Maechler getAcop <- function(family, check=TRUE) { if(is.character(family)) { stopifnot(length(family) == 1) if(nchar(family) <= 2) # it's a short name family <- c_longNames[family] COP <- get(c_objNames[family]) # envir = "package:nacopula" if(check && !is(COP, "acopula")) stop(paste("invalid acopula-family object, family=",family)) COP } else if(is(family, "acopula")) family else stop("'family' must be an \"acopula\" object or family name") } coeffG.methods <- eval(formals(coeffG)$method)# - namespace hidden
## --> accesses formals(dsumSibuya) .. hence at end