%% 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))