Unique solution for probabilistic PCA

@ferrine
I generated some fake data

import numpy as np
import datetime
from pymc3.variational.inference import ADVI
from pymc3.variational.opvi import Group, node_property
from pymc3.variational.approximations import MeanFieldGroup, SingleGroupApproximation
import pymc3 as pm
import theano
N=500
d = np.random.normal(size=(3, N))
d1 = (d[0]>0).astype(int)
d2 = -1*d[1]+np.random.normal(size=(1,N))/2
d3 = 2*d1+np.random.normal(size=(1,N))
d4=.5*d2+np.random.normal(size=(1,N))
yd = 4.*d1 +7.*d2 +5*d3+.3*d4+d[2]
#yd=.5*d1+d[2]
#d1[[1,5,25,14,32,11,125,157, 13,46,81,99,168,101,145,111]]=-999
d1True=d1.copy()
#d1[:450]=-999.
d1

Here is the model code:

import theano.tensor as tt
obs=np.array([d1,d2[0],d3[0],d4[0], yd[0]]).T
with pm.Model() as model1:
    s = pm.HalfCauchy('s', 10., shape=5)
    latent_std=10#HalfCauchy('l_std',100.)
    w = pm.Normal('w',0.,10, shape=(5,1))
    
    latent_ = pm.Normal('latent_',0.,1., shape=(len(d1),1))
    latent = pm.Deterministic('latent', latent_*latent_std)
    p = pm.Deterministic('p',tt.dot(latent, w.T))
    cov = pm.Deterministic('cov', latent_std**2*tt.dot(w, w.T)+(s**2)*np.eye(5))
    sd=pm.Deterministic('sd', tt.sqrt(tt.diag(cov)))
    latent_like = pm.Normal('latent_like', mu=p, sd=s, observed=obs)

with model1:
    approx=pm.fit(method=SADVI())
    trace=approx.sample()

I know for sure that the parameter w exhibits sign-switching, at least. Make sure to sample with NUTS 4 chains to reliably see the sign-switching occur. Or run ADVI to see parameters regularized to 0.

Just FYI, in this model, w is the matrix of factor loadings (5,1)