1、EM算法在高斯混合模型中的應用1 .定義對于一個隨機信號生成器，假設他的模型參數為，我們能觀測到的數據輸出為X,不能觀測到白數據輸出為 Y,且隨機系統模型結構的概率密度函數為p(x,y| )(1)能夠觀測到的一部分數據輸出數據Xi,X2,.,Xn,模型的另一部分輸出數據 未知，模型的參數也未知。EMJ法就是要求我們從觀測數據Xi,X2,.,Xn中估計出參數 。2 .EM算法的描述假設每一對隨機系統的輸出樣本(xn, yn)對于不同的n相互獨立，這樣當 p(x,y, ) , x和y都已知的情況下，概率p(x,y,)也已知。未觀測的輸出y的概率 分布也屬于待求參數。根據獨立性假設有：Np(x, y

2、| )P(Xn,yn| )(2)n 13 .EM算法的基本思路基本問題是求解下面的方程的解：argmax p(x,y | )(3)由于X®確定量，Y是未知的，因此即使給定了，也無法求得p(x,y| )的值，因此我們只能退一步求：argmax p(x| )(4)其中p(x| )y Y p(x,y| ) y Yp(y| ), p(x|y, )(5)表示考慮了未知數據y的所有可能的取值Y后又tp(x|y,)求平均值。最后根據10g函數的單調性得到(4)的等效形式：argmaxlog p(x | )(6)對于(6)給出的最優化問題，考慮用下面的遞推算法解決，即：先給定一個估值k 并計算p(x

3、| k),然后更新k得到k1并且有log p( x | k 1) log p(x| k)log p(x| ) log yYp(y| )p(x| y,)logy Y p(y |x,k p y | p(x| y,)/ I kp(y|x,)(8)y Y p(y |x, k)logB( , k)其中，等號在B( k, k)時成立，即：B( k, k) logp(x|p(y| )p(x|y,)p(y|x, k)k)(9)于是對10gp(x| )的遞推算法(7)可通過B( , k)進行，步驟為：1)令k=0,先給出估值k2)然后找出 k1 滿足 B( k1, k) B( k, k)(10)3) k更新為k

4、+1并返回步驟2)直到收斂令k 1 arg max B( , k)(11)處理后k1 argmaxB( , k)argmax 丫丫 p(y|x,k)1ogp(y| )P(x|y,)p(y|x, k)p(y|x, k)10gp(y|x, k) (12)kargmax 丫丫 P(y|x, )1og p(y| )p(x| y,)kargmax 、 p(y|x, )1og p(y| )p(x| y,)argmaxC( , k)其中C( , k) y y p(y |x, k)1og p(y| ) p(x|y, )(13)4 .EMB法與高斯混合模型在隨機系統模型中，假設m是通道m的隨機信號生成器的概率密

5、度函數的參數，p(y m)是選中通道m的概率。記為am。(14)(15)假設M個隨機信號生成器和通道選擇隨機生成器是相互獨立的，從通道 m 輸出的數據x的概率是：am Pm(xI m) 不考慮通信信息，輸出x的概率為： Mp(x| )ampm(x| m)m 1其中：m :是第m個通道隨機信號生成器的參數。:參數集合m,am m12 M。,m., ' I觀測數據為一批隨機產生的輸出信號,并且每個輸出都是相互獨立的，而每個輸出來自哪個通道不可測。于是系統模型參數估計問題就變為通過有限的輸出 樣本xi,X2,., Xn估計 M 個通道參數m,am (m 1,2,., M ).應用(12)求解

6、，其中C( , k)可以簡化為：M NM NC( , k)log(am)p(m|xn, k)m 1 n 1Ci( , k) C2( , k)k Xlog(pm(Xm| m)p(m|Xn,)m 1 n 1(16)其中:Ci( , k)C2( , k)k、10g(am)p(m|Xn,)miniM N k、10g( Pm(Xm | m)P(m|Xn,)mini這樣我們把am和Pm分別放在兩項里面，他們不相關，可以獨立考慮。 在C( , k)中應用約束條件：Mam 1m 1用拉格朗日乘子優化am得到:1 NuNn1P(m|Xn，)上式的含義是，選中m號通道的概率估計am1是每個觀測數據來自于m通道的條

7、件概率(根據上一次估值 得出。估算)的平均。其中的p(m|Xn，k)通過下式p(m|4, k) MamPm(Xn| I amPm'(Xn| Tm 1C2( , k)中的m的優化取決于分布函數的類型，對于Pm(Xm | m)為高斯分布時，m m, m其中m是分布的均值，m是方差。再經過推導，有： k 11 NkamP(m |Xn,)N n 1 NXnP(m|Xn, k) k 1_n_Jm-Nk .1/2|2P(m| Xn,) n 1 N P(m| Xn, k) |Xn k 1_njm =NP(m|Xn, k) n 1 k m通道參數m, m得更新可以看作是對Xn的加權，加權系數P(m|

8、Xm,)可以看成是根據上一次的參數估計k算出來得Xn率屬于m通道的概率。最后，上面的EMT法可能收斂到局部極大點，因此需要選擇多個參數 的初 始值進行迭代計算，并選擇使得P(x| )最大的解，p(x| )最大的解可由下式算 出：p(x| )p(y | )p(x | y, )yYMM MNLm1 1 m2 1mN 1 n 1amn p(xn | mm )Nn1am p(xn | m1m)5 .EM算法在matlab中的實現利用上面推導出的公式，我們以二個一維的高斯分布(Ni, N2)來驗證EMJ法的性能,首先用二個一維的高斯分布來建立一個高斯混合模型Nh假設：22 N1 : N( 1, 12)

9、, N2: N( 2, 22) NH1N12N2其中i與2為混合系數，且有i 2 1 ,我們要用EMT法估計混合系數和 各一維高斯分布的均值和方差(1, 2, 1, 2, 12, 22)。并將利用EMT法估計出的 值與真實值做比較，就可以得到該算法的性能。首先我們取( 1, 2, 1, 2, 12, 22) 的真實值為( 0.4,0.6,1,2,0.25,0.36)這樣我們得到一個混合高斯分布， 他的密度函數為 N H , 然后產生1000個服從Nh的分布的觀測樣本點。接下來要做的就是對這1000t樣本點用EMT法進行處( 1, 2, 1, 2, 12, 22) 的值。在使用EMJ法時，要首先

10、給(1, 2, 1, 2, 12, 22)設定一組初值這里假設初值為 1 =0.3 ，2 =0.7 ， 1 0.8 ， 2 1.8 ，12 0.2 ，220.25Matlab 具體程序如下 ：Y=zeros(1,10000);for i=1:10000if rand(1)>0.3Y(i)=normrnd(2,sqrt(0.36),1,1); elseY(i)=normrnd(1,sqrt(0.25),1,1);%高斯混合模型%設置參數的初值%設置均值 的初值%設置方差2 的初值endendA=0.3 0.7;M=0.8 1.8;S=0.2 0.25;for n=1:1000for j=1

11、:2a3=0;a4=0;a5=0;for k=1:10000a1=0;for t=1:2a1=A(t)*1 /sqrt(2*pi*S(t)*exp(-(Y(k)-M(t)A2 /(2*S(t)+a1;endf= A(j) * 1/sqrt(2*pi*S0)*exp(-(Y(k)-M0)F2 /(2*S(j);a2=f/a1;a3=a2+a3;a4=a2*Y(k)+a4;a5=a2*(Y(k)-M(j)A2+a5;endA(j)=a3/10000;M(j)=a4/a3;S(j)=a5/a3;endendNk%a3對應公式p(m | xn,)n1Nk%a4對應公式p(m| xn,)xnn1Nkk12%a5對應公式p(m|xn,)(xn m )n1%循環更新系數值%循環更新均值值%循環更新方差值運行程序，查看變量A,M,S