data {
int<lower=0> N;
int<lower=1> K;
int<lower=1> J;
int<lower=1> L;
int<lower=1,upper=J> jj[N];
matrix[N,K] x;
row_vector[L] u[J];
vector[N] y;
}
parameters {
corr_matrix[K] Omega;
vector<lower=0>[K] tau;
matrix[L,K] gamma;
vector[K] beta[J];
real<lower=0> sigma;
}
model {
tau ~ cauchy(0,2.5);
Omega ~ lkj_corr(2);
to_vector(gamma) ~ normal(0, 5);
{
row_vector[K] u_gamma[J];
for (j in 1:J)
u_gamma[j] <- u[j] * gamma;
beta ~ multi_normal(u_gamma, quad_form_diag(Omega, tau));
}
{
vector[N] x_beta_jj;
for (n in 1:N)
x_beta_jj[n] <- x[n] * beta[jj[n]];
y ~ normal(x_beta_jj, sigma);
}
}