第6章 概率统计方法模型

第6章概率统计方法模型在对实际问题进行数学建模的过程中，人们经常遇到随机性的不确定问题，用传统的数学建模方法难以解决。此时，就需要基于概率论和数理统计知识，运用概率统计的方法建立数学模型，对实际问题进行求解，揭示事物发展的基本规律。本章详细介绍用概率统计方法建模的基本思路，结合实际的案例，指出如何用随机变量和概率分布来描述随机不确定事件，说明求解概率统计类模型的一般过程，并指出该类数学模型在社会调查、影响因素分析、发展趋势模拟等方面的广泛应用。§6.1概率模型与MonteCarlo模拟6.1.1概率模型（1）传染病随机模型在各种传染病的流行过程中，无论健康人还是病人，任何两个人之间接触的机会都是随机的，而且当健康人与病人接触时，健康人是否被传染也是一个随机的事件。下面通过建立传染病随机模型来分析这些随机规律。假设人群总的规模为n，在总人群中，病人的数量为m，健康人的数量为s，即满足n=m+s。在人们的日常生活中，任意两人之间（包括健康人和病人）接触的概率相同，每人平均与k个人接触。当健康人和病人接触时，被传染的概率为p。在以上假设的参数中，m和s通常是已知的，k和p可以通过专家的经验和统计数据获得。我们分析的目的是寻找健康人群中每天平均被感染的人数与已知参数之间的关系，以及初始参数对传染病的扩散速度和流行趋势的影响。我们首先以每一名健康人为研究对象，探讨其每天被感染的概率，而每一名健康人被一名指定病人接触并传染的概率等于每名健康人与指定传染者接触的概率乘以接触时感染的概率。记人群中任意两人接触的概率为q，则对每一名健康人来说，其每天接触的人数服从二项分布，分布函数为,（6.1.1）这个分布的期望为k，即，进而。这样，一名健康人被一名指定病人接触并感染的概率为.进一步，对人群中的每一名健康人来说，其每天不被感染的概率为，被感染的概率为.（6.1.2）所以，对人群中的所有健康人来说，每天被感染的人数服从二项分布，分布函数为,（6.1.3）每天被感染的人数期望为，标准差为为了得到简明的结果，对进行近似计算，由于通常人群的总数，且根据Talyor展开，得，因此，.（6.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，平均每天被感染人数与相对误差的变化趋势如图6.1.1和图6.1.2所示。可见被感染人数随着病人数量的增大而增大，直到病人数量占总人群数量的一半时达到最大，随后呈下降趋势。随着病人人口的增加每天被感染人数的相对误差一直呈减少趋势，尤为明显的是病人数量增长的前期，相对误差急剧减少。图6.1.1平均每天被感染人数的趋势图图6.1.2平均每天被感染人数的相对误差趋势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="平均每天被的传染人数",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）企鹅繁殖模型企鹅的繁殖过程是一个典型的随机不确定模型。首先，每只母企鹅下蛋的数量是随机的，服从泊松分布，其次，每个企鹅蛋是否可以成功孵化也是不确定的。针对这一问题，我们在合理假设的基础上，建立概率模型，求企鹅后代个数的期望值。据统计，企鹅生蛋的个数是服从参数为的泊松分布，即（6.1.6）而每个生蛋能发育成企鹅的概率为p，且每个生蛋能否发育成企鹅是彼此独立的随机事件。令代表企鹅后代的个数，且有。可见取非负的整数值0,1,2,…，对于的概率，我们利用全概率公式(6.1.7)注意到每个生蛋发育成小企鹅是相互独立的，且发育成小企鹅的概率为p，因此，实际上反映了有k个生蛋，每个生蛋独立发育，恰好发育成r个企鹅的概率。显然，它是一个伯努利试验，因而(6.1.8)于是有得到企鹅后代的个数服从参数为的泊松分布，从而企鹅后代个数的期望值为。也说明了企鹅后代的个数与生蛋的个数以及发育成功的概率成正比。6.1.2MonteCarlo模拟MonteCarlo（蒙特卡洛）模拟，也称为统计模拟方法。该方法是上世纪40年代，由JohnvonNeumann（冯·诺依曼），StanislawUlam和NicholasMetropolis在洛斯阿拉莫斯国家实验室进行核武器计划的工作时发明的，后来该方法的得名是由于Ulam的叔叔常在驰名世界的赌城（摩纳哥的MonteCarlo）输钱。事实上，MonteCarlo模拟是由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的非常重要的数值计算方法。该方法是一种使用随机数来解决很多计算问题的方法。目前，蒙特卡罗方法在金融工程学，宏观经济学，生物医学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。通常蒙特卡罗方法可以粗略地分成两类：一类是所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接模拟这种随机的过程。另一种类型是所求解问题可以转化为某种随机分布的特征数，比如随机事件出现的概率，或者随机变量的期望值。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于求解复杂的多维积分问题。例如，我们要计算一个不规则图形的面积，蒙特卡罗方法基于这样的思想：假想你有一袋豆子，把豆子均匀地朝这个图形上撒，然后数这个图形之中有多少颗豆子，这个豆子的数目就是图形的面积。当你的豆子越小，撒的越多的时候，结果就越精确。借助计算机程序可以生成大量的均匀分布坐标点，然后统计出图形内的点数，通过它们占总点数的比例和坐标点生成范围的面积就可以求出图形面积。可以看出，MonteCarlo得到概率模型的解是通过试验得到的，而不是计算出来的。也正是由于这个原因，对于那些由于计算过于复杂而难以得到解析解，或者根本没有解析解的问题，MonteCarlo方法是一种有效的求出数值解的方法。我们利用MonteCarlo模拟的方法实现对圆周率的估计。考虑边长为1的正方形，1为半径的四分之一圆弧，如图6.1.1所示。图6.1.1MonteCarlo对的估计在边长为1的正方形内，等概率产生n个随机点，。这样和就是（0,1）上均匀分布的随机数。当n个点中有k个点落在四分之一圆内，即有k个点满足关系式：，则当时，有如下关系：此时，圆周率的估计值为。通过R语言编程如下：monte<-function(n){k<-0x<-runif(n)#runif()函数的作用是产生均匀分布的随机数y<-runif(n)for(iin1:n){if(x[i]^2+y[i]^2<=1)k<-k+1}pi<-4*k/n}其中，runif(n,a,b)的意义是在(a,b)区间上，产生n个均匀分布的随机数。runif(n)的意义是在(0，1)区间上产生n个均匀分布的随机数，调用monte函数，当n取不同值时，得到不同的的估计值见表6.1.1。表6.1.1的估计值列表n取值300100030005000100003.05333.10803.16803.15683.1404例6.1.1求定积分，被积函数如下图。图6.1.2被积函数图像解：因为该积分不能直接求解，这里采用mentocarlo模拟来求解，见下图图6.1.3Montecarlo模拟求解定积分在MonteCarlo模拟中，经常需要用到随机数。实际上随机数产生的方法有很多，现以R软件为例，介绍用计算机软件产生随机数的方法：（1）在（a,b）范围内产生n个均匀分布的随机数:runif(n,a,b)。当a,b默认时，为（0,1）区间上的随机数。例如：>runif(10,2,3)#在(2,3)范围内产生10个随机数[1]2.5971942.9984072.2032092.8972732.4036392.8735412.5089252.878656[9]2.0031512.096483（2）产生n个均值为，标准差为的正态分布随机数：rnorm(n,a,b)。当a,b默认时，为标准正态分布N(0,1)的随机数。例如：>rnorm(30,1,1)[1]0.96233460.11937530.9032579-0.66492201.23465020.7596317[7]-0.20194210.70260930.2871598-0.38244900.79133980.7948092[13]0.78835492.20743630.55595563.02397531.68790512.1590337[19]-0.33323471.45901871.2827723-0.93002920.18665621.1634624[25]-1.49307231.89482051.08554200.15516741.66714240.6467984（3）产生n个参数为的随机数：rpois(n,)。例如：>rpois(30,3)[1]325451222052414416463427421962下面通过一个例子说明如何产生具有一定分布律的离散型随机变量的随机数。例6.1.2产生具有分布律的离散型随机变量X的随机数。解设是(0,1)上均匀分布的随机数，令.则()是具有随机变量X分布律的随机数。例如产生20个随机数，编程如下：n=100;r<-runif(n);x<-array(0,dim=c(n));for(iin1:n){if(r[i]<=0.3)x[i]<-0elseif(r[i]<=0.6)x[i]<-1elsex[i]<-2}x[1]21210122210022020222212201201[30]10002210100122110122012202211[59]22021200101220002022202200220[88]1221110211001例6.1.3一列火车从A站开往B站，某人每天赶往B站上火车，他已经了解到火车从A站到B站的运行时间是服从均值为30min，标准差为2min的正态随机变量。火车大约下午13:00离开A站，此人大约13:30达到B站。火车离开A站的时刻及概率及此人到达B站的时刻及概率如下表所示。表6.1.2火车离开A站和此人到达B站的时刻概率分布表火车离开时刻13:0013:0513:1012:55概率0.700.200.100人到站时刻13:2813:3013:3213:34概率0.1问他能赶上火车的概率是多少。利用MonteCarlo方法进行分析。以下是求解过程的R程序：MC<-function(n){r1<-runif(n);r2<-runif(n);t2<-rnorm(n,30,2)t1<-array(0,dim=c(1,n));t3<-t1;for(iin1:n){if(r1[i]<=0.7){t1[i]<-0}elseif(r1[i]<=0.9){t1[i]<-5}elset1[i]<-10}for(iin1:n){if(r2[i]<=0.3){t3[i]<-28}elseif(r2[i]<=0.7){t3[i]<-30}elseif(r2[i]<=0.9){t3[i]<-32}elset3[i]<-34}k<-0for(iin1:n)if(t1[i]+t2[i]>t3[i])k<-k+1k/n}做一万次模拟，得到：MC(10000)[1]0.6341此人能赶上火车的概率大约是0.6341。§6.2报童问题与随机库存模型本节应用概率统计知识，首先介绍基本的Newsboy模型，然后将Newsboy模型拓展，讨论随机库存模型。6.2.1Newsboy问题中，报童每天清晨从报社购进报纸，通过一天的零售后，晚上将没有卖掉的报纸以低于购进价的价格退回。设进价为c，零售价为s，剩余退回的价格为a，问其如何确定每天购进的数量，使其期望获益最大。这里满足。从过上述假设，报童每正常卖掉一份报纸利润为，退回一份赔，由于需求量事先无法确定，是随机的。若通过以往销售的经验了解到需求量的随机规律，销售份的概率为，。我们根据以及报纸的进价、零售价和剩余退回价格来建立优化模型，求解最优的订购量。假设报童早晨购进报纸的量为n，则或，所以每天的收入也是不确定的。这里考虑报童在不同销售情况下，建立每天销售收入的期望函数，则（6.2.1）接下来求当n为何值时，达到最大。由于r为离散的，这里用差分的方法来求式（6.2.1）的极值。令令，且，则（6.2.2）也就是说，当，a，s和c具体确定时，n即可确定。例6.2.1某服装店出售某款夏季时装。该款衣服成本100元，售价200元。如整个夏季不能售出，则必须降价为70元。设降价后一定可以售出，已知售货量r服从泊松分布为平均出售数，根据以往经验，平均出售数为120件。问该店的订货量应该为多少单位？解：由题意知：s=200，a=70，c=100代入式（6.2.2），可得编程如下：poisson<-function(r)#泊松分布{lamda<-120#期望y<-((exp(-lamda))*(lamda^r))/(factorial(r));y}###寻找n值f<-0f[1]<-poisson(1)for(iin1:1000){f[i+1]<-f[i]+poisson(i+1)if(f[i+1]>=(10/13))break}f;i可得且，所以更接近于10/13。故，最佳订购量应该为126件。6.2.2由于市场对于商品的需求是随机变量，事前难以知道需求的准确数值。此时，无论工厂或商店均无法决定存贮策略，从存贮的角度来考虑，假设在一个阶段开始的时刻原有的库存为I，如供应不足则须承担缺货费，如供应有余，则多余的部分仍须存贮起来。由于存在这种不确定性，就需要计算随机变量的期望值，从而定出最佳的存贮量。我们考虑一个时间段落。做下列符号假设：原有的存贮量为I；存贮货物的单价为k；订购一次的订购费为C1，如订货量为Q时，所需要的订货费为；单位货物的存贮费为，缺货费为；需求量为r的概率为。当本阶段开始时，订货量为Q，存储量达到I+Q。则本阶段所需要的各种费用由订货费、存贮费和缺货费构成。订货费：；存贮费：当需求时，未能售出的存贮部分必须付存贮费；时，不需要付存贮费。因此，所需要存贮费的期望值为：。当时，不付存贮费及缺货费。缺货费：当需求时，则会发生缺货现象，必须付缺货费缺货费用的期望值为：综上，在整个阶段所需的订货费、缺货费及存贮费的期望之和为：（6.2.3）为简便起见，记，则式（6.2.3）即为（6.2.4）求S值使C(S)达到最小。将需求r的随机值按大小顺序排列为：，其中，，（）。S只从中取值。当S取值为时，记为，则。与newsboy模型中求极值的方法类似，我们求的最小值。（6.2.5）记，则(6.2.6)令，由于，所以，我们有（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需求量（吨）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线性回归模型回归分析是应用数理统计学研究问题的一种重要的方法模型，它的目的是研究变量之间的相互关系，建立变量之间的经验公式，以便达到预测和控制的目的。6.3.1在一元线性回归分析里，我们要考察的是：随机变量y与普通变量x之间的联系。对于x和y，通过观测或实验，得到若干对数据，，…，例如，从总体上考查中国居民收入与消费支出的关系。研究中国居民人均国民生产总值X与中国居民人均消费支出Y的关系。观测值见表6.3.1，X与Y观测值的散点图如图6.3.1。表6.3.1GDP与人均消费支出年份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.2635.619962889.11322.819851185.271619973111.91380.919861269.6746.519983323.11460.619871393.6788.319993529.31564.419881527836.420003789.71690.819891565.9779.7从散点图上发现，观测点基本在一条线附近，从而可以认为Y与X的关系是线性的，既因变量Y主要受自变量X的影响，而这些观测点与直线的偏离都是由其它一些不确定因素造成的。因此，我们作如下假定：（6.3.1）其中，式（6.3.1）被称为总体回归函数，和是未知参数，称为回归系数；表示Y随X的变化而线性变化的部分，是随机误差，称为随机干扰项，反应了未列入方程式的其他一切不确定影响因素对Y影响的总和，通常假定，且随机误差项与自变量X线性无关；称函数为一元线性回归函数，X为自变量，Y为因变量。考虑到，，…，图6.3.1中国人均GDP与人均消费支出是(X,Y)的一组观测值，则一元线性回归模型可表示为：（6.3.2）式中，，，，。求回归参数的一种思路是要求图6.3.1中的点与直线上的点偏离越小越好，若和是未知参数和的估计值，则被为回归值或拟合值。和的最小二乘估计是指使（6.3.3）成立，经计算可得：，（6.3.4）式中，，，由此可得回归方程为。通常取（6.3.5）为参数的估计量，进一步可证明，为的无偏估计，即。关于和估计的标准差分别为，（6.3.6）因此，根据数理统计学中区间估计的原理，，。两个回归系数的区间估计值为：，从回归参数的估计式（6.3.4）可以知道，在回归系数的估计中，不一定要知道Y与X是否有线性关系，但如果不存在这种关系，那么回归方程便毫无意义。因此需要对回归方程进行检验。在统计意义上是E(Y)随X线性变化的变化率，若，则E(Y)实际上并不随X作线性变化，仅当时，一元线性回归方程才有意义。因此假设检验为：.通常用三种检验方法：（1）t检验。当成立时，统计量对于给定的显著性水平，检验的拒绝域为。（2）F检验。当成立时，统计量.对于给定的显著性水平，检验的拒绝域为。（3）相关系数检验。记，则称R为样本相关系数，对于给定的显著性水平，查相关系数临界值表可得，则检验的拒绝域为。当拒绝时，认为线性回归方程式显著的。这里，进一步分析回归方程对样本拟合程度的评价，所谓拟合程度，是指样本观测值聚集在样本回归线周围的紧密程度。判断回归模型拟合程度优劣最常用的数量尺度为样本决定系数，它是建立在对总离差平方和进行分解的基础之上的。样本决定系数的公式为：的取值范围为。由公式可以看出当所有的样本点都位于回归直线上时，，说明总离差可以完全由所估计的样本回归直线来解释。根据图6.3.1中，人均消费支出与人均GDP的数据，建立回归方程模型，进行参数估计并做相应的检验，编程如下：x<-c(675.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);regression<-lm(y~1+x)#做线性回归summary(regression)运行结果见图6.3.2。图6.3.2回归结果第一部分（call）列出了相应的回归模型公式，其中y~1+x表示，第二部分（Residuals）列出的是残差的最小值点，1/4分位点，中位数点，3/4分位点和最大值点。在计算结果的第三部分（Coefficients）中，Estimate表示回归方程参数的估计，即和；Std.Error表示回归参数的标准差，即和；tvalue为t值，即，表示p值，即，并且有显著性程度的标记。在计算结果的第四部分，Residualstandarderror表示残差的标准差，即式（6.3.5）中的，自由度为n-2。MultipleR-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,new,interval="prediction",level=0.95)lm.predfitlwrupr1972.2582893.05521051.461可知，当X=2000时，得到相应的预测值为972.2582，预测区间为[893.0552,1051.461]。6.3.2一元线性回归是一个主要影响因素作为自变量来解释因变量的变化。然而，在现实问题研究中，一种现象常常是与多个因素相联系的，此时就需要用两个或两个以上的影响因素作为自变量来解释因变量的变化。由多个自变量的最优组合共同来预测或估计因变量，比只用一个自变量进行预测或估计更有效，更符合实际。在回归分析中，如果有两个或两个以上的自变量，就称为多元回归或多重回归。在实际的应用中，多元线性回归比一元线性回归用途更广且实用意义更大。在建立多元线性回归模型时，随机变量与一般变量的多元线性回归模型为：其中是个未知参数，称为回归常数项，称为回归系数；称为被解释变量（因变量），是个可以精确测量并可控制的一般变量，称为解释变量（自变量）。为随机扰动项，代表主观或客观原因造成的不可观测的随机误差，它是一个随机变量通常假定满足。多元线性总体回归的方程为：系数表示在其它自变量不变的情况下，自变量变动一个单位时引起的因变量的平均变动单位，其它回归系数的含义类似。（1）样本回归模型的建立设，是随机变量与一般变量的n次独立观测值，则此时多元线性模型可表示为：（6.3.7）其中，独立同分布。多元线性回归样本方程为：,式中为的估计值。为方便起见，令，，，，，则式（6.3.7）可改写为：（6.3.8）且满足，。回归方程可改写为：（6.3.9）多元线性回归方程中回归系数的估计采用最小二乘法。记残差平方和为，根据微积分中求极小值原理，可知残差平方和存在最小值，欲使达到最小，用对求偏导数，使其值等于零。加以整理后可以得到个方程式：通过求解这一方程组便可求出的估计值。得（6.3.10）则为残差向量，取（6.3.11）为的估计，也称为的最小二乘估计。可以证明：.进一步可以证明的方差估计为：.相应的的标准差为其中是对角线上第i个元素。（2）显著性检验与一元线性回归分析不同，在多元线性回归分析中，很难用图形来判断E(y)是否随作线性变化，因而显著性检验尤为重要。对多元线性回归方程的拟合程度进行测定、检验回归方程和回归系数的显著性。①拟合优度检验测定多元线性回归的拟合程度，使用多重判定系数，其定义为：式中SSR为回归平方和，SSE为残差平方和，SST为总离差平方和。当的值范围为，越接近1，回归平面拟合程度越高；反之越接近0，回归平面拟合程度越低。②回归方程的显著性检验（F检验）所谓回归方程的显著性检验就是检验假设：所有回归系数都等于零，即检验：；不全为0。多元线性回归方程的显著性检验一般采用F检验。F统计量定义为回归平方和的平均与残差平方和的平均（均方误差）之比，对于多元线性回归方程，在成立的条件下：式中，SSR为回归平方和，SSE为残差平方和，为样本，为自变量个数。F统计量服从的是第一自由度为，第二自由度为的F分布。从F统计量的定义式可看出，如果F值较大，则说明自变量造成的因变量的变动远远大于随机因素对因变量造成的影响。另外，从另一个角度来看，F统计量也可以反映回归方程的拟合优度。将F统计量的公式与的公式作结合转换，可得：可见，如果回归方程的拟合优度高，F统计量就越显著；F统计量越显著，回归方程的拟合优度就越高。利用F统计量进行回归方程显著性检验的步骤总结如下：Step1提出假设：,不全为0Step2在成立条件下，计算F统计量由样本观测值计算F值。Step3根据给定的显著性水平确定临界值，或者计算F值所对应的相伴概率值p。如果（或者），就拒绝原假设，接受备择假设，认为所有回归系数同时与零有显著性差异，自变量与应变量之间存在显著性的线性关系，自变量的变化确实能够反映因变量的线性变化，回归方程显著。如果（或者），则接受原假设，认为所有回归系数同时与零无显著性差异，自变量与应变量之间不存在显著性的线性关系，自变量的变化无法反映因变量的线性变化，回归方程不显著。③回归系数显著性检验（t检验）回归方程的显著性检验是对线性回归方程的一个整体性检验。如果我们检验的结果是拒绝原假设，则意味着因变量Y线性地依赖于自变量，这个回归自变量的整体。但是，这并不排除Y并不依赖于其中某些自变量。因此，我们还要对每个自变量逐一做显著性检验，即回归系数的显著性检验。回归系数的显著性检验是检验各自变量对因变量的影响是否显著，从而找出哪些自变量对的影响是重要的，哪些是不重要的。对于多元回归方程，回归系数的显著性检验，即检验假设，在假设成立的条件下，T统计量式中为的对角线上第j个元素。t检验步骤如下：Step1提出假设；式中，表示零假设，表示备择假设。如果零假设成立，则说明对没有显著性的影响，反之，则说明对有显著性的影响；Step2在成立的前提下，计算回归系数的T统计量；Step3给定的显著性水平，确定临界值,或者计算t值所对应的相伴率值p的大小。应注意的是，t检验的临界值是由显著性水平和自由度决定的，这里进行的检验是双侧检验，所以临界值为。如果（或者）,就拒绝原假设，接受备择假设，认为回归系数与零有显著性差异，该自变量和应变量之间存在显著的线性关系，它的变动较好地解释说明应变量的变动，应保留在回归方程中；反之，如果（或者）,就接受原假设，认为回归系数与零无显著性差异，该自变量和应变量之间不存在显著的线性关系，它的变动无法较好地解释说明应变量的变动，应剔除回归方程。例6.3.1近年来，高等学校招生规模急剧扩大，在教育产业化的背景下，选取1985-2003年相关数据，对影响我国高校招生人数的各因素及其影响程度的大小进行定量分析。表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.3199029649433.9686.320001284842085.6792253.4199129679482.2708.620011651972582.382366.4199233439564.978420022026113114.2382475.6199342145644.4921.620032689253453.862622.21994508648841221数据EXCEL文件如图6.3.3所示，首先将数据文件存为csv格式。然后编程如下：X<-read.csv("d:\\programFiles\\R\\chengxu\\data2.csv",header=TRUE)y<-X[,2];x1<-X[,3];x2<-X[,4];lm.sol<-lm(y~x1+x2)summary(lm.sol)运行结果见图6.3.4。从运行结果可以看出，回归系数和回归方程的检验都是显著的，因此，回归方程为：.图6.3.3数据文件（3）回归系数的区间估计由参数的统计性质可知：（6.3.12）因此，的区间估计为：图6.3.4运行结果（4）预测当多元线性回归方程经过检验是显著的，且其中每个系数均显著不为0，则说明回归的结果是合理的，在此基础上可用回归方程作预测。当时，代入回归方程可得.当观测值为时，的置信度为的预测区间为：.如求例6.3.1回归方程中，当x=(900,1300)T时的预测区间。这里即为求x=(900,1300)T时的预测值，和置信程度为0.95的置信区间。编程如下：new<-data.frame(x1<-900,x2<-1300)lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95)lm.predfitlwrupr151668.328980.9574355.66因此，当x=(900,1300)T时，得到相应的预测值为51668.3，预测区间为[28980.95,74355.66]。6.3.3逐步回归的实质是在建立多元回归方程的过程中，首先按偏相关系数的大小次序，将自变量逐个引入方程，并且对引入方程中的每个自变量偏相关系数进行统计检验，效应显著的自变量留在回归方程内。继续遴选下一个自变量，如果效应不显著，停止引入新自变量。由于新自变量的引入，原已引入方程中的自变量由于变量之间的相互作用其效应有可能变的不显著，经统计检验后要随时从方程中剔除，只保留效应显著的自变量，直到不再引入和剔除自变量为止，从而得到最优的回归方程。一般来说，如果在一个回归方程中忽略了对Y有显著影响的自变量，那么所建立的方程必与实际有较大的偏离，但变量选的过多，使用就不方便，特别当方程中含有对Y影响不大的变量时，可能由于残差平方和自由度的减小而使的估计增大，从而影响使用回归方程作预测的精度。因此，适当地选择变量以建立一个“最优”的回归方程是十分重要的。在多元线性逐步回归中，“最优”的含义是指从可供选择的所有变量中选出对Y有显著影响的变量建立方程，且在方程中不含对Y无显著影响的变量。R软件提供了较为方便的“逐步回归”计算函数step()，它是以信息统计量为准则，通过选择最小的AIC信息统计量，来达到删除或增加变量的目的。我们通过一个例子来说明如何通过R软件来实现逐步回归的过程。例6.3.2某水泥在凝固时放出的热量y（单位：卡/克）与水泥中下列四种化学成分有关：：3CaO.Al2O3的成分(%);：3CaO.Si2O3的成分(%);：4CaO.Al2O3.Fe2O3的成分(%);：2CaO.SiO2的成分(%)。所测定的数据如表6.3.3所示。用回归分析建立y与四种化学成分的线性回归模型。表6.3.3水泥中所含化学成分表试验序号172666078.52129155274.331156820104.34113184787.6575263395.961155922109.27371176102.78131224472.59254182293.1102147426115.911140233483.8121166912113.3131068812109.4首先做多元线性回归如下：x<-read.csv("d:\\programFiles\\R\\chengxu\\zhubuhuigui.csv",header=F)y<-x[,5];x1<-x[,1];x2<-x[,2];x3<-x[,3];x4<-x[,4];lms<-lm(y~x1+x2+x3+x4)summary(lms)运行结果如下：Call:lm(formula=y~x1+x2+x3+x4)Residuals:Min1QMedian3QMax-3.1750-1.67090.25081.37833.9254Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)62.405470.07100.8910.3991x11.55110.74482.0830.0708.x20.51020.72380.7050.5009x30.10190.75470.1350.8959x4-0.14410.7091-0.2030.8441---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:2.446on8degreesoffreedomMultipleR-squared:0.9824,AdjustedR-squared:0.9736F-statistic:111.5on4and8DF,p-value:4.756e-07从回归方程的计算可以看到，如果选择全部变量作回归方程，效果是不好的，因为回归方程的系数只有一个通过了检验。在R软件中，step()函数时逐步回归函数，它是以AIC信息统计量为准则，通过删除或增加变量，使能得到的“最优”回归方程的AIC值达到最小。lm2<-step(lms)显示如下结果：Start:AIC=26.94y~x1+x2+x3+x4DfSumofSqRSSAIC-x310.109147.97324.974-x410.247048.11125.011-x212.972550.83625.728<none>47.86426.944-x1125.950973.81530.576Step:AIC=24.97y~x1+x2+x4DfSumofSqRSSAIC<none>47.9724.974-x419.9357.9025.420-x2126.7974.7628.742-x11820.91868.8860.629从结果可以看出，用全部变量做回归时，AIC值为26.94；去掉变量时，AIC变为24.974；去掉变量时，AIC值变为25.011；去掉变量时，AIC值为25.728；去掉变量时，AIC值为30.576。因此，R软件自动去掉变量，进行下一轮计算。下一轮计算中，无论去掉哪个变量后，AIC值均会升高，因此，R软件终止计算，得到“最优”回归方程。下面分析一下计算过程，见图6.3.5。图6.3.5逐步回归后的结果分析由逐步回归后的结果可以看出，回归系数检验的显著性水平有很大提高，但变量的回归系数检验显著性水平仍不理想。下面该如何处理呢？在R软件中，还有两个函数可以用来作逐步回归，这两个函数是add1()和drop1()。>drop1(lm2)SingletermdeletionsModel:y~x1+x2+x4DfSumofSqRSSAIC<none>47.9724.974x11820.91868.8860.629x2126.7974.7628.742x419.9357.9025.420从运算变量来看，如果删去变量，AIC的值会从24.97增加到25.42，是增加的最少的。另外，除AIC准则外，残差的平方和也是逐步回归的重要指标之一，一般来说，拟合越好的方程，残差的平方和越小。去掉变量，残差平方和上升9.93，也是最少的。因此，从这两项指标来看，应该再去掉变量。图6.3.6逐步回归后的最终最后的回归结果见图6.3.6，从最后的回归结果可以看出，该结果的回归方程和回归系数都通过了检验。因此，逐步回归的结果为：.§6.4非线性回归模型非线性回归模型按变量个数也可以分为一元非线性回归模型和多元非线性回归模型。曲线的形式也因实际情况不同而有多种形式，如指数曲线、双曲线、S形曲线等。下面我们列出几类典型的非线性回归模型的函数形式：（1）双曲线模型：（6.4.1）（2）多项式模型：（6.4.2）（3）对数模型：（6.4.3）（4）三角函数模型：（6.4.4）（5）指数模型：（6.4.5）（6.4.6）（6）幂函数模型：（6.4.7）我们将上述非线性回归模型分为两类来处理：第一类：直接换元型。这类非线性回归模型通过简单的变量代换可直接转化为线性回归模型，如式（6.4.1）、式（6.4.2）、式（6.4.3）和式（6.4.4）。第二类：间接代换型。这类非线性回归模型通过对数变形代换可间接地转化为线性回归模型，如：式（6.4.5）、式（6.46）和式（6.对于式（6.1.1）、式（6.1.2）、式（6.1.3）和式（6.1.4）所示的非线性回归模型，虽然包含有非线性变量，但因变量与待估计系数之间的关系却是线性的。对于此类模型，可以直接通过变量代换将其化为线性模型，具体代换方法见表6.4.1。表6.4.1变量代换表原模型模型代换代换后模型参数估计一元线性回归OLS法多元线性回归OLS法一元线性回归OLS法一元线性回归OLS法对于式（6.4.5）、式（6.4.6）和式（6.4.7）所示的非线性回归模型，因变量与待估计参数之间的关系也是非线性的。因此不能通过直接换元化为线性模型。对此类模型，可通过对回归方程两边取对数转换为可以直接换元的形式。这种先取对数再进行变量代换的方法称为间接换元法。为使取对数后回归方程的形式更为简捷，我们不妨将式（6.4.5）和式（6.4.7）中随机扰动项的形式进行变换，将式（6.4.5）和式（6.4（6.4.5‘）（6.4.7‘）对（6.1.5‘）、式（6.1.6）和式（6.1.7’）两边取对数，得（6.4.8）（6.4.9）（6.4.10）式（6.4.8）、式（6.4.9）和式（6.4.10）皆可经过适当的换元直接转化为线性回归方程，通过线性回归的方法来进行参数估计。下面，我们来研究不能通过上述两种方法来处理的非线性回归模型。设非线性回归模型具有如下形式：（6.4.11）其中，。设()是的n次独立观测值，则多元非线性模型(6.4.11)可表示为（6.4.12）其中，且独立同分布。为方便起见，将式（6.4.12）简写为，其中，。为求参数的估计值，转化为求解最小二乘问题（6.4.13）式（6.4.13）的解作为参数的估计值。可以证明，的最小二乘估计也是其最大似然估计。在R软件中，一般通过函数nls()求解非线性最小二乘问题，下面通过例子来说明求解过程。例6.4.1在化学工业的可靠性研究中，对象是某种产品A。在对产品进行制造的过程中，单位产品中必须含有0.50的有效氯气，产品中的氯气随着时间的增加而减少，在产品到达用户之前的最初8周内，氯气含量衰减到0.49。但由于随后出现了许多无法控制的因素，因而在后8周理论的计算对有效氯气的进一步预报是不可靠的。为有利于管理，需要决定产品中所含的有效氯气随时间的变化规律。在一段时间中观测若干盒产品得到的数据如表6.4.1。假定非线性模型：利用非线性最小二乘法进行参数估计。表6.4.1单位产品有效氯气百分数序号生产后时间有效氯气序号生产后时间有效氯气180.4910260.412100.4811270.403120.4612280.404140.4313300.415160.4414320.406180.4615340.407200.4216360.428220.4117380.389240.4218400.39R编程如下：>data<-data.frame(+x=c(8,10,12,14,16,18,20,22,24,26,27,28,30,32,34,36,38,40),+y=c(0.49,0.48,0.46,0.43,0.44,0.46,0.42,0.41,0.42,0.41,0.40,0.40,0.41,0.40,0.40,0.42,0.38,0.39)+)>nls.sol<-nls(y~a+(0.49-a)*exp(-b*(x-8)),data=data,start=list(a=0.1,b=0.01))>summary(nls.sol)Formula:y~a+(0.49-a)*exp(-b*(x-8))Parameters:EstimateStd.ErrortvaluePr(>|t|)a0.386680.0109935.172<2e-16***b0.082470.021393.8550.0014**---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:0.01196on16degreesoffreedomNumberofiterationstoconvergence:19Achievedconvergencetolerance:1.996e-06从计算结果看出，参数回归结果通过了显著性检验，得函数表达式如下：用R函数编程，画拟合函数与样本点的图像，如图6.4.1所示。>xfit<-seq(8,44,len=200)>yfit<-predict(nls.sol,data.frame(x=xfit))>plot(data$x,data$y)>lines(xfit,yfit)图6.4.1函数的拟合曲线§6.5方差分析模型在社会实践中，影响一件事的因素是很多的，人们总是希望通过随机试验来观察各种因素对实验结果的影响。方差分析是基于一定的统计数据，定量地分析一个或多个因素对某个（些）响应变量影响和作用的显著性，这种显著性是基于一定概率条件下而言的，其前提是各因素的作用下，响应变量的分布具有正态性和等方差性。因此，本节先给出样本的正态性检验方法，然后分别介绍单因素和双因素的方差分析。6.5.1样本分布的设是来自总体X的样本，我们先通过直方图，核密度估计曲线和经验分布来描述样本数据的分布，然后对其进行正态性检验。(1)直方图对于数据分布，经常用直方图进行描述，首先将数据的取值范围分成若干区间。在等间隔的情况下，每个区间的长度称为组距。考察数据落入每一区间的频数与频率，在每个区间上画一个矩形，它的宽度是组距，高度是频数或者频率。我们应该注意的是，组距对直方图的形态有很大的影响，组距太小，每组的频数较少，由于随机性的影响，临近区间上的频数可能相差很大；组距太大，直方图所反映的数据形态就不够灵敏。步骤如下：Step1找出数据的最大和最小值，即和；Step2确定数据分布的区间；Step3将区间m等分，即:;Step4统计数据落入区间中的频数，；Step5画图。（2）核密度估计与直方图相配套的是核密度估计（kernaldensityestimate）函数，其目的是用已知样本去估计其密度。核密度估计是一种从数据样本本身出发研究数据分布特征的方法，在R软件中，利用非参数方法进行核密度估计的函数为density()。下面举例说明其用法。例6.5.1某班有50名学生，随机抽取20人进行语文测验，考试成绩如下：7566848092747088909565837277696481777379用hist函数，画20个样本的直方图和核密度估计图。编程如下：x<-c(75,66,84,80,92,74,70,88,90,95,65,83,72,77,69,64,81,77,73,79)hist(x,breaks=c(5*12:19),freq=FALSE,col="yellow",border="red")#频率直方图lines(density(x),col="black",lty=3,lwd=2)#核密度估计a<-65:95lines(a,dnorm(a,mean(x),sd(x)),col="blue",lty=2,lwd=2)#正态分布概率密度函数得到如图6.5.1所示。图6.5.1考试成绩直方图、密度估计曲线(虚线)与正态分布概率密度曲线（实线）通过图6.5.1可以看出，成绩的直方图和核密度估计曲线均反映了20名学生考试成绩的分布特点，注意到密度估计曲线与正态分布的概率密度曲线还有一定的区别。（3）经验分布直方图适合于总体为连续性分布的场合。对于更加一般的总体分布，若要估计它的总体分布函数，可用经验分布函数作估计。设是来自总体X的样本，则称（6.5.1）为经验分布函数(expiricaldistribution)，其中表示中不大于的个数。经验分布函数也可以表示成，（6.5.2）是一个跳跃函数，其跳跃点是样本观测值，在每个跳跃点处跳跃度均为。在R软件中，用函数ecdf()绘制样本的经验分布函数。绘制例6.5.1中20名同学考试成绩的经验分布和相应的正态分布图。编程如下：plot(ecdf(x),verticals=T,do.p=F)a<-64:95lines(a,pnorm(a,mean(x),sd(x)))#正态分布函数曲线其中，verticals是逻辑变量，当verticals=T时，表示画竖线，否则不画竖线。do.p是逻辑变量，当do.p=FALSE时，表示不画点处的记号；否则画记号。运行程序，得到图6.5.2。可见正态分布曲线与经验分布函数具有一致性。图6.5.2考试成绩的经验分布图和正态分布曲线（4）Q-Q图Q-Q图可以帮助我们鉴别样本的分布是否近似某种类型的分布。对于正态Q-Q图检验来说，若为随机变量X的n个观测样本，将其由小到大排序后的顺序统计量为.根据经验分布函数，若，则应该有（6.5.3）又因，则由（这里将修正为为了避免出现1的情况）确定的分位数（）应近似满足。因此，与具有线性相关关系。作正态Q-Q图的步骤如下：Step1将样本的观测值排序：；Step2计算样本分位数对应的概率值，，；Step3计算标准正态分布对应的分位数满足，;Step4将数对，画在直角坐标系中，若呈直线状，则认为是正态的；否则认为是非正态数据。下面我们通过R软件中的qqnorm()和qqline()函数画正态Q-Q图和相应的直线。根据正态Q-Q图，对例6.5.1中的数据进行正态性检验。编程如下：x<-c(75,66,84,80,92,74,70,88,90,95,65,83,72,77,69,64,81,77,73,79)qqnorm(x)qqline(x)得到正态Q-Q图，如图6.5.3所示。从正态Q-Q图来看，样本数据基本上可以看成来自正态总体。可以使用相关系数法，记与的相关系数为，对任意的显著性水平，若，则拒绝正态性假设，即认为数据不是来自正态总体；反之，当则接受原假设，认为数据来自正态总体。图6.5.3学生成绩样本数据正态Q-Q图（4）Shapiro-WilkW统计量检验利用Shapiro-WilkW统计量作正态性检验，因此称这种检验方法为正态W检验方法。该方法是Shapiro和Wilk于1965年提出的一种灵敏度高，计算简单，需要的样本容量较小的正态性检验方法，这一方法是由样本的顺序统计量所构成的统计量W。（6.5.4）其中（6.5.5）这里当n为偶数时，取l=n/2；当n为奇数时，取l=(n+1)/2。可以证明对于任何分布，W的值都介于0和1之间，越接近正态分布，W的值就越接近于1，的取值可通过查表取得，进一步可查得样本容量为n时W的下侧临界值，当时，接受总体为正态分布的假设，否则拒绝正态性假设。在R软件中，函数shapiro.test()提供W统计量和相应的p值，当p值小于某个显著性水平时（例如0.05），则认为样本不是来自正态分布的总体；否则认为样本是来自正态分布的总体。对于例6.5.1中的数据，利用Shapiro-WilkW统计量进行正态性检验，编程如下：x<-c(75,66,84,80,92,74,70,88,90,95,65,83,72,77,69,64,81,77,73,79)Shapiro.test(x)运行程序有如下结果：Shapiro-Wilknormalitytestdata:xW=0.9696,p-value=0.7463得出W统计量为0.9696，相应的p值为0.7463大于0.05的显著性水平，因此可以认为来自正态分布的总体，与QQ图得到结论相同。或者，通过查W统计量的分位数表可得当样本量为n=20时，，即，因此接受总体为正态分布的假设，也可以得到相同的结论。6.5.2若只定量地分析一个因素对响应变量的影响，并且用方差进行分析，这种方法被称为单因素方差分析。方差分析的目的是通过试验数据分析推断因素A对响应变量的影响是否显著，即当因素A取不同水平的响应变量有无显著差异。若因素A有s个水平分别为，在其它条件都不变的情况下，各进行次重复独立试验，得到观测数据列表，如表6.5.1所示。不同因素水平下的观测值之间相互独立，将水平下的试验结果看成来自总体的样本观测值，由于，因此数据可以分解为：（6.5.6）表6.5.1单因素方差分析数据结构表因素水平样本（观测值）总体………其中，为第i个总体的均值，为相应的观测或试验误差。比较因素A的s个水平的差异，归结为比较这s个总体的均值。即检验假设：；不全相等.（6.5.7）我们记则为总体的均值，为水平对指标的效应。则式（6.5.6）也可以写成：（6.5.8）该模型为单因素方差分析的数学模型。假设（6.5.7）等价于；不全相等.（6.5.9）为了检验假设（6.5.9），导出检验统计量，引入总离差平方和统计量，即.它是所有样本观测数据与总平均值差的平方和，描述了所有观测数据的离散程度。经计算可以证明，其中这里的为组内平方和，是由随机误差产生的组内差异的度量；称为组间平方和，表示各影响因素水平下样本均值与总平均值之间的差异之和，反映了s个总体均值之间的差异。可以证明，当假设成立时，统计量（6.5.10）给定显著性水平，用表示F分布上分位点。若，则拒绝原假设，认为因素A的s个水平有显著性差异。也可以通过计算p值的方法来决定是接受还是拒绝原假设。p值为表示服从自由度为(s-1,n-s)的F分布的随机变量大于F的概率，因此，p值小于等价于，表示小概率事件发生了，意味着要拒绝原假设。通常将单因素方差分析的结果列成表6.5.2的形式，称为方差分析表。表6.5.2单因素方差分析表方差来源自由度平方和均方F比p值因素As-1F值p误差n-s总和n-1下面通过例子说明方差分析的过程和求解。例6.5.2有三台机器，用来生产厚度为0.25厘米的铝合金板，现在要了解各机器产品的平均厚度是否相同，取样测量精确至千分之一厘米，得结果如下：机器1：0.236,0.238,0.248,0.245,0.243机器2：0.257,0.253,0.255,0.254,0.261,0,258机器3：0.264,0.259,0.267,0.262下面我们用R软件提供的aov()函数来进行方差分析。解：用数据框的格式输入数据，调用aov()函数计算方差分析，用summary()提取方法分析的信息>df<-data.frame(+X=c(0.236,0.238,0.248,0.245,0.243,0.257,0.253,0.255,+0.254,0.261,0.258,0.264,0.259,0.267,0.262),+A=factor(c(rep(1,5),rep(2,6),rep(3,4)))+)>dfXA10.236120.238130.248140.245150.243160.257270.253280.255290.2542100.2612110.2582120.2643130.2593140.2673150.2623>df.aov<-aov(X~A,data=df)>summary(df.aov)DfSumSqMeanSqFvaluePr(>F)A20.00107005.35e-0436.627.79e-06***Residuals120.00017531.46e-05---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1从最后的计算结果可以看出，F统计量为36.62，p值为，查表可知，，且p值也小于0.05的显著性水平，说明拒绝原假设，即各台及其生产的薄板厚度有显著性差异。6.5.3在大量的实际问题中，经常需要考虑影响试验数据的因素多于一个的情形。本节我们考虑影响因素为两个的情况。例如，影响火箭射程有两个重要因素，燃料因素A和推进器因素B，这两个因素都可能影响结果，研究人员很想知道不同水平组合对火箭射程是否有显著影响？在对火箭进行模拟的试验中，取出了四种燃料和三种推进器，每种燃料和每种推进器的组合各发射火箭三次，得到数据如表6.5.3所示。表6.5.3火箭射程试验数据（单位：km）63.264.762.156.255.757.060.359.761.659.159.760.054.153.854.541.642.740.668.167.968.070.968.971.039.238.940.575.875.476.258.259.157.640.740.641.0事实上，不考虑交互作用的方差分析有一个前提假设，即假定A，B两个因素对指标的效应是可以叠加的，而且认为因素A各种效应的比较与因素B处在什么水平无关。而这种影响在大部分实际问题中也应该给予足够重视，这里，我们考虑因素的各种水平组合的不同带来的影响，即建立考虑交互作用的双因素方差分析模型。对于两个因素A和B，设因素A有r个水平分别为，因素B有s个水平分别为，在其它条件都不变的情况下，每种水平组合下重复试验t次，得到观测数据列表，如表6.5.4所示。表6.5.4双………将水平组合下的试验结果看成来自总体的样本观测值，由于且相互独立，则数据可以分解为：（6.5.11）其中，为总平均，为因素A第i个水平的效应，为因素B第j个水平的效应，表示和的交互效应，因此有：判断因素A，B的及交互效应的影响是否显著等价于检验下列假设：（6.5.12）为了检验假设（6.5.12），导出检验统计量，其中，,,其中为总离差平方和，为误差平方和，为因素的平方和,为因素的平方和，为交互效应平方和。可以证明：①当成立时，；②当成立时，；③当成立时，.分别以、和作为、和检验统计量，将计算结果列成方差分析表，如表6.5.5所示。与单因素方差分析类似，给定显著性水平，若检验统计量的F值大于相应显著性水平下对应的分位点，则拒绝原假设，反之则接受原假设；也可以通过p值来检验，即p值小于显著性水平，则表示小概率事件发生了，应拒绝原假设，反之则接受原假设。表6.5.5双方差来源自由度平方和均方F比p值因素Ar-1因素Bs-1交互效应(r-1)(s-1)误差rs(t-1)总和n-1对于表6.5.3中所给出的数据，对实验结果进行双因素方差分析，确定火箭的燃料和推进器因素是否对火箭射程有显著性影响。通过R软件编程，用数据框的形式输入数据，调用aov()函数计算，用summary()，提取方差分析的信息。程序和计算结果见如图6.5.4所示。从计算结果可以看出，因素A和B以及交互因素都对火箭射程有显著影响。图6.5.4火箭射程双因素方差分析§6.6主成分分析和因子分子模型在对实际问题的建模过程中，为了全面的分析问题，往往涉及众多有关的变量。但是变量太多不但会增加计算的复杂性，而且也给合理地分析和解释问题带来困难。一般来说，虽然每个变量都提供了一定的信息，但其重要性有所不同。实际上，在很多情况下，众多变量之间有一定的相关关系，人们希望利用这种相关性对这些变量加以“改造”，用维数较少的新变量来反映原变量所提供的大部分信息，通过对新变量的分析达到解决问题的目的，主成分分析和因子分析便