Ppca

Modeling
1

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.

2

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:

3

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

4

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

5

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.

data{
    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);
    }
}
model{
	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 {
  to_vector(z)~normal(0,1);
  for(d in 1:D) w[d]~normal(0,sigma*alpha);
  sigma~lognormal(0,1);
  alpha~inv_gamma(1,1);
  for(n in 1:N) x[n]~normal(w*col(z,n),sigma);
}