在 numpy 数组中查找模式的最有效方法

Question

我有一个包含整数（正数或负数）的二维数组。每行代表特定空间站点随时间变化的值，而每列代表给定时间的各种空间站点的值。

所以如果数组是这样的：

 1 3 4 2 2 7
5 2 2 1 4 1
3 3 2 2 1 1

结果应该是

1 3 2 2 2 1

请注意，当模式有多个值时，可以将任何一个（随机选择）设置为模式。

我可以一次遍历列查找模式，但我希望 numpy 可能有一些内置函数来做到这一点。或者如果有一个技巧可以在不循环的情况下有效地找到它。

原文由 Nik 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 1

检查 scipy.stats.mode() （灵感来自@tom10 的评论）：

 import numpy as np
from scipy import stats

a = np.array([[1, 3, 4, 2, 2, 7],
              [5, 2, 2, 1, 4, 1],
              [3, 3, 2, 2, 1, 1]])

m = stats.mode(a)
print(m)

输出：

 ModeResult(mode=array([[1, 3, 2, 2, 1, 1]]), count=array([[1, 2, 2, 2, 1, 2]]))

如您所见，它返回模式和计数。您可以直接通过 m[0] 选择模式：

 print(m[0])

输出：

 [[1 3 2 2 1 1]]

原文由 fgb 发布，翻译遵循 CC BY-SA 3.0 许可协议

Answer 2

更新

scipy.stats.mode 函数自这篇文章以来已得到显着优化，并且将是推荐的方法

旧答案

这是一个棘手的问题，因为没有多少东西可以计算沿轴的模式。 The solution is straight forward for 1-D arrays, where numpy.bincount is handy, along with numpy.unique with the return_counts arg as True .我看到的最常见的 n 维函数是 scipy.stats.mode，尽管它非常慢——尤其是对于具有许多唯一值的大型数组。作为解决方案，我开发了这个功能，并大量使用它：

 import numpy

def mode(ndarray, axis=0):
    # Check inputs
    ndarray = numpy.asarray(ndarray)
    ndim = ndarray.ndim
    if ndarray.size == 1:
        return (ndarray[0], 1)
    elif ndarray.size == 0:
        raise Exception('Cannot compute mode on empty array')
    try:
        axis = range(ndarray.ndim)[axis]
    except:
        raise Exception('Axis "{}" incompatible with the {}-dimension array'.format(axis, ndim))

    # If array is 1-D and numpy version is > 1.9 numpy.unique will suffice
    if all([ndim == 1,
            int(numpy.__version__.split('.')[0]) >= 1,
            int(numpy.__version__.split('.')[1]) >= 9]):
        modals, counts = numpy.unique(ndarray, return_counts=True)
        index = numpy.argmax(counts)
        return modals[index], counts[index]

    # Sort array
    sort = numpy.sort(ndarray, axis=axis)
    # Create array to transpose along the axis and get padding shape
    transpose = numpy.roll(numpy.arange(ndim)[::-1], axis)
    shape = list(sort.shape)
    shape[axis] = 1
    # Create a boolean array along strides of unique values
    strides = numpy.concatenate([numpy.zeros(shape=shape, dtype='bool'),
                                 numpy.diff(sort, axis=axis) == 0,
                                 numpy.zeros(shape=shape, dtype='bool')],
                                axis=axis).transpose(transpose).ravel()
    # Count the stride lengths
    counts = numpy.cumsum(strides)
    counts[~strides] = numpy.concatenate([[0], numpy.diff(counts[~strides])])
    counts[strides] = 0
    # Get shape of padded counts and slice to return to the original shape
    shape = numpy.array(sort.shape)
    shape[axis] += 1
    shape = shape[transpose]
    slices = [slice(None)] * ndim
    slices[axis] = slice(1, None)
    # Reshape and compute final counts
    counts = counts.reshape(shape).transpose(transpose)[slices] + 1

    # Find maximum counts and return modals/counts
    slices = [slice(None, i) for i in sort.shape]
    del slices[axis]
    index = numpy.ogrid[slices]
    index.insert(axis, numpy.argmax(counts, axis=axis))
    return sort[index], counts[index]

结果：

 In [2]: a = numpy.array([[1, 3, 4, 2, 2, 7],
                         [5, 2, 2, 1, 4, 1],
                         [3, 3, 2, 2, 1, 1]])

In [3]: mode(a)
Out[3]: (array([1, 3, 2, 2, 1, 1]), array([1, 2, 2, 2, 1, 2]))

一些基准：

 In [4]: import scipy.stats

In [5]: a = numpy.random.randint(1,10,(1000,1000))

In [6]: %timeit scipy.stats.mode(a)
10 loops, best of 3: 41.6 ms per loop

In [7]: %timeit mode(a)
10 loops, best of 3: 46.7 ms per loop

In [8]: a = numpy.random.randint(1,500,(1000,1000))

In [9]: %timeit scipy.stats.mode(a)
1 loops, best of 3: 1.01 s per loop

In [10]: %timeit mode(a)
10 loops, best of 3: 80 ms per loop

In [11]: a = numpy.random.random((200,200))

In [12]: %timeit scipy.stats.mode(a)
1 loops, best of 3: 3.26 s per loop

In [13]: %timeit mode(a)
1000 loops, best of 3: 1.75 ms per loop

编辑：提供更多背景知识并修改方法以提高内存效率

原文由 Devin Cairns 发布，翻译遵循 CC BY-SA 4.0 许可协议