简介

线性回归可以解决回归问题，同时也是一种思想简单，实现容易的算法。线性回归算法前提是可以预见所要解决的问题基本呈一个线性分布，比如房屋的面积和房屋的价格。这里使用到的Anaconda中的Jupyter工具编写代码，会使用到 numpy, sklearn两个Python库。

简单线性回归

简单的线性回归就是解决 $y = ax + b$ 给定x和y的值，求a和b的值。求解线性回归，如果所示

图中蓝色的点就是实际采集的数据点，红色的线就是需要得到的线性回归方程 $y = ax + b$ 的直线，就是要找到一条直线，所有的 $x$ 在直线上的点对应 $\hat{y}$ 与该点实际的 $y$ 的差值 $\Delta{y}$ 最小，即所有绿色的线对应的残差平方和 $\sum{i=1}^m\Delta{y}$ 最小。
这里不使用 $\sum{i=1}^m\Delta{y}$ 的原因在于 $\Delta{y}$ 的值有正负，求和之后正负抵消就不能体现出回归方程的结果的正确性，至于不使用绝对值的方式，是由于绝对值不方便数学求导，所以最后使用 $\sum_{i=1}^m(\Delta{y})^2$。

最小二乘法

数学推导，可以跳过
上面的说明最后得到的式子就是最小二乘法描述的公式，即求损失函数$J$的最小值：

$\sum_{i=1}^m(y_i - \hat{y_i})^2$ $\sum_{i=1}^m(y_i - \hat{y_i})^2 = \sum_{i=1}^m(y_i - ax_i - b)^2 = J(a,b)$

对它在 $a$，$b$ 两个方向上求偏导数，令偏导数等于0。由于 $b$ 求解较方便，对 $b$ 求偏导：

$\frac{\delta{J}}{\delta{b}} = 2\sum_{i=1}^m(y_i - ax_i - b)(-1) = 0$ $\sum_{i=1}^m{y_i} - a\sum_{i=1}^m{x_i} = \sum_{i=1}^m{b} = mb$

$a,b$ 是常量每项除以 m，得到 $b$的等式:

$b = \overline{y} - a\overline{x}$

现在对a求偏导：

$\frac{\delta{J}}{\delta{a}} = 2\sum_{i=1}^m(y_i - ax_i - b)(-x_i) = 0$

带入 $b = \overline{y} - a\overline{x}$ 移项化简得到

$a = \frac{\sum_{i=1}^m(x_iy_i - x_i\overline{y})}{\sum_{i=1}^m(x_i^2 - x_i\overline{x})}$

由于Python处理矩阵十分有优势，所以这里做个向量化的转换。由于

$\sum_{i=1}^mx_i\overline{y}=\overline{y}\sum_{i=1}^mx_i=\overline{y}*m\overline{x}$ $= \overline{x}\sum_{i=1}^my_i = \sum_{i=1}^my_i\overline{x} = \sum_{i=1}^m\overline{x}\overline{y}$

变换上面 $a$ 的等式

$a = \frac{\sum_{i=1}^m(x_iy_i - x_i\overline{y} - \overline{x}y_i + \overline{x}\overline{y})}{\sum_{i=1}^m(x_i^2 - x_i\overline{x} - x_i\overline{x} + \overline{x}^2)}$ $= \frac{\sum_{i=1}^m(x_i - \overline{x})(y_i - \overline{y})}{\sum_{i=1}^m(x_i - \overline{x})^2}$

代码实现

下面进行代码上的操作，这里用到了 numpy 这个库

1
2
3

# 模拟数据
  x = np.array([1., 2., 3., 4., 5.])
  y = np.array([3., 5.5, 4.7, 8., 9.])

$x,y$ 模拟数据在图上的显示：
初始化数据
现在根据公式

$a = \frac{\sum_{i=1}^m(x_i- \overline{x})(y_i - \overline{y})}{\sum_{i=1}^m(x_i - \overline{x})^2}$ $b = \overline{y} - a\overline{x}$

# 首先求出平均值，直接使用numpy的函数
  x_mean = np.mean(x)
  y_mean = np.mean(y)
  # 分别求和分子和分母，然后相除
  num = 0.0
  den = 0.0
  for i in range(len(x)):
      num += (x[i] - x_mean)*(y[i] - y_mean)
      den += (x[i] - x_mean)**2
  a = num / den
  b = y_mean - a * x_mean
# 这组数据求得 a = 1.45, b = 1.69

原始点和回归后的直线作图
回归结果
补充向量的方法，$a$ 的表达式分子分母都是两个元素相乘然后求和，和矩阵的点乘定义相同，所以求 $a$ 的代码可以修改为

1	a = (x - x_mean).dot(y - y_mean) / (x - x_mean).dot(x - x_mean)

同样可以算出 $a = 1.45$
完整代码

import numpy as np

x = np.array([1., 2., 3., 4., 5.])
x_mean = np.mean(x)
y_mean = np.mean(y)
# num = 0.0
# den = 0.0
# for i in range(len(x)):
#     num += (x[i] - x_mean)*(y[i] - y_mean)
#     den += (x[i] - x_mean)**2
# a = num / den
a = (x - x_mean).dot(y - y_mean) / (x - x_mean).dot(x - x_mean)
b = y_mean - a * x_mean

多元线性回归

最小二乘法

多元线性回归需要解决的问题比简单线性回归复杂一些，从解决 $y=ax + b$ 的问题变成了

$y = \theta_0 + \theta_1x_1 + \theta_2x_2 + \theta_3x_3 + \dots + \theta_nx_n$

要求 $n+1$ 个系数的值，将解看成一个向量，即求 $\theta_0,\theta_1,\dots\theta_n$的值，由于便于下面的计算，使用列向量

$\theta = （\theta_0, \theta_1, \theta_2\dots \theta_n)^T$

同样，我们需要求损失函数 $J$ 的最小值，$m$表示有多少组数据：

$J = \sum_{i=1}^m(y_i - \hat{y_i})^2$

由于 $\hat{y} = \theta_0 + \theta_1x_1 + \theta_2x_2 + \theta_3x_3 + \dots + \theta_mx_m$ 给每个元素添加一个可以添加一个 $x_0$，其中 $x_0 = 1$,可以获得$X_b$的矩阵：

$设X_b = \begin{pmatrix} 1 & x^{(1)}_1 & x^{(1)}_2 & \cdots & x^{(1)}_n \\ 1 & x^{(2)}_1 & x^{(2)}_2 & \cdots & x^{(2)}_n \\ \cdots \\ 1 & x^{(m)}_1 & x^{(m)}_2 & \cdots & x^{(m)}_n \\ \end{pmatrix}$ $\theta = \begin{pmatrix} \theta_0 \\ \theta_1 \\ \cdots\\ \theta_n \\ \end{pmatrix}$

可以变成矩阵点乘的形式

$\hat{y} = X_b\cdot{\theta}$

则 $\hat{y}$ 就是 $m$ 个样本的计算值，损失函数也可以修改成点乘的形式，由于要使用矩阵的点乘，进行一下转置：

$J = (y - X_b\cdot{\theta})^T\cdot(y - X_b\cdot{\theta})$

到这里，数值的运算转换成了矩阵的运算，这个矩阵的最小值通过最小二乘法的计算（矩阵求导之类的知识我也就不了解了，有兴趣的可以自己去看看），最终计算出的使得损失函数的结果最小的的$\theta$如下：

$\theta = (X_b^T\cdot{X_b})^{-1}\cdot{X_b^T}\cdot{y}$

直接通过公式求出最小解的方法时间复杂度为 $O(n^3)$, 优化之后 $O(n^{2.4})$，同时，现实生活中能够直接计算出公式的模型是比较少的，通常多元线性回归的计算会通过其他的方式来解决，比如梯度下降法

代码实现

与简单线性回归的流程相同，准备数据，这里使用sklearn里的Boston房价的数据集

from sklearn import datasets
boston = datasets.load_boston()
X = boston.data
y = boston.target

其中 X 是一个 $506*13$ 的一个矩阵，表示包含506条数据，每条数据有13个特征值，包括城镇人均犯罪率、一氧化氮浓度、是否靠近查尔斯河流、房屋数量等，y则是一个列向量，包含506个数据，表示506个房屋样本对应的房屋价格。

# 首先将X转化成X_b: 添加一列到第一列，所有元素都是 1
X_b = np.hstack((np.ones((len(X), 1)), X))
# 根据公式
theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

对于结果我做了处理，将结果 3.22611069e+01 转换成 32.3，方便直接看

1	theta = （32.3, -0.11, 0.04, -0.04, 0.45, -12.4, 3.75, -0.02, -1.21, 0.25, -0.01.3, -0.84, 0.0079, -0.35）

回归系数解读

使用线性回归算法计算的结果都是有可以解释性的，这是一个线性回归比较重要的特点。比如上面Boston房价的问题，从sklearn的数据集中了解到，样本数据 X 的每一个特征（列）代表的内容为：

1	names = ('CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT')

$\theta$ (theta)中第一维是常数项，不讨论，所以对应的系数值为

1	theta = (-0.11, 0.04, -0.04, 0.45, -12.4, 3.75, -0.02, -1.21, 0.25, -0.01.3, -0.84, 0.0079, -0.35）

抽取几个维度来查看，比如

‘CRIM’ （per capita crime rate by town, 城镇人均犯罪率）对应的系数为 -0.11，表示，随着犯罪率的升高，房价是呈下降的趋势；
‘CHAS’ (Charles River dummy variable (= 1 if tract bounds river; 0 otherwise))，是否傍河 Charles River，对应的系数为0.45，表示临河的房子贵
‘NOX’ ，一氧化氮浓度，是系数最小的项，表示浓度越高，房屋价格下降，同时也是下降的最厉害的
‘RM’，平均房间数，是系数最大的项，房间数量越多，房价越高，也很符合实际意义