计数不带循环的“较大”事件

Question

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

如何使A中的元素计数大于B中的每个值？我期望一个(10,)的形状，其中第一个元素是“多少in A大于B[0]”，第二个是“多少in 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

我想你可以用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.

然后向后计算总和

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

输出

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

Answer 2

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

解释

需要了解broadcasting并为此进行重塑。通过将B重塑为shape（len（B），1），它可以与A一起广播以生成包含所有比较的shape（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]])

而且{}是

array([8, 7, 2])

Answer 3

缓慢但简单

如果您不理解错误消息，则错误消息很神秘，但它会告诉您该怎么做。数组维度从右边缘开始排列在一起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)

快速优雅

采用完全不同（但可能更有效）的方法，您可以使用^{}：

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

默认情况下，searchsorted返回左索引，B的每个元素将被插入Aat。这会立即告诉您A中有多少元素大于这个值

基准

此处，ALGO标记如下：

• 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高很多，因为它可以承担更多的数据方面的事情