From c115e924c11d8cd40aec834d36afa443c1372cb4 Mon Sep 17 00:00:00 2001
From: Pierre Santagostini <pierre.santagostini@agrocampus-ouest.fr>
Date: Wed, 24 Jan 2024 20:10:01 +0100
Subject: [PATCH] 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 | 175 ++++++++++++++++++++++++++++++++++++-------------
1 file changed, 130 insertions(+), 45 deletions(-)
diff --git a/R/lauricella.R b/R/lauricella.R
index e56e4e0..8a5d3df 100644
--- a/R/lauricella.R
+++ b/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)
--
GitLab