R/MBCSec_class.R

In rdmatheus/mbcsec: Multivariate Box-Cox Symmetric Distributions Generated From Elliptical Copula

#' @name MBCSec
#'
#' @title Multivariate Box-Cox Symmetric Class of Distributions Generating
#'     by Elliptical Copula
#'
#' @description These functions provide the joint probability density function
#'     and a random generator for the multivariate Box-Cox symmetric class of
#'     distributions generated by elliptical copula.
#'
#' @param x Vector or matrix of non-negative quantiles. If \code{x} is a matrix,
#'     each row is taken to be a quantile. For methods, \code{x} is
#'     a \code{"mbcsec.data"} object.
#' @param n Number of random values to return.
#' @param d Dimension of the data to be generated. Argument used only in
#'     \code{rmbcsec()}.
#' @param param A list that has as components the vectors of marginal
#'     parameters \code{mu}, \code{sigma}, \code{lambda}, and \code{nu}. In case
#'     the marginal distribution is not indexed by \code{nu}, it receives
#'     \code{NULL} by default.
#' @param P Association matrix. If it is specified with \code{NULL} (default),
#'     then the identity matrix is used.
#' @param delta Extra parameters associated with the copula. For example,
#'     the degrees of freedom of the t copula.
#' @param copula Character; informs which distribution in the class of
#'     elliptical distributions should be used to generate the elliptical
#'     copula. Currently, the copulas available are:
#'     Gaussian (\code{"gaussian"}), Student's t (\code{"t"}), Cauchy
#'     (\code{"cauchy"}), and double exponential (\code{"dexponential"}).
#' @param gen A character or a vector character; specifies the distributions
#'     generating the marginal BCS distributions. If all BCS margins are
#'     generated by the same generating density, it is sufficient to enter only
#'     one character. A table with the current available generating
#'     distributions for the BCS class can be seen in \code{\link{BCSgen}}.
#' @param log.p logical; if \code{TRUE}, probabilities \code{p} are given
#'     as \code{log(p)}.
#' @param method A character string passed to the plot method indicating which
#'     correlation coefficient is to be computed. One of \code{"pearson"}
#'     (default), \code{"kendall"}, or \code{"spearman"}: can be abbreviated.
#' @param ... Further arguments for other specific methods.
#'
#' @return \code{dmbcsec} returns the evaluated joint density function and
#'  \code{rmbcsec} returns a \code{"mbcsec.data"} object that contains the
#'      observations generated and the specifications for the multivariate
#'      distribution. More specifically, it returns a list with the
#'      following components:
#'
#'  \itemize{
#'    \item{data:}{ Random observations generated from the multivariate
#'      distribution specified via \code{rmbcsec()}.}
#'    \item{spec:}{ A list with the distribution specifications:
#'      marginal parameters (\code{param}), association matrix (\code{P}),
#'      extra parameters of the copula (\code{delta}), copula (\code{copula}), and
#'      generating functions for margins (\code{gen}).}
#'  }
#'
#' @references
#'  Ferrari, S. L., & Fumes, G. (2017). Box-Cox symmetric distributions and
#'      applications to nutritional data. AStA Advances in Statistical Analysis,
#'      101, 321--344.
#'
#'  Vanegas, L. H., & Paula, G. A. (2016). Log-symmetric distributions:
#'      statistical properties and parameter estimation. Brazilian Journal of
#'      Probability and Statistics, 30, 196--220.
#'
#' @examples
#' \dontrun{
#' ### Sample size and dimension
#' n <- 500
#' d <- 3
#'
#' ### Association matrix
#' P <- matrix(c(1, 0.7, -0.5, 0.7, 1, 0.2, -0.5, 0.2, 1), 3, 3)
#'
#' ### Marginals specifications
#'
#' # Marginals
#' gen <- c("NO", "DE", "ST")
#'
#' # Marginal parameters
#' param <- list(mu = c(10, 6, 1), sigma = c(0.15, 0.15, 0.1),
#'               lambda = c(-2, 0, 2), nu = c(NA, NA, 2))
#'
#' ### Copula
#' copula <- "gaussian"
#'
#' ### Generating observations
#' y <- rmbcsec(n, param, P, d, copula = copula, gen = gen)
#'
#' y
#' plot(y)
#' }
#'
#' @author Rodrigo M. R. Medeiros <\email{rodrigo.matheus@live.com}>
#'

# Joint density
#' @rdname MBCSec
#' @export
dmbcsec <- function(x, param, P = NULL, delta = NULL,
                    copula = c("gaussian", "t", "cauchy", "dexponential"),
                    gen = "NO", log.p = FALSE){

  ### Reading the copula generating density
  copula <- match.arg(copula, c("gaussian", "t", "cauchy", "dexponential"))

  ### Setting dimensions
  if (is.vector(x))  x <- matrix(x, ncol = length(x))

  n <- dim(x)[1]
  d <- dim(x)[2]

  ### Default association matrix
  if (is.null(P)) P <- diag(d)

  ### Marginal parameters
  mu <- matrix(param$mu, ncol = d)
  sigma <- matrix(param$sigma, ncol = d)
  lambda <- matrix(param$lambda, ncol = d)
  nu <- param$nu

  mu <- do.call(rbind, replicate(n/dim(mu)[1], mu, simplify=FALSE))
  sigma <- do.call(rbind, replicate(n/dim(sigma)[1], sigma, simplify=FALSE))
  lambda <- do.call(rbind, replicate(n/dim(lambda)[1], lambda, simplify=FALSE))

  if (!is.null(nu)){
    nu <- matrix(nu, n, d, byrow = TRUE)
  } else{
    nu <- matrix(NA, n, d)
  }

  gen.out <- gen

  if (length(gen) == 1){
    gen <- matrix(gen, n, d); gen.out <- gen[1]
  } else{
    gen <- matrix(gen, n, d, byrow = TRUE)
  }

  if ((ell(copula)$nextraP > 0 & is.null(delta)))
    delta <- rep(4, ell(copula)$nextraP)

  ### Copula density argument
  #w <- matrix(mapply(pBCS, x, mu, sigma, lambda, nu, gen), ncol = d)

  ### Log-likelihood
  #copulaGEN <- ell(copula)$gen
  #epsilon <- BCSgen(copulaGEN)$q(w, delta)

  #den <- log(ell(copula)$Md(epsilon, P, delta)) +
  #  apply(log(matrix(mapply(dBCS, x, mu, sigma, lambda, nu, gen), ncol = d)), 1, sum) -
  #  apply(log(BCSgen(copulaGEN)$d(epsilon^2, delta)), 1, sum)

  copulaGEN <- ell(copula)$gen
  w <- epsilon <- aux <- matrix(NA, n, d)
  for(i in 1:d){
    w[, i] <- do.call(pBCS, args = list(x[, i], mu[, i], sigma[, i], lambda[, i], nu[, i], gen[, i]))
    epsilon[, i] <- BCSgen(copulaGEN)$q(w[, i], delta)
    aux[, i] <- log(do.call(dBCS, args = list(x[, i], mu[, i], sigma[, i], lambda[, i], nu[, i], gen[, i])))
  }

  den <- log(ell(copula)$Md(epsilon, P, delta)) +
         apply(aux, 1, sum) -
         apply(log(BCSgen(copulaGEN)$d(epsilon^2, delta)), 1, sum)


  if (log.p){
    den
  } else{
    pmax(exp(den), .Machine$double.eps)
  }

}

# Random generator
#' @rdname MBCSec
#' @export
rmbcsec <- function(n, param, P = NULL, d = 2L, delta = NULL,
                    copula = c("gaussian", "t", "cauchy", "dexponential"),
                    gen = "NO"){

  ### Reading the copula generating density
  copula <- match.arg(copula, c("gaussian", "t", "cauchy", "dexponential"))

  ### Default association matrix
  if (is.null(P)) P <- diag(d)

  ### Marginal parameters
  mu <- matrix(param$mu, ncol = d)
  sigma <- matrix(param$sigma, ncol = d)
  lambda <- matrix(param$lambda, ncol = d)
  nu <- param$nu

  mu <- do.call(rbind, replicate(n/dim(mu)[1], mu, simplify=FALSE))
  sigma <- do.call(rbind, replicate(n/dim(sigma)[1], sigma, simplify=FALSE))
  lambda <- do.call(rbind, replicate(n/dim(lambda)[1], lambda, simplify=FALSE))

  if (!is.null(nu)){
    nu <- matrix(nu, n, d, byrow = TRUE)
  } else{
    nu <- matrix(NA, n, d)
  }

  gen.out <- gen

  if (length(gen) == 1){
    gen <- matrix(gen, n, d); gen.out <- gen[1]
  } else{
    gen <- matrix(gen, n, d, byrow = TRUE)
  }

  if ((ell(copula)$nextraP > 0 & is.null(delta)))
    delta <- rep(4, ell(copula)$nextraP)

  ## Sampling for elliptical copula
  u <- BCSgen(ell(copula)$gen)$p(ell(copula)$Mr(n, d, P, delta), delta)

  ## Sampling for the multivariate distribution
  y <- matrix(NA, n, d)
  for(i in 1:d){
    y[, i] <- do.call(qBCS, args = list(u[, i], mu[, i], sigma[, i],
                                        lambda[, i], nu[, i], gen[, i]))
  }
  #y <- matrix(mapply(qBCS, u, mu, sigma, lambda, nu, gen), ncol = d)

  if (n == 1) y <- as.numeric(y)

  out <- list(data = y,
              spec = list(param = param,
                          P = P,
                          delta = delta,
                          copula = copula,
                          gen = gen.out))

  class(out) <- "mbcsec.data"
  out

}

#' @rdname MBCSec
#' @export
print.mbcsec.data <- function(x, ...){

  n <- dim(x$data)[1]
  d <- dim(x$data)[2]
  colnames(x$data) <- paste0("y", 1:d)
  row.names(x$data) <- 1:n
  colnames(x$spec$P) <- row.names(x$spec$P) <- paste0("y", 1:d)

  cat("------------",
      "\n   MBCSec",
      "\n------------",
      "\nCopula:", ell(x$spec$copula)$name)
  if (!is.null(x$spec$delta)) cat(" (delta = ", x$spec$delta, ")", sep = "")
  cat("\nMarginals:", paste(x$spec$gen, sep = ","),
      "\nAssociation matrix:\n\n")
  print(x$spec$P)
  cat("\n\nFirst observations:\n")

  print(utils::head(x$data))

  cat("---")
}

#' @rdname MBCSec
#' @export
plot.mbcsec.data <- function(x, method = c("pearson", "kendall", "spearman"),
                             ...){

  y <- x$data
  d <- dim(y)[2]
  param <- x$spec$param
  P <- x$spec$P
  copula <- x$spec$copula
  gen <- x$spec$gen

  if (length(gen) == 1) gen <- rep(gen, d)

  op <- graphics::par(mfrow = c(d, d), mar = c(2, 2, 1, 1) + 0.1)
  for(i in 1:d){
    for(j in 1:d){

      if (j <= i){
        if (i == j){
          graphics::hist(y[, i], prob = TRUE, main = " ", xlab = quote(y[i]))
          graphics::curve(dBCS(x, param$mu[i], param$sigma[i],
                               param$lambda[i], param$nu[i], gen = gen[i]),
                          add = TRUE, col = 2, lwd = 2)
        } else {
          Paux <- matrix(c(1, P[j, i], P[j, i], 1), 2, 2)
          param_aux <- list(mu = param$mu[c(j, i)],
                            sigma = param$sigma[c(j, i)],
                            lambda = param$lambda[c(j, i)],
                            nu = param$nu[c(j, i)])

          faux <- function(x, y){
            dmbcsec(cbind(x, y), param_aux, Paux,
                    copula = copula, gen = gen[c(j, i)])
          }

          z <- outer(seq(min(y[, j]), max(y[, j]), length.out = 30),
                     seq(min(y[, i]), max(y[, i]), length.out = 30),
                     faux)
          graphics::smoothScatter(y[, c(j, i)],
                                  xlab = quote(y[j]), ylab = quote(y[i]))
          graphics::contour(seq(min(y[, j]), max(y[, j]), length.out = 30),
                            seq(min(y[, i]), max(y[, i]), length.out = 30),
                            z, add = TRUE)
        }
      } else {
        plot(1:10, 1:10, type = "n", axes = FALSE)
        graphics::text(5, 5,
                       paste("Corr:\n", round(stats::cor(y[, j], y[, i]), 2)), cex = 1.2)
      }


    }
  }
  graphics::par(op)

}