Modifications to reduce the execution time.

R/dlauricella.R

-#' @importFrom data.table rbindlist
-
 dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
-
+  
   p <- length(lambda)
-
+  
   buildMlist <- function(j, isupp, k, k1, p = p) {
     jsupp <- isupp[, j]
     Ml <- rep(list(0:k), p)
@@ -11,7 +9,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       Ml[[j]] <- k1
     return(expand.grid(Ml))
   }
-
+  
   buildxlist <- function(j, isupp, x, xsupp, p = p) {
     jsupp <- isupp[, j]
     xl <- rep(list(x), p)
@@ -19,29 +17,29 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
     xl[-jsupp] <- x[-jsupp]
     return(expand.grid(xl))
   }
-
+  
   prodlambda <- prod(lambda)
-
+  
   lambdanu <- lambda*nu1/nu2
   prodlambdanu <- prod(lambdanu)
-
+  
   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]
-
+  
   Munique <- 0:k
-
+  
   # Sum of the indices
   Msum <- rowSums(M)
-  # Msumunique <- 0:max(Msum)
-
+  Msumunique <- 0:max(Msum)
+  
   kstep <- 5
-
+  
   if (lambdanu[p] <= 1) {  # lambda[1] < ... < lambda[p] < 1
-
+    
     if (p > 1 & lambda[p] == 1) {
       lambda <- lambda[1:(p-1)]
       p1 <- p-1
@@ -52,25 +50,25 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
     } else {
       p1 <- p
     }
-
+    
     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 - lambdanu[1:p1]), 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(matpoch + lambdaM - matfact)
     # matcommun <- as.data.frame(apply(matcommun, 2, Re))
     matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
     colnames(matcommunsupp) <- colnames(matcommun)
-
+    
     gridcommun <- expand.grid(matcommun)
     gridcommun <- gridcommun[-1, , drop = FALSE]
     # commun <- apply(gridcommun, 1, prod)
     commun <- rowSums(gridcommun)
-
+    
     # d <- sum(
     #   commun * sapply(Msum, function(i) {
     #     pochhammer(1, i) / ( pochhammer((nu1 + p)/2, i) * i )
@@ -81,7 +79,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
         lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
       })
     )))
-
+    
     # Next elements of the sum, until the expected precision
     k1 <- 1:k
     derive <- 0
@@ -90,7 +88,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       k <- k1[length(k1)]
       k1 <- k + (1:kstep)
       derive <- derive + d
-
+      
       # M: data.frame of the indices for the nested sums
       # M <- expand.grid(rep(list(k1), p1))
       # if (p1 > 1) {
@@ -117,12 +115,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
         }
         M <- data.frame(rbindlist(Mlist))
       }
-
+      
       # Sum of the indices
       Msum <- rowSums(M)
-
+      
       Munique <- (max(Munique)+1):max(M)
-      # Msumunique <- (max(Msumunique)+1):max(Msum)
+      Msumunique <- unique(Msum)
+      names(Msumunique) <- as.character(Msumunique)
 
       Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
       matpoch <- matrix(rep(Mpoch, p1), ncol = p1, byrow = FALSE)
@@ -130,13 +129,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       matlabmda <- matrix(rep(log(1 - lambdanu[1:p1]), 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(matpoch + lambdaM - matfact)
       # matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
       colnames(matcommunsupp) <- colnames(matcommun)
-
+      
       communlist <- list(expand.grid(matcommunsupp))
       if (p1 == 1) {
         gridcommun <- communlist[[1]]
@@ -152,49 +151,53 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       }
       # commun <- apply(gridcommun, 1, prod)
       commun <- rowSums(gridcommun)
-
+      
       # d <- sum(
       #   commun * sapply(Msum, function(i) {
       #     pochhammer(1, i) / ( pochhammer((nu1 + p)/2, i) * i )
       #   })
       # )
-      d <- Re(sum(exp(
-        commun + sapply(Msum, function(i) {
-          lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
-        })
-      )))
-
+      # d <- Re(sum(exp(
+      #   commun + sapply(Msum, function(i) {
+      #     lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
+      #   })
+      # )))
+      pochMsum <- sapply(Msumunique, function(i) {
+        lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
+      })
+      names(pochMsum) <- as.character(Msumunique)
+      d <- Re(sum(exp(commun + pochMsum[as.character(Msum)])))
     }
-
+    
     derive <- as.numeric(-0.5 * log(prodlambda) - (nu2 + p)/2 * derive)
-
+    
   } else if (lambdanu[1] > 1) { # 1 < lambda[1] < ... < lambda[p]
-
+    
     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/lambdanu), 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(matpoch + lambdaM - matfact)
     # matcommun <- as.data.frame(apply(matcommun, 2, Re))
     matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
     colnames(matcommunsupp) <- colnames(matcommun)
-
+    
     gridcommun <- expand.grid(matcommun)
     gridcommun <- gridcommun[-1, , drop = FALSE]
     # commun <- apply(gridcommun, 1, prod)
     commun <- rowSums(gridcommun)
-
+    
     d <- 0
     for (i in 1:length(Msum)) {
       A <- sum(1/(0:(Msum[i]-1) + (nu1+p)/2))
       d <- d - exp(commun[i]) * A
     }
     d <- Re(d)
-
+    
     # Next elements of the sum, until the expected precision
     k1 <- 1:k
     derive <- 0
@@ -205,7 +208,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       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) {
@@ -219,7 +222,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       #   }
       # }
       # 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) {
@@ -246,12 +249,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
         }
         M <- data.frame(rbindlist(Mlist))
       }
-
+      
       Msum <- rowSums(M)
-
+      
       Munique <- (max(Munique) + 1):max(M)
-      # Msumunique <- (max(Msumunique) + 1):max(Msum)
-
+      Msumunique <- unique(Msum)
+      names(Msumunique) <- as.character(Msumunique)
+      
       # d <- 0
       # for (i in 1:length(Msum)) {
       #   commun <- prod(
@@ -262,20 +266,20 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       #   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/lambdanu), 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(matpoch + lambdaM - matfact)
       # matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
       colnames(matcommunsupp) <- colnames(matcommun)
-
+      
       communlist <- list(expand.grid(matcommunsupp))
       if (p == 1) {
         gridcommun <- communlist[[1]]
@@ -291,22 +295,29 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       }
       # commun <- apply(gridcommun, 1, prod)
       commun <- rowSums(gridcommun)
-
+      
       # 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) + (nu1+p)/2)))
+      # d <- commun + pochMsum[Msum]
+      # A <- sapply(Msum, function (i) sum(1/(0:(i-1) + (nu1+p)/2)))
+      pochMsum <- sapply(Msumunique, function(i) {
+        sum(1/(0:(i-1) + (nu1+p)/2))
+      })
+      names(pochMsum) <- as.character(Msumunique)
 
+      A <- pochMsum[as.character(Msum)]
+      
       d <- Re(-sum(exp(commun)*A))
-
+      
     }
-
+    
     derive <- as.numeric(-0.5 * log(prodlambda) - (nu2 + p)/2 * prod(1/sqrt(lambdanu)) * derive)
-
+    
   } else { # lambda[1] < ... < 1 < ... < lambda[p]
-
+    
     Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
     matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), 0)
     # Mpochp <- sapply(Munique, function(j) lnpochhammer(nu1/2, j))
@@ -315,18 +326,18 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
     matlabmda <- matrix(rep(log(1 - c(lambdanu[-p], 1)/lambdanu[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(matpoch + lambdaM - matfact)
     # matcommun <- as.data.frame(apply(matcommun, 2, Re))
     matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
     colnames(matcommunsupp) <- colnames(matcommun)
-
+    
     gridcommun <- expand.grid(matcommun)
     gridcommun <- gridcommun[-1, , drop = FALSE]
     # commun <- apply(gridcommun, 1, prod)
     commun <- rowSums(gridcommun)
-
+    
     # d <- sum(
     #   commun * sapply(1:length(Msum), function(i) {
     #     pochhammer(0.5, M[i, p]) * pochhammer(1, Msum[i]) / ( pochhammer((nu1 + p)/2, Msum[i]) * Msum[i] )
@@ -338,7 +349,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
           ( lnpochhammer((nu1 + p)/2, Msum[i]) + log(Msum[i]) )
       })
     )))
-
+    
     # Next elements of the sum, until the expected precision
     k1 <- 1:k
     derive <- 0
@@ -347,7 +358,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       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) {
@@ -361,7 +372,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       #   }
       # }
       # M <- rbind( M, expand.grid(rep(list(k1), p)) )
-
+      
       # M <- expand.grid(rep(list(k1), p))
       # if (p > 1) {
       #   for (i in 1:(p-1)) {
@@ -375,7 +386,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       #     }
       #   }
       # }
-
+      
       # M: data.frame of the indices for the nested sums
       Mlist <- list(expand.grid(rep(list(k1), p)))
       if (p == 1) {
@@ -389,12 +400,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
         }
         M <- data.frame(rbindlist(Mlist))
       }
-
+      
       Msum <- rowSums(M)
-
+      
       Munique <- (max(Munique)+1):max(M)
-      # Msumunique <- (max(Msumunique)+1):max(Msum)
-
+      Msumunique <- unique(Msum)
+      names(Msumunique) <- as.character(Msumunique)
+      
       # d <- 0
       # for (i in 1:length(Msum)) {
       #   commun <- prod(
@@ -405,7 +417,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       #   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)
       # Mpochp <- sapply(Munique, function(j) lnpochhammer(nu1/2, j))
@@ -414,13 +426,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       matlabmda <- matrix(rep(log(1 - c(lambdanu[-p], 1)/lambdanu[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(matpoch + lambdaM - matfact)
       # matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
       colnames(matcommunsupp) <- colnames(matcommun)
-
+      
       communlist <- list(expand.grid(matcommunsupp))
       if (p == 1) {
         gridcommun <- communlist[[1]]
@@ -436,27 +448,38 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
       }
       # commun <- apply(gridcommun, 1, prod)
       commun <- rowSums(gridcommun)
-
+      
       # d <- sum(
       #   commun * sapply(1:length(Msum), function(i) {
       #     pochhammer(0.5, M[i, p]) * pochhammer(1, Msum[i]) / ( pochhammer((nu1 + p)/2, Msum[i]) * Msum[i] )
       #   })
       # )
+      # d <- Re(sum(exp(
+      #   commun + sapply(1:length(Msum), function(i) {
+      #     lnpochhammer(nu1/2, M[i, p]) + lnpochhammer(1, Msum[i]) -
+      #       ( lnpochhammer((nu1 + p)/2, Msum[i]) + log(Msum[i]) )
+      #   })
+      # )))
+      pochnup <- sapply(0:max(M), function(i)  {
+        lnpochhammer(nu1/2, i)
+      })
+      names(pochnup) <- as.character(0:max(M))
+      pochMsum <- sapply(Msumunique, function(i) {
+        lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
+      })
+      names(pochMsum) <- as.character(Msumunique)
       d <- Re(sum(exp(
-        commun + sapply(1:length(Msum), function(i) {
-          lnpochhammer(nu1/2, M[i, p]) + lnpochhammer(1, Msum[i]) -
-            ( lnpochhammer((nu1 + p)/2, Msum[i]) + log(Msum[i]) )
-        })
+        commun + pochMsum[as.character(Msum)] + pochnup[as.character(M[, p])]
       )))
-
+      
     }
-
+    
     derive <- as.numeric(-0.5 * log(prodlambda) + (nu2 + p)/2 * (log(lambdanu[p]) - derive))
-
+    
   }
-
+  
   attr(derive, "epsilon") <- eps
   attr(derive, "k") <- k
-
+  
   return(derive)
 }