史海芳, 姬永剛

(1. 吉林大學(xué) 數(shù)學(xué)研究所, 長(zhǎng)春 130012; 2. 中國(guó)民航大學(xué) 理學(xué)院, 天津 300300)

0 引 言

考慮單向方差分析模型：

yij=μi+ij,j=1,2,…,ni,i=1,2,…,k.

(1)

其中: 響應(yīng)變量yij表示第i種試驗(yàn)第j個(gè)個(gè)體的值;μi為第i種試驗(yàn)的平均值; 隨機(jī)誤差ij為相互獨(dú)立的服從正態(tài)分布N(0,σ2)的隨機(jī)變量.

在很多實(shí)際問(wèn)題中, 通過(guò)一些先驗(yàn)信息可知總體均值滿(mǎn)足一定的序約束. 如簡(jiǎn)單半序μ1≤…≤μk, 簡(jiǎn)單樹(shù)半序μ1≤μi(i=2,3,…,k), 傘型半序μ1≤…≤μp≥μp+1≥…≥μk等. 文獻(xiàn)[1-2]給了在這些先驗(yàn)信息下一些感興趣參數(shù)的估計(jì)和檢驗(yàn). 此外, 文獻(xiàn)[3-5]分別考慮了構(gòu)造總體均值置信區(qū)間的問(wèn)題. 進(jìn)一步, 文獻(xiàn)[6]介紹了用U-I(union-intersection)檢驗(yàn)法考慮似然比檢驗(yàn)統(tǒng)計(jì)量的精確分布, 通過(guò)計(jì)算p值, 解決了簡(jiǎn)單半序約束下的檢驗(yàn)問(wèn)題. 但上述的頻率方法在理論上一般較復(fù)雜, 實(shí)際操作不太方便. 因此, 近年人們?cè)噲D從Bayes角度考慮序約束下的估計(jì)和檢驗(yàn)問(wèn)題, 如文獻(xiàn)[7-8]分別用Bayes方法考慮了如下等值檢驗(yàn)問(wèn)題:

H0:μ1=…=μkv.s.H1:μ1≤…≤μk,μ1<μk;

(2)

的估計(jì), 同時(shí)考慮了如下等值檢驗(yàn)問(wèn)題:

(3)

最后用Gibbs抽樣和Metropolis-Hastings方法給出了數(shù)值模擬. 與非Bayes方法相比, Bayes方法在處理中、 小樣本問(wèn)題時(shí)具有很大優(yōu)勢(shì). 特別地, 如果拒絕原假設(shè), 根據(jù)Bayes方法還可以知道哪些總體均值和標(biāo)準(zhǔn)差之比是不同的.

1 Bayes分層模型

[δhρh,θh]=ρhI(δh=0)+(1-ρh)g(δhθh)I(δh>0),h=1,2,…,k-1,

(4)

(5)

綜上, Bayes分層模型為:

2 滿(mǎn)條件后驗(yàn)分布

由Bayes分層模型知, 所有未知參數(shù)包括δj,ρj,θj,ξ1,σi(j=1,2,…,k-1;i=1,2,…,k).

2.1 δj的滿(mǎn)條件后驗(yàn)分布

由Bayes公式知

因此, 給定y,{σi},ρj,θj,ξ1,δ-j下,δj的滿(mǎn)條件后驗(yàn)分布是一個(gè)混合分布:

其中

(7)

2.2 ξ1的滿(mǎn)條件后驗(yàn)分布

由Bayes公式知

2.3 ρi的滿(mǎn)條件后驗(yàn)分布

由Bayes公式知

ρi的滿(mǎn)條件后驗(yàn)分布可以記作

當(dāng)α0=β0=1時(shí), 有

2.4 θi的滿(mǎn)條件后驗(yàn)分布

由Bayes公式知

θi的滿(mǎn)條件后驗(yàn)分布可以記作

2.5 σi的滿(mǎn)條件后驗(yàn)分布

由Bayes公式, 有

顯然,σi的滿(mǎn)條件后驗(yàn)不是一個(gè)常見(jiàn)分布, 所以本文用Metropolis-Hastings方法給出σi的數(shù)值模擬.

3 數(shù)值模擬

3.1 Gibbs抽樣

Gibbs抽樣步驟如下：

2) ① 先通過(guò)下式計(jì)算fj和ej:

(8)

解方程組(8)得

② 計(jì)算r.

③ 從(0,1)上的均勻分布抽取一隨機(jī)數(shù)u:

3.2 數(shù)值模擬

下面通過(guò)數(shù)值模擬檢驗(yàn)本文的方法. 考慮總體個(gè)數(shù)k=4時(shí)的情況, 已知δj=μj+1/σj+1-μj/σj,j=1,2,3. 共模擬4組總體. 為方便, 令每組總體中均值和標(biāo)準(zhǔn)差比的間隔δj都分別相等, 分別為0.5,0.75,1,1.25. 因此產(chǎn)生如下4組數(shù)據(jù):

表1 不同樣本數(shù)下檢驗(yàn)問(wèn)題(3)的勢(shì)

表2～表4分別列出了一次抽樣中參數(shù)的估計(jì)值和δj=0的后驗(yàn)概率, 其中Bayes估計(jì)是本文方法給出的估計(jì), Frequentist估計(jì)是利用文獻(xiàn)[10]中頻率方法得到的參數(shù)估計(jì)值. 由表2可見(jiàn), 對(duì)于參數(shù)δ1,δ2和δ3, 當(dāng)n=30,100時(shí), 本文方法都給出了比頻率方法更精確的估計(jì). 盡管當(dāng)n=30時(shí), 頻率方法給出的參數(shù)ξ估計(jì)值比本文方法更接近真值, 但當(dāng)n=100 時(shí), 本文方法給出的估計(jì)值更接近真值. 由表3可見(jiàn), 本文方法給出的參數(shù)ξ,δ1,δ3估計(jì)值優(yōu)于頻率方法給出的估計(jì)值, 但對(duì)于參數(shù)δ2, 頻率方法給出了比本文方法更精確的估計(jì). 由表4可見(jiàn), 只有參數(shù)δ1當(dāng)n=30時(shí)和參數(shù)δ3當(dāng)n=100時(shí), 頻率方法優(yōu)于本文方法, 其他情況都是本文方法優(yōu)于頻率方法. 總之, 兩種方法給出的估計(jì)值類(lèi)似, 都與真值很接近. 對(duì)于檢驗(yàn)問(wèn)題(3), 從表中數(shù)據(jù)可見(jiàn)應(yīng)用本文給出的Bayes方法得到的結(jié)果較理想. 其中在表3中, 盡管當(dāng)間隔δj較小且樣本量不大(n=30)時(shí)出現(xiàn)了與實(shí)際不符的檢驗(yàn)結(jié)果(P(δ2=0Y)=0.611 7), 但隨著樣本的增多或間隔的增大(表4), 從后驗(yàn)概率都小于0.5易見(jiàn)4個(gè)正態(tài)總體均值與標(biāo)準(zhǔn)差比都不同. 這主要是由于間隔較小時(shí)很難正確識(shí)別出所有不同的均值和標(biāo)準(zhǔn)差的比.

表2 當(dāng)μ=(1,1,1,1), σ=(1,1,1,1)時(shí)一次抽樣中參數(shù)的估計(jì)值和δj=0的后驗(yàn)概率

表3 當(dāng)μ=(0,1,2,3), σ=(1,1,1,1)時(shí)一次抽樣中參數(shù)的估計(jì)值和δj=0的后驗(yàn)概率

表4 當(dāng)μ=(0,2,4,6), σ=(1,1,1,1)時(shí)一次抽樣中參數(shù)的估計(jì)值和δj=0的后驗(yàn)概率

進(jìn)一步, 當(dāng)原假設(shè)成立時(shí), 本文對(duì)犯第一類(lèi)錯(cuò)誤的概率進(jìn)行了模擬, 結(jié)果列于表5. 考慮4組總體, 數(shù)據(jù)按如下總體產(chǎn)生:

其中σ=(σ1,σ2,σ3,σ4)分別取(1,1,1,1),(3,3,3,3),(6,6,6,6),(10,10,10,10).

表5 不同樣本數(shù)和方差下犯第一類(lèi)錯(cuò)誤的概率

由表5可見(jiàn), 本文方法能很好地控制犯第一類(lèi)錯(cuò)誤的概率, 甚至當(dāng)標(biāo)準(zhǔn)差很大(σ=10)時(shí), 犯第一類(lèi)錯(cuò)誤的概率仍然很小, 因此本文方法是一個(gè)相對(duì)保守的方法. 此外, 在表5中, 當(dāng)σ=(1,1,1,1)時(shí)似乎隨著樣本量的增加犯第一類(lèi)錯(cuò)誤的概率也增加, 為了檢驗(yàn)這一結(jié)論本文做了很多模擬, 結(jié)果顯示當(dāng)n=49時(shí), 犯第一類(lèi)錯(cuò)誤的概率為0.005 0, 當(dāng)n=60,70,80,90,110,120時(shí), 犯第一類(lèi)錯(cuò)誤的概率均為0, 因此, 并不存在犯第一類(lèi)錯(cuò)誤的概率隨樣本量增加而增大的趨勢(shì). 當(dāng)間隔δ不是太小時(shí), 本文的Bayes方法在處理檢驗(yàn)問(wèn)題(3)時(shí), 在控制犯第一類(lèi)錯(cuò)誤很小的同時(shí), 也能達(dá)到一個(gè)合理的檢驗(yàn)勢(shì).

[1] Robertson T, Wright F T, Dykstra R L. Order Restricted Statistical Inference [M]. New Jersey: John Wiley & Sons Inc, 1988.

[2] Silvapulle M J, Sen P K. Constrained Statistical Inference: Inequality, Order and Shape Restrictions [M]. New Jersey: John Wiley & Sons Inc, 2005.

[3] Hayter A J. A One-Sided Studentized Range Test for Testing against a Simple Order Alternative [J]. Journal of the American Statistical Association, 1990, 85(411): 778-785.

[4] LIU Lin. Simultaneous Statistical Inference for Monotone Dose-Response Means [D]. Newfoundlan: Memorial Universtiy of Newfoundland, 2001.

[5] LIU Lin, Lee C C, PENG Jian-an. Max-Min Multiple Comparison Procedure for Isotonic Dose-Response Curves [J]. Journal of Statistical Planning and Inference, 2002, 107(1/2): 133-141.

[6] 史寧中. 統(tǒng)計(jì)檢驗(yàn)的理論與方法 [M]. 北京: 科學(xué)出版社, 2008.

[7] SHANG Jun-feng, Cavanaugh J E, Wright F T. A Bayesian Multiple Comparison Procedure for Order-Restricted Mixed Models [J]. International Statistical Review, 2008, 76(2): 268-284.

[8] Oh M S, Shin D W. A Unified Bayesian Inference on Treatment Means with Order Constraints [J]. Computational Statistics & Data Analysis, 2011, 55(1): 924-934.

[9] Chou Y M, Owen D B. A Likelihood Ratio Test for the Equality of Proportions of Two Normal Populations [J]. Communications in Statistics: Theory and Methods, 1991, 20(8): 2357-2374.

[10] LI Shu-you, SHI Ning-zhong, ZHANG Bao-xue. Testing Ratios of Means to Standard Deviations from Normal Populations under Order Restrictions with Generalizedp-Values [J]. Chinese Journal of Applied Probability and Statistics, 2009, 25(1): 77-85. (李樹(shù)有, 史寧中, 張寶學(xué). 正態(tài)總體均值與標(biāo)準(zhǔn)差比在序約束下的廣義p-值檢驗(yàn) [J]. 應(yīng)用概率統(tǒng)計(jì), 2009, 25(1): 77-85.)