! Last change: DOS 30 Jul 2000 1:50 pm ! *** copyright 2000 *** ! *** filename augrain1.f95 *** John F. Monahan ** ! ********************** ! ! SOME COMPILERS REQUIRE MODULE TO PREVIOUSLY EXIST OR BE ! AT THE BEGINNING MODULE passinfo real(KIND=8) sx,slx integer, parameter :: n = 45 ! number of years END MODULE passinfo ! program augrain1 ! augrain1 -- mle via concentrated likelihood ! for fitting two parameter gamma ! data are August rainfall data for Raleigh, NC ! declare functions passed as arguments ! pass data to functions using module use passinfo implicit none real(KIND=8), external :: dlc,dll real(KIND=8), dimension(2) :: v,typx,del real(KIND=8), dimension(3) :: h real(KIND=8) x,sxx,a0,b0,al,ar,f integer i ! formats 20 format(' Fit two param gamma to rainfall data'/ & & ' get mle via concentrated likelihood') 21 format(2x,f5.2) 22 format(' sum, variance, sum of logs',3f12.6) 23 format(' initial estimates from method of moments, '/ & & ' alpha=',f12.6,' beta=',f12.6,' initial ll=',f12.6) 24 format(/' final estimates from golden section, '/ & & ' alpha=',f12.6,' beta=',f12.6,' final ll=',f12.6) 25 format(' gradient vector ',2f12.6/ & & ' hessian matrix (aa, ab, bb)',3f12.6) 26 format(/' covariance matrix of estimates '/' alpha ',2e16.4/ & & ' beta ',2e16.4) 27 format(' std err(alpha)= ',e16.6,9x,'std err(beta)= ',e16.6) ! open( unit=6, file='augrain1.out') ! read August rainfall data -- 45 years open( unit=8, file='augrain.dat') ! write(6,20) ! initialize sums sx = 0.d0 sxx = 0.d0 slx = 0.d0 ! read data and compute sums do i = 1,n read(8,21) x sx = sx + x if( i .gt. 1 ) & & sxx = sxx + (real(i,8)*x-sx)**2 / real(i*(i-1),8) !KIND slx = slx + log(x) end do ! loop on i ! get sample variance sxx = sxx / real(n-1,8) !KIND write(6,22) sx,sxx,slx ! get starting values from method of moments b0 = real(n,8)*sxx / sx !KIND a0 = sx*sx / (real(n*n,8)*sxx) !KIND ! how good is prelim est v(1) = a0 v(2) = b0 f = dll(v) write(6,23) a0,b0,f ! max concentrated likelihood via golden al = a0/5.d0 ar = a0*5.d0 call golden(dlc,25,al,ar) v(1) = al v(2) = sx / (n*al) f = dll(v) write(6,24) v(1),v(2),f ! now get gradient and hessian typx(1) = v(1)*1.d-4 typx(2) = v(2)*1.d-4 call del12f(dll,v,2,typx,del,h) write(6,25) del,h ! get inverse x = h(1)*h(3) - h(2)*h(2) h(1) = - h(1)/x h(3) = - h(3)/x h(2) = h(2)/x write(6,26) h(3),h(2),h(2),h(1) ! standard errors al = sqrt(h(1)) ar = sqrt(h(3)) write(6,27) ar,al stop end program augrain1 real(KIND=8) function dlc(a) ! compute the concentrated log-likelihood function for the ! two parameter gamma family use passinfo implicit none real(KIND=8) a,dll,v(2) ! evaluate regular log-likelihood v(1) = a v(2) = sx / (n*a) dlc = dll(v) return end real(KIND=8) function dll(v) ! likelihood function for two parameter gamma use passinfo implicit none real(KIND=8), DIMENSION(2), INTENT(IN) :: v real(KIND=8) a,b,dlgama a = v(1) b = v(2) dll = - n*a*log(b) - n*dlgama(a) + (a-1.d0)*slx - sx/b return end REAL(KIND=8) FUNCTION DLGAMA(X) ! COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS ! X .LE. 0. RETURNS -1.0 IMPLICIT NONE REAL(KIND=8), INTENT(IN) :: X REAL(KIND=8), DIMENSION(7) :: P ! COEFS FOR RATIONAL FUNCTION REAL(KIND=8), DIMENSION(6) :: Q REAL(KIND=8), DIMENSION(3) :: PT ! COEFS FOR STIRLING CORRECTION REAL(KIND=8), DIMENSION(2) :: QT REAL(KIND=8) XN,A,V REAL(KIND=8), PARAMETER :: LSQR2PI = .918938533204672742D0 DATA P/ .378601050348257245d4, .207745979389418732d4, & & .893581804523749814d3, .222112396168011795d3, & & .489543462279099381d2, .61260674503360843d1, & & .778079585613300576d0 / DATA Q/ .378601050348257187d4, .476793860503687915d3, & & -.867230987531102994d3, .835500586679197696d2, & & .507884753288954097d2, -.134004147857813483d2 / DATA PT/ .288119283935546015d0, .498030766924499634d0, & & .691561607375687d-1 / DATA QT/ .345743140722674507d1, .609161691641660296d1 / ! USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS' ! ED. BY JF HART, SIAM SERIES, WILEY(1968) ! 27 MAR 74 JFM W/RWS,RM ! ** MODIFICATIONS JUNE 1988, 1993 ** ! RECODED OCTOBER 1999 FOR FORTRAN 95 IF( X .LE. 0.D0 ) THEN DLGAMA = -1.D0 ! NOT DESIGNED FOR NEGATIVE ARGUMENTS RETURN ENDIF ! FOR LARGE X (>8) USE STIRLING ROUTE IF( X .GT. 8.D0 ) THEN ! LARGE X -- USE STIRLING WITH CORRECTION XN = X * X ! CORRECTION FOR STIRLING V = (PT(1)+(PT(2)+PT(3)/XN)/XN) / (QT(1)+(QT(2)+1.D0/XN)/XN) DLGAMA = (X-0.5D0)*LOG(X) - X + V/X + LSQR2PI ELSE ! FOR SMALL X, USE REPRODUCTION FORMULA ! TO GET NUMBER BETWEEN 2 AND 3 V = 1.D0 XN = X DO WHILE( XN .LT. 2.D0 ) V = V/XN XN = XN + 1.D0 END DO ! WHILE( XN .LT. 2.D0 ) DO WHILE( XN .GT. 3.D0 ) XN = XN - 1.D0 V = V * XN END DO ! WHILE( XN .GT. 3.D0 ) A = XN-2.D0 ! RATIONAL FUNCTION APPROXIMATION V =(P(1)+A*(P(2)+A*(P(3)+A*(P(4)+A*(P(5)+A*(P(6)+A*P(7))))))) & & * V / (Q(1)+A*(Q(2)+A*(Q(3)+A*(Q(4)+A*(Q(5)+A*(Q(6)+A)))))) DLGAMA = LOG(V) ENDIF RETURN END FUNCTION DLGAMA SUBROUTINE DEL12F(F,X,N,DX,D1F,D2F) ! COMPUTES NUMERICAL FIRST AND SECOND DERIVATIVES FOR FUNCTION F ! GRADIENT VECTOR AND HESSIAN MATRIX ! ! F FUNCTION OF INTEREST ! X REAL V POINT AT WHICH GRADIENT AND HESSIAN ARE TO BE COMPUTED ! N INTEGER DIMENSION OF X ! DX REAL V CHANGES IN EACH COMPONENT OF X ! D1F REAL V ON OUTPUT, THE GRADIENT VECTOR ! D2F REAL V ON OUTPUT, THE HESSIAN MATRIX IN SYMMETRIC STORAGE ! IMPLICIT NONE REAL(KIND=8) F INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( N ), INTENT(IN OUT) :: X REAL(KIND=8), DIMENSION( N ), INTENT(IN) :: DX REAL(KIND=8), DIMENSION( N ), INTENT(OUT) :: D1F REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(OUT) :: D2F REAL(KIND=8) V0,VP,VN,XH,XG INTEGER I,J,K,KI,KJ ! CENTER POINT V0 = F(X) ! GET CENTERED FIRST DIFFERENCES DO I = 1,N XH = X(I) X(I) = XH + DX(I) D1F(I) = F(X) ! STORE F(X+H*EI) HERE X(I) = XH - DX(I) KI = (I*(I+1))/2 D2F(KI) = F(X) ! STORE F(X-H*EI) HERE X(I) = XH END DO ! LOOP ON I ! ! HAVE AXIS POINTS STORED, NOW DO THE HESSIAN DO I = 1,N XG = X(I) KI = (I*(I+1))/2 ! COMMON FOR EACH I DO J = I,N ! SKIP THE SECOND PARTIALS FOR LAST IF( I .NE. J ) THEN XH = X(J) X(I) = XG + DX(I) X(J) = XH + DX(J) ! EVALUATED AT THIS CORNER ON THE IJ FACE VP = F(X) ! X(I) = XG - DX(I) X(J) = XH - DX(J) ! EVALUATE AT THE NEGATIVE CORNER ON THE IJ FACE VN = F(X) ! X(J) = XH KJ = (J*(J+1))/2 K = KJ - J + I D2F(K) = ( ( (VP-D1F(I)) - (D1F(J)-V0) ) + ( (VN-D2F(KI)) - & & (D2F(KJ)-V0) ) ) / (2.D0*DX(I)*DX(J)) ! ENDIF END DO ! LOOP ON J ! WHEN DONE WITH ALL I CASES, GET GRADIENT AND HESSIAN VP = D1F(I) VN = D2F(KI) D1F(I) = (VP-VN)/(2.D0*DX(I)) D2F(KI) = ( (VP-V0) - (V0-VN) ) / (DX(I)*DX(I)) X(I) = XG ! END DO ! LOOP ON I ! RETURN END SUBROUTINE DEL12F SUBROUTINE GOLDEN(F,MIT,XL,XR) ! GOLDEN SECTION SEARCH FOR THE MAX OF F ON THE INTERVAL (XL,XR) ! F MUST BE UNIMODAL; ON RETURN XL=MODE, XR=VALUE OF F THERE IMPLICIT NONE REAL(KIND=8) F REAL(KIND=8), INTENT(IN OUT) :: XL,XR INTEGER, INTENT(IN) :: MIT REAL(KIND=8), PARAMETER :: TAU = .618034 ! GOLDEN RATIO REAL(KIND=8), PARAMETER :: TAUSQ = .381966 REAL(KIND=8) XLM,XRM,FXLM,FXRM INTEGER I,IT ! KEEP TRACK OF BETTER INSIDE POINT ! I=1 LEFT, I=-1 RIGHT I=0 DO IT = 1,MIT ! INTERVAL OF UNCERTAINTY IS (XL,XR) XL LE XLM LE XRM LE XR XLM = XL*TAU + XR*TAUSQ XRM = XL*TAUSQ + XR*TAU ! ARE ANY OF THESE ALREADY EVALUATED? IF(I.NE.-1) FXLM=F(XLM) IF(I.NE.1) FXRM=F(XRM) ! MAIN COMPARISON IF( FXLM .GT. FXRM ) THEN ! F(XLM) .GT. F(XRM) XR = XRM FXRM = FXLM I = 1 ELSE ! F(XLM) .LT. F(XRM) XL = XLM FXLM = FXRM I = -1 ENDIF END DO ! LOOP ON IT ! ALL DONE, TAKE THE BEST POINT SO FAR IF( I .EQ. -1 ) THEN XL = XRM XR = FXRM ELSE XL = XLM XR = FXLM ENDIF RETURN END SUBROUTINE GOLDEN ! *** end of filename augrain1.f95 *********************