EM简单实现

format long;clc;clear all;X=load('datasetoriginal.txt');[num,dim]=size(X);%disp(X);%参数初始化mu1=[170.0 60.0];mu2=[160.0 50.0];sigma1=[5.0 2.0;2.0 5.0];sigma2=[5.0 2.0;2.0 5.0];fai=0.5;Qb=zeros(num,1);Qg=zeros(num,1);L1=zeros(50,1);L2=zeros(50,1);prel=-3000;for flag=1:50%%%% E步 %%%%    for i=1:num    %男生    a=mvnpdf(X(i,1:2),mu1,sigma1)*fai;    b=mvnpdf(X(i,1:2),mu2,sigma2)*(1-fai);    Qb(i)=a/(a+b);    %女生    Qg(i)=b/(a+b);    end%%%% M步 %%%%%计算Fai值fai=sum(Qb)/num;%计算Mu1和Mu2mu1=Qb'*X(:,1:2)/sum(Qb);mu2=Qg'*X(:,1:2)/sum(Qg);%计算Sigma1和sigma2Ssum_1=[0 0;0 0];Ssum_2=[0 0;0 0];        for i=1:num    Ssum_1=Ssum_1+( (X(i,1:2)-mu1)' * (X(i,1:2)-mu1) )* Qb(i);    Ssum_2=Ssum_2+( (X(i,1:2)-mu2)' * (X(i,1:2)-mu2) )* Qg(i);    endsigma1=Ssum_1/sum(Qb);sigma2=Ssum_2/sum(Qg);%判断收敛,并画出L函数的增长图像    for i=1:num    L1(flag)=L1(flag) + Qb(i) * log(mvnpdf(X(i,1:2),mu1,sigma1)*fai/Qb(i)) + Qg(i) * log(mvnpdf(X(i,1:2),mu2,sigma2)*(1-fai)/Qg(i));    end    L2(flag)=flag;    if L1(flag)-prel<0.0001        x=L2(1:flag);        y=L1(1:flag);        plot(x,y);        break;    end    prel=L1(flag);endresult=zeros(num,1);for i=1:num    if Qb(i)>Qg(i)        result(i)=1;    else        result(i)=0;    endendcornum=0;for i=1:num    if result(i)==X(i,3)        cornum=cornum+1;    endend