% LMex030203_4th
% Probability of winning a series in 4, 5, 6 or 7 games
% if the probability of winning each game is 0.5
% Nahin's Dueling idiots also solves the problem 
% Written by E. Gallagher, revised 10/19/10
% Eugene.Gallagher@umb.edu
% Dept. of Environmental, Earth & Ocean Sciences
% Reference: Larsen & Marx (2006) Introduction to Mathematical
% Statistics, 4th edition, page 133-134.

% Call Matlab's stats toolbox binomial pdf function:
% p4 is the probability of 1 team winning 3 of 3 games * p winning a 
% 4th game; the total probability of a series ending is the probabilty
% of the 2nd team winning too, so since p=0.5, just multiply by 2
% p4=binopdf(4,4,0.5)*2 or 
p4=binopdf(3,3,0.5)*0.5; p4=2*p4
p5=binopdf(3,4,0.5)*0.5; p5=2*p5
p6=binopdf(3,5,0.5)*0.5; p6=2*p6
p7=binopdf(3,6,0.5)*0.5; p7=2*p7

% These 4 calculations can be handled in 1 Matlab line
% pdfseries=2*binopdf(3,3:6,0.5)*0.5
pdfseries=binopdf(3,3:6,0.5)
 % The observed and theoretical probabilities are shown in Table 3.2.3
 SL=[4:7]';
 ObservedP=[19/59 15/59 15/59 10/59];
 Table030203=[SL ObservedP' pdfseries'];
 fprintf('Table 3.2.3\n\n')
 fprintf('Series  Observed    Theoretical\n')
 fprintf('Length  Proportion  Probability\n')
 % Note that a loop structure can be avoided by using fprintf with
 % a matrix. The fprintf comman reads entries row-wise, hence the transpose
 fprintf('  %1.0f       %5.3f        %5.3f\n',Table030203')
 fprintf('Series  Observed    Expected\n')

% Optional produce a bar plot, with each bar of width=1;
Bardata=[ObservedP' pdfseries'];
bar(4:7,Bardata,1,'grouped')
% ylim([0 23])
legend('Observed P','Theoretical P (Pwin=0.5)',...
     'Location','NorthEast','FontSize',18)
xlabel('Series Length','FontSize',16)
ylabel('Proportion of Series','FontSize',16)
title('Example 3.2.3','FontSize',20);figure(gcf);
  
 
 % This next section uses some advance fitting algorithms to fit the model.
 % Optional fitting algorithms [Not required at all for EEOS601]
 % The observed proportions had more short games, what is the 
 % series length if the odds of winning differ from 0.5? More short
 % games would result if one team had a higher probability of winning.
 % This can be solved iteratively by plugging different values into
 % equations and visually examining the output. This was done and 0.75
 % seemed to provide a good fit.
 % Matlab has a suberb set of fitting tools, and fminbnd can be used
 % to find the value of pwin that provides the best fit. There are
 % several different ways of examining the best fit. One is minimizing
 % the sum of squared residuals sum((Observed-Expected).^2). The other,
 % more classic approach, is to minimize the chi-square statistic
 % Chi square = sum((Observed-Expected).^2./Expected). Both approaches
 % will be used with these Stanley cup data.
 
 pwin=0.75;
 % Can no longer multiply the prob of one team winning by 2
 pdfseries2=binopdf(3,3:6,pwin)*pwin + binopdf(3,3:6,1-pwin)*(1-pwin);
 FittedTable030203=[SL ObservedP' pdfseries2'];
 fprintf('\n\n\nFitted Table if pwin=%5.3f\n\n',pwin)
 fprintf('Series  Observed    Theoretical\n')
 fprintf('Length  Proportion  Probability\n')
 % Note that a loop structure can be avoided by using fprintf with
 % a matrix. The fprintf comman reads entries row-wise, hence the transpose
 fprintf('  %1.0f       %5.3f        %5.3f\n',FittedTable030203')

 % find the optimum value of pwin using Matlab's bounded search algorithm,
 % fminbnd
 % First, the minimization of squared residuals will be used, then 
 % using chi square stat: sum((O-E).^2./E)
 % I'm following the help file here and using an implicit function that 
 % I'll call f that will take as input the observed proportions and pwin
 % and return the sum of the squared residuals
 f=@(pwin,ObservedP) ...
 sum(  ( (binopdf(3,3:6,pwin)*   pwin   + ...
         binopdf(3,3:6,1-pwin)*(1-pwin))-ObservedP).^2);
         % Call the function f to find pwin that minimizes the function
         % with bounds 0.49 <= p <= 1
         pwino = fminbnd(@(pwin) f(pwin,ObservedP),0.49,1);
         % Note that a search, with 0,1 bounds returned 0.272
         % So the bounds should be set to find the p for the best team
         % with an allowance for a tie, with 0.49 as lower bound,
         % fminbind returns 0.728
 % Use this 'best' estimate to find the pdf        
 pdfseries3=binopdf(3,3:6,pwino)*pwino + binopdf(3,3:6,1-pwino)*(1-pwino);
 FittedTable=[SL ObservedP' pdfseries3'];
 fprintf('\n\n\nFitted Table if pwin=%5.3f\n\n',pwino)
 fprintf('Series  Observed    Theoretical\n')
 fprintf('Length  Proportion  Probability\n')
 fprintf('  %1.0f       %5.3f        %5.3f\n',FittedTable')
 
 % Now fit the problem minimizing the chi square statistic. Divide
 % each squared residual by the expected value. Observed & Expected values
 % should be in the form of counts, so although unimportant, I'll redefine
 % the function to use the observed & expected games (multiply by 59)
 f_cs=@(pwin,ObservedP) ...
 sum(  ((ObservedP-(binopdf(3,3:6,pwin)*   pwin   + ...
         binopdf(3,3:6,1-pwin)*(1-pwin)))*59).^2./ ...
        (binopdf(3,3:6,pwin)*   pwin   +     ...
         binopdf(3,3:6,1-pwin)*(1-pwin))*59);
         % 
         pwino_cs = fminbnd(@(pwin) f_cs(pwin,ObservedP),0.49,1);
 pdfseries3=binopdf(3,3:6,pwino_cs)*pwino_cs + ...
            binopdf(3,3:6,1-pwino_cs)*(1-pwino_cs);
        Games=59;
        CS=sum((Games*(ObservedP-pdfseries3)).^2./(Games*pdfseries3));
 % What is the probability of observing such a Chi Square statistic by
 % chance? This is covered in Chapter 10 (page 615 in the text). 
 % The degrees of freedom is equal to 1-the number of cells in 
 % the data, since if the total number of games played is known
 % and the number that ended in less than 7 games is known, 
 % then the number of series taking 7 games can
 % be calculated by difference. Only seriesnumber-1 degrees of freedom
 % are present. Also, an addtional degree of freedom is lost for every
 % parameter fit
 % with the data. In the chi-square fitting routine, I fit the best
 % estimate of the winning percentage, or P=0.7265, so another degree
 % of freedom is lost
 df=4-1-1;
 P = 1- chi2cdf(CS,df);
 FittedTable=[SL Games*ObservedP' Games*pdfseries3'];
 fprintf(...
     '\n\n\n Pwin=%6.4f, Chi Square=%5.2f, df=%1.0f, P=%4.2f\n\n',...
      pwino_cs, CS,df, P)
 fprintf('Series  Observed    Expected\n')
 fprintf('Length  Games       Games\n')
 % Note that a loop structure can be avoided by using fprintf with
 % a matrix. The fprintf comman reads entries row-wise, hence the transpose
 fprintf('  %1.0f    %5.0f         %4.1f\n',FittedTable')
 
 % What was the chi-square distribution and p value for the initial
 % fit which assumed that each team was equally likely to win? Note
 % that this fit requires only the loss of a single degree of freedom
 % since the 0.5 did not have to be estimated from the data.
CS2=sum((Games*(ObservedP-pdfseries)).^2./(Games*pdfseries));
df2=4-1;
 P2 = 1- chi2cdf(CS2,df2);
 FittedTable2=[SL Games*ObservedP' Games*pdfseries'];
 fprintf(...
     '\n\n\n Pwin=%6.4f, Chi Square=%5.2f, df=%1.0f, P=%6.2g\n\n',...
      0.5, CS2,df2,P2)

 fprintf('Series  Observed    Expected\n')
 fprintf('Length  Games       Games\n')
 fprintf('  %1.0f    %5.0f         %4.1f\n',FittedTable2')
Bardata=[FittedTable(:,2:3) FittedTable2(:,3)]; 
bar(4:7,Bardata,1,'grouped')
ylim([0 23])
legend('Observed','Expected (Pwin=0.7265)','Expected (Pwin=0.5)',...
     'Location','BestOutside','FontSize',18)
xlabel('Series Length','FontSize',16)
ylabel('Number of Series','FontSize',16)
title('Example 3.2.3','FontSize',20);figure(gcf);
