clear
rng(0)
%First, let's make the data
    p = [0.9 , 0.1 ; 0.05 , 0.95];

    mu_1 = 5;
    mu_2 = 10;
    sigma_1 = 2;
    sigma_2 = 10;
    
    T = 100;
    x_sim(1) = 1;
    for t = 1:T
        temp = rand;
        if x_sim(t) == 1 & temp < p(1,1)
            x_sim(t+1) = 1;
        elseif x_sim(t) == 1 & temp >= p(1,1)
            x_sim(t+1) = 2;
        end

        if x_sim(t) == 2 & temp < p(2,2)
            x_sim(t+1) = 2;
        elseif x_sim(t) == 2 & temp >= p(2,2)
            x_sim(t+1) = 1;
        end
        
        if x_sim(t) == 1
            y_sim(t) = mu_1+sqrt(sigma_1)*randn;
        elseif x_sim(t) == 2
            y_sim(t) = mu_2+sqrt(sigma_2)*randn;
        end
    end;
    
%Plot the data

%Estimate the likelihoods, given beliefs
    x_belief = [1;0];
    %pr(X|Y) = pr(X)*pr(Y|X)/pr(Y)

for t = 2:T
        %Pr in state 1:
            pr1 = x_belief(1,t-1).*p(1,1)+x_belief(2,t-1).*p(2,1);
        %Pr in state 2:
            pr2 = x_belief(1,t-1).*p(1,2)+x_belief(2,t-1).*p(2,2);
            
        %Bayes Rule
            x_belief(1,t) = pr1*pdf('norm',y_sim(t),mu_1,sigma_1)./(pr1*pdf('norm',y_sim(t),mu_1,sigma_1)+pr2*pdf('norm',y_sim(t),mu_2,sigma_2));
            x_belief(2,t) = pr2*pdf('norm',y_sim(t),mu_2,sigma_2)./(pr1*pdf('norm',y_sim(t),mu_1,sigma_1)+pr2*pdf('norm',y_sim(t),mu_2,sigma_2));
end
    
    subplot(2,1,1)
    plot(x_sim(1:T)-1)
    hold on
    plot(x_belief(2,:),'--r')
    ylim([-0.1,1.1])
    subplot(2,1,2)
    plot(y_sim)
    