From d1f584245f1164dbdbe3b00f6bc798d1865dc769 Mon Sep 17 00:00:00 2001
From: Pierre Santagostini <pierre.santagostini@agrocampus-ouest.fr>
Date: Wed, 17 Apr 2024 19:20:00 +0200
Subject: [PATCH] There was a problem if lambda[p]==1: Problem when computing
log(1-lambda[p])
---
R/kldcauchy.R | 48 +++++++++++++++++++++++++++++-------------------
1 file changed, 29 insertions(+), 19 deletions(-)
diff --git a/R/kldcauchy.R b/R/kldcauchy.R
index d2a40b2..7194a23 100644
--- a/R/kldcauchy.R
+++ b/R/kldcauchy.R
@@ -115,6 +115,16 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
if (lambda[p] <= 1) { # equations 109 & 102
+ if (p > 1 & lambda[p] == 1) {
+ lambda <- lambda[1:(p-1)]
+ p1 <- p-1
+ M <- expand.grid(rep(list(0:k), p1))
+ M <- M[-1, , drop = FALSE]
+ # Sum of the indices
+ Msum <- rowSums(M)
+ } else {
+ p1 <- p
+ }
# The first 5 elements of the sum
@@ -129,11 +139,11 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
# }
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 - lambda), length(Munique)), ncol = p, byrow = TRUE)
+ matpoch <- matrix(rep(Mpoch, p1), ncol = p1, byrow = FALSE)
+ matM <- matrix(rep(Munique, p1), ncol = p1, byrow = FALSE)
+ matlabmda <- matrix(rep(log(1 - lambda), length(Munique)), ncol = p1, byrow = TRUE)
lambdaM <- matlabmda*matM
- matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
+ matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
matcommun <- exp(matpoch + lambdaM - matfact)
matcommun <- as.data.frame(apply(matcommun, 2, Re))
@@ -177,15 +187,15 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
# 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) {
- for (i in 1:(p-1)) {
- indsupp <- combn(p, i)
+ M <- expand.grid(rep(list(k1), p1))
+ if (p1 > 1) {
+ for (i in 1:(p1-1)) {
+ indsupp <- combn(p1, i)
for (j in 1:ncol(indsupp)) {
jsupp <- indsupp[, j]
- Mlist <- vector("list", p)
+ Mlist <- vector("list", p1)
for (l in jsupp) Mlist[[l]] <- k1
- for (l in (1:p)[-jsupp]) Mlist[[l]] <- 0:k
+ for (l in (1:p1)[-jsupp]) Mlist[[l]] <- 0:k
M <- rbind(M, expand.grid(Mlist))
}
}
@@ -208,26 +218,26 @@ kldcauchy <- function(Sigma1, Sigma2, eps = 1e-06) {
# }
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 - lambda), length(Munique)), ncol = p, byrow = TRUE)
+ matpoch <- matrix(rep(Mpoch, p1), ncol = p1, byrow = FALSE)
+ matM <- matrix(rep(Munique, p1), ncol = p1, byrow = FALSE)
+ matlabmda <- matrix(rep(log(1 - lambda), length(Munique)), ncol = p1, byrow = TRUE)
lambdaM <- matlabmda*matM
- matfact <- matrix(rep(lfactorial(Munique), p), ncol = p, byrow = FALSE)
+ matfact <- matrix(rep(lfactorial(Munique), p1), ncol = p1, byrow = FALSE)
matcommun <- rbind(matcommun, matcommunsupp)
matcommunsupp <- exp(matpoch + lambdaM - matfact)
matcommunsupp <- as.data.frame(apply(matcommunsupp, 2, Re))
gridcommun <- expand.grid(matcommunsupp)
- if (p > 1) {
- for (i in 1:(p-1)) {
- indsupp <- combn(p, i)
+ if (p1 > 1) {
+ for (i in 1:(p1-1)) {
+ indsupp <- combn(p1, i)
for (j in 1:ncol(indsupp)) {
jsupp <- indsupp[, j]
- communlist <- vector("list", p)
+ communlist <- vector("list", p1)
names(communlist) <- names(gridcommun)
for (l in jsupp) communlist[[l]] <- matcommunsupp[, l]
- for (l in (1:p)[-jsupp]) communlist[[l]] <- matcommun[, l]
+ for (l in (1:p1)[-jsupp]) communlist[[l]] <- matcommun[, l]
gridcommun <- rbind(gridcommun, expand.grid(communlist))
}
}
--
GitLab