R/MBCSec_class.R

In rdmatheus/mbcsec2: 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 df Degrees of freedom associated with 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 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}),
#'      degrees of freedom (\code{df}), 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", "ST", "DE")
#'
#' # 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, 2, NA))
#'
#' ### 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, df = 4,
                    copula = c("gaussian", "t", "cauchy", "dexponential"),
                    gen = "NO"){

  ### 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 <- param$mu
  sigma <- param$sigma
  lambda <- param$lambda
  nu <- param$nu

  if (is.null(nu)) nu <- rep(NA, d)

  ### Copula density argument
  if (length(gen) == 1){
    gen <- rep(gen, d)
    w <- pBCS(x, mu, sigma, lambda, nu, gen)
  }else{
    w <- matrix(mapply(pBCS, t(x), mu, sigma, lambda, nu, gen),
                byrow = TRUE, ncol = d)
    #w <- matrix(apply(matrix(gen, ncol = 1), 1,
    #       function(gen) pBCS(x, mu, sigma, lambda, nu, gen))[, 1], ncol = d)
  }

  if (copula != "t"){
    q <- BCSgen(ell(copula)$gen)$q(w)
    den <- log(ell(copula)$Md(q, P)) +
      apply(log(matrix(mapply(dBCS, t(x), mu, sigma, lambda, nu, gen),
            byrow = TRUE, ncol = d)), 1, sum) -
      apply(matrix(log(BCSgen(ell(copula)$gen)$d(q^2)), ncol = d), 1, sum)

      #mvtnorm::dmvnorm(q, sigma = P, log = TRUE) +
      #apply(matrix(log(dBCS(x, mu, sigma, lambda, nu, gen[1])), ncol = d), 1, sum) -
      #apply(matrix(log(stats::dnorm(q)), ncol = d), 1, sum)
  } else {
    q <- BCSgen(ell(copula)$gen)$q(w, df)
    den <- log(ell(copula)$Md(q, P, df)) +
      apply(log(matrix(mapply(dBCS, t(x), mu, sigma, lambda, nu, gen),
            byrow = TRUE, ncol = d)), 1, sum) -
      apply(matrix(log(BCSgen(ell(copula)$gen)$d(q^2, df)), ncol = d), 1, sum)

    #q <- stats::qt(w, df = df)
    #den <- mvtnorm::dmvt(q, delta = rep.int(0, d), sigma = P,
    #                     df = df, log = TRUE) +
    #apply(matrix(log(dBCS(x, mu, sigma, lambda, nu, gen[1])), ncol = d), 1, sum) -
    #apply(matrix(log(stats::dt(q, df = df)), ncol = d), 1, sum)

  }

  pmax(exp(den), .Machine$double.eps)

}

# Random generator
#' @rdname MBCSec
#' @export
rmbcsec <- function(n, param, P = NULL, d = 2L, df = 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 <- param$mu
  sigma <- param$sigma
  lambda <- param$lambda
  nu <- param$nu

  if (length(gen) == 1) gen <- rep(gen, d)
  if (is.null(nu)) nu <- rep(NA, d)


  if (copula != "t"){
    x <- BCSgen(ell(copula)$gen)$p(ell(copula)$Mr(n, P))
    y <- matrix(mapply(qBCS, t(x), mu, sigma, lambda, nu, gen),
                byrow = TRUE, ncol = d)

    #x <- mvtnorm::rmvnorm(n, sigma = P)
    #y <- qBCS(stats::pnorm(x), mu, sigma, lambda, nu, gen)
  } else {

    if (is.null(df)) df <- 4

    x <- BCSgen(ell(copula)$gen)$p(ell(copula)$Mr(n, P, df), df)
    y <- matrix(mapply(qBCS, t(x), mu, sigma, lambda, nu, gen),
                byrow = TRUE, ncol = d)

    #x <- mvtnorm::rmvt(n, sigma = P, df = df)
    #y <- qBCS(stats::pt(x, df), mu, sigma, lambda, nu, gen)
  }

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

  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$df)) cat(" (Df = ", x$spec$df, ")", 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)

}