线性回归

课程与数据集

本文对应 Andrew Ng（吴恩达）《Machine Learning》单变量 / 多变量线性回归部分。代码使用的 ex1data1.txt（城市人口 vs 连锁店利润）为课程官方编程练习数据。

数据与参考代码可在黄海广教授整理的中文笔记仓库下载：fengdu78/Coursera-ML-AndrewNg-Notes（见 code/ex1/ 目录）。课程主页：Coursera Machine Learning。

---

一、数学原理

1.1 模型表达

我们统计了某个连锁店在不同人口城市的年利润值，例如在 6.1101 万人的城市，年利润可以达到 17.5920 万人民币。现在我要在一个新的城市开这么一家连锁店，人口可以查到，希望预估一下我的年利润值。

这里直接使用吴恩达课程的数据表。你问 Fitten 它也是差不多的例子，你也可以直接去看 GitHub 上其他人的笔记。只不过我这里几乎不再使用 Numpy 中的 matrix 类型，而是依赖于 array 来实现线性代数的工作。下载 Watt Toolkit 打开 GitHub 加速是最直接的加速 GitHub 的方法。

我们先来回忆高中数学，好像确实学过这个东西，只不过今天我们要利用高等数学工具。

模型自然是一元线性方程 $y=wx+b$，我们要做的就是调节 $w$ 和 $b$，让它适合我们的数据点。

1.2 让线适合数据点

我们在高中天天用那个公式算（现在看来高中简直是我们最快乐的时光），现在我们要利用高等数学的知识找到我们的 $w$ 和 $b$。

首先我们来评估方程的误差，很显然用高中的方差就可以评估拟合的效果。

此处我们讨论更一般的情形：

$$f(\mathit{x})={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_n{x}_n$$

表示有 $n$ 个变量决定函数最终的值，${\theta }_0$ 为截距，$\mathit{x}$ 是一个 $n$ 维向量，里面是各个自变量。每一个数据对应的误差为 $f(x)-y$。

如果我们已经收集到了 $m$ 组数据，用上角标来标记它们，那么第 $i$ 组数据的误差为 $f({\mathit{x}}^{(i)})-y^{(i)}$。现在我们要反过来以向量

$$\mathit{\theta }=({\theta }_0{\theta }_1{\theta }_2\cdots {\theta }_n{)}^T$$

为自变量，讨论什么时候误差最小。我们定义代价函数：

$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(f(x^{(i)})-y^{(i)}{)}^2$$

这似乎就是方差除以了 2。虽然除以 2 并不影响误差的评估效果，但为什么不直接使用方差呢？这其实是因为我们要求导，平方项会有一个 2 放下来。

1.3 梯度下降

这是我们要利用的新知识。但让我们先看看梯度下降的原理：在纸面上画一个二次函数，任取一个点 $(x_0,y_0)$，计算它在这一个点的导数值，并作出其切线，取一个微小的步长 $\alpha >0$，让这个点按照

$$x_0^{\prime }=x_0-\alpha f^{\prime }(x_0)$$

移动。当它慢慢移动，一旦移动到最低点，导数值为 0，就会停止移动，这样就让误差停留在了最小值。

你会发现它总会朝着最低点方向移动，这就是梯度下降的原理。前提是步长 $\alpha$ 不能太大，一旦太大，就不是趋近于最小值了。我们把步长 $\alpha$ 称为学习率。

敏锐的同学肯定注意到，我们可以以任何一个点为起点进行梯度下降，但是我们取到的永远是局部最小值，也就是极小值。只不过对于那些能够很好拟合成直线的数据，其三维图都是碗状的，就像二维里面的二次函数，会趋于同一个最小值。

我们将式子 $x_0^{\prime }=x_0-\alpha f^{\prime }(x_0)$ 推广到一般形式，这里就直接用计算机赋值的表达方式。每一次操作称为一次迭代：

$${\theta }_j:={\theta }_j-\alpha \frac{\partial J(\theta )}{\partial {\theta }_j}$$

其中：

$$\frac{\partial J(\theta )}{\partial {\theta }_j}=\frac{1}{m}\sum _{i=1}^m(f(x^{(i)})-y^{(i)})x_j^{(i)}$$

这个线性函数求偏导当然不在话下。但这里有一个代码编写中需要注意的点：这个公式毫无疑问要用到原来的 $\theta$ 向量，所有的 ${\theta }_j$ 用的都是同一个 $\theta$，但是当你改变了 ${\theta }_0$ 后，原来的 $\theta$ 向量已经发生了改变，这就导致了计算不同步。所以在循环体内部要定义一个 temp 来储存原向量，用 temp 去进行运算。这个在代码编写过程中，会很容易 get 到。

---

二、代码编写

2.1 导入库与读取数据

首先我们按照惯例导入要用的库：

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sympy.abc import theta

在 GitHub、Gitee 或者课程上下载第一周的文件，请注意要将 ex1data1.txt 移动到代码文件的同一个文件夹里。

# 用 pandas 读取数据并保存在变量 data 中
path = 'ex1data1.txt'
# 文件并没有给列命名，所以 header=None，我们将数据命名为 Population, Profit
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
# 先来看一下散点图 (scatter)
data.plot(kind='scatter', x='Population', y='Profit', figsize=(12, 8))

2.2 代价函数

接下来就是编写代价函数了。在这个函数里，我们需要传入自变量数组 X、真实值数组 Y，以及我们要计算的 $\theta$ 数组。我们后续把 theta 定义为一维的行数组，所以矩阵乘法前要用到转置（.T 方法）。

def cost_function(X, Y, theta):
    inner = np.power((X @ theta.T) - Y, 2)  # 对应元素相乘，此处 inner 也是数组
    return np.sum(inner) / (2 * len(X))

2.3 准备训练数据

接下来我们把 txt 文件中的数据拆分，使之成为我们需要的几个数组。

# 现在最左端添加一列 1，这样就无需进行显式的加法操作表达截距
data.insert(0, 'ones', 1)
col = data.shape[1]  # 计算有多少列，shape[] 中 0 和 1 分别表示统计行数和统计列数
X = data.iloc[:, 0:col - 1]  # 用 iloc 方法取所有行，舍弃最后一列
Y = data.iloc[:, col - 1:col]  # 取最后一列
theta_begin = np.array([0, 0]).reshape(1, 2)  # 此段代码已经放弃 matrix，所以注意用 reshape 匹配原矩阵形状

# 检查数组形状
print(X.shape)
print(theta_begin.shape)
print(Y.shape)

2.4 梯度下降

下面我们来实现梯度下降函数。我们不可能让偏导数为 0 的时候才停下来，计算机很难算出标准的 0，这里我们指定一个迭代次数 iters。显然迭代 iters 次有一个循环，我们要改变 $\theta$ 数组里所有的 $n+1$ 个元素又需要一个循环。

def gradientDescent(x, y, theta, alpha, iters):
    temp = np.zeros(theta.shape)  # 我们不能在第二个循环里直接改变 theta 的值，创建一个 temp 储存迭代后的数据
    parameters = theta.ravel().shape[0]  # .ravel 平摊数组统计 theta 参数个数，这就是第二个循环的循环次数
    # 我们需要第一个循环，来进行 iters 次迭代
    for i in range(iters):
        error = (x @ theta.T) - y  # 每一次的误差都是不一样的，误差写在第二个循环前
        for j in range(parameters):
            term = np.multiply(error, x[:, j].reshape(-1, 1))  # 第二个循环每一次都在改变角标为 j 的 theta 的值
            temp[0, j] = theta[0, j] - (alpha / len(x)) * np.sum(term)  # 这里 temp 就发挥了存储作用
        theta = temp.copy()  # 跳出第二个循环后将 theta 值统一改变
    return theta  # 返回 theta 数组

alpha = 0.01
iters = 1000
g = gradientDescent(X, Y, theta_begin, alpha, iters)
print(g)

2.5 可视化拟合结果

现在我们来可视化数据：

x = np.linspace(data.Population.min(), data.Population.max(), 100)  # 横坐标
f = g[0, 0] + (g[0, 1] * x)  # 拟合的直线

fig, ax = plt.subplots(figsize=(12, 8))  # 创建图形对象 fig，轴对象 ax
ax.plot(x, f, 'r', label='Prediction')  # 直线
ax.scatter(data.Population, data.Profit, label='Training Data')  # 散点
ax.legend(loc=2)  # 添加图例，loc=2 左上方
ax.set_xlabel('Population')  # 横轴标签
ax.set_ylabel('Profit')  # 纵轴标签
ax.set_title('Predicted Profit vs. Population Size')  # 标题
plt.show()  # 显示

2.6 记录并可视化代价函数变化

我们想看看代价函数的变化怎么弄呢？

# 在全局定义一个 cost 变量，用以记录代价函数的值，注意放在 iters = 1000 的后面
cost = np.zeros(iters)

将梯度下降函数修改为：

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

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

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

        theta = temp.copy()
        cost[i] = cost_function(X, Y, theta)

    return theta

接下来我们可视化代价函数图像：

fig2, ax2 = plt.subplots(figsize=(12, 8))
ax2.plot(np.arange(iters), cost, 'r')
ax2.set_xlabel('Iterations')
ax2.set_ylabel('Cost')
ax2.set_title('costchange')
plt.show()

---

下一篇：多变量情形下的特征缩放见《特征缩放》笔记（数据集 ex1data2.txt，同样来自上述仓库）。