sunmog.m

高斯混合模型参数估计,EM算法,sunMOG.m为函数,testMOG4.m为测试程序

% 高斯混合模型参数估计：EM算法
%  重庆大学 学生：孙世军，导师：彭承琳   2008年3月8日  
function [omega mu nu  t  delta]=sunMOG(x, J, maxit, pre)
% function [omega mu nu  t  delta]=sunMOG(x, J, maxit, pre)
% x为一维信号， J为高斯函数个数，maxit为最大迭代次数，pre为精度
% p(x)=g(x,mu(1),nu(1))*w(1) + ... + g(x,mu(J),nu(J))*w(J)
% g为高斯函数，[w(1) ... w(J)]为线性组合系数
% [mu(1)... mu(J)]为各个高斯函数的均值，[nu(1)..nu(J)]为各个高斯函数的均值sigma^2
if nargin<4
    pre=1e-3;
end
if nargin<3
    maxit=388;
end
if nargin<2
    J=3;
end

N=size(x,2);
omega=rand(1,J);
omega=omega/sum(omega);
mu=rand(1,J);
nu=rand(1,J);
delta=1;
for t=1:maxit
    oldomega=omega;
    oldmu=mu;
    oldnu=nu;
    % calculate p_x
    p_x=zeros(1,N);
    for i=1:N
        p_x(i)=px(x(i),oldomega,oldmu,oldnu);
    end
    % calculate omegam, mu, nu
    omega=zeros(1,J);
    mu=zeros(1,J);
    nu=zeros(1,J);
    for j=1:J
        for n=1:N
            omega(j)=omega(j)+oldomega(j)*g(x(n),oldmu(j),oldnu(j))/p_x(n);
            mu(j)=mu(j)+oldomega(j)*g(x(n),oldmu(j),oldnu(j))*x(n)/p_x(n);
        end
    end
    mu=mu./omega;
    for j=1:J
        for n=1:N
            nu(j)=nu(j)+oldomega(j)*g(x(n),oldmu(j),oldnu(j))*(x(n)-mu(j))^2/p_x(n);
        end
    end
    nu=nu./omega;
    omega=omega/N;
    delta=norm(mu-oldmu)/norm(oldmu)+norm(nu-oldnu)/norm(oldnu);
    if delta<pre        
        break;        
    end    
end
[t delta]
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 高斯函数
function p=g(x,mu,nu)
    p=exp(-(x-mu)^2/nu/2) / sqrt(2*pi*nu);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%一维信号x的pdf
function p=px(x,omega,mu,nu)
    J=size(omega,2);
    p=0;
    for j=1:J
        p=p + g(x,mu(j),nu(j))*omega(j);
    end
end