% Monty_Crawl_Sol computes the winning probabilities in the Monty Crawl 
% game. Three player strategies are considered.

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_Smart_S = zeros(1,N);
X_Smart_NS = zeros(1,N);
H_ChooseL1 = zeros(1,N);
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 orginal choice of the player was correct,...
        X_H(k) = A(1); % Host always chooses the lower number door that is left
    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 Player's original choice
    if X_H(k)==A(2) % Host choose the higher number. We then know that the car must be in the lower number.
        % H_ChooseL1(k) = 0; % not needed
        X_Smart_S(k) = A(1);
        X_Smart_NS(k) = A(1);
        X_S(k) = A(1);
    else
        H_ChooseL1(k) = 1;
        A = A(A~=X_H(k)); % Remove the host choice
        X_Smart_NS(k) = X1(k);
        X_Smart_S(k) = A; % Switch to the one remaining door
        X_S(k) = A;
    end
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)

SmartNotSwitchCorrect = (car == X_Smart_NS);
p_smartnotswitch = cumsum(SmartNotSwitchCorrect)./(1:N);
plot(1:N,p_smartnotswitch,'g','LineWidth',1.5)

legend('(i) Switch','(ii) Not Switch','(iii) Smart, Not Switch')
ylabel('Relative Frequency of Winning')
xlabel('Number of Trials')
grid on
title('Winning Probabilities for the Monty Crawl Game')