获得没有零值的矩阵方差 numpy

get variance of matrix without zero values numpy

如何计算没有零元素的方差?

例如

np.var([[1, 1], [1, 2]], axis=1) -> [0, 0.25]

我需要:

var([[1, 1, 0], [1, 2, 0]], axis=1) -> [0, 0.25]

是您要找的吗?您可以过滤掉所有值为 0(或至少一个值不为 0)的列。

m = np.array([[1, 1, 0], [1, 2, 0]])
np.var(m[:, np.any(m != 0, axis=0)], axis=1)

# Output
array([0.  , 0.25])

V1

您可以使用掩码数组:

data = np.array([[1, 1, 0], [1, 2, 0]])
np.ma.array(data, mask=(data == 0)).var(axis=1)

结果是

masked_array(data=[0.  , 0.25],
             mask=False,
       fill_value=1e+20)

原始 numpy 数组是生成的屏蔽数组的 data 属性:

>>> np.ma.array(data, mask=(data == 0)).var(axis=1).data
array([0.  , 0.25])

V2

如果没有掩码数组,在每行中删除可变数量的元素的操作有点棘手。根据公式 sum(x**2) / N - (sum(x) / N)**2 和 ufuncs.

的部分缩减来实现方差会更简单

首先我们需要找到拆分索引和段长度。在一般情况下,看起来像

lens = np.count_nonzero(data, axis=1)
inds = np.r_[0, lens[:-1].cumsum()]

现在您可以对混淆后的屏蔽数据进行操作了:

mdata = data[data != 0]
mdata2 = mdata**2
var = np.add.reduceat(mdata2, inds) / lens - (np.add.reduceat(mdata, inds) / lens)**2

这为您提供了与 var 相同的结果(顺便说一句,可能比屏蔽版本更有效):

array([0.  , 0.25])

V3

var 函数似乎使用了更传统的公式 (x - x.mean()).mean()。您可以使用上面的数量来实现它,只需多做一点工作:

means = (np.add.reduceat(mdata, inds) / lens).repeat(lens)
var = np.add.reduceat((mdata - means)**2, inds) / lens

比较

这是两种方法的快速基准:

def nzvar_v1(data):
    return np.ma.array(data, mask=(data == 0)).var(axis=1).data

def nzvar_v2(data):
    lens = np.count_nonzero(data, axis=1)
    inds = np.r_[0, lens[:-1].cumsum()]
    mdata = data[data != 0]
    return np.add.reduceat(mdata**2, inds) / lens - (np.add.reduceat(mdata, inds) / lens)**2

def nzvar_v3(data):
    lens = np.count_nonzero(data, axis=1)
    inds = np.r_[0, lens[:-1].cumsum()]
    mdata = data[data != 0]
    return np.add.reduceat((mdata - (np.add.reduceat(mdata, inds) / lens).repeat(lens))**2, inds) / lens

np.random.seed(100)
data = np.random.randint(10, size=(1000, 1000))

%timeit nzvar_v1(data)
18.3 ms ± 278 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit nzvar_v2(data)
5.89 ms ± 69.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit nzvar_v3(data)
11.8 ms ± 62.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

因此,对于大型数据集,第二种方法虽然需要更多代码,但似乎比屏蔽数组快约 3 倍,比使用传统公式快约 2 倍。