EM算法MATLAB实现

EM算法.rar EM算法 MATLAB 实现

close all;clear;clc;
%   参考书籍Pattern.Recognition.and.Machine.Learning.pdf
%   http://www.pr-ml.cn
%   lwm@pr-ml.cn
%   2009/10/15
%% 
M=2;          % number of Gaussian
N=200;        % total number of data samples
th=0.000001;  % convergent threshold
K=2;          % demention of output signal

% 待生成数据的参数
a_real =[4/5;1/5];
mu_real=[3 4;
         5 3];
cov_real（::1）=[5 0;
                 0 0.2];
cov_real（::2）=[0.1 0;
                 0 0.1];
% generate the data
x=[ mvnrnd（ mu_real（:1）  cov_real（::1）  round（N*a_real（1）） ）‘    mvnrnd（mu_real（:2）cov_real（::2）N-round（N*a_real（1）））‘];

% for i=1:round（N*a_real（1））
%     while （~（（x（1i）>0）&&（x（2i）>0）&&（x（1i）<10）&&（x（2i）<10）））
%         x（:i）=mvnrnd（mu_real（:1）cov_real（::1）1）‘;
%     end
% end
% 
% for i=round（N*a_real（1））+1:N
%     while （~（（x（1i）>0）&&（x（2i）>0）&&（x（1i）<10）&&（x（2i）<10）））
%         x（:i）=mvnrnd（mu_real（:1）cov_real（::1）1）‘;
%     end
% end

figure（1）plot（x（1:）x（2:）‘.‘）

%这里生成的数据全部符合标准
%% %%%%%%%%%%%%%%%% 参数初始化
a=[1/32/3];
mu=[1 2;2 1];%均值初始化完毕
cov（::1）=[1 0;
            0 1];
cov（::2）=[1 0;
            0 1];%协方差初始化
%% EM Algorothm
% loop
count=0;
figure（2）hold on
while 1
    a_old  = a;
    mu_old = mu;
    cov_old= cov;      
    rznk_p=zeros（MN）;
    for cm=1:M
        mu_cm=mu（:cm）;
        cov_cm=cov（::cm）;
        for cn=1:N
            p_cm=exp（-0.5*（x（:cn）-mu_cm）‘/cov_cm*（x（:cn）-mu_cm））;
            rznk_p（cmcn）=p_cm;
        end
        rznk_p（cm:）=rznk_p（cm:）/sqrt（det（cov_cm））;
    end
    rznk_p=rznk_p*（2*pi）^（-K/2）;
%E step
    %开始求rznk
    rznk=zeros（MN）;%r（Z
    pikn=zeros（1M）;%r（Z
    pikn_sum=0;
    for cn=1:N
        for cm=1:M
            pikn（1cm）=a（cm）*rznk_p（cmcn）;
%           pikn_sum=pikn_sum+pikn（1cm）;
        end
        for cm=1:M
            rznk（cmcn）=pikn（1cm）/sum（pikn）;
        end
    end
        %求rank结束
% M step
    nk=zeros（1M）;
    for cm=1:M
        for cn=1:N
            nk（1cm）=nk（1cm）+rznk（cmcn）;
        end
    end
    a=nk/N;
    rznk_sum_mu=zeros（M1）;
        
    % 求均值MU
    for cm=1:M
        rznk_sum_mu=0;%开始的时候就是错在这里，这里要置零。
        for cn=1:N
            rznk_sum_mu=rznk_sum_mu+rznk（cmcn）*x（:cn）;
        end
        mu（:cm）=rznk_sum_mu/nk（cm）;
    end
    
    % 求协方差COV   
    for cm=1:M
        rznk_sum_cov=zeros（KM）;
        for cn=1:N
            rznk_sum_cov=rznk_sum_cov+rznk（cmcn）*（x（:cn）-mu（:cm））*（x（:cn）-mu（:cm））‘;
        end
        cov（::cm）=rznk_sum_cov/nk（cm）;
    end
         
    t=max（[norm（a_old（:）-a（:））/norm（a_old（:））;norm（mu_old（:）-mu（:））/norm（mu_old（:））;norm（cov_old（:）-cov（:））/norm（cov_old（:））]）;
 
    temp_f=sum（log（sum（pikn）））;
    plot（counttemp_f‘r+‘）
    count=count+1;        
        
    if t        break;
    end
 
end  %while 1
   
hold off
f=sum（log（sum（pikn）））;
    
a_best=a;
mu_best=mu;
cov_best=cov;
f_best=f;

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

     文件    9028817  2009-12-08 18:43  PatternRecognitionAndMachineLearning.pdf

     文件       3765  2009-12-09 15:54  em_my3.m

     文件       3329  2009-12-08 22:43  em_my1.m

     文件       2927  2009-12-08 23:19  em_my2.m

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

              9038838                    4