线性回归

1 线性回归理论

1.1 问题描述

对于只有一个自变量x和一个因变量y的数据集，我们假设(x,y)之间的关系可以用

$$h(x)=\theta x + b$$ 进行映射，那么我们就把求解这个函数的过程称为线性回归求解。可以在二维坐标系拟合出一条直线来代表(x,y)的趋势。很显然，不可能所有的点刚好在这条直线上，所以求解最优的标准就是，使点集(x,y)与拟合函数间的误差（平方误差）最小。

推广到多个自变量，假设一个数据集有n个自变量，也可以说成n个属性，$\textbf{x}=(x_1;x_2;...;x_n)$，如果使用线性模型去拟合则有

$$h_{\theta}(\textbf{x})=\theta_1 x_1 +\theta_2 x_2 +...+\theta_n x_n +b$$ 这种方式称为多元线性回归（multivariate linear regression）.需要求解的未知量有 $\theta_1, \theta_2,...,\theta_n,b$. 为了便于讨论，令 $b=\theta_0 x_0$，其中 $x_0=1$，那么，原来的目标函数就可以表示为 $$h_{\theta}(\textbf{x}_i)=\theta_0 x_0 + \theta_1 x_1 +\theta_2 x_2 +...+\theta_n x_n$$ 用矩阵的方式表示为 $$h_{\theta}(\textbf{x}_i)=\boldsymbol{\theta} ^T\textbf{x}_i, i=0,1,...,n$$ 即 $$\textbf{y}=\boldsymbol{\theta} ^T\textbf{X}$$ 其中 $\textbf{X}$ 是一个m*(n+1)的矩阵，前n维数据是数据集的数据，最后一维值为1，代表 $x_0$ $$\textbf{X}= \begin{pmatrix} x_{11} & x_{12} & \cdots & x_{1n} & 1 \\ x_{21} & x_{22} & \cdots & x_{2n} & 1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} & 1 \\ \end{pmatrix}$$ 所以说，求解线性模型的解，就是求解 $\boldsymbol{\theta}$。当 $\sum_{i=1}^m (h_{\theta}(x^{(i)}) - y^{(i)})^2$最小，即平方误差最小时，有最优解。 求解方法有梯度下降法（Gradient Descent）、Normal Equation等等。梯度下降有如下特点：需要预先选定步长a、需要多次迭代、特征值需要Scaling（统一到同一个尺度范围）。因此比较复杂，还有一种不需要迭代的求解方式-Normal Equation，简单、方便、不需要Feature Scaling。但是Normal Equation最小二乘法需要计算X的转置与逆矩阵，计算量很大，因此特征个数多时计算会很慢，只适用于特征个数小于100000时使用；当特征数量大于100000时使用梯度法。另外，当X不可逆时就有岭回归算法的用武之地了。

1.2 Normal Equation

Normal Equation最常用的就是最小二乘法。

对$\theta$矩阵求导，得到其最优解为

$$\hat{\boldsymbol{\theta}}^* = (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T \textbf{y}$$

1.3 Gradient Descent

梯度下降法（Gradient Descent）需要根据平方误差，定义该线性回归模型的损失函数（Cost Function）：

$$J(\theta)=J(\theta_0,\theta_1,...,\theta_n)=\frac{1}{2m} \sum_{i=1}^m (h_{\theta}(x^{(i)}) - y^{(i)})^2$$

线性回归的损耗函数的值与回归系数θ的关系是碗状的，只有一个最小点。线性回归的求解过程如同Logistic回归，区别在于学习模型函数hθ(x)不同.

1.4 局部加权线性回归

线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有最小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得最好的预测效果。所以有些方法允许在估计中引人一些偏差，从而降低预测的均方误差。其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression, LWLR )。在该算法中，我们给待预测点附近的每个点赋予一定的权重.于是公式变为:

$$\hat{\boldsymbol{\theta}} = (\textbf{X}^T \textbf{W} \textbf{X})^{-1} \textbf{X}^T \textbf{W} \textbf{y}$$
W是（m,m）矩阵，m表示样本数。

LWLR使用 “核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：

$$w(i,i)=exp(\,\frac{|x^{(i)}-x|}{-2k^2})\,$$
k需要优化选择.

局部加权线性回归也存在一个问题，即增加了计算量，因为它对每个点做预测时都必须使用整个数据集，而不是计算出回归系数得到回归方程后代入计算即可。因此该算法不被推荐。

2 简单的Python实例

根据线性回归的理论的Normal Equation方法，使用Python进行简单的求解，并把结果可视化。

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline

# 使用Pandas读取数据，并查看几条样例数据。
# 第一列全是0，根据上面的理论知识，这是x0的占位符，不要太在意。我们需要拟合的是x1与y之间的线性关系。
ex0 = pd.read_csv("ex0.txt", sep="\t", names=["x0", "x1", "y"])
ex0.head()

	x0	x1	y
0	1.0	0.067732	3.176513
1	1.0	0.427810	3.816464
2	1.0	0.995731	4.550095
3	1.0	0.738336	4.256571
4	1.0	0.981083	4.560815

2.1 整体拟合效果预览

# seaborn可以直接拟合出直线，我们先看一下最终的可视化效果。
# 下图的圆点就是原始数据(x1,y),直线就是数据集的回归线
sns.lmplot(x = 'x1', y='y', data=ex0)

<seaborn.axisgrid.FacetGrid at 0x7fa34ec26438>

2.2 逐步求解线性回归参数

下面根据理论知识，自己动手，逐步求解回归线的参数$\boldsymbol{\theta}$.

根据上面推导的多元线性回归最优解公式

$$\hat{\boldsymbol{\theta}}^* = (\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T \textbf{y}$$

其中矩阵$\textbf{X}$和矩阵$\textbf{y}$可从刚刚读取的ex0变量中获得。

# 获取矩阵X和矩阵y
# 为了使用NumPy中矩阵计算的方法，将Pandas的DataFrame对象转换成NumPy中的ndarray对象
# 如果读取的数据文件没有这里的x0这个数据，需要手动添加一列全为1的数据到xMat对象中
xMat = ex0.as_matrix(['x1', 'x0'])
yMat = ex0.as_matrix(['y'])

# 计算XTX，注意矩阵乘法是np.dot()，而不是np.mul()
xTx = np.dot(xMat.T, xMat)
# 由于有矩阵的逆的计算，所以需要判断行列式不为零
np.linalg.det(xTx) == 0

False

# 完成最后的theta值的计算
ws = np.dot(np.linalg.inv(xTx), xMat.T)
ws = np.dot(np.linalg.inv(xTx), np.dot(xMat.T, yMat))
ws

array([[ 1.69532264],
       [ 3.00774324]])

于是，我们得到解

$$f(x)=1.69532264 x_1 + 3.00774324$$ 注意，这里ws变量中哪个是x0的系数，根据矩阵运算的过程，这个系数在ws中的位置和x0在xMat变量的位置一致，即最后一位。

根据这个解，画出$\textbf{X}$对应的最佳拟合直线图

# 直线上的y值，即预测值
yHat = np.dot(xMat, ws)
# 画原始值的散点图，即(x1,y)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(ex0['x1'], yMat)
# 画拟合的直线
ax.plot(xMat[:,0], yHat)
plt.show()

2.3 评估拟合效果

可以发现，不论什么样的数据集，都可以使用线性回归来拟合，那么拟合的效果如何呢。可以使用皮尔森相关系数（Pearson correlation coefficient）也称皮尔森积矩相关系数(Pearson product-moment correlation coefficient)来衡量，值在[-1,1],表示数据的线性相关程度。

# 计算y的估计值和实际值的相关系数，可以反应这套数据使用线性回归模型拟合的效果。
# 两个向量都必须是行向量
np.corrcoef(yHat.T, yMat.T)

array([[ 1.        ,  0.98647356],
       [ 0.98647356,  1.        ]])

从结果可以看出，变量和自己的相关度都是1，表示完全呈线性相关，这是显然的。估计值和实际值线性相关度为0.98，说明(x1,y)呈线性正相关，线性拟合效果不错。

2.4 线性回归工具

在开发中，我们往往不必自己从零实现线性回归求解过程，可以使用SciPy中的linear_model包，或者StatsModels包，下面以Scipy为例。

from sklearn.linear_model import LinearRegression
regr = LinearRegression()
# 线性拟合，fix(x,y)，xy都必须是列向量
regr.fit(xMat[:, 0].reshape(200,1), yMat)
predictions={}
# 根据拟合的模型计算预测值
predictions['value'] = regr.predict(xMat[:, 0].reshape(200,1))
predictions['coefficient'] = regr.coef_
predictions['intercept'] = regr.intercept_
predictions

{'coefficient': array([[ 1.69532264]]),
 'intercept': array([ 3.00774324]),
 'value': array([[ 3.12257084],
        [ 3.73301922],
        [ 4.69582855],
        ...
        [ 3.20467701]])}

可以看出coefficient就是上面我们自己求解的系数x1，intercept就是$\theta_0$，value是使用求解的函数预测的y值。结果和上面自己实现的方法一致。