# The parameters of model can not be recovered

I simulated a set of data

library(tidyverse)

df <- tibble(
x    = seq(from = -13,  to = 13, by = 1),
z    = rnorm(n = 27, mean =  4.8 + 0.35 * x, sd = 2),
y    = rpois(n = 27, lambda = exp(0.18 - 0.13 * z))
)

# A tibble: 5 × 3
#      x      z     y
#  <dbl>  <dbl> <int>
#1   -13  1.57      1
#2   -12 -0.705     1
#3   -11  4.47      0
#4   -10  0.123     0
#5    -9 -0.471     3


I want to recover the parameters \beta s,

\begin{aligned} y_{i} &\sim \text{Poisson}\left( \lambda_i\right) \\ \lambda_i &=\exp \left(\beta_1 + \beta_2 z_i\right) \\ z_i & \sim \text{Normal}\left( \beta_3 + \beta_4 x_i, \sigma\right) \\ \end{aligned}
data {
int<lower=1> N;
vector[N] x;
array[N] int y;
}

parameters {
vector[4] beta;
real<lower=0> sigma;
}
model {
vector[N] z;

beta ~ normal(0, 20);
sigma ~ exponential(1);

z ~ normal(beta[3] + beta[4] * x, sigma);
y ~ poisson_log(beta[1] + beta[2] * z);
}



but the model failed, pop up a message without samping

Running MCMC with 4 sequential chains...


I don’t know how to solve it. Thank you very much for your help.

code.Rmd (1.4 KB)

You need to move z to the parameters block.

To see why, recall that z ~ normal(beta[3] + beta[4] * x, sigma) is just a shorthand for target += normal_lupdf(z | beta[3] + beta[4] * x, sigma). But z is just an empty vector at this point; it doesn’t contain anything. You’re expecting z ~ normal(beta[3] + beta[4] * x, sigma) to define z to be a sample from the distribution on the right-hand side of the ~, but that’s not how Stan works.

The lack of an informative error message here is definitely not ideal, though. @WardBrian @stevebronder @rok_cesnovar , on my system if I run the rmd provided in the OP I can compile the model and run the line

fit <- model\$sample(data = stan_data, refresh = 1)


with no warning messages or errors. Clearly it doesn’t sample; the only output to the console is

Running MCMC with 4 sequential chains...


But I only see an error if I then run, e.g.

> fit
Error: Fitting failed. Unable to print.


I have range checks turned off; but I don’t think that should matter?

> cmdstan_make_local()

@Wang_MinJie , as an aside, even once you fix this I don’t think you can expect this model to do a good job of parameter recovery. The parameters are not identified except via the prior (which is very weak). In particular, I can multiply beta[3], beta[4], and sigma by a constant and divide beta[2] by the same constant and the model will be completely unchanged.