cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Program: Log_Power.for c c This program computes the power of the generalized c test for testing the mean of a lognormal distribution with parameter c u_l and sigma_l; u_l = population mean of the logged data and c sigma_l = population standard deviation of the logged data c For a given sample size n, u_l and sigma_l, this program computes the c power for testing c c H_0: M >= M_0 vs H_a: M < M_0 (*) c c where M = u_l + sigma_l^2/2 and M_0 is a specified value of M. c Note that exp(M) is the mean of the lognormal distribution. Testing c (*) is equivalent to testing c c H_0: exp(M) >= exp(M_0) vs H_a: exp(M) < exp(M_0) c c INPUT: u_l, sigma_l, mean_0 = M_0 and n c OUTPUT: power c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real n, mean, mean_0 read *, u_l, sigma_l, mean_0, n c nominal level of the test alpha = 0.05 iseed = 123457 df = n-1.0 mm1 = 2500 do i = 1, mm1 var = sigma_l**2*2.0*gamran(df/2.0)/df s = sqrt(var) xbar = u_l + anorran()*sigma_l/sqrt(1.0*n) const1 = s*sqrt(df)/sqrt(n) const2 = var*df*0.5 pvalue = 0.0 mm2 = 5000 do j = 1, mm2 z = anorran() v = 2.0*gamran(df/2.0) stat = xbar + z*const1/sqrt(v)+const2/v if(stat .gt. mean_0) pvalue = pvalue + 1.0 end do pvalue = pvalue/mm2 if(pvalue .lt. alpha) power = power + 1.0 end do power = power/mm1 print *, power end c Function subroutine to generate random numbers from a gamma c distribution with shape parameter a 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 c The function f(x) is called by anorran() 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 SUBROUTINE SHELL(N,ARR) parameter(ALN2I=1./0.69314718,TINY=1.E-5) real ARR(N) LOGNB2=INT(ALOG(FLOAT(N))*ALN2I+TINY) m=n do nn=1,LOGNB2 m=m/2 k=n-m do j=1,k i=j 10 continue l=i+m if(ARR(l).LT.ARR(i)) then t=ARR(i) ARR(i)=ARR(l) ARR(l)=t i=i-m if(i.GE.1) GOTO 10 end if end do end do return END