close all; clear all;

%% Probability Simulation

n = 1e3;
Number_Dice = 2;
D = randi([1,6],Number_Dice,n);
S = sum(D); % sum along each column

RF11 = cumsum(S==11)./(1:n);
plot(1:n,RF11)
hold on
RF12 = cumsum(S==12)./(1:n);
plot(1:n,RF12,'r')

xlabel('Number of Rolls')
ylabel('Relative Frequency')
legend('Sum = 11', 'Sum = 12')
grid on

figure
S_Support = (1*Number_Dice):(6*Number_Dice);
hist(S,S_Support)

N_S_Sim = hist(S,S_Support);


%% Probability Calculation

Dice_Support = 1:6;
S = [];
for k1 = Dice_Support
    for k2 = Dice_Support
            S = [S k1+k2];
    end
end
Size_SampleSpace = length(S);

Number_11 = sum(S==11)
Number_12 = sum(S==12)

% Count all possible cases at once
N_S = hist(S,S_Support) 

P = sym(N_S)/Size_SampleSpace

figure
stem(S_Support,P)

%%
hold on
RF = N_S_Sim./n;
plot(S_Support,RF,'rx')
legend('probability','relative frequency')
xlabel('sum of two dice')
grid on

