How to construct correlation matrix from pymc.LKJCorr?

Question

I want to construct a correlation matrix explicitly from using the pymc.LKJCorr distribution class, but I don't trust my understanding of the pymc.expand_packed_triangular. Here is a minimal working example.

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm


with pm.Model() as model:
    R_upper = pm.LKJCorr('LKJCorr', n=2, eta=2)
    R = pm.expand_packed_triangular(n=2, packed=R_upper, lower=False)
    corr = pm.Deterministic('RCorr', var=R + np.eye(2))

with model:
    idata = pm.sample()
    az.plot_trace(idata, combined=True)
    plt.tight_layout()
    plt.show()

Here is an example plot:

What doesn't make sense to me is the fact that I am seeing one of the RCorr parameters centered around one. That shouldn't be the case, and looks like it is just a translation of "LKJCorr" by unity. I am concerned that pm.expand_packed_triangular is assuming there is a diagonal when there isn't, or something like that.

How can I reconstruct the correlation matrix from using an instance of pymc.LKJCorr in this toy example?

Galen · Accepted Answer · 2023-06-25 01:09:59Z

0

In this simple case of a 2x2 correlation matrix you can do the following

with pm.Model() as model:
    R_upper = pm.LKJCorr("LKJCorr", eta=1, n=2)


    corr = pm.Deterministic(
        "RCorr", 
        at.fill_diagonal(R_upper[np.zeros((2, 2), dtype=np.int64)], 1)
    )

It works by creating a 2x2 array whoses elements are all the correlation score, which is fine here since there is only one, and then filling the diagonal with ones.

Not sure how to do the n-dimensional case yet.

answered Jun 25, 2023 at 1:09

Galen

1,3942 gold badges16 silver badges42 bronze badges

Sign up to request clarification or add additional context in comments.

Comments

velochy · Accepted Answer · 2024-01-09 14:08:18Z

0

The full general case can be achieved with some help from numpy triu_indices.

The following function returns a distribution over the matrices:

def LKJCorrMat(name,n,eta=1):
    c_arr = pm.LKJCorr(name+'_raw',eta=eta,n=n)
    tri = at.zeros( (n,n) )
    tri = at.subtensor.set_subtensor(tri[np.triu_indices(n,1)],c_arr)
    return pm.Deterministic(name,tri + tri.T + at.diag(np.ones(n)))

edited Jan 9, 2024 at 14:08

answered Jan 9, 2024 at 12:15

velochy

4433 silver badges14 bronze badges

Collectives™ on Stack Overflow

How to construct correlation matrix from pymc.LKJCorr?

2 Answers 2

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related