用Scikit-learn和TensorFlow进行机器学习（四）

训练模型概述

作用：快速找到恰当的机器学习模型、训练算法、好的额假设集、有助于调试参数和有效的误差分析。

以简单线性回归为例，讨论不同训练方法的模型最优解问题：

• 直接采用封闭方程进行求解，得到当前训练集的最优参数（最小化损失函数）
• 迭代优化（梯度下降，GD）：通过逐渐调整模型参数以获得最小的损失函数。批量梯度下降（Batch GD）、小批量梯度下降（Mini-batch GD）、随机梯度下降（Stochastic GD）

多项式回归模型：通过曲线判断模型是否过拟合，并介绍正则化方法减少过拟合风险。

Logistic回归和Softmax回归

一、线性回归

重要：理解思想

1、模型

线性回归预测模型：
$$\hat{y}={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+...++{\theta }_n{x}_n$$
其中，$\hat{y}$ 表示预测结果，$n$ 表示特征的个数，$x_i$ 表达第 $i$ 个特征的值，${\theta }_j$ 表达第 $j$ 个参数（包括偏置项${\theta }_0$和特征权重值${\theta }_1,{\theta }_2,...,{\theta }_n$）

线性回归预测模型（向量形式）
$$\hat{y}=h_{\theta }(X)={\theta }^{T}X$$
其中，$\theta$ 表示模型的参数向量包括偏置项${\theta }_0$和特征权重值${\theta }_1,{\theta }_2,...,{\theta }_n$；${\theta }^T$ 表示向量 $\theta$ 的转置（行向量转换为列向量）；$X$ 为每个样本中特征值的向量形式，包括 $x_1$ 到 $x_n$，且 $x_0$ 恒为1；$h_{\theta }$ 表示假设函数。

2、评估方法

评估方法：均方根误差（RMSE），常用于回归模型。最小化均方根差等价于最小化均方根误差且更加简单，两者可得到相同的参数 $\theta$。

损失函数：$h_{\theta }$ 的均方差（MSE）
$$MSE(X,h_{\theta })=\frac{1}{m}\sum _{i=1}^m({\theta }^T\cdot x^{(i)}-y^{(i)}{)}^2$$

3、正规方程（The Normal Equation）

解方程得到最优解
$$\hat{\theta }=(X^T\cdot X{)}^{-1}\cdot X^T\cdot y$$
其中，$\hat{\theta }$ 表示最小化损失 $\theta$ 的值；$y$ 是一个向量，其包含了 $y^{(1)}$ 到 $y^{(m)}$ 的值。

import numpy as np
import matplotlib.pyplot as plt

# 生成数据
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

## 参数求解：inv() 计算矩阵的逆，dot() 计算矩阵的乘法
X_b = np.c_[np.ones((100, 1)), X]
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print(theta_best)

## 利用theta_best预测
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
y_predict = X_new_b.dot(theta_best)
print(y_predict)

## 绘制图像
plt.plot(X_new, y_predict,"r--",label="Predict")
plt.plot(X, y, "b.", label="Train")
plt.axis([0, 2, 0, 15])
plt.show()

输出结果

[[4.1001622 ]
 [2.89052294]]
[[4.1001622 ]
 [9.88120808]]

等价方法

from sklearn.linear_model import LinearRegression
line_reg = LinearRegression()
line_reg.fit(X, y)
print(line_reg.intercept_, line_reg.coef_)
print(line_reg.predict(X_new))

输出结果

[4.1001622] [[2.89052294]]
[[4.1001622 ]
 [9.88120808]]

4、计算复杂度

一个矩阵求逆的运算复杂度大约在 $O(n^{2.4})$ 到 $O(n^3)$ 之间，具体值取决于计算方式。当特征个数较大时，正规方程求解将会非常慢，但是线性回归模型预测是非常快的。

5、梯度下降

（1）概述

思路：通过迭代来逐渐调整参数使得损失函数达到最小值。

具体来说：开始选定一个随机的 $\theta$（称为随机初始值），然后计算参数 $\theta$ 的局部梯度，同时沿着梯度下降得方向进行下一次迭代，直到算法收敛到一个最小值。

参数：步长，超参数学习率的值决定了步长的大小。若学习率太小，迭代次数多，耗时；若学习率太大，可能不收敛、产生跌荡、发散。

两个主要挑战：

• 随机初始值可能收敛到局部最小值，而非全局最小值。（初始值在左侧）
• 需要很长时间才能收敛，若过早结束训练，则可能不会到全局最小值。（初始值在右侧）
  

（2）线性回归的梯度下降

梯度下降可以无线接近全局最小值：

• 线性回归模型的均方差损失函数是一个凸函数，表示损失函数没有局部最小值，只有一个全局最小值。
• 且它也是一个斜率不能突变的连续函数。

注意：在使用梯度下降的时候，应保证所有的特征有着相近的尺度范围（eg：StandardScaler类）
，否则它需要很长时间才能收敛。

（3）批量梯度下降

计算每一个 ${\theta }_j$ 下损失函数的梯度（求偏导数），记为 $\frac{\partial }{\partial {\theta }_j}MSE(\theta )$

损失函数的偏导数
$$\frac{\partial }{\partial {\theta }_j}MSE(\theta )=\frac{2}{m}\sum _{i=1}^m({\theta }^T\cdot x^{(i)}-y^{(i)})x_j^{(i)}$$

梯度向量计算 ${\nabla }_{\theta }MSE(\theta )$，其包含损失函数所有的偏导数
$${\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}{X}^T\cdot (X\cdot \theta -y)$$

该算法每一次训练过程都是用所有的训练数据。

参数更新：
$${\theta }^{nextstep}=\theta -\eta {\nabla }_{\theta }MSE(\theta )$$
其中，学习率 $\eta$ 和梯度向量的乘积决定了每一次更新的步长。

eta = 0.1  			 # learning rate
n_iterations = 60000 # 迭代次数
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

print(theta)

输出结果

[[4.19096745]
 [2.91710735]]

学习率的选取：网格搜索

迭代次数选取：通常设置一个非常大的迭代次数，但是当梯度向量变得非常小的时候，结束迭代。即，梯度向量小于一个值 $\epsilon$ （称为容差）

收敛速率：损失函数是凸函数，且斜率不能突变，则学习率固定后，它的收敛速率为 $O(\frac{1}{iterations})$

缺点：在较大规模数据集上计算每一步的梯度时非常慢。

（4）随机梯度下降

关键点：每次迭代，随机在训练集中选取训练集中的一个样本，所以随机梯度下降算法可以在大规模数据记上使用。
==》由于随机性，所以呈现更多的不规律性：它达到最小值不是平缓的下降，损失函数会忽高忽低，只是在大体上呈下降趋势。当算法停止时，最后的参数还不错，但非最优值。

优点：

• 随机梯度下降在寻找全局最小值上比批量梯度下降表现好，它可以跳过局部最优值。

可以跳过局部最优值，但是它却不能达到最小值。
==》解决方法：模拟退火（逐渐降低学习率：开始时，学习率较大，然后变得越来越小）
==》决定每次迭代的学习率的函数：learning schedule

n_epochs = 50
t0, t1 = 5, 50

def learning_schedule(t):
    return t0 / (t + t1)
theta = np.random.randn(2, 1)
# 每一轮迭代：一代，epoch
for epoch in range(n_epochs): 
    for i in range(m):
        random_index = np.random.randint(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 + 1)
        theta = theta - eta * gradients
print(theta)

输出结果

[[4.07764462]
 [2.83037765]]

由于每个实例的选择是随机的，所以有的实例可能未被选到，若要遍历所有实例：

• 打乱训练集，然后依次选择实例，但收敛速度会慢；
• sklearn中的 SGDRegressor 类，该类默认优化的是均方差损失函数。

迭代50代，学习率为0.1，使用默认的 learning schedule，没有正则项（penalty=None）

from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(n_iter=50, penalty=None, eta0=0.1)
sgd_reg.fit(X, y.ravel())
## 输出参数
print(sgd_reg.intercept_, sgd_reg.coef_)

输出结果

[3.96146609] [3.15715714]

（5）小批量梯度下降

使用一个随机的小型实例集进行每次的迭代，可以通过矩阵运算的硬件优化得到一个较好的训练表现，尤其使用GPU时。

不同方法的参数空间的梯度下降路径

如果你使用了一个较好的 learning schedule ， 随机梯度和小批量梯度也可以得到最小值。

6、优化算法比较

算法	样本个数m	Out-of-core support	特征个数 n	超参数个数	scaling	sklearn库
正规方程	Fast	No	Slow	0	No	LinearRegression
Batch GD	slow	No	Fast	2	Yes	N/A
Stochastic GD	Fast	Yes	Fast	>=2	Yes	SGDRegressor
Mini-batch GD	Fast	Yes	Fast	>=2	Yes	N/A

二、多项式回归

使用线性模型拟合非线性数据。

多项式回归：对每个特征进行加权后作为新的特征，然后训练一个线性模型在这个扩展的特征集。

1、数据生成

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.show()

输出结果

2、拟合

sklearn 中 PolynomialFeatures 类进行训练集的转换，让训练集中的每个特侦的平方（2次多项式）作为新特征。

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

## degree定义当前阶数下特征的所有组合
poly_features = PolynomialFeatures(degree=2, include_bias=False)
# X_ploy 包含原始特征X并加上这个特征的平方 X^2
X_ploy = poly_features.fit_transform(X)
print(X[0])
print(X_ploy[0])

## 训练
line_reg = LinearRegression()
line_reg.fit(X_ploy, y)
print(line_reg.intercept_, line_reg.coef_)

输出结果

[-0.54569402]
[-0.54569402  0.29778196]
[2.03735568] [[0.94998295 0.51120366]]

==》
模型预测函数为 $\hat{y}=0.51120366x_1^2+0.94998295x_1+2.03735568$，
原始函数为 $y=0.5x_1^2+1.0x_1+2.0$ 再加上一些高斯噪声。

注意：PolynomialFeatures(degree=d) 将一个包含 $n$ 个特征的数组转换为一个包含 $\frac{(n+d)!}{d!n!}$ 特征的数组。

3、学习曲线

对于数据生成函数未知的数据，如何决定模型的复杂度？如何检验模型是过拟合还是欠拟合？

• 交叉验证：
• 观察学习曲线

（1）欠拟合

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.3)
    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")

line_reg = LinearRegression()
plot_learning_curves(line_reg, X, y)
plt.legend()
plt.show()

欠拟合：两条曲线都达到高原地带并趋于稳定，且卒后两条曲线非常接近，同时误差值非常大。
==》更复杂的模型 或 更好的特征

（2）过拟合

相同数据上10阶多项式模型拟合的学习曲线

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)
plt.axis([0, 70, 0, 3])
plt.legend()
plt.show()

区别：

• 在训练集上，误差比线性模型低很多；
• 两条曲线之间有间隔，表示模型在训练集表现比验证集上好。

==》改善方法：增加训练数据，直到训练误差和验证误差相等、正则化

（3）偏差和方差

一个模型的泛化误差由三个不同误差的和决定：偏差、方差和不可约误差

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

三、正则化

限制模型的复杂度。

对于线性模型，正则化的典型实现就是约束模型中参数的权重。三种不同约束权重的方法：Ridge回归、Lasso回归和Elastic Net

1、岭（Ridge）回归

（1）概述

岭回归（也称Tikhonov正则化）：在损失函数上直接加上一个正则项 $\alpha {\sum }_{i=1}^n{\theta }_i^2$。其中超参数 $\alpha$ 决定正则化模型的强度。
==》拟合数据+模型参数小
==》正则项仅在训练过程中有，在训练完之后，使用没有正则化的测量方法来评价模型的表现。

一般情况下，训练过程中使用的损失函数和测试过程中使用的评价函数是不一样的。此外，训练时的损失函数应该在优化过程中易于求导，而在测试过程中，评价函数更应该接近最后的客观表现。eg：在分类训练中使用对数损失作为损失函数，而使用精确率/召回率来作为它的评价函数。

（2）岭回归的损失函数

$$J(\theta )=MSE(\theta )+\alpha \frac{1}{2}\sum _{i=1}^n{\theta }_i^2$$
注意，偏差 ${\theta }_0$ 没有被正则化的。
定义 $W$ 表示特征的权重向量（${\theta }_1$到${\theta }_n$），则正则化项可简写为 $\frac{1}{2}(||W|{|}_2{)}^2$，其中 $||\cdot |{|}_2$ 表示权重的 $l_2$ 范数。

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

下图为相同线性数据上使用不同的 $\alpha$ 值的岭回归模型最后的表现。其中，左图为简单的岭回归模型，右图为10阶的 PolynomialFeatures 进行扩展。

分析：当 $\alpha$ 增大时，导致预测曲线变得扁平（即少了极端值，多了一般值），这样减少了模型的方差，却增加了模型的偏差。

（3）岭回归的封闭方程

岭回归的封闭方程的解：
$$\hat{\theta }=(X^T\cdot X+\alpha A{)}^{-1}\cdot X^T\cdot y$$
矩阵 $A$ 是一个除了左上角有一个 0 的 $n\times n$ 的单位矩阵，这个 0 代表偏差项。注意：偏差 ${\theta }_0$ 不被正则化的。

from sklearn.linear_model import Ridge

## 封闭方程求解——cholesky法：进行矩阵分解
ridge_reg = Ridge(alpha=1, solver="cholesky")
ridge_reg.fit(X, y)
predict1 = ridge_reg.predict([[1.5]])
print("The predict is ", predict1)

## 梯度下降法求解
# penalty参数：正则项的惩罚类型
sgd_reg = SGDRegressor(penalty="l2")
sgd_reg.fit(X, y.ravel())
predict2 = sgd_reg.predict([[1.5]])
print("The predict is ", predict2)

输出结果

The predict is  [[4.92157355]]
The predict is  [4.11126873]

2、LASSO回归

LASSO回归（也称 Least Absolute Shrinkage，或 Selection Operator Regression）：其正则化项是权重向量的 $l1$ 范数

（1）损失函数

$$J(\theta )=MSE(\theta )+\alpha \sum _{i=1}^n|{\theta }_i|$$

下图为相同线性数据上使用不同的 $\alpha$ 值的LASSO回归模型最后的表现。其中，左图为简单的LASSO回归模型，右图为10阶的 PolynomialFeatures 进行扩展。

分析：倾向于完全消除最不重要的特征的权重（置为0）
特性：自动进行特征选择同时输出一个稀疏模型。

（2）对比

其中，$L_1$ 范数正则化、$L_2$ 范数正则化都有助于降低过拟合风险，$L_2$ 范数通过对参数向量各元素平方和求平方根，使得范数最小，从而使得参数的各个元素接近0 ，但不等于0。 而 $L_1$ 范数正则化比范数更易获得“稀疏”解，即 $L_1$ 范数正则化求得的会有更少的非零分量，所以 $L_1$ 范数可用于特征选择，而 $L_2$ 范数在参数规则化时经常用到（事实上， $L_0$ 范数得到的“稀疏”解最多，但 $L_0$ 范数$||x||=\#(i|x_i̸=0)$ 是 $x$ 中非零元素的个数，不连续，难以优化求解。因此常用 $L_1$ 范数来近似代替）。

为什么 $L_1$ 正则化更易获得“稀疏”解呢？

假设仅有两个属性，$\omega$ 只有两个参数${\omega }_1$, ${\omega }_2$，绘制不带正则项的目标函数-平方误差项等值线，再绘制 $L_1$，$L_2$ 范数等值线，如下图所示正则化后优化目标的解要在平方误差项和正则化项之间折中，即出现在图中等值线相交处采用。$L_1$ 范数时，交点常出现在坐标轴上，即 ${\omega }_1$ 或 ${\omega }_2$为 0 ;而采用 $L_2$ 范数时，交点常出现在某个象限中，即 ${\omega }_1$, ${\omega }_2$均非0。也就是说，$L_1$ 范数比 $L_2$ 范数更易获得“稀疏”解。

注意：在Lasso损失函数中，批量梯度下降得路径趋向在低谷有一个反弹。这是因为在 ${\omega }_1$ 时，斜率会有一个突变。为了最后收敛到全局最小值，需要降低学习率。

Lasso损失函数在 ${\theta }_i=0(i=1,2,...,n)$ 处无法进行微分运算，但是使用子梯度向量 $g$ 后，它可以在任何 ${\theta }_i=0$ 的情况下进行计算。如下式所示。
Lasso回归子梯度向量：
$$g(\theta ,J)={\nabla }_{\theta }MSE(\theta )+\alpha \left(\begin{array}{c}sign({\theta }_1)\\ sign({\theta }_2)\\ ⋮\\ sign({\theta }_n)\end{array}\right)wheresign({\theta }_i)=\{\begin{array}{l}-1,{\theta }_i<0\\ 0,{\theta }_i=0\\ +1,{\theta }_i>0\end{array}$$

from sklearn.linear_model import Lasso

lasso_reg = Lasso(alpha=0.1)
lasso_reg.fit(X, y)
predict3 = lasso_reg.predict([[1.5]])
print("The predict of Lasso is ", predict3)

输出结果

The predict of Lasso is  [4.89631518]

或可以使用 SGDRegressor(penalty="l1")

3、ElasticNet（弹性网络）

正则项是 Ridge回归 和 Lasso回归正则项的简单混合，同时可以控制它们的混合率 $r$。当 $r=0$ 时，ElasticNet就是Ridge回归，当 $r=1$ 时，其就是 Lasso 回归。

（1）损失函数

$$J(\theta )=MSE(\theta )+r\alpha \sum _{i=1}^n|{\theta }_i|+\frac{1-r}{2}\alpha \sum _{i=1}^n{\theta }_i^2$$

（2）实现

l1_ratio 就是混合率 $r$

from sklearn.linear_model import ElasticNet

elastic_net = ElasticNet(alpha=0.1, l1_ratio=0.5)
elastic_net.fit(X, y)
predict4 = elastic_net.predict([[1.5]])
print("The predict of ElasticNet is ", predict4)

结果输出

The predict of ElasticNet is  [4.895499]

4、四种模型的选择

三种模型：线性回归、岭回归、Lasso回归、弹性网络。

• 避免使用简单的线性回归；
• 首选：岭回归；
• 若特征仅有少数是真正起作用的，则选择：Lasso和弹性网络；==》可以将无用特征的权重降为零
  • 一般来说，弹性网络的表现要比Lasso好。
  • 原因：当特征数量比样本的数量大的时候，或特征之间有很强的相关性时，Lasso可能会表现得不规律。

5、早期停止法（ Early Stopping）

（1）概述

随着迭代次数（epoch）的增加，算法在训练集上的预测误差（RMSE）逐渐下降，而对于验证误差开始下降之后开始上升。这意味着模型在训练集上出现过拟合。
==》一旦验证错误达到最小值，便提早停止训练。
==》“完美的免费午餐”（Geoffrey Hinton）

注意：随机梯度和mini-batch梯度下降不是平滑曲线，如何找到最小值。==》解决方案： 只有在验证误差高于最小值一段时间后才停止，之后将模型参数回滚到验证误差最小值。

（2）实现

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

from sklearn.base import clone

sgd_reg = SGDRegressor(n_iter=1, warm_start=True, penalty=None,
                       learning_rate="constant", eta0=0.0005)
mininum_val_error = float("inf")
best_epoch = None
best_model = None
for epoch in range(1000):
    sgd_reg.fit(X_train_ploy_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 < mininum_val_error:
        mininum_val_error = val_error
        best_epoch = epoch
        best_model = clone(sgd_reg)

四、逻辑回归

Logistic回归（也称 Logit 回归），通常用于估计一个实例属于某个特定类别的概率，并根据sigoid函数划分类别（二分类）。

1、概率估计

Logistic回归模型计算输入特征的加权和（加上偏差项），之后将中间结果输入 logistic() 函数进行二次加工后输出。

逻辑回归模型的概率估计（向量形式）
$$\hat{p}=h_{\theta }(x)=\sigma ({\theta }^T\cdot x)$$
其中，$\sigma ()$ 表示 logistic 函数（也称 logit 函数）
$$\sigma (t)=\frac{1}{1+e^{-t}}$$

一旦 logistic 回归模型估计得到了 $X$ 属于正类的概率 $\hat{p}=h_{\theta }(x)$，则很容易得到预测结果 $\hat{y}$。

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

2、训练和损失函数

注：该部分以手写笔记为主

逻辑回归的损失函数（对数损失）
$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m[y^{(i)}log(h_{\theta }(x_i))+(1-y^{(i)})log(1-h_{\theta }(x_i))]$$

损失函数关于第 $j$ 个模型参数 ${\theta }_j$ 的偏导数
$$\frac{\partial }{\partial {\theta }_j}J({\theta }_j)=\frac{1}{m}\sum _{i=1}^m(\sigma ({\theta }^T\cdot X^{(i)})-y^{(i)})x_j^{(i)}$$

3、决策边界

对鸢尾花数据进行 logistic 回归。
鸢尾花数据：150朵三种不同的鸢尾花的萼片和花瓣的长度和宽度，类别为Setosa、Versicolor、Virginica

目标：建立一个分类器，仅使用花瓣的宽度特征识别Virginica。

from sklearn import datasets
from sklearn.linear_model import LogisticRegression

## 建立仅使用花瓣的宽度特征来识别
iris = datasets.load_iris()
print(list(iris.keys()))
X = iris["data"][:, 3:]  ## 只读最后一个特征
y = (iris["target"] == 2).astype(np.int)  # 取出判断是否为第3类的label

## 训练
log_reg = LogisticRegression()
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)
plt.plot(X_new, y_proba[:, 1], "g-", label="Iris-Virginica")
plt.plot(X_new, y_proba[:, 0], "b--", label="Not Iris-Virginica")
plt.axis([0, 3, 0, 1])
plt.legend()
plt.show()

五、Softmax回归

Logistic回归扩展到多分类==》Softmax回归

1、原理

思路：给定实例 $x$ 时，Softmax回归模型先计算 $k$ 类的分数 $s_k(x)$，然后将分数应用到 Softmax函数，估计出每类的概率。将估计概率最高（它只是得分最高的类）的类别作为预测结果。

softmax函数：估计样本属于第 $k$ 类的概率 $${\hat{p}}_k=\sigma (s(x){)}_k=\frac{exp(s_k(x))}{{\sum }_{j=1}^{K}exp(s_j(x))}$$
其中：
$$s_k(x)={\theta }_k^T\cdot x$$
其中，$K$ 表示类别的个数；需要注意的是 ${\theta }_k$ 说明每个类别对应有自己的 $\theta$，所有 ${\theta }_k$ 组合起来是全部的参数。$s(x)$ 表示包含样本 $x$ 每一类得分的向量；$\sigma (s(x){)}_k$表示给定每一类分数之后，实例 $x$ 属于第 $k$ 类的概率。

注意：softmax回归分类器一次只能预测一个类（智能用于互斥的类别，不能用于多标签）

2、训练

目标：建立一个模型在目标类别上有较高的概率。
==》损失函数：交叉熵（衡量待测类别与目标类别的匹配程度）
$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m\sum _{k=1}^K{y}_k^{(i)}log({\hat{p}}_k^{(i)})$$
若对于第 $i$ 个实例的目标类是 $k$，那么 $y_k^{(i)}=1$， 反之 $y_k^{(i)}=0$。

计算每一类的梯度向量，然后采用梯度下降（或其他优化算法）找到使损失函数达到最小值的参数矩阵 $\theta$。
$k$ 类交叉熵关于 ${\theta }_k$ 的梯度向量：
$${\nabla }_{{\theta }_k}J(\theta )=\frac{1}{m}\sum _{i=1}^m({\hat{p}}_k^{(i)}-y_k^{(i)})x^{(i)}$$

from sklearn.linear_model import LogisticRegression
X = iris["data"][:,(2, 3)] # petal length, petal width
y = iris["target"]

# multinomial设置为softmax回归，lbfgs求解器
softmax_reg = LogisticRegression(multi_class="multinomial", solver="lbfgs", C=10)
softmax_reg.fit(X, y)
predict = softmax_reg.predict([[5, 2]])
predict_proba  = softmax_reg.predict_proba([[5, 2]])
print("The predict is ", predict)
print("The predict_proba is ", predict_proba)

输出结果

The predict is  [2]
The predict_proba is  [[6.33134077e-07 5.75276067e-02 9.42471760e-01]]

回答 94.2% 是 Virginica 花（ 第二类） ， 或者 5.8% 是其他鸢尾花。