Open
Description
Describe the bug
Running loo(moment_match = TRUE) on a fit object that has a parameter of type cholesky_factor_cov[K] yields the error "cholesky_factor_free: y is not lower triangular; y[1,3]=x (found before start of program)" where x may vary from one run to the next; in one case it was 6.65e271, in another case it was 1.26e-312.
To Reproduce
create_model <- function() {
code <- '
data {
int<lower=1> K;
int<lower=1> N;
array[N] vector[K] X;
real<lower = K - 1 + 1e-3> nu;
cov_matrix[K] scale_matrix;
}
transformed data {
cholesky_factor_cov[K] L_scale_matrix = cholesky_decompose(scale_matrix);
array[N] vector[K] zero = rep_array(rep_vector(0.0, K), N);
}
parameters {
cholesky_factor_cov[K] LX;
}
model {
X ~ multi_normal_cholesky(zero, LX);
LX ~ inv_wishart_cholesky(nu, L_scale_matrix);
}
generated quantities {
vector[N] log_lik;
real lp;
for (i in 1:N)
log_lik[i] = multi_normal_cholesky_lpdf(X[i] | rep_vector(0, K), LX);
lp = inv_wishart_cholesky_lpdf(LX | nu, L_scale_matrix) + // prior
sum(log(diagonal(LX))) + // Jacobian
sum(log_lik);
}
'
cmdstanr::cmdstan_model(cmdstanr::write_stan_file(code))
}
create_data <- function() {
N <- 2000
K <- 5
if (is.null(df))
df <- K
Sigma <- diag(rep(1+1/K, K)) - 1/K
LX <- t(chol(Sigma))
scale_matrix <- diag(rep(1, K))
nu <- K
set.seed(23478)
X <- lapply(1:N, function(i){ as.vector(mvtnorm::rmvnorm(1, sigma=Sigma)) })
N1 <- N %/% 100
mu <- 1.0
X[1:N1] <- lapply(1:N1, function(i){ mu + as.vector(mvtnorm::rmvnorm(1, sigma=16*Sigma))})
list(
K = K,
N = N,
X = X,
nu = nu,
scale_matrix = scale_matrix
)
}
dummy_model <- create_model()
dummy_data <- create_data()
dummy_fit <- dummy_model$sample(data = dummy_data)
dummy_loo <- dummy_fit$loo(moment_match = TRUE)
Expected behavior
I expect the loo computation to proceed without error.
Operating system
MacOS Ventura 13.4
CmdStanR version number
0.6.0
Additional context
I have seen problems like this with other models that had cholesky_factor_cov parameters. This is just the simplest example I could create.