我试图使用回归线的第5和95个百分位数来识别数据集中的异常值，所以我在Python中使用了分位数回归和statsmodel、matplotlib和pandas。基于blokeley的this answer，我可以创建数据的散点图，并基于分位数回归显示最佳拟合线以及第5百分位数和第95百分位数的线。但是，我如何识别那些位于这些行之上和之下的点，然后将它们保存到pandas数据框中呢？
我的数据如下（总共有95个值）：
Month Year LST NDVI
0 June 1984 310.550975 0.344335
1 June 1985 310.495331 0.320504
2 June 1986 306.820900 0.369494
3 June 1987 308.945602 0.369946
4 June 1988 308.694022 0.31863
二
到目前为止我的剧本是这样的：
import pandas as pd
excel = my_excel
df = pd.read_excel(excel)
df.head()
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
model = smf.quantreg('NDVI ~ LST',df)
quantiles = [0.05,0.95]
fits = [model.fit(q=q) for q in quantiles]
figure,axes = plt.subplots()
x = df['LST']
y = df['NDVI']
axes.scatter(x,df['NDVI'],c='green',alpha=0.3,label='data point')
fit = np.polyfit(x, y, deg=1)
axes.plot(x, fit[0] * x + fit[1], color='grey',label='best fit')
_x = np.linspace(x.min(),x.max())
for index, quantile in enumerate(quantiles):
_y = fits[index].params['LST'] * _x + fits[index].params['Intercept']
axes.plot(_x, _y, label=quantile)
title = 'LST/NDVI Jun-Aug'
plt.title(title)
axes.legend()
axes.set_xticks(np.arange(298,320,4))
axes.set_yticks(np.arange(0.25,0.5,.05))
axes.set_xlabel('LST')
axes.set_ylabel('NDVI');
我得到的图表是：
因此，我可以明确地看到第95行以上和第5行以下的数据点，我会将它们归类为异常值，但我想识别原始数据框中的数据点，并可能将它们绘制在购物车上，或以某种方式突出显示为“异常值”。
我正在寻找一种方法，但却一无所获，需要一些帮助。
最佳答案
你需要弄清楚某个点是高于95%分位数线还是低于5%分位数线。这可以使用交叉产品来实现，请参见this answer以获得简单的实现。
在您的示例中，您需要将分位数线上方和下方的点组合在一起，可能是在遮罩中。
下面是一个例子：
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
df = pd.DataFrame(np.random.normal(0, 1, (100, 2)))
df.columns = ['LST', 'NDVI']
model = smf.quantreg('NDVI ~ LST', df)
quantiles = [0.05, 0.95]
fits = [model.fit(q=q) for q in quantiles]
figure, axes = plt.subplots()
x = df['LST']
y = df['NDVI']
fit = np.polyfit(x, y, deg=1)
_x = np.linspace(x.min(), x.max(), num=len(y))
# fit lines
_y_005 = fits[0].params['LST'] * _x + fits[0].params['Intercept']
_y_095 = fits[1].params['LST'] * _x + fits[1].params['Intercept']
# start and end coordinates of fit lines
p = np.column_stack((x, y))
a = np.array([_x[0], _y_005[0]]) #first point of 0.05 quantile fit line
b = np.array([_x[-1], _y_005[-1]]) #last point of 0.05 quantile fit line
a_ = np.array([_x[0], _y_095[0]])
b_ = np.array([_x[-1], _y_095[-1]])
#mask based on if coordinates are above 0.95 or below 0.05 quantile fitlines using cross product
mask = lambda p, a, b, a_, b_: (np.cross(p-a, b-a) > 0) | (np.cross(p-a_, b_-a_) < 0)
mask = mask(p, a, b, a_, b_)
axes.scatter(x[mask], df['NDVI'][mask], facecolor='r', edgecolor='none', alpha=0.3, label='data point')
axes.scatter(x[~mask], df['NDVI'][~mask], facecolor='g', edgecolor='none', alpha=0.3, label='data point')
axes.plot(x, fit[0] * x + fit[1], label='best fit', c='lightgrey')
axes.plot(_x, _y_095, label=quantiles[1], c='orange')
axes.plot(_x, _y_005, label=quantiles[0], c='lightblue')
axes.legend()
axes.set_xlabel('LST')
axes.set_ylabel('NDVI')
plt.show()
关于python - 用分位数回归和Python识别异常值，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/51791790/