Method for the computation of the elements of the sum changed: it now runs...

Method for the computation of the elements of the sum changed: it now runs faster and, in the "while" loop, the table M for the combinations of the indices is corrected (there was an error for dim > 3)

R/lauricella.R

@@ -33,6 +33,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   #' @references N. Bouhlel, A. Dziri, Kullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions.
   #' IEEE Signal Processing Letters, vol. 26 no. 7, July 2019.
   #' \doi{10.1109/LSP.2019.2915000}
+  #' @importFrom utils combn
   #' @export
   
   # # Internal functions
@@ -56,7 +57,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   #     return(sum(log(abs(y)) + (y < 0)*pi*1i))
   #   }
   # }
-  # 
+  
   # Number of variables
   n <- length(x)
   
@@ -68,12 +69,16 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   if (any(abs(x) >= 1))
     stop("The series does not converge for these x values.")
   
+  access <- function(ind, tab) {
+    sapply(1:n, function(ii) tab[ind[ii] + 1, ii])
+  }
+  
   k <- 5
   
-  # M: data.frame of the indices for the nested sums
-  # (i.e. all arrangements of n elements from {0:k})
+  # # 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), n))
-  
+
   Msum <- apply(M, 1, sum)
   
   Munique <- 0:max(M)
@@ -82,27 +87,36 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   # Product x^{m_1}/m_1! * ... * x^{m_n}/m_n! for m_1 = 0...k, ...,  m_n = 0...k
   # xfact <- apply(M, 1, function(Mi) prod( x^Mi ))
   # lnfact <- apply(M, 1, function(Mi) sum(sapply(Mi, lnfactorial)))
-  xfact <- matrix(nrow = max(M) + 1, ncol = n, dimnames = list(Munique, 1:n))
+  xfact <- as.data.frame(
+    matrix(nrow = max(M) + 1, ncol = n, dimnames = list(Munique, 1:n)))
   for (i in Munique) for (j in 1:n) {
     # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
     xfact[i+1, j] <- x[j]^i
   }
-  prodxfact <- function(ind) {
-    prod(mapply(function(i, j) xfact[i, j], ind+1, 1:n))
-  }
+  # prodxfact <- function(ind) {
+  #   # prod(mapply(function(i, j) xfact[i, j], ind+1, 1:n))
+  #   prod(access(ind, xfact))
+  # }
+  gridxfact <- expand.grid(xfact)
+  prodxfact <- apply(gridxfact, 1, prod)
   
-  lnfact <- matrix(sapply(Munique, lnfactorial), nrow = length(Munique),
-                   ncol = n, byrow = FALSE, dimnames = list(Munique, 1:n))
-  sumlnfact <- function(ind) {
-    sum(mapply(function(i, j) lnfact[i, j], ind+1, 1:n))
-  }
+  lnfact <- as.data.frame(
+    matrix(sapply(Munique, lnfactorial), nrow = length(Munique),
+           ncol = n, byrow = FALSE, dimnames = list(Munique, 1:n)))
+  # sumlnfact <- function(ind) {
+  #   # sum(mapply(function(i, j) lnfact[i, j], ind+1, 1:n))
+  #   sum(access(ind, lnfact))
+  # }
+  gridlnfact <- expand.grid(as.data.frame(lnfact))
+  sumlnfact <- rowSums(gridlnfact)
   
   # pochhammer(a, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
   # lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
   lnapoch <- sapply(Msumunique, function(j) lnpochhammer(a, j))
   names(lnapoch) <- Msumunique
 
-  lnbpoch <- matrix(nrow = max(M) + 1, ncol = n, dimnames = list(Munique, 1:n))
+  lnbpoch <- as.data.frame(
+    matrix(nrow = max(M) + 1, ncol = n, dimnames = list(Munique, 1:n)))
   for (i in Munique) for (j in 1:n) {
     # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
     lnbpoch[i+1, j] <- lnpochhammer(b[j], i)
@@ -113,9 +127,12 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   #   matlnbpoch[i+1, j] <- lnpochhammer(b[j], Munique[i+1])
   # }
   # lnbpoch <- data.frame(stack(matlnbpoch), M = 0:maxM)
-  sumlnpochb <- function(ind) {
-    sum(mapply(function(i, j) lnbpoch[i, j], ind+1, 1:n))
-  }
+  # sumlnpochb <- function(ind) {
+  #   # sum(mapply(function(i, j) lnbpoch[i, j], ind+1, 1:n))
+  #   sum(access(ind, lnbpoch))
+  # }
+  gridlnbpoch <- expand.grid(as.data.frame(lnbpoch))
+  sumlnpochb <- rowSums(gridlnbpoch)
   
   # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
   # lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
@@ -124,15 +141,30 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
   
   # res1 <- lnapoch + lnbpoch - lngpoch - lnfact
   # res2 <- sum(xfact * exp(res1))
-  res1 <- lnapoch[Msum+1] + apply(M, 1, sumlnpochb) - lngpoch[Msum+1] -
-    apply(M, 1, sumlnfact)
-  res2 <- sum( apply(M, 1, prodxfact) * exp(res1) )
+  res1 <- lnapoch[Msum+1] + sumlnpochb - lngpoch[Msum+1] - sumlnfact
+  # res1 <- lnapoch[Msum+1] + sumbpoch - lngpoch[Msum+1] -
+  #   apply(M, 1, sumlnfact)
+  res2 <- sum( prodxfact * exp(res1) )
+  # res2 <- sum( apply(M, 1, prodxfact) * exp(res1) )
   
   kstep <- 5
   
   k1 <- 1:k
   result <- 0
   
+  # prodxfact <- function(ind) {
+  #   # prod(mapply(function(i, j) xfact[i, j], ind+1, 1:n))
+  #   prod(access(ind, xfact))
+  # }
+  # sumlnfact <- function(ind) {
+  #   # sum(mapply(function(i, j) lnfact[i, j], ind+1, 1:n))
+  #   sum(access(ind, lnfact))
+  # }
+  # sumlnpochb <- function(ind) {
+  #   # sum(mapply(function(i, j) lnbpoch[i, j], ind+1, 1:n))
+  #   sum(access(ind, lnbpoch))
+  # }
+  
   while (abs(res2) > eps/10 & !is.nan(res2)) {
     
     epsret <- signif(abs(res2), 1)*10
@@ -142,21 +174,35 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
     k1 <- k + (1:kstep)
     result <- result + res2
     
-    # M: data.frame of the indices for the nested sums
-    M <- as.data.frame(matrix(nrow = 0, ncol = n))
+    # # M: data.frame of the indices for the nested sums
+    # M <- expand.grid(rep(list(k1), n))
+    # if (n > 1) {
+    #   for (i in 1:(n-1)) {
+    #     Mlist <- c( rep(list(0:k), n-i), rep(list(k1), i) )
+    #     M <- rbind( M, expand.grid(Mlist) )
+    #     for (j in 1:(n-1)) {
+    #       Mlist <- Mlist[c(n, 1:(n-1))]
+    #       M <- rbind(M, expand.grid(Mlist))
+    #     }
+    #   }
+    # }
+    # M1 <- M
+    M <- expand.grid(rep(list(k1), n))
     if (n > 1) {
       for (i in 1:(n-1)) {
-        Mlist <- c( rep(list(0:k), n-i), rep(list(k1), i) )
-        M <- rbind( M, expand.grid(Mlist) )
-        for (j in 1:(n-1)) {
-          Mlist <- Mlist[c(n, 1:(n-1))]
+        indsupp <- combn(n, i)
+        for (j in 1:ncol(indsupp)) {
+          jsupp <- indsupp[, j]
+          Mlist <- vector("list", n)
+          for (l in jsupp) Mlist[[l]] <- k1
+          for (l in (1:n)[-jsupp]) Mlist[[l]] <- 0:k
           M <- rbind(M, expand.grid(Mlist))
         }
       }
     }
-    M <- rbind( M, expand.grid(rep(list(k1), n)) )
+    # M2 <- M
     
-    Msum <- apply(M, 1, sum)
+    Msum <- rowSums(M)
 
     Munique <- (max(Munique)+1):max(M)
     Msumunique <- (max(Msumunique)+1):max(Msum)
@@ -164,22 +210,44 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
     # Product x^{m_1}/m_1! * ... * x^{m_n}/m_n! for m_1, ...,  m_n given by the rows of M
     # xfact <- apply(M, 1, function(Mi) prod( x^Mi ))
     # lnfact <- apply(M, 1, function(Mi) sum(sapply(Mi, lnfactorial)))
-    xfactsupp <- matrix(nrow = length(Munique), ncol = n, dimnames = list(Munique, 1:n))
+    xfactsupp <- as.data.frame(
+      matrix(nrow = length(Munique), ncol = n, dimnames = list(Munique, 1:n)))
     for (i in 1:length(Munique)) for (j in 1:n) {
       # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
       xfactsupp[i, j] <- x[j]^Munique[i]
     }
-    xfact <- rbind(xfact, xfactsupp)
-    # prodxfact <- function(ind) {
-    #   prod(mapply(function(i, j) xfact[i, j], ind+1, 1:n))
-    # }
     
-    lnfactsupp <- matrix(sapply(Munique, lnfactorial), nrow = length(Munique),
-                         ncol = n, dimnames = list(Munique, 1:n))
-    lnfact <- rbind(lnfact, lnfactsupp)
-    # sumlnfact <- function(ind) {
-    #   sum(mapply(function(i, j) lnfact[i, j], ind+1, 1:n))
-    # }
+    gridxfact <- expand.grid(xfactsupp)
+    for (i in 1:(n-1)) {
+      indsupp <- combn(n, i)
+      for (j in 1:ncol(indsupp)) {
+        jsupp <- indsupp[, j]
+        xfactlist <- vector("list", n)
+        names(xfactlist) <- names(gridxfact)
+        xfactlist[jsupp] <- xfactsupp[jsupp]
+        xfactlist[-jsupp] <- xfact[-jsupp]
+        gridxfact <- rbind(gridxfact, expand.grid(xfactlist))
+      }
+    }
+    prodxfact <- apply(gridxfact, 1, prod)
+    
+    lnfactsupp <- as.data.frame(
+      matrix(sapply(Munique, lnfactorial), nrow = length(Munique),
+             ncol = n, dimnames = list(Munique, 1:n)))
+    
+    gridlnfact <- expand.grid(lnfactsupp)
+    for (i in 1:(n-1)) {
+      indsupp <- combn(n, i)
+      for (j in 1:ncol(indsupp)) {
+        jsupp <- indsupp[, j]
+        lnfactlist <- vector("list", n)
+        names(lnfactlist) <- names(gridlnfact)
+        lnfactlist[jsupp] <- lnfactsupp[jsupp]
+        lnfactlist[-jsupp] <- lnfact[-jsupp]
+        gridlnfact <- rbind(gridlnfact, expand.grid(lnfactlist))
+      }
+    }
+    sumlnfact <- rowSums(gridlnfact)
     
     # pochhammer(a, m_1+...+m_n) for m_1, ...,  m_n given by the rows of M
     # lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
@@ -192,15 +260,26 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
     #   # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
     #   lnbpoch[i] <- sum( sapply(1:n, function(j) lnpochhammer(b[j], M[i, j])) )
     # }
-    lnbpochsupp <- matrix(nrow = length(Munique), ncol = n, dimnames = list(Munique, 1:n))
+    lnbpochsupp <- as.data.frame(
+      matrix(nrow = length(Munique), ncol = n, dimnames = list(Munique, 1:n)))
     for (i in 1:length(Munique)) for (j in 1:n) {
       # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
       lnbpochsupp[i, j] <- lnpochhammer(b[j], Munique[i])
     }
-    lnbpoch <- rbind(lnbpoch, lnbpochsupp)
-    sumlnpochb <- function(ind) {
-      sum(mapply(function(i, j) lnbpoch[i, j], ind+1, 1:n))
+    
+    gridlnbpoch <- expand.grid(lnbpochsupp)
+    for (i in 1:(n-1)) {
+      indsupp <- combn(n, i)
+      for (j in 1:ncol(indsupp)) {
+        jsupp <- indsupp[, j]
+        lnbpochlist <- vector("list", n)
+        names(lnbpochlist) <- names(gridlnbpoch)
+        lnbpochlist[jsupp] <- lnbpochsupp[jsupp]
+        lnbpochlist[-jsupp] <- lnbpoch[-jsupp]
+        gridlnbpoch <- rbind(gridlnbpoch, expand.grid(lnbpochlist))
+      }
     }
+    sumlnpochb <- rowSums(gridlnbpoch)
     
     # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ...,  m_p = 0...k
     # lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
@@ -210,8 +289,14 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
     
     # res1 <- lnapoch + lnbpoch - lngpoch - lnfact
     # res2 <- sum(xfact * exp(res1))
-    res1 <- lnapoch[Msum] + apply(M, 1, sumlnpochb) - lngpoch[Msum] - apply(M, 1, sumlnfact)
-    res2 <- sum( apply(M, 1, prodxfact) * exp(res1) )
+    # res1 <- lnapoch[Msum+1] + apply(M, 1, sumlnpochb) - lngpoch[Msum+1] - apply(M, 1, sumlnfact)
+    # res2 <- sum( apply(M, 1, prodxfact) * exp(res1) )
+    res1 <- lnapoch[Msum+1] + sumlnpochb - lngpoch[Msum+1] - sumlnfact
+    res2 <- sum( prodxfact * exp(res1) )
+    
+    xfact <- rbind(xfact, xfactsupp)
+    lnfact <- rbind(lnfact, lnfactsupp)
+    lnbpoch <- rbind(lnbpoch, lnbpochsupp)
   }
   
   result <- Re(result)