TRANS = [.95 .05; .1 .9;];

EMIS = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6;...
1/2, 1/10, 1/10, 1/10, 1/10, 1/10];

%To generate a random sequence of states and emissions from the model, use hmmgenerate:
%[seq,states] = hmmgenerate(1000,TRANS,EMIS,'Statenames',{'F','L'});
[seq,states] = hmmgenerate(1000,TRANS,EMIS);

char(seq+'0')
char(states+'0')

figure(1);
subplot(3,1,1);
index1 = find(states==1);
index2 = find(states==2);
plot(index1, ones(size(index1)), 'b.', index2, ones(size(index2)), 'r.');
set(gca, 'box','off');


%likelystates is a sequence the same length as seq.
likelystates = hmmviterbi(seq, TRANS, EMIS);
subplot(3,1,2);
index1 = find(likelystates==1);
index2 = find(likelystates==2);
plot(index1, ones(size(index1)), 'b.', index2, ones(size(index2)), 'r.');
set(gca, 'box','off');



% To test the accuracy of hmmviterbi, compute the percentage of the 
%actual sequence states that agrees with the sequence likelystates.
sum(states==likelystates)/1000



%Estimating Posterior State Probabilities
[PSTATES,logpseq] = hmmdecode(seq,TRANS,EMIS);
subplot(3,1,3);
plot(PSTATES');









%% Estimating transiton and Emission Matrices



% Using hmmestimate.   
%To function hmmestimate requires that you know the sequence of states states that the model went through to generate seq.
%The following takes the emission and state sequences and returns estimates of the transition and emission matrices:

[TRANS_EST, EMIS_EST] = hmmestimate(seq, states)

%Using hmmtrain.   If you do not know the sequence of states
%states, but you have initial guesses for TRANS and EMIS, you can still estimate TRANS and EMIS using hmmtrain.

TRANS_GUESS = [.1 .9; .9 .1];
EMIS_GUESS = [.17 .16 .17 .16 .17 .17;.6 .08 .08 .08 .08 08];

%You estimate TRANS and EMIS as follows:
[TRANS_EST2, EMIS_EST2] = hmmtrain(seq, TRANS_GUESS, EMIS_GUESS)