! Last change: DOS 26 Jul 2000 8:40 pm ! *** copyright 2000 *** ! *** filename flyreg.f95 *** John F. Monahan ** ! ********************** program pflyreg ! flyreg -- regression algorithm using Givens transformations ! updates "on the fly" implicit none INTEGER, parameter :: n = 22 INTEGER, parameter :: p = 4 REAL, DIMENSION(p) :: x REAL, DIMENSION(p+2) :: b REAL, DIMENSION( (p*(p+1))/2 ) :: r real s,y ! integer i,j,k ! 21 format(14x,f5.3,f5.0,f4.0,f6.2) 22 format(i4,f12.6,2x,5f12.6) 25 format(//' regression coefficients '/4x,6f12.6) 26 format(/' SSE =',f14.4,4x,'Nobs ',F5.0,' Variance est',f12.6) 27 format(/'Covariance matrix of coefficients'/) 28 format(i4,2x,6f12.6) ! data b/6*0./, r/10*0./ ! ! use problem in Brownlee, p.463 on factors related to heart ! disease in various countries open( unit=5, file='brown463.dat' ) OPEN( UNIT=6, FILE='flyreg.out' ) ! ! read in each new observation and update matrices ! do i = 1,n read(5,21) x(2),x(3),x(4),y x(1) = 1. ! write it back out to verify write(6,22) i,y,(x(j),j=1,p) call flyreg(y,x,p,r,b) ! write it again -- y is now a z, x's made to zero write(6,22) i,y,(x(j),j=1,p) END DO ! loop on i ! ! now solve system of equations R b = z(1) for b call chlzih(r,p,b) write(6,25) (b(j),j=1,p) ! ! compute estimate of variance -- error mean square s = b(p+1)/(b(p+2)-p) write(6,26) b(5),b(6),s ! ! r is Cholesky factor is symmetric storage mode ! compute the inverse of (XTX) from it and then ! multiply to get covariance matrix of regression coefficients call chlzoi(r,p) ! write(6,27) k = 0 do i = 1,p do j = 1,i k = k + 1 r(k) = r(k)*s END DO ! loop on j write(6,28) i,(r(k-i+j),j=1,i) END DO ! loop on i stop end program pflyreg subroutine flyreg(y,x,p,r,z1) ! subroutine to compute regression using Givens' transformations ! on the fly -- takes one new row of x and new y at a time ! ! y new dependent observation ! x new vector of explanatory variables for this obs ! p dimension of x (integer) ! r updated cholesky factor of XX matrix ! - in symmetric storage ! z1 transformed previous y's for z1(1)...z1(p) ! z1(p+1) updated error sum of squares ! z1(p+2) number of observations ! ! *** both r and z1 must be initialized to zero to start *** ! dimension of vectors: ! x has length p = number of explanatory vectors ! r has length p(p+1)/2 ! z1 has length p+2 ! ! *** How to use: ! 0) initialize r and z1 to zero ! 1) for each new observation, input y, x and p ; and old r, z1 ! 2) when done solve R b = z1 for b using chlzih ! 3) sse stored in z1(p+1), observations counted in z(p+2) ! 4) invert (XTX) from Cholesky factor r using chlzoi ! ! J F Monahan (Sept 93) Dept of Statistics, NCSU, Raleigh, NC USA ! revised Sept 1999, April 2000 for Fortran 95 ! implicit none integer, intent(in) :: p real, intent(in out) :: y real, dimension(p), intent(in out) :: x real, dimension( (p*(p+1))/2 ), intent(in out) :: r real, dimension(p+2) :: z1 real a,b,s integer i,j,jh,ji ! do j = 1,p ! are we still starting out? jh = j a = r( (j*(j+1))/2 ) ! find Givens transformation b = x(j) s = a*a + b*b if( s .gt. 0. ) then s = sqrt( s ) a = a/s b = b/s ! apply transform to y s = a*z1(j) + b*y y = a*y - b*z1(j) z1(j) = s ! apply transform to x do i = j,p ji = (i*(i-1))/2 + j s = a*r( ji ) + b*x(i) x(i) = a*x(i) - b*r( ji ) r(ji) = s END DO ! loop on i ! END if END DO ! loop on j ! done with most updating, finish and go z1(p+1) = z1(p+1) + y*y z1(p+2) = z1(p+2) + 1. return end subroutine flyreg 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 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 SUBROUTINE CHLZOI(A,N) ! ! OVERWRITES L WITH INVERSE(L')*INVERSE(L) WHERE L IS A LOWER ! TRIANGULAR (N BY N) MATRIX STORED IN A ! ! TO COMPUTE THE INVERSE OF A POSITIVE DEFINITE MATRIX A ! FIRST: FACTOR A = L * L' WHERE L IS LOWER TRIANGULAR ! USING CHLSKY(A,N,D,IDET) ! SECOND: COMPUTE THE INVERSE(A) = INVERSE(L') * INVERSE(L) ! USING CHLSOI(A,N) ! ! J F MONAHAN (JULY,1984) DEPT OF STATISTICS, N C S U, RALEIGH, N C 2769 ! 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 S INTEGER I,J,K,IM1,KP1 ! FIRST INVERT L IN PLACE DO K = 1,N ! SOLVE L * X = E(K) A( (K*(K+1))/2 ) = 1./A( (K*(K+1))/2 ) ! BRANCH AROUND IF K EQUALS N IF( K .LT. N ) THEN KP1 = K + 1 DO I = KP1,N S = 0. IM1 = I - 1 DO J = K,IM1 S = S - A( (I*(I-1))/2 + J )*A( (J*(J-1))/2 + K ) ENDDO ! LOOP ON J A( (I*(I-1))/2 + K ) = S/A( (I*(I+1))/2 ) ENDDO ! LOOP ON I ENDIF ENDDO ! LOOP ON K ! DONE WITH INVERSION, NOW COMPUTE INNER PRODUCT DO J = 1,N DO I = J,N ! GET (I,J) ELEMENT OF MATRIX IS THE INNER PRODUCT ! OF COLUMN I AND COLUMN J S = 0. DO K = I,N S = S + A( (K*(K-1))/2 + I )*A( (K*(K-1))/2 + J ) ENDDO ! LOOP ON K ! STORE THE INNER PRODUCT A( (I*(I-1))/2 + J ) = S ENDDO ! LOOP ON I ENDDO ! LOOP ON J RETURN END SUBROUTINE CHLZOI ! *** end of filename flyreg.f95 *********************