rdstudent replaced with diststudent

R/rdstudent.R

-rdstudent <- function(nu1, Sigma1, nu2, Sigma2, bet, eps = 1e-06) {
-  #' Renyi Divergence between Centered Multivariate $t$ Distributions
-  #'
-  #' Computes the Renyi divergence between two random vectors distributed
-  #' according to multivariate $t$ distributions (MTD) with zero mean vector.
-  #'
-  #' @aliases rdstudent
-  #'
-  #' @usage rdstudent(nu1, Sigma1, nu2, Sigma2, bet, eps = 1e-06)
-  #' @param nu1 numéric. The degrees of freedom of the first distribution.
-  #' @param Sigma1 symmetric, positive-definite matrix. The correlation matrix of the first distribution.
-  #' @param nu2 numéric. The degrees of freedom of the second distribution.
-  #' @param Sigma2 symmetric, positive-definite matrix. The correlation matrix of the second distribution.
-  #' @param bet numeric., positive and not equal to 1. Order of the Renyi divergence.
-  #' @param eps numeric. Precision for the computation of the partial derivative of the Lauricella \eqn{D}-hypergeometric function (see Details). Default: 1e-06.
-  #' @return A numeric value: the Renyi divergence between the two distributions,
-  #' with two attributes \code{attr(, "epsilon")} (precision of the result of the Lauricella \eqn{D}-hypergeometric function,see Details)
-  #' and \code{attr(, "k")} (number of iterations).
-  #'
-  #' @details Given \eqn{X_1}, a random vector of \eqn{R^p} distributed according to the MTD
-  #' with parameters \eqn{(\nu_1, \mathbf{0}, \Sigma_1)}
-  #' and \eqn{X_2}, a random vector of \eqn{R^p} distributed according to the MTD
-  #' with parameters \eqn{(\nu_2, \mathbf{0}, \Sigma_2)}.
-  #'
-  #' Let \eqn{\delta_1 = \frac{\nu_1 + p}{2} \beta}, \eqn{\delta_2 = \frac{\nu_2 + p}{2} (1 - \beta)}
-  #' and \eqn{\lambda_1, \dots, \lambda_p} the eigenvalues of the square matrix \eqn{\Sigma_1 \Sigma_2^{-1}}
-  #' sorted in increasing order: \deqn{\lambda_1 < \dots < \lambda_{p-1} < \lambda_p}
-  #' The Renyi divergence between \eqn{X_1} and \eqn{X_2} is:
-  #' \deqn{D_R^\beta(\mathbd{X}_1||\mathbd{X}_1) = \frac{1}{\beta - 1} \bigg[ \beta \ln\left(\frac{\Gamma\left(\frac{\nu_1+p}{2}\right) \Gamma\left(\frac{\nu_2}{2}\right) \nu_2^{\frac{p}{2}}}{\Gamma\left(\frac{\nu_2+p}{2}\right) \Gamma\left(\frac{\nu_1}{2}\right) \nu_1^{\frac{p}{2}}}\right) + \ln\left(\frac{\Gamma\left(\frac{\nu_2+p}{2}\right)}{\Gamma\left(\frac{\nu_2}{2}\right)}\right) + \ln\left(\frac{\Gamma\left(\delta_1 + \delta_2 - \frac{p}{2})}{\Gamma(\delta_1 + \delta_2)}\right) - \frac{\beta}{2} \sum_{i=1}^p{\ln\lambda_i} + \ln F_D^{(p)}{\bigg( \delta_1, \frac{1}{2}, \dots, \frac{1}{2}; \delta_1+\delta_2; 1-\frac{\nu_2}{\nu_1 \lambda_1}, \dots, 1-\frac{\nu_2}{\nu_1 \lambda_p} \bigg)} \bigg]}
-  #'
-  #' where \eqn{F_D^{(p)}} is the Lauricella \eqn{D}-hypergeometric function defined for \eqn{p} variables:
-  #' \deqn{ \displaystyle{ F_D^{(p)}\left(a; b_1, ..., b_p; g; x_1, ..., x_p\right) = \sum\limits_{m_1 \geq 0} ... \sum\limits_{m_p \geq 0}{ \frac{ (a)_{m_1+...+m_p}(b_1)_{m_1} ... (b_p)_{m_p} }{ (g)_{m_1+...+m_p} } \frac{x_1^{m_1}}{m_1!} ... \frac{x_p^{m_p}}{m_p!} } } }
-  #' Its computation uses the \code{\link{lauricella}} function.
-  #'
-  #' @author Pierre Santagostini, Nizar Bouhlel
-  #' @references N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions, IEEE Signal Processing Letters.
-  #' \doi{10.1109/LSP.2023.3324594}
-  #'
-  #' @examples
-  #' \donttest{
-  #' nu1 <- 2
-  #' Sigma1 <- matrix(c(2, 1.2, 0.4, 1.2, 2, 0.6, 0.4, 0.6, 2), nrow = 3)
-  #' nu2 <- 4
-  #' Sigma2 <- matrix(c(1, 0.3, 0.1, 0.3, 1, 0.4, 0.1, 0.4, 1), nrow = 3)
-  #' rdstudent(nu1, Sigma1, nu2, Sigma2, bet = 0.5)
-  #' rdstudent(nu2, Sigma2, nu1, Sigma1, bet = 0.5)
-  #' }
-  #'
-  #' @export
-
-  # We must have: bet>1 and bet!=1
-  if (bet <= 0 | bet == 1)
-    stop("bet must be positive, different to 1.")
-
-  # Sigma1 and Sigma2 must be matrices
-  if (is.numeric(Sigma1) & !is.matrix(Sigma1))
-    Sigma1 <- matrix(Sigma1)
-  if (is.numeric(Sigma2) & !is.matrix(Sigma2))
-    Sigma2 <- matrix(Sigma2)
-
-  # Number of variables
-  p <- nrow(Sigma1)
-
-  # Sigma1 and Sigma2 must be square matrices with the same size
-  if (ncol(Sigma1) != p | nrow(Sigma2) != p | ncol(Sigma2) != p)
-    stop("Sigma1 et Sigma2 must be square matrices with rank p.")
-
-  # IS Sigma1 symmetric, positive-definite?
-  if (!isSymmetric(Sigma1))
-    stop("Sigma1 must be a symmetric, positive-definite matrix.")
-  lambda1 <- eigen(Sigma1, only.values = TRUE)$values
-  if (any(lambda1 < .Machine$double.eps))
-    stop("Sigma1 must be a symmetric, positive-definite matrix.")
-
-  # Is Sigma2 symmetric, positive-definite?
-  if (!isSymmetric(Sigma2))
-    stop("Sigma2 must be a symmetric, positive-definite matrix.")
-  lambda2 <- eigen(Sigma2, only.values = TRUE)$values
-  if (any(lambda2 < .Machine$double.eps))
-    stop("Sigma2 must be a symmetric, positive-definite matrix.")
-
-  # Eigenvalues of Sigma1 %*% inv(Sigma2)
-  lambda <- sort(eigen(Sigma1 %*% solve(Sigma2), only.values = TRUE)$values, decreasing = FALSE)
-  loglambda <- log(lambda)
-  nulambda <- lambda*nu1/nu2
-
-  delta1 <- bet*(nu1 + p)/2; delta2 <- (1 - bet)*(nu2 + p)/2
-
-  if (nulambda[1] > 1) {
-    lauric <- lauricella(delta1, rep(0.5, p), delta1 + delta2, 1 - 1/nulambda, eps = eps)
-  } else if (nulambda[p] < 1) {
-    lauric <- lauricella(delta2, rep(0.5, p), delta1 + delta2, 1 - nulambda, eps = eps)
-    lauric <- prod(sqrt(nulambda)) * lauric
-  } else {
-    lauric <- lauricella(delta2, c(rep(0.5, p-1), delta1 + delta2 - p/2), delta1 + delta2,
-                         c(1 - lambda[1:(p-1)]/lambda[p], 1/nulambda[p]), eps = eps)
-    lauric <- (nulambda[p])^(-delta2) * prod(sqrt(nulambda)) * lauric
-  }
-
-  gamma1p <- gamma((nu1 + p)/2); gamma2p <- gamma((nu2 + p)/2)
-  gamma1 <- gamma(nu1/2); gamma2 <- gamma(nu2/2)
-
-  result <- (
-    bet*log((gamma1p*gamma2*nu2^(p/2))/(gamma2p*gamma1*nu1^(p/2))) +
-    log(gamma2p/gamma2) + log(gamma(delta1 + delta2 - p/2)/gamma(delta1 + delta2)) -
-    bet/2 * sum(loglambda) + log(lauric)
-  )/(bet - 1)
-
-  return(result)
-}