diff --git a/R/lauricella.R b/R/lauricella.R
index 459631bfcf0174732a8007eb72031a68897e0762..9c45b46d8aad0f6b1c42bf054ef13388a0316808 100644
--- a/R/lauricella.R
+++ b/R/lauricella.R
@@ -53,8 +53,8 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
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))
# Sum of the indices
@@ -88,14 +88,17 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# # sum(mapply(function(i, j) lnfact[i, j], ind+1, 1:n))
# sum(access(ind, lnfact))
# }
+ # Table of all combinations of the log(m_i)
gridlnfact <- expand.grid(as.data.frame(lnfact))
+ # Sum of the logarithms
sumlnfact <- rowSums(gridlnfact)
- # pochhammer(a, m_1+...+m_n) for m_1 = 0...k, ..., m_p = 0...k
+ # Logarithms of pochhammer(a, m_1+...+m_n) for m_1 = 0...k, ..., m_n = 0...k
# lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
lnapoch <- sapply(Msumunique, function(j) lnpochhammer(a, j))
names(lnapoch) <- Msumunique
+ # Logarithms of pochhammer(b_i, m_i) for m_1 = 0...k, ..., m_n = 0...k
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) {
@@ -112,10 +115,12 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# # sum(mapply(function(i, j) lnbpoch[i, j], ind+1, 1:n))
# sum(access(ind, lnbpoch))
# }
+ # Table of all combinations of the log(pochhammer(b_i, m_i))
gridlnbpoch <- expand.grid(as.data.frame(lnbpoch))
+ # Sum of the logarithms
sumlnpochb <- rowSums(gridlnbpoch)
- # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ..., m_p = 0...k
+ # Logarithms of pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ..., m_n = 0...k
# lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
lngpoch <- sapply(Msumunique, function(j) lnpochhammer(g, j))
names(lngpoch) <- Msumunique
@@ -168,6 +173,8 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# }
# }
# M1 <- M
+
+ # 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)) {
@@ -183,12 +190,13 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
}
# M2 <- M
+ # Sum of the indices
Msum <- rowSums(M)
Munique <- (max(Munique)+1):max(M)
Msumunique <- (max(Msumunique)+1):max(Msum)
- # Product x^{m_1}/m_1! * ... * x^{m_n}/m_n! for m_1, ..., m_n given by the rows of M
+ # Product x^{m_1} * ... * x^{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 <- as.data.frame(
@@ -197,7 +205,6 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
xfactsupp[i, j] <- x[j]^Munique[i]
}
-
gridxfact <- expand.grid(xfactsupp)
for (i in 1:(n-1)) {
indsupp <- combn(n, i)
@@ -212,10 +219,11 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
}
prodxfact <- apply(gridxfact, 1, prod)
+ # Logarithm of the product m_1! * ... * m_n! for m_1, ..., m_n given by the rows of M
+ # i.e. \sum_{i=0}^n{\log{m_i!}}
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)
@@ -230,7 +238,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
}
sumlnfact <- rowSums(gridlnfact)
- # pochhammer(a, m_1+...+m_n) for m_1, ..., m_n given by the rows of M
+ # Logarithms of 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))
lnapochsupp <- sapply(Msumunique, function(j) lnpochhammer(a, j))
names(lnapochsupp) <- Msumunique
@@ -241,13 +249,14 @@ 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])) )
# }
+ # Logarithms of pochhammer(b_i, m_i) for m_1, ..., m_n given by the rows of M
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])
}
-
+ # Table of all combinations of the log(pochhammer(b_i, m_i))
gridlnbpoch <- expand.grid(lnbpochsupp)
for (i in 1:(n-1)) {
indsupp <- combn(n, i)
@@ -260,9 +269,10 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
gridlnbpoch <- rbind(gridlnbpoch, expand.grid(lnbpochlist))
}
}
+ # Sum of the logarithms
sumlnpochb <- rowSums(gridlnbpoch)
- # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ..., m_p = 0...k
+ # Logarithms of pochhammer(g, m_1+...+m_n) for m_1, ..., m_n given by the rows of M
# lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
lngpochsupp <- sapply(Msumunique, function(j) lnpochhammer(g, j))
names(lngpochsupp) <- Msumunique
@@ -275,6 +285,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
res1 <- lnapoch[Msum+1] + sumlnpochb - lngpoch[Msum+1] - sumlnfact
res2 <- sum( prodxfact * exp(res1) )
+ # Add the new calculated values to the tables of former values
xfact <- rbind(xfact, xfactsupp)
lnfact <- rbind(lnfact, lnfactsupp)
lnbpoch <- rbind(lnbpoch, lnbpochsupp)