Python: 根据这些分组将一个坐标进行分组并对另一个坐标进行平均

Question

我有两个向量rev_count和stars。它们的元素形成一对（假设rev_count是x坐标，stars是y坐标）。

我想按rev_count分组数据，然后计算一个rev_count bin中stars的平均值（我想沿x轴分组，并计算该分组中y坐标的平均值）。

这是我尝试使用的代码（受到我的Matlab背景启发）：

import matplotlib.pyplot as plt
import numpy

binwidth = numpy.max(rev_count)/10
revbin = range(0, numpy.max(rev_count), binwidth)
revbinnedstars = [None]*len(revbin)

for i in range(0, len(revbin)-1):
    revbinnedstars[i] = numpy.mean(stars[numpy.argwhere((revbin[i]-binwidth/2) < rev_count < (revbin[i]+binwidth/2))])

print('Plotting binned stars with count')
plt.figure(3)
plt.plot(revbin, revbinnedstars, '.')
plt.show()

然而，这种方法似乎非常缓慢/低效。在Python中有更自然的方法吗？

- Ilya

Answer 1

Scipy有一个用于此目的的函数：

from scipy.stats import binned_statistic

revbinnedstars, edges, _ = binned_statistic(rev_count, stars, 'mean', bins=10)
revbin = edges[:-1]

如果您不想使用scipy，还可以使用numpy中的直方图函数：

sums, edges = numpy.histogram(rev_count, bins=10, weights=stars)
counts, _ = numpy.histogram(rev_count, bins=10)
revbinnedstars = sums / counts

- user2379410

---

会尝试，看起来很有前途，我已经在代码的另一部分中使用了scipy。 - Ilya

Answer 2

我猜您正在使用Python 2，但如果不是，计算步长时应将除法更改为//（地板除法），否则numpy将无法将浮点数解释为步长而感到烦恼。

binwidth = numpy.max(rev_count)//10 # Changed this to floor division
revbin = range(0, numpy.max(rev_count), binwidth)
revbinnedstars = [None]*len(revbin)

for i in range(0, len(revbin)-1):
    # I actually don't know what you wanted to do but I guess you wanted the
    # "logical and" combination in that bin (you don't need to use np.where here)
    # You can put that all in one statement but it gets crowded so I'll split it:
    index1 = revbin[i]-binwidth/2 < rev_count
    index2 = rev_count < revbin[i]+binwidth/2)
    revbinnedstars[i] = numpy.mean(stars[np.logical_and(index1, index2)])

至少这个方法可以工作并给出正确的结果。如果您有大型数据集并且想要超过10个bin，那么这将非常低效。

一个非常重要的收获:

• 如果要索引数组，请不要使用np.argwhere。该结果仅应为“人类可读”。如果您确实想要坐标，请使用np.where。它可以用作索引，但如果输入是多维的，则不太容易阅读。

numpy documentation在这一点上支持我：

argwhere的输出不适合用于索引数组。为此，请改用where(a)。

这也是您的代码如此缓慢的原因。它试图做一些您不想要的事情，并且在内存和CPU使用方面可能非常昂贵。而且还无法给您正确的结果。

我在这里所做的被称为布尔掩码。它比np.where(condition)更短，并且少了一个计算步骤。

---

完全矢量化的方法可以通过定义一个网格来使用，该网格知道哪些星星在哪个箱子中：

bins = 10
binwidth = numpy.max(rev_count)//bins
revbin = np.arange(0, np.max(rev_count)+binwidth+1, binwidth)

更好的定义区间的方法是，要注意将最大值加一，因为你想要包括它，并且将区间数加一，因为你关心的是区间的起始和结束点而不是区间的中心：

number_of_bins = 10
revbin = np.linspace(np.min(rev_count), np.max(rev_count)+1, number_of_bins+1)

然后你可以设置网格：

grid = np.logical_and(rev_count[None, :] >= revbin[:-1, None], rev_count[None, :] < revbin[1:, None])

网格大小为bins x rev_count（由于广播的缘故，我增加了每个数组的维度，但并不相同）。这本质上检查一个点是否大于下限范围并且小于上限范围（因此使用了[:-1]和[1:]索引）。这是在多维中完成的，其中计数在第二维（numpy轴=1）中，而bin在第一维（numpy轴=0）中。

因此，我们可以通过将这些与该网格相乘来获取星星在适当bin中的Y坐标。

stars * grid

为了计算平均值，我们需要该区间坐标的总和并将其除以该区间中星星的数量（区间沿axis=1，不在该区间内的星星在此轴上仅具有零值）。

revbinnedstars = np.sum(stars * grid, axis=1) / np.sum(grid, axis=1)

我其实不知道那样是否更有效率。它会占用更多的内存，但或许在 CPU 上会稍微节省一些。

- MSeifert

---

这是Python 3，numpy没有抱怨，但我会改为floor division。我没有意识到Python支持布尔掩码，我现在会尝试。虽然代码仍然很慢，但我会在第一个执行完成后尝试您的第二种方法。谢谢你的帮助！编辑：哦，我读到a<x<b语法更喜欢用于链接比较（而不是使用and），我应该使用numpy的logical_and吗？ - Ilya

@Ilya - 第二种方法仍然存在一个小错误。我已经更新了答案。根据您的样本大小和箱数，这些方法在执行时间和内存使用方面完全不同。您有这些大小的一些数字吗？ - MSeifert

有数十亿行代码。我会尝试另一篇帖子中提到的 scipy.stats.binned_statistics 函数。 - Ilya

Answer 3

我用于将(x,y)数据分组并确定这些分组中的平均值等汇总统计信息的函数是基于scipy.stats.statistic()函数的。由于我经常使用它，因此我为其编写了一个包装器。您可能会发现这很有用...

def binXY(x,y,statistic='mean',xbins=10,xrange=None):
    """
    Finds statistical value of x and y values in each x bin. 
    Returns the same type of statistic for both x and y.
    See scipy.stats.binned_statistic() for options.
    
    Parameters
    ----------
    x : array
        x values.
    y : array
        y values.
    statistic : string or callable, optional
        See documentation for scipy.stats.binned_statistic(). Default is mean.
    xbins : int or sequence of scalars, optional
        If xbins is an integer, it is the number of equal bins within xrange.
        If xbins is an array, then it is the location of xbin edges, similar
        to definitions used by np.histogram. Default is 10 bins.
        All but the last (righthand-most) bin is half-open. In other words, if 
        bins is [1, 2, 3, 4], then the first bin is [1, 2) (including 1, but 
        excluding 2) and the second [2, 3). The last bin, however, is [3, 4], 
        which includes 4.    
        
    xrange : (float, float) or [(float, float)], optional
        The lower and upper range of the bins. If not provided, range is 
        simply (x.min(), x.max()). Values outside the range are ignored.
    
    Returns
    -------
    x_stat : array
        The x statistic (e.g. mean) in each bin. 
    y_stat : array
        The y statistic (e.g. mean) in each bin.       
    n : array of dtype int
        The count of y values in each bin.
        """
    x_stat, xbin_edges, binnumber = stats.binned_statistic(x, x, 
                                 statistic=statistic, bins=xbins, range=xrange)
    
    y_stat, xbin_edges, binnumber = stats.binned_statistic(x, y, 
                                 statistic=statistic, bins=xbins, range=xrange)
    
    n, xbin_edges, binnumber = stats.binned_statistic(x, y, 
                                 statistic='count', bins=xbins, range=xrange)
            
    return x_stat, y_stat, n

- Rupert