cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Name: Log2Power.for c c Consider two lognormal distributions: lognormal(u1, sigma1^2) and c lognormal(u2, sigma2^2). Let meani denote the mean of lognormal(ui,sigmai^2). c That is, meani = exp(ui+sigmai^2/2). For given sample sizes n1 and n2, and c parameters (u1,u2) and (sigma1^2,sigma2^2), this program computes the power c for testing c c H_0: mean1 <= mean2 vs H_a: mean1 > mean2 c c The test is based on a generalized varaible approach, and is given in c Krishnamoorthy, K. and Mathew, T. (2003). Inferences on the means of c lognormal distributions using generalized p-values and generalized confidence c intervals. To appear in Journal of Statistical Planning and Inference. c c Input: n1, n2, u1, u2, sigma1^2, sigma2^2 c alpha = level of the test c Output: Power c c To compute the power of a two-tail test, define p1 = proprtion of g_var > 0 c and p2 = proportion of g_var < 0. The p-value = 2 x min{p1, p2}. c Power = proprtion of the p-values < alpha. c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real n1, n2 print *,'Enter the values of n1 and n2:' read *, n1, n2 print *,' ' print *,'Enter the values of u1 and u2:' read *, u1, u2 print *,' ' print *,'Enter the values of sig1sq and sig2sq:' read *, sig1sq, sig2sq print *,' ' print *,'Enter the level of the test:' read *, 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 mm1 = 2500 mm2 = 5000 do i = 1, mm1 xb1 = u1 + anorran()*sigxb1 xb2 = u2 + anorran()*sigxb2 chi1 = 2.0*gamran(shape1) v1 = sig1**2*chi1/df1 s1 = sqrt(v1) chi2 = 2.0*gamran(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, mm2 z1 = anorran() chi1 = 2.0*gamran(shape1) v1 = chi1 z2 = anorran() chi2 = 2.0*gamran(shape2) v2 = chi2 g_var = (xb1 + z1*const11*sqrtdf1/sqrt(v1)+const12/v1) & - (xb2 + z2*const21*sqrtdf2/sqrt(v2)+const22/v2) if(g_var .lt. 0.0) prvalue = prvalue + 1.0 end do if(prvalue/mm2 .lt. alpha) power = power + 1.0 end do power = power/mm1 print *,'----------------------------------' print *, 'Power = ',power end c Fuction subroutine GAMRAN for gamma variate generator real function gamran(a) if(a .eq. 1.0) then gamran = -alog(grnd()) return end if x3 = a-1.0 d = sqrt(x3) ek =1.0 x1 = 0 x2 = 0 f2 = 0 if (d .ge. x3) go to 2 x2 = x3 -d ek = 1.0- x3/x2 x1 = x2 + 1.0/ek f2 = fgam(x2, a) 2 x4 = x3 + d r = 1.0 - x3/x4 x5 = x4 + 1.0/r f4 = fgam(x4, a) p1 = x4 - x2 p2 = p1 - f2/ek p3 = p2 + f4/r 3 u = grnd() v = grnd() u = u*p3 If (u .gt. p1) go to 4 x = x2 + u If(x .gt. x3 .and. v .le.f4 + (x4 - x)*(1.0 - f4)/(x4 - x3))goto 7 If(x .lt. x3 .and. v .le.f2 + (x - x2)*(1.0 - f2)/(x3 - x2))goto 7 go to 6 4 If (u .gt. p2) go to 5 u = (u - p1)/(p2 - p1) x = x2 - alog(u)/ek If (x .lt. 0.0) go to 3 v = v*f2*u If (v .le. f2*(x - x1)/(x2 - x1))go to 7 go to 6 5 u = (u - p2)/(p3 - p2) x = x4 - alog(u)/r v = v*f4*u If (v .le. f4*(x5 - x)/(x5 - x4))go to 7 6 If (alog(v) .le. x3*alog(x/x3) + x3 - x) then go to 7 else go to 3 end if 7 gamran = x end c The following subroutine FGAM is called by GAMRAN real function fgam(x, a) x3 = a - 1.0 fgam = exp(x3*alog(x/x3)+x3-x) end c Function subroutine for generating N(0,1) variates real function anorran() implicit double precision (a-h, o-z) real u(1), v(1), w(1) xi = 2.216035867166471d0 c u(1) = grnd() c if(u(1) .lt. 0.8840704022 98758d0) then v(1) = grnd() x = xi*(1.131131635444180*u(1)+v(1)-1.) go to 12 end if c if(u(1) .lt. 0.973310954173898d0) goto 4 2 v(1) = grnd() w(1) = grnd() t = xi**2/2.0-alog(w(1)) if(v(1)**2*t .gt. xi**2/2) goto 2 if(u(1) .lt.0.986655477086949d0) then x = sqrt(2.*t) goto 12 else x = -sqrt(2.*t) goto 12 end if 4 if(u(1) .lt. 0.958720824790463d0) goto 6 5 v(1) = grnd() w(1) = grnd() z = v(1)-w(1) t = xi - 0.630834801921960d0*min(v(1),w(1)) if(max(v(1),w(1)) .le. 0.75559131667601d0) goto 9 if(0.034240503750111d0*abs(z) .le. f(t)) goto 9 goto 5 6 if(u(1) .lt. 0.911312780288703d0) goto 8 7 v(1) = grnd() w(1) = grnd() z = v(1)-w(1) t = 0.479727404222441d0+1.105473661022070d0*min(v(1),w(1)) if(max(v(1),w(1)) .le. 0.872834976671790d0) goto 9 if(0.049264496373128d0*abs(z) .le. f(t)) goto 9 goto 7 8 v(1) = grnd() w(1) = grnd() z = v(1)-w(1) t = 0.479727404222441d0 -.595507138015940d0*min(v(1),w(1)) if(max(v(1),w(1)) .lt. 0.805577924423817d0) goto 9 goto 8 c 9 if(z .lt. 0.0) then x = t else x = -t end if c 12 anorran = x end double precision function f(x) implicit double precision (a-h, o-z) xi = 2.216035867166471d0 if(abs(x) .lt. xi) then f = exp(-x**2/2)/2.50663-0.1800251910*(xi-abs(x)) else f = 0.0 end if end *********************************************************************** * A C-program for MT19937: Real number version * genrand() generates one pseudorandom real number (double) * which is uniformly distributed on [0,1]-interval, for each * call. sgenrand(seed) set initial values to the working area * of 624 words. Before genrand(), sgenrand(seed) must be * called once. (seed is any 32-bit integer except for 0). * Integer generator is obtained by modifying two lines. * Coded by Takuji Nishimura, considering the suggestions by * Topher Cooper and Marc Rieffel in July-Aug. 1997. * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later * version. * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. * See the GNU Library General Public License for more details. * You should have received a copy of the GNU Library General * Public License along with this library; if not, write to the * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA * 02111-1307 USA * * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. * When you use this, send an email to: matumoto@math.keio.ac.jp * with an appropriate reference to your work. * ************************************************************************ * Fortran translation by Hiroshi Takano. Jan. 13, 1999. * * genrand() -> double precision function grnd() * sgenrand(seed) -> subroutine sgrnd(seed) * integer seed * * This program uses the following non-standard intrinsics. * ishft(i,n): If n>0, shifts bits in i by n positions to left. * If n<0, shifts bits in i by n positions to right. * iand (i,j): Performs logical AND on corresponding bits of i and j. * ior (i,j): Performs inclusive OR on corresponding bits of i and j. * ieor (i,j): Performs exclusive OR on corresponding bits of i and j. * ************************************************************************ ************************************************************************ subroutine sgrnd(seed) * implicit integer(a-z) * * Period parameters parameter(N = 624) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ * * setting initial seeds to mt[N] using * the generator Line 25 of Table 1 in * [KNUTH 1981, The Art of Computer Programming * Vol. 2 (2nd Ed.), pp102] * mt(0)= iand(seed,-1) do 1000 mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) 1000 continue * return end ************************************************************************ double precision function grnd() * implicit integer(a-z) * * Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) * constant vector a parameter(UMASK = -2147483648) * most significant w-r bits parameter(LMASK = 2147483647) * least significant r bits * Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ data mti/N1/ * mti==N+1 means mt[N] is not initialized * dimension mag01(0:1) data mag01/0, MATA/ save mag01 * mag01(x) = x * MATA for x=0,1 * TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) * if(mti.ge.N) then * generate N words at one time if(mti.eq.N+1) then * if sgrnd() has not been called, call sgrnd(4357) * a default initial seed is used endif * do 1000 kk=0,N-M-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 1000 continue do 1100 kk=N-M,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 1100 continue y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif * y=mt(mti) mti=mti+1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) * if(y.lt.0) then grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0) else grnd=dble(y)/(2.0d0**32-1.0d0) endif * return end