pairedT.stan

data{
  int<lower=0> N;
  vector[2] X[N];
}
parameters{
  vector[2] mu;
  vector<lower=0>[2] sigma;
  real<lower=-1,upper=1> rho;
}
transformed parameters{
  cov_matrix[2] Sigma;
  Sigma[1,1] = sigma[1]^2;
  Sigma[2,2] = sigma[2]^2;
  Sigma[1,2] = sigma[1] * sigma[2] * rho;
  Sigma[2,1] = sigma[2] * sigma[1] * rho;
}
model{
  X ~ multi_normal(mu,Sigma);
}