在没有循环的情况下计算“更大”的出现次数

Question

假设我们有一个数组A，如果形状为(100,) 和B，形状为(10,)。两者都包含 [0,1] 中的值。

我们如何使A 中的元素计数大于B 中的每个值？我期望形状为(10,)，其中第一个元素是“A 中有多少大于B[0]”，第二个是“A 中有多少大于B[1]”，等等。 .

不使用循环。

我尝试了以下方法，但没有成功：

import numpy as np
import numpy.random as rdm

A = rdm.rand(100)
B = np.linspace(0,1,10)
def occ(z: float) ->float:
     return np.count_nonzero(A > z)
occ(B)

Python 不会将我的函数用作 B 上的标量函数，这就是我得到的原因：

operands could not be broadcast together with shapes (10,) (100,) 

我也尝试过np.greater，但我遇到了同样的问题......

Answer 1

缓慢但简单

如果您不理解错误消息，它会很神秘，但它会告诉您该怎么做。数组维度是broadcast，从右边缘开始排列。如果您将操作分成两部分，这将特别有用：

1. 创建一个(100, 10) 掩码，显示A 的哪些元素大于B 的哪些元素：

    mask = A[:, None] > B
   

2. 将上一次运算的结果沿A对应的轴求和：

    result = np.count_nonzero(mask, axis=0)
   

   或

    result = np.sum(mask, axis=0)
   

这可以写成单行：

(A[:, None] > B).sum(0)

或

np.count_nonzero(A[:, None] > B, axis=0)

您可以切换尺寸并将B放在第一个轴上以获得相同的结果：

(A > B[:, None]).sum(1)

快速而优雅

采用完全不同（但可能更有效）的方法，您可以使用np.searchsorted：

A.sort()
result = A.size - np.searchsorted(A, B)

默认情况下，searchsorted 返回将 B 的每个元素插入到 A 处的左索引。这几乎可以立即告诉您A 有多少元素大于此值。

基准测试

这里，算法标记如下：

• B0: (A[:, None] > B).sum(0)
• B1: (A > B[:, None]).sum(1)
• HH: np.cumsum(np.histogram(A, bins=B)[0][::-1])[::-1]
• SS: A.sort(); A.size - np.searchsorted(A, B)

+--------+--------+----------------------------------------+
| A.size | B.size |        Time (B0 / B1 / HH / SS)        |
+--------+--------+----------------------------------------+
|    100 |     10 |  20.9 µs / 15.7 µs / 68.3 µs / 8.87 µs |
+--------+--------+----------------------------------------+
|   1000 |     10 |   118 µs / 57.2 µs /  139 µs / 17.8 µs |
+--------+--------+----------------------------------------+
|  10000 |     10 |   987 µs /  288 µs / 1.23 ms /  131 µs |
+--------+--------+----------------------------------------+
| 100000 |     10 |  9.48 ms / 2.77 ms / 13.4 ms / 1.42 ms |
+--------+--------+----------------------------------------+
|    100 |    100 |  70.7 µs / 63.8 µs /   71 µs / 11.4 µs |
+--------+--------+----------------------------------------+
|   1000 |    100 |   518 µs /  388 µs /  148 µs / 21.6 µs |
+--------+--------+----------------------------------------+
|  10000 |    100 |  4.91 ms / 2.67 ms / 1.22 ms /  137 µs |
+--------+--------+----------------------------------------+
| 100000 |    100 |  57.4 ms / 35.6 ms / 13.5 ms / 1.42 ms |
+--------+--------+----------------------------------------+

内存布局很重要。 B1 总是比 B0 快。发生这种情况是因为对连续（缓存）元素（沿 C 顺序中的最后一个轴）求和总是比跳过行以获得下一个元素要快。广播对于B 的小值表现良好。请记住，B0 和B1 的时间和空间复杂度都是O(A.size * B.size)。两种直方图解决方案的复杂度应该在O(A.size * log(A.size)) 左右，但SS 的实现效率要比HH 高得多，因为它可以对数据进行更多假设。

Answer 2

我认为你可以使用np.histogram 来完成这项工作

A = rdm.rand(100)
B = np.linspace(0,1,10)
np.histogram(A, bins=B)[0]

给出输出

array([10,  9,  8, 11,  9, 14, 10, 12, 17])

B[9] 将始终为空，因为没有大于 1 的值。

并向后计算 cumsum

np.cumsum(np.histogram(A, bins=B)[0][::-1])[::-1]

输出

array([100,  90,  81,  73,  62,  53,  39,  29,  17])

Answer 3

np.sum(A>B.reshape((-1,1)), axis=1)

说明

需要了解broadcasting 并为此重塑。通过将 B 重新整形为形状 (len(B), 1)，它可以与 A 一起广播以生成包含所有比较的形状为 (len(B), len(A)) 的数组。然后对轴 1（沿 A）求和。

换句话说，A < B 不起作用，因为A 有 100 个条目，而 B 有 10 个。如果您阅读 broadcasting rules，您会看到 numpy 将从最后一个维度开始，如果它们是大小相同，则可以一对一比较。如果两者之一是1，则该维度被拉伸或“复制”以匹配另一个。如果它们不相等且都不等于1，则失败。

用一个更短的例子：

A = np.array([0.5112744 , 0.21696187, 0.14710105, 0.98581087, 0.50053359,
              0.54954654, 0.81217522, 0.50009166, 0.42990167, 0.56078499])
B = array([0.25, 0.5 , 0.75])

(A>B.reshape((-1,1))) 的转置（为了便于阅读）

np.array([[ True,  True, False],
          [False, False, False],
          [False, False, False],
          [ True,  True,  True],
          [ True,  True, False],
          [ True,  True, False],
          [ True,  True,  True],
          [ True,  True, False],
          [ True, False, False],
          [ True,  True, False]])

而np.sum(A>B.reshape((-1,1)), axis=1) 是

array([8, 7, 2])