model4.stan

data {
  int<lower=0> n;
  vector[n] y;
  vector[n] x2;
  vector[n] x3;
  real<lower= -1, upper= 1> rho;
  real g; // slope for u
  vector[2] prior;
}
transformed data {
  matrix[2,2] sigx;
  vector[2] mux;
  sigx[1,1] = variance(x2); 
  sigx[1,2] = sd(x2)*rho;
  sigx[2,1] = sigx[1,2]; 
  sigx[2,2] = 1.0;
  mux[1] = mean(x2);
  mux[2] = 0.0;
}
parameters {
  real b0;
  real<lower=0> b1;
  real<lower=0> b2;
  real<lower=0> sig;
  vector[n] u;
}
transformed parameters {
  // multi_normal needs vector
  vector[2] xx[n];
  for(i in 1:n) {
    xx[i,1] = x2[i];
    xx[i,2] = u[i];
  }
}
model {
  real eta[n];
  b0 ~ normal(100, 100);
  b1 ~ normal(0, prior[1]);
  b2 ~ normal(0, prior[2]);
  sig ~ gamma(1, 1);
  for (i in 1:n) {
    eta[i] = b0 + b1*xx[i][1] + g*xx[i][2] + b2*x3[i];
  }
  xx ~ multi_normal(mux, sigx);
  y ~ normal( eta, sig);
}