Numpy协方差矩阵实现

Question

对于一个项目，我正在实现CSP算法，其中我需要一个协方差矩阵。在多个文献来源中，我发现了计算协方差矩阵的下列表达式：(X乘以X^{T}) /跟踪(X乘以X^{T})，用于我实现的矩阵X：

covariance = np.matmul(A, np.transpose(A)) / (np.matmul(np.dot(A, np.transpose(A))))

从数学的角度来看，这在我的计算中是有意义的，但这个表达式给出的结果与标准的numpy实现np.cov(A)不一样。我的问题是，我找到的表达式是否实际上是协方差矩阵的有效计算？

https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0125039#pone.0125039.ref006给出了截面材料和方法下的数学。

Answer 1

你的公式有很多问题。让我一步一步地做

首先定义一个包含两个变量和五个观测值的矩阵A。使用随机数进行演示

np.random.seed(23423)
n_obs = 5
n_var = 2
A = np.random.normal(0.0,1.0, size=(n_var,n_obs))
A

输出

array([[ 1.51978937,  0.56200696, -0.17770035,  0.01628744,  2.45856104],
       [-1.25514739, -0.17201069, -0.42486057, -0.09215201, -0.21674946]])

第二，我们需要对变量进行中心化(减去它们的平均值)。A0将是A的中心版本

A0 = A - A.mean(axis=1).reshape(n_var,-1)
A0

输出

array([[ 0.64400048, -0.31378193, -1.05348924, -0.85950145,  1.58277214],
       [-0.82296337,  0.26017334,  0.00732345,  0.34003201,  0.21543456]])

现在，我们可以使用公式的一个版本计算协方差(只有分子，但除以观测值-1)。

(A0 @ A0.T)/(n_obs-1)

输出

array([[ 1.21673642, -0.14265396],
       [-0.14265396,  0.22676158]])

让我们将其与np.cov进行比较

np.cov(A)

输出

array([[ 1.21673642, -0.14265396],
       [-0.14265396,  0.22676158]])

两场比赛就像承诺的那样。

您的原始公式看起来更像是相关性而不是协方差；但它也存在多个问题(只有一个参数到matmul，需要两个参数，缺少对中的步骤等等)。