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

阅读 289
收藏 5
2018-01-15
原文链接:zhuanlan.zhihu.com

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

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

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

一、最小二乘法

设有 n 组样本点: (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 是与横轴 x_i 对应的直线上的点的纵坐标(称为线性模型预测值),它与样本点 x_i 对应的真实值 y_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) 的极小值,在其一阶偏导 =0 处取到,为了对 \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

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

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

注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}

即可,对应 m 个自变量, n 组样本,第 i 组样本为 (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 . 即将 y 的数据取对数作为因变量,再与自变量 x 数据做线性回归,得到回归系数 \beta_0, \beta_1 ,再由 \beta_0 = \ln a, \beta_1 = b 可得到 a= e^{\beta_0}, b= \beta_1.

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

再例如,自变量 x_1, x_2 , 因变量 y ,想做如下非线性回归:

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

X_1 = x_1, X_2 = x_2, X_3 = x_1^2, X_4 = x_1 x_2, X_5= x_2^2 ,做 5 元线性回归即可。

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


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

那么,什么时候拟合效果最好,当然是预测值 \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 为独立同分布, 服从均值为 0 ,方差为 \sigma^2 的正态分布。

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

\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

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

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. 拟合优度检验

计算 R^2 ,反映了自变量所能解释的方差占总方差的百分比:

总平方和: 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}

R^2 值越大说明模型拟合效果越好。通常可以认为当 R^2 > 0.9 时,所得到的回归直线拟合得很好,而当 R^2 < 0.5 时,所得到的回归直线很难说明变量之间的依赖关系。

2. 回归方程参数的检验

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

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

F 检验

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

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

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

t 检验

H_0: \beta_1 = 0t 检验与 F 检验是等价的( t^2 = F )。


3. 用回归方程做预测

得到回归方程 \hat{y} = \beta_0 + \beta_1 x 后,预测 x=x_0 处的 y 值为 \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} 的自由度为 n-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)

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

参考文献:

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

版权所有,转载请注明

评论