[数据结构与算法]【机器学习】DS的基础学习笔记1：线性回归

线性回归

1.1 专栏介绍

??本专栏学习材料为吴恩达《机器学习》课程，希望自学吴恩达的课程后通过自己的语言对算法进行描述且加强自己对模型的理解，后续在本科课程系统学习机器学习课程后，会对自学的笔记进行修改和补充，同时希望大神们及时批评与指正。

1.2 单变量线性回归

1.2.1 模型表示

课程中吴恩达老师以房屋交易问题为例，训练集如图所示：

我们希望通过这些数据来预测房价，我们需要对一些符号进行定义。

因为本例中的特征只有房屋大小，价格暂时只考虑与房屋大小有关，所以称此模型为单变量线性回归，而我们对数据的处理方法大概如下所示。

在此例中，我们该如何选择$h$呢？
因为该例只有一个特征，即房屋尺寸大小，所以我们可将假设函数设为
$$h_{\theta }(x)={\theta }_0+{\theta }_{1}x$$
$h$表示一个函数，输入的是房屋的尺寸大小，$h$ 根据输入的 $x$ 值来得出 $y$ 值，$y$ 值对应房屋出售价格，因此，$h$ 是一个从 $x$ 到 $y$ 的函数映射。

1.2.2 代价函数

??在设计假设函数 $h_{\theta }(x)={\theta }_0+{\theta }_{1}x$ 的过程也是直线拟合散点图的过程。确定假设函数也就是我们需要选择合适的 ${\theta }_0$ 与 ${\theta }_1$ ，也就是图像上直线的截距与斜率。我们选择的参数决定了我们得到的直线相对于我们的训练集的准确程度，模型预测值与实际值的差距也就是建模误差。

??我们的目标便是选择出可以使得建模误差的平方和能够最小的模型参数。在这里我们命名一种代价函数来表示模型误差。
$$minimizeJ({\theta }_0,{\theta }_1)=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})?y^{(i)}{)}^2$$
??代价函数又叫平方误差函数，解决回归问题，平方误差函数是一种常用手段。我们可以编程将参数${\theta }_0$ 、 ${\theta }_1$与代价函数$J$的关系可视化。

??可以看出，会有一组参数 $({\theta }_0,{\theta }_1)$ 使得 $J$ 处于最低点，也就是我们希望的误差最小的参数对。幸运的是，可以通过数学证明线性回归的代价函数总为凸函数，需要用到引理：二阶可微函数为严格凸函数的充要条件为黑塞矩阵为正定矩阵，因此我们总能在线性回归的代价函数找到唯一的极小值点。接下来我们需要探索一种可以自动找到最佳参数的方法。

1.2.3 梯度下降法

??梯度下降是一个用来求函数最小值的算法，我们会用该算法对本例假设函数进行最小化处理。该算法的思想大致为：首先随机选择一个参数组合$({\theta }_0,{\theta }_1,\dots ,{\theta }_n)$，计算代价函数，接着寻找下一个能让代价函数下降最多的参数组合，直到找到局部最低点为止。而我们不能确定该局部最低点是否为全局最低点，且不同的初始参数组合可能会有不同的局部最低点。

批量梯度下降(batch gradient descent)算法的公式为：
$$\begin{gathered} repeatuntilconvergence: \\ {\theta }_j={\theta }_j?\alpha ?\frac{?}{?{\theta }_j}J({\theta }_0,{\theta }_1)(forj=0andj=1) \end{gathered}$$
??其中 $\alpha$ 为学习率，它决定代价函数下降程度最大方向向下迈出的步子有多大。在批量梯度下降中，我们每一次迭代让所有参数减去学习率乘代价函数的导数，直至 $\theta$ 收敛。
注意： 我们需要同时更新所有参数 $\theta$

$temp0={\theta }_0?\alpha ?\frac{?}{?{\theta }_0}J({\theta }_0,{\theta }_1)$
$temp1={\theta }_1?\alpha ?\frac{?}{?{\theta }_1}J({\theta }_0,{\theta }_1)$
${\theta }_0=temp0$
${\theta }_1=temp1$

??我们对 $\alpha$ 有一定要求，若 $\alpha$ 过大，可能导致其无法收敛甚至发散，若 $\alpha$ 过小，则收敛速度会很慢。若已到达局部极值点，$\frac{?}{?\theta }J({\theta }_0,{\theta }_1)=0$，$\theta$不再改变。
??让我们把梯度下降法利用在线性回归上，则有如下公式：
$$\begin{gathered} repeatuntilconvergence: \\ {\theta }_0={\theta }_0?\alpha ?\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})?y^{(i)}) \\ {\theta }_1={\theta }_1?\alpha ?\frac{1}{m}\sum _{i=1}^m[(h_{\theta }(x^{(i)})?y^{(i)})?x^{(i)}] \\ update{\theta }_{0}and{\theta }_{1}simultaneously \end{gathered}$$

1.3 多变量线性回归

1.3.1 多维特征

??目前为止，我们探讨了单变量特征房屋大小对房价的影响，然而实际上有更多的特征会对房价造成影响，我们需要构建一个含有多个特征变量的模型，模型特征为$(x_1,x_2,x_3,\dots ,x_n)$,我们有新的训练集如下：

??同样地，我们需要对新的符号进行定义

??此时支持多变量的假设函数$h$为：
$$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\dots +{\theta }_n{x}_n$$
??这个公式有n+1个参数和n个变量，为了使公式更加简单，我们可以引入新变量 $x_0=1$ ，则公式转化为：
$$h_{\theta }(x)={\theta }_0{x}_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\dots +{\theta }_n{x}_n$$
??此时模型中的参数 $\theta$ 是一个n+1维向量，而训练实例 $x$ 也是n+1维向量，有如下表示：
$$\theta =?\begin{array}{c}{\theta }_0\\ {\theta }_1\\ {\theta }_2\\ \dots \\ {\theta }_n\end{array}?$$
$$X=?\begin{array}{c}{x}_0\\ x_1\\ x_2\\ \dots \\ x_m\end{array}?$$
??故假设函数 $h$ 可简化为：
$$h_{\theta }(x)={\theta }^{T}X$$
??其中上标T表示矩阵转置.

1.3.2 多元梯度下降

??与单变量线性回归类似，在多变量线性回归中，我们同样构建一个代价函数，则这个代价函数是所有建模误差的平方和，即：
$$minimizeJ({\theta }_0,{\theta }_1,\dots ,{\theta }_n)=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})?y^{(i)}{)}^2$$
其中 $h_{\theta }(x)={\theta }^{T}X={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\dots +{\theta }_n{x}_n$.
??我们的目标与单变量线性回归时一样，找到使得代价函数最小的参数。多变量线性回归的算法为：
$$\begin{gathered} repeatuntilconvergence: \\ {\theta }_j={\theta }_j?\alpha ?\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})?y^{(i)})?x^{(i)} \\ simultaneouslyupdate{\theta }_{j}forj=0,1,\dots ,n \end{gathered}$$
??注意上式的 $j$ 从0开始，因为我们引入了 $x_0=1$ ,使得算法对每一个参数 $\theta$ 都成立.

1.3.3 特征缩放

??在我们面对多维特征问题的时候，我们要保证这些特征都具有相近的尺度，这将帮助梯度下降算法更快地收敛。
??以房价问题为例，假设我们使用两个特征，房屋的尺寸和房间的数量，尺寸的值为 0-2000 平方英尺，而房间数量的值则是 0-5，两个特征绝对值取值差异非常的大，以两个参数分别为横纵坐标，绘制代价函数的等高线图能，看出图像会显得很扁，梯度下降算法需要非常多次的迭代才能收敛。
??我们解决的方法为尽量使得每个特征的取值范围为-1到1之间.

??更一般的情况，我们会利用标准差标准化公式：
$$x_n=\frac{x_n?{\mu }_n}{s_n}$$
??其中 ${\mu }_n$ 为平均值，$s_n$为标准差.

1.3.4 关于学习率

??我们不能提前预知梯度下降法的迭代次数，但可以通过可视化，找到代价函数 $J(\theta )$ 与 迭代次数的关系，帮助我们判断是否接近收敛。我们希望代价函数 $J(\theta )$ 在每一次迭代后都会下降。
??与此同时，有一种可自动判断收敛的方法，即如果代价函数 $J(\theta )$ 在一次迭代后下降少于$\epsilon$，则其收敛。这种方法的 $\epsilon$ 较小，比如$10^{?3}$，但因 $\epsilon$ 不好取值，因此可视化判断收敛程度更佳。
??在单变量线性回归中，我们讨论过学习率对于梯度下降收敛的影响，如果 $\alpha$ 过小，则达到收敛所需要的迭代次数会非常高；若 $\alpha$ 过大，每次迭代可能不会减小代价函数，可能会越过局部最小值导致无法收敛。
??通常可以尝试一些学习率：0.01，0.03，0.1，0.3，1，3……

1.3.5 特征与多项式回归

??有时候训练集给出的特征并不能很好地表现出与目标变量的影响，需要我们自己处理过后得到新的特征变量来描述关系。比如房价预测问题中，训练集给出房子的长度和宽度两个特征，不如直接用其乘积得到面积来体现与房价的关系，从而减少参数也更加直接。使得原本的假设函数$$h_{\theta }(x)={\theta }_0+{\theta }_1?width+{\theta }_2?length$$
变为

$$h_{\theta }(x)={\theta }_0+{\theta }_1?area$$
??有时线性回归并不适用于所有数据，有时更需要曲线来进行拟合，此时可以利用多项式回归这一思想，比如用二次函数来拟合。通常我们需要先可视化观察散点图后再判断准备用什么模型去尝试，之后对该模型进行转换，可化为线性模型。

??同样地，我们仍以房屋大小为特征变量，房价为目标变量为例。当我们想用三次函数模型时，可以令$x_0=1,x_1=x,x_2=x^2,x_3=x^3$,这样就又化为我们的线性回归模型。
??根据训练集所构建的散点图，我们还可以使：
$$h_{\theta }(x)={\theta }_0+{\theta }_1?size+{\theta }_2?\sqrt{size}$$
注意：在多项式回归时，我们利用平方关系、立方关系作为新的特征变量，此时变量取值差异十分明显，应更注重特征缩放这一环节。

1.3.6 正规方程

??我们之前在寻找最合适的参数 $\theta$ 时，所用的算法是梯度下降法，但在有的时候有一种叫做正规方程的方法可以不需要迭代更快地解出 $\theta$.
??正规方程法通过求解$\frac{?}{?\theta }J(\theta )=0$，得到合适的 $\theta$.假设我们的训练集特征矩阵为 $X$ (包含了$x_0=1$)，目标向量为 $y$，则利用正规方程解出向量
$$\theta =(X^{T}X{)}^{?1}{X}^{T}y$$
??我们对此进行以下简单推导：
$$\begin{gathered} 由已知可得，J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})?y^{(i)}{)}^2 \\ 其中，h_{\theta }(x)={\theta }^{T}X={\theta }_0{x}_0+{\theta }_1{x}_1+\dots +{\theta }_n{x}_n \end{gathered}$$
??将标量形式转化为矩阵形式，有
$$\begin{array}{rl}J(\theta )& =\frac{1}{2m}(X\theta ?y{)}^2\\ & =\frac{1}{2m}(X\theta ?y{)}^T(X\theta ?y)\\ & =\frac{1}{2m}({\theta }^T{X}^T?y^T)(X\theta ?y)\\ & =\frac{1}{2m}({\theta }^T{X}^{T}X\theta ?{\theta }^T{X}^{T}y?y^{T}X\theta +y^{T}y)\end{array}$$
??接下来对 $J(\theta )$ 求偏导，需要用到以下矩阵求导法则
$$\begin{array}{rl}\frac{dAB}{dB}& =A^T\\ \frac{d{B}^{T}A}{dB}& =A\\ \frac{d{X}^{T}AX}{dX}& =2AX\end{array}$$
??利用以上法则，继续对 $J(\theta )$ 求偏导,且令其为0
$$\begin{array}{rl}\frac{?J(\theta )}{?\theta }& =\frac{1}{2m}(2X^{T}X\theta ?X^{T}y?X^{T}y?0)\\ & =\frac{1}{m}(X^{T}X\theta ?X^{T}y)\\ & =0\end{array}$$
??解得
$$\theta =(X^{T}X{)}^{?1}{X}^{T}y$$
??关于矩阵求导公式及原理，可浏览其他博主的文章，本文便不再详述。可以看出，正规方程法对小数据集处理起来是非常快捷有效的。让我们将梯度下降与正规方程二者做个对比。

??除此之外，还需要注意一点：如果 $X^{T}X$ 不可逆，这种情况一般很少发生，那么有可能是因为训练集内有相似特征变量或者训练集内特征变量过多（比如$n>m$)，对于前者我们需要保证特征变量的线性无关性，后者可以删除一些特征或者使用正则化方法来找到参数。而且在Octave中的pinv()可以求解矩阵的伪逆，一样可以求解参数 $\theta$.

1.4 配套作业的Python实现

1.4.1 单变量线性回归

1. 问题背景：
??假如你是一个连锁餐厅的CEO，已有不同城市的人口与利润数据集(ex1data1.txt)，第一列为人口，第二列为利润，请你使用一个变量实施线性回归来预测连锁餐厅利润。
2. 代码实现：

??导入所需要的包

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

??导入数据集，一定要把数据文件ex1data1.txt放在和程序同一个文件夹里，否则需要使用绝对路径访问文件

path = 'ex1data1.txt'
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
print(data.head())  # 预览数据

print(data.describe())

??数据可视化，绘制散点图
??可利用matplotlib库将Dataframe转为数组后plt.scatter

data_x = data["Population"].values
data_y = data["Profit"].values
plt.figure(figsize=(12,8))
plt.scatter(data_x, data_y)
plt.xlabel("Population")
plt.ylabel("Profit")
plt.show()

??也可以直接用Pandas直接作图

data.plot(kind="scatter", x="Population", y="Profit", figsize=(12, 8))
plt.show()

??梯度下降法，以实现最小化代价函数。首先需要创建代价函数$$J({\theta }_0,{\theta }_1)=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})?y^{(i)}{)}^2$$

def compute_cost(X, y, theta):
    inner = np.power((X * theta.T)-y, 2)
    cost = np.sum(inner)/(2*len(X))
    return cost

??函数设计思路：因为$h_{\theta }(x)={\theta }_0+{\theta }_{1}x$，可利用矩阵形式批量计算误差后平方累和最后除以2倍样本数。这里的矩阵形式是因为theta的维度是一个（1，2）的行向量，而y是（97，1）的列向量，因此需要把theta取转置。
??由前文可知，引入$x_0=1$后将更方便利用矩阵处理代价函数与梯度，因此我们需要在DataFrame中引入该列。

data.insert(0, 'Ones', 1)

??创建好代价函数之后，我们需要处理DataFrame中的数据，使其可以直接传入compute_cost中，也就是得到所需的X，y，theta。

cols = data.shape[1]
X = data.iloc[:, 0:cols-1]  # X是所有行，去掉最后一列
y = data.iloc[:, cols-1:cols]  # y是所有行，最后一列
print(X.head())
print(y.head())

??我们利用切片处理数据时记得检查是否该数据是我们所需要的。

??检查完毕后，将这些Pandas数据类型转化为可传入函数的矩阵类型。

X = np.matrix(X.values)
y = np.matrix(y.values)
theta = np.matrix(np.array([0, 0]))

??转化完毕后可查看各矩阵的维度是否满足传入函数的条件。

print(X.shape, y.shape, theta.shape)

(97, 2) (97, 1) (1, 2)

??可以看出，果然与前文描述的一样，y为列向量，theta为行向量需要在函数内转置才可运算。接下来便可执行compute_cost函数，观察在
${\theta }_0=0,{\theta }_1=0$ 的条件下的代价函数值。

print(compute_cost(X, y, theta))

32.072733877455676

??接下来是批量梯度下降的代码实现，是本节最重要的知识之一。
$$\begin{gathered} repeatuntilconvergence: \\ {\theta }_j={\theta }_j?\alpha ?\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})?y^{(i)})?x^{(i)} \\ simultaneouslyupdate{\theta }_{j}forj=0,1,\dots ,n \end{gathered}$$

def gradientDescent(X, y, theta, alpha, iters):
    temp = np.matrix(np.zeros(theta.shape))
    parameters = int(theta.ravel().shape[1])
    cost = np.zeros(iters)

    for i in range(iters):
        error = (X * theta.T) - y

        for j in range(parameters):
            term = np.multiply(error, X[:, j])
            temp[0, j] = theta[0, j] - alpha / len(X) * np.sum(term)

        theta = temp
        cost[i] = compute_cost(X, y, theta)
    return theta, cost

??这里运用到的temp变量是为了保证每个 $\theta$ 同时更新，先把算出来的 ${\theta }_j$存起来，再把它们一起赋值给theta向量。细心的读者可以发现，这个批量梯度下降算法实际上可以直接处理多变量线性回归问题。接着，我们需要初始化学习率与迭代次数。

alpha = 0.01
iters = 1000

??最后让我们用梯度下降算法得到最新的参数 $\theta$.

new_theta, cost = gradientDescent(X, y, theta, alpha, iters)
print(new_theta)

[[-3.24140214 1.1272942 ]]

3. 可视化
??得到迭代出的参数 $\theta$ 后，便可以将线性模型可视化在散点图中观察拟合效果。

x = np.linspace(data.Population.min(), data.Population.max(), 100)
f = new_theta[0, 0] + (new_theta[0, 1] * x)

fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(x, f, 'r', label='Prediction')
ax.scatter(data.Population, data.Profit, label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
plt.show()

??有了结果之后，自然还可以做很多有意思的检验，做很多有趣的可视化，比如代价函数值随迭代次数的变化，观察学习率 $\alpha$ 是否设置合理，是否每一次迭代后都有所下降以及其收敛程度。

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

??还可以作出不同参数 $\theta$ 对代价函数的影响的三维图与等高线图。

theta_0 = np.linspace(-10, 10, 50)
theta_1 = np.linspace(-1, 4, 50)
Theta_0, Theta_1 = np.meshgrid(theta_0, theta_1)
Cost = np.zeros((50, 50))
for i in range(len(theta_0)):
    for j in range(len(theta_1)):
        Theta = np.matrix([theta_0[i], theta_1[j]])
        Cost[j][i] = compute_cost(X, y, Theta)
ax = plt.axes(projection='3d')
ax.plot_surface(Theta_0, Theta_1, Cost, cmap="viridis")
ax.set_xlabel("$\\theta_0$")
ax.set_ylabel("$\\theta_1$")
ax.set_zlabel("$J(\\theta)$")
plt.show()

theta_0 = np.linspace(-10, 10, 100)
theta_1 = np.linspace(-1, 4, 100)
[Theta_0, Theta_1] = np.meshgrid(theta_0, theta_1)
print(Theta_1)
Cost = np.zeros((100, 100))
for i in range(len(theta_0)):
    for j in range(len(theta_1)):
        Theta = np.matrix([theta_0[i], theta_1[j]])
        Cost[j, i] = compute_cost(X, y, Theta)
plt.contour(Theta_0, Theta_1, Cost, np.linspace(5, 20, 5))
plt.xlabel("$\\theta_0$")
plt.ylabel("$\\theta_1$")
plt.scatter(new_theta[0, 0], new_theta[0, 1], marker="x")
plt.show()

1.4.2 多变量线性回归

1. 问题背景

??你想要卖你的房子，你想知道什么是好的市场价格。我们收集了最近售出房屋的信息，实现一个多变量线性回归的模型。数据集为ex1data2.txt，里面包含了俄勒冈州波特兰市的房价信息，其中第一列为房子的大小（以平方英尺为单位），第二列为卧室数，第三列为房价。

2. 代码实现

??此问与第一问不同的是，该数据集的特征有两个，为多变量线性回归问题。我们依然选择使用批量梯度下降算法来解决该问题。首先仍是导入数据集ex1data2.txt.

path = 'ex1data2.txt'
data2 = pd.read_csv(path, header=None, names=['Size', 'Bedrooms', 'Price'])
print(data2.head())

??观察该数据集，我们很容易发现上文提到的特征缩放问题，卧室数与房价的范围差距太大了，因此我们需要对该数据集进行数据预处理——特征归一化。

data2 = (data2 - data2.mean()) / data2.std()
print(data2.head())

??现在我们重复第一问单变量线性回归的步骤，对数据进行预处理。

data2.insert(0, 'Ones', 1)
cols = data2.shape[1]
X2 = data2.iloc[:, 0:cols-1]
y2 = data2.iloc[:, cols-1:cols]
X2 = np.matrix(X2.values)
y2 = np.matrix(y2.values)
theta2 = np.matrix(np.array([0, 0, 0]))

??同样地应用批量梯度下降算法，得到合适的 $\theta$ 与 $J(\theta {)}_{min}$.

alpha = 0.01
iters = 1000
new_theta2, cost2 = gradientDescent(X2, y2, theta2, alpha, iters)
print(compute_cost(X2, y2, new_theta2))
print(new_theta2)

0.1307033696077189
[[-1.10910099e-16 8.78503652e-01 -4.69166570e-02]]

3. 可视化

??同样地，我们也可以将训练进程可视化。

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost2, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

1.4.3 正规方程

??对此问题我们依然使用第二问数据集ex1data2.txt.正规方程为
$$\theta =(X^{T}X{)}^{?1}{X}^{T}y$$

def normalEqn(X, y):
    theta = np.linalg.inv(X.T@X)@X.T@y
    return theta

final_theta = normalEqn(X2, y2)
print(final_theta.T)

[[-1.10910099e-16 8.78503652e-01 -4.69166570e-02]] # 梯度下降
[[-1.11022302e-16 8.84765988e-01 -5.31788197e-02]] # 正规方程

??可以发现，梯度下降的结果与正规方程有一定差距，原因可能是迭代次数不够等等。

1.4.4 线性回归的scikit-learn实现

??我们利用第一问的数据集ex1data1.txt，使用scikit-learn库实现。

from sklearn import linear_model
model = linear_model.LinearRegression()
model.fit(X, y)

x = np.array(X[:, 1].A1)
f = model.predict(X).flatten()

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x, f, 'r', label='Prediction')
ax.scatter(data.Population, data.Profit, label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
plt.show()

??以上，便是DS的machine learning学习笔记的第一节：线性回归的全部内容，由于是自学吴恩达课程，会有不严谨甚至错误的地方，欢迎大家指出。在本科系统学习《机器学习》课程后会对该节进行修改更新。下一节内容为逻辑回归。