應用|使用正則化線性模型和XGboost對價格建模

編者按：房屋需求是每個家庭剛性需求，關乎到家庭的幸福生活。如何購買到合適的房子，瞭解影響房子價格的因素以及基於歷史資料對房屋價格做一些預測性探索與分析，既有趣味又有價值。本文利用線性回歸正則化技術和xgboost技術對房屋價格進行建模。我建立了Python語言微信群，定位：Python語言學習與實踐，想入群的朋友，請新增我的微信：luqin360。

我們想要建模房子的價格，我們知道價格取決於房子的位置，房子的面積，建成年限，翻新的年限，臥室的數量，車庫的數量等等。因此，這些因素促成了這種樣式——優質地段通常會導致更高的價格。然而，在同一區域內，面積相同的所有房子的價格並不完全相同。價格的變化就是噪聲。我們對價格建模的標的是對樣式進行建模，忽略噪聲。同樣的概念也適用於酒店房價的建模。

因此，首先，我們將對房價資料進行線性回歸的正則化技術。

資料集

一個很好的房屋價格資料集可以在這裡找到。

import warnings
def ignore_warn(*args, **kwargs):
    pass
warnings.warn = ignore_warn
import numpy as np 
import pandas as pd 
%matplotlib inline
import matplotlib.pyplot as plt 
import seaborn as sns
from scipy import stats
from scipy.stats import norm, skew
from sklearn import preprocessing
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.linear_model import ElasticNetCV, ElasticNet
from xgboost import XGBRegressor, plot_importance 
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import StratifiedKFold
pd.set_option('display.float_format', lambda x: '{:.3f}'.format(x))
df = pd.read_csv('house_train.csv')
df.shape

(df.isnull().sum() / len(df)).sort_values(ascending=False)[:20]

好訊息是我們有很多特徵可以使用(81)，壞訊息是有19個特徵有缺失值，其中4個特徵缺失值超過80%。對於任何一個特徵，如果它缺失了80%的值，那麼它就沒有那麼重要了，因此，我決定刪除這4個特徵。

df.drop(['PoolQC', 'MiscFeature', 'Alley', 'Fence', 'Id'], axis=1, inplace=True)

探索特徵

標的特徵分佈

sns.distplot(df['SalePrice'] , fit=norm);

# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(df['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

# Now plot the distribution
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('Sale Price distribution')

# Get also the QQ-plot
fig = plt.figure()
res = stats.probplot(df['SalePrice'], plot=plt)
plt.show();

標的特徵-銷售價格是右傾斜。作為正態分佈資料的線性模型，我們將對銷售價格進行變換，使其更加正態分佈。

sns.distplot(np.log1p(df['SalePrice']) , fit=norm);

# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(np.log1p(df['SalePrice']))
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

# Now plot the distribution
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('log(Sale Price+1) distribution')

# Get also the QQ-plot
fig = plt.figure()
res = stats.probplot(np.log1p(df['SalePrice']), plot=plt)
plt.show();

數值特徵之間的相關性

pd.set_option('precision',2)
plt.figure(figsize=(10, 8))
sns.heatmap(df.drop(['SalePrice'],axis=1).corr(), square=True)
plt.suptitle("Pearson Correlation Heatmap")
plt.show();

這些特徵之間存在著很強的相關性。例如，GarageYrBlt和YearBuilt、TotRmsAbvGrd和GrLivArea、GarageArea和GarageCars是強相關的。它們實際上表達的是差不多一樣的東西。稍後我將讓ElasticNetCV幫助減少冗餘。

銷售價格與其他數字特徵之間的相關性

corr_with_sale_price = df.corr()["SalePrice"].sort_values(ascending=False)
plt.figure(figsize=(14,6))
corr_with_sale_price.drop("SalePrice").plot.bar()
plt.show();

銷售價格與總體質量的相關性最大(約0.8)。GrLivArea的相關係數超過0.7,GarageCars的相關係數超過0.6。讓我們更詳細地看看這4個特徵。

sns.pairplot(df[['SalePrice', 'OverallQual', 'GrLivArea', 'GarageCars']])
plt.show();

特徵工程

• 具有高度傾斜分佈(傾斜> 0.75)的對數變換特性
• 偽編碼分類特徵
• 用列的平均值填充NaN
• 訓練和測試集劃分

df["SalePrice"] = np.log1p(df["SalePrice"])

#log transform skewed numeric features:
numeric_feats = df.dtypes[df.dtypes != "object"].index

skewed_feats = df[numeric_feats].apply(lambda x: skew(x.dropna())) #compute skewness
skewed_feats = skewed_feats[skewed_feats > 0.75]
skewed_feats = skewed_feats.index

df[skewed_feats] = np.log1p(df[skewed_feats])
df = pd.get_dummies(df)
df = df.fillna(df.mean())

X, y = df.drop(['SalePrice'], axis = 1), df['SalePrice']
X_train, X_test, y_train, y_test  = train_test_split(X, y, test_size = 0.2, random_state = 0)

ElasticNet模型

• Ridge回歸和Lasso回歸是正則化線性回歸模型。
• ElasticNet本質上是一個Lasso/Ridge混合結構，它需要最小化一個包含L1 (Lasso)和L2(Ridge)規範的標的函式。
• 當有多個特徵具有相關性時，ElasticNet是有用的。
• 類ElasticNetCV透過交叉驗證可設定引數α(α)和l1_ratio(ρ)。
• ElasticNetCV: ElasticNet模型透過交叉驗證選擇最佳模型。

讓我們看看ElasticNetCV會為我們選擇什麼。

cv_model = ElasticNetCV(l1_ratio=[.1, .5, .7, .9, .95, .99, 1], eps=1e-3, n_alphas=100, fit_intercept=True, 
                        normalize=True, precompute='auto', max_iter=2000, tol=0.0001, cv=6, 
                        copy_X=True, verbose=0, n_jobs=-1, positive=False, random_state=0)

cv_model.fit(X_train, y_train)
print('Optimal alpha: %.8f'%cv_model.alpha_)
print('Optimal l1_ratio: %.3f'%cv_model.l1_ratio_)
print('Number of iterations %d'%cv_model.n_iter_)

0

模型評價

y_train_pred = cv_model.predict(X_train)
y_pred = cv_model.predict(X_test)
print('Train r2 score: ', r2_score(y_train_pred, y_train))
print('Test r2 score: ', r2_score(y_test, y_pred))
train_mse = mean_squared_error(y_train_pred, y_train)
test_mse = mean_squared_error(y_pred, y_test)
train_rmse = np.sqrt(train_mse)
test_rmse = np.sqrt(test_mse)
print('Train RMSE: %.4f' % train_rmse)
print('Test RMSE: %.4f' % test_rmse)

這裡的RMSE實際上是RMSLE(均方根對數誤差)。因為我們取了實際值的對數。這裡有一篇很好的文章來解釋RMSE和RMSLE之間的區別。

特徵的重要性

feature_importance = pd.Series(index = X_train.columns, data = np.abs(cv_model.coef_))

n_selected_features = (feature_importance>0).sum()
print('{0:d} features, reduction of {1:2.2f}%'.format(
    n_selected_features,(1-n_selected_features/len(feature_importance))*100))

feature_importance.sort_values().tail(30).plot(kind = 'bar', figsize = (12,5));

減少到58.91%的特徵看起來很有效。ElasticNetCV選擇的4個最重要的特性是Condition2_PosN、MSZoning_C(all)、Exterior1st_BrkComm & GrLivArea。我們將看到這些特徵如何與Xgboost所選擇的特徵進行比較。

xgboost

第一個Xgboost模型，我們從預設引數開始。

xgb_model1 = XGBRegressor()
xgb_model1.fit(X_train, y_train, verbose=False)
y_train_pred1 = xgb_model1.predict(X_train)
y_pred1 = xgb_model1.predict(X_test)

print('Train r2 score: ', r2_score(y_train_pred1, y_train))
print('Test r2 score: ', r2_score(y_test, y_pred1))
train_mse1 = mean_squared_error(y_train_pred1, y_train)
test_mse1 = mean_squared_error(y_pred1, y_test)
train_rmse1 = np.sqrt(train_mse1)
test_rmse1 = np.sqrt(test_mse1)
print('Train RMSE: %.4f' % train_rmse1)
print('Test RMSE: %.4f' % test_rmse1)

它已經比ElasticNetCV選擇的模型好得多!
在第二個Xgboost模型中，我們逐步添加了一些引數，這些引數假定可以增加模型的精度。

xgb_model2 = XGBRegressor(n_estimators=1000)
xgb_model2.fit(X_train, y_train, early_stopping_rounds=5, 
             eval_set=[(X_test, y_test)], verbose=False)
y_train_pred2 = xgb_model2.predict(X_train)
y_pred2 = xgb_model2.predict(X_test)

print('Train r2 score: ', r2_score(y_train_pred2, y_train))
print('Test r2 score: ', r2_score(y_test, y_pred2))
train_mse2 = mean_squared_error(y_train_pred2, y_train)
test_mse2 = mean_squared_error(y_pred2, y_test)
train_rmse2 = np.sqrt(train_mse2)
test_rmse2 = np.sqrt(test_mse2)
print('Train RMSE: %.4f' % train_rmse2)
print('Test RMSE: %.4f' % test_rmse2)

又有了進步!
第三個Xgboost模型，我們增加了一個學習率，希望它能產生一個更精確的模型。

xgb_model3 = XGBRegressor(n_estimators=1000, learning_rate=0.05)
xgb_model3.fit(X_train, y_train, early_stopping_rounds=5, 
             eval_set=[(X_test, y_test)], verbose=False)
y_train_pred3 = xgb_model3.predict(X_train)
y_pred3 = xgb_model3.predict(X_test)

print('Train r2 score: ', r2_score(y_train_pred3, y_train))
print('Test r2 score: ', r2_score(y_test, y_pred3))
train_mse3 = mean_squared_error(y_train_pred3, y_train)
test_mse3 = mean_squared_error(y_pred3, y_test)
train_rmse3 = np.sqrt(train_mse3)
test_rmse3 = np.sqrt(test_mse3)
print('Train RMSE: %.4f' % train_rmse3)
print('Test RMSE: %.4f' % test_rmse3)

遺憾的是，沒有任何改善。我的結論是xgb_model2是最好的模型。

特徵的重要性

from collections import OrderedDict
OrderedDict(sorted(xgb_model2.get_booster().get_fscore().items(), key=lambda t: t[1], reverse=True))

Xgboost選擇的最重要的4個特性是LotArea、GrLivArea、OverallQual和TotalBsmtSF。
只有一個特徵GrLivArea被ElasticNetCV和Xgboost選擇。
現在我們要選擇一些相關的特徵並再次擬合Xgboost。

most_relevant_features= list( dict((k, v) for k, v in xgb_model2.get_booster().get_fscore().items() if v >= 4).keys())
train_x=df[most_relevant_features]
train_y=df['SalePrice']
X_train, X_test, y_train, y_test  = train_test_split(train_x, train_y, test_size = 0.2, random_state = 0)
xgb_model5 = XGBRegressor(n_estimators=1000)
xgb_model5.fit(X_train, y_train, early_stopping_rounds=5, 
             eval_set=[(X_test, y_test)], verbose=False)
y_train_pred5 = xgb_model5.predict(X_train)
y_pred5 = xgb_model5.predict(X_test)

print('Train r2 score: ', r2_score(y_train_pred5, y_train))
print('Test r2 score: ', r2_score(y_test, y_pred5))
train_mse5 = mean_squared_error(y_train_pred5, y_train)
test_mse5 = mean_squared_error(y_pred5, y_test)
train_rmse5 = np.sqrt(train_mse5)
test_rmse5 = np.sqrt(test_mse5)
print('Train RMSE: %.4f' % train_rmse5)
print('Test RMSE: %.4f' % test_rmse5)

另一個小的進步!
Jupyter notebook可以在Github上找到。

您若是覺得有幫助，請點贊或者分享給朋友。
文中連結，請點選閱讀原文。您有任何問題，請留言。

原文連結：
https://www.kdnuggets.com/2019/05/modeling-price-regularized-linear-model-xgboost.html