gqs() does not work with transformed parameters

[maxmantei]

Summary:

Standalone generated quantities via qgs() don't work with transformed parameters.

Description:

See this discourse topic: Using transformed parameters leads to error Wrong number of parameter values in draws from fitted model. Expecting 2 columns, found 3 columns. which, I guess, comes from stan/src/stan/services/sample/standalone_gqs.hpp.

Reproducible Steps:

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

sc <- "
transformed data{
  vector[10] y = [4.65, 6.02, 4.92, 4.77, 8.12, 6.93, 7.39, 7.6, 5.68, 2.14]';
}
parameters{
  real mu;
  real log_sigma;
}
transformed parameters{
  real<lower=0> sigma = exp(log_sigma);
}
model{
  y ~ normal(mu, sigma);
  mu ~ normal(0, 10);
  log_sigma ~ normal(0, 2.5);
}

"

sm <- stan_model(model_code = sc)

post <- sampling(sm)

sc_gq <- "
parameters{
  real mu;
  real log_sigma;
}
transformed parameters{
  real<lower=0> sigma = exp(log_sigma);
}
generated quantities{
  vector[10] y_rep;
  
  for (i in 1:10)
    y_rep[i] = normal_rng(mu, sigma);
}

"

sm_gq <- stan_model(model_code = sc_gq)

rep <- gqs(sm_gq, draws = as.matrix(post))

rep

Current Output:

Wrong number of parameter values in draws from fitted model.  Expecting 2 columns, found 3 columns.

Inference for Stan model: 9a372375694ef14366f13bf8c7378254.
1 chains, each with iter=4000; warmup=0; thin=1; 
post-warmup draws per chain=4000, total post-warmup draws=4000.

          mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
y_rep[1]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[2]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[3]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[4]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[5]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[6]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[7]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[8]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[9]     0     NaN  0    0   0   0   0     0   NaN  NaN
y_rep[10]    0     NaN  0    0   0   0   0     0   NaN  NaN

Samples were drawn using  at Tue Nov 05 15:31:10 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Expected Output:

Inference for Stan model: bdf3c0b7af993ed90b29cc81706a90ea.
1 chains, each with iter=4000; warmup=0; thin=1; 
post-warmup draws per chain=4000, total post-warmup draws=4000.

          mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
y_rep[1]  5.82    0.04 2.16 1.54 4.51 5.82 7.08 10.20  3487    1
y_rep[2]  5.83    0.03 2.14 1.54 4.51 5.85 7.13  9.99  3851    1
y_rep[3]  5.82    0.04 2.16 1.64 4.44 5.83 7.13 10.17  3743    1
y_rep[4]  5.75    0.04 2.19 1.36 4.48 5.75 7.10 10.09  3354    1
y_rep[5]  5.85    0.03 2.14 1.77 4.48 5.81 7.17 10.21  3879    1
y_rep[6]  5.82    0.03 2.14 1.52 4.52 5.83 7.15 10.07  3830    1
y_rep[7]  5.83    0.04 2.17 1.60 4.52 5.82 7.13 10.07  3509    1
y_rep[8]  5.86    0.03 2.08 1.64 4.52 5.85 7.17 10.09  3700    1
y_rep[9]  5.83    0.04 2.13 1.54 4.53 5.82 7.14 10.11  3362    1
y_rep[10] 5.77    0.03 2.11 1.72 4.47 5.76 7.08 10.01  3740    1

Samples were drawn using  at Tue Nov 05 15:34:33 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

RStan Version:

> packageVersion("rstan")
[1] ‘2.19.9’

R Version:

> R.version.string
[1] "R version 3.6.1 (2019-07-05)"

Operating System:

Win7 Pro SP1

[bob-carpenter]

Thanks for filing a careful issue. To work properly, the parameters have to be the same, but everything else could vary in a model used for standalone generated quantities. So I think you're right that we could allow arbitrary transformed parameters. We should at least allow the transformed parameters from the original Stan program.

To work around for now, you can move that transformed parameter definition into the generated quantities block and it should be OK.

[maxmantei]

To work around for now, you can move that transformed parameter definition into the generated quantities block and it should be OK.

This is fine for simple models, but a bit annoying for more complex ones. (Maybe I should just do that with different #includes for now...) I think having a workflow where you just have to "swap out" the model block with a generated quantities block would be really great. :)

[bob-carpenter]

Right---that was the original intent---just swap in a new generated quantities block. It should be possible to run standalone generated quantities with only the original parameters declaration and a generated quantities block. If the data block is included in the standalone gq program, the data should be read and available for generated quantities. It doesn't need to match the original data. That's how we should be able to do predictive inference---read in new x_tilde and generate new y_tilde and don't read in original data (which might be large). The transformed parameters declarations and definitions could conceivably change from the original program. It'd probably be better to make that kind of change directly in the generated quantities block. So the question is whether we should flag it as an error, warning, or just allow it without comment as things run. I don't think the model needs to be executed in standalone generated quantities (even if we did execute it, the constants change and target() wouldn't be the same). So it should be optional in the standalone generated quantities programs and perhaps raise a warning. I think we should probably allow lp__ to be read in so that target() is well defined in the gq block. If any of this isn't happening, I'd consider it a bug (other than lp__, which I don't think we ever discussed as a feature).

…

On Nov 7, 2019, at 11:47 AM, Max Mantei ***@***.***> wrote: To work around for now, you can move that transformed parameter definition into the generated quantities block and it should be OK. This is fine for simple models, but a bit annoying for more complex ones. (Maybe I should just do that with different #includes for now...) I think having a workflow where you just have to "swap out" the model block with a generated quantities block would be really great. :) — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[bnicenboim]

Hi,
I was having the same issue and I found this bug.
An easier workaround is to do this:

post <- sampling(sm)

sc_gq <- "
parameters{
  real mu;
  real log_sigma;
  real<lower=0> sigma;
}
generated quantities{
  vector[10] y_rep;
  
  for (i in 1:10)
    y_rep[i] = normal_rng(mu, sigma);
}

"

Moving transformed data to the generated quantities block is problematic, because the output saves more variables. I was using this feature because the output was already huge and I wanted to limit the iterations.

[mguzmann]

No progress on this?