前言
本次項目是就某聯合循環發電廠的數據，運用線性回歸模型進行預測電能輸出，若文中出現錯誤的地方，還望指正，謝謝！
目錄
1. 數據來源及背景
2. 數據探索分析
3. 相關分析
4. 回歸分析
5. 多重共線性
6. 模型應用
正文
1. 數據來源及背景
數據來源: http://archive.ics.uci.edu/ml/machine-learning-databases/00294/
該數據集是收集於聯合循環發電廠的9568個數據點, 共包含5個特征: 每小時平均環境變量溫度（AT），環境壓力（AP），相對濕度（RH），排氣真空（V）和凈每小時電能輸出（PE）, 其中電能輸出PE是我們要預測的變量.
2. 數據探索分析
由於我們的數據是excel格式, 而pandas處理excel文件依賴xlrd, 因此我們首先需要安裝它.
1. 安裝xlrd
pip install xlrd
2. 讀取數據
excel文件中共有5個sheet, 我們讀取最后一個
import pandas as pd #第二參數代表要讀取的sheet, 0表示第一個, 1表示第二個..., pandas默認讀取第一個 df = pd.read_excel(r'D:\Data\CCPP.xlsx', 4)
3. 查看數據前3行/后3行
pd.set_option('display.max_rows', 6) df
數據維度9568行X5列, 均是數值型數據
4. 查看數據整體信息
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 9568 entries, 0 to 9567 Data columns (total 5 columns): AT 9568 non-null float64 V 9568 non-null float64 AP 9568 non-null float64 RH 9568 non-null float64 PE 9568 non-null float64 dtypes: float64(5) memory usage: 373.8 KB
無缺失數據, 且均為浮點型.
5. 描述性統計
pd.set_option('display.max_rows', None) df.describe()
溫度(AT): 范圍1.81--37.11, 均值19.65, 中位數20.35, 左偏分布
排氣真空(V): 范圍25.36--81.56, 均值54.31, 中位數52.08, 右偏分布
環境壓力(AP): 范圍992.89--1033.30, 均值1013.26, 中位數1012.94, 右偏分布
相對濕度(RH): 范圍25.56--100.16, 均值73.31, 中位數74.975, 左偏分布
電能輸出(PE): 范圍420.26--495.76, 均值454.37, 中位數451.55, 右偏分布
通過中位數和均值可以大致看出分布狀況, 但是對於具體的偏斜程度, 就需要用偏態系數去衡量
6. 偏態系數
偏態系數反映的是數據的偏斜程度, 若偏態系數大於1或小於-1, 則稱為高度偏態分布, 若在0.5~1或-1~-0.5之間, 稱為中度偏態分布, 否則, 稱為低度偏態分布
for i in df.columns: print('{}偏態系數: '.format(i), df[i].skew())
AT偏態系數為: -0.1363930494749227 V偏態系數為: 0.19852101136676173 AP偏態系數為: 0.2654446935825862 RH偏態系數為: -0.4318387491833329 PE偏態系數為: 0.3065094354204023
可以看出以上均為低度偏態分布. 我們來看看分布形狀的另一個度量--峰態系數
7. 峰態系數
峰態系數反映的是數據的扁平程度, 峰態常常是與標准正態分布對比, 若峰態系數為0, 則說明服從標准正態分布, 若大於0, 則說明比標准正態分布更尖, 稱為尖峰分布, 若小於0, 則稱為平峰分布.
for i in df.columns: print('{}峰態系數為: '.format(i), df[i].kurt())
3. 相關分析
import seaborn as sns import matplotlib.pyplot as plt #初始化 sns.set() #繪制分布矩陣 sns.pairplot(df) #保存圖片 plt.savefig('ccpp.png') plt.show()
#計算相關系數, 默認為皮爾遜相關系數 correlation = df.corr() correlation
皮爾遜相關系數反映的是線性相關程度, 其取值介於-1~1之間. 若取0, 說明兩變量無關, 若取+1, 說明完全正相關, 若取-1, 說明完全負相關.
可以將相關系數繪制成熱力圖
3. 熱力圖
#繪制熱力圖, 設置線寬, 最大值, 最小值, 線條白色, 顯示數值, 方形 sns.heatmap(correlation, linewidths=0.2, vmax=1, vmin=-1, linecolor='w', annot=True,square=True) plt.savefig('correlation.png') plt.show()
可以根據熱力圖的顏色來判斷關系強度: 自變量中AT和V與因變量PE為高度相關(0.8~1或-1~-0.8), 自變量AP與因變量PE為中度相關(0.5~0.8或-0.8~-0.5), 自變量RH與因變量PE為低度相關(0.3~0.5或-0.5~-0.3). 另外, 自變量中V和AT高度相關.
接下來我們開始回歸分析
4. 回歸分析
相關分析是確定是否存在關系以及關系強度, 而回歸分析則是確定它們到底是什么樣的關系, 建立具體的關系表達式.
由於我們的自變量有4個且各自變量與因變量均有不同程度的線性相關關系, 因此我們假設全部自變量與因變量滿足多元線性回歸模型:
回歸分析的目的就是確定這
, 
, 
, 
, 
5個模型參數, 我們這里直接采用sklearn庫的線性回歸模型(基於最小二乘法)去求解這個5參數.
在求解之前, 將數據集划分為訓練集和測試集.
1. 划分數據集
from sklearn.model_selection import train_test_split X, y = df[['AT', 'V', 'AP', 'RH']], df['PE'] #按照8:2的比例划分為訓練數據集和測試數據集 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2)
2. 構建模型
from sklearn.linear_model import LinearRegression #模型實例化 LR = LinearRegression() #訓練模型 LR.fit(X_train, y_train) print("截距: ", LR.intercept_) print("回歸系數: ", LR.coef_ )
截距: 445.62291957190826 回歸系數: [-1.97153297 -0.23188593 0.07065982 -0.15715584]
因此, 我們的線性回歸模型為: PE=455.62-1.97153297AT-0.23188593V+0.07065982AP -0.15715584RH
3. 模型評估
對於回歸模型, 常用的評估方法有: 平均絕對誤差(MAE), 均方誤差(MSE), 均方根誤差(RMSE), 多重判定系數(R2).
1) 平均絕對誤差MAE(mean_absolute_error)
其中, 
表示實際值, 
表示估計值. 平均絕對誤差可以有效避免誤差相互抵消的問題, 可以更好反映預測誤差的實際情況
from sklearn import metrics #分別對訓練集和測試集進行預測 y_train_pred = LR.predict(X_train) y_test_pred = LR.predict(X_test) #分別計算訓練集和測試集的平均絕對誤差 y_train_mae = metrics.mean_absolute_error(y_train, y_train_pred) y_test_mae = metrics.mean_absolute_error(y_test, y_test_pred) print('訓練集MAE: ', y_train_mae) print("測試集MAE: ", y_test_mae)
2) 均方誤差MSE(mean_square_error)
#分別計算訓練集和測試集的均方誤差 y_train_mse = metrics.mean_squared_error(y_train, y_train_pred) y_test_mse = metrics.mean_squared_error(y_test, y_test_pred) print('訓練集MSE: ', y_train_mse) print("測試集MSE: ", y_test_mse)
訓練集MSE: 20.684348049573618 測試集MSE: 21.115806888804244
同樣是測試集比訓練集誤差大, 由於其數值的單位也是平方級別, 為了更好地描述對其開根號, 也就是均方根誤差.
3) 均方根誤差RMSE(root_mean_square_error)
from math import sqrt y_train_rmse = sqrt(y_train_mse) y_test_rmse = sqrt(y_test_mse) print('訓練集RMSE: ', y_train_rmse) print("測試集RMSE: ", y_test_rmse)
訓練集RMSE: 4.548004842738584 測試集RMSE: 4.5951938902296865
這樣就方便去描述了, 即用我們所建立的回歸模型來預測每小時凈電能輸出時, 平均的預測誤差為4.6MW, 與訓練集相比略大.
4) 多重判定系數R2 (r2_score)
其中, SST則為總誤差平方和, 即SST=SSR+SSE, SSR為回歸平方和, 反映的是能被回歸方程解釋的誤差; SSE為殘差平方和, 反映的是不能被回歸方程解釋的誤差. 也就是說, 多重判定系數反映的是被回歸方程解釋的誤差占總誤差的比例, 回歸平方和所占比重越大, 說明回歸方程擬合程度越好
#分別計算訓練集和測試集的多重判定系數 y_train_r2 = metrics.r2_score(y_train, y_train_pred) y_test_r2 = metrics.r2_score(y_test, y_test_pred) print('訓練集R2: ', y_train_r2) print("測試集R2: ", y_test_r2)
訓練集R2: 0.928600622449398 測試集R2: 0.9289130339796056
即在訓練集和測試集上總變差中能被回歸方程所解釋的比例分別為92.86%, 92.89%, 在測試集上比訓練集表現略好.
在sklearn中還有自帶的評估器, 不妨來看看.
#直接用訓練好的模型去評分 y_train_score = LR.score(X_train, y_train) y_test_score = LR.score(X_test, y_test) print('訓練集score: ', y_train_score) print("測試集score: ", y_test_score)
訓練集score: 0.928600622449398 測試集score: 0.9289130339796056
可以看得出與多重判定系數的結果是一致的, 說明sklearn自帶的評估器采用的是多重判定系數.
由於我們是假設全部自變量與因變量服從線性關系, 那么到底是否真正服從呢? 這就需要進行檢驗才能得知, 即進入下一步: 模型檢驗
4. 模型檢驗
首先, 我們需要檢驗模型是否服從線性關系,
1) 線性關系檢驗
由於該檢驗統計量是以殘差平方和(SSE)和回歸平方和(SSR)為基礎, 故除以它們各自的自由度來構造F檢驗統計量,
即服從分子自由度為k, 分母自由度為n-k-1的F分布, 其中k為自變量個數, n為樣本量. 進行假設檢驗主要分為5個步驟:
① 提出原假設和備擇假設: H0: 
, H1: 不全為0
② 確定顯著性水平: 
③ 選擇檢驗統計量: 
④ 建立決策准則:有兩種方式一種是拒絕域, 一種是P值.
a. 拒絕域: 即拒絕原假設的區域, 該拒絕域通過臨界點
(由於F統計量為二次型, 因此取值非負, 僅有一個拒絕域, 故為α) 來確定, 如果F檢驗統計量落入拒絕域, 則拒絕原假設從而接受備擇假設, 否則不能拒絕原假設.
b. P值: P值表示原假設發生的可能性(概率), 若可能性小於顯著性水平, 即小概率事件, 那么就拒絕原假設從而接受備擇假設, 否則, 不能拒絕原假設.
⑤ 下結論: 根據決策准則, 做出是否拒絕原假設的結論.
from pandas import Series from scipy.stats import f #將array轉為series格式 y_train_pred = Series(y_train_pred, index=y_train.index) #分別計算訓練數據上的SSR和SSE y_train_ssr = y_train_pred.apply(lambda x: (x - y_train.mean())**2).sum() y_train_sse = y_train.sub(y_train_pred).apply(lambda x: x**2).sum() #dn是SSR的自由度(自變量個數), df則是SSE的自由度 dn, df = 4, y_train.shape[0]-dn-1 #計算F值 y_train_f = (y_train_ssr/dn) / (y_train_sse/df) #計算p值 p = f.sf(y_train_f, dn, df) #計算0.05顯著性水平下臨界值 cr_value = f.isf(0.05, dn, df) print('訓練數據集F值: ', y_train_f) print('0.05顯著性水平下臨界值: ', cr_value) print('訓練數據集P值: %.20f'% p)
訓練數據集F值: 24870.196368594174 0.05顯著性水平下臨界值: 2.3730935370191646 訓練數據集P值: 0.00000000000000000000
F值大於臨界值, 拒絕原假設H0 ( P值小於顯著性水平, 同樣可以拒絕原假設), 表明全部自變量與因變量的線性關系是顯著的, 但這並不意味着每個自變量都和因變量的線性關系顯著, 因此還需要進行回歸系數(自變量)檢驗.
2) 回歸系數檢驗
對每個回歸系數構造t檢驗統計量:
其中, 
表示第i個回歸系數的估計值(因為是根據最小二乘法估計的, 而不是真實值), 
表示第i個回歸系數的標准誤差. 則:
其中, MSE表示無偏殘差平方和, 即殘差平方和SSE除以它的自由度(n-k-1), Cii表示(XXT)-1逆矩陣i行第i列的元素. t檢驗過程同樣是五步驟:
① 提出原假設和備擇假設: H0: 
, H1: 不為0
② 確定顯著性水平:
③ 選擇檢驗統計量: 
④ 建立決策准則: 同樣有兩種方式, 不過t檢驗的拒絕域有兩個.
⑤ 下結論: 根據決策准則做出是否拒絕原假設的結論
from scipy.stats import t def get_tvalue(sse, df, matr, beta, i): '''計算t值''' mse = sse / df sbeta = sqrt(matr[i+1, i+1]* mse) t = beta / sbeta return t limit = t.isf(0.025, df) print('0.05顯著性水平下的臨界值: ', limit) X_train['B'] = 1 X_train = X_train.reindex(columns=['B', 'AT', 'V', 'AP', 'RH']) #轉成矩陣 xm = np.mat(X_train) #計算(X'X)的逆矩陣 xmi = np.dot(xm.T, xm).I index, betas = range(4), LR.coef for i, beta in zip(index, betas): tvalue = get_tvalue(y_train_sse, df, xmi, beta, i) pvalue = t.sf(abs(tvalue), df)*2 print('beta{0}的t值: '.format(i+1), tvalue) print('beta{0}的p值: '.format(i+1), pvalue)
0.05顯著性水平下的臨界值: 2.085963447265837 beta1的t值: -5.898890146924707 beta1的p值: 9.049074948455422e-06 beta2的t值: -1.4562620735727096 beta2的p值: 0.16084038317451793 beta3的t值: 0.34176894237398453 beta3的p值: 0.7360897563076247 beta4的t值: -1.734158669159414 beta4的p值: 0.09827823435663764
從以上結果中可以看到只有beta1(自變量AT)的t值的絕對值大於臨界值, 即只有beta1拒絕原假設, 也就是說4個自變量中只有beta1影響是顯著的, 那么為什么其他自變量不呢? 造成這種的原因可能是自變量之間高度相關, 這種問題稱為多重共線性.
5. 多重共線性
引用百度百科的定義
多重共線性是指線性回歸模型中的解釋變量之間由於存在精確相關關系或高度相關關系而使模型估計失真或難以估計准確。
根據定義可知, 多重共線性的判別方法是相關系數矩陣, 查看之前的相關系數矩陣發現在AT和V這兩個解釋變量(自變量)之間存在高度相關關系: 數值為0.84. 因此可以認為變量間存在多重共線性.
對於多重共線性的處理方法有: 刪除共線變量, 轉換模型的形式, 逐步回歸法, 嶺回歸法, 主成分分析法. 在這里我采用前兩種方法.
1. 刪除共線變量: 通常是刪除高度相關的變量, 保留重要的變量. 重要的變量是指通過t檢驗的變量, 於是我們刪除V這個解釋變量.
X1, y1 = df[['AT', 'AP', 'RH']], df['PE'] X1_train, X1_test, y1_train, y1_test = train_test_split(X1, y1, test_size = 0.2) #訓練模型, 預測以及計算均方根誤差和多重判定系數 LR1 = LinearRegression() LR1.fit(X1_train, y1_train) inter, co = LR1.intercept_, LR1.coef_ y1_train_pred = LR1.predict(X1_train) y1_train_rmse = sqrt(metrics.mean_squared_error(y1_train, y1_train_pred)) y1_train_score = LR1.score(X1_train, y1_train) print("回歸模型: PE={0}+{1}AT+{2}AP+{3}RH".format(inter,co[0], co[1], co[2])) print('訓練集RMSE: ', y1_train_rmse) print('訓練集擬合優度: ', y1_train_score) #計算F檢驗中0.05顯著性水平下的P值 y1_train_pred = Series(y1_train_pred, index=y1_train.index) y1_train_ssr = y1_train_pred.apply(lambda x: (x - y1_train.mean())**2).sum() y1_train_sse = y1_train.sub(y1_train_pred).apply(lambda x: x**2).sum() dn1, df1 = 3, y1_train.shape[0]-3-1 y1_train_f = (y1_train_ssr/dn1) / (y1_train_sse/df1) y1_p = f.sf(y1_train_f, dn1, df1) # cr_value = f.isf(0.05, dn, df) print('F檢驗 0.05顯著性水平下訓練數據的P值: %.20f'% y1_p) #計算t檢驗在0.05顯著性水平下的P值 def get_t1value(sse, df, matr, beta, i): mse = sse / df sbeta = sqrt(matr[i+1, i+1]* mse) t = beta / sbeta return t X1_train['B'] = 1 X1_train = X1_train.reindex(columns=['B', 'AT', 'AP', 'RH']) xm1 = np.mat(X1_train) xmi1 = np.dot(xm1.T, xm1).I index, betas = range(3), LR1.coef_ for i, beta in zip(index, betas): tvalue = get_t1value(y1_train_sse, df1, xmi1, beta, i) pvalue = t.sf(abs(tvalue), df1)*2 print('t檢驗 0.05顯著性水平下beta{0}的p值: '.format(i+1), pvalue)
回歸模型: PE=485.39627734198206+-2.374833922164532AT+0.030214741968456443AP+-0.20456036084843196RH 訓練集RMSE: 4.839764102773314 訓練集擬合優度: 0.9195738379930241 F檢驗 0.05顯著性水平下訓練數據的P值: 0.00000000000000000000 t檢驗 0.05顯著性水平下beta1的p值: 0.0 t檢驗 0.05顯著性水平下beta2的p值: 0.006605096641032901 t檢驗 0.05顯著性水平下beta3的p值: 0.0
均方根誤差比之前增大0.3, 擬合優度相比之前減少0.009, F檢驗通過, 且回歸系數均通過, 可見刪除共線變量還是可以有效解決多重共線性問題, 但是誤差相比之前增大0.2, 因此這種方法還是有一定局限性. 換另外一種方法試試.
2. 轉換模型形式: 將數據進行對數轉換.
df_log = np.log(df)
回歸模型: PE=2.8675586960921957+-0.05125235913516803AT+-0.05608345239278856V+0.5268863140584524AP+-0.0059335463238909605RH 訓練集RMSE: 0.010378041999902672 訓練集擬合優度: 0.9229456762426529 F檢驗 0.05顯著性水平下訓練數據的P值: 0.00000000000000000000 t檢驗 0.05顯著性水平下beta1的p值: 0.0 t檢驗 0.05顯著性水平下beta2的p值: 0.0 t檢驗 0.05顯著性水平下beta3的p值: 8.936059858590797e-108 t檢驗 0.05顯著性水平下beta4的p值: 1.7118836102396494e-20
同樣也都通過檢驗, 而且經過對數轉換后的均方根誤差更小. 對比可以發現進行對數轉換比刪除變量的擬合優度更好些. 我們接下來便是用轉換后的模型去預測了.
6. 模型應用
import matplotlib.pyplot as plt #用訓練好的模型去預測 y2_test_pred = LR2.predict(X2_test) y2_test_rmse = sqrt(metrics.mean_squared_error(y2_test, y2_test_pred)) # y1_test_rmse = sqrt(metrics.mean_squared_error(y1_test, y1_test_pred)) y2_test_score = LR2.score(X2_test, y2_test) print('測試集RMSE: ', y2_test_rmse) print('測試集擬合優度: ', y2_test_score) #繪制曲線 plt.figure(figsize=(12,8)) plt.plot(range(len(y2_test)), y2_test, 'g', label='test data') plt.plot(range(len(y2_test_pred)), y2_test_pred, 'r', label='predict data', linewidth=2, alpha=0.8) plt.legend() plt.savefig('tp.png') plt.show()
測試集RMSE: 0.01038585831215862 測試集擬合優度: 0.9223780289364194
可以看到在測試集上表現還是相當不錯的.
以上便是對線性回歸模型學習過程的總結, 若存在問題, 希望指出, 謝謝!
參考
網易雲課堂《吳恩達機器學習》