%% Problem Sheet 7
% Stochastic Simulation

close all; clear all; clc;

%% Exercise 3

N = 10^5; % sample size

% a)

% transition matrix + initial distribution
P = [1/6, 2/6, 2/6, 1/6;
     0, 1/2, 1/2, 0;
     1/6, 2/6, 2/6, 1/6;
     1/2, 0, 0, 1/2];
lambda = [1/4, 1/4, 1/4, 1/4];

% money won/lost in each state
money = [5, 10, -5, -10];
 
% simulate the Markov chain
X = zeros(N,1);
X(1) = find(rand < cumsum(lambda), 1);
for i = 2:N
    X(i) = find(rand < cumsum(P(X(i-1),:)), 1);
end

% compute money won/lost per match
money_per_match = mean(money(X));
fprintf('On average, you have won %f Euro per match.\n', money_per_match);

% b)

% compute and plot running average
money_avg = zeros(N/10,1);
for i = 1:(N/10)
    money_avg(i) = mean(money(X(1:i)));
end
figure(1);
plot(money_avg);
title('Running average of money won or lost per match');


%% Exercise 4

n = 6; % number of vertices
m = 5; % number of colors
N = 10^4; % sample size

% adjacency matrix
E = [0 1 1 0 0 0
     1 0 0 1 1 0
     1 0 0 1 0 1
     0 1 1 0 1 0
     0 1 0 1 0 1
     0 0 1 0 1 0];

% initial coloring and vector for estimation
X = [1 4 4 3 5 2]; 
atmost4 = zeros(N,1);

% Gibb's sampler
for i=1:N
    Y = X;
    % choose vertex at random
    v = ceil(n*rand());
    % collect colors used by neighbors
    used = [];
    for j=1:n
        if E(v,j) == 1
            used = [used, X(j)];
        end
    end
    % draw uniformly from remaining colors
    not_used = setdiff(1:m, used);
    index = ceil(length(not_used)*rand());
    Y(v) = not_used(index);
    % update X and remember if 4 or less colors have been used
    X = Y;
    atmost4(i) = (sum(hist(X,1:m)>0) <= 4);
end

fprintf('The fraction with 4 colors or less is %f\n', mean(atmost4));