先驗(yàn)信息對(duì)MCMC方法估計(jì)概化理論方差分量變異量的影響

黎光明，張敏強(qiáng)

(1.廣州大學(xué) 教育學(xué)院，廣州 510006；2.華南師范大學(xué) 心理應(yīng)用研究中心，廣州 510631)

0 引言

MCMC(Markov Chain Monte Carlo)方法源于物理學(xué)研究，在統(tǒng)計(jì)物理中得到廣泛應(yīng)用已有四十多年的歷史[1]。MCMC方法是一種動(dòng)態(tài)的計(jì)算機(jī)模擬技術(shù)，根據(jù)貝葉斯（Bayes）推斷為中心的多元后驗(yàn)分布，來模擬隨機(jī)樣本的一種方法[2]。上世紀(jì)50年代，統(tǒng)計(jì)物理學(xué)家Metropolis等人引入MCMC，它是一種動(dòng)態(tài)的蒙特卡洛方法[3]。大規(guī)模使用MCMC是在上世紀(jì)90年代以后，與傳統(tǒng)的Monte Carlo方法相比，這種方法可以解決觀測分布是復(fù)雜的、高維的、非標(biāo)準(zhǔn)化形式的分布[1]。MCMC方法是通過模擬服從某一分布也即平穩(wěn)分布（Stationary Distribution）(一般是待估參數(shù)的聯(lián)合后驗(yàn)分布)的馬爾可夫鏈，然后根據(jù)模擬的馬爾可夫鏈上的樣本點(diǎn)對(duì)待估參數(shù)進(jìn)行估計(jì)。在給定數(shù)據(jù)的條件下，使用特定的轉(zhuǎn)移核（Markov Kernels）來模擬未知參數(shù)的分布，MCMC鏈就產(chǎn)生在這個(gè)迭代過程中。進(jìn)行MCMC可解決高維積分問題，特別可應(yīng)用于貝葉斯公式計(jì)算，貝葉斯公式[4]可表示如下：

其中，θ=(θ1,θ2,…)，p(θ|X)表示后驗(yàn)條件概率，X為原始分?jǐn)?shù)（本研究原始分?jǐn)?shù)表示模擬生成的數(shù)據(jù)），θ表示待估的統(tǒng)計(jì)量，p(θ)為先驗(yàn)的概率，∫? p(x|θ)p(θ)dθ 為需要高維積分的部分。很顯然，如果待估參數(shù)較少且形式簡單，傳統(tǒng)的Monte Carlo就可以解決。但是，情況往往并非如此，很多情況下會(huì)采用MCMC方法，因?yàn)楣剑?）高維積分較為困難，就不得不使用其它方法，如MCMC方法。用MCMC方法生成多條鏈，通過每條鏈分別構(gòu)造轉(zhuǎn)移核以使最后能達(dá)到平穩(wěn)分布，然后求取平穩(wěn)分布上統(tǒng)計(jì)量的平均值來估計(jì)相應(yīng)的統(tǒng)計(jì)量。不管最初值如何，如果MCMC多條鏈能夠收斂，最后都能得到一個(gè)常數(shù)值。因?yàn)榻?jīng)過高維積分后貝葉斯公式分母變成一個(gè)常數(shù)（用k表示），那么貝葉斯公式就可表示為：

從公式（2）可知，后驗(yàn)概率∝似然函數(shù)（likelihood）×先驗(yàn)概率，從此關(guān)系可以看出，為了得到后驗(yàn)概率，MCMC方法應(yīng)該具備三個(gè)必要條件：初始值、似然函數(shù)和先驗(yàn)概率。如果先驗(yàn)概率分布已知，則MCMC方法稱為有先驗(yàn)信息的 MCMC方法（MCMC with informative priors，MCMC inf），否則為無先驗(yàn)信息的MCMC方法（MCMC with non-informative priors，MCMC non-inf）。

Geman和Geman(1984)引入Gibbs算法[5]，現(xiàn)已成為一種標(biāo)準(zhǔn)化的統(tǒng)計(jì)計(jì)算工具。進(jìn)行MCMC的Gibbs算法典型軟件是WinBUGS軟件(Windows operating system version of BUGS:Bayesian Analysis Using Gibbs Sampling)。MCMC方法可以估計(jì)模型的各種統(tǒng)計(jì)量，如平均數(shù)、標(biāo)準(zhǔn)差、百分位數(shù)、方差分量等。方差分量（variance component）是指矩陣中的各元素都由若干成分組成，每個(gè)分量有各自含義，某一分量占總分量的比例可以用來說明該分量在總方差或協(xié)方差中所起的作用。

一般地，僅根據(jù)一個(gè)樣本的統(tǒng)計(jì)量來估計(jì)總體參數(shù)，可能存在偏差。在樣本統(tǒng)計(jì)量研究中，僅用一個(gè)（次）樣本平均數(shù)來估計(jì)總體均值，存在較大的風(fēng)險(xiǎn)，因?yàn)闃颖酒骄鶖?shù)容易受抽樣的影響。例如，某班某次考試平均分為70分，用這個(gè)成績估計(jì)這個(gè)班的能力水平，必然存在較大風(fēng)險(xiǎn)，這個(gè)班下次考試成績的平均分有可能變?yōu)?5分。為了降低這種風(fēng)險(xiǎn)，通常的做法是用標(biāo)準(zhǔn)誤或置信區(qū)間等變異量來估計(jì)這個(gè)班的實(shí)際水平，如70±10或[60,80]。與平均數(shù)做法類似，MCMC方法所估計(jì)的方差分量，也受到抽樣的影響，用某個(gè)（次）樣本方差分量來估計(jì)總體方差分量（參數(shù)），存在一定誤差。為了降低這種誤差帶來的風(fēng)險(xiǎn)，需要報(bào)告方差分量對(duì)應(yīng)的變異量（如標(biāo)準(zhǔn)誤或置信區(qū)間），來反映可能存在的變化程度。

概化理論是現(xiàn)代心理與教育測量理論之一[6]。方差分量估計(jì)是概化理論的必用技術(shù)，是進(jìn)行概化理論分析的關(guān)鍵。與其它統(tǒng)計(jì)量一樣，依據(jù)概化理論模型所估計(jì)出的方差分量受限于抽樣，不同的抽樣樣本，所估計(jì)的方差分量可能不一樣，這就要求進(jìn)行方差分量估計(jì)時(shí)需要對(duì)其變異量進(jìn)行探討。Tong和Brennan(2006)[7]認(rèn)為，對(duì)方差分量變異量的估計(jì)可以使用馬爾可夫鏈蒙特卡洛(MCMC)方法。探討方差分量變異量具有重要意義，因?yàn)閳?bào)告這些變異量可以在一定程度上說明方差分量測量的可靠性[8,9]。

國內(nèi)一些學(xué)者[10]～[12]對(duì)使用MCMC方法估計(jì)模型的方差分量進(jìn)行過研究，另有一些學(xué)者在概化理論中使用MCMC方法估計(jì)方差分量變異量[13,14]。但鮮有學(xué)者關(guān)注有無先驗(yàn)信息對(duì)MCMC方法估計(jì)概化理論方差分量變異量的影響。先驗(yàn)概率是MCMC方法估計(jì)各種統(tǒng)計(jì)量的必要條件之一，但在現(xiàn)實(shí)生活中，我們常常是沒有相關(guān)先驗(yàn)信息的，這個(gè)時(shí)候，有先驗(yàn)信息的MCMC方法優(yōu)勢何在，先驗(yàn)信息的有無設(shè)定，是否對(duì)MCMC方法估計(jì)概化理論方差分量變異量有影響等，都需要進(jìn)行探討。

1 方法

1.1 數(shù)據(jù)模擬

運(yùn)用Monte Carlo數(shù)據(jù)模擬技術(shù)，產(chǎn)生概化理論 p×i（person×item）矩陣數(shù)據(jù)。數(shù)據(jù)模擬所使用的軟件為R軟件，產(chǎn)生模擬數(shù)據(jù)的過程如下[15]：

（1）定義 p×i數(shù)學(xué)模型

在式（3）和式（4）中，μp-μ=πp表示 p的效應(yīng)，μi-μ=βi表示i的效應(yīng)，Xpi-μp-μi+μ=εpi表示 pi的效應(yīng)。

（2）將公式（4）轉(zhuǎn)換成 Xpi=μ+σpzp+σizi+σpizpi，在參數(shù)σp、σi和σpi已知的情況下（μ通常設(shè)定為0），調(diào)用R軟件中的rnorm函數(shù)隨機(jī)生成三個(gè)服從正態(tài)分布的zp、zi和zpi。

（3）連乘相加獲得模擬數(shù)據(jù)Xpi。假定為參數(shù)，其值分別設(shè)定為4、16和64。分別構(gòu)建 p×i設(shè)計(jì)不同樣本容量的三種情況，即100×20、100×40、100×80。生成的模擬數(shù)據(jù)為矩陣數(shù)據(jù)（p×i），模擬次數(shù)為1000。

1.2 先驗(yàn)信息

依據(jù)概化理論模型 Xpi=μ+πp+βi+εpi，為了獲得后驗(yàn)概率，參考Mao,Shin和Brennan（2005）的做法[16]，定義似然函數(shù)如下：

通過式（6）～（8）設(shè)定模型的先驗(yàn)分布，根據(jù)共軛分布性質(zhì)，方差分量對(duì)應(yīng)的分布為逆τ分布，且τ=1/σ2。為定義無先驗(yàn)信息的分布，σ常取0.001，那么τ=106，即認(rèn)為τ為無窮大。如果為有先驗(yàn)信息的分布，則依具體的先驗(yàn) 值 而 定 ，本 研 究 設(shè) 定 p～ τ(2,4) ，i～ τ(2,16) ，pi～τ(2,64)，初始值均為0.001。

1.3 比較標(biāo)準(zhǔn)

設(shè)定兩種比較標(biāo)準(zhǔn)[14]：一是方差分量標(biāo)準(zhǔn)誤估計(jì)的比較標(biāo)準(zhǔn)，為“偏差”（bias），計(jì)算方法是 bias=(θ^i-θ)，θ^i表示統(tǒng)計(jì)量的估計(jì)值，θ為參數(shù)，偏差的絕對(duì)值（稱為“絕對(duì)偏差”）越大，表明估計(jì)值與參數(shù)的差異越大，結(jié)果越不可靠；二是方差分量置信區(qū)間估計(jì)的比較標(biāo)準(zhǔn)，為“80%置信區(qū)間包含率”，包含率越接近0.80，表明結(jié)果越可靠。

1.4 分析工具

R軟件、WinBUGS軟件、R2WinBUGS軟件包和Coda軟件包。借助這些軟件或軟件包，自編完成研究程序。MCMC方法對(duì)方差分量及其變異量的估計(jì)是通過自編的R程序觸發(fā)WinBUGS軟件“間接”地實(shí)現(xiàn)。這個(gè)過程要求R軟件先調(diào)用R2WinBUGS和Coda兩個(gè)軟件包，R2Win-BUGS軟件包的作用在于成為R軟件和WinBUGS軟件的“橋梁”，Coda軟件包的作用在于輸出WinBUGS軟件生成的MCMC結(jié)果，如10%和90%兩分位點(diǎn)對(duì)應(yīng)的方差分量。

2 結(jié)果

2.1 MCMC方法估計(jì)的方差分量及其標(biāo)準(zhǔn)誤

利用編制的程序估計(jì)三種不同樣本容量數(shù)據(jù)的方差分量及其標(biāo)準(zhǔn)誤，所得結(jié)果如表1所示。

在表1中，np表示人的樣本容量，ni表示題目的樣本容量，parameters表示參數(shù)，vc.p、vc.i和vc.pi分別表示人的方差分量、題目的方差分量以及人與題目交互（包括殘差）的方差分量，SE(vc.p)、SE(vc.i)和SE(vc.pi)分別表示人的方差分量標(biāo)準(zhǔn)誤、題目的方差分量標(biāo)準(zhǔn)誤、人與題目交互（包括殘差）的方差分量標(biāo)準(zhǔn)誤。MCMC inf和MCMC non-inf分別表示有先驗(yàn)信息的MCMC方法和無先驗(yàn)信息的MCMC方法。

MCMC inf和MCMC non-inf對(duì)應(yīng)兩行數(shù)據(jù)，上一行表示估計(jì)的方差分量及其標(biāo)準(zhǔn)誤，下一行表示估計(jì)的方差分量及其標(biāo)準(zhǔn)誤與真值（參數(shù)）的偏差，即bias。

表1 方差分量及標(biāo)準(zhǔn)誤

2.2 MCMC方法估計(jì)的方差分量置信區(qū)間包含率

利用編制的程序估計(jì)三種不同樣本容量數(shù)據(jù)的方差分量置信區(qū)間包含率，所得結(jié)果如表2所示。

表2 MCMC方法不同樣本容量下估計(jì)的方差分量置信區(qū)間包含率

在表2中，np表示人的樣本容量，ni表示題目的樣本容量，coverage（80%）表示“80%置信區(qū)間包含率”。計(jì)算80%置信區(qū)間包含率的方法是：通過判斷參數(shù)是否落入10%到90%兩分位點(diǎn)對(duì)應(yīng)的方差分量之間，如果某次成功，則包含次數(shù)加1，最后計(jì)算落入的總次數(shù)，并除以1000，即為最后的包含率。MCMC inf和MCMC non-inf分別表示有先驗(yàn)信息的MCMC方法和無先驗(yàn)信息的MCMC方法。

在表2中，括號(hào)中的數(shù)值表示估計(jì)的方差分量置信區(qū)間包含率與0.800（參數(shù)）的差值。

3 分析與討論

3.1 MCMC方法估計(jì)的方差分量標(biāo)準(zhǔn)誤偏差分析

從表1可以看出，p的樣本容量固定為100，i的樣本容量逐漸增大（20、40、80）。對(duì)于方差分量標(biāo)準(zhǔn)誤參數(shù)，其值隨著i的樣本容量增大而趨于減小。例如，i的樣本容量為20、40、80時(shí)，三個(gè)方差分量標(biāo)準(zhǔn)誤參數(shù)分別為1.0287、0.7968、0.6824，這些值逐漸減小。MCMC inf和MCMC non-inf估計(jì)的方差分量標(biāo)準(zhǔn)誤隨著i的樣本容量增大也趨于減小，例如，i的樣本容量為20、40、80時(shí)，MCMC non-inf估計(jì)p的方差分量標(biāo)準(zhǔn)誤分別為1.0443、0.8369、0.7149，這些值逐漸減小。方差分量標(biāo)準(zhǔn)誤隨著i的樣本容量增大而趨于減小，符合標(biāo)準(zhǔn)誤與樣本容量之間的關(guān)系。

進(jìn)一步分析表1可以看出，MCMC兩種方法估計(jì)的方差分量標(biāo)準(zhǔn)誤偏差主要表現(xiàn)在（i題目）上。例如，當(dāng)np=100，ni=20時(shí)，MCMC inf估計(jì)的SE(vc.p)、SE(vc.i)、SE(vc.pi)的偏差分別-0.0340、0.0487、-0.0181，SE(vc.i)的偏差最大。再如，當(dāng)np=100，ni=40時(shí)，MCMC non-inf估計(jì)的 SE(vc.p)、SE(vc.i)、SE(vc.pi)的偏差分別0.0401、0.2378、-0.0110，SE(vc.i)的偏差最大。因此，考察兩種MCMC方法估計(jì)的題目的方差分量標(biāo)準(zhǔn)誤偏差，更為重要。

比較表1中MCMC inf和MCMC non-inf估計(jì)的題目的方差分量標(biāo)準(zhǔn)誤偏差，發(fā)現(xiàn)前者小于后者，表明有先驗(yàn)信息的MCMC方法較無先驗(yàn)信息的MCMC方法對(duì)方差分量標(biāo)準(zhǔn)誤的估計(jì)要精確些。例如，當(dāng)題目樣本容量增大時(shí)（20、40、80），MCMC inf估計(jì)SE(vc.i)的偏差分別為0.0487、-0.1090、0.0859（已用虛線方框括起），MCMC non-inf估計(jì)SE(vc.i)的偏差分別為1.0434、0.2378、0.1059（已用實(shí)線方框括起），明顯可以看出，MCMC non-inf估計(jì)的SE(vc.i)的偏差要大于MCMC inf估計(jì)的SE(vc.i)的偏差。

但是，隨著i的樣本容量增大，MCMC non-inf和MCMC inf估計(jì)的題目的方差分量標(biāo)準(zhǔn)誤偏差，趨于接近。例如，當(dāng)np=100，ni=20時(shí)，兩種方法估計(jì)的SE(vc.i)的偏差分別為1.0434和0.0487，偏差差值為0.9947；當(dāng) np=100，ni=40時(shí)，兩種方法估計(jì)的SE(vc.i)的偏差分別為0.2378和-0.1090，偏差差值為0.3468；當(dāng)np=100，ni=80時(shí)，兩種方法估計(jì)的SE(vc.i)的偏差分別為0.1059和0.0859，偏差差值為0.0200。

3.2 MCMC方法估計(jì)的方差分量置信區(qū)間包含率分析

從表2可知，有先驗(yàn)信息的MCMC方法和無先驗(yàn)信息的MCMC方法隨著i的樣本容量增大，估計(jì)的方差分量置信區(qū)間包含率并沒有表現(xiàn)出一定的規(guī)律性。但是，從兩種MCMC方法估計(jì)的置信區(qū)間包含率看，兩種方法都接近于0.800，并且包含率與0.800的差值的絕對(duì)值在0.000～0.031之間，差值的絕對(duì)值較小。因此，可以認(rèn)為兩種MCMC方法都較好地估計(jì)了方差分量的置信區(qū)間。

進(jìn)一步分析表2可以看出，隨著i的樣本容量增大（20、40、80），無先驗(yàn)信息的MCMC方法并不是在所有方差分量置信區(qū)間估計(jì)上都優(yōu)于有先驗(yàn)信息的MCMC方法，例如，當(dāng)np=100,ni=80時(shí)，有先驗(yàn)信息的MCMC方法估計(jì)的p的方差分量置信區(qū)間包含率為0.802(0.002)（已用方框括起），而無先驗(yàn)信息的MCMC方法為0.811(0.011)，有先驗(yàn)信息的MCMC方法好些。對(duì)于估計(jì)的方差分量置信區(qū)間包含率與0.800的差值，從總體上看，有先驗(yàn)信息的MCMC方法與無先驗(yàn)信息的MCMC方法相距較小，表明兩種方法估計(jì)方差分量置信區(qū)間精確度相當(dāng)。這就是說，在變化樣本容量的條件下，有時(shí)無先驗(yàn)信息的MCMC方法估計(jì)某些方差分量置信區(qū)間的結(jié)果好些，有時(shí)卻是有先驗(yàn)信息的MCMC方法好些。但從總體上看，兩種方法精確度相當(dāng)。

4 結(jié)論

（1）在方差分量標(biāo)準(zhǔn)誤這個(gè)變異量估計(jì)上，有先驗(yàn)信息的MCMC方法較無先驗(yàn)信息的MCMC方法要精確些，但隨著i的樣本容量增大，這種趨勢減小。

（2）在方差分量置信區(qū)間這個(gè)變異量估計(jì)上，隨著i的樣本容量增大，有先驗(yàn)信息的MCMC方法和無先驗(yàn)信息的MCMC方法的精確度相當(dāng)。

[1]茆詩松,王靜龍,濮曉龍.高等數(shù)理統(tǒng)計(jì)[M].北京:高等教育出版社,1998.

[2]王權(quán)編譯.馬爾可夫鏈蒙特卡洛(MCMC)方法在估計(jì)IRT模型參數(shù)中的應(yīng)用[J].考試研究,2006,(4).

[3]Metropolis,N.,Rosenbluth,A.W.,Rosenbluth,M.N.,Teller,A.H.Teller,E.Equations of State Calculations by Fast Computing Machines[J].Journal of Chemical Physics,1953,(21).

[4]魏宗舒.概率論與數(shù)理統(tǒng)計(jì)教程[M].北京:高等教育出版社,1983.

[5]German,S.,German,D.Stochastic Relaxation,Gibbs Distribution and the Bayesian Restoration of Images[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1984,(6).

[6]楊志明,張雷.測評(píng)的概化理論及其應(yīng)用[M].北京:教育科學(xué)出版社,2003.

[7]Tong,Y.,Brennan,R.L.Bootstrap Techniques for Estimating in Gener?alizability Theory(CASMA Research Report No 15).Iowa City,IA:Center for Advanced Studies in Measurement and Assessment,Uni?versity of Iowa[EB/OL].http://www.education.uiowa.edu/casma,2006.

[8]American Educational Research Association.American Psychological Association,National Council on Measurement in Education.Stan?dards for Educational and Psychological Testing[Z].Washington,DC,1985.

[9]American Educational Research Association.American Psychological Association,National Council on Measurement in Education.Stan?dards for Educational and Psychological Testing(Rev.ed.)[Z].Wash?ington,DC,1999.

[10]馬躍淵,徐勇勇,奚苗苗,遇蘇寧.有缺失數(shù)據(jù)的生物等效性評(píng)價(jià)的MCMC方法[J].中國衛(wèi)生統(tǒng)計(jì),2004,21(4).

[11]馬躍淵.醫(yī)學(xué)數(shù)據(jù)統(tǒng)計(jì)分析中MCMC算法的實(shí)現(xiàn)與應(yīng)用[D].第四軍醫(yī)大學(xué)碩士論文,2004.

[12]郜艷暉,姜慶五,孟煒,陳啟明,趙耐青,沈福民.REML法和MCMC法在數(shù)量性狀核心家系遺傳方差分量模型中的參數(shù)估計(jì)的比較[J].復(fù)旦學(xué)報(bào),2003,30(4).

[13]黎光明,張敏強(qiáng).一項(xiàng)跨分布研究:基于概化理論的方差分量變異量估計(jì)[D].北京師范大學(xué)首屆全國心理學(xué)博士學(xué)術(shù)論壇論文集,2009a.

[14]黎光明,張敏強(qiáng).基于概化理論的方差分量變異量估計(jì)[N].心理學(xué)報(bào),2009b,41(9).

[15]Brennan,R.L.Generalizability Theory[M].New York:Springer-Ver?lag,2001.

[16]Mao,X.,Shin,D.,Brennan,R.L.Estimating the Variability of Esti?mated Variance Components and Related Statistics Using the MC?MC Procedure:An Exploratory Study[C].Paper Presented at the An?nual Meeting of the National Council on Measurement in Education,Montreal,2005.