EM算法（期望最大算法）实现

该EM算法是本人自己实现。Mablab编程，懂算法原理的话，很容易翻译成C/C++实现。

% Written by Chunjian Ren
% 2010.3.12-14
% all rights reserved

clear

%%-------------------------------------------------------------------------
% 随机产生符合GMM的随机数
n=10000;
pw1=0.6;
u1=-2;
std1=2;
pw2=0.4;
u2=3;
std2=1;

x=zeros（1n）;

y1=randn（1floor（n*pw1））*std1+u1;
y2=randn（1floor（n*pw2））*std2+u2;

x（11:floor（n*pw1））=y1;
x（1（floor（n*pw1）+1）:n）=y2;

z=-10:0.2:10;
r1=hist（xz）;
r1=r1./sum（r1）;

figureplot（zr1‘b-‘）;hold on;

clear pw1 pw2 u1 u2 std1 std2 y1 y2 r1
% %%-------------------------------------------------------------------------
% % GMM的概率密度函数
% r2=pw1*normpdf（xu1std1）+pw2*normpdf（xu2std2）;
% r2=r2./sum（r2）;
% plot（xr2‘r-‘）;hold on;
% 
% clear y1 y2 z r1 r2;
%%-------------------------------------------------------------------------
% EM算法
% n -- 样本数
% z -- 样本值范围
% x -- 样本集

% threshold=1e-15;
iter=0;
% rep=10; % 找10个局部最大值

% [0.6 0.4]  [-2 3]  [2 1]
alpha=[0.6 0.4];  % 混合高斯权重
mu=[-2 3];         % 高斯函数均值
sigma=[4 1];        % 高斯函数方差

class=zeros（n2）;
class（1:n/21）=1;
class（n/2+1:end2）=1;
class_num=sum（class）;
        
% f_best=-inf;
% for crep=1:rep
%     fprintf（‘第 %02d 个局部最大值\n‘crep）;
    while（1）
        M=2;N=n;
        gamma=zeros（MN）;
        for i=1:M
            for j=1:N
                gamma（ij） = 1 / （ sqrt（sigma（i）*2*pi） ） * exp（ -（x（j）-mu（i））^2 / （2*sigma（i）） ）;
            end
        end

        % E step
        pikn=zeros（1M）;
        for j=1:N
            for i=1:M
                pikn（i） = alpha（i） * gamma（ij）;
            end
            for i=1:M
                gamma（ij）=pikn（i）/sum（pikn）;
            end
        end

        % M step
        nk=zeros（1M）; % alpha
        for i=1:M
            for j=1:N
                nk（i）=nk（i）+gamma（ij）;
            end
        end
        alpha=nk/N;

        sum_mu=zeros（M1）; % mu
        for i=1:M
            for j=1:N
                sum_mu（i）=sum_mu（i）+class（ji）*x（j）;
%               sum_mu（i）=sum_mu（i）+gamma（ij）*x（j）;
            end
                mu（i）=sum_mu（i）/class_num（i）;
%               mu（i）=sum_mu（i）/nk（i）;
        end

        sum_sigma=zeros（M1）; % sigma
        for i=1:M
            for j=1:N
                sum_sigma（i）=sum_sigma（i）+class（ji）*（ x（j）-mu（i） ）^2;
%               sum_sigma（i）=sum_sigma（i）+gamma（ij）*（ x（j）-mu（i） ）^2;
            end
                sigma（i）=sum_sigma（i）/class_num（i）;
%               sigma（i）=sum_sigma（i）/nk（i）;
        end
        
        % 更新class
        for j=1:N
            [zhixia] = max（gamma（:j））;
            class（j:）=0;
            class（jxia）=1;
        end
        
        class_num=sum（class）;
        
        iter=iter+1;
        if（iter>100）
            break;
        end
    end % while
    
%     f=sum（log（sum（pikn）））;
%     
%     if f>f_best
%         alpha_best=alpha;
%         mu_best=mu;
%         sigma_best=sigma;
%         f_best=f;
%     end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3442  2010-03-14 23:54  EM_Algorithm\EM.m

     文件       3969  2010-03-12 15:23  EM_Algorithm\EMDownLoad.m

     文件        387  2010-03-14 16:39  EM_Algorithm\MLE.m

     目录          0  2010-03-14 23:56  EM_Algorithm

----------- ---------  ---------- -----  ----

                 7798                    4