Does anyone have some plain code for PPCA implemented in Stan that I could use for a more complicated hierarchical problem? I’m aware of this code, but something with less fancy footwork would be preferable because it’s complicated already.

Looks like page 34 of the ADVI paper has some super simple Stan code for PPCA:

Although I’ll give this one a go too:

We have code for PPCA, either plain or reparameterized via Housholder transformations to remove the unidentifiable rotation on the latent space

Thank you! I actually found this already. Do you think it’d be simple enough for me to extend your code to MPPCA?

This variant seems nice, although I understand it won’t be identifiable due to rotational symmetry. It has a massive advantage over the code in the ADVI paper (bottom of response) namely that it doesn’t have more parameters than data points… This makes it near impossible to generate a posterior sample for a moderately large dataset (e.g. 200k rows) without lots of RAM.

It’s also pretty straightforward to add a dimension to mu, sigma_noise and W in order to create a mixture of PPCAs although it remains to be seen how well it’ll work.

    int<lower=0> N;
    int<lower=1> D;
    int<lower=1> Q;
    vector[D] Y[N];
transformed data{
parameters {
    matrix[D, Q] W;
    real<lower=0> sigma_noise;
	vector[D] mu;
transformed parameters {
    cholesky_factor_cov[D] L;
        //matrix[D, D] K = W*W'; // tcrossprod(matrix x)
        matrix[D, D] K = tcrossprod(W);
        for (d in 1:D)
            K[d,d] += square(sigma_noise) + 1e-14;
        L = cholesky_decompose(K);
	mu ~ normal(0, 10);
    sigma_noise ~ normal(0,1);
    to_vector(W) ~ normal(0,1);
    Y ~ multi_normal_cholesky(mu, L);

PPCA from the ADVI paper:

data {
  int <lower=0> N;
  int <lower=0> D;
  int <lower=0> M;
  vector[D] x[N];
parameters {
  matrix[M,N] z;
  matrix[D,M] w;
  real <lower=0> sigma;
  vector<lower=0>[M] alpha;
model {
  for(d in 1:D) w[d]~normal(0,sigma*alpha);
  for(n in 1:N) x[n]~normal(w*col(z,n),sigma);