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.