! Last change: DOS 29 Jul 2000 12:24 pm ! *** copyright 2000 *** ! *** filename splrgm.f95 *** John F. Monahan ** ! ********************** PROGRAM PSPLRGM ! TEST PROGRAM FOR SPLRGM -- INTERPOLATORY CUBIC SPLINE REGRESSION ! A LA -- D. J. POIRIER, PIECEWISE REGRESSION USING CUBIC SPLINES ! (1973) JASA 68, 515-524. ! IMPLICIT NONE INTEGER, PARAMETER :: NOBS = 55 ! NUMBER OF OBSERVATIONS INTEGER, PARAMETER :: N = 4 ! NUMBER OF KNOTS REAL, DIMENSION(NOBS) :: X,Y REAL, DIMENSION(N) :: XK,M,H,Z REAL, DIMENSION( (N*(N+1))/2 ) :: COVZ REAL XI,YI,SSE REAL SPLEV INTEGER I,J,IX ! 20 FORMAT(1X,I4,3X,F7.3) 21 FORMAT(1X,F5.0,2X,F7.3) 22 FORMAT(1X,4F12.6) 25 FORMAT(' KNOTS, FITTED VALUES, SECOND DERIVS THERE') 27 FORMAT(/'COVARIANCE MATRIX OF FITTED VALUES') 28 FORMAT(/' SSE FIRST FROM SPLRGM, THEN USING FITTED VALUES') ! ! N=4 KNOTS OR JOIN POINTS -- REGIME CHANGES DATA XK/ 1., 7.5, 33.5, 61. / ! DATA ARE IN FILE 'indy500.dat' OPEN( UNIT=5, FILE='indy500.dat') ! PUT OUTPUT HERE OPEN( UNIT=6, FILE='splrgm.out') ! DO I = 1,NOBS READ(5,20) IX,Y(I) X(I) = IX - 1910 WRITE(6,21) X(I),Y(I) END DO ! LOOP ON I ! NOTICE THAT POIRIER USES THESE AS END CONDITIONS M(1) = 2. M(N) = 2. CALL SPLRGM(X,Y,NOBS,XK,Z,H,M,N,COVZ,SSE) ! WRITE(6,25) DO J = 1,N WRITE(6,22) XK(J),Z(J),M(J) END DO ! LOOP ON J WRITE(6,27) DO I = 1,N WRITE(6,22) ( COVZ( (I*(I-1))/2 + J ),J=1,I) END DO ! LOOP ON I WRITE(6,28) ! THIS SSE FROM SPLRGM WRITE(6,22) SSE SSE = 0. DO I = 1,NOBS XI = X(I) YI = SPLEV(XI,XK,Z,M,H,N) SSE = SSE + ( YI - Y(I) )*( YI - Y(I) ) END DO ! LOOP ON I ! THIS SSE DIRECTLY USING SPLEV'S FITTED VALUES WRITE(6,22) SSE STOP END PROGRAM PSPLRGM SUBROUTINE SPLRGM(X,Y,NOBS,XK,Z,H,M,N,COVZ,SSE) ! COMPUTES THE FITTED VALUES AND SLOPES FOR INTERPOLATORY CUBIC ! SPLINE REGRESSION MODEL ! ! INPUT ! X REAL VECTOR DESIGN POINTS ! Y REAL VECTOR ENDOGENOUS VARIABLES (RESPONSES) ! NOBS INTEGER NUMBER OF OBSERVATIONS *** AT LEAST N *** ! XK REAL VECTOR KNOTS (JOIN POINTS) ! N INTEGER NUMBER OF KNOTS *** AT LEAST 3 *** ! M(1),M(N) REALS PI(FIRST) AND PI(LAST) FOR END CONDITION ! ! OUTPUT ! Z REAL VECTOR FITTED VALUES AT KNOTS ! M REAL VECTOR SECOND DERIVATIVES OF FITTED SPLINE AT KNOTS ! H REAL VECTOR INTERVAL LENGTHS XK(J)-XK(J-1) J=2,...,N ! COVZ REAL VECTOR COVARIANCE MATRIX OF FITTED VALUES Z ! SSE REAL ERROR SUM OF SQUARES AFTER FITTING SPLINE ! ! EXTERNAL SUBPROGRAMS ! ROT734 SUBROUTINE TO COMPUTE GIVENS TRANSFORMATIONS ! CHLZOI SUBROUTINE TO COMPUTE INVERSE OF MATRIX FROM CHOLESKY FACTOR ! ! J F MONAHAN (FEBRUARY, 1988) DEPT OF STATISTICS, NCSU, RALEIGH NC ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! ! REFERENCE ! DALE J POIRIER, PIECEWISE REGRESSION USING CUBIC SPLINES, ! J AM STATISTICAL ASSN 68, PP. 515-524. ! IMPLICIT NONE INTEGER, INTENT(IN) :: NOBS,N REAL, DIMENSION( NOBS ), INTENT(IN) :: X,Y REAL, DIMENSION( N ), INTENT(IN OUT) :: XK,Z,M,H REAL, DIMENSION( (N*(N+1))/2 ), INTENT(OUT) :: COVZ REAL, INTENT(OUT) :: SSE REAL, DIMENSION(3*N) :: W REAL S,SC,SL,YC,XC,DXKC,ALP,BET,GAM,SIG,RNU ! INTEGER NM1,N2,I,J,JX,JJ,JP1,K,L INTEGER JFIND ! NM1 = N - 1 N2 = 2*N ! DO J = 2,N H(J) = XK(J) - XK(J-1) END DO ! LOOP ON J ! DO J = 2,NM1 S = H(J) + H(J+1) W( N2 + J ) = H(J+1)/S W( N + J ) = H(J)/S W(J) = 2. END DO ! LOOP ON J W(1) = 2. W(N) = 2. ! THIS IS PI(FIRST) RATIO OF FIRST TWO SECOND DERIVS W(N2+1) = -2.*M(1) ! THIS IS PI(LAST) RATIO OF LAST TWO SECOND DERIVS W(N2) = - 2.*M(N) ! ZEROS FOR M(1) AND M(N) MEAN NATURAL SPLINE ! ! IN W IS POIRIER'S LAMDA MATRIX -- NOW FACTOR DO J = 2,N W(N+J) = W(N+J) / W(J-1) W(J) = W(J) - W(N+J)*W(N2+J-1) END DO ! LOOP ON J ! NOW FACTORED INTO LU BOTH BIDIAG, L IS UNIT LOWER ! ! INITIALIZE Z AND COVZ FOR R OF REGRESSION ! Z = 0. COVZ = 0. SSE = 0. ! ! FOR OBS I SET UP ROW VECTOR OF DESIGN MATRIX W ! DO I = 1,NOBS ! INITIALIZE DO J = 1,N M(J) = 0. END DO ! LOOP ON J ! YC = Y(I) XC = X(I) ! JX = JFIND(XC,XK,N) ! DXKC = XK(JX) - XC M(JX-1) = (DXKC/(6.*H(JX))) * (DXKC*DXKC - H(JX)*H(JX) ) DXKC = XC - XK(JX-1) M(JX) = (DXKC/(6.*H(JX))) * (DXKC*DXKC - H(JX)*H(JX) ) ! HAVE A ROW OF P MULT IN BACK BY INVERSE OF LAMDA M(1) = M(1) / W(1) DO J = 2,N M(J) = ( M(J) - W(N2+J-1)*M(J-1) ) / W(J) END DO ! LOOP ON J DO JJ = 2,N J = N+1-JJ M(J) = M(J) - W(N+J+1)*M(J+1) END DO ! LOOP ON JJ ! MULT BY THETA NOW ON BACK ! THETA IS TRIDIAGONAL SL = 0. SC = 0. ! LOOP INDEX J IS ROW OF THETA DO J = 2,NM1 SL = SL + M(J) / ( H(J) * (H(J)+H(J+1)) ) SC = SC - M(J) / ( H(J) * H(J+1) ) ! STORE PREVIOUS ONE -- SL M(J-1) = 6.*SL ! SHIFT ROLES FOR NEXT J SL = SC SC = M(J) / ( H(J+1) * (H(J) + H(J+1)) ) END DO ! LOOP ON J ! FINISH THE UNCOMPLETED ONES M(NM1) = 6.*SL M(N) = 6.*SC ! NOW ADD Q ON TO THIS ROW M(JX-1) = M(JX-1) + ( XK(JX) - XC )/H(JX) M(JX) = M(JX) + ( XC - XK(JX-1) )/H(JX) ! NOW HAVE A ROW OF W IN M FOR OBSERVATION I ! ! DO GIVENS TRANSFORMS ! DO J = 1,N ALP = COVZ( (J*(J+1))/2 ) BET = M(J) CALL ROT734(ALP,BET,GAM,SIG,RNU) ! HAVE ROTATION FOR THIS COMPONENT J ! APPLY TO ALL OF ROW DO JJ = J,N K = (JJ*(JJ-1))/2 + J ! SAVE R IN TRANSPOSED FORM ALP = COVZ(K) BET = M(JJ) COVZ(K) = GAM*ALP + SIG*BET M(JJ) = GAM*BET - SIG*ALP END DO ! LOOP ON JJ ! APPLY ROTATION TO THE Y ALP = Z(J) BET = YC Z(J) = GAM*ALP + SIG*BET YC = GAM*BET - SIG*ALP END DO ! LOOP ON J SSE = SSE + YC*YC ! COMPUTATIONS COMPLETED FOR OBSERVATION I END DO ! LOOP ON I ! ! NOW SOLVE FOR FITTED VALUES Z -- HAVE UPPER TRIANGULAR R IN COVZ ! BUT SAVED IN TRANSPOSED FORM ! (I,J) IN (J*(J-1))/2 + I Z(N) = Z(N) / COVZ( (N*(N+1))/2 ) DO J = N-1,1,-1 JP1 = J + 1 DO L = JP1,N Z(J) = Z(J) - COVZ( (L*(L-1))/2 + J ) * Z(L) END DO ! LOOP ON L Z(J) = Z(J) / COVZ( (J*(J+1))/2 ) END DO ! LOOP ON JJ ! ! HAVE FITTED VALUES AT KNOTS IN Z -- FIRST GET SLOPES ! M(1) = 0. M(N) = 0. DO J = 2,NM1 M(J) = 6.* ( Z(J-1)/(H(J)*(H(J)+H(J+1))) & & - Z(J)/(H(J)*H(J+1)) + Z(J+1)/(H(J+1)*(H(J)+H(J+1))) ) END DO ! LOOP ON J ! ! HAVE LAMDA TIMES FITTED VALUES, NOW SOLVE FOR SLOPES ! DO J = 2,N M(J) = M(J) - W(N+J)*M(J-1) END DO ! LOOP ON J M(N) = M(N) / W(N) DO JJ = 2,N J = N + 1 - JJ M(J) = ( M(J) - W(N2+J)*M(J+1) ) / W(J) END DO ! LOOP ON JJ ! ! -1 ! COMPUTE (X'X) USING CHLZOI CALL CHLZOI(COVZ,N) ! ! NOW JUST PUT ON OUR ESTIMATE OF VARIANCE TO SCALE ! JJ = (N*(N+1))/2 S = SSE / REAL( NOBS - N ) DO J = 1,JJ COVZ(J) = COVZ(J) * S END DO ! LOOP ON J ! END OF SPLRGM COMPUTATIONS RETURN END SUBROUTINE SPLRGM REAL FUNCTION SPLEV(X,XK,Y,M,H,N) ! EVALUATES SPLINE AT X XK(JC-1) LT X LE XK(JC) ! CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS ! M(J) GIVES SECOND DERIV AT XKNOT(J) J=1,...,N ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION( N ), INTENT(IN) :: XK,Y,M,H REAL DA,DB INTEGER, SAVE :: J = 0 INTEGER JFIND ! INITIAL VALUE -- RETURN THIS IF OUT OF RANGE SPLEV = 0. ! IF X OUT OF RANGE, SPLEV=0 AND QUIT IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN ! IF LAST CALL IN SAME INTERVAL, NO NEED TO SEARCH IF( ( J .LE. 1 ) .OR. ( J .GT. N ) ) THEN J = JFIND(X,XK,N) ! FIND THE INTERVAL ELSE IF( (X .LE. XK(J-1)) .OR. (X .GT. XK(J)) ) THEN J = JFIND(X,XK,N) ! FIND THE INTERVAL ENDIF ENDIF ! COMPUTE DIFFERENCES FROM NEARBY KNOTS DB = X-XK(J-1) DA = XK(J)-X ! EVALUATE CUBIC SPLINE SPLEV = ( (M(J-1)*DA**3 + M(J)*DB**3)/6. & & + (Y(J-1)-M(J-1)*H(J)*H(J)/6.)*DA & & + (Y(J)-M(J)*H(J)*H(J)/6.)*DB ) / H(J) RETURN END FUNCTION SPLEV INTEGER FUNCTION JFIND(X,XK,N) ! FIND J SUCH THAT XK(J-1) LT X LE XK(J) ! JFIND=2 IF X EQ XK(1) JFIND=0 IF X IS OUT OF RANGE ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION( N ), INTENT(IN) :: XK INTEGER JU,JL,I JFIND = 0 ! IF X OUT OF RANGE, JFIND=0 AND QUIT IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN ! CHECK SPECIAL CASE IF( X .EQ. XK(1) ) THEN JFIND = 2 RETURN END IF ! INITIALIZE POINTERS JU = N JL = 1 DO WHILE( JU-JL .GT. 1 ) ! TEST THIS ONE I = ( JU + JL ) / 2 IF( X .LE. XK(I) ) THEN JU = I ELSE JL = I END IF END DO ! WHILE( JU-JL .GT. 1 ) JFIND = JU RETURN END FUNCTION JFIND 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 SUBROUTINE ROT734(ALPHA,BETA,GAMA,SIGMA,RNU) ! GIVEN ALPHA AND BETA, PRODUCES GIVENS TRANSFORMATION VALUES ! ! ( GAMA SIGMA ) ( ALPHA ) ( RNU ) ! ( ) ( ) = ( ) ! ( -SIGMA GAMA ) ( BETA ) ( ZERO ) ! ! ALGORITHM 7.3.4 FROM STEWART ! IMPLICIT NONE REAL, INTENT(IN) :: ALPHA,BETA REAL, INTENT(OUT) :: GAMA,SIGMA,RNU REAL ETA,DELTA ! ETA = MAX( ABS(ALPHA), ABS(BETA) ) ! IF BOTH ALPHA AND BETA ARE ZERO THEN ! USE IDENTITY TRANSFORMATION IF( ETA .EQ. 0. ) then ! SET UP IDENTITY TRANSFORMATION IF NO WORK DONE RNU = 0. GAMA = 1. SIGMA = 0. RETURN else DELTA = (ALPHA/ETA)*(ALPHA/ETA) + (BETA/ETA)*(BETA/ETA) DELTA = SQRT(DELTA) RNU = ETA * DELTA GAMA = ALPHA / RNU SIGMA = BETA / RNU END if RETURN END SUBROUTINE ROT734 ! *** end of filename splrgm.f95 *********************