clear; %Clear Workspace
close all; %Close all open figures
clc; %Clear the Command Window
N = 1e6; %number of trials for my tests

%% Problem 3.11
%Part a
    %Calculate probabilities for X
    Carlos = sum(round(rand(2, N)));%number of heads for Carlos
    Michael = sum(round(rand(2, N)));%number of heads for Michael
    p = max([Carlos;Michael]);%Chooses the max number of heads from each
    prob_0 = length(find(p==0))/N;
    prob_1 = length(find(p==1))/N;
    prob_2 = length(find(p==2))/N;
    %Plot pmf for X
    figure('Name','pmf X, Prob 3.11 (a)','NumberTitle','off');
    prob_vect = [prob_0 prob_1 prob_2];
    bar(0:(length(prob_vect)-1), prob_vect);
    axis([-.25, 2.25, -.25, .75]);
    fprintf('Problem 3.11, Part a: \n\tP(x=0)=%f\n\tP(x=1)=%f\n\tP(x=2)=%f\n\n', [prob_0, prob_1, prob_2])
    
%Part b
    %Calculate probabilities for X
    Carlos = sum(round(rand(2, N) + 0.25));%number of heads for Carlos w/ P(H) = 3/4
    Michael = sum(round(rand(2, N)));%number of heads for Michael still w/ P(H) = 1/2
    p = max([Carlos;Michael]);%Chooses the max number of heads from each
    prob_0 = length(find(p==0))/N;
    prob_1 = length(find(p==1))/N;
    prob_2 = length(find(p==2))/N;
    %Plot pmf for X
    figure('Name','pmf X, Prob 3.11 (b)','NumberTitle','off');
    prob_vect = [prob_0 prob_1 prob_2];
    bar(0:(length(prob_vect)-1), prob_vect);
    axis([-.25, 2.25, -.25, .75]);
    fprintf('Problem 3.11, Part b: \n\tP(x=0)=%f\n\tP(x=1)=%f\n\tP(x=2)=%f\n\n', [prob_0, prob_1, prob_2])
    
%% Problem 3.20
%Part a (assuming difference > 0)
    %Calculate pmf for X
    diffdice = abs(ceil(6*rand(1,N)) - ceil(6*rand(1,N))); %Calculate the abs diff of two dice tosses
    prob_0 = length(find(diffdice==0))/N;
    prob_1 = length(find(diffdice==1))/N;
    prob_2 = length(find(diffdice==2))/N;
    prob_3 = length(find(diffdice==3))/N;
    prob_4 = length(find(diffdice==4))/N;
    prob_5 = length(find(diffdice==5))/N;
    fprintf('Problem 3.20, Part a (assuming difference > 0): \n\tP(x=0)=%f\n\tP(x=1)=%f\n\tP(x=2)=%f\n\tP(x=3)=%f\n\tP(x=4)=%f\n\tP(x=5)=%f\n\n', [prob_0, prob_1, prob_2, prob_3, prob_4, prob_5])
    %Plot pmf for X
    figure('Name','pmf X, Prob 3.20 (a)','NumberTitle','off');
    prob_vect = [prob_0 prob_1 prob_2, prob_3, prob_4, prob_5];
    bar(0:(length(prob_vect)-1), prob_vect);
    axis([-.25, 5.25, -.25, .5]);
%Part a (assuming difference can be negative)
    %Calculate pmf for X
    diffdice = ceil(6*rand(1,N)) - ceil(6*rand(1,N)); %Calculate the abs diff of two dice tosses
    prob_neg1 = length(find(diffdice==-1))/N;
    prob_neg2 = length(find(diffdice==-2))/N;
    prob_neg3 = length(find(diffdice==-3))/N;
    prob_neg4 = length(find(diffdice==-4))/N;
    prob_neg5 = length(find(diffdice==-5))/N;
    prob_0 = length(find(diffdice==0))/N;
    prob_1 = length(find(diffdice==1))/N;
    prob_2 = length(find(diffdice==2))/N;
    prob_3 = length(find(diffdice==3))/N;
    prob_4 = length(find(diffdice==4))/N;
    prob_5 = length(find(diffdice==5))/N;
    fprintf('Problem 3.20, Part a (assuming difference can be negative): \n\tP(x=-5)=%f\n\tP(x=-4)=%f\n\tP(x=-3)=%f\n\tP(x=-2)=%f\n\tP(x=-1)=%f\n\tP(x=0)=%f\n\tP(x=1)=%f\n\tP(x=2)=%f\n\tP(x=3)=%f\n\tP(x=4)=%f\n\tP(x=5)=%f\n\n', [prob_neg5, prob_neg4, prob_neg3, prob_neg2, prob_neg1, prob_0, prob_1, prob_2, prob_3, prob_4, prob_5])
    %Plot pmf for X
    figure('Name','pmf X, Prob 3.20 (a)','NumberTitle','off');
    prob_vect = [prob_neg5 prob_neg4 prob_neg3 prob_neg2 prob_neg1 prob_0 prob_1 prob_2 prob_3 prob_4 prob_5];
    bar(-5:5, prob_vect);
    axis([-5.25, 5.25, -.25, .5]);
%Part b
    %Calculate pmf for |X| <= k
    diffdice = abs(ceil(6*rand(1,N)) - ceil(6*rand(1,N))); %Calculate the abs diff of two dice tosses
    prob_0 = length(find(diffdice<=0))/N;
    prob_1 = length(find(diffdice<=1))/N;
    prob_2 = length(find(diffdice<=2))/N;
    prob_3 = length(find(diffdice<=3))/N;
    prob_4 = length(find(diffdice<=4))/N;
    prob_5 = length(find(diffdice<=5))/N;
    fprintf('Problem 3.20, Part b: \n\tP(|X|<=0)=%f\n\tP(|X|<=1)=%f\n\tP(|X|<=2)=%f\n\tP(|X|<=3)=%f\n\tP(|X|<=4)=%f\n\tP(|X|<=5)=%f\n\n', [prob_0, prob_1, prob_2, prob_3, prob_4, prob_5])

%% Problem 3.21
%Part a
    %Calculate probabilities for Y
    Y = sum(round(rand(2, N)));%number of heads for two tosses
    ExpectedY = mean(Y);%Calculate the expected value of X
    %Calculate probabilities for X
    Carlos = sum(round(rand(2, N)));%number of heads for Carlos
    Michael = sum(round(rand(2, N)));%number of heads for Michael
    X = max([Carlos;Michael]);%Chooses the max number of heads from each
    ExpectedX = mean(X);%Calculate the expected value of X
    fprintf('Problem 3.21, Part a: \n\tE[Y]=%f\n\tE[X]=%f\n\n',[ExpectedY,ExpectedX]);
%Part b
    %Calculate the variance for X and Y
    VarX = mean(X.^2) - ExpectedX^2;% Variance of X, also VarX = var(X)
    VarY = mean(Y.^2) - ExpectedY^2;% Variance of Y, also VarY = var(Y)
    fprintf('Problem 3.21, Part b: \n\tVAR[X]=%f\n\tVAR[Y]=%f\n\n',[VarX,VarY]);