Nothing
## Copyright (C) 2015-2022 Leo Belzile
##
## This file is part of the "lcopula" package for R. 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 2 of the License, or (at your
## option) any later version.
##
## This library 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, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston,
## MA 02111-1307 USA or look up the web page
## http://www.gnu.org/copyleft/gpl.html.
##########################################################################
## Archimedean Liouville generation and estimation ##
## Version of July 05th, 2016 ##
## Built - 0.206 ##
##########################################################################
# Code adapted from GumbelLiouville.R (Neslehova and Genest - 2013)
# and from LiouvilleFunction.R (McNeil and Neslehova - 2010)
# using functions implemented by Hofert, Machler and McNeil (2012)
# in the package copula
# Warning: different definition used by latter for the Clayton
# The implementation below uses the generator definition found in QRM
# The Archimedean families implemented are Clayton, Gumbel, Frank,
# Ali-Mikhail-Haq (abbreviated AMH) and Joe
# require(copula)
##########################################################################
## Necessary global imports for the methods for NAMESPACE
## Note that stats cannot be imported because it is already loaded by copula
#' @useDynLib lcopula
#' @exportPattern "^[[:alpha:]]+"
#' @importFrom graphics "axis" "lines" "mtext" "plot.window" "points"
#' @importFrom stats "integrate" "optimise" "pbeta" "quantile" "rgamma" "runif" "uniroot" "cor"
#' @import copula
#' @importFrom utils combn
#' @importFrom Rcpp evalCpp
NULL
#' Archimedean copula sampler
#'
#' Sampler based on the Marshall-Olkin algorithm
#'
#' @param n sample size
#' @param family family of the Archimedean copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param d dimension of sample
#' @param theta parameter of the Archimedean copula
#' @return a sample of dimension \code{n} by \code{d} from the Archimedean copula
#' @export
#' @examples
#' #Sample from a Gumbel Archimedean copula
#' rarchi(n = 100, "gumbel", d = 4, theta = 2)
#' #Sample from the independence copula
#' rarchi(n = 100, "gumbel", d = 4, theta = 1)
rarchi <- function(n, family, d, theta){
family <- match.arg(family, c("clayton","gumbel","frank","AMH","joe"))
#Marshall-Olkin algorithm
illegalpar <- switch(family,
clayton = copClayton@paraConstr(theta),
gumbel = copGumbel@paraConstr(theta),
frank = copFrank@paraConstr(theta),
AMH = copAMH@paraConstr(theta),
joe = copJoe@paraConstr(theta))
if(!illegalpar){
stop("Illegal parameter value")
}
independence <- switch(family,
clayton = (theta == 0),
gumbel = (theta == 1),
frank = (theta == 0),
AMH = (theta == 0),
joe = (theta == 1))
U <- runif(n * d)
U <- matrix(U, nrow = n, ncol = d)
if(independence)
return(U)
Y <- switch(family,
clayton = rgamma(n, shape = 1/theta, scale = theta),
gumbel = copGumbel@V0(n, theta),
frank = copFrank@V0(n, theta),
AMH = copAMH@V0(n, theta),
joe = copJoe@V0(n, theta))
Y <- matrix(Y, nrow = n, ncol = d)
phi.inverse <- switch(family,
clayton = function(t, theta){(1 + theta*t)^(-1/theta)},
gumbel = copGumbel@psi,
frank = copFrank@psi,
AMH = copAMH@psi,
joe = copJoe@psi)
return(phi.inverse( - log(U)/Y, theta))
}
##########################################################################
#' Liouville copulas
#'
#' Multivariate density, survival copula
#' and random generation for the Liouville copulas.
#'
#' \code{rliouv} generates draws from the Liouville copula. \code{dliouv} evaluates the density of an \code{n} by \code{d} matrix of observations. \code{pliouv} is the (survival) copula associated with the Liouville vector and is as such the multivariate distribution function for uniform observations.
#'
#' Liouville copulas were introduced in McNeil and Neslehova (2010), generalizing Archimedean copulas. Like the latter, they
#' are survival copulas, which means that the copula is evaluated using the (multivariate) survival function of Liouville vectors. See also \code{\link{sliouv}} for the latter.
#' @references McNeil A.J. and Neslehova, J.G. (2010) From Archimedean to Liouville Copulas. \emph{J. Multivar. Anal.}, \strong{101}(8): 1772--1790.
#' @name Liouville
#' @aliases rliouv pliouv dliouv
#' @param n sample size
#' @param x matrix of quantiles from a Liouville copula
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers) Specifies (implictly) the dimension of sample
#' @param theta parameter of the corresponding Archimedean copula
#' @param reverse if \code{TRUE}, return sample from the corresponding survival copula
#' @param is.log if \code{TRUE}, will return the log-likelihood value
#' @return either a matrix of dimension \code{n} by \code{length(alphavec)} with the corresponding quantile, probability, survival probability or sample from the Liouville vector
#' @export
#' @seealso \code{\link{Liouville_marginal}}
#' @examples
#' \dontrun{
#' #Multivariate density of Clayton Liouville copula
#' x <- rliouv(n = 100, family = "clayton", alphavec <- c(2,3), theta = 2)
#' dliouv(x=x, family="clayton", alphavec=c(2,3), theta=2, TRUE)
#' #Distribution function, multivariate sample
#' x <- rliouv(n=100, family="frank", theta=1.5, alphavec=c(2,3))
#' pliouv(theta=1.5, x=x,family="frank", alphavec=c(2,3))
#' }
rliouv <- function(n = 100, family, alphavec, theta, reverse = FALSE){
family <- match.arg(family, c("clayton","gumbel","frank","AMH","joe"))
alphavec <- as.integer(alphavec)
if(isTRUE(any(alphavec==0))){stop("Invalid parameters")}
alpha <- sum(alphavec)
d <- length(alphavec)
illegalpar <- switch(family,
clayton = copClayton@paraConstr(theta),
gumbel = copGumbel@paraConstr(theta),
frank = copFrank@paraConstr(theta),
AMH = copAMH@paraConstr(theta),
joe = copJoe@paraConstr(theta))
if(!illegalpar)
stop("Illegal parameter value")
t <- rarchi(n, family, alpha, theta)
Xtildedata <- switch(family,
clayton = (t^(-theta)-1)/theta,
gumbel = copGumbel@iPsi(t, theta),
frank = copFrank@iPsi(t, theta),
AMH = copAMH@iPsi(t, theta),
joe = copJoe@iPsi(t, theta))
groups <- c(0, cumsum(alphavec))
Xdata <- matrix(0, nrow = n, ncol = d)
Udata <- Xdata
for (j in 1:d){
class <- ((groups[j]+1):groups[j+1])
for (i in class)
Xdata[, j] <- Xdata[, j]+Xtildedata[, i]
}
for (j in 1:d){
Udata[, j] <- sliouvm(Xdata[, j], family, alphavec[j], theta)
}
if(reverse == FALSE){
return(Udata)
} else{
return(1-Udata)
}
}
##########################################################################
## Marginal functions for Liouville vectors
# Default behavior as in original GumbelLiouville.R
#' @name Liouville_marginal
#' @aliases pliouvm dliouvm isliouvm sliouvm
#' @title Liouville vectors marginal functions
#' @description Marginal density, distribution, survival and inverse survival functions for Liouville copulas or Liouville vectors.
#' The inverse survival function of Liouville vectors is not available in closed-form and is obtained numerically by root-finding.
#' As such, Monte-Carlo approximation have been considered for dealing with inference to avoid computational bottlenecks.
#' Note: the arguments of \code{sliouv} are reversed since they are meant to be called inside \code{optim}. The functions borrow
#' \eqn{psi} functions and their derivatives from the \code{\link[copula]{copula-package}}.
#'
#' @param u vector of quantiles or survival probabilities, (pseudo)-uniform variates
#' @param x vector of quantiles from a Liouville copula (or a Liouville vector for the survival function , with support on the positive real line)
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param alpha integer Dirichlet parameter
#' @param theta parameter of the corresponding Archimedean copula
#' @return a vector with the corresponding quantile, probability, survival probabilities
#' @export
#' @examples
#' \dontrun{
#' #Marginal density
#' samp <- rliouv(n = 100, family = "clayton", alphavec <- c(2,3), theta = 2)
#' dliouvm(x=samp[,1], family="clayton", alpha=2, theta=2)
#' sum(log(dliouvm(x=samp[,1], family="clayton", alpha=2, theta=2)))
#' #Marginal distribution and (inverse) survival function
#' x <- rliouv(n = 100, family = "gumbel", alphavec <- c(2,3), theta = 2)
#' pliouvm(x[,1], family="gumbel", alpha=alphavec[1], theta=2)
#' su <- sliouvm(1-x[,1], family="gumbel", alpha=alphavec[1], theta=2)
#' isliouvm(u=su, family="clayton", alpha=2, theta=2)
#' #pliouv is the same as sliouv(isliouvm)
#' }
#' @rdname Liouville_marginal
#' @export
sliouvm <- function(x, family, alpha, theta){
family <- match.arg(family, c("clayton","gumbel","frank","AMH","joe"))
alpha <- as.integer(alpha)
if(alpha==0){stop("Invalid parameter alpha")}
out <- switch(family,
gumbel = copGumbel@psi(t = x, theta),
clayton = copClayton@psi(t = x*theta, theta),
frank = copFrank@psi(t = x, theta),
AMH = copAMH@psi(t = x, theta),
joe = copJoe@psi(t = x, theta))
if(alpha>1){
for(i in 1:(alpha-1)){
out <- out + exp(i*log(ifelse(is.infinite(x),1e40,x))-lgamma(i+1)+
switch(family,
gumbel = copGumbel@absdPsi(x, theta = theta, degree = i, log = TRUE),
clayton = copClayton@absdPsi(x*theta, theta = theta, degree = i, log = TRUE)+i*log(theta),
frank = copFrank@absdPsi(x, theta = theta, degree = i, log = TRUE),
AMH = copAMH@absdPsi(x, theta = theta, degree = i, log = TRUE),
joe = copJoe@absdPsi(x, theta = theta, degree = i, log = TRUE))
)
}
}
return(out)
}
##########################################################################
## Marginal distribution function for Liouville vectors
# Calls \code{sliouvm} to return cdf
#'@rdname Liouville_marginal
#'@export
pliouvm <- function(x, family, alpha, theta){
family <- match.arg(family, c("clayton","gumbel","frank","AMH","joe"))
return(1-sliouvm(x, family, alpha, theta))
}
##########################################################################
#' Multiple marginal survival function for Liouville vectors
#'
#' This function is a wrapper around \code{\link{sliouv}}; it allows the user to treat all the data matrix simultaneously
#' by applying different parameters to each margin.
#'
#' @param x sample from copula
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#' @param theta parameter of the corresponding Archimedean copula
#' @examples
#' x <- rliouv(n = 100, family = "gumbel", alphavec <- c(2,3), theta = 2)
#' sliouv_m(x, family="gumbel", alphavec=c(2,3), theta=2)
#' all(sliouv_m(x, family="gumbel", alphavec=c(2,3), theta=2)[,1]-
#' sliouvm(x[,1], family="gumbel", alpha=2, theta=2)==0)
#' @return a matrix of same length as \code{x} with the survival probabilities
sliouv_m <- function(x, family, alphavec, theta){
alphavec <- as.integer(alphavec)
if(any(alphavec)==0){stop("Invalid parameters")}
family <- match.arg(family, c("clayton","gumbel","frank","AMH","joe"))
eval.Surv <- function(x, alpha){
der <- matrix(nrow = alpha, ncol = length(x))
der[1, ] <- switch(family,
gumbel = copGumbel@psi(t = x, theta = theta),
clayton = copClayton@psi(t = x*theta, theta = theta),
frank = copFrank@psi(t = x, theta = theta),
AMH = copAMH@psi(t = x, theta = theta),
joe = copJoe@psi(t = x, theta = theta))
if(alpha > 1){
for(i in 1:(alpha-1)){
der[i+1, ] <- exp(i*log(x)-lgamma(i+1)+switch(family,
gumbel = copGumbel@absdPsi(x, theta = theta, degree = i, log = TRUE),
clayton = copClayton@absdPsi(x*theta, theta = theta, degree = i, log = TRUE)+i*log(theta),
frank = copFrank@absdPsi(x, theta = theta, degree = i, log = TRUE),
AMH = copAMH@absdPsi(x, theta = theta, degree = i, log = TRUE),
joe = copJoe@absdPsi(x, theta = theta, degree = i, log = TRUE)))
}
}
return(colSums(der))
}
if(is.matrix(x)){
if(dim(x)[2] != length(alphavec))
stop("Length of the parameter vector does not match input")
result <- list()
for(j in 1:(dim(x)[2])){
result[[j]] <- eval.Surv(x[, j], alphavec[j])
}
result <- unlist(result)
} else{
result <- eval.Surv(x, alphavec[1])
}
return(matrix(result, nrow = nrow(x), ncol = ncol(x)))
}
##########################################################################
## Marginal inverse survival function for Liouville vectors
# Survival analog to \code{'q'} in R for functions
#'@rdname Liouville_marginal
#'@export
isliouvm <- function(u, family, alpha, theta){
alpha <- as.integer(alpha)
if(length(alpha)>1){warning("`alpha' must be an integer. Using only the first argument")}
alpha <- alpha[1]
if(alpha==0){stop("Invalid parameter alpha")}
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
out <- rep(0, length(u))
rootfunc <- function(z, u, family, theta, alpha) {
sliouvm(exp(z), family, alpha, theta)-u
}
for (i in 1:length(u)) {
#Handling of 0 for the Gumbel
tmp <- uniroot(rootfunc, lower = -1e60, upper = 1e+60, u = u[i], f.lower=1-u[i], f.upper=-u[i],
#f.lower = 1, f.upper = -1e+60,
extendInt="downX", family = family, theta = theta, alpha = alpha)
if(tmp$root > 1e58){#Beyond numerical precision, probably not needed
tmp$root <- Inf
} else if(tmp$root < -1e40){
tmp$root <- -Inf
}
out[i] <- exp(tmp$root)
}
out
}
##########################################################################
#' Multiple marginal inverse survival function of Liouville vectors
#'
#' This function is a wrapper around \code{\link{isliouvm}}; it allows the user to treat all the data matrix simultaneously
#' by applying different parameters to each margin.
#' @param u vector of survival probabilities
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#' @param theta parameter of the corresponding Archimedean copula
#' @return a vector of same length as \code{u} with the quantile at 1-u
#' @examples
#' u <- rliouv(n = 10, family = "clayton", alphavec <- c(2,3), theta = 2)
#' isliouv_m(u=u, family="clayton", alphavec=c(2,3), theta=2)
isliouv_m <- function(u, family, alphavec, theta){
alphavec <- as.integer(alphavec)
if(isTRUE(any(alphavec==0))){stop("Invalid parameters")}
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
out <- matrix(nrow = nrow(u), ncol = ncol(u))
rootfunc <- function(z, u, family, theta, alpha) {
sliouvm(exp(z), family, alpha, theta)-u
}
for (i in 1:nrow(u)){
for(j in 1:ncol(u)){
tmp <- uniroot(rootfunc, lower = -1e60, upper = 1e+60, u = u[i, j], f.lower=1-u[i,j], f.upper=-u[i,j],
#f.lower = 1, f.upper = -1e+60,
family = family, theta = theta, alpha = alphavec[j])
if(tmp$root > 1e58){#Beyond numerical precision
tmp$root <- Inf
} else if(tmp$root < -1e40){
tmp$root <- -Inf
}
out[i, j] <- exp(tmp$root)
}
}
out
}
##########################################################################
#' @title Joint survival function of Liouville vectors
#' @description \code{sliouv} returns the survival function of a Liouville vector. For the survival copula
#' and the associated probability, see \code{\link{pliouv}}.
#' @param x an \code{n} by \code{d} matrix of observations, each on the positive real line
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param alphavec \code{d}-vector of Dirichlet allocations (must be a vector of integers)
#' @param theta parameter of the corresponding Archimedean copula
#' @export
sliouv <- function(theta, x, family, alphavec){
alphavec <- as.integer(alphavec)
if(isTRUE(c(any(alphavec==0), any(x<0)))){stop("Invalid parameters")}
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
combo = unique(combn(unlist(sapply(alphavec, function(y){seq(1:y)})), length(alphavec)), MARGIN = 2)-1
surv <- function(x, alphasubvec){
if(length(alphasubvec)!=length(x)){stop("Invalid call to `surv' in function `sliouv'")}
if(sum(alphasubvec) != 0){
a <- alphasubvec*log(x)
if(isTRUE(any(alphasubvec==0, is.infinite(x)))){#Dealing with limit cases on the boundary,
#for which the formula ain't valid - setting the function to zero b/c the generator only matters in the limit
a[union(which(alphasubvec==0),which(is.infinite(x)))] <- 0
}
#Yields 0*Inf or similar problematic cases
exp(sum(a)-sum(lgamma(alphasubvec+1)) +
switch(family,
gumbel = copGumbel@absdPsi(sum(x), theta = theta, degree = sum(alphasubvec), log = TRUE),
clayton = copClayton@absdPsi(t = theta*sum(x), theta = theta, degree = sum(alphasubvec), log = TRUE) +
sum(alphasubvec)*log(theta),
frank = copFrank@absdPsi(sum(x), theta = theta, degree = sum(alphasubvec), log = TRUE),
AMH = copAMH@absdPsi(sum(x), theta = theta, degree = sum(alphasubvec), log = TRUE),
joe = copJoe@absdPsi(sum(x), theta = theta, degree = sum(alphasubvec), log = TRUE)))
} else {
switch(family,
gumbel = copGumbel@psi(sum(x), theta = theta),
clayton = copClayton@psi(t = theta*sum(x), theta = theta),
frank = copFrank@psi(sum(x), theta = theta),
AMH = copAMH@psi(sum(x), theta = theta),
joe = copJoe@psi(sum(x), theta = theta))
}
}
result = apply(combo, 2, function(alpha_combo){
if(is.matrix(x)){
apply(x,1, function(x){surv(x, alphasubvec=alpha_combo)})
} else{
surv(x)
}
}
)
if(is.matrix(x)){
return(rowSums(result))
} else{
return(sum(result))
}
}
##########################################################################
#' \bold{Liouville copula function}
#'
#' @details The Liouville copula is by definition a survival copula. The function thus
#' maps marginally observations from the unit interval to the positive half-line using the marginal inverse
#' survival function \code{\link{isliouvm}} of the Liouville vector, and then evaluating the survival
#' distribution at the resulting Liouville vector.
#' @rdname Liouville
#' @export
pliouv <- function(x, theta, family, alphavec){
alphavec <- as.integer(alphavec)
if(isTRUE(any(alphavec==0))){stop("Invalid parameters")}
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
sliouv(theta, isliouv_m(x, family, alphavec, theta), family, alphavec)
}
##########################################################################
## Liouville copula density function
#' @rdname Liouville
#' @export
dliouv <- function(x, family, alphavec, theta, is.log = FALSE){
alphavec <- as.integer(alphavec)
if(isTRUE(any(alphavec==0))){stop("Invalid parameters")}
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
if(!dim(x)[2] == length(alphavec))
{
stop("Invalid argument: alpha and x are not of same length")
}
#Check for validity of the returned parameter in light of family
alpha <- sum(alphavec)
x.norm <- rowSums(x)
psi.alpha <- switch(family,
gumbel = copGumbel@absdPsi(x.norm, theta = theta, degree = alpha, log = TRUE),
clayton = copClayton@absdPsi(x.norm*theta, theta = theta, degree = alpha, log = TRUE)+alpha*log(theta),
frank = copFrank@absdPsi(x.norm, theta = theta, degree = alpha, log = TRUE),
AMH = copAMH@absdPsi(x.norm, theta = theta, degree = alpha, log = TRUE),
joe = copJoe@absdPsi(x.norm, theta = theta, degree = alpha, log = TRUE))
arg1 <- t(apply(x, 1, function(x){(alphavec-1)*log(x)}))
return_val <- sum(arg1)+sum(psi.alpha)-dim(x)[1]*sum(lgamma(alphavec))
if(is.log){return(return_val)
} else{
return(exp(return_val))
}
}
##########################################################################
## Marginal density function of Liouville copulas
#' @rdname Liouville_marginal
#' @export
dliouvm <- function(x, family, alpha, theta){
if(ncol(as.matrix(x))>1) stop("x should be a vector or a n by 1 matrix")
alpha <- as.integer(alpha)
if(alpha==0){stop("Invalid parameter alpha")}
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
val <- switch(family,
gumbel = x^(alpha-1)/gamma(alpha)*copGumbel@absdPsi(x, theta = theta, degree = alpha),
clayton = x^(alpha-1)/gamma(alpha)*copClayton@absdPsi(t = theta*x, theta, degree = alpha)*(theta^alpha),
frank = x^(alpha-1)/gamma(alpha)*copFrank@absdPsi(x, theta = theta, degree = alpha),
AMH = x^(alpha-1)/gamma(alpha)*copAMH@absdPsi(x, theta = theta, degree = alpha),
joe = x^(alpha-1)/gamma(alpha)*copJoe@absdPsi(x, theta = theta, degree = alpha))
return(pmax(val, 0))
}
##########################################################################
#' Old code for the copula likelihood function for Liouville copulas
#'
#' The function is used internally for optimization.
#'
#' @param theta parameter of the corresponding Archimedean copula
#' @param data sample matrix from a Liouville copula
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#' @param MC.approx whether to use Monte-Carlo approximation for the inverse survival function (default is \code{TRUE})
#'
#' @return value of marginal density
.pliouv.opt_old <- function(theta, data, family, alphavec, MC.approx = TRUE){
alphavec <- as.integer(alphavec)
if(isTRUE(any(alphavec==0))){stop("Invalid parameters")}
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
#changed order of parameters for optimization functions
illegalpar <- switch(family,
clayton = copClayton@paraConstr(theta),
gumbel = copGumbel@paraConstr(theta),
frank = copFrank@paraConstr(theta),
AMH = copAMH@paraConstr(theta),
joe = copJoe@paraConstr(theta))
if(!illegalpar)
stop("Illegal parameter value")
n <- dim(data)[1]
d <- dim(data)[2]
Xdata = data
if(MC.approx){
Xdata <- H_inv(data, alphavec = alphavec, family = family, theta = theta)
}
Xmargin <- data
for (j in 1:d)
{
if(!MC.approx){
Xdata[, j] <- isliouvm(data[, j], family, alphavec[j], theta)
}
Xmargin[, j] <- log(dliouvm(Xdata[, j], family, alphavec[j], theta))
}
numerator.log <- dliouv(Xdata, family, alphavec, theta, is.log = TRUE)
denominator.log <- sum(Xmargin)
numerator.log - denominator.log
}
##########################################################################
#' Internal code for the copula likelihood function of Liouville copulas
#'
#' The function is used internally for optimization.
#'
#' @param theta parameter of the corresponding Archimedean copula
#' @param data sample matrix from a Liouville copula
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#' @param MC.approx whether to use Monte-Carlo approximation for the inverse survival function (default is \code{TRUE})
#'
#' @return value of marginal density
pliouv.opt <- function(theta, data, family, alphavec, MC.approx = TRUE){
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
#changed order of parameters for optimization functions
illegalpar <- switch(family,
clayton = copClayton@paraConstr(theta),
gumbel = copGumbel@paraConstr(theta),
frank = copFrank@paraConstr(theta),
AMH = copAMH@paraConstr(theta),
joe = copJoe@paraConstr(theta))
if(!illegalpar)
stop("Illegal parameter value")
n <- dim(data)[1]
d <- dim(data)[2]
Xdata = data
if(MC.approx == TRUE){
Xdata <- H_inv(data, alphavec = alphavec, family = family, theta = theta)
}
Xmargin <- data
for (j in 1:d){
if(MC.approx == FALSE){
Xdata[, j] <- isliouvm(data[, j], family, alphavec[j], theta)
}
Xmargin[, j] <- switch(family,
gumbel = copGumbel@absdPsi(Xdata[, j], theta = theta, degree = alphavec[j], log = TRUE),
clayton = copClayton@absdPsi(t = theta*Xdata[, j], theta, degree = alphavec[j], log = TRUE),
frank = copFrank@absdPsi(Xdata[, j], theta = theta, degree = alphavec[j], log = TRUE),
AMH = copAMH@absdPsi(Xdata[, j], theta = theta, degree = alphavec[j], log = TRUE),
joe = copJoe@absdPsi(Xdata[, j], theta = theta, degree = alphavec[j], log = TRUE)
)
}
alpha = sum(alphavec)
numerator.log <- sum(switch(family,
gumbel = copGumbel@absdPsi(rowSums(Xdata), theta = theta, degree = alpha, log = TRUE),
clayton = copClayton@absdPsi(t = theta*rowSums(Xdata), theta, degree = alpha, log = TRUE),
frank = copFrank@absdPsi(rowSums(Xdata), theta = theta, degree = alpha, log = TRUE),
AMH = copAMH@absdPsi(rowSums(Xdata), theta = theta, degree = alpha, log = TRUE),
joe = copJoe@absdPsi(rowSums(Xdata), theta = theta, degree = alpha, log = TRUE)
))
denominator.log <- sum(Xmargin)
numerator.log - denominator.log
}
##########################################################################
#' Maximization of Liouville copula likelihood function
#'
#' Two methods, either numerical optimization or method-of-moments
#'
#' A wrapper to \code{optim} using the Nelder-Mead algorithm or using the methods of moments,
#' to maxime pointwise given every \code{alphavec} over a grid.
#' Returns the maximum for \code{alphavec} and \code{theta}.
#' @aliases liouv.maxim liouv.maxim.mm
#' @param data sample matrix from a Liouville copula
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param interval interval over which to look for \code{theta} (bounds for Nelder-Mead)
#' @param boundary vector of endpoints for search of Dirichlet allocation parameters. Either \code{boundary} or \code{lattice.mat} can be supplied
#' @param lattice.mat matrix of tuples of Dirichlet allocation parameters at which to evaluate the likelihood
#' @param return_all should all results (as list) or only maximum value be returned. Defaults to \code{FALSE}
#' @param MC.approx whether to use Monte-Carlo approximation for the inverse survival function (default is \code{TRUE})
#'
#' @return a list with values of \code{theta} and Dirichlet parameter along with maximum found. Gives index of maximum amongst models fitted.
#' @export
#' @examples
#' \dontrun{
#' data <- rliouv(n=100, family="joe", alphavec=c(1,2), theta=2)
#' liouv.maxim(data=data, family="j", interval=c(1.25,3), boundary=c(2,2),return_all=TRUE)
#' lattice.mat <- t(combn(1:3,2))
#' liouv.maxim(data=data, family="j", interval=c(1.25,3), lattice.mat=lattice.mat, return_all=FALSE)
#' #data <- rliouv(n=1000, family="gumbel", alphavec=c(1,2), theta=2)
#' liouv.maxim.mm(data=data, family="gumbel", boundary=c(3,3),return_all=TRUE)
#' lattice.mat <- t(combn(1:3,2))
#' liouv.maxim.mm(data=data, family="gumbel", lattice.mat=lattice.mat, return_all=FALSE)
#' }
liouv.maxim <- function(data, family, interval, boundary = NULL, lattice.mat = NULL, return_all = FALSE, MC.approx = TRUE){
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
stopifnot(!is.null(lattice.mat) || !is.null(boundary))
if(!is.null(boundary)){stopifnot(length(boundary) == dim(data)[2]) }
if(!is.null(lattice.mat)){stopifnot(ncol(lattice.mat) == ncol(data)) }
stopifnot(length(interval) == 2, interval[1]<interval[2],
switch(family, clayton = copClayton@paraConstr(interval[1]),
gumbel = copGumbel@paraConstr(interval[1]),
frank = copFrank@paraConstr(interval[1]),
AMH = copAMH@paraConstr(interval[1]),
joe = copJoe@paraConstr(interval[1])),
switch(family, clayton = copClayton@paraConstr(interval[2]),
gumbel = copGumbel@paraConstr(interval[2]),
frank = copFrank@paraConstr(interval[2]),
AMH = copAMH@paraConstr(interval[2]),
joe = copJoe@paraConstr(interval[2])))
if(is.null(lattice.mat)){
term <- list()
for(i in 1:length(boundary)){
term[[i]] = 1:boundary[i]
}
lattice.mat <- matrix(unlist(expand.grid(term)),ncol=length(boundary))
}
#Define containers
like_array <- array(dim=apply(lattice.mat,2,max))
theta_array <- array(dim=apply(lattice.mat,2,max))
#Optimization over theta
for(i in 1:nrow(lattice.mat)){
lattice.vec <- lattice.mat[i,]
# Choice of maximization function, see http://cran.r-project.org/web/views/Optimization.html
res <- optimise(pliouv.opt, data = data, family = family, alphavec = lattice.vec,
interval = interval, maximum = TRUE, MC.approx = MC.approx)
theta_array[matrix(lattice.vec, 1)] <- res$maximum
like_array[matrix(lattice.vec, 1)] <- res$objective
}
indexed <- which.max(like_array)
pos <- which( like_array==max(like_array,na.rm=TRUE) , arr.ind = TRUE )
if(return_all){
return(list(index = pos,
theta_max = theta_array[indexed],
theta_val = theta_array,
log.lik_max = like_array[indexed],
log.lik = like_array))
} else{
return(list(index = pos,
theta_max = theta_array[indexed], log.lik = like_array[indexed]))
}
}
##########################################################################
## Maximize copula likelihood function for Liouville copulas via methods of moments
#' @export
liouv.maxim.mm <- function(data, family, boundary = NULL, lattice.mat = NULL, return_all = FALSE, MC.approx = TRUE){
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
stopifnot(!is.null(lattice.mat) || !is.null(boundary))
if(is.null(lattice.mat)){
stopifnot(length(boundary) == dim(data)[2] && dim(data)[2] == 2)
term <- list()
for(i in 1:length(boundary)){
term[[i]] = 1:as.integer(boundary[i])
}
lattice.mat <- matrix(unlist(expand.grid(term)),ncol=length(boundary))
}
if(!is.null(lattice.mat)){
stopifnot(ncol(lattice.mat) == ncol(data), ncol(data) == 2)
}
if (!requireNamespace("wdm", quietly = TRUE)) {
# warning("Install package \"wdm\" to speed up calculations.",
# call. = FALSE)
tau_hat <- cor(x = data[,1],
y = data[,2],
method = "kendall")
} else{
tau_hat <- wdm::wdm(x = data[,1],
y = data[,2],
method = "kendall")
}
#Create containers for values
like_array <- array(dim=apply(lattice.mat,2,max))
theta_array <- array(dim=apply(lattice.mat,2,max))
#Optimization over theta
for(i in 1:nrow(lattice.mat)){
lattice.vec <- lattice.mat[i,]
theta_hat<- liouv.iTau(tau_hat, family, alphavec = lattice.vec)
theta_array[matrix(lattice.vec, 1)] <- theta_hat
like_array[matrix(lattice.vec, 1)] <- pliouv.opt(theta_hat, data, family, lattice.vec, MC.approx = MC.approx)
}
#Which is the max
indexed <- which.max(like_array)
pos <- which( like_array==max(like_array,na.rm=TRUE) , arr.ind = TRUE )
if(return_all){
return(list(index = pos,
theta_max = theta_array[indexed],
theta_val = theta_array,
log.lik_max = like_array[indexed],
log.lik = like_array))
} else{
return(list(index = pos,
theta_max = theta_array[indexed], log.lik = like_array[indexed]))
}
}
##########################################################################
#' Computes Kendall's \eqn{\tau}{tau} for Clayton or Gumbel Liouville copula
#'
#' The function computes Kendall's \eqn{\tau}{tau} for the given model, given \code{alphavec}
#'
#' @param theta parameter of the corresponding Archimedean copula
#' @param family family of the Liouville copula. Either \code{"clayton"} or \code{"gumbel"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#'
#' @return value of \eqn{\tau}{tau}
.liouv.Tau_s <- function(theta, family, alphavec){
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
if(length(alphavec) != 2)
stop("Illegal parameter length - only implemented for the bivariate case")
a1 <- alphavec[1]; a2 <- alphavec[2]; a <- a1+a2
if(family == "clayton"){
express_coef_clayton <- function(i, j, theta, a){
return(exp(lgamma(1/theta+i+j)+lgamma(1/theta+a)+lgamma(2/theta)+lgamma(a+i+j)-
2*log(theta)-2*lgamma(1/theta+1)-lgamma(2/theta+i+j+a)-lgamma(a)))
}
}
if(family == "gumbel"){
express_coef_gumb <- function(index_k, index_l, theta, a){
if(index_l == 0){
d1 <- .coeffG(d = index_k, alpha = 1/theta, method = "direct", log = TRUE)
intern <- seq(1:index_k)
d2 <- lgamma(intern)-intern*log(2)
return(theta*sum(exp(d1+d2))/gamma(a))
} else{
d1 <- outer(.coeffG(d = index_k, alpha = 1/theta, method = "direct", log = TRUE),
.coeffG(d = index_l, alpha = 1/theta, method = "direct", log = TRUE), "+")
intern <- outer(seq(1:index_k), seq(1:index_l), "+")
d2 <- lgamma(intern)-intern*log(2)
return(theta*sum(exp(d1+d2))/gamma(a))
}
}
}
count <- sum(sapply(0:(a1-1), function(i){
sum(sapply(0:(a2-1), function(j){
exp(lbeta(a1+i, a2+j)-lbeta(a1, a2)-lgamma(i+1)-lgamma(j+1))*
switch(family,
gumbel = express_coef_gumb(index_k = a, index_l = i+j, theta, a = a),
clayton = express_coef_clayton(i = i, j = j, theta, a = a))
}))
}))
return(count*4-1)
}
#' Computes Kendall's tau for Clayton or Gumbel Liouville copula
#'
#' The function computes Kendall's \eqn{\tau}{tau} for the given model, given \code{alphavec}
#'
#' @param theta parameter of the corresponding Archimedean copula
#' @param family family of the Liouville copula. Either \code{"clayton"} or \code{"gumbel"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#' @examples
#' liouv.Tau(theta=2, family="gumbel", alphavec=c(1,2))
#' liouv.Tau(theta=1, family="clayton", alphavec=c(2,1))
#' @return vector of \eqn{\tau}{tau}
#' @export
liouv.Tau = Vectorize(FUN = .liouv.Tau_s, vectorize.args = "theta")
##########################################################################
#' Moment estimate for \code{theta} derived from Kendall's \eqn{\tau}{tau}
#'
#' The moment estimates are based on inversion of the formula Kendall's \eqn{\tau}{tau} for Clayton or Gumbel Liouville copula
#'
#' @param tau_hat estimated Kendall's \eqn{\tau}{tau} value from data
#' @param family family of the Liouville copula. Either \code{"clayton"} or \code{"gumbel"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#'
#' @return Value of theta
.liouv.iTau_s <- function(tau_hat, family, alphavec){
alphavec <- as.integer(alphavec)
if(isTRUE(any(alphavec==0))){stop("Invalid parameters")}
family <- match.arg(family, c("clayton", "gumbel"))
if(family != "gumbel" && family != "clayton")
stop("Not yet implemented for this family")
invtau <- function(theta, tau, family, alphavec)
{
liouv.Tau(theta, family, alphavec)-tau
}
lb <- switch(family, gumbel = 1+1e-60, clayton = 1e-60)
#Not implemented for negative theta values (Clayton family)
out <- uniroot(invtau, lower = lb, upper = 100, f.lower = -1, f.upper = 1,
alphavec = alphavec, tau = tau_hat, family = family)
out$root
}
#' Moment estimate for theta derived from Kendall's tau formula
#'
#' The moment estimates are based on inversion of the formula Kendall's tau for Clayton or Gumbel Liouville copula
#'
#' @param tau_hat estimated vector of Kendall's tau values
#' @param family family of the Liouville copula. Either \code{"clayton"} or \code{"gumbel"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#' @examples
#' liouv.iTau(0.5,family="gumbel", c(1,2))
#' liouv.iTau(0.5,family="clayton", c(3,2))
#' @return Vector of theta
#' @export
liouv.iTau = Vectorize(.liouv.iTau_s, "tau_hat")
##########################################################################
#' Inverse survival function if Monte-Carlo approximation is set to \code{TRUE} in \code{liouv.maxim}
#'
#' The function is used internally for optimization.
#'
#' @keywords internal
#'
#' @param u data at which to compute the survival inverse
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param theta parameter of the corresponding Archimedean copula
#' @param MC number of Monte-Carlo points for evaluation
#' @param trunc whether to truncate at low quantile. Will be based on numerical root finding for the lower 0.025 fraction of the data
#'
#' @return Inverse survival function values
#' @examples
#' \dontrun{
#' u <- rliouv(n = 100, family = "frank", alphavec <- c(2,3), theta = 1)
#' H_inv(u=u, family="frank", alphavec=c(2,3), theta=2)
#' #Difference between true value and approximation (can be large depending on family)
#' sum(abs(H_inv(u=u, family="frank", alphavec=c(2,3), theta=2)-
#' isliouv_m(u=u, family="frank", alphavec=c(2,3), theta=2)))
#' }
H_inv <- function(u, alphavec, family, theta, MC = 100000, TRUNC = FALSE){
alphavec <- as.integer(alphavec)
if(isTRUE(any(alphavec==0))){stop("Invalid parameters")}
family <- match.arg(family, c("clayton", "gumbel", "frank", "AMH", "joe"))
d <- sum(alphavec)
#Parametrization for Clayton is different
if(family == "clayton"){
mc <- (rCopula(n = MC, claytonCopula(param = theta, dim = d))^(-theta)-1)/theta
} else if(family=="AMH") {
mc <- iPsi(get(paste(tolower(family), "Copula", sep = ""))(param = theta),
rarchi(n = MC, family="AMH", theta = theta, d = d))
} else{
mc <- iPsi(get(paste(tolower(family), "Copula", sep = ""))(param = theta),
rCopula(n = MC, get(paste(tolower(family), "Copula", sep = ""))(param = theta, dim = d)))
}
mc <- .marginCombo(alphavec, mc)
inverse = numeric()
for(i in 1:length(alphavec)){
inverse = cbind(inverse, as.numeric(quantile(mc[, i], probs = (1-u[, i]))))
}
if(TRUNC){
hquant = which(u<0.025, arr.ind = TRUE)
for(i in 1:(dim(hquant)[1])){
inverse[hquant[i, 1], hquant[i, 2]] = isliouvm(u[hquant[i, 1], hquant[i, 2]], family, alphavec[hquant[i, 2]], theta)
}
hquant = which(u>0.975, arr.ind = TRUE)
for(i in 1:(dim(hquant)[1])){
inverse[hquant[i, 1], hquant[i, 2]] = isliouvm(u[hquant[i, 1], hquant[i, 2]], family, alphavec[hquant[i, 2]], theta)
}
}
return(inverse)
}
##########################################################################
### Rcpp function to sum margins - i.e. rowSums
#marginCombo() in Rcpp file
##########################################################################
#' Clayton coefficient for Kendall's tau formula
#'
#' The function is used internally in \cite{liouv.Tau}.
#'
#' @keywords internal
#'
#' @param i first index
#' @param j second index
#' @param theta parameter of the corresponding Archimedean copula
#' @param a \eqn{L_1}{L1}-norm of Dirichlet allocation vector
#'
#' @return Clayton coefficient
express_coef_clayton <- function(i, j, theta, a){
return(exp(lgamma(1/theta+i+j)+lgamma(1/theta+a)+lgamma(2/theta)+lgamma(a+i+j)-
2*log(theta)-2*lgamma(1/theta+1)-lgamma(2/theta+i+j+a)-lgamma(a)))
}
# Note: function coeffG is a copy of that found in the `copula' package
# the latter is unexported, and CRAN won't allow access via :::
# All copyrights and credits to copula authors
# (c) Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan
.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)
method <- match.arg(method)
if (method == "dsumSibuya") {
message("method 'dsumSibuya' is deprecated.\n ", "use method = 'dsSib.log' instead")
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)))
lS <- lapply(1:d, function(n) log(Stirling2.all(n)))
wrong.sign <- integer()
for (k in 1:d) {
j <- k:d
b <- j * log(alpha) + ls[j] + unlist(lapply(j, function(i) lS[[i]][k]))
b.max <- max(b)
exponents <- b - b.max
exps <- exp(exponents)
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 = {
k <- 1:d
ck <- 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)
pol <- vapply(k, function(k.) {
j <- 0:(d - 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)
k <- seq_len(d)
S <- lapply(k, Stirling2.all)
vapply(k, function(k.) {
j <- k.:d
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(gettextf("unsupported method '%s' in coeffG", method),
domain = NA))
}
##########################################################################
#' Gumbel coefficient for Kendall's tau formula
#'
#' The function is used internally in \cite{liouv.Tau}.
#'
#' @keywords internal
#'
#' @param index_k first index
#' @param index_l second index
#' @param theta parameter of the corresponding Archimedean copula
#' @param a sum of alphavec
#'
#' @return Gumbel coefficient
express_coef_gumb <- function(index_k, index_l, theta, a){
if(index_l == 0){
d1 <- .coeffG(d = index_k, alpha = 1/theta, method = "direct", log = TRUE)
intern <- seq_len(index_k)
d2 <- lgamma(intern)-intern*log(2)
return(theta*sum(exp(d1+d2))/gamma(a))
} else{
d1 <- outer(.coeffG(d = index_k, alpha = 1/theta, method = "direct", log = TRUE),
.coeffG(d = index_l, alpha = 1/theta, method = "direct", log = TRUE), "+")
intern <- outer(seq(1:index_k), seq(1:index_l), "+")
d2 <- lgamma(intern)-intern*log(2)
return(theta*sum(exp(d1+d2))/gamma(a))
}
}
##########################################################################
#' Parametric bootstrap confidence interval for the parameter \code{theta} for Liouville copula
#'
#' The parametric bootstrap provides confidence intervals by repeatedly sampling datasets from the postulated
#' Liouvilla copula model. If \eqn{d=2} and the model is either \code{gumbel} or \code{clayton}, the value of
#' Kendall's \eqn{\tau}{tau} is calculated from the sample, and the confidence interval or the quantiles correspond
#' to the inverse \eqn{\tau^{-1}(\tau(\theta))}{tau} for the bootstrap quantile values of \eqn{\tau}{tau} (using monotonicity).
#'
#' Install package \code{wdm} to speed up calculation of Kendall's tau.
#'
#' Since no closed-form formulas exist for the other models or in higher dimension,
#' the method is extremely slow since it relies on maximization
#' of a new sample from the model and look up the corresponding parameters.
#'
#' @param B number of bootstrap replicates
#' @param family family of the Liouville copula. Either \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"AMH"} or \code{"joe"}
#' @param alphavec vector of Dirichlet allocations (must be a vector of integers)
#' @param n sample size
#' @param theta.hat estimate of theta
#' @param quant if the vector of probability is specified, the function will return the corresponding bootstrap quantiles
#' @param silent boolean for output progress. Default is \code{FALSE}, which means iterations are printed if \eqn{d>2}.
#' @examples
#' \dontrun{
#' theta.bci(B=99, family="gumbel", alphavec=c(2,3), n=100, theta.hat=2)
#' theta.bci(B=19, family="AMH", alphavec=c(1,2), n=100, theta.hat=0.5, quant=c(0.05,0.95))
#' theta.bci(B=19, family="frank", alphavec=c(1,2,3), n=100, theta.hat=0.5, quant=c(0.05,0.95))
#' }
#' @return a list with a 95% confidence inteval unless selected quantiles \code{quant} are supplied
#' and the bootstrap values of Kendall's tau in \code{boot_tau} if \eqn{d=2} and the model is either \code{gumbel} or \code{clayton}.
#' Otherwise, the list contains \code{boot_theta}.
#' @export
theta.bci <- function(B = 1999,
family,
alphavec,
n,
theta.hat,
quant = c(0.025,0.975),
silent = FALSE){
family <- match.arg(arg = family,
several.ok = FALSE,
choices = c("clayton",
"gumbel",
"frank",
"AMH",
"joe"))
alphavec <- as.integer(alphavec)
if(isTRUE(any(alphavec==0))){
stop("Invalid parameters")
}
illegalpar <-
switch(family,
clayton = copClayton@paraConstr(theta.hat),
gumbel = copGumbel@paraConstr(theta.hat),
frank = copFrank@paraConstr(theta.hat),
AMH = copAMH@paraConstr(theta.hat),
joe = copJoe@paraConstr(theta.hat))
if(!illegalpar){
stop("Illegal parameter value")
}
d <- length(alphavec)
if(d == 2 && family %in% c("gumbel","clayton")){
boot.rep <- sapply(seq_len(B), function(x){
if (!requireNamespace("wdm", quietly = TRUE)) {
cor(x = rliouv(n = n,
family = family,
alphavec = alphavec,
theta = theta.hat),
method = "kendall")
} else{
wdm::wdm(x = rliouv(n = n,
family = family,
alphavec = alphavec,
theta = theta.hat),
method = "kendall")
} })
quantiles = NaN
if(is.vector(quant)){
quantiles = liouv.iTau(
quantile(boot.rep, probs = quant),
family, alphavec)
}
return(list(quantiles = quantiles, boot_kendall = boot.rep))
} else{
boot.rep <- sapply(1:B, function(x){
if(silent == FALSE){
print(paste("Step", x))
}
boot.samp <- rliouv(n = n, family, alphavec, theta.hat)
boot.opt <- liouv.maxim(boot.samp, family, lattice.mat = as.matrix(t(alphavec)), MC.approx = FALSE,
interval = switch(family, gumbel = c(1.001, 10),
clayton = c(0.001, 10),
AMH = c(0.001, 0.999),
joe = c(1.001, 10),
frank = c(0.001, 10)))
boot.opt$theta_max})
quantiles = NaN
if(is.vector(quant)){
quantiles = quantile(boot.rep, probs = quant)
}
return(list(quantiles = quantiles, boot_theta = boot.rep))
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.