机器学习-正则化线性模型实例

原文：Regularized Linear Models

原文尝试了 Regularized Linear Regression（正则化线性回归）的使用，且作者认为此模型在较少的 feature engineering 的情况下表现出色，关键在于将偏离的数字变量通过 log 进行转换，因为大部分数字变量都是偏离的。

本流程结构非常简单：

1. 导入
2. 数据清洗及补全，包含纠偏（待解释）、变量转换、数据补充
3. 尝试 sklearn 的两个 model，根据误差选择表现更好的 model
4. 补充 xgboost 模型
5. 将 sklearn 中的 Lasso 与 xgboost 加权进行预测得到最终学习结果

其中：

1. 数据清洗及补全的所有步骤都可以直接使用
2. 筛选 Lasso 和 Ridge 模型时的 alpha 值选择可以借鉴，重点参考语法的使用
3. 将 Lasso 与 Xgboost 加权结合的参数选择可以借鉴

具体算法流程

1. 导入工具及数据准备

1.1 导入所有工具

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt
from scipy.stats import skew
from scipy.stats.stats import pearsonr
%config InlineBackend.figure_format = 'retina' # set 'png' here when working on notebook
%matplotlib inline #用法搞准

1.2 导入数据

数据下载地址 https://www.kaggle.com/c/house-prices-advanced-regression-techniques/data

train = pd.read_csv("../input/train.csv")
test = pd.read_csv("../input/test.csv")

2. 数据处理

数据处理的整体思路较为简单，关键三点：

• 把偏离的数字型 feature 通过 log(feature+1) 进行转换，可以使 feature 更正常（朴素贝叶斯用到过）
• 为分类型 feature 创建假/模拟变量 （用 pandas 的 get_dummies 函数）
• 将数字值中的空值替换为该列平均数

2.1 拼接 train 和 test 数据

train.head()
all_data = pd.concat((train.loc[:,'MSSubClass':'SaleCondition'],
                    test.loc[:,'MSSubClass':'SaleCondition']))

截取两批数据的 MSSubClass 到 SaleCondition 的列合并为 all_data 所有数据。

2.2 数据的 log 转换（纠偏）

2.2.1 以 price 值的 log 转换示意：

matplotlib.rcParams['figure.figsize'] = (12.0, 6.0)
prices = pd.DataFrame({"price":train["SalePrice"], "log(price + 1)":np.log1p(train["SalePrice"])})
prices.hist()#直方图

outcome

左：log(feature+1) 转换后，右：转换前

2.2.2 整体数据的 log 转换

#log transform the target:
train["SalePrice"] = np.log1p(train["SalePrice"])

#log transform skewed numeric features:
numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index#object 指 categorical features 的数据类型

skewed_feats = train[numeric_feats].apply(lambda x: skew(x.dropna())) #compute skewness 计算偏离度
skewed_feats = skewed_feats[skewed_feats > 0.75] #筛选高偏离数据
skewed_feats = skewed_feats.index

all_data[skewed_feats] = np.log1p(all_data[skewed_feats])

2.3 为分类型 feature 创建假/模拟变量

all_data = pd.get_dummies(all_data)

2.4 将数字值中的空值替换为该列平均数

#filling NA's with the mean of the column:
all_data = all_data.fillna(all_data.mean())

2.5 为模型创建矩阵

#creating matrices for sklearn:
X_train = all_data[:train.shape[0]]#shape 语法:读取矩阵的长度，比如 shape[0] 就是读取矩阵第一维度的长度
X_test = all_data[train.shape[0]:]
y = train.SalePrice

3. 套用模型

现在准备使用 regularized linear regression （正则化线性回归）模型。打算尝试两种 regularization（正则）方式：l_1(Lasso) and l_2(Ridge). 并且还会写一个函数，返回 cross-validation rmse error（交叉验证均方根误差）以便评估模型，选择最好的 tuning par（调节参数）。par for parameter

3.1 导入模块

导入 sklearn 模块里的 Lasso 和 Ridge，并创建返回 cross-validation rmse error 的函数 rmse_cv

from sklearn.linear_model import Ridge, RidgeCV, ElasticNet, LassoCV, LassoLarsCV
from sklearn.model_selection import cross_val_score

def rmse_cv(model):
    rmse= np.sqrt(-cross_val_score(model, X_train, y, scoring="neg_mean_squared_error", cv = 5))
    return(rmse)

3.2 选用 Ridge 并选取主参数

对于 Ridge 模型来说，主要调节参数是 alpha —— 一个衡量模型灵活度的正则参数。正则度越高，越不可能 overfit。但是它也会导致模型灵活度的降低，可能无法捕捉数据中的所有信号。

关于 alpha 的选择，示意如下：

model_ridge = Ridge()
alphas = [0.05, 0.1, 0.3, 1, 3, 5, 10, 15, 30, 50, 75]#通常值
cv_ridge = [rmse_cv(Ridge(alpha = alpha)).mean() 
            for alpha in alphas]
cv_ridge = pd.Series(cv_ridge, index = alphas)
cv_ridge.plot(title = "Validation - Just Do It")
plt.xlabel("alpha")
plt.ylabel("rmse")

nike

上图可见对勾形状的曲线（耐克），横轴为 alpha 值，纵轴为误差。当 alpha 过大时，正则过强，模型就无法捕捉到数据中的复杂性（即误差较高）。但若 alpha 过小（模型灵活度太高），则模型开始 overfit（误差也高）。根据上图，alpha = 10 是比较合适的一个参数选择（误差最低）。

cv_ridge.min()
0.12733734668670765

得出 Ridge 回归下的误差大约是 0.127。

3.3 选用 Lasso

然后再尝试一下 Lasso 模型。这里的方法与之前稍有不同，直接使用 Lasso CV 内置设置挑选最好的 alpha。==出于某种原因 Lasso CV 里的 alpha 是 Ridge 里的反数==。

model_lasso = LassoCV(alphas = [1, 0.1, 0.001, 0.0005]).fit(X_train, y) #试出来的 #fit 把数据套进模型里跑
rmse_cv(model_lasso).mean()
0.12314421090977427

3.5 使用 Lasso 选择 feature

可以看出 Lasso 表现更好（平均误差更低：0.123<0.127），因此就选择这个模型做预测了。Lasso 还有个好处是它会帮你选 feature——把它认为不重要的 feature 的系数设为 0. 下面来看一下参数：

coef = pd.Series(model_lasso.coef_, index = X_train.columns)
print("Lasso picked " + str(sum(coef != 0)) + " variables and eliminated the other " +  str(sum(coef == 0)) + " variables")

Lasso 选了 110 个变量，删了其他的 178 个。

Lasso 做得很好。但这里有一点要注意：被选出来的 features 不一定是“正确”的 features——尤其是在这个数据集里有很多 features 是同线的。有一种做法是在 boostrapped samples 上多跑几次 Lasso，看 feature 选择的稳定程度。也可以直接看一下最重要的系数都是哪些：（所以前面挑出来的 feature 只是看看？）

imp_coef = pd.concat([coef.sort_values().head(10),
                     coef.sort_values().tail(10)])#选头尾各10条
matplotlib.rcParams['figure.figsize'] = (8.0, 10.0)
imp_coef.plot(kind = "barh")
plt.title("Coefficients in the Lasso Model")   

重要系数

最重要的正向 feature 是 GrLivArea——地上部分的面积。这个肯定说得通。还有其他的地点和质量 feature 也是正向的。一部分负向的 feature 不太说得通，需要进一步观察——它们可能是由不平衡的分类型变量引起的。

还需要注意的是，不像 feature importance 可以直接从 random forest 里得到，这些参数是模型里==真正的==参数——所以你可以很明确地知道为什么预测出的会是这个价格。唯一的问题就是我们把目标和数字 feature 都 log 转换了，所以实际的数量级有点难判断。

3.6 检查 residuals（预测误差，预测值-实际值）

#let's look at the residuals as well:
matplotlib.rcParams['figure.figsize'] = (6.0, 6.0)#图片大小 6*6
preds = pd.DataFrame({"preds":model_lasso.predict(X_train), "true":y}) # predict：预测结果，y：实际结果，二者之差为 residual
preds["residuals"] = preds["true"] - preds["preds"]
preds.plot(x = "preds", y = "residuals",kind = "scatter")

residual plot

预测误差的图看着也不错，大部分基本集中在 0 附近。最后预测一下 test 数据集，提交到 leaderboard。

3.7 添加一个 xgboost 模型

在目前的线性模型上再加一个 xgboost 模型，看能否进一步提高分数（这里的分数应该是 kaggle 上对于机器学习模型的一个得分？）。

import xgboost as xgb
dtrain = xgb.DMatrix(X_train, label = y)#DMatrix 是 xgb 存储信息的单位，本步骤把数据放进这里面去
dtest = xgb.DMatrix(X_test)
params = {"max_depth":2, "eta":0.1}#max_depth：最大深度。eta 和 gradiant boosting 中的 learning rate 参数类似。通过减少每一步的权重，可以提高模型的稳定性。 典型值为 0.01-0.2。
model = xgb.cv(params, dtrain,  num_boost_round=500, early_stopping_rounds=100)# CV 用法及参数见下文
model.loc[30:,["test-rmse-mean", "train-rmse-mean"]].plot()

Early_stopping_rounds: 提前终止程序

如果有评价数据，可以提前终止程序，这样可以找到最优的迭代次数。如果要提前终止程序必须至少有一个评价数据在参数evals中。 超过一个则使用最后一个。

train(..., evals=evals, early_stopping_rounds=10)

此模型下，机器会一直学习到 validation score 不再增长。每经过 early_stopping_rounds 轮，误差应该都有所下降，否则不应该继续学习。

如果出现了提前终止，模型会出现两个额外情况：bst.best_score 和 bst.best_iteration. 注意此处 train() 会返回最后一次循环的模型，而非最好的循环的模型。

此方法适用于各类最低 (RMSE, log loss, etc.) 或最高 (MAP, NDCG, AUC) 误差计算。

观察误差情况

横轴 boost round，最大 500；纵轴是平均误差。

model_xgb = xgb.XGBRegressor(n_estimators=360, max_depth=2, learning_rate=0.1) #用了xgb.cv 调参。n_estimators：训练的轮数；max_depth：最大深度。
model_xgb.fit(X_train, y)

#结果：
XGBRegressor(base_score=0.5, colsample_bylevel=1, colsample_bytree=1, gamma=0,
       learning_rate=0.1, max_delta_step=0, max_depth=2,
       min_child_weight=1, missing=None, n_estimators=360, nthread=-1,
       objective='reg:linear', reg_alpha=0, reg_lambda=1,
       scale_pos_weight=1, seed=0, silent=True, subsample=1)
# 所有参数详解见后文

xgb_preds = np.expm1(model_xgb.predict(X_test))#expm1：对数列所有对象执行 exp(x) - 1 计算，这里是对 test 预测值进行调整
lasso_preds = np.expm1(model_lasso.predict(X_test))#lasso也一样
predictions = pd.DataFrame({"xgb":xgb_preds, "lasso":lasso_preds})
predictions.plot(x = "xgb", y = "lasso", kind = "scatter")

横轴：xgb 预测值，纵轴：lasso 预测值。散点图可以看出强线性相关，两个预测结果大部分基本一致。

很多情况下把不相关的结果进行加权平均是有用的，通常能够优化结果，虽然在这个案例里帮助不大。

preds = 0.7*lasso_preds + 0.3*xgb_preds#把 lasso 和 xgb 加权后得到最终预测值，为啥是 0.7 和 0.3？经验吧
solution = pd.DataFrame({"id":test.Id, "SalePrice":preds})
solution.to_csv("ridge_sol.csv", index = False)

XGBoost 的参数

1. eta [默认 0.3]

和 GBM 中的 learning rate 参数类似。 通过减少每一步的权重，可以提高模型的稳定性。 典型值为 0.01-0.2。

2. min_child_weight [默认 1]

决定最小叶子节点样本权重和。和 GBM 的 min_child_leaf 参数类似，但不完全一样。XGBoost 的这个参数是最小样本权重的和，而 GBM 参数是最小样本总数。这个参数用于避免过拟合。当它的值较大时，可以避免模型学习到局部的特殊样本。但是如果这个值过高，会导致欠拟合。这个参数需要使用 CV 来调整。

3. max_depth [默认 6]

和 GBM 中的参数相同，这个值为树的最大深度。这个值也是用来避免过拟合的。max_depth 越大，模型会学到更具体更局部的样本。需要使用 CV 函数来进行调优。 典型值：3-10

4. max_leaf_nodes

树上最大的节点或叶子的数量。 可以替代 max_depth 的作用。因为如果生成的是二叉树，一个深度为 n 的树最多生成 n2 个叶子。 如果定义了这个参数，GBM 会忽略 max_depth 参数。

5. gamma [默认 0]

在节点分裂时，只有分裂后损失函数的值下降了，才会分裂这个节点。Gamma 指定了节点分裂所需的最小损失函数下降值。 这个参数的值越大，算法越保守。这个参数的值和损失函数息息相关，所以是需要调整的。

6、max_delta_step[默认 0]

这参数限制每棵树权重改变的最大步长。如果这个参数的值为 0，那就意味着没有约束。如果它被赋予了某个正值，那么它会让这个算法更加保守。 通常，这个参数不需要设置。但是当各类别的样本十分不平衡时，它对逻辑回归是很有帮助的。 这个参数一般用不到，但是你可以挖掘出来它更多的用处。

7. subsample [默认 1]

和 GBM 中的 subsample 参数一模一样。这个参数控制对于每棵树，随机采样的比例。 减小这个参数的值，算法会更加保守，避免过拟合。但是，如果这个值设置得过小，它可能会导致欠拟合。 典型值：0.5-1

8. colsample_bytree [默认 1]

和 GBM 里面的 max_features 参数类似。用来控制每棵随机采样的列数的占比 (每一列是一个特征)。 典型值：0.5-1

9. colsample_bylevel [默认 1]

用来控制树的每一级的每一次分裂，对列数的采样的占比。 我个人一般不太用这个参数，因为 subsample 参数和 colsample_bytree 参数可以起到相同的作用。但是如果感兴趣，可以挖掘这个参数更多的用处。

10. lambda [默认 1]

权重的 L2 正则化项。(和 Ridge regression 类似)。 这个参数是用来控制 XGBoost 的正则化部分的。虽然大部分数据科学家很少用到这个参数，但是这个参数在减少过拟合上还是可以挖掘出更多用处的。

11. alpha [默认 1]

权重的 L1 正则化项。(和 Lasso regression 类似)。 可以应用在很高维度的情况下，使得算法的速度更快。

12. scale_pos_weight [默认 1]

在各类别样本十分不平衡时，把这个参数设定为一个正值，可以使算法更快收敛。

学习目标参数

这个参数用来控制理想的优化目标和每一步结果的度量方法。

1. objective [默认 reg:linear]

这个参数定义需要被最小化的损失函数。最常用的值有：

binary:logistic 二分类的逻辑回归，返回预测的概率 (不是类别)。 multi:softmax 使用 softmax 的多分类器，返回预测的类别 (不是概率)。

在这种情况下，你还需要多设一个参数：num_class(类别数目)。 multi:softprob 和 multi:softmax 参数一样，但是返回的是每个数据属于各个类别的概率。

2. eval_metric [默认值取决于 objective 参数的取值]

对于有效数据的度量方法。对于回归问题，默认值是 rmse，对于分类问题，默认值是 error。 典型值有：

rmse 均方根误差、mae 平均绝对误差、logloss 负对数似然函数值、error 二分类错误率 (阈值为 0.5)、merror 多分类错误率、mlogloss 多分类 logloss 损失函数、auc 曲线下面积

3. seed [默认 0]

随机数的种子设置它可以复现随机数据的结果，也可以用于调整参数。