Drawing a Cholesky-decomposed correlation matrix

Questions
1

I would like to draw a Cholesky-decomposed correlation matrix using the LKJ prior. I have tried to use the LKJCholeskyCov with an sd_dist=pm.Beta.dist(10**3,10**-3), as indicated in Uses of LKJCholeskyCov and LKJCorr - #2 by chartl . However, this leads to a warning for “Bad initial energy”.

Can you tell what the problem is? Is there any other way to draw the decomposed correaltion matrix?

2

Have you tried an alternative parameterization (halfnormal, then rescaling)? If the initial energy is still bad, then the error is occurring somewhere else in the model. If not, then the spike at 1 in the beta dist is wreaking havoc (it might if the matrix is somewhat large).

3 
sd_dist = pm.HalfNormal.dist(sd=1)
packed_chol = pm.LKJCholeskyCov('chol_cov', eta=2, n=n, sd_dist=sd_dist)
chol = pm.expand_packed_triangular(n, packed_chol, lower=True)
cov = pm.Deterministic('cov', tt.dot(chol, chol.T))
sd = tt.sqrt(tt.diag(cov))
corr = cov / sd[:, None] / sd[None, :]
4

Thanks a lot for your reply.

I would like to avoid using the piece of script you provided above, as it increases the runtime. Just for context, I am running a Hierarchical Bayesian Model with several units and I need a correlation matrix per unit. Using the beta spike appeared as a computationally cheaper approach.

At the moment, the correlation matrix per unit is 5x5. The beta spike method didnt work even when running the model for just 1 unit. I should also mention that testing LKJ_sd_dist=pm.Beta.dist(2,2) also lead to “Bad initial energy”. I am posting the model below (please note that my data has been standardized):

LKJ_sd_dist=pm.Beta.dist(10** 3,10**-3)
eta=2

with pm.Model() as train_model:
chol_packed = pm.LKJCholeskyCov(‘chol_packed’, n=data.shape[1], eta=eta, sd_dist=LKJ_sd_dist)
chol = pm.expand_packed_triangular(n=data.shape[1], packed=chol_packed)
obs = pm.MvNormal(‘obs’, mu=0, chol=chol, observed=data)

5

I assume that it works, then. Is your data normalized to have unit standard deviation? If not, then you need a more flexible sd distribution like pm.HalfNormal.

6

The following works only when using LKJ_sd_dist=pm.Exponential.dist(1):

The data is normalized to unit standard deviation. Do you think that the beta spike method can be applied to avoid using the script above?

7

That seems particularly strange, since the standard deviation is normalized out. Have you tried a normal spike like pm.Normal.dist(mean=1, sd=1e-3) (maybe bound it to be more careful)?

8

Using a spike on a Normal distribution worked! Thanks a lot.