% hmm_example.m
%
% Neal Patwari, 4 Nov 2010
%

possible_s   = 3;    % s is in the set {1, 2, ..., possible_s}
possible_x   = 2;    % X is in the set {1, 2, ..., possible_x}

% Pr{s_n | x_n}.  Each row corresponds to a possible value of X_n, each column corresponds to
% a possible values of s_n.
pr_s_given_x = [0.90 0.09 0.01; 
                0.33 0.34 0.33];  

% Markov chain (X_n)  Transition Probability Matrix 
P            = [0.95 0.05; 
                0.07 0.93];

% Initial state probabilities.
pi_0         = [0.5 0.5];

% Simulate the HMC
duration     = 20;
[S,X]        = SimHMC(pr_s_given_x, P, pi_0, duration);  % vector of measurements

% Solve forward approach to determine
% F_n(j) = Pr{S^n = s_n, X_n = j}
% For n=1,2,... and j=1...N
[F, pr_x_given_S]  = ForwardSolveHMC(S, pr_s_given_x, P, pi_0);

Expected_X_given_S = pr_x_given_S * (1:possible_x)' 


% Plot the results.
figure(1)
subplot(2,1,1)
g = plot(1:duration, S,'ko')
set(g,'MarkerSize',9);
set(g,'LineWidth',2);
set(gca,'FontSize',20);
ylabel('Observation s_n');
grid
set(gca,'ylim',[0.5, possible_s + 0.5])
subplot(2,1,2)
h = plot(1:duration, X,'ko', 1:duration, Expected_X_given_S, 'r-*');
set(gca,'FontSize',20);
set(h,'MarkerSize',9);
set(h,'LineWidth',2);
grid
set(gca,'ylim',[0.5, possible_x + 0.5])
legend('Actual','E[X_n | S^n]',0)
ylabel('State X_n')
xlabel('Time n')
