2

我有一个二维数组 x，每一行都有不同数量的 nan 值：

array([[   nan, -0.355, -0.036, ...,    nan,    nan],
       [   nan, -0.341, -0.047, ...,    nan,  0.654],
       [  .016, -1.147, -0.667, ...,    nan,    nan],
       ..., 
       [   nan,  0.294, -0.235, ...,    0.65,   nan]])

给定这个数组，对于每一行，我想计算前 25 个百分位数内所有值的平均值。我正在执行以下操作：

limit = np.nanpercentile(x, 25, axis=1) # output 1D array
ans = np.nanmean(x * (x < limit[:,None]), axis=1)

但这给出了错误的结果——特别是计数（np.nansum/np.nanmean）无论我选择什么百分位数都保持不变，因为比较在不正确的地方产生零，并被视为平均值的有效值。我不能简单地使用x[x>limit[:,None]]，因为它给出了一个 1D 数组，我需要一个 2D 结果。

我通过以下方式解决了它：

f = x.copy()
f[f > limit[:,None]] = np.nan
ans = np.nanmean(f, axis=1)

有更好的方法吗？

4

1 回答 1

2

方法#1：您可以创建一个无效掩码，它NaNs来自原始数组，掩码来自f > limit[:,None]. np.nanmean然后，使用此掩码通过仅考虑有效掩码来执行等效方法masking。使用的好处masks/boolean arrays在于内存，因为它占用的内存比浮动 pt 数组少 8 倍。因此，我们会有这样的实现 -

# Create mask of non-NaNs and thresholded ones
mask = ~np.isnan(x) & (x <= limit[:,None])

# Get the row, col indices. Use the row indices for bin-based summing and
# finally averaging by using those indices to get the group lengths.
r,c = np.where(mask)
out = np.bincount(r,x[mask])/np.bincount(r)

方法#2：我们也np.add.reduceat可以在这里使用 which 是有帮助的，因为 bin 已经根据掩码进行了排序。所以，更高效一点就是这样 -

# Get the valid mask as before
mask = ~np.isnan(x) & (x <= limit[:,None])

# Get valid row count. Use np.add.reduceat to perform grouped summations
# at intervals separated by row indices.
rowc = mask.sum(1)
out = np.add.reduceat(x[mask],np.append(0,rowc[:-1].cumsum()))/rowc

基准测试

函数定义——

def original_app(x, limit):
    f = x.copy()
    f[f > limit[:,None]] = np.nan
    ans = np.nanmean(f, axis=1) 
    return ans

def proposed1_app(x, limit):
    mask = ~np.isnan(x) & (x <= limit[:,None])
    r,c = np.where(mask)
    out = np.bincount(r,x[mask])/np.bincount(r)
    return out

def proposed2_app(x, limit):
    mask = ~np.isnan(x) & (x <= limit[:,None])
    rowc = mask.sum(1)
    out = np.add.reduceat(x[mask],np.append(0,rowc[:-1].cumsum()))/rowc
    return out

计时和验证 -

In [402]: # Setup inputs
     ...: x = np.random.randn(400,500)
     ...: x.ravel()[np.random.randint(0,x.size,x.size//4)] = np.nan # Half as NaNs
     ...: limit = np.nanpercentile(x, 25, axis=1)
     ...: 

In [403]: np.allclose(original_app(x, limit),proposed1_app(x, limit))
Out[403]: True

In [404]: np.allclose(original_app(x, limit),proposed2_app(x, limit))
Out[404]: True

In [405]: %timeit original_app(x, limit)
100 loops, best of 3: 5 ms per loop

In [406]: %timeit proposed1_app(x, limit)
100 loops, best of 3: 4.02 ms per loop

In [407]: %timeit proposed2_app(x, limit)
100 loops, best of 3: 2.18 ms per loop
于 2016-11-21T10:11:22.733 回答