一、算法目的
对于线性模型: Y = β 0 + β 1 X + ε (1) Y=\beta_{0}+\beta_{1}X+\varepsilon \tag{1} Y=β0+β1X+ε(1) 通过最小一乘法进行中位数回归。给出参数 β 0 \beta_{0} β0和 β 1 \beta_{1} β1的估计 β 0 ^ \widehat{\beta_{0}} β0 和 β 1 ^ \widehat{\beta_{1}} β1 .
二、算法推导
我们的目标函数为:
min
Q
(
β
0
,
β
1
)
=
∑
i
=
1
n
∣
y
i
−
β
0
−
β
1
x
i
∣
(2)
\min Q(\beta_{0},\beta_{1})=\sum_{i=1}^{n}|y_{i}-\beta_{0}-\beta_{1}x_{i}| \tag{2}
minQ(β0,β1)=i=1∑n∣yi−β0−β1xi∣(2)
首先对
β
0
,
β
1
\beta_{0},\beta_{1}
β0,β1加以约束,使得回归直线
y
=
β
0
+
β
1
x
y=\beta_{0}+\beta_{1}x
y=β0+β1x过点
(
x
k
,
y
k
)
(x_{k},y_{k})
(xk,yk),其中
y
k
y_{k}
yk为数据
y
i
,
i
=
1
,
2
,
⋯
,
n
y_{i},i=1,2,\cdots,n
yi,i=1,2,⋯,n的中位数。
作变换如下:
{
x
i
′
=
x
i
−
x
k
y
i
′
=
y
i
−
y
k
,
i
=
1
,
2
,
⋯
,
n
(3)
\begin{cases} x_{i}'=x_{i}-x_{k}\\ y_{i}'=y_{i}-y_{k} \end{cases},i=1,2,\cdots,n \tag{3}
{xi′=xi−xkyi′=yi−yk,i=1,2,⋯,n(3)
这样,(2)就转化成了(4).
min
Q
(
β
1
)
=
∑
i
=
1
n
∣
y
i
′
−
β
1
x
i
′
∣
(4)
\min Q(\beta_{1})=\sum_{i=1}^{n}|y_{i}'-\beta_{1}x_{i}'| \tag{4}
minQ(β1)=i=1∑n∣yi′−β1xi′∣(4)
为了后面的书写运算方便,我们仍用
(
x
i
,
y
i
)
(x_{i},y_{i})
(xi,yi)来表示经过变换(3)得到的数据
(
x
i
′
,
y
i
′
)
(x_{i}',y_{i}')
(xi′,yi′),这样(4)式就可以写成下式:
min
Q
(
β
1
)
=
∑
i
=
1
n
∣
y
i
−
β
1
x
i
∣
(5)
\min Q(\beta_{1})=\sum_{i=1}^{n}|y_{i}-\beta_{1}x_{i}| \tag{5}
minQ(β1)=i=1∑n∣yi−β1xi∣(5)
在求解(5)之前我们先来考虑对于任意的
i
i
i,
Q
i
(
β
1
)
=
∣
y
i
−
β
i
x
i
∣
Q_{i}(\beta_{1})=|y_{i}-\beta_{i}x_{i}|
Qi(β1)=∣yi−βixi∣的最小值。如图2.1,在
β
=
y
i
/
x
i
\beta=y_{i}/x_{i}
β=yi/xi时,
Q
i
(
β
)
Q_{i}(\beta)
Qi(β)取最小值,图中两条直线的斜率互为相反数分别为
−
∣
x
i
∣
-|x_{i}|
−∣xi∣和
∣
x
i
∣
|x_{i}|
∣xi∣。
现在考虑如下数据表:
序号 | x i x_{i} xi | y i y_{i} yi |
---|---|---|
1 | 1 | 3 |
2 | 1 | 1 |
3 | 2 | 4 |
分别画出
Q
1
(
β
1
)
=
∣
3
−
β
1
∣
,
Q
2
(
β
1
)
=
∣
1
−
β
1
∣
,
Q
3
(
β
1
)
=
∣
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}|
Q1(β1)=∣3−β1∣,Q2(β1)=∣1−β1∣,Q3(β1)=∣4−2β1∣和
∑
i
=
1
3
Q
i
(
β
1
)
\sum_{i=1}^{3}Q_{i}(\beta_{1})
∑i=13Qi(β1)的图形,如图2.2。
可以看出
∑
i
=
1
3
Q
i
(
β
1
)
\sum_{i=1}^{3}Q_{i}(\beta_{1})
∑i=13Qi(β1)是一条折线凸函数。该结论对
∑
i
=
1
n
Q
i
(
β
1
)
\sum_{i=1}^{n}Q_{i}(\beta_{1})
∑i=1nQi(β1)也是成立的,即
∑
i
=
1
n
Q
i
(
β
1
)
\sum_{i=1}^{n}Q_{i}(\beta_{1})
∑i=1nQi(β1)是一折线凸函数。
所以,对于 Q ( β 1 ) = ∑ i = 1 n Q i ( β 1 ) Q(\beta_{1})=\sum_{i=1}^{n}Q_{i}(\beta_{1}) Q(β1)=∑i=1nQi(β1), Q ( β 1 ) Q(\beta_{1}) Q(β1)有 n n n个折点。在 β 1 \beta_{1} β1轴上,每个折点的横坐标为 y i / x i , i = 1 , 2 , ⋯ , n y_{i}/x_{i},i=1,2,\cdots,n yi/xi,i=1,2,⋯,n。下面我们按照 y i / x i y_{i}/x_{i} yi/xi的大小来排序,我们令最小的 y i / x i = β ( 1 ) y_{i}/x_{i}=\beta_{(1)} yi/xi=β(1),同时,将 β ( 1 ) \beta_{(1)} β(1)所对应的 Q i ( β 1 ) Q_{i}(\beta_{1}) Qi(β1)重新赋值给 Q 1 ( β 1 ) Q_{1}(\beta_{1}) Q1(β1)以此类推: β ( 1 ) ≤ β ( 2 ) ≤ ⋯ ≤ β ( n ) \beta_{(1)}\leq \beta_{(2)}\leq \cdots\leq \beta_{(n)} β(1)≤β(2)≤⋯≤β(n),且 n n n个 Q i ( β 1 ) Q_{i}(\beta_{1}) Qi(β1)在平面坐标系中从左到右依次为 Q 1 ( β 1 ) , Q 2 ( β 1 ) , ⋯ , Q n ( β 1 ) Q_{1}(\beta_{1}),Q_{2}(\beta_{1}),\cdots,Q_{n}(\beta_{1}) Q1(β1),Q2(β1),⋯,Qn(β1),且它们原来所对应的 x i x_{i} xi分别重新记为 x 1 , x 2 , ⋯ , x n x_{1},x_{2},\cdots,x_{n} x1,x2,⋯,xn,则当 β ≤ β ( 1 ) \beta\leq \beta_{(1)} β≤β(1)时,对于所有的 i i i都有 Q i ( β ) Q_{i}(\beta) Qi(β)的斜率为 − ∣ x i ∣ -|x_{i}| −∣xi∣,故此时, Q ( β 1 ) = − ∑ i = 1 n ∣ x i ∣ Q(\beta_{1})=-\sum_{i=1}^{n}|x_{i}| Q(β1)=−∑i=1n∣xi∣,我们记 M = ∑ i = 1 n ∣ x i ∣ M=\sum_{i=1}^{n}|x_{i}| M=∑i=1n∣xi∣。
当 β ( 1 ) ≤ β ≤ β ( 2 ) \beta_{(1)}\leq \beta\leq\beta_{(2)} β(1)≤β≤β(2)时,由于 Q 1 ( β 1 ) Q_{1}(\beta_{1}) Q1(β1)的斜率为 ∣ x i ∣ |x_{i}| ∣xi∣,故 Q ( β 1 ) = − ∑ i = 1 n ∣ x i ∣ + 2 ∣ x 1 ∣ Q(\beta_{1})=-\sum_{i=1}^{n}|x_{i}|+2|x_{1}| Q(β1)=−∑i=1n∣xi∣+2∣x1∣。
由此可见,折线 Q ( β 1 ) Q(\beta_{1}) Q(β1)的斜率只在 β 1 = β ( i ) , i = 1 , 2 , ⋯ , n \beta_{1}=\beta_{(i)},i=1,2,\cdots,n β1=β(i),i=1,2,⋯,n时改变。所以 β 1 = β ( i ) , i = 1 , 2 , ⋯ , n \beta_{1}=\beta_{(i)},i=1,2,\cdots,n β1=β(i),i=1,2,⋯,n是折线 Q ( β 1 ) Q(\beta_{1}) Q(β1)的折点的横坐标,并且只有这些折点。从左到右经过每个折点,斜率就会增加 2 ∣ x i ∣ 2|x_{i}| 2∣xi∣。根据这些,我们就可以找到求出 Q ( β 1 ) Q(\beta_{1}) Q(β1)最小值的方法。
设最小值点在 β ( r ) \beta_{(r)} β(r)处达到,那么对于 r r r有
{ − ∑ i = 1 n ∣ x i ∣ + 2 ∑ j = 1 r − 1 ∣ x 1 ∣ < 0 , − ∑ i = 1 n ∣ x i ∣ + 2 ∑ j = 1 r ∣ x 1 ∣ ≥ 0 (6) \begin{cases} -\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}|\ge0 \end{cases} \tag{6} {−∑i=1n∣xi∣+2∑j=1r−1∣x1∣<0,−∑i=1n∣xi∣+2∑j=1r∣x1∣≥0(6)或者等价表示为
{ ∑ j = 1 r − 1 ∣ x 1 ∣ < M 2 , ∑ j = 1 r ∣ x 1 ∣ ≥ M 2 (7) \begin{cases} \sum_{j=1}^{r-1}|x_{1}|<\frac{M}{2},\\ \sum_{j=1}^{r}|x_{1}|\ge\frac{M}{2} \end{cases} \tag{7} {∑j=1r−1∣x1∣<2M,∑j=1r∣x1∣≥2M(7)
下面我们分两种情况进行讨论:
(1)若 ∑ j = 1 r ∣ x 1 ∣ > M 2 \sum_{j=1}^{r}|x_{1}|>\frac{M}{2} ∑j=1r∣x1∣>2M,此时 β 1 ^ = β ( r ) \widehat{\beta_{1}}=\beta_{(r)} β1 =β(r),由 y k = β 0 ^ + β 1 ^ x k y_{k}=\widehat{\beta_{0}}+\widehat{\beta_{1}}x_{k} yk=β0 +β1 xk得到, β 0 ^ = y k − β 1 ^ x k = y k − β ( r ) x k \widehat{\beta_{0}}=y_{k}-\widehat{\beta_{1}}x_{k}=y_{k}-\beta_{(r)}x_{k} β0 =yk−β1 xk=yk−β(r)xk。此时中位数回归直线方程为: y = y k − β ( r ) x k + β ( r ) x (8) y=y_{k}-\beta_{(r)}x_{k}+\beta_{(r)}x \tag{8} y=yk−β(r)xk+β(r)x(8)
(2)若 ∑ j = 1 r ∣ x 1 ∣ = M 2 \sum_{j=1}^{r}|x_{1}|=\frac{M}{2} ∑j=1r∣x1∣=2M,此时 [ β ( r ) , β ( r + 1 ) ] [\beta_{(r)},\beta_{(r+1)}] [β(r),β(r+1)]上所有的点都是 β 1 \beta_{1} β1的最优估计,我们不妨就取 β 1 ^ = β ( r ) \widehat{\beta_{1}}=\beta_{(r)} β1 =β(r)。最后中位数回归直线方程于式(8)相同。
三、实际案例与python编程计算
3.1中位数回归方程的计算
我们来考虑如下这样一组数据:
y y y | x x x |
---|---|
220 | 4 |
146 | 3 |
438 | 7 |
49 | 1 |
95 | 2 |
303 | 6 |
261 | 5 |
下面给出计算中位数回归方程完整 p y t h o n python 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数据可视化
下面画出回归方程与数据点的图:
可以看出中位数回归直线一定过点 ( x k , y k ) (x_{k},y_{k}) (xk,yk),而且还会经过数据点中的另一个点。
数据可视化部分的代码如下:
#数据可视化
#画出数据散点图
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回归参数检验
下面我们对回归的情况进行一些必要的检验:
参数检验代码如下:
#回归参数检验
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.