Probabilistic principan component analysis using Stan

PCA.stan

data {
  int<lower=0> N;
  int<lower=0> D;
  int<lower=0> R;
  matrix[N,D] Y;
}
parameters {
  row_vector[D] mu;
  matrix[D,R] W;
  real<lower=0> sig2;
}
model {
  for(i in 1:N){
    Y[i,:] ~ multi_normal(mu, add_diag(W*W', sig2));    
  }
  to_vector(W) ~ normal(0, 1);
  sig2 ~ student_t(3, 0, 1);
}
generated quantities {
  matrix[R,N] Z;
  matrix[R,R] M = inverse(add_diag(W'*W, sig2));
  for(i in 1:N){
    Z[:,i] = multi_normal_rng((Y[i,:]-mu)*W*M, sig2*M);
  }
}