kldcauchy.R

kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
  #' Kullback-Leibler Divergence between Centered Multivariate Cauchy Distributions
  #'
  #' Computes the Kullback-Leibler divergence between two random vectors distributed
  #' according to multivariate Cauchy distributions (MCD) with zero location vector.
  #'
  #' @aliases kldcauchy
  #'
  #' @usage kldcauchy(Sigma1, Sigma2, eps = 1e-06)
  #' @param Sigma1 symmetric, positive-definite matrix. The scatter matrix of the first 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)
  #' and \code{attr(, "k")} (number of iterations).
  #'
  #' @details Given \eqn{X_1}, a random vector of \eqn{R^p} distributed according to the MCD
  #' with parameters \eqn{(0, \Sigma_1)}
  #' and \eqn{X_2}, a random vector of \eqn{R^p} distributed according to the MCD
  #' with parameters \eqn{(0, \Sigma_2)}.
  #' 
  #' Let \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}
  #' Depending on the values of these eigenvalues,
  #' the computation of the Kullback-Leibler divergence of \eqn{X_1} from \eqn{X_2}
  #' is given by:
  #' \itemize{
  #' \item if \eqn{\lambda_1 < 1} and \eqn{\lambda_p > 1}:\cr
  #' \eqn{ \displaystyle{ KL(X_1||X_2) = -\frac{1}{2} \ln{ \prod_{i=1}^p{\lambda_i}} + \frac{1+p}{2} \bigg( \ln{\lambda_p} } } \cr
  #' \eqn{ \displaystyle{ - \frac{\partial}{\partial a} \bigg\{ F_D^{(p)} \bigg( a, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}, a + \frac{1}{2}}_p ; a + \frac{1+p}{2} ; 1 - \frac{\lambda_1}{\lambda_p}, \dots, 1 - \frac{\lambda_{p-1}}{\lambda_p}, 1 - \frac{1}{\lambda_p} \bigg) \bigg\}\bigg|_{a=0} \bigg) } }
  #' 
  #' \item if \eqn{\lambda_p < 1}:\cr
  #' \eqn{ \displaystyle{ KL(X_1||X_2) = -\frac{1}{2} \ln{ \prod_{i=1}^p{\lambda_i}} } }
  #' \eqn{ \displaystyle{ - \frac{1+p}{2} \frac{\partial}{\partial a} \bigg\{ F_D^{(p)} \bigg( a, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p ; a + \frac{1+p}{2} ; 1 - \lambda_1, \dots, 1 - \lambda_p \bigg) \bigg\}\bigg|_{a=0} } }
  #' 
  #' \item if \eqn{\lambda_1 > 1}:\cr
  #' \eqn{ \displaystyle{ KL(X_1||X_2) = -\frac{1}{2} \ln{ \prod_{i=1}^p{\lambda_i}} - \frac{1+p}{2} \prod_{i=1}^p\frac{1}{\sqrt{\lambda_i}} } } \cr
  #' \eqn{ \displaystyle{ \times \frac{\partial}{\partial a} \bigg\{ F_D^{(p)} \bigg( \frac{1+p}{2}, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p ; a + \frac{1+p}{2} ; 1 - \frac{1}{\lambda_1}, \dots, 1 - \frac{1}{\lambda_p} \bigg) \bigg\}\bigg|_{a=0} } }
  #' }
  #' 
  #' 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!} } } }
  #' 
  #' The computation of the partial derivative uses the \code{\link{pochhammer}} function.
  #' 
  #' @author Pierre Santagostini, Nizar Bouhlel
  #' @references N. Bouhlel, D. Rousseau, A Generic Formula and Some Special Cases for the Kullback–Leibler Divergence between Central Multivariate Cauchy Distributions.
  #' Entropy, 24, 838, July 2022.
  #' \doi{10.3390/e24060838}
  #'
  #' @examples
  #' \donttest{
  #' Sigma1 <- matrix(c(1, 0.6, 0.2, 0.6, 1, 0.3, 0.2, 0.3, 1), nrow = 3)
  #' Sigma2 <- matrix(c(1, 0.3, 0.1, 0.3, 1, 0.4, 0.1, 0.4, 1), nrow = 3)
  #' kldcauchy(Sigma1, Sigma2)
  #' kldcauchy(Sigma2, Sigma1)
  #' 
  #' 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)
  #' }
  #' 
  #' @importFrom utils combn
  #' @export
  
  # 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)
  
  prodlambda <- prod(lambda)
  
  k <- 5
  
  # M: data.frame of the indices for the nested sums
  # (i.e. all arrangements of n elements from {0:k})
  M <- expand.grid(rep(list(0:k), p))
  M <- M[-1, , drop = FALSE]
  
  # Sum of the indices
  Msum <- rowSums(M)
  
  Munique <- 0:max(M)
  Msumunique <- 0:max(Msum)
  
  kstep <- 5
  
  if (lambda[p] <= 1) { # equations 109 & 102
    
    if (p > 1 & lambda[p] == 1) {
      lambda <- lambda[1:(p-1)]
      p1 <- p-1
      M <- expand.grid(rep(list(0:k), p1))
      M <- M[-1, , drop = FALSE]
      # Sum of the indices
      Msum <- rowSums(M)
    } else {
      p1 <- p
    }
    
    # The first 5 elements of the sum
    
    # d <- 0
    # for (i in 1:length(Msum)) {
    #   commun <- prod(
    #     sapply(1:p, function(j) {
    #       pochhammer(0.5, M[i, j])*(1 - lambda[j])^M[i, j]/factorial(M[i, j])
    #     })
    #   )
    #   d <- d + commun * pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
    # }
    
    Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
    matpoch <- matrix(rep(Mpoch, p1), ncol = p1, byrow = FALSE)
    matM <- matrix(rep(Munique, p1), ncol = p1, byrow = FALSE)
    matlabmda <- matrix(rep(log(1 - lambda), length(Munique)), ncol = p1, byrow = TRUE)
    lambdaM <- matlabmda*matM
    matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
    
    matcommun  <- exp(matpoch + lambdaM - matfact)
    matcommun <- as.data.frame(apply(matcommun, 2, Re))
    matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
    colnames(matcommunsupp) <- colnames(matcommunsupp)
    
    gridcommun <- expand.grid(matcommun)
    gridcommun <- gridcommun[-1, , drop = FALSE]
    commun <- apply(gridcommun, 1, prod)
    
    d <- sum(
      commun * sapply(Msum, function(i) {
        pochhammer(1, i) / ( pochhammer((1 + p)/2, i) * i )
      })
    )
    
    # Next elements of the sum, until the expected precision
    k1 <- 1:k
    derive <- 0
    while (abs(d) > eps/10 & !is.nan(d)) {
      
      # cat(tail(k1, 1), "\t"); print(d)
      
      epsret <- signif(abs(d), 1)*10
      k <- k1[length(k1)]
      k1 <- k + (1:kstep)
      derive <- derive + d
      
      # # M: data.frame of the indices for the nested sums
      # M <- as.data.frame(matrix(nrow = 0, ncol = p))
      # if (p > 1) {
      #   for (i in 1:(p-1)) {
      #     Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
      #     M <- rbind( M, expand.grid(Mlist) )
      #     for (j in 1:(p-1)) {
      #       Mlist <- Mlist[c(p, 1:(p-1))]
      #       M <- rbind(M, expand.grid(Mlist))
      #     }
      #   }
      # }
      # M <- rbind( M, expand.grid(rep(list(k1), p)) )
      
      # M: data.frame of the indices for the nested sums
      M <- expand.grid(rep(list(k1), p1))
      if (p1 > 1) {
        for (i in 1:(p1-1)) {
          indsupp <- combn(p1, i)
          for (j in 1:ncol(indsupp)) {
            jsupp <- indsupp[, j]
            Mlist <- vector("list", p1)
            for (l in jsupp) Mlist[[l]] <- k1
            for (l in (1:p1)[-jsupp]) Mlist[[l]] <- 0:k
            M <- rbind(M, expand.grid(Mlist))
          }
        }
      }
      
      # Sum of the indices
      Msum <- rowSums(M)
      
      Munique <- (max(Munique)+1):max(M)
      Msumunique <- (max(Msumunique)+1):max(Msum)
      
      # d <- 0
      # for (i in 1:length(Msum)) {
      #   commun <- prod(
      #     sapply(1:p, function(j) {
      #       pochhammer(0.5, M[i, j])*(1 - lambda[j])^M[i, j]/factorial(M[i, j])
      #     })
      #   )
      #   d <- d + commun * pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
      # }
      
      Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
      matpoch <- matrix(rep(Mpoch, p1), ncol = p1, byrow = FALSE)
      matM <- matrix(rep(Munique, p1), ncol = p1, byrow = FALSE)
      matlabmda <- matrix(rep(log(1 - lambda), length(Munique)), ncol = p1, byrow = TRUE)
      lambdaM <- matlabmda*matM
      matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
      
      matcommun <- rbind(matcommun, matcommunsupp)
      matcommunsupp  <- exp(matpoch + lambdaM - matfact)
      matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
      
      gridcommun <- expand.grid(matcommunsupp)
      if (p1 > 1) {
        for (i in 1:(p1-1)) {
          indsupp <- combn(p1, i)
          for (j in 1:ncol(indsupp)) {
            jsupp <- indsupp[, j]
            communlist <- vector("list", p1)
            names(communlist) <- names(gridcommun)
            for (l in jsupp) communlist[[l]] <- matcommunsupp[, l]
            for (l in (1:p1)[-jsupp]) communlist[[l]] <- matcommun[, l]
            gridcommun <- rbind(gridcommun, expand.grid(communlist))
          }
        }
      }
      commun <- apply(gridcommun, 1, prod)
      
      d <- sum(
        commun * sapply(Msum, function(i) {
          pochhammer(1, i) / ( pochhammer((1 + p)/2, i) * i )
        })
      )
      
    }
    
    # Next elements of the sum, with step=1, while not NaN
    # if (is.nan(d)) {
    #   k1 <- k
    #   d <- 0
    #   while (!is.nan(d)) {
    #     if (d > 0)
    #       epsret <- signif(abs(d), 1)*10
    #     k <- k1
    #     k1 <- k + 1
    #     derive <- derive + d
    # 
    #     # # M: data.frame of the indices for the nested sums
    #     # M <- as.data.frame(matrix(nrow = 0, ncol = p))
    #     # if (p > 1) {
    #     #   for (i in 1:(p-1)) {
    #     #     Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
    #     #     M <- rbind( M, expand.grid(Mlist) )
    #     #     for (j in 1:(p-1)) {
    #     #       Mlist <- Mlist[c(p, 1:(p-1))]
    #     #       M <- rbind(M, expand.grid(Mlist))
    #     #     }
    #     #   }
    #     # }
    #     # M <- rbind( M, rep(k1, p) )
    #     
    #     # M: data.frame of the indices for the nested sums
    #     M <- expand.grid(rep(list(k1), p))
    #     if (p > 1) {
    #       for (i in 1:(p-1)) {
    #         indsupp <- combn(p, i)
    #         for (j in 1:ncol(indsupp)) {
    #           jsupp <- indsupp[, j]
    #           Mlist <- vector("list", p)
    #           for (l in jsupp) Mlist[[l]] <- k1
    #           for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
    #           M <- rbind(M, expand.grid(Mlist))
    #         }
    #       }
    #     }
    #     
    #     # Sum of the indices
    #     Msum <- rowSums(M)
    #     
    #     Munique <- (max(Munique)+1):max(M)
    #     Msumunique <- (max(Msumunique)+1):max(Msum)
    #     
    #     # d <- 0
    #     # for (i in 1:length(Msum)) {
    #     #   commun <- prod(
    #     #     sapply(1:p, function(j) {
    #     #       pochhammer(0.5, M[i, j])*(1 - lambda[j])^M[i, j]/factorial(M[i, j])
    #     #     })
    #     #   )
    #     #   d <- d + commun * pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
    #     # }
    #     
    #     Mpoch <- sapply(Munique, function(j) pochhammer(0.5, j))
    #     matpoch <- matrix(rep(Mpoch, p), ncol = p, byrow = FALSE)
    #     matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
    #     matlabmda <- matrix(rep(1 - lambda, length(Munique)), ncol = p, byrow = TRUE)
    #     lambdaM <- matlabmda^matM
    #     matfact <- matrix(rep(factorial(Munique), p), ncol = p, byrow = FALSE)
    #     
    #     matcommun  <- as.data.frame(matpoch*lambdaM/matfact)
    #     
    #     gridcommun <- expand.grid(matcommun)
    #     gridcommun <- gridcommun[-1, , drop = FALSE]
    #     commun <- apply(gridcommun, 1, prod)
    #     
    #     d <- sum(
    #       commun * sapply(Msum, function(i) {
    #         pochhammer(1, i) / ( pochhammer((1 + p)/2, i) * i )
    #       })
    #     )
    # 
    #   }
    # }
    
    result <- as.numeric(-0.5 * log(prodlambda) - (1 + p)/2 * derive)
    
  } else if (lambda[1] > 1) { # equations 110
    
    # The first 5 elements of the sum
    
    # d <- 0
    # for (i in 1:length(Msum)) {
    #   commun <- prod(
    #     sapply(1:p, function(j) {
    #       pochhammer(0.5, M[i, j])*(1 - 1/lambda[j])^M[i, j]/factorial(M[i, j])
    #     })
    #   )
    #   A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
    #   d <- d - commun * A # / pochhammer((1 + p)/2, Msum[i])
    # }
    
    Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
    matpoch <- matrix(rep(Mpoch, p), ncol = p, byrow = FALSE)
    matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
    matlabmda <- matrix(rep(log(1 - 1/lambda), length(Munique)), ncol = p, byrow = TRUE)
    lambdaM <- matlabmda*matM
    matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
    
    matcommun  <- exp(matpoch + lambdaM - matfact)
    matcommun <- as.data.frame(apply(matcommun, 2, Re))
    matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
    colnames(matcommunsupp) <- colnames(matcommunsupp)
    
    gridcommun <- expand.grid(matcommun)
    gridcommun <- gridcommun[-1, , drop = FALSE]
    commun <- apply(gridcommun, 1, prod)
    
    d <- 0
    for (i in 1:length(Msum)) {
      A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
      d <- d - commun[i] * A
    }
    
    # Next elements of the sum, until the expected precision
    k1 <- 1:k
    derive <- 0
    # vd <- vderive <- numeric()
    while (abs(d) > eps/10 & !is.nan(d)) {
      epsret <- signif(abs(d), 1)*10
      k <- k1[length(k1)]
      k1 <- k + (1:kstep)
      derive <- derive + d
      # vd <- c(vd, d); vderive <- c(vderive, derive)
      
      # # M: data.frame of the indices for the nested sums
      # M <- as.data.frame(matrix(nrow = 0, ncol = p))
      # if (p > 1) {
      #   for (i in 1:(p-1)) {
      #     Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
      #     M <- rbind( M, expand.grid(Mlist) )
      #     for (j in 1:(p-1)) {
      #       Mlist <- Mlist[c(p, 1:(p-1))]
      #       M <- rbind(M, expand.grid(Mlist))
      #     }
      #   }
      # }
      # M <- rbind( M, expand.grid(rep(list(k1), p)) )
      
      # M: data.frame of the indices for the nested sums
      M <- expand.grid(rep(list(k1), p))
      if (p > 1) {
        for (i in 1:(p-1)) {
          indsupp <- combn(p, i)
          for (j in 1:ncol(indsupp)) {
            jsupp <- indsupp[, j]
            Mlist <- vector("list", p)
            for (l in jsupp) Mlist[[l]] <- k1
            for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
            M <- rbind(M, expand.grid(Mlist))
          }
        }
      }
      
      Msum <- rowSums(M)
      
      Munique <- (max(Munique) + 1):max(M)
      Msumunique <- (max(Msumunique) + 1):max(Msum)
      
      # d <- 0
      # for (i in 1:length(Msum)) {
      #   commun <- prod(
      #     sapply(1:p, function(j) {
      #       pochhammer(0.5, M[i, j])*(1 - 1/lambda[j])^M[i, j]/factorial(M[i, j])
      #     })
      #   )
      #   A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
      #   d <- d - commun * A # / pochhammer((1 + p)/2, Msum[i])
      # }
      
      Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
      matpoch <- matrix(rep(Mpoch, p), ncol = p, byrow = FALSE)
      matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
      matlabmda <- matrix(rep(log(1 - 1/lambda), length(Munique)), ncol = p, byrow = TRUE)
      lambdaM <- matlabmda*matM
      matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
      
      matcommun <- rbind(matcommun, matcommunsupp)
      matcommunsupp  <- exp(matpoch + lambdaM - matfact)
      matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
      
      gridcommun <- expand.grid(matcommunsupp)
      if (p > 1) {
        for (i in 1:(p-1)) {
          indsupp <- combn(p, i)
          for (j in 1:ncol(indsupp)) {
            jsupp <- indsupp[, j]
            communlist <- vector("list", p)
            names(communlist) <- names(gridcommun)
            for (l in jsupp) communlist[[l]] <- matcommunsupp[, l]
            for (l in (1:p)[-jsupp]) communlist[[l]] <- matcommun[, l]
            gridcommun <- rbind(gridcommun, expand.grid(communlist))
          }
        }
      }
      commun <- apply(gridcommun, 1, prod)
      
      # d <- 0
      # for (i in 1:length(Msum)) {
      #   A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
      #   d <- d - commun[i] * A
      # }
      A <- sapply(Msum, function (i) sum(1/(0:(i-1) + (1+p)/2)))
      
      d <- -sum(commun*A)
      
    }
    
    # Next elements of the sum, with step=1, while not NaN
    # if (is.nan(d)) {
    #   k1 <- k
    #   d <- 0
    #   while (!is.nan(d)) {
    #     if (d > 0)
    #       epsret <- signif(abs(d), 1)*10
    #     k <- k1
    #     k1 <- k + 1
    #     derive <- derive + d
    #     
    #     # # M: data.frame of the indices for the nested sums
    #     # M <- as.data.frame(matrix(nrow = 0, ncol = p))
    #     # if (p > 1) {
    #     #   for (i in 1:(p-1)) {
    #     #     Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
    #     #     M <- rbind( M, expand.grid(Mlist) )
    #     #     for (j in 1:(p-1)) {
    #     #       Mlist <- Mlist[c(p, 1:(p-1))]
    #     #       M <- rbind(M, expand.grid(Mlist))
    #     #     }
    #     #   }
    #     # }
    #     # M <- rbind( M, rep(k1, p) )
    #     
    #     # M: data.frame of the indices for the nested sums
    #     M <- expand.grid(rep(list(k1), p))
    #     if (p > 1) {
    #       for (i in 1:(p-1)) {
    #         indsupp <- combn(p, i)
    #         for (j in 1:ncol(indsupp)) {
    #           jsupp <- indsupp[, j]
    #           Mlist <- vector("list", p)
    #           for (l in jsupp) Mlist[[l]] <- k1
    #           for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
    #           M <- rbind(M, expand.grid(Mlist))
    #         }
    #       }
    #     }
    #     
    #     Msum <- rowSums(M)
    #     
    #     d <- 0
    #     for (i in 1:length(Msum)) {
    #       commun <- prod(
    #         sapply(1:p, function(j) {
    #           pochhammer(0.5, M[i, j])*(1 - 1/lambda[j])^M[i, j]/factorial(M[i, j])
    #         })
    #       )
    #       A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
    #       d <- d - commun * A # / pochhammer((1 + p)/2, Msum[i])
    #     }
    #     
    #   }
    # }
    
    result <- as.numeric(-0.5 * log(prodlambda) - (1 + p)/2 * prod(lambda)^-0.5 * derive)
    
  } else { # equations 106 & 101
    
    # The first 5 elements of the sum
    
    # d <- 0
    # for (i in 1:length(Msum)) {
    #   commun <- prod(
    #     sapply(1:(p-1), function(j) {
    #       pochhammer(0.5, M[i, j])*(1 - lambda[j]/lambda[p])^M[i, j]/factorial(M[i, j])
    #     })
    #   )
    #   commun <- commun*(1 - 1/lambda[p])^M[i, p]/factorial(M[i, p])
    #   d <- d + commun * pochhammer(0.5, M[i, p])*pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
    # }
    
    Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
    matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), 0)
    matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
    matlabmda <- matrix(rep(log(1 - c(lambda[-p], 1)/lambda[p]), length(Munique)), ncol = p, byrow = TRUE)
    lambdaM <- matlabmda*matM
    matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
    
    matcommun  <- exp(matpoch + lambdaM - matfact)
    matcommun <- as.data.frame(apply(matcommun, 2, Re))
    matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
    colnames(matcommunsupp) <- colnames(matcommunsupp)
    
    gridcommun <- expand.grid(matcommun)
    gridcommun <- gridcommun[-1, , drop = FALSE]
    commun <- apply(gridcommun, 1, prod)
    
    d <- sum(
      commun * sapply(Msum, function(i) {
        pochhammer(0.5, i) * pochhammer(1, i) / ( pochhammer((1 + p)/2, i) * i )
      })
    )
    
    # Next elements of the sum, until the expected precision
    k1 <- 1:k
    derive <- 0
    while (abs(d) > eps/10 & !is.nan(d)) {
      epsret <- signif(abs(d), 1)*10
      k <- k1[length(k1)]
      k1 <- k + (1:kstep)
      derive <- derive + d
      
      # # M: data.frame of the indices for the nested sums
      # M <- as.data.frame(matrix(nrow = 0, ncol = p))
      # if (p > 1) {
      #   for (i in 1:(p-1)) {
      #     Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
      #     M <- rbind( M, expand.grid(Mlist) )
      #     for (j in 1:(p-1)) {
      #       Mlist <- Mlist[c(p, 1:(p-1))]
      #       M <- rbind(M, expand.grid(Mlist))
      #     }
      #   }
      # }
      # M <- rbind( M, expand.grid(rep(list(k1), p)) )
      
      # M: data.frame of the indices for the nested sums
      M <- expand.grid(rep(list(k1), p))
      if (p > 1) {
        for (i in 1:(p-1)) {
          indsupp <- combn(p, i)
          for (j in 1:ncol(indsupp)) {
            jsupp <- indsupp[, j]
            Mlist <- vector("list", p)
            for (l in jsupp) Mlist[[l]] <- k1
            for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
            M <- rbind(M, expand.grid(Mlist))
          }
        }
      }
      
      Msum <- rowSums(M)
      
      Munique <- (max(Munique)+1):max(M)
      Msumunique <- (max(Msumunique)+1):max(Msum)
      
      # d <- 0
      # for (i in 1:length(Msum)) {
      #   commun <- prod(
      #     sapply(1:(p-1), function(j) {
      #       pochhammer(0.5, M[i, j])*(1 - lambda[j]/lambda[p])^M[i, j]/factorial(M[i, j])
      #     })
      #   )
      #   commun <- commun*(1 - 1/lambda[p])^M[i, p]/factorial(M[i, p])
      #   d <- d + commun * pochhammer(0.5, M[i, p])*pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
      # }
      
      Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
      matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), 0)
      matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
      matlabmda <- matrix(rep(log(1 - c(lambda[-p], 1)/lambda[p]), length(Munique)), ncol = p, byrow = TRUE)
      lambdaM <- matlabmda*matM
      matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
      
      matcommun <- rbind(matcommun, matcommunsupp)
      matcommunsupp  <- exp(matpoch + lambdaM - matfact)
      matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
      
      gridcommun <- expand.grid(matcommunsupp)
      if (p > 1) {
        for (i in 1:(p-1)) {
          indsupp <- combn(p, i)
          for (j in 1:ncol(indsupp)) {
            jsupp <- indsupp[, j]
            communlist <- vector("list", p)
            names(communlist) <- names(gridcommun)
            for (l in jsupp) communlist[[l]] <- matcommunsupp[, l]
            for (l in (1:p)[-jsupp]) communlist[[l]] <- matcommun[, l]
            gridcommun <- rbind(gridcommun, expand.grid(communlist))
          }
        }
      }
      commun <- apply(gridcommun, 1, prod)
      
      d <- sum(
        commun * sapply(Msum, function(i) {
          pochhammer(0.5, i) * pochhammer(1, i) / ( pochhammer((1 + p)/2, i) * i )
        })
      )
      
    }
    
    # # Next elements of the sum, with step=1, while not NaN
    # if (is.nan(d)) {
    #   k1 <- k
    #   d <- 0
    #   while (!is.nan(d)) {
    #     if (d > 0)
    #       epsret <- signif(abs(d), 1)*10
    #     k <- k1
    #     k1 <- k + 1
    #     derive <- derive + d
    #     
    #     # # M: data.frame of the indices for the nested sums
    #     # M <- as.data.frame(matrix(nrow = 0, ncol = p))
    #     # if (p > 1) {
    #     #   for (i in 1:(p-1)) {
    #     #     Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
    #     #     M <- rbind( M, expand.grid(Mlist) )
    #     #     for (j in 1:(p-1)) {
    #     #       Mlist <- Mlist[c(p, 1:(p-1))]
    #     #       M <- rbind(M, expand.grid(Mlist))
    #     #     }
    #     #   }
    #     # }
    #     # M <- rbind( M, rep(k1, p) )
    #     
    #     # M: data.frame of the indices for the nested sums
    #     M <- expand.grid(rep(list(k1), p))
    #     if (p > 1) {
    #       for (i in 1:(p-1)) {
    #         indsupp <- combn(p, i)
    #         for (j in 1:ncol(indsupp)) {
    #           jsupp <- indsupp[, j]
    #           Mlist <- vector("list", p)
    #           for (l in jsupp) Mlist[[l]] <- k1
    #           for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
    #           M <- rbind(M, expand.grid(Mlist))
    #         }
    #       }
    #     }
    #     
    #     Msum <- rowSums(M)
    #     
    #     d <- 0
    #     for (i in 1:length(Msum)) {
    #       commun <- prod(
    #         sapply(1:(p-1), function(j) {
    #           pochhammer(0.5, M[i, j])*(1 - lambda[j]/lambda[p])^M[i, j]/factorial(M[i, j])
    #         })
    #       )
    #       commun <- commun*(1 - 1/lambda[p])^M[i, p]/factorial(M[i, p])
    #       d <- d + commun * pochhammer(0.5, M[i, p])*pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
    #     }
    #     
    #   }
    # }
    
    result <- as.numeric(-0.5 * log(prodlambda) + (1 + p)/2 * (log(lambda[p]) - derive))
    
  }
  
  # if (is.nan(d)) {
  #   epsret <- signif(epsret, 1)
  #   warning("Cannot reach the precision ", eps, " due to NaN\n",
  #           "Number of iteration: ", k, "\n",
  #           "Precision reached:", epsret)
  #   attr(result, "epsilon") <- epsret
  # } else {
  attr(result, "epsilon") <- eps
  # }
  attr(result, "k") <- k
  
  return(result)
}