%generatehmmmodel.m 
%p’s are state transition probabilities %and h are generation probabilities 
%p11->transition from 1 to 1 
%p12->transition from 1 to 2 
%p21->transition from 2 to 1 
%p22->transition from 2 to 2 
%h11->being in the state 1, 1 is generated 
%h12->being in the state 1, 2 is generated 
%h21->being in the state 2, 1 is generated 
%h22->being in the state 2, 2 is generated 
prior=1; p11=0.4; p12=1-p11; p21=0.6;
p22=1-p21; h11=0.3; h12=1-h11; h21=0.5; h22=1-h21; 
REF=[prior p11 p21 h11 h21];
[state,data]=hmmgenerateseq(p11,p21,h11,h21,prior);
%Given the sequence data with state as hidden, 
%HMM is obtained as follows. 
%Initialize prior,p11,p21,h11,h21 and 
%the intial state is set to 1 %h11 and h21 are the generating probabilities 
%p11 and p21 are the state transition probabilities 

%prior=1; p11=rand; p12=1-p11; p21=rand; p22=1-p21; 
%h11=rand; h12=1-h11; h21=rand; h22=1-h21; 

%Expectation stage %In this r(n,i) e(n,i,n-1,j) and r(1,i) are obtained. 
%Given the sequence, r(n,i) gives the responsibility factor that 
%nth frame being generated with ith state 
%e(n,i,n-1,j) gives the probability that nth frame generated 
%from the ith state and (n-1)th frame generated from the jth state. 
%Maximization stage gives the actual state transition probability 
%and the generating probability 
%alpha(n,k) is the probability p(x1 x2 ... xn,znk) 
%(i.e) nth frame generated from kth state 
PX=[]; 
COL=[]; 
for iter=1:1:50 
    if(data(1)==1) 
        alpha(1,1)=prior*h11;
        alpha(1,2)=(1-prior)*h21;
    else
        alpha(1,1)=prior*h12;
        alpha(1,2)=(1-prior)*h22;
end
for n=2:1:length(data)
    if(data(n)==1)
        alpha(n,1)=h11*(alpha(n-1,1)*p11+alpha(n-1,2)*p21);
        alpha(n,2)=h21*(alpha(n-1,1)*p12+alpha(n-1,2)*p22);
    else
        alpha(n,1)=h12*(alpha(n-1,1)*p11+alpha(n-1,2)*p21); 
        alpha(n,2)=h22*(alpha(n-1,1)*p12+alpha(n-1,2)*p22);
    end
end
%beta(n,k) is the probability p(xn+1 xn+2 ... xN/Znk) 
%beta(N,1) and beta(N,2) are considered as 1 
beta(length(data),1)=1; 
beta(length(data),2)=1; 
for m=length(data)-1:-1:1 
    if(data(m+1)==1) 
        beta(m,1)=(h11*beta(m+1,1)*p11+h21*beta(m+1,2)*p12);
        beta(m,2)=(h11*beta(m+1,1)*p21+h21*beta(m+1,2)*p22);
    else
        beta(m,1)=(h12*beta(m+1,1)*p11+h22*beta(m+1,2)*p12);
        beta(m,2)=(h12*beta(m+1,1)*p21+h22*beta(m+1,2)*p22);
    end
end
%probability of data generation is represented as px 
%computation of gamma 
px=alpha(length(data),1)*beta(length(data),1)+alpha(length(data),2)*beta(length(data),2);
for n=1:1:length(data) 
    gamma(n,1)=(alpha(n,1)*beta(n,1))/(px); 
    gamma(n,2)=(alpha(n,2)*beta(n,2))/(px);
end
%computation of e(n,i,j) (nth data generated from ith state 
%and (n-1) the data generted from jth state. 
for n=2:1:length(data)
if(data(n)==1)
    e(n,1,1)=alpha(n-1,1)*h11*beta(n,1)*p11;
    e(n,1,2)=alpha(n-1,2)*h11*beta(n,1)*p21;
    e(n,2,1)=alpha(n-1,1)*h21*beta(n,2)*p12;
    e(n,2,2)=alpha(n-1,2)*h21*beta(n,2)*p22;
else
    e(n,1,1)=alpha(n-1,1)*h12*beta(n,1)*p11;
    e(n,1,2)=alpha(n-1,2)*h12*beta(n,1)*p21;
    e(n,2,1)=alpha(n-1,1)*h22*beta(n,2)*p12;
    e(n,2,2)=alpha(n-1,2)*h22*beta(n,2)*p22;
end
end
%Maximization stage (Estimate p11,p12,h11,h12, prior... 
%using e,gamma 
prior=gamma(1,1)/(gamma(1,1)+gamma(1,2)); 
num1=0; 
den1=0;
num2=0;
den2=0; 
for n=2:1:length(data)
    num1=num1+e(n,1,1); 
    num2=num2+e(n,1,2);
    for i=1:2 
        den1=den1+e(n,i,1); 
        den2=den2+e(n,i,2); 
    end
end
p11=(num1)/(den1);
p21=(num2)/(den2);
[u v]=find(data==1);
h11=sum(gamma(v,1))/sum(gamma(:,1));
h12=1-h11; 
h21=sum(gamma(v,2))/sum(gamma(:,2));
h22=1-h21;
COL{iter}=[prior p11 p21 h11 h21];
PX=[PX px];
end
[p,q]=max(gamma');
length(find((data-q)==0))
display('ORIGINAL probabilities of the HMM Model')
REF 
display('Probabilities obtained using EM algorithm')
COL{iter} 
figure
plot(PX) 
xlabel('Iterations') 
ylabel('Generating probability')
