在x范围内的y值运行中位数

Question

下面是我用两个numpy数组制作的散点图。

散点图示例

我想在这个图上添加一个y值的移动中位数，范围是x值的某个区间。我在Photoshop里做了个示例：

修改后的散点图

具体来说，我需要在x轴上每1个单位的区间内，计算出数据点的中位数，这个区间的值会在多个图之间变化，但我可以手动调整。我很感谢任何能给我指明方向的建议。

Answer 1

我在C#里写过类似的东西。我不太懂Python，所以这里是伪代码：

• 创建一个List，用来存放我们要计算中位数的数据。
• 根据x值对散点图的点进行排序。
• 按x值的顺序逐个处理这些排序后的点。
• 对于每个点，把它的Y值插入到中位数列表中，这样中位数列表就会保持有序。也就是说，要插入Y值时，确保它上面的值都比它小，下面的值都比它大。你可以看看这里的内容：在Python中将值插入列表的特定位置。
• 每次添加Y值后，中位数的值就是当前列表中间位置的值，也就是List(List.Length/2)。

希望这对你有帮助！

Answer 2

你可以创建一个函数，利用 numpy.median() 来计算给定区间的中位数值：

import numpy as np

def medians(x, y, intervals):
    out = []
    for xmin, xmax in intervals:
        mask = (x >= xmin) & (x < xmax)
        out.append(np.median(y[mask]))
    return np.array(out)

然后用这个函数来处理你想要的区间：

import matplotlib.pyplot as plt

intervals = ((18, 19), (19, 20), (20, 21), (21, 22))
centers = [(xmin+xmax)/2. for xmin, xmax in intervals]

plt.plot(centers, medians(x, y, intervals)

Answer 3

这个问题也可以通过 python pandas（Python数据分析库）来高效解决，它提供了方便的数据切割和分析方法。

考虑一下这个例子：

（感谢@Hooked提供的例子，我借用了其中的X和Y数据）

 import pandas as pd
 df = pd.DataFrame({'X' : X, 'Y' : Y})  #we build a dataframe from the data

 data_cut = pd.cut(df.X,bins)           #we cut the data following the bins
 grp = df.groupby(by = data_cut)        #we group the data by the cut

 ret = grp.aggregate(np.median)         #we produce an aggregate representation (median) of each bin

 #plotting

 plt.scatter(df.X,df.Y,color='k',alpha=.2,s=2)
 plt.plot(ret.X,ret.Y,'r--',lw=4,alpha=.8)
 plt.show()

备注：这里红色曲线的x值是按区间计算的中位数（可以使用区间的中点）。

Answer 4

我会使用 np.digitize 来帮你进行分组排序。这样你就可以轻松地应用任何函数，并设置你感兴趣的范围。

import numpy as np
import pylab as plt

N = 2000
total_bins = 10

# Sample data
X = np.random.random(size=N)*10
Y = X**2 + np.random.random(size=N)*X*10

bins = np.linspace(X.min(),X.max(), total_bins)
delta = bins[1]-bins[0]
idx  = np.digitize(X,bins)
running_median = [np.median(Y[idx==k]) for k in range(total_bins)]

plt.scatter(X,Y,color='k',alpha=.2,s=2)
plt.plot(bins-delta/2,running_median,'r--',lw=4,alpha=.8)
plt.axis('tight')
plt.show()

为了展示这种方法的灵活性，我们来加上每个分组的标准差作为误差条：

running_std    = [Y[idx==k].std() for k in range(total_bins)]
plt.errorbar(bins-delta/2,running_median,
              running_std,fmt=None)