如何快速获取numpy数组的众数？

Question

我需要找出从hdf5文件中读取的NumPy数组的众数。该NumPy数组是一维的，包含浮点值。

my_array=f1[ds_name].value    
mod_value=scipy.stats.mode(my_array)

我的数组是一维的，包含大约100万个值。我的脚本需要大约15分钟才能返回众数值。有没有什么方法可以使这个过程更快？

另一个问题是为什么scipy.stats.median(my_array)无法工作，而mode可以？

AttributeError: module 'scipy.stats' has no attribute 'median'

- Heli

---

听起来是IO限制，因为我认为剩下的代码已经很优化了。所以请检查您的hdf缓冲区、压缩等。另外：scipy.stats没有一个名为median的函数，可以通过阅读文档轻松确认。您可以使用numpy的中位数函数。 - sascha

@sascha 这个文件在0.02秒内被读取。15分钟的时间花费在计算这行代码中的众数 "scipy.stats.median(my_array)" 上。 - Heli

也许你应该展示更多的代码，因为你的计时结果与给出答案中显示的这些合成示例差距很大（这些示例还表明我错了；可以实现不那么难的加速）。 - sascha

@sascha，请查看完整代码的回复。如果您需要我的输入文件进行测试，请告诉我，谢谢。 - Heli

Answer 1

scipy.stats.mode的实现在处理多维数组时使用了Python循环以处理axis参数。对于仅限一维数组的情况，以下简单实现更快：

该函数的实现在处理多维数组时需要使用Python循环来处理axis参数。对于仅限一维数组的情况，下面这个简单的实现更快：

def mode1(x):
    values, counts = np.unique(x, return_counts=True)
    m = counts.argmax()
    return values[m], counts[m]

以下是一个例子。首先，创建一个长度为1000000的整数数组。

In [40]: x = np.random.randint(0, 1000, size=(2, 1000000)).sum(axis=0)

In [41]: x.shape
Out[41]: (1000000,)

检查scipy.stats.mode和mode1是否给出相同的结果。

In [42]: from scipy.stats import mode

In [43]: mode(x)
Out[43]: ModeResult(mode=array([1009]), count=array([1066]))

In [44]: mode1(x)
Out[44]: (1009, 1066)

现在检查性能。

In [45]: %timeit mode(x)
2.91 s ± 18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [46]: %timeit mode1(x)
39.6 ms ± 83.8 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

mode(x)需要2.91秒，而mode1(x)仅需要39.6毫秒。

- Warren Weckesser

---

为了好玩，我尝试使用“Counter”进行计算。运行时间大约是“scs.mode”的一半，但是您的解决方案轻松胜出。 - Brad Solomon

Answer 2

这里是基于排序的一种方法 -

def mode1d(ar_sorted):
    ar_sorted.sort()
    idx = np.flatnonzero(ar_sorted[1:] != ar_sorted[:-1])
    count = np.empty(idx.size+1,dtype=int)
    count[1:-1] = idx[1:] - idx[:-1]
    count[0] = idx[0] + 1
    count[-1] = ar_sorted.size - idx[-1] - 1
    argmax_idx = count.argmax()

    if argmax_idx==len(idx):
        modeval = ar_sorted[-1]
    else:
        modeval = ar_sorted[idx[argmax_idx]]
    modecount = count[argmax_idx]
    return modeval, modecount

请注意，这会改变输入数组并对其进行排序。因此，如果您想保持输入数组未被更改或者不希望输入数组被排序，请传递一个副本。

1百万元素的样本运行 -

In [65]: x = np.random.randint(0, 1000, size=(1000000)).astype(float)

In [66]: from scipy.stats import mode

In [67]: mode(x)
Out[67]: ModeResult(mode=array([ 295.]), count=array([1098]))

In [68]: mode1d(x)
Out[68]: (295.0, 1098)

运行时测试

In [75]: x = np.random.randint(0, 1000, size=(1000000)).astype(float)

# Scipy's mode
In [76]: %timeit mode(x)
1 loop, best of 3: 1.64 s per loop

# @Warren Weckesser's soln
In [77]: %timeit mode1(x)
10 loops, best of 3: 52.7 ms per loop

# Proposed in this post
In [78]: %timeit mode1d(x)
100 loops, best of 3: 12.8 ms per loop

使用副本，mode1d 的时间与 mode1 相当。

- Divakar

Answer 3

我将上面回复中的两个函数mode1和mode1d添加到我的脚本中，并尝试与scipy.stats.mode进行比较。

dir_name="C:/Users/test_mode"
file_name="myfile2.h5"
ds_name="myds"
f_in=os.path.join(dir_name,file_name)

def mode1(x):
    values, counts = np.unique(x, return_counts=True)
    m = counts.argmax()
    return values[m], counts[m]

def mode1d(ar_sorted):
    ar_sorted.sort()
    idx = np.flatnonzero(ar_sorted[1:] != ar_sorted[:-1])
    count = np.empty(idx.size+1,dtype=int)
    count[1:-1] = idx[1:] - idx[:-1]
    count[0] = idx[0] + 1
    count[-1] = ar_sorted.size - idx[-1] - 1
    argmax_idx = count.argmax()

    if argmax_idx==len(idx):
        modeval = ar_sorted[-1]
    else:
        modeval = ar_sorted[idx[argmax_idx]]
    modecount = count[argmax_idx]
    return modeval, modecount


startTime=time.time()
with h5py.File(f_in, "a") as f1:

        myds=f1[ds_name].value
        time1=time.time()
        file_read_time=time1-startTime
        print(str(file_read_time)+"\t"+"s"+"\t"+str((file_read_time)/60)+"\t"+"min")

        print("mode_scipy=")
        mode_scipy=scipy.stats.mode(myds)
        print(mode_scipy)
        time2=time.time()
        mode_scipy_time=time2-time1
        print(str(mode_scipy_time)+"\t"+"s"+"\t"+str((mode_scipy_time)/60)+"\t"+"min")

        print("mode1=")
        mode1=mode1(myds)
        print(mode1)
        time3=time.time()
        mode1_time=time3-time2
        print(str(mode1_time)+"\t"+"s"+"\t"+str((mode1_time)/60)+"\t"+"min")

        print("mode1d=")
        mode1d=mode1d(myds)
        print(mode1d)
        time4=time.time()
        mode1d_time=time4-time3
        print(str(mode1d_time)+"\t"+"s"+"\t"+str((mode1d_time)/60)+"\t"+"min")

运行numpy数组约为1M的脚本的结果是：

mode_scipy= ModeResult(mode=array([ 1.11903353e-06], dtype=float32), count=array([304909])) 938.8368742465973 秒
15.647281237443288 分钟

mode1=(1.1190335e-06, 304909)

0.06500649452209473 秒
0.0010834415753682455 分钟

mode1d=(1.1190335e-06, 304909)

0.06200599670410156 秒
0.0010334332784016928 分钟

- Heli