周國瓊, 蔣文江

(1.昭通學院 數(shù)學與統(tǒng)計學院，云南 昭通 657000 ；2.云南師范大學 泛亞商學院，云南 昆明 650092)

在傳統(tǒng)的因子分析模型中，因子被假設為服從標準正態(tài)分布，可取實數(shù)域上的任何值；而在日常生活中，往往存在一些非負或非正的數(shù)據(jù)，如學生成績和壽命數(shù)據(jù)等，當對這類數(shù)據(jù)做因子分析時，假設因子服從一個非負的概率分布(如指數(shù)分布)是較為合理的選擇.

傳統(tǒng)的因子分析模型假設模型中各個公共因子及特殊因子之間相互獨立且公共因子服從標準正態(tài)分布，因而可以直接利用極大似然估計[1-2]對模型中的參數(shù)做估計；但當因子分析模型的因子服從指數(shù)分布時，會導致模型中的似然函數(shù)沒有顯式表達，不能直接使用極大似然估計法估計參數(shù)[3]；針對這一情況，本文采用蒙特卡洛方法來解決極大似然估計中似然函數(shù)沒有顯式表達的問題，即利用基于EM算法的極大似然估計法對模型中的參數(shù)做估計[4-6]，其中EM算法中的E步采用馬爾科夫鏈蒙特卡洛方法中的M-H算法[7]從一個非常規(guī)的復雜分布中通過抽樣來完成積分計算[8].

1 模型及待估參數(shù)

設y=(y1,y2,…,yp)T是p維可觀測的隨機向量，μ=(μ1，μ2，…，μp)T是p維截距向量，Λ=(λij)p×q是p×q(p>q)維的因子載荷矩陣，x=(x1,x2,…,xq)T是q維潛在因子的隨機向量，ε=(ε1，ε2，…，εp)T是p維誤差隨機向量，本文所研究的模型為

y=μ+Λx+ε

(1)

該模型與傳統(tǒng)的因子分析模型[1]的不同之處是本文假設誤差隨機向量ε～N(0,∑)，其中∑=(σij)p×p為對角矩陣；潛在因子的隨機向量x中每一個隨機變量xk～exp(βk),k=1,2,…,q；誤差隨機向量的分量εj(j=1,2,…,p)與潛在因子隨機向量的分量xk之間相互獨立.模型中的待估參數(shù)有μ,Λ,Σ,βk,本文的核心工作就是對這些參數(shù)進行估計.

2 參數(shù)估計研究

2.1 模型參數(shù)極大似然估計的方法與原理

把模型(1)中的待估參數(shù)記為參數(shù)向量θ=(μ，Λ，∑，βk)，記Y=(y1,y2,…,yn)為觀測到的數(shù)據(jù)矩陣，X=(x1,x2,…,xn)是潛在的因子矩陣.

2.1.1 模型的似然函數(shù)

由于潛在因子隨機向量x的每一個分量之間相互獨立且xk～exp(βk),所以隨機向量x的聯(lián)合密度函數(shù)

(2)

由于誤差向量ε～N(0,∑)，所以當潛在的因子隨機向量已知時，可觀測隨機向量y～N(μ+Λx,∑)，于是根據(jù)條件概率的定義可得y的條件概率密度函數(shù)[9]

(3)

根據(jù)式(2)和式(3)，基于觀測數(shù)據(jù)Y的似然函數(shù)

(4)

對應的對數(shù)似然函數(shù)

(5)

(5)中的多重積分沒有顯式表達，所以很難通過極大化(5)來獲得參數(shù)θ的極大似然估計[10].

為解決上述積分問題，考慮使用數(shù)據(jù)添加算法[11]，把潛在的因子矩陣X=(x1,x2,…,xn)看作缺失的數(shù)據(jù)添加到觀測到的數(shù)據(jù)矩陣Y=(y1,y2,…,yn)中，從而得到完全數(shù)據(jù)矩陣Z=(X,Y)，然后用EM算法來獲得參數(shù)θ的極大似然估計.

2.1.2 EM算法估計參數(shù)

記Z=(X,Y)為添加數(shù)據(jù)后的完全數(shù)據(jù)矩陣，由于把潛在的不可直接觀測到的因子矩陣當作是已知的能夠觀測到的數(shù)據(jù)矩陣添加到觀測數(shù)據(jù)Y中，式(4)和式(5)中對X的積分運算失效，則添加數(shù)據(jù)后的對數(shù)似然函數(shù)

(6)

通過EM算法可求得式(6)中參數(shù)θ的極大似然估計.

根據(jù)EM算法的原理[12]，基于θ的第t次迭代值對對數(shù)似然函數(shù)(6)中的X求期望，得到Q函數(shù)

(7)

需要計算的積分有E[xi|θ(t),Y]和E[xixiT|θ(t),Y].

利用蒙特卡洛方法，通過從密度函數(shù)p(x|y,θ)中抽取M個x的樣本x(l),l=1,2,…,M，然后根據(jù)樣本觀測值對x求期望，得

由貝葉斯公式[13]可得密度函數(shù)

(8)

該密度函數(shù)并不是某一常見分布[14]的密度函數(shù)，所以想要用常規(guī)的方法從(8)中抽取x的樣本非常困難[15].本文采用馬爾科夫鏈蒙特卡洛(MCMC)方法中的Metropolis-Hastings(M-H)算法來對x進行抽樣.

根據(jù)Metropolis-Hastings算法的原理[16]，視密度函數(shù)p(x|y,θ)為產(chǎn)生樣本x的目標密度，然后找一個容易產(chǎn)生樣本的分布作為建議分布，通常選取正態(tài)分布N(x(l-1),σ2Ω)作為建議分布[17]，其中x(l-1)是產(chǎn)生第l個x樣本的上一步第l-1步所產(chǎn)生的x樣本，Ω為目標密度關于x的Fisher信息矩陣，σ2是任一給定的調(diào)節(jié)參數(shù)，通過調(diào)整不同的σ2值來得到樣本不同的接受概率，通常需要把樣本的接受概率控制在0.25到0.50之間[18]，而目標分布與建議分布之間通過信息矩陣Ω進行聯(lián)系.

根據(jù)極大似然估計原理，M步要對E步中所得的Q函數(shù)進行極大化，從而求得參數(shù)θ的極大似然估計，即解決優(yōu)化問題

也即求解如下方程

用矩陣微商的求導公式[19]分別對函數(shù)Q(θ|θ(t))中的參數(shù)μ、Λ、∑和βk求一階偏導,得

分別令上述一階偏導為零，即可得每一個參數(shù)的估計式

2.2 EM算法的收斂準則

對于該模型參數(shù)的極大似然估計，本文在EM算法里采用的停止準則[20]是當估計值θ前后兩次迭代值的絕對差值|θ(t)-θ(t-1)|<10-2時停止迭代，重復進行100次后計算參數(shù)極大似然估計的平均值；在EM算法的E步中，本文用M-H方法產(chǎn)生樣本量為15 000的樣本，然后使用后5 000個樣本來進行統(tǒng)計推斷，以保證所使用的樣本來自目標分布.

2.3 模擬研究

基于所述的參數(shù)估計的方法及原理，下面對該模型中的參數(shù)估計問題進行模擬研究，從下述定義的因子服從指數(shù)分布的因子分析模型中產(chǎn)生模擬數(shù)據(jù).

首先，假設模型中潛在因子隨機向量中包含四個潛變量，即x=(x1,x2,x3,x4)T；可觀測的隨機向量中包含8個可觀測變量，即y=(y1,y2,y3,y4,y5,y6,y7,y8)T，則模型為

y=μ+Λx+ε

待估參數(shù)有

上述矩陣中的常數(shù)表示不進行參數(shù)估計.根據(jù)模型的定義，共有31個待估參數(shù)，參數(shù)取如下真值進行模擬：

分別在樣本量n為200、300和400的情況下進行100次重復估計并取其均值作為參數(shù)估計值(結(jié)果見表1).

從表1的模擬結(jié)果可知，每種樣本量下參數(shù)估計值與真實值之間的偏差都較小，這表明用MCECM算法對該模型的參數(shù)做估計具有良好的效果，即估計值與真實值非常接近.特別地，當樣本量為200時，真實值與估計值之間偏差絕對值的最大值約為0.06，而最小值約為0.001，所以在樣本量為200的情況下，參數(shù)的估計結(jié)果較為準確；另一方面，隨著樣本量的增加，參數(shù)估計值的標準誤差整體在逐漸減小，可見隨著樣本量的增加，參數(shù)的估計結(jié)果越來越準確，雖然其標準誤減小的幅度較小，但可認為這種參數(shù)估計方法對樣本量的要求不高，故MCECM算法對于參數(shù)估計而言具有較強的實用性，針對需要解決的參數(shù)估計問題，樣本量的要求不是很高.

表1 參數(shù)模擬結(jié)果

n代表樣本量；Bias代表真實值與估計值之間的偏差；SD代表標準差.

3 結(jié)語

針對諸多應用領域中相關數(shù)據(jù)為非負的情形，提出了一個傳統(tǒng)因子模型的替代模型；用發(fā)展成熟的MCECM算法[21-22]來對模型中的參數(shù)進行估計.研究結(jié)果表明，用MCECM算法對因子分析模型中的參數(shù)進行估計是一種有效的方法.