Implement random method for LKJCorr

[junpenglao]

using the algorithm in LKJ 2009(vine method based on a C-vine)

• implement test.

[ferrine]

Is it possible to vectorize it?

[junpenglao]

I can not see an easy way to do it, the step below converting the partial correlation to raw correlation is especially tricky.

[junpenglao]

Any suggestion of implementing a test? I am stuck. @ferrine @aseyboldt

[ferrine]

What about comparing kde with theoretical density?

[ferrine]

You'd better create these samples beforehand

[twiecki]

Or you could run a KS-test between true and sampled density.

[junpenglao]

I see - is there any similar test I can reference to in our test? @ferrine @twiecki

[aseyboldt]

This is similar I think: https://github.com/pymc-devs/pymc3/blob/master/pymc3/tests/sampler_fixtures.py#L28

[junpenglao]

Thanks everyone! I think I finally figure it out ;-)

[twiecki]

pep8 wants white spaces around math operators

[twiecki]

white-spaces

[junpenglao]

what i am doing here is slow in a for loop, should I change it to reduce @ColCarroll?

[aseyboldt]

For stuff like this we could also consider to use numba if it is available. (check if it is available and create a no-op replacement if it is not).

[ferrine]

should be static I suppose

[ferrine]

seems like size argument is ignored

[junpenglao]

Yep. I ignored the size as the _random method here can only generate 1 slide of the random matrix.

[junpenglao]

Thanks for the suggestion - I manage to use the size and vectorized the _random method.

[junpenglao]

OK the implementation should be correct as the distribution of the matrix elements match beta() distribution. However I am getting none positive-definite matrix when the dimension n>3:

import numpy as np
from scipy import stats

# cherry picked LKJ random function
def lkj_random(n, eta, size=None):
    beta0 = eta - 1 + n/2
    shape = n * (n-1) // 2
    triu_ind = np.triu_indices(n, 1)
    beta = np.array([beta0 - k/2 for k in triu_ind[0]])
    # partial correlations sampled from beta dist.
    P = np.ones((n, n) + (size,))
    P[triu_ind] = stats.beta.rvs(a=beta, b=beta, size=(size,) + (shape,)).T
    # scale partial correlation matrix to [-1, 1]
    P = (P - .5) * 2
    
    for k, i in zip(triu_ind[0], triu_ind[1]):
        p = P[k, i]
        for l in range(k-1, -1, -1):  # convert partial correlation to raw correlation
            p = p * np.sqrt((1 - P[l, i]**2) *
                            (1 - P[l, k]**2)) + P[l, i] * P[l, k]
        P[k, i] = p
        P[i, k] = p

    return np.transpose(P, (2, 0 ,1))

def is_pos_def(A):
    if np.array_equal(A, A.T):
        try:
            np.linalg.cholesky(A)
            return 1
        except np.linalg.linalg.LinAlgError:
            return 0
    else:
        return 0

P = lkj_random(4, 1., 1000)
k=0
for i, p in enumerate(P):
    k+=is_pos_def(p)
print(k)

Thoughts? @aseyboldt

[aseyboldt]

@junpenglao I haven't look at it in detail, but should there really be that many 1. in the partial correlations P? Maybe something like P = np.ones((n, n) + (size,)) / 2?

[junpenglao]

@aseyboldt at the end only the diagonal is 1. Currently LKJ random returns the upper triangular elements (same as the distribution) but I was just doing some test in the code above.

[aseyboldt]

Sorry, you are right.
The partial correlations are just a scaled version of the precision matrix, right? So maybe you could use that to check the result?

[junpenglao]

@aseyboldt not sure I understand what you mean.

I also compare with a Julia implementation, both have the same output. So it might be a numerical stability issue.

[EDIT], trying to figure out whether I have the same output as Julia or Stan.

[aseyboldt]

hm. But shouldn't A_{ij}/np.sqrt(A_ii * A_jj) equal the original partial correlations, were A = \sigma^{-1}? It isn't. And sometimes the eigenvalues are definitely negative, much more so than could be explained by some simple rounding.

[aseyboldt]

I'm using this to compute the partial correlations:

tau = linalg.inv(P[0])

partial = tau.copy()
partial /= np.sqrt(np.diag(tau))[None, :]
partial /= np.sqrt(np.diag(tau))[:, None]

(And just added a print statement in the function to get the original values)

[junpenglao]

I will also compare with the implementation in R from @rmcelreath https://github.com/rmcelreath/rethinking/blob/master/R/distributions.r#L165-L184, I cannot really figure it out in Stan and Julia

[aseyboldt]

Sorry, my last comment was wrong, I looked at the wrong array. I just deleted it.

[junpenglao]

Turns out the R implementation from @rmcelreath is the most stable - test pass locally and samples are mostly positive definite (failed sometimes when n is large and eta<<1, but much better than the previous implementation nonetheless).

[junpenglao]

@aseyboldt The current random method could potentially be extended for generating LKJCholskyCov as it produces a triangular matrix as an intermedia step.