原文:
annas-archive.org/md5/7039549ff2b32d189f96a3420dc66360
译者:飞龙
第七章:交叉验证和后模型工作流
在本章中,我们将涵盖以下内容:
-
使用交叉验证选择模型
-
K 折交叉验证
-
平衡交叉验证
-
使用 ShuffleSplit 的交叉验证
-
时间序列交叉验证
-
使用 scikit-learn 进行网格搜索
-
使用 scikit-learn 进行随机搜索
-
分类度量
-
回归度量
-
聚类度量
-
使用虚拟估计器比较结果
-
特征选择
-
L1 范数上的特征选择
-
使用 joblib 或 pickle 持久化模型
介绍
这也许是最重要的章节。本章所探讨的基本问题如下:
- 我们如何选择一个预测良好的模型?
这就是交叉验证的目的,不管模型是什么。这与传统统计略有不同,传统统计更关心我们如何更好地理解现象。(为什么要限制我对理解的追求?好吧,因为有越来越多的数据,我们不一定能够全部查看、反思并创建理论模型。)
机器学习关注预测以及机器学习算法如何处理新的未见数据并得出预测。即使它看起来不像传统统计,你可以使用解释和领域理解来创建新列(特征)并做出更好的预测。你可以使用传统统计来创建新列。
书的早期,我们从训练/测试拆分开始。交叉验证是许多关键训练和测试拆分的迭代,以最大化预测性能。
本章探讨以下内容:
-
交叉验证方案
-
网格搜索——在估计器中找到最佳参数是什么?
-
评估指标比较
y_test
与y_pred
——真实目标集与预测目标集
下面一行包含交叉验证方案 cv = 10
,用于 neg_log_lost
评分机制,该机制由 log_loss
指标构建:
cross_val_score(SVC(), X_train, y_train, cv = 10, scoring='neg_log_loss')
Scikit-learn 的一部分力量在于在一行代码中包含如此多的信息。此外,我们还将看到一个虚拟估计器,查看特征选择,并保存训练好的模型。这些方法真正使得机器学习成为它所是的东西。
使用交叉验证选择模型
我们看到了自动交叉验证,在 第一章,高性能机器学习 – NumPy 中的 cross_val_score
函数。这将非常相似,除了我们将使用鸢尾花数据集的最后两列作为数据。本节的目的是选择我们可以选择的最佳模型。
在开始之前,我们将定义最佳模型为得分最高的模型。如果出现并列,我们将选择波动最小的模型。
准备就绪
在这个配方中,我们将执行以下操作:
-
加载鸢尾花数据集的最后两个特征(列)
-
将数据拆分为训练数据和测试数据
-
实例化两个k 近邻(KNN)算法,分别设置为三个和五个邻居。
-
对两个算法进行评分
-
选择得分最好的模型
从加载数据集开始:
from sklearn import datasets
iris = datasets.load_iris()
X = iris.data[:,2:]
y = iris.target
将数据拆分为训练集和测试集。样本是分层抽样的,书中默认使用这种方法。分层抽样意味着目标变量在训练集和测试集中的比例相同(此外,random_state
被设置为7
):
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, stratify = y,random_state = 7)
如何操作……
- 首先,实例化两个最近邻算法:
from sklearn.neighbors import KNeighborsClassifier
kn_3 = KNeighborsClassifier(n_neighbors = 3)
kn_5 = KNeighborsClassifier(n_neighbors = 5)
- 现在,使用
cross_val_score
对两个算法进行评分。查看kn_3_scores
,这是得分的列表:
from sklearn.model_selection import cross_val_score
kn_3_scores = cross_val_score(kn_3, X_train, y_train, cv=4)
kn_5_scores = cross_val_score(kn_5, X_train, y_train, cv=4)
kn_3_scores
array([ 0.9 , 0.92857143, 0.92592593, 1\. ])
- 查看
kn_5_scores
,这是另一个得分列表:
kn_5_scores
array([ 0.96666667, 0.96428571, 0.88888889, 1\. ])
- 查看两个列表的基本统计信息。查看均值:
print "Mean of kn_3: ",kn_3_scores.mean()
print "Mean of kn_5: ",kn_5_scores.mean()
Mean of kn_3: 0.938624338624
Mean of kn_5: 0.95496031746
- 查看分布,查看标准差:
print "Std of kn_3: ",kn_3_scores.std()
print "Std of kn_5: ",kn_5_scores.std()
Std of kn_3: 0.037152126551
Std of kn_5: 0.0406755710299
总体来说,当算法设置为五个邻居时,kn_5
的表现比三个邻居稍好,但它的稳定性较差(它的得分有点分散)。
- 现在我们进行最后一步:选择得分最高的模型。我们选择
kn_5
,因为它的得分最高。(该模型在交叉验证下得分最高。请注意,涉及的得分是最近邻的默认准确率得分:正确分类的比例除以所有尝试分类的数量。)
它是如何工作的……
这是一个 4 折交叉验证的示例,因为在cross_val_score
函数中,cv = 4
。我们将训练数据,或交叉验证集(X_train
),拆分为四个部分,或称折叠。我们通过轮流将每个折叠作为测试集来迭代。首先,折叠 1 是测试集,而折叠 2、3 和 4 一起构成训练集。接下来,折叠 2 是测试集,而折叠 1、3 和 4 是训练集。我们还对折叠 3 和折叠 4 进行类似的操作:
一旦我们将数据集拆分为折叠,我们就对算法进行四次评分:
-
我们在折叠 2、3 和 4 上训练其中一个最近邻算法。
-
然后我们对折叠 1 进行预测,即测试折。
-
我们衡量分类准确性:将测试折与该折的预测结果进行比较。这是列表中第一个分类分数。
该过程执行了四次。最终输出是一个包含四个分数的列表。
总体来说,我们进行了整个过程两次,一次用于kn_3
,一次用于kn_5
,并生成了两个列表以选择最佳模型。我们从中导入的模块叫做model_selection
,因为它帮助我们选择最佳模型。
K 折交叉验证
在寻找最佳模型的过程中,你可以查看交叉验证折叠的索引,看看每个折叠中有哪些数据。
准备工作
创建一个非常小的玩具数据集:
import numpy as np
X = np.array([[1, 2], [3, 4], [5, 6], [7, 8],[1, 2], [3, 4], [5, 6], [7, 8]])
y = np.array([1, 2, 1, 2, 1, 2, 1, 2])
如何操作……
- 导入
KFold
并选择拆分的数量:
from sklearn.model_selection import KFold
kf= KFold(n_splits = 4)
- 你可以遍历生成器并打印出索引:
cc = 1
for train_index, test_index in kf.split(X):
print "Round : ",cc,": ",
print "Training indices :", train_index,
print "Testing indices :", test_index
cc += 1
Round 1 : Training indices : [2 3 4 5 6 7] Testing indices : [0 1]
Round 2 : Training indices : [0 1 4 5 6 7] Testing indices : [2 3]
Round 3 : Training indices : [0 1 2 3 6 7] Testing indices : [4 5]
Round 4 : Training indices : [0 1 2 3 4 5] Testing indices : [6 7]
你可以看到,例如,在第一轮中有两个测试索引,0
和1
。[0 1]
构成了第一个折叠。[2 3 4 5 6 7]
是折叠 2、3 和 4 的组合。
- 你还可以查看拆分的次数:
kf.get_n_splits()
4
分割数为 4
,这是我们实例化 KFold
类时设置的。
还有更多…
如果愿意,可以查看折叠数据本身。将生成器存储为列表:
indices_list = list(kf.split(X))
现在,indices_list
是一个元组的列表。查看第四个折叠的信息:
indices_list[3] #the list is indexed from 0 to 3
(array([0, 1, 2, 3, 4, 5], dtype=int64), array([6, 7], dtype=int64))
此信息与前面的打印输出信息相匹配,但它以两个 NumPy 数组的元组形式给出。查看第四个折叠的实际数据。查看第四个折叠的训练数据:
train_indices, test_indices = indices_list[3]
X[train_indices]
array([[1, 2],
[3, 4],
[5, 6],
[7, 8],
[1, 2],
[3, 4]])
y[train_indices]
array([1, 2, 1, 2, 1, 2])
查看测试数据:
X[test_indices]
array([[5, 6],
[7, 8]])
y[test_indices]
array([1, 2])
平衡交叉验证
在将不同折叠中的不同数据集分割时,您可能会想知道:k 折交叉验证中每个折叠中的不同集合可能会非常不同吗?每个折叠中的分布可能会非常不同,这些差异可能导致得分的波动。
对此有一个解决方案,使用分层交叉验证。数据集的子集看起来像整个数据集的较小版本(至少在目标变量方面)。
准备工作
创建一个玩具数据集如下:
import numpy as np
X = np.array([[1, 2], [3, 4], [5, 6], [7, 8],[1, 2], [3, 4], [5, 6], [7, 8]])
y = np.array([1, 1, 1, 1, 2, 2, 2, 2])
如何做…
- 如果我们在这个小型玩具数据集上执行 4 折交叉验证,每个测试折叠将只有一个目标值。可以使用
StratifiedKFold
来解决这个问题:
from sklearn.model_selection import StratifiedKFold
skf = StratifiedKFold(n_splits = 4)
- 打印出折叠的索引:
cc = 1
for train_index, test_index in skf.split(X,y):
print "Round",cc,":",
print "Training indices :", train_index,
print "Testing indices :", test_index
cc += 1
Round 1 : Training indices : [1 2 3 5 6 7] Testing indices : [0 4]
Round 2 : Training indices : [0 2 3 4 6 7] Testing indices : [1 5]
Round 3 : Training indices : [0 1 3 4 5 7] Testing indices : [2 6]
Round 4 : Training indices : [0 1 2 4 5 6] Testing indices : [3 7]
观察 skf
类的 split
方法,即分层 k 折叠分割,具有两个参数 X
和 y
。它试图在每个折叠集中以相同的分布分配目标 y
。在这种情况下,每个子集都有 50% 的 1
和 50% 的 2
,就像整个目标集 y
一样。
还有更多…
您可以使用 StratifiedShuffleSplit
重新洗牌分层折叠。请注意,这并不会尝试制作具有互斥测试集的四个折叠:
from sklearn.model_selection import StratifiedShuffleSplit
sss = StratifiedShuffleSplit(n_splits = 5,test_size=0.25)
cc = 1
for train_index, test_index in sss.split(X,y):
print "Round",cc,":",
print "Training indices :", train_index,
print "Testing indices :", test_index
cc += 1
Round 1 : Training indices : [1 6 5 7 0 2] Testing indices : [4 3]
Round 2 : Training indices : [3 2 6 7 5 0] Testing indices : [1 4]
Round 3 : Training indices : [2 1 4 7 0 6] Testing indices : [3 5]
Round 4 : Training indices : [4 2 7 6 0 1] Testing indices : [5 3]
Round 5 : Training indices : [1 2 0 5 4 7] Testing indices : [6 3]
Round 6 : Training indices : [0 6 5 1 7 3] Testing indices : [2 4]
Round 7 : Training indices : [1 7 3 6 2 5] Testing indices : [0 4]
这些分割不是数据集的分割,而是随机过程的迭代,每个迭代的训练集大小为整个数据集的 75%,测试集大小为 25%。所有迭代都是分层的。
使用 ShuffleSplit 的交叉验证
ShuffleSplit 是最简单的交叉验证技术之一。使用这种交叉验证技术只需取数据的样本,指定的迭代次数。
准备工作
ShuffleSplit 是一种简单的验证技术。我们将指定数据集中的总元素数量,其余由它来处理。我们将通过估计单变量数据集的均值来示例化。这类似于重新采样,但它将说明为什么我们要在展示交叉验证时使用交叉验证。
如何做…
- 首先,我们需要创建数据集。我们将使用 NumPy 创建一个数据集,其中我们知道底层均值。我们将对数据集的一半进行采样以估计均值,并查看它与底层均值的接近程度。生成一个均值为 1000,标准差为 10 的正态分布随机样本:
%matplotlib inline
import numpy as np
true_mean = 1000
true_std = 10
N = 1000
dataset = np.random.normal(loc= true_mean, scale = true_std, size=N)
import matplotlib.pyplot as plt
f, ax = plt.subplots(figsize=(10, 7))
ax.hist(dataset, color='k', alpha=.65, histtype='stepfilled',bins=50)
ax.set_title("Histogram of dataset")
- 估计数据集的一半的平均值:
holdout_set = dataset[:500]
fitting_set = dataset[500:]
estimate = fitting_set[:N/2].mean()
estimate
999.69789261486721
- 你也可以获取整个数据集的均值:
data_mean = dataset.mean()
data_mean
999.55177343767843
- 它不是 1,000,因为随机选择了点来创建数据集。为了观察
ShuffleSplit
的行为,写出以下代码并绘制图表:
from sklearn.model_selection import ShuffleSplit
shuffle_split = ShuffleSplit(n_splits=100, test_size=.5, random_state=0)
mean_p = []
estimate_closeness = []
for train_index, not_used_index in shuffle_split.split(fitting_set):
mean_p.append(fitting_set[train_index].mean())
shuf_estimate = np.mean(mean_p)
estimate_closeness.append(np.abs(shuf_estimate - dataset.mean()))
plt.figure(figsize=(10,5))
plt.plot(estimate_closeness)
估计的均值不断接近数据的均值 999.55177343767843,并在距离数据均值 0.1 时停滞。它比仅用一半数据估算出的均值更接近数据的均值。
时间序列交叉验证
scikit-learn 可以对时间序列数据(例如股市数据)进行交叉验证。我们将使用时间序列拆分,因为我们希望模型能够预测未来,而不是从未来泄漏信息。
准备工作
我们将为时间序列拆分创建索引。首先创建一个小的玩具数据集:
from sklearn.model_selection import TimeSeriesSplit
import numpy as np
X = np.array([[1, 2], [3, 4], [1, 2], [3, 4],[1, 2], [3, 4], [1, 2], [3, 4]])
y = np.array([1, 2, 3, 4, 1, 2, 3, 4])
如何做…
- 现在创建一个时间序列拆分对象:
tscv = TimeSeriesSplit(n_splits=7)
- 遍历它:
for train_index, test_index in tscv.split(X):
X_train, X_test = X[train_index], X[test_index]
y_train, y_test = y[train_index], y[test_index]
print "Training indices:", train_index, "Testing indices:", test_index
Training indices: [0] Testing indices: [1]
Training indices: [0 1] Testing indices: [2]
Training indices: [0 1 2] Testing indices: [3]
Training indices: [0 1 2 3] Testing indices: [4]
Training indices: [0 1 2 3 4] Testing indices: [5]
Training indices: [0 1 2 3 4 5] Testing indices: [6]
Training indices: [0 1 2 3 4 5 6] Testing indices: [7]
- 你也可以通过创建一个包含元组的列表来保存索引:
tscv_list = list(tscv.split(X))
还有更多…
你也可以使用 NumPy 或 pandas 创建滚动窗口。时间序列交叉验证的主要要求是测试集必须出现在训练集之后;否则,你就会从未来预测过去。
时间序列交叉验证很有趣,因为根据数据集的不同,时间的影响是不同的。有时,你不需要将数据行按时间顺序排列,但你永远不能假设你知道过去的未来。
使用 scikit-learn 进行网格搜索
在模型选择和交叉验证章节的开头,我们尝试为鸢尾花数据集的最后两个特征选择最佳的最近邻模型。现在,我们将使用 GridSearchCV
在 scikit-learn 中重新聚焦这一点。
准备工作
首先,加载鸢尾花数据集的最后两个特征。将数据拆分为训练集和测试集:
from sklearn import datasets
iris = datasets.load_iris()
X = iris.data[:,2:]
y = iris.target
from sklearn.model_selection import train_test_split, cross_val_score
X_train, X_test, y_train, y_test = train_test_split(X, y, stratify = y,random_state = 7)
如何做…
- 实例化一个最近邻分类器:
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
- 准备一个参数网格,这是网格搜索所必需的。参数网格是一个字典,包含你希望尝试的参数设置:
param_grid = {'n_neighbors': list(range(3,9,1))}
-
实例化一个网格搜索,传入以下参数:
-
估计器
-
参数网格
-
一种交叉验证方法,
cv=10
-
from sklearn.model_selection import GridSearchCV
gs = GridSearchCV(knn_clf,param_grid,cv=10)
- 拟合网格搜索估计器:
gs.fit(X_train, y_train)
- 查看结果:
gs.best_params_
{'n_neighbors': 3}
gs.cv_results_['mean_test_score']
zip(gs.cv_results_['params'],gs.cv_results_['mean_test_score'])
[({'n_neighbors': 3}, 0.9553571428571429),
({'n_neighbors': 4}, 0.9375),
({'n_neighbors': 5}, 0.9553571428571429),
({'n_neighbors': 6}, 0.9553571428571429),
({'n_neighbors': 7}, 0.9553571428571429),
({'n_neighbors': 8}, 0.9553571428571429)]
它是如何工作的…
在第一章中,我们尝试了蛮力法,即使用 Python 扫描最佳得分:
all_scores = []
for n_neighbors in range(3,9,1):
knn_clf = KNeighborsClassifier(n_neighbors = n_neighbors)
all_scores.append((n_neighbors, cross_val_score(knn_clf, X_train, y_train, cv=10).mean()))
sorted(all_scores, key = lambda x:x[1], reverse = True)
[(3, 0.95666666666666667),
(5, 0.95666666666666667),
(6, 0.95666666666666667),
(7, 0.95666666666666667),
(8, 0.95666666666666667),
(4, 0.94000000000000006)]
这种方法的问题是,它更加耗时且容易出错,尤其是当涉及更多参数或额外的转换(如使用管道时)时。
请注意,网格搜索和蛮力法方法都会扫描所有可能的参数值。
使用 scikit-learn 进行随机搜索
从实际角度来看,RandomizedSearchCV
比常规网格搜索更为重要。因为对于中等大小的数据,或者涉及少量参数的模型,进行完整网格搜索的所有参数组合计算开销太大。
计算资源最好用于非常好地分层采样,或者改进随机化过程。
准备就绪
如前所述,加载鸢尾花数据集的最后两个特征。将数据拆分为训练集和测试集:
from sklearn import datasets
iris = datasets.load_iris()
X = iris.data[:,2:]
y = iris.target
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, stratify = y,random_state = 7)
如何操作…
- 实例化一个最近邻分类器:
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
- 准备一个参数分布,这是进行随机网格搜索时必需的。参数分布是一个字典,包含你希望随机尝试的参数设置:
param_dist = {'n_neighbors': list(range(3,9,1))}
-
实例化一个随机网格搜索并传入以下参数:
-
估算器
-
参数分布
-
一种交叉验证类型,
cv=10
-
运行此过程的次数,
n_iter
-
from sklearn.model_selection import RandomizedSearchCV
rs = RandomizedSearchCV(knn_clf,param_dist,cv=10,n_iter=6)
- 拟合随机网格搜索估算器:
rs.fit(X_train, y_train)
- 查看结果:
rs.best_params_
{'n_neighbors': 3}
zip(rs.cv_results_['params'],rs.cv_results_['mean_test_score'])
[({'n_neighbors': 3}, 0.9553571428571429),
({'n_neighbors': 4}, 0.9375),
({'n_neighbors': 5}, 0.9553571428571429),
({'n_neighbors': 6}, 0.9553571428571429),
({'n_neighbors': 7}, 0.9553571428571429),
({'n_neighbors': 8}, 0.9553571428571429)]
- 在这种情况下,我们实际上对所有六个参数进行了网格搜索。你本可以扫描更大的参数空间:
param_dist = {'n_neighbors': list(range(3,50,1))}
rs = RandomizedSearchCV(knn_clf,param_dist,cv=10,n_iter=15)
rs.fit(X_train,y_train)
rs.best_params_
{'n_neighbors': 16}
- 尝试使用 IPython 计时此过程:
%timeit rs.fit(X_train,y_train)
1 loop, best of 3: 1.06 s per loop
- 计时网格搜索过程:
from sklearn.model_selection import GridSearchCV
param_grid = {'n_neighbors': list(range(3,50,1))}
gs = GridSearchCV(knn_clf,param_grid,cv=10)
%timeit gs.fit(X_train,y_train)
1 loop, best of 3: 3.24 s per loop
- 查看网格搜索的最佳参数:
gs.best_params_
{'n_neighbors': 3}
- 结果表明,3-最近邻的得分与 16-最近邻相同:
zip(gs.cv_results_['params'],gs.cv_results_['mean_test_score'])
[({'n_neighbors': 3}, 0.9553571428571429),
({'n_neighbors': 4}, 0.9375),
...
({'n_neighbors': 14}, 0.9553571428571429),
({'n_neighbors': 15}, 0.9553571428571429),
({'n_neighbors': 16}, 0.9553571428571429),
({'n_neighbors': 17}, 0.9464285714285714),
...
因此,我们在三分之一的时间内得到了相同的分数。
是否使用随机搜索,这是你需要根据具体情况做出的决定。你应该使用随机搜索来尝试了解某个算法的表现。可能无论参数如何,算法的表现都很差,这时你可以换一个算法。如果算法表现非常好,可以使用完整的网格搜索来寻找最佳参数。
此外,除了专注于穷举搜索,你还可以通过集成、堆叠或混合一组合理表现良好的算法来进行尝试。
分类指标
在本章早些时候,我们探讨了基于邻居数量 n_neighbors
参数选择几个最近邻实例的最佳方法。这是最近邻分类中的主要参数:基于 KNN 的标签对一个点进行分类。所以,对于 3-最近邻,根据三个最近点的标签对一个点进行分类。对这三个最近点进行多数投票。
该分类指标在此案例中是内部指标 accuracy_score
,定义为正确分类的数量除以分类总数。还有其他指标,我们将在这里进行探讨。
准备就绪
- 首先,从 UCI 数据库加载 Pima 糖尿病数据集:
import pandas as pd
data_web_address = "https://archive.ics.uci.edu/ml/machine-learning-databases/pima-indians-diabetes/pima-indians-diabetes.data"
column_names = ['pregnancy_x',
'plasma_con',
'blood_pressure',
'skin_mm',
'insulin',
'bmi',
'pedigree_func',
'age',
'target']
feature_names = column_names[:-1]
all_data = pd.read_csv(data_web_address , names=column_names)
- 将数据拆分为训练集和测试集:
import numpy as np
import pandas as pd
X = all_data[feature_names]
y = all_data['target']
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7,stratify=y)
- 回顾上一部分,使用 KNN 算法运行随机搜索:
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import RandomizedSearchCV
knn_clf = KNeighborsClassifier()
param_dist = {'n_neighbors': list(range(3,20,1))}
rs = RandomizedSearchCV(knn_clf,param_dist,cv=10,n_iter=17)
rs.fit(X_train, y_train)
- 然后显示最佳准确率得分:
rs.best_score_
0.75407166123778502
- 此外,查看测试集上的混淆矩阵:
y_pred = rs.predict(X_test)
from sklearn.metrics import confusion_matrix
confusion_matrix(y_test, y_pred)
array([[84, 16],
[27, 27]])
混淆矩阵提供了更具体的关于模型表现的信息。有 27 次模型预测某人没有糖尿病,尽管他们实际上有。这比 16 个被认为有糖尿病的人实际上没有糖尿病的错误更为严重。
在这种情况下,我们希望最大化灵敏度或召回率。在检查线性模型时,我们查看了召回率或灵敏度的定义:
因此,在这种情况下,灵敏度得分为 27/ (27 + 27) = 0.5。使用 scikit-learn,我们可以方便地按如下方式计算此值。
如何操作…
- 从 metrics 模块导入
recall_score
。使用y_test
和y_pred
测量集合的灵敏度:
from sklearn.metrics import recall_score
recall_score(y_test, y_pred)
0.5
我们恢复了之前手动计算的召回得分。在随机搜索中,我们本可以使用 recall_score
来找到具有最高召回率的最近邻实例。
- 导入
make_scorer
并使用带有两个参数的函数recall_score
和greater_is_better
:
from sklearn.metrics import make_scorer
recall_scorer = make_scorer(recall_score, greater_is_better=True)
- 现在执行随机网格搜索:
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import RandomizedSearchCV
knn_clf = KNeighborsClassifier()
param_dist = {'n_neighbors': list(range(3,20,1))}
rs = RandomizedSearchCV(knn_clf,param_dist,cv=10,n_iter=17,scoring=recall_scorer)
rs.fit(X_train, y_train)
- 现在查看最高得分:
rs.best_score_
0.5649632669176643
- 查看召回得分:
y_pred = rs.predict(X_test)
recall_score(y_test,y_pred)
0.5
- 结果与之前相同。在随机搜索中,你本可以尝试
roc_auc_score
,即曲线下面积(ROC AUC):
from sklearn.metrics import roc_auc_score
rs = RandomizedSearchCV(knn_clf,param_dist,cv=10,n_iter=17,scoring=make_scorer(roc_auc_score,greater_is_better=True))
rs.fit(X_train, y_train)
rs.best_score_
0.7100264217324479
还有更多…
你可以为分类设计自己的评分器。假设你是一个保险公司,并且你为混淆矩阵中的每个单元格分配了成本。相对成本如下:
我们正在查看的混淆矩阵的成本可以按如下方式计算:
costs_array = confusion_matrix(y_test, y_pred) * np.array([[1,2],
[100,20]])
costs_array
array([[ 84, 32],
[2700, 540]])
现在加总总成本:
costs_array.sum()
3356
现在将其放入评分器中并运行网格搜索。评分器中的参数 greater_is_better
设置为 False
,因为成本应尽可能低:
def costs_total(y_test, y_pred):
return (confusion_matrix(y_test, y_pred) * np.array([[1,2],
[100,20]])).sum()
costs_scorer = make_scorer(costs_total, greater_is_better=False)
rs = RandomizedSearchCV(knn_clf,param_dist,cv=10,n_iter=17,scoring=costs_scorer)
rs.fit(X_train, y_train)
rs.best_score_
-1217.5879478827362
得分为负,因为当 make_scorer
函数中的 greater_is_better
参数为 false 时,得分会乘以 -1
。网格搜索试图最大化该得分,从而最小化得分的绝对值。
测试集的成本如下:
costs_total(y_test,rs.predict(X_test))
3356
在查看这个数字时,别忘了查看测试集中涉及的个体数量,共有 154 人。每个人的平均成本约为 21.8 美元。
回归指标
使用回归指标的交叉验证在 scikit-learn 中非常简单。可以从 sklearn.metrics
导入评分函数并将其放入 make_scorer
函数中,或者你可以为特定的数据科学问题创建自定义评分器。
准备就绪
加载一个使用回归指标的数据集。我们将加载波士顿房价数据集并将其拆分为训练集和测试集:
from sklearn.datasets import load_boston
boston = load_boston()
X = boston.data
y = boston.target
from sklearn.model_selection import train_test_split, cross_val_score
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7)
我们对数据集了解不多。我们可以尝试使用高方差算法进行快速网格搜索:
from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import RandomizedSearchCV
knn_reg = KNeighborsRegressor()
param_dist = {'n_neighbors': list(range(3,20,1))}
rs = RandomizedSearchCV(knn_reg,param_dist,cv=10,n_iter=17)
rs.fit(X_train, y_train)
rs.best_score_
0.46455839325055914
尝试一个不同的模型,这次是一个线性模型:
from sklearn.linear_model import Ridge
cross_val_score(Ridge(),X_train,y_train,cv=10).mean()
0.7439511908709866
默认情况下,两个回归器都衡量 r2_score
,即 R 平方,因此线性模型更好。尝试一个不同的复杂模型,一个树的集成:
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
cross_val_score(GradientBoostingRegressor(max_depth=7),X_train,y_train,cv=10).mean()
0.83082671732165492
集成模型的表现更好。你也可以尝试随机森林:
cross_val_score(RandomForestRegressor(),X_train,y_train,cv=10).mean()
0.82474734196711685
现在我们可以通过最大化内部 R-squared 梯度提升评分器,专注于利用当前评分机制来改进梯度提升。尝试进行一两次随机搜索。这是第二次搜索:
param_dist = {'n_estimators': [4000], 'learning_rate': [0.01], 'max_depth':[1,2,3,5,7]}
rs_inst_a = RandomizedSearchCV(GradientBoostingRegressor(), param_dist, n_iter = 5, n_jobs=-1)
rs_inst_a.fit(X_train, y_train)
为 R-squared 优化返回了以下结果:
rs_inst_a.best_params_
{'learning_rate': 0.01, 'max_depth': 3, 'n_estimators': 4000}
rs_inst_a.best_score_
0.88548410382780185
梯度提升中的树的深度为三。
如何做到…
现在我们将执行以下操作:
-
创建一个评分函数。
-
使用该函数创建一个 scikit-scorer。
-
运行网格搜索以找到最佳的梯度提升参数,最小化误差函数。
让我们开始:
- 使用 Numba 即时编译(JIT)编译器创建平均百分比误差评分函数。原始的 NumPy 函数如下:
def mape_score(y_test, y_pred):
return (np.abs(y_test - y_pred)/y_test).mean()
- 让我们使用 Numba JIT 编译器重写这个函数,稍微加速一些。你可以用类似 C 的代码,通过 Numba 按位置索引数组:
from numba import autojit
@autojit
def mape_score(y_test, y_pred):
sum_total = 0
y_vec_length = len(y_test)
for index in range(y_vec_length):
sum_total += (1 - (y_pred[index]/y_test[index]))
return sum_total/y_vec_length
- 现在创建一个评分器。得分越低越好,不像 R-squared,那里的得分越高越好:
from sklearn.metrics import make_scorer
mape_scorer = make_scorer(mape_score, greater_is_better=False)
- 现在进行网格搜索:
param_dist = {'n_estimators': [4000], 'learning_rate': [0.01], 'max_depth':[1,2,3,4,5]}
rs_inst_b = RandomizedSearchCV(GradientBoostingRegressor(), param_dist, n_iter = 3, n_jobs=-1,scoring = mape_scorer)
rs_inst_b.fit(X_train, y_train)
rs_inst_b.best_score_
0.021086502313661441
rs_inst_b.best_params_
{'learning_rate': 0.01, 'max_depth': 1, 'n_estimators': 4000}
使用此度量,最佳得分对应于深度为 1 的梯度提升树。
聚类度量
衡量聚类算法的性能比分类或回归要复杂一些,因为聚类是无监督机器学习。幸运的是,scikit-learn 已经非常直接地为我们提供了帮助。
准备开始
为了衡量聚类性能,首先加载鸢尾花数据集。我们将鸢尾花重新标记为两种类型:当目标是 0 时为类型 0,当目标是 1 或 2 时为类型 1:
from sklearn.datasets import load_iris
import numpy as np
iris = load_iris()
X = iris.data
y = np.where(iris.target == 0,0,1)
如何做到…
- 实例化一个 k-means 算法并训练它。由于该算法是聚类算法,因此在训练时不要使用目标值:
from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=2,random_state=0)
kmeans.fit(X)
- 现在导入所有必要的库,通过交叉验证对 k-means 进行评分。我们将使用
adjusted_rand_score
聚类性能指标:
from sklearn.metrics.cluster import adjusted_rand_score
from sklearn.metrics import make_scorer
from sklearn.model_selection import cross_val_score
cross_val_score(kmeans,X,y,cv=10,scoring=make_scorer(adjusted_rand_score)).mean()
0.8733695652173914
评估聚类算法与评估分类算法非常相似。
使用虚拟估算器进行结果比较
本食谱是关于创建虚拟估算器的;这不是很华丽或令人兴奋的部分,但它为你最终构建的模型提供了一个参考点,值得一做。
准备开始
在这个食谱中,我们将执行以下任务:
-
创建一些随机数据。
-
拟合各种虚拟估算器。
我们将对回归数据和分类数据执行这两步操作。
如何做到…
- 首先,我们将创建随机数据:
from sklearn.datasets import make_regression, make_classification
X, y = make_regression()
from sklearn import dummy
dumdum = dummy.DummyRegressor()
dumdum.fit(X, y)
DummyRegressor(constant=None, quantile=None, strategy='mean')
- 默认情况下,估算器将通过取值的均值并多次输出它来进行预测:
dumdum.predict(X)[:5]
>array([-25.0450033, -25.0450033, -25.0450033, -25.0450033, -25.0450033])
还有两种其他策略可以尝试。我们可以预测一个提供的常量(参考此食谱中的第一条命令块中的 constant=None
)。我们还可以预测中位数值。仅当策略为常量时才会考虑提供常量。
- 我们来看看:
predictors = [("mean", None),
("median", None),
("constant", 10)]
for strategy, constant in predictors:
dumdum = dummy.DummyRegressor(strategy=strategy,
constant=constant)
dumdum.fit(X, y)
print "strategy: {}".format(strategy), ",".join(map(str, dumdum.predict(X)[:5]))
strategy: mean -25.0450032962,-25.0450032962,-25.0450032962,-25.0450032962,-25.0450032962
strategy: median -37.734448002,-37.734448002,-37.734448002,-37.734448002,-37.734448002
strategy: constant 10.0,10.0,10.0,10.0,10.0
- 我们实际上有四种分类器的选择。这些策略与连续情况类似,只不过更加倾向于分类问题:
predictors = [("constant", 0),("stratified", None),("uniform", None),("most_frequent", None)]
#We'll also need to create some classification data:
X, y = make_classification()
for strategy, constant in predictors:
dumdum = dummy.DummyClassifier(strategy=strategy,
constant=constant)
dumdum.fit(X, y)
print "strategy: {}".format(strategy), ",".join(map(str,dumdum.predict(X)[:5]))
strategy: constant 0,0,0,0,0
strategy: stratified 1,0,1,1,1
strategy: uniform 1,0,1,0,1
strategy: most_frequent 0,0,0,0,0
它是如何工作的…
测试你的模型是否能够比最简单的模型表现更好总是个不错的做法,这正是虚拟估计器所能提供的。例如,假设你有一个欺诈检测模型。在这个模型中,数据集中只有 5%是欺诈行为。因此,我们很可能通过仅仅不猜测数据是欺诈的,就能拟合出一个相当不错的模型。
我们可以通过使用分层策略并执行以下命令来创建此模型。我们还可以得到一个很好的例子,说明类别不平衡是如何导致问题的:
X, y = make_classification(20000, weights=[.95, .05])
dumdum = dummy.DummyClassifier(strategy='most_frequent')
dumdum.fit(X, y)
DummyClassifier(constant=None, random_state=None, strategy='most_frequent')
from sklearn.metrics import accuracy_score
print accuracy_score(y, dumdum.predict(X))
0.94615
我们实际上经常是正确的,但这并不是重点。重点是这是我们的基准。如果我们不能创建一个比这个更准确的欺诈检测模型,那么就不值得花费时间。
特征选择
本食谱以及接下来的两个将围绕自动特征选择展开。我喜欢将其视为参数调优的特征类比。就像我们通过交叉验证来寻找合适的参数一样,我们也可以找到一个合适的特征子集。这将涉及几种不同的方法。
最简单的想法是单变量选择。其他方法则涉及特征的组合使用。
特征选择的一个附加好处是,它可以减轻数据收集的负担。假设你已经基于一个非常小的数据子集建立了一个模型。如果一切顺利,你可能想扩大规模,在整个数据子集上预测模型。如果是这种情况,你可以在这个规模上减轻数据收集的工程负担。
准备工作
使用单变量特征选择时,评分函数将再次成为焦点。这一次,它们将定义我们可以用来消除特征的可比度量。
在本食谱中,我们将拟合一个包含大约 10,000 个特征的回归模型,但只有 1,000 个数据点。我们将逐步了解各种单变量特征选择方法:
from sklearn import datasets
X, y = datasets.make_regression(1000, 10000)
现在我们已经有了数据,我们将比较通过各种方法包含的特征。这实际上是你在处理文本分析或某些生物信息学领域时非常常见的情况。
如何操作…
- 首先,我们需要导入
feature_selection
模块:
from sklearn import feature_selection
f, p = feature_selection.f_regression(X, y)
- 这里,
f
是与每个线性模型拟合相关的f
分数,该模型仅使用一个特征。然后我们可以比较这些特征,并根据这种比较来剔除特征。p
是与f
值相关的p
值。在统计学中,p
值是指在给定的检验统计量值下,出现比当前值更极端的结果的概率。在这里,f
值是检验统计量:
f[:5]
array([ 1.23494617, 0.70831694, 0.09219176, 0.14583189, 0.78776466])
p[:5]
array([ 0.26671492, 0.40020473, 0.76147235, 0.7026321 , 0.37499074])
- 正如我们所见,许多
p
值相当大。我们希望p
值尽可能小。因此,我们可以从工具箱中取出 NumPy,选择所有小于.05
的p
值。这些将是我们用于分析的特征:
import numpy as np
idx = np.arange(0, X.shape[1])
features_to_keep = idx[p < .05]
len(features_to_keep)
496
如你所见,我们实际上保留了相对较多的特征。根据模型的上下文,我们可以缩小这个p
值。这将减少保留的特征数量。
另一个选择是使用VarianceThreshold
对象。我们已经了解了一些它的内容,但需要明白的是,我们拟合模型的能力在很大程度上依赖于特征所产生的方差。如果没有方差,那么我们的特征就无法描述因变量的变化。根据文档的说法,它的一个优点是由于它不使用结果变量,因此可以用于无监督的情况。
- 我们需要设定一个阈值,以决定去除哪些特征。为此,我们只需要取特征方差的中位数并提供它:
var_threshold = feature_selection.VarianceThreshold(np.median(np.var(X, axis=1)))
var_threshold.fit_transform(X).shape
(1000L, 4888L)
如我们所见,我们已经去除了大约一半的特征,这也大致符合我们的预期。
它是如何工作的…
一般来说,这些方法都是通过拟合一个只有单一特征的基本模型来工作的。根据我们是分类问题还是回归问题,我们可以使用相应的评分函数。
让我们看看一个较小的问题,并可视化特征选择如何去除某些特征。我们将使用第一个示例中的相同评分函数,但只使用 20 个特征:
X, y = datasets.make_regression(10000, 20)
f, p = feature_selection.f_regression(X, y)
现在让我们绘制特征的 p 值。我们可以看到哪些特征将被去除,哪些特征将被保留:
%matplotlib inline
from matplotlib import pyplot as plt
f, ax = plt.subplots(figsize=(7, 5))
ax.bar(np.arange(20), p, color='k')
ax.set_title("Feature p values")
如我们所见,许多特征将不会被保留,但有一些特征会被保留。
基于 L1 范数的特征选择
我们将使用一些与 LASSO 回归配方中看到的类似的思想。在那个配方中,我们查看了具有零系数的特征数量。现在我们将更进一步,利用与 L1 范数相关的稀疏性来预处理特征。
准备工作
我们将使用糖尿病数据集来进行回归拟合。首先,我们将使用 ShuffleSplit 交叉验证拟合一个基本的线性回归模型。完成后,我们将使用 LASSO 回归来找到系数为零的特征,这些特征在使用 L1 惩罚时会被去除。这有助于我们避免过拟合(即模型过于专门化,无法适应它未训练过的数据)。换句话说,如果模型过拟合,它对外部数据的泛化能力较差。
我们将执行以下步骤:
-
加载数据集。
-
拟合一个基本的线性回归模型。
-
使用特征选择去除不具信息量的特征。
-
重新拟合线性回归模型,并检查它与完全特征模型相比的拟合效果。
如何操作…
- 首先,让我们获取数据集:
import sklearn.datasets as ds
diabetes = ds.load_diabetes()
X = diabetes.data
y = diabetes.target
- 让我们创建
LinearRegression
对象:
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
- 让我们从 metrics 模块导入
mean_squared_error
函数和make_scorer
包装器。从model_selection
模块,导入ShuffleSplit
交叉验证方案和cross_val_score
交叉验证评分器。接下来,使用mean_squared_error
度量来评分该函数:
from sklearn.metrics import make_scorer, mean_squared_error
from sklearn.model_selection import cross_val_score,ShuffleSplit
shuff = ShuffleSplit(n_splits=10, test_size=0.25, random_state=0)
score_before = cross_val_score(lr,X,y,cv=shuff,scoring=make_scorer(mean_squared_error,greater_is_better=False)).mean()
score_before
-3053.393446308266
- 现在我们已经有了常规拟合,让我们在去除系数为零的特征后检查一下。让我们拟合 LASSO 回归:
from sklearn.linear_model import LassoCV
lasso_cv = LassoCV()
lasso_cv.fit(X,y)
lasso_cv.coef_
array([ -0\. , -226.2375274 , 526.85738059, 314.44026013,
-196.92164002, 1.48742026, -151.78054083, 106.52846989,
530.58541123, 64.50588257])
- 我们将删除第一个特征。我将使用 NumPy 数组来表示要包含在模型中的列:
import numpy as np
columns = np.arange(X.shape[1])[lasso_cv.coef_ != 0]
columns
- 好的,现在我们将使用特定的特征来拟合模型(请参见以下代码块中的列):
score_afterwards = cross_val_score(lr,X[:,columns],y,cv=shuff, scoring=make_scorer(mean_squared_error,greater_is_better=False)).mean()
score_afterwards
-3033.5012859289677
之后的得分并没有比之前好多少,尽管我们已经消除了一个无信息特征。在*还有更多内容…*部分,我们将看到一个额外的示例。
还有更多内容…
- 首先,我们将创建一个具有许多无信息特征的回归数据集:
X, y = ds.make_regression(noise=5)
- 创建一个
ShuffleSplit
实例,进行 10 次迭代,n_splits=10
。测量普通线性回归的交叉验证得分:
shuff = ShuffleSplit(n_splits=10, test_size=0.25, random_state=0)
score_before = cross_val_score(lr,X,y,cv=shuff, scoring=make_scorer(mean_squared_error,greater_is_better=False)).mean()
- 实例化
LassoCV
来消除无信息的列:
lasso_cv = LassoCV()
lasso_cv.fit(X,y)
- 消除无信息的列。查看最终得分:
columns = np.arange(X.shape[1])[lasso_cv.coef_ != 0]
score_afterwards = cross_val_score(lr,X[:,columns],y,cv=shuff, scoring=make_scorer(mean_squared_error,greater_is_better=False)).mean()
print "Score before:",score_before
print "Score after: ",score_afterwards
Score before: -8891.35368845
Score after: -22.3488585347
在我们移除无信息特征后,最后的拟合效果要好得多。
使用 joblib 或 pickle 持久化模型
在这个教程中,我们将展示如何将模型保留以供以后使用。例如,你可能希望使用一个模型来预测结果并自动做出决策。
做好准备
创建数据集并训练分类器:
from sklearn.datasets import make_classification
from sklearn.tree import DecisionTreeClassifier
X, y = make_classification()
dt = DecisionTreeClassifier()
dt.fit(X, y)
如何做…
- 使用 joblib 保存分类器所做的训练工作:
from sklearn.externals import joblib
joblib.dump(dt, "dtree.clf")
['dtree.clf']
打开已保存的模型
- 使用 joblib 加载模型。使用一组输入进行预测:
from sklearn.externals import joblib
pulled_model = joblib.load("dtree.clf")
y_pred = pulled_model.predict(X)
我们不需要重新训练模型,并且节省了大量训练时间。我们只是使用 joblib 重新加载了模型并进行了预测。
还有更多内容…
你也可以在 Python 2.x 中使用cPickle
模块,或者在 Python 3.x 中使用pickle
模块。就个人而言,我使用这个模块处理几种类型的 Python 类和对象:
- 首先导入
pickle
:
import cPickle as pickle #Python 2.x
# import pickle # Python 3.x
- 使用
dump()
模块方法。它有三个参数:要保存的数据、保存目标文件和 pickle 协议。以下代码将训练好的树保存到dtree.save
文件:
f = open("dtree.save", 'wb')
pickle.dump(dt,f, protocol = pickle.HIGHEST_PROTOCOL)
f.close()
- 按如下方式打开
dtree.save
文件:
f = open("dtree.save", 'rb')
return_tree = pickle.load(f)
f.close()
- 查看树:
return_tree
DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=None,
max_features=None, max_leaf_nodes=None,
min_impurity_split=1e-07, min_samples_leaf=1,
min_samples_split=2, min_weight_fraction_leaf=0.0,
presort=False, random_state=None, splitter='best'
第八章:支持向量机
在本章中,我们将涵盖以下内容:
-
使用线性 SVM 进行数据分类
-
优化 SVM
-
使用 SVM 进行多类分类
-
支持向量回归
介绍
在本章中,我们将首先使用支持向量机(SVM)与线性核,以大致了解 SVM 的工作原理。它们创建一个超平面,或在多个维度中的线性面,最佳地分隔数据。
在二维空间中,这很容易看出:超平面是分隔数据的直线。我们将看到 SVM 的系数和截距数组。它们一起唯一地描述了一个scikit-learn
线性 SVC 预测器。
在本章的其余部分,SVM 使用径向基函数(RBF)核。它们是非线性的,但具有平滑的分隔面。在实际应用中,SVM 在许多数据集上表现良好,因此是scikit-learn
库的一个重要组成部分。
使用线性 SVM 进行数据分类
在第一章中,我们看到了一些使用 SVM 进行分类的示例。我们重点讨论了 SVM 在分类性能上略优于逻辑回归,但大部分时间我们并未深入探讨 SVM。
在这里,我们将更仔细地关注它们。虽然 SVM 没有容易的概率解释,但它们有一个直观的几何解释。线性 SVM 的主要思想是通过最佳的平面分隔两个类。
让我们使用 SVM 对两个类进行线性分隔。
准备工作
让我们从加载并可视化scikit-learn
中提供的鸢尾花数据集开始:
加载数据
加载部分鸢尾花数据集。这将使我们能够与第一章进行轻松的比较:
#load the libraries we have been using
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt #Library for visualization
from sklearn import datasets
iris = datasets.load_iris()
X_w = iris.data[:, :2] #load the first two features of the iris data
y_w = iris.target #load the target of the iris data
现在,我们将使用 NumPy 掩码来关注前两个类:
#select only the first two classes for both the feature set and target set
#the first two classes of the iris dataset: Setosa (0), Versicolour (1)
X = X_w[y_w < 2]
y = y_w[y_w < 2]
可视化这两个类
使用 matplotlib 绘制0
和1
类。请记住,X_0[:,0]
表示 NumPy 数组的第一列。
在以下代码中,X_0
表示与目标y
为0
相对应的输入子集,而X_1
是目标值为1
的匹配子集:
X_0 = X[y == 0]
X_1 = X[y == 1]
#to visualize within IPython:
%matplotlib inline
plt.figure(figsize=(10,7)) #change figure-size for easier viewing
plt.scatter(X_0[:,0],X_0[:,1], color = 'red')
plt.scatter(X_1[:,0],X_1[:,1], color = 'blue')
从图表中可以清楚地看出,我们可以找到一条直线来分隔这两个类。
如何操作…
找到 SVM 直线的过程很简单。这与任何scikit-learn
的监督学习估计器的过程相同:
-
创建训练集和测试集。
-
创建 SVM 模型实例。
-
将 SVM 模型拟合到加载的数据。
-
使用 SVM 模型进行预测,并在准备好对未见数据进行预测之前,衡量模型的性能。
让我们开始吧:
- 将前两个类的前两个特征的数据集进行划分。对目标集进行分层:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7,stratify=y)
- 创建 SVM 模型实例。将核设置为线性,因为我们希望有一条线来分隔这个例子中涉及的两个类:
from sklearn.svm import SVC
svm_inst = SVC(kernel='linear')
- 拟合模型(训练模型):
svm_inst.fit(X_train,y_train)
- 使用测试集进行预测:
y_pred = svm_inst.predict(X_test)
- 测量 SVM 在测试集上的表现:
from sklearn.metrics import accuracy_score
accuracy_score(y_test, y_pred)
1.0
它在测试集上表现得非常好。这并不奇怪,因为当我们可视化每个类别时,它们很容易被视觉上分开。
- 通过在二维网格上使用估算器,来可视化决策边界,即分隔类别的直线:
from itertools import product
#Minima and maxima of both features
xmin, xmax = np.percentile(X[:, 0], [0, 100])
ymin, ymax = np.percentile(X[:, 1], [0, 100])
#Grid/Cartesian product with itertools.product
test_points = np.array([[xx, yy] for xx, yy in product(np.linspace(xmin, xmax), np.linspace(ymin, ymax))])
#Predictions on the grid
test_preds = svm_inst.predict(test_points)
- 通过为预测着色来绘制网格。请注意,我们已经修改了之前的可视化图像,加入了 SVM 的预测:
X_0 = X[y == 0]
X_1 = X[y == 1]
%matplotlib inline
plt.figure(figsize=(10,7)) #change figure-size for easier viewing
plt.scatter(X_0[:,0],X_0[:,1], color = 'red')
plt.scatter(X_1[:,0],X_1[:,1], color = 'blue')
colors = np.array(['r', 'b'])
plt.scatter(test_points[:, 0], test_points[:, 1], color=colors[test_preds], alpha=0.25)
plt.scatter(X[:, 0], X[:, 1], color=colors[y])
plt.title("Linearly-separated classes")
我们通过在二维网格上进行预测,详细描述了 SVM 线性决策边界。
它是如何工作的……
有时,在整个网格上进行预测计算可能非常昂贵,尤其是当 SVM 预测许多类别且维度较高时。在这种情况下,您将需要访问 SVM 决策边界的几何信息。
线性决策边界,一个超平面,是由一个垂直于超平面的向量和一个截距唯一确定的。法向量包含在 SVM 实例的coef_ data
属性中。截距包含在 SVM 实例的intercept_ data
属性中。查看这两个属性:
svm_inst.coef_
array([[ 2.22246001, -2.2213921 ]])
svm_inst.intercept_
array([-5.00384439])
您可能会很快发现,coef_[0]
向量垂直于我们绘制的分隔两个鸢尾花类别的直线。
每次,这两个 NumPy 数组 svm_inst.coef_
和 svm_inst.intercept_
的行数是相同的。每一行对应一个分隔相关类别的平面。在这个例子中,两个类别通过一个超平面线性分开。特定的 SVM 类型,SVC 在这种情况下实现了一个一对一分类器:它会绘制一个唯一的平面来分隔每一对类别。
如果我们尝试分离三个类别,那么有三种可能的组合,3 x 2/2 = 3。对于 n 个类别,SVC 提供的平面数如下:
coef_ data
属性中的列数是数据中特征的数量,在本例中是两个。
要找到关于空间中某一点的决策,求解以下方程为零:
如果您只关心平面的唯一性,可以存储元组(coef_, intercept_)
。
还有更多……
此外,您还可以查看实例的参数以了解更多信息:
svm_inst
SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
decision_function_shape=None, degree=3, gamma='auto', kernel='linear',
max_iter=-1, probability=False, random_state=None, shrinking=True,
tol=0.001, verbose=False)
传统上,SVC 预测性能通过以下参数进行优化:C、gamma 和核的形状。C 描述了 SVM 的边距,默认设置为 1。边距是超平面两侧没有类别示例的空白区域。如果数据集有许多噪声观察值,可以尝试使用交叉验证来提高 C 的值。C 与边距上的错误成正比,随着 C 值的增大,SVM 将尝试使边距更小。
关于 SVM 的最后一点是,我们可以重新缩放数据,并通过交叉验证测试该缩放效果。方便的是,鸢尾花数据集中的所有输入单位都是厘米,所以不需要重新缩放,但对于任意数据集,你应该考虑这个问题。
优化 SVM
在本示例中,我们将继续使用鸢尾花数据集,但使用两种难以区分的品种——变色鸢尾和维吉尼卡鸢尾。
在本节中,我们将重点关注以下内容:
-
设置 scikit-learn 管道:一系列变换,最后是一个预测模型
-
网格搜索:对多个版本的支持向量机(SVM)进行性能扫描,并改变其参数
准备工作
加载鸢尾花数据集中的两个类别和两个特征:
#load the libraries we have been using
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import datasets
iris = datasets.load_iris()
X_w = iris.data[:, :2] #load the first two features of the iris data
y_w = iris.target #load the target of the iris data
X = X_w[y_w != 0]
y = y_w[y_w != 0]
X_1 = X[y == 1]
X_2 = X[y == 2]
如何操作…
- 首先将数据分为训练集和测试集:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7,stratify=y)
构建一个管道
- 然后构建一个包含两个步骤的管道:一个缩放步骤和一个 SVM 步骤。在将数据传递给 SVM 之前,最好先进行缩放:
from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
svm_est = Pipeline([('scaler',StandardScaler()),('svc',SVC())])
请注意,在管道中,缩放步骤的名称是scaler
,SVM 的名称是svc
。这些名称在接下来的步骤中将非常关键。还要注意,默认的 SVM 是 RBF SVM,它是非线性的。
为管道构建参数网格
- 以对数方式改变相关的 RBF 参数 C 和 gamma,每次改变一个数量级:
Cs = [0.001, 0.01, 0.1, 1, 10]
gammas = [0.001, 0.01, 0.1, 1, 10]
- 最后,将参数网格构建成字典。SVM 参数字典的键名以
svc__
开头,取管道 SVM 的名称并加上两个下划线。然后是 SVM 估计器内的参数名称,C
和gamma
:
param_grid = dict(svc__gamma=gammas, svc__C=Cs)
提供交叉验证方案
- 以下是一个分层且经过洗牌的拆分。
n_splits
参数指的是数据集将被拆分成的折叠数或尝试次数。test_size
参数则指每个折叠中留出来用于测试的数据量。估计器将在每个折叠中使用测试集评分:
from sklearn.model_selection import StratifiedShuffleSplit
cv = StratifiedShuffleSplit(n_splits=5, test_size=0.2, random_state=7)
分层洗牌的最重要元素是,每个折叠都保持每个类别样本的比例。
- 对于普通的交叉验证方案,将
cv
设置为一个整数,表示折叠的数量:
cv = 10
执行网格搜索
网格搜索需要三个必需的元素:
-
估计器
-
参数网格
-
一个交叉验证方案
- 我们有这三项元素。设置网格搜索,并在训练集上运行:
from sklearn.model_selection import GridSearchCV
grid_cv = GridSearchCV(svm_est, param_grid=param_grid, cv=cv)
grid_cv.fit(X_train, y_train)
- 查找通过网格搜索找到的最佳参数:
grid_cv.best_params_
{'svc__C': 10, 'svc__gamma': 0.1}
- 查找与最佳估计器相关的最佳得分:
grid_cv.best_score_
0.71250000000000002
还有更多内容…
让我们从其他角度看一下 SVM 分类。
随机网格搜索替代方案
scikit-learn 的GridSearchCV
会执行一个完整的扫描,以寻找估计器的最佳参数集。在此情况下,它会搜索由param_grid
参数指定的 5 x 5 = 25(C,gamma)对。
另一种选择是使用RandomizedSearchCV
,通过使用以下这一行代替GridSearchCV
所用的那一行:
from sklearn.model_selection import RandomizedSearchCV
rand_grid = RandomizedSearchCV(svm_est, param_distributions=param_grid, cv=cv,n_iter=10)
rand_grid.fit(X_train, y_train)
它得到了相同的C
和gamma
:
rand_grid.best_params_
{'svc__C': 10, 'svc__gamma': 0.001}
可视化非线性 RBF 决策边界
使用类似于之前配方的代码可视化 RBF 决策边界。首先,创建一个网格并预测网格上每个点对应的类别:
from itertools import product
#Minima and maxima of both features
xmin, xmax = np.percentile(X[:, 0], [0, 100])
ymin, ymax = np.percentile(X[:, 1], [0, 100])
#Grid/Cartesian product with itertools.product
test_points = np.array([[xx, yy] for xx, yy in product(np.linspace(xmin, xmax), np.linspace(ymin, ymax))])
#Predictions on the grid
test_preds = grid_cv.predict(test_points)
现在可视化网格:
X_1 = X[y == 1]
X_2 = X[y == 2]
%matplotlib inline
plt.figure(figsize=(10,7)) #change figure-size for easier viewing
plt.scatter(X_2[:,0],X_2[:,1], color = 'red')
plt.scatter(X_1[:,0],X_1[:,1], color = 'blue')
colors = np.array(['r', 'b'])
plt.scatter(test_points[:, 0], test_points[:, 1], color=colors[test_preds-1], alpha=0.25)
plt.scatter(X[:, 0], X[:, 1], color=colors[y-1])
plt.title("RBF-separated classes")
请注意,在结果图中,RBF 曲线看起来相当直,但实际上它对应的是一条轻微的曲线。这是一个 gamma = 0.1 和 C = 0.001 的 SVM:
C 和 gamma 的更多含义
更直观地说,gamma 参数决定了单个样本对每单位距离的影响程度。如果 gamma 较低,则样本在较长距离处具有影响。如果 gamma 较高,则它们的影响仅限于较短的距离。SVM 在其实现中选择支持向量,gamma 与这些向量的影响半径成反比。
关于 C,较低的 C 会使决策面更加平滑,而较高的 C 会使 SVM 尝试正确分类所有样本,导致不太平滑的决策面。
使用 SVM 进行多类别分类
我们将扩展前面的配方,通过两个特征对所有鸢尾花类型进行分类。这不是二分类问题,而是多分类问题。这些步骤是在之前配方的基础上扩展的。
准备就绪
对于多分类问题,SVC 分类器(scikit 的 SVC)可以稍作修改。为此,我们将使用鸢尾数据集中的所有三个类别。
为每个类别加载两个特征:
#load the libraries we have been using
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import datasets
iris = datasets.load_iris()
X = iris.data[:, :2] #load the first two features of the iris data
y = iris.target #load the target of the iris data
X_0 = X[y == 0]
X_1 = X[y == 1]
X_2 = X[y == 2]
将数据拆分为训练集和测试集:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7,stratify=y)
如何操作…
OneVsRestClassifier
- 在管道中加载
OneVsRestClassifier
:
from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.multiclass import OneVsRestClassifier
svm_est = Pipeline([('scaler',StandardScaler()),('svc',OneVsRestClassifier(SVC()))])
- 设置参数网格:
Cs = [0.001, 0.01, 0.1, 1, 10]
gammas = [0.001, 0.01, 0.1, 1, 10]
- 构建参数网格。注意,表示
OneVsRestClassifier
SVC 的语法非常特殊。当在管道中命名为svc
时,字典中的参数键名以svc__estimator__
开头:
param_grid = dict(svc__estimator__gamma=gammas, svc__estimator__C=Cs)
- 加载一个随机化的超参数搜索。拟合它:
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import StratifiedShuffleSplit
cv = StratifiedShuffleSplit(n_splits=5, test_size=0.2, random_state=7)
rand_grid = RandomizedSearchCV(svm_est, param_distributions=param_grid, cv=cv,n_iter=10)
rand_grid.fit(X_train, y_train)
- 查找最佳参数:
rand_grid.best_params_
{'svc__estimator__C': 10, 'svc__estimator__gamma': 0.1}
可视化它
我们将通过调用训练好的 SVM 来预测二维网格中每个点的类别:
%matplotlib inline
from itertools import product
#Minima and maxima of both features
xmin, xmax = np.percentile(X[:, 0], [0, 100])
ymin, ymax = np.percentile(X[:, 1], [0, 100])
#Grid/Cartesian product with itertools.product
test_points = np.array([[xx, yy] for xx, yy in product(np.linspace(xmin, xmax,100), np.linspace(ymin, ymax,100))])
#Predictions on the grid
test_preds = rand_grid.predict(test_points)
plt.figure(figsize=(15,9)) #change figure-size for easier viewing
plt.scatter(X_0[:,0],X_0[:,1], color = 'green')
plt.scatter(X_1[:,0],X_1[:,1], color = 'blue')
plt.scatter(X_2[:,0],X_2[:,1], color = 'red')
colors = np.array(['g', 'b', 'r'])
plt.tight_layout()
plt.scatter(test_points[:, 0], test_points[:, 1], color=colors[test_preds], alpha=0.25)
plt.scatter(X[:, 0], X[:, 1], color=colors[y])
SVM 生成的边界通常是平滑曲线,这与我们将在接下来的章节中看到的基于树的边界非常不同。
如何操作…
OneVsRestClassifier
创建许多二元 SVC 分类器:每个类别与其他类别进行对比。在这种情况下,将计算三个决策边界,因为有三个类别。这种类型的分类器很容易理解,因为决策边界和面较少。
如果有 10 个类别,使用默认的 OneVsOneClassifier
(SVC)会有 10 x 9/2 = 45 个边界。另一方面,使用 OneVsAllClassifier
会有 10 个边界。
支持向量回归
我们将利用 SVM 分类的配方,在 scikit-learn 的糖尿病数据集上执行支持向量回归。
准备就绪
加载糖尿病数据集:
#load the libraries we have been using
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import datasets
diabetes = datasets.load_diabetes()
X = diabetes.data
y = diabetes.target
将数据划分为训练集和测试集。此情况下回归问题没有分层:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7)
如何操作…
- 在管道中创建一个
OneVsRestClassifier
,并从sklearn.svm
导入支持向量回归(SVR):
from sklearn.svm import SVR
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.multiclass import OneVsRestClassifier
svm_est = Pipeline([('scaler',StandardScaler()),('svc',OneVsRestClassifier(SVR()))])
- 创建一个参数网格:
Cs = [0.001, 0.01, 0.1, 1]
gammas = [0.001, 0.01, 0.1]
param_grid = dict(svc__estimator__gamma=gammas, svc__estimator__C=Cs)
- 执行随机搜索以寻找最佳超参数,C 和 gamma:
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import StratifiedShuffleSplit
rand_grid = RandomizedSearchCV(svm_est, param_distributions=param_grid, cv=5,n_iter=5,scoring='neg_mean_absolute_error')
rand_grid.fit(X_train, y_train)
- 查看最佳参数:
rand_grid.best_params_
{'svc__estimator__C': 10, 'svc__estimator__gamma': 0.1}
- 查看最佳分数:
rand_grid.best_score_
-58.059490084985839
分数似乎不是很好。尝试不同的算法和不同的评分设置,看看哪个表现最好。
第九章:树算法与集成方法
本章将涵盖以下内容:
-
使用决策树进行基本分类
-
使用 pydot 可视化决策树
-
调整决策树
-
使用决策树进行回归
-
使用交叉验证减少过拟合
-
实现随机森林回归
-
使用最近邻法进行包外回归
-
调整梯度提升树
-
调整 AdaBoost 回归器
-
使用 scikit-learn 编写堆叠聚合器
介绍
本章重点讨论决策树和集成算法。决策算法容易解释和可视化,因为它们是我们熟悉的决策过程的概述。集成方法可以部分解释和可视化,但它们包含许多部分(基础估计器),因此我们不能总是轻松地读取它们。
集成学习的目标是多个估计器比单个估计器表现更好。scikit-learn 中实现了两种集成方法:平均方法和提升方法。平均方法(如随机森林、包外法、额外树)通过平均多个估计器的预测来减少方差。提升方法(如梯度提升和 AdaBoost)通过依次构建基础估计器来减少偏差,从而减少整个集成方法的偏差。
许多集成方法的共同特点是使用随机性来构建预测器。例如,随机森林使用随机性(正如其名字所示),我们也将通过许多模型参数的搜索来利用随机性。本章中的随机性思路可以帮助你在工作中降低计算成本,并生成更高分的算法。
本章最后介绍了一个堆叠聚合器,它是一个可能非常不同模型的集成方法。堆叠中的数据分析部分是将多个机器学习算法的预测作为输入。
很多数据科学任务计算量很大。如果可能的话,使用多核计算机。在整个过程中,有一个名为 n_jobs
的参数设置为 -1
,它会利用计算机的所有核心。
使用决策树进行基本分类
在这里,我们使用决策树进行基本分类。决策树用于分类时,是一系列决策,用于确定分类结果或类别结果。此外,决策树可以通过 SQL 由同一公司内的其他人员查看数据来进行检查。
准备工作
重新加载鸢尾花数据集并将数据分为训练集和测试集:
from sklearn.datasets import load_iris
iris = load_iris()
X = iris.data
y = iris.target
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y)
如何实现…
- 导入决策树分类器并在训练集上训练它:
from sklearn.tree import DecisionTreeClassifier
dtc = DecisionTreeClassifier() #Instantiate tree class
dtc.fit(X_train, y_train)
- 然后在测试集上测量准确度:
from sklearn.metrics import accuracy_score
y_pred = dtc.predict(X_test)
accuracy_score(y_test, y_pred)
0.91111111111111109
决策树看起来很准确。让我们进一步检查它。
使用 pydot 可视化决策树
如果您想生成图表,请安装 pydot
库。不幸的是,对于 Windows 来说,这个安装可能会比较复杂。如果您在安装 pydot
时遇到困难,请专注于查看图表,而不是重现它们。
如何操作…
- 在 IPython Notebook 中,执行多个导入并键入以下脚本:
import numpy as np
from sklearn import tree
from sklearn.externals.six import StringIO
import pydot
from IPython.display import Image
dot_iris = StringIO()
tree.export_graphviz(dtc, out_file = dot_iris, feature_names = iris.feature_names)
graph = pydot.graph_from_dot_data(dot_iris.getvalue())
Image(graph.create_png())
它是如何工作的…
这是通过训练生成的决策树;通过在 X_train
和 y_train
上调用 fit
方法来得到的。仔细查看它,从树的顶部开始。您在训练集中有 105 个样本。训练集被分成三组,每组 35 个:value = [35, 35, 35]。具体来说,这些是 35 个 Setosa、35 个 Versicolor 和 35 个 Virginica 花:
第一个决策是花瓣长度是否小于或等于 2.45。如果答案为真,或者是,花朵将被分类为第一类,value = [35, 0, 0]。该花被分类为 Setosa 类别。在鸢尾花数据集分类的多个示例中,这是最容易分类的一个。
否则,如果花瓣长度大于 2.45,第一个决策将导致一个较小的决策树。这个较小的决策树仅包含最后两类花:Versicolour 和 Virginica,value = [0, 35, 35]。
算法继续生成一个四层的完整树,深度为 4(注意,最上层节点不算在层数之内)。用正式语言来说,图中表现决策的三个节点被称为分裂。
还有更多…
你可能会想知道在决策树的可视化中,gini 参考是什么。Gini 指的是 gini 函数,它衡量分裂的质量,三个节点表示一个决策。当算法运行时,考虑了优化 gini 函数的几个分裂。选择产生最佳 gini 不纯度度量的分裂。
另一个选择是测量熵来确定如何分割树。您可以尝试这两种选择,并通过交叉验证确定哪种效果最好。按如下方式更改决策树中的标准:
from sklearn.tree import DecisionTreeClassifier
dtc = DecisionTreeClassifier(criterion='entropy')
dtc.fit(X_train, y_train)
这将生成如下决策树图:
您可以使用 GridSearchCV
检查此标准在交叉验证下的表现,并在参数网格中更改标准参数。我们将在下一节中进行此操作。
调整决策树
我们将继续深入探索鸢尾花数据集,重点关注前两个特征(花萼长度和花萼宽度),优化决策树并创建一些可视化图表。
准备工作
- 加载鸢尾花数据集,专注于前两个特征。并将数据分为训练集和测试集:
from sklearn.datasets import load_iris
iris = load_iris()
X = iris.data[:,:2]
y = iris.target
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y)
- 使用 pandas 查看数据:
import pandas as pd
pd.DataFrame(X,columns=iris.feature_names[:2])
- 在优化决策树之前,我们先试试一个默认参数的单一决策树。实例化并训练一个决策树:
from sklearn.tree import DecisionTreeClassifier
dtc = DecisionTreeClassifier() #Instantiate tree with default parameters
dtc.fit(X_train, y_train)
- 测量准确度分数:
from sklearn.metrics import accuracy_score
y_pred = dtc.predict(X_test)
accuracy_score(y_test, y_pred)
0.66666666666666663
使用graphviz
可视化树,揭示了一棵非常复杂的树,包含许多节点和层级(这张图片仅供参考:如果你看不懂也没关系!它是一棵非常深的树,存在很多过拟合!):
这是过拟合的一个例子。决策树非常复杂。整个鸢尾数据集由 150 个样本组成,而一个非常复杂的树是不可取的。回想一下,在之前的章节中,我们使用了线性 SVM,它通过几条直线简单地划分空间。
在继续之前,使用 matplotlib 可视化训练数据点:
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=((12,6)))
plt.xlabel(iris.feature_names[0])
plt.ylabel(iris.feature_names[1])
plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train)
如何实现…
- 为了优化决策树的性能,使用
GridSearchCV
。首先实例化一个决策树:
from sklearn.tree import DecisionTreeClassifier
dtc = DecisionTreeClassifier()
- 然后,实例化并训练
GridSearchCV
:
from sklearn.model_selection import GridSearchCV, cross_val_score
param_grid = {'criterion':['gini','entropy'], 'max_depth' : [3,5,7,20]}
gs_inst = GridSearchCV(dtc,param_grid=param_grid,cv=5)
gs_inst.fit(X_train, y_train)
请注意,在参数网格param_grid
中,我们将分割评分标准在gini
和entropy
之间变化,并调整树的max_depth
。
- 现在尝试在测试集上评分准确度:
from sklearn.metrics import accuracy_score
y_pred_gs = gs_inst.predict(X_test)
accuracy_score(y_test, y_pred_gs)
0.68888888888888888
准确度略有提高。让我们更详细地看看GridSearchCV
。
- 查看网格搜索中尝试的所有决策树的分数:
gs_inst.grid_scores_
[mean: 0.78095, std: 0.09331, params: {'criterion': 'gini', 'max_depth': 3},
mean: 0.68571, std: 0.08832, params: {'criterion': 'gini', 'max_depth': 5},
mean: 0.70476, std: 0.08193, params: {'criterion': 'gini', 'max_depth': 7},
mean: 0.66667, std: 0.09035, params: {'criterion': 'gini', 'max_depth': 20},
mean: 0.78095, std: 0.09331, params: {'criterion': 'entropy', 'max_depth': 3},
mean: 0.69524, std: 0.11508, params: {'criterion': 'entropy', 'max_depth': 5},
mean: 0.72381, std: 0.09712, params: {'criterion': 'entropy', 'max_depth': 7},
mean: 0.67619, std: 0.09712, params: {'criterion': 'entropy', 'max_depth': 20}]
请注意,这种方法将在未来版本的 scikit-learn 中不可用。你可以使用zip(gs_inst.cv_results_['mean_test_score'],gs_inst.cv_results_['params'])
来产生类似的结果。
从这个分数列表中,你可以看到较深的树表现得比浅层的树差。详细来说,训练集中的数据被分为五部分,训练发生在四部分中,而测试发生在五部分中的一部分。非常深的树会过拟合:它们在训练集上表现很好,但在交叉验证的五个测试集上表现很差。
- 使用
best_estimator_
属性选择表现最好的树:
gs_inst.best_estimator_
DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=3,
max_features=None, max_leaf_nodes=None,
min_impurity_split=1e-07, min_samples_leaf=1,
min_samples_split=2, min_weight_fraction_leaf=0.0,
presort=False, random_state=None, splitter='best')
- 使用
graphviz
可视化树:
import numpy as np
from sklearn import tree
from sklearn.externals.six import StringIO
import pydot
from IPython.display import Image
dot_iris = StringIO()
tree.export_graphviz(gs_inst.best_estimator_, out_file = dot_iris, feature_names = iris.feature_names[:2])
graph = pydot.graph_from_dot_data(dot_iris.getvalue())
Image(graph.create_png())
还有更多…
- 为了获得更多的洞察,我们将创建一个额外的可视化。首先按如下方式创建一个 NumPy 网格:
grid_interval = 0.02
x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
xmin, xmax = np.percentile(X[:, 0], [0, 100])
ymin, ymax = np.percentile(X[:, 1], [0, 100])
xmin_plot, xmax_plot = xmin - .5, xmax + .5
ymin_plot, ymax_plot = ymin - .5, ymax + .5
xx, yy = np.meshgrid(np.arange(xmin_plot, xmax_plot, grid_interval),
np.arange(ymin_plot, ymax_plot, grid_interval))
- 使用网格搜索中的
best_estimator_
属性,预测刚刚创建的 NumPy 网格上的场景:
test_preds = gs_inst.best_estimator_.predict(np.array(zip(xx.ravel(), yy.ravel())))
- 看一下可视化结果:
import matplotlib.pyplot as plt
%matplotlib inline
X_0 = X[y == 0]
X_1 = X[y == 1]
X_2 = X[y == 2]
plt.figure(figsize=(15,8)) #change figure-size for easier viewing
plt.scatter(X_0[:,0],X_0[:,1], color = 'red')
plt.scatter(X_1[:,0],X_1[:,1], color = 'blue')
plt.scatter(X_2[:,0],X_2[:,1], color = 'green')
colors = np.array(['r', 'b','g'])
plt.scatter(xx.ravel(), yy.ravel(), color=colors[test_preds], alpha=0.15)
plt.scatter(X[:, 0], X[:, 1], color=colors[y])
plt.title("Decision Tree Visualization")
plt.xlabel(iris.feature_names[0])
plt.ylabel(iris.feature_names[1])
- 使用这种类型的可视化,你可以看到决策树尝试构建矩形来分类鸢尾花的种类。每个分割都会创建一条与某个特征垂直的线。在下面的图中,有一条垂直线表示第一个决策,是否花萼长度大于(线的右边)或小于(线的左边)5.45。键入
plt.axvline(x = 5.45, color='black')
与前面的代码一起,会得到如下结果:
- 可视化前三行:
plt.axvline(x = 5.45, color='black')
plt.axvline(x = 6.2, color='black')
plt.plot((xmin_plot, 5.45), (2.8, 2.8), color='black')
水平线 sepal_width = 2.8
比较短,并且结束于 x = 5.45
,因为它不适用于 sepal_length >= 5.45 的情况。最终,多个矩形区域被创建。
- 下图展示了应用于过拟合的非常大的决策树的相同类型的可视化。决策树分类器试图将一个矩形框住鸢尾花数据集中许多特定样本,这显示了它在面对新样本时的泛化能力差:
- 最后,你还可以绘制最大深度如何影响交叉验证得分。编写一个网格搜索脚本,设置最大深度范围从 2 到 51:
from sklearn.tree import DecisionTreeClassifier
dtc = DecisionTreeClassifier()
from sklearn.model_selection import GridSearchCV, cross_val_score
max_depths = range(2,51)
param_grid = {'max_depth' : max_depths}
gs_inst = GridSearchCV(dtc, param_grid=param_grid,cv=5)
gs_inst.fit(X_train, y_train)
plt.plot(max_depths,gs_inst.cv_results_['mean_test_score'])
plt.xlabel('Max Depth')
plt.ylabel("Cross-validation Score")
该图从另一个角度展示了,较高的最大深度倾向于降低交叉验证得分。
使用决策树进行回归
回归树和分类树非常相似。开发回归模型的过程包括四个部分:
-
加载数据集
-
将数据集拆分为训练集/测试集
-
实例化一个决策树回归器并训练它
-
在测试子集上评分模型
准备工作
对于这个例子,加载 scikit-learn 的糖尿病数据集:
#Use within an Jupyter notebook
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
X = diabetes.data
y = diabetes.target
X_feature_names = ['age', 'gender', 'body mass index', 'average blood pressure','bl_0','bl_1','bl_2','bl_3','bl_4','bl_5']
现在我们已经加载了数据集,必须将数据分为训练集和测试集。在此之前,使用 pandas 可视化目标变量:
pd.Series(y).hist(bins=50)
这是一个回归示例,我们在拆分数据集时不能使用 stratify=y
。相反,我们将对目标变量进行分箱:我们将记录目标变量是否小于 50,或在 50 到 100 之间,等等。
创建宽度为 50 的区间:
bins = 50*np.arange(8)
bins
array([ 0, 50, 100, 150, 200, 250, 300, 350])
使用 np.digitize
对目标变量进行分箱:
binned_y = np.digitize(y, bins)
使用 pandas 可视化 binned_y
变量:
pd.Series(binned_y).hist(bins=50)
NumPy 数组 binned_y
记录了每个 y
元素所属的区间。现在,将数据集拆分为训练集和测试集,并对 binned_y
数组进行分层:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,stratify=binned_y)
如何操作…
- 为了创建一个决策树回归器,实例化决策树并训练它:
from sklearn.tree import DecisionTreeRegressor
dtr = DecisionTreeRegressor()
dtr.fit(X_train, y_train)
- 为了衡量模型的准确性,使用测试集对目标变量进行预测:
y_pred = dtr.predict(X_test)
- 使用误差度量比较
y_test
(真实值)和y_pred
(模型预测)。这里使用mean_absolute_error
,它是y_test
和y_pred
之间差异的绝对值的平均值:
from sklearn.metrics import mean_absolute_error
mean_absolute_error(y_test, y_pred)
58.49438202247191
- 作为替代方案,衡量平均绝对百分比误差,它是绝对值差异的平均值,差异被除以真实值元素的大小。这衡量了相对于真实值元素大小的误差幅度:
(np.abs(y_test - y_pred)/(y_test)).mean()
0.4665997687095611
因此,我们已经建立了关于糖尿病数据集的性能基准。模型的任何变化都可能影响误差度量。
还有更多…
- 使用 pandas,你可以快速可视化误差的分布。将真实值
y_test
和预测值y_pred
之间的差异转换为直方图:
pd.Series((y_test - y_pred)).hist(bins=50)
- 你也可以对百分比误差进行同样的操作:
pd.Series((y_test - y_pred)/(y_test)).hist(bins=50)
- 最后,使用前面部分的代码,查看决策树本身。注意,我们并未对最大深度进行优化:
这棵树非常复杂,很可能会发生过拟合。
使用交叉验证减少过拟合
在这里,我们将使用前面食谱中的糖尿病数据集进行交叉验证,以提高性能。首先加载数据集,如前面的食谱中所示:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
X = diabetes.data
y = diabetes.target
X_feature_names = ['age', 'gender', 'body mass index', 'average blood pressure','bl_0','bl_1','bl_2','bl_3','bl_4','bl_5']
bins = 50*np.arange(8)
binned_y = np.digitize(y, bins)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,stratify=binned_y)
如何操作…
- 使用网格搜索来减少过拟合。导入决策树并实例化它:
from sklearn.tree import DecisionTreeRegressor
dtr = DecisionTreeRegressor()
- 然后,导入
GridSearchCV
并实例化该类:
from sklearn.model_selection import GridSearchCV
gs_inst = GridSearchCV(dtr, param_grid = {'max_depth': [3,5,7,9,20]},cv=10)
gs_inst.fit(X_train, y_train)
- 查看使用
best_estimator_
属性的最佳估计器:
gs_inst.best_estimator_
DecisionTreeRegressor(criterion='mse', max_depth=3, max_features=None,
max_leaf_nodes=None, min_impurity_split=1e-07,
min_samples_leaf=1, min_samples_split=2,
min_weight_fraction_leaf=0.0, presort=False, random_state=None,
splitter='best')
- 最佳估计器的
max_depth
为3
。现在检查误差指标:
y_pred = gs_inst.predict(X_test)
from sklearn.metrics import mean_absolute_error
mean_absolute_error(y_test, y_pred)
54.299263338774338
- 检查平均百分比误差:
(np.abs(y_test - y_pred)/(y_test)).mean()
0.4672742120960478
还有更多…
最后,使用graphviz
可视化最佳回归树:
import numpy as np
from sklearn import tree
from sklearn.externals.six import StringIO
import pydot
from IPython.display import Image
dot_diabetes = StringIO()
tree.export_graphviz(gs_inst.best_estimator_, out_file = dot_diabetes, feature_names = X_feature_names)
graph = pydot.graph_from_dot_data(dot_diabetes.getvalue())
Image(graph.create_png())
该树具有更好的准确度指标,并且已经通过交叉验证以最小化过拟合。
实现随机森林回归
随机森林是一种集成算法。集成算法将多种算法结合使用以提高预测准确性。Scikit-learn 有几种集成算法,大多数都使用树来预测。让我们从扩展决策树回归开始,使用多棵决策树在随机森林中共同工作。
随机森林是由多个决策树组成的混合体,每棵树都为最终预测提供一个投票。最终的随机森林通过平均所有树的结果来计算最终输出。
准备就绪
像我们之前使用决策树一样加载糖尿病回归数据集。将所有数据分成训练集和测试集:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
X = diabetes.data
y = diabetes.target
X_feature_names = ['age', 'gender', 'body mass index', 'average blood pressure','bl_0','bl_1','bl_2','bl_3','bl_4','bl_5']
#bin target variable for better sampling
bins = 50*np.arange(8)
binned_y = np.digitize(y, bins)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,stratify=binned_y)
如何操作…
- 让我们深入了解并导入并实例化一个随机森林。训练这个随机森林:
from sklearn.ensemble import RandomForestRegressor
rft = RandomForestRegressor()
rft.fit(X_train, y_train)
- 测量预测误差。尝试在测试集上运行随机森林:
y_pred = rft.predict(X_test)
from sklearn.metrics import mean_absolute_error
mean_absolute_error(y_test, y_pred)
48.539325842696627
(np.abs(y_test - y_pred)/(y_test)).mean()
0.42821508503434541
与单棵决策树相比,误差略有下降。
- 要访问构成随机森林的任何一棵树,请使用
estimators_
属性:
rft.estimators_
[DecisionTreeRegressor(criterion='mse', max_depth=None, max_features='auto',
max_leaf_nodes=None, min_impurity_split=1e-07,
min_samples_leaf=1, min_samples_split=2,
min_weight_fraction_leaf=0.0, presort=False,
random_state=492413116, splitter='best')
...
- 要在
graphviz
中查看列表中的第一棵树,请参考列表中的第一个元素rft.estimators_[0]
:
import numpy as np
from sklearn import tree
from sklearn.externals.six import StringIO
import pydot
from IPython.display import Image
dot_diabetes = StringIO()
tree.export_graphviz(rft.estimators_[0], out_file = dot_diabetes, feature_names = X_feature_names)
graph = pydot.graph_from_dot_data(dot_diabetes.getvalue())
Image(graph.create_png())
-
要查看第二棵树,请使用
best_rft.estimators_[1]
。要查看最后一棵树,请使用best_rft.estimators_[9]
,因为默认情况下,有 10 棵树,索引从 0 到 9,这些树构成了随机森林。 -
随机森林的一个附加特性是通过
feature_importances_
属性来确定特征的重要性:
rft.feature_importances_
array([ 0.06103037, 0.00969354, 0.34865274, 0.09091215, 0.04331388,
0.04376602, 0.04827391, 0.02430837, 0.23251334, 0.09753567])
- 你也可以可视化特征的重要性:
fig, ax = plt.subplots(figsize=(10,5))
bar_rects = ax.bar(np.arange(10), rft.feature_importances_,color='r',align='center')
ax.xaxis.set_ticks(np.arange(10))
ax.set_xticklabels(X_feature_names, rotation='vertical')
最具影响力的特征是体重指数(BMI),其次是bl_4
(六项血清测量中的第四项),然后是平均血压。
最近邻的 Bagging 回归
Bagging 是一种附加的集成方法,有趣的是,它不一定涉及决策树。它在第一个训练集的随机子集上构建多个基础估计器实例。在本节中,我们尝试将k-最近邻(KNN)作为基础估计器。
实际上,bagging 估计器对于减少复杂基础估计器(例如,具有多个层次的决策树)的方差非常有效。另一方面,boosting 通过减少弱模型的偏差(例如,层次较少的决策树或线性模型)来提升性能。
为了尝试 bagging,我们将使用 scikit-learn 的随机网格搜索来寻找最佳参数,即超参数搜索。像之前一样,我们将进行以下流程:
-
确定在算法中需要优化的参数(这些是研究人员在文献中认为最值得优化的参数)。
-
创建一个参数分布,其中最重要的参数会发生变化。
-
执行随机网格搜索。如果你使用的是集成方法,开始时保持较低的估计器数量。
-
使用上一阶段的最佳参数和多个估计器。
准备工作
再次加载上一节中使用的糖尿病数据集:
import numpy as np
import pandas as pd
from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
X = diabetes.data
y = diabetes.target
X_feature_names = ['age', 'gender', 'body mass index', 'average blood pressure','bl_0','bl_1','bl_2','bl_3','bl_4','bl_5']
#bin target variable for better sampling
bins = 50*np.arange(8)
binned_y = np.digitize(y, bins)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,stratify=binned_y)
如何操作…
- 首先,导入
BaggingRegressor
和KNeighborsRegressor
。另外,还需要导入RandomizedSearchCV
:
from sklearn.ensemble import BaggingRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import RandomizedSearchCV
-
然后,为网格搜索设置一个参数分布。对于 bagging 元估计器,一些要变化的参数包括
max_samples
、max_features
、oob_score
和估计器数量n_estimators
。估计器的数量初始设置为较低的 100,以便在尝试大量估计器之前优化其他参数。 -
此外,还有一个 KNN 算法的参数列表。它被命名为
base_estimator__n_neighbors
,其中n_neighbors
是 KNN 类中的内部名称。base_estimator
是BaggingRegressor
类中基础估计器的名称。base_estimator__n_neighbors
列表包含数字3
和5
,它们表示最近邻算法中的邻居数量:
param_dist = {
'max_samples': [0.5,1.0],
'max_features' : [0.5,1.0],
'oob_score' : [True, False],
'base_estimator__n_neighbors': [3,5],
'n_estimators': [100]
}
- 实例化
KNeighboursRegressor
类,并将其作为BaggingRegressor
中的base_estimator
:
single_estimator = KNeighborsRegressor()
ensemble_estimator = BaggingRegressor(base_estimator = single_estimator)
- 最后,实例化并运行一个随机搜索。进行几次迭代,
n_iter = 5
,因为这可能会耗时较长:
pre_gs_inst_bag = RandomizedSearchCV(ensemble_estimator,
param_distributions = param_dist,
cv=3,
n_iter = 5,
n_jobs=-1)
pre_gs_inst_bag.fit(X_train, y_train)
- 查看随机搜索运行中的最佳参数:
pre_gs_inst_bag.best_params_
{'base_estimator__n_neighbors': 5,
'max_features': 1.0,
'max_samples': 0.5,
'n_estimators': 100,
'oob_score': True}
- 使用最佳参数训练
BaggingRegressor
,除了n_estimators
,你可以增加它。在这种情况下,我们将估计器的数量增加到 1,000:
rs_bag = BaggingRegressor(**{'max_features': 1.0,
'max_samples': 0.5,
'n_estimators': 1000,
'oob_score': True,
'base_estimator': KNeighborsRegressor(n_neighbors=5)})
rs_bag.fit(X_train, y_train)
- 最后,在测试集上评估算法的表现。虽然该算法的表现不如其他算法,但我们可以在后续的堆叠聚合器中使用它:
y_pred = rs_bag.predict(X_test)
from sklearn.metrics import r2_score, mean_absolute_error
print "R-squared",r2_score(y_test, y_pred)
print "MAE : ",mean_absolute_error(y_test, y_pred)
print "MAPE : ",(np.abs(y_test - y_pred)/y_test).mean()
R-squared 0.498096653258
MAE : 44.3642741573
MAPE : 0.419361955306
如果仔细观察,你会发现,在上一节中,袋装回归的表现稍微优于随机森林,因为无论是平均绝对误差还是平均绝对百分误差都更好。始终记住,你不必将集成学习仅限于树形结构——在这里,你可以使用 KNN 算法构建一个集成回归器。
调整梯度提升树
我们将使用梯度提升树来分析加利福尼亚住房数据集。我们整体的方法与之前相同:
-
关注梯度提升算法中的重要参数:
-
max_features
-
max_depth
-
min_samples_leaf
-
learning_rate
-
loss
-
-
创建一个参数分布,其中最重要的参数会有所变化。
-
执行随机网格搜索。如果使用集成方法,开始时保持估计器的数量较低。
-
使用上一阶段的最佳参数和更多估计器进行训练。
准备就绪
加载加利福尼亚住房数据集,并将加载的数据集分为训练集和测试集:
%matplotlib inline
from __future__ import division #Load within Python 2.7 for regular division
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing
cali_housing = fetch_california_housing()
X = cali_housing.data
y = cali_housing.target
#bin output variable to split training and testing sets into two similar sets
bins = np.arange(6)
binned_y = np.digitize(y, bins)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,stratify=binned_y)
如何操作…
- 加载梯度提升算法和随机网格搜索:
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import RandomizedSearchCV
- 为梯度提升树创建一个参数分布:
param_dist = {'max_features' : ['log2',1.0],
'max_depth' : [3, 5, 7, 10],
'min_samples_leaf' : [2, 3, 5, 10],
'n_estimators': [50, 100],
'learning_rate' : [0.0001,0.001,0.01,0.05,0.1,0.3],
'loss' : ['ls','huber']
}
- 运行网格搜索以找到最佳参数。进行 30 次迭代的随机搜索:
pre_gs_inst = RandomizedSearchCV(GradientBoostingRegressor(warm_start=True),
param_distributions = param_dist,
cv=3,
n_iter = 30, n_jobs=-1)
pre_gs_inst.fit(X_train, y_train)
- 现在以数据框形式查看报告。查看报告的函数已经封装,可以多次使用:
import numpy as np
import pandas as pd
def get_grid_df(fitted_gs_estimator):
res_dict = fitted_gs_estimator.cv_results_
results_df = pd.DataFrame()
for key in res_dict.keys():
results_df[key] = res_dict[key]
return results_df
def group_report(results_df):
param_cols = [x for x in results_df.columns if 'param' in x and x is not 'params']
focus_cols = param_cols + ['mean_test_score']
print "Grid CV Report \n"
output_df = pd.DataFrame(columns = ['param_type','param_set',
'mean_score','mean_std'])
cc = 0
for param in param_cols:
for key,group in results_df.groupby(param):
output_df.loc[cc] = (param, key, group['mean_test_score'].mean(), group['mean_test_score'].std())
cc += 1
return output_df
- 查看展示梯度提升树在不同参数设置下表现的数据框:
results_df = get_grid_df(pre_gs_inst)
group_report(results_df)
从这个数据框来看;ls
作为损失函数明显优于huber
,3
是最佳的min_samples_leaf
(但4
也可能表现不错),3
是最佳的max_depth
(尽管1
或2
也能有效),0.3
作为学习率效果很好(0.2
或0.4
也可以),max_features
为1.0
效果良好,但也可以是其他数字(如特征的一半:0.5
)。
- 有了这些信息,尝试另一次随机搜索:
param_dist = {'max_features' : ['sqrt',0.5,1.0],
'max_depth' : [2,3,4],
'min_samples_leaf' : [3, 4],
'n_estimators': [50, 100],
'learning_rate' : [0.2,0.25, 0.3, 0.4],
'loss' : ['ls','huber']
}
pre_gs_inst = RandomizedSearchCV(GradientBoostingRegressor(warm_start=True),
param_distributions = param_dist,
cv=3,
n_iter = 30, n_jobs=-1)
pre_gs_inst.fit(X_train, y_train)
- 查看生成的新报告:
results_df = get_grid_df(pre_gs_inst)
group_report(results_df)
- 有了这些信息,你可以使用以下参数分布再进行一次随机搜索:
param_dist = {'max_features' : [0.4, 0.5, 0.6],
'max_depth' : [5,6],
'min_samples_leaf' : [4,5],
'n_estimators': [300],
'learning_rate' : [0.3],
'loss' : ['ls','huber']
}
- 将结果存储在
rs_gbt
中,并使用 4,000 个估计器最后一次进行训练:
rs_gbt = GradientBoostingRegressor(warm_start=True,
max_features = 0.5,
min_samples_leaf = 4,
learning_rate=0.3,
max_depth = 6,
n_estimators = 4000,loss = 'huber')
rs_gbt.fit(X_train, y_train)
- 使用 scikit-learn 的
metrics
模块描述测试集上的误差:
y_pred = rs_gbt.predict(X_test)
from sklearn.metrics import r2_score, mean_absolute_error
print "R-squared",r2_score(y_test, y_pred)
print "MAE : ",mean_absolute_error(y_test, y_pred)
print "MAPE : ",(np.abs(y_test - y_pred)/y_test).mean()
R-squared 0.84490423214
MAE : 0.302125381378
MAPE : 0.169831775387
如果你还记得,随机森林的 R 平方值略低,为 0.8252。这种算法稍微好一些。对于两者,我们都进行了随机化搜索。请注意,如果你频繁地进行树的超参数优化,可以自动化多个随机化参数搜索。
还有更多内容…
现在,我们将优化一个梯度提升分类器,而不是回归器。过程非常相似。
寻找梯度提升分类器的最佳参数
使用梯度提升树进行分类与我们之前做的回归非常相似。我们将再次执行以下操作:
-
查找梯度提升分类器的最佳参数。这些参数与梯度提升回归器相同,不同之处在于损失函数选项有所不同。参数名称相同,具体如下:
-
max_features
-
max_depth
-
min_samples_leaf
-
learning_rate
-
loss
-
-
使用最佳参数运行估计器,但在估计器中使用更多的树。在下面的代码中,请注意损失函数(称为 deviance)的变化。为了进行分类,我们将使用一个二元变量。回顾目标集
y
的可视化:
pd.Series(y).hist(bins=50)
在最右侧,似乎存在异常:分布中的许多值等于五。也许我们想将该数据集分开并单独分析。作为这个过程的一部分,我们可能希望能够预先确定一个点是否应该属于异常集。我们将构建一个分类器,将 y
大于或等于五的点分离出来:
- 首先,将数据集分割为训练集和测试集。对分箱变量
binned_y
进行分层:
bins = np.arange(6)
binned_y = np.digitize(y, bins)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,stratify=binned_y)
- 创建一个二元变量,当目标变量
y
大于或等于 5 时其值为1
,小于 5 时为0
。注意,如果二元变量为1
,则表示它属于异常集:
y_binary = np.where(y >= 5, 1,0)
- 现在,使用
X_train
的形状将二元变量分割为y_train_binned
和y_test_binned
:
train_shape = X_train.shape[0]
y_train_binned = y_binary[:train_shape]
y_test_binned = y_binary[train_shape:]
- 执行随机网格搜索:
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import RandomizedSearchCV
param_dist = {'max_features' : ['log2',0.5,1.0],
'max_depth' : [2,3,6],
'min_samples_leaf' : [1,2,3,10],
'n_estimators': [100],
'learning_rate' : [0.1,0.2,0.3,1],
'loss' : ['deviance']
}
pre_gs_inst = RandomizedSearchCV(GradientBoostingClassifier(warm_start=True),
param_distributions = param_dist,
cv=3,
n_iter = 10, n_jobs=-1)
pre_gs_inst.fit(X_train, y_train_binned)
- 查看最佳参数:
pre_gs_inst.best_params_
{'learning_rate': 0.2,
'loss': 'deviance',
'max_depth': 2,
'max_features': 1.0,
'min_samples_leaf': 2,
'n_estimators': 50}
- 增加估计器数量并训练最终的估计器:
gbc = GradientBoostingClassifier(**{'learning_rate': 0.2,
'loss': 'deviance',
'max_depth': 2,
'max_features': 1.0,
'min_samples_leaf': 2,
'n_estimators': 1000, 'warm_start':True}).fit(X_train, y_train_binned)
- 查看算法的性能:
y_pred = gbc.predict(X_test)
from sklearn.metrics import accuracy_score
accuracy_score(y_test_binned, y_pred)
0.93580426356589153
该算法是一个二分类器,准确率大约为 94%,能够判断房屋是否属于异常集。梯度提升分类器的超参数优化非常相似,具有与梯度提升回归相同的重要参数。
调整 AdaBoost 回归器
在 AdaBoost 回归器中,重要的调节参数是 learning_rate
和 loss
。与之前的算法一样,我们将执行随机参数搜索,以找到该算法能达到的最佳分数。
如何操作…
- 导入算法和随机网格搜索。尝试一个随机的参数分布:
from sklearn.ensemble import AdaBoostRegressor
from sklearn.model_selection import RandomizedSearchCV
param_dist = {
'n_estimators': [50, 100],
'learning_rate' : [0.01,0.05,0.1,0.3,1],
'loss' : ['linear', 'square', 'exponential']
}
pre_gs_inst = RandomizedSearchCV(AdaBoostRegressor(),
param_distributions = param_dist,
cv=3,
n_iter = 10,
n_jobs=-1)
pre_gs_inst.fit(X_train, y_train)
- 查看最佳参数:
pre_gs_inst.best_params_
{'learning_rate': 0.05, 'loss': 'linear', 'n_estimators': 100}
- 这些表明进行另一次带有参数分布的随机搜索:
param_dist = {
'n_estimators': [100],
'learning_rate' : [0.04,0.045,0.05,0.055,0.06],
'loss' : ['linear']
}
- 复制包含最佳参数的字典。将副本中的估计器数量增加到 3,000:
import copy
ada_best = copy.deepcopy(pre_gs_inst.best_params_)
ada_best['n_estimators'] = 3000
- 训练最终的 AdaBoost 模型:
rs_ada = AdaBoostRegressor(**ada_best)
rs_ada.fit(X_train, y_train)
- 在测试集上衡量模型性能:
y_pred = rs_ada.predict(X_test)
from sklearn.metrics import r2_score, mean_absolute_error
print "R-squared",r2_score(y_test, y_pred)
print "MAE : ",mean_absolute_error(y_test, y_pred)
print "MAPE : ",(np.abs(y_test - y_pred)/y_test).mean()
R-squared 0.485619387823
MAE : 0.708716094846
MAPE : 0.524923208329
不幸的是,这个模型明显不如其他树模型表现得好。我们将暂时搁置它,不再进一步优化,因为这样做会增加更多的训练时间和 Python 开发时间。
还有更多…
我们已经找到了几个算法的最佳参数。下面是一个表格,总结了每个算法在交叉验证下需要优化的参数。建议您从优化这些参数开始:
使用 scikit-learn 编写堆叠聚合器
在本节中,我们将使用 scikit-learn 编写一个堆叠聚合器。堆叠聚合器将可能非常不同类型的模型进行混合。我们看到的许多集成算法混合的是同类型的模型,通常是决策树。
堆叠聚合器中的基本过程是,我们使用几个机器学习算法的预测作为训练另一个机器学习算法的输入。
更详细地说,我们使用一对 X
和 y
集合(X_1
,y_1
)训练两个或多个机器学习算法。然后我们在第二个 X
集(X_stack
)上做出预测,得到 y_pred_1
、y_pred_2
等。
这些预测,y_pred_1
和 y_pred_2
,成为一个机器学习算法的输入,训练输出为 y_stack
。最后,可以在第三个输入集 X_3
和真实标签集 y_3
上衡量错误。
在一个示例中会更容易理解。
如何操作…
- 再次从加利福尼亚住房数据集加载数据。观察我们如何再次创建分箱,以便对一个连续变量进行分层:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing
cali_housing = fetch_california_housing()
X = cali_housing.data
y = cali_housing.target
bins = np.arange(6)
from __future__ import division
from sklearn.model_selection import train_test_split
binned_y = np.digitize(y, bins)
from sklearn.ensemble import RandomForestRegressor, AdaBoostRegressor, ExtraTreesRegressor, GradientBoostingRegressor
from sklearn.model_selection import GridSearchCV
- 现在通过使用
train_test_split
两次,将一对X
和y
拆分为三个X
和y
对,输入和输出。注意我们在每个阶段如何对连续变量进行分层:
X_train_prin, X_test_prin, y_train_prin, y_test_prin = train_test_split(X, y,
test_size=0.2,
stratify=binned_y)
binned_y_train_prin = np.digitize(y_train_prin, bins)
X_1, X_stack, y_1, y_stack = train_test_split(X_train_prin,
y_train_prin,
test_size=0.33,
stratify=binned_y_train_prin )
- 使用
RandomizedSearchCV
查找堆叠聚合器中第一个算法的最佳参数,在此例中是多个最近邻模型的袋装算法:
from sklearn.ensemble import BaggingRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import RandomizedSearchCV
param_dist = {
'max_samples': [0.5,1.0],
'max_features' : [0.5,1.0],
'oob_score' : [True, False],
'base_estimator__n_neighbors': [3,5],
'n_estimators': [100]
}
single_estimator = KNeighborsRegressor()
ensemble_estimator = BaggingRegressor(base_estimator = single_estimator)
pre_gs_inst_bag = RandomizedSearchCV(ensemble_estimator,
param_distributions = param_dist,
cv=3,
n_iter = 5,
n_jobs=-1)
pre_gs_inst_bag.fit(X_1, y_1)
- 使用最佳参数,训练一个使用多个估算器的袋装回归器,在此例中为 3,000 个估算器:
rs_bag = BaggingRegressor(**{'max_features': 0.5,
'max_samples': 0.5,
'n_estimators': 3000,
'oob_score': False,
'base_estimator': KNeighborsRegressor(n_neighbors=3)})
rs_bag.fit(X_1, y_1)
- 对
X_1
和y_1
集合进行梯度提升算法的相同处理:
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import RandomizedSearchCV
param_dist = {'max_features' : ['log2',0.4,0.5,0.6,1.0],
'max_depth' : [2,3, 4, 5,6, 7, 10],
'min_samples_leaf' : [1,2, 3, 4, 5, 10],
'n_estimators': [50, 100],
'learning_rate' : [0.01,0.05,0.1,0.25,0.275,0.3,0.325],
'loss' : ['ls','huber']
}
pre_gs_inst = RandomizedSearchCV(GradientBoostingRegressor(warm_start=True),
param_distributions = param_dist,
cv=3,
n_iter = 30, n_jobs=-1)
pre_gs_inst.fit(X_1, y_1)
- 使用更多估算器训练最佳参数集:
gbt_inst = GradientBoostingRegressor(**{'learning_rate': 0.05,
'loss': 'huber',
'max_depth': 10,
'max_features': 0.4,
'min_samples_leaf': 5,
'n_estimators': 3000,
'warm_start': True}).fit(X_1, y_1)
- 使用
X_stack
和两个算法预测目标:
y_pred_bag = rs_bag.predict(X_stack)
y_pred_gbt = gbt_inst.predict(X_stack)
- 查看每个算法产生的指标(错误率)。查看袋装回归器的指标:
from sklearn.metrics import r2_score, mean_absolute_error
print "R-squared",r2_score(y_stack, y_pred_bag)
print "MAE : ",mean_absolute_error(y_stack, y_pred_bag)
print "MAPE : ",(np.abs(y_stack- y_pred_bag)/y_stack).mean()
R-squared 0.527045729567
MAE : 0.605868386902
MAPE : 0.397345752723
- 查看梯度提升的指标:
from sklearn.metrics import r2_score, mean_absolute_error
print "R-squared",r2_score(y_stack, y_pred_gbt)
print "MAE : ",mean_absolute_error(y_stack, y_pred_gbt)
print "MAPE : ",(np.abs(y_stack - y_pred_gbt)/y_stack).mean()
R-squared 0.841011059404
MAE : 0.297099247278
MAPE : 0.163956322255
- 创建一个包含两个算法预测的 DataFrame。或者,你也可以创建一个 NumPy 数组来存储数据:
y_pred_bag = rs_bag.predict(X_stack)
y_pred_gbt = gbt_inst.predict(X_stack)
preds_df = pd.DataFrame(columns = ['bag', 'gbt'])
preds_df['bag'] = y_pred_bag
preds_df['gbt'] = y_pred_gbt
- 查看新的预测数据框:
preds_df
- 查看预测列之间的相关性。这些列是相关的,但不是完全相关。理想的情况是算法之间没有完全相关,并且两个算法都表现良好。在这种情况下,袋装回归器的表现远不如梯度提升:
preds_df.corr()
- 现在使用第三个算法进行随机搜索。这个算法将前两个算法的预测作为输入。我们将使用额外的树回归器对其他两个算法的预测进行预测:
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.model_selection import RandomizedSearchCV
param_dist = {'max_features' : ['sqrt','log2',1.0],
'min_samples_leaf' : [1, 2, 3, 7, 11],
'n_estimators': [50, 100],
'oob_score': [True, False]}
pre_gs_inst = RandomizedSearchCV(ExtraTreesRegressor(warm_start=True,bootstrap=True),
param_distributions = param_dist,
cv=3,
n_iter = 15)
pre_gs_inst.fit(preds_df.values, y_stack)
- 复制参数字典,并在复制的字典中增加估算器的数量。如果你想查看最终的字典,可以查看:
import copy
param_dict = copy.deepcopy(pre_gs_inst.best_params_)
param_dict['n_estimators'] = 2000
param_dict['warm_start'] = True
param_dict['bootstrap'] = True
param_dict['n_jobs'] = -1
param_dict
{'bootstrap': True,
'max_features': 1.0,
'min_samples_leaf': 11,
'n_estimators': 2000,
'n_jobs': -1,
'oob_score': False,
'warm_start': True}
- 在预测数据框上训练额外的树回归器,使用
y_stack
作为目标:
final_etr = ExtraTreesRegressor(**param_dict)
final_etr.fit(preds_df.values, y_stack)
- 为了检查堆叠聚合器的整体性能,你需要一个函数,该函数以
X
集合为输入,通过袋装回归器和梯度提升创建一个数据框,并最终在这些预测上进行预测:
def handle_X_set(X_train_set):
y_pred_bag = rs_bag.predict(X_train_set)
y_pred_gbt = gbt_inst.predict(X_train_set)
preds_df = pd.DataFrame(columns = ['bag', 'gbt'])
preds_df['bag'] = y_pred_bag
preds_df['gbt'] = y_pred_gbt
return preds_df.values
def predict_from_X_set(X_train_set):
return final_etr.predict(handle_X_set(X_train_set))
- 使用
X_test_prin
进行预测,这是我们刚刚创建的有用的predict_from_X_set
函数,在留出的X
集合上进行预测:
y_pred = predict_from_X_set(X_test_prin)
- 测量模型的性能:
from sklearn.metrics import r2_score, mean_absolute_error
print "R-squared",r2_score(y_test_prin, y_pred)
print "MAE : ",mean_absolute_error(y_test_prin, y_pred)
print "MAPE : ",(np.abs(y_test_prin- y_pred)/y_test_prin).mean()
R-squared 0.844114615094
MAE : 0.298422222752
MAPE : 0.173901911714
接下来怎么办?R 平方指标略有改善,我们为这一小小的提升付出了很多努力。接下来我们可以编写更健壮、更像生产环境的代码,为堆叠器添加许多不相关的估计器,使其更易于操作。
此外,我们可以进行特征工程——通过数学和/或加利福尼亚房地产行业的领域知识来改进数据列。你还可以尝试对不同的输入使用不同的算法。两个列:纬度和经度,非常适合随机森林,而其他输入可以通过线性算法来很好地建模。
第三,我们可以在数据集上探索不同的算法。对于这个数据集,我们专注于复杂的、高方差的算法。我们也可以尝试一些简单的高偏差算法。这些替代算法可能有助于我们最终使用的堆叠聚合器。
最后,关于堆叠器,你可以通过交叉验证旋转X_stacker
集合,以最大化训练集的利用率。