From 9a48960a7605a5a7cf9fdf3ed69d63cd97c54c84 Mon Sep 17 00:00:00 2001
From: Pierre Santagostini <pierre.santagostini@agrocampus-ouest.fr>
Date: Wed, 13 Mar 2024 16:57:59 +0100
Subject: [PATCH] Fixed: the calculation of the data frame containing the
indices for the nested sums was wrong for dimension p>3.
---
R/kldcauchy.R | 173 +++++++++++++++++++++++++++++++++++++++-----------
1 file changed, 135 insertions(+), 38 deletions(-)
diff --git a/R/kldcauchy.R b/R/kldcauchy.R
index 3654089..48fd5ab 100644
--- a/R/kldcauchy.R
+++ b/R/kldcauchy.R
@@ -63,6 +63,7 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
#' kldcauchy(Sigma2, Sigma1)
#' }
#'
+ #' @importFrom utils combn
#' @export
# Sigma1 and Sigma2 must be matrices
@@ -129,19 +130,35 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
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) {
+ # for (i in 1:(p-1)) {
+ # Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
+ # M <- rbind( M, expand.grid(Mlist) )
+ # for (j in 1:(p-1)) {
+ # Mlist <- Mlist[c(p, 1:(p-1))]
+ # M <- rbind(M, expand.grid(Mlist))
+ # }
+ # }
+ # }
+ # M <- rbind( M, expand.grid(rep(list(k1), p)) )
+
# M: data.frame of the indices for the nested sums
- M <- as.data.frame(matrix(nrow = 0, ncol = p))
+ M <- expand.grid(rep(list(k1), p))
if (p > 1) {
for (i in 1:(p-1)) {
- Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
- M <- rbind( M, expand.grid(Mlist) )
- for (j in 1:(p-1)) {
- Mlist <- Mlist[c(p, 1:(p-1))]
+ indsupp <- combn(p, i)
+ for (j in 1:ncol(indsupp)) {
+ jsupp <- indsupp[, j]
+ Mlist <- vector("list", p)
+ for (l in jsupp) Mlist[[l]] <- k1
+ for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
M <- rbind(M, expand.grid(Mlist))
}
}
}
- M <- rbind( M, expand.grid(rep(list(k1), p)) )
+
Msum <- apply(M, 1, sum)
d <- 0
@@ -166,22 +183,38 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
k <- k1
k1 <- k + 1
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) {
+ # for (i in 1:(p-1)) {
+ # Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
+ # M <- rbind( M, expand.grid(Mlist) )
+ # for (j in 1:(p-1)) {
+ # Mlist <- Mlist[c(p, 1:(p-1))]
+ # M <- rbind(M, expand.grid(Mlist))
+ # }
+ # }
+ # }
+ # M <- rbind( M, rep(k1, p) )
# M: data.frame of the indices for the nested sums
- M <- as.data.frame(matrix(nrow = 0, ncol = p))
+ M <- expand.grid(rep(list(k1), p))
if (p > 1) {
for (i in 1:(p-1)) {
- Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
- M <- rbind( M, expand.grid(Mlist) )
- for (j in 1:(p-1)) {
- Mlist <- Mlist[c(p, 1:(p-1))]
+ indsupp <- combn(p, i)
+ for (j in 1:ncol(indsupp)) {
+ jsupp <- indsupp[, j]
+ Mlist <- vector("list", p)
+ for (l in jsupp) Mlist[[l]] <- k1
+ for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
M <- rbind(M, expand.grid(Mlist))
}
}
}
- M <- rbind( M, rep(k1, p) )
- Msum <- apply(M, 1, sum)
+ Msum <- apply(M, 1, sum)
+
d <- 0
for (i in 1:length(Msum)) {
commun <- prod(
@@ -191,7 +224,7 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
)
d <- d + commun * pochhammer(1, Msum[i]) / ( pochhammer((1 + p)/2, Msum[i]) * Msum[i] )
}
-
+
}
}
@@ -222,19 +255,35 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
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) {
+ # for (i in 1:(p-1)) {
+ # Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
+ # M <- rbind( M, expand.grid(Mlist) )
+ # for (j in 1:(p-1)) {
+ # Mlist <- Mlist[c(p, 1:(p-1))]
+ # M <- rbind(M, expand.grid(Mlist))
+ # }
+ # }
+ # }
+ # M <- rbind( M, expand.grid(rep(list(k1), p)) )
+
# M: data.frame of the indices for the nested sums
- M <- as.data.frame(matrix(nrow = 0, ncol = p))
+ M <- expand.grid(rep(list(k1), p))
if (p > 1) {
for (i in 1:(p-1)) {
- Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
- M <- rbind( M, expand.grid(Mlist) )
- for (j in 1:(p-1)) {
- Mlist <- Mlist[c(p, 1:(p-1))]
+ indsupp <- combn(p, i)
+ for (j in 1:ncol(indsupp)) {
+ jsupp <- indsupp[, j]
+ Mlist <- vector("list", p)
+ for (l in jsupp) Mlist[[l]] <- k1
+ for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
M <- rbind(M, expand.grid(Mlist))
}
}
}
- M <- rbind( M, expand.grid(rep(list(k1), p)) )
+
Msum <- apply(M, 1, sum)
d <- 0
@@ -261,19 +310,35 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
k1 <- k + 1
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) {
+ # for (i in 1:(p-1)) {
+ # Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
+ # M <- rbind( M, expand.grid(Mlist) )
+ # for (j in 1:(p-1)) {
+ # Mlist <- Mlist[c(p, 1:(p-1))]
+ # M <- rbind(M, expand.grid(Mlist))
+ # }
+ # }
+ # }
+ # M <- rbind( M, rep(k1, p) )
+
# M: data.frame of the indices for the nested sums
- M <- as.data.frame(matrix(nrow = 0, ncol = p))
+ M <- expand.grid(rep(list(k1), p))
if (p > 1) {
for (i in 1:(p-1)) {
- Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
- M <- rbind( M, expand.grid(Mlist) )
- for (j in 1:(p-1)) {
- Mlist <- Mlist[c(p, 1:(p-1))]
+ indsupp <- combn(p, i)
+ for (j in 1:ncol(indsupp)) {
+ jsupp <- indsupp[, j]
+ Mlist <- vector("list", p)
+ for (l in jsupp) Mlist[[l]] <- k1
+ for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
M <- rbind(M, expand.grid(Mlist))
}
}
}
- M <- rbind( M, rep(k1, p) )
+
Msum <- apply(M, 1, sum)
d <- 0
@@ -315,19 +380,35 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
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) {
+ # for (i in 1:(p-1)) {
+ # Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
+ # M <- rbind( M, expand.grid(Mlist) )
+ # for (j in 1:(p-1)) {
+ # Mlist <- Mlist[c(p, 1:(p-1))]
+ # M <- rbind(M, expand.grid(Mlist))
+ # }
+ # }
+ # }
+ # M <- rbind( M, expand.grid(rep(list(k1), p)) )
+
# M: data.frame of the indices for the nested sums
- M <- as.data.frame(matrix(nrow = 0, ncol = p))
+ M <- expand.grid(rep(list(k1), p))
if (p > 1) {
for (i in 1:(p-1)) {
- Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
- M <- rbind( M, expand.grid(Mlist) )
- for (j in 1:(p-1)) {
- Mlist <- Mlist[c(p, 1:(p-1))]
+ indsupp <- combn(p, i)
+ for (j in 1:ncol(indsupp)) {
+ jsupp <- indsupp[, j]
+ Mlist <- vector("list", p)
+ for (l in jsupp) Mlist[[l]] <- k1
+ for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
M <- rbind(M, expand.grid(Mlist))
}
}
}
- M <- rbind( M, expand.grid(rep(list(k1), p)) )
+
Msum <- apply(M, 1, sum)
d <- 0
@@ -354,19 +435,35 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
k1 <- k + 1
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) {
+ # for (i in 1:(p-1)) {
+ # Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
+ # M <- rbind( M, expand.grid(Mlist) )
+ # for (j in 1:(p-1)) {
+ # Mlist <- Mlist[c(p, 1:(p-1))]
+ # M <- rbind(M, expand.grid(Mlist))
+ # }
+ # }
+ # }
+ # M <- rbind( M, rep(k1, p) )
+
# M: data.frame of the indices for the nested sums
- M <- as.data.frame(matrix(nrow = 0, ncol = p))
+ M <- expand.grid(rep(list(k1), p))
if (p > 1) {
for (i in 1:(p-1)) {
- Mlist <- c( rep(list(0:k), p-i), rep(list(k1), i) )
- M <- rbind( M, expand.grid(Mlist) )
- for (j in 1:(p-1)) {
- Mlist <- Mlist[c(p, 1:(p-1))]
+ indsupp <- combn(p, i)
+ for (j in 1:ncol(indsupp)) {
+ jsupp <- indsupp[, j]
+ Mlist <- vector("list", p)
+ for (l in jsupp) Mlist[[l]] <- k1
+ for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
M <- rbind(M, expand.grid(Mlist))
}
}
}
- M <- rbind( M, rep(k1, p) )
+
Msum <- apply(M, 1, sum)
d <- 0
--
GitLab