Hi,

I am trying to fit the model below where I propagate uncertainty about the predictors through the model. I’ve a working version in which I fit the models separately and plug in the central tendency for xbar but obviously that’s unsatisfactory. Do you have any advice on how to speed this up?

```
data {
int J;
int G;
matrix[J, G] xbar1;
matrix[J, 3] ybar;
}
parameters {
matrix<lower = 0, upper = 1>[G,3] beta;
simplex[G] xbar2[J];
simplex[G] xbar3[J];
vector<lower = 0>[3] sigma;
}
model {
matrix[J, G] xbar2_mat;
matrix[J, G] xbar3_mat;
for (i in 1:3) target += normal_lpdf(beta[:,1] | 0.5, 0.25);
for (j in 1:J){
xbar2[j,:] ~ dirichlet(softmax(xbar1[j]' .* beta[:,1])); //softmax(xbar1[j]' .* beta[:,1])
xbar3[j,:] ~ dirichlet(softmax(xbar2[j] .* beta[:,2])); //softmax(xbar2[j]' .* beta[:,2])
xbar2_mat[j] = xbar2[j,:]';
xbar3_mat[j] = xbar3[j,:]';
}
target += normal_lpdf(sigma | 0, 2);
target += normal_lpdf(ybar[:,1] | xbar1 * beta[:,1], sigma[1]);
target += normal_lpdf(ybar[:,2] | xbar2_mat * beta[:,2], sigma[2]);
target += normal_lpdf(ybar[:,3] | xbar3_mat * beta[:,3], sigma[3]);
}
generated quantities {
matrix[J, 3] y_bar_hat;
y_bar_hat[,1] = to_vector(normal_rng(xbar1 * beta[:,1], sigma[1]));
for (j in 1:J){
y_bar_hat[j,2] = normal_rng(xbar2[j]' * beta[:,2], sigma[2]);
y_bar_hat[j,3] = normal_rng(xbar3[j]' * beta[:,3], sigma[3]);
}
}
```