From a59d00dc1826f6b06200d381a368cb75def953d9 Mon Sep 17 00:00:00 2001
From: Pierre Santagostini <pierre.santagostini@agrocampus-ouest.fr>
Date: Wed, 17 Jan 2024 17:56:42 +0100
Subject: [PATCH] lauricella now uses lnfactorial and lnpochhammer (instead of
factorial and pochhammer). Typing error corrected in ref.
---
R/lauricella.R | 101 +++++++++++--------------------------------------
1 file changed, 22 insertions(+), 79 deletions(-)
diff --git a/R/lauricella.R b/R/lauricella.R
index a218f1e..5dbe91f 100644
--- a/R/lauricella.R
+++ b/R/lauricella.R
@@ -30,7 +30,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
#' If this happens, the \code{attr(, "epsilon")} attribute is the precision that was really reached.
#'
#' @author Pierre Santagostini, Nizar Bouhlel
- #' @references N. Bouhlel, A. Dziri, Jullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions.
+ #' @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}
#' @export
@@ -55,35 +55,37 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
Msum <- apply(M, 1, sum)
# 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 / factorial(Mi) ))
+ xfact <- apply(M, 1, function(Mi) prod( x^Mi ))
+ lnfact <- apply(M, 1, function(Mi) sum(sapply(Mi, lnfactorial)))
# pochhammer(a, m_1+...+m_n) for m_1 = 0...k, ..., m_p = 0...k
- apoch <- sapply(Msum, function(j) pochhammer(a, j))
+ lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
- bpoch <- numeric(nrow(M))
+ lnbpoch <- numeric(nrow(M))
for (i in 1:nrow(M)) {
# Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
- bpoch[i] <- prod( sapply(1:n, function(j) pochhammer(b[j], M[i, j])) )
+ lnbpoch[i] <- sum( sapply(1:n, function(j) lnpochhammer(b[j], M[i, j])) )
}
# pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ..., m_p = 0...k
- gpoch <- sapply(Msum, function(j) pochhammer(g, j))
+ lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
- res1 <- sum( xfact * apoch * bpoch / gpoch )
+ res1 <- lnapoch + lnbpoch - lngpoch - lnfact
+ res2 <- sum(xfact * exp(res1))
kstep <- 5
k1 <- 1:k
result <- 0
- while (abs(res1) > eps/10 & !is.nan(res1)) {
+ while (abs(res2) > eps/10 & !is.nan(res2)) {
- epsret <- signif(abs(res1), 1)*10
+ epsret <- signif(abs(res2), 1)*10
k <- k1[length(k1)]
k1 <- k + (1:kstep)
- result <- result + res1
+ result <- result + res2
# M: data.frame of the indices for the nested sums
M <- as.data.frame(matrix(nrow = 0, ncol = n))
@@ -102,86 +104,27 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
Msum <- apply(M, 1, sum)
# 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 / factorial(Mi) ))
+ xfact <- apply(M, 1, function(Mi) prod( x^Mi ))
+ lnfact <- apply(M, 1, function(Mi) sum(sapply(Mi, lnfactorial)))
# pochhammer(a, m_1+...+m_n) for m_1, ..., m_n given by the rows of M
- apoch <- sapply(Msum, function(j) pochhammer(a, j))
+ lnapoch <- sapply(Msum, function(j) lnpochhammer(a, j))
- bpoch <- numeric(nrow(M))
+ lnbpoch <- numeric(nrow(M))
for (i in 1:nrow(M)) {
# Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
- bpoch[i] <- prod( sapply(1:n, function(j) pochhammer(b[j], M[i, j])) )
+ lnbpoch[i] <- sum( sapply(1:n, function(j) lnpochhammer(b[j], M[i, j])) )
}
# pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ..., m_p = 0...k
- gpoch <- sapply(Msum, function(j) pochhammer(g, j))
+ lngpoch <- sapply(Msum, function(j) lnpochhammer(g, j))
- res1 <- sum( xfact * apoch * bpoch / gpoch )
-
- }
-
- if (is.nan(res1)) {
-
- res1 <- 0
- k1 <- k
-
- while(!is.nan(res1)) {
-
- if (res1 > 0)
- epsret <- abs(res1)*10
-
- k <- k1
-
- k1 <- k + 1
- result <- result + res1
-
- # M: data.frame of the indices for the nested sums
- M <- as.data.frame(matrix(nrow = 0, ncol = 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))
- }
- }
- }
- M <- rbind( M, rep(k1, n) )
-
- Msum <- apply(M, 1, sum)
-
- # Product x^{m_1} * ... * x^{m_n} for m_1 = 0...k, ..., m_n = 0...k
- xfact <- apply(M, 1, function(Mi) prod( x^Mi / factorial(Mi) ))
-
- # pochhammer(a, m_1+...+m_n) for m_1 = 0...k, ..., m_p = 0...k
- apoch <- sapply(Msum, function(j) pochhammer(a, j))
-
- bpoch <- numeric(nrow(M))
- for (i in 1:nrow(M)) {
- # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
- bpoch[i] <- prod( sapply(1:n, function(j) pochhammer(b[j], M[i, j])) )
- }
-
- # pochhammer(c, m_1+...+m_n) for m_1 = 0...k, ..., m_p = 0...k
- gpoch <- sapply(Msum, function(j) pochhammer(g, j))
-
- res1 <- sum( xfact * apoch * bpoch / gpoch )
-
- }
- }
-
- if (is.nan(res1)) {
- epsret <- signif(epsret, 1)
- warning("\nCannot reach the precision ", eps, " due to NaN\n",
- "Number of iteration: ", k, "\n",
- "Precision reached:", epsret)
- attr(result, "epsilon") <- epsret
- } else {
- # result <- result + res1
- attr(result, "epsilon") <- eps
+ res1 <- lnapoch + lnbpoch - lngpoch - lnfact
+ res2 <- sum(xfact * exp(res1))
}
+ result <- Re(result)
+ attr(result, "epsilon") <- eps
attr(result, "k") <- k
# Returns the result of the nested sums
--
GitLab