Numpy找到两个二维ndarray的协方差答案

【问题标题】：Numpy find covariance of two 2-dimensional ndarrayNumpy找到两个二维ndarray的协方差
【发布时间】：2018-09-06 07:37:41
【问题描述】：

我是 numpy 的新手，遇到了这个问题。我有两个二维 numpy 数组，例如

x = numpy.random.random((10, 5))
y = numpy.random.random((10, 5))

我想使用 numpy cov 函数逐行查找这两个 ndarray 的协方差。即，对于上面的示例，输出数组应由 10 个元素组成，每个元素表示 ndarrays 的相应行的协方差。我知道我可以通过遍历行并找到两个一维数组的协方差来做到这一点，但这不是 pythonic。

Edit1：两个数组的协方差表示0, 1索引处的元素。

Edit2：目前这是我的实现

s = numpy.empty((x.shape[0], 1))
for i in range(x.shape[0]):
    s[i] = numpy.cov(x[i], y[i])[0][1]

【问题讨论】：

标签： python numpy

【解决方案1】：

使用协方差的定义：E(XY) - E(X)E(Y)。

import numpy as np
x = np.random.random((10, 5))
y = np.random.random((10, 5))

n = x.shape[1]
cov_bias = np.mean(x * y, axis=1) - np.mean(x, axis=1) * np.mean(y, axis=1))
cov_bias * n / (n-1)

注意cov_bias对应numpy.cov(bias=True)的结果。

【讨论】：

这是最快的方法，我见过。
确实不错！您可以进一步优化它：np.einsum('ij,ij->i',x,y)/x.shape[1] 替换 np.mean(x * y, axis=1)。

【解决方案2】：

这是使用covariance 的定义并受corr2_coeff_rowwise 启发的一个-

def covariance_rowwise(A,B):
    # Rowwise mean of input arrays & subtract from input arrays themeselves
    A_mA = A - A.mean(-1, keepdims=True)
    B_mB = B - B.mean(-1, keepdims=True)

    # Finally get covariance
    N = A.shape[1]
    return np.einsum('ij,ij->i',A_mA,B_mB)/(N-1)

示例运行 -

In [66]: np.random.seed(0)
    ...: x = np.random.random((10, 5))
    ...: y = np.random.random((10, 5))

In [67]: s = np.empty((x.shape[0]))
    ...: for i in range(x.shape[0]):
    ...:     s[i] = np.cov(x[i], y[i])[0][1]

In [68]: np.allclose(covariance_rowwise(x,y),s)
Out[68]: True

【讨论】：

【解决方案3】：

这可行，但我不确定对于较大的矩阵 x 和 y 是否更快，调用 numpy.cov(x, y) 会计算出许多我们用 numpy.diag 丢弃的条目：

x = numpy.random.random((10, 5))
y = numpy.random.random((10, 5))

# with loop
for (xi, yi) in zip(x, y):
    print(numpy.cov(xi, yi)[0][1])

# vectorized
cov_mat = numpy.cov(x, y)
covariances = numpy.diag(cov_mat, x.shape[0])
print(covariances)

我还为n x n 大小的方阵做了一些计时：

import time
import numpy

def run(n):
    x = numpy.random.random((n, n))
    y = numpy.random.random((n, n))

    started = time.time()
    for (xi, yi) in zip(x, y):
        numpy.cov(xi, yi)[0][1]

    needed_loop = time.time() - started

    started = time.time()
    cov_mat = numpy.cov(x, y)
    covariances = numpy.diag(cov_mat, x.shape[0])
    needed_vectorized = time.time() - started
    print(
        f"n={n:4d} needed_loop={needed_loop:.3f} s "
        f"needed_vectorized={needed_vectorized:.3f} s"
    )

for n in (100, 200, 500, 600, 700, 1000, 2000, 3000):
    run(n)

我的慢速 MacBook Air 上的输出是

n= 100 needed_loop=0.006 s needed_vectorized=0.001 s
n= 200 needed_loop=0.011 s needed_vectorized=0.003 s
n= 500 needed_loop=0.033 s needed_vectorized=0.023 s
n= 600 needed_loop=0.041 s needed_vectorized=0.039 s
n= 700 needed_loop=0.043 s needed_vectorized=0.049 s
n=1000 needed_loop=0.061 s needed_vectorized=0.130 s
n=2000 needed_loop=0.137 s needed_vectorized=0.742 s
n=3000 needed_loop=0.224 s needed_vectorized=2.264 s

所以盈亏平衡点大约在n=600

【讨论】：

【解决方案4】：

选取cov(x,y)的对角向量并展开dims：

numpy.expand_dims(numpy.diag(numpy.cov(x,y),x.shape[0]),1)

【讨论】：