Skip to content
Snippets Groups Projects
Commit ec59dc98 authored by SANTAGOSTINI Pierre's avatar SANTAGOSTINI Pierre
Browse files

Ref corrected

parent 788d0fc2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
R/lauricella.R
+ 36
− 36
View file @ ec59dc98
......@@ -17,52 +17,52 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
#' @details If \eqn{n} is the length of the \eqn{b} and \code{x} vectors,
#' the Lauricella \eqn{D}-hypergeometric Function function is given by:
#' \deqn{\displaystyle{F_D^{(n)}\left(a, b_1, ..., b_n, g; x_1, ..., x_n\right) = \sum_{m_1 \geq 0} ... \sum_{m_n \geq 0}{ \frac{ (a)_{m_1+...+m_n}(b_1)_{m_1} ... (b_n)_{m_n} }{ (g)_{m_1+...+m_n} } \frac{x_1^{m_1}}{m_1!} ... \frac{x_n^{m_n}}{m_n!} } }}
#'
#'
#' where \eqn{(x)_p} is the Pochhammer symbol (see \code{\link{pochhammer}}).
#'
#'
#' If \eqn{|x_i| < 1, i = 1, \dots, n}, this sum converges.
#' Otherwise there is an error.
#'
#'
#' The \code{eps} argument gives the required precision for its computation.
#' It is the \code{attr(, "epsilon")} attribute of the returned value.
#'
#'
#' Sometimes, the convergence is too slow and the required precision cannot be reached.
#' 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, Kullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions.
#' IEEE Signal Processing Letters, vol. 26 no. 7, July 2019.
#' @references N. Bouhlel and D. Rousseau, Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions.
#' IEEE Signal Processing Letters Processing Letters, vol. 26 no. 7, July 2019.
#' \doi{10.1109/LSP.2019.2915000}
#' @importFrom utils combn
#' @export
# Number of variables
n <- length(x)
# Do x and b have the same length?
if (length(b) != n)
stop("x and b must have the same length")
# Condition for the convergence: are all abs(x) < 1 ?
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 <- expand.grid(rep(list(0:k), n))
# Sum of the indices
Msum <- rowSums(M)
Munique <- 0:max(M)
Msumunique <- 0:max(Msum)
# 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 ))
# lnfact <- apply(M, 1, function(Mi) sum(sapply(Mi, lnfactorial)))
......@@ -78,7 +78,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# }
gridxfact <- expand.grid(xfact)
prodxfact <- apply(gridxfact, 1, prod)
# Logarithm of the product m_1! * ... * m_n! for m_1 = 0...k, ..., m_n = 0...k
# i.e. \sum_{i=0}^n{\log{m_i!}}
# lnfact <- as.data.frame(
......@@ -95,7 +95,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
gridlnfact <- expand.grid(lnfact)
# Sum of the logarithms
sumlnfact <- rowSums(gridlnfact)
# 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))
......@@ -122,12 +122,12 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
gridlnbpoch <- expand.grid(as.data.frame(lnbpoch))
# Sum of the logarithms
sumlnpochb <- rowSums(gridlnbpoch)
# 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
# res1 <- lnapoch + lnbpoch - lngpoch - lnfact
# res2 <- sum(xfact * exp(res1))
res1 <- lnapoch[Msum+1] + sumlnpochb - lngpoch[Msum+1] - sumlnfact
......@@ -135,12 +135,12 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# 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))
......@@ -153,16 +153,16 @@ 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))
# }
while (abs(res2) > eps/10 & !is.nan(res2)) {
epsret <- signif(abs(res2), 1)*10
k <- k1[length(k1)]
k1 <- k + (1:kstep)
result <- result + res2
# # M: data.frame of the indices for the nested sums
# M <- expand.grid(rep(list(k1), n))
# if (n > 1) {
......@@ -176,7 +176,7 @@ 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) {
......@@ -192,13 +192,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} * ... * 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)))
......@@ -221,7 +221,7 @@ 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(
......@@ -240,13 +240,13 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
}
}
sumlnfact <- rowSums(gridlnfact)
# 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
lnapoch <- c(lnapoch, lnapochsupp)
# lnbpoch <- numeric(nrow(M))
# for (i in 1:nrow(M)) {
# # Product pochhammer(b_1,m_1) * ...* pochhammer(b_n, m_n)
......@@ -274,30 +274,30 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
}
# Sum of the logarithms
sumlnpochb <- rowSums(gridlnbpoch)
# 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
lngpoch <- c(lngpoch, lngpochsupp)
# 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
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)
}
result <- Re(result)
attr(result, "epsilon") <- eps
attr(result, "k") <- k
# Returns the result of the nested sums
return(result)
}
man/lauricella.Rd
+ 3
− 2
View file @ ec59dc98
......@@ -41,8 +41,9 @@ Sometimes, the convergence is too slow and the required precision cannot be reac
If this happens, the \code{attr(, "epsilon")} attribute is the precision that was really reached.
}
\references{
N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions, IEEE Signal Processing Letters.
\doi{10.1109/LSP.2023.3324594}
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}
}
\author{
Pierre Santagostini, Nizar Bouhlel
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment