em5.m

用matlab实现高斯混合模型Em算法的源程序哪位有呀,帮帮忙啦- 工具箱与

function tralala= em4
%%This programm estimates the parameters of a mixture of two univariate
%%normals using the EM algorithm and the convergence is demonstrated " on-line".
%%It can be used for any set of data that a mixture of normals can be
%%assumed. We used it for a mixture of Normals emerging from interest rate models with Markovian jumps. 
%R.J. Elliot's and R.S.Mamon's 2-factor Vacicek with mean reverting level
%of two states


%Made by Panagiotis Braimakis as a part of a team assignment.
%Finished at Tuesday 4th, 2005 , 2.23am
clf;
clc;
figure(1);
disp('Lets start some parameter estimation, waiting for your click...>');

clear;

format short g
load data
R=D5';           %LOAD YOUR DATA
n=length(R);   %HOW MANY DATA YOU HAVE
subplot(2,1,1);
%Theoretical pdf
syms p 
m=mean(R);
s=std(R);
g=((2*pi*s^2)^(-1/2))*exp(-((p-m)^2)/(2*s^2));    %THE SIMPLE NORMAL DISTRIBUTION
figure(1);
%View some different aspects of our data
[w,x]=ecdf(R);  %empirical Kamplan-Meier cdf
ecdfhist(w,x,25);
hold on

ezplot(g,[-20,20]);
axis([min(R)-5 max(R)+5 0 .5])
title('Histogram of IR''s with 20 bins selected')
ylabel('Probabilities');
xlabel('Interest Rates');
axis manual
%----------------------------------end plots-------------------------------------------
pause;

%Finding the mixtures likelihood

%random allocation of the zeros, ones augmentation
w=binornd(10,0.5,length(R),1)>5;

%The complete dataset
Y=[R w];

%These are the pdf for each observation ri belonging to two different
%states. the parameters are the (a(j),b,s), bu =t we will only estimate the
%mean and the sigma of the mixture.

syms x sigmaa mu1 mu2 p1 p2
tic;
%notes [p.1]

%we will suppose different means only but the same variance!

f1=((2*pi*sigmaa^2)^(-0.5))*exp(-0.5*((x-mu1)/sigmaa)^2);
f2=((2*pi*sigmaa^2)^(-0.5))*exp(-0.5*((x-mu2)/sigmaa)^2);

fi=[f1,f2].';
p=[p1, p2].';
%notes [p.2]
%mixture pdf, not joint! 2 marginals 
%Relation (1.1)
f=sum(p.*fi);

q=toc;
fprintf('Found the marginal pdf in %f seconds.\n',q)
disp(f);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5UP TO HERE WAS THE SYMBOLIC PART

%__________________________________________________________Something to enter the first loop__________________

ll(1)=-999999999;
ll(2)=-99999999;
tol=10^(-30);
i=1;

while abs((ll(i+1)-ll(i))/ll(i+1))>=tol

%estimates of the mixing probabilities

p1(i)=sum(w)/n;      %THIS VALUE DEPENDS ON THE w MATRIX AND CHANGES IN TIME
P1(i)=double(p1(i));
p2(i)=1-p1(1);
P2(i)=1-P1(i);


%Initial values (they will be re-valued in the M-step)
MU1(1)=0.5;            
MU2(1)=10.4;
Sigma1(1)=0.5;
Sigma2(1)=1.7;
%the likelihood for the complete data set GIVEN THE PARAMETERS MU1,MU2, SIGMA

fun=1/2*P1(i)*2^(1/2)/(pi*Sigma1(i)^2)^(1/2)*exp(-1/2*(R-MU1(i)).^2./Sigma1(i).^2)+1/2*P2(i)*2^(1/2)/(pi*Sigma2(i)^2)^(1/2)*exp(-1/2*(R-MU2(i)).^2./Sigma2(i).^2);

ff=double(fun);

%The log-likelihood
ll(i+2)=sum(log(ff(:)));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%E-step*********************%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%**
%nom=f1
x=R;
nom=((2*pi*Sigma1(i)^2)^(-0.5))*exp(-0.5*((x-MU1(i))./Sigma1(i)).^2);
%den=f
den=ff;
w= P1(i)*nom./den; %the probabilities that ith observation belongs to 1st group (the one with (a1))

%Now we know exactly which ri belongs where!
ww=round(w);

%plot according to the partitioning from the E-step.
subplot(2,1,1);
%****************************************again
[wq,x]=ecdf(R);  %empirical Kamplan-Meier cdf
ecdfhist(wq,x,25);
hold on

ezplot(g,[-20,20]);
axis([min(R)-5 max(R)+5 0 .5])
title('Histogram of IR''s with 20 bins selected')
ylabel('Probabilities');
xlabel('Interest Rates');
axis manual
% %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&7
hold on
pp=(min(R)-3):.1:(max(R)+3);
Dist1=(((2*pi*Sigma1(i)^2)^(-1/2))*exp(-0.5*(pp-MU1(i)).^2./(Sigma1(i).^2)))/(n/round(sum(w)));    %THE SIMPLE NORMAL DISTRIBUTION
Dist2=(((2*pi*Sigma2(i)^2)^(-1/2))*exp(-0.5*(pp-MU2(i)).^2./(Sigma2(i).^2)))/(n/round(sum(1-w)));    %THE SIMPLE NORMAL DISTRIBUTION
plot(pp,Dist1,'c-');
plot(pp,Dist2,'r-');
hold off

%Finally the M-step for the next estimation of our parameters
i=i+1;


MU1(i)=sum(ww.*R)/round(sum(w));
MU2(i)=sum((1-ww).*R)/round(sum(1-w));

wk=1-ww;
%Let''s estimate only one sigma (since we have supposed homoskedastidity in our models)
Sigma1(i)=sqrt(sum(ww.*(R-MU1(i)).^2)/(sum(ww)));
Sigma2(i)=sqrt(sum(wk.*(R-MU2(i)).^2)/(n-sum(ww)));


subplot(2,1,2);
plot(1:1:i,ll(1:i),'k-');
hold on;
title('Log-likelihood evolution until convergence');
ylabel('log-likelihood');
xlabel('steps');
pause(.51);

i;
end %while
fprintf('Convergence achieved in %0.0f iterations! \n\n\nHave a nice day.',i)