Gammasum=zeros(1,K);
Scale=zeros(T,1);
Xi=zeros(T-1,K*K);
for n=1:N
iCov=inv(Cov);
k2=k1/sqrt(det(Cov));
for i=1:T
for l=1:K
d=Mu(l,:)-X((n-1)*T+i,:);
B(i,l)=k2*exp(-0.5*d*iCov*d');
end;
end;
scale=zeros(T,1);
alpha(1,:)=Pi.*B(1,:);
scale(1)=sum(alpha(1,:));
alpha(1,:)=alpha(1,:)/scale(1);
for i=2:T
alpha(i,:)=(alpha(i-1,:)*P).*B(i,:);
scale(i)=sum(alpha(i,:));
alpha(i,:)=alpha(i,:)/scale(i);
end;
beta(T,:)=ones(1,K)/scale(T);
for i=T-1:-1:1
beta(i,:)=(beta(i+1,:).*B(i+1,:))*(P')/scale(i);
end;
gamma=(alpha.*beta);
gamma=rdiv(gamma,rsum(gamma));
gammasum=sum(gamma);
xi=zeros(T-1,K*K);
for i=1:T-1
t=P.*( alpha(i,:)' * (beta(i+1,:).*B(i+1,:)));
xi(i,:)=t(:)'/sum(t(:));
end;
Scale=Scale+log(scale);
Gamma=[Gamma; gamma];
Gammasum=Gammasum+gammasum;
Xi=Xi+xi;
end;
%%%% M STEP
% outputs
Mu=zeros(K,p);
Mu=Gamma'*X;
Mu=rdiv(Mu,Gammasum');
% transition matrix
sxi=rsum(Xi')';
sxi=reshape(sxi,K,K);
P=rdiv(sxi,rsum(sxi));
% priors
Pi=zeros(1,K);
for i=1:N
Pi=Pi+Gamma((i-1)*T+1,:);
end
Pi=Pi/N;
% covariance
Cov=zeros(p,p);
for l=1:K
d=(X-ones(T*N,1)*Mu(l,:));
Cov=Cov+rprod(d,Gamma(:,l))'*d;
end;
Cov=Cov/(sum(Gammasum));
oldlik=lik;
lik=sum(Scale);
LL=[LL lik];
fprintf('cycle %i log likelihood = %f ',cycle,lik);
if (cycle<=2)
likbase=lik;
elseif (lik<oldlik)
fprintf('violation');
elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)|~finite(lik))
fprintf('\n');
break;
end;
fprintf('\n');
end
**************************************************************
echo on;
clc;
% load and plot data on geyser eruption durations and waiting times
load geyser.m;
clf;
subplot(211);
plot(geyser(:,1),geyser(:,2),'o');
xlabel('dur'); ylabel('wait time');
axis('square'); hold on;
% Hit any key to continue
pause;
clc;
% divide up dataset
Xtrain=geyser(1:200,:); Xtest=geyser(201:295,:);
% train HMM with 3 states for 30 cycles of EM or until convergence
[Mu,Cov,P,Pi,LL]=hmm(Xtrain,200,3,30);
% plot log likelihood (log L) per sample
subplot(212);
plot(LL/200);
ylabel('Log likelihood per sample');
xlabel('Iterations of EM');
hold on;
% calculate log L for test data
lik=hmm_cl(Xtest,95,3,Mu,Cov,P,Pi);
plot(length(LL),lik/95,'go');