Thanks to everyone for the comments. I have been pondering this for
the last few days and believe I understand the nature of my problem
now. As several of you said, the issue was that my data generation was
not consistent with my Stan model.
I'll summarise my understanding in the hope that (a) it is correct and
(b) it is of some use to others.
In my examples, I have some latent quantity x sampled from a
distribution of known form, with parameters θ_x, and then some
observed quantity sampled from a known distribution centred on x, with
parameters θ_y. The data are only collected if y is greater than
some threshold y_min.
I'll denote the pdf of x as P_x, and the pdf of y as P_y. In my first
toy model, P_x and P_y were both normal while in my semi-realistic
model, P_x was lognormal and P_y was Poissonian.
I generated my data using something like
x = rnorm(N_tot,mu,sigma_x)
y = rnorm(N_tot,x,sigma_y)
and applied my detection threshold with
y = y[y>ymin]
N_obs = length(y)
Now consider the likelihood for y. First, neglecting the truncation at
y_min, the likelihood of a particular pair of x and y is
P_1(y,x|θ_x,θ_y) = P_x(x|θ_x) * P_y(y|x,θ_y)
and we can marginalise out x to get
P_2(y|θ_x,θ_y) = \int dx P_x(x|θ_x) * P_y(y|x,θ_y)
Now, taking the truncation into account, we set P_2(y|θ_x,θ_y)=0 for y<=y_min and
renormalise the distribution to give
P_trunc(y|θ_x,θ_y,y_min) = P_2(y|θ_x,θ_y) / \int_y_min^inf dy P_2(y|θ_x,θ_y)
The problem I had was that my Stan model for truncated data had the form
x ~ normal(mu,sigma_x);
for (i in 1:N_obs)
y[i] ~ normal(x[i], sigma_y) T[ymin,];
but this does not correspond to the likelihood above. Instead, this
P_trunc(y|θ_x,θ_y,y_min) = \int dx P_x(x|θ_x) * P_y(y|x,θ_y) / \int_y_min^∞ dy P_y(y|x,θ_y)
i.e. in this case it is only the distribution of y that is truncated,
while for the first case [equation 1], the x-marginalised distribution
of y is truncated, as in my generated data.
In the simple case of P_x and P_y being normal, then P_trunc
simplifies to a single normal with the variances added to give
sigma_tot^2 (as pointed out earlier in the thread), and you can just
y[i] ~ normal(mu,sigma_tot) T[ymin,]
For the general case where P_x and P_y are not both normal, I don't
know a way of writing a model corresponding to [equation 1] concisely
in Stan, (i.e. keeping the x's as latent parameters rather than
integrating them out). An inelegant version in which I compute the
integrals numerically in user-defined functions, correctly recovers the
parameters in my lognormal-poisson generated data.
Any thoughts on whether it is possible to write the general case in
Stan would be gratefully received!
Thanks for all of your time and input!