CODEWITHSUNDEEP
pythoncorrelationmontecarlopymcprobabilistic-programming
Galen
Galen

Reputation: 1352

How to construct correlation matrix from pymc.LKJCorr?

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:

enter image description here

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?

Upvotes: 0

Views: 146

Answers (2)

velochy
velochy

Reputation: 425

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)))

Upvotes: 0

Galen
Galen

Reputation: 1352

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.

Upvotes: 0

Related Questions