[数据结构与算法]线性回归-从理论到实践

一元一次线性回归

抽象问题

现有一组数据，共有两列，分别为x和y，如下所示

现将数据用python作散点图，观察其变化趋势

import pandas as pd 
import numpy as np
import matplotlib.pyplot as plt

# 读取数据并转为np.array的格式
train = pd.read_csv('click.csv').values 
# 读取所有行的第0列作为训练数据的x
train_x = train[:, 0] 
# 读取所有行的第1列作为训练数据的y
train_y = train[:, 1]
# 画散点图
plt.plot(train_x, train_y, 'o')
plt.show()

我们的问题是要找到x和y的对应关系，根据此对应关系，我们就可以在给定任意x的情况下，计算得到y的值。

定义模型

观察上面的散点图可以看出来，随着x的增大，y值也在近似的不断增大。据此，我们推测，这组数据可能符合一次函数的关系，即：$y=kx+b$?

此式可改写为：$$f_{\theta }(x)={\theta }_0+{\theta }_1?x$$?

其中$f_{\theta }(x)$ 表示含有参数 $\theta$ ，与变量$x$相关的函数。该函数即为：一元一次线性回归模型。

注：统计学中，常用 $\theta$ 来表示未知数和推测值，采用 $\theta$???? 加数字下标的形式，是为了防止当未知数的数量增加时，表达式中大量出现a、b、c、d…这样的符号。这样不但不易理解，还容易出现符号本身不够用的情况。

只要求得了上述公式中的未知参数${\theta }_0$和${\theta }_1$，即知道了x和y的关系。

现在问题转化为求未知参数${\theta }_0$和${\theta }_1$

求解方法

已知1个自变量$x^{(i)}$????，根据模型公式$f_{\theta }(x)={\theta }_0+{\theta }_1?x$???，可以求得与其对应的因变量$f_{\theta }(x^{(i)})$????。

那么 $y^{(i)}?f_{\theta }(x^{(i)})$?表示实际值与模型预测值之间的误差。

我们要求解的模型函数最理想的情况就是使得$f_{\theta }(x^{(i)})=y^{(i)}$，即误差值$y^{(i)}?f_{\theta }(x^{(i)})=0$?。

实际中，很难做到让所有点的误差都等于0，因此，我们要做的就是让所有点的误差之和尽可能的小。

我们一般使用误差的平方来表征误差，一方面是去除误差正负抵消的影响，一方面是为了使之后求解微分更加方便。

所有点的误差的平方和可以用如下公式来表示：

$$E(\theta )=\frac{1}{2}\sum _{i=1}^n{(y^{(i)}?f_{\theta }(x^{(i)}))}^2$$

这个表达式称为目标函数，$E$????是Error的首字母。

现在，我们的问题转化为：求出使得目标函数$E(\theta )$取最小值的参数${\theta }_0$和${\theta }_1$

梯度下降法(最速下降法)

现有函数$g(x)$根据导数的意义：

1. 当$\frac{d}{dx}g(x)>0$????，?即导数为正，$g(x)$随x的增大而增大，因此，要使得$g(x)$减小，需要x负移
2. 当$\frac{d}{dx}g(x)<0$??，即导数为负，$g(x)$随x的增大而减小，因此，要使得$g(x)$减小，需要x正移

据此可知，要使得$g(x)$不断减小，只需要将x向与导数符号相反的方向移动即可。用数学表达式表述为：

$$x:=x?\eta \frac{d}{dx}g(x)$$

1. $A:=B$，表示通过B来定义A
2. $\eta$??，表示学习率。
   1. $\eta$越小，更新次数越多，速度越慢
   2. $\eta$越大，更新次数越少，速度越快，但有可能导致不收敛。

梯度下降法的核心：根据导数的符号来决定x移动的方向。

根据梯度下降法的原理，我们知道了怎样更新参数${\theta }_0$和${\theta }_1$，就能使目标函数$E(\theta )$?的值达到最小。

参数更新方法如下

$$\{\begin{array}{l}{\theta }_0:={\theta }_0?\eta \frac{?}{?{\theta }_0}E(\theta )\\ {\theta }_1:={\theta }_1?\eta \frac{?}{?{\theta }_1}E(\theta )\end{array}$$

现在只要计算出$\frac{?}{?{\theta }_0}E(\theta )$?和$\frac{?}{?{\theta }_1}E(\theta )$?，代入上式，就可不断更新参数，求得最佳的${\theta }_0$?和${\theta }_1$?。

计算上式中的两个偏导数过程如下：

$$\frac{?}{?{\theta }_0}E(\theta )=\frac{1}{2}\stackrel{n}{\sum _{i=1}}2?(y^{(i)}?f_{\theta }(x^{(i)}))?\frac{?(y^{(i)}?f_{\theta }(x^{(i)}))}{?{\theta }_0}$$
$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?\frac{?(f_{\theta }(x^{(i)}))}{?{\theta }_0}$$??
$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?\frac{?({\theta }_0+{\theta }_1?x^{(i)})}{?{\theta }_0}$$?
$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?1$$
$$=\stackrel{n}{\sum _{i=1}}(f_{\theta }(x^{(i)})?y^{(i)})$$

同理：

$$\frac{?}{?{\theta }_1}E(\theta )=\frac{1}{2}\stackrel{n}{\sum _{i=1}}2?(y^{(i)}?f_{\theta }(x^{(i)}))?\frac{?(y^{(i)}?f_{\theta }(x^{(i)}))}{?{\theta }_1}$$?
$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?\frac{?(f_{\theta }(x^{(i)}))}{?{\theta }_1}$$???
$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?\frac{?({\theta }_0+{\theta }_1?x^{(i)})}{?{\theta }_1}$$??
$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?x^{(i)}$$
$$=\stackrel{n}{\sum _{i=1}}(f_{\theta }(x^{(i)})?y^{(i)})?x^{(i)}$$

因此，参数更新表达式如下：

$$\{\begin{array}{l}{\theta }_0:={\theta }_0?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})\\ {\theta }_1:={\theta }_1?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})?x^{(i)}\end{array}$$

根据该式，不断迭代，直到求得最佳的参数${\theta }_0$和${\theta }_1$，即确定了模型表达式$f_{\theta }(x)={\theta }_0+{\theta }_1?x$

python实现

根据上述的理论及推导，我们知道：

1. 模型表达式为：$f_{\theta }(x)={\theta }_0+{\theta }_1?x$
2. 目标函数为：$E(\theta )=\frac{1}{2}{\sum }_{i=1}^n{(y^{(i)}?f_{\theta }(x^{(i)}))}^2$
3. 参数更新表达式为：$\{\begin{array}{l}{\theta }_0:={\theta }_0?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})\\ {\theta }_1:={\theta }_1?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})?x^{(i)}\end{array}$

在进行参数更新之前，应该将数据标准化（也称z-score规范化）。这个步骤非必须，但是做了之后，参数的收敛会更快

$$z^{(i)}=\frac{x^{(i)}?\mu }{\sigma }$$

# 读取数据并转为np.array的格式
train = pd.read_csv('click.csv').values 
train_x = train[:, 0] 
train_y = train[:, 1]

# 对变量train_x进行标准化
mu = train_x.mean()
sigma = train_x.std()
def standardize(x):
    return (x-mu) / sigma
train_x_std = standardize(train_x)

python逐步实现

# 1,随机生成初始参数theta0 和 theta1
theta0, theta1 = np.random.rand(2)

# 2,定义模型函数f(x)
def f(x):
    y = theta0 + theta1 * x
    return y

# 3,定义目标函数E(x, y)
def E(x, y):
    e = 0.5 * np.sum((y-f(x)) ** 2)
    return e 

# 4,初始化学习率η
ETA = 0.001

# 5,初始化误差变化量diff
diff = 1

# 6,初始化更新次数count
count = 0 

# 7,迭代学习
error = E(train_x_std, train_y)    # 计算初始误差
# 开始迭代更新参数，直到误差的变化小于0.01
while diff > 0.01:
    # 更新参数theta0和theta1
    theta0 = theta0 - ETA * np.sum((f(train_x_std) - train_y))
    theta1 = theta1 - ETA * np.sum((f(train_x_std) - train_y) * train_x_std)
    # 计算更新参数后的当前误差
    error_current = E(train_x_std, train_y)
    # 计算原误差与当前误差的差值
    diff = error - error_current
    # 更新误差值
    error = error_current
    # 输出日志
    count += 1
    print(f'第{count}次：theta0 = {theta0:.3f}, theta1 = {theta1:.3f}, 差值 = {diff:.4f}')

部分输出日志如下

注：若执行多次，会发现迭代次数和误差的差值在每次执行时都不一样，这是因为随机初始化的参数${\theta }_0$和${\theta }_1$?的不同而导致的。

此时，最终参数为：

print(f'theta0 = {theta0:.3f}')
print(f'theta1 = {theta1:.3f}')

作图查看拟合结果如下

plt.plot(train_x_std, train_y, 'o')     # 标准化后散点图
x = np.linspace(-3, 3, 100)
plt.plot(x, f(x))                       # 拟合直线
plt.show()

sklearn方法

在实际使用中，机器学习库sklearn为我们提供了更模块化的方式来进行线性回归，如下所示

from sklearn.linear_model import LinearRegression
# 1, 定义线性回归模型
lr = LinearRegression()
# 2, 训练模型
lr.fit(train_x_std.reshape(-1,1), train_y)

# 查看tehta0 和 theta1
print(f'theta0 = {lr.intercept_}')
print(f'theta1 = {lr.coef_}')

结果如下所示

和使用python逐步计算得到的参数一致。

np.polyfit(X, y, n)方法

Numpy库也提供了线性拟合方法：Numpy.polyfit(X, y, n)，使用方法如下所示：

# 一次线性拟合
z1 = np.polyfit(train_x_std, train_y, 1)
# 线性表示
p1 = np.poly1d(z1)
print(z1)
print(p1)

结果如下所示

和使用前两种方法得到的结果一致。

一元多次线性回归-多项式回归

观察前面做出的图像，发现一元一次线性模型$f_{\theta }(x)={\theta }_0+{\theta }_1?x$?的拟合效果并不好。

实际上，对于给定的数据，曲线比直线拟合地更好。因此，我们重新定义模型函数为：$$f_{\theta }(x)={\theta }_0+{\theta }_1?x+{\theta }_2?x^2$$?
该式为一元多次线性回归，即多项式拟合模型表达式。

如果使用更大次数的表达式，如下所示，就能表示更复杂的曲线了

$$f_{\theta }(x)={\theta }_0+{\theta }_1?x+{\theta }_2?x^2+{\theta }_3?x^3+?+{\theta }_n?x^n$$

次数越高，拟合的越好，但可能会出现过拟合问题

对于要解决的问题，在找出合适的表达式之前，需要不断地去尝试。

这里以二次拟合函数为例，我们增加了参数${\theta }_2$，此时：

1. 模型表达式为：$f_{\theta }(x)={\theta }_0+{\theta }_1?x+{\theta }_2?x^2$??
2. 目标函数为：$E(\theta )=\frac{1}{2}{\sum }_{i=1}^n{(y^{(i)}?f_{\theta }(x^{(i)}))}^2$?

使用同样的方式，对$E(\theta )$进行微分，得到参数更新表达式如下：

$$?\begin{array}{l}{\theta }_0:={\theta }_0?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})\\ {\theta }_1:={\theta }_1?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})?x^{(i)}\\ {\theta }_2:={\theta }_2?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})?x^{(i{)}^2}\end{array}$$

即使增加参数（增加次数），比如有${\theta }_3$??，${\theta }_4$?等，依然可以用同样的方法来求出参数更新表达式。

像这样增加函数中多项式的次数，然后再使用函数的分析方法被称为多项式回归

多元线性回归

前述的线性回归都是一元回归，即只有一个变量。

但是，实际中要解决的很多问题是变量超过2个的复杂问题。

例如：有三个变量，分别为$x_1$?，$x_2$?，$x_3$???。此时：

模型表达式为：$$f_{\theta }(x_1,x_2,x_3)={\theta }_0+{\theta }_1?x_1+{\theta }_2?x_2+{\theta }_3?x_3$$

将此式推广到n个变量的情况，此时

模型表达式为：$$f_{\theta }(x_1,\dots ,x_n)={\theta }_0+{\theta }_1?x_1+{\theta }_2?x_2+?+{\theta }_n?x_n$$

我们可以把参数$\theta$和变量$x$看作向量，用黑体表示

$$\theta =?\begin{array}{c}{\theta }_0\\ {\theta }_1\\ {\theta }_2\\ ?\\ {\theta }_n\end{array}?,x=?\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\\ x_2\end{array}\\ ?\\ x_n\end{array}?$$

其中$x_0=1$

则：${\theta }^{T}x={\theta }_0?x_0+{\theta }_1?x_1+{\theta }_2?x_2+?+{\theta }_n?x_n$?

因此：$$f_{\theta }(x)={\theta }^{T}x$$

对目标函数$E(\theta )=\frac{1}{2}{\sum }_{i=1}^n{(y^{(i)}?f_{\theta }(x^{(i)}))}^2$ 求 ${\theta }_j$ 的偏导如下：

$$\frac{?}{?{\theta }_j}E(\theta )=\frac{1}{2}\stackrel{n}{\sum _{i=1}}2?(y^{(i)}?f_{\theta }(x^{(i)}))?\frac{?(y^{(i)}?f_{\theta }(x^{(i)}))}{?{\theta }_j}$$??

$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?\frac{?(f_{\theta }(x^{(i)}))}{?{\theta }_j}$$????
$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?\frac{?({\theta }_0?x_0^{(i)}+{\theta }_1?x_1^{(i)}+{\theta }_2?x_2^{(i)}+?+{\theta }_n?x_n^{(i)})}{?{\theta }_j}$$?
$$=\stackrel{n}{\sum _{i=1}}(y^{(i)}?f_{\theta }(x^{(i)}))?(?1)?x_j^{(i)}$$
$$=\stackrel{n}{\sum _{i=1}}(f_{\theta }(x^{(i)})?y^{(i)})?x_j^{(i)}$$
综上：

1. 模型表达式为：

   $f_{\theta }(x)={\mathbf{\theta }}^{\mathbf{T}}\mathbf{x}={\theta }_0?x_0+{\theta }_1?x_1+{\theta }_2?x_2+?+{\theta }_n?x_n$

2. 目标函数为：$E(\theta )=\frac{1}{2}{\sum }_{i=1}^n{(y^{(i)}?f_{\theta }(x^{(i)}))}^2$

3. 参数更新表达式为：${\theta }_j:={\theta }_j?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})?x_j^{(i)}$

像这样包含了多个变量的回归称为多重回归。

采用矩阵求解参数

由于训练集数据有很多，所以我们把1行数据当作一个训练数据，以矩阵的形式来处理更方便。

n个变量($x_0$?????，$x_1$?????，$x_2$?????，…，$x_n$?????)，n个参数(${\theta }_0$?????，${\theta }_1$?????，${\theta }_2$?????，…，${\theta }_n$???????)，n个训练数据的模型表达式如下：

$$f_{\theta }(x)=X?\theta =\left[\begin{array}{cc}{\mathbf{x}}_0& {\mathbf{x}}_1{\mathbf{x}}_2?{\mathbf{x}}_{\mathbf{n}}\end{array}\right]??\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ?\\ {\theta }_n\end{array}?=?\begin{array}{cccc}{x}_0^{(1)}& x_1^{(1)}& ?& x_n^{(1)}\\ x_0^{(2)}& x_1^{(2)}& ?& x_n^{(2)}\\ x_0^{(3)}& x_1^{(3)}& ?& x_n^{(3)}\\ & & ?& \\ x_0^{(n)}& x_1^{(n)}& ?& x_n^{(n)}\end{array}???\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ?\\ {\theta }_n\end{array}?$$

? $$=?\begin{array}{c}{\theta }_0?x_0^{(1)}+{\theta }_1?x_1^{(1)}+?+{\theta }_n?x_n^{(1)}\\ {\theta }_0?x_0^{(2)}+{\theta }_1?x_1^{(2)}+?+{\theta }_n?x_n^{(2)}\\ {\theta }_0?x_0^{(3)}+{\theta }_1?x_1^{(3)}+?+{\theta }_n?x_n^{(3)}\\ ?\\ {\theta }_0?x_0^{(n)}+{\theta }_1?x_1^{(n)}+?+{\theta }_n?x_n^{(n)}\end{array}?$$

对于参数更新表达式：${\theta }_j:={\theta }_j?\eta ?{\sum }_{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})?x_j^{(i)}$，当$j=0$时，求和项可展开为：

$$\sum _{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})?x_0^{(i)}=(f_{\theta }(x^{(1)})?y^{(1)})?x_0^{(1)}+(f_{\theta }(x^{(2)})?y^{(2)})?x_0^{(2)}+?+(f_{\theta }(x^{(n)})?y^{(n)})?x_0^{(n)}$$

令：

$$\mathbf{f}=?\begin{array}{c}{f}_{\theta }(x^{(1)})?y^{(1)}\\ f_{\theta }(x^{(2)})?y^{(2)}\\ ?\\ f_{\theta }(x^{(n)})?y^{(n)}\end{array}?,{\mathbf{x}}_0=?\begin{array}{c}{x}_0^{(1)}\\ x_0^{(2)}\\ ?\\ x_0^{(n)}\end{array}?$$

则：

$$\sum _{i=1}^n(f_{\theta }(x^{(i)})?y^{(i)})?x_0^{(i)}={\mathbf{f}}^T?{\mathbf{x}}_0$$

Python逐步实现

对于多项式回归：

$$x_0=?\begin{array}{c}1\\ 1\\ ?\\ 1\end{array}?,x_1=?\begin{array}{c}{x}^{(i)}\\ x^{(i)}\\ ?\\ x^{(i)}\end{array}?,x_2=?\begin{array}{c}{x}^{(i{)}^2}\\ x^{(i{)}^2}\\ ?\\ x^{(i{)}^2}\end{array}?,X=[{\mathbf{x}}_0{x}_1{\mathbf{x}}_2]=?\begin{array}{ccc}1& x^{(i)}& x^{(i{)}^2}\\ 1& x^{(i)}& x^{(i{)}^2}\\ & ?& \\ 1& x^{(i)}& x^{(i{)}^2}\end{array}?$$

# 1，创建训练数据的矩阵X
def to_matrix(x):
    return np.vstack([np.ones(x.shape[0]), x, x**2]).T
X = to_matrix(train_x_std)

# 2，随机初始化参数Theta
theta = np.random.rand(3)

# 3，定义预测函数
def f(x):
    return np.dot(x, theta)

# 4, 初始化参数
diff = 1
ETA = 0.001

# 5, 迭代更新参数
error = E(X, train_y)
while diff > 0.01:
    theta = theta - ETA * np.dot(f(X)-train_y, X)
    current_error = E(X, train_y)
    diff = error - current_error
    error = current_error

运行结束后，得到参数如下

将结果绘图展示如下：

x = np.linspace(-3, 3, 100)

plt.plot(train_x_std, train_y, 'o') 	# 原数据散点图
plt.plot(x, f(to_matrix(x)))             # 多项式拟合曲线

plt.show()

np.polyfit(X, y, n)方法

z2 = np.polyfit(train_x_std, train_y, 2)
p2 = np.poly1d(z2)
print(z2)
print(p2)

结果如下所示

与python逐步实现结果一致