一、算法目的
???????对于线性模型: Y = β 0 + β 1 X + ε (1) Y=\beta_{0}+\beta_{1}X+\varepsilon \tag{1} Y=β0?+β1?X+ε(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??β1?xi?∣(2)
???????首先对
β
0
,
β
1
\beta_{0},\beta_{1}
β0?,β1?加以约束,使得回归直线
y
=
β
0
+
β
1
x
y=\beta_{0}+\beta_{1}x
y=β0?+β1?x过点
(
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??xk?yi′?=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′??β1?xi′?∣(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??β1?xi?∣(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??βi?xi?∣的最小值。如图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=13?Qi?(β1?)的图形,如图2.2。
可以看出
∑
i
=
1
3
Q
i
(
β
1
)
\sum_{i=1}^{3}Q_{i}(\beta_{1})
∑i=13?Qi?(β1?)是一条折线凸函数。该结论对
∑
i
=
1
n
Q
i
(
β
1
)
\sum_{i=1}^{n}Q_{i}(\beta_{1})
∑i=1n?Qi?(β1?)也是成立的,即
∑
i
=
1
n
Q
i
(
β
1
)
\sum_{i=1}^{n}Q_{i}(\beta_{1})
∑i=1n?Qi?(β1?)是一折线凸函数。
???????所以,对于 Q ( β 1 ) = ∑ i = 1 n Q i ( β 1 ) Q(\beta_{1})=\sum_{i=1}^{n}Q_{i}(\beta_{1}) Q(β1?)=∑i=1n?Qi?(β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.