從零開始學Python數據分析【25】--嶺回歸及LASSO回歸(實戰部分)

作者：劉順祥

個人微信公眾號：每天進步一點點2015

前文傳送門：

從零開始學Python數據分析【1】--數據類型及結構

從零開始學Python數據分析【2】-- 數值計算及正則表達式

從零開始學Python數據分析【3】-- 控制流與自定義函數

從零開始學Python數據分析【4】-- numpy

從零開始學Python數據分析【5】-- pandas(序列部分)

從零開始學Python數據分析【6】-- pandas(數據框部分01)

從零開始學Python數據分析【7】-- pandas(數據框部分02)

從零開始學Python數據分析【8】-- pandas(數據框部分03)

從零開始學Python數據分析【9】-- pandas(數據框部分04)

從零開始學Python數據分析【10】-- matplotlib(條形圖)

從零開始學Python數據分析【11】-- matplotlib(餅圖)

從零開始學Python數據分析【12】-- matplotlib(箱線圖)

從零開始學Python數據分析【13】-- matplotlib(直方圖)

從零開始學Python數據分析【14】-- matplotlib(折線圖)

從零開始學Python數據分析【15】-- matplotlib(散點圖)

從零開始學Python數據分析【16】-- matplotlib(雷達圖)

從零開始學Python數據分析【17】-- matplotlib(面積圖)

從零開始學Python數據分析【18】-- matplotlib(熱力圖)

從零開始學Python數據分析【19】-- matplotlib(樹地圖)

從零開始學Python數據分析【20】--線性回歸（理論部分）

從零開始學Python數據分析【21】--線性回歸（實戰部分）

從零開始學Python數據分析【22】--線性回歸診斷（第一部分）

從零開始學Python數據分析【23】--線性回歸診斷（第二部分）

從零開始學Python數據分析【24】--嶺回歸及LASSO回歸（理論部分）

前言

      在《從零開始學Python【24】--嶺回歸及LASSO回歸（理論部分）》一文中我們詳細介紹了關於嶺回歸和LASSO回歸的理論知識，其實質就是在線性回歸的基礎上添加了2範數和1範數的懲罰項。這兩個模型的關鍵點是找到一個合理的lambda係數，來平衡模型的方差和偏差，從而得到比較符合實際的回歸係數。本期是基於之前討論的理論部分，採用Python和R語言，完成對嶺回歸和LASSO回歸的實戰。文中所利用的數據集來自R語言ISLR包中的Hitters數據集，描述的是關於棒球運動員收入的相關信息，數據集和腳本可以之文末查看連結獲得。

嶺回歸

      原始數據集存在收入的缺失，我們不妨先把這樣的觀測刪除，同時，數據集中含有離散變量，需要將這些變量轉換為啞變量後方可建模，故第一步需要對原始數據集進行清洗。

# ===== Python3 =====

# 導入第三方包
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import Ridge,RidgeCV
from sklearn.linear_model import Lasso,LassoCV
from sklearn.metrics import mean_squared_error

# 讀取數據
df = pd.read_csv('Hitters.csv')
# 啞變量處理
dummies = pd.get_dummies(df[['League', 'Division', 'NewLeague']])
# 將原始數據集與啞變量數據合併起來
mydf = df.join(dummies)
# 缺失值刪除
mydf = mydf.dropna()
# 刪除不必要的變量(字符串變量和各啞變量中的一個變量)
mydf = mydf.drop([ 'League', 'Division', 'NewLeague', 'League_N', 'Division_W', 'NewLeague_N'], axis = 1)# 前5行展示mydf.head()

上面的數據集清洗完畢，展現的是乾淨數據的前5行信息，下面要基於這個數據集進行建模。建模之前還需要將數據集拆分為兩部分，一部分用於建模，另一部分用於模型的測試。

# 將數據集拆分成訓練集和測試集
predictors = list(mydf.columns)predictors.remove('Salary') # Salay變量為因變量，故需要排除

X_train, X_test, y_train, y_test = train_test_split(mydf[predictors],mydf['Salary'], train_size = 0.8, random_state = 1234 )

# 通過不確定的alphas值，生成不同的嶺回歸模型
alphas = 10**np.linspace(-3,3,100)ridge_cofficients = []

for alpha in alphas:    ridge = Ridge(alpha = alpha, normalize=True)    ridge.fit(X_train, y_train)    ridge_cofficients.append(ridge.coef_)

# 繪製alpha的對數與回歸係數的關係
# 中文亂碼和坐標軸負號的處理
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']plt.rcParams['axes.unicode_minus'] = False
# 設置繪圖風格
plt.style.use('ggplot')
plt.plot(alphas, ridge_cofficients)plt.xscale('log')plt.axis('tight')plt.title('alpha係數與嶺回歸係數的關係')plt.xlabel('Log Alpha')plt.ylabel('Cofficients')plt.show()

從上面的圖形結果來看，alpha在10附近時，所有的自變量係數基本趨於穩定（但也不能完全確定是這個值）。接下來，我們採用交叉驗證（CV）方法確定最佳的lambda值。

基於CV選擇lambda值

# 為了找到最佳的lambda值，我們採用交叉驗證方法

# 嶺回歸模型的交叉驗證
ridge_cv = RidgeCV(alphas = alphas, normalize=True, scoring='mean_squared_error', cv = 10)ridge_cv.fit(X_train, y_train)
# 取出最佳的lambda值ridge_best_alpha = ridge_cv.alpha_ridge_best_alpha

      

不出所料，得到的lambda值確實在10附近，這裡最佳的lambda值為10。下面，我們要基於這個最佳的lambda值進入嶺回歸模型的創建和模型驗證的階段。

# 基於最佳的lambda值建模
ridge = Ridge(alpha = ridge_best_alpha, normalize=True)ridge.fit(X_train, y_train)
# 嶺回歸係數
ridge.coef_

# 預測
ridge_predict = ridge.predict(X_test)
# 預測效果驗證
RMSE = np.sqrt(mean_squared_error(y_test,ridge_predict))RMSE

經過模型的驗證，得到的RMSE為319.9。接下來，我們利用同樣的邏輯，對比一下LASSO回歸模型的效果。

LASSO回歸

# 通過不確定的alphas值，生成不同的LASSO回歸模型
alphas = 10**np.linspace(-3,3,100)lasso_cofficients = []

for alpha in alphas:    lasso = Lasso(alpha = alpha, normalize=True, max_iter=10000)    lasso.fit(X_train, y_train)    lasso_cofficients.append(lasso.coef_)

# 繪製alpha的對數與回歸係數的關係# 中文亂碼和坐標軸負號的處理
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']plt.rcParams['axes.unicode_minus'] = False
# 設置繪圖風格
plt.style.use('ggplot')
plt.plot(alphas, lasso_cofficients)plt.xscale('log')plt.axis('tight')plt.title('alpha係數與LASSO回歸係數的關係')plt.xlabel('Log Alpha')plt.ylabel('Cofficients')plt.show()

從圖形結果來看，lambda值應該在1附近，此時LASSO回歸的係數也基本趨於穩定（但也不能完全確定是這個值）。同樣，我們利用CV方法，來尋找最佳的lambda值。

基於CV選擇lambda值

# LASSO回歸模型的交叉驗證
lasso_cv = LassoCV(alphas = alphas, normalize=True, cv = 10, max_iter=10000)lasso_cv.fit(X_train, y_train)
# 取出最佳的lambda值
lasso_best_alpha = lasso_cv.alpha_lasso_best_alpha

      

通過CV方法得到的lambda結果是0.23，這與與我們看圖得到的1這個值還是有一點差異的。下面，我們就基於交叉驗證得到的最佳lambda值重新構造LASSO回歸模型。

# 基於最佳的lambda值建模
lasso = Lasso(alpha = lasso_best_alpha, normalize=True, max_iter=10000)lasso.fit(X_train, y_train)
# 嶺回歸係數lasso.coef_

# 預測
lasso_predict = lasso.predict(X_test)
# 預測效果驗證
RMSE = np.sqrt(mean_squared_error(y_test,lasso_predict))RMSE

      

對LASSO回歸模型進行驗證發現得到的RMSE更小，說明LASSO回歸模型的擬合效果會更貼近於Hitters數據集的原貌。

      上面的內容是基於Python工具對嶺回歸模型和LASSO回歸模型的實戰，接下來，我們再利用R語言對上面的過程再作一次復現，希望對R語言感興趣的朋友能夠有幫助。

R語言對比

      由於上面的邏輯我們已經通過Python進行了一一說明，這裡就不再贅述，只給出R語言代碼僅供參考。

# 加載第三方包
library(caret)
library(glmnet)
library(ISLR)

# 啞變量處理
dummies <- dummyVars(~League+Division+NewLeague, data = Hitters)dummies <- predict(dummies, newdata = Hitters)
# 數據合併
Hitters_dummy <- cbind(Hitters, dummies)

# 刪除缺失值
Hitters_dummy <- na.omit(Hitters_dummy)
# 刪除不必要的變量
Hitters_dummy <- subset(Hitters_dummy,                       select = -c(League,Division,NewLeague,League.N,Division.W,NewLeague.N))head(Hitters_dummy)

# 構建訓練集和測試集set.seed(1)index <- sample(1:nrow(Hitters_dummy), size = 0.8*nrow(Hitters_dummy))train <- Hitters_dummy[index,]test <- Hitters_dummy[-index,] # 繪製lambda值與嶺回歸係數的關係
fit_ridge <-  glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 0)plot(fit_ridge,xvar = 'lambda',label = T)

# 嶺回歸的交叉驗證，確定最佳的lambda值
fit_ridge_cv <- cv.glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 0)best_lambda_ridge <- fit_ridge_cv$lambda.minbest_lambda_ridge

# 根據最佳lambda構建嶺回歸模型
fit_ridge <- glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 0)coeff_ridge <- predict(fit_ridge, s = best_lambda_ridge, type = 'coefficients')coeff_ridge

# 模型評估
pred_ridge <- predict(fit_ridge, s = best_lambda_ridge, newx = as.matrix(test[,-17]))RMSE <- sqrt(mean((test$Salary-pred_ridge)**2))RMSE

# 繪製lambda值與LASSO回歸係數的關係
fit_lasso <-  glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 1)plot(fit_lasso,xvar = 'lambda',label = T)

# 嶺回歸的交叉驗證，確定最佳的lambda值
fit_lasso_cv <- cv.glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 1)best_lambda_lasso <- fit_lasso_cv$lambda.minbest_lambda_lasso

# 根據最佳lambda構建嶺回歸模型
fit_lasso <- glmnet(x = as.matrix(train[,-17]), y = train[,17], alpha = 1)coeff_lasso <- predict(fit_lasso, s = best_lambda_lasso, type = 'coefficients')coeff_lasso

# 模型評估
pred_lasso <- predict(fit_lasso, s = best_lambda_lasso, newx = as.matrix(test[,-17]))RMSE <- sqrt(mean((test$Salary-pred_lasso)**2))RMSE

結語

      OK，今天關於嶺回歸和LASSO回歸的實戰部分就介紹到這裡，希望對數據挖掘或機器學習的朋友有所幫助，同時，也希望讀者能夠靜下心來好好的復現一遍。如果你有任何問題，歡迎在公眾號的留言區域表達你的疑問。同時，也歡迎各位朋友繼續轉發與分享文中的內容，讓更多的朋友學習和進步。

連結: https://pan.baidu.com/s/1qYOGuc8 密碼: xb3r

Python愛好者社區歷史文章大合集：

Python愛好者社區歷史文章列表（每周append更新一次）

福利：文末掃碼立刻關注公眾號，「Python愛好者社區」，開始學習Python課程：

關注後在公眾號內回復「課程」即可獲取：

0.小編的Python入門視頻課程！！！

1.崔老師爬蟲實戰案例免費學習視頻。

2.丘老師數據科學入門指導免費學習視頻。

3.陳老師數據分析報告製作免費學習視頻。

4.玩轉大數據分析！Spark2.X+Python 精華實戰課程免費學習視頻。

5.丘老師Python網絡爬蟲實戰免費學習視頻。