文章目录
-
第4章 训练模型
1 线性回归
1.1 标准方程
1.2 计算复杂度
2 梯度下降
2.1 批量梯度下降
2.2 随机梯度下降
2.3 小批量梯度下降
3 多项式回归
4 学习曲线
5 正则线性模型
5.1 岭回归
5.2 套索回归
5.3 弹性网络
如何选择线性回归、岭回归、lasso回归和弹性网络呢？5.4 早期停止法
5.5 逻辑回归
5.5.1 概率估算
5.5.2 训练和成本函数
5.5.3 决策边界
5.6 Softmax回归
多分类如何选择逻辑回归和Softmax回归？6 练习
本章主要是讲模型原理，即模型是如何工作的？了解这些内容有助于快速选择合适的模型，正确的训练算法(GD等)，以及一套适当的超参数。
1. 线性回归
预测公式（输入x，预测y）
y ^ = θ 0 + θ 1 x 1 + θ 2 x 2 + . . . + θ n x n \hat y = \theta_0 + \theta_1x_1 + \theta_2 x_2 + ... + \theta_n x_n y^=θ0+θ1x1+θ2x2+...+θnxn
y ^ \hat y y^ 是预测值
n n n 是特征的数量
x i x_i xi 是第i个特征值
θ j \theta_j θj 是第j个模型参数（包括偏置项 θ 0 \theta_0 θ0以及特征权重 θ 1 \theta_1 θ1, θ 2 \theta_2 θ2 … θ n \theta_n θn）
预测公式的向量化表示：
y ^ = h θ ( X ) = θ T ⋅ X \hat y = h_\theta(X) = \theta^T\cdot X y^=hθ(X)=θT⋅X
θ T \theta^T θT为转置后的行向量，不再是列向量。X为特征向量。(二维平面中可理解为一元一次的直线方程)。
那么如何使用这个线性回归模型呢？无非就是迭代训练直到找到最适应训练集的那个 θ \theta θ。怎么才算是最适应？此时需要一个衡量的标准，比如MSE，MAE等。
成本函数(也叫损失函数) MSE:
M S E ( X , h θ ) = 1 m ∑ i = 1 m ( θ T ⋅ X ( i ) − y ( i ) ) 2 MSE(X,h_\theta) = \frac{1}{m}\sum_{i=1}^m(\theta^T\cdot X^{(i)} - y^{(i)} )^2 MSE(X,hθ)=m1i=1∑m(θT⋅X(i)−y(i))2
所以，我们要求出损失函数最小时的 θ 值 \theta值 θ值。
1.1 标准方程
通过上面的分析，要求损失最小时的 θ \theta θ，有一个闭式解公式可以直接得出：
θ ^ = ( X T ⋅ X ) − 1 ⋅ X T ⋅ y \hat \theta = (X^T \cdot X)^{-1}\cdot X^T \cdot y θ^=(XT⋅X)−1⋅XT⋅y
验证公式：
import numpy as np
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)#高斯噪声
# 每个实例加入1，原来的X为100行1列，变为100行2列，第一列为1
# np.ones((100, 1)) 生成100行1列的矩阵，值全为1
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的值
theta_best
原函数的参数为： θ 0 = 4 , θ 1 = 3 \theta_0=4,\theta_1=3 θ0=4,θ1=3 ，我们训练出来的参数还比较接近 θ 0 = 4.305 , θ 1 = 2.943 \theta_0=4.305,\theta_1=2.943 θ0=4.305,θ1=2.943。如果把噪声去掉，预测出来的就是4和3。
我们用得出的 θ \theta θ进行预测:
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的实现(注意:sk将截距
θ
0
\theta_0
θ0和系数
θ
1
\theta_1
θ1分开了)。
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)
1.2 计算复杂度
求逆矩阵
X
T
⋅
X
X^T \cdot X
XT⋅X，是n 乘 n矩阵(n是特征数量)。这种矩阵求逆的计算复杂度通常为
O
(
n
2.4
)
O(n^{2.4})
O(n2.4) 到
O
(
n
3
)
O(n^3)
O(n3)。如果特征数量n翻倍，比如由1变为2，那么计算时间大约为
2
2.4
=
5.3
2^{2.4} = 5.3
22.4=5.3 倍到
2
3
=
8
2^3=8
23=8倍之间。
2. 梯度下降
找到最优解 θ \theta θ 的一种优化算法，通过不断的迭代不断调整参数，使损失函数最小化。
具体做法：通过测量参数向量 θ \theta θ的相关误差函数的局部梯度，并不断沿着降低梯度的方向调整，直到梯度降为0，到到最小值。
梯度下降中一个重要参数:每一步的步长，取决于超参数学习率。学习率太低，需要大量迭代，学习率太高，容易扯着裆，直接跳过最低点，可能无法找到最小值。
贴张大神吴恩达的帮助图理解下：
注意：应用梯度下降时，需保证所有特征值的大小比例都差不多(比如使用StandardScaler类)，否则收敛时间会长很多。
训练模型就是搜寻使成本函数(在训练集上)最小化的参数组合。
2.1 批量梯度下降
要实现梯度下降，需要计算关于参数
θ
j
\theta_j
θj的成本函数的梯度，即偏导数。(根据复合函数的求导法则得到)
成本函数的偏导数:
∂ ∂ θ j M S E ( θ ) = 2 m ∑ i = 1 m ( θ T ⋅ x ( i ) − y ( i ) ) x j ( i ) \frac{\partial}{\partial\theta_j}MSE(\theta) = \frac{2}{m}\sum_{i=1}^m(\theta^T\cdot x^{(i)} - y^{(i)} )x_j^{(i)} ∂θj∂MSE(θ)=m2i=1∑m(θT⋅x(i)−y(i))xj(i)
成本函数的梯度向量:
∇ θ M S E ( θ ) = 2 m X T ⋅ ( X ⋅ θ − y ) \nabla_\theta MSE(\theta) = \frac{2}{m}X^T \cdot (X\cdot \theta - y) ∇θMSE(θ)=m2XT⋅(X⋅θ−y)
计算梯度下降的每一步时，都是基于完整的训练集X。所以叫批量梯度下降。
得到梯度向量之后，朝反方向下坡。也就是从 θ \theta θ中减去 ∇ θ M S E ( θ ) \nabla_\theta MSE(\theta) ∇θMSE(θ)。梯度下降步长：
θ ( n e x t s t e p ) = θ − η ∇ θ M S E ( θ ) \theta^{(next step)} = \theta - \eta \nabla_\theta MSE(\theta) θ(nextstep)=θ−η∇θMSE(θ)
快速实现
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
if iteration < 3:
print(theta)
theta
3种不同学习率前10步的梯度下降寻找最优解过程：
theta_path_bgd = []
#迭代的同时画出梯度下降寻找最优theta的过程。
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()
要找到合适的学习率，可以使用网格搜索。
2.2 随机梯度下降
每一步在训练集中随机选择一个实例，并且仅基于该单个实例来计算梯度。但是，由于算法的随机性质，SGD比BGD要不规则得多。
优点：可以逃离局部最优，缺点：永远定位不出最小值。解决方法：逐步降低学习率(每次迭代的学习率都不同)。这个过程叫做模拟退火。确定每个迭代学习率的函数叫学习计划。
利用简单的学习计划实现随机梯度下降SGD：
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 < 20: # not shown in the book
y_predict = X_new_b.dot(theta) # not shown
style = "b-" if i > 0 else "r--" # not shown
plt.plot(X_new, y_predict, style) # not shown
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) #求的是单个样本实例的步长，区别于BGD
eta = learning_schedule(epoch * m + i)
theta = theta - eta * gradients
theta_path_sgd.append(theta) # not shown
plt.plot(X, y, "b.") # not shown
plt.xlabel("$x_1$", fontsize=18) # not shown
plt.ylabel("$y$", rotation=0, fontsize=18) # not shown
plt.axis([0, 2, 0, 15]) # not shown
save_fig("sgd_plot") # not shown
plt.show()
sklearn SGD实现：
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=50, tol=-np.infty, penalty=None, eta0=0.1, random_state=42)
sgd_reg.fit(X, y.ravel())
sgd_reg.intercept_, sgd_reg.coef_
2.3 小批量梯度下降
每一步的梯度计算，基于一小部分随机的实例集。
theta_path_mgd = []
n_iterations = 50
minibatch_size = 20
np.random.seed(42)
theta = np.random.randn(2,1) # random initialization
t0, t1 = 200, 1000
def learning_schedule(t):
return t0 / (t + t1)
t = 0
for epoch in range(n_iterations):
shuffled_indices = np.random.permutation(m)
X_b_shuffled = X_b[shuffled_indices]
y_shuffled = y[shuffled_indices]
for i in range(0, m, minibatch_size):
t += 1
xi = X_b_shuffled[i:i+minibatch_size]
yi = y_shuffled[i:i+minibatch_size]
gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)
eta = learning_schedule(t)
theta = theta - eta * gradients
theta_path_mgd.append(theta)
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=2, label="Mini-batch")
plt.plot(theta_path_bgd[:, 0], theta_path_bgd[:, 1], "b-o", linewidth=3, 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])
save_fig("gradient_descent_paths_plot")
plt.show()
3. 多项式回归
如果数据比简单的直线更为复杂怎么办？其实也可以用线性模型来拟合非线性数据。
方法：将每个特征的幂次方添加为一个新特征，然后在扩展的特征集上训练线性模型。这种方法被称为多项式回归。
制造些非线性数据
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.axis([-3, 3, 0, 10])
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.show()
save_fig("quadratic_data_plot")
显然，直线不可能拟合这个数据。使用sk将每个特征的平方扩展为新的特征加入训练集。然后用线性回归去拟合非线性的训练数据。
from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False) #2次，不包含偏置项
X_poly = poly_features.fit_transform(X) #
X[:2]
X_poly[:2]
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) #linspace：返回在指定范围内的均匀间隔的数字（组成的数组），也即返回一个等差数列
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$", rotation=0, fontsize=18)
plt.legend(loc="upper left", fontsize=14)
plt.axis([-3, 3, 0, 10])
plt.show()
多项式回归能发现特征与特征的关系(纯线性回归做不到)。原因：PolynomialFeatures会在给定的阶数下，添加所有特征组合。举例：特征a和b，阶数degree=3，除了a,b，扩展的新特征有 a 2 、 a 3 、 b 2 、 b 3 、 a b 、 a 2 b 以 及 a b 2 a^2、a^3、b^2、b^3、ab、a^2b以及ab^2 a2、a3、b2、b3、ab、a2b以及ab2。
所以，要特别小心:特征组合数量的爆炸。
4. 学习曲线
学习曲线描述的是预测函数和实际函数的对比，及特征X和label的关系。个人认为是预测准确度的体现，还有一种描述性能的学习曲线。
分别画出300,2,1阶的效果。
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", polybig_features),
("std_scaler", std_scaler),
("lin_reg", lin_reg),
])
polynomial_regression.fit(X, y) #训练模型 (此时X已经是扩展过的训练集，并进行了标准化处理)
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()
另一种学习曲线：模型在训练集和验证集上，关于训练集大小的性能函数。
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) # not shown in the book
plt.xlabel("Training set size", fontsize=14) # not shown
plt.ylabel("RMSE", fontsize=14)
lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)
plt.axis([0, 80, 0, 3])
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]) # not shown
save_fig("learning_curves_plot") # not shown
plt.show()
分析：上图可以看到，训练集上表现比验证集上要好很多，说明过拟合了。
******************小结（重要）******************：
① 学习曲线可以用于选择模型，可以观察到是否过拟合等等。
② 偏差/方差的权衡，统计和机器学习领域重要理论结果，：
模型的泛化误差 = 偏差 + 方差 + 不可避免的误差。
偏差：假设函数就没选对，比如假设函数是线性的，实际数据确是2次的，高偏差可能造成欠拟合。
方差：模型对训练数据的微小变化比较敏感，高度自由的模型很可能有高方差，所以很容易对训练数据过拟合。
不可避免的的误差： 数据本身所致，比如异常值，缺失值等，需要数据处理。
总结，增加模型的复杂度会显著提升模型的方差，减小偏差。反之，降低模型的复杂度则会降低模型的方差，提升偏差。所以，一个好的模型要平衡好偏差和方差，才能使得泛化能力更好。
再贴张大佬图帮助理解：
5. 正则线性模型
减少过拟合的一个好办法：对模型正则化 ( 约束他，不要过度自由)。
5.1 岭回归
也叫吉洪诺夫正则化，是线性回归的正则化版。加入L2正则项 α ∑ i = 1 n θ i 2 \alpha \sum_{i=1}^n \theta_i^2 α∑i=1nθi2。
岭回归成本函数：
J ( θ ) = M S E ( θ ) + α 1 2 ∑ i = 1 n θ i 2 J(\theta) = MSE(\theta)+\alpha\frac{1}{2} \sum_{i=1}^n \theta_i^2 J(θ)=MSE(θ)+α21i=1∑nθi2
注意：这里偏置项 θ 0 \theta_0 θ0没有正则化，求和是从 θ 1 \theta_1 θ1 开始的。如果我们将w定义为特征权重的向量( θ 1 到 θ n \theta_1到\theta_n θ1到θn)，那么正则项即等于 1 2 ( ∥ w ∥ 2 ) 2 \frac{1}{2}(\left \|w \right \|_2)^2 21(∥w∥2)2，其中 ∥ w ∥ 2 \left \|w \right \|_2 ∥w∥2为权重向量的l2范数。
注意：在执行岭回归之前，要对特征进行缩放 (例如StandardScaler)，因为它对输入特征的大小非常敏感，大多数正则化模型都是如此。
下面使用多个
α
\alpha
α对某线性数据进行训练的两种岭回归模型。左图为直接使用岭回归模型，右图为岭正则化后的多项式回归。
上图代码：
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) #获取0到3之间的数值，符合等差数列
#绘制模型图
# model_class 模型
# polynomial 是否高阶
# alphas 学习率取值列表
# model_kargs 模型其他参数
def plot_model(model_class, polynomial, alphas, **model_kargs):
for alpha, style in zip(alphas, ("b-", "g--", "r:")):
#如果学习率大于0就用正则化模型，否则就线性回归模型（alpha = 0时其实就是LinearRegression）
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_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)
#直接使用岭回归
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()
岭回归也有闭式解：
θ ^ = ( X T ⋅ X + α A ) − 1 ⋅ X T ⋅ y \hat\theta = (X^T\cdot X+\alpha A)^{-1}\cdot X^T\cdot y θ^=(XT⋅X+αA)−1⋅XT⋅y
使用sklearn执行闭式解的岭回归 ( 利用上面公式的变体：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]])
使用随机梯度下降
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]])
上面2种方式都属于岭回归。
5.2 套索回归
线性回归的另一种正则化，简称Lasso，其中La表示Least Absolute的意思，即最小绝对值。它也是向成本函数增加一个正则项，但是它增加的是权重向量的L1范数。
Lasso回归成本函数:
J ( θ ) = M S E ( θ ) + α ∑ i = 1 n ∣ θ i ∣ J(\theta)=MSE(\theta)+ \alpha\sum_{i=1}^n|\theta_i| J(θ)=MSE(θ)+αi=1∑n∣θi∣
Lasso回归的一个重要特点：倾向于完全消除掉最不重要的特征权重（也就是设置为0）。如下图，当 α = 1 0 − 7 \alpha=10^{-7} α=10−7时，快接近于线性：因为所有高阶多项式的特征权重都等于0。换句话说，Lasso回归会输出一个稀疏模型，只有很少的特征有非零权重。
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()
sklearn的Lasso实现
from sklearn.linear_model import Lasso
lasso_reg = Lasso(alpha=0.1)
lasso_reg.fit(X, y)
lasso_reg.predict([[1.5]])
5.3 弹性网络
顾名思义，所谓弹性，指的就是它位于岭回归和lasso回归之间的中间地带。混合的比例通过r来控制，r=0时，弹性网络即等于岭回归。
弹性网络成本函数：
J ( θ ) = M S E ( θ ) + r α ∑ i = 1 n ∣ θ i ∣ + 1 − r 2 α ∑ i = 1 n θ i 2 J(\theta)=MSE(\theta)+ r\alpha\sum_{i=1}^n|\theta_i|+\frac{1-r}{2} \alpha\sum_{i=1}^n \theta_i^2 J(θ)=MSE(θ)+rαi=1∑n∣θi∣+21−rαi=1∑nθi2
如何选择线性回归、岭回归、lasso回归和弹性网络呢？
① 通常来说，有正则化总比没有要好。所以岭回归可以是默认选择。
② 若实际用到的特征很少，那就应该更倾向于lasso回归或弹性网络，因为它们会将无用特征的权重都降维0.
③ 一般情况下，弹性网络优于lasso，因为当特征数量超过实例数量，或几个特征强相关时，lasso可能非常不稳定。
所以，一般来说，默认用岭回归，第②种情况用弹性网络。
sk种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]])
5.4 早期停止法
针对梯度下降这类迭代学习的算法而言，早停是另一种正则化方法：在验证集误差达到最小时停止迭代训练。
上图代码
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)
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) # mse得分
if val_error < minimum_val_error: #
minimum_val_error = val_error
best_epoch = epoch
best_model = clone(sgd_reg)
注意：当warm_start=True时，调用fit()，会从停下的地方继续开始训练，而不会重新开始。
5.5 逻辑回归
其实就是用线性的方式去解决分类问题，更详细通俗的请参考我之前写的博客逻辑回归。
逻辑回归的成本函数没有闭式方程，不能直接求出 θ \theta θ，但是它是一个凸函数，可以使用梯度下降等算法求出全局最小值。
5.5.1 概率估算
逻辑回归也是计算输入特征的加权和，但是不同的是它输出的是一个概率 （0 - 1之间的数字）。
5.5.2 训练和成本函数
成本函数为log损失函数，具体参考之前博客。
5.5.3 决策边界
略
5.9 Softmax回归
也叫多元逻辑回归，它是逻辑回归的扩展，支持多个类别。
原理步骤：对于一个给定的实例
① Softmax先计算出每个类别k的分数 s k ( x ) s_k(x) sk(x)；
类别k的Softmax分数： s k ( x ) = θ k T ⋅ x s_k(x)=\theta_k^T\cdot x sk(x)=θkT⋅x
② 对这些分数应用Softmax函数（也叫归一化函数），估算出每个类别的概率。
Softmax分数归一化（除以所有指数的总合）即得到 p ^ k \hat p_k p^k ：
p ^ k = h ( s ( x ) ) k = e x p ( s k ( x ) ) ∑ j = 1 k e x p ( s j ( x ) ) \hat p_k = h(s(x))_k = \frac {exp(s_k(x))}{\sum_{j=1}^kexp(s_j(x))} p^k=h(s(x))k=∑j=1kexp(sj(x))exp(sk(x))
K是类别的数量；
s(x)是实例x每个类别的分数的向量；
h ( s ( x ) ) k h(s(x))_k h(s(x))k是实例x属于类别k的概率。
③ 最终的Softmax分类器预测函数为：
y ^ = a r g m a x ( k ) h ( s ( x ) ) k = a r g m a x ( k ) ( s k ( x ) ) = a r g m a x ( k ) ( θ k T ⋅ X ) \hat y = argmax_{(k)} \ h(s(x))_k = argmax_{(k)} \ (s_k(x))=argmax_{(k)} \ (\theta_k^T\cdot X) y^=argmax(k) h(s(x))k=argmax(k) (sk(x))=argmax(k) (θkT⋅X)
argmax返回的是使函数最大化所对应的变量的值。 在上面的等式里，返回的是使估算概率
h
(
x
(
x
)
)
k
h(x(x))_k
h(x(x))k最大时 k的值。
注意：Softmax回归分类器一次只会预测一个类别（它是多类别，但不是多输出，每次只会单输出众多类别中的一个），比如，用它识别照片中的多个人是不行的。
Softmax成本函数（交叉熵），为啥使用交叉熵？因为交叉熵经常用于衡量一组估算类别概率和目标类别的匹配程度。公式：
J ( Θ ) = − 1 m ∑ i = 1 m ∑ k = 1 k y k ( i ) l o g ( p ^ k ( i ) ) J(\Theta) = - \frac{1}{m}\sum_{i=1}^{m}\sum_{k=1}^{k}y_k^{(i)}log(\hat p_k^{(i)}) J(Θ)=−m1i=1∑mk=1∑kyk(i)log(p^k(i))
如果第i个实例的目标类别为k，则 y k ( i ) y_k^{(i)} yk(i) = 1，否则为0。
其中大写的theta是参数矩阵，存储的是每个类别自己的参数向量 θ k \theta_k θk 。
注意：当只有2个类别时(k=2)，它等价于逻辑回归的成本函数(log损失函数)。
类别k的交叉熵梯度向量（偏导）：
∇ θ k J ( Θ ) = 1 m ∑ i = 1 m ( p ^ k ( i ) − y k ( i ) ) x ( i ) \nabla_{\theta_k}J(\Theta)=\frac{1}{m}\sum_{i=1}^m(\hat p_k^{(i)}-y_k^{(i)})x^{(i)} ∇θkJ(Θ)=m1i=1∑m(p^k(i)−yk(i))x(i)
现在就可以计算每个类别的梯度向量，然后使用梯度下降找到最小化损失的参数矩阵 Θ \Theta Θ了。
下面使用Softmax回归将鸢尾花分为三类。默认情况下，LR是一对多的训练方式，将超参数multi_class设置为"multinomial"就可以切换成Softmax回归，还要必须指定一个支持Softmax回归的求解器，比如lbfgs求解器。默认l2正则，可以通过超参数C进行控制。
from sklearn.linear_model import LogisticRegression
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)
# 给你一朵花，花瓣长5 宽2
# Softmax预测
softmax_reg.predict([[5,2]]) #输出为第2类Virginica鸢尾花
softmax_reg.predict_proba([[5,2]]) # 94.2%为第2类Virginica鸢尾花，5.8%概率为Versicolor鸢尾花
决策边界
x0, x1 = np.meshgrid(#坐标矩阵
np.linspace(0, 8, 500).reshape(-1, 1),# 0-8区间生成500个符合等差数列的数字
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()
412
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
