Hello!
This is a follow-up to Faster / better loading of sampler diagnostics in cmdstanr?.
I’ve fit a model where one of the parameters of interest is a hierarchical GP that has a varying effect across a predictor grid d according to a variable typgrp.
parameters {
matrix[n_d, n_typgrp] GP_typgrp_std; // standard mvnorm
real<lower=0> alpha_typgrp; // marginal deviance
real<lower=0> rho_typgrp; // length-scale
...
}
transformed parameters {
...
matrix[n_d, n_t] GP_typgrp; // typgrp GPs
// typgrp-level non-centered Gaussian Processes (mvnorm std * L_cov)
{
matrix[n_d, n_d] cov_typgrp; // typgrp covariance
matrix[n_d, n_d] L_cov_typgrp; // lower tri
cov_typgrp = cov_exp_quad(ds, alpha_typgrp, rho_typgrp)
+ diag_matrix(rep_vector(1e-9, n_d)); // avoid float errors
L_cov_typgrp = cholesky_decompose(cov_typgrp);
GP_typgrp = L_cov_typgrp * GP_typgrp_std; // non-centered
}
...
}
I can successfully fit the model in cmdstanr, and I can recover the parameters and everything looks great - my variable GP_typgrp varies across d by typgrp as expected. Very happy!
As suggested in my previous post, though, for the sake of efficiency, I’ve split out the generated quantities block into a separate .stan program, and this is where things fall apart. My generated quantities program is identical to the model program minus the model block and with the following gq block:
generated quantities {
matrix[n_d, n_t] eta_new;
eta_new = GP_typgrp;
}
(this is a simplified version for illustration)
My call to the gq_program is as follows, where the variable fit is the sampled CmdStanMCMC object:
gq_program <- file.path(..., "gq_typgrp.stan")
mod_gq <- cmdstan_model(gq_program)
gq <- mod_gq$generate_quantities(
fit,
data = data,
parallel_chains = 4
)
When I run the generated quantities program, the resulting matrix eta_new is a matrix with only the first column repeated n_t times. That is to say that the value of eta_new over d is identical across typgrp - different than the GP_typgrp variable recovered as a posterior sample in the model program.
I’m using the following {posterior} + {tidybayes} workflow to look at the draws:
post_draws <- gq$draws()
eta_draws_df <- post_draws %>%
posterior::subset_draws(variable = "eta_new") %>%
posterior::as_draws_df() %>%
tidybayes::gather_variables() %>%
ungroup() %>%
mutate(
d = as.integer(
stringr::str_extract(.variable, "(?<=\\[)[:digit:]+")
),
typgrp = as.integer(
stringr::str_extract(.variable, "(?<=,)[:digit:]+")
),
.variable = "eta"
) %>%
arrange(d, .draw)
I think I must be misunderstanding something basic about how generated quantities handles / accepts matrices, or am missing something else. I can’t post the data / full model, but can try to make a reproducible example if no obvious operator error jumps out. Any advice is appreciated, with thanks!
[edited to add relevant parameters declaration and figures]