%% Problem Sheet 9
% Stochastic Simulation

close all; clear all; clc;

%% Exercise 3

% initialization
N = 10^4;
X = zeros(N+1, 1);
X(1) = 0;
f = @(x) (exp(-(x-2).^2/2) + exp(-(x+2).^2/2))/(2*sqrt(2*pi));

for i = 1:N
    Y = rand() * f(X(i));
    % shift function f such that g is zero if f(x)=Y
    g = @(x) f(x) - Y;
    % compute start/end points of interval(s)
    a = fzero(g, 1.5);
    b = fzero(g, 2.5);
    if abs(b-a)<10^-6 % found the same point twice -> one large interval
        X(i+1) = 2*a*rand() - a;
    else % two smaller intervals
        x = rand();
        if x > 0.5
            X(i+1) = 2*(x-0.5)*(b-a) + a;
        else
            X(i+1) = 2*(x-0.5)*(b-a) - a;
        end
    end
end

% plot
figure(1);
h = hist(X, -10:0.5:10);
bar(-10:0.5:10, h/(sum(h)*0.5), 1); 
hold on;
x = -10:0.1:10;
plot(x, f(x), 'green');
hold off;
axis([-10, 10, 0, 0.3])


%% Exercise 4

clear all;

f = @(x) x.^2 .* exp(-x/2) .* (x >= 0);
N = 10^5;

X = zeros(N,1);
X(1) = 2;
for i = 2:N
    % proposal: current value + std normal
    Y = X(i-1) + randn();
    % acceptance probability
    alpha = f(Y)/f(X(i-1));
    if rand() < alpha
        X(i) = Y;
    else
        X(i) = X(i-1);
    end
end

% plot
figure(2);
h = hist(X, 0:50);
bar(0:50, h/sum(h), 1);
hold on;
x = 0:0.1:50;
C_est = 1/(sum(f(x))*0.1);
plot(x, C_est*f(x), 'green');
hold off;
axis([-10, 40, 0, 0.2]);