数学建模笔记—— 回归分析

回归分析

在统计学中，回归分析( regression analysis)指的是确定俩种或两种以上变量间相互依赖的定量关系的一种统计分析方法。回归分析按照涉及的变量的多少，分为一元回归和多元回归分析；按照因变量的多少，可分为简单回归分析和多重回归分析；按照自变量和因变量之间的关系类型，可分为线性回归分析和非线性回归分析。

在大数据分析中，回归分析是一种预测性的建模技术，它研究的是因变量(目标)和自变量(预测器) 之间的关系。这种技术通常用于预测分析以及发现变量之间的因果关系。例如，司机的鲁莽驾驶与道路交通事故数量之间的关系，最好的研究方法就是回归。

1. 回归分析的一般步骤

1. 确定回归方程中的解释变量和被解释变量
2. 确定回归模型和建立回归方程
3. 对回归方程进行各种检验
4. 利用回归方程进行预测

2. 一元线性回归分析

2.1 具体过程

2.1.1 确定回归方程中的解释变量和被解释变量

• 因变量：被预测或被解释的变量，用$y$表示
• 自变量：预测或解释因变量的一个或者多个变量，用$x$表示
• 对于只有线性关系的两个变量，可以用一个方程来表示它们之间的线性关系
• 描述因变量$y$如何依赖于自变量$x$和误差项$ϵ$的方程称为回归模型

2.1.2 确定回归模型和建立回归方程

对于只涉及一个自变量的一元线性回归模型可表示为
$$y={\beta }_0+{\beta }_{1}x+ϵ$$
在这个模型里：

• $y$叫做因变量或被解释变量
• $x$叫做自变量或解释变量
• ${\beta }_0$表示截距
• ${\beta }_1$表示斜率
• $ϵ$表示误差项，反映除$x$和y之间的线性关系之外的随机因素对$y$的影响，是不可避免的

回归方程：

• 描述因变量$y$的期望值如何依赖于自变量$x$的方程称为回归方程。根据对一元线性回归模型的假设，可以得到它的回归方程为：

$$E(y)={\beta }_0+{\beta }_{1}x$$

• 如果回归方程中的参数已知，对于一个给定的$x$值，利用回归方程就能计算出$y$的期望值

• 用样本统计量代替回归方程中的未知参数，就得到估计的回归方程，简称回归直线

参数的最小二乘法估计：

对于回归直线，关键在于求解参数，常用高斯提出的最小二乘法，它是使因变量的观察值$y$与估计值之间的离差平方和达到最小来求解的
$$Q=\sum (y-\hat{y}{)}^2=\sum (y-{\hat{\beta }}_0-{\hat{\beta }}_{1}x{)}^2$$
展开可得：
$$Q=\sum (y-\hat{y}{)}^2=\sum y^2+n{\hat{\beta }}_0^2+{\hat{\beta }}_1^2\Sigma x^2+2{\hat{\beta }}_0{\hat{\beta }}_1\Sigma x-2{\hat{\beta }}_0\Sigma y-2{\hat{\beta }}_1\Sigma xy$$
为了求该式子的极值，对${\beta }_0$和$\beta 1$求偏导得：
$$\{\begin{array}{l}\sum y=n{\hat{\beta }}_0+{\hat{\beta }}_1\Sigma x\\ \sum xy=2{\hat{\beta }}_0\Sigma x+{\hat{\beta }}_1\Sigma x^2\end{array}$$
求解：
$$\{\begin{array}{l}{\hat{\beta }}_1=\frac{n\sum xy-\sum x\sum y}{n\sum x^2-(\sum x{)}^2}\\ {\hat{\beta }}_0=\bar{y}-{\hat{\beta }}_1\bar{x}\end{array}$$

2.1.3 利用回归直线进行估计和预测

• 点估计：利用估计的回归方程，对于$x$的某一个特定的值，求出$y$的一个估计值就是点估计
• 区间估计：利用估计的回归方程，对于$x$的一个特定值，求出$y$的一个估计值的区间就是区间估计

估计标准误差的计算：

为了度量回归方程的可靠性，通常计算估计标准误差，它度量观察值回绕着回归直线的变化程度或分散程度，估计标准误差定义为：

$$S_e=\sqrt{\frac{\sum (y-\hat{y}{)}^2}{n-2}}$$

• 估计标准误差越大，则数据点围绕回归直线的分散程度就越大，回归方程代表性就越小
• 估计标准误差越小，则数据点围绕回归直线的分散程度越小，回归方程的代表性越大，可靠性越高

置信区间和预测区间的估计：

置信区间:
$$\begin{gathered} {\hat{y}}_0\pm t_{\frac{\alpha }{2}}{S}_e\sqrt{\frac{1}{n}+\frac{(x_0-\bar{x}{)}^2}{\sum (x-\bar{x}{)}^2}} \\ 其中\alpha 为显著性水平，t_{\alpha }可以通过查表得，S_e是估计标准误差 \end{gathered}$$
预测区间：
$$\begin{gathered} {\hat{y}}_0\pm t_{\frac{\alpha }{2}}{s}_e\sqrt{1+\frac{1}{n}+\frac{(x_0-\bar{x}{)}^2}{\sum (x-\bar{x}{)}^2}} \\ 其中\alpha 为显著性水平，t_{\alpha }可以通过查表得，S_e是估计标准误差 \end{gathered}$$

2.1.4 对回归方程进行各种检验(补充)

1. 回归直线的拟合优度

回归直线与各观测点的接近程度称为回归直线对数据的拟合优度

• 总平方和 (TSS)：反映因变量的$n$个观察值与其均值的总离差
  $$TSS=\sum y_i^2=\sum (y_i-{\bar{y}}_i{)}^2$$

• 回归平方和 (ESS)：反映了$y$的总变差中，由于$x$与$y$之间的线性关系引起的$y$的变化部分
  $$ESS=\sum {\hat{y}}_i^2=\sum ({\hat{y}}_i-{\bar{y}}_i{)}^2$$

• 残差平方和(RSS)：反映了除了$x$对$y$的线性影响之外的其他因素对$y$变差的作用，是不能由回归直线来解释的y的变差部分
  $$RSS=\sum e_i^2=\sum (y_i-{\hat{y}}_i{)}^2$$

总平方和可以分解为回归平方和、残差平方和两部分：
$$TSS=ESS+RSS$$
判定系数：

回归平方和占总平方和的比例，用$R$表示，其值在0到1之间
$$R^2=\frac{\mathrm{ESS}}{T\mathrm{SS}}=1-\frac{\mathrm{RSS}}{T\mathrm{SS}}=\frac{\sum ({\hat{y}}_i-{\bar{y}}_i{)}^2}{\sum (y_i-{\bar{y}}_i{)}^2}$$

• $R^2=0$：说明$y$的变化与$x$无关，$x$无助于解释$y$的变差
• $R^2=1$：说明残差平方和为0，拟合是完全的，$y$的变化只与$x$有关

2. 显著性检验

显著性检验的主要目的是根据所建立的估计方程用自变量$x$来估计或预测因变量$y$的取值，当建立了估计方程后，还不能马上进行估计或预测，因为该估计方程是根据样本数据得到的，它是否真实的反映了变量$x$和$y$之间的关系，则需要通过检验后才能证实。

根据样本数据拟合回归方程时，实际上就已经假定变量$x$和$y$之间存在着线性关系，并假定误差项是一个服从正态分布的随机变量，且具有相同的方差，但这些假设是否成立需要检验。

显著性检验包括两方面：

• 线性关系检验

  线性关系检验是检验自变量$x$和因变量y之间的线性关系是否显著，或者说，它们之间能否用一个线性模型来表示

  检验原理：

  将均方回归(MSR)同均方残差(MSE)加以比较，应用F检验来分析二者之间的差别是否显著

  • 均方回归(MSR)：回归平方和ESS除以相应的回归自由度(自变量的个数k)
  • 均方残差(MSE)：残差平方和RSS除以相应的残差自由度(n-k+1)

  检验过程：

  • $H_0$(原假设)：${\beta }_1=0$,回归系数与0无显著差异，$y$与$x$的线性关系不显著

  • $H_1$:${\beta }_1\ne 0$, 回归显著，认为$y$与$x$存在线性关系，所求的线性回归方程有意义

  • 计算检验统计量$F:$

    $H_0$成立时
    $$F=\frac{ESS/1}{RSS/(n-2)}=\frac{MSR}{MSE}\sim F(1,n-2)$$
    若$F>F_{1-\alpha }(1,n-2)$，拒绝$H_0$，否则接受$H_0$

• 回归系数检验

  回归系数显著性检验的目的是通过检验回归系数$\beta$的值与0是否有显著性差异，来判断$Y$与$X$之间是否有显著的线性关系，若$\beta =0$,则总体回归方程中不含$X$项(即$Y$不随$X$变动而变动),因此，变量$Y$与$X$之间并不存在线性关系；若$\beta \ne 0$, 说明变量$Y$与$X$之间存在显著的线性关系。

  ${\hat{\beta }}_1$是根据最小二乘法求出的样本统计量，服从正态分布

  ${\hat{\beta }}_1$的分布具有如下性质：

  • 数学期望：$E({\hat{\beta }}_1)={\beta }_1$

  • 标准差：${\sigma }_{{\hat{\beta }}_1}=\frac{\sigma }{\sqrt{\sum x_i^2-\frac{1}{n}(\sum x_i{)}^2}}$

  由于$\sigma$未知，需用其估计量$S_{\mathrm{e}}$来代替得到${\hat{\beta }}_1$的估计标准差
  $$S_{{\hat{\beta }}_1}=\frac{S_e}{\sqrt{\sum x_i^2-\frac{1}{n}(\sum x_i{)}^2}}{S}_e=\sqrt{\frac{\sum (y_i-{\hat{y}}_i{)}^2}{n-k-1}}=\sqrt{MSE}$$
  $t$检验的统计量：
  $$t=\frac{{\hat{\beta }}_1-{\beta }_1}{s_{{\hat{\beta }}_1}}\sim t(n-2)$$

线性关系检验与回归系数检验的区别：
线性关系的检验是检验自变量与因变量是否可以用线性来表达，而回归系数的检验是对样本数据计算的回归系数检验是否为0

• 在一元线性回归中，自变量只有一个，线性关系检验与回归系数检验是等价的
• 在多元回归分析中，这两种检验的意义是不同的。线性关系检验只能用来检验总体回归关系的显著性，而回归系数检验可以对各个回归系数分别进行检验

2.2 典型例题

某团队测了16名成年女子的身高与腿长所得数据如下

身高	143	145	146	147	49	153	154	155	156	157	158	159	160	162	164
腿长	88	85	88	91	92	93	95	96	98	97	96	98	999	100	102

我们发现，腿长和身高有很强的正相关性，身高越高，腿长越长，那么通过这组数据，我们可以预测身高170的成年女性腿长是多长吗？

2.2.1 进行手算

一元线性表达式为：
$$\hat{y}={\hat{\beta }}_0+{\hat{\beta }}_{1}x$$
对参数进行最小二乘估计：
$$\{\begin{array}{l}{\hat{\beta }}_1=\frac{n\sum xy-\sum x\sum y}{n\sum x^2-(\sum x{)}^2}=\frac{16\times 232566-2458\times 1511}{16\times 378220-2458^2}=0.7194\\ {\hat{\beta }}_0=\bar{y}-{\hat{\beta }}_1\bar{x}=94.4375-0.7194\times 153.625=-16.08\end{array}$$
将估计得到的参数代回一元线性表达式：
$$\hat{y}=0.794x-16.08$$
进行预测：
$$x=170,\hat{y}=106.218$$
预测身高170的成年女性腿长是106.218cm

计算估计标准误差：
$$s_e=\sqrt{\frac{\sum (y-\hat{y}{)}^2}{n-2}}=\sqrt{1.744}=1.32$$
我们令显著性水平$\alpha =0.05$，则置信水平$1-\alpha =0.95$

自由度为16，计算$t_{\frac{\alpha }{t}}$
$$t_{\frac{\alpha }{2}}(n-k)=t_{0.025}(16-2)=t_{0.025}(14)=2.145(\text{查表})$$
置信区间为：
$${\hat{y}}_0\pm t_{\frac{\alpha }{2}}{s}_e\sqrt{\frac{1}{n}+\frac{(x_0-\bar{x}{)}^2}{\sum (x-\bar{x}{)}^2}}=106.218\pm 2.007$$
预测区间为：
$${\hat{y}}_0\pm t_{\frac{\alpha }{2}}{s}_e\sqrt{1+\frac{1}{n}+\frac{(x_0-\bar{x}{)}^2}{\sum (x-\bar{x}{)}^2}}=106.218\pm 3.472$$
判定系数为：
$$R^2=\frac{\mathrm{ESS}}{T\mathrm{SS}}=1-\frac{\mathrm{RSS}}{T\mathrm{SS}}=\frac{\overline{\sum }({\hat{y}}_i-{\bar{y}}_i{)}^2}{\sum (y_i-{\bar{y}}_i{)}^2}=1-\frac{\sum (y_i-{\hat{y}}_i{)}^2}{\sum (y_i-{\bar{y}}_i{)}^2}=0.9282,\text{ 是比较接近与1的}$$
进行线性关系检验：
$$F=\frac{ESS/1}{RSS/(n-2)}=\frac{MSR}{MSE}\sim F(1,n-2)=180.9631>4.6$$
所以线性关系显著

2.2.2 使用matlab进行计算

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

• 输入变量

  • $X,Y$——自变量和因变量的样本值
  • $alpha$——显著性水平，默认为0.05

• 输出变量

  • $b$——回归系数

  • $bint$——回归系数的区间估计

  • $r$——残差

  • $rint$——置信区间

  • $stats$——用于检验回归模型的统计量

    $stats$有四个数值：决定系数$R^2$、$F$值、与$F$对应的概率$P$、无偏估计$\sigma$

具体代码：

% 1. 输入数据
% 输入X的样本值
x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164]';
% 插入\beta_0对应列
X=[ones(16,1) x];
% 输入Y的样本值
Y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]';

% 2. 回归分析及检验
[b,bint,r,rint,stats]=regress(Y,X);
b,bint,stats

% 3. 残差分析，作残差图
rcoplot(r,rint)

% 4. 预测及作图
z=b(1)+b(2)*x;
plot(x,Y,'k+',x,z,'r')

输出：

b =

  -16.0730
    0.7194


bint =

  -33.7071    1.5612
    0.6047    0.8340


stats =

    0.9282  180.9531    0.0000    1.7437

得到的残差图：

得到的回归曲线：

3. 多元线性回归分析

经常会遇到某一现象的发展和变化取决于(几个影响因素的情况,也就是一个因变量和几个自变量有依存关系的情况，这时需用多元线性回归分析，多元线性回归分析预测法，是指通过对两个或两个以上的自变量与一个因变量的相关分析,建立预测模型进行预测和控制的方法

3.1 算法介绍

3.1.1 多元线性回归模型

一般称由$y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_k{x}_k$确定的模型：
$$\{\begin{array}{c}Y=X\beta +ϵ\\ E(ϵ)=0,COV(ϵ,ϵ)={\sigma }^2{I}_n\end{array}$$
为$k$元线性回归模型，并简记为$(Y,X\beta ,{\sigma }^2{I}_n)$

其中
$$Y=\left[\begin{array}{c}{y}_1\\ y_2\\ ...\\ y_n\end{array}\right],X=\left[\begin{array}{ccccc}1& x_{11}& x_{12}& \cdots & x_{1k}\\ 1& x_{21}& x_{22}& \cdots & x_{2k}\\ ...& ...& ...& ...& ...\\ 1& x_{n1}& x_{n2}& \cdots & x_{nk}\end{array}\right],\beta =\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ...\\ {\beta }_k\end{array}\right],ϵ=\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ...\\ {ϵ}_n\end{array}\right]$$
$y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_k{x}_k$称为回归平面方程

3.1.2 线性模型$(Y,X\beta ,{\sigma }^2{I}_n)$考虑的主要问题

• 对参数β和${\sigma }^2$作点估计，建立$y$与$x_1,x_2,\cdot \cdot \cdot ,x_k$之间的数量关系
• 对模型参数、模型结果等做检验
• 对$\gamma$的值作预测，即对$\gamma$作点(区间)估计

3.1.3 多元线性回归模型的参数估计

用最小二乘法求${\beta }_0,...,{\beta }_k$的估计量：作离差平方和
$$Q=\sum _{i=1}^n(y_i-{\beta }_0-{\beta }_1{x}_{i1}-\cdots -{\beta }_k{x}_{ik}{)}^2$$
选择${\beta }_0,...,{\beta }_k$使$Q$达到最小，解得估计值：
$$\hat{\beta }=(X^{T}X{)}^{-1}(X^{T}Y)$$
将得到的$\widehat{{\beta }_i}$代入回归平面方程得：
$$y={\hat{\beta }}_0+{\hat{\beta }}_1{x}_1+\cdots +{\hat{\beta }}_k{x}_k$$
称为经验回归平面方程，$\widehat{{\beta }_i}$称为经验回归系数

3.1.4 多元线性回归模型和回归系数的检验

1. $F$检验法

   当$H_0$成立时
   $$\begin{gathered} F=\frac{U/k}{Q_e/(n-k-1)}\sim F(k,n-k-1) \\ 其中U=\sum _{i=1}^n({\hat{y}}_i-\bar{y}{)}^2\left(\text{回归平方和}\right),Q_e=\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2(残差平方和) \end{gathered}$$
   如果$F>F_{1-\alpha }(k,n-k-1)$，则拒绝$H_0$,认为$y$与$x_1,\cdots ,x_k$之间显著地有线性关系；否则就接受$H_0$，认为$y$与$x_1,\cdots ,x_k$之间线性关系不显著。

2. $r$检验法

   定义：
   $$R=\sqrt{\frac{U}{L_{yy}}}=\sqrt{\frac{U}{U+Q_e}}$$
   称$y$与$x_1,...,x_k$的多元相关系数或复相关系数，由于
   $$F=\frac{n-k-1}{k}\frac{R^2}{1-R^2}$$
   故$F$与$R$的检验是等效的

3.1.5 多元线性回归模型的预测

1. 点预测

   求出回归方程$\hat{y}={\hat{\beta }}_0+{\hat{\beta }}_1{x}_1+\cdots +{\hat{\beta }}_k{x}_k$, 对于给定自变量的值$x_1^{\ast },\cdots ,x_k^{\ast }$, 用${\hat{y}}^{\ast }={\hat{\beta }}_0+{\hat{\beta }}_1{x}_1^{\ast }+\cdots +$
   ${\hat{\beta }}_k{x}_k^{\ast }$来预测$y^{\ast }={\beta }_0+{\beta }_1{x}_1^{\ast }+\cdots +{\beta }_k{x}_k^{\ast }+ϵ$,称${\hat{\gamma }}^{\ast }$为$y^{\ast }$的点预测

2. 区间预测

   $y$的$1-\alpha$的预测区间(置信)区间为(${\hat{y}}_1,{\hat{y}}_2)$, 其中
   $$\{\begin{array}{l}{\hat{y}}_1=\hat{y}-{\hat{\sigma }}_e\sqrt{1+{\sum }_{i=0}^k{\sum }_{j=0}^k{c}_{ij}{x}_i{x}_j{t}_{1-\frac{\alpha }{2}}(n-k-1)}\\ {\hat{y}}_2=\hat{y}+{\hat{\sigma }}_e\sqrt{1+{\sum }_{i=0}^k{\sum }_{j=0}^k{c}_{ij}{x}_i{x}_j{t}_{1-\frac{\alpha }{2}}(n-k-1)}\end{array}$$

   $$C=L^{-1}=\left(c_{ij}\right),L=X^{\prime }X$$

3.2 使用matlab实现

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

• 输入变量

  • $X,Y$——自变量和因变量的样本值
  • $alpha$——显著性水平，默认为0.05

• 输出变量

  • $b$——回归系数

  • $bint$——回归系数的区间估计

  • $r$——残差

  • $rint$——置信区间

  • $stats$——用于检验回归模型的统计量

    $stats$有四个数值：决定系数$R^2$、$F$值、与$F$对应的概率$P$、无偏估计$\sigma$

$$b=\left[\begin{array}{c}{\hat{\beta }}_0\\ {\hat{\beta }}_1\\ ...\\ {\hat{\beta }}_p\end{array}\right],Y=\left[\begin{array}{c}{y}_1\\ y_2\\ ...\\ y_n\end{array}\right],X=\left[\begin{array}{ccccc}1& x_{11}& x_{12}& \cdots & x_{1p}\\ 1& x_{21}& x_{22}& \cdots & x_{2p}\\ ...& ...& ...& ...& ...\\ 1& x_{n1}& x_{n2}& \cdots & x_{np}\end{array}\right]$$

3.3 多元线性回归中的逐步回归

“最优”的回归方程就是包含所有对$Y$有影响的变量，而不包含对Y影响不显者的变量回归方程
逐步回归分析法的思想：

• 从一个自变量开始，根据自变量对$Y$作用的显著程度，从大到小地依次逐个引入回归方程。
• 当引入的自变量由于后引入变量而变得不显著时，要将其剔除掉。
• 引入一个自变量或从回归方程中剔除一个自变量，为逐步回归的一步。
• 对于每一步都要进行$Y$值检验，以确保每次引入新的显著性变量前回归方程中只包含对$Y$作用显著的变量。
• 这个过程反复进行，直至既无不显著的变量从回归方程中剔除，又无显著变量可引入回归方程时为止。

逐步回归的matlab实现：

$stepwise(x,y,inmodel,alpha)$

• $x$——自变量数据，$n\times m$阶矩阵
• $y$——因变量数据，$n\times 1$阶矩阵
• $inmodel$——矩阵的列数的指标，给出初始模型中包含的子集（缺省时设定为全部自变量)
• $alpha$——默认为0.05

3.4 典型例题

例1：

某建材公司对某年20个地区的建材销售量$Y$(千方)、推销开支、实际帐目数、同类商品竞争数和地区销售潜力分别进行了统计。试分析推销开支、实际帐目数、同类商品竞争数和地区销售潜力对建材销售量的影响作用。试建立回归模型，且分析哪些是主要的影响因素。

序号	推销开支	实际账目数	同类商品竞争数	地区销售潜力	建材销售量
1	5.5	31	10	8	79.3
2	2.5	55	8	6	200.1
3	8.0	67	12	9	163.2
4	3.0	50	7	16	200.1
5	3.0	38	8	15	146.0
6	2.9	71	12	17	177.7
7	8.0	30	12	8	30.9
8	9.0	56	5	10	291.9
9	4.0	42	8	4	160.0
10	6.5	73	5	16	339.4
11	5.5	60	11	7	159.6
12	5.0	44	12	12	86.3
13	6.0	50	6	6	237.5
14	5.0	39	10	4	107.2
15	3.5	55	10	4	155.0
16	8.0	70	6	14	201.4
17	6.0	40	11	6	100.2
18	4.0	50	11	8	135.8
19	7.5	62	9	13	223.3
20	7.0	59	9	11	195.0

设：

• 推销开支$--x_1$
• 实际账目数——$x_2$
• 同类商品竞争数——$x_3$
• 地区销售潜力一一$x_4$

clc
clear
% 1. 数据输入
x1 = [5.5, 2.5, 8.0, 3.0, 3.0, 2.9, 8.0, 9.0, 4.0, 6.5, 5.5, 5.0, 6.0, 5.0, 3.5, 8.0, 6.0, 4.0, 7.5, 7.0]';
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.0, 177.7, 30.9, 291.9, 160.0, 339.4, 159.6, 86.3, 237.5, 107.2, 155.0, 201.4, 100.2, 135.8, 223.3, 195.0]';
X=[ones(size(x1)),x1,x2,x3,x4];

% 2. 求结果
[b,bint,r,rint,stats]=regress(y,X)

% 3. 画残差图
rcoplot(r,rint)

输出：

b =

  191.9158
   -0.7719
    3.1725
  -19.6811
   -0.4501


bint =

  103.1071  280.7245
   -7.1445    5.6007
    2.0640    4.2809
  -25.1651  -14.1972
   -3.7284    2.8283


r =

   -6.3045
   -4.2215
    5.1293
   -3.1536
    0.0469
    6.6035
  -10.2401
   32.1803
   -2.8219
   26.5194
    1.2254
    0.2303
   12.3801
   -5.9704
  -10.0875
  -82.0245
    5.2103
    8.4422
   23.4625
    3.3938


rint =

  -56.5085   43.8995
  -51.7960   43.3530
  -40.4734   50.7319
  -49.7055   43.3983
  -44.3616   44.4554
  -33.2319   46.4388
  -53.1505   32.6704
   -9.9247   74.2853
  -52.7697   47.1259
  -17.9495   70.9883
  -50.3051   52.7559
  -50.3940   50.8545
  -37.1038   61.8641
  -57.7043   45.7635
  -58.2382   38.0632
  -94.9983  -69.0506
  -47.1746   57.5953
  -44.4834   61.3677
  -26.9939   73.9190
  -50.3179   57.1054


stats =

    0.9034   35.0509    0.0000  644.6510

得到的残差图：

例2：

水泥凝固时放出的热量y与水泥中4种化学成分$x_1,x_2,x_3,x_4$有关，今测得一组数据如下，试用逐步回归法确定一个线性模型

序号	1	2	3	4	5	6	7	8	9	10	11	12	13
$x_1$	7	1	11	11	7	11	3	1	2	21	1	11	10
$x_2$	26	29	56	31	52	55	71	31	54	47	40	66	68
$x_3$	6	15	8	8	6	9	17	22	18	4	23	9	8
$x_4$	60	52	20	47	33	22	6	44	22	26	34	12	12
$y$	78.5	74.3	104.3	87.6	95.9	109.2	102.7	72.5	93.1	115.9	83.8	113.3	109.4

使用逐步回归法确定模型

clc
clear
% 1. 数据输入
x1 = [7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10]';
x2 = [26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68]';
x3 = [6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8]';
x4 = [60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12]';
y = [78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1, 115.9, 83.8, 113.3, 109.4]';
x=[x1 x2 x3 x4];
% 2. 逐步回归
stepwise(x,y)

1. 先在初始模型中取全部自变量

   

2. 点击直线3和直线4来移除变量$x_3$和变量$x_4$

   

进行线性回归

% 3. 线性回归
X=[ones(13,1) x1 x2];
b=regress(y,X)

输出：

b =

   52.5773
    1.4683
    0.6623

绘图：

% 4. 作图
% 定义自变量的范围
x1=-10:0.1:10;
x2=-10:0.1:10;
% 计算对应的因变量值
[X1,X2]=meshgrid(x1,x2); % 创建网格点坐标矩阵
Z=52.5773+1.4683*X1+0.6623*X2;
% 绘制三维图形
figure
surf(X1,X2,Z);
xlabel('x1')
ylabel('x2')
zlabel('z')
title('Z=52.5773+1.4683*X1+0.6623*X2')

输出：

4. 非线性回归分析

4.1 模型原理

直线关系是两变量向最简单的一种关系，曲线回归分析的基本任务是通过两个相关变量$x$与$y$的实际观测数据建立曲线回归方程，以揭示$x$与$y$间的曲线联系的形式。曲线回归分析最困难和首要的工作是确定自变量与因变量间的曲线关系的类型，曲线回归分析的基本过程：

• 先将$x$与$\nu$进行变量转换
• 对新变量进行直线回归分析、建立直线回归方程并进行显著性检验和区间估计
• 将新变量还原为原变量，由新变量的直线回归方程和置信区间得出原变量的曲线回归方程和置信区间

由于曲线回归模型种类繁多，所以没有通用的回归方程可直接使用，但是对某些特殊的回归模型，可以通过变量代换、取对数等方法将其线性化，然后使用标准方程求解参数，再将参数带回原方程就是所求。

如有一组数据$x,y$数据，画出散点图后显示$x$与$y$的变动关系为一条递减的双曲线
设
$$\hat{y}=\hat{a}+\hat{b}\frac{1}{x}$$
令$\frac{1}{x}=x^{\prime }$,原式变为
$$\hat{y}=\hat{a}+\hat{b}{x}^{\prime }$$
把原始数据转换为$\frac{1}{x},y$， 求解出$\hat{a},\hat{b}$，再把$x^{\prime }=\frac{1}{x}$带入回去得到曲线回归方程

多重共线性：

回归模型中两个或两个以上的自变量彼此相关的现象，多重共线性带来的问题有：

• 回归系数估计值的不稳定增强
• 回归系数假设检验的结果不显著等

多重共线性检验的主要方法：

• 容忍度：
  $$To{l}_i=1-R^2$$

  • $R_i$是解释变量$x_i$与方程中其他解释变量间的负相关系数
  • 容忍度在0~1之间，越接近与0，表示多重共线性越强，越接近于1，表示多重共线性越弱

• 方差膨胀因子(VIF):
  $$VI{F}_i=\frac{1}{1-R_i^2}$$

  • $VI{F}_i$越大，特别是大于等于10，说明解释变量$x_{\mathrm{i}}$与方程中其他解释变量之间有严重的多重共线心
  • $VI{F}_i$越接近1，表明解释变量$x_{\mathrm{i}}$和其他解释变量之间的多重共线性越弱

4.2 配曲线求解

配曲线的一般方法是：

1. 画出散点图
2. 根据散点图确定须配曲线的类型
3. 估计其中的未知参数

4.2.1 常见的六类曲线

1. 双曲函数曲线：变形双曲线

   方程形式：
   $$\{\begin{array}{l}\hat{y}=\frac{x}{a+bx}\\ \hat{y}=\frac{a+bx}{x}\\ \hat{y}=\frac{1}{a+bx}\end{array}$$
   图形：

   
   • $\hat{y}=\frac{x}{a+bx}$该曲线通过原点(0,0)， 当$a>0$、$b>0$时，$y$随$x$的增大而增大，但速率趋小，并向$y=\frac{1}{b}$ 渐进，是凸曲线；当 $a>0$、$b<0$时，$y$随$x$的增大而增大，速率趋大，并向$x=-\frac{a}{b}$渐进，是凹曲线
   • 变换方式：
     • $\hat{y}=\frac{x}{a+bx}$，两边取倒数后，令$y^{\prime }=\frac{x}{\hat{y}}$,得$y^{\prime }=ax+b$
     • $\hat{y}=\frac{a+bx}{x}$，令$y^{\prime }=\hat{y}x$,得$y^{\prime }=ax+b$ ; $\hat{y}=\frac{1}{a+bx}$，两边取倒数后，令$y^{\prime }=\frac{1}{y}$,得${\mathrm{y}}^{\prime }=ax+b$

2. 幂函数曲线

   方程形式：
   $$\hat{y}=a{x}^b$$
   图形：

   
   • 当 $a>0$、$b>1$时，$y$随$x$的增大而增大(增长),是凹曲线；
   • 当$a>0$、$0<b<1$时，$y$随$x$的增大而增大(增长),但变化缓慢，是凸曲线；
   • 当 $a>0$、$b<0$时，$y$随$x$的增大而减小，且以$x,y$轴为渐近线，是凹曲线。
   • 变换方式：两边取对数，令$y^{\prime }=ln\hat{y},x^{\prime }=lnx,a^{\prime }=lna$，得$y^{\prime }=a^{\prime }+b{x}^{\prime }$

3. 指数函数曲线

   方程形式：
   $$\{\begin{array}{c}\hat{y}=a{e}^{bx}\\ \hat{y}=a{b}^x\end{array}$$
   图形：

   
   • $\hat{y}=a{e}^{bx}$
     • 参数$b$一般用来描述增长或衰减的速度。
     • 当 $a>0、b>0$时，$y$随$x$的增大而增大(增长),是凹曲线；
     • 当 $a>0、b<0$时，$y$随$x$的增大而减小(衰减),是凹曲线。
     • 变换方式：两边取对数，令$y^{\prime }=ln\hat{y},a^{\prime }=lna$,得$y^{\prime }=a^{\prime }+bx$

4. 倒指数曲线

   方程形式：
   $$\hat{y}=a{e}^{\frac{b}{x}},\text{其中}a>0$$
   图形：

   
   • 当 $a>0、b>0$时，$y$随$x$的增大而减小(衰减),是凹曲线；
   • 当 $a>0$、$b<0$时，$y$随$x$的增大而增大(增长),是先凹后凸曲线。
   • 变换方式：两边取对数，令$y^{\prime }=ln\hat{y},a^{\prime }=lna,x^{\prime }=\frac{1}{r}$ ,得$y^{\prime }=a^{\prime }+b{x}^{\prime }$

5. 对数函数曲线

   方程形式：
   $$\hat{y}=a+blnx(x>0)$$
   图形：

   
   • 对数函数表示：$x$变数的较大变化可引起$y$变数的较小变化。
     • $b>0$时$,y$随$x$的增大而增大，是凸曲线；
     • $b<0$时，$y$随$x$的增大而减小，是凹曲线。
   • 变换方式：令$x^{\prime }=lnx$,得$\hat{y}=a+b{x}^{\prime }$

6. $S$型曲线

   • 主要描述动、植物的自然生长过程，又称生长曲线，也可以描述传染病的发展趋势。
   • 生长过程的基本特点是开始增长较慢，而在以后的某一范围内迅速增长，达到一定的限度后增长又缓慢下来，曲线呈拉长的’$S^{\prime }$型曲线。著名的${}^{\prime }{S}^{\prime }$型曲线是Logistic生长曲线。

   方程形式：
   $$\begin{gathered} \hat{y}=\frac{k}{1+a{e}^{-bx}}(a,b,k\text{均大于}0) \\ x=0,\hat{y}=\frac{k}{1+a};x\to \infty ,\hat{y}=k \end{gathered}$$
   图形：

   
   • 所以时间为0的起始量为$\frac{1}{1+a}$，时间为无限延长的，终极量为$k$，曲线$x=\frac{lna}{b}$时有一个拐点，这时$\hat{y}=\frac{k}{2}$，恰好是终极量的一半拐点左侧，为凹曲线，速率由小趋大；拐点右侧，为凸曲线速率由大趋小。
   • 变换方式:两边取倒数再取对数后,$y^{\prime }=ln(\frac{k-\hat{y}}{\hat{y}}),a^{\prime }=lna$，得$y^{\prime }=a^{\prime }+bx$

4.2.2 例题求解

出钢时所用的盛钢水的钢包，由于钢水对耐火材料的侵蚀，容积不断增大。我们希望知道使用次数与增大的容积之间的关系，对一钢包作试验，测得的数据列于下表：

使用次数	增大容积	使用次数	增大容积
1	/	9	9.99
2	6.42	10	10.49
3	8.20	11	10.59
4	9.58	12	10.60
5	9.50	13	10.80
6	9.70	14	10.60
7	10.00	15	10.90
8	9.93	16	10.76

画出散点图：

由散点图我们可以选配倒指数曲线
$$y=a{e}^{\frac{b}{x}}$$
取对数后得
$$lny=lna+\frac{b}{x}lne=lna+\frac{b}{x}$$
设
$${\hat{y}}^{\prime }=lny,x^{\prime }=\frac{1}{x},\hat{b}=lna,\hat{A}=b$$
得
$$\hat{y}=\hat{b}+\hat{A}{x}^{\prime }$$
利用线性回归的方法，算得
$$\begin{array}{rl}& \hat{A}=-1.1107,\hat{b}=2.4587\\ & a=e^{\hat{b}}=11.6789\end{array}$$
所以回归模型为：
$$y=11.6789e^{-\frac{1.1107}{x}}$$

4.3 多项式回归求解

当六类曲线都无法很好地拟合模型时，我们可以使用多项式回归。

设变量$x$、$Y$的回归模型为
$$Y={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+\cdots +{\beta }_p{x}^p+ϵ$$
其中$p$是已知的，${\beta }_i(i=1,2,\cdots ,p)$是未知参数，$ϵ$服从正态分布$N(0,\sigma )$
$$Y={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+\cdots +{\beta }_p{x}^p$$
称为回归多项式，上面的回归模型称为多项式回归。

4.3.1 一元多项式回归

matlab求解方法：
$$y=a_1{x}^m+a_2{x}^{m-1}+\cdots +a_{m}x+a_{m+1}$$

1. 回归：

   1. 确定多项式系数的命令 ：
      $$[p,S]=polyfit(x,y,m)$$
      其中$x=(x_1,x_2,...,x_n),y=(y_1,y_2,...,y_n);p=(a_1,a_2,\cdots ,a_{m+1})$是多项式$y=a_1{x}^m+a_2{x}^{m-1}+\cdots +a_{m}x+a_{m+1}$的系数；$S$是一个结构数据，用来估计预测误差

   2. 一元多项式回归命令：$polytool(x,y,m)$

2. 预测和预测误差估计：

   1. $Y=polyval(p,x)$ 求$polyfit$所得的回归多项式在$x$处预测值$Y;$

   2. $[Y,DELTA]=polyconf(p,x,S,alpha)$

      求$polyfit$所得的回归多项式在$x$处的预测值$Y$及预测值的显著性为$alpha$的置信区间
      $Y\pm DELTA,alpha$缺省时为0.05

示例：

观测物体降落的距离s与时间$t$的关系，得到数据如下表，求s关于$t$的回归方程

$$\hat{s}=a+bt+c{t}^2$$

$t(s)$	1/30	2/30	3/30	4/30	5/30	6/30	7/30
$s(cm)$	11.86	15.67	20.60	26.69	33.71	41.93	51.13

$t(s)$	8/30	9/30	10/30	11/30	12/30	13/30	14/30
$s(cm)$	61.49	72.90	85.44	99.08	113.77	129.54	146.48

方法一：直接作多项式回归

clc
clear
t=1/30:1/30:14/30;
s = [11.86, 15.67, 20.60, 26.69, 33.71, 41.93, 51.13, 61.49, 72.90, 85.44, 99.08, 113.77, 129.54, 146.48];
[p,S]=polyfit(t,s,2);% 进行二次多项式拟合

% 预测及作图
Y=polyconf(p,t,S); % 计算拟合曲线的预测值
plot(t,s,'k+',t,Y,'r'); % 绘制原始数据点和拟合曲线
xlabel('t')
ylabel('s')
title('二次多项式回归')
legend("原始数据","拟合曲线")

结果：

方法二：化为多元线性回归

clc
clear
t=1/30:1/30:14/30;
s = [11.86, 15.67, 20.60, 26.69, 33.71, 41.93, 51.13, 61.49, 72.90, 85.44, 99.08, 113.77, 129.54, 146.48];
% 令x1=t,x2=t^2
T=[ones(14,1),t',(t.^2)'];
[b,bint,r,rint,stats]=regress(s',T);
b

输出：

b =

    9.1329
   65.8896
  489.2946

4.3.2 多元二项式回归

matlab求解方法：

命令：$rstool(x,y,mode{l}^{\prime },alpha)$

$x$——$n$(样本个数)$\times m$(自变量个数)矩阵

$y$——$n$(样本个数)维向量

$alpha$——显著性水平（默认0.05）

$model$——由下列4个模型中选择1个（用字符串输入，缺省时为线性模型）
$$\begin{array}{rl}& linear(\text{线性}):y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_m{x}_m\\ & purequadratic(\text{纯二次}):y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_m{x}_m+\sum _{j=1}^n{\beta }_{jj}{x}_j^2\\ & interaction(\text{交叉}):y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_m{x}_m+\sum _{1\le j\ne k\le m}{\beta }_{jk}{x}_j{x}_k\\ & quadratic(\text{完全二次}):y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_m{x}_m+\sum _{1\le j,k\le m}{\beta }_{jk}{x}_j{x}_k\end{array}$$
示例：

设某商品的需求量与消费者的平均收入、商品价格的统计数据如下，建立回归模型，预测平均收入为1000、价格为6时的商品需求量

需求量	100	75	80	70	50	65	90	100	110	60
收入	1000	600	1200	500	300	400	1300	1100	1300	300
价格	5	7	6	6	8	7	5	4	3	9

选择纯二次模型，即$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_{11}{x}_1^2+{\beta }_{22}{x}_2^2$

x1 = [1000, 600, 1200, 500, 300, 400, 1300, 1100, 1300, 300];
x2= [5, 7, 6, 6, 8, 7, 5, 4, 3, 9];
y= [100, 75, 80, 70, 50, 65, 90, 100, 110, 60];
x=[x1',x2'];
rstool(x,y,'purequadratic')
beta,rmse

选择$x1$和$x2$的值进行预测

输出：

beta =

  110.5313
    0.1464
  -26.5709
   -0.0001
    1.8475


rmse =

    4.5362

方法二：化为多元线性回归

x1 = [1000, 600, 1200, 500, 300, 400, 1300, 1100, 1300, 300];
x2= [5, 7, 6, 6, 8, 7, 5, 4, 3, 9];
y= [100, 75, 80, 70, 50, 65, 90, 100, 110, 60]';
% 令x3=x1^2; x4=x2^2
X=[ones(10,1),x1',x2',(x1.^2)',(x2.^2)'];
[b,bint,r,rint,stats]=regress(y,X);
b,stats

结果：

b =

  110.5313
    0.1464
  -26.5709
   -0.0001
    1.8475


stats =

    0.9702   40.6656    0.0005   20.5771

4.4 一般求解

1. 回归命令

   • 确定回归系数的命令：

     $[beta,r,J]=nlinfit(x,y{,}^{\prime }mode{l}^{\prime },beta0)$

     $beta$——估计出的回归系数

     $r$——残差

     $J$——$Jacobian$矩阵
     $x,y$——输入数据$x,y$分别为$n\times m$矩阵和$n$维列向量，对一元非线性回归，$x$为$n$维列向量
     ${}^{\prime }mode{l}^{\prime }$——是事先用$m$文件定义的非线性函数
     $beta0$——回归系数的初值

   • 非线性回归命令：$nlintool(x,y{,}^{\prime }mode{l}^{\prime },beta0,alpha)$

2. 预测和预测误差估计
   $[Y,DELTA]=nlpredci{(}^{\prime }mode{l}^{\prime },x,beta,r,J)$
   求$nlinfit$或$nlintool$所得的回归函数在$x$处的预测值$Y$及预测值的显著性为$alpha$的置信区间$ Y\pm DELTA $

示例：

出钢时所用的盛钢水的钢包，由于钢水对耐火材料的侵蚀，容积不断增大。我们希望知道使用次数与增大的容积之间的关系，对一钢包作试验，测得的数据列于下表：

使用次数	增大容积	使用次数	增大容积
1	/	9	9.99
2	6.42	10	10.49
3	8.20	11	10.59
4	9.58	12	10.60
5	9.50	13	10.80
6	9.70	14	10.60
7	10.00	15	10.90
8	9.93	16	10.76

由散点图我们可以选配倒指数曲线$y=a{e}^{\frac{b}{x}}$。

对将要拟合的非线性模型$y=a{e}^{\frac{b}{x}}$，建立$volum.m$文件

function yhat=volum(beta,x)
yhat=beta(1)*exp(beta(2)./x);
end

主函数如下

% 输入数据
x=2:16;
y = [ 6.42, 8.20, 9.58, 9.50, 9.70, 10.00, 9.93, 9.99, 10.49, 10.59, 10.60, 10.80, 10.60, 10.90, 10.76];
beta0=[8 2]';
% 求回归系数
[beta,r,j]=nlinfit(x',y','volum',beta0);
beta
% 预测及作图
[YY delta]=nlpredci('volum',x',beta,r,j);
plot(x,y,'k+',x,YY,'r')

输出：

beta =

   11.6037
   -1.0641