Using multivariate normal with widely varying means

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?

did you try to set the test val of mu? It might help but I am not sure how much.

That actually helped a lot, I’m embarrassed not to have thought of it. Thanks!

You are welcome! I think if you tune for a really long time it should be able to find it but