第6章概率统计方法模型上

1、第6章概率统计方法模型在对实际问题进行数学建模的过程中，人们经常遇到随机性的不确定问题，用传统的数学建模方法难以解决。此时，就需要基于概率论和数理统计知识，运用概率统计的方法建立数学模型，对实际问题进行求解，揭示事物发展的基本规律。本章详细介绍用概率统计方法建模的基本思路，结合实际的案例，指出如何用随机变量和概率分布来描述随机不确定事件，说明求解概率统计类模型的一般过程，并指出该类数学模型在社会调查、影响因素分析、发展趋势模拟等方面的广泛应用。§6.1概率模型与Monte Carlo模拟概率模型（1）传染病随机模型在各种传染病的流行过程中，无论健康人还是病人，任何两个人之间接触的机会

2、都是随机的，而且当健康人与病人接触时，健康人是否被传染也是一个随机的事件。我们通过建立传染病随机模型来分析这些随机规律。假设人群总的规模为N，在总人群中，病人的数量为m，健康人的数量为s，即满足N=m+s。在人们的日常生活中，任意两人之间（包括健康人和病人）接触的概率相同，每人平均与k个人接触。当健康人和病人接触时，被传染的概率为p。在以上假设的参数中，m和s通常是已知的，k和p可以通过专家的经验和统计数据获得。我们分析的目的是寻找健康人群中每天平均被感染的人数与已知参数之间的关系，以及初始参数对传染病的扩散速度和流行趋势的影响。我们首先以每一名健康人为研究对象，探讨其每天被感染的概率，而每一

3、名健康人被一名指定病人接触并传染的概率等于每名健康人与指定传染者接触的概率乘以接触时感染的概率。记人群中任意两人接触的概率为q，则对每一名健康人来说，其每天接触的人数服从二项分布，分布函数为,（6.1.1）这个分布的期望为k，即，进而。这样，一名健康人被一名指定病人接触并感染的概率为.进一步，对人群中的每一名健康人来说，其每天不被感染的概率为，被感染的概率为. （6.1.2）所以，对人群中的所有健康人来说，每天被感染的人数服从二项分布，分布函数为, （6.1.3）每天被感染的人数期望为，标准差为为了得到简明的结果，对进行近似计算，由于通常人群的总数，且根据Talyor展开，得，因此，. （6.

4、1.4）通过式（6.1.4）可以看出平均每天被感染的人数与s、m、p和k之间的关系。进而可以度量平均每天被感染人数的相对误差即 （6.1.5）由式（6.1.4）可以看出，对于健康人群来说，每天平均被感染的人数与人群中每人每天平均接触的人数k，健康人与病人接触时被感染的概率p成正比。当n，p，k都确定的情况下，时，也就是在整个人群中，病人和健康人的数量各占一半时，每天被感染的人数达到最大。为了对传染病的传染过程有一个直观的了解，假设一个人口总量n=10000的人群，在日常生活中，平均每人每天接触的人数k=18，健康人与病人接触时被感染的概率p=10%，对于不同的m，平均每天被感染人数与相对误差的

5、变化趋势如所示。可见被感染人数随着病人数量的增大而增大，直到病人数量占总人群数量的一半时达到最大，随后呈下降趋势。随着病人人口的增加每天被感染人数的相对误差一直呈减少趋势，尤为明显的是病人数量增长的前期，相对误差急剧减少。图6.1.1 平均每天被感染人数的趋势图 平均每天被感染人数的相对误差趋势R编程如下：crb <- function(m, n=10000, p=0.1, k=18) #函数 u<-(m*(n-m)*p*k)/(n-1);u m<-1:10000plot(1:10000,crb(m),xlab="m",ylab="平均每天被的传

6、染人数",type="l", col="blue") crb1 <- function(s, n=10000, p=0.1, k=18) #相对误差函数 miugama<-(n-1-m*p*k)/(n-m)*m*p*k)0.5;miugama m<-1:6000plot(1:6000,crb1(m),xlab="m",ylab="相对误差",type="l", col="red") （2）企鹅繁殖模型企鹅的繁殖过程是一个典型的随机不确定模型。首先，

7、每只母企鹅下蛋的数量是随机的，服从泊松分布，其次，每个企鹅蛋是否可以成功孵化也是不确定的。针对这一问题，我们在合理假设的基础上，建立概率模型，求企鹅后代个数的期望值。根据人们的统计，企鹅生蛋的个数是服从参数为的泊松分布，即 （6.1.6）而每个生蛋能发育成企鹅的概率为p，且每个生蛋能否发育成企鹅是彼此独立的随机事件。令代表企鹅后代的个数，且有。可见取非负的整数值0, 1, 2, ，对于的概率，我们利用全概率公式 (6.1.7)注意到每个生蛋发育成小企鹅是相互独立的，且发育成小企鹅的概率为p，因此，实际上反映了有k个生蛋，每个生蛋独立发育，恰好发育成k个企鹅的概率。显然，它是一个伯努利试验，因而

8、 (6.1.8)于是有得到企鹅后代的个数服从参数为的泊松分布，从而企鹅后代个数的期望值为。也说明了企鹅后代的个数与生蛋的个数以及发育成功的概率成正比。Monte Carlo模拟Monte Carlo（蒙特卡洛）模拟，也称统计模拟方法。该方法是上世纪40年代，由John von Neumann（冯·诺依曼），Stanislaw Ulam和Nicholas Metropolis在洛斯阿拉莫斯国家实验室进行核武器计划的工作时发明的，后来该方法的得名是由于Ulam的叔叔常在驰名世界的赌城（摩纳哥的Monte Carlo）输钱。事实上，Monte Carlo模拟是由于科学技术的发展和电子计算机

9、的发明，而被提出的一种以概率统计理论为指导的非常重要的数值计算方法。该方法是一种使用随机数来解决很多计算问题的方法。目前，蒙特卡罗方法在金融工程学，宏观经济学，生物医学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。通常蒙特卡罗方法可以粗略地分成两类：一类是所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接模拟这种随机的过程。另一种类型是所求解问题可以转化为某种随机分布的特征数，比如随机事件出现的概率，或者随机变量的期望值。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于

10、求解复杂的多维积分问题。例如，我们要计算一个不规则图形的面积，蒙特卡罗方法基于这样的思想：假想你有一袋豆子，把豆子均匀地朝这个图形上撒，然后数这个图形之中有多少颗豆子，这个豆子的数目就是图形的面积。当你的豆子越小，撒的越多的时候，结果就越精确。借助计算机程序可以生成大量均匀分布坐标点，然后统计出图形内的点数，通过它们占总点数的比例和坐标点生成范围的面积就可以求出图形面积。可以看出，Monte Carlo得到概率模型的解是通过试验得到的，而不是计算出来的。也正是由于这个原因，对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，Monte Carlo方法是一种有效的求出数值解的方法。

11、我们利用Monte Carlo模拟的方法实现对圆周率的估计。考虑边长为1的正方形，1为半径的四分之一圆弧，如图6.1.1所示。 Monte Carlo对的估计在边长为1的正方形内，等概率的产生n个随机点，。这样和就是（0,1）上均匀分布的随机数。当n个点中有k个点落在四分之一圆内，既有k个点满足关系式：，则当时，有如下关系：此时，圆周率的估计值为。通过R语言编程如下：monte<-function(n) k<-0 x<-runif(n) #runif()函数的作用是产生均匀分布的随机数y<-runif(n) for(i in 1:n) if (xi2+yi2<=1

12、)k<-k+1 pi<-4*k/n 其中，runif(n, a, b)的意义是在(a, b)区间上，产生n个均匀分布的随机数。Runif(n) 的意义是在 (0，1) 区间上产生n个均匀分布的随机数，调用monte函数，当n取不同值时，得到不同的的估计值。表6.1.1 的估计值列表n取值300100030005000100003.05333.10803.16803.15683.1404例求定积分，被积函数如下图。图6.1.2 被积函数图像解：因为该积分不能直接求解，这里采用mento carlo模拟来求解，见图图 6.1.3 monte carlo模拟求解定积分在Monte Car

13、lo模拟中，经常需要用到随机数。实际上随机数产生的方法有很多，现以R软件为例，介绍用计算机软件产生随机数的方法：（1）在（a,b）范围内产生n个均匀分布的随机数: runif(n, a, b)。当a,b默认时，为（0, 1）区间上的随机数。例如：>runif(10,2,3) #在(2,3)范围内产生10个随机数1 2.597194 2.998407 2.203209 2.897273 2.403639 2.873541 2.508925 2.8786569 2.003151 2.096483（2）产生n个均值为，标准差为的正态分布随机数：rnorm(n, a, b)。当a,b默认时，为标

14、准正态分布N(0,1)的随机数。例如：> rnorm(30,1,1)1 0.9623346 0.1193753 0.9032579 -0.6649220 1.2346502 0.75963177 -0.2019421 0.7026093 0.2871598 -0.3824490 0.7913398 0.794809213 0.7883549 2.2074363 0.5559556 3.0239753 1.6879051 2.159033719 -0.3332347 1.4590187 1.2827723 -0.9300292 0.1866562 1.163462425 -1.493072

15、3 1.8948205 1.0855420 0.1551674 1.6671424 0.6467984（3）产生n个参数为的随机数：rpois(n, )。例如：> rpois(30,3)1 3 2 5 4 5 1 2 2 2 0 5 2 4 1 4 4 1 6 4 6 3 4 2 7 4 2 1 9 6 2下面通过一个例子说明如何产生具有一定分布律的离散型随机变量的随机数。产生具有分布律0的离散型随机变量X的随机数。解 设是(0,1)上均匀分布的随机数，令.则()是具有随机变量X分布律的随机数。例如产生20个随机数，编程如下：n=100; r<-runif(n

16、); x<-array(0,dim=c(n); for (i in 1:n) if (ri<=0.3) xi<-0 else if (ri<=0.6) xi<-1 else xi<-2 x1 2 1 2 1 0 1 2 2 2 1 0 0 2 2 0 2 0 2 2 2 2 1 2 2 0 1 2 0 1 30 1 0 0 0 2 2 1 0 1 0 0 1 2 2 1 1 0 1 2 2 0 1 2 2 0 2 2 1 1 59 2 2 0 2 1 2 0 0 1 0 1 2 2 0 0 0 2 0 2 2 2 0 2 2 0 0 2 2 0 88 1 2

17、 2 1 1 1 0 2 1 1 0 0 1一列火车从A站开往B站，某人每天赶往B站上火车，他已经了解到火车从A站到B站的运行时间是服从均值为30min，标准差为2min的正态随机变量。火车大约下午13:00离开A站，此人大约13:30达到B站。火车离开A站的时刻及概率及此人到达B站的时刻及概率如下表所示。火车离开时刻13:0013:0513:1012:55概率0.700.200.100人到站时刻13:2813:3013:3213:34概率0.1问他能赶上火车的概率是多少。利用Monte Carlo方法进行分析。以下是求解过程的R程序：MC<-function(n)r1

18、<-runif(n); r2<-runif(n); t2<-rnorm(n,30,2)t1<-array(0,dim=c(1,n); t3<-t1;for(i in 1:n)if (r1i<=0.7)t1i<-0else if (r1i<=0.9)t1i<-5elset1i<-10for(i in 1:n)if (r2i<=0.3)t3i<-28else if (r2i<=0.7)t3i<-30else if (r2i<=0.9)t3i<-32elset3i<-34k<-0for(i i

19、n 1:n) if (t1i+t2i>t3i) k<-k+1 k/n做一万次模拟，得到：MC(10000)1 0.6341此人能赶上火车的概率大约是0.6341。§6.2报童问题与随机库存模型本节应用概率统计知识，首先介绍基本的Newsboy模型，然后将Newsboy模型拓展，讨论随机库存模型。6.2.1报童问题Newsboy问题中，报童每天清晨从报社购进报纸，通过一天的零售后，晚上将没有卖掉的报纸以低于购进价的价格退回。设进价为c，零售价为s，剩余退回的价格为a，问其如何确定每天购进的数量，使其期望获益最大。这里满足。从过上述假设，报童每正常卖掉一份报纸利润为，退回一份

20、赔，由于需求量事先无法确定，是随机的。若通过以往销售的经验了解到需求量的随机规律，销售份的概率为，。我们根据以及报纸的进价、零售价和剩余退回价格来建立优化模型，求解最优的订购量。假设报童早晨购进报纸的量为n，则 或 ，所以每天的收入也是不确定的。这里考虑报童在不同销售情况下，建立每天销售收入的期望函数，则 （6.2.1）接下来求当n为何值时，达到最大？由于r为离散的，这里用差分的方法来求式（6.2.1）的极值。令令，且，则 （6.2.2）也就是说，当，a，s和c具体确定时，n即可确定。例 6.2.1 某服装店出售某款夏季时装。该款衣服成本100元，售价200元。如整个夏季不能售出，则必须降价为

21、70元。设降价后一定可以售出，已知售货量r服从泊松分布为平均出售数，根据以往经验，平均出售数为120件。问该店的订货量应该为多少单位？解：由题意知：s=200，a=70，c=100代入式（6.2.2），可得编程如下：poisson<-function(r) #泊松分布 lamda<-120 #期望 y<-(exp(-lamda)*(lamdar)/(factorial(r); y #寻找n值f<-0f1<-poisson(1)for (i in 1:1000) fi+1<-fi+poisson(i+1) if (fi+1>=(10/13) breakf

22、;i可得且，所以更接近于10/13。故，最佳订购量应该为126件。随机库存模型由于市场对于商品的需求是随机变量，事前难以知道需求的准确数值。此时，无论工厂或商店无法决定存贮策略，从存贮的角度来考虑，假设在一个阶段开始的时刻原有的库存为I，如供应不足则须承担缺货费，如供应有余，则多余的部分仍须存贮起来。有余存在这种不确定性，就需要计算随机变量的期望值，从而定出最佳的存贮量。我们考虑一个时间段落。做下列符号假设：原有的存贮量为I；存贮货物的单价为k；订购一次的订购费为C1，如订货量为Q时，所需要的订货费为；单位货物的存贮费为，缺货费为；需求量为r的概率为。当本阶段开始时，订货量为Q，存储量达到I+

23、Q。则本阶段所需要的各种费用由订货费、存贮费和缺货费构成。订货费：；存贮费：当需求时，未能售出的存贮部分必须付存贮费；时，不需要付存贮费。因此，所需要存贮费的期望值为：。当时，不付存贮费及缺货费。缺货费：当需求时，则会发生缺货现象，必须付缺货费缺货费用的期望值为：综上，在整个阶段所需的订货费、缺货费及存贮费的期望之和为：（6.2.3）为简便起见，记，则式（6.2.3）即为 （6.2.4）求S值使C(S)达到最小。将需求r的随机值按大小顺序排列为：，其中，（）。S只从中取值。当S取值为时，记为，则。与newsboy模型中求极值的方法类似，我们求的最小值。 （6.2.5）记，则(6.2.6) 令，

24、由于，所以，我们有 （6.2.7）式（6.2.7）右端的数值称为临界值，记为。我们选使不等式成立的得最小值为S，则订货量为。模型中还有一个问题需要我们解决，那就是原库存消耗到什么水平时，需要订货？假设这一水平是s，当时，可以不订货，当时要订货，使库存达到S，订货量为。要想确定s，首先需要考察不等式 （6.2.8）因s也只能从中取值，使式（6.2.8）成立的（）值中最小者定为s。当时，式（6.2.8）左端缺货费用的期望值虽然在增加，但订货费及存贮费期望值都在减少。在最不利的情况下，如时，不等式使成立的，因此s值一定存在。例 6.2.2 某汽车零部件生产企业，对某型号钢材的需求量的概率为：需求量（

25、吨）8090100110120P(r)0.30.1已知每吨钢材的购价为k=7500元，订货费为元，存贮费元，缺货费元。求该企业最优的存贮策略。解：（1）临界值另外，且，因此，S=90吨为最优订货量。（2）利用式（6.2.8）求s：由于S=90，式（6.2.8）右端为当s=80时，式（6.2.8）左端为此时，式（6.2.8）成立，故。可知，该企业最优的存贮策略为每当钢材的库存低于80吨，补充存贮使存贮量达到90吨，当存贮量大于80吨时，不需要补充。§6.3线性回归模型回归分析是应用数理统计学研究问题的一种重要的方法模型，它的目的是研究变量之间的相互关系，建立变量之间的

26、经验公式，以便达到预测和控制的目的。6.3.1一元线性回归模型在一元线性回归分析里，我们要考察的是：随机变量y与普通变量x之间的联系。对于x和y，通过观测或实验，得到若干对数据，。表6.3.1 GDP与人均消费支出年份XY年份XY1978675.1359.819901602.3797.11979716.943719911727.2861.41980763.7464.119921949.8966.61981792.4501.919932187.91048.61982851.1533.519942436.11108.71983931.4572.819952663.71213.119841059.2

27、635.619962889.11322.819851185.271619973111.91380.919861269.6746.519983323.11460.619871393.6788.319993529.31564.419881527836.420003789.71690.819891565.9779.7图6.3.1 中国人均GDP与人均消费支出从散点图上发现，观测点基本在一条线附近，从而可以认为Y与X的关系式线性的，既因变量Y主要受自变量X的影响，而这些观测点与直线的偏离都是由其它一些不确定因素造成的。因此，我们做如下假定： （6.3.1）其中，式（6.3.1）被称为总体回归函数，和是

28、未知参数，称为回归系数；表示Y随X的变化而线性变化的部分，是随机误差，称为随机干扰项，反应了未列入方程式的其他一切不确定影响因素对Y影响的总和，通常假定，且随机误差项与自变量X线性无关；称函数为一元线性回归函数，X为自变量，Y为因变量。考虑到，是(X,Y)的一组观测值，则一元线性回归模型可表示为： （6.3.2）式中，。求回归参数的一种思路是要求图6.3.1中的点与直线上的点偏离越小越好，若和是未知参数和的估计值，则被为回归值或拟合值。和的最小二乘估计是指使 （6.3.3）成立，经计算可得：， （6.3.4）式中，由此可得回归方程为。通常取 （6.3.5）为参数的估计量，进一步可证明，为的无偏

29、估计，即。关于和估计的标准差分别为， （6.3.6）因此，根据数理统计学中区间估计的原理，。我们有回归系数的区间估计为：，从回归参数的估计式（6.3.4）可以知道，在回归系数的估计中，不一定要知道Y与X是否有线性关系，但如果不存在这种关系，那么回归方程便毫无意义。因此需要对回归方程进行检验。在统计意义上是E(Y)随X线性变化的变化率，若，则E(Y)实际上并不随X作线性变化，仅当时，一元线性回归方程才有意义。因此假设检验为：.通常用三种检验方法：（1）t检验。当成立时，统计量对于给定的显著性水平，检验的拒绝域为（2）F检验。当成立时，统计量.对于给定的显著性水平，检验的拒绝域为.（3）相关系数检

30、验。记，则称R为样本相关系数，对于给定的显著性水平，查相关系数临界值表可得，则检验的拒绝域为当拒绝时，认为线性回归方程式显著的。另外，进一步介绍回归方程对样本拟合程度的评价，所谓拟合程度，是指样本观测值聚集在样本回归线周围的紧密程度。判断回归模型拟合程度优劣最常用的数量尺度为样本决定系数，它是建立在对总离差平方和进行分解的基础之上的。样本决定系数的公式为：的取值范围为。由公式可以看出当所有的样本点都位于回归直线上时，说明总离差可以完全由所估计的样本回归直线来解释。根据图6.3.1中，人均消费支出与人均GDP的数据，建立回归方程模型，进行参数估计并做相应的检验，编程如下：x<-c(675.

31、1,716.9,763.7,792.4,851.1,931.4,1059.2,1185.2,1269.6,1393.6,1527,1565.9,1602.3,1727.2,1949.8,2187.9,2436.1,2663.7,2889.1,3111.9,3323.1,3529.3,3789.7);y<-c(359.8,437,464.1,501.9,533.5,572.8,635.6,716,746.5,788.3,836.4,779.7,797.1,861.4,966.6,1048.6,1108.7,1213.1,1322.8,1380.9,1460.6,1564.4,1690.8)

32、;regression<-lm(y1+x) #作线性模型summary(regression)运行结果见图6.3.2。图6.3.2 回归结果第一部分（call）列出了相应的回归模型的公式，其中y1+x表示，第二部分（Residuals）列出的是残差的最小值点，1/4分位点，中位数点，3/4分位点和最大值点。在计算结果的第三部分（Coefficients）中，Estimate表示回归方程参数的估计，即和；Std. Error表示回归参数的标准差，即和；t value为t值，即，表示p值，即，并且有显著性程度的标记。在计算结果的第四部分，Residual standard error表示残差

33、的标准差，即式（6.3.5）中的，自由度为n-2。Multiple R-Squared为相关系数的平方，即，F-statistic表示F统计量，即，其自由度为(1, n-2)，p-value为p值，即概率值。从计算结果可以看出，回归方程通过了回归参数的显著性检验和回归方程的检验，因此得到的回归方程为：Y=196+0.3881X.进一步，可以通过R软件求预测值和预测区间。例如求回归方程中X=2000时的预测区间。这里即为求X=2000时的预测值，和置信程度为的置信区间。编程如下：new<-data.frame(x=2000)lm.pred<-predict(regression,ne

34、w,interval="prediction",level=0.95)lm.predfit lwr upr1972.2582 893.0552 1051.461因此当X=2000时，得到相应的预测值为972.2582，预测区间为893.0552, 1051.461。6.3.2多元线性回归模型一元线性回归是一个主要影响因素作为自变量来解释因变量的变化。然而，在现实问题研究中，一种现象常常是与多个因素相联系的，此时就需要用两个或两个以上的影响因素作为自变量来解释因变量的变化，由多个自变量的最优组合共同来预测或估计因变量，比只用一个自变量进行预测或估计更有效，更符合实际。在回归分

35、析中，如果有两个或两个以上的自变量，就称为多元回归或多重回归。在实际的应用中，多元线性回归比一元线性回归用途更广且实用意义更大。在建立多元线性回归模型时，随机变量与一般变量的多元线性回归模型为：其中是个未知参数，称为回归常数项，称为回归系数；称为被解释变量（因变量），是个可以精确测量并可控制的一般变量，称为解释变量（自变量）。为随机扰动项，代表主观或客观原因造成的不可观测的随机误差，它是一个随机变量通常假定满足。则多元线性总体回归方程为：系数表示在其它自变量不变的情况下，自变量变动一个单位时引起的因变量的平均变动单位，其它回归系数的含义类似。（1）样本回归模型建立设，是随机变量与一般变量的n次

36、独立观测值，则此时多元线性模型可表示为： （）其中，独立同分布。多元线性回归样本方程为：,式中为的估计值。为了方便起见，令，则式（6.3.7）可改写为： （6.3.8）且满足，。回归方程可改写为： （6.3.9）多元线性回归方程中同归系数的估计采用最小二乘法。若对于变量记残差平方和为，根据微积分中求极小值原理，可知残差平方和存在最小值，欲使达到最小，用对的偏导数使其值等于零。加以整理后可以得到个方程式：通过求解这一方程组便可求出的估计值。得 （）取为残差向量，取 （6.3.11）为的估计，也称为的最小二乘估计。可以证明：.进一步可以证明的方差估计为：.相应的的标准差为其中是对角线上第i个元素。

37、（2）显著性检验与一元线性回归分析不同，在多元线性回归分析中，很难用图形来判断E(y)是否随作线性变化，因而显著性检验尤为重要。对多元线性回归方程的拟合程度进行测定、检验回归方程和回归系数的显著性。拟合优度检验测定多元线性回归的拟合程度，使用多重判定系数，其定义为：式中SSR为回归平方和，SSE为残差平方和，SST为总离差平方和。当的值范围为，越接近1，回归平面拟合程度越高；反之越接近0，回归平面拟合程度越低。回归方程的显著性检验（F检验）所谓回归方程的显著性检验就是检验假设：所有回归系数都等于零，即检验： ；不全为0。多元线性回归方程的显著性检验一般采用F检验。F统计量的定义为，平均的回归平

38、方和与平均的残差平方和（均方误差）之比，对于多元线性回归方程，在成立的条件下：式中，SSR为回归平方和，SSE为残差平方和，为样本，为自变量个数。F统计量服从的是第一自由度为，第二自由度为的F分布。从F统计量的定义式可看出，如果F值较大，则说明自变量造成的因变量的变动远远大于随机因素对因变量造成的影响。另外，从另一个角度来看，F统计量也可以反映回归方程的拟合优度。将F统计量的公式与的公式作结合转换，可得：可见，如果回归方程的拟合优度高，F统计量就越显著；F统计量越显著，回归方程的拟合优度就越高。利用F统计量进行回归方程显著性检验的步骤总结如下：Step 1提出假设：,不全为0Step 2在成立

39、条件下，计算F统计量由样本观测值计算F值。Step 3根据给定的显著性水平确定临界值，或者计算F值所对应的相伴概率值p。如果（或者），就拒绝原假设，接受备择假设，认为所有回归系数同时与零有显著性差异，自变量与应变量之间存在显著性的线性关系，自变量的变化确实能够反映因变量的线性变化，回归方程显著。如果（或者），则接受原假设，认为所有回归系数同时与零无显著性差异，自变量与应变量之间不存在显著性的线性关系，自变量的变化无法反映因变量的线性变化，回归方程不显著。回归系数显著性检验（t检验）回归方程的显著性检验是对线性回归方程的一个整体性检验。如果我们检验的结果是拒绝原假设，则意味着因变量Y线性地依赖于

40、自变量，这个回归自变量的整体。但是，这并不排除Y并不依赖于其中某些自变量。因此，我们还要对每个自变量逐一做显著性检验，即回归系数的显著性检验。回归系数的显著性检验是检验各自变量对因变量的影响是否显著，从而找出哪些自变量对的影响是重要的，哪些是不重要的。对于多元回归方程，回归系数的显著性检验，即检验假设，在假设成立的条件下，T统计量式中为的对角线上第j个元素。t检验步骤如下：Step 1提出假设；式中，表示零假设，表示备择假设。如果零假设成立，则说明对没有显著性的影响，反之，则说明对有显著性的影响；Step 2 在成立的前提下，计算回归系数的t统计量；Step 3 给定的显著性水平，确定临界值,

41、或者计算t值所对应的相伴率值p的大小。应注意的是，t检验的临界值是由显著性水平和自由度决定的，这里进行的检验是双侧检验，所以临界值为。如果（或者）,就拒绝原假设，接受备择假设，认为回归系数与零有显著性差异，该自变量和应变量之间存在显著的线性关系，它的变动较好地解释说明应变量的变动，应保留在回归方程中；反之，如果（或者）,就接受原假设，认为回归系数与零无显著性差异，该自变量和应变量之间不存在显著的线性关系，它的变动无法较好地解释说明应变量的变动，应剔除回归方程。例近年来，高校招生规模急剧扩大，在教育产业化的背景下，选取1985-2003年相关数据，对影响我国高校招生人数的各因素及其影响程度的大小

42、进行定量分析。表 6.3.2 高校招生人数影响因素的多元线性回归分析原始数据年份高校招生数y国家财政教育经费x1农村家庭平均收入x2年份高校招生数y国家财政教育经费x1农村家庭平均收入x2198546871227.9397.61995510531028.41577.7198641310270.4423.81996593981211.911926.1198739017285.9462.61997637491357.732090.1198835645340.7544.91998725081565.5922162198928569397.7601.51999922251815.762210.3199

43、029649433.9686.320001284842085.6792253.4199129679482.2708.620011651972582.382366.4199233439564.978420022026113114.2382475.6199342145644.4921.620032689253453.862622.21994508648841221数据EXCEL文件如图6.3.3，首先将数据文件存为csv格式。然后编程如下：X<-read.csv("d:program FilesRchengxudata2.csv",header=TRUE)y<-X,

44、2;x1<-X,3;x2<-X,4;lm.sol<-lm(yx1+x2)summary(lm.sol)运行结果见图6.3.4。从运行结果可以看出，回归系数和回归方程的检验都是显著的，因此，回归方程为：.图数据文件图 6.3.4 运行结果（3）回归系数的区间估计由参数的统计性质可知： （）因此，的区间估计为：（4）预测当多元线性回归方程经过检验是显著的，且其中每个系数均显著不为0，则说明回归的结果是合理的，在此基础上可用回归方程作预测。当时，代入回归方程可得.当观测值为时，的置信度为的预测区间为：.例如求回归方程中，当x=(900, 1300)T时的预测区间。这里即为求x=(9

45、00, 1300)T时的预测值，和置信程度为的置信区间。编程如下：new<-data.frame(x1<-900,x2<-1300)lm.pred<-predict( lm.sol, new, interval="prediction", level=0.95)lm.predfit lwr upr1 51668.3 28980.95 74355.66因此，当x=(900, 1300)T时，得到相应的预测值为51668.3，预测区间为28980.95, 74355.66。6.3.3逐步回归模型逐步回归的实质是在建立多元回归方程的过程中，首先按偏相关系数

46、大小次序，将自变量逐个引入方程，并且对引入方程中的每个自变量偏相关系数进行统计检验，效应显著的自变量留在回归方程内。继续遴选下一个自变量，如果效应不显著，停止引入新自变量。由于新自变量的引入，原已引入方程中的自变量由于变量之间的相互作用其效应有可能变的不显著，经统计检验后要随时从方程中剔除，只保留效应显著的自变量，直到不再引入和剔除自变量为止，从而得到最优的回归方程。一般来说，如果在一个回归方程中忽略了对Y有显著影响的自变量，那么所建立的方程必与实际有较大的偏离，但变量选的过多，使用就不方便，特别当方程中含有对Y影响不大的变量时，可能由于残差平方和自由度的减小而使的估计增大，从而影响使用回归方

47、程作预测的精度。因此，适当地选择变量以建立一个“最优”的回归方程是十分重要的。在多元线性逐步回归中，“最优”的含义是指从可供选择的所有变量中选出对Y有显著影响的变量建立方程，且在方程中不含对Y无显著影响的变量。R软件提供了较为方便的“逐步回归”计算函数step( )，它是以信息统计量为准则，通过选择最小的AIC信息统计量，来达到删除或增加变量的目的。我们通过一个例子来说明如何通过R软件来实现逐步回归的过程。某水泥在凝固时放出的热量y（单位：卡/克）与水泥中下列四种化学成分有关：3CaO.Al2O3的成分(%);：3CaO.Si2O3的成分(%);：4CaO.Al2O3.Fe2O3的成分(%);

48、：2CaO.SiO2的成分(%)。所测定的数据如表6.3.3所示。用回归分析建立y与四种化学成分的线性回归模型。表 水泥中所含化学成分表试验序号172666078.52129155274.331156820104.34113184787.6575263395.961155922109.27371176102.78131224472.59254182293.1102147426115.911140233483.8121166912113.3131068812109.4数据来源：首先做多元线性回归如下：x<-read.csv("d:program FilesRchengxuzhubuhuigui.csv",header=F)y<-x,5;x1<-x,1;x2<-x,2;x3<-x,3;x4<-x,4;lms<-lm(yx1+x2+x3+x4)summary(lms)运行结果如下：Call:lm(formula = y x1 + x2 + x3 + x4)Residuals: Min 1Q Median 3Q Max -3.1750 -1.67