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?