概率统计Python计算：单因素试验参数的区间估计

对单因素试验模型 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)$的样本数据。调用平方和分解函数sfeDecompose(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$。利用其中的$n$，$s$，$S_A$和$S_E$，调用假设检验函数sfeTest可以算得假设$H_0:{\mu }_1={\mu }_2=\cdots ={\mu }_s$在显著水平$\alpha$下的检验。无论检验结果是接受还是拒绝假设$H_0$，我们都用数轴量$\frac{S_E}{{\sigma }^2}$~${\chi }^2(n-s)$可算得${\sigma }^2$的置信水平为$1-\alpha$的置信区间。若接受$H_0$，则可利用枢轴量$\frac{\overline{X}-\mu }{\sqrt{\frac{S_T}{(n-1)n}}}$~$t(n-1)$算得$\mu$的置信区间。若拒绝假设$H_0$，则可利用枢轴量$\frac{{\overline{X}}_i-{\overline{X}}_j-({\mu }_i-{\mu }_j)}{\sqrt{\frac{S_E}{n-s}\left(\frac{1}{n_i}+\frac{1}{n_j}\right)}}$~$t(n-s)$算得${\mu }_i-{\mu }_j$的置信区间，$1\le i<j\le s$。下列代码按上述思想定义计算单因素试验参数区间估计的函数。

import numpy as np
def sfeEstimat(accept, n, s, X_bar, Xt_bar, ST, SE, alpha):
    ans=[]                                              #初始化返回值
    nt=n.sum()                                          #数据总容量
    (a, b)=sigma2Bounds(SE, nt-s, 1-alpha)             #sigma^2的置信区间
    ans.append((a, b))
    if accept:                                          #若H0为真
        d=np.sqrt(ST/(nt-1)/nt)
        (a, b)=muBounds(Xt_bar, d, 1-alpha, nt-1)       #计算mu的置信区间
        ans.append((a, b))
    else:                                               #若H0为假
        for i in range(s):                              #对每个i
            for j in range(i+1, s):                     #对每个j>i
                mean=X_bar[i]-X_bar[j]                  #差mui-muj
                S_E=SE/(nt-s)                           #sigma^2估计值
                d=np.sqrt(S_E*(1/n[i]+1/n[j]))          #置信区间增量因子
                (a, b)=muBounds(mean, d, 1-alpha, nt-s) #置信区间
                ans.append((a, b))                      #置信区间
    return np.array(ans)

函数sfeEstimat的参数accept表示是否接受假设$H_0$，除此之外的其他参数均与由调用函数sfeDecompose算得的同名变量的意义相同，此不赘述。第3行将返回值ans初始化为空的list。第4行计算数据总容量nt。第5行调用计算正态总体$N(\mu ,{\sigma }^2)$的参数${\sigma }^2$的函数sigma2Bounds，计算参数${\sigma }^2$的置信水平为$1-\alpha$的置信区间。第6行将区间数据(a,b)加入ans。第7~18行的if-else语句分别就假设$H_0$为真或假计算$\mu$的置信区间或诸${\mu }_i-{\mu }_j$，$1\le i<j\le s$的置信区间。其中第8~10行调用计算正态总体参数$\mu$的置信区间的函数muBounds，计算参数$\mu$的置信水平为$1-\alpha$的置信区间(a,b)。第12~18行的双重for语句，计算诸${\mu }_i-{\mu }_j$的置信水平为$1-\alpha$的置信区间。第14行计算差${\overline{X}}_i-{\overline{X}}_j$为mean，第15行计算${\sigma }^2$的无偏估计值$\frac{S_E}{n-s}$为S_E，第16行计算置信区间增量因子$\sqrt{\frac{S_E}{n-s}\left(\frac{1}{n_i}+\frac{1}{n_j}\right)}$为d。第17行调用函数muBounds计算${\mu }_i-{\mu }_j$的置信区间。
例1制造某型号计算器要用到某种类型的电路板。电路板由四家工厂提供，分别随机选取使用来自各厂家电路板的计算器，其响应时间（以毫秒计）列表如下：
$$\begin{gathered} \text{厂家I：}19,22,20,18,15 \\ \text{厂家II：}20,21,33,27,40 \\ \text{厂家III：}16,15,18,26,17 \\ \text{厂家IV：}18,22,19 \end{gathered}$$
判断不同厂家的电路是否显著影响计算器的计算响应时间。
解： 本例中，试验指标为计算响应时间。可变因素为使用的不同厂家的电路板，该因素有4个水平。这也是一个单因素试验，设用第$i$个厂家生产的电路计算响应时间为随机变量$X_i$~$N({\mu }_i,{\sigma }^2)$，$i=1,2,3,4$。为判断不同的厂家的电路是否显著影响计算器的计算响应时间，利用试验数据检验假设：
$$H_0:{\mu }_1={\mu }_2={\mu }_3={\mu }_4(H_1:{\mu }_1,{\mu }_2,{\mu }_3,{\mu }_4\text{不全相等}).$$
下列代码完成本例计算。

import numpy as np											#导入numpy
X=np.array([np.array([19, 22, 20, 18, 15]),					#试验数据
           np.array([20, 21, 33, 27, 40]),
           np.array([16, 15, 18, 26, 17]),
           np.array([18, 22, 19])])
alpha=0.05													#显著水平
(n, s, X_bar, Xt_bar, ST, SA, SE)=sfeDecompose(X)			#方差分解
accept=sfeTest(n, s, SA, SE, alpha)							#假设检验
ans=sfeEstimat(accept, n, s, X_bar, Xt_bar, ST, SE, alpha)	#参数估计
print(‘H0 is %s’%accept)
for a,b in ans:
   print('(%.3f, %.3f)'%(a,b))

第2~6行按题面设置数据。第7行调用平方和分解函数sfeDecompose(X)，算得数据项$n$，$s$，$S_A$和$S_E$等，第8行调用假设检验函数sfeTest传递这些数据和显著水平$\alpha$，计算对假设$H_0:{\mu }_1={\mu }_2={\mu }_3={\mu }_4$的检验，结果存储于accept。第9行调用参数区间估计计算函数sfeEstimat根据假设检验的结果accept，利用数据$n$，$s$，$S_T$和$S_E$和显著水平$\alpha$计算各参数的置信区间。运行程序，输出

H0 is False
(15.141, 70.259)
(-16.609, -2.191)
(-6.809, 7.609)
(-9.191, 7.458)
(2.591, 17.009)
(0.209, 16.858)
(-9.591, 7.058)

其中，第1行表示拒绝假设$H_0$，第2行表示${\sigma }^2$的置信区间，由于拒绝假设$H_0:{\mu }_1={\mu }_2={\mu }_3={\mu }_4$，故第3~8行显示${\mu }_1-{\mu }_2$，${\mu }_1-{\mu }_3$，${\mu }_1-{\mu }_4$，${\mu }_2-{\mu }_3$，${\mu }_2-{\mu }_4$和${\mu }_3-{\mu }_4$在0.95的置信水平下的置信区间。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
代码诚可贵，原理价更高。若为AI学，读正版书好。
返回《导引》