data nullpop;
 alpha=.05;
 tails=2;
 mu0=6;
 mu1=5;
 sigma2=2.5;
 n=10;
 xbarzc2= mu0 + probit(1-alpha/tails) * sqrt(sigma2/n);
 xbarzc1= mu0 + probit(alpha/tails)   * sqrt(sigma2/n);
 xbartc2= mu0 + tinv((1-alpha/tails),(n-1)) * sqrt(sigma2/n);
 xbartc1= mu0 + tinv((alpha/tails),(n-1))   * sqrt(sigma2/n);
 tal= (xbartc1-mu1) / sqrt(sigma2/n);
 tau= (xbartc2-mu1) / sqrt(sigma2/n);
 beta=probt(tau,n-1) - probt(tal,n-1);
 power=1-beta;
proc print;
run;

data one;
set nullpop;
samples=10000;
do sample=1 to samples;
  do mu=mu0, mu1;
    sigmax=0;
    do i=1 to n;
      array r {*} r1-r50;
      x=mu + sqrt(sigma2)*normal(1);
      r{i}=x;
      sigmax+x;
      xbar=sigmax/n;
      z=(xbar-mu)/(sqrt(sigma2/n));
      ss=0;
      if i=n then do s=1 to n;
        d=(r{s}-xbar)**2;
        ss+d;
        var=ss/(n-1);
        t=(xbar-mu)/(sqrt(var/n));
        drop alpha tails samples xbarzc1 xbarzc2 xbartc1 xbartc2
             x i s ss d sigmax sigma2 n r1-r50;
        xbar=round(xbar,.01);
        z=round(z,.01);
        t=round(t,.01);
        if s=n then output;
      end;
    end;
  end;
end;

proc sort data=one;
  by mu;
proc freq noprint data=one;
  by mu;
  tables xbar / out=two;
  *tables z / out=two;
  *tables t / out=two;
run;

proc gplot data=two;
  where mu=6;
  *symbol v=plus i=sm50s w=3;
  symbol i=sm50s w=3;
  plot percent*xbar=mu
        / vaxis=0 to 1 by .1 vminor=0
          haxis=4 to 10 by 1 hminor=10
          href=7.13;
run;