如何从 pymc.LKJCorr 构造相关矩阵？

Question

我想使用

pymc.LKJCorr

分布类显式构造一个相关矩阵，但我不相信我对

pymc.expand_packed_triangular

的理解。这是一个最小的工作示例。

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

这是一个示例图：

对我来说没有意义的是，我看到

RCorr

参数之一以 1 为中心。应该不是这样的，看起来这只是统一翻译的

"LKJCorr"

。我担心

pm.expand_packed_triangular

假设有对角线，但实际上并不存在，或者类似的东西。

如何在这个玩具示例中使用

pymc.LKJCorr

的实例重建相关矩阵？

Answer 1

在 2x2 相关矩阵的简单情况下，您可以执行以下操作

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

它的工作原理是创建一个 2x2 数组，其元素都是相关分数，这在这里很好，因为只有一个，然后用 1 填充对角线。

还不知道如何处理n维情况。

Answer 2

在 numpy triu_indices 的帮助下可以实现完整的一般情况。

以下函数返回矩阵的分布：

def LKJCorrMat(name,n,eta=1):
    c_arr = pm.LKJCorr.dist(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)))