第四章 训练模型

这篇文章是本人学习 《Hands-On-Machine-Learning-with-Scikit-Learn-and-TensorFlow》的读书笔记第三篇。整理出来是希望在巩固自己的学习效果的同时，希望能够帮助到同样想学习的人。本人也是小白，可能很多地方理解和翻译不是很到位，希望大家多多谅解和提意见。

Setup

# Common imports
import numpy as np
import os

# to make this notebook's output stable across runs
np.random.seed(42)

# To plot pretty figures
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)

# Where to save the figures
PROJECT_ROOT_DIR = "."
CHAPTER_ID = "training_linear_models"

def save_fig(fig_id, tight_layout=True):
    path = os.path.join(PROJECT_ROOT_DIR, "images", CHAPTER_ID, fig_id + ".png")
    print("Saving figure", fig_id)
    if tight_layout:
        plt.tight_layout()
    plt.savefig(path, format='png', dpi=300)

# Ignore useless warnings (see SciPy issue #5998)
import warnings
warnings.filterwarnings(action="ignore", message="^internal gelsd")

线性回归

生成数据：

import numpy as np

X = 2 * np.random.rand(100,1)
y = 4 + 3*X +np.random.randn(100,1)

可视化数据，观察下基本特征

plt.plot(X,y,'b.')
plt.xlabel('$x_1$',fontsize=18)
plt.ylabel('$y$', rotation=0, fontsize=18)
plt.axis([0,2,0,15])
save_fig('generated_data_plot')
plt.show()

使用正规方程来求解线性回归系数：$\hat{\theta }={\left({\mathbf{X}}^T\cdot \mathbf{X}\right)}^{-1}\cdot {\mathbf{X}}^T\cdot \mathbf{y}$

X_b = np.c_[np.ones((100,1)), X] # add x0=1 to each instance
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

使用求解的系数做个预测：

X_new = np.array([[0],[2]])
X_new_b = np.c_[np.ones((2,1)), X_new]  # add x0 = 1 to each instance
y_predict = X_new_b.dot(theta_best)
y_predict

使用 Scikit-Learn 中的 LinearRegression 来求解线性回归的系数。

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X,y)
lin_reg.intercept_, lin_reg.coef_

使用 lin_reg.intercept_ 显示出截距项，lin_reg.coef_ 显示出系数。

LinearRegression 类求解实际是基于 scipy.linalg.lstsq()函数，可以直接使用这个函数来求解。

theta_best_svd, residuals, rank, s = np.linalg.lstsq(X_b, y, rcond=1e-6)
theta_best_svd

这个函数通过计算 ${\mathbf{X}}^{+}\mathbf{y}$的值来求解，其中 ${\mathbf{X}}^{+}$ 是$\mathbf{X}$的伪逆。可以使用 np.linalg.pinv()函数来直接计算伪逆。

np.linalg.pinv(X_b).dot(y)

时间复杂度

${\mathbf{X}}^T\cdot \mathbf{X}$ 是一个 $n\ast n$ 的矩阵，计算它的逆需要的时间为 $O(n^3)$。使用正规方程求解，时间复杂度随着训练集的大小成正比增加，当数据量较大且能够全部被读入内存时，这还是有优势的。另外当模型的特征过多时，算法的复杂度增加太快。所以当特征数较多，数据不能一次读入内存时，应该考虑其他的求解方法。

梯度下降

批量梯度下降

损失函数的偏导数：$\frac{\partial }{\partial {\theta }_j}\mathrm{MSE}(\theta )=\frac{2}{m}{\sum }_{i=1}^m\left({\theta }^T\cdot {\mathbf{x}}^{(i)}-y^{(i)}\right)x_j^{(i)}$
损失函数的梯度向量：${\nabla }_{\theta }MSE(\theta )=\left(\begin{array}{c}\frac{\partial }{\partial {\theta }_0}MSE(\theta )\\ \frac{\partial }{\partial {\theta }_1}MSE(\theta )\\ ⋮\\ \frac{\partial }{\partial {\theta }_n}MSE(\theta )\end{array}\right)=\frac{2}{m}{\mathbf{X}}^T\cdot (\mathbf{X}\cdot \theta -y)$
在这个方程中每一步计算都包含了整个数据集，因此在大数据集上速度会变得很慢。但是梯度下降的运算规模和特征数量成正比，因此在数千特征的数据集上，梯度下降比正规方程快。

梯度下降的步长：${\theta }^{(\text{next step})}=\theta -\eta {\nabla }_{\theta }MSE(\theta )$

eta = 0.1
n_iterations = 1000
m = 100
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

当使用不同的学习率是会发生什么，下图展示了三个使用不同的学习率求解的前10部运算（虚线代表初始位置）。

theta_path_bgd = []

def plot_gradient_descent(theta, eta, theta_path=None):
    m = len(X_b) # 样本数
    plt.plot(X,y,'b.')
    n_iterations = 1000
    for iteration in range(n_iterations):
        if iteration < 10:
            y_predict = X_new_b.dot(theta)
            style = 'b-' if iteration > 0 else 'r--'
            plt.plot(X_new, y_predict, style)
        gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
        theta = theta - eta * gradients
        if theta_path is not None:
            theta_path.append(theta)
    plt.xlabel('$x_1$', fontsize=18)
    plt.axis([0,2,0,15])
    plt.title(r'$\eta = {}$'.format(eta), fontsize=16)


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

plt.figure(figsize=(10,4))
plt.subplot(131); plot_gradient_descent(theta, eta=0.02)
plt.ylabel('$y$', rotation=0, fontsize=18)
plt.subplot(132); plot_gradient_descent(theta, eta=0.1, theta_path = theta_path_bgd)
plt.subplot(133); plot_gradient_descent(theta, eta=0.5)

save_fig("gradient_descent_plot")
plt.show()

随机梯度下降

随机梯度下降在每一步计算时，随机选取训练集中的一个数据计算梯度。这使得它在大规模数据上训练更快。当损失函数很不规则时， 随机梯度下降算法能够跳过局部最小值。 因此， 随机梯度下降在寻找全局最小值上比批量梯度下降表现要好。

theta_path_sgd = []
m = len(X_b)
np.random.seed(42)

n_epochs = 50
t0, t1 = 5, 50 # learning schedule hyperparameters

def learning_schedule(t):
    return t0 / (t + t1)

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

for epoch in range(n_epochs):
    for i in range(m):
        if epoch == 0 and i < 10:
            y_predict = X_new_b.dot(theta)
            style = 'b-' if i > 0 else 'r--'
            plt.plot(X_new, y_predict, style)
        random_index = np.random.randint(m)
        xi = X_b[random_index:random_index+1]  # 每次只选取一个点
        yi = y[random_index:random_index+1]
        gradients = xi * xi.T.dot(xi.dot(theta) - yi)  # 针对该点计算梯度
        eta = learning_schedule(epoch *m + i)
        theta = theta -eta * gradients
        theta_path_sgd.append(theta)
        
plt.plot(X, y, "b.")                               
plt.xlabel("$x_1$", fontsize=18)                    
plt.ylabel("$y$", rotation=0, fontsize=18)          
plt.axis([0, 2, 0, 15])                            
save_fig("sgd_plot")                  
plt.show()

如果想保证每一代的迭代过程，算法能够遍历所有实例，可以先将训练数据打乱，然后逐个遍历数据。下一代训练时再重复这个过程，但是这样会使得训练的速度变慢。

使用 Scikit-Learn 的 SGDRegressor 类完成线性回归的梯度下降，这个类默认优化的是均方差损失函数。

from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=50, tol=-np.infty, penalty=None, 
                       eta0 = 0.1, random_state=42)  # 迭代50代，学习率为0.1
sgd_reg.fit(X,y.ravel())

使用 sgd_reg.intercept_, sgd_reg.coef_ 分别获取截距和系数。

小批量梯度下降

小批量梯度下降中， 它则使用一个随机的小型数据集。 它比随机梯度的主要优点在于你可以通过矩阵运算的硬件优化得到一个较好的训练表现， 尤其当你使用 GPU 进行运算的时候。

小批量梯度下降在参数空间上的表现比随机梯度下降要好的多， 尤其在有大量的小型实例集时。 因此，小批量梯度下降会比随机梯度更靠近最小值。 但是， 另一方面， 它有可能陷在局部最小值中。下图显示了训练期间三种梯度下降算法在参数空间中所采用的路径。 他们都接近最小值，但批量梯度的路径最后停在了最小值， 而随机梯度和小批量梯度最后都在最小值附近摆动。

theta_path_mgd = []

n_iterations = 50
minibatch_size = 20

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

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)

多项式回归

通过加入每个特征的阶数来构造新的特征，在新的特征上进行线性回归。

数据集：

m = 100
X = 6 * np.random.rand(m,1) - 3
y = 0.5 * X**2 +X + 2 + np.random.randn(m,1)

plt.plot(X, y, 'b.')
plt.xlabel('$x_1$',fontsize=18)
plt.ylabel('$y$',rotation=0, fontsize=18)
plt.axis([-3,3,0,10])
save_fig("quadratic_data_plot")
plt.show()

我们使用 Scikit-Learning 的 PolynomialFeatures 类进行训练数据集的转换， 加入训练集中每个特征的平方（ 2 次多项式）作为新特征。

from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)

lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
lin_reg.intercept_, lin_reg.coef_

学习曲线

第二章中我们知道了可以利用交叉验证的方法来评估模型的泛化能力。另一种方法就是观察学习曲线。在学习曲线上我们画出模型在训练集和验证集上的表现，作为训练集规模的函数。为了产生这些图形，只需要在不同大小的训练集上多训练几次即可。

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, random_state=10)
    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[:m], y_train_predict))
        val_errors.append(mean_squared_error(y_val, y_val_predict))

    plt.plot(np.sqrt(train_errors), "r-+", linewidth=2, label="train")
    plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="val")
    plt.legend(loc="upper right", fontsize=14)   
    plt.xlabel("Training set size", fontsize=14) 
    plt.ylabel("RMSE", fontsize=14) 

lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)
plt.axis([0, 80, 0, 3])                         # not shown in the book
save_fig("underfitting_learning_curves_plot")   # not shown
plt.show()  

这是个典型的欠拟合模型的学习曲线，两条线都趋于平缓，很接近但是都在高位。

接下来我们看下在相同数据上的十阶多项式模型的学习曲线。

from sklearn.pipeline import Pipeline

polynomial_regression = Pipeline([
        ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
        ("lin_reg", LinearRegression()),
    ])

plot_learning_curves(polynomial_regression, X, y)
plt.axis([0, 80, 0, 3])          
save_fig("learning_curves_plot")  
plt.show()       

这幅图中两条线都处于相对低位，且两者存在较大差距。这意味着模型在训练集上的表现比验证集上要好很多，这是过拟合的特征。

偏差与方差均衡
偏差：泛化误差是由于错误的假设，高偏差的模型更可能是欠拟合。
方差：这部分误差主要是模型对训练数据中的小变动异常敏感造成的，是过拟合的结果。

正则化的线性模型

岭回归

岭回归是线性回归的正则化版本，在其损失函数后加入了 正则化项。这迫使模型在拟合数据的同时，使得模型的参数权重尽可能小。注意这个正则化项只有在训练时才会被加入到损失函数中。当模型训练好以后，需要用没有正则化的测量方法去评估模型。

岭回归损失函数：$J(\theta )=MSE(\theta )+\alpha \frac{1}{2}{\sum }_{i=1}^n{\theta }_i^2$
岭回归的解： $\hat{\theta }={\left({\mathbf{X}}^T\cdot \mathbf{X}+\alpha \mathbf{A}\right)}^{-1}\cdot {\mathbf{X}}^T\cdot \mathbf{y}$
在使用岭回归以前进行数据缩放是很有必要的，因为它对输入特征的尺度很敏感。

使用 Scikit-Learn 求解岭回归：

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

使用随机梯度下降法求解：

sgd_reg = SGDRegressor(max_iter=50, tol=-np.infty, penalty="l2", random_state=42)
sgd_reg.fit(X, y.ravel())
sgd_reg.predict([[1.5]])

Lasso 回归

像岭回归一样，Lasso回归也在损失函数上添加了一个正则化项，但是它使用的是 $l_1$ 范数而不是$l_2$ 范数。
Lasso 回归的损失函数：$J(\theta )=MSE(\theta )+\alpha {\sum }_{i=1}^n\left|{\theta }_i\right|$
Lasso回归子梯度向量：
$g(\theta ,J)={\nabla }_{\theta }MSE(\theta )+\alpha \left(\begin{array}{c}\mathrm{sign}\left({\theta }_1\right)\\ \mathrm{sign}\left({\theta }_2\right)\\ ⋮\\ \mathrm{sign}\left({\theta }_n\right)\end{array}\right)$
where sign $\left({\theta }_i\right)=\{\begin{array}{ll}-1,& {\theta }_i<0\\ 0,& {\theta }_i=0\\ +1,& {\theta }_i>0\end{array}$

使用 Scikit-Learn 求解Lasso回归：

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

Lasso回归的一个重要特点是它倾向于完全忽略不重要特征的权重。换句话说，Lasso回归会自动进行特征选择并产生一个稀疏的模型。

弹性网络(ElasticNet)

弹性网络介于 Ridge 回归和 Lasso 回归之间，它的正则项是 Ridge 回归和 Lasso 回归正则项的简单混合。

弹性网络损失函数： $J(\theta )=MSE(\theta )+r\alpha {\sum }_{i=1}^n\left|{\theta }_i\right|+\frac{1-r}{2}\alpha {\sum }_{i=1}^n{\theta }_i^2$
使用 Scikit-Learn 求解ElasticNet：

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

如何在线性回归、岭回归、Lasso回归、Elastic Net中做选择？通常来说应该避免直接使用简单的线性回归。岭回归通常是个不错的选择，如果你担心只有部分特征是有用的，可以考虑 Lasso 或者 Elastic Net，因为它们会倾向于将不重要的参数权重降到0。通常来说，Elastic Net优于Lasso因为当特征的数量大于训练集数量或者几个特征存在相关性时Lasso会很不稳定。

早期停止法

对于迭代学习算法，有一种非常特殊的正则化方法，当验证误差达到最小时便停止训练，这种方法被称为早期停止法。

np.random.seed(42)
m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 2 + X + 0.5 * X**2 + np.random.randn(m, 1)

X_train, X_val, y_train, y_val = train_test_split(X[:50], y[:50].ravel(), test_size=0.5, random_state=10)

poly_scaler = Pipeline([
        ("poly_features", PolynomialFeatures(degree=90, include_bias=False)),
        ("std_scaler", StandardScaler()),
    ])

X_train_poly_scaled = poly_scaler.fit_transform(X_train)
X_val_poly_scaled = poly_scaler.transform(X_val)

sgd_reg = SGDRegressor(max_iter=1,
                       tol=-np.infty,
                       penalty=None,
                       eta0=0.0005,
                       warm_start=True,
                       learning_rate="constant",
                       random_state=42)

n_epochs = 500
train_errors, val_errors = [], []
for epoch in range(n_epochs):
    sgd_reg.fit(X_train_poly_scaled, y_train)
    y_train_predict = sgd_reg.predict(X_train_poly_scaled)
    y_val_predict = sgd_reg.predict(X_val_poly_scaled)
    train_errors.append(mean_squared_error(y_train, y_train_predict))
    val_errors.append(mean_squared_error(y_val, y_val_predict))

best_epoch = np.argmin(val_errors)
best_val_rmse = np.sqrt(val_errors[best_epoch])

plt.annotate('Best model',
             xy=(best_epoch, best_val_rmse),
             xytext=(best_epoch, best_val_rmse + 1),
             ha="center",
             arrowprops=dict(facecolor='black', shrink=0.05),
             fontsize=16,
            )

best_val_rmse -= 0.03  # just to make the graph look better
plt.plot([0, n_epochs], [best_val_rmse, best_val_rmse], "k:", linewidth=2)
plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="Validation set")
plt.plot(np.sqrt(train_errors), "r--", linewidth=2, label="Training set")
plt.legend(loc="upper right", fontsize=14)
plt.xlabel("Epoch", fontsize=14)
plt.ylabel("RMSE", fontsize=14)
save_fig("early_stopping_plot")
plt.show()

下面是一个早期停止法的基本运用。

from sklearn.base import clone
sgd_reg = SGDRegressor(max_iter=1, tol=-np.infty, warm_start=True, penalty=None,
                       learning_rate="constant", eta0=0.0005, random_state=42)

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)  # continues where it left off
    y_val_predict = sgd_reg.predict(X_val_poly_scaled)
    val_error = mean_squared_error(y_val, y_val_predict)
    if val_error < minimum_val_error:
        minimum_val_error = val_error
        best_epoch = epoch
        best_model = clone(sgd_reg)

best_epoch, best_model

逻辑回归

逻辑回归被用来估计一个实例属于某个特定类的概率，如果估计的概率大于50%，那么模型预测这个实例属于当前类，否则不属于该类。，Logistic 回归模型计算输入特征的加权和（加上偏差项），但它不像线性回归模型那样直接输出结果，而是把结果输入 logistic() 函数进行二次加工后进行输出。

逻辑回归模型的概率估计（向量形式）：$\hat{p}=h_{\theta }(\mathbf{x})=\sigma \left({\theta }^T\cdot \mathbf{x}\right)$
逻辑函数： $\sigma (t)=\frac{1}{1+\exp (-t)}$
逻辑函数的图像：

逻辑回归预测模型：$\hat{y}=\{\begin{array}{ll}0,& \hat{p}<0.5\\ 1,& \hat{p}\ge 0.5\end{array}$

逻辑回归的损失函数：
$J(\theta )=-\frac{1}{m}{\sum }_{i=1}^m\left[y^{(i)}\log \left({\hat{p}}^{(i)}\right)+\left(1-y^{(i)}\right)\log \left(1-{\hat{p}}^{(i)}\right)\right]$

这个函数对于求解最小化损失函数的 $\theta$ 是没有公式解的（没有等价的正规方程）。我们只能用梯度下降的方法找到全局最小值。

逻辑回归损失函数的偏导数：
$\frac{\partial }{\partial {\theta }_j}J\left({\theta }_j\right)=\frac{1}{m}{\sum }_{i=1}^m\left(\sigma \left({\theta }^T\cdot {\mathbf{x}}^{(i)}\right)-y^{(i)}\right)x_j^{(i)}$

决策边界

在Iris数据集中，我们只基于花瓣宽度的特征构造一个分类器去检测 Iris-Virginica 类型。

from sklearn import datasets
iris = datasets.load_iris()
list(iris.keys())

特征抽取

X = iris['data'][:,3:] # petal width
y = (iris['target'] == 2).astype(np.int)  # 1 if Iris-Virginica, else 0

训练一个逻辑回归模型

from sklearn.linear_model import LogisticRegression
log_reg = LogisticRegression(solver="liblinear", random_state=42)
log_reg.fit(X, y)

观察下模型估计的花瓣宽度在0到3里面的概率估计。

X_new = np.linspace(0, 3, 1000).reshape(-1, 1)
y_proba = log_reg.predict_proba(X_new)
decision_boundary = X_new[y_proba[:, 1] >= 0.5][0]

plt.figure(figsize=(8, 3))
plt.plot(X[y==0], y[y==0], "bs")
plt.plot(X[y==1], y[y==1], "g^")
plt.plot([decision_boundary, decision_boundary], [-1, 2], "k:", linewidth=2)
plt.plot(X_new, y_proba[:, 1], "g-", linewidth=2, label="Iris-Virginica")
plt.plot(X_new, y_proba[:, 0], "b--", linewidth=2, label="Not Iris-Virginica")
plt.text(decision_boundary+0.02, 0.15, "Decision  boundary", fontsize=14, color="k", ha="center")
plt.arrow(decision_boundary, 0.08, -0.3, 0, head_width=0.05, head_length=0.1, fc='b', ec='b')
plt.arrow(decision_boundary, 0.92, 0.3, 0, head_width=0.05, head_length=0.1, fc='g', ec='g')
plt.xlabel("Petal width (cm)", fontsize=14)
plt.ylabel("Probability", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 3, -0.02, 1.02])
save_fig("logistic_regression_plot")
plt.show()

使用 decision_boundary 可得出模型的决策边界，当特征值大于该边界是为正类，否则为负类。

下图使用相同的数据集，但是采用花瓣的宽度和长度两个特征。一旦训练完毕，Logistic 回归分类器就可以根据这两个特征来估计一朵花是 Virginica 的可能性。 虚线表示这时两类情况出现的概率都等于 50%：这是模型的决策边界。

from sklearn.linear_model import LogisticRegression

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

log_reg = LogisticRegression(solver="liblinear", C=10**10, random_state=42)
log_reg.fit(X, y)

x0, x1 = np.meshgrid(
        np.linspace(2.9, 7, 500).reshape(-1, 1),
        np.linspace(0.8, 2.7, 200).reshape(-1, 1),
    )
X_new = np.c_[x0.ravel(), x1.ravel()]

y_proba = log_reg.predict_proba(X_new)

plt.figure(figsize=(10, 4))
plt.plot(X[y==0, 0], X[y==0, 1], "bs")
plt.plot(X[y==1, 0], X[y==1, 1], "g^")

zz = y_proba[:, 1].reshape(x0.shape)
contour = plt.contour(x0, x1, zz, cmap=plt.cm.brg)


left_right = np.array([2.9, 7])
boundary = -(log_reg.coef_[0][0] * left_right + log_reg.intercept_[0]) / log_reg.coef_[0][1]

plt.clabel(contour, inline=1, fontsize=12)
plt.plot(left_right, boundary, "k--", linewidth=3)
plt.text(3.5, 1.5, "Not Iris-Virginica", fontsize=14, color="b", ha="center")
plt.text(6.5, 2.3, "Iris-Virginica", fontsize=14, color="g", ha="center")
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.axis([2.9, 7, 0.8, 2.7])
save_fig("logistic_regression_contour_plot")
plt.show()

Softmax回归

Logistic回归可以被推广到支持多类别的分类而不用组合和训练多个二分类器。基本思路是当给定一个实例 $x$ 时，Softmax回归模型首先计算每个类别 $K$ 的分数 $s_k(x)$，然后通过 Softmax函数计算出每类的概率。

K 类的Softmax 得分：$s_k(\mathbf{x})={\theta }^T\cdot \mathbf{x}$。

计算出样本属于每一类的得分后，可以通过 Softmax函数估计出样本属于第$k$类的概率 ${\hat{p}}_k$。

Softmax 函数：${\hat{p}}_k=\sigma (\mathbf{s}(\mathbf{x}){)}_k=\frac{\exp \left(s_k(\mathbf{x})\right)}{{\sum }_{j=1}^K\exp \left(s_j(\mathbf{x})\right)}$

• $K$ 表示有多少类
• $s(x)$是一个向量，里面包含着$x$属于每一类的得分
• $\sigma (\mathbf{s}(\mathbf{x}){)}_k$表示给定每一类的分数后，实例 $x$属于第$k$类的概率。

和 Logistic 回归分类器一样，Softmax 回归分类器将估计概率最高的那类作为预测结果。注意，Softmax回归分类器每次只能预测一个类（它是多类的，但不是多结果的），因此它只能被用来判断互斥的类别。不能用它来识别一张照片中的几个人。

Softmax 回归模型分类器预测结果：
$\hat{y}=\mathrm{argmax}\sigma (\mathbf{s}(\mathbf{x}){)}_k=\mathrm{argmax}{s}_k(\mathbf{x})=\mathrm{argmax}\left({\theta }_k^T\cdot \mathbf{x}\right)$

损失函数，也成为交叉熵，通常用来衡量待测类别与目标类别的匹配程度。$J(\Theta )=-\frac{1}{m}{\sum }_{i=1}^m{\sum }_{k=1}^K{y}_k^{(i)}\log \left({\hat{p}}_k^{(i)}\right)$。如果对于第 $i$个实例的目标类是 $k$，那么 $y_k^{(i)}=1$，反之 $y_k^{(i)}=0$。

关于损失函数的梯度向量为：
${\nabla }_{{\theta }_k}J(\Theta )=\frac{1}{m}{\sum }_{i=1}^m\left({\hat{p}}_k^{(i)}-y_k^{(i)}\right){\mathbf{x}}^{(i)}$

现在可以计算每一类的梯度向量，然后使用梯度下降或其他的优化算法找到使损失函数最小的参数矩阵 $\Theta$。当我们使用 Scikit-Learn 的 LogisticRegression 类来进行多分类时，默认采用1对所有的方式。但是当我们设置 multi_class = multinomial 时，则切换到 Softmax 回归。

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

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

画出分类决策边界

x0, x1 = np.meshgrid(
        np.linspace(0, 8, 500).reshape(-1, 1),
        np.linspace(0, 3.5, 200).reshape(-1, 1),
    )
X_new = np.c_[x0.ravel(), x1.ravel()]


y_proba = softmax_reg.predict_proba(X_new)
y_predict = softmax_reg.predict(X_new)

zz1 = y_proba[:, 1].reshape(x0.shape)
zz = y_predict.reshape(x0.shape)

plt.figure(figsize=(10, 4))
plt.plot(X[y==2, 0], X[y==2, 1], "g^", label="Iris-Virginica")
plt.plot(X[y==1, 0], X[y==1, 1], "bs", label="Iris-Versicolor")
plt.plot(X[y==0, 0], X[y==0, 1], "yo", label="Iris-Setosa")

from matplotlib.colors import ListedColormap
custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])

plt.contourf(x0, x1, zz, cmap=custom_cmap)
contour = plt.contour(x0, x1, zz1, cmap=plt.cm.brg)
plt.clabel(contour, inline=1, fontsize=12)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 7, 0, 3.5])
save_fig("softmax_regression_contour_plot")
plt.show()

程序

我把书中的程序都用 Python 3 运行了一遍，确保没有Bug并且都加了注释，方便大家理解。原书的数据集和代码在这个网站上，我自己运行的程序在我的GitHub上。