劉鶴飛，王 坤，宗鳳喜

（曲靖師范學院 數(shù)學與統(tǒng)計學院，云南 曲靖 655011）

0 引言

混合模型是一種常見的統(tǒng)計學模型，在工業(yè)統(tǒng)計，生物醫(yī)學統(tǒng)計、經(jīng)濟管理、人工智能等領(lǐng)域都有著非常廣泛的應用。張香云和宋紅鳳[1]研究了基于Gibss抽樣的正態(tài)混合模型的參數(shù)估計。蔡鳴晶[2]研究了正態(tài)混合模型的Bayes分析，利用貝葉斯方法對正態(tài)混合模型進行參數(shù)估計。尤芳[3]研究了EM算法、MCEM算法、Gibbs抽樣這三種算法對混合模型參數(shù)估計的效果差異。但是以上的方法都是在混合個數(shù)明確已知的情形下，直接對模型中的未知參數(shù)進行估計。Richardson[4]研究了混合個數(shù)未知時有限正態(tài)混合模型的模型選擇的問題，利用可逆跳MCMC算法計算混合個數(shù)的后驗概率，選出后驗概率最大的那個作為模型的混合個數(shù)。

本文在已有研究的基礎(chǔ)上，將觀測變量從一維擴展到了多維，研究混合個數(shù)未知的多元正態(tài)混合模型。先采用可逆跳MCMC算法對模型中未知的混合個數(shù)進行選擇，然后在混合個數(shù)確定的情形下，用傳統(tǒng)的MCMC算法對混合模型中的未知參數(shù)進行貝葉斯估計。

其中，在利用可逆跳MCMC算法時，需要對模型中所有的待估參數(shù)進行分裂與合并，而對于多維觀測變量來說，還要求協(xié)方差矩陣都必須是正定矩陣，即利用可逆跳MCMC算法分裂與合并協(xié)方差矩陣時，要保證每一個協(xié)方差矩陣Σ都是正定矩陣，這是一維觀測變量模型使用可逆跳MCMC算法時所不曾遇到的問題。

1 模型描述

假設yip是模型的觀測變量，其中i=1，2，3…，N是觀測對象，P是多元正態(tài)分布的維數(shù)。則多元正態(tài)混合模型的定義為：

其中，yi是觀測向量，且y1，y2，…，yN相互獨立，K是模型的混合個數(shù)，表示模型是由K個獨立的多維正態(tài)分布混合而成，θj就是第j個混合部分的正態(tài)分布的參數(shù)，即θj=(μj，Σj)。wj是第j個混合部分的權(quán)重，0

為了研究的方便，本文定義一個隱狀態(tài)Zi，其中：

當隱狀態(tài)Z給定的條件下，觀測向量來自于與其狀態(tài)對應的多元正態(tài)分布，即：

上式就是K個混合部分的多元正態(tài)混合模型的數(shù)學定義。

2 貝葉斯推斷原理

將模型中的觀測變量記為Y，隱狀態(tài)向量集合記為Z，權(quán)重向量記為W，多元正態(tài)分布的參數(shù)（μ,Σ）記為θ，則該模型的貝葉斯推斷問題為相應的MCMC算法的具體執(zhí)行步驟為：

步驟1：更新隱狀態(tài)Z；

步驟2：更新權(quán)重W;

步驟3：更新多維正態(tài)分布的參數(shù)θ;

步驟4：將一個混合成分分裂成兩個混合成分或者將兩個混合成分合并成一個混合成分。

步驟1至步驟3是傳統(tǒng)的MCMC算法，使用Gibbs抽樣和MH算法即可，其中步驟3更新多維正態(tài)分布的參數(shù)θ又可以分解為兩小步，更新多元正態(tài)分布的均值向量μ和更新多元正態(tài)分布的協(xié)方差矩陣Σ，步驟4是可逆跳步，這一步是進行模型選擇的關(guān)鍵，在這一步中，可能增加一個混合成分或者減少一個混合成分，并以一定的概率接受這種改變。

3 先驗分布

貝葉斯學派最基本的一個觀點是，任何一個未知參數(shù)θ都可以看作一個隨機變量，應該用一個概率分布去刻畫θ的未知狀態(tài)，這個概率分布通常被稱為未知參數(shù)的先驗分布[5]。本文參考已有文獻的選取經(jīng)驗，為模型中的各個參數(shù)選擇的先驗分布為：

其中，α是狄尼克雷分布D中的超參數(shù)；μ0j是第j個多維正態(tài)分布的先驗分布均值，Λ0j是第j個多維正態(tài)分布的均值的先驗協(xié)方差矩陣，μ0j和Λ0j都是超參數(shù)；Wishartp表示P維的威希爾分布，R0j和ρ0j是威希爾分布中的兩個超參數(shù)。在本文的研究中，混合個數(shù)也是未知的，在貝葉斯推斷中也應該看成是一個未知參數(shù)，在這里選擇離散的均勻分布作為該參數(shù)的先驗分布。

4 可逆跳躍MCMC過程

1995年，Green[6]提出了可逆跳躍MCMC方法，該方法通過在可變維參數(shù)空間中跳躍式抽樣，實現(xiàn)模型選擇的目的。1997年，Richardson和Green[4]正是利用這種可逆跳MCMC方法，對于正態(tài)混合模型中未知的混合個數(shù)進行選擇。本文就是在此基礎(chǔ)上，將一元正態(tài)混合模型擴展到多元正態(tài)混合模型，以增加一個混合部分為例，該算法的具體執(zhí)行過程描述如下：

（1）從當前的K個混合部分中等概率隨機的挑選一個；

（2）從事先給定的分布中抽取一個隨機數(shù)向量；

（3）利用所抽取的隨機數(shù)向量，將選中的混合部分的參數(shù)分裂成兩個混合部分；

（4）以概率min(1 ，A)接受這個跳躍，其中：

A=似然比×跳躍比×先驗比×建議比×雅克比

假設當前的混合個數(shù)為K，則對于P維觀測變量的多維正態(tài)混合模型來說，此時的權(quán)重向量是K維的，并且有K個P維的均值向量μ，有K個P維正定的協(xié)方差矩陣Σ。根據(jù)可逆跳MCMC算法的原理，參數(shù)分裂之后，K維的權(quán)重向量變成K+1維的，且增加一個P維均值向量μ和一個P維的協(xié)方差矩陣Σ，并且協(xié)方差矩陣必須是正定矩陣。

將選中的待分裂的混合部分標記為j*，將j*分裂后的兩個混合部分分別標記為j1，j2。則生成隨機數(shù)分裂權(quán)重W的方式為：

生成隨機數(shù)分裂均值向量μ的方式為:

生成隨機數(shù)分裂協(xié)方差矩陣Σ的方式非常特殊，根據(jù)可逆跳MCMC算法參數(shù)分裂原理的要求，需要保證分裂和合并的過程是可逆的，即分裂和合并的過程要一一對應而且能完全覆蓋整個參數(shù)空間，另外，還要保證協(xié)方差矩陣Σ是正定矩陣這個特殊的要求。所以，協(xié)方差矩陣Σ的分裂和合并也是將一元正態(tài)混合模型的可逆跳MCMC擴展到多維正態(tài)混合模型的可逆跳MCMC時的困難之處。本文根據(jù)正定矩陣的Cholesky分裂原理，提出了一種分裂和合并的方法，符合以上所有的條件，。

因為正定矩陣都是對稱矩陣，所以P維協(xié)方差矩陣Σ中有個自由元素。根據(jù)正定矩陣的分裂原理，從標準正態(tài)分布中生成/2個隨機數(shù)，將這些隨機數(shù)按次序排列成一個標準的下三角矩陣L，對協(xié)方差矩陣Φj*進行分裂的時候，令Φj1=Φj*，Φj2=L?LT即可。對協(xié)方差矩陣Φj1和Φj2進行合并的時候，令Φj*=Φj1，然后對Φj2進行Cholesky分裂，就可以得到相應的隨機數(shù)，這個分裂和合并的方法是可逆的，而且理論上能取到參數(shù)空間的所有值，更重要的是可以保證分裂和合并的過程中每一個協(xié)方差矩陣都是正定矩陣。

5 后驗推斷

對本文來說，貝葉斯后驗即P(Z,W,θ|Y)，尤其是對模型中的未知參數(shù)W和θ進行估計，其中θ包括均值向量μ和協(xié)方差矩陣Σ。由于這個后驗分布非常的復雜，本文采用MCMC方法來對后驗分布進行后驗模擬，下面來推導后驗模擬時需要用到的滿條件后驗分布：

因為：

f是在第j個混合部分的參數(shù)(μj,Σj) 下，觀測向量yi的似然函數(shù)。

由于各混合部分的權(quán)重只與觀測變量的隱狀態(tài)的取值有關(guān)，所以：

又因為權(quán)重向量W的共軛先驗分布是對稱的狄尼克雷分布所以權(quán)重向量W的滿條件后驗分布[7]為：

其中，nj是隱狀態(tài)為j的觀測變量的個數(shù)。

進一步，推導多元正態(tài)分布均值向量μ的滿條件后驗分布P(μ|Y,Z,W,Σ)，因為一旦所有觀測變量的隱狀態(tài)確定的話，權(quán)重向量W即可隨之確定。此時第j個多維正態(tài)分布的均值向量μj只與觀測變量Yj及協(xié)方差矩陣Σj有關(guān)，也即其中，Yj是隱狀態(tài)取值為j的觀測變量的集合，Σj是第j個多維正態(tài)分布的協(xié)方差矩陣。

因為多維正態(tài)分布N(μ,Σ)的聯(lián)合概率密度函數(shù)為：

由于μj的共軛先驗分布為N(μ0j,Λ0j) ，其中μ0j和Λ0j是超參數(shù)，都可以看成是已知的，所以μj的后驗密度的形式為：

忽略上式中的常數(shù)因子，將其簡化為：

這表明，μ的后驗密度仍然是P維正態(tài)分布，所以為均值向量μ選擇共軛多元正態(tài)分布是合理的，且μk后驗均值向量協(xié)μjp協(xié)方差矩陣Λjp分別為：

也即，μj的滿條件后驗分布[8]為：

其中，nj是觀測變量中隱狀態(tài)為j的總個數(shù)是隱狀態(tài)為j的所有觀測向量的均值向量。

最后來推導多維正態(tài)分布N(μj,Σj) 的協(xié)方差陣∑j的滿條件后驗分布[9]。

∑j的先驗分布是一個威希爾分布Wishartp(R0j,ρ0j) ，其概率密度函數(shù)為：

記Yj為隱狀態(tài)為j的觀測變量的集合，則Σj的滿條件后驗分布的密度函數(shù)有以下形式：

上式表明，Σj的滿條件后驗分布為：

其中，nj是隱狀態(tài)為j的觀測變量的總個數(shù)，IWp表示P階的倒威希爾分布。

上面推導出了所有未知參數(shù)的滿條件后驗分布，但是在混合模型的研究中，還有一個非常重要的問題就是標簽轉(zhuǎn)換。因為如果不處理好標簽轉(zhuǎn)換問題，在迭代的過程中，不同的混合部分就有可能發(fā)生混亂，肯定會影響估計的結(jié)果。交換抽樣技術(shù)[10]是現(xiàn)在最常見的解決混合模型標簽轉(zhuǎn)換的方法之一，但是也有很多文獻是通過對模型中某個特定的參數(shù)進行控制來解決標簽轉(zhuǎn)換這個問題的。比如，Richardson和Green[4]通過限制均值μ1<μ2<…<μn來解決就是一個非常好的例子。本文為了防止不同混合部分下的多維正態(tài)分布的參數(shù)在迭代的過程中發(fā)生混亂，采用限制均值向量的第一維分量的大小的方法來解決這一問題。在每一次迭代更新參數(shù)之后，比較μ1,1,μ2,1,…,μk,1的值，如果μ1,1＜μ2,1＜…＜μk,1，則進行下一輪的迭代更新，否則按照不同混合部分均值向量的第一維分量從小到大的順序重新排序，并且重新分配隱狀態(tài)的標簽值，然后再進行下一輪的迭代。通過這種方法就可以保證在迭代的過程中不同隱狀態(tài)下的參數(shù)不會發(fā)生混亂。

6 數(shù)值模擬

在真實的模型中，本文取混合個數(shù)K=2，觀測對象的個數(shù)N=500，多元正態(tài)分布的維數(shù)P=2，兩個多元正態(tài)分布的權(quán)重分別為0.5和0.5，兩個多元正態(tài)分布N(μ,Σ)的均值向量和協(xié)方差矩陣分別為

取狄尼克來分布中的超參數(shù)α=1，多元正態(tài)分布均值向量中的超參數(shù)為2階單位矩陣，j=1,2;威希爾分布中的超參數(shù)ρ0j=6,j=1,2;R0j=I，I為2階單位矩陣，j=1,2。

利用可逆跳MCMC算法選擇模型中的未知混合個數(shù)時，迭代的總次數(shù)為5萬次，刪去前面的3萬次，用后面的2萬次的迭代結(jié)果來計算未知混合個數(shù)K的后驗概率，圖1是混合個數(shù)的迭代路徑圖，表1是混合個數(shù)的后驗結(jié)果。

圖1 隱狀態(tài)個數(shù)的迭代路徑圖

表1 隱狀態(tài)個數(shù)的后驗結(jié)果

根據(jù)表1混合個數(shù)K的后驗概率，本文選擇模型的混合個數(shù)為2，接下來再用馬爾科夫鏈蒙特卡羅算法對混合個數(shù)確定時模型中的所有參數(shù)進行貝葉斯估計，超參數(shù)的取值不變，取MCMC算法的迭代次數(shù)為1萬，去掉前面的5千次，用后面5千次迭代的后驗均值作為參數(shù)的估計值，估計結(jié)果見表2。

表2 各參數(shù)的MCMC迭代真值和估計值結(jié)果

7 結(jié)論

本文研究了混合個數(shù)未知的多元正態(tài)混合模型后的貝葉斯推斷，利用可逆跳MCMC算法對未知的混合個數(shù)進行了貝葉斯模型選擇，給出了算法的具體執(zhí)行步驟，在分裂和合并多元正態(tài)分布的協(xié)方差矩陣Σ時，除了要滿足可逆跳躍算法的要求，還必須保證分裂和合并過程中每一個協(xié)方差矩陣都必須是正定矩陣，本文根據(jù)正定矩陣的Cholesky分裂原理，提出了一種符合以上所有要求的協(xié)方差矩陣的分裂和合并的方法，這種方法可以推廣到其他多維分布協(xié)方差矩陣的分裂和合并中去。

確定出模型真實的混合個數(shù)以后，即當混合個數(shù)已經(jīng)確定時，利用MCMC算法對模型的其他未知參數(shù)進行了貝葉斯估計，將表2中的各未知參數(shù)的估計值與事先確定的模型的真實值進行對比，顯示該方法的估計效果比較好。