统计回归与Matlab软件实现上(一元多元线性回归模型)

引言

关于数学建模的基本方法

1. 机理驱动
   由于客观事物内部规律的复杂及人们认识程度的限制，无法得到内在因果关系，建立合乎机理规律的数学模型
2. 数据驱动
   直接从数据出发，找到隐含在数据背后的最佳模型，是数学模型建立的另一大思路
   统计回归方法是基于统计理论建立的最基本的一类数据驱动建模方法

学习目的

1. 用统计回归方法估计数据中隐含的模型
2. 对模型参数和模型结果的可靠性做必要检验
3. 分析建模过程中的一些细节问题：异常数据的辨识与处理，变量的筛选
4. 用MATLAB软件实现

一元线性回归模型的概念

一元线性回归模型基本概念

一般地，称由
$$y={\beta }_0+{\beta }_{1}x+\epsilon$$
确定的模型为一元线性回归模型

• ${\beta }_0,{\beta }_1$为回归系数
• $x$为自变量、回归变量或解释变量
• $y$为因变量或被解释变量
• $\epsilon$为随机误差

随机误差$\epsilon$的基本假设

1. 高斯-马尔科夫条件
   $$\{\begin{array}{c}E(\epsilon )=0\\ cov(\epsilon ,\epsilon )={\sigma }^2\end{array}$$

• 随机误差项必须是0均值的
• 方差等于${\sigma }^2$，是恒定的，即与$x$的取值无关

2. 正太分布假设
   $$\epsilon \sim N(0,{\sigma }^2)$$
   随机误差项要服从0均值的正太分布，并且方差同样是恒定的，与$x$无关

一元线性回归分析的主要任务

1. 基于样本数据，对参数${\beta }_0,{\beta }_1,\sigma$做参数估计
2. 对模型参数${\beta }_0,{\beta }_1$以及模型显著性作假设检验分析
3. 对$y$的值作预测，即对$y$作点(区间)估计

Matlab实现

[b, bint, r, rint, stats]=regress(Y,X,alpha)

• bint，回归系数的区间估计
• r，残差
• rint，残差的置信区间
• stats，检验回归模型的统计量：决定系数$r^2$，F值，与F值对应的概率p
• alpha，显著性水平，缺省时为0.05

模型的参数估计与软件实现

回归系数的最小二乘估计

有$n$组独立样本：$(x_1,y_1),(x_2,y_2),\dots ,(x_n,y(n))$，带入回归方程可得
$$\{\begin{array}{c}{y}_i={\beta }_0+{\beta }_1{x}_i+{\epsilon }_i,i=1,2,\dots ,n\\ E({\epsilon }_i)=0,var({\epsilon }_i)={\sigma }^2\end{array}$$
其中，${\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_n$相互独立

• 拟合误差或残差：$r_i=y_i-y_i^{\prime }$
• 最好直线：使残差平方和最小的直线
  $$Q({\beta }_0,{\beta }_1)=\sum _{i=1}^n(y_i-y_i^{\prime }{)}^2=\sum _{i=1}^n(y_i-{\beta }_0-{\beta }_i{x}_i{)}^2$$
  最小化的参数值${\beta }_0^{\prime },{\beta }_1^{\prime }$称为${\beta }_0,{\beta }_1$的最小二乘估计

该优化问题的求解，可以基于极值原理实现
通过残差平方和，分别对${\beta }_0,{\beta }_1$求偏导数，令偏导数等于0
$$\{\begin{array}{c}\frac{\partial Q}{\partial {\beta }_0}=0\\ \frac{\partial Q}{\partial {\beta }_1}=0\end{array}$$
得到的是二元一次线性方程组

相应的最小二乘估计为
$$\{\begin{array}{c}{\hat{\beta }}_0=\bar{y}-{\hat{\beta }}_1\bar{x}\\ {\hat{\beta }}_1=\frac{\bar{x}\bar{y}-\bar{x}\bar{y}}{\overline{x^2}-{\bar{x}}^2}\end{array}$$
其中
$$\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i,\bar{y}=\frac{1}{n}\sum _{i=1}^n{y}_i,$$
$${\bar{x}}^2=\frac{1}{n}\sum _{i=1}^n{x}_i^2,\bar{x}\bar{y}=\frac{1}{n}\sum _{i=1}^n{x}_i{y}_i$$

Matlab实现

regress命令

b=regress(Y, X)

待求解的线性方程组
$$y_i={\beta }_0+{\beta }_1{x}_i+{\epsilon }_i,i=1,2,\dots ,n$$
$$\left[\begin{array}{c}{y}_1\\ y_2\\ \dots \\ y_n\end{array}\right]=\left[\begin{array}{ccc}1& & x_1\\ 1& & x_2\\ \dots & & \dots \\ 1& & x_n\end{array}\right]\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\end{array}\right]$$

• Y指的是y的矩阵向量
• X指的是系数矩阵x
• 等式左边的b指的是参数$\beta$的最小二乘估计

1. 输入数据

x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164]';
X=[ones(16,1)x];
Y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]';

x表示所有的自变量，16各人的身高数据，是列向量，‘是转秩
X是两列的矩阵，第一个是16个1组成的列向量，第二个是自变量x的列向量
Y是16个腿长数据，也是列向量
2. 参数估计

b=regress(Y, X)

得b

b=[-16.0730 0.7194]'

$${\beta }_0^{\prime }=-16.0730;{\beta }_1^{\prime }=0.7194$$
经验回归方程
$$y=-16.0739+0.7194x$$

回归系数的置信区间估计

在正太假设的前提下
$${\hat{\beta }}_0\sim N\left({\beta }_0,\left(\frac{1}{n}+\frac{{\hat{x}}^2}{L_{xx}}\right){\sigma }^2\right)$$
$$\widehat{{\beta }_1}\sim N\left({\beta }_1,\frac{{\sigma }^2}{L_{xx}}\right)$$
其中
$$L_{xx}=\sum _{i=1}^n(x_i-\bar{x}{)}^2$$

由于$\sigma$未知，可以构造t统计量来进行区间估计
$$t=\frac{{\beta }_1^{\prime }-{\beta }_1}{\sqrt{\frac{({\sigma }^{\prime }{)}^2}{L_{xx}}}}\sim t(n-2)$$
其中
$$L_{xx}=\sum _{i=1}^n(x_i-\bar{x}{)}^2$$
$${\hat{\sigma }}^2=\frac{1}{n-2}\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$$

$$P\left(-t_{\frac{a}{2}}(n-2)<\frac{\widehat{{\beta }_1}-{\beta }_2}{\sqrt{\frac{{\hat{\sigma }}^2}{L_{xx}}}}<t_{\frac{a}{2}}(n-2)\right)=1-\alpha$$
t统计量落在这两个值之间的概率是$1-\alpha$

故${\beta }_1$的置信水平为$1-\alpha$的置信区间估计为
$$\left[{\beta }_1^{\prime }-t_{\frac{a}{2}}(n-2)\sqrt{\frac{({\sigma }^{\prime }{)}^2}{L_{xx}}},{\beta }_1^{\prime }+t_{\frac{a}{2}}(n-2)\sqrt{\frac{({\sigma }^{\prime }{)}^2}{L_{xx}}}\right]$$
同理也可以得到${\beta }_0$的置信区间估计

Matlab实现

[b, bint]=regress(Y, X, 0.05)

bint就是区间估计结果
Y，隐变量的取值向量
X，系数矩阵
0.05，置信水平=1-0.05=95%

用上面的数据得

b = -16.0730 0.7194bint = -33.7071  1.56120.6047   0.8340

bint第一行是${\beta }_0^{\prime }$的置信区间估计结果
第二行是${\beta }_1^{\prime }$的置信区间估计结果
置信水平是95%的区级估计
$$[-33.71,1.56]和[0.60,0.83]$$

模型的残差分析

残差分析的基本概念

残差：样本的观测值与样本的预测值之差
$$r_i=y_i-y_i^{\prime }$$
残差向量：所有样本的拟合误差，组成的列向量
$$r=\left[\begin{array}{c}{r}_1\\ r_2\\ \dots \\ r_n\end{array}\right]$$
残差应该满足的一些基本性质
0均值
$$E({\epsilon }_i)=0$$
残差与残差之间是不相关的
$$cov({\epsilon }_i,{\epsilon }_j)=0$$
残差的方差
$$var({\epsilon }_i)=(1-h_{ii}){\sigma }^2$$

残差图分析

横坐标是自变量x，纵坐标是残差
残差是在0附近随机波动，残差与残差之间不存在明显的关联性

异方差现象，与x有关系，不符合条件

不是0均值，残差与残差之间有联系

前后之间有关联，不是随机波动

软件实现

简单残差图命令

plot(r, '*')
hold on
ezplot('0',[1,length(r)])

用plot命令画出残差，用星号表示
hold on，表示前面的不要擦除，继续画图
ezplot，画出0的基准线

Matlab残差图作图命令

rcoplot(r, rint)

r表示残差向量
rint表示残差的置信区间

中间的圆圈，表示残差
每个残差都有区间线段，表示置信区间

• 一般认为，正常的样本，残差的置信区间，应该是要跨越0的
• 如果置信区间明显远离0，表示这个样本是异常的

模型的检验与软件实现

模型检验之决定系数

总体平方和
$$TSS=\sum _{i=1}^n(y_i-\bar{y}{)}^2$$
隐变量的观测值减去观测值得平均值的平方和
能够反应样本观测值与中心的偏离程度

能够近似衡量样本观测值序列所包含的信息的多少

TSS的分解
$$TSS=\sum _{i=1}^n(y_i-\bar{y}{)}^2=\sum _{i=1}^n(y_i-\widehat{y_i}+\widehat{y_i}-\bar{y}{)}^2$$
$$=\sum _{i=1}^n({\hat{y}}_i-\bar{y}{)}^2+(y_i-{\hat{y}}_i{)}^2+2(y_i-{\hat{y}}_i)({\hat{y}}_i-\bar{y})$$
$$=\sum _{i=1}^n({\hat{y}}_i-\bar{y}{)}^2+\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2+\sum _{i=1}^{n}2(y_i-{\hat{y}}_i)({\hat{y}}_i-\bar{y})$$
交叉项的和严格等于0
$$\sum _{i=1}^n(y_i-\bar{y}{)}^2=\sum _{i=1}^n({\hat{y}}_i-\bar{y}{)}^2+\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$$
总体平方和=回归平方和(ESS)+残差平方和(RSS)

• 回归平方和表示的是模型能够解释的那一部分平方和的信息，反应的是回归模型能够解释的观测值中的信息的多少
• 残差平方和表示模型没有学习到的信息的多少

决定模型$R^2$统计量：
$$R^2=\frac{ESS}{TSS}=1-\frac{RSS}{TSS}$$
$R^2$也被称为判定系数或拟合优度

• 取值范围一定在01之间
• 越接近1，样本数据拟合效果越好

Matlab实现

[b, bint, r, rint, stats] = regress(Y, X, 0.05)

得

stats:0.9282 180.9531 0.0000 1.7437

因此$R^2=0.9282$

模型检验之F统计量检验

原假设$H_0$：回归方程$y={\beta }_0+{\beta }_{1}x$不显著成立
也就是线性项可有可无，即${\beta }_1=0$
备择假设$H_1$回归方程$y={\beta }_0+{\beta }_{1}x$显著成立
即${\beta }_1\ne 0$

在$H_0$成立的假定下，构造统计量
$$F=\frac{\frac{ESS}{1}}{\frac{RSS}{n-2}}\sim F(1,n-2)$$
其中
$$ESS=\sum _{i=1}^n({\hat{y}}_i-\bar{y}{)}^2$$
自由度是1
$$RSS=\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$$
自由度是n-2

概率密度曲线
$F_{\alpha }$：上$\alpha$分位点，临界值点
如果F值，大于临界值，就拒绝原假设，即线性回归模型显著
如果F值。小于临界值，接受原假设，即线性回归模型不显著

Matlab实现

[b, bint, r, rint, stats] = regress(Y, X, 0.05)

得

stats:0.9282 180.9531 0.0000 1.7437

因此F值=180.9531
因为数据是16个人的数据，临界值$F_{\alpha }(1,n-2)$，就是$F_{\alpha }(1,14)$
$\alpha$取0.05
可以查询F分布表，查到分位点
或

x_a = finv(0.95, 1, 14)

0.95表示落在临界值左侧的概率
1和14分别表示F分布的两个自由度
返回值就是临界值
得：$F_{0.05}(1,14)=4.6001$
有$F值\gg F_{0.05}(1,14)$，可以得出拒绝原假设得结论，所以线性回归关系是显著成立的

与F值对应的p值

F值对应的右侧的这一块面积，就是p值
是分布落在F值右边的概率
当原假设成立的前提下，自由度是1和n-2的随机变量落在F值右侧的概率
$$p=P(F(1,n-2)>F值|H_0成立)$$

• p值可以理解为接受回归模型的风险，即犯错的概率

Matlab实现

[b, bint, r, rint, stats] = regress(Y, X, 0.05)

得

stats:0.9282 180.9531 0.0000 1.7437

p值是0.000

>> stats(3)

输入以上语句来得到p值

ans=2.1312e-09

$$2.1312\times 10^{-9}$$
因此，接受回归模型正确的风险为2.1312e-09

${\sigma }^2$的无偏估计

stats的第四个统计指标
${\sigma }^2$是模型的随机误差项的方差
$${\hat{\sigma }}^2=\frac{1}{n-2}\sum _{i=1}^n{\epsilon }_i^2=\frac{1}{n-2}\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$$
残差平方和除以自由度n-2

模型预测

1. 点预测
   将对$x_0$代入经验回归方程，得点预测结果
   $${\hat{y}}_0={\hat{\beta }}_0+{\hat{\beta }}_1{x}_0$$
2. 区间预测
   置信水平$1-\alpha$下，对$y_0$进行区间估计
   $$[\hat{y}-\delta (x_0),\hat{y}+\delta (x_0)]$$
   其中
   $$\delta (x_0)=\hat{\sigma }\sqrt{1+\frac{1}{n}+\frac{(x_0-\bar{x}{)}^2}{L_{xx}}}{t}_{\frac{a}{2}}(n-2)$$
   $$X^{T}X=(C_{ij})$$

Matalab实现

y_hat = b(1)+b(2)*x
plot(x, Y, 'k+', x, y_hat, 'r')

把x的取值向量，直接代入到经验回归模型当中
b(1)就是${\beta }_0$
b(2)就是${\beta }_1$
y_hat是预测值
Y表示原始的样本观测值，用黑色加号表示
预测值用红色的实线表示

多元线性回归模型与软件实现

基本概念

一般地，称由
$$y={\beta }_0+{\beta }_1{x}_2+\cdots +{\beta }_m{x}_m+ϵ$$
确定的模型，为m元线性回归模型，也可表示为矩阵形式
$$\{\begin{array}{c}Y=X\beta +ϵ\\ E(\epsilon )=0,cov(\epsilon ,\epsilon )={\sigma }^2{I}_n\end{array}$$
其中
$$Y=\left(\begin{array}{c}{y}_1\\ y_2\\ \dots \\ y_n\end{array}\right)X=\left(\begin{array}{ccccccccc}1& & x_{11}& & x_{12}& & \dots & & x_{1m}\\ 1& & x_{12}& & x_{22}& & \dots & & x_{2m}\\ \dots & & \dots & & \dots & & \dots & & \dots \\ 1& & x_{1n}& & x_{n2}& & \dots & & x_{nm}\end{array}\right)$$
$$\beta =\left(\begin{array}{c}{\beta }_0\\ {\beta }_1\\ \dots \\ {\beta }_n\end{array}\right)\epsilon =\left(\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ \dots \\ {\epsilon }_n\end{array}\right)$$

主要任务

1. 对参数$\beta$和${\sigma }^2$作点估计
2. 对模型参数、模型显著性作检验分析
3. 对$y$的值作预测，即对$y$作点(区间)估计

模型参数的估计

用最小二乘法对${\beta }_0,{\beta }_1\dots {\beta }_m$进行参数估计
$$mi{n}_{{\beta }_0,{\beta }_1,\dots ,{\beta }_m}Q=\sum _{i=1}^n(y-{\beta }_0-{\beta }_1{x}_1-\cdots -{\beta }_m{x}_m{)}^2$$
解得最小二乘估计为
$$\hat{\beta }=(X^{T}X{)}^{-1}(X^{T}Y)$$

模型的检验

类似于一元线性回归情形

1. 拟合优度检验
2. 方程显著性的F检验
3. 变量显著性的t检验
   基于t统计量：
   $$t=\frac{{\hat{\beta }}_i}{\hat{\sigma }\sqrt{c_{ii}}}\sim t(n-k-1)$$
   对参数${\beta }_i$进行显著性检验$(H_0:{\beta }_i=0)$
   其中，$X^{T}X=(c_{jj})$

模型的预测

1. 点预测
   将对$(x_1^{\ast },x_2^{\ast },\dots ,x_m^{\ast })$代入经验回归方程，得点预测结果
   $${\hat{y}}^{\ast }=\widehat{{\beta }_0}+\widehat{{\beta }_1}{x}_1^{\ast }+\widehat{{\beta }_2}{x}_2^{\ast }+\cdots +\widehat{{\beta }_m}{x}_m^{\ast }$$
2. 区间预测
   $$\left[\begin{array}{c}\hat{y}-\widehat{\sigma t_{\frac{a}{2}}(n-k-1)\sqrt{1+{\sum }_{i,j}{c}_{ij}{x}_i{x}_j},\hat{y}+\hat{\sigma }{t}_{\frac{a}{2}}(n-k-1)\sqrt{1+{\sum }_{i,j}{c}_{ij}{x}_i{x}_j}}\end{array}\right]$$
   其中$X^{T}X=(c_{ij})$

Matlab实现

[b, bint, r, rint, stats]=regress(Y, X, alpha)

• 拟合优度$r^2$越接近1，说明回归方程越显著
• $F>F_{\alpha }(m,n-m-1)$时拒绝H0，$F$越大，说明回归方程越显著
• 与$F$对应得概率$p<\alpha$时拒绝H0，回归模型成立

建材销售量的回归模型

求解

1. 建立建材销售量$y$与推销开支$x_1$、实际账目数$x_2$、同类商品竞争数$x_3$和地区销售潜力$x_4$的线性回归模型
   $$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+{\beta }_4{x}_4+ϵ$$
2. 用regress命令进行线性回归模型求解

x1=[5.5 2.5 8 3 3 2.9 8 9 4 6.5 5.5 5 6 5 3.5 8 6 4 7.5 7]';
x2=[31 55 67 50 38 71 30 56 42 73 60 44 50 39 55 70 40 50 62 59]';
x3=[10 8 12 7 8 12 12 5 8 5 11 12 6 10 10 6 11 11 9 9]';
x4=[8 6 9 16 15 17 8 10 4 16 7 12 6 4 4 14 6 8 13 11]';
y=[79.3 200.1 163.2 200.1 146 177.7 30.9 291.9 160 339.4 159.6 86.3 237.5 ... 107.2 155 201.4 100.2 135.8 223.3 195]';
x=[ones(size(x1)), x1, x2, x3, x4];
[b, bint, r, rint, stats]=regress(y, X);

3. 程序求解结果分析
   

• y的90.34%可由模型确定
• F值远超过临界值$F_{0.05}(4.15)=3.0556$
• p值远小于$\sigma =0.05$
• 模型整体上成立
• ${\beta }_1和{\beta }_4$置信区间包含零点，$x_1,x_4$对y的影响不太显著

4. 模型的残差分析
   

• 第16个样本为异常样本

5. 模型的改进
   首先，剔除异常样本，并重新回归计算

y(16)=[];
x(16,:)=[];
[b, bint, r, rint, stats]=regress(y, X);

$R^2。F$都有较大改进，但回归系数${\beta }_4$的置信区间包含零点
剔除不显著的变量$x_4$
$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+ϵ$$

[b, bint, r, rint, stats]=regress(y, X(:,1:end-1))

• 置信区间越短了
  

• 残差图基本正常

• 最终模型具有较好的显著性