Hey y'all, I think I may have found a potential off by one indexing bug in fit_aitchison.R:
|
## Next: take average of EM samples past burn. Included init, so have EMiter+1 total |
|
b0_EM <- colMeans(b0_list[(EMburn + 1):(EMiter + 1), ]) |
|
b_list_reduced <- b_list[, , (EMburn + 1):(EMiter + 1)] |
|
sigma_em <- apply(sigma_list[, , (EMburn + 1):(EMiter + 1)], c(1, 2), mean) |
I believe these three should all be (EMburn + 2):(EMiter + 1).
E.g., for EMburn <- 3 and EMiter <- 6, the first item is the init, next three are EM iters in the burn in stage, next three are EM iters past the burn in stage. But as written these would all be indexing (3+1):(6+1) -> 4:7 -> c(4, 5, 6, 7). That is the last EM iter in the burn stage, and the three after the burn stage. Doing (EMburn + 2):(EMiter + 1) instead gives you the 5, 6, and 7, which are only the EM iters past the burn in which is what it should be if I'm not mistaken.
Here's where the things are initialized for reference:
In all three cases the first item is the "init". Here is b_list and b0_list being initialized before the EM loop starts:
|
b0_list <- matrix(0, nrow = EMiter + 1, ncol = length(b0)) |
|
b0_list[1,] <- b0 |
|
if (!no_covariates) { |
|
if (is.matrix(b)) { |
|
b_list <- array(0, dim = c(dim(b), EMiter + 1)) |
|
b_list[,, 1] <- b |
|
} else { |
|
b_list <- matrix(0, nrow = length(b), ncol = EMiter + 1) |
|
b_list[, 1] <- b |
|
} |
|
} |
|
|
And here is sigma_list:
|
sigma_list <- array(0, dim = c(dim(sigma), EMiter + 1)) |
|
sigma_list[, , 1] <- sigma |
Hey y'all, I think I may have found a potential off by one indexing bug in fit_aitchison.R:
DivNet/R/fit_aitchison.R
Lines 176 to 177 in 63908cf
DivNet/R/fit_aitchison.R
Line 185 in 63908cf
DivNet/R/fit_aitchison.R
Line 199 in 63908cf
I believe these three should all be
(EMburn + 2):(EMiter + 1).E.g., for
EMburn <- 3andEMiter <- 6, the first item is the init, next three are EM iters in the burn in stage, next three are EM iters past the burn in stage. But as written these would all be indexing(3+1):(6+1)->4:7->c(4, 5, 6, 7). That is the last EM iter in the burn stage, and the three after the burn stage. Doing(EMburn + 2):(EMiter + 1)instead gives you the 5, 6, and 7, which are only the EM iters past the burn in which is what it should be if I'm not mistaken.Here's where the things are initialized for reference:
In all three cases the first item is the "init". Here is
b_listandb0_listbeing initialized before the EM loop starts:DivNet/R/fit_aitchison.R
Lines 103 to 114 in 63908cf
And here is
sigma_list:DivNet/R/fit_aitchison.R
Lines 117 to 118 in 63908cf