Consider a standard factor model such as the one taken from this blog:
data {
int<lower=1> n;
int<lower=1> nX;
vector [n] y;
matrix [n,nX] X;
}
parameters {
vector[nX] beta;
real<lower=0> sigma;
}
transformed parameters {
vector[n] mu;
mu = X*beta;
}
model {
// Likelihood
y~normal(mu,sigma);
// Priors
beta ~ normal(0,100);
sigma~cauchy(0,5);
}
generated quantities {
vector[n] log_lik;
for (i in 1:n) {
log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);
}
}
Now consider that one wants to add correlated residuals (error terms) between items in this model. What would be the best way to do so? Do you need to literally calculate the residuals in the generated quantities
block and then assume a multivariate normal distribution for them in the model
block? Or should this be specified in some other way?