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.
matrix[K, J] beta;
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.
For animations of algorithms like the ones in Stan, check out the NUTS samplers here: https://chi-feng.github.io/mcmc-demo/app.html
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).
I like Radford Neal’s paper a lot too: https://arxiv.org/abs/1206.1901
The Stan manual has some info on the samplers, but the above references are probably better: https://mc-stan.org/docs/2_24/reference-manual/hmc-chapter.html