! Last change: DOS 3 Aug 2000 5:22 pm ! *** copyright 2000 *** ! *** filename gilks.f95 *** John F. Monahan ** ! ********************** program pgilks ! Gilks' adaptive acceptance-rejection method ! test by generating from standard normal distribution ! and doing the same goodness of fit tests as before implicit none real, external :: h integer, parameter :: n = 4096 ! sample size real, dimension(n) :: x real XS,X2S,FN,PHIXI,DKSP,DKSM,DKS,AD real ran,gilks,alnphi,cdfn integer i ! 21 FORMAT(' Sample Mean and Variance for ',I8,' Normal Deviates' & & /' Mean is',F16.8,4x,'Variance is',F16.8) 22 FORMAT(' Goodness of Fit Statistics: '/4x,'Anderson-Darling is',& & F12.6/4x,'Compare to .05, .01 critical values of 2.492, 3.857'//& & 4x,'Kolmogorov-Smirnov D+',F9.6,4x,' D-',F9.6/ & & 4x,'SQRT(N)*MAX(D+,D-) is',F9.4/4x,'Compare to .05, .01', & & ' critical values of 1.36, 1.63') 25 format(i4,f8.4) ! open output file open( unit=6, file='gilks.out' ) ! let's check out 2**14 = 16384 FN = REAL(N) ! initialize uniform generator XS = RAN(5151917) ! initialize moment sums X2S = 0. XS = 0. DO I = 1,N X(I) = gilks(h, -2., -1., 1. ) ! use these three x's to start !** write(6,25) i,x(i) XS = XS + X(I) IF( I .GT. 1 ) X2S = X2S + (XS - I*X(I))*(XS - I*X(I))/(I*(I-1)) END DO ! LOOP ON I ! get sample mean and variance XS = XS/FN X2S = X2S/(FN-1.) WRITE(6,21) N,XS,X2S ! sort the data -- get order statistics CALL HSORT(X,N) ! now compute K-S and A-D statistics DKSP = 0. DKSM = 0. AD = 0. DO I = 1,N AD = AD - (2*I-1)*( ALNPHI(X(I)) + ALNPHI(-X(N+1-I)) )/FN - 1. PHIXI = CDFN( X(I) ) DKSP = MAX( DKSP, REAL(I)/FN - PHIXI ) DKSM = MAX( DKSM, PHIXI - REAL(I-1)/FN ) END DO ! LOOP ON I ! usual standardization of K-S DKS = SQRT(FN)*MAX( DKSP, DKSM ) WRITE(6,22) AD,DKSP,DKSM,DKS STOP END PROGRAM PGILKS real function h(x) ! unscaled log of density function ! a concave function real x h = - x*x / 2. ! normal distribution ! return end function h real function gilks(flog,x1,x2,x3) ! Gilks' secant-style adaptive acceptance/rejcection algorithm ! for generating from log-concave density ! ! flog real function to compute (unscaled) log of density ! x1 real one of three starting values ! x2 real second ! x3 real third ! ! J F Monahan (July, 1999) ! Recoded March 2000 for Fortran 95 implicit none real, intent(in) :: x1,x2,x3 real flog integer, parameter :: nmax = 100 ! max number of trials real, dimension(nmax) :: x,h,psum,c,d integer i,ii,n,nm1,nint real hmax,hmaxsum,xm,p,z,u real ran integer ifind ! initial setup hmaxsum = 0. x(1) = x1 h(1) = flog(x1) x(2) = x2 h(2) = flog(x2) x(3) = x3 h(3) = flog(x3) n = 3 call getcd ! check to see if left-most slope is postive do while( d(1) .le. 0. ) ! add a point on left n = n + 1 if( n .gt. nmax ) then write(*) 'failure in subroutine gilks' stop end if ! ( n .gt. nmax ) x(n) = ( 3.*x(1) - x(2) ) / 4. h(n) = flog( x(n) ) - hmaxsum call getcd end do ! while( d(1) .le. 0. ) ! left-most slope is positive -- check right do while( d(nm1) .ge. 0. ) ! add a point on right n = n + 1 if( n .ge. nmax ) then write(*) 'failure in subroutine gilks' stop end if ! ( n .gt. nmax ) x(n) = ( 5.*x(n-1) - x(n-2) ) / 4. h(n) = flog( x(n) ) - hmaxsum call getcd end do ! while( d(nm1) .ge. 0. ) ! got both ends tied down ! ! finally ready to generate ! do while( n .le. nmax ) !WW write(6,20) n !WW 20 format(i4) ! sum probabilities ! first two are leftmost intervals psum(1) = exp( c(1) + d(1)*x(1) ) / d(1) z = x(2) - x(1) ! beware zero slopes if( d(2) .ne. 0. ) z = ( 1. - exp( -d(2)*z ) ) / d(2) psum(2) = psum(1) + exp( c(2) + d(2)*x(2) ) * z ! next two are rightmost intervals z = x(n) - x(n-1) if( d(n-2) .ne. 0. ) z = ( 1. - exp( -d(n-2)*z ) ) / d(n-2) psum(3) = psum(2) + exp( c(n-2) + d(n-2)*x(n) ) * z psum(4) = psum(3) - exp( c(nm1) + d(nm1)*x(n) ) / d(nm1) if( n .gt. 3 ) then ! ! now do the others in pairs do ii = 3,nm1 ! find point in interval (x(ii-1),x(ii)) xm = ( c(ii) - c(ii-2) ) / ( d(ii-2) - d(ii) ) i = 2*ii - 1 ! left part of interval z = xm - x(ii-1) if( d(ii-2) .ne. 0. ) z = ( 1. - exp( -d(ii-2)*z ) ) / d(ii-2) psum(i) = psum(i-1) + exp( c(ii-2) + d(ii-2)*xm ) * z i = i + 1 ! right part of interval z = x(ii) - xm if( d(ii) .ne. 0. ) z = ( 1. - exp( -d(ii)*z ) ) / d(ii) psum(i) = psum(i-1) + exp( c(ii) + d(ii)*x(ii) ) * z end do ! loop on ii ! end if ! ( n .gt. 3 ) ! nint = number of intervals nint = 2*(n-3) + 4 !WW write(6,22) (psum(i),i=1,nint) !WW 22 format(10f8.2) ! acceptance/rejection p = ran(1) p = p*psum(nint) i = ifind(p,psum,nint) + 1 ! now generate in interval i ! first special cases select case(i) case(1) gilks = x(1) + log( ran(2) ) /d(1) z = log( ran(3) ) + c(1) + d(1)*gilks case(2) u = ran(2) xm = u*( x(2) - x(1) ) if( d(2) .ne. 0. ) xm = & & log( 1.+ u*(exp( d(2)*(x(2)-x(1)) )-1.) )/d(2) gilks = x(1) + xm z = log( ran(3) ) + c(2) + d(2)*gilks case(3) u = ran(2) xm = u*( x(n) - x(nm1) ) if( d(n-2) .ne. 0. ) xm = & & log(1.+u*(exp(d(n-2)*(x(n)-x(nm1)))-1.))/d(n-2) gilks = x(nm1) + xm z = log( ran(3) ) + c(n-2) + d(n-2)*gilks case(4) gilks = x(n) + log( ran(2) ) / d(nm1) z = log( ran(3) ) + c(nm1) + d(nm1)*gilks case default ! which interval? ii = ( i + 1) / 2 ! find point in interval (x(ii-1),x(ii)) xm = ( c(ii) - c(ii-2) ) / ( d(ii-2) - d(ii) ) if( 2*ii .ne. i ) then ! left part of interval ( x(ii-1), x(ii) ) u = ran(2) if( d(ii-2) .ne. 0. ) then gilks = x(ii-1) + & & log( 1. + u*(exp( d(ii-2)*(xm-x(ii-1)) )-1.) ) / d(ii-2) else gilks = xm + u*(xm-x(ii-1)) end if ! ( d(ii-2) .ne. 0. ) z = log( ran(3) ) + c(ii-2) + d(ii-2)*gilks ! quick accept check if( z .le. c(ii-1) + d(ii-1)*gilks ) return else ! right part of interval u = ran(2) if( d(ii) .ne. 0. ) then gilks = xm + & & log( 1. + u*(exp( d(ii)*(x(ii)-xm) )-1.) ) / d(ii) else gilks = xm + u*(x(ii)-xm) end if ! ( d(ii) .ne. 0. ) z = log( ran(3) ) + c(ii) + d(ii)*gilks ! quick accept check if( z .le. c(ii-1) + d(ii-1)*gilks ) return end if ! ( 2*ii .gt. i ) end select ! main acceptance rejection check h(n+1) = flog(gilks) - hmaxsum if( z .le. h(n+1) ) return n = n + 1 x(n) = gilks call getcd end do ! while ( n .le. nmax ) ! only get here if nmax failures in accept/reject write(*) 'failure in subroutine gilks' stop contains subroutine getcd ! computes slopes and intercepts (d's and c's) of secant lines ! through points (h(i),x(i)), i=1,...,n call hpsort(x,h,n) ! find and subtract max h -- rescale f hmax = h(1) do i = 2,n hmax = max( h(i), hmax ) end do ! loop on i do i = 1,n h(i) = h(i) - hmax end do ! loop on i hmaxsum = hmaxsum + hmax ! ! compute slopes and intercepts of secants nm1 = n - 1 do i = 1,nm1 d(i) = (h(i+1) - h(i)) / ( x(i+1) - x(i) ) c(i) = h(i) - d(i)*x(i) end do ! loop on i !WW 21 format(10f8.3) !WW write(6,21) (x(i),i=1,n) !WW write(6,21) (h(i),i=1,n) !WW write(6,21) (c(i),i=1,nm1) !WW write(6,21) (d(i),i=1,nm1) end subroutine getcd end function gilks SUBROUTINE HPSORT(K,M,N) ! HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N ! ! ARGUMENTS ! K REAL VECTOR OF KEYS TO BE SORTED ! M INTEGER VECTOR TO BE MOVED IN PARALLEL TO K ! N NUMBER OF ITEMS TO BE SORTED ! ! TO SORT A PARALLEL VECTOR OF RECORDS, USE HPSORT ! ! J F MONAHAN (DEC, 1999) FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN OUT) :: K,M REAL KK,MM INTEGER I,L,NCUR ! ! DO NOTHING IF THERE'S NOTHING TO DO IF( N .LE. 1 ) RETURN ! INITIALIZE TO BUILDHEAP PART (LOOP ON L) L = N/2 + 1 NCUR = N DO I = L,1,-1 CALL HEAPIFY(I) END DO ! LOOP ON I DO I = 2,N ! SWITCH CURRENT LARGEST WITH BOTTOM KK = K(1) K(1) = K(NCUR) K(NCUR) = KK ! ALSO SWITCH PARALLEL VECTOR MM = M(1) M(1) = M(NCUR) M(NCUR) = MM ! REHEAP WITH ONE SHORTER NCUR = NCUR - 1 CALL HEAPIFY(1) END DO ! LOOP ON NCUR RETURN CONTAINS SUBROUTINE HEAPIFY(II) INTEGER, INTENT(IN) :: II INTEGER I,J I = II DO J = 2*I ! IS IT A LEAF OR ARE THERE SONS? IF( J > NCUR ) EXIT ! A LEAF IF( J < NCUR ) THEN ! ANOTHER SON OF I IF( K(J+1) > K(J) ) J = J+1 ! LARGER SON IS K END IF ! A LEAF IF( K(J) > K(I) ) THEN ! EXCHANGE KK = K(J) K(J) = K(I) K(I) = KK ! ALSO SWITCH PARALLEL VECTOR MM = M(J) M(J) = M(I) M(I) = MM I = J ELSE ! EXIT -- HEAP PROPERTY EXIT END IF END DO ! WHILE NOT A LEAF END SUBROUTINE HEAPIFY END SUBROUTINE HPSORT INTEGER FUNCTION IFIND(X,XK,N) ! FIND I SUCH THAT XK(I) LE X LT XK(I+1) ! IFIND=0 IF X LT XK(1) ! IFIND=N IF X GT XK(N) ! J F MONAHAN JAN 1982 DEPT OF STAT, N C S U, RALEIGH, NC 27650 ! RECODED FEBRUARY, 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION(N), INTENT(IN) :: XK INTEGER IL,IU ! SPECIAL CASES IFIND = 0 IF( X .LT. XK(1) ) RETURN IFIND = N IF( X .GE. XK(N) ) RETURN IL = 1 IU = N DO WHILE( IU-IL .GT. 1 ) IFIND = ( IU + IL ) / 2 IF( X .EQ. XK(IFIND) ) RETURN IF( X .LT. XK(IFIND) ) THEN IU = IFIND ELSE IL = IFIND END IF ! ( X .LT. XK(IFIND) ) END DO ! WHILE ( IU-IL .GT. 1 ) IFIND = IL RETURN END FUNCTION IFIND SUBROUTINE HSORT(K,N) ! HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N ! ! ARGUMENTS ! K REAL VECTOR OF KEYS TO BE SORTED ! N NUMBER OF ITEMS TO BE SORTED ! ! TO SORT A PARALLEL VECTOR OF RECORDS, USE HKSORT ! ! J F MONAHAN (DEC, 1999) FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN OUT) :: K REAL KK INTEGER I,L,NCUR ! ! DO NOTHING IF THERE'S NOTHING TO DO IF( N .LE. 1 ) RETURN ! INITIALIZE TO BUILDHEAP PART (LOOP ON L) L = N/2 + 1 NCUR = N DO I = L,1,-1 CALL HEAPIFY(I) END DO ! LOOP ON I DO I = 2,N ! SWITCH CURRENT LARGEST WITH BOTTOM KK = K(1) K(1) = K(NCUR) K(NCUR) = KK ! REHEAP WITH ONE SHORTER NCUR = NCUR - 1 CALL HEAPIFY(1) END DO ! LOOP ON NCUR RETURN CONTAINS SUBROUTINE HEAPIFY(II) INTEGER, INTENT(IN) :: II INTEGER I,J I = II DO J = 2*I ! IS IT A LEAF OR ARE THERE SONS? IF( J > NCUR ) EXIT ! A LEAF IF( J < NCUR ) THEN ! ANOTHER SON OF I IF( K(J+1) > K(J) ) J = J+1 ! LARGER SON IS K END IF ! A LEAF IF( K(J) > K(I) ) THEN ! EXCHANGE KK = K(J) K(J) = K(I) K(I) = KK I = J ELSE ! EXIT -- HEAP PROPERTY EXIT END IF END DO ! WHILE NOT A LEAF END SUBROUTINE HEAPIFY END SUBROUTINE HSORT REAL FUNCTION RAN(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RAN = REAL(IX)*RP RETURN END FUNCTION RAN REAL FUNCTION CDFN(X) ! CDF FOR NORMAL DIST 22 MAY 74 JFM ! ** RECODED AS SEPARATE SUBPROGRAMS ** 14 JULY 1988 ** ! ** RECODED FOR READABILITY ** 25 JUNE 1996 ** ! ** RECODED FOR FORTRAN 95 ** 30 JAN 2000 ** ! ! CDFN=CDF FOR NORMAL DIST CDCN=1-CDFN ! CDFN=0.5*(1.0+ERF()/SQRT(2))) ! ! RANGE FOR X METHOD SIZE NO. ACCURACY ! (0.0,0.14) CONTINUED FRACTION 3 6.7.9 7 ! (0.14,3.54) RATIONAL FUNCTION 3,4 5663 6,MOSTLY 7 ! (3.54,11.3) RATIONAL FUNCTION 4,5 5704 5,MOSTLY 6 ! (11.3,13.15) RATIONAL FUNCTION 1,2 5721 5 ! NO. REFERS TO EQN NO. OR METHOD FROM 'COMPUTER APPROXIMATIONS' ! ED BY J F HART WILEY 1968 ! REFERENCE FOR ACCURACY 'TABLES OF ERROR FN AND ITS DERIVATIVE' ! NATL BUREAU OF STANDARDS AMS 41 1958 ! ACCURACY IS IN SIG FIGS FOR UNIVAC 1108 ! IMPLICIT NONE REAL, INTENT(IN) :: X REAL, DIMENSION(4) :: P2 REAL, DIMENSION(5) :: Q2,P3 REAL, DIMENSION(6) :: Q3 REAL, DIMENSION(2) :: P4 REAL, DIMENSION(3) :: Q4 REAL XA,XS,V DATA P2/ .100046412E+2, .842655286E+1, .346025933E+1, & & .562353612/ DATA Q2/ .100046412E+2, .197155807E+2, .157022881E+2, & & .609074879E+1, 1.0 / DATA P3/ .183983860E+2, .224003952E+2, .130613857E+2, & & .402838962E+1, .564201377/ DATA Q3/ .183983860E+2, .431607573E+2, .433645871E+2, & & .236403038E+2, .714083720E+1, 1.0/ DATA P4/ .148459208, .564187743/ DATA Q4/ .511437285, .262770594, 1.0/ ! DEFINE SOME STATEMENT FUNCTIONS ! RESCALE SINCE ORIGINAL CODING WAS FOR ERROR FUNCTION XA = ABS(X/1.414214) ! ABS( X / SQRT(2) ) XS = X*X/2. IF( XA .LE. 0.1 ) THEN ! FIRST PIECE FOR SMALLEST X -- USE CONTINUED FRACTION V = 1.12837917*XA / ( EXP(XS) * (1.0-2.0*XS/(3.0+0.8*XS)) ) CDFN = ( 1. + SIGN( V, X) )/2. RETURN END IF ! IF( XA .LE. 0.1 ) IF( XA .LE. 2.5 ) THEN ! SECOND PIECE -- USE STRAIGHT RATIONAL FUNCTION (3,4) CDFN = (P2(1)+XA*(P2(2)+XA*(P2(3)+XA*P2(4)))) / & & (2.*EXP(XS)*(Q2(1)+XA*(Q2(2)+XA*(Q2(3)+XA*(Q2(4)+XA*Q2(5)))))) IF( X .GT. 0. ) CDFN = 1. - CDFN RETURN END IF ! IF( XA .LE. 2.5 ) IF( X .GT. 5. ) THEN ! TOO FAR IN RIGHT TAIL -- RETURN ONE CDFN = 1. RETURN END IF ! IF( X .GT. 5. ) IF( XA .LE. 8. ) THEN ! THIRD PIECE -- RATIONAL FUNCTION (4,5) FOR NEAR TAIL CDFN = (P3(1)+XA*(P3(2)+XA*(P3(3)+XA*(P3(4)+XA*P3(5))))) / (2.* & & (Q3(1)+XA*(Q3(2)+XA*(Q3(3)+XA*(Q3(4)+XA*(Q3(5)+XA*Q3(6))))))) CDFN = CDFN / EXP(XS) IF( X .GT. 0. ) CDFN = 1. - CDFN RETURN END IF ! IF( XA .LE. 8. ) IF( X .GE. -9.3 ) THEN ! FOURTH (TAIL) PIECE -- ANOTHER RATIONAL FUNCTION (1,2) FOR TAIL CDFN = (P4(1)+XA*P4(2))*EXP(-XS)/(2.*(Q4(1)+XA*(Q4(2)+XA*Q4(3)))) RETURN END IF ! IF( XA .LE. 9.3 ) ! TOO FAR IN LEFT TAIL -- RETURN ZERO CDFN = 0. RETURN END FUNCTION CDFN REAL FUNCTION ALNPHI(X) ! COMPUTES THE NATURAL LOGARITHM OF THE NORMAL DF ! J F MONAHAN JAN 1981 DEPT. OF STAT., N C S U, RALEIGH, N C 27650 ! RELATIVE ACCURACY 4E-6 FOR -4 14 IF( X .GT. 14.) RETURN AX = ABS(X) IF( AX .GT. 2. ) THEN ! DO THE TAILS FIRST XX = X*X XR32 = (((PR3/XX+PR2)/XX+PR1)/AX+AX)/((QR2/XX+QR1)/XX+1.) IF( X .GT. 0. ) THEN OMPHI = EXP(-X*X/2.)*RSR2PI/XR32 ! OMPHI=1-PHI NEXT USE PADE APPROX TO GET LN(1-OMPHI) ALNPHI = -OMPHI*(6-OMPHI)/(6-4*OMPHI) RETURN END IF ! IF( X .GT. 0. ) ALNPHI = LOG(RSR2PI/XR32) - X*(0.5*X) RETURN END IF ! IF( AX .GT. 2. ) IF( X .LE. 0.60 ) THEN ! NOW DO THE RANGE -2 < X < 0.6 ALNPHI = (PS0+X*(PS1+X*(PS2+X*(PS3+X*PS4))))/(1.0+X*(QS1+X*QS2)) RETURN END IF ! IF( X .LT. 0.60 ) ! LAST TAKE CARE OF 0.6 < X < 2 ALNPHI = (PT0+X*(PT1+X*(PT2+X*(PT3+X*PT4)))) & & *EXP(-0.5*X*X)/(1.0+X*(QT1+X*QT2)) RETURN END FUNCTION ALNPHI ! *** end of filename gilks.f95 *********************