从零开始实现一个端到端的机器学习项目[6]

基于书籍《Hands-on Machine Learning with Scikit-Learn, Keras & TensorFlow》的笔记

3 选择和训练模型

至此，我们以及提出了问题、获得了数据、进行了EDA、对数据集进行抽样分为了训练集和测试集，还编写了转换流水线能够实现自动化清洗和准备机器学习算法的数据，接下来，让我们一起来训练模型吧。

3.1 训练和评估训练集

事实上，数据处理步骤往往是最复杂的，“data is all you need！”，而现在训练模型的步骤在框架的帮助下比想象中容易很多。

线性模型是最简单的，让我们先尝试一下其效果。

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels)

ok，训练事实上只需要三行代码，现在让我们看看其效果如何：

>>> some_data = housing.iloc[:5]  #DataFrame取索引专用iloc
>>> some_labels = housing_labels.iloc[:5]
>>> some_data_prepared = full_pipeline.transform(some_data)
>>> print("Predictions:", lin_reg.predict(some_data_prepared))
Predictions: [ 210644.6045 317768.8069 210956.4333 59218.9888 189747.5584]
>>> print("Labels:", list(some_labels))
Labels: [286600.0, 340600.0, 196900.0, 46300.0, 254500.0]

可以看到，预测误差还是较大，特别是第一次预测误差接近40%。可以使用sklearn的mean_squared_error()函数来衡量整个训练集上回归模型的RMSE（root mean squre error）。

>>> from sklearn.metrics import mean_squared_error
>>> housing_predictions = lin_reg.predict(housing_prepared)
>>> lin_mse = mean_squared_error(housing_labels, housing_predictions)
>>> lin_rmse = np.sqrt(lin_mse)
>>> lin_rmse
68628.19819848922

大多数区域的median_housing_values分布在12000~265000美元之间，所以典型的预测误差达到68628美元确实是差强人意的。

这是一个典型的欠拟合的案例，对于这种情况，一个较大的可能是目前的特征不足以支撑我们去预测出更加准确的房价。

解决的办法是

1. 选择更强大的模型（more 参数）但这也需要更多数据，否则易过拟合。
2. 为算法提高更多更好的特征
3. 减少对模型的限制（regulation）-此例被排除，我们并没有正则化

我们来尝试一下训练一个更复杂的模型，决策树回归器。其善于从数据中找到复杂的非线性关系。

from sklearn.tree import DecisionTreeRegressor
tree_reg = DecisionTreeRegressor()
tree_reg.fit(housing_prepared, housing_labels)
>>> housing_predictions = tree_reg.predict(housing_prepared)
>>> tree_mse = mean_squared_error(housing_labels, housing_predictions)
>>> tree_rmse = np.sqrt(tree_mse)
>>> tree_rmse
0.0

ok，误差居然为0，这明显是不可能的，说明其发生了严重的过拟合。我们使用交叉验证集来check一下，测试集不到最后准备上线时都是不会用的。

3.2 使用交叉验证来更好地进行评估

使用前面提到的train_test_split函数可以将训练集拆分为训练集和交叉验证集（cross-validation），然后重新训练并在交叉验证集上进行验证。

另外一个常用的优秀选择是使用sklearn的K-折(K-fold)交叉验证功能。以下是执行K-折交叉验证的代码：它将训练集随机分割成10个不同的子集，每个子集成为一个fold，然后对决策树模型进行10次训练和评估-每次挑选1个fold进行评估，使用另外的9个fold进行训练。产生的结果是一个包含10次评估分数的数组：

from sklearn.model_selection import cross_val_score
scores = cross_val_score(tree_reg, housing_prepared, housing_labels,
 scoring="neg_mean_squared_error", cv=10)
tree_rmse_scores = np.sqrt(-scores)  #这里之所以会是负数，是因为交叉验证功能更倾向于使用正向的评估函数（越大越好）而不是成本函数（越小越好）所以上面用的scoring是“neg_mean_squared_error",下面我们得到rmse的时候需要先加入一个负号。

>>> def display_scores(scores):
... print("Scores:", scores)
... print("Mean:", scores.mean())
... print("Standard deviation:", scores.std())
...
>>> display_scores(tree_rmse_scores)
Scores: [70194.33680785 66855.16363941 72432.58244769 70758.73896782
 71115.88230639 75585.14172901 70262.86139133 70273.6325285
 75366.87952553 71231.65726027]
Mean: 71407.68766037929
Standard deviation: 2439.4345041191004
    
    
>>> lin_scores = cross_val_score(lin_reg, housing_prepared, housing_labels,
... scoring="neg_mean_squared_error", cv=10)
...
>>> lin_rmse_scores = np.sqrt(-lin_scores)
>>> display_scores(lin_rmse_scores)
Scores: [66782.73843989 66960.118071 70347.95244419 74739.57052552
 68031.13388938 71193.84183426 64969.63056405 68281.61137997
 71552.91566558 67665.10082067]
Mean: 69052.46136345083
Standard deviation: 2731.674001798348

重新审视整个结果，我们会发现决策树模型的表现甚至不如线性模型好，其确实出现了严重的过拟合k-fold的一个典型优势是我们还可以得出你预测性能的标准差来检查模型预测的robust，其代价也就是需要多次训练。-需要大量算力（特别是深度学习）。

最后，我们来尝试一个更加强大的模型-随机森林：其主要原理是：**通过对特征的随机子集进行多个决策树的训练，然后对其预测取平均。**在多个模型的基础上建立模型，称为集成学习（ensemble learning）,这是进一步提升机器学习系统性能的良策。

>>> from sklearn.ensemble import RandomForestRegressor
>>> forest_reg = RandomForestRegressor()
>>> forest_reg.fit(housing_prepared, housing_labels)
>>> [...] #和前面代码一样，省略了部分
>>> forest_rmse
18603.515021376355
>>> display_scores(forest_rmse_scores)
Scores: [49519.80364233 47461.9115823 50029.02762854 52325.28068953
 49308.39426421 53446.37892622 48634.8036574 47585.73832311
 53490.10699751 50021.5852922 ]
Mean: 50182.303100336096
Standard deviation: 2097.0810550985693

效果好多了，但是训练集上的分数仍然远低于验证集，这意味着该模型仍然对训练集过拟合。过拟合的可能解决方案包括：

1. 简化模型
2. 约束模型（加入正则化）
3. 获得更多训练数据

在一般的工程中，我们会再尝试几个模型（如不同内核的SVM和DNN），从而迅速筛选出2~5个有效的模型，再对具体的各个模型进行深入的探索，最后才是调整超参数。

我们应该对每一个尝试过的模型进行妥善的保存，以便随时可以进行回滚（rollback）。同时还需要保存超参数和训练过的参数，以及交叉验证的评分和实际预测的结果，可以方便最后对不同模型的评分或者不同模型造成的错误类型进行对比（无论是自用还是展示）。

import joblib
joblib.dump(my_model, "my_model.pkl")
# and later...
my_model_loaded = joblib.load("my_model.pkl")

3.3 微调（fine-tune) 模型

假设我们现在已经有了一个有效模型的候选列表，现在需要对它们进行微调。

3.3.1 网格搜索（grid search）

一种微调的方法是手动调整超参数，直到找到一组很好的超参数值组合。当然，这样太慢也不利于探索更多的组合，我们得想点更加聪明点的办法。

超参数的网格搜索是在实践中常用的方法，使用sklearn的GridSearchCV来替你实现探索。
其主要参数有：你需要试验的超参数是什么、以及需要尝试的值，它将会使用交叉验证来评估参数值的所有可能组合。

from sklearn.model_selection import GridSearchCV
param_grid = [
 {'n_estimators': [3, 10, 30], 'max_features': [2, 4, 6, 8]},
 {'bootstrap': [False], 'n_estimators': [3, 10], 'max_features': [2, 3, 4]},
 ]
forest_reg = RandomForestRegressor()
grid_search = GridSearchCV(forest_reg, param_grid, cv=5,
 scoring='neg_mean_squared_error',
return_train_score=True)
grid_search.fit(housing_prepared, housing_labels)

一般常用的超参数赋值是尝试10的连续幂次方，如果想要更细粒度的搜索，可以使用更小的数，例如这个示例中的n_estimators超参数。

上面的代码告诉sklearn，首先评估第一个dict的3*4=12种组合，然后评估第二个超参数组合，其加入了一个bootstrap为超参数，其默认值为True。

总而言之，网格搜索将探索RandomForestRegressor超参数值的18种组合，并对每个组合进行5折交叉验证，也就是最终会进行90次训练！这可能需要一些时间，但是完成后你就可以获得最佳的参数组合：

>>> grid_search.best_params_
{'max_features': 8, 'n_estimators': 30}

因为被评估的两个超参数的最大值和最好值都是8和30，所以你还可以试试更高的值，其评分或许还会继续被改善。

GridSearchCV如果被初始化为refit=True（其为默认值），那么一旦通过交叉验证找到了最佳估算器，它将在整个训练集上重新训练。这通常是个好办法，因为提供更多数据可能提升其性能。

gird_search为我们提供了两个属性（成员变量）可以访问最佳模型以及评估分数。

>>> grid_search.best_estimator_
RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
 max_features=8, max_leaf_nodes=None, min_impurity_decrease=0.0,
 min_impurity_split=None, min_samples_leaf=1,
 min_samples_split=2, min_weight_fraction_leaf=0.0,
 n_estimators=30, n_jobs=None, oob_score=False, random_state=None,
 verbose=0, warm_start=False)




>>> cvres = grid_search.cv_results_
>>> for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
... print(np.sqrt(-mean_score), params)
...
63669.05791727153 {'max_features': 2, 'n_estimators': 3}
55627.16171305252 {'max_features': 2, 'n_estimators': 10}
53384.57867637289 {'max_features': 2, 'n_estimators': 30}
60965.99185930139 {'max_features': 4, 'n_estimators': 3}
52740.98248528835 {'max_features': 4, 'n_estimators': 10}
50377.344409590376 {'max_features': 4, 'n_estimators': 30}
58663.84733372485 {'max_features': 6, 'n_estimators': 3}
52006.15355973719 {'max_features': 6, 'n_estimators': 10}
50146.465964159885 {'max_features': 6, 'n_estimators': 30}
57869.25504027614 {'max_features': 8, 'n_estimators': 3}
51711.09443660957 {'max_features': 8, 'n_estimators': 10}
49682.25345942335 {'max_features': 8, 'n_estimators': 30}
62895.088889905004 {'bootstrap': False, 'max_features': 2, 'n_estimators': 3}
54658.14484390074 {'bootstrap': False, 'max_features': 2, 'n_estimators': 10}
59470.399594730654 {'bootstrap': False, 'max_features': 3, 'n_estimators': 3}
52725.01091081235 {'bootstrap': False, 'max_features': 3, 'n_estimators': 10}
57490.612956065226 {'bootstrap': False, 'max_features': 4, 'n_estimators': 3}
51009.51445842374 {'bootstrap': False, 'max_features': 4, 'n_estimators': 10}

在本例中，我们得到的最佳解决方案是max_features=8、n_estimators=30. 其RMSE=49682,略低于之前使用默认超参数的分数50182.现在我们已经将模型调整到了最佳模式！

事实上，有一些数据准备的步骤也可以被当成超参数（例如是否使用转换器Combined-AttributesAdder的超参数add_bedrooms_per_room)使用网格搜索来进超参数搜索。同样，还可以用它来自动寻找处理问题的最佳方法，例如处理异常值、缺失特征，以及特征选择等。

3.3.2 随机搜索

如果探索的组合数量较少（例如上一个示例），那么网络搜索是一个不错的选择。但当超参数的搜索范围较大的时候，RandomizedSearchCV通常会是更好的选择，其使用起来与GridSearchCV类大致相同，但它不会尝试所有可能的组合，而是在每次迭代中为每个超参数选择一个随机值（在参数空间中随机采样，然后对给定数量的随机组合进行评估。这种方法的显著好处是：

1. 如果你让其跑1000次迭代，其就会对每一个超参数尝试1000个值而不是像网格搜索那样智能尝试几次。
2. 通过简单地设置迭代数量，你可以更好地控制在超参数搜索上所花费的计算时间。

https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.RandomizedSearchCV.html#sklearn.model_selection.RandomizedSearchCV-其官方API介绍：

我们选取两个重要的参数看看，其传入字典或者字典列表，以参数名称为key，值的分布或者列表作为要尝试的值。

如果给了分布，就进行随机采样，如果是list，就均匀采样，如果给了一个字典列表，则首先均匀采样字典列表中的参数。

n_iter参数控制总的采样数，其常常被用于权衡解决方案的质量和运行时间。

下面给一个示例代码：

param_dist = {'C': np.linspace(0.1, 10, 10),  #SVC的超参数
              'gamma': np.linspace(1, 0.01, 10)
              }
 
grid = RandomizedSearchCV(SVC(), param_distributions=param_dist, cv=4)

3.3.3 集成方法（模型集成-融合）

还有一种微调系统的方法是将最优的模型组合起来。组合（或“集成”）方法通常比最佳的单一模型更好（就像随机森林比其所依赖的任何单个决策树模型更好一样），特别是单一模型会产生不同类型误差时候，即如果不同的模型刚好学习到了数据的不同模式，那么综合起来将会是一个更加强大的系统。

3.3.4 分析最佳模型及其误差

通过检查最佳模型，我们可以获得一些好的洞见，例如在进行accurate 预测时，RandomForestRegressor 可以指出每个属性的相对重要程度：

>>> feature_importances = grid_search.best_estimator_.feature_importances_
>>> feature_importances
array([7.33442355e-02, 6.29090705e-02, 4.11437985e-02, 1.46726854e-02,
	   1.41064835e-02, 1.48742809e-02, 1.42575993e-02, 3.66158981e-01,
 	   5.64191792e-02, 1.08792957e-01, 5.33510773e-02, 1.03114883e-02,
 	   1.64780994e-01, 6.02803867e-05, 1.96041560e-03, 2.85647464e-03])

将其更好的可视化:

>>> extra_attribs = ["rooms_per_hhold", "pop_per_hhold", "bedrooms_per_room"]
>>> cat_encoder = full_pipeline.named_transformers_["cat"]
>>> cat_one_hot_attribs = list(cat_encoder.categories_[0])-#获取五个类别特征
>>> attributes = num_attribs + extra_attribs + cat_one_hot_attribs-#8+3（extra）+5=16个特征
>>> sorted(zip(feature_importances, attributes), reverse=True)
[(0.3661589806181342, 'median_income'),
 (0.1647809935615905, 'INLAND'),
 (0.10879295677551573, 'pop_per_hhold'),
 (0.07334423551601242, 'longitude'),
 (0.0629090704826203, 'latitude'),
 (0.05641917918195401, 'rooms_per_hhold'),
 (0.05335107734767581, 'bedrooms_per_room'),
 (0.041143798478729635, 'housing_median_age'),
 (0.014874280890402767, 'population'),
 (0.014672685420543237, 'total_rooms'),
 (0.014257599323407807, 'households'),
 (0.014106483453584102, 'total_bedrooms'),
 (0.010311488326303787, '<1H OCEAN'),
 (0.002856474637320158, 'NEAR OCEAN'),
 (0.00196041559947807, 'NEAR BAY'),
 (6.028038672736599e-05, 'ISLAND')]

通过对特征重要性的观察，我们可以尝试删除一些不太重要的特征(例如，本例中只有一个ocean-proximity（类别属性）是有用的，我们可以尝试删除其他所有特征)。

这时候应该去看一下系统的预测误差具体是什么，然后尝试去理解为什么会这样并尝试修复它。（例如添加额外特征、删除没有信息的特征、清除异常值等（我们前面提了但没做））。

3.3.5 通过测试集评估系统是否能够上线

在模型上线的前一步，我们需要最后使用测试集来评估模型的性能-值得注意的是，这里的full_pipeline只需要transform()即可，因为其已经在训练集上拟合过。

final_model = grid_search.best_estimator_

X_test = strat_test_set.drop("median_house_value", axis=1)
y_test = strat_test_set["median_house_value"].copy()

X_test_prepared = full_pipeline.transform(X_test)
final_predictions = final_model.predict(X_test_prepared)

final_mse = mean_squared_error(y_test, final_predictions)
final_rmse = np.sqrt(final_mse) # => evaluates to 47,730.2

我们可以通过scipy.stats.t.interval()来计算泛发误差的95%置信区间-（这将有利于我们决策其是否比当前已经部署的解决方案更好）:

>>> from scipy import stats
>>> confidence = 0.95
>>> squared_errors = (final_predictions - y_test) ** 2
>>> np.sqrt(stats.t.interval(confidence, len(squared_errors) - 1,
... loc=squared_errors.mean(),
... scale=stats.sem(squared_errors)))
...
array([45685.10470776, 49691.25001878])

如果之前进行过大量的超参数调整，这时的评估结果通常会略逊于你之前使用交叉验证时的表现结果（因为通过不断调整，系统在验证数据上终于表现良好，但是在未知数据集上通常可能达不到那么好的效果）。

OK，我们的模型从数据集洞察到最终模型的训练与观察就到此为止了，最后我们需要进行的是PPT展示，系统上线，然后对其进行监控和维护，这是一个冗长的过程，下面我们提几个重点吧。

3.3.6 上线、监控与维护机器学习系统

1.监控模型的实时性能，一般从下游任务推断模型的性能指标，例如当你的模型是推荐系统的一部分时，例如电商系统，如果电商系统的每日商品销售量突然下滑，就可能说明为用户的推荐不再那么精准了，那么主要的嫌疑就是模型，其可能原因是数据流水线被损坏了，你可能需要对新数据重新训练一个模型（实时性下降）。

2.人工评估是不可以被完全抹去的，AI距离超过人类还有很长一段路要走。

3.建立一个完整的故障解决方案是工程中必要的，其工作量往往比构建和训练模型更大。

4.当数据不断发现，则需要更新数据集并定期重新训练模型，整个流程应该自动化：

1. 定期收集模型并标注（如众包人工标注）
2. 编写脚本训练模型并自动微调超参数，该脚本能够根据你设定的时间定期运行
3. 编写另一个脚本，该脚本在更新的测试集上评估新模型和旧模型，如果性能良好，则将新模型部署到生产环境中（如果性能反而降低，则查找具体原因）

5.监控模型的输入与输出同样重要，因为模型的输入出现问题后反馈到输出往往需要一段时间。常见问题有：（越来越多的输入缺少了某种特征，其平均值或标准差离训练集太远（异常值），分类特征中出现了新类别）。

6.保存每个模型的备份，方便新模型出现故障时快速回滚，防止整个系统直接崩溃。

7.保留每个版本的数据集，作用和上面一样，且备份的数据集还可以用于评估新的模型（与新数据训练的模型比较）-如果发现新添入的数据都是离群值，方便快速回滚。

8.创建几个测试集的子集，用于评估模型在特定部分数据的效果如何（例如房价预测在某一个特定的沿海区域的预测误差如何）,

用于观察模型鲁棒性。

总结：

到这里，我已经带亲爱的读者您完整的过了一遍整个机器学习项目的流程，需要更加详细的流程可以查阅我发布的机器学习流程清单，此篇博客基于
《Hands-on Machine Learning with Scikit-Learn, Keras & TensorFlow》，在翻译的基础上加入更多自己的理解。欢迎您阅读原书，加深自己的理解。你还可以根据我的教程去kaggle动手试试看，看一下您学的怎样，如果您觉得学到很多，希望您能给一个一键三连。如果有机会，后面我将会在B站发布相关的讲解视频-带您手把手操作。
最后，我放一个覆盖完整的数据准备和最终预测的流水线的源码在此-所有代码都已经开源在-https://github.com/ageron/handson-ml2/blob/master/02_end_to_end_machine_learning_project.ipynb

Question: Try creating a single pipeline that does the full data preparation plus the final prediction.

prepare_select_and_predict_pipeline = Pipeline([
    ('preparation', full_pipeline),
    ('feature_selection', TopFeatureSelector(feature_importances, k)),
    ('svm_reg', SVR(**rnd_search.best_params_))
])

摘取一个原作者提出的有趣的问题和答案：

#Question: Try adding a transformer in the preparation pipeline to select only the most important attributes.(尝试在数据准备的流水线中加入一个转换器，从而只选出最重要的属性)

from sklearn.base import BaseEstimator, TransformerMixin

def indices_of_top_k(arr, k):
    return np.sort(np.argpartition(np.array(arr), -k)[-k:])

class TopFeatureSelector(BaseEstimator, TransformerMixin):
    def __init__(self, feature_importances, k):
        self.feature_importances = feature_importances
        self.k = k
    def fit(self, X, y=None):
        self.feature_indices_ = indices_of_top_k(self.feature_importances, self.k)
        return self
    def transform(self, X):
        return X[:, self.feature_indices_]
#Note: this feature selector assumes that you have already computed the feature importances somehow (for example using a RandomForestRegressor). You may be tempted to compute them directly in the TopFeatureSelector's fit() method, however this would likely slow down grid/randomized search since the feature importances would have to be computed for every hyperparameter combination (unless you implement some sort of cache).