%% Problem Sheet 3
% Methods of Monte Carlo Simulation

clear all; close all; clc;

%% Exercise 1 (discrete acceptance-rejection)

% a) and b)

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

N = 10^4; % sample size

% (i) (uniform distribution)

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

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

Y1 = [];
while length(Y1) < N
    X = floor(11*rand()); % proposal
    U = rand();
    if U <= f(X+1)/(C_i*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 (i)');

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

% (ii) (geometric distribution)

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

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

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

Y2 = [];
while length(Y2) < N
    % generate exponential using inverse transform
    X = floor(-log(rand())/lambda);
    U = rand();
    if (X <= 10) && (U <= f(X+1)/(C_ii*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 (ii)');

% c)

fprintf('The acceptance probabilities are %f for (i) and %f for (ii).\n', 1/C_i, 1/C_ii);
fprintf('=> The geometric proposals should be preferred.\n');


%% Exercise 3 (sampling from a shape)

close all; clear all;

R = 5;
r = 2;

N = 10^4;

rejections = 0;

X = zeros(N,3);
for i = 1:N
    Y(1:2) = rand(1,2)*2*(r+R) - (r+R);
    Y(3) = rand()*2*r - r;
    while (R - sqrt(Y(1)^2 + Y(2)^2))^2 + Y(3)^2 > r^2
        Y(1:2) = rand(1,2)*2*(r+R) - (r+R);
        Y(3) = rand()*2*r - r;
        rejections = rejections+1;
    end
    X(i,:) = Y;
end

figure(1);
plot3(X(:,1), X(:,2), X(:,3), 'o');
axis([-10 10 -10 10 -10 10]);
title('Points uniformly distributed in a torus');

fprintf('The estimated acceptance probability is %f\n', N/(N+rejections));


%% Exercise 6 (Box-Muller)

close all; clear all;

% a)

% see function box_muller

N = 10^5;
X = box_muller(N);

% histogram (normalized, so the area of bars is 1)
figure(1);
w = 0.5; % width of bins
h = hist(X, -20:w:20);
bar(-20:w:20, h/sum(h)/w, 1);
axis([-7, 7, 0, 0.6]);
ylab('relative frequency ()');

% add density
hold on;
x = linspace(-10, 10, 100);
plot(x, (1/sqrt(2*pi))*exp(-x.^2/2), 'g');
hold off;

% b)

mu = 20;
sigma_sq = 100;

X = sqrt(sigma_sq)*box_muller(N) + mu;

figure(2);
hist(X, -40:5:80);
axis([-30, 70, 0, 3e4]);
ylab('absolute frequency');
title(sprintf('Histogram of N(%d,%d)-distributed random numbers', mu, sigma_sq));