一、分位数回归概念

分位数回归是估计一组回归变量X与被解释变量Y的分位数之间线性关系的建模方法。

以往的回归模型实际上是研究被解释变量的条件期望。而人们也关心解释变量与被解释变量分布的中位数，分位数回归是Koenker和Bassett针对最小二乘的不足提出来的一种新的估计方法，它不仅能够度量回归变量对分布中心的影响，而且能度量回归变量对分布上尾和下尾的影响，在不同的分位数下进行预测，得到的信息更为全面和精确。时间序列分析是对时间序列数据建立模型,分析数据的内部结构和自身规律,从而对未来的发展进行预测。用分位数回归方法来估计时间序列模型时,对随机误差的分布不做要求,能更加全面的刻画分布的特征。

OLS回归估计量的计算是基于最小化残差平方。分位数回归估计量的计算也是基于一种非对称形式的绝对值残差最小化。其中，中位数回归运用的是最小绝对值离差估计(LAD，least absolute deviations estimator)。

分位数回归的优点
(1)能够更加全面的描述被解释变量条件分布的全貌，而不是仅仅分析被解释变量的条件期望（均值），也可以分析解释变量如何影响被解释变量的中位数、分位数等。
(2)不同分位数下的回归系数估计量常常不同，即解释变量对不同水平被解释变量的影响不同。
(3)中位数回归的估计方法与最小二乘法相比，估计结果对离群值则表现的更加稳健，而且，分位数回归对误差项并不要求很强的假设条件，因此对于非正态分布而言，分位数回归系数估计量则更加稳健。

二、普通回归优化为分位数回归的过程：

在一般线性回归中，我们估计的是一些变量y的平均值，条件是自变量x的值。
当我们在数据上拟合一般最小二乘回归模型时，我们对线性模型中的随机误差项做了一个关键假设。我们假设误差项在自变量x的值上有一个常数方差。

• 当这个假设不再成立时会发生什么？
• 另外，除了估计自变量的平均值，我们还能估计自变量的中位数、0.3分位数或0.8分位数吗？

这就是分位数回归发挥作用的地方。
接下来将编写一些代码来更好地理解这一点。创建一些数据并绘制出来。

import numpy as np
import matplotlib.pyplot as plt

# 生成一些具有恒定方差/噪声的数据
x = np.arange(100).reshape(100,1)
intercept_ = 6
slope_ = 0.1
# 非常数误差
error_ = np.random.normal(size = (100,1), loc = 0.0, scale = 1)
# 回归方程
y = intercept_ + slope_ * x + error_

plt.figure(1)
plt.scatter(x, y)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Data with constant variance")

自变量为x，因变量为y，噪声，误差_u是高斯单位方差。

Figure_1.png

如图所示，当我们沿着x轴从左向右移动时，我们看不到y值有很大的变化。
一般的最小二乘回归是建立数据模型的理想候选。下面对数据进行最小二乘回归。

# 对上述数据集进行最小二乘回归
from sklearn.linear_model import LinearRegression
model1 = LinearRegression(fit_intercept = True, normalize = False)
model1.fit(x, y)
y_pred1 = model1.predict(x)
print("Mean squared error: {0:.2f}".format(np.mean((y_pred1 - y) ** 2)))
print('Variance score: {0:.2f}'.format(model1.score(x, y)))

# 绘制回归图
plt.figure(2)
plt.scatter(x, y,  color='black')
plt.plot(x, y_pred1, color='blue', linewidth=3)

plt.xticks(())
plt.yticks(())
plt.xlabel("x")
plt.ylabel("y and predicted y")
plt.title("Linear regression")

Figure_2.png

方差得分为1.0，我们对数据进行了完美的建模。我们的回归线图也证实了这一点。
现在在数据中引入一些可变噪声。我们的噪声根据x值的范围而变化。

# 生成一些具有非常量方差的数据
x_ = np.arange(100).reshape(100,1)
intercept_ = 6
slope_ = 0.1
# 非常数方差
var_ = 0.1 + 0.05 * x_
# 非常数误差
error_ = np.random.normal(size = (100,1), loc = 0.0, scale = var_)
# 回归方程
y_ = intercept_ + slope_ * x + error_

plt.figure(3)
plt.scatter(x_, y_)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Data with non-constant variance")

Figure_3.png

误差计算的比例参数不再是前一种情况下的1。比例是x值的线性函数。
当y的变异性在x值的范围内不相等时，这种现象称为异方差。如图所示，它呈圆锥形。Y变量随着x值的增加而变宽。
接下来尝试将线性回归拟合到此数据集。

# 尝试拟合线性回归
model2 = LinearRegression(fit_intercept = True, normalize = False)
model2.fit(x_, y_)

y_pred2 = model2.predict(x_)

print("Mean squared error: {0:.2f}"
      .format(np.mean((y_pred2 - y_) ** 2)))
print('Variance score: {0:.2f}'.format(model1.score(x_, y_)))

# 绘制回归图
plt.figure(4)
plt.scatter(x_, y_,  color='black')
plt.plot(x_, y_pred2, color='blue', linewidth=3)

plt.xticks(())
plt.yticks(())
plt.xlabel("x")
plt.ylabel("y and predicted y")
plt.title("Linear regression on data with non-constant variance")

Figure_4.png

方差得分为0.43的线性回归总体效果不好。当x值接近0时，线性回归可以很好地估计y，
但是我们接近x值的末尾时，预测的y与实际值相差很远，因此变得完全没有意义。

这就是分位数回归的救星。我使用了python包statsmodels 0.8.0进行分位数回归。
让我们从找到条件中位数0.5分位数的回归系数开始。

# 中位数的分位数回归，0.5分位数
import pandas as pd
data = pd.DataFrame(data = np.hstack([x_, y_]), columns = ["x", "y"])
print(data.head())

import statsmodels.formula.api as smf
mod = smf.quantreg('y ~ x', data)
res = mod.fit(q=.5)
print(res.summary())

首先，我们将数据放入一个pandas数据框架中，这样我们就可以更容易地使用StatsModel接口。
我们的数据帧数据有两列：“x”和“y”。然后，我们继续为中位数0.5分位数建立分位数回归模型。

模型的总结： 屏幕快照 2019-05-07 上午11.00.20.png

截距是6.0849，斜率或x的系数是0.1027。这些是Y的0.5分位数的参数。同样，可以为其他分位数建立模型。

# 为其他分位数建立模型
quantiles = np.arange(0.1,1,0.1)
print(quantiles)
models = []
params = []

for qt in quantiles:
    print(qt)
    res = mod.fit(q = qt )
    models.append(res)
    params.append([qt, res.params['Intercept'], res.params['x']] + res.conf_int().ix['x'].tolist())

params = pd.DataFrame(data = params, columns = ['qt','intercept','x_coef','cf_lower_bound','cf_upper_bound'])
print(params)

在for循环的一侧，为列表中的每个分位数构建模型。在构建这些模型时，还将模型参数存储在一个名为params的列表中。制作了一个同名的数据框架，这样我们就可以查看不同的模型。

数据框架的截图 屏幕快照 2019-05-07 上午11.00.32.png

正如在上面的输出中看到的，0.1th分位数的截距值是5.863，斜率是0.042，还有下限和上限都在结果中输出了，也就是x截距值的间隔。

# 根据原始数据绘制0.1th、0.5和0.9分位数模型
plt.figure(5)
plt.scatter(x_, y_,  color='black')
plt.plot(x_, y_pred2, color='blue',
         linewidth=3, label='Lin Reg')

y_pred3 = models[0].params['Intercept'] + models[0].params['x'] * x_
plt.plot(x_, y_pred3, color='red', linewidth=3, label='Q Reg : 0.1')

y_pred4 = models[4].params['Intercept'] + models[4].params['x'] * x_
plt.plot(x_, y_pred4, color='green', linewidth=3, label='Q Reg : 0.5')

y_pred5 = models[8].params['Intercept'] + models[8].params['x'] * x_
plt.plot(x_, y_pred5, color='cyan',
         linewidth=3, label='Q Reg : 0.9')

plt.xticks(())
plt.yticks(())
plt.xlabel("x")
plt.ylabel("y and predicted y")
plt.title("Quantile regression on data with non-constant variance")
plt.legend()

普通线性回归模型用蓝色线绘制。您可以将该模型与其他分位数模型进行比较。
另一种有趣的可视化方法是不同分位数的坡度值及其上/下限。

# 绘制分位数系数的变化图
plt.figure(6)
params.plot(x = 'qt', y=['x_coef', 'cf_lower_bound', 'cf_upper_bound'],
            title = 'Slope for different quantiles', kind ='line', style = ['b-','r--','g--'])
plt.show()

Figure_7-1.png

可以看到不同分位数的坡度值是如何变化的。与所有分位数的普通最小二乘回归相比，分位数回归允许我们调查数据的不同区域并对其进行适当的建模。