概率统计Python计算：单因素试验总偏差平方和的分解

单因素试验模型$X$是一个数组的数组： X = { X 1 , X 2 , ⋯   , X s } X=\{X_1,X_2,\cdots,X_s\} X={X1​,X2​,⋯,Xs​}，其中 X i = { X i 1 , X i 2 , ⋯   , X i n i } X_i=\{X_{i1},X_{i2},\cdots,X_{in_i}\} Xi​={Xi1​,Xi2​,⋯,Xini​​}，$i=1,2\cdots ,s$（诸$n_i$未必相等），表示来自对应水平$A_i$的试验指标$N({\mu }_i,{\sigma }^2)$的样本数据。记${\overline{X}}_i=\frac{1}{n_i}\sum _{j=1}^{n_i}{X}_{ij}$，$n=n_1+n_2+\cdots +n_s$，$\overline{X}=\frac{1}{n}\sum _{i=1}^s\sum _{j=1}^{n_i}{X}_{ij}$。$S_T=\sum _{i=1}^s\sum _{j=1}^{n_i}(X_{ij}-\overline{X}{)}^2$，则$S_T$可分解成误差平方和$S_E=\sum _{i=1}^s\sum _{j=1}^{n_i}(X_{ij}-{\overline{X}}_i{)}^2$与效应平方和$S_A=\sum _{i=1}^s\sum _{j=1}^{n_i}({\overline{X}}_i-\overline{X}{)}^2$之和，即
$$S_T=S_E+S_A.$$
$X$，$n$，$s$， { n 1 , n 2 , ⋯   , n s } \{n_1,n_2,\cdots,n_s\} {n1​,n2​,⋯,ns​}， { X ‾ 1 , X ‾ 2 , ⋯   , X ‾ s } \{\overline{X}_1,\overline{X}_2,\cdots,\overline{X}_s\} {X1​,X2​,⋯,Xs​}，$\overline{X}$，$S_T$，$S_E$和$S_A$是单因素试验方差分析的基础数据。下列代码，定义了根据单因素试验模型数据$X$，计算并分解$S_T$的函数。

import numpy as np                                          #导入numpy
def sfeDecompose(X):                                        #X为试验样本数据
    s=X.shape[0]                                            #水平数s
    n=np.array([X[i].size for i in range(s)])               #各水平样本容量
    nt=n.sum()                                              #样本数据总容量
    X_bar=np.array([X[i].mean() for i in range(s)])         #各水平样本均值
    Xt_bar=(X_bar*n).sum()/nt                               #样本数据总均值
    ST=np.sum([((X[i]-Xt_bar)**2).sum() for i in range(s)]) #总平方和ST
    SA=(n*(X_bar**2)).sum()-nt*(Xt_bar**2)                  #效应平方和SA
    SE=ST-SA                                                #误差平方和
    return (n, s, X_bar, Xt_bar, ST, SA, SE)

函数sfeDecompose（sfe是单因素试验single-factor experiment的缩写）的参数X表示单因素试验的样本数据。第3行计算因素个数$s$——它是X的行数。第4行计算对应每个水平的样本容量$n_i$，存于数组n。第5行计算数据总容量$n={\sum }_{i=1}^s{n}_i$为nt。第6行计算对应每个水平的样本均值${\overline{X}}_i$，存于数组X_bar。第7行计算数据总均值$\overline{X}=\frac{1}{n}{\sum }_{i=1}^2{n}_i{\overline{X}}_i$为Xt_bar。第8行计算总平方和$S_T={\sum }_{i=1}^s{\sum }_{j=1}^{n_i}(X_{ij}-\overline{X}{)}^2$为ST。第9行计算效应平方和$S_A={\sum }_{i=1}^s{n}_i{\overline{X}}_i^2-n{\overline{X}}^2$为SA。第10行计算误差平方和$S_E=S_T-S_A$为SE。第11行将所有计算结果作为一个元组返回。
利用函数sfeDecompose算得的数据$n$，$s$，$S_A$和$S_E$，可以计算显著水平$\alpha$下假设
$$H_0:{\mu }_1={\mu }_2=\cdots ={\mu }_s.$$
的检验，下列代码定义计算该假设检验。

def sfeTest(n, s, SA, SE, alpha):
    nt=n.sum()                          #数据总容量
    F=(nt-s)/(s-1)*SA/SE                #检验统计量值
    accept=ftestR(F, s-1, nt-s, alpha)  #F分布的分位点
    return accept

函数sfeTest的参数n，s，SA，SE和alpha中除了alpha表示的是显著水平$\alpha$外，其余的均与函数sfeDecompose函数中所算得的同名变量意义相同，此不赘述。第2行计算数据总容量$n={\sum }_{i=1}^s{n}_i$，第3行计算检验统计量值$\frac{S_A/(s-1)}{S_E/(n-s)}$为F。第4行调用ftestR函数用p值法计算假设$H_0:{\mu }_1={\mu }_2=\cdots ={\mu }_s=\mu$的检验，结果为accept。第6行的返回值accept若为True，则接受假设$H_0$，否则拒绝。
例1设有三台机器，用来生产规格相同铝合金薄板。取样，测量薄板的厚度（精确至千分之一厘米）。得到如下结果：
$$\begin{gathered} \text{机器I：}0.236,0.238,0.248,0.245,0.243 \\ \text{机器II：}0.257,0.253,0.255,0.254,0.261 \\ \text{机器III：}0.258,0.264,0.259,0.267,0.262 \end{gathered}$$
对于第$i$台机器，所生产的薄板厚度为随机变量$X_i$~$N({\mu }_i,{\sigma }^2)$。目标是，利用试验数据检验假设
$$H_0:{\mu }_1={\mu }_2={\mu }_3(H_1:{\mu }_1,{\mu }_2,{\mu }_3\text{不全相等}).$$
即判断不同的机器是否显著影响所生产的铝合金薄板的厚度。
解： 本例中，试验指标为薄板厚度。试验中诸如材料、操作人员等条件均视为不变的因素。唯一可变因素是所用的机器。不同的三台机器就是该因素的三个不同水平。因此，这是一个单因素试验。下列代码完成本例计算。

import numpy as np                                          #导入numpy
alpha=0.05                                                  #显著水平
X=np.array([np.array([0.236, 0.238, 0.248, 0.245, 0.243]),  #试验数据
            np.array([0.257, 0.253, 0.255, 0.254, 0.261]),
            np.array([0.258, 0.264, 0.259, 0.267, 0.262])])
(n, s, X_bar, Xt_bar, ST, SA, SE)=sfeDecompose(X)           #平方和分解
accept=sfeTest(n, s, SA, SE, alpha)							#假设检验
print('显著水平%.2f下H0为%s'%(alpha,accept))

第2~5行根据题面设置试验数据X和显著水平alpha。第6行调用函数sfeDecompose(X)计算平方和分解。第7行调用函数sfeTest(n, s, SA, SE, alpha)计算假设$H_0:{\mu }_1={\mu }_2={\mu }_3$的检验。运行程序，输出：

显著水平0.05下H0为False

表示拒绝假设$H_0$，即认为不同的机器显著影响所生产的铝合金薄板的厚度。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
代码诚可贵，原理价更高。若为AI学，读正版书好。
返回《导引》