Fixed: 2 errors in the computing. Order of the arguments. The example is no more 'donttest'.

R/kldstudent.R

-kldstudent <- function(Sigma1, nu1, Sigma2, nu2, eps = 1e-06) {
+kldstudent <- function(nu1, Sigma1, nu2, Sigma2, eps = 1e-06) {
   #' Kullback-Leibler Divergence between Centered Multivariate \eqn{t} Distributions
   #'
   #' Computes the Kullback-Leibler divergence between two random vectors distributed
@@ -6,11 +6,11 @@ kldstudent <- function(Sigma1, nu1, Sigma2, nu2, eps = 1e-06) {
   #'
   #' @aliases kldstudent
   #'
-  #' @usage kldstudent(Sigma1, nu1, Sigma2, nu2, eps = 1e-06)
-  #' @param Sigma1 symmetric, positive-definite matrix. The scatter matrix of the first distribution.
+  #' @usage kldstudent(nu1, Sigma1, nu2, Sigma2, eps = 1e-06)
   #' @param nu1 numeric. The degrees of freedom of the first distribution.
-  #' @param Sigma2 symmetric, positive-definite matrix. The scatter matrix of the second distribution.
+  #' @param Sigma1 symmetric, positive-definite matrix. The scatter matrix of the first distribution.
   #' @param nu2 numeric. The degrees of freedom of the second distribution.
+  #' @param Sigma2 symmetric, positive-definite matrix. The scatter matrix of the second distribution.
   #' @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 Kullback-Leibler divergence between the two distributions,
   #' with two attributes \code{attr(, "epsilon")} (precision of the partial derivative of the Lauricella \eqn{D}-hypergeometric function,see Details)
@@ -38,19 +38,13 @@ kldstudent <- function(Sigma1, nu1, Sigma2, nu2, eps = 1e-06) {
   #' \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)
   #'
-  #' Sigma1 <- matrix(c(0.5, 0, 0, 0, 0.4, 0, 0, 0, 0.3), nrow = 3)
-  #' Sigma2 <- diag(1, 3)
-  #' # Case when all eigenvalues of Sigma1 %*% solve(Sigma2) are < 1
-  #' kldcauchy(Sigma1, Sigma2)
-  #' # Case when all eigenvalues of Sigma1 %*% solve(Sigma2) are > 1
-  #' kldcauchy(Sigma2, Sigma1)
-  #' }
+  #' kldstudent(nu1, Sigma1, nu2, Sigma2)
+  #' kldstudent(nu2, Sigma2, nu1, Sigma1)
   #'
   #' @importFrom utils combn
   #' @export
@@ -284,7 +278,7 @@ kldstudent <- function(Sigma1, nu1, Sigma2, nu2, eps = 1e-06) {
 
       }
     }
-    derive <- prod(nulambda)^(-0.5) * derive
+    derive <- prod(1/sqrt(lambdanu)) * derive
 
   } else { # lambda[1] < ... < 1 < ... < lambda[p]
 
@@ -394,12 +388,12 @@ kldstudent <- function(Sigma1, nu1, Sigma2, nu2, eps = 1e-06) {
 
       }
     }
-    derive <- -log(lambdanu[p])*derive
+    derive <- -log(lambdanu[p]) + derive
 
   }
 
   result <- log(gamma((nu1+p)/2)*gamma(nu2/2)*nu2^(p/2) / (gamma((nu2+p)/2)*gamma(nu1/2)*nu1^(p/2)))
-  result <- result + (nu1-nu2)/2 * (digamma((nu1+1)/2) - digamma(nu1/2) - 0.5 * sum(log(lambda)))
+  result <- result + (nu2-nu1)/2 * (digamma((nu1+p)/2) - digamma(nu1/2)) - 0.5 * sum(log(lambda))
   result <- result - (nu2 + p)/2 * derive
 
   if (is.nan(d)) {