chapter4 training_linear_models

设置

首先，确保这个笔记本在python 2和python 3中都能正常工作，导入一些公共模块，确保MatplotLib内联绘制图形，并准备一个函数来保存这些图形:

from __future__ import division,print_function,unicode_iterals
 
import numpy as np
np.random.seed(42)
 
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
 
import os
PROJECT_ROOT_DIR = "."
CHAPTER_ID = "training_linear_models"
def save_fig(fig_id,tight_layout = True):
    path = os.join.path(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)
 
import warnings
warnings.filterwarnings(action = "ignore",module = "scipy",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()

用来生成数据的函数是y = 4 + 3x0 + 高斯噪声。看一下最后的计算结果。 

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)
theta_best

希望最后得到的参数为θ0 = 4，θ1 = 3而不是 θ0 = 4.215，θ1 = 2.770。这已经足够了，由于存在噪声，参数不可能达到到原始函数的值。

接下来进行预测

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)
y_predict

plt.plot(X_new,y_predict,"r-")
Plt.plot(X,y,"b.")
plt.axis([0,2,0,15])
plt.show()

plt.plot(X_new,y_predict,"r-",linewidth = 2,labei = "Predictions")
PLT.Plot(X,y,"b.")
plt.xlabel("$x_1$",fontsize = 18)
plt.ylabel("$y$,rotation = 0,fontsize = 18")
plt.legend(loc = "upper left",fontsize = 14)
plt.axis([0,2,0,15])
save_fig("linear_model_predictions")
plt.show()

使用下面的 Scikit-Learn 代码可以达到相同的效果：

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

lin_reg.predict(X_new)

LinearRegression类是基于scipy.linalg.lstsq()函数(名称代表“最小二乘”)，可以直接调用它:

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

这个函数计算,是X的伪逆(特别是Moore-Penrose逆)。可以使用np.linalg.pinv()来直接计算伪逆:

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

注意:这本书的第一版暗示LinearRegressor类是基于标准方程。这是一个错误。:如上所述,它是基于伪逆,最终依赖于计算矩阵分解𝐗(有关详细信息,请参阅第八章的计算分解)。其时间复杂度𝑂(𝑛2)和它的工作原理,即使𝑚<𝑛或者当某个功能是其他特性的线性组合(在这种情况下,不是可逆的正规方程失败),查看更多细节问题# 184。然而,这并不改变其他LinearRegression类的描述,特别是,它是基于一个解析解,它不规模与数量的特性,它扩展线性的实例数量,所有的数据必须装入内存,它不需要功能扩展和实例的训练集的顺序无关紧要。

用批量梯度下降处理线性回归

损失函数的梯度向量：

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

梯度下降步长：

eta = 0.1
n_iterations = 1000
m = 100
theta = np.random.randn(2,1)
for iteration in range(n_iterations):
    gradient = 2/m * X_b.T.dot(X_b.dot(theta) - y)
    theta = theta - eta * gradients
theta

X_new_b.dot(theta)

看！正态方程的表现非常好。完美地求出了梯度下降的参数。但是当换一个学习率会发生什么？

np.random.seed(42)
theta = np.random.randn(2,1)
 
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(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)    
 
 
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()

上图展示了使用了三个不同的学习率进行梯度下降的前 10 步运算（虚线代表起始位置）。

在左面的那副图中，学习率是最小的，算法几乎不能求出最后的结果，而且还会花费大量时间。在中间的这幅图中，学习率的表现看起来不错，仅仅几次迭代后，它就收敛到了最后的结果。在右面的那副图中，学习率太大了，算法是发散的，跳过了所有的训练样本，同时每一步都离正确的结果越来越远。

用随机梯度下降处理线性回归

随机梯度下降，在每一步的梯度计算上只随机选取训练集中的一个样本。很明显，由于每一次的操作都使用了非常少的数据，这样使得算法变得非常快。由于每一次迭代，只需要在内存中有一个实例，这使随机梯度算法可以在大规模训练集上使用。

另一方面，由于它的随机性，与批量梯度下降相比，其呈现出更多的不规律性：它到达最小值不是平缓的下降，损失函数会忽高忽低，只是在大体上呈下降趋势。随着时间的推移，它会非常的靠近最小值，但是它不会停止在一个值上，它会一直在这个值附近摆动。因此，当算法停止的时候，最后的参数还不错，但不是最优值。

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

theta_path_sgd = []
m = len(X_b)
np.random.seed(42)
 
 
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
 
 
n_epochs = 50
t0,t1 = 5,50
def learning_schedule(t):
    return t0 / (t + t1)
theta = np.random.randn(2,1)
 
for epoch in range(n_epochs):
    for i in range(m):
        if epoch == 0 and i < 20:
            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]
        gradient = 2 * xi.T.dot(xi.dot(theta) - yi)
        theta = theta - eta * gradients
        theta_path_sgd.append(theta)
theta)

按习惯来讲，进行m轮的迭代，每一轮迭代被称为一代。 

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

上图展示了前 20 次的训练过程（注意每一步的不规则程度）。

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

from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter = 50,penalty = None,eta0 = 0.1,random_state = 42)
sgd_reg.fit(X,y.ravel())

sgd_reg.intercept_,sgd_reg.coef_

可以再一次发现，这个结果非常的接近正态方程的解。

用小批量梯度下降处理线性回归

在迭代的每一步，批量梯度使用整个训练集，随机梯度时候用仅仅一个实例，在小批量梯度下降中，它则使用一个随机的小型实例集。它比随机梯度的主要优点在于你可以通过矩阵运算的硬件优化得到一个较好的训练效果，尤其当你使用 GPU 进行运算的时候。

小批量梯度下降在参数空间上的表现比随机梯度下降要好的多，尤其在有大量的小型实例集时。作为结果，小批量梯度下降会比随机梯度更靠近最小值。但是，另一方面，它有可能陷在局部最小值中（在遇到局部最小值问题的情况下，和我们之前看到的线性回归不一样）。

EPOCH, BATCH, INTERATION

CIFAR10 数据集有 50000 张训练图片，10000 张测试图片。现在选择 Batch Size = 256 对模型进行训练。

每个 Epoch 要训练的图片数量： 50000

训练集具有的 Batch 个数： 50000 / 256 = 195 + 1 = 196

每个 Epoch 需要完成的 Batch 个数： 196

每个 Epoch 具有的 Iteration 个数： 196

每个 Epoch 中发生模型权重更新的次数： 196

训练 10 代后，模型权重更新的次数： 196 * 10 = 1960

不同代的训练，其实用的是同一个训练集的数据。第 1 代和第 10 代虽然用的都是训练集的五万张图片，但是对模型的权重更新值却是完全不同的。因为不同代的模型处于代价函数空间上的不同位置，模型的训练代越靠后，越接近谷底，其代价越小

EPOCH:训练样本全部跑一遍 就是一个EPOCH

BATCH SIZE:就一个BATCH有多少个样本

iteration:迭代，只需要知道乘法表或者一个计算器就可以了。迭代是 batch 需要完成一个 epoch 的次数。记住：在一个 epoch 中，batch 数和迭代数是相等的。

比如对于一个有 2000 个训练样本的数据集。将 2000 个样本分成大小为 500 的 batch，那么完成一个 epoch 需要 4 个 iteration。

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_interations):
    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 = t + 1
        xi = X_b_shuffled[i:i + minibatch_size]
        yi = y_shuffled[i:i + minibatch_size]
        gradient = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)
        theta = theta - eta * gradients
        theta_path_mgd.append(theta)
theta

 

theta_path_bgd = np.array(theta_path_bgd)
theta_path_sgd = np.array(theta_path_sgd)
theta_path_mgd = np.array(theta_path_mgd)
 
 
plt.figure(figsize = (7,4))
#plt.plot(X,y,"b.")
plt.plot(theta_path_sgd[:,0],theta_path_sgd[:,1],"r-s",linewidth = 1,label = "Stochastic")
plt.plot(theta_path_mgd[:,0],theta_path_mgd[:,1],"g-+",linewidth = 1,label = "Mini_batch")
plt.plot(theta_path_bgd[:,0],theta_path_bgd[:,1],"b-o",linewidth = 1,label = "Batch")
plt.legend(loc = "upper left",fontsize = 16)
 
plt.xlabel(r"$\theta_0$",fontsize = 20)
plt.ylabel(r"$\theta_1$",rotation = 0,fontsize = 20)
plt.axis([2.5,4.5,2.3,3.9])
save_fig("gradient_descent_paths_plot")
plt.show()

上图显示了训练期间三种梯度下降算法在参数空间中所采用的路径。 它们都接近最小值，但批量梯度的路径最后停在了最小值，而随机梯度和小批量梯度最后都在最小值附近摆动。 但是，不要忘记，批量梯度需要花费大量时间来完成每一步，但是，如果使用了一个较好的learning schedule，随机梯度和小批量梯度也可以得到最小值。

多项式回归

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

看一个例子。 首先，根据一个简单的二次方程（并加上一些噪声）来生成一些非线性数据：

import numpy as np
import numpy.random as rnd
np.random.seed(42)
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 = ploy_features.fit_transform(X)
X[0],X_poly[0]

,

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

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

X_new = np.linspace(-3,3,100).reshape(100,1)
X_new_poly = poly.feature.transform(X_new)
y_new = lin_reg.predict(X_new_poly)
 
plt.plot(X,y,"b.")
plt.plot(X_new,y_new,"r-",linewidth = 2,label = "Predictions")
plt.xlabel("$x_1$",fontsize = 18)
plt.ylabel("$y$",rotation = 0,fontsize = 18)
plt.legend(loc = "upper left",fontsize = 14)
plt.axis([-3,3,0,10])
save_fig("quadratic_predictions_plot")
plt.show()

还是不错的，模型预测函数是，实际的原始函数是。

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

提示
PolynomialFeatures(degree=d)把一个包含n个特征的数组转换为一个包含(n+d)!/d!n!个特征的数组，n!是n的阶乘，等于1 × 2 × 3 × ⋯ × n。小心大量特征的组合爆炸！

如果使用一个高阶的多项式回归，可能发现它的拟合程度要比普通的线性回归要好的多。

from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
for style,width degree in (("g-",1,300),("b--",2,2),("r-+",2,1)) :
    polybig_features = PolynomialFeatures(degree = degree,include_bias = False)
    std_scaler = StandardScaler()
    lin_reg = LinearRegression()
    polynomial_regression = Pipeline([
        ("poly_features",poly_features),
        ("std_sccaler",std_scaler),
        ("lin_reg",lin_reg),
])
    polynomial_regression.fit(X,y)
    y_newbig = polynomial_regression.predict(X_new)
    plt.plot(X_new,y_newbig,style,label = str(degree),linewidth = width)
 
plt.plot(X,y,"b.",linewidth = 3)
plt.legend(loc = "upper left")
plt.xlabel("$x_1$",fontsize = 18)
plt.ylabel("$y$",rotation = 0,fontsize = 18)
plt.axis([-3,3,0,10])
save_fig("high_degree_polynomials_plot")
plt.show()

上图使用一个 300 阶的多项式模型去拟合之前的数据集，并同简单线性回归、2 阶的多项式回归进行比较。注意 300 阶的多项式模型如何摆动以尽可能接近训练实例。

当然，这种高阶多项式回归模型在这个训练集上严重过拟合了，线性模型则欠拟合。在这个训练集上，二次模型有着较好的泛化能力。那是因为在生成数据时使用了二次模型，但是一般我们不知道这个数据生成函数是什么，那该如何决定我们模型的复杂度呢？如何区分模型是过拟合还是欠拟合？

在第二章，可以使用交叉验证来估计一个模型的泛化能力。如果一个模型在训练集上表现良好，通过交叉验证指标却得出其泛化能力很差，那么你的模型就是过拟合了。如果在这两方面都表现不好，那么它就是欠拟合了。

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

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_ratio = 0.2,random_state = 42)
    train_errors,test_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_train_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 = LinearRegresson()
plot_learning_curves(lin_reg,X,y)
plt.axis([0,80,0,3])
save_fig("underfitting_learning_curves_plot")
plt.show()

这幅图值得深究。首先，观察在训练集上的效果：当训练集只有一两个样本的时候，模型能够非常好的拟合它们，这也是为什么曲线是从零开始的原因。但是当加入了一些新的样本的时候，训练集上的拟合程度变得难以接受，出现这种情况有两个原因，一是因为数据中含有噪声，另一个是数据根本不是线性的。因此随着数据规模的增大，误差也会一直增大，直到达到高原地带并趋于稳定，在之后，继续加入新的样本，模型的平均误差不会变得更好或者更差。继续来看模型在验证集上的表现，当以非常少的样本去训练时，模型不能恰当的泛化，也就是为什么验证误差一开始是非常大的。当训练样本变多的到时候，模型学习的东西变多，验证误差开始缓慢的下降。但是一条直线不可能很好的拟合这些数据，因此最后误差会到达在一个高原地带并趋于稳定，最后和训练集的曲线非常接近。

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

提示

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

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

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

 

这幅图像和之前的有一点点像，但是其有两个非常重要的不同点：

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

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

正则化模型（线性模型的正则化）

降低模型的过拟合的好方法是正则化这个模型（即限制它）：模型有越少的自由度，就越难以拟合数据。例如，正则化一个多项式模型，一个简单的方法就是减少多项式的阶数。

对于一个线性模型，正则化的典型实现就是约束模型中参数的权重。 接下来我们将介绍三种不同约束权重的方法：岭回归，Lasso回归和弹性网络（Elastic Net）。

岭（Ridge）回归

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

超参数α决定了你想正则化这个模型的强度。如果α=0那此时的岭回归便变为了线性回归。如果α非常的大，所有的权重最后都接近于零，最后结果将是一条穿过数据平均值的水平直线。公式 4-8 是岭回归的损失函数：

岭回归损失函数：

值得注意的是偏差项是没有被正则化的（累加运算的开始是i=1而不是i=0。

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

from sklearn.linear_model import Ridge
 
np.random.seed(42)
m = 20
X = 3 *np.random.rand(m,1)
y = 1 + 0.5 * X + np.random.randn(m,1) / 1.5
X_new = np.linspace(0,3,100).reshape(100,1)
def plot_model(model_class,polynomial,alphas,**model_kargs):
    for alpha,style in zip(alphas,("b-","g--","r:")):
        model = model_class(alpha,**model_kargs) if alpha > 0 else LinearRegression()
        if poltnomial:
            model = Pipeline([
                ("poly_features",PolynomialFeatures(degree = 10,include_bias = False)),
                ("std_scaler",StandardScaler()),
                ("regular_reg",model),
])
        model.fit(X,y)
        y_new_regul = model.predict(X_new)
        lw = 2 if alpha > 0 else 1
        plt.plot(X_new,y_new_regul,style,linewidth = lw,label = r"$\alpha = {}$".format(alpha))
    plt.plot(X,y,"b.",linewidth = 3)
    plt.legend(loc = "upper left",fontsize = 15)
    plt.xlabel("$x_1$",fontsize = 18)
    plt.axis([0,3,0,4])
 
plt.figure(figsize = (8,4))
plt.subplot(121)
plt_model(Ridge,polynomial = False,alphas = (0,10,100),random_state = 42)
plt.ylabel("$y$",rotation = 0,fontsize = 18)
plt.subplot(122)
plt_model(Ridge,polynomial = True,alphas = (0,10**-5,1),random_state = 42)
 
save_fig("ridge_regression_plot")
plt.show()

上图展示了在相同线性数据上使用不同α值的岭回归模型最后的表现。左图中，使用简单的岭回归模型，最后得到了线性的预测。右图中的数据首先使用 10 阶的PolynomialFearures进行扩展，然后使用StandardScaler进行缩放，最后将岭模型应用在处理过后的特征上。这就是带有岭正则项的多项式回归。注意当α增大的时候，导致预测曲线变得扁平（即少了极端值，多了一般值），这样减少了模型的方差，却增加了模型的偏差。

对线性回归来说，可以使用封闭方程去计算，也可以使用梯度下降来做岭回归。它们的缺点和优点是一样的。下面的公式表示封闭方程的解（矩阵A 是一个除了左上角有一个0的n × n 单位矩阵，这个0代表偏差项。）

下面是如何使用 Scikit-Learn 来进行封闭方程的求解（使用 Cholesky 法）:

from sklearn.linear_model import Ridge
ridge_reg = Ridge(alpha = 1,solver = "cholesky",random_state = 42)
ridge_reg.fit(X,y)
rideg_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]])

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

ridge_reg = Ridge(alpha = 1,solver = "sag",random_state = 42)
ridge_reg.fit(X,y)
ridge_reg.predict([[1.5]])

Lasso 回归

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

Lasso 回归的损失函数：

from sklearn.linear_model import Lasso
 
plt.figure(figsize = (8,4))
plt.subplot(121)
plt_model(Lasso,polynomial = False,alphas = (0,0.1,1),random_state = 42)
plt.ylabel("$y$",rotation = 0,fontsize = 18)
plt.subplot(122)
plt_model(Lasso,polynomial = True,alphas = (0,10**-7,1),tol = 1,random_state = 42)
 
save_fig("lasso_regression_plot")
plt.show()

 

上图展示了和前边的图相同的事情，仅仅是用 Lasso 模型代替了 Ridge 模型，同时调小了α的值。

Lasso 回归的一个重要特征是它倾向于完全消除最不重要的特征的权重（即将它们设置为零）。例如，右图中的虚线所示（），曲线看起来像一条二次曲线，而且几乎是线性的，这是因为所有的高阶多项特征都被设置为零。换句话说，Lasso回归自动的进行特征选择同时输出一个稀疏模型（即，具有很少的非零权重）。

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

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

弹性网络（ElasticNet）

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

弹性网络损失函数：

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

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

早期停止法（Early Stopping）

对于迭代学习算法，有一种非常特殊的正则化方法，就像梯度下降在验证错误达到最小值时立即停止训练那样。我们称为早期停止法。

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 = 42)
 
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.fit_transform(X_val)
 
sgd_reg = SGDRegressor(max_iter = 1,tol = -np.infty,penalty = Npne,eta0 = 0.0005,
                        warm_start = None,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),)
 
best_val_rmse -= 0.03
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()

上图表示使用随机梯度下降来训练一个非常复杂的模型（一个高阶多项式回归模型）。随着训练的进行，算法一直学习，它在训练集上的预测误差（RMSE）自然而然的下降。然而一段时间后，验证误差停止下降，并开始上升。这意味着模型在训练集上开始出现过拟合。一旦验证错误达到最小值，便提早停止训练。

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

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

from sklearn.base import clone
sgd_reg = SGDRegressor(max_iter = 1,tol = -np.infty,warm_start = True,penalty = None,
            learning_rate = "constant",tea0 = 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)
    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

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

下面为Ridge 回归和 Lasso 回归对比

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
 
 
t1a,t1b,t2a,t2b = -1,3,-1.5,1.5
 
t1s = np.linspace(t1a,t1b,500)
t2s = np.linspace(t2a,t2b,500)
t1,t2 = np.meshgrid(t1s,t2s)
T = np.c_[t1.ravel(),t2.ravel()]
Xr = np.array([[-1,1],[-0.3,-1],[1,0.1]])
yr = 2 * Xr[:,:1] + 0.5 * Xr[:,1:]
 
J = (1/len(Xr)) * np.sum((T.dot(Xr.T) - yr.T)**2,axis = 1)
 
N1 = np.linalg.norm(T,ord = 1,axis = 1).reshape(t1.shape)
N2 = np.linalg.norm(T,ord = 2,axis = 1).reshape(t1.shape)
 
t_min_idx = np.unravel_index(np.argmin(J), J.shape)
t1_min, t2_min = t1[t_min_idx], t2[t_min_idx]
 
t_init = np.array([[0.25], [-1]])
 
 
 
def bgd_path(theta, X, y, l1, l2, core = 1, eta = 0.1, n_iterations = 50):
    path = [theta]
    for iteration in range(n_iterations):
        gradients = core * 2/len(X) * X.T.dot(X.dot(theta) - y) + l1 * np.sign(theta) + 2 * l2 * theta
 
        theta = theta - eta * gradients
        path.append(theta)
    return np.array(path)
 
plt.figure(figsize=(12, 8))
for i, N, l1, l2, title in ((0, N1, 0.5, 0, "Lasso"), (1, N2, 0,  0.1, "Ridge")):
    JR = J + l1 * N1 + l2 * N2**2
    
    tr_min_idx = np.unravel_index(np.argmin(JR), JR.shape)
    t1r_min, t2r_min = t1[tr_min_idx], t2[tr_min_idx]
 
    levelsJ=(np.exp(np.linspace(0, 1, 20)) - 1) * (np.max(J) - np.min(J)) + np.min(J)
    levelsJR=(np.exp(np.linspace(0, 1, 20)) - 1) * (np.max(JR) - np.min(JR)) + np.min(JR)
    levelsN=np.linspace(0, np.max(N), 10)
    
    path_J = bgd_path(t_init, Xr, yr, l1=0, l2=0)
    path_JR = bgd_path(t_init, Xr, yr, l1, l2)
    path_N = bgd_path(t_init, Xr, yr, np.sign(l1)/3, np.sign(l2), core=0)
 
    plt.subplot(221 + i * 2)
    plt.grid(True)
    plt.axhline(y=0, color='k')
    plt.axvline(x=0, color='k')
    plt.contourf(t1, t2, J, levels=levelsJ, alpha=0.9)
    plt.contour(t1, t2, N, levels=levelsN)
    plt.plot(path_J[:, 0], path_J[:, 1], "w-o")
    plt.plot(path_N[:, 0], path_N[:, 1], "y-^")
    plt.plot(t1_min, t2_min, "rs")
    plt.title(r"$\ell_{}$ penalty".format(i + 1), fontsize=16)
    plt.axis([t1a, t1b, t2a, t2b])
    if i == 1:
        plt.xlabel(r"$\theta_1$", fontsize=20)
    plt.ylabel(r"$\theta_2$", fontsize=20, rotation=0)
 
    plt.subplot(222 + i * 2)
    plt.grid(True)
    plt.axhline(y=0, color='k')
    plt.axvline(x=0, color='k')
    plt.contourf(t1, t2, JR, levels=levelsJR, alpha=0.9)
    plt.plot(path_JR[:, 0], path_JR[:, 1], "w-o")
    plt.plot(t1r_min, t2r_min, "rs")
    plt.title(title, fontsize=16)
    plt.axis([t1a, t1b, t2a, t2b])
    if i == 1:
        plt.xlabel(r"$\theta_1$", fontsize=20)
 
save_fig("lasso_vs_ridge_plot")
plt.show()

可以从上图知道为什么会出现这种情况：在左上角图中，后背景的等高线（椭圆）表示了没有正则化的均方差损失函数，白色的小圆圈表示在当前损失函数上批量梯度下降的路径。前背景的等高线（菱形）表示惩罚，黄色的三角形表示了仅在这个惩罚下批量梯度下降的路径（无穷）。注意路径第一次是如何到达，然后向下滚动直到它到达。在右上角图中，等高线表示的是相同损失函数再加上一个的惩罚。左下角的图中，它的全局最小值在这根轴上。批量梯度下降首先到达，然后向下滚动直到达到全局最小值。 两个底部图显示了相同的情况，只是使用了惩罚。 规则化的最小值比非规范化的最小值更接近于，但权重不能完全消除。

逻辑回归

正如之前讨论的那样，一些回归算法也可以用于分类（反之亦然）。 Logistic 回归（也称为 Logit 回归）通常用于估计一个实例属于某个特定类别的概率（例如，这电子邮件是垃圾邮件的概率是多少？）。 如果估计的概率大于 50%，那么模型预测这个实例属于当前类（称为正类，标记为“1”），反之预测它不属于当前类（即它属于负类 ，标记为“0”）。 这样便成为了一个二元分类器。

概率估计

那么它是怎样工作的？ 就像线性回归模型一样，Logistic 回归模型计算输入特征的加权和（加上偏差项），但它不像线性回归模型那样直接输出结果，而是把结果输入logistic()函数进行二次加工后进行输出（详见公式 4-13）。

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

逻辑函数：

训练和损失函数

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

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

t = np.linspace(-10,10,100)
sig = 1 / (1 + np.exp(-t))
plt.figure(fontsize = (9,3))
plt.plot([-10,10],[0,0],"k-")
plt.plot([-10,10],[0.5,0.5],"k:")
plt.plot([-10,10],[1,1],"k:")
plt.plot([0,0],[-1.1,1.1],"k-")
plt.plot(t,sig,"b-",linewwidth = 2,label = r"$\sigma(t) = \frac{1 + e ^ {-t}}$")
plt.xlabel("t")
plt.legend(loc = "upper left",fontsize = 20)
lpt.axis([-10,10,-0.1,1.1])
save_fig("logistic_function_plot")
plt.show()

决策边界

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

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

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

print(iris.DESR)

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

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

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

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.axis([0,3,-0.02,1.02])
save_fig("Logistic_regression_plot")
plt.show()
 

decision_boundary

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

log_reg.predict([1.7],[1.5])

下面用相同的数据集，但是这次使用了两个特征进行判断：花瓣的宽度和长度。 一旦训练完毕，Logistic 回归分类器就可以根据这两个特征来估计一朵花是 Virginica 的可能性。

from sklearn.linear_model import LogisticRegression
 
X = iris["data"][:, (2, 3)] 
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()
plt.plot(X_new,y_proba[:,0],"b--",linewidth = 2,label = "Not Iris-Virginica")

虚线表示这时两类情况出现的概率都等于 50%：这是模型的决策边界。 请注意，它是一个线性边界。每条平行线都代表一个分类标准下的两两个不同类的概率，从 15%（左下角）到 90%（右上角）。越过右上角分界线的点都有超过 90% 的概率是 Virginica 花。

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

 

注意

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

Softmax 回归

Logistic 回归模型可以直接推广到支持多类别分类，不必组合和训练多个二分类器（如第 3 章所述）， 其称为 Softmax 回归或多类别 Logistic 回归。

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

Softmax函数：

• 是类的个数
• 是实例的每个类的得分的向量
• 表示给定每一类分数之后，实例属于第类的概率

Softmax 回归模型分类器预测结果:

argmax运算返回一个函数取到最大值的变量值。 在这个等式，它返回使最大时的的值

提示

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

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

X = iris["data"][:,(2,3)]
y = iris["target"]

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

softmax_reg.predict([[5,2]])

softmax_reg.predict_proba([[5,2]])

可以看到一个花瓣长为 5 厘米，宽为 2 厘米的鸢尾花时，可以问模型它是哪一类鸢尾花，它会回答 94.2% 是 Virginica 花（第二类），或者 5.8% 是其他鸢尾花。

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.contuorf(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("sofemax_regression_contour_plot")
plt.show()

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

练习题
在 Softmax 回归上应用批量梯度下降的早期停止法（不使用 Scikit-Learn）

从加载数据开始。将重新使用之前加载的Iris数据集。

X = iris["data"][:,(2,3)]
y = iris["target"]
X[:5]

需要为每个实例添加偏差项():

X_with_bias = np.c_[np.ones([len(X),1]),X]
X_with_bias[:5]

设置随机种子，这样这个练习的结果是可复制的。

np.random.seed(2042)

将数据集分割为训练集、验证集和测试集的最简单方法是使用Scikit-Learn的train_test_split()函数，但是本练习的目的是通过手动实现来尝试理解算法。这里有一个可能的实现:

test_ratio = 0.2
validation_ratio = 0.2
total_size = len(X_with_bias)

test_size = int(total_size * test_ratio)
validation_size = int(total_size * validation_ratio)
train_size = total_size - test_size - validation_size

rnd_indices = np.random.permutation(total_size)

X_train = X_with_bias[rnd_indices[:train_size]]
y_train = y[rnd_indices[:train_size]]
X_valid = X_with_bias[rnd_indices[train_size:-test_size]]
y_valid = y[rnd_indices[train_size:-test_size]]
X_test = X_with_bias[rnd_indices[-test_size:]]
y_test = y[rnd_indices[-test_size:]]

目前的目标是类指标(0、1或2)，但需要目标类概率来训练Softmax回归模型。每个实例的目标类概率对于所有类来说都是0.0，除了目标类的概率是1.0(换句话说，给定实例的类概率的向量是一个热向量)。写一个小函数转换成一个矩阵的类，该矩阵指标的向量对于每个实例包含一个热向量:

def to_one_hot(y):
    n_classes = y.max() + 1 #2+1
    m = len(y) #3
    Y_one_hot = np.zeros((m,n_classes))
    Y_one_hot[np.arange(m),y] = 1
    return y_one_hot

让我们在前10个实例上测试这个函数:

y_train[:10]

to_one_hot(y_train[:10])

看起来不错，接下来为训练集和测试集创建目标类概率矩阵:

Y_train_one_hot = to_one_hot(y_train)
Y_valid_one_hot = to_one_hot(y_valid)
Y_test_one_hot = to_one_hot(y_test)

现在让我们实现Softmax函数。回忆一下，它是由以下方程定义的:

def softmax(logists):
    exps = np.exp(logists)
    exp_sums = np.sum(exps,axis = 1,keepdims = True)
    return exps / exp_sums

几乎准备好开始训练了。定义输入和输出的数量: 

n_inputs = X_train.shape[1]
n_outputs = len(np.unique(y_train))

现在最困难的部分来了:训练!从理论上讲，这很简单:只需将数学方程转换成Python代码即可。但在实践中，这可能相当棘手:尤其是，很容易混淆各项或指标的顺序。甚至可以最终得到看起来工作正常但实际上并没有正确计算的代码。当不确定时，应该将等式中每个项的形状写下来，并确保代码中相应的项非常匹配。它还可以帮助独立地评估每个术语并将它们打印出来。好消息是，你不必每天都这样做，因为Scikit-Learn很好地实现了这一切，但它将帮助你了解在引线下发生了什么。

所以我们需要的方程是代价函数:

和梯度方程为：

注意,如果= 0，可能不能计算。所以我们将添加一个小值𝜖记录来避免nan值。

eta = 0.01
n_iterations = 5001
m = len(X_train)
epsilon = 1e - 7

Theta = np.random.randn(n_inputs,n_outputs)

for iteration in range(n_iterations):
    logits = X_train.dot(Theta)
    Y_proba = softmax(logits)
    loss = -np.mean(np.sum(Y_train_one_hot * np.log(Y_proba + epsilon),axis = 1))

    error = Y_proba - Y_train_one_hot

    if iteration % 500 == 0:
        print(iteration,loss)

    gradients = 1 /m * X_train.T.dot(error)
    Theta = Theta - eta * gradients

看看模型参数:

Theta 

对验证集进行预测，并检查准确性分数:

logits = X_valid.dot(Theta)

Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba,axis = 1)

y_predict

accuracy_score = np.mean(y_predict == y_valid)
accuracy_score

嗯，这个模型看起来很不错。

添加一点ℓ2正规化。下面的培训代码类似于上面,但是现在有一个额外的损失ℓ2惩罚,和梯度有适当的附加项(注意,我们不规范Theta的第一个元素因为对应偏差项θ)。同时，让我们试着增加学习率eta。

eta = 0.1
n_iterations = 5001
m = len(X_train)
epsilon = 1e - 7
alpha = 0.1

Theta = np.random.randn(n_inputs,n_outputs)

for iteration in range(n_iterations):
    logits = X_train.dot(Theta)
    Y_proba = softmax(logits)
    xentropy_loss = -np.mean(np.sum(Y_train_one_hot * np.log(Y_proba + epsilon), axis=1))    
    l2_loss = 1/2 * np.sum(np.square(Theta[1:]))
    loss = xentropy_loss + alpha * l2_loss

    error = Y_proba - Y_train_one_hot

    if iteration % 500 == 0:
        print(iteration,loss)

    gradients = 1 /m * X_train.T.dot(error) + np.r_[np.zeros([1, n_outputs]), alpha * Theta[1:]]
    Theta = Theta - eta * gradients

因为ℓ2正则化,似乎大于损失,但也许这个模型将表现得更好?来看看:

logits = X_valid.dot(Theta)

Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba,axis = 1)

y_predict

accuracy_score = np.mean(y_predict == y_valid)
accuracy_score

酷,完美的准确性!

现在加上提前停止。为此，只需要在每次迭代时度量验证集的损失，并在错误开始增长时停止。

eta = 0.1
n_iterations = 5001
m = len(X_train)
epsilon = 1e - 7
alpha = 0.1
best_loss = np.infty

Theta = np.random.randn(n_inputs,n_outputs)

for iteration in range(n_iterations):
    logits = X_train.dot(Theta)
    Y_proba = softmax(logits)
    xentropy_loss = -np.mean(np.sum(Y_train_one_hot * np.log(Y_proba + epsilon), axis=1))    
    l2_loss = 1/2 * np.sum(np.square(Theta[1:]))
    loss = xentropy_loss + alpha * l2_loss

    error = Y_proba - Y_train_one_hot

    gradients = 1 /m * X_train.T.dot(error) + np.r_[np.zeros([1, n_outputs]), alpha * Theta[1:]]
    Theta = Theta - eta * gradients


    logits = X_valid.dot(Theta)
    Y_proba = softmax(logits)
    xentropy_loss = -np.mean(np.sum(Y_valid_one_hot * np.log(Y_proba + epsilon), axis=1))    
    l2_loss = 1/2 * np.sum(np.square(Theta[1:]))
    loss = xentropy_loss + alpha * l2_loss


    if iteration % 500 == 0:
        print(iteration,loss)

    if loss < best_boss:
        best_loss = loss
    else:
        print(iteration -1,best_loss)
        print(iteration,loss,"early stopping!")
        break

logits = X_valid.dot(Theta)

Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba,axis = 1)

y_predict

accuracy_score = np.mean(y_predict == y_valid)
accuracy_score

仍然完美，但是更快。

现在在整个数据集上绘制模型的预测

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()]
X_new_with_bias = np.c_[np.ones[len(X_new),1],X_new]
 
logits = X_new_with_bias.dot(Theta)
Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba,axis = 1)

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

现在在测试集中测量最终模型的准确性:

logits = X_test.dot(Theta)

Y_proba = softmax(logits)
y_predict = np.argmax(Y_proba,axis = 1)

y_predict

accuracy_score = np.mean(y_predict == y_test)
accuracy_score

完美模型原来有轻微的缺陷。这种可变性可能是由于数据集非常小:取决于如何采样训练集、验证集和测试集，可以得到非常不同的结果。尝试更改随机种子并再次运行代码几次，将看到结果会有所不同。