Closed khaledbouguila closed 3 weeks ago
Thanks for the report.
I am on it, and it seems that there is a sign error in the objective of the VE-step
Since we are optimizing something related to the negative log-likelihood, the trace term should be
+ 0.5 * trace(Omega * nSigma);
right?
If anybody e.g. (@mahendra-mariadassou) can double check/confirm this...
Ok, now the following piece of code
library("PLNmodels")
# choose data
n <- 100
p <- 2
counts <- matrix(rpois(n*p, c(5,11)), n, p)
covariates <- matrix(1, n, 1)
# Fit PLN on data
data <- prepare_data(counts, covariates)
model <- PLN(Abundance ~ 1 + offset(log(Offset)), data = data)
# take training data
# get parameters for the VE-step
new_n <- n
ind <- 1:new_n
new_data <- data[ind, , drop=FALSE]
new_responses <- new_data$Abundance
new_covariates <- matrix(new_data$V1, ncol=1)
new_offsets <- matrix(rep(log(new_data$Offset), p), ncol = p)
new_weights <- rep(1, new_n)
B <- model$model_par$B
Omega <- model$model_par$Omega
M_init <- model$var_par$M[ind, , drop=FALSE]
S_init <- model$var_par$S[ind, , drop=FALSE]
M_init <- matrix(0 , new_n, p)
S_init <- matrix(0.1, new_n, p)
args <- list(data = list(Y = new_responses, X = new_covariates, O = new_offsets, w = new_weights),
## Initialize the variational parameters with the new dimension of the data
params = list(M = M_init, S = S_init),
B = as.matrix(B),
Omega = as.matrix(Omega) ,
config = PLN_param()$config_optim)
VE <- do.call(PLNmodels:::nlopt_optimize_vestep, args)
mse <- function(a, b) sum((a-b)^2)
print(mse(VE$S**2, S_init**2))
print(mse(VE$M, M_init))
print(mse(VE$Ji, model$loglik_vec[ind]))
is returning, as expected
[1] 0.01935975
[1] 1.701245e-05
[1] 4.389937e-10
Two problems related to the VE-Step optimization:
self$var_par$M
andself$var_par$S
, which are the optimum parameters found by the VEM optimization, the VE-Step optimization changes their values and finds new optimums $M$ and $S$ (the other parameters $B$ and $\Omega$ are fixed).This is an example on trichoptera data that illustrates the problems: