【预测算法】01. 线性回归原理及R语言实现

回归分析是统计学的核心算法，是机器学习最基本算法，也是数学建模最常用的算法之一。

简单来说，回归分析就是用一个或多个自变量来预测因变量的方法，具体是通过多组自变量和因变量的样本数据，拟合出最佳的函数关系。

本篇由前入深将线性回归的原理讲清楚，并用案例演示实际操作。

一、最小二乘法

设有组样本点： $(x_i, y_i), \quad i=1, \cdots, n$

例1，现有10期的广告费用与销售额的数据：

先画散点图观察一下：

cost<-c(30,40,40,50,60,70,70,70,80,90)
sale<-c(143.5,192.2,204.7,266,318.2,457,333.8,312.1,386.4,503.9)
dat<-as.data.frame(cbind(cost,sale))
plot(dat)

可见，这些散点大致在一条直线上，一元线性回归就是寻找一条直线，使得与这些散点拟合程度最好（越接近直线越好）。

比如画这样一条直线，方程可写为： $y=\beta_0 + \beta_1 x$ （线性模型），其中 $\beta_0, \beta_1$ 是待定系数，目标是选取与样本点最接近的直线对应的 $\beta_0, \beta_1$ .

那么，怎么刻画这种“最接近”？

$\hat{y}_i = \beta_0 + \beta_1 x_i$ 是与横轴对应的直线上的点的纵坐标（称为线性模型预测值），它与样本点对应的真实值之差，就是预测误差（红线长度）：

$\varepsilon_i = |y_i - \hat{y}_i|, \quad i = 1, \cdots, n$

适合描述散点到直线的“接近程度”。

但绝对值不容易计算，改用：

$\varepsilon_i^2 = (y_i - \hat{y}_i)^2, \quad i = 1, \cdots, n$

我们需要让所有散点总体上最接近该直线，故需要让总的预测误差

$J(\beta_0, \beta_1) = \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{i=1}^n [y_i - (\beta_0+\beta_1 x_i)]^2$

最小。

于是问题转化为优化问题，选取 $\beta_0, \beta_1$ 使得

$\min \, J(\beta_0, \beta_1) = \sum_{i=1}^n [y_i - (\beta_0+\beta_1 x_i)]^2$ (1)

这就是“最小二乘法”，有着很直观的几何解释。

二、问题(1)求解

这是个求二元函数极小值问题。

根据微积分知识，二元函数极值是在一阶偏导等于0点处取到：

$\left\{ \begin{array}{l} \dfrac{\partial J}{\partial \beta_0} = -2 \sum\limits_{i=1}^n \big[ y_i - \beta_0 - \beta_1 x_i ) \big] =0 \\[0.1cm] \dfrac{\partial J}{\partial \beta_1} = -2 \sum\limits_{i=1}^n \big[ y_i - \beta_0 - \beta_1 x_i ) \big] x_i =0 \end{array} \right.$

解关于 $\beta_0, \beta_1$ 的二元一次方程组得

$\left\{ \begin{array}{l} \beta_0 = \bar{y} - \beta_1 \bar{x} \\[0.2cm] \beta_1 = \dfrac{\sum\limits_{i=1}^n x_i y_i - \bar{y} \sum_\limits{i=1}^n x_i}{\sum\limits_{i=1}^n x^2_i - \bar{x} \sum\limits_{i=1}^n x_i} = \dfrac{\sum\limits_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum\limits_{i=1}^n (x_i - \bar{x})^2} \end{array} \right.$

其中，

$\left\{ \begin{array}{l} \bar{x} = \dfrac{1}{n} \sum\limits_{i=1}^n x_i \\[0.2cm] \bar{y} = \dfrac{1}{n} \sum\limits_{i=1}^n y_i \end{array} \right.$

三、提升：矩阵形式推导

将线性模型的全部预测值，用矩阵来写：

$\begin{bmatrix} 1 & x_1 \\ \vdots & \vdots \\ 1 & x_n \end{bmatrix}_{n \times 2} \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}_{2 \times 1} =\begin{bmatrix} \hat{y}_1 \\ \vdots \\ \hat{y}_n \end{bmatrix}_{n \times 1}$

记

$X = \begin{bmatrix} 1 & x_1 \\ \vdots & \vdots \\ 1 & x_n \end{bmatrix}_{n \times 2} , \quad \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}_{2 \times 1} , \quad \hat{Y} =\begin{bmatrix} \hat{y}_1 \\ \vdots \\ \hat{y}_n \end{bmatrix}_{n \times 1}$

则矩阵表示为

$\hat{Y}= X \beta$

于是，让预测误差最小的“最小二乘法”优化问题就表示为

$\min\limits_{\beta} J(\beta) = \|Y - \hat{Y} \|^2 = \|Y - X \beta \|^2$ (2)

这里 $\| \cdot\|$ 即向量的长度，计算一下：

$\begin{eqnarray*} \|Y - X \beta \|^2 &=& \langle Y-X \beta, Y-X \beta \rangle \\ &=& (Y - X \beta)^T (Y - X \beta) \\ &=& (Y^T - \beta^T X^T)(Y - X \beta) \\ &=& Y^T Y - Y^T X \beta - \beta^T X^T Y + \beta^T X^T X \beta \end{eqnarray*}$

同样 $J(\beta)$ 的极小值，在其一阶偏导处取到，为了对 $\beta$ 求导，引入矩阵微积分知识：

$\frac{\partial x^T A}{\partial x} = \frac{\partial A^T x}{\partial x} = A$

$\frac{\partial x^T A x}{\partial x} = A x + A^T x$

于是，令

$\frac{\partial \|Y-X \beta\|^2}{\partial \beta}= \frac{\partial \big( Y^T Y - Y^T X \beta - \beta^T X^T Y+\beta^T X^T X \beta \big)}{\partial \beta} = 2X^T X \beta - 2X^T Y = 0$

若满秩，则可逆，从而可得 $\beta= (X^T X)^{-1} X^T Y$ .

注1：最小二乘法求解需要矩阵满秩，否则只能采用梯度下降法求解.

注2：矩阵形式非常便于Matlab求解，同时也可以很容易地推广到多元线性回归：取

$X = \begin{bmatrix} 1 & x_1^1 & \cdots & x_m^1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_1^n & \cdots & x_m^n \end{bmatrix}$

即可，对应个自变量，组样本，第组样本为 $(x_1^i, x_2^i, \cdots, x_m^i, y_i)$ ，结果形式不变： $\beta= (X^T X)^{-1} X^T Y$ , 此时 $\beta = [\beta_0, \beta_1, \cdots, \beta_m]^T$ .

注3：一些非线性回归也可以转化为线性回归来做：

例如，人口指数增长模型 $y=a e^{bx}$ ，做对数变换： $\ln y = \ln a + bx$ . 即将的数据取对数作为因变量，再与自变量数据做线性回归，得到回归系数 $\beta_0, \beta_1$ ，再由 $\beta_0 = \ln a, \beta_1 = b$ 可得到 $a= e^{\beta_0}, b= \beta_1$ .

其它可变换为线性回归的函数形式：

再例如，自变量 , 因变量，想做如下非线性回归：

$y=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1^2 + \beta_4 x_1 x_2 + \beta_5 x_2^2$

令，做元线性回归即可。

四、最小二乘法的概率解释

真实数据中，一个值可能对应多个值，因为实际值可能受多种因素的影响，故可以假设任意一个值对应的的真实值是服从正态分布的随机变量。

那么，什么时候拟合效果最好，当然是预测值 $\hat{Y}=X \beta$ 取值概率最大的时候，即最大似然原理。

记

$Y=X\beta + \varepsilon$

其中，
$Y = \begin{bmatrix} y_1 \\ \vdots \\ y_n \end{bmatrix}$ 为真实值， $\varepsilon = \begin{bmatrix} \varepsilon_1 \\ \vdots \\ \varepsilon_n \end{bmatrix}$ 为预测误差（即不能被线性模型刻画的部分）。

假设 $\varepsilon_i$ 为独立同分布，服从均值为，方差为 $\sigma^2$ 的正态分布。

该假设由中心极限定理可以保证，顺便说一句，理想的误差都得服从均值，等方差的正态分布，否则说明建立的模型不充分，误差中还有趋势没有被提取出来。

则 $\varepsilon_i$ 的概率密度为

$p(\varepsilon_i)=\frac{1}{\sqrt{2 \pi } \sigma} \exp \Big( - \frac{\varepsilon_i^2}{2 \sigma^2} \Big), \quad i = 1, \cdots, n$

这就意味着，在给定和参数 $\beta$ 情况下，因变量也服从正态分布，即

$p(y_i|x_i; \beta)=\frac{1}{\sqrt{2 \pi } \sigma} \exp \Big( - \frac{(y_i - x_i \beta)^2}{2 \sigma^2} \Big), \quad i=1, \cdots, n$

其中， $x_i = (x_1^i, \cdots, x_m^i)$

定义所有样本数据关于参数 $\beta$ 的似然函数：

$L(\beta) = p(y|x;\beta)=\prod_{i=1}^n p(y_i | x_i ; \beta) = \prod_{i=1}^n \frac{1}{\sqrt{2 \pi } \sigma} \exp \Big( - \frac{(y_i - x_i \beta)^2}{2 \sigma^2} \Big)$

为了便于最大化似然函数，先取对数：

$\ln \big( L(\beta) \big) = \ln \bigg(\prod_{i=1}^n \frac{1}{\sqrt{2 \pi } \sigma} \exp \Big( - \frac{(y_i - x_i \beta)^2}{2 \sigma^2} \Big) \bigg)=n \ln \Big( \frac{1}{\sqrt{2 \pi} \sigma} \Big) - \sum_{i=1}^n \frac{ (y_i - x_i \beta)^2}{2 \sigma^2}$

于是最大化对数似然函数，就相当于最小化：

$\sum_{i=1}^n (y_i - x_i \beta)^2 = \|Y-X \beta\|^2 = J(\beta)$

这就从概率学角度完美地解释了最小二乘法的合理性。

五、模型检验

1. 拟合优度检验

计算，反映了自变量所能解释的方差占总方差的百分比：

总平方和： $SST = \sum_{i=1}^n (y_i - \bar{y})^2$

解释平方和： $SSE = \sum_{i=1}^n (\hat{y}_i - \bar{y})^2$

残差平方和： $SSR = \sum_{i=1}^n (y_i - \hat{y})^2 = \sum_{i=1}^n (y_i - x_i \beta)^2$

则 $R^2=\frac{SSE}{SST}=1-\frac{SSR}{SST}$

值越大说明模型拟合效果越好。通常可以认为当 0.9″ loading=”lazy”> 时，所得到的回归直线拟合得很好，而当时，所得到的回归直线很难说明变量之间的依赖关系。

2. 回归方程参数的检验

回归方程反应了因变量随自变量变化而变化的规律，若 $\beta_1 = 0$ ，则不随变化，此时回归方程无意义。所以，要做如下假设检验：

$H_0: \beta_1 = 0, \quad H_1: \beta_1 \neq 0$

① 检验

若 $\beta_1 = 0$ 为真，则回归平方和与残差平方和 $\frac{ESS}{n-2}$ 都是 $\sigma^2$ 的无偏估计，因而采用统计量：

$F=\frac{RSS/\sigma^2 / 1}{ESS / \sigma^2 / (n-2)}=\frac{RSS}{ESS/(n-2)} \sim F(1,n-2)$

来检验原假设 $\beta_1 = 0$ 是否为真。

② 检验

对 $H_0: \beta_1 = 0$ 的检验与检验是等价的（）。

3. 用回归方程做预测

得到回归方程 $\hat{y} = \beta_0 + \beta_1 x$ 后，预测处的值为 $\hat{y}_0 = \beta_0 + \beta_1 x_0$ ，其置信区间为：

$\Big(\hat{y}_0 - t_{\alpha/2} \sqrt{h_0 \hat{\sigma}^2}, \hat{y}_0 + t_{\alpha/2} \sqrt{h_0 \hat{\sigma}^2} \Big)$

其中 $t_{\alpha/2}$ 的自由度为 , $h_0 = \frac{1}{n}+\frac{(x_0 - \bar{x})^2}{\sum_{i=1} ^n (x_i - \bar{x})^2}$ 称为杠杆率， $\hat{\sigma}^2=\frac{ESS}{n-2}$ .

六、回到案例——例1

1. 做一元线性回归

lreg<-lm(sale~cost,data=dat)
summary(lreg)

结果说明：

(1) 回归方程为

sale= 5.58*cost – 23

(2) 回归系数显著不为0，其p值=5.94e-5 < 0.05

(3) 拟合优度R2=0.88，说明拟合效果很好

2. 回归系数的置信区间

confint(lreg,parm="cost",level=0.95)

回归系数 $\beta_1$ 的95%置信区间为[3.9, 7.26]

3. 回归诊断

par(mfrow=c(2,2))
plot(lreg)  #绘制回归诊断图

图1是残差拟合，越没有趋势越好，有趋势说明可能需要二次项；

图2是残差正态性检验，越落在虚线上越好（理想的残差服从0附件的正态分布，否则说明模型不够充分还有趋势没有提取出来）；

图3检验残差是否等方差；

图4检验离群点，第6个样本点偏离较远，应该剔除掉重新做回归。

4. 剔除离群点，重新做一元线性回归

dat2<-dat[-6,]  
lreg2<-lm(sale~cost,data=dat2)
summary(lreg2)

新回归模型拟合优度R2=0.946有了明显改进。新回归方程为：

sale = 5.28*cost – 15.15

plot(dat2)
abline(lreg2)
points(dat[6,],pch=8,col="red")

5. 用回归模型做预测

x<-data.frame(cost=c(100,110,120))
predict(lreg2,x,interval="prediction",level=0.95)

注意：做预测的自变量数据需要存为数据框格式。

参考文献：

《R语言实战》,Robort I. Kabacoff
矩阵简介与最小二乘法, 耿修瑞,中国科学院电子学研究所
最小二乘法解的矩阵形式推导
最小二乘法的极大似然解释
《数学建模与数学实验》，吴刚，张敬信等，中国商业出版社

版权所有，转载请注明

今天的文章【预测算法】01. 线性回归原理及R语言实现分享到此就结束了，感谢您的阅读。

版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容，请发送邮件至举报，一经查实，本站将立刻删除。
如需转载请保留出处：https://bianchenghao.cn/15453.html