R函数实现正态总体均值、方差的区间估计及假设检验的设计

张应应，魏 毅

(重庆大学 数学与统计学院，重庆 401331)

0 引言

正态总体均值、方差的区间估计与假设检验是数理统计中的经典内容。数理统计的教材[1～6]一般都会讲到。针对摘要中提到的R软件[7]内置程序t.test()、var.test()函数的缺陷，参考文献[1]中为实现单个、两个正态总体均值、方差的区间估计、假设检验时自编了12个函数interval_estimate1()、 interval_estimate2()、 interval_estimate4()、 interval_estimate5()、interval_var1()、interval_var2()、interval_var3()、interval_var4()、mean.test1()、mean.test2()、var.test1()、var.test2()，这些函数可以实现R软件的内置函数t.test()、var.test()的全部功能，并能有效弥补t.test()和var.test()的缺陷。但是要记住并灵活掌握这么多函数是一件非常麻烦的事情，并且我们也常常同时需要区间估计和假设检验的结果。本文在文献[1]的启发下创造了一个R函数Interval-Estimate_TestOfHypothesis()，它可以实现 t.test()和 var.test()的所有功能及它们不能完成的上述功能，只用一个R函数便能实现单个、两个正态总体均值、方差的所有区间估计及假设检验。

1 程序设计

1.1 P值计算

在软件计算中，通常计算随机变量X大于或小于某个指定值的概率，称为P值。

以正态分布为例，在给定z值后，计算原理[1]如下：

图1 正态总体双边检验(H1: μ≠μ0)

图2 正态总体单边检验(H1: μ＞μ0)

考虑到假设检验中多处需要计算P值，编写R函数p_value.R来实现不同分布P值的计算。相应的，当给定概率值α时，可计算出对应的上分位数q值，编写R函数q_value.R来实现不同分布分位数的计算。

1.2 值的区间估计和假设检验

1.2.1 单个正态总体

正态总体 X～N(μ,σ2)，X1,X2,…,Xn为来自总体 X的一个样本，1-α为置信度，为样本均值，S2为样本方差。在作总体X均值的区间估计时，需分别讨论方差σ2已知和未知两种情形；作假设检验时，在单边、双边检验情况同样需要区分方差σ2已知和未知[1]。程序设计见下图3。

图3

编写R函数命名为single_mean.R

1.2.2 两个正态总体

图4

编写R函数命名为double_mean.R

1.3 方差的区间估计和假设检验

1.3.1 单个正态总体

作方差的区间估计时，分总体X的均值μ已知和未知两种情形讨论，作假设检验时，又须在μ已知和未知的情形下分单边、双边检验。程序设计类似于单个正态总体均值的区间估计和假设检验，故省略。编写R函数命名为single_var.R

1.3.2 两个正态总体

2 R函数IntervalEstimate_TestOfHypothesis()

在以上程序基础上，还需编写主函数来调用子程序以实现不同的功能，主函数使用格式如下：

IntervalEstimate_TestOfHypothesis(x,y=NULL,test=c(“mean”,“variance”),mu=c(Inf,Inf),sigma=c(-1,-1),var.equal=FALSE,ratio=1,side=c(“two.sided”,“less”,“greater”),alpha=0.05)

其中x，y是由样本数据构成的向量；

y默认值为NULL，即默认为对单总体进行操作；

test为检验的类型，默认值为“mean”，代表作均值的区间估计和假设检验，test=“variance”或”var”代表作方差的区间估计和假设检验；

mu为总体的均值向量，在方差的区间估计和假设检验以及单总体均值的假设检验中会用到，默认值为Inf(即未知)；

sigma为总体的标准差向量，默认值均为-1(即未知)，当程序用于作单总体方差假设检验时，默认为检验σ2=1；

var.equal判断两总体方差是否相等，默认为FALSE，此参数在两总体均值检验中用到；

ratio为两总体方差比率，默认为1，此参数在两总体方差检验中用到；

side判断求置信区间和作假设检验的类型，默认值为”two.sided”，即作双边检验并求双侧置信区间；side=”less”或“l”，表示求置信区间上限并作单边检验(H1: μ＜μ0)；side=”greater”或“g”，表示求置信区间下限并作单边检验(H1: μ＞μ0)；

alpha为一个取值为[0,1]的实数，默认为0.05，1-alpha为置信度。

由于本程序是对正态总体进行操作，因此，在使用本函数前须先确认样本数据服从正态分布，为此，编写R程序testNormal_plot.R来对样本做正态性检验，其调用格式为：

表1 正态总体区间估计及假设检验函数IntervalEstimate_TestOfHypothesis()的使用方法表

testNormal_plot(x,alpha)

其中x为待测样本向量，alpha意义同上。

表1中，side=”two.sided”（程序默认值，可不必输入）表示求双侧置信区间并作双边检验，alpha=0.05（默认值，可不必输入）表示显著性水平为0.05，test=”var”（效果与test=”variance”相同）表示对输入变量作方差的区间估计和假设检验，ratio=2（ratio=默认值为1），代表作两总体方差检验的原假设为

接下来，将举例来测试函数IntervalEstimate_TestOf-Hypothesis()（以下简称待测函数），以验证其正确性。

testNormal_plot()用于测试数据x是否来自正态总体。

程序结果解释set.seed(1)x=rnorm(100)testNormal_plot(x,alpha=0.05)Shapiro-Wilk normality test data:x W=0.9956,p-value=0.9876 The data comes from a normal distribution.p-value=0.9876 ＞ 0.05=alpha，接受假设H0:数据x来自正态总体。set.seed(1)x=rt(100,5)testNormal_plot(x,alpha=0.05)Shapiro-Wilk normality test data:x W=0.9232,p-value=2.088e-05 The data does not come from a normal distribution!p-value=2.088e-05＜0.05=alpha，拒绝假设H0:数据x来自正态总体。

正态随机数检验图像

产生样本x，此x将用于例1～例4的测试

程序结果解释set.seed(1)x=rnorm(10,mean=1,sd=0.2);x[1]0.8747092 1.0367287 0.8328743 1.3190562 1.0659016[6]0.8359063 1.0974858 1.1476649 1.1151563 0.9389223产生一个来自正态总体的均值为1，标准差为0.2的容量为10的样本，为便于重复结果，我们使用set.seed(1)。

例1:单总，均值，sigma2(即 σ2，下同)已知

程序结果解释testNormal_plot(x)Shapiro-Wilk normality test data:x W=0.9383,p-value=0.534 The data comes from a normal distribution.p-value=0.534＞0.05=alpha，接受假设H0:数据x来自正态总体。

程序结果解释IntervalEstimate_TestOfHypothesis(x,sigma=0.2)One sample mean alternative hypothesis:true mean is not equal to 0$statistic_df_Pvalue z df p-value 1 16.22945 10 0$OneMinusAlpha_confidence_interval low up 1 0.9024815 1.1504$mean_of_x[1]1.026441待测函数与文献中函数结果一致，95%置信区间为[0.9024815,1.1504]，p-value=0＜0.05=alpha，拒绝原假设，认为均值不为0。interval_estimate4(x,0.2)mean.test1(x,sigma=0.2)＞interval_estimate4(x,0.2)mean df a b 1 1.026441 10 0.9024815 1.1504＞mean.test1(x,sigma=0.2)mean df Z P_value 1 1.026441 10 16.22945 0

例2:单总，均值，sigma2未知

程序结果解释IntervalEstimate_TestOfHypothesis(x,mu=0.95)One sample mean alternative hypothesis:true mean is not equal to 0.95$statistic_df_Pvalue t df p-value 1 1.548364 9 0.1559421$OneMinusAlpha_confidence_interval low up 1 0.914761 1.13812$mean_of_x[1]1.026441待测函数与t.test()结果相同，95%置信区间为[0.914761,1.13812]，p-value=0.1559421＞0.05=alpha，接受原假设，认为均值为0.95。t.test(x,mu=0.95)One Sample t-test data:x t=1.5484,df=9,p-value=0.1559 alternative hypothesis:true mean is not equal to 0.95 95 percent confidence interval:0.914761 1.138120 sample estimates:mean of x 1.026441

例3:单总，方差，mu已知，sigma不输入时，默认检测σ2=1

程序结果解释IntervalEstimate_TestOfHypothesis(x,test="variance",mu=1)one sample variance alternative hypothesis:true variance is not equal to 1$statistic_df_Pvalue chi2 df p-value 1 0.2263442 10 2.816068e-07$OneMinusAlpha_percent_confidence_interval low up 1 0.01105025 0.06970931待测函数与参考文献中结果相同，95%置信区间为[0.01105025,0.06970931]，p-value=2.816068e-07＜0.05=alpha，拒绝原假设，认为方差不为1。interval_var3(x,mu=1)var.test1(x,mu=1)＞interval_var3(x,mu=1)var df a b 1 0.0226344 10 0.011050 0.069709＞var.test1(x,mu=1)var df chisq2 P_value 1 0.022634 10 0.2263442 2.816068e-07

例4:单总，方差，mu未知

程序结果解释IntervalEstimate_TestOfHypothesis(x,test="variance",sigma=0.18)one sample variance alternative hypothesis:true variance is not equal to 0.0324$statistic_df_Pvalue chi2 df p-value 1 6.77016 9 0.6779301$OneMinusAlpha_percent_confidence_interval low up 1 0.01153109 0.08123021待测函数与文献中函数结果相同，95%置信区间为[0.01153109,0.08123021]，p-value=0.6779301＞0.05=alpha，接受原假设，认为标准差为0.18。interval_var3(x)var.test1(x,sigma2=0.18^2)＞interval_var3(x)var df a b 1 0.024372 9 0.011531 0.08123021＞var.test1(x,sigma2=0.18^2)var df chisq2 P_value 1 0.02437258 9 6.77016 0.6779301

产生正态样本，此x,y1,y2将用于例5～例9的测试

set.seed(1);x=rnorm(10,mean=1,sd=0.2);x

set.seed(2);y1=rnorm(15,mean=2,sd=0.3);y1

set.seed(2);y2=rnorm(15,mean=2,sd=0.2);y2

testNormal_plot(y1)

testNormal_plot(y2)

以下的数据正态性检验与之前相似，检验结果不再列出。

例5：双总，均值，方差已知

程序结果解释IntervalEstimate_TestOfHypothesis(x,y1,sigma=c(0.2,0.3))double sample mean alternative hypothesis:true difference in means is not equal to 0$statistic_df_Pvalue z df p-value 1-10.50775 25 7.956809e-26$OneMinusAlpha_confidence_interval low up 1-1.246772-0.8547787$mean_of_x_and_y xb yb 1 1.026441 2.077216待测函数与文献中函数结果相同，95%置信区间为[1.026441,2.077216]，p-value=7.956809e-26＜0.05=alpha，拒绝原假设，认为两总体均值差不为0。interval_estimate5(x,y1,sigma=c(0.2,0.3))mean.test2(x,y1,sigma=c(0.2,0.3))interval_estimate5(x,y1,sigma=c(0.2,0.3))mean df a b 1-1.050775 25-1.246772-0.8547787＞mean.test2(x,y1,sigma=c(0.2,0.3))mean df Z P_value 1-1.050775 25-10.50775 7.956809e-26

例6：双总，均值，方差未知但相等

程序结果解释IntervalEstimate_TestOfHypothesis(x,y2,var.equal=T)double sample mean alternative hypothesis:true difference in means is not equal to 0$statistic_df_Pvalue t df p-value 1-13.73411 23 1.429158e-12$OneMinusAlpha_confidence_interval low up 1-1.17943-0.8706436$mean_of_x_and_y xb yb 1 1.026441 2.051477待测函数与t.test()结果相同，95%置信区间为[-1.17943,-0.8706436]，p-value=1.429158e-12＜0.05=alpha，拒绝原假设，认为两总体均值差不为0。t.test(x,y2,var.equal=T)Two Sample t-test data:x and y2 t=-13.7341,df=23,p-value=1.429e-12 alternative hypothesis:true difference in means is not equal to 0 95 percent confidence interval:-1.1794296-0.8706436 sample estimates:mean of x mean of y 1.026441 2.051477

例7：双总，均值，方差未知且不等

程序结果解释IntervalEstimate_TestOfHypothesis(x,y1,var.equal=F,side="less")double sample mean alternative hypothesis:true difference in means is less than 0$statistic_df_Pvalue t df p-value 1-11.51773 22.10023 4.112604e-11$OneMinusAlpha_confidence_interval low up 1-Inf-0.8941493$mean_of_x_and_y xb yb 1 1.026441 2.077216待测函数与t.test()结果相同，95%置信区间为[-Inf,-0.8941493]，p-value=4.112604e-11＜0.05=alpha，拒绝原假设，认为两总体均值差小于0。t.test(x,y1,var.equal=F,al='less')Welch Two Sample t-test data:x and y1 t=-11.5177,df=22.1,p-value=4.113e-11 alternative hypothesis:true difference in means is less than 0 95 percent confidence interval:-Inf-0.8941493 sample estimates:mean of x mean of y 1.026441 2.077216

说明：例1～例9基本覆盖了单个、两个正态总体均值、方差的区间估计和假设检验的所有情况，IntervalEstimate_TestOfHypothesis()运行的结果与R内置函数t.test()，var.test()所得的结果完全一致，对于t.test()，var.test()所不能解决的例子(包括例1，例3，例4，例5，例8)，笔者使用参考文献中的函数interval_estimate4(),interval_estimate5(),interval_var3(),interval_var4(),mean.test1(),mean.test2(),var.test1(),var.test2()进行了检测，结果也完全吻合。

例8：双总，方差，均值已知

程序结果解释IntervalEstimate_TestOfHypothesis(x,y1,test="var",mu=c(1,2))double sample variance alternative hypothesis:true ratio of the variances is not equal to 1$statistic_df_Pvalue f df1 df2 p-value 1 0.2561489 10 15 0.03500855$OneMinusAlpha_percent_confidence_interval low up 1 0.0837034 0.9020726$ratio_of_variances[1]0.2561489待测函数与文献中函数结果相同，95%置信区间为[0.0837034,0.9020726]，p-value=0.03500855＜0.05=alpha，拒绝原假设，认为方差比率不为1。interval_var4(x,y1,mu=c(1,2))var.test2(x,y1,mu=c(1,2))＞interval_var4(x,y1,mu=c(1,2))rate df1 df2 a b 1 0.2561489 10 15 0.0837034 0.9020726＞var.test2(x,y1,mu=c(1,2))Rate df1 df2 F P_value 1 0.2561489 10 15 0.2561489 0.03500855

例9：双总，方差，均值未知

程序结果解释IntervalEstimate_TestOfHypothesis(x,y1,test="var",ratio=2,side="g")double sample variance alternative hypothesis:true ratio of the variances is greater than 2$statistic_df_Pvalue f df1 df2 p-value 1 0.1380289 9 14 0.9973642$OneMinusAlpha_percent_confidence_interval low up 1 0.1043385 Inf$ratio_of_variances[1]0.2760579待测函数与var.test()结果相同，95%单侧置信区间为[0.1043385,Inf]，p-value=0.9973642＞0.05=alpha，接受原假设，认为方差比率小于或等于2。var.test(x,y1,ratio=2,al="g")F test to compare two variances data:x and y1 F=0.138,num df=9,denom df=14,p-value=0.9974 alternative hypothesis:true ratio of variances is greater than 2 95 percent confidence interval:0.1043385 Inf sample estimates:ratio of variances 0.2760579

3 结论

本文通过简要的理论分析，编写的函数IntervalEstimate_TestOfHypothesis()很好地解决了R软件内置函数t.test()、var.test()的缺陷，同时其参数设计也简单明了，将为需要作相关正态总体区间估计和假设检验的使用者提供方便。

[1]薛毅,陈丽萍.统计建模与R软件[M].北京:清华大学出版社,2007.

[2]陈希孺，倪国熙，陈长顺.数理统计学教程[M].安徽:中国科学技术大学出版社,2009.

[3]杨虎,刘琼荪,钟波.数理统计[M].北京:高等教育出版社,2004.

[4]王学民.多元应用分析[M].上海:上海财经大学出版社,2009.

[5]David Freedman等著,魏宗舒等译.统计学[M].北京:中国统计出版社,1997.

[6]郑忠国.高等统计学[M].北京:北京大学出版社,2012.

[7]R Core Team.R:A Language and Environment for Statistical Computing.R Foundation for Statistical Computing,Vienna,Austria[EB/OL].URL http://www.R-project.org/,2013.