1、西南交通人學(xué)，信息學(xué)院，xgongsvHucn資料整理高斯混合模型(GMM)參數(shù)優(yōu)化及實(shí)現(xiàn)龔勛(2010-11-13)1高斯混合模型概述高斯密度函數(shù)估計(jì)是一種參數(shù)化模型。有單高斯模型(SingleGaussianModelSGM)和高斯混合模型(Gaussianmixturemodel,GMM)兩類。類似于聚類,根據(jù)高斯概率密度函數(shù)(PDF,見(jiàn)公式1)參數(shù)的不同，每一個(gè)高斯模型可以看作一種類別，輸入一個(gè)樣本x,即可通過(guò)PDF計(jì)算其值，然后通過(guò)一個(gè)閾值來(lái)判斷該樣本是否屬于高斯模型。很明顯，SGM適合于僅有兩類別問(wèn)題的劃分，而GMM由于具有多個(gè)模型，劃分更為精細(xì)，適用于多類別的劃分，可以應(yīng)用于復(fù)雜

2、對(duì)彖建模。下面以視頻前景分割應(yīng)用場(chǎng)景為例，說(shuō)明SGM與GMM在應(yīng)用上的優(yōu)劣比較：ISGM需耍進(jìn)行初始化，如在進(jìn)行視頻背景分割時(shí)，這意味著如果人體在前幾幀就出現(xiàn)在攝像頭前，人體將會(huì)被初始化為背景，而使模型無(wú)法使用；ISGM只能進(jìn)行微小性漸變，而不可突變。如戶外亮度隨時(shí)間的漸變是可以適應(yīng)的，如果在明亮的室內(nèi)突然關(guān)燈，單高斯模型就會(huì)將整個(gè)室內(nèi)全部判斷為前景。又如，若在監(jiān)控范圍內(nèi)開了一輛車，并在攝像頭卜開始停留。由于與模型無(wú)法匹配，車會(huì)一直被視為前景。當(dāng)車過(guò)很長(zhǎng)時(shí)間離去時(shí)，由于車停留點(diǎn)的亮度發(fā)生了很人的變化，因此已經(jīng)無(wú)法與先前的背景模型相匹配；ISGM無(wú)法適應(yīng)背景冇多個(gè)狀態(tài)，如窗簾，風(fēng)吹的樹葉。單高斯

3、模型無(wú)法表示這種情況，而使得前背景檢測(cè)混亂，而GMM能夠很好地描述不同狀態(tài)；I相對(duì)于單高斯模型的自適應(yīng)變化，混合高斯模型的自適應(yīng)變化要健壯的多。它能解決單高斯模型很多不能解決的問(wèn)題。如無(wú)法解決同一樣本點(diǎn)的多種狀態(tài)，無(wú)法進(jìn)行模型狀態(tài)轉(zhuǎn)化等。1-1單高斯模型(1)“是模型期望,是模型方差。多維高斯(正態(tài))分布概率密度函數(shù)PDF定義如卜一:注意與一維高斯分布不同，其中X是維數(shù)為d的樣本向杲(列向最),對(duì)丁單高斯模型，由于可以明確訓(xùn)練樣本是否屬于該高斯模型(如訓(xùn)練人臉膚色模型時(shí)，將人臉圖像膚色部分分割出來(lái)，形成訓(xùn)練集)，故“通常由訓(xùn)練樣本均值代替，由樣本方差代替。為了將高斯分布用于模式分類，假設(shè)訓(xùn)練樣

4、本屬于類別C,那么，式(1)可以改為如卜形式：N(x/C)=-j=exp_dx_“)T7(x_“)(2)J(2利斗2厶式表明樣本屬丁類別C的概率人小。從而將任意測(cè)試樣本旺輸入式(2),均可以得到一個(gè)標(biāo)量N(兀:“Q,然后根據(jù)閾值/來(lái)確定該樣本是否屬于該類別。2閾值/的確定：可以為經(jīng)驗(yàn)值，也可以通過(guò)實(shí)驗(yàn)確定。另外也有一些策略可以參考：如文獻(xiàn)1中采用搜索策略：首先令/=0.7,以0.05為步長(zhǎng)一直減到0.1左右，選擇使樣本變化最小的那個(gè)閾值做為最終f值，也就是意味著所選/值所構(gòu)造的分類模型最穩(wěn)定。2幾何意義理解：根據(jù)單高斯分布PDF的含義我們可以知道，符合SGM分布的二維點(diǎn)在平面上應(yīng)該近似橢圓形；

5、相應(yīng)地，三維點(diǎn)在空間中則近似于橢球狀酋南交通人學(xué)，信息學(xué)院，xgongsvHucn資料整理酋南交通人學(xué)，信息學(xué)院，xgongsvHucn資料整理(5) #(5) 1.2高斯混合模型高斯混合模型是單一高斯機(jī)率密度函數(shù)的延伸，由J*GMM能夠平滑地近似任意形狀的密度分布，因此近年來(lái)常被用在語(yǔ)音、圖像識(shí)別等方面，得到不錯(cuò)的效果。例：有一批觀察數(shù)據(jù)X=x.，“,數(shù)據(jù)個(gè)數(shù)為，在d維空間中的分布不是橢球狀（如圖1（a）,那么就不適合以一個(gè)單一的高斯密度函數(shù)來(lái)描述這些數(shù)據(jù)點(diǎn)的機(jī)率密度函數(shù)。此時(shí)我們采用一個(gè)變通方案，假設(shè)每個(gè)點(diǎn)均由一個(gè)單高斯分布生成（如圖1（b）,具體參數(shù)“廠刁未知），而這一批數(shù)據(jù)共由M（明確

6、）個(gè)單高斯模型生成，具體某個(gè)數(shù)據(jù)旺屬丁哪個(gè)單高斯模型未知，且每個(gè)單高斯模型在混合模型中占的比例未知，將所有來(lái)自不同分布的數(shù)據(jù)點(diǎn)混在一起，該分布稱為高斯混合分布。()0酋南交通人學(xué)，信息學(xué)院，xgongsvHucn資料整理酋南交通人學(xué)，信息學(xué)院，xgongsvHucn資料整理(5) #(5) #酋南交通人學(xué)，信息學(xué)院，xgongsvHucn資料整理酋南交通人學(xué)，信息學(xué)院，xgongsvHucn資料整理(5) #(5) #00.51（a）（b）圖1高斯混合模型圖示，（a）表示所有樣本數(shù)據(jù)：（b）表示已經(jīng)明確了樣本的分類從數(shù)學(xué)上講，我們認(rèn)為這些數(shù)據(jù)的概率分布密度函數(shù)可以通過(guò)加權(quán)函數(shù)表示：冃上式即稱為

7、GMM,fa產(chǎn)1，其中Nj（x;“（爲(wèi)）=.1=exp-l（x-“丿丿（4）J（2兀）書|2表示第j個(gè)SGM的PDF。令巧=（“,為），GMM共何M個(gè)SGM模型，現(xiàn)在，我們就需要通過(guò)樣本集X來(lái)估計(jì)GMM的所有參數(shù)：0=（%你）丁。樣本X的概率公式為：/（x|0）=iogHEaA（.）1=1;=1NM=1。遼巧y（x；“尹為）i=i；=1西南交通人學(xué)，信息學(xué)院，資料整理西南交通人學(xué)，信息學(xué)院，xgonMsvHucn資料整理 #西南交通人學(xué)，信息學(xué)院，資料整理西南交通人學(xué)，信息學(xué)院，xgonMsvHucn資料整理 2釆用EM估計(jì)GMM的參數(shù)通常采用EM算法(期塑值最人，ExpectationMax

8、imum)根據(jù)式對(duì)GMM參數(shù)進(jìn)行估計(jì)。具體推導(dǎo)過(guò)程在很多文獻(xiàn)可以參閱，這里不詳述，我們主要介紹采用EM估計(jì)GMM參數(shù)的方法。描述估計(jì)高斯混合分布的EM算法(參考“!EM算法.ppt”)：2初始值：方案1：協(xié)方差矩陣與。設(shè)為單位矩陣；每個(gè)模型比例的先驗(yàn)概率：均值的。設(shè)為隨機(jī)數(shù)；方案2：由k均值(k-mean)聚類算法對(duì)樣本進(jìn)行聚類，利用各類的均值作為的并計(jì)算丹，取各類樣本占樣本總數(shù)的比例。2算法流程：估計(jì)步驟(E-step)：令勺的后驗(yàn)概率為：aN(x0)卩產(chǎn)E(aK；0)=,l/n,jM(6)fI注意：寫代碼實(shí)現(xiàn)公式(6)時(shí)，對(duì)丁每個(gè)SGM分別用公式計(jì)算每個(gè)樣本點(diǎn)兀在該模型卜的概率密度值對(duì)于所

9、苗樣本，得到一個(gè)的向量，計(jì)算M次，得到八M的矩陣，每一列為所有點(diǎn)在該模型卜的概率密度值(PDF)；實(shí)現(xiàn)時(shí)，需要針對(duì)每個(gè)點(diǎn)計(jì)算在各個(gè)SGM的概率值總和。公I(xiàn)式(6)可以通過(guò)如Fmatlab代碼描述：functionprob=GaussPDF(Data,Mu,Sigma)dim,N=size(Data);Data=Data*-repmat(MuNJ);prob=sum(Data*inv(Sigma).*Data,2);prob=exp(-0.5*prob)/sqrt(2*pi)Adim*(abs(det(Sigma)+realmin);最人化步驟(M-step)：(7)0更新權(quán)值：6=幾衛(wèi)N0更新

10、均值：叮封產(chǎn)0更新方差矩陣：1嚴(yán)星也筆巴竺注意，有些文獻(xiàn)把方差矩陣當(dāng)作對(duì)角陣處理，即認(rèn)為2；,=diag52),故采用如卜公式估計(jì)L：(10)經(jīng)本人實(shí)驗(yàn)驗(yàn)證，公式(10)會(huì)給最終結(jié)果帶來(lái)誤差，故在運(yùn)算速度能保證的情況卜建議采用公式(9)估計(jì)方差矩陣。收斂條件方案1：不斷地迭代E和M步驟，重復(fù)更新上而的三個(gè)值，直到|/(X|0)-/(X|0)je,通常E=10-5,/(X|0)通過(guò)公式(5)計(jì)算，/(X|0j表示更新參數(shù)后計(jì)算的值；方案2：不斷地迭代E和M步驟，重復(fù)更新上面的三個(gè)值，直到參數(shù)的變化不顯著，即|0-0|YCrCbIplImage*imgYCrCb=cvCreateIniage(cv

11、GetSize(pImg),IPL_DEPTH_8U,3)y/YCrCb圖像cvCvtColor(pImg,imgYCrCb,CV_BGR2YCrCb);第0,1,2層分別為YCr.CbIpIImage*imgY=cvCreateImage(cvGetSize(pImg)JPL_DEPTH_8UJ);/YCrCb圖像IpIImage*imgCr=cvCreateImage(cvGetSize(pImg)JPL_DEPTH_8UJ)y/YCrCb圖像IpIImage*imgCb=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,1);/YCrCb圖像IpIIm

12、age*imgY32=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32E1)/YCrCb圖像IpIImage*imgCr32=cvCreateImage(cvGetSize(pImg)JPL_DEPTH_32E1);/YCrCb圖像IpIImage*imgCb32=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32F,l);/YCrCb圖像cvSplit(imgYCrCb,imgY,imgCr,imgCb,NULL);cvConvert(imgYimgY32);cvConvert(imgCr,imgCr32);cvConve

13、rt(imgCb,imgCb32);llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllII2根據(jù)SigleGaussianModel計(jì)算顏色模型IpIImage*PCbCr=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32F,l)YCrCb顏色模型IpIImage*tempA=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32El)YCrCb顏色模型IpIImage*tempB=cvCreateImage(cvGetSize(p

14、Img),IPL_DEPTH_32El)y/YCrCb顏色模型cvSubS(imgCb32,cvScalar(M0),imgCb32);/x-mcvSubS(imgCr32,cvScalar(Ml),imgCr32);x-mcvAddWeighted(imgCb32,invC00,imgCr32,invCl0,0,tempA);cvAddWeighted(imgCb32,invC0l,imgCr32,invCll,0,tempB);cvMul(imgCb32,tempA,tempA,-0.5);cvMul(imgCr32,tempB.tempB,-0.5);cvAdd(tempA,tempB,

15、PCbCr);cvExp(PCbCr,PCbCr);doublemax_val=0,min_val=0;cvMinMaxLoc(PCbCr,&min_val,&max_val);IpIImage*proImg=cvCreateImage(cvGetSize(pImg)JPL_DEPTH_8U,1);/YCrCb顏色模型doublea=255/(max_val);cvConvertScaleAbs(PCbCnproImg,a,0);m_proimg=cvCIonelmage(proImg);if(*skinImg)!=NULL)cvReleaselmage(skinlmg);*skinlmg=c

16、vCreateImage(cvGetSize(pImg)JPL_DEPTH_8U,1);/膚色結(jié)呆/釋放內(nèi)存cvReIeaseImage(&prolmg);cvReleaseImage(&imgYCrCb);cvReIeaseImage(&imgY);cvReIeaseImage(&imgC;cvReIeaseImage(&imgCb);cvReIeaseImage(&imgY32);cvReIeaseImage(&imgCr32);cvReIeaseImage(&imgCb32);cvReleaseImage(&PCbCr);西南交通人學(xué)，信息學(xué)院，xgonMsvHucn資料整理西南交通人學(xué)

17、，信息學(xué)院，資料整理 cvReleaseImage(&tempA);cvReleaseImage(&tempB);3.1高斯混合模型以下matlab代碼實(shí)現(xiàn)了高斯混合模型：functionAlpha,Mu,Sigma=GMM_EM(Data,AlphaO,MuO,SigmaO)%EM迭代停止條件loglik.threshold=le-IO;%初始化參數(shù)dim,N=size(Data);M=size(MuO,2);loglik.old=-realmax;nbStep=0;Mu=MuO;Sigma=SigmaO;Alpha=AlphaO;Epsilon=0.0001;while(nbStep120

18、0)nbStep=nbStep+1;%E-步驟%fori=l:M%PDFofeachpointPxi(:,i)=GaussPDF(Data,Mu(:,i),Sigma(:,:,i);end%計(jì)算后驗(yàn)概率beta(ilx)Pix_tmp=repmat(Alpha,Nl).*Pxi;Pix=Pix_tmp./(repmat(sum(Pix_tmp,2),1M)+realmin);Beta=sum(Pix);%M-步驟%fori=l:M%更新權(quán)值A(chǔ)lpha(i)=Beta(i)/N;%更新均值Mu(:,i)=Data*Pix(:,i)/Beta(i);%更新方差Data_tmpl=Datarepma

19、t(Mu(:,i),1,N);Sigma(:;i)=(repmat(Pix(:,i)dim,1).*Data_tmpl*Data_tmpl)/Beta(i);%AddatinyvariancetoavoidnumericalinstabilitySigma(:,：a)=Sigma(:；i)+1E-5.*diag(ones(dim,1);end%Stoppingcriterion1%fori=l:M%Computethenewprobabilityp(xli)%Pxi(:j)=GaussPDF(Data,Mu(:,i),Sigma)；%end%Computetheloglikelihood%F=

20、Pxi*AIpha;%F(find(Frealmin)=realmin;%loglik=mean(log(F);%Stoptheprocessdependingontheincreaseoftheloglikelihood%ifabs(loglik/loglik_old)-l)loglikjhreshold%break;%end%loglik_old=loglik;%Stoppingcriterion2%v=sum(abs(MuMuO),abs(AlphaAlphaO);s=abs(SigmaSigmaO);v2=0;fori=l:Mv2=v2+det(s(:,:,i);endif(sum(v

21、)+v2)Epsilon)break;endMuO=Mu;SigmaO=Sigma;AlphaO=Alpha;endnbStep以卜代碼根據(jù)高斯分布函數(shù)計(jì)算每組數(shù)據(jù)的概率密度，被GMM.EM函數(shù)所調(diào)用functionprob=GaussPDF(Data,Mu,Sigma)%根據(jù)高斯分布函數(shù)計(jì)算每組數(shù)據(jù)的概率密度ProbabilityDensityFunction(PDF)%輸入%oData:DxN,N個(gè)D維數(shù)據(jù)%oMu:DxI,M個(gè)Gauss模型的中心初始值%oSigma:MxM，每個(gè)Gauss模型的方差(假設(shè)每個(gè)方差矩陣都是對(duì)角陣，%即一個(gè)數(shù)和單位矩陣的乘積)%Outputs%oprob:1

22、xNarrayrepresentingtheprobabilitiesforthe%Ndatapointsdim,N=size(Data);Data=Data*-repmat(MuNJ);prob=sum(Data*inv(Sigma).*Data,2);prob=exp(-0.5*prob)/sqrt(2*pi)Adim*(abs(det(Sigma)+realmin);以卜是演示代碼demol.m%高斯混合模型參數(shù)估計(jì)示例(基于EM算法)%2010年11月9日data,mu,vanweight=CreateSample(M,dim,N);/生成測(cè)試數(shù)據(jù)Alpha,Mu,Sigma=GMM_EM(Data,Priors,Mu,Sigma)以下是測(cè)試數(shù)據(jù)生成函數(shù)，為demol.m所調(diào)用：fun