首先以简单的线性回归为例，讨论两种不同的训练方法来得到模型的最优解：
1.直接使用封闭方程进行根运算，得到模型在当前训练集上的最优参数， 即导致损失函数到达最小的模型参数。
2.使用迭代优化方法：梯度下降GD ，在训练集上，可以逐渐调整模型 参数，以获得最小的损失函数，最终，参数会收敛到和第一种相同的值。
梯度下降的变形体：批量梯度下降 Batch GD 、小批量梯度下降 Mini-batch GD 、随机梯度下降 Stochasitic GD 等。
然后讨论更复杂的模型：多项式回归。
最后讨论两种常用的分类模型：Logistic 回归和 Softmax 回归。
线性回归¶
预测模型
$$\widehat{y}=\theta _{0}+\theta _{1}x_{1}+\theta _{2}x_{2}...+\theta _{n}x_{n}$$
- $\widehat{y}$ 表示预测结果
- $n$ 表示特征数量
- $x_{i}$ 表示第 $i$ 个特征的值
- $\theta _{j}$ 表示第 $j$ 个参数
上面的公式可以写成更加简洁的向量形式:
$$\widehat{y}=h_{\theta }(x)=\theta ^{T}\cdot x$$
- $\theta $ 表示模型的参数向量包括偏置项 $\theta _{0}$ 和特征权重 $\theta _{1}$ 到 $\theta _{n}$
- $\theta ^{T}$ 表示向量 $\theta $ 的转置
- $x$ 为每个样本中特征值的向量形式，包括 $x_{1}$ 到 $x_{n}$ ，并且 $x_{0}$ 恒为1
- $\theta ^{T}\cdot x$ 表示 $\theta ^{T}$ 和 $x$ 的点积
- $h_{\theta }$ 表示参数为 $\theta $ 的假设函数
损失函数
$RMSE$ 和 $MSE$ 会得到相同的 $\theta $，$MSE$ :
$$MSE(X,h_{\theta })=\frac{1}{m}\sum_{i=1}^{m}(\theta ^{T}\cdot x^{(i)}-y^{(i)})^{2}$$
正态方程
$$\widehat{\theta }=(X^{T}\cdot X)^{-1}\cdot X^{T}\cdot y$$
- $\widehat{\theta }$ 最小化损失 $\theta $ 的值
- $y$ 是一个向量，其包含 $y^{(1)}$ 到 $y^{(m)}$ 的值
import numpy as np
X = 2*np.random.rand(100, 1)
y = 4+3*X+np.random.randn(100, 1)
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
# 用正态方程来计算，我们将使用 np.linalg 中的 inv 函数来计算矩阵的逆，
# 使用 dot 来计算矩阵的乘法
X_b = np.c_[np.ones((100, 1)), X]
# T 转置， dot 点积， inv 计算矩阵的逆
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
# 事实上 theta_0 = 4, theta_1 = 3
# 现在的 theta_0 = 3.9... theta_1 = 3.0...
# 由于噪声的存在，参数不可能达到原始函数的值
# 现在我们能够使用 theta_hat 来进行预测
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()
# 使用 sklearn 的 LinearRegression 也可以办到正态方程求解
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)
print(lin_reg.intercept_, lin_reg.coef_)
lin_reg.predict(X_new)
计算复杂度¶
正态方程要求我们计算 $(X^{T}\cdot X)^{-1}$ ，它是 $n*n$ 的矩阵，运算复杂度在 0(2.4) 到 O(3) 之间，我们将特征个数翻倍的话，其计算时间大概会变成原来的 5.3 到 8倍。
下面介绍的方法适合在特征个数非常多，训练实例非常多，内存无法满足要求的时候使用。
梯度下降¶
我们迷失在大山中，最好的下山办法是往坡度最陡的地方下山，其实就是梯度下降所做的。计算误差函数关于参数向量theata 的局部梯度。初始，我们要选定一个随机的 theta ，然后去逐渐改进它，每一次变化一小步直到算法收敛到一个最小值。而学习率就是步长。
陷阱：我们可能永远无法到达全局最小值或者只能到达局部最小值。幸运的是线性回归模型的 MSE 是一个凸函数，仅有一个全局最小值，梯度下降可以无限接近全局最小值。
在梯度下降的时候，应该确保所有的特征有着相近的尺度范围，否则将很难收敛。
批量梯度下降¶
在梯度下降的过程中，我们需要计算每一个 theta 下损失函数的梯度。 换而言之，我需要计算的是 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)}$$
我们也可以使用公式，一次计算全部梯度，这被称作梯度向量：
$$\bigtriangledown _{\theta }MSE(\theta )=\bigl(\begin{smallmatrix}\frac{\partial }{\partial \theta _{0}}MSE(\theta ) \\ \frac{\partial }{\partial \theta _{1}}MSE(\theta ) \\ ... \\ \frac{\partial }{\partial \theta _{n}}MSE(\theta ) \end{smallmatrix}\bigr)=\frac{2}{m}X^{T}\cdot (X\cdot \theta -y)$$
每一步计算都包含了整个训练集 X ，这也是为什么这个算法被称为批量梯度下降，因此，在大数据集上，其会变得相当的慢，然而梯度下降的运算规模和特征的数量成正比。训练一个数千数量特征的线性回归模型使用梯度下降要比使用正态方程快得多。
一旦求得了方向是上山的梯度向量，就可以向着相反的方向去下山，同时梯度向量和学习率的积决定了下山每一步的步长：
$$\theta ^{nextstep}=\theta -\eta \bigtriangledown _{\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
theta
# 与上面的内容完全一致，仅仅是将每次迭代画出来
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):
# 只计算前10轮
if iteration < 10:
y_predict = X_new_b.dot(theta)
# iteration = 0 是初始状态
style = 'b-' if iteration > 0 else 'r--'
# 绘画前10轮的梯度下降结果
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)
if iteration >= 10:
break
plt.xlabel('$x_1$', fontsize = 18)
plt.ylabel('$y$', rotation=0, 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)
# width height
plt.figure(figsize = (10, 4))
# 学习率 = 0.02
plt.subplot(131); plot_gradient_descent(theta, eta=0.02)
# 学习率 = 0.1
plt.subplot(132); plot_gradient_descent(theta, eta=0.1,
theta_path=theta_path_bgd)
# 学习率 = 0.5
plt.subplot(133); plot_gradient_descent(theta, eta=0.5)
随机梯度下降¶
# 批量梯度下降的问题是计算每一步梯度时都需要使用整个训练集
# 这导致在规模较大的数据集上，将会变得非常慢，与其相反，随机梯度下降，
# 每一步在梯度计算上只随机选取训练集中的一个样本，由于每一次操作都使用
# 非常少的数据，这样使得算法变得非常快，每一次迭代，只需要在内存中有
# 一个实例，这使得随机梯度算法可以在大规模训练集上使用。
# 另一个方面，因为它的随机性，与批量梯度下降相比，其呈现出更多的不规律性，
# 它到达最小值不是平缓的下降，损失函数会忽高忽低，只是在大体上呈现下降，
# 它会非常靠近最小值，但是它不会停下，它会在这个值附近摆动。结果非最优。
# 当损失函数很不规则的时候，随机梯度下降算法能够跳过局部最小值，因此，
# 随机梯度下降在寻找最小值上比批量梯度下降表现要好。
# 为了解决让随机梯度下降找到最小值点，我们可以尝试减小学习率，开始，每一步
# 都较大，然后越来娱小，从而达到全局最小值，这个过程是模拟退火，
# 因为它类似于熔融金属慢慢冷却的炼金学退火过程，决定学习率的函数是
# learning schedule 。
# 但是如果学习速度降低得过快，那么可能会陷入局部最小，甚至在半路就停止；
# 如果下降过慢，可能在最小值的附近长时间摆动，同时可能出现次优解。
theta_path_sgd = []
np.random.seed(42)
n_epochs = 50
t0, t1 = 5, 50 # learning_schedule 的超参数
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 < 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 = 2/1*xi.T.dot(xi.dot(theta)-yi)
# 更新学习率
eta = learning_schedule(epoch*m+i)
# 改变参数 theta
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])
plt.show()
# 现在每一个实例的选择都是随机的，有的实例可能在每一代中都被选到
# 其他的实例也可能一直不被选到，如果想保证每一代迭代过程，算法可以遍历
# 所有实例，一种方法是将训练集打乱重排，然后选择一个实例，之后再打乱，
# 但是这样收敛会非常慢。
# 通过 sklearn 完成线性回归的随机梯度下降，我们可以使用 SGDRegressor ，
# 它默认优化的是 MSE ，下面的代码迭代了50代，eta = 0.1
# 使用默认的 learning schedule （与之前的不一样）
# penatly = None 没有添加任何正则项
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=50, penalty=None, eta0=0.1)
sgd_reg.fit(X, y.ravel())
# intercept 顾名思义在 y 轴上的拦截
# coef 就是分配的权重
# 合起来我们就可以写出线性方程 coef(x) + intercept
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):
# 随机打乱 100 个数
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)
theta
# 对比三种梯度下降，第一维是 intercept ，第二维是 coef
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(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$ ", fontsize=20, rotation=0)
plt.axis([2.5, 4.5, 2.3, 3.9])
plt.show()
多项式回归¶
# 就算我们的数据实际上币简单的直线更复杂，我们依然可以使用线性模型来
# 拟合非线性数据，一个简单的方法是对每个特征进行加权后作为新的特征，
# 然后训练一个线性模型在这个扩展的特征集，这个方法称为多项式回归。
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])
plt.show()
# 从图中很明显可以发现，直线不能恰当拟合这些数据，于是我们可以使用
# sklearn 的 PolynomiaFeatures 类来进行训练集的转换，让训练
# 集中每个特征的平方作为新特征
from sklearn.preprocessing import PolynomialFeatures
poly_features=PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
X[0]
X_poly[0]
# 现在的 X_poly 包含原始特征 X 和这个特征的平方 X2 ，可以在这个扩展
# 训练集上使用 LinearRegression 模型进行拟合
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
lin_reg.intercept_, lin_reg.coef_
# -3 到 3 的 100 个等间隔点
X_new = np.linspace(-3, 3, 100).reshape(100, 1)
# 获得特征的平方
X_new_poly = poly_features.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$', fontsize = 18, rotation = 0)
plt.legend(loc='upper left', fontsize = 14)
plt.axis([-3, 3, 0, 10])
plt.show()
学习曲线¶
# 如果我们使用一个高阶的多项式回归，可能发现它的拟合程度比普通的
# 线性回归要好。
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_feature', polybig_features),
('std_scalar', 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])
plt.show()
# 当然，高阶模型严重过拟合了，线性模型有欠拟合，在这个训练集上，二次模型
# 有较好的泛化能力。如何知道一个模型是过拟合还是欠拟合？
# 1.如果一个模型在训练集上表现良好，通过交叉验证指标却得出其泛化能力差，
# 说明过拟合了。
# 2.通过观察学习曲线，画出模型在训练集上的表现，同时画出以训练集规模为
# 自变量的训练集函数，为了得到图像，需要在训练集的不同规模子集上进行多次
# 训练。
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)
plt.axis([0, 80, 0, 3])
plt.show()
lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)
# 这幅图值得我们深究，首先，我们观察训练集的表现：
# 1.当训练集只有一两个样本
# 的时候，模型能够非常好的拟合他们，这也是为什么曲线是从零开始的原因。
# 但是当加入了一些新的样本的时候，训练集上的拟合程度变得难以接受，
# 出现这种情况的原因有两个，一是因为数据中含有噪音，另一个原因是数据
# 是非线性的。随着数据规模的增大，误差会一直增大，直到达到高原地带并
# 趋于稳定，在之后，继续加入新的样本，模型的 RMSE 不会变的更好或者
# 更差。
# 2.我们继续来看模型在验证集上的表现，当以非常少的样本去训练时，
# 模型不能恰当的泛化，也就是为什么验证误差一开始是非常大的，当训练样本
# 变多的时候，模型学习的东西变多，验证误差开始缓慢的下降，但是一条直线
# 不可能很好的拟合这些数据，因此最后误差会到达一个高原地带并趋于稳定。
# 这就是一个非常典型的欠拟合模型，两条曲线都到达高原地带并趋于稳定，并且
# 最后两条曲线非常接近，同时误差值非常大。
# 欠拟合情况下，添加更多的样本是没用的，你需要使用一个更复杂的模型或者找到
# 更好的特征。
# 现在我们在一个相同数据上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)
这幅图和之前的有一点点像，但是其有两个非常重要的不同点。
- 在训练集上，误差要比线性回归模型低的多。
- 图中两条曲线之间有间隔，这以为模型在训练集上的表现比验证集好。
改善过拟合的一种方法是提供更多的训练数据，直到训练误差和验证误差相等。 在统计和机器学习领域有个重要的理论，一个模型的泛化误差由是三个不同的 误差和决定：
- 1.偏差，泛化误差的这部分误差是由于错误的假设决定的。
- 2.方差，模型对训练数据的微小变化较为敏感，高阶多项式会导致模型过拟合。
- 3.不可约误差，数据本身的噪声决定的。
线性模型的正则化¶
# 降低模型过拟合的好方法是正则化这个模型，模型自由度越小，越难以拟合数据。
# 正则化一个多项式模型，一个简单的方法是减少多项式的阶数。
# 对于一个线性模型，正则化的电信实现是约束模型中参数的权重，下面有三种不同
# 的约束权重的方法： Ridge 回归、 Lasso 回归、 Elastic Net
岭 (Ridge) 回归¶
岭回归也称为 Tikhonov 正则化，是线性回归的正则化版：在损失函数上直接加上一个正则项 $a\sum_{i=1}^{n}\theta ^{2}$ 。这使得讯息算法不仅能够拟合数据，而且能够使模型的参数权重尽量的小。
注意到这个正则化只有在训练过程中才会被加到损失函数，当得到完成训练的模型后，我们还应该使用没有正则化的测量方法去评价模型的表现。
ps: 一般情况下，训练过程中使用的损失函数和测试过程中使用的评价函数是不一样的。除了正则化，还有一个不同：训练时的损失函数应该在优化过程中易于求导，而在测试过程中，评价函数应该更接近最后的客观表现。一个好的例子：在分类训练中我们使用对数损失作为损失函数，但是我们却使用准确率/召回率来作为它的评价函数。
超参数 $\alpha $ 决定了你想正则化这个模型的强度。如果 $\alpha =0$ 那么此时岭回归便变为了线性回归。如果 $\alpha $ 非常的大，所有的权重最后都接近于零，最后结果将是一条穿过数据平均值的水平直线，岭回归的损失函数：
$$J(\theta )=MSE(\theta )+\alpha \frac{1}{2}\sum_{i=1}^{n}\theta _{i}^{2}$$
值得注意的是偏差 theta_0 是没有被正则化的，累加从 i=1 开始。同时，在使用岭回归之前，对数据进行放缩是非常重要的，算法对于输入特征的数值尺度非常敏感，大多数的正则化模型都是这样的。
import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import PolynomialFeatures
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
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 polynomial:
model = Pipeline([
('poly_features', PolynomialFeatures(degree = 10, include_bias = False)),
('std_scaler', StandardScaler()),
('regul_reg', model),
])
model.fit(X, y)
y_new_regual = model.predict(X_new)
lw = 2 if alpha > 0 else 1
plt.plot(X_new, y_new_regual, style, linewidth=1, 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)
plot_model(Ridge, polynomial=False, alphas=(0, 10, 100), random_state=42)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(122)
plot_model(Ridge, polynomial=True, alphas=(0, 10**-5, 1), random_state=42)
plt.show()
# 上面展示了在相同线性数据上使用不同 alpha 值得岭回归模型最后的表现，左
# 图使用简单的岭回归模型，最后得到了线性的预测；右图使用10阶扩展，然后
# 使用 StandarScaler 进行缩放，最后将岭模型应用到处理过后的特征上。
# 注意将 alpha 加大，导致预测曲线变得扁平，即减少了极端值，多了一般值
# 这样减少了模型的方差，却增加了模型的偏差。
# 对线性回归来说，岭模型我们可以使用封闭方程来计算，也可以使用梯度下降
# 去处理。它们的缺点和优点是一样的。
A: 一个除了左上角有一个 0 的 n 阶单位矩阵，这个 0 代表偏差项。
$$\widehat{\theta }=(X^{T}\cdot X+\alpha A)^{-1}\cdot X^{T}\cdot y$$
# 下面使用封闭方程来求解，使用 cholesky 法进行矩阵分解对公式变形
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]])
# 使用随机梯度法来求解
# penalty 是正则化的惩罚类型，指定 l2 表明你要在损失函数上添加一项：
# 权重向量 l2 范数平方的一般，这就是简单的岭回归
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter = 5, penalty = 'l2', random_state = 42)
sgd_reg.fit(X, y.ravel())
sgd_reg.predict([[1.5]])
# 也可以使用 sovler = sag 来进行随机平均梯度下降
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 Operation Regression ，是另一种正则化版的线性回归，就像岭回归那样，它在损失函数上添加一个正则化项，但是它使用权重向量的 l1 范数而不是权重向量的 l2 范数平方的一半。损失函数：
$$J(\theta )=MSE(\theta)+\alpha \sum_{i=1}^{n}\left | \theta _{i} \right |$$
from sklearn.linear_model import Lasso
plt.figure(figsize=(8, 4))
plt.subplot(121)
plot_model(Lasso, polynomial=False, alphas=(0, 0,1, 1), random_state = 42)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(122)
plot_model(Lasso, polynomial=True, alphas=(0, 10**-7, 1), tol=1, random_state=42)
plt.show()
Lasso 回归的一个重要特征是它倾向于完全消除最不重要的特征的权重，即将他们设置为零，例如 alpha = 10-7 ，看起来像是一条二次曲线，而且几乎是线性的，这是因为所有的高阶多想特征都被设置为0，换句话说， Lasso 回归自动的进行特征选择同时输出一个稀疏模型。
Ridge 和 Lasso 对比： 在 Lasso 损失函数中，批量梯度下降的路径趋向于在低谷有一个反弹，这是因为在 theta2 = 0 时斜率有一个突变，为了最后真正收敛到全局最小值，你需要逐渐降低学习率。
Lasso 损失函数在 theta_i=0 (i=1~n) 处都无法进行微分运算，但是梯度下降如果你使用子梯度向量 g 后，它可以在任何 theta_i=0 的情况下进行计算，下面是 Lasso 损失函数在进行梯度下降的子梯度向量公式：
$$g(\theta ,J)=\bigtriangledown _{\theta }MSE(\theta )+\alpha \begin{pmatrix}sign(\theta _{1}) \\ sign(\theta _{2}) \\ ... \\ sign(\theta _{n}) \end{pmatrix}\quad where\quad sign(\theta {i})=\left\{\begin{matrix}-1,\theta_{i}<0 \\ 0,\theta_{i}=0 \\ +1,\theta_{i}>0 \end{matrix}\right.$$
from sklearn.linear_model import Lasso
lasso_reg = Lasso(alpha=0.1)
lasso_reg.fit(X, y)
lasso_reg.predict([[1.5]])
# 使用随机梯度来求解
# penalty = l1
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=5, penalty='l1', random_state=42)
sgd_reg.fit(X, y.ravel())
sgd_reg.predict([[1.5]])
弹性网络 (ElasticNet)¶
# 弹性网络介于 Ridge 回归和 Lasso 回归之间，它的正则项是 Ridge 回归和 Lasso 回归正则项的简单混合，
# 同时可以控制它们的混合率 r ，当 r = 0 ，弹性网络就是 Ridge 回归，当 r = 1 ，就是 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}$$
我们应该如何选择线性回归？岭回归，Lasso 回归，弹性网络？一般来说有一点正则项的表现更好，因此避免选择简单的线性回归。
- 岭回归是一个很好的首选项
- 如果特征仅有少数是真正有用的，应该选择 Lasso 和弹性网络
- 一般来说，弹性网络的表现比 Lasso 要好，因为当特征数量比样本的数量大的时候，或者特征之间有很强的相关性的时候， Lasso 可能会表现的不规律
# l1_ratio 指的就是混合率 r
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]])
早期停止法 (Early Stopping)¶
# 对于迭代学习算法，有一种非常特殊的正则化方法，就像梯度下降在验证错误达到最小值时立即停止训练那样，
# 我们称之为早期停止法。
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
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,
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)
plt.show()
随机梯度和小批量梯度下降都不是平滑曲线，我们可能很难知道它时候达到了最小值，一种解决方案是，只有在验证误差高于最小值一段时间后才停止，之后将模型参数回滚到验证误差最小值。
# 这是一个早期停止法的基础应用
from sklearn.base import clone
sgd_reg = SGDRegressor(max_iter=1, 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)
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
逻辑回归¶
一些回归算法也可以用于分类，反之亦然。 Logistic 回归通常用于估计一个实例属于某个特定类别的概率，利于这是垃圾邮件的概率是多少？
概率估计¶
就像线性回归模型一样， Logistic 回归模型计算输入特征的加权和 (加上偏差项) ，但它不像线性回归模型那样直接输出结果，而是把结果输入 logistic 函数进行二次加工后进行输出，逻辑回归模型的概率估计：
$$\widehat{p}=h_{\theta }(x)=\sigma (\theta ^{T}\cdot x)$$
Logistic 函数用 sigma() 表示，其是一个 sigmoid 函数，它的输出介于 0~1。
训练和损失函数¶
现在我们知道了 Logistic 回归模型如何估计概率并进行预测，它是如何训练的呢？训练的目的是设置参数向量 theta ，使得正例 (y=1) 的概率增加，负例 (y=0) 的概率减少，其通过在单个训练实例 x 的损失函数来实现。
$$c(\theta )=\left\{\begin{matrix}-log(\widehat{p}),y=1 \\ -log(1-\widehat{p}),y=0 \end{matrix}\right.$$
当 t 接近于 0，-log(t) 将变得非常大，所以如果模型估计一个正例概率接近于0，那么损失函数将会很大； 同时，如果模型估计一个负例的概率接近1，那么损失函数也会很大。反之同理。所以这个损失函数是合理的。
整个训练集的损失函数只是所有训练实例的平均值，可以用一个表达式来统一表示，称之为对数损失。
$$J(\theta )=-\frac{1}{m}\sum_{i=1}^{m}\left [ y^{(i)}log(\widehat{p}^{(i)})+(1-y^{(i)})log(1-\widehat{p}^{(i)}) \right ]$$
现在有一个坏消息和一个好消息：
- 这个损失函数对于最小化损失函数的 theta 没有公式解
- 这是凸函数，一定有全局最小值。
$$\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)}$$
对损失函数求偏导数，可以获得所有偏导数的梯度向量，然后在梯度向量上使用批量梯度下降算法。
决策边界¶
from sklearn import datasets
iris = datasets.load_iris()
# 花瓣的宽度
X = iris['data'][:, 3:]
# 1 if iris-Virginica else 0
y = (iris['target']==2).astype(np.int)
# 逻辑回归
from sklearn.linear_model import LogisticRegression
log_reg = LogisticRegression()
log_reg.fit(X, y)
# 我们看看模型估计的花瓣宽度从 0~3 cm 的概率估计：
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', linewidth=2)
plt.plot(X_new, y_proba[:, 0], "b--", linewidth=2, label="Not Iris-Virginica")
plt.axis([0, 3, -0.02, 1.02])
plt.xlabel("Petal width (cm)", fontsize=14)
plt.ylabel("Probability", fontsize=14)
plt.legend(loc="center left", fontsize=14)
Virginica 花的花瓣宽度在1.4到2.5 cm 之间，而其他种类的话通常具有较小的花瓣宽度，范围从0.1到1.8 cm 。他们有重叠。2 cm 以上，分类器肯定这是 Virginica ，而在 1cm 一下，一定不是 Virginica 。在这两个极端之间，分类器是不确定的。但是，如果你使用它进行预测，它将返回一个最可能的结果，因此在 1.6cm 附近存在一个决策边界，这时两类出现的概率等于 50% 。
# 预测结果
log_reg.predict([[1.7], [1.5]])
# 下面表示相同的数据，但是这次使用两个特征来进行判断：花瓣的宽度和长度，。
# 一旦训练完毕， Logisitic 回归分类器就可以根据这两个特征来估计一朵花是 Virginica 的可能性。
from sklearn.linear_model import LogisticRegression
# petal length, petal width
X = iris['data'][:, (2, 3)]
y = (iris['target'] == 2).astype(np.int)
# C 是 alpha 的逆， alpha 是正则化强度的超参， C 越大，正则化强度越低，逻辑回归默认使用了 l2 惩罚
log_reg = LogisticRegression(C=10**10, random_state=42)
log_reg.fit(X, y)
# 画网格， x0 是 x 轴上的分割点， x1 是 全部一个数字的数字，比如 0.8 、 2.7
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))
# 挑选出是否是 Virginica 花
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])
plt.show()
y_proba.shape
# 这是一条线性边界，每条平行线代表一个分类标准下两个不同类的概率，从 15% 到 90% 。
Softmax 回归¶
逻辑回归模型可以直接推广到支持多类别分类，不必组合和训练多个二分类器，其称为 Softmax 回归或多类别 Logistic 回归。
当给定一个实例 x 时， Softmax 回归模型首先计算 k 类的分数 s_{k}(x) ，然后将分数应用在 Softmax 函数上，估计每类的概率。
$$分数 -> 概率 -> argmax(概率) -> 类别$$
k 类的 sofmax 得分：
$$S_{k}(x)=\theta ^{T}\cdot x$$
每个类都有自己的参数向量 theta，所有这些向量通常作为行放在参数矩阵中。 一旦我们获得每一类的得分，就可以通过 Softmax 函数来估计样本属于第 k 类的概率 pk ：
$$\widehat{p_{k}}=\sigma (s(x))k=\frac{exp(s_{k}(x))}{\sum_{j=1}^{K}exp(s_{j}(x))}$$
- K 表示有多少类
- s(x) 表示包含样本 x 每一类得分的向量
- sigma(s(x))_k 表示给定每一类分数后，实例 x 属于第 k 类的概率
和 logistic 回归一样， softmax 回归也将估计概率最高的那类作为预测结果。
$$\widehat{y}=argmax\sigma (s(x))_{k}=argmaxs_{k}(x)=argmax(\theta _{k}^{T}\cdot x)$$
现在知道如何估计概率并进行预测，接下来介绍如何训练。目标是建立一个模型，在目标类别上有较高的概率，最小化公式可以达到这个目的，其表示当前模型的损失函数，称为交叉熵，当模型对目标类得出一个概率，会惩罚这个模型。交叉熵：
$$J(\Theta )=-\frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{K}y_{k}^{(i)}log(\widehat{p}_{k}^{(i)})$$
如果对于第 i 个实例的目标类是 k ，那么 y_k^i=1 ，反之 y_k^i=0 。
可以发现，当只有两个类 k=2 时，此损失函数等同于 Logistic 回归的损失函数，
k类交叉熵的梯度向量：
$$\bigtriangledown _{\theta _{k}}J(\Theta )=\frac{1}{m}\sum_{i=1}^{m}(\widehat{p}_{k}^{(i)}-y_{k}^{(i)})x^{(i)}$$
现在可以计算每一类的梯度向量，然后使用梯度下降找到使得损失函数达到最小值的参数矩阵 Theta
现在使用 Softmax 回归对三种鸢尾花进行分类，在使用 LogisticRegression 训练时，默认是一对多的模型，我们可以设置 multi_class 参数为 multinomial 来改变它的 softmax 回归，还需要指定一个求解器，比如 lbfgs 求解器，它默认使用 l2 正则化。
# 花瓣的长和宽
X = iris['data'][:, (2, 3)]
y = iris['target']
softmax_reg = LogisticRegression(multi_class='multinomial', solver='lbfgs', C=10)
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])
plt.show()
softmax_reg.predict([[5, 2]])
softmax_reg.predict_proba([[5, 2]])
思考问题
- ridge 和 lasso 的对比
- 逻辑回归不同分类标准下不同的概率