Suppose y_i \sim N(\mu, \sigma) (i = 1, \ldots, N) but we observe y^*_i = \lfloor y_i \rfloor.
A model for inference about \mu, \sigma can be based on integration, i.e.:
p(y^*_i \mid \mu, \sigma) = \int^{y^*_i + 1}_{y^*_i} \mathrm{normal}(y \mid \mu, \sigma) dy
In Stan code:
data {
int<lower=1> N;
array[N] int y_star;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
for (n in 1:N) {
target += log_diff_exp(normal_lcdf(y_star[n] + 1 | mu, sigma),
normal_lcdf(y_star[n] | mu, sigma));
}
}
However, because I want to eventually adapt this idea to a different problem, I would like to “estimate” the unobserved y_i, similar to the approach taken in the Stan user’s guide.
In that example there is a fixed censoring value, so it uses a simple constraint in the parameter declaration. I note that it is possible to set varying constraints, but again because of my planned extension to another case, I would like to do it in a more complicated way: basically using an inverse transform sample:
y_i = \mu + \sigma \Phi^{-1}\left( \Phi(\frac{y^*_i - \mu}{\sigma}) + u_i (\Phi(\frac{y^*_i + 1 - \mu}{\sigma}) - \Phi(\frac{y^*_i - \mu}{\sigma})) \right) where u_i \sim U(0,1).
If I want y_i \sim N(\mu, \sigma), I need to compute the Jacobian:
\frac{dy_i}{du_i} = \sigma \frac{\Phi\left(\frac{y^*_i + 1 - \mu}{\sigma}\right) - \Phi\left(\frac{y^*_i - \mu}{\sigma}\right)}{\phi\left(\frac{y_i - \mu}{\sigma}\right)}
The other parameters aren’t transformed and the Jacobian is triangular, so its log-determinant is just the log of the value above. In Stan:
data {
int<lower=1> N;
array[N] int y_star;
}
parameters {
real mu;
real<lower=0> sigma;
array[N] real<lower=0,upper=1> u;
}
transformed parameters {
array[N] real y;
for (n in 1:N) {
y[n] = mu + sigma * inv_Phi(normal_cdf(y_star[n] | mu, sigma) +
u[n] * (normal_cdf(y_star[n] + 1 | mu, sigma) - normal_cdf(y_star[n] | mu, sigma)));
}
}
model {
y ~ normal(mu, sigma); // (0)
for (n in 1:N) {
// Jacobian adjustment
target += log(sigma); // (1)
target += log_diff_exp(normal_lcdf(y_star[n] + 1 | mu, sigma),
normal_lcdf(y_star[n] | mu, sigma)); // (2)
target += -normal_lpdf(y[n] | mu, sigma); // (3)
}
}
This model doesn’t work well. Line (2)
is the same as in the integrated case. Line (3)
cancels out Line (0)
. So the problem is Line (1)
(the model works if I remove it), but I am not sure why.
Does my logic make sense? Where have I gone wrong?