%
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

xbar    = 47.2;
mu_0    = 45; 
sigma_0 = 15; 

figure; grid on; hold on; 

nArray = [ 50, 100, 200 ]; allPs = []; 
for ni=1:length(nArray),
  % compute the power (1-probability of type II error) in each case:
  % -- this is simple enough that we'll just code this in the loop here ... 
  %
  
  % our alternative H_1 hypothesis:
  mualt = 30:60;

  % our "critical value" used in determining whether to reject H_0 or not:
  cv  = 1.645; 
  % in terms of physical variables of H_0 this translates into 
  sig = sigma_0/sqrt(nArray(ni)); 
  pcv = mu_0 + cv*sig;

  % the probability of a type II error (as a function of \mu) is then: 
  pII = normcdf( pcv, mualt, sig );
  allPs = [allPs, pII(:)]; 
end

hs=plot( mualt, 1-allPs, '-x' ); axis tight; 
legend( hs, { 'n=50', 'n=100', 'n=200' } ); 
title( 'power' ); xlabel( '\mu' ); ylabel( 'probability' ); 
saveas( gcf, 'ex_6_4_p_type_II_error_fn_n', 'epsc' );