LionKing数据科学专栏

购买普通会员高级会员可以解锁网站精华内容且享受VIP服务的优惠

想要查看更多数据科学相关的内容请关注我们的微信公众号知乎专栏

线性回归(linear regression)原理及算法实现

线性回归的原理

假设你有一个数据集,希望通过一些特征$x_1, \ldots, x_p$预测目标变量$y$。最简单的模型是假设目标变量$y$是这些特征的某个线性组合:

$$y \approx a_0 + a_1x_1 + \ldots + a_px_p$$

记第$i$组观测为$(x_1^{(i)}, \ldots, x_n^{(i)}, y^{(i)})$。总共有$n$组观测。

第$i$组观测的预测值为$\hat{y}^{(i)} = a_0 + a_1x_1^{(i)} + \ldots + a_px_p^{(i)}$

我们将$a_0, a_1, \ldots, a_p$视为参数,最小化均方误差(Mean Squared Error):

$$L(a_0, a_1, \ldots, a_p) = \frac{1}{n}[(y^{(1)} - \hat{y}^{(1)})^2 + \ldots + (y^{(n)} - \hat{y}^{(n)})^2]$$

算法

假设$x_0 = 1, x_1, \ldots, x_p$组成的数据矩阵$X_{n \times (p + 1)}$列满秩(full column rank),即秩为$p + 1$。则$X^TX$是秩为$p + 1$的方阵,因此可逆。并且最优的参数$(a_0, a_1, \ldots, a_p)$满足正规方程(normal equation):

$$X^TX\begin{bmatrix}a_0\\a_1\\\ldots\\a_p\end{bmatrix} = X^Ty$$

需要购买普通会员高级会员登录后刷新该页面查看

因为$X^TX$可逆,可以通过求解该线性系统得到待求参数。

当$p$过大时,求解线性系统的效率很低,因此我们转而使用梯度下降(gradient descent)近似地估计参数。

Python实现

Python中的scikit-learn包可以进行线性回归模型的训练:scikit-learn线性回归模型文档

import numpy as np
from sklearn import linear_model

n_train = 1000
n_test = 200
p = 10

# 真实参数
intercept_true = np.random.rand()
coef_true = np.random.rand(p, 1)

# 生成X
X = np.random.rand(n_train + n_test, p)
# 生成y = intercept + X * coef + noise
y = intercept_true + X.dot(coef_true) + np.random.randn(n_train + n_test, 1)

# 拆分训练和测试数据集
X_train = X[:n_train, ]
y_train = y[:n_train, ]
X_test = X[n_train:, ]
y_test = y[n_train:, ]

# 创建线性回归训练器
linear_regressor = linear_model.LinearRegression()
# 训练模型
linear_regressor.fit(X_train, y_train)
# 训练集预测
yhat_train = linear_regressor.predict(X_train)
# 训练均方误差
residual_train = yhat_train - y_train
mse_train = np.mean(residual_train ** 2)
print('训练均方误差为%f' % (mse_train, ))
# 测试集预测
yhat_test = linear_regressor.predict(X_test)
# 测试均方误差
residual_test = yhat_test - y_test
mse_test = np.mean(residual_test ** 2)
print('测试均方误差为%f' % (mse_test, ))

# 模型参数
intercept_fit = linear_regressor.intercept_
coef_fit = linear_regressor.coef_
# 模型参数均方误差
mse_param = (np.sum((coef_fit.flatten() - coef_true.flatten()) ** 2) + (intercept_fit[0] - intercept_true) ** 2) / (p + 1)
print('模型参数均方误差为%f' % (mse_param, ))
              
输出:

训练均方误差为0.984766
测试均方误差为0.943605
模型参数均方误差为0.014939
	      

共线性(collinearity)

当数据矩阵$X_{n \times (p + 1)}$不是列满秩时,正规方程中的$X^TX$不可逆,任何满足正规方程的参数都能够最小化均方误差。满足要求的解有无限多个,因此得出的解容易过拟合(overfitting)。这个问题又称数据共线性。

如果$n \leqslant p$,则数据一定共线性。否则,数据很少精确地满足共线性。即使数据不是精确地满足共线性,也可能会非常接近列不满秩。此时$X^TX$虽然可逆,但是数值求解其逆矩阵(inverse matrix)非常不稳定,这样得到的不稳定解也容易产生过拟合。

解决共线性的方法主要有两种。第一种是找出共线性的列并且删除其中一列直到数据列满秩。第二种是使用正则化(regularization)。

岭回归(ridge regression),Lasso回归和弹性网络(Elastic net)回归

对于线性回归,正则化的两种方法是使用岭回归(ridge regression)或Lasso。

岭回归的思想是在原来的均方误差损失函数的基础上加上参数的$\ell_2$范数平方作为惩罚项。新的损失函数是:

$$L_{\lambda}(a_0, a_1, \ldots, a_p) = L(a_0, a_1, \ldots, a_p) + \frac{\lambda}{n}[a_0^2 + a_1^2 + \ldots + a_p^2]$$

类似原来的正规方程,新的正规方程为

$$(X^TX + \lambda I_{p + 1})\begin{bmatrix}a_0\\a_1\\\vdots\\a_p\end{bmatrix} = X^Ty$$

需要购买高级会员登录后刷新该页面查看

由于$X^TX$一定是半正定(positive semi-definite)的,只要$\lambda \gt 0$, $X^TX + \lambda I_{p + 1}$的每一个特征值都不小于$\lambda$,因此$X^TX + \lambda I_{p + 1}$是正定(positive definite)的,进而可逆。岭回归的解存在且唯一。

Lasso的思想是在均方损失误差的基础上加上参数的$\ell_1$范数作为惩罚项。新的损失函数是:

$$L_{\lambda}(a_0, a_1, \ldots, a_p) = L(a_0, a_1, \ldots, a_p) + \frac{\lambda}{n}[|a_0| + |a_1| + \ldots + |a_p|]$$

Lasso的求解相对于岭回归复杂一些,主流的算法是Least Angle Regression(LARS)。

Lasso倾向于得出稀疏(sparse)解,岭回归倾向于得出稠密(dense)解。

弹性网络(Elastic net)结合了两者,同时使用$\ell_2$范数和$\ell_1$范数作为惩罚项,其形式如下:

$$L_{\lambda}(a_0, a_1, \ldots, a_p) = L(a_0, a_1, \ldots, a_p) + \frac{\lambda}{n}(\tau\|\beta\|_2^2 + (1 - \tau)\|\beta\|_1)$$

其中$0 \leqslant \tau \leqslant 1$调节两种惩罚的比例,若$\tau = 0$,则弹性网络回归化为Lasso回归;若$\tau = 1$,则弹性网络回归化为岭回归。

弹性网络回归的好处是即保留了稀疏性,又可以把相关性高的变量同时找出来。

均方误差的概率意义

线性回归采用最小化残差平方和主要有三个原因:第一,数学上残差平方和的解,在数据列满秩的前提下,存在且唯一;第二,平方和的解有显式数学解可以精确求得;第三,最小化残差平方和等价于正态(Gaussian)假设下的最大似然估计(Maximum likelihood estimation)。

对于第三点,假设数据来自$y = X\beta + \varepsilon$,其中$\varepsilon \sim N(0, \sigma^2)$。

给定数据$(X^{(1)}, y^{(1)}), \ldots, (X^{(n)}, y^{(n)})$,第$i$组数据的似然即$y^{(i)} - X^{(i)}\beta$在正态分布中的概率密度:

$$L_i = \frac{1}{\sqrt{2\pi}\sigma}\exp{\left[-\frac{(y^{(i)} - X^{(i)}\beta)^2}{2\sigma^2}\right]}$$

总体的似然为$L = L_1L_2\ldots L_n$。最大化似然即最大化每个数据的对数似然(log likelihood)之和:

$$l = \log{L_1} + \ldots + \log{L_n}$$

$$\log{L_i} = -\log{\left(\sqrt{2\pi}\sigma\right)} - \frac{1}{2\sigma^2}(y^{(i)} - X^{(i)}\beta)^2$$

$$l = -n\log{\left(\sqrt{2\pi}\sigma\right)} - \frac{1}{2\sigma^2}[(y^{(1)} - \hat{y}^{(1)})^2 + \ldots + (y^{(n)} - \hat{y}^{(n)})^2]$$

因此最大化$l$等价于最小化$[(y^{(1)} - \hat{y}^{(1)})^2 + \ldots + (y^{(n)} - \hat{y}^{(n)})^2]$,即残差平方和。

岭回归的概率意义

岭回归从概率意义上可以看作在数据服从高斯噪声假设下,参数的先验分布(prior distribution)为以0为中心的正态分布。

若假设参数$\beta$服从先验分布$N(\vec{0}, \tau^2I_{p + 1})$,则$\beta$的先验概率为

$$prior(\beta) = C\exp{-\frac{\|\beta\|_2^2}{2\tau^2}}$$

其中$C$为只与$p$和$\tau$有关的常数。

根据贝叶斯定理(Bayes' theorem),数据的后验似然(posterior likelihood)与先验和$L$的乘积成正比。

$$posterior(\beta; X, y) \propto prior(\beta)L$$

$$\log{posterior(\beta; X, y)} = C - \frac{1}{2\tau^2}\|\beta\|_2^2 - \frac{1}{2\sigma^2}[(y^{(1)} - \hat{y}^{(1)})^2 + \ldots + (y^{(n)} - \hat{y}^{(n)})^2]$$

最大化后验似然等价于最小化

$$[(y^{(1)} - \hat{y}^{(1)})^2 + \ldots + (y^{(n)} - \hat{y}^{(n)})^2] + \frac{\sigma^2}{\tau^2}\|\beta\|_2^2$$

取$\lambda = \frac{\sigma^2}{\tau^2}$,最大化后验似然等价于岭回归的损失函数。

类似地,Lasso可以视作给参数一个拉普拉斯分布(Laplace distribution)作为先验。

常见面试问题

Q:为什么需要正则化?

需要购买普通会员高级会员登录后刷新该页面查看

Q:为什么$\ell_1$更容易给出稀疏解(sparse solution)?

需要购买普通会员高级会员登录后刷新该页面查看

Q:$R^2$是什么?应当如何计算?

需要购买普通会员高级会员登录后刷新该页面查看

Q:共线性应如何检测?

需要购买普通会员高级会员登录后刷新该页面查看

更多机器学习相关问题见本网站论坛机器学习理论版面机器学习实践版面

更多面试问题见面试真题汇总

想要查看更多数据科学相关的内容请关注我们的微信公众号知乎专栏