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

Commentaries added/updated

parent c1f4f20d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 21 additions and 10 deletions
R/lauricella.R
+ 21
− 10
View file @ d58981ba
......@@ -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)
......
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