python scipy.stats 分位数_分位数回归及其Python源码导读

天朗气清，惠风和畅。赋闲在家，正宜读书。前人文章，不得其解。代码开源，无人注释。你们不来，我行我上。废话少说，直入主题。o(*￣︶￣*)o

我们要探测自变量

与因变量

的关系，最简单的方法是线性回归，即假设：

我们通过最小二乘方法 (OLS: ordinary least squares)^[1]得到

的无偏估计

，

。为了解决

，

的可靠性问题，我们同时对残差

做了假设，即：

为均值为0，方差恒定的独立随机变量。

即为给定自变量

下，因变量

的条件均值。

假如残差

不满足我们的假设，或者更重要地，我们不仅仅想要知道

的在给定

下的条件均值，而且想知道条件中位数（更一般地，条件分位数），那么OLS下的线性回归就不能满足我们的需求。分位数回归(Quantile Regression)

^[2]解决了这些问题，下面我先给出一个分位数回归的实际应用例子，再简述其原理，最后再分析其在Python实现的源代码。

1. 一个例子：收入与食品消费

这个例子出自statasmodels:Quantile Regression.^[3] 我们想探索家庭收入与食品消费的关系，数据出自working class Belgian households in 1857 (the Engel data).我们用Python包statsmodels实现分位数回归。

1.1 预处理

%

	income	        foodexp
0	420.157651	255.839425
1	541.411707	310.958667
2	901.157457	485.680014
3	639.080229	402.997356
4	750.875606	495.560775

1.2 中位数回归 （分位数回归的特例，q=0.5）

mod 

QuantReg 

由结果可以知道

，如何得到回归系数的估计？结果中的std err, t, Pseudo R-squared等是什么？我会在稍后解释。

1.3 数据可视化

我们先拟合10个分位数回归，分位数q分别在0.05到0.95之间。

quantiles 

      q           a         b        lb        ub
0  0.05  124.880096  0.343361  0.268632  0.418090
1  0.15  111.693660  0.423708  0.382780  0.464636
2  0.25   95.483539  0.474103  0.439900  0.508306
3  0.35  105.841294  0.488901  0.457759  0.520043
4  0.45   81.083647  0.552428  0.525021  0.579835
5  0.55   89.661370  0.565601  0.540955  0.590247
6  0.65   74.033435  0.604576  0.582169  0.626982
7  0.75   62.396584  0.644014  0.622411  0.665617
8  0.85   52.272216  0.677603  0.657383  0.697823
9  0.95   64.103964  0.709069  0.687831  0.730306
{'a': 147.47538852370562, 'b': 0.48517842367692354, 'lb': 0.4568738130184233,

这里拟合了10个回归，其中q是对应的分位数，a是斜率，b是回归系数。lb和ub分别是b的95%置信区间的下界与上界。

现在来画出这10条回归线：

x 

上图中虚线是分位数回归线，红线是线性最小二乘（OLS）的回归线。通过观察，我们可以发现3个现象:

1. 随着收入提高，食品消费也在提高。
2. 随着收入提高，家庭间食品消费的差别拉大。穷人别无选择，富人能选择生活方式，有喜欢吃贵的，也有喜欢吃便宜的。然而我们无法通过OLS发现这个现象，因为它只给了我们一个均值。
3. 对于穷人来说，OLS预测值过高。这是因为少数的富人拉高了整体的均值，可见OLS对异常点敏感，不是一个稳健的模型。

2.分位数回归的原理

这部分是数理统计的内容，只关心如何实现的朋友可以略过。我们要解决以下这几个问题：

1. 什么是分位数？
2. 如何求分位数？
3. 什么是分位数回归？
4. 分位数回归的回归系数如何求得？
5. 回归系数的检验如何进行？
6. 如何评估回归拟合优度？

2.1 分位数的定义^[4]

是随机变量，

的累积密度函数是

.

的

分位数为：

,

假设有100个人，95%的人身高少于1.9m, 1.9m就是身高的95%分位数。

2.2 分位数的求法^[4]

通过选择不同的

值，使

最小，对应的

值即为

的

分位数的估计

.

2.3 分位数回归^[4]

对于OLS, 我们有：

为

最小化所对应的

，类比地，对于

分位数回归，我们有：

为最小化：

即最小化

所对应的

2.4 系数估计^[5][6]

由于

不能直接对

求导，我们只能用迭代的方法来逼近

最小时对应的

值。statsmodels采用了Iteratively reweighted least squares （IRLS）的方法。

假设我们要求

最小化形如下的

范数：

则第t+1步迭代的

值为：

是对角矩阵且初始值为

第t次迭代

以中位数回归为例子(q=0.5),我们求：

即

即最小化形如上的

范数，

为避免分母为0，我们取

,

是一个很小的数，例如0.0001.

2.5 回归系数的检验^[7]

我们通过2.4，多次迭代得出

的估计值，为了得到假设检验的t统计量，我们只需得到

的方差的估计。

分位数回归

的协方差矩阵的渐近估计为：

其中

是对角矩阵，

当

,

, 当

，

的估计为

其中

为核函数（Kernel），可取Epa,Gaussian等.

为根据Stata 12所选的窗宽(bandwidth)

^[5]

回归结果中的std err即由

获得，t统计量等于

。

2.6 拟合优度

对于OLS,我们用

来衡量拟合优度。对于

分位数回归，我们类比得到：

,其中

为所有

观察值的

分位数。

即为回归结果中的Pseudo R-squared。

3.Python源码分析

实现分位数回归的完整源码在^[5] ，里面主要含有两个类QuantReg 和 QuantRegResults. 其中QuantReg是核心，包含了回归系数估计，协方差计算等过程。QuantRegResults计算拟合优度并组织回归结果。

3.1 QuantReg类

#QuantReg是包中RegressionModel的一个子类
  

3.2 QuantRegResults类

这里我只给出计算拟合优度的代码。

class 

4.总结

上文我先给出了一个分位数回归的应用例子，进而叙述了分位数回归的原理，最后再分析了Python实现的源码。

分位数回归对比起OLS回归，虽然较为复杂，但它有三个主要优势：

1. 能反映因变量分位数与自变量的关系，而不仅仅反映因变量均值与自变量的关系。
2. 分位数回归对残差不作任何假设。
3. 分位数回归受异常点的影响较小。^[8]

（欢迎转载，但请注明出处）

参考

1. ^https://en.wikipedia.org/wiki/Ordinary_least_squares
2. ^QUANTILE REGRESSION http://www.econ.uiuc.edu/~roger/research/rq/rq.pdf
3. ^https://www.statsmodels.org/dev/examples/notebooks/generated/quantile_regression.html
4. ^^abchttps://en.wikipedia.org/wiki/Quantile_regression
5. ^^abchttps://www.statsmodels.org/devel/_modules/statsmodels/regression/quantile_regression.html
6. ^https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares
7. ^Green,W. H. (2008). Econometric Analysis. Sixth Edition. International Student Edition.
8. ^https://www.ibm.com/support/knowledgecenter/en/SSLVMB_sub/statistics_mainhelp_ddita/spss/regression/idh_quantile.html