集成学习和随机森林方法

集成

假设已经为某一特定问题选中了最佳的模型，想进一步提升其准确率，就需要应用一些更高级的机器学习技术：集成（Ensemble）。集成是使用一系列学习器进行学习，并使用某种规则把各个学习结果进行整合从而获得比单个学习器更好的学习效果的一种机器学习方法。在集成中，最终的整体输出比任何单个部分的表现更重要。
某种意义上，孔多塞陪审团定理形象的描述了上面提到的集成概念。该定理的内容为：如果评审团的每个成员做出独立判断，并且每个陪审员做出正确决策的概率高于 0.5，那么整个评审团做出正确的总体决策的概率随着陪审员数量的增加而增加，并趋向于一。另一方面，如果每个陪审员判断正确的概率小于 0.5，那么整个陪审团做出正确的总体决策的概率随着陪审员数量的增加而减少，并趋向于零。
该定理的公式为：

image.png

让我们看看另一个集成的例子：群体的智慧。1906 年，Francis Galton 访问了普利茅斯的一个农村集市，在那里他看到一个比赛。800 个参与者尝试估计一头屠宰的牛的重量。所有参与者的预测的平均值为 1197 磅，与牛的真实重量 1198 磅十分接近。

Bootstrapping

Bagging 基于统计学中的 Bootstraping（自助法），该方法令复杂模型的统计评估变得更加可行。
Bootstrap 方法的流程如下：假设有尺寸为 N 的样本 X，从该样本中有放回地随机抽取 N 个样本，以创建一个新样本。换句话说，从尺寸为 N 的原样本中随机选择一个元素，并重复此过程 N 次。选中所有元素的可能性是一样的，因此每个元素被抽中的概率均为 1/N。
假设用 Bootstrap 方法从一个袋子中抽球，每次抽一个。在每一步中，将选中的球放回袋子，这样下一次抽取是等概率的，即，从同样数量的 N 个球中抽取。注意，因为我们把球放回了，新样本中可能有重复的球。把这个新样本称为X 1。重复这一过程 M 次，创建 M 个 Bootstrap 样本X1,…,XM。最后，我们的样本数量就从原先的 1 个扩充到了 M 个，就有充足的样本来计算原始分布的多种统计数据。

image.png

看一个例子，这个例子使用之前的 telecom_churn 数据集。我们曾讨论过这一数据集的特征重要性，其中最重要的特征之一是呼叫客服次数（Customer service calls）。可视化「呼叫客服次数」这一特征，看看其分布。

import warnings
import seaborn as sns
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 10, 6
%matplotlib inline
warnings.filterwarnings('ignore')

telecom_data = pd.read_csv(
    'https://labfile.oss.aliyuncs.com/courses/1283/telecom_churn.csv')

fig = sns.kdeplot(telecom_data[telecom_data['Churn']
                               == False]['Customer service calls'], label='Loyal')
fig = sns.kdeplot(telecom_data[telecom_data['Churn']
                               == True]['Customer service calls'], label='Churn')
fig.set(xlabel='Number of calls', ylabel='Density')
plt.show()

上图表明，相比那些逐渐离网的客户，忠实客户呼叫客服的次数更少。估计每组客户的平均呼叫客服数可能很有用，但由于该数据集很小，如果直接计算原样本的均值，得到的估计可能较差。因此可以应用 Bootstrap 方法，基于原样本生成 1000 个新 Bootstrap 样本，然后计算均值的区间估计。
定义 get_bootstrap_samples() 函数产生 Bootstrap 样本，定义 stat_intervals() 函数产生一个区间估计。

def get_bootstrap_samples(data, n_samples):
    indices = np.random.randint(0, len(data), (n_samples, len(data)))
    samples = data[indices]
    return samples


def stat_intervals(stat, alpha):
    boundaries = np.percentile(
        stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
    return boundaries

分割数据集，分组为忠实客户和离网客户。

loyal_calls = telecom_data.loc[telecom_data['Churn'] == False,
                               'Customer service calls'].values
churn_calls = telecom_data.loc[telecom_data['Churn'] == True,
                               'Customer service calls'].values

固定随机数种子，以得到可重现的结果。

np.random.seed(0)

使用 Bootstrap 生成样本，计算各自的均值。

loyal_mean_scores = [np.mean(sample)
                     for sample in get_bootstrap_samples(loyal_calls, 1000)]
churn_mean_scores = [np.mean(sample)
                     for sample in get_bootstrap_samples(churn_calls, 1000)]

打印区间估计值，使用 stat_intervals() 函数将区间定义为 95%。

print("Service calls from loyal: mean interval",
      stat_intervals(loyal_mean_scores, 0.05))
print("Service calls from churn: mean interval",
      stat_intervals(churn_mean_scores, 0.05))

上述结果表明，有 95% 的概率，忠实客户平均呼叫客服的次数在 1.4 到 1.49 之间，而离网客户平均呼叫客服的次数在 2.06 到 2.40 之间。另外，可以注意到忠实客户的区间更窄，这是合理的，因为相比多次呼叫客服的离网客户，忠实客户呼叫客服的次数更少（0、1、2 次）。

Bagging

假设我们有一个训练集 X。我们使用 Bootstrap 生成样本X1,…,XM。现在，我们在每个Bootstrap样本上分别训练分类器ai(x)，最终的分类器将对所有这些单独的分类器的输出取均值。在分类情形下，这种技术即投票（voting）：

image.png

下图形象化解释了上面的公式：

image.png

image.png

Bagging 通过在不同数据集上训练模型的方式来降低分类器的方差。换句话说，Bagging 可以预防过拟合。由于来自不同训练数据集上的模型各不相同，它们的误差在投票过程中会相互抵消，所以 Bagging 是有效的。此外，某些 Bootstrap 训练样本很可能会忽略异常值。
Bagging 的实际效果，并与决策树比较下，如下图所示：

image.png

image.png

袋外误差（OOB Error）

随机森林不需要使用交叉验证或留置样本，因为它内置了误差估计。随机森林中的决策树基于原始数据集中不同的 Bootstrap 样本构建。对第 K 棵树而言，其特定 Bootstrap 样本大约留置了 37% 的输入。
设数据集中有ℓ个样本。在每一步，每个数据点最终出现在有放回的 Bootstrap 样本中的概率均为1/ℓ。Bootstrap 样本最终不包含特定数据集元素的概率（即该元素在 ℓ 次抽取中都没抽中）等于 (1 - 1/ℓ)^ℓ。当ℓ→+∞ 时，这一概率等于e/1。因此，选中某一特定样本的概率≈1− 1/e≈0.63。

image.png

示意图上方为原始数据集。将其分为训练集（左图）和测试集（右图）。在测试集上，绘制一副网格，完美地实现了分类。现在，我们应用同一副网格到测试集上，以估计分类的正确率。右图可见，分类器在 4 个未曾在训练集中出现的数据点上给出了错误的答案。而测试集中共有 15 个数据点，这 15 个数据点未在训练集中出现。因此，该分类器的准确率为11/15*100% = 73.33% 。
总结一下，每个基础算法在约 63% 的原始样本上训练。该算法可以在剩下的约 37% 的样本上验证。袋外估计仅仅是基础算法在 37% 样本上的平均估计。

随机森林

在 Bagging 中，主要使用的基础分类器为决策树，因为决策树相当复杂，并能在任何样本上达到零分类误差。随机子空间方法可以降低树的相关性，从而避免过拟合。基于 Bagging，基础算法会在不同的原始特征集的随机子集上训练。

image.png

image.png

from sklearn.model_selection import train_test_split
from sklearn.datasets import make_circles
from sklearn.tree import DecisionTreeRegressor, DecisionTreeClassifier
from sklearn.ensemble import BaggingClassifier, BaggingRegressor
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
import seaborn as sns
from matplotlib import pyplot as plt

# Disable warnings in Anaconda
import warnings
import numpy as np
warnings.filterwarnings('ignore')

%matplotlib inline
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 10, 6

n_train = 150
n_test = 1000
noise = 0.1

生成数据。

def f(x):
    x = x.ravel()
    return np.exp(-x ** 2) + 1.5 * np.exp(-(x - 2) ** 2)


def generate(n_samples, noise):
    X = np.random.rand(n_samples) * 10 - 5
    X = np.sort(X).ravel()
    y = np.exp(-X ** 2) + 1.5 * np.exp(-(X - 2) ** 2)\
        + np.random.normal(0.0, noise, n_samples)
    X = X.reshape((n_samples, 1))

    return X, y


X_train, y_train = generate(n_samples=n_train, noise=noise)
X_test, y_test = generate(n_samples=n_test, noise=noise)

使用单棵决策树回归。

dtree = DecisionTreeRegressor().fit(X_train, y_train)
d_predict = dtree.predict(X_test)

plt.figure(figsize=(10, 6))
plt.plot(X_test, f(X_test), "b")
plt.scatter(X_train, y_train, c="b", s=20)
plt.plot(X_test, d_predict, "g", lw=2)
plt.xlim([-5, 5])
plt.title("Decision tree, MSE = %.2f"
          % np.sum((y_test - d_predict) ** 2))

使用 Bagging 决策树回归。

bdt = BaggingRegressor(DecisionTreeRegressor()).fit(X_train, y_train)
bdt_predict = bdt.predict(X_test)

plt.figure(figsize=(10, 6))
plt.plot(X_test, f(X_test), "b")
plt.scatter(X_train, y_train, c="b", s=20)
plt.plot(X_test, bdt_predict, "y", lw=2)
plt.xlim([-5, 5])
plt.title("Bagging for decision trees, MSE = %.2f" %
          np.sum((y_test - bdt_predict) ** 2))

使用随机森林。

rf = RandomForestRegressor(n_estimators=10).fit(X_train, y_train)
rf_predict = rf.predict(X_test)

plt.figure(figsize=(10, 6))
plt.plot(X_test, f(X_test), "b")
plt.scatter(X_train, y_train, c="b", s=20)
plt.plot(X_test, rf_predict, "r", lw=2)
plt.xlim([-5, 5])
plt.title("Random forest, MSE = %.2f" % np.sum((y_test - rf_predict) ** 2))

从上面 3 张图像和 MSE 值可以看到，10 树随机森林比单棵决策树和 10 树 Bagging 的表现要好，当然随机森林在这个例子上不稳定，多次运行结果表明随机森林和 Bagging 决策树互有胜负。随机森林和 Bagging 的主要差别在于，前者分割的最佳特征是从一个随机特征子空间中选取的，而后者在分割时将考虑所有特征。
接下来，我们将查看随机森林和 Bagging 在分类问题上的表现。

np.random.seed(42)
X, y = make_circles(n_samples=500, factor=0.1, noise=0.35, random_state=42)
X_train_circles, X_test_circles, y_train_circles, y_test_circles = train_test_split(
    X, y, test_size=0.2)

dtree = DecisionTreeClassifier(random_state=42)
dtree.fit(X_train_circles, y_train_circles)

x_range = np.linspace(X.min(), X.max(), 100)
xx1, xx2 = np.meshgrid(x_range, x_range)
y_hat = dtree.predict(np.c_[xx1.ravel(), xx2.ravel()])
y_hat = y_hat.reshape(xx1.shape)
plt.contourf(xx1, xx2, y_hat, alpha=0.2)
plt.scatter(X[:, 0], X[:, 1], c=y, cmap='autumn')
plt.title("Decision tree")
plt.show()

b_dtree = BaggingClassifier(
    DecisionTreeClassifier(), n_estimators=300, random_state=42)
b_dtree.fit(X_train_circles, y_train_circles)

x_range = np.linspace(X.min(), X.max(), 100)
xx1, xx2 = np.meshgrid(x_range, x_range)
y_hat = b_dtree.predict(np.c_[xx1.ravel(), xx2.ravel()])
y_hat = y_hat.reshape(xx1.shape)
plt.contourf(xx1, xx2, y_hat, alpha=0.2)
plt.scatter(X[:, 0], X[:, 1], c=y, cmap='autumn')
plt.title("Bagging (decision trees)")
plt.show()

rf = RandomForestClassifier(n_estimators=300, random_state=42)
rf.fit(X_train_circles, y_train_circles)

x_range = np.linspace(X.min(), X.max(), 100)
xx1, xx2 = np.meshgrid(x_range, x_range)
y_hat = rf.predict(np.c_[xx1.ravel(), xx2.ravel()])
y_hat = y_hat.reshape(xx1.shape)
plt.contourf(xx1, xx2, y_hat, alpha=0.2)
plt.scatter(X[:, 0], X[:, 1], c=y, cmap='autumn')
plt.title("Random forest")
plt.show()

上图显示了决策树形成的边界凹凸不平，有大量锐角，这暗示出现过拟合现象。相反，随机森林和 Bagging 的边界相当平滑，没有明显的过拟合的迹象。
随机森林的参数
scikit-learn 库提供了 BaggingRegressor 和 BaggingClassifier 类，其中在创建新模型时需要注意的一些参数如下：
n_estimators 是随机森林中树的数量；
criterion 是衡量分割质量的函数；
max_features 是查找最佳分割时考虑的特征数；
min_samples_leaf 是叶节点的最小样本数；
max_depth 是树的最大深度。
随机森林在现实问题中的应用
使用之前运营商的客户离网数据集作为例子。这是一个分类问题，使用准确率评估模型。
首先，创建一个简单的分类器作为基线。出于简单的目的，只使用数值特征构建分类器。

import pandas as pd
from sklearn.model_selection import cross_val_score, StratifiedKFold, GridSearchCV
from sklearn.metrics import accuracy_score

加载数据集。

df = pd.read_csv(
    "https://labfile.oss.aliyuncs.com/courses/1283/telecom_churn.csv")

选择数值特性。

cols = []
for i in df.columns:
    if (df[i].dtype == "float64") or (df[i].dtype == 'int64'):
        cols.append(i)

将数据集分离为输入变量和目标变量。

X, y = df[cols].copy(), np.asarray(df["Churn"], dtype='int8')

为验证过程进行分层分割。

skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

基于默认参数初始化分类器。

rfc = RandomForestClassifier(random_state=42, n_jobs=-1, oob_score=True)

在训练集上进行训练。

results = cross_val_score(rfc, X, y, cv=skf)

在测试集上评估准确率。

print("CV accuracy score: {:.2f}%".format(results.mean()*100))

得到了 91.48% 的准确率作为基线，现在，尝试改进结果，同时查看一下基本参数变化时的学习曲线。首先，从改进树的数量开始，初始化 K 折交叉验证。

skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

创建列表储存训练集和测试集上的准确率。

train_acc = []
test_acc = []
temp_train_acc = []
temp_test_acc = []

进行网格搜索。下方列表中候选数值较多，会增加线上环境执行时间，可以自行减少数量以加快运行。

trees_grid = [5, 10, 15, 20, 30, 50, 75, 100]

在训练集上训练模型。

for ntrees in trees_grid:
    rfc = RandomForestClassifier(
        n_estimators=ntrees, random_state=42, n_jobs=-1, oob_score=True)
    temp_train_acc = []
    temp_test_acc = []
    for train_index, test_index in skf.split(X, y):
        X_train, X_test = X.iloc[train_index], X.iloc[test_index]
        y_train, y_test = y[train_index], y[test_index]
        rfc.fit(X_train, y_train)
        temp_train_acc.append(rfc.score(X_train, y_train))
        temp_test_acc.append(rfc.score(X_test, y_test))
    train_acc.append(temp_train_acc)
    test_acc.append(temp_test_acc)

打印结果。

train_acc, test_acc = np.asarray(train_acc), np.asarray(test_acc)
print("Best accuracy on CV is {:.2f}% with {} trees".format(max(test_acc.mean(axis=1))*100,
                                                            trees_grid[np.argmax(test_acc.mean(axis=1))]))

更改树的数量为 50 时，得到了 92.44% 的准确率。接下来，绘制相应的学习曲线。

plt.style.use('ggplot')

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(trees_grid, train_acc.mean(axis=1),
        alpha=0.5, color='blue', label='train')
ax.plot(trees_grid, test_acc.mean(axis=1), alpha=0.5, color='red', label='cv')
ax.fill_between(trees_grid, test_acc.mean(axis=1) - test_acc.std(axis=1),
                test_acc.mean(axis=1) + test_acc.std(axis=1), color='#888888', alpha=0.4)
ax.fill_between(trees_grid, test_acc.mean(axis=1) - 2*test_acc.std(axis=1),
                test_acc.mean(axis=1) + 2*test_acc.std(axis=1), color='#888888', alpha=0.2)
ax.legend(loc='best')
ax.set_ylim([0.88, 1.02])
ax.set_ylabel("Accuracy")
ax.set_xlabel("N_estimators")

上图表明，当树的数量增加到一定值时，测试集上的准确率不再提升，同时训练集上的准确率达到了 100%，这意味着出现过拟合。
为了避免过拟合，需要给模型加上正则化参数。下面将树的数目固定为 100，然后看看不同的 max_depth 效果如何，首先创建列表以保存训练和测试集上的准确率。

train_acc = []
test_acc = []
temp_train_acc = []
temp_test_acc = []
max_depth_grid = [3, 5, 7, 9, 11, 13, 15, 17, 20, 22, 24]

在训练集上训练模型。

for max_depth in max_depth_grid:
    rfc = RandomForestClassifier(
        n_estimators=100, random_state=42, n_jobs=-1, oob_score=True, max_depth=max_depth)
    temp_train_acc = []
    temp_test_acc = []
    for train_index, test_index in skf.split(X, y):
        X_train, X_test = X.iloc[train_index], X.iloc[test_index]
        y_train, y_test = y[train_index], y[test_index]
        rfc.fit(X_train, y_train)
        temp_train_acc.append(rfc.score(X_train, y_train))
        temp_test_acc.append(rfc.score(X_test, y_test))
    train_acc.append(temp_train_acc)
    test_acc.append(temp_test_acc)

打印结果。

train_acc, test_acc = np.asarray(train_acc), np.asarray(test_acc)
print("Best accuracy on CV is {:.2f}% with {} max_depth".format(max(test_acc.mean(axis=1))*100,
                                                                max_depth_grid[np.argmax(test_acc.mean(axis=1))]))

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(max_depth_grid, train_acc.mean(axis=1),
        alpha=0.5, color='blue', label='train')
ax.plot(max_depth_grid, test_acc.mean(axis=1),
        alpha=0.5, color='red', label='cv')
ax.fill_between(max_depth_grid, test_acc.mean(axis=1) - test_acc.std(axis=1),
                test_acc.mean(axis=1) + test_acc.std(axis=1), color='#888888', alpha=0.4)
ax.fill_between(max_depth_grid, test_acc.mean(axis=1) - 2*test_acc.std(axis=1),
                test_acc.mean(axis=1) + 2*test_acc.std(axis=1), color='#888888', alpha=0.2)
ax.legend(loc='best')
ax.set_ylim([0.88, 1.02])
ax.set_ylabel("Accuracy")
ax.set_xlabel("Max_depth")

max_depth 在模型中起到了正则化的作用，过拟合现象有所缓解，模型准确率提高了一些。
另一个值得调整的重要参数是 min_samples_leaf，它也能起到正则化作用。首先，创建列表以保存训练和测试集上的准确率。

train_acc = []
test_acc = []
temp_train_acc = []
temp_test_acc = []
min_samples_leaf_grid = [1, 3, 5, 7, 9, 11, 13, 15, 17, 20, 22, 24]

在训练集上训练模型。

for min_samples_leaf in min_samples_leaf_grid:
    rfc = RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1,
                                 oob_score=True, min_samples_leaf=min_samples_leaf)
    temp_train_acc = []
    temp_test_acc = []
    for train_index, test_index in skf.split(X, y):
        X_train, X_test = X.iloc[train_index], X.iloc[test_index]
        y_train, y_test = y[train_index], y[test_index]
        rfc.fit(X_train, y_train)
        temp_train_acc.append(rfc.score(X_train, y_train))
        temp_test_acc.append(rfc.score(X_test, y_test))
    train_acc.append(temp_train_acc)
    test_acc.append(temp_test_acc)

打印结果。

train_acc, test_acc = np.asarray(train_acc), np.asarray(test_acc)
print("Best accuracy on CV is {:.2f}% with {} min_samples_leaf".format(max(test_acc.mean(axis=1))*100,
                                                                       min_samples_leaf_grid[np.argmax(test_acc.mean(axis=1))]))

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(min_samples_leaf_grid, train_acc.mean(
    axis=1), alpha=0.5, color='blue', label='train')
ax.plot(min_samples_leaf_grid, test_acc.mean(
    axis=1), alpha=0.5, color='red', label='cv')
ax.fill_between(min_samples_leaf_grid, test_acc.mean(axis=1) - test_acc.std(axis=1),
                test_acc.mean(axis=1) + test_acc.std(axis=1), color='#888888', alpha=0.4)
ax.fill_between(min_samples_leaf_grid, test_acc.mean(axis=1) - 2*test_acc.std(axis=1),
                test_acc.mean(axis=1) + 2*test_acc.std(axis=1), color='#888888', alpha=0.2)
ax.legend(loc='best')
ax.set_ylim([0.88, 1.02])
ax.set_ylabel("Accuracy")
ax.set_xlabel("Min_samples_leaf")

在这一情形下，验证集上的准确率没有提高，但当验证集的准确率保持 92%以上时，过拟合的情况降低了 2%。
下面考虑 max_features 这一参数。在分类问题中，默认特征为「所有特征数的平方根」。首先，检查在这个例子中选择 4 个特征是否是最佳的，创建列表以保存训练和测试集的准确率。

train_acc = []
test_acc = []
temp_train_acc = []
temp_test_acc = []
max_features_grid = [2, 4, 6, 8, 10, 12, 14, 16]

在训练集上训练模型。

for max_features in max_features_grid:
    rfc = RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1,
                                 oob_score=True, max_features=max_features)
    temp_train_acc = []
    temp_test_acc = []
    for train_index, test_index in skf.split(X, y):
        X_train, X_test = X.iloc[train_index], X.iloc[test_index]
        y_train, y_test = y[train_index], y[test_index]
        rfc.fit(X_train, y_train)
        temp_train_acc.append(rfc.score(X_train, y_train))
        temp_test_acc.append(rfc.score(X_test, y_test))
    train_acc.append(temp_train_acc)
    test_acc.append(temp_test_acc)

打印结果。

train_acc, test_acc = np.asarray(train_acc), np.asarray(test_acc)
print("Best accuracy on CV is {:.2f}% with {} max_features".format(max(test_acc.mean(axis=1))*100,
                                                                   max_features_grid[np.argmax(test_acc.mean(axis=1))]))

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(max_features_grid, train_acc.mean(axis=1),
        alpha=0.5, color='blue', label='train')
ax.plot(max_features_grid, test_acc.mean(axis=1),
        alpha=0.5, color='red', label='cv')
ax.fill_between(max_features_grid, test_acc.mean(axis=1) - test_acc.std(axis=1),
                test_acc.mean(axis=1) + test_acc.std(axis=1), color='#888888', alpha=0.4)
ax.fill_between(max_features_grid, test_acc.mean(axis=1) - 2*test_acc.std(axis=1),
                test_acc.mean(axis=1) + 2*test_acc.std(axis=1), color='#888888', alpha=0.2)
ax.legend(loc='best')
ax.set_ylim([0.88, 1.02])
ax.set_ylabel("Accuracy")
ax.set_xlabel("Max_features")

上图表明，在这个例子中，max_features=12 时模型表现最好。
我们已经看到学习曲线如何随着基本参数的变化而变化。现在，直接使用 GridSearch() 来查找最佳参数。首先，初始化参数集以进行穷举搜索和拟合。

# 原实验搜索参数
# parameters = {'max_features': [4, 7, 10, 13], 'min_samples_leaf': [
#     1, 3, 5, 7], 'max_depth': [5, 10, 15, 20]}
# 为加速线上环境执行，优化搜索参数
parameters = {'max_features': [10, 13], 'min_samples_leaf': [1, 3], 'max_depth': [5, 10]}
rfc = RandomForestClassifier(n_estimators=100, random_state=42,
                             n_jobs=-1, oob_score=True)
gcv = GridSearchCV(rfc, parameters, n_jobs=-1, cv=skf, verbose=1)
gcv.fit(X, y)

打印结果。

gcv.best_estimator_, gcv.best_score_

调参之后，该模型的准确率达到了 92.7%。随机森林最重要的一点是它的准确率不会随着树的增加而下降，这意味着可以使用 10 棵树调整超参数，接着把树的数量从 10 增至 500 来提升准确率。
随机森林的方差和去相关

image.png

image.png

image.png

T1(x),…,Tn(x) 特征是非线性空间分割的结果，它们提供了关于样本之间的相似性的信息。《The Elements of Statistical Learning》一书中有一个例子，演示了随机森林和 k-近邻技术的相似性，如下图所示：

image.png

随机森林的优势和劣势
优势：
高预测准确率。在大多数问题上表现优于线性算法，准确率与 Boosting 相当；
随机取样导致其对离散值的鲁棒性较好；
随机选取子空间导致其对特征缩放及其他单调转换不敏感；
不需要精细的参数调整。
在具有大量特征和类别的数据集上很高效；
既可处理连续值，也可处理离散值；
不容易出现过拟合。在实践中，增加树的数量几乎总是能提升总体表现。不过，当树的数量增大到一定值后，学习曲线将趋于平稳；
有理论方法可以估计特征的重要性；
能够很好地处理数据缺失。即使有大部分的数据缺失现象，仍能保持较好的准确率；
支持整个数据集及单棵树样本上的加权分类；
决策树底层使用的相似性计算可以用于后续的聚类、离散值检测或兴趣数据表示；
以上功能和性质可以扩展到未标注数据，以支持无监督聚类，数据可视化和离散值检测；
易于并行化，伸缩性强。

劣势：
相比单棵决策树，随机森林的输出更难解释；
特征重要性估计没有形式化的 p 值；
在数据稀疏时（比如，文本输入、词袋），表现不如线性模型好；
和线性回归不同，随机森林无法外推。不过，这样也就不会因为离散值的存在导致极端值的出现；
在某些问题上容易过拟合，特别是处理高噪声数据时；
处理数量级不同的类别数据时，随机森林偏重数量级较高的数据，因为这能更明显的提高准确率；
所得模型较大，需要大容量的 RAM 来支撑。

特征的重要性

我们常常需要对算法的输出结果做出解释，在不能完全理解算法的情况下，至少希望能找出哪个输入特征对结果的贡献最大。基于随机森林，我们可以相当容易地获取这类信息。
方法精要
下图很直观地表明，在信用评分问题中，年龄比收入更重要。基于信息增益这一概念，我们可以形式化地解释这一点。

image.png

image.png

from sklearn.ensemble.forest import RandomForestRegressor
import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib import pyplot as plt

%matplotlib inline

导入数据集。

hostel_data = pd.read_csv(
    "https://labfile.oss.aliyuncs.com/courses/1283/hostel_factors.csv")
features = {"f1": u"Staff",
            "f2": u"Hostel booking",
            "f3": u"Check-in and check-out",
            "f4": u"Room condition",
            "f5": u"Shared kitchen condition",
            "f6": u"Shared space condition",
            "f7": u"Extra services",
            "f8": u"General conditions & conveniences",
            "f9": u"Value for money",
            "f10": u"Customer Co-creation"}

使用随机森林训练模型。

forest = RandomForestRegressor(n_estimators=1000, max_features=10,
                               random_state=0)

forest.fit(hostel_data.drop(['hostel', 'rating'], axis=1),
           hostel_data['rating'])
importances = forest.feature_importances_

indices = np.argsort(importances)[::-1]

标出随机森林模型中的特征重要性。

num_to_plot = 10
feature_indices = [ind+1 for ind in indices[:num_to_plot]]

打印特征重要性排名（重要性从高到低排列）。

print("Feature ranking:")

for f in range(num_to_plot):
    print("%d. %s %f " % (f + 1,
                          features["f"+str(feature_indices[f])],
                          importances[indices[f]]))
plt.figure(figsize=(15, 5))
plt.title(u"Feature Importance")
bars = plt.bar(range(num_to_plot),
               importances[indices[:num_to_plot]],
               color=([str(i/float(num_to_plot+1))
                       for i in range(num_to_plot)]),
               align="center")
ticks = plt.xticks(range(num_to_plot),
                   feature_indices)
plt.xlim([-1, num_to_plot])
plt.legend(bars, [u''.join(features["f"+str(i)])
                  for i in feature_indices])

上图表明，消费者常常更关心服务人员素质（Staff）和性价比（Value for money)，这两个因子对最终评分的影响最大。然而，这两项特征和其他特征的差别不是非常大，因此排除任何特征都会导致模型准确率的下降。
梯度提升介绍
在机器学习中，有一些算法通过常被称为弱模型，这里所说的“弱模型”指的是像决策树这样简单的基本模型，也指的是那些精度相对较差的模型，也就是比随机模型好那么一点点的模型。
先来思考一个问题：是否可以从大量相对较弱和简单的模型中通过某种手段而得到一个强大的模型?
这个问题的答案显然是肯定的。集成方法就是其中的一类，集成学习的算法有很多，先从最简单的开始：AdaBoost 。
AdaBoost 采用的是一种贪婪的学习方法。在算法学习时，会给数据样本的每个都会附一个权重值，每一个基分类器开始进行分类前，都会根据前一个分类器的分类误差来调节样本的权重值。然后通线性组合的方式把所有的基分类器集成起来得到一个强分类器。整个过程可以描述为下图所示。

image.png

image.png

image.png

image.png

GBM 算法

定义损失函数 L(y,f(x))L(y,f(x)) 来表示预测值f(x)f(x) 与真实值 yy 之间的距离。通过最小化损失函数，就可以得到：

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

损失函数

如果想要求解一个分类问题，而不是回归问题，那应该怎么做呢？首先选择一个合适的损失函数L(y,f) 很重要。因为损失函数决定了模型将如何优化，以及最终会得到什么样的模型。
一般来说， 并不需要自己去发明一个损失函数，这些工作研究人员已经替我们做了。现在来探索用于最常见的两种目标的损失函数：回归和二分类。
回归损失函数
对于回归问题，通常对于不同的问题，不同的数据，选择的损失函数也不同。最常见的回归损失函数是L2损失或高斯损失函数。该损失函数假定数据服从高斯分布。如果没有额外信息，也不要求模型的鲁棒性，那么可以使用该损失函数。公式如下：

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

构建信用评分预测分类模型

信用评分预测任务

问题
信用评分预测任务是一个二元分类问题，主要是预测用户是否会在 90 天内偿还信用借款。
数据集结构示例如下：

image.png

import warnings
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
%matplotlib inline
sns.set()
warnings.filterwarnings('ignore')

from matplotlib import rcParams
rcParams['figure.figsize'] = 11, 8

我们编写一个函数，使用每列数据的中位数替换 NaN 空值。

def fill_nan(table):
    for col in table.columns:
        table[col] = table[col].fillna(table[col].median())
    return table

现在，读取并预览数据集。

data = pd.read_csv(
    'https://labfile.oss.aliyuncs.com/courses/1283/credit_scoring_sample.csv', sep=';')
data.head()

查看各列变量类型：

data.dtypes

通过绘图检查类别是否均衡：

ax = data['SeriousDlqin2yrs'].hist(orientation='horizontal', color='red')
ax.set_xlabel("number_of_observations")
ax.set_ylabel("unique_value")
ax.set_title("Target distribution")

print('Distribution of the target:')
data['SeriousDlqin2yrs'].value_counts()/data.shape[0]

这里，分离出数据特征列名称：

independent_columns_names = [x for x in data if x != 'SeriousDlqin2yrs']
independent_columns_names

使用预先定义好的缺失值处理函数来填充缺失特征数据。

table = fill_nan(data)

最后，将特征数据和目标值分离并赋值：

X = table[independent_columns_names]
y = table['SeriousDlqin2yrs']
X.shape, y.shape

区间估计

对延迟还款的客户的平均年龄进行间隔估计，置信水平设为 90％。区间估计结果是多少？

def get_bootstrap_samples(data, n_samples):
    """Generate samples using bootstrapping."""
    indices = np.random.randint(0, len(data), (n_samples, len(data)))
    samples = data[indices]
    return samples


def stat_intervals(stat, alpha):
    """Make an interval estimate."""
    boundaries = np.percentile(
        stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
    return boundaries


# Save the ages of those who let a delay
churn = data[data['SeriousDlqin2yrs'] == 1]['age'].values

# Set the random seed for reproducibility
np.random.seed(0)

# Generate bootstrap samples and calculate the mean for each sample
churn_mean_scores = [np.mean(sample)
                     for sample in get_bootstrap_samples(churn, 1000)]

# Print the interval estimate for the sample means
print("Mean interval", stat_intervals(churn_mean_scores, 0.1))

逻辑回归分类器

建立逻辑回归分类预测模型。

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV, StratifiedKFold

现在，我们将创建一个 LogisticRegression 模型，并使用 class_weight ='balanced' 来弥补数据中不平衡类。

lr = LogisticRegression(random_state=5, class_weight='balanced')

让我们试着找到模型最佳正则化系数，该系数在逻辑回归中由 C 系数指定，其作用是得到一个不过度拟合的最优模型。下面给出几个 C 参数的备选项。

parameters = {'C': (0.0001, 0.001, 0.01)}

为了找到最佳 C 值，让我们使用 StratifiedKFold 分层采样交叉验证，以确保训练集、测试集中各类别样本的比例与原始数据集中相同。

skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=5)

本次模型质量评估选用曲线下面积（AUC）指标。AUC 值介于 0 到 1 之间，接近 1 则表示分类模型的质量越好。
请使用 scoring='roc_auc' 作为评估指标，并运用 Grid Search 方法来搜索最佳 C 参数，最优值为多少？

grid_search = GridSearchCV(lr, parameters, n_jobs=-1,
                           scoring='roc_auc', cv=skf)
grid_search = grid_search.fit(X, y)
grid_search.best_estimator_

模型最佳 ROC AUC 值为：

grid_search.best_score_

如果模型测试集上的标准差小于 0.005，则认为模型是稳定的。上面得到的最佳 C 值能否保证模型稳定？

grid_search.cv_results_['std_test_score'][1]

特征重要性

特征重要性由其相应系数的绝对值定义。首先，您需要使用 StandardScaler 规范化所有特征值以便于正确比较。那么，上面得到的最佳逻辑回归模型中，哪个特征最为重要？

from sklearn.preprocessing import StandardScaler
lr = LogisticRegression(C=0.001, random_state=5, class_weight='balanced')
scal = StandardScaler()  # 规范化
lr.fit(scal.fit_transform(X), y)

pd.DataFrame({'feat': independent_columns_names,
              'coef': lr.coef_.flatten().tolist()}).sort_values(by='coef', ascending=False)

随机森林

接下来，我们构建随机森林分类预测模型，首先导入所需类。

from sklearn.ensemble import RandomForestClassifier

初始化 100 颗决策树并同样使用类别平衡。

rf = RandomForestClassifier(n_estimators=100, n_jobs=-1, random_state=42,
                            class_weight='balanced')

接下来，从以下候选参数中进行网格搜索。

# parameters = {'max_features': [1, 2, 4], 'min_samples_leaf': [3, 5, 7, 9], 'max_depth': [5, 10, 15]}
# 为了加快搜索速度，这里仅保留原参数中的最佳参数
parameters = {'max_features': [2], 'min_samples_leaf': [7], 'max_depth': [10]}

同样，这里为了找到最佳参数值，使用 StratifiedKFold 分层采样交叉验证。
使用随机森林方法得到的最佳参数比上文中逻辑回归最佳参数在 ROC AUC 指标上高多少？

rf_grid_search = GridSearchCV(
    rf, parameters, n_jobs=-1, scoring='roc_auc', cv=skf, verbose=True)
rf_grid_search = rf_grid_search.fit(X, y)

rf_grid_search.best_score_ - grid_search.best_score_

rf_grid_search.best_estimator_

随机森林方法最优模型中，特征重要性最弱的是拿一项？

independent_columns_names[np.argmin(
    rf_grid_search.best_estimator_.feature_importances_)]

pd.DataFrame({'feat': independent_columns_names,
              'coef': rf_grid_search.best_estimator_.feature_importances_}).sort_values(by='coef', ascending=False)

对于挑战数据集，Logistic 回归相比于随机森林最显着优势是什么？
特征可解释性。
一方面，随机森林模型更适合信用评分问题。它的性能比逻辑回归提高了 4％。另一方面，Logistic 回归的主要优点是我们可以解释特征对模型结果的影响。