%% Problem Sheet 4
% Stochastic Simulation

close all; clear all; clc;


%% Exercise 3 (acceptance-rejection)

n = 10; % parameters of the binomial distribution
p = 0.2; 

N = 10^4; % sample size

% a) (uniform distribution)

f = binopdf(0:10, n, p); % binomial pdf
g = ones(1, 11)/11; % uniform pdf

C = max(f./g); % optimal scaling constant

Y1 = [];
while length(Y1) < N
    X = floor(11*rand()); % proposal
    U = rand();
    if U <= f(X+1)/(C*g(X+1)) % condition for accepting
        Y1 = [Y1, X];
    end
end

figure(1);
counts = hist(Y1, 0:10);
bar(0:10, counts/sum(counts), 1);
title('Histogram of values simulated in a)');

figure(2);
bar(0:10, f, 1);
title('Binomial distribution (n=10, p=0.2)');


% b)

p2 = 0.3; % parameter for geometric distribution
lambda = -log(1-p2); % corresponding parameter for exponential distribution

f = binopdf(0:10, n, p); % binomial pdf
g = geopdf(0:10, p2); % geometric pdf (where needed)

C = max(f./g); % optimal scaling constant

Y2 = [];
while length(Y2) < N
    X = floor(-log(rand())/lambda);
    U = rand();
    if (X <= 10) && (U <= f(X+1)/(C*g(X+1)))
        Y2 = [Y2, X];
    end
end

figure(3);
counts = hist(Y2, 0:10);
bar(0:10, counts/sum(counts), 1);
title('Histogram of values simulated in b)');


%% Exercise 4

% a)

close all; clear all;

n = 5; % (3, 5, 10)
N = 10^3;

success = zeros(N,1);
rejections = 0;

for i=1:N
    X = ceil(n*rand(1,n));
    while max(hist(X, 1:n)) > 1
       X = ceil(n*rand(1,n));
       % count rejections
       rejections = rejections+1;
    end
    success(i) = (abs(find(X==1) - find(X==n)) == 1);
end
    
fprintf('The estimated acceptance probability is %f\n', N/(rejections+N));
fprintf('The probability of the cars standing next to each other is %f\n', mean(success));

% theoretical acceptance probability
factorial(n)/(n^n)
% theoretical probability for cars standing next to each other
2*(n-1)*factorial(n-2)/factorial(n)

% b)

close all; clear all;

N = 10^3;
n = 3; % (3, 5, 10)
p = 0.1;

success = zeros(N,1);
rejections = 0;

for i=1:N
    % determine which employees are ill
    ill = (rand(1,n) <= p);
    % acceptance-rejection for the positions
    X = floor((n+1)*rand(1,n));
    % check validity
    h = hist(X,0:n);
    accept = (h(1) == sum(ill)) && (max(h(2:end)) <= 1);
    for j = 1:length(ill)
        if ill(j) && (h(j+1) > 0)
            accept = false;
            break;
        end
    end
    % if not, repeat simulation and check
    while ~accept
        X = floor((n+1)*rand(1,n));
        h = hist(X,0:n);
        accept = (h(1) == sum(ill)) && (max(h(2:end)) <= 1);
        for j = 1:length(ill)
            if ill(j) && (h(j+1) > 0)
                accept = false;
                break;
            end
        end
        % count rejections
        rejections = rejections+1;
    end
    if (h(2)>0) && (h(n+1) > 0)
        success(i) = (abs(find(X==1) - find(X==n)) == 1);
    end
end

fprintf('The estimated acceptance probability is %f\n', N/(rejections+N));
fprintf('The probability of the cars standing next to each other is %f\n', mean(success));


%% Exercise 10

close all; clear all;

% a)

% State space: {0,1}^(5x5)
% Proposal distribution: q_X,Y = 1/(n^2) if X and Y differ in exactly
%                        1 element (zero otherwise)
% Acceptance probability: 1 if the configuration is legal, zero otherwise

% b)

N = 10^5; % sample size
n = 5; % dimension
X = zeros(n,n); % initial state

for i=1:N
    % generate a proposal
    index = ceil(n*rand(1,2));
    Y = X;
    Y(index(1),index(2)) = 1-Y(index(1),index(2));
    % check if it is legal
    legal = true;
    for j=1:n
        for k=1:(n-1)
            if Y(j,k)==1 && Y(j,k+1)==1
                legal = false;
                break;
            end
            if Y(k,j)==1 && Y(k+1,j)==1
                legal = false;
                break;
            end
        end
        if ~legal
            break;
        end
    end
    % if yes, accept; otherwise reject and stay in X.
    if legal
        X = Y;
    end
end
X