线性模型(linear model) 试图学得一个通过属性的线性组合来进行预测的函数，即：其中，当确定和之后，整个模型即可确定。
一、线性回归
线性回归是其中最经典的一种回归模型。
1线性回归的求解
求解可分为梯度下降法和最小二乘法。
1.1 梯度下降法
假设函数：为了统一，令，即可得到假设函数。
损失函数即样本值与实际值的差值的平方（欧氏距离）：
代价函数是所有样本点的损失求和，通过均方误差（mean-square error, MSE）来定义，在线性回归中，最小二乘法就是试图找到一条直线，使所有样本到直线上的欧氏距离之和最小。分母中的2用于后续求导的方便：
显然我们希望损失函数最小，即模型预测的结果与实际对应结果的误差最小，因此对求导：为样本数量，表示第个特征的第条数据（矩阵的第行第列），根据上式有迭代公式：其中，为步长（学习率）代表每次沿梯度下降的幅度，当迭代次数足够大时，损失函数最终会收敛到一个局部的极值点。
对于上述的迭代公式，分别代入不同的，即可得到对应的。
举例，考虑两个特征的情况，训练数据集如下：
此时有：给定初始值，步长，第一次迭代：同理可以更新，经多次迭代，最终可以得到一组使损失函数最小的。
训练得到模型参数以后，便可以使用模型进行预测，下面是具体的代码实现：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
| import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model
def load_data():
"""
加载数据集
:return:
"""
diabetes_X, diabetes_y = datasets.load_diabetes(return_X_y=True)
diabetes_y = np.reshape(diabetes_y, (-1, 1))
diabetes_X = np.reshape(diabetes_X[:, 2], (-1, 1))
X_train = diabetes_X[:-20]
X_test = diabetes_X[-20:]
y_train = diabetes_y[:-20]
y_test = diabetes_y[-20:]
return X_train, X_test, y_train, y_test
def cost_compute(X, y, theta):
"""
计算损失值
:param X: X
:param y: y
:param theta: 特征参数
:return:
"""
cost = np.sum(np.power(np.dot(X, theta) - y, 2)) / 2 / X.shape[0]
return cost
def gradient_descent(X, y, theta, iter_num, eta):
"""
线性回归梯度下降算法
:param X: X
:param y: y
:param theta: 特征参数
:param iter_num: 迭代次数
:param eta: 学习率（步长）
:return:
"""
new_col = np.ones(X.shape[0]).T
X = np.insert(X, 0, new_col, axis=1)
n = X.shape[0]
m = X.shape[1]
costs = np.zeros(iter_num)
for i in range(iter_num):
costs[i] = cost_compute(X, y, theta)
for j in range(m):
sum_ = np.sum((y - np.dot(X, theta)) * X[:, np.newaxis, j])
theta[j] = theta[j] + eta / n * sum_
return theta, costs
def norm_X(X):
return (X - np.mean(X)) / np.std(X, ddof=1)
def predict(X, theta):
"""
预测函数
:param X: 测试数据
:param theta: 参数
:return:
"""
new_col = np.ones(X.shape[0]).T
X = np.insert(X, 0, new_col, axis=1)
return np.dot(X, theta)
X_train, X_test, y_train, y_test = load_data()
X_train = norm_X(X_train)
iter_num = 100
eta = 0.05
theta = np.zeros(X_train.shape[1] + 1).reshape((-1, 1))
theta, costs = gradient_descent(X_train, y_train, theta, iter_num, eta)
X_test = norm_X(X_test)
predict_y = predict(X_test,theta)
plt.figure(1)
plt.scatter(X_train, y_train, color="black")
plt.plot(X_train, theta[0] + theta[1] * X_train, color="blue", linewidth=2)
plt.show()
plt.figure(1)
x_axis = np.linspace(1, iter_num, iter_num)
plt.plot(x_axis, costs[0:iter_num], color="red", linewidth=2)
plt.show()
|
1.2 梯度下降的几种方式
梯度下降可以分为：
- 批量梯度下降（Batch Gradient Descent,BGD）
- 随机梯度下降（Stochastic Gradient Descent,SGD）
- 小批量梯度下降（Mini-Batch Gradient Descent,MBGD）
上一节所述的公式为批量梯度下降：每次更新时，使用全部样本进行计算，在小样本的情况下，这种方式更新的幅度相对较大；在样本量多的时候，计算性能上会比较慢。
随机梯度下降是指每次更新时只用使用一个样本：即每次更新时不使用全量样本，而是随机选择一个样本近似所有的样本。这种方式的优点是计算相比BGD快，但带来的问题是每次迭代得到的损失函数并不总是朝着全局最优方向（总体仍然向最优方向）下降，可能会出现摇摆的情况。
小批量梯度下降是前两者的折中，
其中的是全量样本中的一小部分，比如每次更新时选择30个样本，用这些样本近似全部的样本，既考虑了运行效率问题，同时兼顾准确性的问题。
1.3 一元线性回归最小二乘法
另外一种求解方式是最小二乘法。
首先从简单开始，假设特征只有一个，则称为一元线性回归，方程如下：下面的任务就是确定和，对于这个方程来说，就是二维平面上的一条直线。我们通过的一系列取值确定参数，使得拟合的直线尽可能地趋近于真实值。问题的关键是如何衡量和的差别。
由于均方误差是回归任务中最常用的性能度量，它反映了估计量与被估计量之间差异程度的一种度量。自然，我们希望估计的均方误差越小越好，因此希望求得使此函数最小化的。这个过程称为线性回归模型的最小二乘参数估计，同样分别对和求导，得到：和分别令上式为零，即可解得和的最优闭式解。计算过程如下，令（1）式为0，得：令（2）式为0得：又因为样本均值的定义：，，则，代入上式可得：由于，，代入上式即可得：
1.4 多元线性回归最小二乘法
更一般地，考虑个样本，个属性的情形，此时我们试图学得：类似的，可利用最小二乘法来对和进行估计。对于多元线性回归来说，我们可以类似得到为了便于讨论，将和吸入向量形式：将写成向量形式：此时（3）式可以简化为：由向量内积的定义，上式可以写成：
其中：（4）式中，我们将数据集表示为一个大小的矩阵，前个元素对应属性值，最后一列恒置为1：同样将吸入向量形式：所以有：令，对求导，首先将式子展开：求导：这里根据矩阵的求导公式（推导过程见矩阵的导数运算）：得到：令上式为零可得最优解的闭式解，但这里涉及到矩阵的逆运算，需要额外讨论一下。
当为满秩矩阵或正定矩阵时，令上式为零得：其中，为的逆矩阵，令，则最终学得多元线性回归模型为：
当非满秩时，例如在许多任务中我们会遇到大量的特征，其数目甚至超过样本数，导致的列数多于行数，此时矩阵显然不满秩，根据线性方程组的解的情况可知，方程组有无穷多解。即存在多个，它们都能使均方误差最小化，此时选择哪个解作为输出将由学习算法的归纳偏好决定，常见的做法是为求解函数引入惩罚项。例如引入正则化项，在线性回归中，引入L1范数的模型称为LASSO回归，引入L2范数的模型称为岭回归。
依据公式，核心代码实现：
1
2
3
4
5
6
7
8
9
10
11
12
| def least_square_method(X, y):
"""
最小二乘法
:param X: 训练样本
:param y: y
:return:
"""
new_col = np.ones(X.shape[0]).T
X = np.insert(X, 0, new_col, axis=1)
theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
return theta
|
1.5 最小二乘法和梯度下降法的比较
梯度下降法是一种数值优化方法，需要设置初始值，并进行迭代，如果初始值或迭代次数设置不合理，会导致无法收敛的情况发生，类似的方法还有牛顿法等。
最小二乘法不需要设置初始值，直接通过公式即可计算出参数（直接得到解析解）。但是，注意到公式中存在矩阵的逆运算，当特征较多时，计算一个矩阵的逆非常耗费时间，另外，当非满秩（奇异矩阵）时，需要进行特殊处理才可求逆，常用的做法是加上一个对角阵（推导见下文岭回归）：其中为与同大小的单位阵。
在具体的求解中，需要具体问题具体分析，采用不同的求解方法。
2 岭回归求解
2.1 欠拟合和过拟合
当学习器把训练样本学得太好了的时候，就会出现过拟合现象，这样就会导致泛化性能下降。与过拟合相对应的是欠拟合，通常是由于学习能力低下而造成的。
欠拟合比较容易克服，而过拟合则很麻烦。过拟合是机器学习面临的关键障碍，各类学习算法都必然带有一些针对过拟合的措施，然而必须认识到，过拟合是无法彻底避免的，我们能做的只是“缓解”。
2.2 正则化项缓解过拟合
通过在损失函数中添加正则化项，可以减少求出过拟合解的可能性。
对于来说，系数表示了与之相关的特征在整个式子中的权重。显然所有的特征权重均相同，考虑到模型训练中的很多特征是没必要的特征，这些对于模型来说，这些特征就是噪声，因此希望通过调整每个对应的所占权重，只保留重要特征的从而减少模型复杂度，进而达到规避过拟合的目的。
添加正则化项后，可以看做在控制损失函数不变的情况时令正则项最小化，相当于为损失函数添加了一个约束条件，从图像上直观理解就是参数不仅要满足损失函数取得最小值，同时要满足在菱形或圆形内部（分别对应1范数和2范数）。
上面左图，当损失函数取到最小值时，说明这一项被抵消掉；上面右图，当损失函数取到最小值时变得很小，接近于0，而还比较大，说明这一项的影响权重大大降低。
添加正则化项的损失函数如下：其中，为损失函数，为正则化系数，为正则化函数，其一般是一个过原点的凸函数。
范数是机器学习中常用的正则项，它的一般形式如下：在图形上，随着越大，等高线图越接近正方形（正无穷范数）；越小，曲线弯曲越接近原点（负无穷范数）。
当越小，则最优解越容易趋近于某个坐标轴，意味正交于该坐标轴的权值趋向于0，即参数稀疏化；越大则约束项对各个维度的权值约束力趋于相同，最终的效果为等比例地缩小各个权值。
2.3 岭回归推导
岭回归就是在原有的损失函数上添加L2范数正则化：对对求导，前面的结果和普通线性回归求导结果相同，只是后半部分添加了正则项的求导：岭回归的迭代公式：核心代码实现：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
| def ridge_gradient_descent(X, y, theta, iter_num, eta, lambda_):
"""
岭回归梯度下降算法
:param X: X
:param y: y
:param theta: 特征参数
:param iter_num: 迭代次数
:param eta: 学习率（步长）
:param lambda_: 岭回归系数
:return:
"""
new_col = np.ones(X.shape[0]).T
X = np.insert(X, 0, new_col, axis=1)
n = X.shape[0]
m = X.shape[1]
costs = np.zeros(iter_num)
for i in range(iter_num):
costs[i] = cost_compute(X, y, theta)
for j in range(m):
sum_ = np.sum((y - np.dot(X, theta)) * X[:, np.newaxis, j])
theta[j] = theta[j] + eta / n * (sum_ - 2 * lambda_ * theta[j])
return theta, costs
|
最小二乘法的岭回归推导类似，我们为添加L2范数：求导得：令上式为零，得到：
4 使用sklearn求解线性回归
sklearn中提供了普通最小二乘法求解线性回归方法，LinearRegression 定义模型，通过fit方法接收X和y进行训练，并将参数保存在intercept_和coef_中，分别对应和。
代码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
| import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model
def load_data():
"""
加载数据集
:return:
"""
diabetes_X, diabetes_y = datasets.load_diabetes(return_X_y=True)
diabetes_y = np.reshape(diabetes_y, (-1, 1))
diabetes_X = np.reshape(diabetes_X[:, 2], (-1, 1))
X_train = diabetes_X[:-20]
X_test = diabetes_X[-20:]
y_train = diabetes_y[:-20]
y_test = diabetes_y[-20:]
return X_train, X_test, y_train, y_test
X_train, X_test, y_train, y_test = load_data()
model = linear_model.LinearRegression()
model.fit(X_train,y_train)
print(model.coef_)
print(model.intercept_)
plt.figure(1)
plt.scatter(X_train, y_train, color="black")
plt.plot(X_train, model.intercept_ + model.coef_ * X_train, color="blue", linewidth=2)
plt.show()
|
对于岭回归，sklearn提供了Ridge模型，其余部分不变，只需要将模型类替换为Ridge，并添加与正则项相乘的常数alpha即可，例如：
1
| model = linear_model.Ridge(alpha=.5)
|
对于LASSO，将其替换为Lasso即可：
1
| model = linear_model.Lasso(alpha=0.1)
|
对于弹性网模型，使用ElasticNet，其中alpha参数不变，增加了另外一个混合参数l1_ratio代表L1的权重，取值范围为0 <= l1_ratio <= 1 。对于 l1_ratio = 0 ，惩罚是L2惩罚。 For l1_ratio = 1 是L1惩罚。对于 0 < l1_ratio < 1 ，惩罚是 L1 和 L2 的组合。
1
| model = linear_model.ElasticNet(alpha=0.01, l1_ratio=0.7)
|
更多参数可查看sklearn官方文档。
二、逻辑回归
1 广义线性模型与对数线性回归
对于线性回归：可否令模型预测值逼近的衍生物呢，假设我们认为示例所对应的输出标记是在指数尺度上变化，那就可将输出标记的对数作为线性模型逼近的目标，即：这就是对数线性回归（逻辑回归）。上式在形式上虽然是线性回归，但实质上已是在求取输入空间到输出空间的非线性函数映射，它实际上是在试图让逼近。
更一般地，对于单调可微函数，令这样的模型称为广义线性模型，起到了将线性回归模型的预测值与真实标记联系起来的作用，被称为“联系函数”，显然对数几率回归是的特例。
对于广义线性模型的定义，在后面还会继续介绍。
2 二分类任务
对于分类任务，根据（6）式，只需要找到一个单调可微函数将分类任务的真实标记与线性回归模型的预测值联系起来。
首先考虑二分类，输出标记为，即正例与反例。于是需要将实值转换为0或1，最理想的函数是“单位阶跃函数”（unit-step function）：若预测值小于零则判别为反例，大于零判别为正例，等于零则任意判别。
但是由于单位阶跃函数并不连续，不能直接用作（6）式中的，于是我们希望找到能在一定程度上近似单位阶跃函数的替代函数，并希望它单调可微。对数几率函数正是常用的替代函数：
如图，对数几率函数是一种“Sigmoid函数“，它将值转化为一个接近0或1的值，并且其输出值在附近变化很陡峭，将其代入到公式（6）中的，得到：将输出标记的对数作为线性模型的目标函数：将视为正例的可能性，视为反例的可能性，则称为几率，反映了作为正例的相对可能性，取对数则得到对数几率：将这种模型对数几率回归（logistic regression，或logit regression），有些文献又称其为逻辑回归。此函数是任意阶可导的凸函数，有很好的数学性质，现有的许多数值优化算法都可直接用于求取最优解。
将视为类后验概率估计，则（7）式可以重写为：此时有：
3 通过广义线性模型建模
前面对于线性回归和逻辑回归都进行了建模，但是为什么要使用Sigmoid仍然没有明确的解释，对于这两种回归模型，实际上可以统一到一种更加抽象的模型中，也就是广义线性模型。
下面通过广义线性模型对逻辑回归进行建模，通过这种方式，可以更好地理解为什么逻辑回归要用sigmoid函数。
首先是指数分布族的定义，指数分布族是这样一类概率分布，它们的概率密度或分布律满足如下形式：其中，是分布的自然参数，为充分统计量，通常情况下，为对数配分函数，是关于随机变量的函数。
0-1分布（伯努利分布）、二项分布、泊松分布、高斯分布（正态分布）等都属于指数分布族。
广义线性模型建立在三个定义的基础上，分别为：
- 假设变量服从指数分布族的某个分布
- 模型的目标是，即得到一个模型可以预测出的期望值
- 指数分布族中的自然参数，即和呈线性关系
显然，对于二分类问题，可以假设数据服从伯努利分布，伯努利分布的分布律为：
可以证明伯努利分布属于指数分布族。对上式进行变换：对比（9）式，当，伯努利分布符合指数分布族形式。
其中：
依据广义线性模型的定义，伯努利分布属于指数分布。根据广义线性模型的第二条，目标函数：服从伯努利分布，根据期望的定义，，以及（10）式，得到：这就是Sigmoid函数的由来，下面根据广义线性模型第三条定义，，最终得到：这与（8）式建模的模型相同。
4 极大似然法参数估计、梯度下降求解逻辑回归
4.1 极大似然法
为了方便，将（7）式中的吸入向量，令，则。
（9）式可以变为：根据极大似然估计法，写出随机变量的分布律表达式：则对数似然函数为：考虑到可能的取值只有0和1，所以：可以看出在分别取0和取1时，仅相差一个，因此，综合上式，我们为其乘以得到，此时对数似然函数变为：这里求对数似然函数的最大值，相当于求的最小值，因此可将定义为新的损失函数。
以上就是损失函数的推导过程。实际上还有另外一种推导方式，随机变量的分布律表达式不止可以使用（11）式描述，下面是另外的表述方法，同样满足伯努利分布率：同理，此时的对数似然函数为：显然，这种方式更容易推导出式（13）。
4.2 梯度下降求解逻辑回归
下面对损失函数求导：因此有迭代公式：代码实现：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
| import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets
def gradient_descent(X, y, theta, iter_num, eta):
new_col = np.ones(X.shape[0]).T
X = np.insert(X, 0, new_col, axis=1)
n = X.shape[0]
m = X.shape[1]
y = np.reshape(y, (-1, 1))
for i in range(iter_num):
for j in range(m):
sig = 1 / (1 + np.exp(-np.dot(X, theta)))
sum_ = np.sum(X[:, np.newaxis, j] * (y - sig))
theta[j] = theta[j] + eta / n * (sum_)
return theta
def plot(X, theta):
plt.scatter(X[:, 0], X[:, 1], c=Y, edgecolors="k", cmap=plt.cm.Paired)
x1 = np.arange(min(X[:, 0]), max(X[:, 0]), 0.1)
x2 = -(theta[1] * x1 + theta[0]) / theta[2]
plt.plot(x1, x2)
plt.show()
iris = datasets.load_iris()
X = iris.data[:, :2]
Y = iris.target
X = X[Y != 2]
Y = Y[Y != 2]
iter_num = 10000
eta = 0.05
theta = np.zeros(X.shape[1] + 1).reshape((-1, 1))
theta = gradient_descent(X, Y, theta, iter_num, eta)
plot(X, theta)
|
5 逻辑回归的正则化项
同样可以为逻辑回归损失函数中添加正则化项。下面是加入2范数的正则化项后的迭代公式：代码上的实现也很简单，仅添加一项即可：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
| def gradient_descent(X, y, theta, iter_num, eta, lambda_):
new_col = np.ones(X.shape[0]).T
X = np.insert(X, 0, new_col, axis=1)
n = X.shape[0]
m = X.shape[1]
y = np.reshape(y, (-1, 1))
for i in range(iter_num):
for j in range(m):
sig = 1 / (1 + np.exp(-np.dot(X, theta)))
sum_ = np.sum(X[:, np.newaxis, j] * (y - sig))
theta[j] = theta[j] + eta / n *(sum_ + 2 * lambda_ * theta[j])
return theta
|
6 多分类任务
考虑个类别，多分类学习的基本思路是拆解法，即将多分类任务拆为若干个二分类任务求解。具体来说，先对问题进行拆分，然后为拆出的每个二分类任务训练一个分类器。在测试时，对这些分类器的预测结果进行集成以获得最终的多分类结果。
这里的关键是如何对多分类任务进行拆分，以及如何对多个分类器进行集成。西瓜书提到的经典的拆分策略有三种，分别是一对一（OvO）、一对其余（OvR）和多对多（MvM）。
一对一策略
给定数据集和类别，OvO将这个类别两两配对，即进行组合，从而产生个二分类任务，每个分类器都会产生一个正例和一个反例。在预测阶段，新样本将同时提交给所有分类器，通过这些分类器最终得到个结果，最终的结果通过投票产生：将被预测最多的结果当做最终结果。
一对其余策略
OvR则是每次将一个类的样例作为正例、所有其他类的样例作为反例来训练个分类器。在预测时，若仅有一个分类器预测为正类，则对应的类别标记作为最终分类结果；若有多个分类器预测为正类，则通常考虑各分类器的预测置信度，选择置信度最大的类别标记作为分类结果。容易看出，OvR训练的分类器通常少于OvO，因此训练所需的开销也更小。至于预测性能则取决于具体的数据分布， 在多数情形下两者差不多。
多对多策略
MvM是每次将若干个类作为正类，若干个其他类作为反类，这种正反类的构造通常需要进行特殊编码，称为纠错输出码，比如ECOC编码。首先，类别划分通过“编码矩阵“确定，编码矩阵有多种形式， 常见的主要有二元码、三元码。前者将每个类别分别指定为正类和反类，后者在正、反类之外，还可指定停用类。参考西瓜书中图：
对于左图中的二元码，分类器将类作为正例，其余作为反例。在解码阶段，各分类器的预测结果联合起来形成了测试示例的编码，该编码与各类所对应的编码进行比较，将距离最小的编码所对应的类别作为预测结果。距离的计算可以通过欧氏距离，也可以通过海明距离得出。比如图中(a)的解码若根据欧氏距离，则输出作为预测结果。
ECOC 编码对分类器的错误有一定的容忍和修正能力， 上图(a)中对测试示例的正确预测编码是，假设在预测时某个分类器出错了，例如出错从而导致了错误编码，但基于这个编码仍能产生正确的最终分类结果。一般来说，对同一个学习任务， ECOC编码越长，纠错能力越强。然而，编码越长，意味着所需训练的分类器越多，计算、存储开销都会增大；另一方面，对有限类别数可能的组合数目是有限的，码长超过一定范围后就失去了意义。
7 sklearn实现逻辑回归
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
| import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.inspection import DecisionBoundaryDisplay
from sklearn.linear_model import LogisticRegression
iris = datasets.load_iris()
X = iris.data[:, :2]
Y = iris.target
logreg = LogisticRegression(C=1e5,max_iter=150,multi_class="auto")
logreg.fit(X, Y)
_, ax = plt.subplots(figsize=(4, 3))
DecisionBoundaryDisplay.from_estimator(
logreg,
X,
cmap=plt.cm.Paired,
ax=ax,
response_method="predict",
plot_method="pcolormesh",
shading="auto",
xlabel="Sepal length",
ylabel="Sepal width",
eps=0.5,
)
plt.scatter(X[:, 0], X[:, 1], c=Y, edgecolors="k", cmap=plt.cm.Paired)
plt.xticks(())
plt.yticks(())
plt.show()
|
alipay
wechat