高斯混合模型之代码实现

高斯混合模型的代码实现，总体的思路是比较简单的。但涉及到具体的优化，如多维高斯概率分布协方差矩阵的逆矩阵，就是一个很头疼的奇异矩阵问题。这里我只想讲下代码实现的流程。具体的代码可以参照：http://blog.csdn.net/crzy_sparrow/article/details/7413019。（注意他的代码没有考虑协方差逆矩阵的问题）。

高斯混合模型代码实现流程：

（1）·首先是初始化，高斯混合模型的效果很大程度上依赖于初始点的设定。一般我们用K-means聚类生成K个中心节点。对于属于同一节点的数据，我们求其均值，方差以及该节点的概率。这里所谓的均值就是中心节点，协方差矩阵按照定义求解，该节点概率（选择该个高斯模型的概率）= 属于该节点的数据个数 / 总数据个数，这样初始化完成。

（2）·E-STEP:求得Q(j)，这里要将上次得到的均值u，协方差sigma，模型概率pj，带入Q(j)的定义式（见“高斯混合模型之理解”），注意p(x|j)是j类高斯概率分布；

（3）·M-STEP:按照我们推导的公式，更新均值u，协方差sigma和模型概率pj；

（4）·将（3）中更新的参数带入（2）中更新Q(j)；

（5）·最后要设定阀值，使迭代结束。按照定义，我们要将u,sigma,pj，带入L(theta)（最大似然值）公式中，如果t+1时刻的L与t时刻的L的比值接近于1，即可停止。具体的阀值还要应对实际的数据进行调整；

我的代码（MATLAB）：

·初始化函数：

function [ mu,m_sigma,mp ] = GMM_ini( data,n_center )

[m,n]=size(data);
[data_id,centers]=kmeans(data,n_center);
mu=centers;
mp=zeros(1,n_center);
m_sigma=zeros(n,n,n_center);

for i=1:n_center
    tem_id=(data_id==i);
    m_sigma(:,:,i)=sigma(data(tem_id,:));
    mp(i)=sum(tem_id)/m;
end

end

function sig=sigma(data)//计算初始化的方差

[m,n]=size(data);
u=mean(data,1);
tem_data=data-repmat(u,m,1);

sig=zeros(n,n);
for k1=1:m
%     for k2=1:m 
    sig=sig+tem_data(k1,:)'*tem_data(k1,:);
%     end
end
sig=(sig+ 1E-5.*diag(ones(n,1)))/m;
end

·高斯概率分布函数

function gp=GaussianPDF(data,u,sigma)

[m,n]=size(data);

pre_item=1/sqrt(((2*pi)^n)*abs(det(sigma)+realmin));
nxt_item(1:m)=0;
tem_data=data-repmat(u,m,1);
for i=1:m
tem_data_t=tem_data(i,:)';
nxt_item(i)=exp(-0.5*(tem_data(i,:)*(inv(sigma))*tem_data_t));
end

gp=pre_item*nxt_item;

end

·EM算法函数

function [mu,msigma,mp]=GMM(data,n_center,loglik_threshold)

[ mu,msigma,mp ] = GMM_ini( data,n_center );
disp('GMM_Ini Completed ! ');

Qt=E_step(data,mu,msigma,mp);
loglik_pre=loglike(data,mu,msigma,mp);
step=0;

while 1
[mu,msigma,mp]=M_step(Qt,data);
loglik_nxt=loglike(data,mu,msigma,mp);
if abs((loglik_nxt/loglik_pre)-1) < loglik_threshold
break;
end

if step>4
break;
end

step=step+1;
step
loglik_pre=loglik_nxt;
Qt=E_step(data,mu,msigma,mp);

end

end

function Qt=E_step(data,mu,m_sigma,mp)//E_STEP

n_model=length(mp);
m=size(data,1);
pxj(m,n_model)=0;

for j=1:n_model
pxj(:,j)=GaussianPDF(data,mu(j,:),m_sigma(:,:,j));
end

px=pxj.*repmat(mp,m,1);
sp=sum(px,2);
Qt=px./repmat(sp,1,n_model);

end

function [lu,lsigma,lp]=M_step(Qt,data)//M_STEP

[m,n_model]=size(Qt);
n=size(data,2);

lu=zeros(n_model,n);
lsigma=zeros(n,n,n_model);
lp=zeros(1,n_model);

mul_data=zeros(n,n);

for j=1:n_model
lu(j,:)=sum(data.*repmat(Qt(:,j),1,n))/sum(Qt(:,j));
tem_data=data-repmat(lu(j,:),m,1);
for k=1:m
mul_data=mul_data+tem_data(k,:)'*tem_data(k,:)*Qt(k,j);
end
lsigma(:,:,j)=realmin+mul_data/sum(Qt(:,j));
lp(j)=sum(Qt(:,j))/m;
end

end

function loglik=loglike(data,mu,msigma,mp)//似然值

n_center=size(mu,1);
pxj=zeros(size(data,1),n_center);
for j=1:n_center
pxj(:,j) = GaussianPDF(data, mu(j,:), msigma(:,:,j));
end
F = pxj*mp';
F(F<realmin) = realmin;
loglik = log(sum(F));

end

·测试函数

clear all;
clc;
data=rand(1000,128);//1000个128维的数据样本
n_center=4;
thresh=0.0005;
[u,sigma,p]=GMM(data,n_center,thresh);

disp('Test Completed !');

注意，模型数在3-5左右，阀值要在0.0005-0.001，否则容易得到奇异方差矩阵。