% LMcs020402_4th.m
% NOT WORKING - just try montyhall.m from user-contributed m.files
% Case Study 2.4.2 p 50 from 
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Written by E. Gallagher & revised 1/14/11
% Written by Eugene.Gallagher@umb.edu
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% Monty Hall problem
% Two goats (G) and a car (C) are hidden behind 3 doors, producing the
% following sample space S
S=['CGG';'GCG';'GGC'];
% You the contestant have just selected door 2;
% Monty hall will then randomly select either door 1 or 3 to reveal a goat
% Let us state that Monty Hall chose door 3.
% Your choice is whether to switch or stay with your original wager.
% Solve using Bayes theorem:
% P winning if one stays with door 2=P(Car behind door 2|Monty Hall chose
% door 3)=P(A|B), where A=car IS behind door 2; and B=Monty Hall chose door
% 3;
% P(A|B)= P(B|A)*P(A)/[P(B|A)*P(A)+P(B|A)*(P(~A))]
%         P(B|A)=1/2;
PBgivenA=.5
%         P(A)=1/3;
PA=1/3;
%      P(~A)=1-P(A)=2/3;
PnA=2/3;
% So the probability of winning if one didn't switch is
PAgivenB=1/2*1/3/(1/2*1/3+1/2*2/3)
% So the probability of winning if one switches MUST be 2/3, but we can
% work that out.
%--------------------

% Brute enumeration
S=zeros(3^Doors,3); 
i=1;
for Actual=1:3
    for ContestantGuess=1:3
        for MontyOpen=1:3
            S(i,:)=[Actual ContestantGuess MontyOpen];
            i=i+1;
        end
    end
end
disp('The sample space, S:')
S
% But some of these outcomes are incompatible
i=find(S(:,1)==S(:,3)) % Monty wouldn't open a door with a car;
j=find(S(:,2)==S(:,3)) % Monty wouldn't open the door picked by a contestant
k=union(i,j);
S(k,:)=[]
% Assume contestant picks door 2;
[r,c]=size(S);
% Event A=the prize is behind door 1; B is event contestant picked door 2;
% P(A|B)= P(B|A)*P(A)/[P(B|A)*P(A)+P(B|~A)*(P(~A))]
% find PA
i=find(S(:,1)==1);
PA=length(i)/r;
% P(B|A)
j=find(S(i,2)==2);
PBgivenA=length(j)/length(i);
% P(B|~A)
i=find(S(:,1)~=1);
j=find(S(i,2)==2);
PBgivennotA=length(j)/length(i);
PnotA=1-PA;
PAgivenB=PBgivenA*PA/(PBgivenA*PA+PBgivennotA*PnotA)

        