我终于想出了一个解决方案,它看起来可能不是最优雅的,但它工作得相当好:
要估计二维分布的分位数,可以使用scipy 函数binned_statistics,它允许将数据装箱
其中一个并计算另一个中的一些统计数据。
这是此类功能的文档:
https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.stats.binned_statistic.html
哪个语法是:
scipy.stats.binned_statistic(x, values, statistic='mean', bins=10, range=None)
首先,可以选择要使用的垃圾箱数量,例如Nbins=100。
接下来,可以定义一个用户函数作为输入
(这里是一个如何做到这一点的例子:
How to make user defined functions for binned_statistic),我的例子是一个函数,它估计该 bin 中数据的第 n 个百分位数(我称之为 myperc)。最后定义一个函数,例如它接受x、y、Nbins 和nth(所需的百分位数)并返回binned_statistics 给出3 个输出:statistic(其中所需统计的值) bin),bin_edges,binnumber(您的数据点在哪个 bin),还有位于 bin 中心的 x 的值(bin_center)
def quantile2d(x,y,Nbins,nth):
from numpy import percentile
from scipy.stats import binned_statistic
def myperc(x,n=nth):
return(percentile(x,n))
t=binned_statistic(x,y,statistic=myperc,bins=Nbins)
v=[]
for i in range(len(t[0])): v.append((t[1][i+1]+t[1][i])/2.)
v=np.array(v)
return(t,v)
所以v 和t.statistic 将分别给出定义所需百分位数的曲线的 x 和 y 值。
Nbins=100
nth=30.
t,v=me.quantile2d(x,y,Nbins,nth)
ii=[]
for i in range(Nbins):
ii=ii+np.argwhere(((t.binnumber==i) & (y<t.statistic[i]))).flatten().tolist()
ii=np.array(ii,dtype=int)
最后,这给出了以下情节:
plt.plot(x,y,'o',color='gray',ms=1,zorder=1)
plt.plot(v,t.statistic,'r-',zorder=3)
plt.plot(x[ii],y[ii],'o',color='blue',ms=1,zorder=2)
其中第 30 个百分位数的线以红色显示,该百分位数以下的数据以蓝色显示。