%?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