ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Name: Log_Two.for c This program compute the p-value for testing the difference between means of c two lognormal distributions. The hypotheses are c c H_0: mean1 <= mean2 vs H_a: mean1 > mean2 c c Note: Consider two lognormal populations with parameters (u1, sigma1^2) and c (u2,sigma2^2). Then mean1 = exp(u1+sigma1^2/2) and mean2 = exp(u2+sigma2^2/2) c Let xbi = mean of the logged data from the ith population, i=1,2. c si_sq = variance of the logged data from the ith population, i=1,2. c For a given sample sizes n1, n2, xb1, xb2, s1_sq and s2_sq, this program c computes the p-value for testing the hypotheses given above and confidence c interval for the difference mean1 - mean2. c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real n1, n2, g_var(100000) print *,'Enter the sample sizes: ' read *, n1, n2 print *, ' ' print *,'Enter the sample means of the logged data: ' read *, xb1, xb2 print *, ' ' print *,'Enter the sample varainces of the logged data: ' read *, s1_sq, s2_sq print *, ' ' print *,'Enter the confidence level: ' read *, cl print *, ' ' alpha = 1.0-cl mm = 100000 L1 = int(alpha*mm/2.0) L2 = int((1.0-alpha/2.0)*mm) df1 = n1 - 1.0 sqrtdf1 = sqrt(df1) df2 = n2 - 1.0 sqrtdf2 = sqrt(df2) s1 = sqrt(s1_sq) s2 = sqrt(s2_sq) const11 = s1/sqrt(n1) const12 = s1_sq*df1*0.5 const21 = s2/sqrt(n2) const22 = s2_sq*df2*0.5 pvalue = 0.0 iseed = 123457 do j = 1, mm z1 = anorran() v1 = 2.0*gamran(df1*0.5) z2 = anorran() v2 = 2.0*gamran(df2*0.5) T_31 = xb1 + z1*const11*sqrtdf1/sqrt(v1)+const12/v1 T_32 = xb2 + z2*const21*sqrtdf2/sqrt(v2)+const22/v2 g_var(j) = exp(T_31) - exp(T_32) if (g_var(j) .lt. 0.0) pvalue = pvalue + 1.0 end do pvalue = pvalue/mm call hpsort(mm, g_var) print *,'-----------------------------------------------------' print *,'p-value for testing mean_1 <= mean2 vs & mean1 > mean2 is ',pvalue print *,' ' print *,'-----------------------------------------------------' print *,cl,'Confidence Interval for mean1 - mean2' print *,' ' print *,'(',g_var(L1),', ',g_var(L2),')' end c Fuction subroutine GAMRAN for gamma varaite 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 SUBROUTINE HPSORT(N,RA) real RA(N) L=N/2+1 IR=N c !The index L will be decremented from its initial value during the c !"hiring" (heap creation) phase. Once it reaches 1, the index IR c !will be decremented from its initial value down to 1 during the c !"retirement-and-promotion" (heap selection) phase. 10 continue if(L > 1)then L=L-1 RRA=RA(L) else RRA=RA(IR) RA(IR)=RA(1) IR=IR-1 if(IR.eq.1)then RA(1)=RRA return end if end if I=L J=L+L 20 if(J.le.IR)then if(J < IR)then if(RA(J) < RA(J+1)) J=J+1 end if if(RRA < RA(J))then RA(I)=RA(J) I=J; J=J+J else J=IR+1 end if goto 20 end if RA(I)=RRA goto 10 END