% LMex110302_4th.
% See LMex110303_4th.m for a more complete analysis, with prediction
% intervals
% Example 11.3.2 pp. 694 in
% Larsen & Marx (2006)Introduction to Mathematical Statisticcs, 4th edition
% Written by Eugene D. Gallagher, Eugene.Gallagher@umb.edu
% Written November 2001, updated 2/24/2011
% 95% Confidence intervals for a new value of X
% DATA from case study 11.3.2.
DATA=...
[1765 407;1942 466;2132 489;2431 545;2642 610;2895 659;2931 686;3217 724];
% solve the regression equation
[r,c]=size(DATA);
x=[repmat(1,r,1) DATA(:,1)]; % Create a leading column of all 1's
y=DATA(:,2);
[B,BINT,R,RINT,STATS] = REGRESS(y,x);
fprintf('From regress.m, the Y intercept & slope =[%5.2f %4.2f]\n',B);
% This will duplicate case study 11.3.2
% I need Matlab's polynomial fitting algorithms to easily fit intervals
% resolve the regression using polyfit
[p,S] = polyfit(DATA(:,1),y,1);
fprintf('From polyfit.m, the Y intercept & slope =[%5.2f %4.2f]\n',fliplr(p));
% Calculate the 95% confidence interval for a value of 3500
X=3500;
[Y1,DELTA1] = polyconf(p,X,S); % prediction intervals for observations;
[Y2,DELTA2] = polyconf(p,X,S,'predopt','curve'); % confidence intervals
fprintf('The expected expenses for $%7.0f sales = $%6.0f\n',X*1000,Y1*1000);
fprintf('The 95%% Confidence Interval is $%6.0f to $%6.0f\n',...
    (Y2-DELTA2)*1000,(Y2+DELTA2)*1000);
fprintf('The 95%% Prediction Interval is $%6.0f to $%6.0f\n',...
    (Y1-DELTA1)*1000,(Y1+DELTA1)*1000);
fprintf(...
'The prediction interval is %2.0f%% wider than the confidence interval.\n',...
DELTA1/DELTA2*100-100) 
% b[Y3,DELTA3]= polyconf(p,X,S,'simopt','on'); % Hotelling-Working interval                                   
                                           % for observations
% [Y4,DELTA4]= polyconf(p,X,S,'simopt','on','predopt','curve'); % Hotelling-
            % Working interval for mean
            
% Plot figure 11.3.4   
X=[1760:10:3220]';
[Y5,DELTA5] = polyconf(p,X,S,'predopt','curve'); % confidence intervals
LL=Y5-DELTA5;UL=Y5+DELTA5;
% [Y6,DELTA6] = polyconf(p,X,S); % prediction intervals
% LLi=Y6-DELTA6;ULi=Y6+DELTA6;
yl=sprintf('Expenses (thousands)');
plot(X,Y5,'m',X,LL,'--b',X,UL,'--b',...
                       DATA(:,1),y,'o','LineWidth',2,...
                       'MarkerEdgeColor','k',...
                       'MarkerFaceColor','b',...
                       'MarkerSize',8)
                   axis([1450 3550 340 790])
                   ax1=gca;
set(ax1,'xtick',1500:500:3500,'FontSize',18)
set(ax1,'xticklabel',{'$1500','$2000', '$2500','$3000','$3500'},'FontSize',18)
set(ax1,'ytick',350:100:750,'FontSize',18)
set(ax1,'yticklabel',{'$350', '$450','$550','$650','$750'},...
      'FontSize',18)
xlabel('Sales (thousands)','Color','black','FontSize',20)
ylabel(yl,'Color','black','FontSize',20)
title('Figure 11.3.4','FontSize',22)                   
figure(gcf)
pause
