close all 
clear all

%% Part a
trials = 10000;
vector = rand(1,trials);
P = zeros(2,2);
for val = 1:trials;
    if vector(val) < (1/8)
        P(1,1) = P(1,1) + 1;
    elseif vector(val) < (2/8)
        P(2,1) = P(2,1) + 1;
    elseif vector(val) < (1/2)
        P(1,2) = P(1,2) + 1;
    elseif vector(val) <= (1)
        P(2,2) = P(2,2) + 1;
    end
end
P = P./trials;
P
%% Part b
trials = 10000;
vector = rand(1,trials);
Px = zeros(2);
Py = zeros(2);
for val = 1:trials;
    if vector(val) < (3/8)
        Px(1) = Px(1) + 1;
    elseif vector(val) <= (1)
        Px(2) = Px(2) + 1;
    end
end
for val = 1:trials;
    if vector(val) < (1/4)
        Py(1) = Py(1) + 1;
    elseif vector(val) <= (1)
        Py(2) = Py(2) + 1;
    end
end
Px = Px./trials;
Py = Py./trials;
P(1,1) = Px(1)*Py(1);
P(2,1) = Px(2)*Py(1);
P(1,2) = Px(1)*Py(2);
P(2,2) = Px(2)*Py(2);
P
%% Part c
Excalc = 0*(3/8) + 1*(5/8);
Ex = P(2,1) + P(2,2);
Ex
Excalc
%The predicted match the simulated values.
%% Part d
Varxcalc = Excalc - Excalc^2; %E[X^2] = Excalc
Varxcalc
Varx = Ex - Ex^2;
Varx