先驗信息對MCMC方法估計概化理論方差分量變異量的影響

黎光明，張敏強

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

0 引言

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

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

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

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

一般地，僅根據一個樣本的統計量來估計總體參數，可能存在偏差。在樣本統計量研究中，僅用一個（次）樣本平均數來估計總體均值，存在較大的風險，因為樣本平均數容易受抽樣的影響。例如，某班某次考試平均分為70分，用這個成績估計這個班的能力水平，必然存在較大風險，這個班下次考試成績的平均分有可能變為75分。為了降低這種風險，通常的做法是用標準誤或置信區間等變異量來估計這個班的實際水平，如70±10或[60,80]。與平均數做法類似，MCMC方法所估計的方差分量，也受到抽樣的影響，用某個（次）樣本方差分量來估計總體方差分量（參數），存在一定誤差。為了降低這種誤差帶來的風險，需要報告方差分量對應的變異量（如標準誤或置信區間），來反映可能存在的變化程度。

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

國內一些學者[10]～[12]對使用MCMC方法估計模型的方差分量進行過研究，另有一些學者在概化理論中使用MCMC方法估計方差分量變異量[13,14]。但鮮有學者關注有無先驗信息對MCMC方法估計概化理論方差分量變異量的影響。先驗概率是MCMC方法估計各種統計量的必要條件之一，但在現實生活中，我們常常是沒有相關先驗信息的，這個時候，有先驗信息的MCMC方法優勢何在，先驗信息的有無設定，是否對MCMC方法估計概化理論方差分量變異量有影響等，都需要進行探討。

1 方法

1.1 數據模擬

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

（1）定義 p×i數學模型

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

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

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

1.2 先驗信息

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

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

1.3 比較標準

設定兩種比較標準[14]：一是方差分量標準誤估計的比較標準，為“偏差”（bias），計算方法是 bias=(θ^i-θ)，θ^i表示統計量的估計值，θ為參數，偏差的絕對值（稱為“絕對偏差”）越大，表明估計值與參數的差異越大，結果越不可靠；二是方差分量置信區間估計的比較標準，為“80%置信區間包含率”，包含率越接近0.80，表明結果越可靠。

1.4 分析工具

R軟件、WinBUGS軟件、R2WinBUGS軟件包和Coda軟件包。借助這些軟件或軟件包，自編完成研究程序。MCMC方法對方差分量及其變異量的估計是通過自編的R程序觸發WinBUGS軟件“間接”地實現。這個過程要求R軟件先調用R2WinBUGS和Coda兩個軟件包，R2Win-BUGS軟件包的作用在于成為R軟件和WinBUGS軟件的“橋梁”，Coda軟件包的作用在于輸出WinBUGS軟件生成的MCMC結果，如10%和90%兩分位點對應的方差分量。

2 結果

2.1 MCMC方法估計的方差分量及其標準誤

利用編制的程序估計三種不同樣本容量數據的方差分量及其標準誤，所得結果如表1所示。

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

MCMC inf和MCMC non-inf對應兩行數據，上一行表示估計的方差分量及其標準誤，下一行表示估計的方差分量及其標準誤與真值（參數）的偏差，即bias。

表1 方差分量及標準誤

2.2 MCMC方法估計的方差分量置信區間包含率

利用編制的程序估計三種不同樣本容量數據的方差分量置信區間包含率，所得結果如表2所示。

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

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

在表2中，括號中的數值表示估計的方差分量置信區間包含率與0.800（參數）的差值。

3 分析與討論

3.1 MCMC方法估計的方差分量標準誤偏差分析

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

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

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

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

3.2 MCMC方法估計的方差分量置信區間包含率分析

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

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

4 結論

（1）在方差分量標準誤這個變異量估計上，有先驗信息的MCMC方法較無先驗信息的MCMC方法要精確些，但隨著i的樣本容量增大，這種趨勢減小。

（2）在方差分量置信區間這個變異量估計上，隨著i的樣本容量增大，有先驗信息的MCMC方法和無先驗信息的MCMC方法的精確度相當。

[1]茆詩松,王靜龍,濮曉龍.高等數理統計[M].北京:高等教育出版社,1998.

[2]王權編譯.馬爾可夫鏈蒙特卡洛(MCMC)方法在估計IRT模型參數中的應用[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]魏宗舒.概率論與數理統計教程[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]楊志明,張雷.測評的概化理論及其應用[M].北京:教育科學出版社,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]馬躍淵,徐勇勇,奚苗苗,遇蘇寧.有缺失數據的生物等效性評價的MCMC方法[J].中國衛生統計,2004,21(4).

[11]馬躍淵.醫學數據統計分析中MCMC算法的實現與應用[D].第四軍醫大學碩士論文,2004.

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

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

[14]黎光明,張敏強.基于概化理論的方差分量變異量估計[N].心理學報,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.