统计方法5-回归分析

统计方法5回归分析

前面我们讲过曲线拟合问题。曲线拟合问题的特点是，根据得到的假设干有关变量的一组数据，寻

找因变量与（一个或几个）自变审之间的一个函数，使这个函数对那组数据拟合得最好。通常，函数的

形式可以由经验、先验知识或对数据的直观观察决定，要作的工作是由数据用最小二乘法计算函数中的

待定系数。从计算的角度看•，问题似乎已经完全解决了，还有进一步研究的必要吗？

从数理统计的观点看，这里涉及的都是随机变量，我们根据一个样本计算出的那些系数，只是它们

的一个（点）估计，应该对它们作区间估计或假设检验，如果置信区间太大，甚至包含了零点，那么系

数的估计值是没有多大意义的。另外也应该对模型的误差进行分析，对拟合的优劣给出评价。从建模的

角度说，回归分析就是对拟合问题作的统计分析。

具体地说，回归分析在一组数据的根底上研究这样几个问题：

（i）建立因变量y与自变量凡,当，…，4之间的回归模型（经验公式）；

（ii）对回归模型的可信度进行检验；

（iii）判断每个自变量为&=1,2,对），的影响是否显著；

（iv）诊断回归模型是否适合这组数据；

（v）利用回归模型对y进行预报或控制。

§1多元线性回归

回归分析中最简单的形式是》二用+男］，乂），均为标量，用,四为回归系数，称一元线性回归。

它的一个自然推广是人为多元变量，形如

y=A）+夕内+…+4/⑴

m>2,或者更一般地

J=00+4/（X）+…+PJ,n⑴⑵

其中X=（项,…,（），4（j=1,…,"2）是函数。这里），对回归系数4二（夕0,4,…,4）是线性的，称

为多元线性回归。不难看出，对自变量X作变量代换，就可将（2）化为（1）的形式，所以下面以（1）

为多元线性回归的标准型。

1.1模型

在回归分析中自变量工=（8,占，…，七”）是影响因变量）的主要因素，是人们能控制或能观察的，

而），还受到随机因素的干扰，可以合理地假设这种干扰服从零均值的正态分布，于是模型记作

'）'=&+夕内+..•+»〃川+£一

⑶

卜〜刈。,/）

其中未知。现得到〃个独立观测数据（州,孙，…,与〃），i=由（3）得

y,=&+回/+…+笈/油+巧

(4)

弓~阳0面），i=1,•••,/?

记

Y=i(5)

儿

£=困...%]"p=(A)A…4J

(4)表为

[Y=X3+£

0(6)

£~N(O,b)

1.2参数估计

用最小二乘法估计模型(3)中的参数夕。

由(4)式这组数据的误差平方和为

Q(P)=力婷=(y-Xfi)T(Y-X/3)(7)

1=1

求6使。(6)最小，得到£的最小二乘估计，记作方，可以推出

^=(XTX)-]XTY(8)

将分代回原模型得到y的估计值

§=Ba+Bx\+…+瓦/⑼

而这组数据的拟合值为｝;二乂坂，拟合误差c=y-声称为残差，可作为随机误差£的估计，而

1=1/=|

为残差平方和(或剩余平方和)，即。(历。

1.3统计分析

不加证明地给出以下结果：

⑴6是夕的线性无偏最小方差估计。指的是6是y的线性函数；分的期望等于一；在夕的线性

无偏估计中，/H勺方差最小。

(ii)力服从正态分布

/〜N(p,b2(x7"X广)(11)

(iii)对残差平方和。，EQ=(n-m-\)a2,且

乌~/2(〃一，〃一1)(12)

b

由此得到。'的无偏估计

(13)

n一〃l1

是剩余方差(残差的方差)，S称为剩余标准差。

(iv)对总平方和S=£(K-»)2进行分解，有

»=|

S=Q+U,U=£(少一JOz(14)

r»l

其中。是由(10)定义的残差平方和，反映随机误差对)，的影响，U称为回归平方和，反映自变量对y

的影响。

1.4回归模型的假设检验

因变量y与自变量内,…,/之间是否存在如模型(1)所示的线性关系是需要检验的，显然，如果

所有的1夕1"=1,•••,，〃)都很小，》与西，…,4的线性关系就不明显，所以可令原假设为

当Ho成立时由分解式(14)定义的U,Q满足

厂Ulm~,、

F------------------r(m,n-m-\)(15)

QKn-m-\)

在显著性水平a下有1一。分位数几—假设尸—1),接受名；否则，

拒绝。

注意拒绝"°只说明y与玉,…一%的线性关系不明显，可能存在非线性关系，如平方关系。

还有一些衡量y与玉，…,.二相关程度的指标，如用回归平方和在总平方和中的比值定义

R'=—(16)

S

Ae[0,l]称为相关系数，R越大，y与阳，…,七”相关关系越密切，通常，R大于0.81或0.9)才认为

相关关系成立。

1.5回归系数的假设检验和区间估计

当上面的“0被拒绝时，乩・不全为零，但是不排除其中假没干个等于零。所以应进一步作如下切个

检验()=1,：

"产血一。

由⑴）式，呢是（X,X）T对角线上的元素，用一代替,由⑴）~（13）

式，当”『成立时

t.=/JN”~«〃一m一1）（17）

JQ/（〃一〃?_1）

对给定的。，假设a（〃-加一1）,接受“P；否则，拒绝。

1——

2

（⑺式也可用于对乩.作区间估计（/=（）』,…,〃1）,在置信水平1-aF,乩•的置信区间为

\Pj（18）

22

其中s=J—2—。

\n-in-\

1.6利用回归模型进行预测

当回归模型和系数通过检验后，可由给定的…预测光,凡是随机的，显然其预测

值（点估计）为

九=A+A%、+•••+瓦/〃，（⑼

给定a可以算出方的预测区间（区间估计），结果较复杂，但当〃较大且均接近平均值用时，儿的预

测区间可简化为

[丸-〃（20）

I——1——

22

其中"a是标准正态分布的1-q分位数。

1-32

24

对方的区间估计方法可■用于给出数据残差6=x-y,（z=1,­••,/?）的置信区间，号服从均值为零的

正态分布，所以假设某个q的置信区间不包含零点，那么认为这个数据是异常的，可予以剔除。

1.7Matlab实现

Matlab统计工具箱用命令regress实现多元线性回归，用的方法是最小二乘法，用法是：

b=regress（Y,X）

其中Y.X为按（5）式排列的数据，b为回归系数估计值氐,瓦…,瓦。

[b,bint,r,rint,stats]=regress（Y,X,alpha）

这里Y,X同上，alpha为显著性水平（缺省时设定为0.05）,b,bint为回归系数估计值和它们的置信区间，

rjint为残差（向量）及其置信区间，stats是用于检验回归模型的统计量，有三个数值，第一个是（见

（16）式），第二个是尸（见（⑶式），第3个是与产对应的概率〃，〃＜。拒绝乩），回归模型成立。

注关于R2的说明：一般，R?越接近1,回归方程越显著。

残差及其置信区间可以用rcoplot(r,rint)画图。

1.一元线形回归

例1合金的强度y与其中的碳含量x有比拟密切的关系，今从生产中收集了一批数据如下表:

A|0.100.110.120.130.140.150.160.170.18

I42.041.545.045.545.047.549.055.050.0

试先拟合一个函数y(x)，再用回归分析对它进行检验。

解先画出散点图：

x=0.1:0.01:0.18;

y=(42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0];

plot(x,y,'+1)

可知y与x大致上为线性关系。

设回归模型为

)'=Bo+脏(21)

用regress和rcoplot编程如b：

clc,clear

xl=[0.1:0.01:0.18],;

7=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]

x=[ones省原始数据左边加一列1,表示模型包含常数项3

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

b,bint,stats,rcoplot(rzrint)

得到

b=27.4722137.5000

bint=18.685136.2594

75.7755199.2245

stats=0.798527.74690.0012

即3=27.4722，=137.5(X)0,氐的置信区间是[18.6851,36.25941,8的置信区间是

[75.7755,199.2245]；R2=0.7985,F=27.7469,”=0.0012“

可知模型(21)成立。

Ei回归直线

yhat=x*b;%得到y的预测值％

plot(x1,yhat,'linewidth',3);

ResidualCaseOrderPlot

6

4

2

s

rno

p0

s-

a

y

-2

-4

-6

123456789

CaseNumber

观察命令”外1。1任,由［1）所画的残差分布，假设残差的置信区间不包括。点，那么该点可视为异常点。

除第8个数据外其余残差的置信区间均包含零点，第8个点应视为异常点，将其剔除后重新计算，可得

b=30.7820109.3985

bint=26.280535.2834

76.9014141.8955

stats=0.918867.85340.0002

应该用修改后的这个结果。

两条回归直线的比拟

其中红色虚线为调整之后敢回归直线

2.多元线性回归

例2某厂生产的一种电器的销售量），与竞争对手的价格阳和本厂的价格/有关。下表是该商品

在10，、城市的销售记录。

xI元120140190130155175125145180150

元10011090150210150250270300250

y个10210012077469326696585

试根据这些数据建立y与》和々的关系式，对得到的模型和系数进行检验。假设某市本厂产品售价160

〔元），竞争对手售价170（元），预测商品在该市的销售量。

解分别画出y关于M和y关于9的散点图，可以看出了与占有较明显的线性关系，而）'与阳之

间的关系那么难以确定，我们将作几种尝试，用统计分析决定优劣。

%在MATLAB里进行旋转，可以观察2维的情况％

设回归模型为

y=A）+P\x\+P1X2122）

编写如下程序：

xl=[120140190130155175125145180150]

x2=[10011090150210150250270300250]

y=[10210012077469326696585]

x=[ones(10,1)zxl,x2];

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

b,bint,stats

得到

b=66.51760.4139-0.2698

bint=-32.5060165.5411

-0.20181.0296

-0.4611-0.0785

stats=0.65276.57860.0247

可以看出结果不是太好：p=0.0247,取a=0.05时回归模型(22)可用，但取a=0.01那么模

型不能用；肥=06527较小；片,6的置信区间包含了零点。下面将试图用为的二次函数改良它。

1.8多项式回归

如果从数据的散点图上发现)，与x呈较明显的二次(或高次)函数关系，或者用线性模型(1)的

效果不太好，就可以选用多项式回归。

1.8.1一元多项式回归

一元多项式回归可用命令polyfit实现。

p=polyfit(x,y,n)

[P/S]=polyfit(x,y,n)

[p,S,mu]=polyfit(x,y,n),这里n是多项式阶数，p是多项式系数

例3将17至29岁的运发动每两岁一组分为7组，每组两人测量其旋转定向能力，以考察年龄对

这种运动能力的影响。现得到一组数据如下表：

年龄17192123252729

第一人20.4825.1326.1530.026.120.319.35

第二人24.3528.1126.331.426.9225.721.3

试建立二者之间的关系。

解数据的散点图明显地呈现两端低中间高的形状，所以应拟合一条二次曲线。

选用二次模型

2

y=a2x+tz1x+6zo(23)

编写如下程序：

x0=17:2:29;x0=[xO,x0];

y0=[20.4825.1326.1530.026.120.319.35...

24.3528.1126.331.426.9225.721.3];

[p,s]=polyfit(x0,y0z2);p

得到

p=-0.20038.9782-72.2150

即。2=-0.2003,ax=8.9782,a0=-72.2150。

上面的s是一个数据结构，用于计算函数值，如

[y,delta]=polyconf(p,x0,s);y

得至ij)，的拟合值，及预测值)，的置信区间半径delta。

Dogrco

VoJjBS

XValues

用polylool（x0,y0,2）,可以得到一个如上图的交互式画面，在画面中绿色曲线为拟合曲线，它两侧的

红线是）的置信区间。你可以用鼠标移动图中的十字线来改变图下方的x值，也可以在窗口内输入，左

边就给出y的预测值及其置信区间。通过左下方的Export下拉式菜单，可以输出回归系数等。这个命令

的用法与下面将介绍的rstool相似。

多元二项式回归

统计工具箱提供了一个作多元二项式回归的命令rstool,它也产生一个交互式画面，并输出有关信息,

用法是

rstool（x,y,model,alpha）

其中输入数据x,y分别为〃x〃?矩阵和〃维向量，alpha为显著性水平。（缺省时设定为0.05）,model由以

下4个模型中选择1个（用字符串输入，缺省时设定为线性模型）：

linear（线性）：y=戊）+夕西+…+Pmxm

purequadratic（纯二次）：），=&+2内+…+/3mxm+2夕力引

7=|

interaction（交叉）：丁=4+夕内+…+以£”+Z%/汽

14jwk4m

quadratic（完全二次）：），二4+尸内+…++^Jkxjxk

ISj,k"i

我们再作一遍例2商品销售量与价格问题，选择纯二次模型，即

y=Bo+Ax+P1X2+A1X12+P11A（24）

编程如下：

xl=[120140190130155175125145180150]';

x2=[10011090150210150250270300250]1;

y=[10210012077469326696585]1;

x=[xlx2];

rstool（x,y,*purequadratic'）

300

200

PfAriirtpriY

7n2R11

100

♦/-

51msfi

0

-100

Export▼

PureQuad▼

Close

得到一个如下图的交互式画面，左边是芮（=151）固定时的曲线y（再）及其置信区间,右边是当

（=188）固定时的曲线），（々）及其置信区间。用鼠标移动图中的十字线，或在图下方窗口内输入，可改

变用“2。图左边给出y的预测值及其置信区间，就用这种画面可以答复例2提出的“假设某市本厂产品

售价160（元），竞争对手售价170（元），预测该市的销售量”问题。

图的左下方有两个下拉式菜单，一个菜单Export用以向Mailab工作区传送数据，包括beia（回归系数）,

rmse（剩余标准差），residuals（残差）。模型（24）的回归系数却剩余标准差为

beta=-312.58717.2701-1.7337-0.02280.0037

rinse=16.6436

另一个菜单model用以在上述4个模型中选择，你可以比拟一下它们的剩余标准差，会发现以模型

（24）的rmse=16.6436最小。

§2非线性回归和逐步回归

本节介绍怎样用Matlab统计工具箱实现非线性回归和逐步回归。

2.1非线性回归

非线性回归是指因变量y对回归系数四，…，凡，（而不是自变量）是非线性的。Matlab统计工具箱

中的nlinfit,nlparci*nlpredci,nlintool,不仅给出拟合的回归系数，而且可以给出它的置信

区间，及预测值和置信区间等。

b=nlinfit（X,y,funzb0）

fun是函数的形式，bO是叵归系数初值，b为返回的回归系数估计值。

下面通过例题说明这些命令的用法。

例4在研究化学动力学反响过程中，建立了一个反响速度和反响物含量的数学模型，形式为

一&

y=

1+P\X\+P1X2+B3X3

其中4,…,凡是未知的参数，X”%，当是三种反响物（氢，“戊烷，异构戊烷）的含量，）是反响速

度。今测得一组数据如下表，试由此确定参数口，…，风，并给出其置信区间。笈，…,网的参考值为（°」，

0.05,0.02,1,2）。

序号反响速度y氢修n戊烷x2异构戊烷.

18.5547030010

23.792858010

34.8247030()120

40.0247080120

52.754708010

614.3910019010

72.541008065

84.3547019065

913.0010030054

1()8.5()1(X)30()120

110.0510080120

12113228530010

133/p>

解首先，以回归系数和自变量为输入变量，将要拟合的模型写成函数文件huaxue.m：

functionyhat=huaxue(beta,x);

yhat=(beta(4)*x(:z2)-x(:,3)/beta(5))./(1+beta(1)*x(:,1)+...

beLa(2)*x(:,2)+beLa(3)*x(:,3));

然后，用nlinfit计算回归系数,用nlparci计算回归系数的置信区间,用nlpredci计算预测值及其

置信区间，编程如下：

clc,clear

x0=[18.5547330010

23.792858010

34.82479300120

40.0247380120

52.754738010

614.3910019010

72.541038065

84.3547319065

913.0010030054

108.50100300120

110.0510080120

1211.3228530010

133;

x=x0(:,3:5);

y=x0(:,2)/•

beta=[0.1,0.05,0.02,1,2],；考回归系数的初值

[betahat,f,j]=nlinfit(xz,‘huaxue',beta);当f,j是下面命令用的信息

b2taci=nlparci(betahat,rj)；

betaa=[betahat,betaci]学回归系数及其置信区间

[yhat,delta]=nlpredci('huaxue',x,betahat,f,j)

考y的预测值及其置信区间的半径,置信区间为yhat士delta。

用nlintool得到一个交互式画面，左下方的Export可向工作区传送数据，如剩余标准差等。使用命

令

nlintool(x,y,'huaxue',beta)

可看到画面,并传出剩余标准差rmso=0.1933,

2.2逐步回归

实际问题中影响因变量的因素可能很多，我们希望从中习诞出影响显著的自变量来建立回归模型，

这就涉及到变量选择的问题，逐步回归是一种从众多变量中有效地选择重要变量的方法。以下只讨论线

性回归模型(1)式的情况。

变量选择的标准，简单地说就是所有对因变量影响显著的变量都应选入模型，而影响不显著的变量

都不应选入模型，从便于应用的角度应使模型中变量个数尽可能少。

假设候选的自变量集合为S={2,…从中选出一个子集S|uS,设，中有/个自变量

(/=1,•••,〃?)，由S1和因变量)，构造的回归模型的误差平方和为Q,那么模型的剩余标准差的平方

s2=Q,〃为数据样本容量。所选子集S应使s尽量小，通常回归模型中包含的自变量越多，误

n-l-\

差平方和。越小，但假设模型中包含有对y影响很小的变量，那么。不会由于包含这些变量在内而减

少多少，却因/的增加可能使s反而增大，同时这些对y影响不显著的变量也会影响模型的稳定性，因

此可将剩余标准差s最小作为衡量变量选择的一个数量标准。

逐步回归是实现变鼠选择的一种方法，根本思路为，先确定一初始子集，然后每次从子集外影响显

著的变量中引入一个对)，影响最大的，再对原来子集中的变量进行检验，从变得不显著的变量中剔除一

个影响最小的，直到不能引入和剔除为止。使用逐步回归有两点值得注意，一是要适当地选定引入变量

的显著性水平6”和剔除变量的显著性水平。⑼，显然，。讪越大，引入的变量越多；1H越大，剔除的

变量越少。二是由于各个变量之间的相关性，一个新的变量引入后，会使原来认为显著的某个变量变得

不显著，从而被剔除，所以在最初选择变量时应尽量选择相互独立性强的那些。

在Uatlab统计工具箱中用作逐步回归的是命令stepwise,它提供了一个交互式画面，通过这个工具

你可以自由地选择变量，进行统计分析，其通常用法是：

stepwise(x,y,inmodel,alpha)

其中x是自变量数据，y是因变量数据,分别为和，x1矩阵，inmodel是矩阵x的列数的指标，给出

初始模型中包括的子集(缺省时设定为空)，alpha为显著性水平。

StepwiseRegression窗口，显示回归系数及其置信区间.和其它一些统计量的信息。绿色说明在

模型中的变量，红色说明从模型中移去的变量。在这个窗口中有Export按钮，点击Export产生一个菜单,

说明了要传送给Matlab工作区的参数，它们给出了统计计算的一些结果。

下面通过一个例子说明stepwise的用法。

例5水泥凝固时放出的热量)，与水泥中4种化学成分王,工2，与，匕有关，今测得一组数据如下，试

用逐步回归来确定一个线性模型

序号内x?七％y

172666078.5

2129155274.3

31156820104.3

4113184787.6

575263395.9

61155922109.2

7371176102.7

9254182293.1

102147426115.9

11140233483.8

121166912113.3

131068812109.4

编写程序如下：

clc,clear

x0=[l72666078.5

2129155274.3

31156820104.3

4113184787.6

575263395.9

61155922109.2

7371176102.7

8131224472.5

9254182293.1

102147426115.9

11140233483.8

121166912113.3

131068812109.4];

x=x0(:,2:5);

y=x0(:,6);

stepwise(x,y)

得到如以卜.图形界面：

CoefficientswithErrorBarsCoeff.t-statp-val

Nextstep.

XI2.08270.0708

1.5511MoveX3out

X20.5101680.70490.5009

X30.1019090.13500.8959

X4-0.144061-0.20320.8441

-2-10123

Intercept=62.4054R-square=0.982376F=111.479

RMSE=2.44601AdjR-sq=0.97136p=4.75618e-007

ModelHistory

LLSI

W

H

可以看出,

七,匕不显著，移去这两个变量后的统计结果如下:

CoefficientswithErrorBarsCoeff.t-statp-val

Nextstep.

1•

♦1.4683112.10470.0000Movenoterms

NextStep

一0.6622514.44240.0000

AllSteps

0.2500181.35360.2089