! Last change: DOS 29 Jul 2000 12:26 pm ! *** copyright 2000 *** ! *** filename splsmu.f95 *** John F. Monahan ** ! ********************** PROGRAM PSPLSMU ! TEST PROGRAM FOR SPLSMU -- COMPUTING SMOOTHING SPLINES ! ! TEST PROBLEM FROM R. L. EUBANK (1988) EXERCISE 21, PP.271-2 ! SPLINE SMOOTHING AND NONPARAMETRIC REGRESSION ! DATA ARE FROM EPPRIGHT, ET AL. (1972) ! IMPLICIT NONE REAL, DIMENSION(100) :: XK,Y,WE,H,Z,ZZ,M,MM,SCY,R REAL T,DET,A11 ! INTEGER, PARAMETER :: N = 72 INTEGER I,NM2,IDET ! 20 FORMAT(' LEVINSON-DURBIN ROUTE, DET=',F9.6,' IDET=',I9) 21 FORMAT(F5.1,F5.2) 22 FORMAT(6F13.6) 24 FORMAT(//'SMOOTHING SPLINES, X,Y, Z,Z M,M ') ! ! DATA ARE IN eppright.dat OPEN( UNIT=8, FILE='eppright.dat') ! PUT OUTPUT HERE OPEN( UNIT=9, FILE='splsmu.out') ! INITIALIZE TO ZERO R = 0. ! 72 OBSERVATIONS IN THIS PROBLEM NM2 = N - 2 ! EUBANK FINDS BEST GCV LAMBDA IS 1.38 ! T IS SMOOTHING PARAMETER HERE, T = N * LAMBDA T = 72*1.38 R(1) = (1.-24.*T)/(4.+36.*T) R(2) = (6.*T)/(4.+36.*T) ! READ IN THE DATA DO I = 1,N READ(8,21) XK(I),Y(I) WE(I) = 1.0 END DO ! LOOP ON I ! CALL SPLINE SMOOTHING PROGRAM, GET BACK ! FITTED VALUES IN Z ! SECOND DERIVATIVES IN M ! DIFFERENCES IN H ! CALL SPLSMU(XK,Y,WE,T,Z,M,H,N) ! ! NOW DO THIS ANOTHER WAY - MORE DIRECTLY ! EPPRIGHT PROBLEM HAS ALL EQUALLY SPACED X'S ! FIRST SET UP 6*C*Y - CAREFUL OF CANCELLATION DO I = 1,NM2 SCY(I) = 6.*( (Y(I+2)-Y(I+1)) - (Y(I+1)-Y(I)) ) END DO ! LOOP ON I ! EQUAL SPACING MAKES MATRIX WITH SIMPLE BANDS ! DIAGONAL AS IF VARIANCE R(0) ON BOTTOM ! OFF DIAGONALS BECOMES NUMERATORS FOR CORRELS R(1) = (1.-24.*T)/(4.+36.*T) R(2) = (6.*T)/(4.+36.*T) ! SOLVE EQUATION BY LEVINSON-DURBIN CALL LEVDRB(NM2,R,WE,SCY,A11,DET,IDET) WRITE(9,20) DET,IDET ! HAVE 2ND DERIVS IN SCY, NOW GET FITTED VALUES DO I = 1,N ZZ(I) = Y(I) END DO ! LOOP ON I ! MM(1) = 0. DO I = 1,NM2 SCY(I) = SCY(I) / (4.+36.*T) ! INDEXING DIFFERENT HERE, COPY INTO MM ! SO THAT OUTPUT IS EASIER MM(I+1) = SCY(I) ZZ(I) = ZZ(I) - SCY(I)*T ZZ(I+1) = ZZ(I+1) + 2.*SCY(I)*T ZZ(I+2) = ZZ(I+2) - SCY(I)*T END DO ! LOOP ON I ! WRITE OUT OUTPUT WRITE(9,24) ! X'S, Y, FITTED BOTH, SLOPES BOTH ! SPLSMU FIRST, LEVDRB SECOND DO I = 1,N WRITE(9,22) XK(I),Y(I),Z(I),ZZ(I),M(I),MM(I) END DO ! LOOP ON I STOP END PROGRAM PSPLSMU SUBROUTINE SPLSMU(XK,Y,WE,T,Z,M,H,N) ! CONSTRUCT SMOOTHING SPLINES TO MINIMIZE ! 2 2 ! MINIMIZE SUM WE(J)*(Y(J)-Z(J)) + T*INT | S"(X) | DX ! ! WHERE S(X) IS THE CUBIC SPLINE THAT INTERPOLATES (X(I),Z(I)) ! ! XK REAL VECTOR KNOTS ** IN INCREASING ORDER ** NO TIES ** ! Y REAL VECTOR ORDINATES ! WE REAL VECTOR WEIGHTS * IF K Y'S AT ONE X THEN WE = K ! T REAL SMOOTHING PARAMETER ! Z REAL VECTOR FITTED VALUES OF SMOOTHING SPLINE ! M REAL VECTOR M(J) IS SECOND DERIV OF SPLINE AT KNOT XK(J) ! H REAL VECTOR INTERVAL LENGTHS H(J) = XK(J)-XK(J-1) J=2,N ! N INTEGER NUMBER OF POINTS/KNOTS ** N MUST BE AT LEAST 5 ** ! ! INPUT XK,Y,WE,T,N ! OUTPUT Z,M,H VALUE OF SMOOTHNESS PENALITY IN H(1) ! ! J F MONAHAN (JAN,1988) DEPT OF STATISTICS, NCSU, RALEIGH, NC 27695 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( N ), INTENT(IN) :: XK,Y,WE REAL, DIMENSION( N ), INTENT(IN OUT) :: Z,M,H REAL, INTENT(IN) :: T ! REAL T6 REAL, DIMENSION(3*N) :: W INTEGER NM1,NM2,N2,J,JJ ! BEGIN WITH SOME CONSTANTS NM1 = N - 1 NM2 = N - 2 N2 = N + N ! DO J = 2,N ! INTERVAL LENGTHS H(J) = XK(J) - XK(J-1) ! HOLD RECIPROCALS IN Z FOR NOW --WILL CHANGE LATER Z(J) = 1./H(J) ! FIRST DIFFERENCES OF Y'S IN M M(J) = ( Y(J) - Y(J-1) ) / H(J) END DO ! LOOP ON J ! SET UP RIGHT HAND SIDE OF EQUATIONS TO BE SOLVED DO J = 2,NM1 M(J) = 6.*( M(J+1) - M(J) ) END DO ! LOOP ON J ! NOW CREATE BANDED MATRIX AND PUT IT IN W ! J,J ELEMENT IN W(J) ! J+1,J ELEMENT IN W(N+J) ! J+2,J ELEMENT IN W(N2+J) J = 2,...,NM1 T6 = 6.*T DO J = 2,NM1 W(J) = 2.*(H(J)+H(J+1)) + T6*( Z(J)*Z(J)/WE(J-1) & & + (Z(J)+Z(J+1))*(Z(J)+Z(J+1))/WE(J) + Z(J+1)*Z(J+1)/WE(J+1)) IF( J .LT. NM1 ) THEN W(N+J) = H(J+1) - T6*( (Z(J)+Z(J+1))*Z(J+1)/WE(J) & & + Z(J+1)*(Z(J+1)+Z(J+2))/WE(J+1) ) W(N2+J) = T6 * Z(J+1)*Z(J+2)/WE(J+1) ENDIF END DO ! LOOP ON J ! NOW DO CHOLESKY FACTORIZATION OF BANDED MATRIX ! AND DO THE FIRST HALF OF SOLVING ! START WITH THE FIRST CASES -- THEY'RE DIFFERENT W(2) = SQRT( W(2) ) M(2) = M(2)/W(2) W(N+2) = W(N+2) / W(2) W(3) = SQRT( W(3) - W(N+2)*W(N+2) ) M(3) = ( M(3) - W(N+2)*M(2) ) / W(3) ! DO J = 4,NM1 W(N2+J-2) = W(N2+J-2) / W(J-2) W(N+J-1) = ( W(N+J-1) - W(N2+J-2)*W(N+J-2) ) / W(J-1) W(J) = SQRT( W(J) - W(N+J-1)*W(N+J-1) - W(N2+J-2)*W(N2+J-2) ) M(J) = ( M(J) - W(N2+J-2)*M(J-2) - W(N+J-1)*M(J-1) ) / W(J) END DO ! LOOP ON J ! CHOLESKY FACTORING IS DONE -- NOW SOLVE BACKWARDS M(NM1) = M(NM1)/W(NM1) M(NM2) = ( M(NM2) - W(N+NM2)*M(NM1) ) / W(NM2) ! DO JJ = 4,NM1 J = N - JJ + 1 M(J) = ( M(J) - M(J+1)*W(N+J) - M(J+2)*W(N2+J) ) / W(J) END DO ! LOOP ON JJ ! NOW FINISH M -- WE HAVE A NATURAL SPLINE M(1) = 0. M(N) = 0. ! REMAINING TASK IS TO GET THE FITTED VALUES Z(J) ! START WITH THE FIRST TWO Z(1) = Y(1) - (T/WE(1))*Z(2)*M(2) Z(2) = Y(2) + ( (Z(2)+Z(3))*M(2) - Z(3)*M(3) ) * (T/WE(2)) ! OLD Z(J) = 1/H(J) REPLACED BY FITTED VALUES DO J = 3,NM2 Z(J) = Y(J) - (T/WE(J)) * ( Z(J)*M(J-1) - (Z(J)+Z(J+1))*M(J) & & + Z(J+1)*M(J+1) ) END DO ! LOOP ON J ! NOW LAST TWO Z(NM1) = Y(NM1) - (T/WE(NM1))*( Z(NM1)*M(NM2) & & - (Z(NM1)+Z(N))*M(NM1) ) Z(N) = Y(N) - (T/WE(N))*( Z(N)*M(NM1) ) ! FINISH BY COMPUTING SMOOTHNESS PENALTY H(1) = 0. DO J = 2,N H(1) = H(1) + ( M(J-1)*M(J-1) + M(J-1)*M(J) + M(J)*M(J) )*H(J) END DO ! LOOP ON J H(1) = H(1) / 3. ! FINISHED RETURN END SUBROUTINE SPLSMU subroutine levdrb(n,r,b,y,a11,det,idet) ! solve Toeplitz systems A*b = -r and A*x = y using Levinson-Durbin ! algorithm -- first one is Yule-Walker equations ! A is correlation matrix -- A(i,j) = r( |i-j| ) -- Toeplitz ! ! n in integer size of problem/matrix A ! r in real vector correlations r(1),...,r(n) r(0)=1 ! b out real vector solution to Yule-Walker equations A*b = -r ! y in/out real vector right hand side (in) / solution (out) A*x=y ! a11 out real quadratic form in inverse of A with vector of ones ! det out real determinant in det*(2**idet) form ! idet out integer ! ! J F Monahan (June,1992) Dept of Statistics, NCSU, Raleigh, 27695-8203 ! (Sept 1999, April 2000) rewritten for Fortran 95 ! implicit none INTEGER, INTENT(IN) :: n REAL, DIMENSION( n ), INTENT(IN) :: r REAL, DIMENSION( n ), INTENT(IN OUT) :: b,y REAL, INTENT(OUT) :: a11,det INTEGER, INTENT(OUT) :: idet REAL d,e,f,bi,bic integer k,kh,km1,i,ic ! det = 1. idet = 0 ! begin with k=1 b(1) = - r(1) a11 = 1. if( n .eq. 1 ) return d = 1. - b(1)*b(1) ! loop on k do k = 2,n km1 = k - 1 e = r(k) f = y(k) bi = 1. do i = 1,km1 ic = k-i e = e + r(i)*b(ic) f = f - y(i)*r(ic) bi = bi + b(i) end do ! loop on i ! compute the new elements at the bottom det = d * det call adjust(det,idet) b(k)= - e/d y(k) = f/d a11 = a11 + bi*bi/d d = d * (1. - b(k)*b(k)) ! update remaining elements - but treat b special do i = 1,km1 ic = k-i y(i) = y(i) + y(k)*b(ic) end do ! loop on i ! careful with b so do not overwrite incorrectly kh = k / 2 do i = 1,kh ic = k - i bi = b(i) + b(k)*b(ic) bic = b(ic) + b(k)*b(i) b(i) = bi b(ic) = bic END do ! loop on i ! end of loop on k END DO ! loop on k return end subroutine levdrb 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 ! *** end of filename splsmu.f95 *********************