【数学建模算法】(31)方差分析（上）

2019-08-27 本文已影响0人热爱学习的高老板

我们已经作过两个总体均值的假设检验，如两台机床生产的零件尺寸是否相等，病人和正常人的某个生理指标是否一样。如果把这类问题推广一下，要检验两个以上总体的均值彼此是否相等，仍然用以前介绍的方法是很难做到的。而你在实际生产和生活中可以举出许多这样的问题：从用几种不同工艺制成的灯泡中，各抽取了若干个测量其寿命，要推断这几种工艺制成的灯泡寿命是否有显著差异；用几种化肥和几个小麦品种在若干块试验田里种植小麦，要推断不同的化肥和品种对产量有无显著影响。

可以看到，为了使生产过程稳定，达到优质、高产，需要对影响产品质量的因素进行分析，找出有显著影响的那些因素，除了从机理方面进行研究外，常常要作许多试验，对结果作分析、比较，寻求规律。用数理统计分析试验结果、鉴别各因素对结果影响程度的方法称为方差分析（Analysis Of Variance），记作 ANOVA。

人们关心的试验结果称为指标，试验中需要考察、可以控制的条件称为因素或因子，因素所处的状态称为水平。上面提到的灯泡寿命问题是单因素试验，小麦产量问题是双因素试验。处理这些试验结果的统计方法就称为单因素方差分析和双因素方差分析。

1.单因素方差分析

只考虑一个因素 $A$ 对所关心的指标的影响， $A$ 取几个水平，在每个水平上作若干个试验，试验过程中除 $A$ 外其它影响指标的因素都保持不变（只有随机因素存在），我们的任务是从试验结果推断，因素 $A$ 对指标有无显著影响，即当 $A$ 取不同水平时指标有无显著差别。

$A$ 取某个水平下的指标视为随机变量，判断 $A$ 取不同水平时指标有无显著差别，相当于检验若干总体的均值是否相等。

1.1.数学模型

设 $A$ 取 $r$ 个水平 $A_{1}, A_{2}, \cdots, A_{r}$ ，在水平 $A_{i}$ 下总体 $x_{i}$ 服从正态分布 $N\left(\mu_{i}, \sigma^{2}\right)$ ， $i=1, \cdots, r$ 这里 $\mu_{i}, \sigma^{2}$ 未知， $\mu_{i}$ 可以互不相同，但假定 $x_{i}$ 有相同的方差。又设在每个水平 $A_{i}$ 下作了 $n_{i}$ 次独立试验，即从中抽取容量为 $n_{i}$ 的样本，记作 $x_{i j}, j=1, \cdots, n_{i}, x_{i j}$ 服从 $N\left(\mu_{i}, \sigma^{2}\right), \quad i=1, \cdots, r, j=1, \cdots, n_{i}$ 且相互独立。将这些数据列成表单因素试验数据表）的形式。

单因素试验数据表

$A_{1}$	$x_{11}$	$x_{12}$	...	$x_{1n_{1}}$
$A_{2}$	$x_{21}$	$x_{22}$	...	$x_{2n_{2}}$
...	...	...	...	...
$A_{r}$	$x_{r1}$	$x_{r2}$	...	$x_{rn_{r}}$

将第 $i$ 行称为第 $i$ 组数据。判断 $A$ 的 $r$ 个水平对指标有无显著影响，相当于要作以下的假设检验:
$H_{0} : \mu_{1}=\mu_{2}=\cdots=\mu_{r} ; \quad H_{1} : \mu_{1}, \mu_{2}, \cdots, \mu_{r}$ 不全相等。

由于 $x_{i j}$ 的取值既受不同水平 $A_{i}$ 的影响，又受 $A_{i}$ 固定下随机因素的影响，所以将它分解为：
$x_{i j}=\mu_{i}+\varepsilon_{i j}, i=1, \cdots, r, j=1, \cdots, n_{i}（1）$
其中 $\varepsilon_{i j} \sim N\left(0, \sigma^{2}\right)$ ，且相互独立。

记：
$\mu=\frac{1}{n} \sum_{i=1}^{r} n_{i} \mu_{i}, \quad n=\sum_{i=1}^{r} n_{i}, \quad \alpha_{i}=\mu_{i}-\mu, \quad i=1, \cdots, r（2）$

$\mu$ 是总均值， $\alpha_{i}$ 是水平 $A_{i}$ 对指标的效应。由（1）、（2）模型可表为：
$\left\{\begin{array}{l}{x_{i j}=\mu+\alpha_{i}+\varepsilon_{i j}} \\ {\sum_{i=1}^{r} \alpha_{i}=0} \\ {\varepsilon_{i j} \sim N\left(0, \sigma^{2}\right), i=1, \cdots, r, j=1, \cdots, n_{i}}\end{array}\right.（3）$

原假设为（以后略去备选假设）：
$H_{0} : \alpha_{1}=\alpha_{2}=\cdots=\alpha_{r}=0（4）$

1.2.统计分析

记：
$\overline{x}_{i \bullet}=\frac{1}{n_{i}} \sum_{j=1}^{n_{i}} x_{i j}, \overline{x}=\frac{1}{n} \sum_{i=1}^{r} \sum_{j=1}^{n_{i}} x_{i j}（5）$

$\overline{x}_{i}$ 是第 $i$ 组数据的组平均值， $\overline{x}$ 是总平均值。考察全体数据对 $\overline{x}$ 的偏差平方和：
$S_{T}=\sum_{i=1}^{r} \sum_{j=1}^{n_{i}}\left(x_{i j}-\overline{x}\right)^{2}（6）$

经分解可得：
$S_{T}=\sum_{i=1}^{r} n_{i}\left(\overline{x}_{i \cdot}-\overline{x}\right)^{2}+\sum_{i=1}^{r} \sum_{j=1}^{n_{i}}\left(x_{i j}-\overline{x}_{i \bullet}\right)^{2}$

记：
$S_{A}=\sum_{i=1}^{r} n_{i}\left(\overline{x}_{i \cdot}-\overline{x}\right)^{2}（7）$
$S_{E}=\sum_{i=1}^{r} \sum_{j=1}^{n_{i}}\left(x_{i j}-\overline{x}_{i .}\right)^{2}（8）$

则：
$S_{T}=S_{A}+S_{E}（9）$

$S_{A}$ 是各组均值对总方差的偏差平方和，称为组间平方和； $S_{E}$ 是各组内的数据对均值偏差平方和的总和。 $S_{A}$ 反映 $A$ 不同水平间的差异， $S_{E}$ 则表示在同一水平下随机误差的
大小。

注意到 $\sum_{j=1}^{n_{i}}\left(x_{i j}-\overline{x}_{i \bullet}\right)^{2}$ 是总体 $N\left(\mu_{i}, \sigma^{2}\right)$ 的样本方差的 $n_{i}-1$ 倍，于是有：
$\sum_{j=1}^{n_{i}}\left(x_{i j}-\overline{x}_{i .}\right)^{2} / \sigma^{2} \sim \chi^{2}\left(n_{i}-1\right)$
由 $\chi^{2}$ 分布的可加性知：
$S_{E} / \sigma^{2} \sim \chi^{2}\left(\sum_{i=1}^{r}\left(n_{i}-1\right)\right)$
即：
$S_{E} / \sigma^{2} \sim \chi^{2}(n-r)$
且有：
$E S_{E}=(n-r) \sigma^{2}（10）$
对 $S_{A}$ 作进一步分析可得：
$E S_{A}=(r-1) \sigma^{2}+\sum_{i=1}^{r} n_{i} \alpha_{i}^{2}（11）$
若 $H_{0}$ 成立，则：
$E S_{A}=(r-1) \sigma^{2}（12）$

可知若 $H_{0}$ 成立， $S_{A}$ 只反映随机波动，而若 $H_{0}$ 不成立，那它就还反映了 $A$ 的不同水平的效应 $\alpha_{i}$ 。单从数值上看，当 $H_{0}$ 成立时，由（10）、（12）对于一次试验应有：
$\frac{S_{A} /(r-1)}{S_{E} /(n-r)} \approx 1$

而当 $H_{0}$ 不成立时这个比值将远大于 1。当 $H_{0}$ 成立时，该比值服从自由度 $n_{1}=r-1$ ， $n_{2}=(n-r)$ 的 $F$ 分布，即：
$F=\frac{S_{A} /(r-1)}{S_{E} /(n-r)} \sim F(r-1, n-r)（13）$

为检验 $H_{0}$ ，给定显著性水平 $\alpha$ ，记 $F$ 分布的 $1-\alpha$ 分位数为 $F_{1-\alpha}(r-1,(n-r))$ ，检验规则为：
$F<F_{1-\alpha}(r-1,(n-r))$ 时接受 $H_{0}$ ，否则拒绝。

以上对 $S_{A}, S_{E}, S_{T}$ 的分析相当于对组间、组内等方差的分析，所以这种假设检验方法称方差分析。

1.3.方差分析表

将试验数据按上述分析、计算的结果排成表 2 的形式，称为单因素方差分析表（Matlab 中给出的方差分析表）。

单因素方差分析表

方差来源	平方和	自由度	均方	$1-p_{r}$ 分位数	概率
因素 $A$	$S_{A}$	$r-1$	$\overline{S}_{A}=\frac{S_{A}}{r-1}$	$F_{1-p_{r}}(r-1, n-r)$	$p_{r}$
误差	$S_{E}$	$n-r$	$\overline{S}_{E}=\frac{S_{E}}{n-r}$
总和	$S_{T}$	$n-1$

最后一列给出大于 $F$ 值的概率 $p_{r}, \quad F_{1-p_{r}}<F_{1-\alpha}(r-1,(n-r))$ 相当于 $p_{r}>\alpha$ 。

方差分析一般用的显著性水平是：取 $\alpha=0.01$ ，拒绝 $H_{0}$ ，称因素 $A$ 的影响（或 $A$ 各水平的差异）非常显著；取 $\alpha=0.01$ ，不拒绝 $H_{0}$ ，但取 $\alpha=0.05$ ，拒绝 $H_{0}$ ，称因素 $A$ 的影响显著；取 $\alpha=0.05$ ，不拒绝 $H_{0}$ ，称因素 $A$ 无显著影响。

1.4.Matlab实现

Matlab 统计工具箱中单因素方差分析的命令是 anoval。

若各组数据个数相等，称为均衡数据。若各组数据个数不等，称非均衡数据。

1.4.1.均衡数据

处理均衡数据的用法为：

p=anoval(x)

返回值 $p$ 是一个概率，当 $p>\alpha$ 时接受 $H_{0}$ ， $x$ 为 $m \times r$ 的数据矩阵， $x$ 的每一列是一个水平的数据（这里各个水平上的样本容量 $n_{i}=m$ ）。另外，还输出一个方差表和一个
Box 图。

例1 为考察 5 名工人的劳动生产率是否相同，记录了每人 4 天的产量，并算出其平均值，如下表。你能从这些数据推断出他们的生产率有无显著差别吗？

	$A_{1}$	$A_{2}$	$A_{3}$	$A_{4}$	$A_{5}$
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

编写程序如下：

clear,clc
x=[256 254 250 248 236
242 330 277 280 252
280 290 230 305 220
298 295 302 289 252];
p=anova1(x)

Matlab会生成方差分析表：

方差分析表

求得 $p=0.1109>\alpha=0.05$ ，故接受 $H_{0}$ ，即 5 名工人的生产率没有显著差异。方差表对应于上面的单因素方差分析表的 1 ~ 4列， $F=2.26$ 是 $F(4,15)$ 分布的 $1-p$ 分位数，可以验证：

fcdf(2.262,4,15)=0.8891=1-p

同时程序会生成箱式图。

箱式图

注：接受 $H_{0}$ ，是将 5 名工人的生产率作为一个整体进行假设检验的结果，并不表明取其中 2 个工人的生产率作两总体的均值检验时，也一定接受均值相等的假设。实际上，读者可以用 ttest2 对本题作 $H_{0} : \mu_{2}=\mu_{5}$ 的检验，看看会得到什么结果。

1.4.2.非均衡数据

处理非均匀数据的用法为：

p=anova1(x,group)

x为向量，从第1组到第r组数据依次排列；group为与x同长度的向量，标志x中数据的组别（在与x第 $i$ 组数据对应的位置输入整数 $i(i=1,2, \cdots, r)$ ）。

例2 用4种工艺生产灯泡，从各种工艺制成的灯泡中各抽出了若干个测量其寿命，结果如下表，试推断这几种工艺制成的灯泡寿命是否有显著差异。

	$A_{1}$	$A_{2}$	$A_{3}$	$A_{4}$
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<\mathrm{p}=0.0331<0.05$ ，所以几种工艺制成的灯泡寿命有显著差异。

1.5.多重比较

在灯泡寿命问题中，为了确定哪几种工艺制成的灯泡寿命有显著差异，我们先算出各组数据的均值：

工艺	$A_{1}$	$A_{2}$	$A_{3}$	$A_{4}$
均值	1708	1635	1540	1585

虽然 $A_{1}$ 均值最大，但要判断它与其它几种有显著差异，还需做多重比较。一般多重比较要对所有 $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)]