資料科學必備基礎之線性迴歸

挖地兔訂閱號發表於2022-12-08

翻譯整理 | One    本期編輯 | 一隻小綠怪獸 

譯者簡介:西南財經大學應用數學本科,英國曼徹斯特大學金融數學碩士,金融分析師,專注於利用資料建立金融模型,發掘潛在投資價值。

作者:Mirko Stojiljkovic

線性迴歸是統計學和機器學習的基礎之一。無論是做統計,機器學習還是科學計算,都可能用上線性迴歸。因此,建議先掌握線性迴歸,然後再學習更復雜的方法。

透過閱讀本文,你將瞭解到:

· 什麼是線性迴歸

· 線性迴歸的用途

· 線性迴歸的原理

· 在Python中實現線性迴歸


01 

什麼是迴歸

迴歸主要用於研究變數之間的關係。例如,可以透過觀察一些公司的職員,並試著瞭解他們的特徵如何影響他們的薪水,比如經驗,教育水平,職位,工作地點等等。

這就個迴歸問題,與每個職員相關的資料都是代表一次觀察。其中職員的經驗,教育水平,職位,工作地點都是獨立的特徵,這些特徵決定了職員的薪水水平。

類似的,也可以試著建立房屋價格與地域,臥室數量,到市中心距離的數學依賴關係。

在迴歸分析中,經常會發現一些有趣的現象和一些觀察,每次觀察都會有兩個或者更多的特徵。假設其中一個特徵依賴其他幾個特徵,可以試著去建立迴歸發現特徵之間的關係。

換種說法,就是需要找到一個能夠很好地將某些特徵對映到其他特徵的函式。依賴特徵稱為因變數,獨立特徵稱為自變數

迴歸模型中通常有一個連續且無界的因變數。但是,自變數可以是連續的,離散的,甚至是像性別,國籍,品牌等分類資料。

通常用y表示因變數,x表示自變數。如果有兩個或者多個自變數,可以用向量x=(x1, …, xr)表示。


02 

迴歸的用途

通常,迴歸被用來回答某些現象是否會對其他現象產生影響,或者說,幾個變數之間是否存在關聯。例如,可以用迴歸的方法確定工作經驗或者性別是否在某種程度上影響到工資水平。

用新的自變數資料集去預測因變數時,迴歸也是非常有用的。舉個例子,已知室外氣溫,一天的某個時間,以及家裡居住的人數,可以試著預測下一個小時的家庭用電量。

迴歸可用於許多領域:經濟,IT,社會科學等等。隨著大量資料的可獲取性提升以及人們對於資料的實用價值認知的提高,迴歸變得越來越重要了。


03

線性迴歸的原理

線性迴歸是被廣泛使用的迴歸方法中比較重要的一種。它是最基礎也是最簡單的迴歸方法。它的主要優點是易懂。

原理

在對一組自變數x=(x1,…,xr)和因變數y進行迴歸時,我們假設y與x的關係:y=β01x1+···+βrxr。這個方程叫做迴歸方程β0, β1,... ,βr迴歸係數, ε是隨機誤差

線性迴歸計算迴歸係數的估計值,或者說計算以b1,b2,…,br表示並用於預測的權重。其中迴歸函式被定義為f(x)=b0+b1x1+b2x2+···+brxr。迴歸函式可以很好地捕獲自變數與因變數的依賴關係。

對於每一次觀測i=1,…,n迴歸函式f(xi)應該要足夠接近實際觀測的因變數yi值,其中yi與f(xi)的差值稱為殘差。迴歸是為了找到最優權重並用於預測,而殘差最小的時候說明迴歸函式f(xi)與實際觀測的yi是最接近的,此時對應的權重便是最優權重b=(b1,b2,…,br)

通常使用最小化殘差平方和的方法來計算最優權重,即最小化SSR=∑i(yi-f(xi))2。這種方法叫做最小二乘法OLS

評價迴歸模型

實際觀測值yi的變化一部分由於自變數xi的變化導致的。但是實際觀測值yi仍然會有無法用迴歸函式f(xi)解釋的部分。

決定係數R2,用來說明實際觀測值yi中有多大程度可以由迴歸函式f(x)解釋。R2越大,兩者擬合程度越高,意味著觀測值yi-xi資料組可以很好的契合迴歸函式f(x)

R2=1時,即SSR=0,稱作完全擬合,這意味著所有yi-xi可以完全匹配迴歸函式f(x)

單變數線性迴歸

單變數線性迴歸是線性迴歸中最簡單的一種情形,即迴歸模型中只有一個自變數x=x1

資料科學必備基礎之線性迴歸

這些資料集一般來自實際業務中的資料採集。上圖中,最左邊綠點x=5,y=5就是觀測值。

迴歸函式以圖中黑線表示,其表示式為f(x)=b0+b1x。我們需要用最小化殘差平方和SSR的方法來計算迴歸函式中的權重b0與b1。其中b0截距,即迴歸函式f(x)與y軸的交點。b1為迴歸函式的斜率,也可以理解為黑線的傾斜程度。

圖中紅色方塊表示的是xi-f(xi)的資料集,該資料集都落在直線上。例如x=5,f(5)=8.33

有了以上兩組資料後,便可以計算殘差,即yi-f(xi)=yi-b0-b1xi,i=1,…,n。在圖中殘差可表示為綠點到紅方塊的距離,用虛線表示。在執行線性迴歸時,事實上是在試著使圖中綠點到紅方塊的虛線總長度最小,即殘差最小。

多變數線性迴歸

假如僅有兩個自變數,迴歸函式可以寫成f(x1,x2)=b0+b1x1+b2x2。這個函式在三維空間中是一個平面。

這種迴歸的目標也是確定權重b0,b1,b2,即使得這個平面f(x1,x2)與實際觀測的數值間的距離最小,即最小化SSR。

兩個自變數以上的多變數回歸也與上述的方法類似。例如f(x1,…,xr)=b0+b1x1+···+brxr,這個方程中有r+1個權重bi需要估計。

多項式迴歸

可以將多項式迴歸視為線性迴歸的一般情形。自變數與因變數的多項式依賴關係可以透過迴歸得到多項式迴歸函式。

換句話說,迴歸函式f中,除了包含線性部分如b1x1,還可以加入非線性部分如b2x12,b3x13甚至b4x1x2,b5x12x2等等。

單變數多項式迴歸是其中最簡單的一種形式,例如2階單變數多項式迴歸函式:f(x)=b0+b1x+b2x2

現在只需要透過最小二乘法計算b0,b1,b2即可。與上一節的線性迴歸函式f(x1,x2)=b0+b1x1+b2x2相比較,兩種迴歸非常相似,並且都需要估計係數b0,b1,b2。因此,只需要將高階項如x2視為因變數,多項式迴歸與普通的線性迴歸是一樣的。

在2階多變數的情形中,迴歸函式可以像這樣:f(x1,x2)=b0+b1x1+b2x2+b3x12+b4x1x2+b5x22,迴歸的方法和前面敘述的一樣。對五個自變數應用線性迴歸:x1,x2,x12x1x2,x22。迴歸的結果得到6個在SSR最小化時的權重:b0,b1,b2,b3,b4,b5

欠擬合和過擬合

在執行多項式迴歸中,一個需要引起高度重視的問題是:關於多項式迴歸函式最優階數的選擇

階數的選擇並沒有明確的規則。它需要視情況而定。但是,在選擇階數的時候,需要關注兩個問題:欠擬合與過擬合。

當模型無法準確的捕獲資料之間的依賴關係時導致欠擬合,這通常是由於模型過於簡單導致的。欠擬合的模型在用現有資料進行迴歸時會有較低的決定係數R2,同時它的預測能力也不足

當資料和隨機波動性都擬合到模型中時會導致過擬合。換種說法,就是模型和現有資料契合程度過高了。高度複雜的模型一般都有很多自變數,這通常容易導致模型過擬合。將現有資料擬合這種模型的時候一般會有很高的決定係數R2。但是在使用新資料時,模型可能會表現出很弱的預測能力和低決定係數R2

下圖分別展示了欠擬合,好的擬合,和過擬合的模型:


資料科學必備基礎之線性迴歸


左上圖是個低決定係數R2的線性迴歸。這條直線沒有考慮當x變化時觀測值yi的劇烈變動。這也許是一個欠擬合的例子。

右上圖是個2階多項式迴歸。在這個迴歸模型中,2階也許是擬合該組資料的最優階數。這個迴歸模型的R2還不錯,並且可以很好地勾勒出了資料的趨勢。

左下圖是個3階的多項式迴歸。迴歸模型的R2比右上圖的要高。對於現有的資料,這個模型比右上圖的模型要好。但是,它可能有一些過擬合的跡象,尤其是在x接近60的時候,現有的觀測資料並沒有顯示出迴歸模型中給出的下降趨勢。

最後,在右下圖中,這是一個完全擬合,6個觀測值都落在迴歸模型的曲線上。這可能是我們想要的迴歸模型,但是在多數情況下,這是一種過擬合的模型。它可能會有糟糕的預測能力。在這例中,可以看到曲線在x大於50時候毫無徵兆的急速下降,並且在x=60附近歸零。這種曲線可能是過度學習或者過度擬合現有資料導致的。


04

在Python中實現線性迴歸

關於線性迴歸Python庫

numpy是一款基礎Python包,它可以在單維和多維陣列上進行高效的操作。它還提供很多數值計算方法,同時也是開源庫。

如果不熟悉numpy,可以閱讀官方的numpy使用指南,此外還有numpy相關的文章會在附錄中提供。網路上還有大量的numpy資料可供搜尋。

scikit-learn是一個被大量使用者用於機器學習的Python庫,該庫是基於numpy和其他庫建立的。它可以執行資料預處理,降維,迴歸,分類和聚類等功能。和numpy一樣,scikit-learn也是開源的。

可以在scikit-learn官網上檢視Generalized Linear Models頁面瞭解更多關於線性模型的內容以進一步瞭解該庫的原理。

statsmodels庫可以實現一些在做線性迴歸時scikit-learn中沒有的功能。它也是一個不錯的Python庫,主要用於統計模型的估計,模型的測試等等。同時它也是開源的。

基於sklearn的單變數線性迴歸

這裡從最簡單的單變數線性迴歸開始,手把手教你在Python中做線性迴歸模型,整個過程大致分五個基本步驟:

① 匯入相關庫和類

② 讀取相關資料,對資料進行適當處理

③ 建立迴歸模型並對現有資料樣本進行擬合

④ 檢視模型的擬合結果,判斷模型是否合理

⑤ 用模型進行預測

這幾個步驟在大多數迴歸方法中基本是通用的。

Step 1:匯入相關庫和類

首先匯入pandastushare庫和sklearn.linear_model中的LinearRegression類:

import pandas as pd
import numpy as np
import tushare as ts
from sklearn.linear_model import LinearRegression

當然,也可以用numpy庫。pandas和numpy的資料型別可以與sklearn無縫對接,並且高效,可讀性強而類sklearn.linear_model.LinearRegression提供了執行線性迴歸的相關功能。

Step 2:讀取資料

這一步主要是讀取資料,並進行必要的預處理。這裡採用了tushare.pro的個股和指數日收益率資料,並對資料進行了適當的處理。

# 使用tushare獲取資料
ts.set_token('your token')
pro = ts.pro_api()


# 獲取滬深300與中國平安股票2018年全年以及2019年前7個交易日的日線資料
index_hs300 = pro.index_daily(ts_code='399300.SZ', start_date='20180101', end_date='20190110')
stock_000001 = pro.query('daily', ts_code='000001.SZ', start_date='20180101', end_date='20190110', fields='ts_code, trade_date ,pct_chg')

# 保留用於迴歸的資料,資料對齊
# join的inner引數用於剔除指數日線資料中中國平安無交易日的資料,比如停牌。
# 同時保留中國平安和指數的日收益率並分別命名兩組列名y_stock, x_index,確定因變數和自變數
index_pct_chg = index_hs300.set_index('trade_date')['pct_chg']
stock_pct_chg = stock_000001.set_index('trade_date')['pct_chg']
df = pd.concat([stock_pct_chg, index_pct_chg], keys=['y_stock''x_index'], join='inner', axis=1sort=True)


# 選中2018年的x,y資料作為現有資料進行線性迴歸
df_existing_data = df[df.index < '20190101']

# 注意:自變數x為pandas.DataFrame型別,其為二維資料結構,注意與因變數的pandas.Series資料結構的區別。也可以用numpy.array的.reshape((-1, 1))將資料結構改為二維陣列。
x = df_existing_data[['x_index']]
y = df_existing_data['y_stock']

經過簡單處理,現在有了自變數x和因變數y。需要注意的是,sklearn.linear_model.LinearRegression類需要傳入的自變數是一組二維陣列結構,因此可以透過numpy的.reshape()方法將一維陣列變更為二維陣列,另一種方法就是使用pandas的DataFrame格式。

現在自變數和因變數看起來像這樣:

print(x.head())
print(y.head())

            x_index
trade_date         
20180102     1.4028
20180103     0.5870
20180104     0.4237
20180105     0.2407
20180108     0.5173


trade_date
20180102    3.01
20180103   -2.70
20180104   -0.60
20180105    0.38
20180108   -2.56
Namey_stockdtypefloat64

可以透過.shape引數看到x有兩個維度(243, 1),而y是一維陣列,(243, )。

Step 3:建立模型並擬合資料

第三步,建立線性迴歸模型,並用中國平安和滬深300的2018年全年日線收益率資料進行擬合。先建立一個LinearRegression的例項用於建立迴歸模型:

model = LinearRegression()

這裡建立了一個預設引數的LinearRegression類例項。可以為LinearRegression傳入以下幾個可選引數:

fit_intercept布林值,預設為True。True為需要計算截距b0,False為不需要截距,即b0=0。

normalize:布林值,預設為False。引數作用是將資料標準化,False表示不做標準化處理。

copy_x:布林值,預設為True。True為保留x陣列,False為覆蓋原有的x陣列。

n_jobs:整型或者None,預設為None。這個引數用於決定平行計算執行緒數。None表示單執行緒,-1表示使用所有執行緒。

模型建立以後,首先需要在模型上呼叫.fit()方法:

model.fit(xy)

用現有的資料集x-y,透過.fit()方法,可以計算出線性迴歸方程的最優權重b0和b1。換句話說,資料集x-y透過.fit()方法擬合了迴歸模型。.fit()方法輸出self,也就是輸出model模型本身。可以將以上兩個語句簡化為:

model = LinearRegression().fit(x, y)

Step 4:輸出結果

資料擬合模型後,可以輸出結果看看模型是否滿意,以及是否可以合理的解釋資料的趨勢或者依賴關係。

可以透過呼叫.score()方法來獲取model的R2

model = LinearRegression().fit(x, y)
r_sq = model.score(x, y)
print('coefficient of determination: ', r_sq)

coefficient of determination:  0.5674052646582468

同時,也可以透過調取model的.intercept_.coef_屬性獲取截距b0和斜率b1

print('intercept:', model.intercept_)
print('slope:', model.coef_)

intercept: 0.01762302497463801
slope: [1.18891167]

需要注意的是,.intercept_是一個numpy.float64數值,而.coef_是一個numpy.array陣列。

b0=0.017,表示自變數滬深300的日收益率為0時,因變數中國平安股票的日收益率大約為 0.017%。而b1=1.189,表示滬深300的日收益率增加1%時,中國平安股票的日收益率增加1.189%。反之亦然。

在上面的資料載入中,因變數y也可以是二維陣列,這種情形下的輸出結果和以上類似:

y = pd.DataFrame(y)
new_model = LinearRegression().fit(x, y)
print('intercept:', new_model.intercept_)
print('slope:', new_model.coef_)

intercept: [0.01762302]
slope: [[1.18891167]]

可以看到此次輸出結果中,b0變為一維numpy.array陣列,而不是之前的float64數值。而b1輸出一組二維numpy.array陣列。

Step 5:模型預測

獲得了比較合理的模型後,可以用現有資料或者新資料來進行預測。

本例中使用了滬深300指數2019年前7個交易日的收益率資料對中國平安股票進行預測,由於本文主旨是介紹線性迴歸的原理與相關使用方法,模型的預測效果本文中不做探討。當然,也可以用2018年的現有的自變數資料獲取迴歸模型的資料f(x)。

要獲得預測結果,可以使用.predict()方法:

# 選2019年前7個交易日的資料作為新資料進行預測
df_new_data = df[df.index > '20190101']

# 將自變數x,資料型別為DataFrame
new_x = df_new_data[['x_index']]

# 預測
y_pred = model.predict(new_x)
print('predicted response:', y_pred, sep='\n')

predicted response:
[-1.60619253 -0.17022502  2.86601759  0.73929241 -0.23930079  1.21806713
 -0.20601126]

.predict()方法需要將新的自變數傳入該方法中以獲得預測值。也可以用下面這種方法:

b0, b1 = model.intercept_, model.coef_
y_pred = b0 + b1 * new_x
y_pred.columns = ['y_pred']

print('predicted response:', y_pred, sep='\n')

predicted response:
              y_pred
trade_date          
20190102   -1.606193
20190103   -0.170225
20190104    2.866018
20190107    0.739292
20190108   -0.239301
20190109    1.218067
20190110   -0.206011

這種方式可以將模型透過函式形式直觀的表達出來。與.predict()的方法相比,這種方式的輸出結果是一個DataFrame型別,其為二維陣列結構。而.predict()僅輸出一組一維numpy.array陣列。

如果將new_x改為一維陣列,y_pred的輸出結果將與.predict()一致。

實踐中,迴歸模型通常可以用來做預測。這意味著可以將一些新的自變數資料x傳入擬合好的模型中來計算結果。比如假設一組滬深300日收益率資料[0, 1, 2, 3, 4],透過model獲得中國平安股票的日收益率預測結果:

x_new = np.arange(5).reshape((-1, 1))
y_new = model.predict(x_new)
print(y_new)

[0.01762302 1.20653469 2.39544636 3.58435802 4.77326969]

基於sklearn的多變數線性迴歸

多元線性迴歸和基礎線性迴歸實現步驟基本一致。

Step 1:匯入相關庫和類

Step 2:讀取資料

本例使用了tushare.pro中的滬深300作為因變數,中金的細分金融指數和細分有色指數作為兩個自變數。

import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
import tushare as ts


# 使用tushare獲取資料
ts.set_token('your token')
pro = ts.pro_api()

# 獲取滬深300,中金細分有色指數,中金細分金融指數2018年全年日線資料
index_hs300 = pro.index_daily(ts_code='399300.SZ', start_date='20180101', end_date='20190110')
index_metal = ts.pro_bar(ts_code='000811.CSI', freq='D', asset='I', start_date='20180101', end_date='20190110')
index_finance = ts.pro_bar(ts_code='000818.CSI', freq='D', asset='I', start_date='20180101', end_date='20190110')

# 保留用於迴歸的資料,資料對齊
# join的inner引數用於剔除指數日線資料中交易日不匹配的資料
# 同時保留三個指數的日收益率並分別命名兩組列名y_zz300, x1_metal,x2_finance,確定因變數和自變數
index_hs300_pct_chg = index_hs300.set_index('trade_date')['pct_chg']
index_metal_pct_chg = index_metal.set_index('trade_date')['pct_chg']
index_finance_pct_chg = index_finance.set_index('trade_date')['pct_chg']
df = pd.concat([index_hs300_pct_chg, index_metal_pct_chg, index_finance_pct_chg], keys=['y_hs300''x1_metal''x2_finance'], join='inner', axis=1, sort=True)


# 選中2018年的x1,x2,y資料作為現有資料進行多變數線性迴歸
df_existing_data = df[df.index < '20190101']

# 提取多變數x1,x2,其資料型別為DataFrame。另外numpy需要的資料結構可以透過np.array(x)檢視,也是二維陣列。
x = df_existing_data[['x1_metal''x2_finance']]
y = df_existing_data['y_hs300']

因變數和自變數如下

print(x.head())
print(y.head())

            x1_metal  x2_finance
trade_date                      
20180102      0.8574      1.6890
20180103      1.0870     -0.0458
20180104      0.9429     -0.3276
20180105     -0.7645      0.1391
20180108      2.2010      0.0529

trade_date
20180102    1.4028
20180103    0.5870
20180104    0.4237
20180105    0.2407
20180108    0.5173
Namey_hs300dtypefloat64

DataFrame型別可以非常直觀的列出多個自變數,語法也相當簡潔,最重要的是該資料型別可以直接應用於LinearRegression類。當然,這裡也可以使用numpy.array的資料型別。

Step 3:建立模型並擬合資料

和前一個例子一樣,建立模型並資料擬合:

model = LinearRegression().fit(x, y)

Step 4:輸出結果

一樣的方式,透過model的.score().intercept_.coef_屬性獲取迴歸模型的決定係數,截距b0以及權重b1,b2

model = LinearRegression().fit(x, y)
r_sq = model.score(x, y)
print('coefficient of determination:', r_sq)
print('intercept:', model.intercept_)
print('slope:', model.coef_)


coefficient of determination: 0.8923741832489497
intercept: 0.0005333314305123599
slope: [0.30529472 0.64221632]

輸出結果可以看出,截距b0大致為0,即當x1=x2=0時,f(x1, x2)=0。更多地,當x1增加1時,相應的f(x1, x2)增加大約0.305,而當x2增加1時,相應的f(x1, x2)增加大約0.642。反之亦然。

Step 5:模型預測

預測方式與前面也是一致的:

# 選中2019年前7個交易日的資料作為新資料進行預測
df_new_data = df[df.index > '20190101']
# 取自變數x1,x2
new_x = df_new_data[['x1_metal''x2_finance']]
# 預測
y_pred = model.predict(new_x)
print('predicted response:', y_pred, sep='\n')

predicted response:
[-1.4716842   1.17891556  2.55034238  0.17701977 -0.81283822  0.57159138
 -0.19729323]

也可以用下面這種方式:

y_pred = model.intercept_ + (model.coef_ * new_x).sum(axis=1)
print('predicted response:', y_pred, sep='\n')

predicted response:
trade_date
20190102   -1.471684
20190103    1.178916
20190104    2.550342
20190107    0.177020
20190108   -0.812838
20190109    0.571591
20190110   -0.197293
dtype: float64

這種方式就是將每列自變數x1,x2乘上其對應的權重b1,b2,再加上一個固定的截距b0,得到預測的因變數。

基於sklearn的多項式迴歸

執行多項式迴歸步驟與前面也相似,只是需要增加一組經過轉換的自變數作為非線性項,如x2

Step 1:匯入相關庫和類

除了需要匯入pandas和tushare,sklearn.linear_model.LinearRegression外,還需要匯入sklearn.preprocessing中的PolynomialFeatures類用來處理非線性項。

import pandas as pd
import tushare as ts
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

Step 2:讀取資料

本例就採用tushare.pro上面的滬深300和中金細分金融指數資料。

# 獲取滬深300,中金細分金融指數2018年全年日線資料及2019年前7個交易日資料
index_hs300 = pro.index_daily(ts_code='399300.SZ', start_date='20180101', end_date='20190110')
index_finance = ts.pro_bar(ts_code='000818.CSI', freq='D', asset='I', start_date='20180101', end_date='20190110')

index_hs300_pct_chg = index_hs300.set_index('trade_date')['pct_chg']
index_finance_pct_chg = index_finance.set_index('trade_date')['pct_chg']
df = pd.concat([index_hs300_pct_chg, index_finance_pct_chg], keys=['y_hs300''x_finance'], join='inner', axis=1, sort=True)

df_existing_data = df[df.index < '20190101']

x = df_existing_data[['x_finance']]
y = df_existing_data['y_hs300']

Step 3:自變數轉換

這個步驟是多項式迴歸需要執行的步驟。由於多項式迴歸中有x2項,因此需要將x陣列轉換成x2並作為新的列。

有很多種轉換方式,比如使用numpy中的insert()方法,或者pandas的DataFrame新增x2列。本例中使用的是PolynomialFeatures類:

transformer = PolynomialFeatures(degree=2, include_bias=False)

transformer是PolynomialFeatures類的例項,用於對自變數x進行轉換。

PolynomialFeatures類中有以下幾個可選引數:

degree:整型,預設為2。用於決定線性迴歸模型的階數。

interaction_only:布林值,預設為False。如果指定為True,那麼就不會有xi2項,只有如xixj的交叉項。

include_bias:布林值,預設為True。此引數決定是否將截距項新增為迴歸模型中的一項,即增加一列值為1的列。False表示不新增。

先將自變數x資料傳入轉換器transformer中:

transformer.fit(x)

傳入transformer之後,就可以使用.transform()方法獲得轉換後的資料x和x2

x_ = transformer.transform(x)

當然也可以用.fit_transform()來替代以上三個語句:

x_ = PolynomialFeatures(degree=2, include_bias=False).fit_transform(x)

轉換後的陣列效果如下:

print(x_[:5])

[[ 1.6890000e+00  2.8527210e+00]
 [-4.5800000e-02  2.0976400e-03]
 [-3.2760000e-01  1.0732176e-01]
 [ 1.3910000e-01  1.9348810e-02]
 [ 5.2900000e-02  2.7984100e-03]]

可以看到新的資料中有兩列:一列是原始的自變數x的資料,另一列是x2項。

Step 4:建立模型並擬合資料

資料經過轉換後,模型的建立和擬合的方式與前面一樣:

model = LinearRegression().fit(x_, y)

【注】模型中第一個引數變為經過處理後的資料x_,而不是原始的自變數x陣列。

Step 5:輸出結果

r_sq = model.score(x_, y)
print('coefficient of determination:', r_sq)
print('intercept:', model.intercept_)
print('coefficients:', model.coef_)

.score()返回R2=0.813,截距項b0大約為-0.014,x項的權重b1約為0.862,x2項的權重b2約為-0.014。

另外一種迴歸方法,也可以透過向轉換器中傳入不同的引數獲得相同的迴歸結果:

x_ = PolynomialFeatures(degree=2, include_bias=True).fit_transform(x)

這裡將include_bias引數設定為True,其他的依然是預設引數。這種方法和step3中的相比,x_中會多出值為1的列,這列對應的是迴歸模型中的截距項:

print(x_[:5])

[[ 1.0000000e+00  1.6890000e+00  2.8527210e+00]
 [ 1.0000000e+00 -4.5800000e-02  2.0976400e-03]
 [ 1.0000000e+00 -3.2760000e-01  1.0732176e-01]
 [ 1.0000000e+00  1.3910000e-01  1.9348810e-02]
 [ 1.0000000e+00  5.2900000e-02  2.7984100e-03]]

這組資料中,第一列全為1,第二列為x,第三列為x2由於截距項已經包含在新的資料x_中了,因此在建立LinearRegression類時應該忽略截距項,即傳入引數fit_intercept=False

model = LinearRegression(fit_intercept=False).fit(x_, y)

至此,模型和資料x_是匹配的。模型結果如下:

r_sq = model.score(x_, y)
print('coefficient of determination:', r_sq)
print('intercept:', model.intercept_)
print('coefficients:', model.coef_)

coefficient of determination: 0.8133705974408868
intercept: 0.0
coefficients: [-0.01395755  0.86215933 -0.01444605]

可以看到.intercept_為0,但實際上截距b0已經包含在.coef_中了。其他結果都是一樣的。

Step 6:模型預測

同樣的,模型預測只需要呼叫.predict()方法,但是需要注意的是,傳入的新資料也需要做轉換:

# 選中2019年前7個交易日的資料作為新資料進行預測
df_new_data = df[df.index>'20190101']

# 取自變數x
new_x = df_new_data[['x_finance']]

# 自變數轉換
new_x_ = PolynomialFeatures(degree=2, include_bias=True).fit_transform(new_x)

# 預測
y_pred = model.predict(new_x_)
print('predicted response:', y_pred, sep='\n')

predicted response:
[-1.32738808  0.98948547  2.38535216 -0.30448932 -0.40485458  0.78721444
 -0.23448668]

已經瞭解到了,多項式迴歸模型的預測方式與前面兩例基本一樣,只需要將用於預測的資料做下轉換就好了。

如果有多個自變數,執行的過程也是一樣的。這裡使用了和前面提到的多變數回歸中相同的三組指數日收益率資料:

x = df_existing_data[['x1_metal', 'x2_finance']]
y = df_existing_data['y_hs300']


# step3 資料轉換
x_ = PolynomialFeatures(degree=2, include_bias=False).fit_transform(x)

# step4 建立模型並擬合資料
model = LinearRegression().fit(x_, y)

r_sq = model.score(x_, y)
intercept, coefficients = model.intercept_, model.coef_

# step5:模型預測

# 選中2019年前7個交易日的資料作為新資料進行預測
df_new_data = df[df.index > '20190101']
new_x = df_new_data[['x1_metal', 'x2_finance']]

# 自變數轉換
new_x_ = PolynomialFeatures(degree=2, include_bias=False).fit_transform(new_x)

# 預測
y_pred = model.predict(new_x_)

模型結果和預測如下:

print('coefficient of determination:', r_sq)
print('intercept:', intercept)
print('predicted response:', y_pred, sep='\n')

coefficient of determination: 0.8960573740272783
intercept: 0.020518382683722233
predicted response:
[-0.18180669  0.55726287 -0.79432939  0.17404118  2.57114362  1.23813068
 -1.45582482]

此結果的迴歸函式為f(x1, x2)=b0+b1x1+b2x2+b3x12+b4x1x2+b5x22。同時可以看到多項式迴歸的決定係數會高於多變數線性迴歸。決定係數R2高可能說明這是個不錯的迴歸模型。

但是,在實踐中,如此複雜的模型產生的R2,也可能是過擬合導致的。如果想要檢驗迴歸模型的表現,需要使用新的資料測試模型,即用樣本外的資料進行模型訓練。

statsmodels庫的線性迴歸應用

我們也可以使用statsmodels庫實現線性迴歸,這個庫的優點是可以輸出迴歸模型更詳細的結果。使用statsmodels的方法也和scikit-learn類似。

Step 1:匯入相關庫和類

import pandas as pd
import tushare as ts
import statsmodels.api as sm

除了pandas,還需要匯入statsmodels.api模組。

Step 2:輸出結果

x = df_existing_data[['x1_metal''x2_finance']]
y = df_existing_data['y_hs300']

如果用statsmodels計算截距,需要在x中新增全為1的列。statsmodels.api中有新增截距項的方法.add_constant(),該方法在陣列中的第一列新增了值為1的列:

x = sm.add_constant(x)

print(x.head())

            const  x1_metal  x2_finance
trade_date                             
20180102      1.0    0.8574      1.6890
20180103      1.0    1.0870     -0.0458
20180104      1.0    0.9429     -0.3276
20180105      1.0   -0.7645      0.1391
20180108      1.0    2.2010      0.0529

修改後的x現在有三列,第一列是新的截距項,另外兩列是原始資料。

Step 3:建立模型並擬合資料

迴歸模型是在statsmodels.regression.linear_model.OLS類建立的一個例項,OLS為最小二乘法:

model = sm.OLS(y, x)

【注】.OLS()中第一個引數是因變數,後面一個是自變數。其中還有幾個可選引數,可以參考statsmodels的官方文件。

模型建立後進行擬合:

results = model.fit()

透過呼叫.fit(),可以獲得模型結果results,results是statsmodels.regression.linear_model.RegressionResultsWrapper類的例項。模型結果中包含了大量關於這次迴歸模型的資訊。

Step 4:獲取結果

results中包含了迴歸模型的詳細資訊。提取結果的方法可以呼叫.summary():

print(results.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                y_zz300   R-squared:                       0.892
Model:                            OLS   Adj. R-squared:                  0.891
Method:                 Least Squares   F-statistic:                     995.0
Date:                Sat, 20 Jul 2019   Prob (F-statistic):          6.76e-117
Time:                        17:35:37   Log-Likelihood:                -146.32
No. Observations:                 243   AIC:                             298.6
Df Residuals:                     240   BIC:                             309.1
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.0005      0.029      0.019      0.985      -0.056       0.057
x1_metal       0.3053      0.023     13.385      0.000       0.260       0.350
x2_finance     0.6422      0.026     24.357      0.000       0.590       0.694
==============================================================================
Omnibus:                        0.419   Durbin-Watson:                   1.890
Prob(Omnibus):                  0.811   Jarque-Bera (JB):                0.224
Skew:                          -0.055   Prob(JB):                        0.894
Kurtosis:                       3.101   Cond. No.                         2.19
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

這張表提供的非常全面的迴歸結果。裡面可以找到包括R2,b0,b1和b2在內的資訊。同時也可以從上表中提取想要的值:

print('coefficient of determination:', results.rsquared)
print('adjusted coefficient of determination:', results.rsquared_adj)
print('regression coefficients:', results.params)

coefficient of determination: 0.8923741832489497
adjusted coefficient of determination: 0.8914773014426909
regression coefficients: 
const         0.000533
x1_metal      0.305295
x2_finance    0.642216
dtype: float64

以下是獲取線性迴歸的屬性的方法:

.rsquared:代表R2

.rsquared_adj:代表修正R2(經過自變數個數的調整後的R2)

.params:返回迴歸模型引數

可以對比scikit-learn中同類的迴歸模型,得到的模型結果都是一樣的。

Step 5:預測模型

可以使用.fittedvalues或者.predict()方法獲取模型產生的因變數f(x1,x2):

print('predicted response:', results.fittedvalues.head(), sep='\n')

predicted response:
trade_date
20180102    1.346996
20180103    0.302975
20180104    0.078006
20180105   -0.143532
20180108    0.706460
dtype: float64

print('predicted response:', results.predict(x.head()), sep='\n')

predicted response:
trade_date
20180102    1.346996
20180103    0.302975
20180104    0.078006
20180105   -0.143532
20180108    0.706460
dtype: float64

或者使用新資料進行預測:

print('predicted response:', results.fittedvalues.head(), sep='\n')

# 選中2019年前7個交易日的資料作為新資料進行預測
df_new_data = df[df.index>'20190101']

# 取自變數x1,x2
new_x = df_new_data[['x1_metal', 'x2_finance']]

# 新增截距項
new_x = sm.add_constant(new_x)

# 預測
new_y = results.predict(new_x)
print(new_y)

trade_date
20190102   -1.471684
20190103    1.178916
20190104    2.550342
20190107    0.177020
20190108   -0.812838
20190109    0.571591
20190110   -0.197293
dtype: float64

關於線性迴歸的不適用性情形

可以對比scikit-learn中同類的迴歸模型,得到的模型結果都是一樣的。

在資料分析的過程中,有時線性迴歸並不適用,尤其是面對一些比較複雜的非線性模型時。線性迴歸不適用的時候,還有很多其他方法可以作為備選方案,比如支援向量機,決策樹,隨機森林以及神經網路。

這些方法都有Python庫可用。而且大多數是免費開源的。這也是Python為什麼是機器學習主要語言之一的原因。

而scikit-learn庫也提供了除了線性迴歸以外的其他方法,而且用法和本文中的使用方法非常類似。該庫包含了支援向量機,決策樹,隨機森林等等演算法,這裡不過多展開了。這些演算法都可以透過.fit(),.predict(),.score()進行呼叫。


05

總結

本文介紹了什麼是線性迴歸已經如何在Python中實現,可以利用pandas對資料進行處理,pandas直觀高效的處理資料,並且可以與scikit-learn, statsmodels庫實現無縫銜接。

線性迴歸可以透過兩種方式實現:① scikit-learn:如果不需要回歸模型的詳細結果,用sklearn庫是比較合適的。② statsmodels:用於獲取迴歸模型詳細統計結果。兩個庫都可以做進一步的探索。更詳細的使用方法可以參考相關官方文件。

在Python中執行線性迴歸時,可以按照以下步驟執行:

① 匯入需要的庫和類。

② 提供資料,並且進行適當的預處理。

③ 建立迴歸模型,用現有資料進行擬合。

④ 檢驗模型是否合適,比如:欠擬合或者過擬合。

⑤ 應用模型進行預測。

【參考連結】

【1】

【2】

【3】

【4】

來自 “ ITPUB部落格 ” ,連結:http://blog.itpub.net/69908612/viewspace-2675167/,如需轉載,請註明出處,否則將追究法律責任。

相關文章