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

Internal functions removed. Calculation of the sums: with rowSums(M) and not apply(M, 1, sum)

parent e821feae
No related branches found
No related tags found
No related merge requests found
......@@ -35,29 +35,7 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
#' \doi{10.1109/LSP.2019.2915000}
#' @importFrom utils combn
#' @export
# # Internal functions
# lnfactorial <- function(n) {
# # Logarithm of the factorial of a positive integer
# if (n == 0) {
# return(0)
# }
#
# if (n > 0) {
# return(sum(log(1:n)))
# }
# }
# lnpochhammer <- function(x, n) {
# # Logarithm of the Pochhammer Symbol
# if (n == 0) {
# return(0)
# }
# if (n > 0) {
# y <- x + (1:n) - 1
# return(sum(log(abs(y)) + (y < 0)*pi*1i))
# }
# }
# Number of variables
n <- length(x)
......@@ -79,12 +57,13 @@ lauricella <- function(a, b, g, x, eps = 1e-06) {
# # (i.e. all arrangements of n elements from {0:k})
M <- expand.grid(rep(list(0:k), n))
Msum <- apply(M, 1, sum)
# Sum of the indices
Msum <- rowSums(M)
Munique <- 0:max(M)
Msumunique <- 0:max(Msum)
# Product x^{m_1}/m_1! * ... * x^{m_n}/m_n! for m_1 = 0...k, ..., m_n = 0...k
# 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)))
xfact <- as.data.frame(
......@@ -100,6 +79,8 @@ 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(
matrix(sapply(Munique, lnfactorial), nrow = length(Munique),
ncol = n, byrow = FALSE, dimnames = list(Munique, 1:n)))
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment