There was a problem if lambda[p]==1: Problem when computing log(1-lambda[p])

R/kldcauchy.R

@@ -115,6 +115,16 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
   
   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
     
@@ -129,11 +139,11 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
     # }
     
     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 - lambda), length(Munique)), ncol = p, byrow = TRUE)
+    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), p), ncol = p, byrow = FALSE)
+    matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
     
     matcommun  <- exp(matpoch + lambdaM - matfact)
     matcommun <- as.data.frame(apply(matcommun, 2, Re))
@@ -177,15 +187,15 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
       # 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)
+      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", p)
+            Mlist <- vector("list", p1)
             for (l in jsupp) Mlist[[l]] <- k1
-            for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
+            for (l in (1:p1)[-jsupp]) Mlist[[l]] <- 0:k
             M <- rbind(M, expand.grid(Mlist))
           }
         }
@@ -208,26 +218,26 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
       # }
       
       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 - lambda), length(Munique)), ncol = p, byrow = TRUE)
+      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), p), ncol = p, byrow = FALSE)
+      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 (p > 1) {
-        for (i in 1:(p-1)) {
-          indsupp <- combn(p, i)
+      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", p)
+            communlist <- vector("list", p1)
             names(communlist) <- names(gridcommun)
             for (l in jsupp) communlist[[l]] <- matcommunsupp[, l]
-            for (l in (1:p)[-jsupp]) communlist[[l]] <- matcommun[, l]
+            for (l in (1:p1)[-jsupp]) communlist[[l]] <- matcommun[, l]
             gridcommun <- rbind(gridcommun, expand.grid(communlist))
           }
         }