一、算法目的

???????对于线性模型：$$Y={\beta }_0+{\beta }_{1}X+\epsilon \text{\hspace{1em}(1)}$$???????通过最小一乘法进行中位数回归。给出参数${\beta }_0$和${\beta }_1$的估计$\widehat{{\beta }_0}$和$\widehat{{\beta }_1}$.

二、算法推导

???????我们的目标函数为：
$$\min Q({\beta }_0,{\beta }_1)=\sum _{i=1}^n|y_i?{\beta }_0?{\beta }_1{x}_i|\text{\hspace{1em}(2)}$$
???????首先对${\beta }_0,{\beta }_1$加以约束，使得回归直线$y={\beta }_0+{\beta }_{1}x$过点$(x_k,y_k)$，其中$y_k$为数据$y_i,i=1,2,?,n$的中位数。
作变换如下：
$$\{\begin{array}{l}{x}_i^{\prime }=x_i?x_k\\ y_i^{\prime }=y_i?y_k\end{array},i=1,2,?,n\text{\hspace{1em}(3)}$$
???????这样，(2)就转化成了(4).
$$\min Q({\beta }_1)=\sum _{i=1}^n|y_i^{\prime }?{\beta }_1{x}_i^{\prime }|\text{\hspace{1em}(4)}$$
???????为了后面的书写运算方便，我们仍用$(x_i,y_i)$来表示经过变换(3)得到的数据$(x_i^{\prime },y_i^{\prime })$,这样(4)式就可以写成下式：
$$\min Q({\beta }_1)=\sum _{i=1}^n|y_i?{\beta }_1{x}_i|\text{\hspace{1em}(5)}$$
???????在求解(5)之前我们先来考虑对于任意的$i$，$Q_i({\beta }_1)=|y_i?{\beta }_i{x}_i|$的最小值。如图2.1，在$\beta =y_i/x_i$时，$Q_i(\beta )$取最小值，图中两条直线的斜率互为相反数分别为$?|x_i|$和$|x_i|$。

图2.1

???????现在考虑如下数据表：

序号	$x_i$	$y_i$
1	1	3
2	1	1
3	2	4

表2.1

分别画出$Q_1({\beta }_1)=|3?{\beta }_1|,Q_2({\beta }_1)=|1?{\beta }_1|,Q_3({\beta }_1)=|4?2{\beta }_1|$和${\sum }_{i=1}^3{Q}_i({\beta }_1)$的图形，如图2.2。
可以看出${\sum }_{i=1}^3{Q}_i({\beta }_1)$是一条折线凸函数。该结论对${\sum }_{i=1}^n{Q}_i({\beta }_1)$也是成立的，即${\sum }_{i=1}^n{Q}_i({\beta }_1)$是一折线凸函数。

图2.2
???????所以，对于 $Q({\beta }_1)={\sum }_{i=1}^n{Q}_i({\beta }_1)$， $Q({\beta }_1)$有 $n$个折点。在 ${\beta }_1$轴上，每个折点的横坐标为 $y_i/x_i,i=1,2,?,n$。下面我们按照 $y_i/x_i$的大小来排序,我们令最小的 $y_i/x_i={\beta }_{(1)}$，同时，将 ${\beta }_{(1)}$所对应的 $Q_i({\beta }_1)$重新赋值给 $Q_1({\beta }_1)$以此类推： ${\beta }_{(1)}\le {\beta }_{(2)}\le ?\le {\beta }_{(n)}$，且 $n$个 $Q_i({\beta }_1)$在平面坐标系中从左到右依次为 $Q_1({\beta }_1),Q_2({\beta }_1),?,Q_n({\beta }_1)$，且它们原来所对应的 $x_i$分别重新记为 $x_1,x_2,?,x_n$，则当 $\beta \le {\beta }_{(1)}$时，对于所有的 $i$都有 $Q_i(\beta )$的斜率为 $?|x_i|$,故此时， $Q({\beta }_1)=?{\sum }_{i=1}^n|x_i|$，我们记 $M={\sum }_{i=1}^n|x_i|$。
???????当 ${\beta }_{(1)}\le \beta \le {\beta }_{(2)}$时，由于 $Q_1({\beta }_1)$的斜率为 $|x_i|$，故 $Q({\beta }_1)=?{\sum }_{i=1}^n|x_i|+2|x_1|$。
???????由此可见，折线 $Q({\beta }_1)$的斜率只在 ${\beta }_1={\beta }_{(i)},i=1,2,?,n$时改变。所以 ${\beta }_1={\beta }_{(i)},i=1,2,?,n$是折线 $Q({\beta }_1)$的折点的横坐标，并且只有这些折点。从左到右经过每个折点，斜率就会增加 $2|x_i|$。根据这些，我们就可以找到求出 $Q({\beta }_1)$最小值的方法。
???????设最小值点在 ${\beta }_{(r)}$处达到，那么对于 $r$有
$$\{\begin{array}{l}?{\sum }_{i=1}^n|x_i|+2{\sum }_{j=1}^{r?1}|x_1|<0,\\ ?{\sum }_{i=1}^n|x_i|+2{\sum }_{j=1}^r|x_1|\ge 0\end{array}\text{\hspace{1em}(6)}$$或者等价表示为
$$\{\begin{array}{l}{\sum }_{j=1}^{r?1}|x_1|<\frac{M}{2},\\ {\sum }_{j=1}^r|x_1|\ge \frac{M}{2}\end{array}\text{\hspace{1em}(7)}$$
???????下面我们分两种情况进行讨论：
（1）若 ${\sum }_{j=1}^r|x_1|>\frac{M}{2}$,此时 $\widehat{{\beta }_1}={\beta }_{(r)}$，由 $y_k=\widehat{{\beta }_0}+\widehat{{\beta }_1}{x}_k$得到， $\widehat{{\beta }_0}=y_k?\widehat{{\beta }_1}{x}_k=y_k?{\beta }_{(r)}{x}_k$。此时中位数回归直线方程为： $$y=y_k?{\beta }_{(r)}{x}_k+{\beta }_{(r)}x\text{\hspace{1em}(8)}$$
（2）若 ${\sum }_{j=1}^r|x_1|=\frac{M}{2}$,此时 $[{\beta }_{(r)},{\beta }_{(r+1)}]$上所有的点都是 ${\beta }_1$的最优估计，我们不妨就取 $\widehat{{\beta }_1}={\beta }_{(r)}$。最后中位数回归直线方程于式(8)相同。

三、实际案例与python编程计算

3.1中位数回归方程的计算

???????我们来考虑如下这样一组数据：

$y$	$x$
220	4
146	3
438	7
49	1
95	2
303	6
261	5

表3.1.1

???????下面给出计算中位数回归方程完整$python$源代码：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math
#一元线性模型的中位数回归
#导入数据
dataset=pd.read_excel('LAD for dimension of one.xlsx')

#寻找y的中位数,因为取中位数函数np.median()有时会把中间两项求平均值，这与我们的目的不和，故自己写一个程序
dataset=dataset.sort_values(by='y',ascending=True)
n=len(dataset)#数据量大小n
k=math.ceil(n/2)#中位数下标k
yk=dataset.iloc[k-1,0]
xk=dataset.iloc[k-1,1]

#作以中位数为中心的变换
reversef_x=lambda x:x-xk
reversef_y=lambda y:y-yk
dataset['x']=reversef_x(dataset['x'])
dataset['y']=reversef_y(dataset['y'])

#删除xi=0的数据
dataset=dataset.drop(dataset[dataset['x'].isin([0])].index)

#计算βi=yi/xi,并以βi的大小排序
dataset['beta']=dataset['y']/dataset['x']
dataset=dataset.sort_values(by='beta',ascending=True)

#计算M
fabs=lambda x:abs(x)
dataset['|x|']=fabs(dataset['x'])
M=dataset['|x|'].sum()

#寻找最优的βr
a=0
for i in range(len(dataset)):
   a+=dataset.iloc[i,3]
   if a<M/2:
       continue
   else:
       beta_r=dataset.iloc[i,2]
       break
#输出结果
print("中位数回归曲线方程为：y={}+{}x".format(yk-beta_r*xk,beta_r))

???????下面给出程序运行结果：

3.2数据可视化

???????下面画出回归方程与数据点的图：

图3.2.1
???????可以看出中位数回归直线一定过点 $(x_k,y_k)$，而且还会经过数据点中的另一个点。
???????数据可视化部分的代码如下：

#数据可视化
#画出数据散点图
dataset=pd.read_excel('LAD for dimension of one.xlsx')
plt.scatter(dataset['x'],dataset['y'],c='red',label='原始数据点')
#画出回归直线图
x=np.arange(0,9,1)
y=[]
for i in x:
    yi=yk-beta_r*xk+beta_r*i
    y.append(yi)
plt.plot(x,y,label='回归直线')
#设置图例
plt.legend(loc='upper left',frameon=True)
plt.title('中位数回归直线',size=15)
plt.xlabel('x',size=15)
plt.ylabel('y',size=15)
plt.show()

3.3回归参数检验

???????下面我们对回归的情况进行一些必要的检验：

图3.3.1
???????参数检验代码如下：

#回归参数检验
def get_lr_stats(x, y):
    n = len(x)
    y_prd = yk-beta_r*xk+beta_r*x
    Regression = sum((y_prd - np.mean(y)) ** 2)  # 回归平方和
    Residual = sum((y - y_prd) ** 2)  # 残差平方和
    total = sum((y - np.mean(y)) ** 2)  # 总体平方和
    R_square = 1 - Residual / total  # 相关性系数R^2
    message1 = ('相关系数(R^2)： ' + str(R_square) + '；' + '\n' + '总体平方和(TSS)： ' + str(total) + '；' + '\n')
    message2 = ('回归平方和(RSS)： ' + str(Regression) + '；' + '\n残差平方和(ESS)： ' + str(Residual) + '；' + '\n')
    return print('\n' + message1 + message2)
get_lr_stats(dataset['x'], dataset['y'])

四、参考文献

[1]李仲来.最小一乘法介绍[J].数学通报,1992(02):42-45.