/***************************************************************************
 Name: logtwoPower.sas

 This program computes the power of testing two lognormal means.

 Consider two lognormal distributions: lognormal(u1, sigma1^2) and
 lognormal(u2, sigma2^2). Let meani denote the mean of lognormal(ui,sigmai^2).
 That is, meani = exp(ui+sigmai^2/2). For given sample sizes n1 and n2, and
 parameters (u1,u2) and (sigma1^2,sigma2^2), this program computes the power
 for testing

       H_0: mean1 <= mean2   vs   H_a: mean1 > mean2

 The test is based on a generalized varaible approach, and is given in
 Krishnamoorthy, K. and Mathew, T. (2003). Inferences on the means of
 lognormal distributions using generalized p-values and generalized confidence
 intervals. To appear in Journal of Statistical Planning and Inference.

 To compute the power of a two-tail test, define p1 = proprtion of g_var > 0
 and p2 = proportion of g_var < 0. The p-value = 2 x min{p1, p2}.

 Input:
                    n1, n2 = sample sizes
                    u1, u2 = means of two lognormal populations
     variance1, variance 2 = sigma1^1, sigma2^2
                     alpha = level of the test
                        m1 =
                        m2 =
 Output:
       Power = proprtion of the p-values < alpha.
*****************************************************************************/
%macro logtwoPower( n1=, n2=, mean1=, mean2=, variance1=,variance2 =,
                    alpha=0.05, m1=2500, m2=5000);

options ls = 64 ps = 45 nodate nonumber;
title 'Output of power of test of two lognormal means';
proc iml;

n1 = &n1; n2 = &n2;
u1 = &mean1; u2 = &mean2;
sig1sq = &variance1; sig2sq = &variance2;
alpha =&alpha;

sig1 = sqrt(sig1sq);
sig2 = sqrt(sig2sq);

df1 = n1 - 1.0;
sqrtdf1 = sqrt(df1);
df2 = n2 - 1.0;
sqrtdf2 = sqrt(df2);

sqrtn1 = sqrt(n1);
sigxb1 = sig1/sqrtn1;

sqrtn2 = sqrt(n2);
sigxb2 = sig2/sqrtn2;

shape1 = df1/2.0;
shape2 = df2/2.0;

seed1 = int(time());
seed2 = int(time()+12345);
seed3 = int(time()+13456);
seed4 = int(time()+19867);

reset name=noname ;
power = 0.0;
do i = 1 to m1;

   xb1 = u1 + rannor(seed1)*sigxb1;
   xb2 = u2 + rannor(seed2)*sigxb2;

   chi1 = 2.0*rangam(seed1,shape1);
   v1 = sig1**2*chi1/df1;
   s1 = sqrt(v1);

   chi2 = 2.0*rangam(seed2,shape2);
   v2 = sig2**2*chi2/df2;
   s2 = sqrt(v2);

   const11 = s1/sqrtn1;
   const12 = s1**2*df1*0.5;
   const21 = s2/sqrtn2;
   const22 = s2**2*df2*0.5;

   prvalue = 0.0;

      do j = 1 to m2;
         z1 = rannor(seed3);
         chi1 = 2.0*rangam(seed3,shape1);
         v1 = chi1;
         z2 = rannor(seed4);
         chi2 = 2.0*rangam(seed4,shape2);
         v2 = chi2;

         g_var = (xb1 + z1*const11*sqrtdf1/sqrt(v1)+const12/v1)
          - (xb2 + z2*const21*sqrtdf2/sqrt(v2)+const22/v2);
         if g_var < 0.0 then prvalue = prvalue + 1.0;

      end;
      if prvalue/m2 < alpha then power = power + 1.0;

end;
power = power/m1;

print 'n1=' n1[format=2.0] ',' 'u1=' u1[format=6.3] ','
   'sigma1^2=' sig1sq[format=6.3] ;
print 'n2=' n2[format=2.0] ',' 'u2=' u2[format=6.3] ','
   'sigma2^2=' sig2sq[format=6.3] ;
print 'The power for right-tail test is: ' power;

quit;

 %mend logtwoPower;