我们想要建模房子的价格，我们知道价格取决于房子的位置，房子的面积，建成年限，翻新的年限，卧室的数量，车库的数量等等。因此，这些因素促成了这种模式——优质地段通常会导致更高的价格。然而，在同一区域内，面积相同的所有房子的价格并不完全相同。价格的变化就是噪声。我们对价格建模的目标是对模式进行建模，忽略噪声。同样的概念也适用于酒店房价的建模。
因此，首先，我们将对房价数据进行线性回归的正则化技术。
数据集
一个很好的房屋价格数据集可以在这里找到。
代码语言：javascript复制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个特征。
代码语言：javascript复制df.drop(['PoolQC', 'MiscFeature', 'Alley', 'Fence', 'Id'], axis=1, inplace=True)探索特征
目标特征分布
代码语言：javascript复制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();
目标特征-销售价格是右倾斜。作为正态分布数据的线性模型，我们将对销售价格进行变换，使其更加正态分布。
代码语言：javascript复制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();
数值特征之间的相关性
代码语言：javascript复制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帮助减少冗余。
销售价格与其他数字特征之间的相关性
代码语言：javascript复制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个特征。
代码语言：javascript复制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会为我们选择什么。
代码语言：javascript复制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<最优l1_ratio <1，说明惩罚是L1和L2的组合，即Lasso和Ridge的组合。
模型评价
代码语言：javascript复制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之间的区别。
特征的重要性
代码语言：javascript复制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模型，我们从默认参数开始。
代码语言：javascript复制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模型中，我们逐步添加了一些参数，这些参数假定可以增加模型的精度。
代码语言：javascript复制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模型，我们增加了一个学习率，希望它能产生一个更精确的模型。
代码语言：javascript复制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是最好的模型。
特征的重要性
代码语言：javascript复制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。
代码语言：javascript复制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)
