I have a scenario where I have multiple correlated outcomes and I think use of Bayesian multivariate regression (in the seemingly unrelated regression framework) is appropriate. I’ve implemented the following code after encountering it in the guide. My question, though, is about the underlying mechanics of what this code is achieving, especially the ‘model’ component. I’d appreciate if someone could walk me through that.

data {
int<lower=1> K;
int<lower=1> J;
int<lower=0> N;
vector[J] x[N];
vector[K] y[N];
}
parameters {
matrix[K, J] beta;
cov_matrix[K] Sigma;
}
model {
vector[K] mu[N];
for (n in 1:N)
mu[n] = beta * x[n];
y ~ multi_normal(mu, Sigma);
}

The Stan model is basically an elaborate way of writing \log p(x | \theta) p(\theta) where \theta are the parameters you want to sample and x is your data. Once you have that abstraction, it’s MCMC time.

Michael Betancourt’s HMC paper here is the most up to date documentation on the Stan algorithms https://arxiv.org/abs/1701.02434 (other than the code itself).