From 277382b6f40b9c93a771a7c9de4fa648699638ed Mon Sep 17 00:00:00 2001
From: Pierre Santagostini <pierre.santagostini@agrocampus-ouest.fr>
Date: Mon, 8 Jul 2024 11:57:01 +0200
Subject: [PATCH] Fixed: there was an error (line 445)
---
R/dlauricella.R | 140 ++++++++++++++++++++++++------------------------
1 file changed, 71 insertions(+), 69 deletions(-)
diff --git a/R/dlauricella.R b/R/dlauricella.R
index 834b172..80e9277 100644
--- a/R/dlauricella.R
+++ b/R/dlauricella.R
@@ -1,9 +1,7 @@
-#' @importFrom data.table rbindlist
-
dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
-
+
p <- length(lambda)
-
+
buildMlist <- function(j, isupp, k, k1, p = p) {
jsupp <- isupp[, j]
Ml <- rep(list(0:k), p)
@@ -11,7 +9,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
Ml[[j]] <- k1
return(expand.grid(Ml))
}
-
+
buildxlist <- function(j, isupp, x, xsupp, p = p) {
jsupp <- isupp[, j]
xl <- rep(list(x), p)
@@ -19,29 +17,29 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
xl[-jsupp] <- x[-jsupp]
return(expand.grid(xl))
}
-
+
prodlambda <- prod(lambda)
-
+
lambdanu <- lambda*nu1/nu2
prodlambdanu <- prod(lambdanu)
-
+
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), p))
M <- M[-1, , drop = FALSE]
-
+
Munique <- 0:k
-
+
# Sum of the indices
Msum <- rowSums(M)
# Msumunique <- 0:max(Msum)
-
+
kstep <- 5
-
+
if (lambdanu[p] <= 1) { # lambda[1] < ... < lambda[p] < 1
-
+
if (p > 1 & lambda[p] == 1) {
lambda <- lambda[1:(p-1)]
p1 <- p-1
@@ -52,25 +50,25 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
} else {
p1 <- p
}
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- matrix(rep(Mpoch, p1), ncol = p1, byrow = FALSE)
matM <- matrix(rep(Munique, p1), ncol = p1, byrow = FALSE)
matlabmda <- matrix(rep(log(1 - lambdanu[1:p1]), length(Munique)), ncol = p1, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
-
+
# matcommun <- exp(matpoch + lambdaM - matfact)
matcommun <- as.data.frame(matpoch + lambdaM - matfact)
# matcommun <- as.data.frame(apply(matcommun, 2, Re))
matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
gridcommun <- expand.grid(matcommun)
gridcommun <- gridcommun[-1, , drop = FALSE]
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- sum(
# commun * sapply(Msum, function(i) {
# pochhammer(1, i) / ( pochhammer((nu1 + p)/2, i) * i )
@@ -81,7 +79,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
})
)))
-
+
# Next elements of the sum, until the expected precision
k1 <- 1:k
derive <- 0
@@ -90,7 +88,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
k <- k1[length(k1)]
k1 <- k + (1:kstep)
derive <- derive + d
-
+
# M: data.frame of the indices for the nested sums
# M <- expand.grid(rep(list(k1), p1))
# if (p1 > 1) {
@@ -117,10 +115,10 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
M <- data.frame(rbindlist(Mlist))
}
-
+
# Sum of the indices
Msum <- rowSums(M)
-
+
Munique <- (max(Munique)+1):max(M)
# Msumunique <- (max(Msumunique)+1):max(Msum)
@@ -130,13 +128,13 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
matlabmda <- matrix(rep(log(1 - lambdanu[1:p1]), length(Munique)), ncol = p1, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
-
+
matcommun <- rbind(matcommun, matcommunsupp)
# matcommunsupp <- exp(matpoch + lambdaM - matfact)
matcommunsupp <- as.data.frame(matpoch + lambdaM - matfact)
# matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
communlist <- list(expand.grid(matcommunsupp))
if (p1 == 1) {
gridcommun <- communlist[[1]]
@@ -152,7 +150,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- sum(
# commun * sapply(Msum, function(i) {
# pochhammer(1, i) / ( pochhammer((nu1 + p)/2, i) * i )
@@ -163,38 +161,38 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
lnpochhammer(1, i) - ( lnpochhammer((nu1 + p)/2, i) + log(i) )
})
)))
-
+
}
-
+
derive <- as.numeric(-0.5 * log(prodlambda) - (nu2 + p)/2 * derive)
-
+
} else if (lambdanu[1] > 1) { # 1 < lambda[1] < ... < lambda[p]
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- matrix(rep(Mpoch, p), ncol = p, byrow = FALSE)
matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
matlabmda <- matrix(rep(log(1 - 1/lambdanu), length(Munique)), ncol = p, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
-
+
# matcommun <- exp(matpoch + lambdaM - matfact)
matcommun <- as.data.frame(matpoch + lambdaM - matfact)
# matcommun <- as.data.frame(apply(matcommun, 2, Re))
matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
gridcommun <- expand.grid(matcommun)
gridcommun <- gridcommun[-1, , drop = FALSE]
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
d <- 0
for (i in 1:length(Msum)) {
A <- sum(1/(0:(Msum[i]-1) + (nu1+p)/2))
d <- d - exp(commun[i]) * A
}
d <- Re(d)
-
+
# Next elements of the sum, until the expected precision
k1 <- 1:k
derive <- 0
@@ -205,7 +203,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
k1 <- k + (1:kstep)
derive <- derive + d
# vd <- c(vd, d); vderive <- c(vderive, derive)
-
+
# # M: data.frame of the indices for the nested sums
# M <- as.data.frame(matrix(nrow = 0, ncol = p))
# if (p > 1) {
@@ -219,7 +217,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# }
# }
# M <- rbind( M, expand.grid(rep(list(k1), p)) )
-
+
# M: data.frame of the indices for the nested sums
# M <- expand.grid(rep(list(k1), p))
# if (p > 1) {
@@ -246,12 +244,12 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
M <- data.frame(rbindlist(Mlist))
}
-
+
Msum <- rowSums(M)
-
+
Munique <- (max(Munique) + 1):max(M)
# Msumunique <- (max(Msumunique) + 1):max(Msum)
-
+
# d <- 0
# for (i in 1:length(Msum)) {
# commun <- prod(
@@ -262,20 +260,20 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
# d <- d - commun * A # / pochhammer((1 + p)/2, Msum[i])
# }
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- matrix(rep(Mpoch, p), ncol = p, byrow = FALSE)
matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
matlabmda <- matrix(rep(log(1 - 1/lambdanu), length(Munique)), ncol = p, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
-
+
matcommun <- rbind(matcommun, matcommunsupp)
# matcommunsupp <- exp(matpoch + lambdaM - matfact)
matcommunsupp <- as.data.frame(matpoch + lambdaM - matfact)
# matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
communlist <- list(expand.grid(matcommunsupp))
if (p == 1) {
gridcommun <- communlist[[1]]
@@ -291,40 +289,42 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- 0
# for (i in 1:length(Msum)) {
# A <- sum(1/(0:(Msum[i]-1) + (1+p)/2))
# d <- d - commun[i] * A
# }
A <- sapply(Msum, function (i) sum(1/(0:(i-1) + (nu1+p)/2)))
-
+
d <- Re(-sum(exp(commun)*A))
-
+
}
-
+
derive <- as.numeric(-0.5 * log(prodlambda) - (nu2 + p)/2 * prod(1/sqrt(lambdanu)) * derive)
-
+
} else { # lambda[1] < ... < 1 < ... < lambda[p]
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), 0)
+ # Mpochp <- sapply(Munique, function(j) lnpochhammer(nu1/2, j))
+ # matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), Mpochp)
matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
matlabmda <- matrix(rep(log(1 - c(lambdanu[-p], 1)/lambdanu[p]), length(Munique)), ncol = p, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
-
+
# matcommun <- exp(matpoch + lambdaM - matfact)
matcommun <- as.data.frame(matpoch + lambdaM - matfact)
# matcommun <- as.data.frame(apply(matcommun, 2, Re))
matcommunsupp <- as.data.frame(matrix(nrow = 0, ncol = ncol(matcommun)))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
gridcommun <- expand.grid(matcommun)
gridcommun <- gridcommun[-1, , drop = FALSE]
# commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- sum(
# commun * sapply(1:length(Msum), function(i) {
# pochhammer(0.5, M[i, p]) * pochhammer(1, Msum[i]) / ( pochhammer((nu1 + p)/2, Msum[i]) * Msum[i] )
@@ -332,11 +332,11 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# )
d <- Re(sum(exp(
commun + sapply(1:length(Msum), function(i) {
- lnpochhammer(0.5, M[i, p]) + lnpochhammer(1, Msum[i]) -
+ lnpochhammer(nu1/2, M[i, p]) + lnpochhammer(1, Msum[i]) -
( lnpochhammer((nu1 + p)/2, Msum[i]) + log(Msum[i]) )
})
)))
-
+
# Next elements of the sum, until the expected precision
k1 <- 1:k
derive <- 0
@@ -345,7 +345,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
k <- k1[length(k1)]
k1 <- k + (1:kstep)
derive <- derive + d
-
+
# # M: data.frame of the indices for the nested sums
# M <- as.data.frame(matrix(nrow = 0, ncol = p))
# if (p > 1) {
@@ -359,7 +359,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# }
# }
# M <- rbind( M, expand.grid(rep(list(k1), p)) )
-
+
# M <- expand.grid(rep(list(k1), p))
# if (p > 1) {
# for (i in 1:(p-1)) {
@@ -373,7 +373,7 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# }
# }
# }
-
+
# M: data.frame of the indices for the nested sums
Mlist <- list(expand.grid(rep(list(k1), p)))
if (p == 1) {
@@ -387,12 +387,12 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
}
M <- data.frame(rbindlist(Mlist))
}
-
+
Msum <- rowSums(M)
-
+
Munique <- (max(Munique)+1):max(M)
# Msumunique <- (max(Msumunique)+1):max(Msum)
-
+
# d <- 0
# for (i in 1:length(Msum)) {
# commun <- prod(
@@ -403,20 +403,22 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# commun <- commun*(1 - 1/lambda[p])^M[i, p]/factorial(M[i, p])
# d <- d + commun * pochhammer(0.5, M[i, p])*pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
# }
-
+
Mpoch <- sapply(Munique, function(j) lnpochhammer(0.5, j))
matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), 0)
+ # Mpochp <- sapply(Munique, function(j) lnpochhammer(nu1/2, j))
+ # matpoch <- cbind(matrix(rep(Mpoch, p-1), ncol = p-1, byrow = FALSE), Mpochp)
matM <- matrix(rep(Munique, p), ncol = p, byrow = FALSE)
matlabmda <- matrix(rep(log(1 - c(lambdanu[-p], 1)/lambdanu[p]), length(Munique)), ncol = p, byrow = TRUE)
lambdaM <- matlabmda*matM
matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
-
+
matcommun <- rbind(matcommun, matcommunsupp)
# matcommunsupp <- exp(matpoch + lambdaM - matfact)
matcommunsupp <- as.data.frame(matpoch + lambdaM - matfact)
# matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
colnames(matcommunsupp) <- colnames(matcommun)
-
+
communlist <- list(expand.grid(matcommunsupp))
if (p == 1) {
gridcommun <- communlist[[1]]
@@ -430,9 +432,9 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
names(communlist[[1]]) <- names(communlist[[2]])
gridcommun <- data.frame(rbindlist(communlist))
}
- commun <- apply(gridcommun, 1, prod)
+ # commun <- apply(gridcommun, 1, prod)
commun <- rowSums(gridcommun)
-
+
# d <- sum(
# commun * sapply(1:length(Msum), function(i) {
# pochhammer(0.5, M[i, p]) * pochhammer(1, Msum[i]) / ( pochhammer((nu1 + p)/2, Msum[i]) * Msum[i] )
@@ -440,19 +442,19 @@ dlauricella <- function(nu1, nu2, lambda, eps = 1e-6) {
# )
d <- Re(sum(exp(
commun + sapply(1:length(Msum), function(i) {
- lnpochhammer(0.5, M[i, p]) + lnpochhammer(1, Msum[i]) -
+ lnpochhammer(nu1/2, M[i, p]) + lnpochhammer(1, Msum[i]) -
( lnpochhammer((nu1 + p)/2, Msum[i]) + log(Msum[i]) )
})
)))
-
+
}
-
+
derive <- as.numeric(-0.5 * log(prodlambda) + (nu2 + p)/2 * (log(lambdanu[p]) - derive))
-
+
}
-
+
attr(derive, "epsilon") <- eps
attr(derive, "k") <- k
-
+
return(derive)
}
--
GitLab