08-第八章 方差分析与回归分析
依赖于
被以下题目直接调用
正文部分
§8.1 方差分析
- 单因子方差分析
(1) 问题与数据
设因子有 r 个水平,记为 A1,A2,⋯,Ar,在每一水平下各做 m 次独立重复试验,若记第 i 个水平下第 j 次重复的试验结果为 yij,所有试验的结果可列表如下:
因子水平A1A2⋮Ar试验数据y11, y12, ⋯, y1my21, y22, ⋯, y2m⋮yr1, yr2, ⋯, yrm和T1T2⋮TrT平均yˉ1⋅yˉ2⋅⋮yˉr⋅yˉ
对这个试验要研究的问题是:r 个水平 A1,A2,⋯,Ar 间有无显著差异。
(2) 基本假定
A1: 正态性:第 i 个水平下的数据 yi1,yi2,⋯,yim 是来自正态总体 N(μi,σi2) 的一个样本,i=1,2,⋯,r;
A2: 等方差性:r 个方差相同,即
σ12=σ22=⋯=σr2=σ2;
A3: 独立性:诸数据 yij 都相互独立。
在这三个基本假定下,要检验的假设是
H0:μ1=μ2=⋯=μrvsH1:μ1,μ2,⋯,μr 不全相等。
方差分析就是在上述三个基本假定下,对若干个正态均值是否相等作检验。
(3) 平方和分解式
ST=SA+Se,fT=fA+fe,
若记
yˉi⋅=m1j=1∑myij,yˉ=rm1i=1∑rj=1∑myij=r1i=1∑ryˉi⋅,
上述诸平方和分别为
1.
ST=i=1∑rj=1∑m(yij−yˉ)2
称为总偏差平方和,其自由度
fT=rm−1=n−1;
-
SA=mi=1∑r(yˉi⋅−yˉ)2
称为组间偏差平方和或因子 A 的偏差平方和,其自由度
fA=r−1;
-
Se=i=1∑rj=1∑m(yij−yˉi⋅)2
称为组内偏差平方和或误差偏差平方和,其自由度
fe=r(m−1)=n−r.
注:数据 yij 的平移
y′=yij−a
不会改变其平方和的值。用此性质可简化计算。
(4) 方差分析表
来源因子误差总和平方和SA=m1∑i=1rTi2−rmT2Se=ST−SAST=∑i=1r∑j=1myij2−rmT2自由度fA=r−1fe=r(m−1)fT=rm−1均方MSA=SA/fAMSe=Se/feF 比F=MSA/MSe
(5) 判断
在 H0 成立下,
F=MSA/MSe∼F(fA,fe),
对给定的显著性水平 α (0<α<1),其拒绝域为
W={F≥F1−α(fA,fe)},
其中 F1−α(fA,fe) 可从附表 5 中查得。
- 若
F≥F1−α(fA,fe),
则认为因子 A 显著,即诸正态均值间有显著差异;
- 若
F<F1−α(fA,fe),
则说明因子 A 不显著,即接受原假设 H0。
给出检验的 p 值是更常用的,若以 X 记服从 F(fA,fe) 的随机变量,F 为统计量,F=MSA/MSe 的观测值为 F0,则
p=P(X≥F0),
可用软件计算,如在 MATLAB 中使用如下命令:
p=1−fcdf(F0,fA,fe).
- 数据结构式及其参数估计
(1) 数据结构式
yij=μ+ai+εij,i=1,2,⋯,r, j=1,2,⋯,m,
其中 μ 为总均值,ai 为第 i 个水平的效应,且
i=1∑rai=0,
εij 为试验误差,所有 εij 可作为来自 N(0,σ2) 的一个样本,在上述数据结构式下,
yij∼N(μ+ai,σ2).
要检验的假设可改写为
H0:a1=a2=⋯=ar=0vsH1:a1,a2,⋯,ar 不全为 0.
(2) 点估计
- 总均值 μ 的估计
μ^=yˉ;
- 各水平均值 μi 的估计
μ^i=yˉi⋅,i=1,2,⋯,r;
- 主效应 ai 的估计
a^i=yˉi⋅−yˉ,i=1,2,⋯,r;
- 误差方差 σ2 的估计
σ^2=MSe=Se/fe.
(3) 1−α 置信区间(0<α<1)
- μi 的 1−α 置信区间为
[yˉi⋅−σ^⋅t1−α/2(fe)/m, yˉi⋅+σ^⋅t1−α/2(fe)/m].
- 单因子试验的统计分析可给出如下三个结果:
- 因子 A 是否显著;
- 试验误差方差 σ2 的估计;
- 诸水平均值 μi 的点估计与区间估计(此项在因子 A 不显著时无需进行)。
- 重复数不等情形下的方差分析
(1) 数据略有不同
设因子 A 有 r 个水平 A1,A2,⋯,Ar,并且第 i 个水平 Ai 下重复进行 mi 次试验,获得如下数据:
水平A1A2⋮Ar合计重复数m1m2⋮mrn数据y11,y12,⋯,y1m1y21,y22,⋯,y2m2⋮yr1,yr2,⋯,yrmrT和T1T2⋮TrT平均yˉ1⋅yˉ2⋅⋮yˉr⋅yˉ
(2) 基本假定、平方和分解、方差分析及判断准则都和前面一样,只是因子 A 的平方和 SA 的计算公式略有不同:记
n=i=1∑rmi,
则
SA=m1T12+m2T22+⋯+mrTr2−nT2.
(3) 数据结构式及其参数估计基本同前,但要注意以下两点:
- 总均值
μ=n1i=1∑rmiμi;
- 主效应的约束条件为
i=1∑rmiai=0.
习题与解答 8.1
在一个单因子试验中,因子 A 有三个水平,每个水平下各重复 4 次,具体数据如下:
水平一水平二水平三水平\multicolumn4c数据86051017125492
试计算误差平方和 Se、因子 A 的平方和 SA 与总平方和 ST,并指出它们各自的自由度。
解
此处因子水平数 r=3,每个水平下的重复次数 m=4,总试验次数为
n=mr=12.
首先,算出每个水平下的数据和以及总数据和:
T1=8+5+7+4=24,
T2=6+10+12+9=37,
T3=0+1+5+2=8,
T=T1+T2+T3=24+37+8=69.
误差平方和 Se 由三个平方和组成:
Se1=(82+52+72+42)−4242=10,fe1=3,
Se2=(62+102+122+92)−4372=18.75,fe2=3,
Se3=(02+12+52+22)−482=14,fe3=3.
于是
Se=Se1+Se2+Se3=10+18.75+14=42.75,
fe=fe1+fe2+fe3=9.
而
SA=i=1∑rmTi2−nT2=4242+4372+482−12692=105.5,fA=3−1=2,
ST=i=1∑rj=1∑myij2−12T2=(82+52+⋯+22)−12692=148.25,fT=12−1=11.
注:在所有的计算公式中可以用数据和也可以用数据均值,但在实际计算中应尽可能用和而不用均值,这样可避免不必要的误差(如四舍五入造成的误差等)。
在一个单因子试验中,因子 A 有 4 个水平,每个水平下重复次数分别为 5,7,6,8。那么误差平方和、A 的平方和及总平方和的自由度各是多少?
解
此处因子水平数 r=4,总试验的次数
n=5+7+6+8=26,
因而有
误差平方和的自由度 fe=n−r=26−4=22,
因子 A 的平方和的自由度 fA=r−1=3,
总平方和的自由度 fT=n−1=25.
在单因子试验中,因子 A 有 4 个水平,每个水平下各重复 3 次试验,现已求得每个水平下试验结果的样本标准差分别为 1.5,2.0,1.6,1.2,则其误差平方和为多少?误差的方差 σ2 的估计值是多少?
解
此处因子水平数 r=4,每个水平下的试验次数 m=3,误差平方和 Se 由四个平方和组成,它们分别为
Se1=(3−1)×1.52=4.5,Se2=(3−1)×22=8,
Se3=(3−1)×1.62=5.12,Se4=(3−1)×1.22=2.88.
于是
Se=Se1+Se2+Se3+Se4=20.5,
其自由度为
fe=mr−r=8,
误差方差 σ2 的估计值为
σ^2=feSe=820.5=2.5625.
在单因子方差分析中,因子 A 有 3 个水平,每个水平均各做 4 次重复试验。请完成下列方差分析表,并在显著性水平 α=0.05 下对因子 A 是否显著作出检验:
来源因子 A误差 e总和 T平方和4.22.56.7自由度均方F 比p 值
解
补充的方差分析表如下:
来源因子 A误差 e总和 T平方和4.22.56.7自由度2911均方2.10.28F 比7.5p 值0.0121
对于给定的显著性水平 α=0.05,查表得
F0.95(2,9)=4.26,
故拒绝域为
W={F≥4.26}.
由于
F=7.5>4.26,
因而认为因子 A 是显著的。此外,
p=1−fcdf(7.5,2,9)=0.0121.
用 4 种安眠药在兔子身上进行试验,挑选 24 只健康兔子,随机把它们均分为 4 组,每组各服一种安眠药,安眠时间(单位:h)如下:
安眠药A1A2A3A4\multicolumn6c安眠时间6.26.36.85.46.16.57.16.46.06.76.66.26.36.66.86.36.17.16.96.05.96.46.65.9
在显著性水平 α=0.05 下对其进行方差分析,可以得到什么结果?
解
这是一个单因子方差分析的问题。根据样本数据有
安眠药A1A2A3A4和m666624Ti36.639.640.836.2153.2Ti21339.561568.161664.641310.445882.8∑j=1myij2223.36261.76277.62219.06981.8
于是
ST=i=1∑rj=1∑myij2−nT2=981.8−24153.22=3.87,fT=23,
SA=m1i=1∑rTi2−nT2=65882.8−24153.22=2.54,fA=3,
Se=ST−SA=1.33,fe=20.
方差分析表为
来源因子 A误差 e总计平方和2.541.333.87自由度32023均方0.84670.0665F 比12.7323
在显著性水平 α=0.05 下,查表得
F0.95(3,20)=3.10.
由于 F=12.7323>3.10,故认为因子 A(安眠药)是显著的,即四种安眠药对兔子的安眠作用有明显差别。检验的 p 值为
p=1−fcdf(12.7323,3,20)=7.048×10−5.
为研究单因子 A 中咖啡因剂量对人的影响,随机选择 30 名体质大致相同的健康男大学生进行手指叩击试验。咖啡因有三个水平:
A1=0mg,A2=100mg,A3=200mg.
每个水平下各安排 10 人,试验结果如下:
咖啡因剂量A1A2A3\multicolumn10c叩击次数242248246245246248244245250248247252247248248248250250242247246244246248246243245242244250
对于上述数据进行方差分析,从中可得到什么结论?
解
对数据同时减去 240 不改变方差分析结果。将计算结果列入下表:
水平A1A2A3和\multicolumn10c∣数据(原始数据−240)286Ti568Ti24510∑j=1myij28712788810102764686352410486483195230440966889132892824487331463
于是
ST=1463−301952=195.5,fT=29,
SA=1013289−301952=61.4,fA=2,
Se=ST−SA=134.1,fe=27.
方差分析表为
来源因子 A误差 e总计平方和61.4134.1195.5自由度22729均方30.74.9667F 比6.1812
当 α=0.05 时,F0.95(2,27)=3.362。由于 F=6.1812>3.362,故认为因子 A(咖啡因剂量)是显著的,即三种不同剂量对人的作用有明显差别。检验的 p 值为
p=1−fcdf(6.1812,2,27)=0.0062.
某粮食加工厂试验三种不同的储藏方法对粮食含水率有无显著影响。现取一批粮食分成若干份,分别用三种不同的方法储藏,过一段时间后测得的含水率(单位:%)如下:
储藏方法A1A2A3\multicolumn5c含水率数据7.35.47.98.37.49.57.67.110.08.46.89.88.35.38.4
- 假定各种方法储藏的粮食含水率服从正态分布,且方差相等,试在 α=0.05 下检验这三种方法对含水率的平均值有无显著影响;
- 对每种方法的平均含水率给出置信水平为 0.95 的置信区间。
解
(1) 这是一个单因子方差分析问题。由数据计算如下:
储藏方法A1A2A3和Ti39.932.045.6117.5Ti21592.011024.002079.364695.37∑j=1myij2319.39208.66419.26947.31
ST=947.31−15117.52=26.893,fT=14,
SA=54695.37−15117.52=18.657,fA=2,
Se=ST−SA=8.236,fe=12.
于是方差分析表为
来源因子 A误差 e总计平方和18.6578.23626.893自由度21214均方9.3290.686F 比13.599
在显著性水平 α=0.05 下,查表得 F0.95(2,12)=3.89。由于 F=13.599>3.89,故认为因子 A(储藏方法)显著,即三种不同储藏方法对粮食含水率有显著影响。检验的 p 值为
p=1−fcdf(13.599,2,12)=8.2320×10−4.
(2) 各水平均值估计分别为
μ^1=7.98,μ^2=6.40,μ^3=9.12.
误差方差的无偏估计为
σ^2=MSe=0.686,σ^=0.686=0.828.
若取 α=0.05,则
t0.975(12)=2.178,σ^t0.975(12)/5=0.807,
于是三种方法平均含水率的 0.95 置信区间分别为
μ1:[7.98±0.807]=[7.173,8.787],
μ2:[6.40±0.807]=[5.593,7.207],
μ3:[9.12±0.807]=[8.313,9.927].
在人员推销中有五种方法,某大公司想比较这五种方法有无显著的效果差异。设计一项实验:从应聘的且无推销经验的人员中随机挑选一部分人,将他们随机分为五组,每组用一种推销方法进行培训。培训相同时间后观察他们在一个月内的推销额(单位:千元),数据如下:
组别第一组第二组第三组第四组第五组\multicolumn7c推销额20.024.916.017.525.216.821.320.118.226.217.922.617.320.226.921.230.220.917.729.323.929.922.019.130.426.822.526.818.429.722.420.720.816.528.2
- 假定数据满足进行方差分析的条件,对数据进行分析,在 α=0.05 下这五种方法在平均月推销额上有无显著差异;
- 哪种推销方法的效果最好?试对该种方法一个月的平均推销额求置信水平为 0.95 的置信区间。
解
(1) 由数据计算得
推销方法A1A2A3A4A5和Ti149.0172.1143.9127.6195.9788.5Ti222201.0029618.4120707.2116281.7638376.81127185.19∑j=1myij23243.34325.253030.992334.445505.0718439.05
ST=18439.05−35788.52=675.2714,fT=34,
SA=7127185.19−35788.52=405.5343,fA=4,
Se=ST−SA=269.7371,fe=30.
方差分析表为
来源因子 A误差 e总计平方和405.5343269.7371675.2714自由度43034均方101.38368.9912F 比11.2759
在显著性水平 α=0.05 下,查表得 F0.95(4,30)=2.69。由于 F=11.2759>2.69,故认为五种不同推销方法在平均月推销额上有显著差异。检验的 p 值为
p=1−fcdf(11.2759,4,30)=1.0526×10−5.
(2) 各种方法平均月推销额的估计值分别为
μ^1=21.2857,μ^2=24.5857,μ^3=20.5571,μ^4=18.2286,μ^5=27.9857.
从点估计看,第 5 种方法最好。又有
σ^2=MSe=8.9912,σ^=8.9912=2.9985.
查表得
t0.975(30)=2.0423,σ^t0.975(30)/7=2.3146,
于是第 5 种方法下均值的 0.95 置信区间为
μ5:[27.9857±2.3146]=[25.6711,30.3003].
某组装产品内有部分噪声很大的次品,产生次品的原因怀疑是由于这种组装品某个部位的间隙过大引起的。为检验这一认识是否正确,特从正品 A1 和次品 A2 中各抽出 8 个,对其间隙进行测量,数据(单位:μm)如下:
A1A2478122133155841061279
在正态分布假设下,请对 A1 与 A2 中间隙的均值是否存在显著差异进行检验(取 α=0.05)。
解
这是单因子二水平等重复试验,其均值比较可用两种方法。
解法一 方差分析法
先计算
T1=39,T2=86,T=T1+T2=125.
再得
ST=1195−161252=218.4375,
SA=8392+8862−161252=138.0625,
Se=ST−SA=80.375.
方差分析表为
来源因子 A误差 e总计平方和138.062580.375218.4375自由度11415均方138.06255.7411F 比24.0481
当 α=0.05 时,F0.95(1,14)=4.60。由于 F>4.60,故因子 A 显著,即正品与次品该部位的平均间隙有显著差异。
解法二 双样本 t 检验
在正态总体方差相等的条件下,检验统计量为
t=swm11+m21yˉ1−yˉ2∼t(m1+m2−2),
其中
sw2=m1+m2−21[i=1∑m1(y1i−yˉ1)2+i=1∑m2(y2i−yˉ2)2].
由样本可算得
m1=m2=8,yˉ1=4.875,yˉ2=10.75,sw=2.3961,
t=2.396181+814.875−10.75=−4.9039.
对给定显著性水平 α=0.05,拒绝域为
W={∣t∣>t1−α/2(m1+m2−2)}.
查表得
t0.975(14)=2.1448.
由于 ∣t∣>2.1448,故应拒绝两均值相等的假设,此结论与方差分析相同。这并非偶然,因为自由度为 f2 的 t 变量的平方服从 F(1,f2) 分布,本题中
(−4.9039)2=24.0481.
某乳制品公司有四个车间生产同一种酸乳酪,为考察四个车间产品中脂肪含量是否一致,特在每个车间生产的产品中各抽取 8 个样品送往实验室测定,结果(单位:%)如下:
A1A2A3A43.112.943.182.843.363.043.262.953.243.133.482.732.962.863.353.183.153.273.303.043.183.193.062.903.303.103.243.083.062.993.412.98
试比较各车间生产的酸乳酪中脂肪含量均值有无显著差异(取 α=0.05)。
解
为简化运算,令
zij=(yij−3)×100,
则
车间A1A2A3A4\multicolumn9c∣T=386\multicolumn8c∣zij11−618−16\multicolumn2c∑∑zij2=14342Ti36426−5∑j=18zij2241348−27−4−143518152730418196−1030102486−141−213652228−303494160877221518
于是
ST=14342−323862=9685.875,
SA=81i=1∑4Ti2−323862=4604.375,
Se=ST−SA=5081.500.
方差分析表为
来源因子 A误差 e总计平方和4604.3755081.5009685.875自由度32831均方1534.7917181.4821F 比8.46
查表得 F0.95(3,28)≈2.92。由于 F=8.46>2.92,故因子 A 显著,即四个车间生产的酸乳酪中脂肪含量均值有显著差异。
由逆变换
yij=3+100zij
可知
yˉ1=3.17,yˉ2=3.065,yˉ3=3.285,yˉ4=2.963,
σ^y2=10000MSe=0.01814821,σ^y=0.1347.
§8.2 多重比较
- 问题 在单因子方差分析中,当因子 A 显著时,就要进一步比较多个水平均值中任意两个水平之间有无显著差异。这类同时比较多个均值的问题称为多重比较。
若因子 A 有 r 个水平,则同时检验
H0ij:μi=μj,1≤i<j≤r,
其拒绝域可写为
W=1≤i<j≤r⋃{∣yˉi−yˉj∣≥cij}.
给定显著性水平 α (0<α<1) 后,诸临界值 cij 由 P(W)=α 决定。
- 重复数相等场合的 T 法 当各水平试验次数相同时,诸临界值相同:
c=q1−α(r,fe)σ^/m,σ^=MSe,
其中 q1−α(r,fe) 是分布 q(r,fe) 的 1−α 分位数,可由附表 8 查得。
- 重复数不等场合的 S 法 当各水平试验次数不同时,临界值 cij 不同:
cij=(r−1)F1−α(r−1,fe)(mi1+mj1)σ^2,
其中 F1−α(r−1,fe) 为分布 F(r−1,fe) 的 1−α 分位数。
习题与解答 8.2
采用习题 8.1 第 7 题的数据,对三种储藏方法的平均含水率在 α=0.05 下作多重比较。
解
由于习题 8.1 第 7 题中已经指出储藏方法因子是显著的,因此可以作多重比较。各水平下试验次数相同,均为 5,可采用重复数相等场合的 T 法。若取 α=0.05,查表得
q1−0.05(3,12)=3.77,σ^=0.686=0.828,
所以
c=3.77×0.828/5=1.396.
比较各均值差得
∣yˉ1−yˉ2∣=∣7.98−6.4∣=1.58>1.396,
认为 μ1 与 μ2 有显著差别;
∣yˉ1−yˉ3∣=∣7.98−9.12∣=1.14<1.396,
认为 μ1 与 μ3 无显著差别;
∣yˉ2−yˉ3∣=∣6.4−9.12∣=2.72>1.396,
认为 μ2 与 μ3 有显著差别。
采用习题 8.1 第 8 题的数据,对五种推销方法在 α=0.05 下作多重比较。
解
这里各水平下试验次数相同,均为 7,在推销因子显著的前提下,采用重复数相等场合的 T 法。取 α=0.05,查表得
q1−0.05(5,30)=4.10,σ^=8.9912=2.9985,
所以
c=4.10×2.9985/7=4.6466.
于是
∣yˉ1−yˉ2∣=∣21.2857−24.5857∣=3.300<4.6466,
∣yˉ1−yˉ3∣=∣21.2857−20.5571∣=0.7286<4.6466,
∣yˉ1−yˉ4∣=∣21.2857−18.2286∣=3.0571<4.6466,
∣yˉ1−yˉ5∣=∣21.2857−27.9857∣=6.700>4.6466,
∣yˉ2−yˉ3∣=∣24.5857−20.5571∣=4.0286<4.6466,
∣yˉ2−yˉ4∣=∣24.5857−18.2286∣=6.3571>4.6466,
∣yˉ2−yˉ5∣=∣24.5857−27.9857∣=3.400<4.6466,
∣yˉ3−yˉ4∣=∣20.5571−18.2286∣=2.3285<4.6466,
∣yˉ3−yˉ5∣=∣20.5571−27.9857∣=7.4286>4.6466,
∣yˉ4−yˉ5∣=∣18.2286−27.9857∣=9.7571>4.6466.
因此,在显著性水平 0.05 下,第 1,3,4 种方法与第 5 种有显著差别,第 2 种与第 4 种也有显著差别,其余 6 组均无显著差别。
有七种人造纤维,每种抽 4 根测其强度,得每种纤维的平均强度及标准差如下:
iyˉisi16.30.8126.20.9236.71.2246.80.7456.50.8867.00.5877.11.05
假定各种纤维的强度服从等方差的正态分布。
- 试问七种纤维强度间有无显著差异(取 α=0.05);
- 若各种纤维强度间无显著差异,则给出平均强度的置信水平为 0.95 的置信区间;若有显著差异,请进一步在 α=0.05 下进行多重比较,并指出哪种纤维的平均强度最大,同时给出该种纤维平均强度的置信水平为 0.95 的置信区间。
解
(1) 这是一个方差分析问题。由题设可算得
yˉ=71(6.3+6.2+⋯+7.1)=6.6571.
于是
SA=4(6.32+6.22+⋯+7.12)−28×6.65712=2.8045,fA=6,
Se=3(0.812+0.922+⋯+1.052)=17.2554,fe=21.
从而
MSA=62.8045=0.4674,MSe=2117.2554=0.8217,
F=MSeMSA=0.82170.4674=0.5688.
检验的 p 值为
p=1−fcdf(0.5688,6,21)=0.7505.
故认为七种纤维强度间无显著差异。
(2) 由于方差分析结论不显著,应把全部数据看作来自同一个总体。这里
ST=SA+Se=2.8045+17.2554≈20.0608,
从而误差方差的无偏估计为
σ^2=28−1ST=2720.0608=0.743,σ^=0.743=0.862.
平均强度的估计为
μ^=yˉ=6.6571.
若取 α=0.05,则
t0.975(27)=2.0518,σ^t0.975(27)/28=0.3342,
于是平均强度的 0.95 置信区间为
[6.6571±0.3342]=[6.3229,6.9913].
一位经济学家对生产电子计算机设备的企业收集了在一年内生产力提高指数(用 0 到 100 内的数表示),并按过去三年间在科研和开发上的平均花费分为三类:
A1:花费少,A2:花费中等,A3:花费多.
生产力提高指数如下表所示:
水平A1A2A3\multicolumn12c生产力提高指数7.66.78.58.28.19.76.89.410.15.88.67.86.97.89.66.67.79.56.38.97.77.96.08.38.77.18.4
试列出方差分析表,并进行多重比较(取 α=0.05)。
解
由数据计算如下:
水平A1A2A3和mi9126n=27Ti61.997.655.2T=214.7Ti2/mi425.734793.813507.84i=1∑rmiTi2=1727.387∑j=1miyij2431.03800.12511.60i=1∑rj=1∑miyij2=1742.75
于是
ST=1742.75−27214.72=35.487,fT=26,
SA=1727.387−27214.72=20.124,fA=2,
Se=35.487−20.124=15.363,fe=24.
方差分析表为
来源因子 A误差 e总和 T平方和20.12415.36335.487自由度22426均方10.0620.64F 比15.72
若取 α=0.05,查表得 F0.95(2,24)=3.41。由于 F=15.72>3.41,故认为各水平间有显著差异。检验的 p 值为
p=1−fcdf(15.72,2,24)=4.3317×10−5.
由于各水平试验次数不等,可采用重复数不等场合的 S 法。此时
F0.95(2,24)=3.41,σ^2=0.64,
故
c12=2×3.41(91+121)×0.64=0.9213,
c13=2×3.41(91+61)×0.64=1.1011,
c23=2×3.41(121+61)×0.64=1.0446.
又因为
∣yˉ1−yˉ2∣=∣6.88−8.13∣=1.25>0.9213,
∣yˉ1−yˉ3∣=∣6.88−9.2∣=2.32>1.1011,
∣yˉ2−yˉ3∣=∣8.13−9.2∣=1.07>1.0446,
故在显著性水平 0.05 下,各个水平间均有显著差异,第 3 个水平(花费多)对生产力提高最有帮助。
某工厂为了研究本厂生产的垫片与国内外同类产品在耐磨性能上的差别,特选国外一家产品、国内两家产品与本厂产品进行磨损试验。试验数据用磨损率表示,它是愈小愈好,磨损率的计算公式为
γ=初始值初始值−300 小时后的值.
具体数据如下:
水平A1(国外产品)A2(国内甲厂产品)A3(国内乙厂产品)A4(本厂产品)\multicolumn11c∣T=513\multicolumn10c∣磨损率9262420\multicolumn2cn=26Ti12192518mi1426191915282217232715252416131822175014714117546610
试在正态分布假设下对比四家同类产品的磨损率均值有无显著差异;若有显著差异,再作多重比较(取 α=0.05)。
解
先计算各平方和:
ST=i=1∑4j=1∑miyij2−nT2=10769−265132=647.12,
SA=i=1∑4miTi2−nT2=4502+61472+61412+101752−265132=480.62,
Se=ST−SA=166.50.
方差分析表为
来源因子 A误差 e总和 T平方和480.62166.50647.12自由度32225均方160.217.57F 比21.16
在显著性水平 α=0.05 下,F0.95(3,22)=3.05。由于 F>3.05,故因子 A 显著,即四个工厂的磨损率均值间存在显著差异。由于各水平重复数不同,选用 S 法进行多重比较:
cij=[(r−1)F1−α(r−1,fe)MSe(mi1+mj1)]1/2=8.32mi1+mj1.
进一步有
c12=c13=8.3241+61=5.37,
c14=8.3241+101=4.92,
c23=8.3261+61=4.80,
c24=c34=8.3261+101=4.30.
最后比较可得
∣yˉ1−yˉ2∣=12>5.37,∣yˉ1−yˉ3∣=11>5.37,∣yˉ1−yˉ4∣=5>4.92,
∣yˉ2−yˉ3∣=1<4.80,∣yˉ2−yˉ4∣=7>4.30,∣yˉ3−yˉ4∣=6>4.30.
可见,除 A2 与 A3 间无显著差异外,其他水平间都有显著差异,可分为三组:
{A1}, {A4}, {A2,A3}.
组内无显著差异,组间都有显著差异。磨损率愈小愈好,而四个水平的均值估计分别为
yˉ1=12.5,yˉ2=24.5,yˉ3=23.5,yˉ4=17.5,
因此 A1(国外产品)最好,其次是 A4(本厂产品),A2 与 A3 较差。
某厂有四条生产线生产同一种垫片,为了比较它们的断裂强度有无显著差异,特从每条生产线上随机抽取 5 个垫片,测其断裂强度,数据如下:
生产线A1A2A3A4\multicolumn5c断裂强度 yij86.593.488.694.392.087.993.293.385.290.688.892.087.985.592.789.286.088.490.992.5
试在正态分布假设下比较四条生产线上的产品断裂强度;若有显著差异,再作多重比较(取 α=0.05)。
解
为便于计算,把各数据 yij 均减去 85,记
zij=yij−85,
则有
生产线A1A2A3A4\multicolumn6c∣T=98.9\multicolumn5c∣zij1.58.43.69.3Ti7.02.98.28.30.25.63.87.02.90.57.74.21.03.45.97.512.620.829.236.3
利用这些数据可算得
ST=i=1∑4j=1∑5zij2−20T2=160.7895,
SA=i=1∑45Ti2−20T2=63.2855,
Se=ST−SA=97.5040.
于是方差分析表为
来源因子 A误差 e总和 T平方和63.285597.5040160.7895自由度31619均方21.09526.0940F 比3.46
对给定显著性水平 α=0.05,查表得 F0.95(3,16)=3.24。由于 F=3.46>3.24,故因子 A 显著,即四条生产线上垫片断裂强度均值间有显著差异。其均值估计为
yˉ1=87.52,yˉ2=89.16,yˉ3=90.84,yˉ4=92.26.
由于各水平重复数相等,可采用 T 法进行多重比较。对显著性水平 α=0.05,临界值为
c=q0.95(4,16)mMSe=4.05×56.094=4.47.
比较各均值差,仅有
∣yˉ1−yˉ4∣=4.74>4.47,
其余各差值均小于 4.47。因此只有第 1 条和第 4 条生产线上的产品均值存在显著差异。断裂强度愈大愈好,所以第 4 条生产线上的产品质量最好,第 1 条最差。
§8.3 方差齐性检验
- 问题 方差齐性即诸方差相等,是方差分析的基本假定之一。方差齐性检验就是检验这一假定是否成立。其一对假设为
H0:σ12=σ22=⋯=σr2vsH1:诸 σi2 不全相等.
- 哈特利检验(在重复数相等场合使用) 在重复次数均为 m 时,检验统计量为
H=min{s12,s22,…,sr2}max{s12,s22,…,sr2},
其中 si2 为第 i 个水平下样本方差,拒绝域为
W={H≥H1−α(r,f)},f=m−1.
- 巴特利特检验(可在重复数不等场合使用,但样本量不得低于 5) 当重复次数不等且各水平下试验次数均不低于 5 时,可采用巴特利特检验。统计量为
B=C1(felnMSe−i=1∑rfilnsi2),
其中 fi=mi−1 为 si2 的自由度,
MSe=fe1i=1∑rfisi2,fe=i=1∑rfi,
而
C=1+3(r−1)1(i=1∑rfi1−fe1).
拒绝域为
W={B≥χ1−α2(r−1)}.
- 修正的巴特利特检验(在样本量相等或不等、样本量较小或较大均可使用) 一般场合可采用修正的巴特利特检验,其统计量为
B′=f1(A−BC)f2BC,
其中 B,C 同前,
f1=r−1,f2=(C−1)2r+1,A=2−C+2/f2f2.
拒绝域为
W={B′≥F1−α(f1,f2)}.
若 f2 不是整数,可通过对 F 分布分位数表作线性内插获得近似值。
习题与解答 8.3
采用教材中例 8.1.1 的数据,在显著性水平 α=0.05 下用哈特利检验考察三个总体方差是否彼此相等。
解
这是一个检验方差是否相等的问题。各数据计算如下:
水平A1A2A3\multicolumn8c∣数据(原始数据−1000)7310793Ti99229均值 yˉi60−1080∑j=18(yij−yˉi)2110921290221274329122292814819458535424.2573.12544.255319.517576.95319.5
于是
s12=75319.5=759.93,s22=717576.9=2510.99,s32=75319.5=759.93.
从而
H=759.932510.99=3.304.
当 α=0.05 时,由附表 10 查得 H0.95(3,7)=6.94。由于 H<6.94,故接受原假设 H0,认为三个总体方差间无显著差异。
在安眠药试验(见习题 8.1 第 5 题)中已求得四个样本方差
s12=0.02,s22=0.08,s32=0.036,s42=0.1307.
请用哈特利检验在显著性水平 α=0.05 下考察四个总体方差是否彼此相等。
解
这是关于方差齐性的检验问题。这里 r=4,m=6,已知统计量
H=0.020.1307=6.535.
当 α=0.05 时,查附表 10 知 H0.95(4,5)=13.7。由于 H<13.7,故应接受原假设 H0,认为四个总体方差间无显著差异。
在生产力提高的指数研究中(见习题 8.2 第 4 题)可求得三个样本方差:
s12=0.662,s22=0.573,s32=0.752.
请用巴特利特检验在显著性水平 α=0.05 下考察三个总体方差是否彼此相等。
解
由题设知各组样本量分别为 9,12,6,最小样本量大于 5,可采用巴特利特检验。此时
fe=8+11+5=24,
MSe=241(8×0.662+11×0.573+5×0.752)=0.64,
C=1+3(3−1)1(81+111+51−241)=1.0624.
因此
B=1.06241[24ln0.64−(8ln0.662+11ln0.573+5ln0.752)]=0.1315.
对显著性水平 α=0.05,查表知 χ0.952(3−1)=5.9915,拒绝域为
W={B≥5.9915}.
由于 B<5.9915,故应接受原假设 H0,认为三个总体的方差无显著差异。
在人员推销效果研究中(见习题 8.1 第 8 题),分别用哈特利检验和巴特利特检验在显著性水平 α=0.05 下对五个总体作方差齐性检验。
解
先用哈特利检验判断等方差性。通过习题 8.1 第 8 题的解答可知各组内平方和
Q1=71.7286,Q2=94.0486,Q3=72.8171,Q4=8.4743,Q5=22.6686.
利用公式 si2=Qi/fi 可得各组样本方差
s12=11.9548,s22=15.6748,s32=12.1362,s42=1.4124,s52=3.7781.
于是
H=1.412415.6748=11.098.
对显著性水平 α=0.05,由附表 10 查得 H0.95(5,6)=12.1。由于 H<12.1,所以应接受原假设 H0,即认为各个总体方差相等。
接下来计算巴特利特检验统计量。由习题 8.1 第 8 题已知 MSe=8.9912,且
C=1+3(5−1)1(65−301)=1.0667.
于是
B=1.06671[30ln8.9912−6(ln11.9548+ln15.6748+ln12.1362+ln1.4124+ln3.7781)]=8.8723.
对显著性水平 α=0.05,查表知 χ0.952(5−1)=9.4877。由于 B<9.4877,故也应接受原假设 H0。两种检验的结果是一致的。
在对粮食含水率的研究中(见习题 8.1 第 7 题)已求得 3 个水平下的组内平方和
Q1=0.988,Q2=3.860,Q3=3.388.
请用修正的巴特利特检验在显著性水平 α=0.05 下考察三个总体方差是否彼此相等。
解
由已给条件及每组样本量均为 5,可得三个样本方差
s12=0.2470,s22=0.965,s32=0.847.
且
C=1+3(3−1)1(43−121)=1.1111.
从而
B=1.11111[12ln0.686−4(ln0.247+ln0.965+ln0.847)]=1.6899.
进一步,
f1=3−1=2,f2=(1.1111−1)23+1=324.0648,
A=2−1.1111+2/324.0648324.0648=362.0546.
因此修正的巴特利特检验统计量为
B′=2(362.0546−1.6899×1.1111)324.0648×1.6899×1.1111=0.845.
对显著性水平 α=0.05,F0.95(2,324.0648)≈F0.95(2,∞)=3。由于 B′<3,故接受原假设 H0,认为三个水平下的方差间无显著差异。
针对食品包装研究的数据(见教材中例 8.1.4),请用修正的巴特利特检验在显著性水平 α=0.05 下考察四个总体是否满足方差齐性假定。
解
在例 8.1.4 中,r=4,各组样本量不相等,且样本量分别为 2,3,3,2,都不大,因此只能用修正的巴特利特检验。由例 8.1.4 数据可得各水平下样本方差
s12=18,s22=1,s32=4,s42=18,
并且已知 MSe=7.67。于是
C=1+3(4−1)1(1+21+21+1−61)=1.3148,
B=1.31481[6ln7.67−(ln18+2ln1+2ln4+ln18)]=2.79.
进一步,
f1′=4−1=3,f2′=(1.3148−1)24+1=50.45,
A=2−1.3148+2/50.4550.45=69.60.
因而修正的巴特利特检验统计量为
B′=3(69.6−2.79×1.3148)50.45×2.79×1.3148=0.9356.
若取显著性水平 α=0.05,则
F0.95(3,50.45)≈F0.95(3,60)=2.76.
由于 B′<2.76,故接受原假设,可以认为四个总体方差彼此相等。
在垫片的耐磨试验(见习题 8.2 补充习题 5)中,关于磨损率有四个样本,它们的样本方差 si2、样本量 ni、误差均方和 MSe 与其自由度 fe 分别为
n1=4, s12=7,n2=6, s22=9.9,n3=6, s32=7.5,n4=10, s42=6.5,
fe=22,MSe=7.57.
现要对“四个总体方差彼此相等”的假设作出判断。
解
由于四个样本量不全相等,其中有一个样本量小于 5,故选用修正的巴特利特检验。先计算中间结果:
C=1+3(4−1)1(31+51+51+91−221)=1.0888,
B=C1(felnMSe−i=1∑4filnsi2)=1.088822ln7.57−(3ln7+5ln9.9+5ln7.5+9ln6.5)=0.2857,
f1′=r−1=3,f2′=(C−1)2r+1=0.088825=634.08,
A=2−1.0888+2/634.08634.08=693.47.
于是
B′=f1′(A−BC)f2′BC=3(693.47−0.2857×1.0888)634.08×0.2857×1.0888=0.0949.
又因
F0.95(3,634.08)≈F0.95(3,∞)=2.61,
故 B′<2.61,应接受原假设,可以认为四个总体方差彼此相等。
杀伤性武器迫击炮弹的出口速度与炮口面积有关,现选择四个炮口面积作为因子 A 的四个水平,安排一个重复数为 8 的单因子试验,所测出口速度(已排序)如下:
j12345678Tiyˉisi2A1=0.01643.650.155.356.056.358.262.964.9447.355.9145.82A2=0.03045.746.346.948.151.054.554.955.7403.150.3917.42A3=0.04417.719.624.227.627.928.528.630.9205.025.6322.17A4=0.0587.410.614.314.616.016.022.122.2123.215.4025.82
- 对各水平下的数据作正态性检验;
- 对各水平下的数据作方差齐性检验;
- 计算各平方和并作方差分析;
- 若因子 A 显著,再作多重比较。
解
(1) 由于各水平下样本量均为 8,故选用 W 检验作正态性检验。其统计量定义为
ui=j=1∑4aj[yi(9−j)−yij],Qi=j=1∑8(yij−yˉi)2.
\begingroup
\footnotesize
\setlength{\tabcolsep}{4pt}
\renewcommand{\arraystretch}{1.2}
j1234uiQiWi=ui2/Qiyi(9−j)−yijyi8−yi1yi7−yi2yi6−yi3yi5−yi4A121.312.82.90.317.4630320.750.9508A210.08.67.62.910.2604121.950.8633A313.29.04.30.311.6026155.160.8676A414.811.51.71.412.9704180.740.9308aj0.60520.31640.17430.0561
\endgroup
对显著性水平 α=0.05,查附表 7 在 n=8 时 W0.05=0.818。可见诸 Wi 都超过 0.818,可以认为这四个样本都来自不同正态分布。
(2) 由于各水平下重复数相等,故可选用哈特利检验进行方差齐性检验,其统计量
H=min{s12,s22,s32,s42}max{s12,s22,s32,s42}=17.4245.82=2.63.
在显著性水平 α=0.05 下,查附表 10 得 H0.95(4,7)=8.44。由于 H<8.44,故可以认为四个水平下的正态总体方差彼此相等。
(3) 为进行方差分析,需计算各平方和:
ST=i=1∑4j=1∑8yij2−32T2=53249.86−321178.62=9840.55,
SA=81i=1∑4Ti2−32T2=9061.96,
Se=ST−SA=778.59.
把这些平方和移入方差分析表,继续计算得
来源因子 A误差 e总和 T平方和9061.96778.599840.55自由度32831均方3020.6527.81F 比108.62
对给定的显著性水平 α=0.05,查表得 F0.95(3,28)≈2.92。由于 F>2.92,故因子 A 显
著,即四个平均出口速度间有显著差异。
(4)由于各水平的重复数相等,故用 T 法进行多重比较。对给定显著性水平 α=0.05,查附表 8 得
q0.95(4,28)≈3.85,
临界值为
c=q0.95(4,28)mMSe≈3.85×827.81=7.18.
把 c 与各均值差的绝对值
∣yˉi−yˉj∣
分别进行比较,只有
∣yˉ1−yˉ2∣=5.52<c=7.18,
其他 5 个绝对值都大于 7.18,这表明 A1 与 A2 间无显著差异,其他 5 对水平间都有显著差异。
§8.4 一元线性回归
- 问题 考察两个变量 x 与 y 之间是否存在线性相关关系,其中 x 是一般(可控)变量,y 是随机变量,其线性回归关系可表示如下(可用散点图显示):
y=β0+β1x+ε,
其中 β0 为截距,β1 为斜率,ε 为随机误差,常假设
ε∼N(0,σ2).
这里 β0,β1,σ2 是三个待估参数.上式表明,y 与 x 之间有线性关系,但受到随机误差的干扰.
- 数据 对 x 与 y 通过试验或观察可得 n 对数据(注:数据是成对的,不允许错位).在 y 与 x 之间存在线性关系的假设下,有如下统计模型:
{yi=β0+β1xi+εi,各 εi 独立同分布,其分布为 N(0,σ2).i=1,2,…,n,
利用成对数据可获得 β0 与 β1 的估计,设估计分别为 β^0,β^1,则称
y^=β^0+β^1x
为回归方程,其图形称为回归直线.
- 参数估计 用最小二乘法可得 β0 与 β1 的无偏估计
{β^1=lxy/lxx,β^0=yˉ−β^1xˉ,
其中
xˉ=n1∑xi,yˉ=n1∑yi,
lxy=∑(xi−xˉ)(yi−yˉ)=∑xiyi−nxˉyˉ=∑xiyi−n1∑xi∑yi,
lxx=∑(xi−xˉ)2=∑xi2−nxˉ2=∑xi2−n1(∑xi)2,
lyy=∑(yi−yˉ)2=∑yi2−nyˉ2=∑yi2−n1(∑yi)2.
- 回归方程的显著性检验 回归方程的显著性检验就是要对如下一个假设作出判断:
H0:β1=0vsH1:β1=0.
对此可采用如下两种等价的检验方法:
(1)F 检验
如下的平方和分解式是非常重要的,它在许多统计领域都有应用:
ST=SR+Se,fT=fR+fe,
其中
ST=∑(yi−yˉ)2=lyy
是总偏差平方和,其自由度 fT=n−1;
SR=∑(y^i−yˉ)2=∑(β^0+β^1xi−yˉ)2=β^1lxy=β^12lxx
是回归平方和,其自由度 fR=1;
Se=∑(yi−y^i)2=∑(yi−β^0−β^1xi)2
是残差平方和,其自由度 fe=n−2.
而 y^i=β^0+β^1xi 是在 x=xi 时的回归值(拟合值),它与实测值 yi 通常是不相等的.
在原假设 H0 成立的条件下,检验统计量
F=Se/(n−2)SR∼F(1,n−2),
拒绝域为
W={F≥F1−α(1,n−2)}.
上述检验过程一般用如下方差分析表列出:
来源回归残差总计平方和SRSeST自由度fR=1fe=n−2fT=n−1均方MSR=SRMSe=n−2SeF 比F=MSeMSR
(2)t 检验
检验统计量为
t=σ^/lxxβ^1,σ^=Se/(n−2).
在原假设成立下,
t∼t(n−2),
因此拒绝域为
W={∣t∣≥t1−α/2(n−2)}.
注意到
t2=F,
因此 t 检验与 F 检验是等价的,选其中之一使用即可.
- 相关系数及其检验
(1)相关系数
对容量为 n 的二维样本
{(xi,yi), i=1,2,…,n}
的线性相关程度可用如下(样本)相关系数量
r=∑(xi−xˉ)2∑(yi−yˉ)2∑(xi−xˉ)(yi−yˉ)=lxxlyylxy
来衡量.
- r=±1,n 个点完全在一条直线上,此时两者之间可能是确定性关系;
- r>0,当 x 增加时,y 有线性增加趋势,此时称正相关;
- r<0,当 x 增加时,y 反而有线性减少趋势,此时称负相关;
- r=0,n 个点可能杂乱无章,也可能呈某种曲线趋势,此时称不(线性)相关.
(2)相关系数的检验
记 ρ 为二维总体的相关系数,于是可建立如下假设:
H0:ρ=0vsH1:ρ=0.
对此,采用检验统计量 r=lxy/lxxlyy,拒绝域为
W={∣r∣≥r1−α(n−2)},
其中 r1−α(n−1) 是 ∣r∣ 分布的 1−α 分位数,可查附表 9.
(3)检验统计量 r 与 F 统计量之间关系
r2=F+(n−2)F.
这表明 ∣r∣ 是 F 的严格增函数,所以相关系数检验与前面的 F 检验也是等价的.
- 估计与预测 回归方程的应用
- 当 x=x0 时,y^0=β^0+β^1x0 是 E(y0)=β0+β1x0 的点估计;
- 当 x=x0 时,E(y0)=β0+β1x0 的置信水平为 1−α 的置信区间是
[y^0−δ0, y^0+δ0],
其中
δ0=t1−α/2(n−2)σ^n1+lxx(x0−xˉ)2,σ^=MSe;
- 当 x=x0 时,y0=β0+β1x0+ε0 的 1−α 预测区间是
[y^0−δ, y^0+δ],
其中
δ=δ(x0)=t1−α/2(n−2)σ^1+n1+lxx(x0−xˉ)2.
注:E(y0) 是未知参数,而 y0 是随机变量.对 E(y0) 谈论的是置信区间,对 y0 谈论的是预测区间,两者是不同的,显然,预测区间要比置信区间宽很多.
要提高预测区间(置信区间也一样)的精度,即要使 δ(或 δ0)较小,这要求:
- 增大样本量 n;
- 增大 lxx,即要求 x1,x2,…,xn 较为分散;
- 使 x0 靠近 xˉ.
习题与解答 8.4
假设回归直线过原点,即一元线性回归模型为
yi=βxi+εi,i=1,2,…,n,
E(εi)=0,Var(εi)=σ2,
诸观测值相互独立.
- 写出 β 的最小二乘估计和 σ2 的无偏估计;
- 对给定的 x0,其对应的因变量均值的估计为 y^0,求 Var(y^0).
解
由最小二乘法原理,令
Q=i=1∑n(yi−βxi)2,
则正则方程为
∂β∂Qβ^=−2i=1∑n(yi−β^xi)xi=0.
从中解得 β 的最小二乘估计为
β^=∑i=1nxi2∑i=1nxiyi.
不难看出
E(β^)=β,Var(β^)=∑i=1nxi2σ2.
于是,由
Se=∑(yi−y^i)2=∑(βxi+εi−β^xi)2=∑[xi2(β^−β)2+εi2−2(β^−β)xiεi],
可得
E(Se)=∑xi2Var(β^)+nVar(ε)−2∑xiE(β^εi).
将 β^ 写成 y1,y2,…,yn 的线性组合,利用 yj 与 εi(i=j)间的独立性,有
E(β^εi)=∑j=1nxj2xiσ2.
由此即有
∑xiE(β^εi)=σ2,
从而
E(Se)=(n−1)σ2.
这给出 σ2 的无偏估计为
σ^2=n−11Se.
对给定的 x0,对应的因变量均值的估计为
y^0=β^x0,
于是
Var(y^0)=x02Var(β^)=∑i=1nxi2x02σ2.
设回归模型为
{yi=β0+β1xi+εi,各 εi 独立同分布,其分布为 N(0,σ2),i=1,2,…,n,
试求 β0,β1 的最大似然估计,它们与其最小二乘估计一致吗?
解
似然函数为
L=(2πσ1)nexp{−2σ2∑i=1n(yi−β0−β1xi)2},
其对数似然函数为
l=−2nlnσ2−2σ21i=1∑n(yi−β0−β1xi)2
(忽略常数项).将其分别对 β0,β1 求导,并令导函数为 0,得到如下似然方程组:
⎩⎨⎧∂β0∂lβ^0=∑i=1n(yi−β^0−β^1xi)=0,∂β1∂lβ^1=∑i=1n(yi−β^0−β^1xi)xi=0.
经过整理可以解出
{β^1=lxy/lxx,β^0=yˉ−β^1xˉ.
可以看到 β0,β1 的最大似然估计与其最小二乘估计是一致的.
在回归分析计算中,常对数据进行变换
y~i=d1yi−c1,x~i=d2xi−c2,i=1,2,…,n,
其中 c1,c2,d1(d1>0),d2(d2>0) 是适当选取的常数.
- 试建立由原始数据和变换后数据得到的最小二乘估计、总平方和、回归平方和以及残差平方和之间的关系;
- 证明:由原始数据和变换后数据得到的 F 检验统计量的值保持不变.
解
经变换后,各平方和的表达式如下:
xˉ~=n1∑x~i=d21(xˉ−c2),yˉ~=n1∑y~i=d11(yˉ−c1),
l~x~y~=∑(x~i−xˉ~)(y~i−yˉ~)=d1d21lxy,
l~x~x~=∑(x~i−xˉ~)2=d221lxx,l~y~y~=∑(y~i−yˉ~)2=d121lyy.
所以由原始数据和变换后数据得到的最小二乘估计间的关系为
β~1=l~x~x~l~x~y~=d1d2β^1,β~0=yˉ~−β~1xˉ~=d11β^0−d11(c1−β^1c2).
在实际应用中,人们往往在先由变换后的数据求出 β~1,β~0,然后再据此给出 β^1,β^0,它们的关系为
β^1=d2d1β~1,β^0=d1β~0+c1(1−d2/c2d1/c1β~1).
总平方和、回归平方和以及残差平方和分别为
ST=lyy=d12l~y~y~=d12S~T,
SR=β^12lxx=d22d12β~12⋅d22lx~x~=d12S~R,
Se=d12S~e.
由此知道
F=Se/(n−2)SR=S~e/(n−2)S~R=F~,
即说明了由原始数据和变换后数据得到的 F 检验统计量的值保持不变.
对给定的 n 组数据 (xi,yi),i=1,2,⋯,n,若我们关心的是 y 如何依赖 x 的取值而变动,则可以建立回归方程
y^=a+bx.
反之,若我们关心的是 x 如何依赖 y 的取值而变动,则可以建立另一个回归方程
x^=c+dy.
试问这两条直线在直角坐标系中是否重合?为什么?若不重合,它们有无交点?若有,试给出交点的坐标.
解
一般不重合.因为回归方程 y^=a+bx 可化为
y^−yˉ=lxxlxy(x−xˉ),
而 x^=c+dy 化为
x^−xˉ=lyylxy(y−yˉ).
当且仅当
lxy2=lxxlyy
时两条直线重合.我们知道,lxy2=lxxlyy 表示相关系数的绝对值为 1,即 n 组数据 (xi,yi),i=1,2,⋯,n 在一条直线上,这在实际中极其罕见,所以说“一般不重合”.
注:不重合时,它们一定有交点 (xˉ,yˉ).
为考察某种维尼纶纤维的耐水性能,安排了一组试验,测得其中醇浓度 x 及相应的“缩醇化度” y 数据如下:
xy1826.862028.352228.752428.872629.752830.003030.36
- 作散点图;
- 求样本相关系数;
- 建立一元线性回归方程;
- 对建立的回归方程作显著性检验(α=0.01).
解
(1)散点图如图 8.1,y 有随着 x 增加而增加趋势.
\FigureEightOne
(2)由样本数据可算得
∑xi=168,∑yi=202.94,
lxx=∑(xi−xˉ)2=112,lyy=∑(yi−yˉ)2=8.4931,
lxy=∑(xi−xˉ)(yi−yˉ)=29.6.
因此样本相关系数
r=lxxlyylxy=112×8.493129.6=0.9597.
(3)应用最小二乘估计公式,
β^1=lxxlxy=11229.6=0.2643,β^0=yˉ−β^1xˉ=22.6482,
于是,一元线性回归方程为
y^=22.6482+0.2643x.
(4)首先计算几个平方和
ST=lyy=8.4931,SR=β^12lxx=0.26432×112=7.8237,
Se=ST−SR=0.6694.
将各平方和移入方差分析表,继续计算,可以得到
来源回归残差总计平方和7.82370.66948.4931自由度156均方7.82370.1339F 比58.43
若取 α=0.01,查表知
F0.99(1,5)=16.26<58.43,
拒绝域为
W={F≥16.26},
现检验统计量值落入拒绝域,因此在显著性水平 0.01 下回归方程是显著的.此处,回归方程显著性检验的 p 值为(用 MATLAB 语句表示)
p=1−fcdf(58.43,1,5)=0.0006.
测得一组弹簧形变 x(单位:cm)和相应的外力 y(单位:N)数据如下:
yx13.081.23.761.44.311.65.021.85.512.06.252.26.742.47.402.88.543.09.24
由胡克定律知 y^=kx,试估计 k,并在 x0=2.6 cm 处给出相应的外力 y0 的 0.95 预测区间.
解
由本节的第 1 题可以给出 k 的最小二乘估计为
k^=∑xi2∑xiyi=395.3199128.296=0.3245.
在第 1 题中已经给出 k^ 的均值和方差分别为 k 和 σ2/∑xi2,所以
k^x0∼N(E(y0),∑xi2x02σ2),
其中
E(y0)=kx0,y0∼N(E(y0),σ2),
且两者独立,从而有
y0−y^0∼N(0,(1+∑xi2x02)σ2).
因此 y0 的预测区间为
(y^0−δ, y^0+δ),
其中
δ=t1−α/2(n−1)σ^1+∑xi2x02,σ^=n−1∑(yi−k^xi)2.
由样本数据可计算得到
Se=∑(k^xi−yi)2=k^2∑xi2−2k^∑xiyi+∑yi2
=0.32452×395.3199−2×0.3245×128.296+41.64=0.0032,
从而
σ^=0.0032/(10−1)=0.0189.
而 x0=2.6 cm 对应的外力的预测值为
y^0=0.3245×2.6=0.8437,
当 α=0.05 时,查表知
t0.975(9)=2.2622,
故
δ=2.2622×0.0189×1+395.31992.62=0.0431.
因而得到 y0 的预测区间为
[0.8437−0.0431, 0.8437+0.0431]=[0.8006, 0.8868].
设由 (xi,yi) (i=1,2,⋯,n) 可建立一元线性回归方程,y^i 是由回归方程得到的拟合值,证明:样本相关系数 r 满足关系
r2=∑i=1n(yi−yˉ)2∑i=1n(y^i−yˉ)2,
上式也称为回归方程的决定系数.
解
因为
SR=β^12lxx=lxxlxy2,
将之代入样本相关系数 r 的表达式中,即有
r2=lxxlyylxy2=lxxlyySRlxx=STSR=∑i=1n(yi−yˉ)2∑i=1n(y^i−yˉ)2.
证明完成.
现收集了 16 组合金钢中的碳含量 x 及强度 y 的数据,求得
xˉ=0.125,yˉ=45.7886,lxx=0.3024,lxy=25.5218,lyy=2432.4566.
- 建立 y 关于 x 的一元线性回归方程 y^=β^0+β^1x;
- 写出 β^0 和 β^1 的分布;
- 求 β^0 和 β^1 的相关系数;
- 列出对回归方程作显著性检验的方差分析表(α=0.05);
- 给出 β1 的 0.95 置信区间;
- 在 x0=0.15 时求对应的 y0 的 0.95 预测区间.
解
(1)根据已知数据可以得到回归系数的估计为
β^1=lxxlxy=0.302425.5218=84.3975,β^0=yˉ−β^1xˉ=35.2389.
于是 y 关于 x 的一元线性回归方程为
y^=35.2389+84.3975x.
(2)我们知道
β^0∼N(β0,(n1+lxxxˉ2)σ2),β^1∼N(β1,lxxσ2).
利用已给数据可计算出
lxx1=0.30241=3.3069,n1+lxxxˉ2=161+0.30240.1252=0.1142,
由此可得到 β^0,β^1 的分布分别为
β^0∼N(β0,0.1142σ2),β^1∼N(β1,3.3069σ2).
(3)由于
Cov(β^0,β^1)=−lxxxˉσ2=−0.30240.125σ2=−0.4134σ2,
故 β^0 和 β^1 的相关系数为
rβ^0,β^1=Var(β^0)Var(β^1)Cov(β^0,β^1)=0.1142×3.3069−0.4134=−0.6727.
(4)首先计算三个平方和
ST=lyy=2432.4566,SR=lxxlxy2=2153.9758,Se=ST−SR=278.4808.
于是可建立如下方差分析表:
来源回归残差总计平方和2153.9758278.48082432.4566自由度11415均方2153.975819.8915F 比108.2862
若取显著性水平 α=0.05,查表知
F0.95(1,14)=4.60,
拒绝域为
W={F≥4.60},
此处检验统计量落入拒绝域,因此,在显著性水平 0.05 下回归方程是显著的.此处,回归方程显著性检验的 p 值为
p=1−fcdf(108.2862,1,14)=5.6929×10−8.
(5)由教材中定理 8.4.1 与定理 8.4.3 知,
β^1∼N(β1,lxxσ2),σ2Se∼χ2(n−2),
且 β^1 与
σ^2=n−21Se
相互独立,因此 β1 的置信区间为
(β^1−t1−α/2(n−2)lxxσ^, β^1+t1−α/2(n−2)lxxσ^).
其中
σ^=19.8915=4.46,t0.975(14)=2.1448,
由此可得到 β1 的置信区间为
[84.3975−2.1448×4.46/0.3024, 84.3975+2.1448×4.46/0.3024]=[67.0022, 101.7928].
(6)首先算出 x0=0.15 对应的 y0 的预测值为
y^0=35.2389+84.3975×0.15=47.8985.
而
δ=δ(x0)=t1−α/2(n−2)σ^1+n1+lxx(x0−xˉ)2
=2.1448×4.46×1+161+0.3024(0.15−0.125)2=9.8698,
所以 x0=0.15 时求对应的 y0 的 0.95 预测区间为
[47.8985−9.8698, 47.8985+9.8698]=[38.0287, 57.7683].
设回归模型为
{yi=β0+β1xi+εi,εi∼N(0,σ2),
现收集了 15 组数据,计算有
xˉ=0.85,yˉ=25.60,lxx=19.56,lxy=32.54,lyy=46.74,
后经校对,发现有一组数据记录错误,正确数据为 (1.2,32.6),记录为 (1.5,32.3).
- 求 β0,β1 修正后的 LSE;
- 对回归方程作显著性检验(α=0.05);
- 若 x0=1.1,给出对应响应变量的 0.95 预测区间.
解
由于有一组数据记录错误,应将 xˉ,yˉ,lxx,lyy,lxy 作修正,修正后的量分别记为
xˉ′, yˉ′, lxx′, lyy′, lxy′,
则
xˉ′=xˉ+151(1.2−1.5)=0.83,yˉ′=yˉ+151(32.6−32.3)=25.62,
lxx′=lxx+nxˉ2−nxˉ′2−1.52+1.22=19.254,
lyy′=lyy+nyˉ2−nyˉ′2−32.32+32.62=50.844,
lxy′=lxy+nxˉyˉ−nxˉ′yˉ′−1.5×32.3+1.2×32.6=30.641.
根据修正后的数据可计算得到 β0,β1 的 LSE 为
β^1=lxx′lxy′=1.5914,β^0=yˉ′−β^1xˉ′=24.2991.
利用修正后的数据计算三个平方和为
ST=lyy′=50.844,fT=14,
SR=lxx′(lxy′)2=48.7624,fR=1,
Se=ST−SR=2.0816,fe=13,
从而检验统计量
F=MSeMSR=2.0816/1348.7624=304.5746.
若取显著性水平 α=0.05,查表知
F0.95(1,13)=4.67,
拒绝域为
W={F≥4.67},
由于检验统计量落入拒绝域,因此回归方程是显著的.此处,回归方程显著性检验的 p 值为
p=1−fcdf(304.5746,1,13)=2.1043×10−10.
对于 x0=1.1,其对应响应变量的预测值为
y^0=24.2991+1.5914×1.1=26.0496,
而
σ^=0.1601=0.4001,t0.975(13)=2.1604,
δ(x0)=2.1604×0.40011+151+19.254(1.1−0.83)2=0.8943,
因此响应变量的 0.95 预测区间为
[26.0496−0.8943, 26.0496+0.8943]=[25.1553, 26.9439].
在生产中积累了 32 组某种铸件在不同腐蚀时间 x 下腐蚀深度 y 的数据,求得回归方程为
y^=−0.4441+0.002263x,
且误差方差的无偏估计为
σ^2=0.001452,
总偏差平方和为 0.1246.
- 对回归方程作显著性检验(α=0.05),列出方差分析表;
- 求样本相关系数;
- 若腐蚀时间 x=870,试给出 y 的 0.95 近似预测区间.
解
(1)由已知条件可以得到
ST=0.1246,Se=(n−2)σ^2=0.04356,
因此
SR=0.1246−0.04356=0.08104.
把这些平方和移至如下方差分析表上,继续计算:
来源回归残差总计平方和0.081040.043560.1246自由度13031均方0.081040.001452F 比55.8127
若取显著性水平 α=0.05,则
F0.95(1,30)=4.17<55.8127,
因此回归方程是显著的,此处回归方程检验的 p 值为
p=1−fcdf(55.8127,1,30)=2.5101×10−8.
(2)样本相关系数
r=STSR=0.12460.08104=0.8065.
(3)若腐蚀时间 x=870,则 y 的预测值为
y^=−0.4441+0.002263×870=1.5247.
其 0.95 近似预测区间的半径为
δ≈σ^u0.975=0.001452×1.96=0.0747,
从而 y 的 0.95 近似预测区间为
[1.5247−0.0747, 1.5247+0.0747]=[1.4500, 1.5994].
我们知道营业税税收总额 y 与社会商品零售总额 x 有关.为能从社会商品零售总额去预测税收总额,需要了解两者之间的关系.现收集了如下 9 组数据(单位:亿元):
序号xy1142.083.932177.305.963204.687.854242.689.825316.2412.506341.9915.557332.6915.798389.2916.399453.4018.45
- 画散点图;
- 建立一元线性回归方程,并作显著性检验(取 α=0.05),列出方差分析表;
- 若已知某年社会商品零售总额为 300 亿元,试给出营业税税收总额的概率为 0.95 的预测区间;
- 若已知回归直线过原点,试求回归方程,并在显著性水平 0.05 下作显著性检验.
解
(1)散点图如图 8.2.
\FigureEightTwo
(2)用 MATLAB 统计软件进行回归,得到的回归方程为
y^=−2.26+0.0487x,
其方差分析表如下:
来源回归残差总计平方和203.407.93211.33自由度178均方和203.401.13F 比180p 值0.000
根据以上结果,在显著性水平 α=0.05 下,回归方程是显著的.
(3)按 regression 的 prediction of new observation 功能,当自变量 x 的值取 300 时,可得到 y 的 0.95 预测区间为
[9.688, 14.999].
(4)若拟合不带截距的过原点的回归方程,只要在 options 中不选 Fit intercept 选项,即可得到过原点的回归直线为
y^=0.0417x,
此时检验的 p 值为 0.000,因此在显著性水平 α=0.05 下,过原点的回归方程是显著的.
补充习题及解答
求一回归直线 y=A+Bx,使所有样本点 (x1,y1),(x2,y2),⋯,(xn,yn) 到该直线的垂直距离平方和最小.
解
点 (xi,yi) 到直线 y=A+Bx 的垂直距离的平方为
di2=1+B2(yi−A−Bxi)2.
如今要求 A 与 B,使
φ(A,B)=i=1∑n1+B2(yi−A−Bxi)2
最小,使用微分法,并令其导数为零,可得如下两个方程:
i=1∑n(yi−A^−B^xi)=0,或yˉ−A^−B^xˉ=0,(1)
(1+B^2)i=1∑n(yi−A^−B^xi)xi+B^i=1∑n(yi−A^−B^xi)2=0.(2)
由(1)式可得
A^=yˉ−B^xˉ,
并将其代入(2)式,可得
(1+B^2)i=1∑n[(yi−yˉ)−B^(xi−xˉ)]xi+B^i=1∑n[(yi−yˉ)−B^(xi−xˉ)]2=0.
注意到恒等式
i=1∑n[(yi−yˉ)−B^(xi−xˉ)]xˉ=0,
可将上式化为
(1+B^2)i=1∑n[(yi−yˉ)−B^(xi−xˉ)](xi−xˉ)+B^i=1∑n[(yi−yˉ)−B^(xi−xˉ)]2=0.
使用相同的记号
lxx=∑(xi−xˉ)2,lyy=∑(yi−yˉ)2,lxy=∑(xi−xˉ)(yi−yˉ),
则上式可表示为
(1+B^2)(lxy−B^lxx)+B^(lyy+B^2lxx−2B^lxy)=0.
整理后可得如下 B^ 的二次方程
lxyB^2+(lxx−lyy)B^−lxy=0.
由于判别式
Δ=(lxx−lyy)2+4lxy2≥0,
故此二次方程有实根,
B^=2lxy(lyy−lxx)±(lxx−lyy)2+4lxy2,
这里 B^ 是斜率,根据散点图上的上升趋势或下降趋势选择 B^ 表达式中的正负号.
讨论:这样一来,通过二维样本 (xi,yi)(i=1,2,⋯,n) 建立一元线性回归方程有两个标准:
- 使残差平方和 ∑i=1n(yi−a−bxi)2 达到最小;
- 使回归方程的垂直距离平方和
i=1∑n1+B2(yi−A−Bxi)2
达到最小.
它们导出的回归系数估计也是不同的,但两条回归直线都通过点 (xˉ,yˉ).哪个更好呢?这要看你使用什么标准,实际情况是(1)的理论已很完善,使用的人已很多,效果也很好,而(2)的研究甚少,几乎还无人用在实际中使用.实际研究表明,当相关系数的绝对值 ∣r∣ 接近于 1 时,两个标准下得到的回归系数都很接近,即 a≈A,b≈B.倘若 ∣r∣ 远离 1 时,两个标准下得到的回归系数相差就大了,接下去的两个习题分别说明这个现象.
在用光电比色计检验尿素时,对给定的尿素含量 x(单位:mg/l),消光系数 y 服从正态分布,且方差与 x 无关,观测得如下数据:
尿素含量 x消光系数 y26441386205828510360
试用两个标准分别建立一元回归方程.
解
由这组数据可计算得到
xˉ=6,yˉ=210.4,lxx=40,lyy=54649.2,lxy=1478.
(1)用残差平方和最小的标准,可得两回归系数为
b^=lxxlxy=401478=36.95,a^=yˉ−b^xˉ=210.4−36.95×6=−11.3.
(2)用到回归直线垂直距离平方和最小的标准(见补充习题 12),可得两回归系数为
B^=2×147854649.2−40+(54649.2−40)2+4×14782=36.9751,
A^=210.4−36.9751×6=−11.4506.
比较两个标准下的结果,可见 a^≈A^,b^≈B^,这是因为其相关系数
r=0.99966
很接近 1.
某合金钢的抗拉强度 y 与碳含量 x 有关,现有 92 炉钢样数据,从中算得
xˉ=0.1255,yˉ=45.80,
lxx=0.3018,lyy=2981,lxy=26.70,
试用两个标准分别建立一元回归方程.
解
(1)用残差平方和最小的标准,可得两回归系数为
b^=lxxlxy=0.301826.7=88.47,a^=yˉ−b^xˉ=45.80−88.47×0.1255=34.70.
(2)用到回归直线垂直距离平方和最小的标准(见补充习题 12),可得两回归系数为
B^=2×26.70(2981−0.3018)+(2981−0.3018)2+4×26.72=111.6456,
A^=45.80−111.6456×0.1255=31.7885.
比较两种标准下的结果,可见 a^ 与 A^,b^ 与 B^ 之间相差较大,这是因其相关系数
r=0.8902
与 1 有较大差距.
§8.5 一元非线性回归
- 非线性函数形式 根据二维样本的散点图确定可能的非线性函数形式,部分常见的非线性函数及其线性化变换如下表:
函数名称双曲线函数幂函数指数函数指数函数对数函数S 形曲线函数表达式y1=a+xby=axby=aebxy=aeb/xy=a+blnxy=a+be−x1线性化变换v=y1, u=x1v=lny, u=lnxv=lny, u=xv=lny, u=x1v=y, u=lnxv=y1, u=e−x
- 参数估计 通过适当变换,把非线性函数转化为线性函数形式,然后对未知参数寻求最小二乘估计.譬如由
y1=a+xb
可转化为
v=a+bu,
只要令
v=y1,u=x1
即可.
- 评价标准 常用的曲线回归方程的好坏评价标准有两个:
- 决定系数
R2=1−∑(yi−yˉ)2∑(yi−y^i)2,
愈大愈好;
- 剩余标准差
s=n−2∑(yi−y^i)2,
愈小愈好.
这两个评价标准是一致的,只是从两个侧面作出评价.
习题与解答 8.5
设曲线函数形式为 y=a+blnx,试给出一个变换将之化为一元线性回归的形式.
解
令
u=lnx,v=y,
则原曲线函数化为
v=a+bu,
即为一元线性回归的形式.
设曲线函数形式为 y=a+bx,试给出一个变换将之化为一元线性回归的形式.
解
令
u=x,v=y,
则原函数化为
v=a+bu.
设曲线函数形式为 y−100=ae−x/b (b>0),试给出一个变换将之化为一元线性回归的形式.
解
根据原函数形式,可考虑作如下变换:
u=x,v=ln(y−100).
变换后的线性函数为
v=lna−b1u,
进一步,可将之规范化,令
β0=lna,β1=−1/b,
则最后的回归函数化为
v=β0+β1u.
设曲线函数形式为 y=a+ebx,问能否找到一个变换将之化为一元线性回归的形式?若能,试给出;若不能,说明理由.
解
不能.此处 a 是未知参数,我们不能采用如上题所用的方法,即取
v=ln(y−a),
这样的变换是不通的,因为这样变换后的 v 无法观测.
设曲线函数形式为 y=a+be−x1,问能否找到一个变换将之化为一元线性回归的形式?若能,试给出;若不能,说明理由.
解
能.令
u=e−x,v=y1,
则变换后的函数形式为
v=a+bu.
设曲线函数形式为 y=aeb/x,问能否找到一个变换将之化为一元线性回归的形式?若能,试给出;若不能,说明理由.
解
能.令
u=x1,v=lny,
则变换后的函数形式为
v=lna+bu.
为了检验 X 射线的杀菌作用,用 200 kV 的 X 射线照射杀菌,每次照射 6 min,照射次数为 x,照射后所剩细菌数为 y,下表是一组试验结果:
xy1783262134334431528762517175815491291010311721250134314311528162017161812199207
根据经验知道 y 关于 x 的曲线回归方程形如
y^=aebx,
试给出具体的回归方程,并求其对应的决定系数 R2 和剩余标准差 s.
解
令
v=lny,β0=lna,
则回归方程 y^=aebx 化为
v^=β0+bx.
由数据可算得
lxx=2870−2102/20=665,lxv=751.389−210×87.2250/20=−164.4735,
从而
b^=−164.4735/665=−0.2473,β^0=87.225/20+0.2473×210/20=6.9579.
于是就得到了 lny 关于 x 的线性回归方程
lny^=6.9579−0.2473x,
所以 y 关于 x 的曲线回归方程为
y^=1051.4232e−0.2473x.
R2=1−∑(yi−yˉ)2∑(yi−y^i)2=1−1611149−36552/209210.66=0.9902,
s=189210.66=22.6209.
{\small\bfseries 续表\par}
{\small
\setlength{\tabcolsep}{6pt}
\renewcommand{\arraystretch}{1.08}
| x | y | v=lny | x2 | xv | y2 | y^ | (y−y^)2 |
|---|
| 11 | 72 | 4.2767 | 121 | 47.043 | 5184 | 69.242 | 7.61 |
| 12 | 50 | 3.9120 | 144 | 46.944 | 2500 | 54.071 | 16.57 |
| 13 | 43 | 3.7612 | 169 | 48.896 | 1849 | 42.224 | 0.60 |
| 14 | 31 | 3.4340 | 196 | 48.076 | 961 | 32.973 | 3.89 |
| 15 | 28 | 3.3322 | 225 | 49.983 | 784 | 25.749 | 5.07 |
| 16 | 20 | 2.9957 | 256 | 47.932 | 400 | 20.108 | 0.01 |
| 17 | 16 | 2.7726 | 289 | 47.134 | 256 | 15.702 | 0.09 |
| 18 | 12 | 2.4849 | 324 | 44.728 | 144 | 12.262 | 0.07 |
| 19 | 9 | 2.1972 | 361 | 41.747 | 81 | 9.575 | 0.33 |
| 20 | 7 | 1.9459 | 400 | 38.918 | 49 | 7.478 | 0.23 |
| 和 | 210 | 3655 | 87.2250 | 2870 | 751.389 | 1611149 | {} | 9210.66 |
}
评论
支持 Markdown 和 LaTeX 数学公式。