线性回归：

线性回归模型

• x 为每个样例中特征值的向量形式，包括 $x_1$ 到 $x_n$ ，而且$x_0$ 恒为1。矩阵点乘，相同规模对应相乘。

损失函数

矩阵形式

一行为一个实例

优化线性回归损失函数的两种方法：1）最小二乘法原理的正态方程 2）梯度下降

正态方程

这里的$X^TX$要满足满秩矩阵，然而，现实大多数任务不会满足这个条件。好像只有线性回归能用正态方程求解。优点是一次计算；缺点是矩阵的逆计算慢，尤其是特征数量很多的情况下就更糟糕了，但是一旦你得到了线性回归模型（通过解正态方程或者其他的算法），进行预测是非常快的。

梯度下降

批量梯度下降

批量梯度下降：使用梯度下降的过程中，你需要计算每一个θj 下代价函数的梯度
代价函数的偏导数: 利用公式2对θj求导，其余 θ看做常数。

更新：

为了避免单独计算每一个梯度，你也可以使用下面的公式来一起计算它们。梯度向量记为$\nabla_{\theta}MSE(\theta) $ ，其包含了代价函数所有的偏导数(每个模型参数只出现一次)。利用正态方程最后的推导即可）

在这个方程中每一步计算时都包含了整个训练集X ，这也是为什么这个算法称为批量梯度下降：每一次训练过程都使用所有的的训练数据。因此，在大数据集上，其会变得相当的慢（但是我们接下来将会介绍更快的梯度下降算法）。然而，梯度下降的运算规模和特征的数量成正比。训练一个数千数量特征的线性回归模型使用梯度下降要比使用正态方程快的多.

更新：

我们来看一下这个算法的应用：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

import numpy as np
eta = 0.03
n_iterations = 15000
m = 100    # 样本数目，

X = 2 * np.random.rand(100, 1)    # 产生100行1列的0~2的数值
y = 4 + 3 * X + np.random.randn(100, 1)
# np.c_表示按列操作拼接，np.r_表示按行操作拼接
X_b = np.c_[np.ones((100, 1)), X]     # x0 = 1
theta = np.random.randn(2, 1)

for iteration in range(n_iterations):
    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
    theta = theta - eta * gradients    # 最后输出的是所有系数矩阵
>>> theta
array([[4.11509616],[2.87011339]])	# 理论θ_0=4, θ_3, 由于噪声会有点误差

梯度下降的一些要点：

• 应该确保所有的特征有着相近的尺度范围（例如：使用Scikit_Learn的 StandardScaler类）
• 学习率$\lambda $ 要自适应

随机梯度下降

批量梯度下降的最要问题是计算每一步的梯度时都需要使用整个训练集，这导致在规模较大的数据集上，其会变得非常的慢。与其完全相反的随机梯度下降，在每一步的梯度计算上只随机选取训练集中的一个样例。

虽然随机性可以很好的跳过局部最优值，但同时它却不能达到最小值。解决这个难题的一个办法是逐渐降低学习率。开始时，走的每一步较大（这有助于快速前进同时跳过局部最小值），然后变得越来越小，从而使算法到达全局最小值。 这个过程被称为模拟退火

下面的代码使用一个简单的learning schedule来实现随机梯度下降：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

n_epochs = 50 
t0, t1 = 5, 50  #learning_schedule的超参数

def learning_schedule(t):
    return t0 / (t + t1)
    
# np.random.randn返回2行1列符合标准正态分布的数；
# np.random.rand返回[0,1）的随机数；
# randint返回范围内的整数
theta = np.random.randn(2,1)

for epoch in range(n_epochs):
    for i in range(m):
        random_index = np.random.randint(m)    # m个样本随机选一个样本
        xi = X_b[random_index:random_index+1]
        yi = y[random_index:random_index+1]
        gradients = 2 * xi.T.dot(xi.dot(theta)-yi)	# 单个个体视为批量
        eta = learning_schedule(epoch * m + i)    # 根据迭代情况调整学习速率
        theta = theta - eta * gradiens
>>> theta
array([[3.96100095],[3.0580351 ]])

通过使用Scikit-Learn完成线性回归的随机梯度下降，你需要使用SGDRegressor类，这个类默认优化的是均方差代价函数。下面的代码迭代了50代，其学习率eta为0.1，使用默认的learning schedule(与前面的不一样)，同时也没有添加任何正则项（penalty = None）：

1
2
3
4
5

# 因为这个函数需要的y是一个行向量，所以压扁;
# 另外，numpy.flatten() 与 numpy.ravel()将多维数组降位一维，前者会进行拷贝处理
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRregressor(n_iter=50, penalty=None, eta=0.1)
sgd_reg.fit(X,y.ravel()) 

结果很接近正态方程的解

1
2

>>> sgd_reg.intercept_, sgd_reg.coef_
(array([4.18380366]),array([2.74205299]))

小批量梯度下降

小批量梯度下降中，它则使用一个随机的小型实例集，小批量梯度下降在参数空间上的表现比随机梯度下降要好的多，尤其在有大量的小型实例集时，主要利用了矩阵运算的硬件优化

也看一下这个算法的应用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

theta_path_mgd = []

n_iterations = 50
minibatch_size = 20

np.random.seed(42)
theta = np.random.randn(2,1)  # random initialization

t0, t1 = 200, 1000
def learning_schedule(t):
    return t0 / (t + t1)

t = 0
for epoch in range(n_iterations):
    shuffled_indices = np.random.permutation(m)
    X_b_shuffled = X_b[shuffled_indices]	# 打乱所有样本顺序
    y_shuffled = y[shuffled_indices]
    for i in range(0, m, minibatch_size):
        t += 1
        xi = X_b_shuffled[i:i+minibatch_size]
        yi = y_shuffled[i:i+minibatch_size]
        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)
        eta = learning_schedule(t)
        theta = theta - eta * gradients
        theta_path_mgd.append(theta)

多项式回归

如果你的数据实际上比简单的直线更复杂呢？ 令人惊讶的是，你依然可以使用线性模型来拟合非线性数据。 一个简单的方法是对每个特征进行加权后作为新的特征，然后训练一个线性模型在这个扩展的特征集。 这种方法称为多项式回归。

于是，我们使用Scikit-Learning的PolynomialFeatures类进行训练数据集的转换，让训练集中每个特征的平方（2次多项式）作为新特征（在这种情况下，仅存在一个特征）：

1
2
3
4
5
6
7

>>> from sklearn.preprocessing import PolynomialFeatures
>>> poly_features = PolynomialFeatures(degree=2,include_bias=False)
>>> X_poly = poly_features.fit_transform(X)    # 转换特征，包含原始特征和二次项特征
>>> X[0]
array([-0.75275929])
>>> X_poly[0]
array([-0.75275929, 0.56664654])

现在包含原始特X并加上了这个特征的平方X^2。现在你可以在这个扩展训练集上使用LinearRegression模型进行拟合

1
2
3
4
5

>>> lin_reg = LinearRegression()
>>> lin_reg.fit(X_poly, y)
>>> lin_reg.intercept_, lin_reg.coef_
(array([ 1.78134581]), array([[ 0.93366893, 0.56456263]]))
# 模型预测函数y=0.56*x_1^2+0.93*x_1+1.78

请注意，当存在多个特征时，多项式回归能够找出特征之间的关系（这是普通线性回归模型无法做到的）。 这是因为LinearRegression会自动添加当前阶数下特征的所有组合。例如，如果有两个特征a,b，使用3阶（degree=3）的LinearRegression时，不仅仅只有a2,a3,b2,同时也会有它们的其他组合项ab,a2b,ab2。

学习曲线

我们可以使用交叉验证来估计一个模型的泛化能力。如果一个模型在训练集上表现良好，通过交叉验证指标却得出其泛化能力很差，那么你的模型就是过拟合了。如果在这两方面都表现不好，那么它就是欠拟合了。这种方法可以告诉我们，你的模型是太复杂了还是太简单了。

另一种方法是观察学习曲线：画出模型在训练集上的表现，同时画出以训练集规模为自变量的训练集函数。为了得到图像，需要在训练集的不同规模子集上进行多次训练。下面的代码定义了一个函数，用来画出给定训练集后的模型学习曲线：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

def plot_learning_curves(model, X, y):
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2)
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):    # 根据样本规模画出模型的表现
        model.fit(X_train[:m], y_train[:m])
        y_train_predict = model.predict(X_train[:m])
        y_val_predict = model.predict(X_val)
        train_errors.append(mean_squared_error(y_train_predict, y_train[:m]))
        val_errors.append(mean_squared_error(y_val_predict, y_val))
plt.plot(np.sqrt(train_errors), "r-+", linewidth=2, label="train")    # 训练损失
plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="val")    # 验证损失

我们一起看一下简单线性回归模型的学习曲线

1
2

lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)

上面的曲线表现了一个典型的欠拟合模型，两条曲线都到达高原地带并趋于稳定，并且最后两条曲线非常接近，同时误差值非常大。

如果你的模型在训练集上是欠拟合的，添加更多的样例是没用的。你需要使用一个更复杂的模型或者找到更好的特征。

现在让我们看一个在相同数据上10阶多项式模型拟合的学习曲线

1
2
3
4
5
6
7
8

from sklearn.pipeline import Pipeline

polynomial_regression = Pipeline((
    ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
    ("sgd_reg", LinearRegression()),
))

plot_learning_curves(polynomial_regression, X, y)

• 在训练集上，误差要比线性回归模型低的多。
• 图中的两条曲线之间有间隔，这意味模型在训练集上的表现要比验证集上好的多，这也是模型过拟合的显著特点。当然，如果你使用了更大的训练数据，这两条曲线最后会非常的接近。

改善模型过拟合的一种方法是提供更多的训练数据，直到训练误差和验证误差相等

在统计和机器学习领域有个重要的理论：一个模型的泛化误差由三个不同误差的和决定：

• 偏差：泛化误差的这部分误差是由于错误的假设决定的。例如实际是一个二次模型，你却假设了一个线性模型。一个高偏差的模型最容易出现欠拟合。
• 方差：这部分误差是由于模型对训练数据的微小变化较为敏感，一个多自由度的模型更容易有高的方差（例如一个高阶多项式模型），因此会导致模型过拟合。
• 不可约误差：这部分误差是由于数据本身的噪声决定的。降低这部分误差的唯一方法就是进行数据清洗（例如：修复数据源，修复坏的传感器，识别和剔除异常值）。

下图依次为欠拟合，过拟合，较合适。

正则化

范数

器学习中有几个常用的范数，分别是：

• $L_1$−范数：$\Vert x\Vert_1 =\sum_{i=1}^n\vert x_i\vert$
• $L_2$−范数：$\Vert x\Vert_ 2=(\sum_{i=1}^n\vert x_i^2\vert)^{\frac{1}{2}}$
• $L_p$−范数：$\Vert x\Vert_p =(\sum_{i=1}^n\vert x_i^p\vert)^{\frac{1}{p}}$
• $L_∞$−范数：$\Vert x\Vert_∞=lim_{p→∞}(\sum_{i=1}^n\vert x_i^p\vert)^{\frac{1}{p}}$

岭回归(Ridge)

岭回归（也称为Tikhonov正则化）是线性回归的正则化版，是L2正则的基础，注意到这个正则项只有在训练过程中才会被加到代价函数。当得到完成训练的模型后，我们应该使用没有正则化的测量方法去评价模型的表现。

一般情况下，训练过程使用的代价函数和测试过程使用的评价函数不一样样的。除了正则化，还有一个不同：训练时的代价函数应该在优化过程中易于求导，而在测试过程中，评价函数更应该接近最后的客观表现。一个好的例子：在分类训练中我们使用对数损失（马上我们会讨论它）作为代价函数，但是我们却使用精确率/召回率来作为它的评价函数。

岭回归代价函数:

超参数α 决定了你想正则化这个模型的强度,正则化强度越大，模型会越简单。如果α=0 那此时的岭回归便变为了线性回归。如果α 非常的大，所有的权重最后都接近与零，最后结果将是一条穿过数据平均值的水平直线

值得注意的是偏差 $\theta_0$是没有被正则化的（累加运算的开始是 i=1而不是i=0）。如我定义$w$作为特征的权重向量($\theta_1$到$\theta_n$)，那么正则项可以简写成$\frac{1}{2} (\Vert w\Vert_2)^2$, 其中$\Vert \cdot \Vert_2$ 表示权重向量的L2范数。对于梯度下降来说仅仅在均方差梯度向量加上一项$\alpha w$ ,加上$\alpha\theta$是$1/2∗\alpha∗\theta^2$求偏导的结果

在使用岭回归前，对数据进行放缩（可以使用StandardScaler）是非常重要的,算法对于输入特征的数值尺度（scale）非常敏感。大多数的正则化模型都是这样的。

对线性回归来说，对于岭回归，我们可以使用封闭方程去计算，也可以使用梯度下降去处理.

岭回归的封闭方程的解

令

求出

矩阵$I$是是一个除了左上角有一个0的n×n的单位矩阵，这个0代表偏差项。偏差$\theta_0$不被正则化的。

下面是如何使用 Scikit-Learn 来进行封闭方程的求解（使用 Cholesky 法进行矩阵分解对上面公式进行变形）:

1
2
3
4
5

>>> from sklearn.linear_model import Ridge
>>> ridge_reg = Ridge(alpha=1, solver="cholesky")
>>> ridge_reg.fit(X, y)
>>> ridge_reg.predict([[1.5]])
array([[ 1.55071465]]

使用随机梯度法进行求解：

1
2
3
4

>>> sgd_reg = SGDRegressor(penalty="l2")
>>> sgd_reg.fit(X, y.ravel())
>>> sgd_reg.predict([[1.5]])
array([[ 1.13500145]])

penalty参数指的是正则项的惩罚类型。指定“l2”表明你要在损失函数上添加一项：权重向量 L2范数平方的一半，这就是简单的岭回归。

Lasso 回归

Lasso 回归（也称 Least Absolute Shrinkage，或者 Selection Operator Regression）是另一种正则化版的线性回归：L1正则的基础，就像岭回归那样，它也在损失函数上添加了一个正则化项，但是它使用权重向量的L1范数而不是权重向量L2范数的一半。

Lasso回归的代价函数:

Lasso回归的一个重要特征是它倾向于完全消除最不重要的特征的权重（即将它们设置为零）

下面是一个使用Lasso类的小Scikit-Learn示例。你也可以使用SGDRegressor(penalty=”l1”)来代替它

1
2
3
4
5

>>> from sklearn.linear_model import Lasso
>>> lasso_reg = Lasso(alpha=0.1)
>>> lasso_reg.fit(X, y)
>>> lasso_reg.predict([[1.5]])
array([ 1.53788174]

弹性网络(ElasticNet)

弹性网络介于Ridge回归和Lasso回归之间。它的正则项是Ridge回归和Lasso回归正则项的简单混合，同时你可以控制它们的混合率r，当r=0时，弹性网络就是Ridge回归，当r=1时，其就是Lasso回归

弹性网络代价函数：

那么我们该如何选择线性回归，岭回归，Lasso回归，弹性网络呢？一般来说有一点正则项的表现更好，因此通常你应该避免使用简单的线性回归。岭回归是一个很好的首选项，但是如果你的特征仅有少数是真正有用的，你应该选择Lasso和弹性网络。就像我们讨论的那样，它两能够将无用特征的权重降为零。一般来说，弹性网络的表现要比Lasso好，因为当特征数量比样例的数量大的时候，或者特征之间有很强的相关性时，Lasso可能会表现的不规律。下面是一个使用Scikit-Learn 弹性网络ElasticNet（l1_ratio指的就是混合率r）的简单样例:

1
2
3
4
5

>>> from sklearn.linear_model import ElasticNet
>>> elastic_net = ElasticNet(alpha=0.1, l1_ratio=0.5)
>>> elastic_net.fit(X, y)
>>> elastic_net.predict([[1.5]])
array([ 1.54333232])

正则化的作用

那为什么正则化能起作用呢？首先L0范数（元素非零个数，严格上来说不能算是范数）和L1范数都可以实现权重稀疏。L1范数和L0范数可以实现稀疏，L1范数是L0范数的最优凸近似，L1因具有比L0更好的优化求解特性而被广泛应用。L1会趋向于产生少量的特征，而其他的特征都是0，而L2会选择更多的特征，这些特征都会接近于0。

L1范数的主要作用的实现稀疏特征，那么L2范数可以起什么样作用呢？

执行 L2 正则化对模型具有以下影响

• 使权重的平均值接近于 0，且呈正态（钟形曲线或高斯曲线）分布。

• 使权重值接近于 0（但并非正好为 0）

L2 正则化可能会导致对于某些信息缺乏的特征，模型会学到适中的权重。L2 正则化降低较大权重的程度高于降低较小权重的程度。随着权重越来越接近于 0.0，L2 将权重“推”向 0.0 的力度越来越弱。L2 正则化会使相似度高(存在噪点)两个特征的权重几乎相同。按照我自己的理解，不同的权重会有不同程度的拟合效果，权重较小，低阶的w控制曲线的整体走势，权重较大，高阶的w控制曲线的局部形态，以此类推。这样看来L2正则项的作用就很明显了，要改变预测曲线的整体细节走势肯地会造成损失函数的不满，但是把曲线的形态熨平似乎并没有什么不妥，会降低过拟合的风险。

L2除了能防止过拟合，提升模型的泛化能力。还有另外的一点好处：优化计算。 从优化或者数值计算的角度来说，L2范数有助于处理 condition number不好的情况下矩阵求逆很困难的问题。conditionnumber是一个矩阵（或者它所描述的线性系统）的稳定性或者敏感度的度量，如果一个矩阵的condition number在1附近，那么它就是well-conditioned的，如果远大于1，那么它就是ill-conditioned的，如果一个系统是ill-conditioned的，它的输出结果就不要太相信了。

然而，如果当我们的样本X的数目比每个样本的维度还要小的时候，矩阵XTX将会不是满秩的，也就是XTX会变得不可逆，所以w*就没办法直接计算出来了。或者更确切地说，将会有无穷多个解（因为我们方程组的个数小于未知数的个数）。也就是说，我们的数据不足以确定一个解，如果我们从所有可行解里随机选一个的话，很可能并不是真正好的解，总而言之，我们过拟合了。

但如果加上L2规则项，就变成了下面这种情况，就可以直接求逆了：

这里面，专业点的描述是：要得到这个解，我们通常并不直接求矩阵的逆，而是通过解线性方程组的方式（例如高斯消元法）来计算。考虑没有规则项的时候，也就是λ=0的情况，如果矩阵XTX的 condition number 很大的话，解线性方程组就会在数值上相当不稳定，而这个规则项的引入则可以改善condition number。

早期停止法（Early Stopping）

随着训练的进行，算法一直学习，它在训练集上的预测误差（RMSE）自然而然的下降。然而一段时间后，验证误差停止下降，并开始上升。这意味着模型在训练集上开始出现过拟合。一旦验证错误达到最小值，便提早停止训练.

随机梯度和小批量梯度下降不是平滑曲线，你可能很难知道它是否达到最小值。 一种解决方案是，只有在验证误差高于最小值一段时间后（你确信该模型不会变得更好了），才停止，之后将模型参数回滚到验证误差最小值。

下面是一个早期停止法的基础应用

1
2
3
4
5
6
7
8
9
10
11
12
13
14

from sklearn.base import clone
sgd_reg = SGDRegressor(n_iter=1, warm_start=True, penalty=None,learning_rate="constant", eta0=0.0005)

minimum_val_error = float("inf")
best_epoch = None
best_model = None
for epoch in range(1000):
    sgd_reg.fit(X_train_poly_scaled, y_train)    # 训练多项式的新特征，拟合非线性
    y_val_predict = sgd_reg.predict(X_val_poly_scaled)
    val_error = mean_squared_error(y_val_predict, y_val)
    if val_error < minimum_val_error:
        minimum_val_error = val_error
        best_epoch = epoch
        best_model = clone(sgd_reg)

注意：当warm_start=True时，调用fit()方法后，训练会从停下来的地方继续，而不是从头重新开始

逻辑回归

逻辑回归会生成一个介于 0 到 1 之间（不包括 0 和 1）的概率值，而不是确切地预测结果是 0 还是 1。以用于检测垃圾邮件的逻辑回归模型为例。如果此模型推断某一特定电子邮件的值为 0.932，则意味着该电子邮件是垃圾邮件的概率为 93.2%。更准确地说，这意味着在无限训练样本的极限情况下，模型预测其值为 0.932 的这组样本实际上有 93.2% 是垃圾邮件，其余的 6.8% 不是垃圾邮件。

逻辑回归模型的概率估计（向量形式）：

Logistic函数（也称为logit），用σ() 表示，其是一个sigmoid函数（图像呈S型），它的输出是一个介于0和1之间的数字
逻辑函数(S函数)

Logistic函数（也称为logit），用σ() 表示，其是一个sigmoid函数（图像呈S型），它的输出是一个介于0和1之间的数字
逻辑函数(S函数)

逻辑回归预测模型(σ() 概率输出以0.5作为二分类门槛):

单个样例的代价函数:

这个代价函数是合理的，因为当t接近0时，-log(t)变得非常大，所以如果模型估计一个正例概率接近于0，那么代价函数将会很大，同时如果模型估计一个负例的概率接近1，那么代价函数同样会很大。 另一方面，当t接近于1时， -log(t)接近0，所以如果模型估计一个正例概率接近于0，那么代价函数接近于0，同时如果模型估计一个负例的概率接近0，那么代价函数同样会接近于0， 这正是我们想的.（简单来说,y=1时，概率p越接近1损失越小；相反y=0时，概率p越接近0时损失越小）

整个训练集的代价函数只是所有训练实例的平均值。可以用一个表达式（你可以很容易证明）来统一表示，称为对数损失

逻辑回归的代价函数（对数损失）：

但是这个代价函数对于求解最小化代价函数的θ 是没有公式解的（没有等价的正态方程）。 但好消息是，这个代价函数是凸的，所以梯度下降（或任何其他优化算法）一定能够找到全局最小值（如果学习速率不是太大，并且你等待足够长的时间）。下面公式给出了代价函数关于第j个模型参数θj 的偏导数。

逻辑回归代价函数的偏导数:

这个公式首先计算每个样例的预测误差，然后误差项乘以第j项特征值，最后求出所有训练样例的平均值。 一旦你有了包含所有的偏导数的梯度向量，你便可以在梯度向量上使用批量梯度下降算法。 也就是说：你已经知道如何训练Logistic回归模型。 对于随机梯度下降，你当然只需要每一次使用一个实例，对于小批量梯度下降，你将每一次使用一个小型实例集。

决策边界

我们使用鸢尾花数据集来分析Logistic回归。 这是一个著名的数据集，其中包含150朵三种不同的鸢尾花的萼片和花瓣的长度和宽度。这三种鸢尾花为：Setosa，Versicolor，Virginica

让我们尝试建立一个分类器，仅仅使用花瓣的宽度特征来**识别Virginica**，首先让我们加载数据：

1
2
3
4
5
6

>>> from sklearn import datasets
>>> iris = datasets.load_iris()
>>> list(iris.keys())
['data', 'target_names', 'feature_names', 'target', 'DESCR']
>>> X = iris["data"][:, 3:] # petal width
>>> y = (iris["target"] == 2).astype(np.int)

接下来，我们训练一个逻辑回归模型：

1
2
3
4

from sklearn.linear_model import LogisticRegression

log_reg = LogisticRegression()
log_reg.fit(X, y) # 训练模型

我们来看看模型估计的花瓣宽度从0到3厘米的概率估计

1
2
3
4

X_new = np.linspace(0, 3, 1000).reshape(-1, 1)    # 构造花瓣宽度从0到3厘米的所有特征
y_proba = log_reg.predict_proba(X_new)    # 预测概率
plt.plot(X_new, y_proba[:, 1], "g-", label="Iris-Virginica")
plt.plot(X_new, y_proba[:, 0], "b--", label="Not Iris-Virginica"

Virginica花的花瓣宽度（用三角形表示）在1.4厘米到2.5厘米之间，而其他种类的花（由正方形表示）通常具有较小的花瓣宽度，范围从0.1厘米到1.8厘米。注意，它们之间会有一些重叠。在大约2厘米以上时，分类器非常肯定这朵花是Virginica花（分类器此时输出一个非常高的概率值），而在1厘米以下时，它非常肯定这朵花不是Virginica花（不是Virginica花有非常高的概率）。在这两个极端之间，分类器是不确定的。但是，如果你使用它进行预测（使用predict()方法而不是predict_proba()方法），它将返回一个最可能的结果。因此，在1.6厘米左右存在一个决策边界，这时两类情况出现的概率都等于50％：如果花瓣宽度大于1.6厘米，则分类器将预测该花是Virginica，否则预测它不是（即使它有可能错了）：

1
2

>>> log_reg.predict([[1.7], [1.5]])
array([1, 0])

下图的线性决策边界表示相同的数据集，但是这次使用了两个特征进行判断：花瓣的宽度和长度。 一旦训练完毕，Logistic回归分类器就可以根据这两个特征来估计一朵花是Virginica的可能性。 虚线表示这时两类情况出现的概率都等于50％：这是模型的决策边界。 请注意，它是一个线性边界。每条平行线都代表一个分类标准下的两两个不同类的概率，从15％（左下角）到90％（右上角）。越过右上角分界线的点都有超过90％的概率是Virginica花

就像其他线性模型，逻辑回归模型也可以ℓ1或者ℓ2 惩罚使用进行正则化。Scikit-Learn默认添加了ℓ2 惩罚

在Scikit-Learn的LogisticRegression模型中控制正则化强度的超参数不是α （与其他线性模型一样），而是是它的逆：C. C的值越大，模型正则化强度越低

Softmax回归

Logistic回归模型可以直接推广到支持多类别分类，不必组合和训练多个二分类器， 其称为Softmax回归或多类别Logistic回归.

这个想法很简单：当给定一个实例x 时，Softmax回归模型首先计算k类的分数sk(x) ，然后将分数应用在Softmax函数（也称为归一化指数）上，估计出每类的概率。 计算sk(x) 的公式看起来很熟悉，因为它就像线性回归预测的公式一样

k类的Softmax得分: $s_k(x)=θ^T⋅x$

注意，每个类都有自己独一无二的参数向量θk 。 所有这些向量通常作为行放在参数矩阵Θ 中

一旦你计算了样例x 的每一类的得分，你便可以通过Softmax函数估计出样例属于第k类的概率p^k ：通过计算e的sk(x) 次方，然后对它们进行归一化（除以所有分子的总和）。

和Logistic回归分类器一样，Softmax回归分类器将估计概率最高（它只是得分最高的类）的那类作为预测结果，如公式4-21所示

Softmax回归分类器一次只能预测一个类（即它是多类的，但不是多输出的），因此它只能用于判断互斥的类别，如不同类型的植物。 你不能用它来识别一张照片中的多个人。

现在我们知道这个模型如何估计概率并进行预测，接下来将介绍如何训练。我们的目标是建立一个模型在目标类别上有着较高的概率（因此其他类别的概率较低），最小化公式4-22可以达到这个目标，其表示了当前模型的代价函数，称为交叉熵，当模型对目标类得出了一个较低的概率，其会惩罚这个模型。 交叉熵通常用于衡量待测类别与目标类别的匹配程度（我们将在后面的章节中多次使用它）

交叉熵

熵的本质是香农信息量$log\frac{1}{p}$的期望。信息熵代表的是随机变量或整个系统的不确定性，熵越大，随机变量或系统的不确定性就越大。在给定的真实分布下，使用非真实分布所指定的策略消除系统的不确定性所需要付出的努力的大小（猜题次数、编码长度等），就是用交叉熵来衡量的。

现有关于样本集的2个概率分布p和q，其中p为真实分布，q非真实分布。按照真实分布p来衡量识别一个样本的所需要的编码长度的期望(即平均编码长度)为$H(p)=\sum \limits_{i=1}^n p(i)\cdot log\frac{1}{p(i)}$ 。如果使用错误分布q来表示来自真实分布p的平均编码长度，则应该是$H(p,q)=\sum\limits_{i=1}^n p(i)\cdot log\frac{1}{q(i)}$ 。因为用q来编码的样本来自分布p，所以期望H(p,q)中概率是p(i)。H(p,q)我们称之为“交叉熵”。当q为真实分布p时，交叉熵达到最小值1，否则将会大于1。我们将由q得到的平均编码长度比由p得到的平均编码长度多出的bit数称为“相对熵”：$D(p\Vert q)=H(p,q)-H(p)=\sum\limits_{i=1}^n p(i)\cdot log\frac{p(i)}{q(i)}$ ,其又被称为KL散度(Kullback–Leibler divergence，KLD)。它表示两个概率分布的差异性：差异越大则相对熵越大，差异越小则相对熵越小，特别地，若2者相同则熵为0。

另外，通常“相对熵”也可称为“交叉熵”，因为真实分布p是固定的，D(p||q)由H(p,q)决定。所以他们得到的相对效果是一样程度的。当然也有特殊情况，彼时两者须区别对待。

上面这个公式由公式4-22求导得到，过程和逻辑回归损失函数一样，只不过将每个类别都纳入计算而已，当k=2则计算正负两类，与逻辑回归一模一样。现在你可以计算每一类的梯度向量，然后使用梯度下降（或者其他的优化算法）找到使得代价函数达到最小值的参数矩阵Θ。

让我们使用Softmax回归对三种鸢尾花进行分类。当你使用LogisticRregression对模型进行训练时，Scikit_Learn默认使用的是一对多模型，但是你可以设置multi_class参数为“multinomial”来把它改变为Softmax回归。你还必须指定一个支持Softmax回归的求解器，例如“lbfgs”求解器（有关更多详细信息，请参阅Scikit-Learn的文档）。其默认使用ℓ12 正则化，你可以使用超参数C控制它。

1
2
3
4
5

X = iris["data"][:, (2, 3)] # petal length, petal width
y = iris["target"]

softmax_reg = LogisticRegression(multi_class="multinomial",solver="lbfgs", C=10)
softmax_reg.fit(X, y)

所以下次你发现一个花瓣长为5厘米，宽为2厘米的鸢尾花时，你可以问你的模型你它是哪一类鸢尾花，它会回答94.2％是Virginica花（第二类），或者5.8％是其他鸢尾花

1
2
3
4

>>> softmax_reg.predict([[5, 2]])
array([2])
>>> softmax_reg.predict_proba([[5, 2]])
array([[ 6.33134078e-07, 5.75276067e-02, 9.42471760e-01]])是

图4-25用不同背景色表示了结果的决策边界。注意，任何两个类之间的决策边界是线性的。 该图的曲线表示Versicolor类的概率（例如，用0.450标记的曲线表示45％的概率边界）。注意模型也可以预测一个概率低于50％的类。 例如，在所有决策边界相遇的地方，所有类的估计概率相等，分别为33％。

练习题

1. 如果你有一个数百万特征的训练集，你应该选择哪种线性回归训练算法？
2. 假设你训练集中特征的数值尺度（scale）有着非常大的差异，哪种算法会受到影响？有多大的影响？对于这些影响你可以做什么？
3. 训练 Logistic 回归模型时，梯度下降是否会陷入局部最低点？
4. 在有足够的训练时间下，是否所有的梯度下降都会得到相同的模型参数？
5. 假设你使用批量梯度下降法，画出每一代的验证误差。当你发现验证误差一直增大，接下来会发生什么？你怎么解决这个问题？
6. 当验证误差升高时，立即停止小批量梯度下降是否是一个好主意？
7. 哪个梯度下降算法（在我们讨论的那些算法中）可以最快到达解的附近？哪个的确实会收敛？怎么使其他算法也收敛？
8. 假设你使用多项式回归，画出学习曲线，在图上发现学习误差和验证误差之间有着很大的间隙。这表示发生了什么？有哪三种方法可以解决这个问题？
9. 假设你使用岭回归，并发现训练误差和验证误差都很高，并且几乎相等。你的模型表现是高偏差还是高方差？这时你应该增大正则化参数$\alpha$ 还是降低它？
10. 你为什么要这样做：

• 使用岭回归代替线性回归？
• Lasso 回归代替岭回归？
• 弹性网络代替 Lasso 回归？

1. 假设你想判断一副图片是室内还是室外，白天还是晚上。你应该选择二个逻辑回归分类器，还是一个 Softmax 分类器？

1、如果您拥有具有数百万个功能的训练集，则可以使用随机梯度下降或小批量梯度下降，如果计算内存足够的话，则可使用批量梯度下降。 但是你不能使用正态方程，因为计算复杂度随着特征数量的增长而快速增长（超过二次方），求矩阵特征的逆非常花时间。

2、如果训练集中的特征具有非常不同的比例，则损失函数将具有细长碗的形状，因此梯度下降优化将花费很长时间来收敛。 要解决此问题，您应该在训练模型之前缩放数据。 另外，正态方程在没有缩放的情况下可以正常工作。

3、在训练Logistic回归模型时，梯度下降不会陷入在局部最小值，因为它的损失函数是凸函数的。

4、如果优化问题是凸函数的（例如线性回归或逻辑回归），并且假设学习速率不是太高，则所有梯度下降算法将接近全局最优并最终产生相当类似的模型。 但是，除非你逐渐降低学习率，否则随机梯度下降和小批量GD将永远不会真正收敛; 相反，他们将继续围绕全局最佳状态来回跳跃。 这意味着即使你让它们运行很长时间，这些Gradient Descent算法也会产生略微不同的模型。

5、如果验证误差在每个时期之后一直上升，则一种可能性是学习速率太高并且算法发散。如果训练误差也会增加，那么这显然是问题，你应该降低学习率。 但是，如果训练错误没有增加，那么您的模型将过度拟合训练集，您应该停止训练。

6、由于随机性，随机梯度下降和小批量梯度下降都不能保证在每次训练迭代中都取得进展。 因此，如果在验证损失增加时立即停止训练，你可能会在达到最佳值之前过早停止。 更好的选择是定期保存模型，当它长时间没有改进时（意味着它可能永远不会超过记录），你可以恢复到最佳保存模型。

7、随机梯度下降具有最快的训练迭代，因为它一次只考虑一个训练实例，因此它通常是第一个到达全局最优值（或具有非常小的小批量大小的Minibatch GD）附近。 但是，如果有足够的训练时间，只有批量梯度下降实际上会收敛。 如上所述，除非你逐渐降低学习速度，否则随机指标GD和小批量GD将在最佳状态下反弹。

8、如果验证误差远远高于训练误差，则可能是因为你的模型过度拟合了训练集。 尝试解决此问题的一种方法是降低多项式度：具有较少自由度的模型不太可能过度拟合。 你可以尝试的另一件事是加入正则项，例如，通过在成本函数中添加ℓ2惩罚（岭）或ℓ1惩罚（Lasso）。 这也会降低模型的自由度。 最后，你还可以尝试增加训练集的大小。

9、如果训练误差和验证误差几乎相等且相当高，则模型可能欠拟合训练集，这意味着它具有高偏差。 你应该尝试减少正则化超参数α。

10、

• 具有一些正则化的模型通常比没有任何正则化的模型表现更好，因此通常应该优先选择岭回归而不是简单的线性回归。
• Lasso回归使用ℓ1惩罚，这往往会将权重降低到恰好为零。 这导致稀疏模型，除了最重要的权重之外，所有权重都为零。 这是一种自动执行特征选择的方法，如果你怀疑只有少数特征真正重要，这是很好的。 当你不确定时，你应该更偏向岭回归。
• 弹性网络常比Lasso更受欢迎，因为Lasso在某些情况下可能表现不稳定（当有些特征强烈相关或者特征数量比训练样本数量还要多）。 但是，它确实添加了一个额外的超参数来调整。 如果你想要具有稳定行为的Lasso，你可以使用弹性网络，并设置比率r接近1。

11、如果你想将图片分类为室外/室内和白天/夜晚，因为这些不是专属类别（即，所有四种组合都是可能的），你应该训练两个Logistic回归分类器。

MNIST：手写数字分类数据集

1
2
3
4
5
6
7
8
9
10
11
12
13
14

>>> from sklearn.datasets import fetch_mldata
>>> mnist = fetch_mldata('MNIST original')
>>> mnist
{'COL_NAMES': ['label', 'data'],
'DESCR': 'mldata.org dataset: mnist-original',   # DESCR键描述数据集
'data': array([[0, 0, 0, ..., 0, 0, 0],    # 数组的一行表示一个样例，一列表示一个特征
                [0, 0, 0, ..., 0, 0, 0],
                [0, 0, 0, ..., 0, 0, 0],
                ...,
                [0, 0, 0, ..., 0, 0, 0],
                [0, 0, 0, ..., 0, 0, 0],
                [0, 0, 0, ..., 0, 0, 0]], dtype=uint8),
'target': array([ 0., 0., 0., ..., 9., 9., 9.])}    # target键存放一个标签数组
X, y = mnist["data"], mnist["target"] # 获取样本或标签

MNIST 有 70000 张图片，每张图片有 784 个特征。这是因为每个图片都是28×28像素的，并且每个像素的值介于 0~255 之间。让我们看一看数据集的某一个数字。你只需要将某个实例的特征向量，reshape为28*28的数组，然后使用 Matplotlib 的imshow函数展示出来

1
2
3
4
5
6
7
8
9
10

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
some_digit = X[36000]    
some_digit_image = some_digit.reshape(28, 28)    # 将样本转为28大小的像素矩阵
# 按‘0’‘1’数值转为灰度图像 
# interpolation当小图像放大时,interpolation ='nearest'效果很好，否则用None。
plt.imshow(some_digit_image, cmap = matplotlib.cm.binary, interpolation="nearest")    
plt.axis("off")
plt.show()

MNIST 数据集已经事先被分成了一个训练集（前 6000 张图片）和一个测试集（最后 10000 张图片）

1

X_train, X_test, y_train, y_test = X[:60000], X[60000:], y[:60000], y[60000:]

打乱数据

1
2
3
4

import numpy as np

shuffle_index = np.random.permutation(60000)
X_train, y_train = X_train[shuffle_index], y_train[shuffle_index]

训练一个二分类器

现在我们简化一下问题，只尝试去识别一个数字，比如说，数字 5。这个“数字 5 检测器”就是一个二分类器，能够识别两类别，“是 5”和“非 5”。让我们为这个分类任务创建目标向量：

1
2
3

# 在训练和测试集上区分是否为5转为0,1标签矩阵
y_train_5 = (y_train == 5) # True for all 5s, False for all other digits.
y_test_5 = (y_test == 5)

采用随机梯度下降分类器

1
2
3

from sklearn.linear_model import SGDClassifier
sgd_clf = SGDClassifier(random_state=42)    #如果你想重现结果，你应该固定参数random_state 
sgd_clf.fit(X_train, y_train_5)

输出预测结果

1
2

>>> sgd_clf.predict([some_digit])
array([ True], dtype=bool)

分类器猜测这个数字代表 5（True）。看起来在这个例子当中，它猜对了。现在让我们评估这个模型的性能。

使用交叉验证测量准确性

评估一个模型的好方法是使用交叉验证，像之前提过一样。但有时为了有更好的控制权，可以写自己版本的交叉验证，以下代码粗略地做了和cross_val_score()相同的事情，并且输出相同的结果。

1
2
3
4
5
6
7
8
9
10
11
12
13

from sklearn.model_selection import StratifiedKFold
from sklearn.base import clone
skfolds = StratifiedKFold(n_splits=3, random_state=42)    # 三组
for train_index, test_index in skfolds.split(X_train, y_train_5):
    clone_clf = clone(sgd_clf)
    X_train_folds = X_train[train_index]
    y_train_folds = (y_train_5[train_index])
    X_test_fold = X_train[test_index]
    y_test_fold = (y_train_5[test_index])
    clone_clf.fit(X_train_folds, y_train_folds)
    y_pred = clone_clf.predict(X_test_fold)
    n_correct = sum(y_pred == y_test_fold)
    print(n_correct / len(y_pred)) # prints 0.9502, 0.96565 and 0.96495

StratifiedKFold类实现了分层采样，生成的折（fold）包含了各类相应比例的样例。在每一次迭代，上述代码生成分类器的一个克隆版本，在训练折（training folds）的克隆版本上进行训，在测试折（test folds）上进行预测。然后它计算出被正确预测的数目和输出正确预测的比例。

这里使用sklearn提供的cross_val_score()函数来评估SGDClassifier模型

1
2
3

>>> from sklearn.model_selection import cross_val_score
>>> cross_val_score(sgd_clf, X_train, y_train_5, cv=3, scoring="accuracy")
array([ 0.9502 , 0.96565, 0.96495]

有大于 95% 的精度（accuracy），特别高！但要注意这是一个有数据偏差的数据集，这是因为只有 10% 的图片是数字 5，所以你总是猜测某张图片不是 5，你也会有90%的可能性是对的。处理这类问题，要回归到之前讲的准确率和召回率和ORC曲线了。

混淆矩阵

对分类器来说，一个好得多的性能评估指标是混淆矩阵，为了计算混淆矩阵，首先你需要有一系列的预测值，这样才能将预测值与真实值做比较。你或许想在测试集上做预测。但是我们现在先不碰它。（记住，只有当你处于项目的尾声，当你准备上线一个分类器的时候，你才应该使用测试集）。相反，你应该使用cross_val_predict()函数

1
2

from sklearn.model_selection import cross_val_predict
y_train_pred = cross_val_predict(sgd_clf, X_train, y_train_5, cv=3)

就像 cross_val_score()，cross_val_predict()也使用 K 折交叉验证。它不是返回一个评估分数，而是返回基于每一个测试折做出的一个预测值。这意味着，对于每一个训练集的样例，你得到一个干净的预测（“干净”是说一个模型在训练过程当中没有用到测试集的数据）。

现在使用 confusion_matrix()函数，你将会得到一个混淆矩阵。传递目标类(y_train_5)和预测类（y_train_pred）给它。

1
2
3
4

>>> from sklearn.metrics import confusion_matrix
>>> confusion_matrix(y_train_5, y_train_pred)
array([[53272, 1307],
        [ 1077, 4344]])

混淆矩阵中的每一行表示一个实际的类, 而每一列表示一个预测的类。该矩阵的第一行认为“非 5”（反例）中的 53272 张被正确归类为 “非 5”（他们被称为真反例，true negatives）, 而其余 1307 被错误归类为”是 5” （假正例，false positives）。第二行认为“是 5” （正例）中的 1077 被错误地归类为“非 5”（假反例，false negatives），其余 4344 正确分类为 “是 5”类（真正例，true positives）。一个完美的分类器将只有真反例和真正例，所以混淆矩阵的非零值仅在其主对角线（左上至右下）。

Scikit-Learn 提供了一些函数去计算分类器的指标，包括精确率和召回率（之前的文章是tensorflow，这里主要讲Scikit-Learn）

1
2
3
4
5

>>> from sklearn.metrics import precision_score, recall_score
>>> precision_score(y_train_5, y_pred) # == 4344 / (4344 + 1307)
0.76871350203503808
>>> recall_score(y_train_5, y_train_pred) # == 4344 / (4344 + 1077)
0.79136690647482011

通常结合精确率和召回率会更加方便，这个指标叫做“F1 值”，特别是当你需要一个简单的方法去比较两个分类器的优劣的时候。F1 值是精确率和召回率的调和平均。普通的平均值平等地看待所有的值，而调和平均会给小的值更大的权重。所以，要想分类器得到一个高的 F1 值，需要召回率和精确率。

为了计算 F1 值，简单调用f1_score()

1
2
3

>>> from sklearn.metrics import f1_score
>>> f1_score(y_train_5, y_pred)
0.78468208092485547

F1 支持那些有着相近精确率和召回率的分类器。这不会总是你想要的。有的场景你会绝大程度地关心精确率，而另外一些场景你会更关心召回率。不幸的是，你不能同时拥有两者。增加精确率会降低召回率，反之亦然。这叫做精确率与召回率之间的折衷. 一般来说，提高分类阈值会减少假正例，从而提高精确率。降低分类阈值会提高召回率。

Scikit-Learn 不让你直接设置阈值，但是它给你提供了设置决策分数的方法，这个决策分数可以用来产生预测。它不是调用分类器的predict()方法，而是调用decision_function()方法。这个方法返回每一个样例的分数值，然后基于这个分数值，使用你想要的任何阈值做出预测。

1
2
3
4
5
6

>>> y_scores = sgd_clf.decision_function([some_digit])
>>> y_scores
array([ 161855.74572176])
>>> threshold = 0
>>> y_some_digit_pred = (y_scores > threshold)
array([ True], dtype=bool)

SGDClassifier用了一个等于 0 的阈值，所以前面的代码返回了跟predict()方法一样的结果（都返回了true）。让我们提高这个阈值：

1
2
3
4

>>> threshold = 200000
>>> y_some_digit_pred = (y_scores > threshold)
>>> y_some_digit_pred
array([False], dtype=bool)

这证明了提高阈值会降调召回率。这个图片实际就是数字 5，当阈值等于 0 的时候，分类器可以探测到这是一个 5，当阈值提高到 20000 的时候，分类器将不能探测到这是数字 5。

那么，你应该如何使用哪个阈值呢？首先，你需要再次使用cross_val_predict()得到每一个样例的分数值，但是这一次指定返回一个决策分数，而不是预测值。

1
2

y_scores = cross_val_predict(sgd_clf, X_train, y_train_5, cv=3, 
                            method="decision_function")

现在有了这些分数值。对于任何可能的阈值，使用precision_recall_curve(),你都可以计算精确率和召回率:

1
2

from sklearn.metrics import precision_recall_curve
precisions, recalls, thresholds = precision_recall_curve(y_train_5, y_scores)

最后，你可以使用 Matplotlib 画出精确率和召回率，这里把精确率和召回率当作是阈值的一个函数。

1
2
3
4
5
6
7
8

def plot_precision_recall_vs_threshold(precisions, recalls, thresholds):
    plt.plot(thresholds, precisions[:-1], "b--", label="Precision")
    plt.plot(thresholds, recalls[:-1], "g-", label="Recall")
    plt.xlabel("Threshold")
    plt.legend(loc="upper left")
    plt.ylim([0, 1])
plot_precision_recall_vs_threshold(precisions, recalls, thresholds)
plt.show()

你也许会好奇为什么精确率曲线比召回率曲线更加起伏不平（右上部分）。原因是精确率有时候会降低，尽管当你提高阈值的时候，通常来说精确率会随之提高。另一方面，当阈值提高时候，召回率只会降低。这也就说明了为什么召回率的曲线更加平滑。

现在你可以选择适合你任务的最佳阈值。另一个选出好的精确率/召回率折衷的方法是直接画出精确率对召回率的曲线(PR曲线)，如图所示。

我们假设你决定达到 90% 的准确率，在 70000 附近找到一个阈值。为了作出预测（目前为止只在训练集上预测），你可以运行以下代码，而不是运行分类器的predict()方法。

1

y_train_pred_90 = (y_scores > 70000)

检查这些预测的准确率和召回率：

1
2
3
4

>>> precision_score(y_train_5, y_train_pred_90)
0.8998702983138781
>>> recall_score(y_train_5, y_train_pred_90)
0.63991883416343853

ROC 曲线

受试者工作特征（ROC）曲线是另一个二分类器常用的工具。它非常类似与准确率/召回率曲线（PR曲线），但不是画出准确率对召回率的曲线，ROC 曲线是真正例率（true positive rate，另一个名字叫做召回率）对假正例率（false positive rate, FPR）的曲线。FPR 是反例被错误分成正例的比率。它等于 1 减去真反例率（true negative rate， TNR）。TNR是反例被正确分类的比率。TNR也叫做特异性。所以 ROC 曲线画出召回率对（1 减特异性）的曲线。

为了画出 ROC 曲线，你首先需要计算各种不同阈值下的 TPR、FPR，使用roc_curve()函数：

1
2

from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(y_train_5, y_scores)

然后你可以使用 matplotlib，画出 FPR 对 TPR 的曲线

1
2
3
4
5
6
7
8

def plot_roc_curve(fpr, tpr, label=None):
    plt.plot(fpr, tpr, linewidth=2, label=label)
    plt.plot([0, 1], [0, 1], 'k--')
    plt.axis([0, 1, 0, 1])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
plot_roc_curve(fpr, tpr)
plt.show()

一个比较分类器之间优劣的方法是：测量ROC曲线下的面积（AUC）**。一个完美的分类器的 ROC AUC 等于 1，而一个纯随机分类器的 ROC AUC 等于 0.5。Scikit-Learn 提供了一个函数来计算 ROC AUC：

1
2
3

>>> from sklearn.metrics import roc_auc_score
>>> roc_auc_score(y_train_5, y_scores)
0.97061072797174941

因为 ROC 曲线跟准确率/召回率曲线（或者叫 PR）很类似，你或许会好奇如何决定使用哪一个曲线呢？一个笨拙的规则是，优先使用 PR 曲线当正例很少，或者当你关注假正例多于假反例的时候。其他情况使用 ROC 曲线。举例子，回顾前面的 ROC 曲线和 ROC AUC 数值，你或许人为这个分类器很棒。但是这几乎全是因为只有少数正例（“是 5”），而大部分是反例（“非 5”）。相反，PR 曲线清楚显示出这个分类器还有很大的改善空间（PR 曲线应该尽可能地靠近右上角）。

我们训练一个RandomForestClassifier，然后拿它的的ROC曲线和ROC AUC数值去跟SGDClassifier的比较。首先你需要得到训练集每个样例的数值。但是由于随机森林分类器的工作方式，RandomForestClassifier不提供decision_function()方法。相反，它提供了predict_proba()方法。Skikit-Learn分类器通常二者中的一个。predict_proba()方法返回一个数组，数组的每一行代表一个样例，每一列代表一个类。数组当中的值的意思是：给定一个样例属于给定类的概率。比如，70%的概率这幅图是数字 5。

1
2
3
4

from sklearn.ensemble import RandomForestClassifier
forest_clf = RandomForestClassifier(random_state=42)
y_probas_forest = cross_val_predict(forest_clf, X_train, y_train_5, cv=3,
                                    method="predict_proba")

但是要画 ROC 曲线，你需要的是样例的分数，而不是概率。一个简单的解决方法是使用正例的概率当作样例的分数。

1
2

y_scores_forest = y_probas_forest[:, 1] # score = proba of positive class 预测为正例概率
fpr_forest, tpr_forest, thresholds_forest = roc_curve(y_train_5,y_scores_forest)

现在你即将得到 ROC 曲线。将前面一个分类器的 ROC 曲线一并画出来是很有用的，可以清楚地进行比较。

1
2
3
4

plt.plot(fpr, tpr, "b:", label="SGD")
plot_roc_curve(fpr_forest, tpr_forest, "Random Forest")
plt.legend(loc="bottom right")
plt.show()

如你所见，RandomForestClassifier的 ROC 曲线比SGDClassifier的好得多：它更靠近左上角。所以，它的 ROC AUC 也会更大。

1
2

>>> roc_auc_score(y_train_5, y_scores_forest)
0.99312433660038291

现在你知道如何训练一个二分类器，选择合适的标准，使用交叉验证去评估你的分类器，选择满足你需要的准确率/召回率折衷方案，和比较不同模型的 ROC 曲线和 ROC AUC 数值。现在让我们检测更多的数字，而不仅仅是一个数字 5。

多类别分类

一些算法（比如随机森林分类器或者朴素贝叶斯分类器）可以直接处理多类分类问题。其他一些算法（比如 SVM 分类器或者线性分类器）则是严格的二分类器。然后，有许多策略可以让你用二分类器去执行多类分类。

• 一个方法是：训练10个二分类器，每一个对应一个数字（探测器 0，探测器 1，探测器 2，以此类推）。然后当你想对某张图片进行分类的时候，让每一个分类器对这个图片进行分类，选出决策分数最高的那个分类器（One vs all 里面分数最高的One）。这叫做“一对所有”（OvA）策略
• 另一个策略是对每一对数字都训练一个二分类器：一个分类器用来处理数字 0 和数字 1，一个用来处理数字 0 和数字 2，一个用来处理数字 1 和 2，以此类推。这叫做“一对一”（OvO）策略。如果有 N 个类。你需要训练N*(N-1)/2个分类器。

一些算法（比如 SVM 分类器）在训练集的大小上很难扩展，所以对于这些算法，OvO 是比较好的，因为它可以在小的数据集上面可以更多地训练，较之于巨大的数据集而言。但是，对于大部分的二分类器来说，OvA 是更好的选择。Scikit-Learn 可以探测出你想使用一个二分类器去完成多分类的任务，它会自动地执行 OvA（除了 SVM 分类器，它使用 OvO）让我们试一下SGDClassifier.

1
2
3

>>> sgd_clf.fit(X_train, y_train) # y_train, not y_train_5
>>> sgd_clf.predict([some_digit])
array([ 5.])

上面的代码在训练集上训练了一个SGDClassifier。这个分类器处理原始的目标class，从 0 到 9（y_train），而不是仅仅探测是否为 5 （y_train_5）。然后它做出一个判断（在这个案例下只有一个正确的数字）。在幕后，Scikit-Learn 实际上训练了 10 个二分类器，每个分类器都产到一张图片的决策数值，选择数值最高的那个类。

为了证明这是真实的，你可以调用decision_function()方法。不是返回每个样例的一个数值，而是返回 10 个数值，一个数值对应于一个类。

1
2
3
4
5
6

>>> some_digit_scores = sgd_clf.decision_function([some_digit])
>>> some_digit_scores
array([[-311402.62954431, -363517.28355739, -446449.5306454 ,
        -183226.61023518, -414337.15339485, 161855.74572176,
        -452576.39616343, -471957.14962573, -518542.33997148,
        -536774.63961222]])

最高数值是对应于类别 5 :

1
2
3
4
5
6

>>> np.argmax(some_digit_scores)    # 找最大值的索引
5
>>> sgd_clf.classes_
array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])
>>> sgd_clf.classes[5]    # 用索引匹配类别
5.0

一个分类器被训练好了之后，它会保存目标类别列表到它的属性classes_ 中去，按照值排序。在本例子当中，在classes_ 数组当中的每个类的索引方便地匹配了类本身，比如，索引为 5 的类恰好是类别 5 本身。但通常不会这么幸运。

如果你想强制 Scikit-Learn 使用 OvO 策略或者 OvA 策略，你可以使用OneVsOneClassifier类或者OneVsRestClassifier类。创建一个样例，传递一个二分类器给它的构造函数。举例子，下面的代码会创建一个多类分类器，使用 OvO 策略，基于SGDClassifier。

1
2
3
4
5
6
7

>>> from sklearn.multiclass import OneVsOneClassifier
>>> ovo_clf = OneVsOneClassifier(SGDClassifier(random_state=42))
>>> ovo_clf.fit(X_train, y_train)
>>> ovo_clf.predict([some_digit])
array([ 5.])
>>> len(ovo_clf.estimators_)
45

训练一个RandomForestClassifier同样简单：

1
2
3

>>> forest_clf.fit(X_train, y_train)
>>> forest_clf.predict([some_digit])
array([ 5.])

这次 Scikit-Learn 没有必要去运行 OvO 或者 OvA，因为随机森林分类器能够直接将一个样例分到多个类别。你可以调用predict_proba()，得到样例对应的类别的概率值的列表：

1
2

>>> forest_clf.predict_proba([some_digit])
array([[ 0.1, 0. , 0. , 0.1, 0. , 0.8, 0. , 0. , 0. , 0. ]])

你可以看到这个分类器相当确信它的预测：在数组的索引 5 上的 0.8，意味着这个模型以 80% 的概率估算这张图片代表数字 5。它也认为这个图片可能是数字 0 或者数字 3，分别都是 10% 的几率。

现在当然你想评估这些分类器。像平常一样，你想使用交叉验证。让我们用cross_val_score()来评估SGDClassifier的精度。

1
2

>>> cross_val_score(sgd_clf, X_train, y_train, cv=3, scoring="accuracy")
array([ 0.84063187, 0.84899245, 0.86652998])

在所有测试折（test fold）上，它有 84% 的精度。如果你是用一个随机的分类器，你将会得到 10% 的正确率。所以这不是一个坏的分数，但是你可以做的更好。举例子，简单将输入正则化，将会提高精度到 90% 以上。

1
2
3
4
5

>>> from sklearn.preprocessing import StandardScaler
>>> scaler = StandardScaler()
>>> X_train_scaled = scaler.fit_transform(X_train.astype(np.float64))    # 特征正则化，没说用哪种
>>> cross_val_score(sgd_clf, X_train_scaled, y_train, cv=3, scoring="accuracy")
array([ 0.91011798, 0.90874544, 0.906636 ])

误差分析：

首先，你可以检查混淆矩阵。你需要使用cross_val_predict()做出预测，然后调用confusion_matrix()函数，像你早前做的那样。

1
2
3
4
5
6
7
8
9
10
11
12
13

>>> y_train_pred = cross_val_predict(sgd_clf, X_train_scaled, y_train, cv=3)
>>> conf_mx = confusion_matrix(y_train, y_train_pred)
>>> conf_mx
array([[5725, 3, 24, 9, 10, 49, 50, 10, 39, 4],
        [ 2, 6493, 43, 25, 7, 40, 5, 10, 109, 8],
        [ 51, 41, 5321, 104, 89, 26, 87, 60, 166, 13],
        [ 47, 46, 141, 5342, 1, 231, 40, 50, 141, 92],
        [ 19, 29, 41, 10, 5366, 9, 56, 37, 86, 189],
        [ 73, 45, 36, 193, 64, 4582, 111, 30, 193, 94],
        [ 29, 34, 44, 2, 42, 85, 5627, 10, 45, 0],
        [ 25, 24, 74, 32, 54, 12, 6, 5787, 15, 236],
        [ 52, 161, 73, 156, 10, 163, 61, 25, 5027, 123],
        [ 43, 35, 26, 92, 178, 28, 2, 223, 82, 5240]])

这里是一对数字。使用 Matplotlib 的matshow()函数，将混淆矩阵以图像的方式呈现，将会更加方便

1
2

plt.matshow(conf_mx, cmap=plt.cm.gray)    # #灰度图,对应位置的值越大色块越亮
plt.show()

这个混淆矩阵看起来相当好，因为大多数的图片在主对角线上。在主对角线上意味着被分类正确。数字 5 对应的格子看起来比其他数字要暗淡许多。这可能是数据集当中数字 5 的图片比较少，又或者是分类器对于数字 5 的表现不如其他数字那么好。你可以验证两种情况.

让我们关注仅包含误差数据的图像呈现。首先你需要将混淆矩阵的每一个值除以相应类别的图片的总数目。这样子，你可以比较错误率，而不是绝对的错误数（这对大的类别不公平）。

1
2

row_sums = conf_mx.sum(axis=1, keepdims=True)
norm_conf_mx = conf_mx / row_sums

现在让我们用 0 来填充对角线。这样子就只保留了被错误分类的数据。让我们画出这个结果。(此时数值为错误率)

1
2
3

np.fill_diagonal(norm_conf_mx, 0)
plt.matshow(norm_conf_mx, cmap=plt.cm.gray)
plt.show()

现在你可以清楚看出分类器制造出来的各类误差。记住：行代表实际类别，列代表预测的类别。第 8、9 列相当亮，这告诉你许多图片被误分成数字 8 或者数字 9。相似的，第 8、9 行也相当亮，告诉你数字 8、数字 9 经常被误以为是其他数字。相反，一些行相当黑，比如第一行：这意味着大部分的数字 1 被正确分类（一些被误分类为数字 8 ）。留意到误差图不是严格对称的。举例子，比起将数字 8 误分类为数字 5 的数量，有更多的数字 5 被误分类为数字 8。

分析混淆矩阵通常可以给你提供深刻的见解去改善你的分类器。回顾这幅图，看样子你应该努力改善分类器在数字 8 和数字 9 上的表现，和纠正 3/5 的混淆。举例子，你可以尝试去收集更多的数据，或者你可以构造新的、有助于分类器的特征。举例子，写一个算法去数闭合的环（比如，数字 8 有两个环，数字 6 有一个， 5 没有）。又或者你可以预处理图片（比如，使用 Scikit-Learn，Pillow， OpenCV）去构造一个模式，比如闭合的环。

分析独特的误差，是获得关于你的分类器是如何工作及其为什么失败的洞见的一个好途径。但是这相对难和耗时。举例子，我们可以画出数字 3 和 5 的例子

1
2
3
4
5
6
7
8
9
10
11

cl_a, cl_b = 3, 5
X_aa = X_train[(y_train == cl_a) & (y_train_pred == cl_a)]
X_ab = X_train[(y_train == cl_a) & (y_train_pred == cl_b)]
X_ba = X_train[(y_train == cl_b) & (y_train_pred == cl_a)]
X_bb = X_train[(y_train == cl_b) & (y_train_pred == cl_b)]
plt.figure(figsize=(8,8))
plt.subplot(221); plot_digits(X_aa[:25], ../images_per_row=5)
plt.subplot(222); plot_digits(X_ab[:25], ../images_per_row=5)
plt.subplot(223); plot_digits(X_ba[:25], ../images_per_row=5)
plt.subplot(224); plot_digits(X_bb[:25], ../images_per_row=5)
plt.show()

左边两个5*5的块将数字识别为 3，右边的将数字识别为 5。一些被分类器错误分类的数字（比如左下角和右上角的块）是书写地相当差，甚至让人类分类都会觉得很困难（比如第 8 行第 1 列的数字 5，看起来非常像数字 3 ）。但是，大部分被误分类的数字，在我们看来都是显而易见的错误。很难明白为什么分类器会分错。原因是我们使用的简单的SGDClassifier，这是一个线性模型。它所做的全部工作就是分配一个类权重给每一个像素，然后当它看到一张新的图片，它就将加权的像素强度相加，每个类得到一个新的值。所以，因为 3 和 5 只有一小部分的像素有差异，这个模型很容易混淆它们。

3 和 5 之间的主要差异是连接顶部的线和底部的线的细线的位置。如果你画一个 3，连接处稍微向左偏移，分类器很可能将它分类成 5。反之亦然。换一个说法，这个分类器对于图片的位移和旋转相当敏感。所以，减轻 3/5 混淆的一个方法是对图片进行预处理，确保它们都很好地中心化和不过度旋转。这同样很可能帮助减轻其他类型的错误。

多标签分类

先看一个简单点的例子，仅仅是为了阐明的目的

1
2
3
4
5
6

from sklearn.neighbors import KNeighborsClassifier
y_train_large = (y_train >= 7)
y_train_odd = (y_train % 2 == 1)
y_multilabel = np.c_[y_train_large, y_train_odd]
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_multilabel)

这段代码创造了一个y_multilabel数组，里面包含两个目标标签。第一个标签指出这个数字是否为大数字（7，8 或者 9），第二个标签指出这个数字是否是奇数。接下来几行代码会创建一个KNeighborsClassifier样例（它支持多标签分类，但不是所有分类器都可以），然后我们使用多目标数组来训练它。现在你可以生成一个预测，然后它输出两个标签：

1
2

>>> knn_clf.predict([some_digit])
array([[False, True]], dtype=bool)

它工作正确。数字 5 不是大数（False），同时是一个奇数（True）

有许多方法去评估一个多标签分类器，和选择正确的量度标准，这取决于你的项目。举个例子，一个方法是对每个个体标签去量度 F1 值（或者前面讨论过的其他任意的二分类器的量度标准），然后计算平均值。下面的代码计算全部标签的平均 F1 值：

1
2
3

>>> y_train_knn_pred = cross_val_predict(knn_clf, X_train, y_train, cv=3)
>>> f1_score(y_train, y_train_knn_pred, average="macro")
0.96845540180280221

这里假设所有标签有着同等的重要性，但可能不是这样。特别是，如果你的 Alice 的照片比 Bob 或者 Charlie 更多的时候，也许你想让分类器在 Alice 的照片上具有更大的权重。一个简单的选项是：给每一个标签的权重等于它的支持度（比如，那个标签的样例的数目）。为了做到这点，简单地在上面代码中设置average=”weighted”。

多输出分类

我们即将讨论的最后一种分类任务被叫做“多输出-多类分类”（或者简称为多输出分类）。它是多标签分类的简单泛化，在这里每一个标签可以是多类别的（比如说，它可以有多于两个可能值）。

为了说明这点，我们建立一个系统，它可以去除图片当中的噪音。它将一张混有噪音的图片作为输入，期待它输出一张干净的数字图片，用一个像素强度的数组表示，就像 MNIST 图片那样。注意到这个分类器的输出是多标签的（一个像素一个标签）和每个标签可以有多个值（像素强度取值范围从 0 到 255）。所以它是一个多输出分类系统的例子。

分类与回归之间的界限是模糊的，比如这个例子。按理说，预测一个像素的强度更类似于一个回归任务，而不是一个分类任务。而且，多输出系统不限于分类任务。你甚至可以让你一个系统给每一个样例都输出多个标签，包括类标签和值标签。

让我们从 MNIST 的图片创建训练集和测试集开始，然后给图片的像素强度添加噪声，这里是用 NumPy 的randint()函数。目标图像是原始图像。

1
2
3
4
5
6

noise = rnd.randint(0, 100, (len(X_train), 784))
noise = rnd.randint(0, 100, (len(X_test), 784))
X_train_mod = X_train + noise
X_test_mod = X_test + noise
y_train_mod = X_train
y_test_mod = X_test

让我们看一下测试集当中的一张图片（是的，我们在窥探测试集，所以你应该马上邹眉）：

左边的加噪声的输入图片。右边是干净的目标图片。现在我们训练分类器，让它清洁这张图片：

1
2
3

knn_clf.fit(X_train_mod, y_train_mod)
clean_digit = knn_clf.predict([X_test_mod[some_index]])
plot_digit(clean_digit)

到这里就讲完分类的内容了，有点混乱对不对，我们来总结梳理一下。

• 要掌握自定义k折交叉验证的方法（≈cross_val_score）

• cross_val_score为验证模型的一个好方法，但是只能得到准确率的评估分数

• 如果正反例数据偏差大，我们需要用到混淆矩阵，这个矩阵要用到预测值而不是评估分数，所以改cross_val_predict，这会返回每个测试折做出的预测值，即y_train_pred

• 利用预测值y_train_pred可以得到混淆矩阵，精确率，召回率，F1

• 有时我们需要阈值来平衡精确率，召回率，而Scikit-Learn 不让你直接设置阈值，它会调用decision_function()方法。返回样例的分数值，然后基于这个分数值，使用你想要的任何阈值做出预测。

  1 2 3 4 5 6

  >>> y_scores = sgd_clf.decision_function([some_digit]) >>> y_scores array([ 161855.74572176]) >>> threshold = 0 >>> y_some_digit_pred = (y_scores > threshold) array([ True], dtype=bool)

• 每次都设定阈值不是一个完美的方法，如何才能找到合适的阈值呢？你需要再次使用cross_val_predict()得到每一个样例的分数值，但是这一次指定返回一个决策分数，而不是预测值。(阈值相关，就要进行打分)

  1 2	y_scores = cross_val_predict(sgd_clf, X_train, y_train_5, cv=3, method="decision_function")

  现在有了这些分数值。对于任何可能的阈值，使用precision_recall_curve(),你都可以计算精确率和召回率；precisions, recalls, thresholds是任何阈值的范围值，可以变化曲线和PR曲线

  1 2 3

  from sklearn.metrics import precision_recall_curve precisions, recalls, thresholds = precision_recall_curve(y_train_5, y_scores)

• 与PR曲线另一个相关的是ROC曲线（TPR/FPR），为了画出 ROC 曲线，你首先需要计算各种不同阈值下的 TPR、FPR，使用roc_curve()函数（还是要打分）；跳过ROC曲线(其实相当于已经做了)，想直接计算出ROC AUC也行。

  1 2 3 4 5

  from sklearn.metrics import roc_curve fpr, tpr, thresholds = roc_curve(y_train_5, y_scores) from sklearn.metrics import roc_auc_score roc_auc_score(y_train_5, y_scores)

  如果想得到RandomForestClassifier的ROC曲线，由于RandomForestClassifier不提供decision_function()方法，相反，它提供了predict_proba()方法（另外一种概率打分），返回概率值，此时用正例概率作为分值。例如70%的概率是垃圾邮件。

  另外，因为 ROC 曲线跟准确率/召回率曲线（或者叫 PR）很类似，你或许会好奇如何决定使用哪一个曲线呢？一个笨拙的规则是，优先使用 PR 曲线当正例很少，或者当你关注假正例多于假反例的时候。其他情况使用 ROC 曲线

• 多类别分类有一对一ovo, 一对多ova两种方法，一般svm由于在训练集的大小上很难扩展，因为它可以在小的数据集上面可以更多地训练，故用ovo，其他大部分用ova。如果Scikit-Learn嗅探出你想做一个多分类任务，它会自动使用ova，svm训练器除外

• 误差分析，将混淆矩阵归一化后用图片色块输出，查看哪些类别经常被错误分类。

写在前面：这个系列打算把「Hands-On Machine Learning with Scikit-Learn and TensorFlow 」重新梳理一遍，这本书在看完机器学习基础知识之后有一个很好的算法实践，对于算法落地有很多帮助。这次写的Sklearn 与 TensorFlow 机器学习实用指南系列，目的是让自己更清楚算法的每个流程处理，加强对一些机器学习模型理解。这本书在github有中文的翻译版本（还在更新）.

拆分数据集

训练集+测试集

1
2
3
4
5
6
7
8
9
10

import numpy as np

def split_train_test(data, test_ratio):
    shuffled_indices = np.random.permutation(len(data))    # 打乱序列
    test_set_size = int(len(data) * test_ratio)    # 拆分比例
    test_indices = shuffled_indices[:test_set_size]
    train_indices = shuffled_indices[test_set_size:]
    return data.iloc[train_indices], data.iloc[test_indices]
    
train_set, test_set = split_train_test(housing, 0.2)    # housing数据二八拆分

或者直接将整体数据打乱，然后按需取量。(california_housing_dataframe为谷歌机器学习教程提供的加州住房数据)

1
2
3
4

california_housing_dataframe = california_housing_dataframe.reindex(	# 整体打乱
    np.random.permutation(california_housing_dataframe.index))
train_set = california_housing_dataframe.head(12000)
test_set = california_housing_dataframe.tail(5000)

以上为训练集+测试集的拆分方式

训练集+验证集+测试集

这样的拆分方式主要有存在一些不足。1、程序多次运行后，测试集的数据有可能会加入到训练集当中，调参时用于改进模型超参数的测试集会造成过拟合。2、不便于新数据的加入

更好的办法是将数据集拆分为训练集+验证集+测试集。

那如何解决新加入数据的问题呢？一个通常的解决办法是使用每个实例的识别码，以判定是否这个实例是否应该放入测试集（假设实例有单一且不变的识别码）。例如，你可以计算出每个实例识别码的哈希值，只保留其最后一个字节，如果值小于等于 51（约为 256 的 20%），就将其放入测试集。这样可以保证在多次运行中，测试集保持不变，即使更新了数据集。新的测试集会包含新实例中的 20%，但不会有之前位于训练集的实例。可能很多数据没有稳定的特征，最简单的办法就是利用索引作为识别码。下面的代码根据识别码按0.7,0.2,0.1比例拆分训练集、验证集和测试集。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# 参数identifier为单一且不变的识别码，可以为索引id
# hash(np.int64(identifier)).digest()[-1]返回识别码的哈希摘要值的最后一个字节
def validate_set_check(identifier, validate_ratio, test_ratio, hash):
    return 256 * test_ratio <= hash(np.int64(identifier)).digest()[-1] < 256 *      (validate_ratio+test_ratio)  

def test_set_check(identifier, test_ratio, hash):
    return hash(np.int64(identifier)).digest()[-1] < 256 * test_ratio    # 记录满足条件的索引

def split_train_test_by_id(data, validate_ratio, test_ratio, id_column, hash=hashlib.md5):
    ids = data[id_column]	# 确定识别码
    in_validate_set = ids.apply(lambda id_: validate_set_check(id_, validate_ratio, test_ratio，hash))
    in_test_set = ids.apply(lambda id_: test_set_check(id_, test_ratio, hash))
    combine_set = np.bitwise_or(in_validate_set, in_test_set)
    return data.loc[~combine_set], data.loc[in_validate_set], data.loc[in_test_set]			                       
housing_with_id = housing.reset_index()   # housing数据增加一个索引列，放在数据的第一列
train_set, validate_set, test_set = split_train_test_by_id(housing_with_id, 0.2, 0.1, "index")

分成采样

另外一种拆分方式：分成采样

将人群分成均匀的子分组，称为分层，从每个分层去除合适数量的实例，以保证测试集对总人数有代表性。例如，美国人口的 51.3% 是女性，48.7% 是男性。所以在美国，严谨的调查需要保证样本也是这个比例：513 名女性，487 名男性作为数据样本。数据集中的每个分层都要有足够的实例位于你的数据中，这点很重要。否则，对分层重要性的评估就会有偏差。这意味着，你不能有过多的分层，且每个分层都要足够大。后面的代码通过将收入中位数除以 1.5（以限制收入分类的数量），创建了一个收入类别属性，用ceil对值舍入（以产生离散的分类），然后将所有大于 5的分类归入到分类5 ：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# 预处理，创建"income_cat"属性 
# 凡是会对原数组作出修改并返回一个新数组的，往往都有一个 inplace可选参数
# inplace=True,原数组名对应的内存值直接改变;inplace=False,原数组名对应的内存值并不改变，新的结果赋给一个新的数组.
housing["income_cat"] = np.ceil(housing["median_income"] / 1.5)
housing["income_cat"].where(housing["income_cat"] < 5, 5.0, inplace=True)

# 现在，就可以根据收入分类，进行分层采样。你可以使用 Scikit-Learn 的StratifiedShuffleSplit类
from sklearn.model_selection import StratifiedShuffleSplit

# random_state为随机种子生成器，可以得到相同的随机结果
# n_splits是将训练数据分成train/test对的组数，这里汇总成一组数据
split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)    

for train_index, test_index in split.split(housing, housing["income_cat"]):
    strat_train_set = housing.loc[train_index]
    strat_test_set = housing.loc[test_index]

# 现在，你需要删除income_cat属性，使数据回到初始状态：    
for set in (strat_train_set, strat_test_set):
    set.drop(["income_cat"], axis=1, inplace=True)

数据预处理

将原始数据映射到特征

我们在进行机器学习的时候，采用的数据样本往往是矢量（特征矢量），而我们的原始数据并不是以矢量的形式呈现给我们的，这是便需要将数据映射到特征

整数和浮点数映射

直接映射便ok（虽然机器学习是根据浮点值进行的训练，但是不需要将整数6转换为6.0，这个过程是默认的）

字符串映射

好多时候，有的特征是字符串，比如此前训练的加利福尼亚房产数据集中的街区名称，机器学习是无法根据字符串来学习规律的，所以需要转换。但是存在一个问题，如果字符特征是’’一环’’ ‘’二环’’ ‘’三环’’…（代表某个城市的地理位置），那么对其进行数值转换的时候，是不可以编码为形如1，2，3，4…这样的数据的，因为其存在数据大小的问题，学习模型会把他们的大小关系作为特征而学习，所以我们需要引入独热编码,（具体解释见链接，解释的很好）.我们需要把这些文本标签转换为数字。Scikit-Learn 为这个任务提供了一个转换器LabelEncoder：

1
2
3
4
5
6
7
8
9
10
11

# 简单来说 LabelEncoder 是对不连续的数字或者文本进行编号
# le.fit([1,5,67,100])
# le.transform([1,1,100,67,5])
# 输出： array([0,0,3,2,1])

>>> from sklearn.preprocessing import LabelEncoder
>>> encoder = LabelEncoder()
>>> housing_cat = housing["ocean_proximity"]
>>> housing_cat_encoded = encoder.fit_transform(housing_cat)	# 装换器
>>> housing_cat_encoded
array([1, 1, 4, ..., 1, 0, 3])

译注:

在原书中使用LabelEncoder转换器来转换文本特征列的方式是错误的，该转换器只能用来转换标签（正如其名）。在这里使用LabelEncoder没有出错的原因是该数据只有一列文本特征值，在有多个文本特征列的时候就会出错。应使用factorize()方法来进行操作：

1 2	housing_cat_encoded, housing_categories = housing_cat.factorize() housing_cat_encoded[:10]

处理离散特征这还不够，Scikit-Learn 提供了一个编码器OneHotEncoder，用于将整书分类值转变为独热向量。注意fit_transform()用于 2D 数组，而housing_cat_encoded是一个 1D 数组，所以需要将其变形：

1
2
3
4
5
6
7
8

# reshape(-1,1)里面的-1代表将数据自动计算有多少行，但是列数明确设置为1
# reshape(-1)则是变形为1行和自动计算有多少列
>>> from sklearn.preprocessing import OneHotEncoder
>>> encoder = OneHotEncoder()
>>> housing_cat_1hot = encoder.fit_transform(housing_cat_encoded.reshape(-1,1))
>>> housing_cat_1hot
<16513x5 sparse matrix of type '<class 'numpy.float64'>'
    with 16513 stored elements in Compressed Sparse Row format>

注意输出结果是一个 SciPy 稀疏矩阵，而不是 NumPy 数组。当类别属性有数千个分类时，这样非常有用。经过独热编码，我们得到了一个有数千列的矩阵，这个矩阵每行只有一个 1，其余都是 0。使用大量内存来存储这些 0 非常浪费，所以稀疏矩阵只存储非零元素的位置。你可以像一个 2D 数据那样进行使用，但是如果你真的想将其转变成一个（密集的）NumPy 数组，只需调用toarray()方法：

1
2
3
4
5
6
7
8

>>> housing_cat_1hot.toarray()
array([[ 0.,  1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.],
       ...,
       [ 0.,  1.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.]])

使用类LabelBinarizer，我们可以用一步执行这两个转换（从文本分类到整数分类，再从整数分类到独热向量）：

1
2
3
4
5
6
7
8
9
10
11

>>> from sklearn.preprocessing import LabelBinarizer
>>> encoder = LabelBinarizer()
>>> housing_cat_1hot = encoder.fit_transform(housing_cat)
>>> housing_cat_1hot
array([[0, 1, 0, 0, 0],
       [0, 1, 0, 0, 0],
       [0, 0, 0, 0, 1],
       ...,
       [0, 1, 0, 0, 0],
       [1, 0, 0, 0, 0],
       [0, 0, 0, 1, 0]])

注意默认返回的结果是一个密集 NumPy 数组。向构造器LabelBinarizer传递sparse_output=True，就可以得到一个稀疏矩阵。

译注:

在原书中使用LabelBinarizer的方式也是错误的，该类也应用于标签列的转换。正确做法是使用sklearn即将提供的CategoricalEncoder类。如果在你阅读此文时sklearn中尚未提供此类，用如下方式代替：（来自Pull Request #9151）

1 2 3 4 5 6

#from sklearn.preprocessing import CategoricalEncoder # in future versions of Scikit-Learn cat_encoder = CategoricalEncoder() housing_cat_reshaped = housing_cat.values.reshape(-1, 1) housing_cat_1hot = cat_encoder.fit_transform(housing_cat_reshaped) housing_cat_1hot

寻找良好特征（的特点）

当得到特征之后，还是要进行筛选的，因为有的特征没有参考价值，就像我们的在做合成特征的时候，正常的特征数据是人均几间房间，而有的人是几十间，这明显没有参考价值
良好特征的几点原则

• 避免很少使用的离散特征值：如果只是出现了一两次的特征几乎是没有意义的

• 最好具有清晰明确的含义：特征的含义不仅仅是让机器学习的模型学习的，人也要知道其具体的含义，不然不利于分析数据（最好将数值很大的秒转换为天数，或者年，让人看起来直观一些）

• 将“神奇”的值与实际数据混为一谈：有些特征中会出现一些”神奇的数据”，当然这些数据并不是很少的特征，而是超出范围的异常值，比如特征应该是介于0——1之间的，但是因为这个数据是空缺的，而采用的默认数值-1，那么这样的数值就是”神奇”，解决办法是，将该特征转换为两个特征：

  • 一个特征只存储质正常范围的值，不含神奇值。
  • 一个特征存储布尔值，表示的信息为是否为空

• 考虑上游不稳定性：由经验可知，特征的定义不应随时间发生变化，代表城市名称的话，那么特征值始终都该是城市的名称，但是有的时候，上游模型将特征值处理完毕后，返还给下游模型的却变成了数值，这样是不好的，因为这种表示在未来运行其他模型时可能轻易发生变化，那么特征就乱套了

  ​

可视化数据寻找规律：

1
2
3
4
5

housing.plot(kind="scatter", x="longitude", y="latitude", alpha=0.4,
    s=housing["population"]/100, label="population",
    c="median_house_value", cmap=plt.get_cmap("jet"), colorbar=True,
)
plt.legend()

每个圈的半径表示街区的人口（选项s），颜色代表价格（选项c）。我们用预先定义的名为jet的颜色图（选项cmap），它的范围是从蓝色（低价）到红色（高价）：

皮尔逊相关系数

1

corr_matrix = housing.corr()

Pandas 的scatter_matrix函数

另一种检测属性间相关系数的方法是使用 Pandas 的scatter_matrix函数,它能画出每个数值属性对每个其它数值属性的图。因为现在共有 11 个数值属性，你可以得到11 ** 2 = 121张图。

1
2
3
4
5

from pandas.tools.plotting import scatter_matrix

attributes = ["median_house_value", "median_income", "total_rooms",
              "housing_median_age"]
scatter_matrix(housing[attributes], figsize=(12, 8))

得到两个属性的散点图

清查数据

截至目前，我们假定用于训练和测试的所有数据都是值得信赖的。在现实生活中，数据集中的很多样本是不可靠的，原因有以下一种或多种：

• 遗漏值。 例如，有人忘记为某个房屋的年龄输入值。(值会为-1，所以要分为两个特征，忘了的看上面)
• 重复样本。 例如，服务器错误地将同一条记录上传了两次。
• 不良标签。 例如，有人错误地将一颗橡树的图片标记为枫树。
• 不良特征值。 例如，有人输入了多余的位数，或者温度计被遗落在太阳底下。

一旦检测到存在这些问题，通常需要将相应样本从数据集中移除，从而“修正”不良样本。要检测遗漏值或重复样本，可以编写一个简单的程序。检测不良特征值或标签可能会比较棘手，可采用可视化数据的方法。

对于处理特征丢失的问题。前面，你应该注意到了属性total_bedrooms有一些缺失值。有三个解决选项：

• 去掉对应的街区；（数据大可用）
• 去掉整个属性；
• 进行赋值（0、平均值、中位数等等）。

用DataFrame的dropna()，drop()，和fillna()方法，可以方便地实现：

1
2
3
4

housing.dropna(subset=["total_bedrooms"])    # 选项1
housing.drop("total_bedrooms", axis=1)       # 选项2    axis=0对行操作，axis=1对列操作
median = housing["total_bedrooms"].median()
housing["total_bedrooms"].fillna(median)     # 选项3

如果选择选项 3，你需要计算训练集的中位数，用中位数填充训练集的缺失值，不要忘记保存该中位数。后面用测试集评估系统时，需要替换测试集中的缺失值，也可以用来实时替换新数据中的缺失值。

Scikit-Learn 提供了一个方便的类来处理缺失值：Imputer。下面是其使用方法：首先，需要创建一个Imputer实例，指定用该属性的中位数替换它的每个缺失值：

1
2
3

from sklearn.preprocessing import Imputer

imputer = Imputer(strategy="median")    # 进行中位数赋值

因为只有数值属性才能算出中位数，我们需要创建一份不包括文本属性ocean_proximity的数据副本：

1

housing_num = housing.drop("ocean_proximity", axis=1) # 去除ocean_proximity不为数值属性的特征

现在，就可以用fit()方法将imputer实例拟合到训练数据：

1

imputer.fit(housing_num)

imputer计算出了每个属性的中位数，并将结果保存在了实例变量statistics_中。只有属性total_bedrooms有缺失值，但是我们确保一旦系统运行起来，新的数据中没有缺失值，所以安全的做法是将imputer应用到每个数值：

1
2
3
4

>>> imputer.statistics_    # 实例变量statistics_和housing_num数值数据得到的中位数是一样的
array([ -118.51 , 34.26 , 29. , 2119. , 433. , 1164. , 408. , 3.5414])
>>> housing_num.median().values
array([ -118.51 , 34.26 , 29. , 2119. , 433. , 1164. , 408. , 3.5414])

现在，你就可以使用这个“训练过的”imputer来对训练集进行转换，通过将缺失值替换为中位数：

1

X = imputer.transform(housing_num)

结果是一个普通的 Numpy 数组，包含有转换后的特征。如果你想将其放回到 PandasDataFrame中，也很简单：

1

housing_tr = pd.DataFrame(X, columns=housing_num.columns) # 得到处理缺失值后的DF数据

整理数据：

数据缩放

有两种常见的方法可以让所有的属性有相同的量度：线性函数归一化（Min-Max scaling）和标准化（standardization）。Scikit-Learn 提供了一个转换器MinMaxScaler来实现这个功能。它有一个超参数feature_range，可以让你改变范围，如果不希望范围是 0 到 1；Scikit-Learn 提供了一个转换器StandardScaler来进行标准化

min-max方式,对应的方法为

1

MinMaxScaler(self, feature_range=(0, 1), copy=True)

standardization 标准化数据,对应的方法为

1

StandardScaler(self, copy=True, with_mean=True, with_std=True)

警告：与所有的转换一样，缩放器只能向训练集拟合，而不是向完整的数据集（包括测试集）。只有这样，你才能用缩放器转换训练集和测试集（和新数据）。

处理极端离群值

还是举加利福尼亚州住房数据集中的人均住房数的例子，有的极端值达到了50
对于这些极端值其实很好处理，无非几个办法

• 对数缩放

1

roomsPerPerson = log((totalRooms / population) + 1)

• 特征值限制到 某个上限或者下限

1

roomsPerPerson = min(totalRooms / population, 4)	# 大于4.0的取4.0

分箱

分箱其实是一个形象化的说法，就是把数据分开来，装在一个个箱子里，这样一个箱子里的数据就是一家人了。
那有什么用呢？下面就举个栗子！

在数据集中，latitude 是一个浮点值。不过，在我们的模型中将 latitude 表示为浮点特征没有意义。这是因为纬度和房屋价值之间不存在线性关系。例如，纬度 35 处的房屋并不比纬度 34 处的房屋贵 35/34（或更便宜）。但是，纬度或许能很好地预测房屋价值。

我们现在拥有 11 个不同的布尔值特征（LatitudeBin1、LatitudeBin2、…、LatitudeBin11），而不是一个浮点特征。拥有 11 个不同的特征有点不方便，因此我们将它们统一成一个 11 元素矢量。这样做之后，我们可以将纬度 37.4 表示为：

[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]

分箱之后，我们的模型现在可以为每个纬度学习完全不同的权重。（是不是觉得有点像独热编码，没错，就是的）

为了简单起见，我们在纬度样本中使用整数作为分箱边界。如果我们需要更精细的解决方案，我们可以每隔 1/10 个纬度拆分一次分箱边界。添加更多箱可让模型从纬度 37.4 处学习和维度 37.5 处不一样的行为，但前提是每 1/10 个纬度均有充足的样本可供学习。

另一种方法是按分位数分箱，这种方法可以确保每个桶内的样本数量是相等的。按分位数分箱完全无需担心离群值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

'''
分桶也称为分箱。
例如，我们可以将 population 分为以下 3 个分桶：
bucket_0 (< 5000)：对应于人口分布较少的街区
bucket_1 (5000 - 25000)：对应于人口分布适中的街区
bucket_2 (> 25000)：对应于人口分布较多的街区
根据前面的分桶定义，以下 population 矢量：
[[10001], [42004], [2500], [18000]]
将变成以下经过分桶的特征矢量：
[[1], [2], [0], [1]]
这些特征值现在是分桶索引。请注意，这些索引被视为离散特征。通常情况下，这些特征将被进一步转换为上述独热表示法，但这是以透明方式实现的。

要为分桶特征定义特征列，我们可以使用 bucketized_column（而不是使用 numeric_column），该列将数字列作为输入，并使用 boundardies 参数中指定的分桶边界将其转换为分桶特征。以下代码为 households 和 longitude 定义了分桶特征列；get_quantile_based_boundaries 函数会根据分位数计算边界，以便每个分桶包含相同数量的元素
'''
def get_quantile_based_boundaries(feature_values, num_buckets):
    boundaries = np.arange(1.0, num_buckets) / num_buckets
    quantiles = feature_values.quantile(boundaries)
    return [quantiles[q] for q in quantiles.keys()]

# Divide households into 7 buckets.
households = tf.feature_column.numeric_column("households")		# 定义数值特征
# 分桶特征bucketized_column第一个参数用数字列 numeric_column得到的households，第二个参数用上面get_quantile_based_boundaries方法得到的分桶数据，返回的bucketized_households为可使用的分桶特征
bucketized_households = tf.feature_column.bucketized_column(       
    households,boundaries=get_quantile_based_boundaries(california_housing_dataframe["households"], 7))

自定义转换器

尽管 Scikit-Learn 提供了许多有用的转换器，你还是需要自己动手写转换器执行任务，比如自定义的清理操作，或属性组合。你需要让自制的转换器与 Scikit-Learn 组件（比如流水线）无缝衔接工作，因为 Scikit-Learn 是依赖鸭子类型的（而不是继承，忽略对象，只要行为像就行），你所需要做的是创建一个类并执行三个方法：fit()（返回self），transform()，和fit_transform()。通过添加TransformerMixin作为基类，可以很容易地得到最后一个。另外，如果你添加BaseEstimator作为基类（且构造器中避免使用args和kargs），你就能得到两个额外的方法（get_params()和set_params()），二者可以方便地进行超参数自动微调。例如，一个小转换器类添加了上面讨论的属性：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

# 添加一个特征组合的装换器
from sklearn.base import BaseEstimator, TransformerMixin
rooms_ix, bedrooms_ix, population_ix, household_ix = 3, 4, 5, 6

# 这里的示例没有定义fit_transform()，可能是因为fit()没有做任何动作（我猜的
class CombinedAttributesAdder(BaseEstimator, TransformerMixin):
    def __init__(self, add_bedrooms_per_room = True): # no *args or **kargs
        self.add_bedrooms_per_room = add_bedrooms_per_room
    def fit(self, X, y=None):
        return self  # nothing else to do
    def transform(self, X, y=None):
        rooms_per_household = X[:, rooms_ix] / X[:, household_ix]  # X[:,3]表示的是第4列所有数据
        population_per_household = X[:, population_ix] / X[:, household_ix]
        if self.add_bedrooms_per_room:
            bedrooms_per_room = X[:, bedrooms_ix] / X[:, rooms_ix]
            return np.c_[X, rooms_per_household, population_per_household, # np.c_表示的是拼接数组。
                         bedrooms_per_room]
        else:
            return np.c_[X, rooms_per_household, population_per_household]

attr_adder = CombinedAttributesAdder(add_bedrooms_per_room=False)
housing_extra_attribs = attr_adder.transform(housing.values)    # 返回一个加入新特征的数据

在这个例子中，转换器有一个超参数add_bedrooms_per_room，默认设为True（提供一个合理的默认值很有帮助）。这个超参数可以让你方便地发现添加了这个属性是否对机器学习算法有帮助。更一般地，你可以为每个不能完全确保的数据准备步骤添加一个超参数。数据准备步骤越自动化，可以自动化的操作组合就越多，越容易发现更好用的组合（并能节省大量时间）。

另外sklearn是不能直接处理DataFrames的，那么我们需要自定义一个处理的方法将之转化为numpy类型

1
2
3
4
5
6
7

class DataFrameSelector(BaseEstimator,TransformerMixin):
    def __init__(self,attribute_names): #可以为列表
        self.attribute_names = attribute_names
    def fit(self,X,y=None):
        return self
    def transform(self,X):
        return X[self.attribute_names].values #返回的为numpy array

转换流水线

目前在数据预处理阶段，我们需要对缺失值进行处理、特征组合和特征缩放。每一步的执行都有着先后顺序，存在许多数据转换步骤，需要按一定的顺序执行。sklearn提供了Pipeline帮助顺序完成转换幸运的是，Scikit-Learn 提供了类Pipeline，来进行这一系列的转换。下面是一个数值属性的小流水线：

1
2
3
4
5
6
7
8
9
10

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

num_pipeline = Pipeline([
        ('imputer', Imputer(strategy="median")),    # 处理缺失值
        ('attribs_adder', CombinedAttributesAdder()),    # 特征组合
        ('std_scaler', StandardScaler()),    # 特征缩放
        ])

housing_num_tr = num_pipeline.fit_transform(housing_num)

Pipeline构造器需要一个定义步骤顺序的名字/估计器对的列表。除了最后一个估计器，其余都要是转换器（即，它们都要有fit_transform()方法）。名字可以随意起。

当你调用流水线的fit()方法，就会对所有转换器顺序调用fit_transform()方法，将每次调用的输出作为参数传递给下一个调用，一直到最后一个估计器，它只执行fit()方法。

估计器（Estimator）：很多时候可以直接理解成分类器，主要包含两个函数：fit()和predict()
转换器（Transformer）：转换器用于数据预处理和数据转换，主要是三个方法：fit（）,transform()和fit_transform()

最后的估计器是一个StandardScaler，它是一个转换器，因此这个流水线有一个transform()方法，可以顺序对数据做所有转换（它还有一个fit_transform方法可以使用，就不必先调用fit()再进行transform()）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

num_attribs = list(housing_num) # 返回的为列名[col1,col2,....]
cat_attribs = ["ocean_proximity"]

num_pipeline = Pipeline([ # 数值类型
        ('selector', DataFrameSelector(num_attribs)),	# DataFrames转为numpy array
        ('imputer', Imputer(strategy="median")),	# 缺失值处理
        ('attribs_adder', CombinedAttributesAdder()),	# 特征组合
        ('std_scaler', StandardScaler()),	# 缩放
    ])

cat_pipeline = Pipeline([ # 标签类型
        ('selector', DataFrameSelector(cat_attribs)), 	# DataFrames转为numpy array
        ('cat_encoder', CategoricalEncoder(encoding="onehot-dense")),	
    ])

上面定义的为分别处理数值类型和标签类型的转换流程，housing_num为DataFrame类型，list(DataFrame)的结果返回的为列名字。上面着两个流程还可以再整合一起。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

from sklearn.pipeline import FeatureUnion
full_pipeline = FeatureUnion(transformer_list=[
        ("num_pipeline", num_pipeline),
        ("cat_pipeline", cat_pipeline),
    ])
housing_prepared = full_pipeline.fit_transform(housing) # 最终的结果
>>> housing_prepared
array([[ 0.73225807, -0.67331551,  0.58426443, ...,  0.        ,
         0.        ,  0.        ],
       [-0.99102923,  1.63234656, -0.92655887, ...,  0.        ,
         0.        ,  0.        ],
       [...]
>>> housing_prepared.shape
(16513, 17)

译注:

如果你在上面代码中的cat_pipeline流水线使用LabelBinarizer转换器会导致执行错误，解决方案是用上文提到的CategoricalEncoder转换器来代替：

1 2 3 4

cat_pipeline = Pipeline([ ('selector', DataFrameSelector(cat_attribs)), ('cat_encoder', CategoricalEncoder(encoding="onehot-dense")), ])

每个子流水线都以一个选择转换器开始：通过选择对应的属性（数值或分类）、丢弃其它的，来转换数据，并将输出DataFrame转变成一个 NumPy 数组。Scikit-Learn 没有工具来处理 PandasDataFrame，因此我们需要写一个简单的自定义转换器来做这项工作：

1
2
3
4
5
6
7
8
9

from sklearn.base import BaseEstimator, TransformerMixin

class DataFrameSelector(BaseEstimator, TransformerMixin):
    def __init__(self, attribute_names):
        self.attribute_names = attribute_names
    def fit(self, X, y=None):
        return self
    def transform(self, X):
        return X[self.attribute_names].values

训练模型

我们先来训练一个线性回归模型：

1
2
3
4

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels)    # 利用预处理好的数据进行训练模型

完毕！你现在就有了一个可用的线性回归模型。用一些训练集中的实例做下验证：

1
2
3
4
5
6
7

>>> some_data = housing.iloc[:5]    # 前五个作为预测数据
>>> some_labels = housing_labels.iloc[:5]
>>> some_data_prepared = full_pipeline.transform(some_data)
>>> print("Predictions:\t", lin_reg.predict(some_data_prepared))    # 预测结果
Predictions:     [ 303104.   44800.  308928.  294208.  368704.]
>>> print("Labels:\t\t", list(some_labels))
Labels:         [359400.0, 69700.0, 302100.0, 301300.0, 351900.0]    # 实际结果

行的通，尽管预测并不怎么准确（比如，第二个预测偏离了 50%！）。让我们使用 Scikit-Learn 的mean_squared_error函数，用全部训练集来计算下这个回归模型的 RMSE：

1
2
3
4
5
6

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

OK，有总比没有强，但显然结果并不好，这是一个模型欠拟合训练数据的例子。当这种情况发生时，意味着特征没有提供足够多的信息来做出一个好的预测，或者模型并不强大，修复欠拟合的主要方法是选择一个更强大的模型，给训练算法提供更好的特征，或去掉模型上的限制，你可以尝试添加更多特征（比如，人口的对数值），但是首先让我们尝试一个更为复杂的模型，看看效果。训练一个决策树模型DecisionTreeRegressor。这是一个强大的模型，可以发现数据中复杂的非线性关系。

1
2
3
4
5
6
7
8
9

from sklearn.tree import DecisionTreeRegressor

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

等一下，发生了什么？没有误差？这个模型可能是绝对完美的吗？当然，更大可能性是这个模型严重过拟合数据。如何确定呢？如前所述，直到你准备运行一个具备足够信心的模型，都不要碰测试集，因此你需要使用训练集的部分数据来做训练，用一部分来做模型验证。

用交叉验证做更佳的评估

使用 Scikit-Learn 的交叉验证功能。下面的代码采用了 K 折交叉验证（K-fold cross-validation）：它随机地将训练集分成十个不同的子集，成为“折”，然后训练评估决策树模型 10 次，每次选一个不用的折来做评估，用其它 9 个来做训练。结果是一个包含 10 个评分的数组：

1
2
3
4

from sklearn.model_selection import cross_val_score
scores = cross_val_score(tree_reg, housing_prepared, housing_labels,
                         scoring="neg_mean_squared_error", cv=10)
rmse_scores = np.sqrt(-scores)

警告：Scikit-Learn 交叉验证功能期望的是效用函数（越大越好）而不是损失函数（越低越好），因此得分函数实际上与 MSE 相反（即负值），这就是为什么前面的代码在计算平方根之前先计算-scores。

来看下结果

1
2
3
4
5
6
7
8
9
10
11

>>> def display_scores(scores):
...     print("Scores:", scores)
...     print("Mean:", scores.mean())
...     print("Standard deviation:", scores.std())
...
>>> display_scores(tree_rmse_scores)
Scores: [ 74678.4916885   64766.2398337   69632.86942005  69166.67693232
          71486.76507766  73321.65695983  71860.04741226  71086.32691692
          76934.2726093   69060.93319262]
Mean: 71199.4280043
Standard deviation: 3202.70522793

现在决策树就不像前面看起来那么好了。实际上，它看起来比线性回归模型还糟！注意到交叉验证不仅可以让你得到模型性能的评估，还能测量评估的准确性（即，它的标准差）。决策树的评分大约是 71200，通常波动有 ±3200。如果只有一个验证集，就得不到这些信息。但是交叉验证的代价是训练了模型多次，不可能总是这样。

让我们计算下线性回归模型的的相同分数，以做确保：

1
2
3
4
5
6
7
8
9
10

>>> lin_scores = cross_val_score(lin_reg, housing_prepared, housing_labels,
...                              scoring="neg_mean_squared_error", cv=10)
...
>>> lin_rmse_scores = np.sqrt(-lin_scores)
>>> display_scores(lin_rmse_scores)
Scores: [ 70423.5893262   65804.84913139  66620.84314068  72510.11362141
          66414.74423281  71958.89083606  67624.90198297  67825.36117664
          72512.36533141  68028.11688067]
Mean: 68972.377566
Standard deviation: 2493.98819069

判断没错：决策树模型过拟合很严重，它的性能比线性回归模型还差

现在再尝试最后一个模型：RandomForestRegressor（随机森林），随机森林是通过用特征的随机子集训练许多决策树。在其它多个模型之上建立模型成为集成学习（Ensemble Learning），它是推进 ML 算法的一种好方法。我们会跳过大部分的代码，因为代码本质上和其它模型一样：

1
2
3
4
5
6
7
8
9
10
11
12

>>> from sklearn.ensemble import RandomForestRegressor
>>> forest_reg = RandomForestRegressor()
>>> forest_reg.fit(housing_prepared, housing_labels)
>>> [...]
>>> forest_rmse
22542.396440343684
>>> display_scores(forest_rmse_scores)
Scores: [ 53789.2879722   50256.19806622  52521.55342602  53237.44937943
          52428.82176158  55854.61222549  52158.02291609  50093.66125649
          53240.80406125  52761.50852822]
Mean: 52634.1919593
Standard deviation: 1576.20472269

现在好多了：随机森林看起来很有希望。但是，训练集的评分仍然比验证集的评分低很多。解决过拟合可以通过简化模型，给模型加限制（即，正则化），或用更多的训练数据。在深入随机森林之前，你应该尝试下机器学习算法的其它类型模型（不同核心的支持向量机，神经网络，等等），不要在调节超参数上花费太多时间。目标是列出一个可能模型的列表（两到五个）。

提示：你要保存每个试验过的模型，以便后续可以再用。要确保有超参数和训练参数，以及交叉验证评分，和实际的预测值。这可以让你比较不同类型模型的评分，还可以比较误差种类。你可以用 Python 的模块pickle，非常方便地保存 Scikit-Learn 模型，或使用sklearn.externals.joblib，后者序列化大 NumPy 数组更有效率：

1
2
3
4
5

from sklearn.externals import joblib

joblib.dump(my_model, "my_model.pkl")
# 然后
my_model_loaded = joblib.load("my_model.pkl")

模型微调

网格搜索：使用 Scikit-Learn 的GridSearchCV来做这项搜索工作。你所需要做的是告诉GridSearchCV要试验有哪些超参数，要试验什么值，GridSearchCV就能用交叉验证试验所有可能超参数值的组合。例如，下面的代码搜索了RandomForestRegressor超参数值的最佳组合（很费时间）：

1
2
3
4
5
6
7
8
9
10
11
12
13

from sklearn.model_selection import GridSearchCV

param_grid = [
    {'n_estimators': [3, 10, 30], 'max_features': [2, 4, 6, 8]},
    {'bootstrap': [False], 'n_estimators': [3, 10], 'max_features': [2, 3, 4]},
  ]

forest_reg = RandomForestRegressor()

grid_search = GridSearchCV(forest_reg, param_grid, cv=5,
                           scoring='neg_mean_squared_error')

grid_search.fit(housing_prepared, housing_labels)

当你不能确定超参数该有什么值，一个简单的方法是尝试连续的 10 的幂（如果想要一个粒度更小的搜寻，可以用更小的数，就像在这个例子中对超参数n_estimators做的）。

param_grid告诉 Scikit-Learn 首先评估所有的列在第一个dict中的n_estimators和max_features的3 × 4 = 12种组合（不用担心这些超参数的含义，会在第 7 章中解释）。然后尝试第二个dict中超参数的2 × 3 = 6种组合，这次会将超参数bootstrap设为False而不是True（后者是该超参数的默认值）。完成后，你就能获得参数的最佳组合，如下所示：

1
2

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

你还能直接得到最佳的估计器：

1
2
3
4
5
6

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

当然，也可以得到评估得分：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

>>> cvres = grid_search.cv_results_
... for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
...     print(np.sqrt(-mean_score), params)
...
64912.0351358 {'max_features': 2, 'n_estimators': 3}
55535.2786524 {'max_features': 2, 'n_estimators': 10}
52940.2696165 {'max_features': 2, 'n_estimators': 30}
60384.0908354 {'max_features': 4, 'n_estimators': 3}
52709.9199934 {'max_features': 4, 'n_estimators': 10}
50503.5985321 {'max_features': 4, 'n_estimators': 30}
59058.1153485 {'max_features': 6, 'n_estimators': 3}
52172.0292957 {'max_features': 6, 'n_estimators': 10}
49958.9555932 {'max_features': 6, 'n_estimators': 30}
59122.260006 {'max_features': 8, 'n_estimators': 3}
52441.5896087 {'max_features': 8, 'n_estimators': 10}
50041.4899416 {'max_features': 8, 'n_estimators': 30}
62371.1221202 {'bootstrap': False, 'max_features': 2, 'n_estimators': 3}
54572.2557534 {'bootstrap': False, 'max_features': 2, 'n_estimators': 10}
59634.0533132 {'bootstrap': False, 'max_features': 3, 'n_estimators': 3}
52456.0883904 {'bootstrap': False, 'max_features': 3, 'n_estimators': 10}
58825.665239 {'bootstrap': False, 'max_features': 4, 'n_estimators': 3}
52012.9945396 {'bootstrap': False, 'max_features': 4, 'n_estimators': 10}

在这个例子中，我们通过设定超参数max_features为 6，n_estimators为 30，得到了最佳方案。对这个组合，RMSE 的值是 49959，这比之前使用默认的超参数的值（52634）要稍微好一些。祝贺你，你成功地微调了最佳模型！

随机搜索：当探索相对较少的组合时，就像前面的例子，网格搜索还可以。但是当超参数的搜索空间很大时，最好使用RandomizedSearchCV。这个类的使用方法和类GridSearchCV很相似，但它不是尝试所有可能的组合，而是通过选择每个超参数的一个随机值的特定数量的随机组合。这个方法有两个优点：

• 如果你让随机搜索运行，比如 1000 次，它会探索每个超参数的 1000 个不同的值（而不是像网格搜索那样，只搜索每个超参数的几个值）
• 你可以方便地通过设定搜索次数，控制超参数搜索的计算量。

分析最佳模型和它们的误差

通过分析最佳模型，常常可以获得对问题更深的了解。比如，RandomForestRegressor可以指出每个属性对于做出准确预测的相对重要性：

1
2
3
4
5
6
7
8

>>> feature_importances = grid_search.best_estimator_.feature_importances_
>>> feature_importances
array([  7.14156423e-02,   6.76139189e-02,   4.44260894e-02,
         1.66308583e-02,   1.66076861e-02,   1.82402545e-02,
         1.63458761e-02,   3.26497987e-01,   6.04365775e-02,
         1.13055290e-01,   7.79324766e-02,   1.12166442e-02,
         1.53344918e-01,   8.41308969e-05,   2.68483884e-03,
         3.46681181e-03])

将重要性分数和属性名放到一起：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

>>> extra_attribs = ["rooms_per_hhold", "pop_per_hhold", "bedrooms_per_room"]
>>> cat_one_hot_attribs = list(encoder.classes_)
>>> attributes = num_attribs + extra_attribs + cat_one_hot_attribs
>>> sorted(zip(feature_importances,attributes), reverse=True)
[(0.32649798665134971, 'median_income'),
 (0.15334491760305854, 'INLAND'),
 (0.11305529021187399, 'pop_per_hhold'),
 (0.07793247662544775, 'bedrooms_per_room'),
 (0.071415642259275158, 'longitude'),
 (0.067613918945568688, 'latitude'),
 (0.060436577499703222, 'rooms_per_hhold'),
 (0.04442608939578685, 'housing_median_age'),
 (0.018240254462909437, 'population'),
 (0.01663085833886218, 'total_rooms'),
 (0.016607686091288865, 'total_bedrooms'),
 (0.016345876147580776, 'households'),
 (0.011216644219017424, '<1H OCEAN'),
 (0.0034668118081117387, 'NEAR OCEAN'),
 (0.0026848388432755429, 'NEAR BAY'),
 (8.4130896890070617e-05, 'ISLAND')]

有了这个信息，你就可以丢弃一些不那么重要的特征（比如，显然只要一个分类ocean_proximity就够了，所以可以丢弃掉其它的）。你还应该看一下系统犯的误差，搞清为什么会有些误差，以及如何改正问题（添加更多的特征，或相反，去掉没有什么信息的特征，清洗异常值等等）。

模型评估

调节完系统之后，你终于有了一个性能足够好的系统。现在就可以用测试集评估最后的模型了。这个过程没有什么特殊的：从测试集得到预测值和标签，运行full_pipeline转换数据（调用transform()，而不是fit_transform()！），再用测试集评估最终模型：

1
2
3
4
5
6
7
8
9
10
11

final_model = grid_search.best_estimator_

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

X_test_prepared = full_pipeline.transform(X_test)

final_predictions = final_model.predict(X_test_prepared)

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

评估结果通常要比交叉验证的效果差一点，如果你之前做过很多超参数微调（因为你的系统在验证集上微调，得到了不错的性能，通常不会在未知的数据集上有同样好的效果）。这个例子不属于这种情况，但是当发生这种情况时，你一定要忍住不要调节超参数，使测试集的效果变好；这样的提升不能推广到新数据上。

兼容Python2+和Python3+

很多时候你可能希望你开发的程序能够同时兼容Python2+和Python3+。

试想你有一个非常出名的Python模块被很多开发者使用着，但并不是所有人都只使用Python2或者Python3。这时候你有两个办法。第一个办法是开发两个模块，针对Python2一个，针对Python3一个。还有一个办法就是调整你现在的代码使其同时兼容Python2和Python3。

本节中，我将介绍一些技巧，让你的脚本同时兼容Python2和Python3。

Future模块导入

第一种也是最重要的方法，就是导入__future__模块。它可以帮你在Python2中导入Python3的功能。这有一组例子：

上下文管理器是Python2.6+引入的新特性，如果你想在Python2.5中使用它可以这样做：

1

from __future__ import with_statement

在Python3中print已经变为一个函数。如果你想在Python2中使用它可以通过__future__导入：

1
2
3
4
5
6
7

print
# Output:

from __future__ import print_function
print(print)
# Output: <built-in function print>

模块重命名

首先，告诉我你是如何在你的脚本中导入模块的。大多时候我们会这样做：

1
2
3

import foo 
# or
from foo import bar

你知道么，其实你也可以这样做：

1

import foo as foo

这样做可以起到和上面代码同样的功能，但最重要的是它能让你的脚本同时兼容Python2和Python3。现在我们来看下面的代码：

1
2
3
4

try:
    import urllib.request as urllib_request  # for Python 3
except ImportError:
    import urllib2 as urllib_request  # for Python 2

过期的Python2内置功能

另一个需要了解的事情就是Python2中有12个内置功能在Python3中已经被移除了。要确保在Python2代码中不要出现这些功能来保证对Python3的兼容。这有一个强制让你放弃12内置功能的方法：

1

from future.builtins.disabled import *

现在，只要你尝试在Python3中使用这些被遗弃的模块时，就会抛出一个NameError异常如下：

1
2
3
4

from future.builtins.disabled import *

apply()
# Output: NameError: obsolete Python 2 builtin apply is disabled

标准库向下兼容的外部支持

有一些包在非官方的支持下为Python2提供了Python3的功能。例如，我们有：

• enum pip install enum34
• singledispatch pip install singledispatch
• pathlib pip install pathlib

想更多了解，在Python文档中有一个全面的指南可以帮助你让你的代码同时兼容Python2和Python3。

函数缓存 (Function caching)

函数缓存允许我们将一个函数对于给定参数的返回值缓存起来。当一个I/O密集的函数被频繁使用相同的参数调用的时候，函数缓存可以节约时间。在Python 3.2版本以前我们只有写一个自定义的实现。在Python 3.2以后版本，有个lru_cache的装饰器，允许我们将一个函数的返回值快速地缓存或取消缓存。

Python 3.2及以后版本

我们来实现一个斐波那契计算器，并使用lru_cache。

1
2
3
4
5
6
7
8
9
10

from functools import lru_cache

@lru_cache(maxsize=32)
def fib(n):
    if n < 2:
        return n
    return fib(n-1) + fib(n-2)

>>> print([fib(n) for n in range(10)])
# Output: [0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

那个maxsize参数是告诉lru_cache，最多缓存最近多少个返回值。

我们也可以轻松地对返回值清空缓存，通过这样：

1

fib.cache_clear()

上下文管理器(Context managers)

上下文管理器允许你在有需要的时候，精确地分配和释放资源。上下文管理器的常用于一些资源的操作，需要在资源的正确获取与释放相关的操作 ，先看一个例子,我们经常会用到 try … catch … finally 语句确保一些系统资源得以正确释放。如:

1
2
3
4
5
6
7
8

try:
    f = open('somefile')
    for line in f:
        print(line)
except Exception as e:
    print(e)
finally:
    f.close()

我们经常用到上面的代码模式，用复用代码的模式来讲，并不够好。于是 with 语句出现了，通过定义一个上下文管理器来封装这个代码块:

1
2
3

with open('somefile') as f:
    for line in f:
        print(line)

使用上下文管理器最广泛的案例就是with语句了。想象下你有两个需要结对执行的相关操作，然后还要在它们中间放置一段代码。 上下文管理器就是专门让你做这种事情的。上面这段代码打开了一个文件，往里面写入了一些数据，然后关闭该文件。如果在往文件写数据时发生异常，它也会尝试去关闭文件。这就是with语句的主要优势，它确保我们的文件会被关闭，而不用关注嵌套代码如何退出。

上下文管理器的一个常见用例，是资源的加锁和解锁，以及关闭已打开的文件（就像我已经展示给你看的）。

实际上，我们可以同时处理多个上下文管理器：

1
2

with A() as a, B() as b:
    suite

上下文管理协议

与迭代器类似，实现了迭代协议的函数/对象即为迭代器。实现了上下文协议的函数/对象即为上下文管理器。迭代器协议是实现了__iter__方法。上下文管理协议则是一个类实现__enter__ (self)和__exit__(self, exc_type, exc_valye, traceback)方法就可以了。实行如下结构：

1
2
3
4
5
6
7
8
9
10
11
12

class Contextor:
    # __enter__返回一个对象，通常是当前类的实例，也可以是其他对象。
    def __enter__(self):
        pass
    
    def __exit__(self, exc_type, exc_val, exc_tb):
        pass

contextor = Contextor()

with contextor [as var]:
    with_body

Contextor 实现了__enter__和__exit__这两个上下文管理器协议，当Contextor调用/实例化的时候，则创建了上下文管理器contextor

通过定义__enter__和__exit__方法的类（包括自己定义的类，只要加上特定的两个方法即可），我们可以在with语句里使用它。我们来看看在底层都发生了什么。

执行步骤：

1. 执行 contextor (实例化具有上下文协议的对象，这里也称为上下文表达式)以获取上下文管理器，上下文表达式就是 with 和 as 之间的代码。
2. 加载上下文管理器对象的 exit()方法，备用。
3. 调用上下文管理器的 enter() 方法
4. 如果有 as var 从句，则将 enter() 方法的返回值赋给 var
5. 执行子代码块 with_body
6. with语句调用上下文管理器之前暂存的 exit() 方法，如果 with_body 的退出是由异常引发的，那么该异常的 type、value 和 traceback 会作为参数传给 exit()，否则传三个 None。然后，exit()需要明确地返回 True 或 False。当返回 True 时，异常不会被向上抛出，当返回 False 时曾会向上抛出。
7. 如果 with_body 的退出由异常引发，它让exit()方法来处理异常，并且 exit() 的返回值等于 False，那么这个异常将被with语句重新引发抛出一次；如果 exit() 的返回值等于 True，那么这个异常就被无视掉，继续执行后面的代码。

上下文管理器工具

通过实现上下文协议定义创建上下文管理器很方便，Python为了更优雅，还专门提供了一个模块用于实现更函数式的上下文管理器用法。Python的contextlib模块专门用于这个目的。

AbstractContextManager ： 此类在 Python3.6中新增，提供了默认的enter()和exit()实现。enter()返回自身，exit()返回 None。

contextmanager： 我们要实现上下文管理器，总是要写一个类。此函数则容许我们通过一个装饰一个生成器函数得到一个上下文管理器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

import contextlib

@contextlib.contextmanager
def database():
    db = Database()
    try:
        if not db.connected:
            db.connect()
        yield db	# 生成器
    except Exception as e:
        db.close()

def handle_query():
    with database() as db:
        print 'handle ---', db.query()

使用contextlib 定义一个上下文管理器函数，通过with语句，database调用生成一个上下文管理器，然后调用函数隐式的__enter__方法，并将结果通yield返回。最后退出上下文环境的时候，在exception代码块中执行了__exit__方法。当然我们可以手动模拟上述代码的执行的细节。注意：yield 只能返回一次，返回的对象 被绑定到 as 后的变量，不需要返回时可以直接 yield，不带返回值。退出时则从 yield 之后执行。由于contextmanager继承自ContextDecorator，所以被contextmanager装饰过的生成器也可以用作装饰器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

In [1]: context = database()    # 创建上下文管理器

In [2]: context
<contextlib.GeneratorContextManager object at 0x107188f10>

In [3]: db = context.__enter__() # 进入with语句

In [4]: db                          # as语句，返回 Database实例
Out[4]: <__main__.Database at 0x107188a10>

In [5]: db.query()       
Out[5]: 'query data'

In [6]: db.connected
Out[6]: True

In [7]: db.__exit__(None, None, None)    # 退出with语句

In [8]: db
Out[8]: <__main__.Database at 0x107188a10>

In [9]: db.connected
Out[9]: False

ContextDecorator： 我们可以实现一个上下文管理器，同时可以用作装饰器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

class AContext(ContextDecorator):

    def __enter__(self):
        print('Starting')
        return self

    def __exit__(self, exc_type, exc_value, traceback):
        print('Finishing')
        return False

# 在 with 中使用
with AContext():
    print('祖国伟大')

# 用作装饰器
@AContext()
def print_sth(sth):
    print(sth)

print_sth('祖国伟大')

#在这两种写法中，有没有发现，第二种写法更好，因为我们减少了一次代码缩进，可读性更强

还有一种好处：当我们已经实现了某个上下文管理器时，只要增加一个继承类，该上下文管理器立刻编程装饰器。

1
2
3
4
5
6
7

from contextlib import ContextDecorator
class mycontext(ContextBaseClass, ContextDecorator):
    def __enter__(self):
        return self

    def __exit__(self, *exc):
        return False

协程

子程序，或者称为函数，在所有语言中都是层级调用，比如A调用B，B在执行过程中又调用了C，C执行完毕返回，B执行完毕返回，最后是A执行完毕。所以子程序调用是通过栈实现的，一个线程就是执行一个子程序。

子程序调用总是一个入口，一次返回，调用顺序是明确的。而协程的调用和子程序不同。协程看上去也是子程序，但执行过程中，在子程序内部可中断，然后转而执行别的子程序，在适当的时候再返回来接着执行。

注意，在一个子程序中中断，去执行其他子程序，不是函数调用，有点类似CPU的中断。比如子程序A、B：

1
2
3
4
5
6
7
8
9

def A():
    print '1'
    print '2'
    print '3'

def B():
    print 'x'
    print 'y'
    print 'z'

假设由协程执行，在执行A的过程中，可以随时中断，去执行B，B也可能在执行过程中中断再去执行A，结果可能是：

1
2
3
4
5
6

1
2
x
y
3
z

多线程比，协程有何优势？

最大的优势就是协程极高的执行效率。因为子程序切换不是线程切换，而是由程序自身控制，因此，没有线程切换的开销，和多线程比，线程数量越多，协程的性能优势就越明显。

第二大优势就是不需要多线程的锁机制，因为只有一个线程，也不存在同时写变量冲突，在协程中控制共享资源不加锁，只需要判断状态就好了，所以执行效率比多线程高很多。

因为协程是一个线程执行，那怎么利用多核CPU呢？最简单的方法是多进程+协程，既充分利用多核，又充分发挥协程的高效率，可获得极高的性能。Python对协程的支持是通过generator实现的。

一个例子：生产者－消费者的协程

现在我们要让生产者发送1,2,3,4,5给消费者，消费者接受数字，返回状态给生产者，而我们的消费者只需要3,4,5就行了，当数字等于3时，会返回一个错误的状态。最终我们需要由主程序来监控生产者－消费者的过程状态，调度结束程序。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

#-*- coding:utf-8
def consumer():
    status = True
    while True:
        n = yield status
        print("我拿到了{}!".format(n))
        if n == 3:
            status = False

def producer(consumer):
    n = 5
    while n > 0:
        # yield给主程序返回消费者的状态
        # consumer.send(n)把n传值给c生成器，同时返回c生成器yield的结果（相当于fetch取一个放一个东西）
        yield consumer.send(n)	
        n -= 1
 
if __name__ == '__main__':
    c = consumer()	# c产生一个生成器(带yield语句)
    c.send(None)	# consumer()程序推进到yield，但yield还未被执行.send()是传值给生成器的语句
    p = producer(c)	# p也产生一个生成器，但传入c生成器，与p进行通信
    for status in p:# 循环获取p生成器yield回来的状态
        if status == False:
            print("我只要3,4,5就行啦")
            break
    print("程序结束")

上面这个例子是典型的生产者－消费者问题，我们用协程的方式来实现它。

第一句c = consumer()，因为consumer函数中存在yield语句，python会把它当成一个generator，因此在运行这条语句后，python并不会像执行函数一样，而是返回了一个generator object。

第二条语句c.send(None)，这条语句的作用是将consumer（即变量c，它是一个generator）中的语句推进到第一个yield语句出现的位置，那么在例子中，consumer中的status = True和while True:都已经被执行了，程序停留在n = yield status的位置（注意：此时这条语句还没有被执行），上面说的send(None)语句十分重要，如果漏写这一句，那么程序直接报错

第三句p = producer(c)，这里则像上面一样定义了producer的生成器，注意的是这里我们传入了消费者的生成器，来让producer跟consumer通信。

第四句for status in p:，这条语句会循环地运行producer和获取它yield回来的状态。

现在程序流进入了producer里面，我们直接看yield consumer.send(n)，生产者调用了消费者的send()方法，把n发送给consumer（即c），在consumer中的n = yield status，n拿到的是消费者发送的数字，同时，consumer用yield的方式把状态（status）返回给消费者，注意：这时producer（即消费者）的consumer.send()调用返回的就是consumer中yield的status！消费者马上将status返回给调度它的主程序，主程序获取状态，判断是否错误，若错误，则终止循环，结束程序。上面看起来有点绕，其实这里面generator.send(n)的作用是：把n发送generator(生成器)中yield的赋值语句中，同时返回generator中yield的变量（结果）。

于是程序便一直运作，直至consumer中获取的n的值变为3！此时consumer把status变为False，最后返回到主程序，主程序中断循环，程序结束。

Coroutine与Generator

有些人会把生成器（generator）和协程（coroutine）的概念混淆，我以前也会这样，不过其实发现，两者的区别还是很大的。

直接上最重要的区别：

• generator总是生成值，一般是迭代的序列
• coroutine关注的是消耗值，是数据(data)的消费者
• coroutine不会与迭代操作关联，而generator会
• coroutine强调协同控制程序流，generator强调保存状态和产生数据

相似的是，它们都是不用return来实现重复调用的函数/对象，都用到了yield(中断/恢复)的方式来实现

asyncio

asyncio是Python 3.4版本引入的标准库，直接内置了对异步IO的支持。

asyncio的编程模型就是一个消息循环。我们从asyncio模块中直接获取一个EventLoop的引用，然后把需要执行的协程扔到EventLoop中执行，就实现了异步IO。用asyncio实现Hello world代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

import asyncio

#@asyncio.coroutine把一个generator标记为coroutine类型，然后，我们就把这个coroutine扔到EventLoop中执行。
@asyncio.coroutine
def hello():
    print("Hello world!")
    # 异步调用asyncio.sleep(1):
    r = yield from asyncio.sleep(1)
    print("Hello again!")

# 获取EventLoop:
loop = asyncio.get_event_loop()
# 执行coroutine
loop.run_until_complete(hello())	# 循环执行EventLoop里要完成的事件
loop.close()

hello()会首先打印出Hello world!，然后，yield from语法可以让我们方便地调用另一个generator。由于asyncio.sleep()也是一个coroutine，所以线程不会等待asyncio.sleep()，而是直接中断并执行下一个消息循环。当asyncio.sleep()返回时，线程就可以从yield from拿到返回值（此处是None），然后接着执行下一行语句。把asyncio.sleep(1)看成是一个耗时1秒的IO操作，在此期间，主线程并未等待，而是去执行EventLoop中其他可以执行的coroutine了，因此可以实现并发执行。

我们用Task封装两个coroutine试试：

1
2
3
4
5
6
7
8
9
10
11
12
13

import threading
import asyncio

@asyncio.coroutine
def hello():
    print('Hello world! (%s)' % threading.currentThread())
    yield from asyncio.sleep(1)
    print('Hello again! (%s)' % threading.currentThread())

loop = asyncio.get_event_loop()
tasks = [hello(), hello()]
loop.run_until_complete(asyncio.wait(tasks))
loop.close()

观察执行过程：

1
2
3
4
5

Hello world! (<_MainThread(MainThread, started 140735195337472)>)
Hello world! (<_MainThread(MainThread, started 140735195337472)>)	# asyncio.sleep(1)会挂起并去执行下一个hello()协程
(暂停约1秒)
Hello again! (<_MainThread(MainThread, started 140735195337472)>)
Hello again! (<_MainThread(MainThread, started 140735195337472)>)

总结：

• asyncio提供了完善的异步IO支持；
• 异步操作需要在coroutine中通过yield from`完成；
• 多个coroutine可以封装成一组Task然后并发执行

async/await

用asyncio提供的@asyncio.coroutine可以把一个generator标记为coroutine类型，然后在coroutine内部用yield from调用另一个coroutine实现异步操作。

为了简化并更好地标识异步IO，从Python 3.5开始引入了新的语法async和await，可以让coroutine的代码更简洁易读。

请注意，async和await是针对coroutine的新语法，要使用新的语法，只需要做两步简单的替换：

1. 把@asyncio.coroutine替换为async；
2. 把yield from替换为await。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

import asyncio

async def compute(x, y):
    print("Compute %s + %s ..." % (x, y))
    await asyncio.sleep(1.0)
    return x + y

async def print_sum(x, y):
    result = await compute(x, y)
    print("%s + %s = %s" % (x, y, result))

loop = asyncio.get_event_loop()
loop = asyncio.get_event_loop()
loop.run_until_complete(print_sum(1, 2))
# tasks = [print_sum(1, 2), print_sum(3, 4)]
# loop.run_until_complete(asyncio.wait(tasks))
loop.close()

当事件循环开始运行时，它会在Task中寻找coroutine来执行调度，因为事件循环注册了print_sum()，因此print_sum()被调用，执行result = await compute(x, y)这条语句（等同于result = yield from compute(x, y)），因为compute()自身就是一个coroutine，因此print_sum()这个协程就会暂时被挂起，compute()被加入到事件循环中，程序流执行compute()中的print语句，打印”Compute %s + %s …”，然后执行了await asyncio.sleep(1.0)，因为asyncio.sleep()也是一个coroutine，接着compute()就会被挂起，等待计时器读秒，在这1秒的过程中，事件循环会在队列中查询可以被调度的coroutine，而因为此前print_sum()与compute()都被挂起了，因此事件循环会停下来等待协程的调度(如果有其他协程task就会在等待时间内去执行并返回)，当计时器读秒结束后，程序流便会返回到compute()中执行return语句，结果会返回到print_sum()中的result中，最后打印result，事件队列中没有可以调度的任务了，此时loop.close()把事件队列关闭，程序结束。

使用C扩展

CPython还为开发者实现了一个有趣的特性，使用Python可以轻松调用C代码

开发者有三种方法可以在自己的Python代码中来调用C编写的函数-ctypes，SWIG，Python/C API。每种方式也都有各自的利弊。

首先，我们要明确为什么要在Python中调用C？

常见原因如下：

• 你要提升代码的运行速度，而且你知道C要比Python快50倍以上
• C语言中有很多传统类库，而且有些正是你想要的，但你又不想用Python去重写它们
• 想对从内存到文件接口这样的底层资源进行访问
• 不需要理由，就是想这样做

CTypes

Python中的ctypes模块可能是Python调用C方法中最简单的一种。ctypes模块提供了和C语言兼容的数据类型和函数来加载dll文件，因此在调用时不需对源文件做任何的修改。也正是如此奠定了这种方法的简单性。

示例如下

实现两数求和的C代码，保存为add.c