clear;
close all;
clc;

%% Monte Carlo P(x^2+2x-1 > 0)

N = 1e6;
x = randn(1,N);
num_of_succes = 0;

for i = 1:N
    if(g1(x(i)) > 0)
        num_of_succes = num_of_succes + 1;
    end
end

num_of_succes / N

%% Monte Carlo P(sum x^2 > 1)

num_of_succes = 0;

% Three independent samples
x1 = randn(1,N);
x2 = randn(1,N);
x3 = randn(1,N);

for i = 1:N
    if(g2(x1(i), x2(i), x3(i)) > 1)
        num_of_succes = num_of_succes + 1;
    end
end

num_of_succes / N

%% Bootstrap

x = [0.5377, -1.3499, 0.6715, 0.8884, -0.1022, -0.8637, -1.0891, -0.6156, 1.4193, 0.7254];

num_of_succes = 0;
for i = 1:N
    % Generate three random indices
    idx = randi(length(x), 1, 3);
    if g2(x(idx(1)), x(idx(2)), x(idx(3))) > 1
        num_of_succes = num_of_succes + 1;
    end
end

num_of_succes / N

%% Bootsrap with variable Threshold

u = rand(1,N);
T = 0:0.1:3;

P = zeros(1, length(T));
num_of_succes = 0;

y_k = zeros(1, N);
y_k(1) = u(1)^2;
y_k(2) = u(2)^2+u(1)^2;

for i=3:N
    for j=0:2
        y_k(i) = y_k(i) + u(i-j)^2;
    end
end

for t = 1:length(T)
    num_of_succes = 0;
    for i = 1:N
       if y_k(i) < T(t)
         num_of_succes = num_of_succes + 1;
       end
    end
    P(t) = num_of_succes / N;
end

plot(T, P, LineWidth=2);
xlabel('T');
ylabel('P');
grid on;
title('Plot of P vs T');

%% Function Definition
function result = g1(x)
    result = x^2+2*x-1;
end

function result = g2(x1, x2, x3)
    result = x1^2+x2^2+x3^2;
end