! Last change: DOS 30 Jul 2000 2:01 pm ! *** copyright 2000 *** ! *** filename chex98.f95 *** John F. Monahan ** ! ********************** program chex98 ! Example 9.8 -- Logistic Regression using Finney's data on ! vaso-constriction of digits in response to air flow ! in rate and volume ! See D.J. Finney (1947) Biometrika, pp. 320-334 ! implicit none integer, parameter :: nobs = 39 integer, parameter :: np = 3 integer, dimension(nobs) :: y integer i,nn,nr,mi,k,krep,idet real, dimension(nobs) :: x2,x3 real, dimension(np) :: beta,dell real, dimension( (np*(np+1))/2 ) :: del2l real gam,d,pi,rll,det,flogit ! ! i = 1,...,39 ! ! mi number in each i category (mi=1 here) ! y(i) number of (positive) responses ! 20 format(2f6.3,i3) 21 format(' number of patients',i4,' positive responses',i4) 22 format(' new estimates',3f14.7) 23 format(' gradient -- dell',3f14.7) 24 format(' new step ',3f14.7) 25 format(/' iteration ',i3,' log-likelihood',f14.7) ! ! read in the data ! open( unit=5, file='finney.dat') open( unit=6, file='chex98.out') ! mi = 1 nn = 0 nr = 0 do i = 1,nobs read(5,20) x2(i),x3(i),y(i) x2(i) = log( x2(i) ) x3(i) = log( x3(i) ) nn = nn + mi nr = nr + y(i) end do ! loop on i (nobs) write(6,21) nn,nr ! set initial rate values beta(1) = log(real(nr)/real(nn-nr)) write(6,22) beta ! start looping do krep = 1,10 ! ! create dell and del2l ! ! initialize to zero rll = 0. dell = 0. del2l = 0. ! go through data do i = 1,nobs gam = beta(1) + x2(i)*beta(2) + x3(i)*beta(3) pi = flogit(gam) ! prob of response for this obs ! first log likelihood rll = rll + y(i)*gam - mi*log(1. + exp(gam) ) ! d = y(i) - mi*pi ! dell dell(1) = dell(1) + d dell(2) = dell(2) + d*x2(i) dell(3) = dell(3) + d*x3(i) ! del2l in symmetric form gam = mi*pi*flogit(-gam) del2l(1) = del2l(1) + gam del2l(2) = del2l(2) + x2(i)*gam del2l(3) = del2l(3) + x2(i)*x2(i)*gam del2l(4) = del2l(4) + x3(i)*gam del2l(5) = del2l(5) + x2(i)*x3(i)*gam del2l(6) = del2l(6) + x3(i)*x3(i)*gam ! updated dell and del2l for this obs end do ! loop on i (obs) ! write rll and dell write(6,25) krep,rll write(6,23) dell ! call chlzky(del2l,np,det,idet) ! solve for Newton step call chlzhi(del2l,np,dell) call chlzih(del2l,np,dell) ! write out next step write(6,24) dell ! move beta - notice sign adjustments do k = 1,np beta(k) = beta(k) + dell(k) end do ! loop on k write(6,22) beta ! end of loop end do ! loop on krep ! stop end program chex98 REAL FUNCTION FLOGIT(S) ! COMPUTE THE DISTRIBUTION FUNCTION OF LOGISTIC IMPLICIT NONE REAL, INTENT(IN) :: S REAL AS,E AS = ABS(S) E = 0.D0 IF(AS .LT. 125. ) E = EXP(-AS) IF( S .LT. 0. ) FLOGIT = E/(1. + E) IF( S .GE. 0. ) FLOGIT = 1./(1. + E) RETURN END FUNCTION FLOGIT SUBROUTINE ADJUST(D,I) ! NORMALIZES D WHILE KEEPING CONSTANT VALUE OF D * (2**I) INTEGER, PARAMETER :: IBIG = -2147483644 REAL, INTENT (IN OUT) :: D INTEGER, INTENT(IN OUT) :: I IF ( D .GT. 0.0 ) THEN DO WHILE ( D .LT. 1.0 ) D = D*16. I = I-4 ENDDO DO WHILE ( D .GT. 16.0 ) D = D/16. I = I+4 ENDDO ELSE I = IBIG ! IF D < 0 THEN I = - BIG ENDIF RETURN END SUBROUTINE ADJUST SUBROUTINE CHLZKY(A,N,DET,IDET) ! CHLSKY COMPUTES THE CHOLESKY (SQUARE-ROOT) FACTORIZATION ! A = L * ( L - TRANSPOSE ) WHERE L IS LOWER TRIANGULAR ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! A IS OVERWRITTEN BY L IN ITS LOWER TRIANGULAR PART ! SUBROUTINE ADJUST KEEPS DETERMINANT FROM EXPLODING USING ! 1 LE DET LE 16 AND DETERMINANT OF A IS DET*2**IDET ! ! J F MONAHAN (DEC,1983) DEPT OF STAT, N C S U, RALEIGH, N C 27650 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL, INTENT(OUT) :: DET INTEGER, INTENT(OUT) :: IDET REAL T,S INTEGER I,J,K,IM1,JM1 ! INITIALIZE FOR DETERMINANT DET = 1. IDET = 0 ! DO I-TH ROW DO I = 1,N T = A( (I*(I+1))/2 ) ! FIRST ROW IS TRIVIAL IF(I .GT. 1) THEN IM1 = I-1 DO J=1,IM1 ! WORK ON (I,J)-TH ELEMENT S = A( (I*(I-1))/2 + J ) IF(J .GT. 1) THEN JM1 = J-1 DO K = 1,JM1 S = S - A( (J*(J-1))/2 + K ) * A( (I*(I-1))/2 + K ) ENDDO ! LOOP ON K ENDIF A( (I*(I-1))/2 + J )= S / A( (J*(J+1))/2 ) ! WORK ON DIAGONAL ELEMENT T = T - A( (I*(I-1))/2 + J ) * A( (I*(I-1))/2 + J ) ENDDO ! LOOP ON J ! FINISHED WITH LOOP ON J ! UPDATE DETERMINANT WITH DIAGONAL ENDIF DET = DET * T ! KEEP DET FROM EXPLODING CALL ADJUST(DET,IDET) ! ABORT IF DET IS NEGATIVE OR ZERO IF(IDET.LT.-2000000000) RETURN ! FINISH A POSITIVE DIAGONAL A( (I*(I+1))/2 ) = SQRT(T) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZKY SUBROUTINE CHLZHI(A,N,Y) ! ! OVERWRITES Y WITH INVERSE( L ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE QUADRATIC FORM IN THE INVERSE OF A POSITIVE DEFINITE ! MATRIX A, THAT IS Y' INVERSE( A ) Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE THE SUM OF SQUARES OF THE ELEMENTS OF Y ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL, DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,J,IM1 ! ! SINCE THE MATRIX IS LOWER TRIANGULAR, SOLVE FROM THE TOP DOWN ! DO I = 1,N ! BRANCH AROUND ON THE FIRST ONE IF( I .GT. 1 ) THEN IM1 = I - 1 DO J = 1,IM1 Y(I) = Y(I) - A( (I*(I-1))/2 + J )*Y(J) ENDDO ENDIF Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZHI SUBROUTINE CHLZIH(A,N,Y) ! *** NOTICE TRANSPOSE *** ! OVERWRITES Y WITH INVERSE( L' ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE THE PRODUCT OF THE INVERSE OF A POSITIVE DEFINITE MATRIX ! AND A VECTOR Y, THAT IS INVERSE( A ) * Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE INVERSE(L') * ( INVERSE(L)*Y ) BY CHLSIH(A,N,Y) ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL, DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,II,IP1,J ! ! SINCE THE MATRIX IS UPPER TRIANGULAR, SOLVE FROM THE BOTTOM UP ! ! DO THE LAST ONE SEPARATELY Y(N) = Y(N) / A( (N*(N+1))/2 ) ! IF N = 1 THEN QUIT ELSE START LOOPING IF( N .EQ. 1 ) RETURN DO II = 2,N ! I BELOW NOW GOES FROM N-1 DOWN TO 1 I = N + 1 - II IP1 = I + 1 DO J = IP1,N Y(I) = Y(I) - A( (J*(J-1))/2 + I )*Y(J) ENDDO ! LOOP ON J Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON II,I RETURN END SUBROUTINE CHLZIH ! *** end of filename chex98.f95 *********************