混合個(gè)數(shù)未知的多元正態(tài)混合模型的貝葉斯推斷

劉鶴飛，王 坤，宗鳳喜

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

0 引言

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

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

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

1 模型描述

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

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

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

當(dāng)隱狀態(tài)Z給定的條件下，觀測(cè)向量來自于與其狀態(tài)對(duì)應(yīng)的多元正態(tài)分布，即：

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

2 貝葉斯推斷原理

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

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

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

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

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

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

3 先驗(yàn)分布

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

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

4 可逆跳躍MCMC過程

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

（1）從當(dāng)前的K個(gè)混合部分中等概率隨機(jī)的挑選一個(gè)；

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

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

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

A=似然比×跳躍比×先驗(yàn)比×建議比×雅克比

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

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

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

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

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

5 后驗(yàn)推斷

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

因?yàn)椋?/p>

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

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

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

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

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

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

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

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

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

也即，μj的滿條件后驗(yàn)分布[8]為：

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

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

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

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

上式表明，Σj的滿條件后驗(yàn)分布為：

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

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

6 數(shù)值模擬

在真實(shí)的模型中，本文取混合個(gè)數(shù)K=2，觀測(cè)對(duì)象的個(gè)數(shù)N=500，多元正態(tài)分布的維數(shù)P=2，兩個(gè)多元正態(tài)分布的權(quán)重分別為0.5和0.5，兩個(gè)多元正態(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算法選擇模型中的未知混合個(gè)數(shù)時(shí)，迭代的總次數(shù)為5萬次，刪去前面的3萬次，用后面的2萬次的迭代結(jié)果來計(jì)算未知混合個(gè)數(shù)K的后驗(yàn)概率，圖1是混合個(gè)數(shù)的迭代路徑圖，表1是混合個(gè)數(shù)的后驗(yàn)結(jié)果。

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

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

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

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

7 結(jié)論

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

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