一元線性模型的中位數回歸
- 一、演算法目的
- 二、演算法推導
- 三、實際案例與python編程計算
- 3.1中位數回歸方程的計算
- 3.2資料可視化
- 3.3回歸引數檢驗
- 四、參考文獻
一、演算法目的
對于線性模型: 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.
轉載請註明出處,本文鏈接:https://www.uj5u.com/houduan/275794.html
標籤:python
