%% Problem Sheet 3
% Stochastic Simulation

close all; clear all; clc;

%% Exercise 5

N = 30; % persons
n = [100, 1000, 10000]; % steps
m = max(n);

P = zeros(N+1, N+1);
for i = 1:N
    P(i, i+1) = (N-(i-1))/N;
end
for i = 2:(N+1);
    P(i, i-1) = (i-1)/N;
end

% simulation of MC
X = zeros(m, 1);
X(1) = 0;
for j = 1:m
    U = rand;
    % states and indices are shifted -> +-1
    X(j+1) = find(cumsum(P(X(j)+1,:)) > U, 1) - 1;
end

% histograms
for i = 1:3
    figure(i);
    hist(X(1:n(i)), 0:30);
end

% %%% not asked in exercise
% 
% % nice to know: no limiting distribution
% % the Markov chain is periodic!
% 
% m_even = 1:2:m;
% m_odd = 2:2:m;
% figure(4);
% hist(X(m_even), 0:30);
% figure(5);
% hist(X(m_odd), 0:30);
% %%%


%% Exercise 6

close all; clear all;

% transition matrix and invariant distribution of Y_n
P = [0, 3/4, 1/4;
     1/2, 0, 1/2;
     1/3, 1/3, 1/3]; 
mu = [12/41, 14/41, 15/41]; % stationary initial distribution

% transition matrix of time-reversal
P_hat = zeros(3,3);
for i=1:3
    for j=1:3
        P_hat(i,j) = P(j,i)*mu(j)/mu(i);
    end
end

% a)

n = 50;
X = zeros(n+1,1);

% simulate X_n, i.e., X_0 of the backwards chain
X(1) = find(rand() < cumsum(mu), 1 );
for i=1:n
    % simulate each step using P_hat
    X(i+1) = find(rand() < cumsum(P_hat(X(i),:)), 1);
end

figure(1);
clf
for i=1:n
    line([i-1, i], [X(i), X(i)]);
end
axis([0, n, -1, 5]);

% b)

n = 10^5;
X = zeros(n+1,1);
P_est = zeros(3,3);

X(1) = find(rand() < cumsum(mu), 1);
for i=1:n
    X(i+1) = find(rand() < cumsum(P_hat(X(i),:)), 1);
    % count (backwards) transitions
    P_est(X(i+1),X(i)) = P_est(X(i+1),X(i))+1;
end
% normalize P_est
for i=1:3
    P_est(i,:) = P_est(i,:)/sum(P_est(i,:));
end

P_est
P