【问题标题】:weighted covariance matrix in numpynumpy中的加权协方差矩阵
【发布时间】:2016-07-19 02:21:03
【问题描述】:

我想计算 n 测量 p 的协方差 C,其中每个单独的数量测量都有自己的权重。也就是说,我的权重数组W 与我的数量数组Q 具有相同的形状(n by p)。原生 np.cov() 函数仅支持赋予单个测量值的权重(即长度为 n 的向量)。

我可以通过p 矩阵初始化p 并进行迭代,但是如果p 很大,那么这是一个非常缓慢的过程。

由于已知Q 的每个数量(Q 的列)的均值为零,因此C 的每个元素的显式公式为

C[i,j] = np.sum(
    Q[:, i] * Q[:, j] * W[:, i] * W[:, j]) / np.sum(W[:, i] * W[:, j])

如果我将分子重新排列为Q[:, i] * W[:, i] * Q[:, j] * W[:, j],似乎我应该能够将Q * W 的列相乘和相加,然后类似地做分母(使用W * W 除外)。

有没有办法用np.einsum() 做到这一点?

为了测试,让我们定义以下内容:

C = array([[ 1.  ,  0.1 ,  0.2 ], # set this beforehand, to test whether 
           [ 0.1 ,  0.5 ,  0.15], # we get the correct result
           [ 0.2 ,  0.15,  0.75]])

Q = array([[-0.6084634 ,  0.16656143, -1.04490324],
           [-1.51164337, -0.96403094, -2.37051952],
           [-0.32781346, -0.19616374, -1.32591578],
           [-0.88371729,  0.20877833, -0.52074272],
           [-0.67987913, -0.84458226,  0.02897935],
           [-2.01924756, -0.51877396, -0.68483981],
           [ 1.64600477,  0.67620595,  1.24559591],
           [ 0.82554885,  0.14884613, -0.15211434],
           [-0.88119527,  0.11663335, -0.31522598],
           [-0.14830668,  1.26906561, -0.49686309]])

W = array([[ 1.01133857,  0.91962164,  1.01897898],
           [ 1.09467975,  0.91191381,  0.90150961],
           [ 0.96334661,  1.00759046,  1.01638749],
           [ 1.04827001,  0.95861001,  1.01248969],
           [ 0.91572506,  1.09388218,  1.03616461],
           [ 0.9418178 ,  1.07210878,  0.90431879],
           [ 1.0093642 ,  1.00408472,  1.07570172],
           [ 0.92203074,  1.00022631,  1.09705542],
           [ 0.99775598,  0.01000000,  0.94996408],
           [ 1.02996389,  1.01224303,  1.00331465]])

【问题讨论】:

    标签: python arrays numpy covariance


    【解决方案1】:

    您可以使用非常有效的矩阵乘法与np.dot -

    QW = Q*W
    C = QW.T.dot(QW)/W.T.dot(W)
    

    【讨论】:

    • 你说得对,这个时钟比我的 100x10 数组方法快约 3 倍,而且这个因子随着 QW 的大小而增长。
    【解决方案2】:

    经过更多的实验,我发现以下工作:

    A = np.einsum('ki,kj->ij', Q*W, Q*W)
    B = np.einsum('ki,kj->ij', W, W)
    C = A/B
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2017-08-27
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多