
 disp('Lead Exposure: Variances')
 lw = 3; 
 set(0, 'DefaultAxesFontSize', 16);
 fs = 15;
 msize = 10;



% In the context of Example Lead Exposure.
% Approximate the power of the two sided, 5\% level test of equality of variances
% against the specific alternative $H_1: \sigma_1^2/\sigma^2_2 = 1.8.$ 
% Recall, $n_1 = 12$ and $n_2 = 15.$ 
% According to (\ref{eq:pow2var}),

  n1 = 12;  n2=15; alpha = 0.05; F = 1.8;
  normcdf(   sqrt(2 * (n1-1)*(n2 -2)/(n1 + n2 -2  )) ...
              * abs(log(F)) - norminv(1-alpha/2) )
% 0.5112
% What is the retrospective power of test of variances in Example  .   
%   
    s12 = 0.004^2;  s22 = 0.006^2;
    F=s12^2/s22^2; %0.1975
    normcdf(   sqrt(2 * (n1-1)*(n2 -2)/(n1 + n2 -2  )) ...
             * abs(log(F)) - norminv(1-0.05/2) )
%  0.9998
%  What sample size would guerantee  that the effect of 
%  size 1.5 would be detected with the power of 90\% in a two sided 
%  5\%  level test.
%   
  alpha = 0.05; beta = 0.1; eff = 1.5;
  n = ( (norminv(1-alpha/2) + norminv(1-beta))/log(eff))^2 + 2
 % n=66 (65.9130)
  

%%
 alpha = 0.05; n=50;
 MM = 10000; sigma  = 1; mean = 0;
 count1 = [];
 count2 = [];
 for i = 1:MM
     sample = normrnd(mean, sigma, [n 1]);
     s2 = var(sample);
LB1 = (n-1) * s2/chi2inv(1-alpha/2, n-1);
UB1 = (n-1) * s2/chi2inv(alpha/2, n-1);
LB2 = s2 - norminv(1-alpha/2)*sqrt(2)* s2/sqrt(n-1);
UB2 = s2 + norminv(1-alpha/2)*sqrt(2)* s2/sqrt(n-1);


if (LB1 > 1 | UB1 < 1)  
    add1 = 0;
else add1 = 1;
end
if (LB2 > 1 | UB2 < 1) 
    add2 = 0;
else add2 = 1;
end
count1 =[count1 add1];
count2 =[count2 add2];
 end
 
 sum(count1)/MM
 
 sum(count2)/MM
