Stan slow with latent-variable model for rounded data

[Andrew Gelman]

Hi, I mentioned this one on the Stan list a few weeks ago, and recently I explored it further. Stan and R code are attached but for simplicity I'm also including the Stan program here:

data {
int N;
vector[2] y[N];
vector[2] mu;
cov_matrix[2] Sigma;
}
parameters {
vector<lower=-0.5, upper=0.5>[2] round_err[N];
}
transformed parameters {
vector[2] z[N];
for (n in 1:N) {
z[n] = y[n] - round_err[n];
}
}
model {
z ~ multi_normal(mu, Sigma);
}

Step 1:

The data are heights of mothers and daughters, modeled as coming from a bivariate normal with mean vector mu and cov matrix Sigma, with all heights rounded to the nearest inch. For simplicity I'm assuming mu and Sigma are known.

I follow the recommended approach in which the rounding errors are bounded in the range (-0.5, +0.5), and then the model above allows us to estimate the underlying continuous heights z. (Because the data are bivariate, we can't easily integrate out z using the normal cdf.)

I fit the model with some simulated data and it's slow. For example, with N=1000 and running the default 1000 warmup + 1000 saved iterations per chain, it took 150-270 seconds per chain (mostly in warmup, of course).

When this came up before, Michael suggested that the problem is at the boundaries of the rounding errors, thus in the tails of these parameters on the transformed scale.

Step 2:

I was curious so I re-ran but this time removing the multivariate normal from the model, thus simply imputing the rounding errors from their uniform prior distributions. This time, it was much faster, only taking 75-100 seconds per chain. At first I was thinking that the small amount of information provided by that multivariate normal in Step 1 above was making the geometry worse. Or maybe it's just that the multivariate normal is expensive for Stan to evaluate.

Question:

Is there any good solution to this problem? It's not a huge deal, but we do have this rounded-data model in the manual so it would be good to have a clean solution. I think that Bob at one point suggested hacking the problem by removing the hard boundaries and letting the transitions be soft. This violates the rounding model but I suppose that, realistically, when there's rounding error there's also measurement error. That's certainly the case with height.

Anyway, I was wondering what people's thoughts are.

See you
A