%% Problem Sheet 1
% Stochastic Simulation

%% Exercise 5, a)

P = [1/2, 1/4, 1/4;
     1/4, 1/2, 1/4;
     1/4, 1/4, 1/2]; 
mu = [1, 0, 0];

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

% draw X_0 from mu
U = rand();
X(1) = find(U < cumsum(mu), 1);
% simulate the Markov chain
for i = 1:n
    U = rand();
    X(i+1) = find(U < cumsum(P(X(i),:)), 1);
end

figure(1);
hist(X, 1:3);
title('Histogram of visited states');
xlabel('state');
ylabel('absolute frequency');

% % for relative frequencies
% figure(2);
% h = hist(X, 1:3);
% bar(1:3, h/sum(h), 1);
% title('Histogram of visited states');
% xlabel('state');
% ylabel('relative frequency');


%% Exercise 5, b)

for k = [5, 10, 365]
    fprintf('Distribution for n=%d: [%f %f %f]\n', k, mu*P^k);
end


%% Exercise 6

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

X(1) = 1;
for i=1:n
    U = rand();
    if U <= (X(i)-1)/X(i)
        X(i+1) = X(i)-1;
    else
        X(i+1) = X(i)+1;
    end
end

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