高斯混合模型、EM参数估计及其代码

高斯混合模型

有时随机变量的分布不太符合一个尖的高斯分布。比如一个班级的学生的身高，它的频率分布直方图应该就有两个尖，一个尖是男生的平均身高，另一个尖是女生的平均身高。对于这样的数据，用混合高斯模型来拟合就比用高斯分布模型更加合适。
混合高斯模型是多个高斯模型的叠加，它的概率密度函数为：

p(X|Y)=∑k=1KπkN(x|μk,Σk)

其中 N()表示高斯分布，Y是表示某个样点x来自于哪一个高斯成分的隐变量，Y也是个随机变量，它的k维中有且只有一位为1，其余都为0，表示样点x来自为1的那一位所代表的高斯成分。且满足

p(Y(k)=1)=πk,∑k=1Kπk=1

EM参数估计

假设一个随机变量服从混合高斯分布，则可以用最大化似然函数L(π,μ,Σ)=p(π,μ,Σ|X,Y)的ML(maximum likelyhood)方法来估计这个分布的参数。但是这要求X和Y都已知，然而实际中往往只能观测到随机变量X的样点x，却不知道各个样点x到底属于混合高斯分布中的哪个高斯成分。因此，用EM方法来估计参数：E，对Y取数学期望，消去这个参数；M，最大化期望似然函数。EM参数估计会收敛到似然函数的局部极大值，而不是全局的，所以选择合适的初始参数很重要。

EM的一般方法：

1. Expectation：在X和当前的参数θt已知的条件下，将log似然函数对Y取期望
   定义Q(θ,θt)为
   
   Q(θ,θt)=E(log[L(θ)]|X,θt)=∫γlog(L(θ))f(Y|X,θt)dy
2. Maximum：极大化期望似然函数，求出最优的t，作为新的t进行下一轮迭代。。。
   
   θt+1=maxθQ(θ,θt)

EM在混合高斯模型中的具体操作

p(k|x)=p(x,k)p(x)=p(x|k)p(k)/p(x)

(当然这是指的在t时刻的参数下的p(k|x)，要用它来迭代求下一次的参数)

αt+1k=1n∑i=1Np(k|xi)

μt+1k=∑Ni=1p(k|xi)xi∑i=1Np(k|xi)

Σt+1k=∑Ni=1p(k|xi)(xi−μt+1k)T(xi−μt+1k)∑Ni=1p(k|xi)

用EM求解混合高斯模型的Matlab代码

%%=======================混合高斯模型==============================%%%kmeans初始化K = 2;D=3;N=spnum;[data_id, centers] = kmeans(seg_vals, K);   mu = centers;Sigma = zeros(3,3,K);alpha = zeros(1,K);    for k = 1:K    Sigma(:,:,k) = cov(seg_vals(data_id==k,:));    alpha(k) = numel(find(data_id==k));endalpha = alpha/sum(alpha);niter = 1;maxiter = 1000;loglik_th = 10^-5;loglik_old = realmin;while 1    %E step    %p(x|k)第k个高斯，每个点在第k个高斯中的概率    pxk = zeros(N,K);    for k = 1:K        pxk(:,k) = mvnpdf(seg_vals, mu(k,:), Sigma(:,:,k));    end    %p(k|x) = p(x,k)/p(x) = p(x|k)p(k)/p(x)    pkx = bsxfun(@times,pxk,alpha);    pkx = bsxfun(@times,pkx,1./sum(pkx,2));    E = sum(pkx);    %M step    %更新参数    for k = 1:K        alpha(k) = E(k) / N;        mu(k,:) = pkx(:,k)' * seg_vals / E(k);        tdata = bsxfun(@minus,seg_vals, mu(k,:));%尼玛居然有这种函数，以前都蠢蠢地repmat呢        Sigma(:,:,k) = bsxfun(@times,tdata,pkx(:,k))'*tdata / E(k);    end    %判断收敛    %k个高斯分布下各个数据的概率密度函数值    for k = 1:K        pxk(:,k) = mvnpdf(seg_vals, mu(k,:), Sigma(:,:,k));    end    %log likelyhood    F = pxk *alpha';    F(F<realmin)=realmin;    loglik = mean(log(F));    if abs(loglik/loglik_old-1) < loglik_th || niter>maxiter        break;    end    loglik_old = loglik;    niter = niter+1;end%画图figure;hold on;scatter(seg_vals(:,1),seg_vals(:,2));drawGaussian(mu(1,[1,2]),Sigma([1,2],[1,2],1));drawGaussian(mu(2,[1,2]),Sigma([1,2],[1,2],2));%%======================莫听穿林打叶声，何妨吟啸且徐行=======================%%%绘图函数（此函数为复制粘贴来自网络)function Z=drawGaussian(u,v)% u,vector,expactation;v,covariance matrixx=linspace(0,1);y=linspace(0,1);[X,Y]=meshgrid(x,y);DX=v(1,1);     %X的方差dx=sqrt(DX);DY=v(2,2);     %Y的方差dy=sqrt(DY);COV=v(1,2);     %X Y的协方差r=COV/(dx*dy);part1=1/(2*pi*dx*dy*sqrt(1-r^2));p1=-1/(2*(1-r^2));px=(X-u(1)).^2./DX;py=(Y-u(2)).^2./DY;pxy=2*r.*(X-u(1)).*(Y-u(2))./(dx*dy);Z=part1*exp(p1*(px-pxy+py));contour(x,y,Z);   

结果

博客～http://www.xyzeng.pub/