《数学建模》算法全收录（算法大全）：第11章 方差分析

第十一章方差分析 我们已经作过两个总体均值的假设检验,如两台机床生产的零件尺寸是否相等,病 人和正常人的某个生理指标是否一样。如果把这类问题推广一下,要检验两个以上总体 的均值彼此是否相等,仍然用以前介绍的方法是很难做到的。而你在实际生产和生活中 可以举出许多这样的问题:从用几种不同工艺制成的灯泡中,各抽取了若干个测量其寿 命,要推断这几种工艺制成的灯泡寿命是否有显著差异:用几种化肥和几个小麦品种在 若干块试验田里种植小麦,要推断不同的化肥和品种对产量有无显著影响 可以看到,为了使生产过程稳定,达到优质、高产,需要对影响产品质量的因素进 行分析,找出有显著影响的那些因素,除了从机理方面进行研究外,常常要作许多试验, 对结果作分析、比较,寻求规律。用数理统计分析试验结果、鉴别各因素对结果影响程 度的方法称为方差分析( Analysis Of Variance),记作 ANOVA。 人们关心的试验结果称为指标,试验中需要考察、可以控制的条件称为因素或因子 因素所处的状态称为水平。上面提到的灯泡寿命问题是单因素试验,小麦产量问题是双 因素试验。处理这些试验结果的统计方法就称为单因素方差分析和双因素方差分析 §1单因素方差分析 只考虑一个因素A对所关心的指标的影响,A取几个水平,在每个水平上作若干 个试验,试验过程中除A外其它影响指标的因素都保持不变(只有随机因素存在),我 们的任务是从试验结果推断,因素A对指标有无显著影响,即当A取不同水平时指标 有无显著差别。 A取某个水平下的指标视为随机变量,判断A取不同水平时指标有无显著差别, 相当于检验若干总体的均值是否相等。 1.1数学模型 设A取r个水平A1,A12…,A,在水平A下总体x服从正态分布N(u1,a2), i=1,…,r,这里12a未知,共可以互不相同,但假定x2有相同的方差。又设在每 个水平A下作了n次独立试验,即从中抽取容量为n的样本,记作x,j=1…,n,x 服从N(H12-),i=1…,r,j=1,…,n且相互独立。将这些数据列成表1(单因素试 验数据表)的形式。 表1单因素试验数据表 A, A 将第i行称为第i组数据。判断A的r个水平对指标有无显著影响,相当于要作以 下的假设检验 H0:1=12=…=山;H1:H1,42,…,H1不全相等 由于x的取值既受不同水平A的影响,又受A固定下随机因素的影响,所以将它 分解为 x=H1+5,i=l…,r,j=1,…,n (1) 213

-213- 第十一章 方差分析 我们已经作过两个总体均值的假设检验，如两台机床生产的零件尺寸是否相等，病 人和正常人的某个生理指标是否一样。如果把这类问题推广一下，要检验两个以上总体 的均值彼此是否相等，仍然用以前介绍的方法是很难做到的。而你在实际生产和生活中 可以举出许多这样的问题：从用几种不同工艺制成的灯泡中，各抽取了若干个测量其寿 命，要推断这几种工艺制成的灯泡寿命是否有显著差异；用几种化肥和几个小麦品种在 若干块试验田里种植小麦，要推断不同的化肥和品种对产量有无显著影响。 可以看到，为了使生产过程稳定，达到优质、高产，需要对影响产品质量的因素进 行分析，找出有显著影响的那些因素，除了从机理方面进行研究外，常常要作许多试验， 对结果作分析、比较，寻求规律。用数理统计分析试验结果、鉴别各因素对结果影响程 度的方法称为方差分析（Analysis Of Variance），记作 ANOVA。 人们关心的试验结果称为指标，试验中需要考察、可以控制的条件称为因素或因子， 因素所处的状态称为水平。上面提到的灯泡寿命问题是单因素试验，小麦产量问题是双 因素试验。处理这些试验结果的统计方法就称为单因素方差分析和双因素方差分析。 §1 单因素方差分析 只考虑一个因素 A 对所关心的指标的影响， A 取几个水平，在每个水平上作若干 个试验，试验过程中除 A 外其它影响指标的因素都保持不变（只有随机因素存在），我 们的任务是从试验结果推断，因素 A 对指标有无显著影响，即当 A 取不同水平时指标 有无显著差别。 A 取某个水平下的指标视为随机变量，判断 A 取不同水平时指标有无显著差别， 相当于检验若干总体的均值是否相等。 1.1 数学模型 设 A 取 r 个水平 A A Ar , , , 1 2 L ，在水平 Ai 下总体 i x 服从正态分布 ( , ) 2 N μi σ ， i = 1,L,r ，这里 2 μi ,σ 未知， μi 可以互不相同，但假定 i x 有相同的方差。又设在每 个水平 Ai 下作了 ni 次独立试验，即从中抽取容量为 ni 的样本，记作 ij ni x , j =1,L, ， ij x 服从 ( , ) 2 N μi σ ， ni i =1,L,r, j =1,L, 且相互独立。将这些数据列成表 1（单因素试 验数据表）的形式。 表 1 单因素试验数据表 A1 11 x 12 x … 1n1 x A2 21 x 22 x … 2n2 x … … … … … Ar r1 x r 2 x … r rn x 将第i 行称为第i 组数据。判断 A 的 r 个水平对指标有无显著影响，相当于要作以 下的假设检验 H0 μ1 = μ2 =L= μ r : ； H μ μ μ r : , , , 1 1 2 L 不全相等 由于 ij x 的取值既受不同水平 Ai 的影响，又受 Ai 固定下随机因素的影响，所以将它 分解为 ij i ij x = μ + ε ，i = 1,L,r ， ni j =1,L, （1）

其中E~N(0.0),且相互独立。记 n1,n=∑n,a1=H-H,i=1,…,r (2) μ是总均值,∝1是水平A对指标的效应。由(1)、(2)模型可表为 x=u+a+++++ta (3) N(0,a2),i=1,…,r,j=1, 原假设为(以后略去备选假设) H0:a1=a2=…=ar=0 (4) 12统计分析 x n 是第i组数据的组平均值,是总平均值。考察全体数据对x的偏差平方和 S=∑∑(x1-x) (6) 经分解可得 n ( x (7) ∑∑(x1-x 则 ST=s+sE (9) SA是各组均值对总方差的偏差平方和,称为组间平方和;Sg是各组内的数据对均值偏 差平方和的总和。S』反映A不同水平间的差异,SE则表示在同一水平下随机误差的 大小。 注意到∑(x1-x)是总体N(A,o2)的样本方差的n-1倍,于是有 ∑(x-x) 由x2分布的可加性知

-214- 其中 ~ (0, ) 2 ε ij N σ ，且相互独立。记 ∑= = r i ni i n 1 1 μ μ ， ∑= = r i n ni 1 ，αi = μi − μ ， i = 1,L,r （2） μ 是总均值，αi 是水平 Ai 对指标的效应。由（1）、（2）模型可表为 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = = = = + + ∑= ij i r i i ij i ij N i r j n x ~ (0, ), 1, , , 1, , 0 2 1 ε σ L L α μ α ε （3） 原假设为（以后略去备选假设） H0 :α1 = α 2 =L= α r = 0 （4） 1.2 统计分析 记 ∑= • = i n j ij i i x n x 1 1 ， ∑∑= = = r i n j ij i x n x 1 1 1 （5） i• x 是第i 组数据的组平均值， x 是总平均值。考察全体数据对 x 的偏差平方和 ∑∑= = = − r i n j T ij i S x x 1 1 2 ( ) （6） 经分解可得 ∑ ∑∑= = • = = • − + − r i n j ij i r i T i i i S n x x x x 1 1 2 1 2 ( ) ( ) 记 ∑= = • − r i A i i S n x x 1 2 ( ) （7） ∑∑= = = − • r i n j E ij i i S x x 1 1 2 ( ) （8） 则 T A E S = S + S （9） A S 是各组均值对总方差的偏差平方和，称为组间平方和； E S 是各组内的数据对均值偏 差平方和的总和。 A S 反映 A 不同水平间的差异， E S 则表示在同一水平下随机误差的 大小。 注意到∑= − • ni j ij i x x 1 2 ( ) 是总体 ( , ) 2 N μi σ 的样本方差的ni −1倍，于是有 ( ) ~ ( 1) 2 1 2 2 ∑ − − = • i n j xij xi n i σ χ 由 2 χ 分布的可加性知

S x2∑(n1-1) S/o2-x(n-r 且有 (10) 对S:作进一步分析可得 ESA=(r-1)a2+∑nax2 (11) 当H。成立时 Es,=(r-D)o (12) 可知若H0成立,S4只反映随机波动,而若H不成立,那它就还反映了A的不同水平 的效应a1。单从数值上看,当H0成立时,由(10)、(12)对于一次试验应有 -1)=1 SE/(n-r) 而当H0不成立时这个比值将远大于1。当H0成立时,该比值服从自由度n1=r-1, n2=(n-r)的F分布,即 F=s w F(r-l, n-r) (13) Se/n-r 为检验H,给定显著性水平a,记F分布的1-a分位数为Fa(r-1、(n-r),检验 规则为 F<Fa(-1,(n-r)时接受Ho,否则拒绝 以上对SA,Sg,S的分析相当于对组间、组内等方差的分析,所以这种假设检验方法称 方差分析 13方差分析表 将试验数据按上述分析、计算的结果排成表2的形式,称为单因素方差分析表 ( Matlab中给出的方差分析表)。 表2单因素方差分析表 方差来源 平方和 自由度 均方 P分位数 概率 因素A SA 1|S SA Fi-p, (r-I, n-r)I p 误差 S 总和 最后一列给出大于F值的概率P,F-n<F(r-1(n-)相当于P,>a 215-

-215- ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∑ − = r i SE ni 1 2 2 σ ~ χ ( 1) 即 ~ ( ) 2 2 S n r E σ χ − 且有 2 ESE = (n − r)σ （10） 对 A S 作进一步分析可得 ∑= = − + r i A ni i ES r 1 2 2 ( 1)σ α （11） 当 H0 成立时 2 ES A = (r −1)σ （12） 可知若 H0 成立， A S 只反映随机波动，而若 H0 不成立，那它就还反映了 A 的不同水平 的效应αi 。单从数值上看，当 H0 成立时，由（10）、（12）对于一次试验应有 1 /( ) /( 1) ≈ − − S n r S r E A 而当 H0 不成立时这个比值将远大于 1。当 H0 成立时，该比值服从自由度 1 n1 = r − ， ( ) 2 n = n − r 的 F 分布，即 ~ ( 1, ) /( ) /( 1) F r n r S n r S r F E A − − − − = （13） 为检验 H0 ，给定显著性水平α ，记 F 分布的1−α 分位数为 ( 1,( )) 1 F r − n − r −α ，检验 规则为 ( 1,( )) 1 F < F r − n − r −α 时接受 H0 ，否则拒绝。 以上对 A E T S , S , S 的分析相当于对组间、组内等方差的分析，所以这种假设检验方法称 方差分析。 1.3 方差分析表 将试验数据按上述分析、计算的结果排成表 2 的形式，称为单因素方差分析表 （Matlab 中给出的方差分析表）。 表 2 单因素方差分析表 方差来源 平方和 自由度 均方 r 1− p 分位数 概率 因素 A A S r −1 −1 = r S S A A ( 1, ) 1 F r n r pr − − − r p 误差 E S n − r n r S S E E − = 总和 T S n −1 最后一列给出大于 F 值的概率 r p ， ( 1,( )) 1 1 F F r n r pr − < −α − − 相当于 pr > α

方差分析一般用的显著性水平是:取a=0.01,拒绝H0,称因素A的影响(或A 各水平的差异)非常显著;取a=0.01,不拒绝H,但取a=0.05,拒绝H0,称因 素A的影响显著;取a=0.05,不拒绝H0,称因素A无显著影响 4 Matlab实现 Matlab统计工具箱中单因素方差分析的命令是 anova 若各组数据个数相等,称为均衡数据。若各组数据个数不等,称非均衡数据。 (1)均衡数据 处理均衡数据的用法为 p=anova(x) 返回值p是一个概率,当p>a时接受H0,x为mxr的数据矩阵,x的每一列是一个 水平的数据(这里各个水平上的样本容量n1=m)。另外,还输出一个方差表和一个 Box图 例1为考察5名工人的劳动生产率是否相同,记录了每人4天的产量,并算出其 平均值,如表3。你能从这些数据推断出他们的生产率有无显著差别吗? 表3 A A2 A, A As 1234 256 50248 330 280252 0 平均产量269 9225264.75280.5240 解编写程序如下 x=[256 254 236 242 330 280252 280 290 230 305220 295 302289252]; p=-anoval(x) 求得P=0.1109>a=0.05,故接受H0,即5名工人的生产率没有显著差异。方 差表对应于上面的单因素方差分析表的1~4列,F=2262是F(415)分布的1-p分 位数,可以验证 fcdf(2.262,4,15)=0.8891=1-p Box图反映了各组数据的特征。 注:接受H0,是将5名工人的生产率作为一个整体进行假设检验的结果,并不表 明取其中2个工人的生产率作两总体的均值检验时,也一定接受均值相等的假设。实际 上,读者可以用 ttest2对本题作H0:42=p5的检验,看看会得到什么结果。 (2)非均衡数据 处理非均衡数据的用法为: p=anova (x, group ⅹ为向量,从第1组到第r组数据依次排列; group为与x同长度的向量,标志x中数 据的组别(在与x第i组数据相对应的位置处输入整数(=1,2,…,r))

-216- 方差分析一般用的显著性水平是：取α = 0.01，拒绝 H0 ，称因素 A 的影响（或 A 各水平的差异）非常显著；取α = 0.01，不拒绝 H0 ，但取α = 0.05 ，拒绝 H0 ，称因 素 A 的影响显著；取α = 0.05 ，不拒绝 H0 ，称因素 A 无显著影响。 1.4 Matlab 实现 Matlab 统计工具箱中单因素方差分析的命令是 anoval。 若各组数据个数相等，称为均衡数据。若各组数据个数不等，称非均衡数据。 （1）均衡数据 处理均衡数据的用法为： p=anoval(x) 返回值 p 是一个概率，当 p > α 时接受 H0 ，x 为m× r 的数据矩阵，x 的每一列是一个 水平的数据（这里各个水平上的样本容量 ni = m ）。另外，还输出一个方差表和一个 Box 图。 例 1 为考察 5 名工人的劳动生产率是否相同，记录了每人 4 天的产量，并算出其 平均值，如表 3。你能从这些数据推断出他们的生产率有无显著差别吗？ 表 3 工人 天 A1 A2 A3 A4 A5 1 256 254 250 248 236 2 242 330 277 280 252 3 280 290 230 305 220 4 298 295 302 289 252 平均产量 269 292.25 264.75 280.5 240 解 编写程序如下： x=[256 254 250 248 236 242 330 277 280 252 280 290 230 305 220 298 295 302 289 252]; p=anova1(x) 求得 p = 0.1109 >α = 0.05，故接受 H0 ，即 5 名工人的生产率没有显著差异。方 差表对应于上面的单因素方差分析表的1 ~ 4 列，F = 2.262 是 F(4,15) 分布的1− p 分 位数，可以验证 fcdf(2.262,4,15)=0.8891=1-p Box 图反映了各组数据的特征。 注：接受 H0 ，是将 5 名工人的生产率作为一个整体进行假设检验的结果，并不表 明取其中 2 个工人的生产率作两总体的均值检验时，也一定接受均值相等的假设。实际 上，读者可以用 ttest2 对本题作 0 2 5 H : μ = μ 的检验，看看会得到什么结果。 （2）非均衡数据 处理非均衡数据的用法为： p=anova1(x,group) x 为向量，从第 1 组到第 r 组数据依次排列；group 为与 x 同长度的向量，标志 x 中数 据的组别（在与 x 第i 组数据相对应的位置处输入整数i(i = 1,2,L,r)）

例2用4种工艺生产灯泡,从各种工艺制成的灯泡中各抽出了若干个测量其寿命, 结果如下表,试推断这几种工艺制成的灯泡寿命是否有显著差异。 表4 工艺 A, A3 44 1620158014601500 1670160015401550 4 1680 1800 解编写程序如下 [1620158014601500 1670160015401550 1700164016201610 1750172016801800] x=[x(1:4),x(16),x(5:8),x(9:11),x(12:15)]; g=[ones(1,5),2*ones(1,4),3*ones(1,3),4+ones(1,4) p=anova (x, g) 求得0.01<p=0.0331<0.05,所以几种工艺制成的灯泡寿命有显著差异。 1.5多重比较 在灯泡寿命问题中,为了确定哪几种工艺制成的灯泡寿命有显著差异,我们先算出 各组数据的均值 工艺 均值17081635 15401585 虽然A的均值最大,但要判断它与其它几种有显著差异,还需做多重比较。一般 多重比较要对所有r个总体作两两对比,分析相互间的差异。根据问题的具体情况可以 减少对比次数。 对于上述问题, Matlab多重比较的程序为 x=[1620158014601500 1670160015401550 1700164016201610 1750172016801800] x=[x(1:4),x(16),x(5:8),x(9:11),x(12:15) g=[ones(1,5),2*ones(1,4),3*ones(1,3),4*ones(1,4)]; Ip, t, st]=anoval(x, g) [c,m,h, nms]= multcompare(st) 11(m)] §2双因素方差分析 如果要考虑两个因素A,B对指标的影响,A,B各划分几个水平,对每一个水平组 合作若干次试验,对所得数据进行方差分析,检验两因素是否分别对指标有显著影响, 或者还要进一步检验两因素是否对指标有显著的交互影响。 2.1数学模型 设A取r个水平A,A2,…,A,B取$个水平B1,B2,…,B4,在水平组合(4,B)

-217- 例 2 用 4 种工艺生产灯泡，从各种工艺制成的灯泡中各抽出了若干个测量其寿命， 结果如下表，试推断这几种工艺制成的灯泡寿命是否有显著差异。 表 4 工艺 A1 A2 A3 A4 序号 1 1620 1580 1460 1500 2 1670 1600 1540 1550 3 1700 1640 1620 1610 4 1750 1720 1680 5 1800 解 编写程序如下： x=[1620 1580 1460 1500 1670 1600 1540 1550 1700 1640 1620 1610 1750 1720 1680 1800]; x=[x(1:4),x(16),x(5:8),x(9:11),x(12:15)]; g=[ones(1,5),2*ones(1,4),3*ones(1,3),4*ones(1,4)]; p=anova1(x,g) 求得 0.01<p=0.0331<0.05，所以几种工艺制成的灯泡寿命有显著差异。 1.5 多重比较 在灯泡寿命问题中，为了确定哪几种工艺制成的灯泡寿命有显著差异，我们先算出 各组数据的均值： 工艺 A1 A2 A3 A4 均值 1708 1635 1540 1585 虽然 A1 的均值最大，但要判断它与其它几种有显著差异，还需做多重比较。一般 多重比较要对所有 r 个总体作两两对比，分析相互间的差异。根据问题的具体情况可以 减少对比次数。 对于上述问题，Matlab 多重比较的程序为 x=[1620 1580 1460 1500 1670 1600 1540 1550 1700 1640 1620 1610 1750 1720 1680 1800]; x=[x(1:4),x(16),x(5:8),x(9:11),x(12:15)]; g=[ones(1,5),2*ones(1,4),3*ones(1,3),4*ones(1,4)]; [p,t,st]=anova1(x,g) [c,m,h,nms] = multcompare(st); [nms num2cell(m)] §2 双因素方差分析 如果要考虑两个因素 A,B 对指标的影响， A,B 各划分几个水平，对每一个水平组 合作若干次试验，对所得数据进行方差分析，检验两因素是否分别对指标有显著影响， 或者还要进一步检验两因素是否对指标有显著的交互影响。 2.1 数学模型 设 A 取 r 个水平 A A Ar , , , 1 2 L ，B 取 s 个水平 B B Bs , , , 1 2 L ，在水平组合( , ) Ai Bj