I want to model a matrix as a multivariate normal with univariate normal priors on the column means and an LKJ prior on the covariance matrix. Something like this works fine, as expected:
df = (df - df.mean()) / df.std() # necessary!
m = len(df.columns)
mu = pm.Normal("mu", 0, sd=1, shape=m)
sd_dist = pm.HalfCauchy.dist(beta=1, shape=m)
chol_packed = pm.LKJCholeskyCov("chol_packed", n=m, eta=1, sd_dist=sd_dist)
chol = pm.expand_packed_triangular(m, chol_packed)
pm.MvNormal("Y", mu=mu, chol=chol, observed=df)
However, if I remove the standardization step (i.e., so columns no longer have zero mean and unit variance) sampling goes haywire, with poor convergence, autocorrelation, etc. I have narrowed down the problem to the means of the columns being highly variable.
At first I thought this was related to funneling but using a non-centered parameterization on the means did not fix the issue. Any ideas?
