% LMFig060203a_4th.
% based almost entirely on 
% LMFig060402_4th.m
% based almost exactly on 
% 90 min; just change figure642=0 to plot LM Fig 6.2.3a
% LMFig060203a_4th.m
% Larsen & Marx (2006, p. 448).
% Normal probability pdf:
% The normal probability equation is provided on p. 264.
% This is for mean 0, and unit standard
% deviation. The more general equation (Legendre & Legendre, 1998 p. 147) is:
% f(y_j)=1/(sqrt(2*pi)*sigma_j)*exp(-1/2*((y_j-mu_j)/sigma_j)^2)
% Written by Eugene.Gallagher@umb.edu
% written for EEOS601
% find the gas mileage, yup, that produces a p value of 0.05
% Use Matlab's equation solver and the anonymous function
n=30;sigma=2.4;mu=25;alpha=0.05
figure642=0;
figure643=0;
onetailed=1;
K = menu('Choose a graph','Figure 6.2.3a, 1-tailed','Figure 6.2.3a, 2-tailed','Figure 6.4.2','Figure 6.4.3');
if K==2
    onetailed=0;
elseif K==3
    figure642=1;
elseif K==4
    figure643=1;
end
if onetailed==1;
    % This is one tailed:
    yup=norminv(1-alpha,mu,sigma/sqrt(n))
    % or for those who'd like to try the optimization toolbox:
    yup2 = fsolve(@(x) (1-alpha)-normcdf((x-mu)./(sigma./sqrt(n))),25,optimset('Display','off'))
    yup-yup2  % <e-9
    fprintf('The critical value is %6.4f\n',yup)
else
    % This is two-tailed you must solve for both tails
    yup=norminv(1-alpha/2,mu,sigma/sqrt(n))
    yup3 = fsolve(@(x) (1-alpha/2)-normcdf((x-mu)./(sigma./sqrt(n))),25,optimset('Display','off'))
    yup-yup3
    ylow=norminv(alpha/2,mu,sigma/sqrt(n))
    ylow3 = fsolve(@(x) (alpha/2)-normcdf((x-mu)./(sigma./sqrt(n))),25,optimset('Display','off'))
    ylow-ylow3
    fprintf('The lower critical value is %6.4f\n',ylow)
    fprintf('The upper critical value is %6.4f\n',yup)
end
p=1-normcdf((yup-mu)/(sigma/sqrt(n)))
muj=mu;
sigmaj=sigma/sqrt(n); % sigmaj is the standard error of the mean
miny=22.9;
if figure643==1
    maxy=28.1;
else
    maxy=27.1;
end
y=miny:.01:maxy;
fyj=normpdf(y,mu,sigmaj); 
% Plot using ax1 handle, saved above,to save this graph
% on top of the previous graph.
plot(y,fyj,'linestyle','--','color','r','linewidth',3)
ylabel('Density','FontSize',20)
xlabel('Gas Mileage','FontSize',22)
axis([miny maxy 0 1.05])
set(gca,'Ytick',[0 0.5 1],'FontSize',18)
set(gca,'Xtick',[23:.5:25 yup 26.5:0.5:28],'FontSize',18)
ax1=gca;% save the handle of the graph
if figure642==1
    title('Figure 6.4.2','FontSize',22)
elseif figure643==1
    title('Figure 6.4.3','FontSize',22)
else
    title('Figure 6.2.3a','FontSize',22)
end
hold on
fyj1=normpdf(yup,mu,sigmaj);
plot([yup yup]',[0 fyj1]','-k','linewidth',1) 
fymu=normpdf(mu,mu,sigmaj);
plot([mu mu]',[0 fymu]','-.r','linewidth',1) 
% Fill in the upper tail with fill
y2=yup:.1:maxy;
fyj2=normpdf(y2,mu,sigmaj);
fyjmax=normpdf(maxy,mu,sigmaj);
fill([yup y2 maxy maxy],[0 fyj2 fyjmax 0],[1 .4 .4])
t=sprintf('Mileage>%5.2f\nArea=%5.3f\n',yup,p);
text(26,.1,t,'Color','k','FontSize',16);
if onetailed ~=1
    % Fill in the lower tail
    y4=miny:.01:ylow;
    fyj4=normpdf(y4,mu,sigmaj);
    fyjmin=normpdf(miny,mu,sigmaj);
    fill([miny miny y4 ylow],[0 fyjmin fyj4 0],[1 .4 .4])
    t=sprintf('Mileage<%5.2f\nArea=%5.3f\n',ylow,p);
    text(miny+.2,.1,t,'Color','k','FontSize',16);
end
t6=sprintf('Sampling\ndistribution of Y\nwhen Ho true');
text(23.3,.9,t6,'Color','r','FontSize',20);
if figure642==1
    mu2=25.750;
    fyj3=normpdf(y,mu2,sigmaj); 
    % Plot using ax1 handle, saved above,to save this graph
    % on top of the previous graph.
    plot(y,fyj3,'linestyle','--','color','b','linewidth',3);
    y5=miny:.01:yup;
    fyj5=normpdf(y5,mu2,sigmaj);
    fyj5min=normpdf(miny,mu2,sigmaj);
    fyj5max=normpdf(yup, mu2,sigmaj);
    fill([miny miny y5 yup],[0 fyj5min fyj5 0],[.8 .8 1])
    t=sprintf('Sampling\ndistribution of Y\nwhen mu=%5.2f',mu2);
    text(26,.9,t,'Color','b','FontSize',20);
    fyj6=normpdf(yup,mu2,sigmaj);
    plot([yup yup]',[0 fyj6]','-b','linewidth',1) 
    h1 = findobj(gca,'Type','patch'); 
    set(h1,'facealpha',0.5);
    Beta=normcdf(yup,mu2,sigmaj);
    t5=sprintf('Beta=%6.4f\n',Beta);
    text(25.05,.05,t5,'Color','b','FontSize',16);
    fymu6=normpdf(mu2,mu2,sigmaj);
    plot([mu2 mu2]',[0 fymu6]','-.b','linewidth',1)
elseif figure643==1
    mu2=26.8;
    fyj3=normpdf(y,mu2,sigmaj); 
    % Plot using ax1 handle, saved above,to save this graph
    % on top of the previous graph.
    plot(y,fyj3,'linestyle','--','color','g','linewidth',3);
    y5=miny:.01:yup;
    fyj5=normpdf(y5,mu2,sigmaj);
    fyj5min=normpdf(miny,mu2,sigmaj);
    fyj5max=normpdf(yup, mu2,sigmaj);
    fill([miny miny y5 yup],[0 fyj5min fyj5 0],[.8 .8 1])
    t=sprintf('Sampling\ndistribution of Y\nwhen mu=%5.2f',mu2);
    text(27,.9,t,'Color','g','FontSize',20);
    fyj6=normpdf(yup,mu2,sigmaj);
    plot([yup yup]',[0 fyj6]','-g','linewidth',1) 
    h1 = findobj(gca,'Type','patch'); 
    set(h1,'facealpha',0.5);
    Beta=normcdf(yup,mu2,sigmaj);
    t5=sprintf('Beta=%6.4f\n',Beta);
    text(24.6,.05,t5,'Color','g','FontSize',16);
    fymu6=normpdf(mu2,mu2,sigmaj);
    plot([mu2 mu2]',[0 fymu6]','-.g','linewidth',1)
end
figure(gcf)
hold off