% S. S. Singh, CUED, Oct 2007
% 4F7, Adaptive Filters
% simulating the dishonest casino problem

clear all

P=[0.95 .05; 0.1 0.9];  % state 1 is fair dice, state 2 is unfair dice
Q=[1/6, 1/6, 1/6, 1/6, 1/6, 1/6;
    1/10 1/10 1/10 1/10 1/10 5/10];

noStates=2;  % declare no of hidden states
noObs=6;     % declare number of observations

T=100; % number of dice rolls
pi0=[0.5 .5]';  % initial distribution of hidden state

state=zeros(1,T+1);
state(1)=2;
Pcumsum=cumsum(P,2);
Qcumsum=cumsum(Q,2);
obs=zeros(1,T);

for i=1:T  % now generate the hidden state and obs sequence
    
[temp, state(i+1)]=max(rand < Pcumsum(state(i),:)) ;
[temp, obs(i)]=max(rand < Qcumsum(state(i+1),:)) ;
end


Pi=zeros(noStates,T+1);
Pi(:,1)=pi0;

% do the forward filtering
for i=1:T
    B=diag(Q(:,obs(i)));
    Pi(:,i+1)=B*P'*Pi(:,i);
    Pi(:,i+1)=Pi(:,i+1)/sum(Pi(:,i+1));
    
end

%now backward filter
beta=ones(noStates,T+1);

for i=1:T
    B=diag(Q(:,obs(T-i+1)));
    beta(:,T-i+1)=P*B*beta(:,T-i+2);
    beta(:,T-i+1)=beta(:,T-i+1)/sum(beta(:,T-i+1));  % normalizing for numerical stability
end

% now combine forward and backward filter to get smoother
smoother=beta(:,2:T+1).*Pi(:,2:T+1);
smoother=smoother./(ones(noStates,1)*sum(smoother,1));

% now plot results
figure
subplot(3,1,1)
plot([obs],'.','MarkerSize',14)
title('Observations','FontSize',14)
grid on
axis([0 T 0 7] )

subplot(3,1,2)
plot(state(2:end),'rs','MarkerSize',10)
title('Hidden state: 1=Fair, 2=Biased','FontSize',14)
xlabel('dice rolls')
grid on


subplot(3,1,3)
hold on
grid on
plot(Pi(2,2:end),'b')
plot(smoother(2,:),'--k')
title('Prob of Biased','FontSize',14)
legend('Filter', 'Smoother')