c NAME: logmean.for c ----------------- c This program computes confidence limits and p-values for testing c about the mean exp(mu + sigma^2/2) of a lognormal distribution c with parameter mu and sigma^2. The methods are based on the c generalized p-value and the generalized limit. Let x1,...,xn be a c sample from a lognormal population. c c Input: n = sample size c xbar = mean of ln(x1),...,ln(xn) c s = std deviation of ln(x1),...,ln(xn) c Output: 1. p-values for left-tail, right-tail and two-tail tests c 2. one-sided and two-sided confidence limits c Reference: Krishnamoorthy, K. and Mathew, T. (2002). Inferences on c the means of lognormal distributions using generalized p-values and c generalized confidence intervals. To appear in the Journal of Statistical c Planning and Inference c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real g_variable(100000), chi(1) real n mm = 100000 c Validate the input values 1 print *,'Enter the sample size = ' read *, n if(n .le. 0) then print *, 'sample size must be positive integer !! ' go to 1 end if print *,'Enter the mean of the logged data = ' read *, xbar 2 print *,'Enter the std deviation of the logged data = ' read *, s if(s .le. 0.0) then print *,'standard deviation must be a positive number !!' goto 2 end if 4 print *,'Enter 1 for Hypothesis Test; 2 for confidence interval = ' read *, check if(check .eq. 1) then print *,'Enter the value of the lognormal mean under H0 = ' read *, pel pel = alog(pel) elseif (check .eq. 2) then 3 print *,'Enter the confidence level = ' read *, conflev if(conflev .lt. 0.5 .or. conflev .ge. 1.0) then print *,'confidence level mut be in (0.5, 1.0)' goto 3 end if else goto 4 end if df = n-1.0 const1 = s*sqrt(df)/sqrt(n) const2 = s**2*df*0.5 c iseed = 123457 c pval_left = 0.0 pval_rigt = 0.0 pval_two = 0.0 c do j = 1, mm z = anorran() v = 2.0*gamran(df/2.0) g_variable(j) = xbar + z*const1/sqrt(v)+const2/v if(g_variable(j) .gt. pel) pval_left = pval_left + 1.0 if(g_variable(j) .lt. pel) pval_rigt = pval_rigt + 1.0 end do c if(check .eq. 2) then l1 = int((1.0-conflev)*mm) l2 = int(conflev*mm) lt1 = int(0.5*(1.0-conflev)*mm) lt2 = int(0.5*(1.0+conflev)*mm) call hpsort(mm, g_variable) print *, ' ' print *,'-------------------------------------------------' print *, ' ' write (*,5) conflev, exp(g_variable(l1)) 5 format(2x, f6.3, ' one-sided lower limit is ',f12.4) write (*,6) conflev, exp(g_variable(l2)) 6 format(2x, f6.3, ' one-sided upper limit is ',f12.4) write (*,7) conflev, exp(g_variable(lt1)), exp(g_variable(lt2)) 7 format(2x, f6.3, ' confidence interval is (',f12.4,',',f12.4,')') print *, ' ' print *,'-------------------------------------------------' end if pval_left = pval_left/mm pval_rigt = pval_rigt/mm pval_two = 2.0*min(pval_left, pval_rigt) if(check .eq. 1) then write(*,*) ' ' write(*,*) ' ' write(*,*) '----------------------------------------' write(*,*) ' ' write(*,8) pval_left 8 format(2x,'p-value for left-tail test is ',f6.4/) write(*,9) pval_rigt 9 format(2x,'p-value for right-tail test is ',f6.4/) write(*,10) pval_two 10 format(2x,'p-value for two-tail test is ',f6.4/) write(*,*) ' ' write(*,*) '----------------------------------------' end if 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