我发现一个特定的稀疏矩阵构造器可以非常有效地达到预期的结果。它有点晦涩，但我们可以为此目的滥用它。下面的函数可以用几乎与scipy.stats.binned_statistic相同的方式使用，但可以快几个数量级

import numpy as np
from scipy.sparse import csr_matrix

def binned_statistic(x, values, func, nbins, range):
    '''The usage is nearly the same as scipy.stats.binned_statistic''' 

    N = len(values)
    r0, r1 = range

    digitized = (float(nbins)/(r1 - r0)*(x - r0)).astype(int)
    S = csr_matrix((values, [digitized, np.arange(N)]), shape=(nbins, N))

    return [func(group) for group in np.split(S.data, S.indptr[1:-1])]

我避免使用np.digitize，因为它没有利用所有bin宽度相等的事实，因此速度很慢，但我使用的方法可能无法完美地处理所有边缘情况。

我假设在digitize的例子中所做的合并是不能改变的，这是一种方法，你可以一次性地完成排序。

vals = np.random.random(1e4)
nbins = 100
bins = np.linspace(0, 1, nbins+1)
ind = np.digitize(vals, bins)

new_order = argsort(ind)
ind = ind[new_order]
ordered_vals = vals[new_order]
# slower way of calculating first_hit (first version of this post)
# _,first_hit = unique(ind,return_index=True)
# faster way:
first_hit = searchsorted(ind,arange(1,nbins-1))
first_hit.sort()

#example of using the data:
for j in range(nbins-1):
    #I am using a plotting function for your f, to show that they cluster
    plot(ordered_vals[first_hit[j]:first_hit[j+1]],'o')

该图显示了条柱实际上是预期的聚类：

通过先对数组排序，然后使用np.searchsorted，可以将计算时间减半。

vals = np.random.random(1e8)
vals.sort()

nbins = 100
bins = np.linspace(0, 1, nbins+1)
ind = np.digitize(vals, bins)

results = [func(vals[np.searchsorted(ind,j,side='left'):
                     np.searchsorted(ind,j,side='right')])
           for j in range(1,nbins)]

使用1e8作为我的测试用例，我的计算时间从34秒减少到17秒。

一种高效的解决方案是使用numpy_indexed包（免责声明：我是它的作者）：

import numpy_indexed as npi
npi.group_by(ind).split(vals)

2023更新：新增Scipy值_指数

为了记录在案，在我最初的问题八年之后，Scipy 1.10 in January 2023引入了一个新函数scipy.ndimage.value_indices，它完全符合我在问题中提出的要求，文档甚至明确提到他们试图模拟IDL功能
IDL用户注意事项：这提供了与IDL的REVERSE_INDICES选项等效的功能（根据HISTOGRAM函数的IDL文档）。
使用新的Scipy函数，可接受答案中建议的函数的等效函数如下所示

import numpy as np
from scipy.ndimage import value_indices

def binned_statistic(x, values, func, nbins, extent):
    '''The usage is nearly the same as scipy.stats.binned_statistic''' 

    r0, r1 = extent
    digitized = (float(nbins)/(r1 - r0)*(x - r0)).astype(int)

    val_indices = value_indices(digitized)
    keys = val_indices.keys()

    return [func(values[val_indices[j]]) for j in keys]

此函数可按如下方式使用

import numpy as np

# This function func() is just a placemark for a much more complicated function.
# I am aware that my problem could be easily sped up in the specific case of
# of the sum() function, but I am looking for a general solution to the problem.
def func(x):
    y = np.sum(x)
    return y

vals = np.random.random(int(1e8))
nbins = 100
extent = [0, 1]
res = binned_statistic(x, vals, func, nbins, extent)

我根据当前公认的答案计算了新函数的时间，发现它在给定的示例中具有相当的速度，但它是1.7倍 * 慢 *。因此，这不明显应该成为公认的答案，因为效率取决于问题的大小。

Pandas有一个非常快的分组代码（我认为它是用C编写的），所以如果你不介意加载库，你可以这样做：

import pandas as pd

pdata=pd.DataFrame({'vals':vals,'ind':ind})
resultsp = pdata.groupby('ind').sum().values

或者更一般地说：

pdata=pd.DataFrame({'vals':vals,'ind':ind})
resultsp = pdata.groupby('ind').agg(func).values

虽然后者对于标准聚合函数（如sum、mean等）来说速度较慢