design.stan

data{
 int<lower=0> N;
 vector[2] X1[N];
 vector[2] X2[N];
}
parameters{
  real mu;
  real<lower=0> sig[2];
  real<lower=-1,upper=1> rho;
  real effectA; //between
  real effectB; //within
  real effectAB;//interaction
}
transformed parameters{
  cov_matrix[2] Sigma;
  vector[2] mu2;
  vector[2] ideal1;
  vector[2] ideal2;
  Sigma[1,1] = sig[1]*sig[1];
  Sigma[2,2] = sig[2]*sig[2];
  Sigma[1,2] = sig[1]*sig[2]*rho;
  Sigma[2,1] = Sigma[1,2];
  mu2[1] = mu;
  mu2[2] = mu + effectB;
  ideal1[1] = mu2[1] + effectA + effectAB;
  ideal1[2] = mu2[1] - effectA - effectAB;
  ideal2[1] = mu2[2] + effectA - effectAB;
  ideal2[2] = mu2[2] - effectA + effectAB;
}
model{
  X1 ~ multi_normal(ideal1,Sigma);
  X2 ~ multi_normal(ideal2,Sigma);
  mu ~ normal(60,100);
  rho ~ uniform(-1,1);
  sig ~ cauchy(0,2.5);
}