%% Power Calculations

% Power is P(reject H0 | H1 true)
% Given H1 true and a normally distributed population with true mean mu1
% and null hypothesis mean mu0

mu0 = 8.4 % from BLS data
mu1 = 7.4 % One hour different


s = 1.3 % standard deviation (from our example)
n = 7 % our sample size

% noncentrality parameter
ncp = (mu1-mu0)/(s/sqrt(n))
df = n-1; % degrees of freedom

% Plot noncentral t
x = linspace(-6,6,100);
nct = nctpdf(x,df,ncp);
fig1 = figure; plot(x,nct); 
title(sprintf('Noncentral t distribution, %i degrees of freedom, ncp=%g',n-1,ncp));

alpha = 0.05
% 
fig2 = figure;
x2 = linspace(-3,3,100);
tprobdf = tpdf(x2,df);
tcumdf = tcdf(x2,df);
plot(x2,tprobdf,x2,tcumdf);
legend('Central t PDF, df=6','Central t CDF, df=6','Location','West');
title('Probability and Cumulative Distribution Functions of the central t distribution with 6 df')
pause

%% Find Critical T values by reading off central t CDF at alpha/2 and 1-alpha/2 percentiles
figure(fig2) %make sure the appropriate figure is current
Tcrit1 = tinv(alpha/2,df);
line([-3,Tcrit1],[alpha/2,alpha/2],'Color','r')
Tcrit2 = tinv(1-alpha/2,df);
line([-3,Tcrit2],[1-alpha/2,1-alpha/2],'Color','r')
pause

%% Plot Critical T values on axis
figure(fig2)
line([Tcrit2,Tcrit2],[0,1-alpha/2],'Color','r')
text(Tcrit2,.5,[sprintf('Critical T value = %0.3g',Tcrit2),'\rightarrow'],'HorizontalAlignment','right')
line([Tcrit1,Tcrit1],[0,0.1],'Color','r')
text(Tcrit1,.1,['\leftarrow ',sprintf('Critical T value = %0.3g',Tcrit1)],'HorizontalAlignment','left')

pause

%% Shade Rejection Region
figure(fig2)

shade_region([-1,1],[Tcrit1,Tcrit2],x2,tprobdf,tpdf([Tcrit1,Tcrit2],df));
pause

%% Back to the noncentral figure
figure(fig1);
pause
% Mark ame critical values
figure(fig1);
line([Tcrit2,Tcrit2],[0,0.3],'Color','r')
text(Tcrit2,.2,[sprintf('Critical T value = %0.3g',Tcrit2),'\rightarrow'],'HorizontalAlignment','right')
line([Tcrit1,Tcrit1],[0,0.1],'Color','r')
text(Tcrit1,.1,['\leftarrow ',sprintf('Critical T value = %0.3g',Tcrit1)],'HorizontalAlignment','left')

shade_region([-1,1],[Tcrit1,Tcrit2],x,nct,nctpdf([Tcrit1,Tcrit2],df,ncp));
% Area of the shaded region is the power. Outside the critical T values, we
% will reject H0.  Under H1, that was the right thing to do!
% Between the critical values, we will not reject H0.  Under H1, that was
% the wrong thing to do.  Not rejecting a false H0 is a Type II error.  The
% area in this region = probability of a Type II error, which is designated
% beta. 
pause

% plot CDF over this
hold on
plot(x,nctcdf(x,df,ncp),'m')
pause
%% Find percentiles to calculate area for power
% lines up from critical values to CDF, then over to Y-axis 
line([Tcrit1,Tcrit1],[0,nctcdf(Tcrit1,df,ncp)],'Color','r');
line([Tcrit2,Tcrit2],[0,nctcdf(Tcrit2,df,ncp)],'Color','r');
pause

% lines from CDF over to Y-axis, with labels
line([Tcrit1,min(x)],[nctcdf(Tcrit1,df,ncp),nctcdf(Tcrit1,df,ncp)],'Color','r');
text((Tcrit1+min(x))/2,nctcdf(Tcrit1,df,ncp),sprintf('CDF value = %0.3g',nctcdf(Tcrit1,df,ncp)));
line([Tcrit2,min(x)],[nctcdf(Tcrit2,df,ncp),nctcdf(Tcrit2,df,ncp)],'Color','r');
text((Tcrit2+min(x))/2,nctcdf(Tcrit2,df,ncp),sprintf('CDF value = %0.5g',nctcdf(Tcrit2,df,ncp)));


power = nctcdf(tinv(alpha/2,df),df,ncp)+(1-nctcdf(tinv(1-alpha/2,df),df,ncp))