%% Problem Sheet 1
% Methods of Monte Carlo Simulation


%% Exercise 1 (surface area estimation)

close all; clear all; clc;

% circle centers and radius
center1 = [0, 0.5];
center2 = [1, 0.5];
r = 1/sqrt(2);

% sample sizes
N = [10^2, 10^3, 10^5];

for n = N
    % initialize vector for counting how many points fell inside
    insideA = zeros(n,1);
    for i = 1:n
        % draw random point in (0,1)^2
        x = rand(1,2);
        % check if it is in the gray region
        insideA(i) = (sum((center1-x).^2) <= r^2) && (sum((center2-x).^2) <= r^2);
    end
    fprintf('The estimated area using %d points is: %f\n', n, mean(insideA));
end
fprintf('The true surface area is: %f\n', pi/4 - 1/2);


%% Exercise 2 (RANDU)

close all; clear all;

% sample size
N = 10^4;
X = zeros(N,1);
S = zeros(N+1,1);

% seed
S(1) = 15;

% RANDU
a = 2^16 + 3;
c = 0;
m = 2^31;

for i = 1:N    
    S(i+1) = mod(a*S(i) + c, m);
    X(i) = S(i+1) / m;
end

% 1-dim plot
figure(1)
plot(X);

% 2-dim plot
figure(2)
scatter(X(1:N-1), X(2:N));

% 3-dim plot
figure(3)
scatter3(X(1:N-2), X(2:N-1), X(3:N));


%% Exercise 3 (Monte Carlo Integration)

close all; clear all;

% sample size
N = 10^4;   
M = 3 * N;  

U = zeros(M, 1);
V = zeros(M+3, 1);
W = zeros(M+3, 1);

% parameters of generator
m1 = 2^32 - 209;
m2 = 2^32 - 22853;
a11 = 1403580;
a12 = 810728;
a21 = 527612;
a22 = 1370589;

% seed
V(1) = 1359;
V(2) = 4976;
V(3) = 1279;
W(1) = 5879;
W(2) = 8265;
W(3) = 2269;

for i = 1:M 
    V(i+3) = mod(a11 * V(i+1) - a12 * V(i), m1);
    W(i+3) = mod(a21 * W(i+2) - a22 * W(i), m2);
    if V(i+3) <= W(i+3)
        U(i) = (V(i+3) - W(i+3) + m1) / (m1 + 1);
    else
        U(i) = (V(i+3) - W(i+3)) / (m1 + 1);
    end
end

% split into 3 vectors
X1 = U(1:3:end);
X2 = U(2:3:end);
X3 = U(3:3:end);

% plug into function
Y = cos(X3.*X2) .* sin(cos(X3.*X1));

int_est = mean(Y);
fprintf('The estimated value of the integral is: %f\n', int_est);