% Monty_Sol use simulation to estimate the winning probabilites in the 
% Monty Hall game.

clear all; close all;
D = 3; % number of doors
N = 1e3; % #times to repeat the game

% Initialization
car = randi(D,1,N); % The winning door
X1 = randi(D,1,N); % Player's original choice

OriginalCorrect = (car == X1);
% Preallocation
X_H = zeros(1,N); % Door open by the Host
X_S = zeros(1,N); % Remaining door that the player can switch to
X_NS = X1;
for k=1:N    
    OC = OriginalCorrect(k);
    % Host's selection
    % A is a temporary variable to keep the choice(s) under consideration 
    A = 1:D; % Start with all the choices for doors
    A = A(A~=X1(k)); % Remove the original choice
    A = A(A~=car(k)); % Remove the car position
    if OC % If original choice of the player was correct,...
        X_H(k) = A(randi(2)); % Host randomly chooses one of the two remaining doors
    else % If orginal choice of the player was incorrect,...
        X_H(k) = A; % Host chooses the one remaining door
    end
    % Player's selection
    A = 1:D; % Start with all the choices for doors
    % X_NS(k) = X1(k); % This line is not needed because we have already
    % equated X_NS and X1 earlier.
    A = A(A~=X1(k)); % Remove the original choice   
    A = A(A~=X_H(k)); % Remove the host choice
    X_S(k) = A; % Swicth to the one remaining door
end

SwitchCorrect = (car == X_S);
p_switch = cumsum(SwitchCorrect)./(1:N);
plot(1:N,p_switch,'r','LineWidth',1.5)

hold on

p_original = cumsum(OriginalCorrect)./(1:N);
plot(1:N,p_original,'LineWidth',1.5)

legend('(i) Switch','(ii) Not Switch')

ylabel('Relative Frequency of Winning')
xlabel('Number of Trials')
grid on
title('Winning Probabilities for the Monty Hall Game')