home *** CD-ROM | disk | FTP | other *** search
/ High Voltage Shareware / high1.zip / high1 / DIR9 / ACMALG01.ZIP / ACM642.FOR < prev    next >
Text File  |  1991-02-18  |  37KB  |  913 lines
  1. C     ALGORITHM 642 COLLECTED ALGORITHMS FROM ACM.
  2. C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.12, NO. 2,
  3. C     JUNE, 1986, PP. 150-153.
  4. C   SUBROUTINE NAME     - CUBGCV
  5. C
  6. C--------------------------------------------------------------------------
  7. C
  8. C   COMPUTER            - VAX/DOUBLE
  9. C
  10. C   AUTHOR              - M.F.HUTCHINSON
  11. C                         CSIRO DIVISION OF MATHEMATICS AND STATISTICS
  12. C                         P.O. BOX 1965
  13. C                         CANBERRA, ACT 2601
  14. C                         AUSTRALIA
  15. C
  16. C   LATEST REVISION     - 15 AUGUST 1985
  17. C
  18. C   PURPOSE             - CUBIC SPLINE DATA SMOOTHER
  19. C
  20. C   USAGE               - CALL CUBGCV (X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
  21. C
  22. C   ARGUMENTS    X      - VECTOR OF LENGTH N CONTAINING THE
  23. C                           ABSCISSAE OF THE N DATA POINTS
  24. C                           (X(I),F(I)) I=1..N. (INPUT) X
  25. C                           MUST BE ORDERED SO THAT
  26. C                           X(I) .LT. X(I+1).
  27. C                F      - VECTOR OF LENGTH N CONTAINING THE
  28. C                           ORDINATES (OR FUNCTION VALUES)
  29. C                           OF THE N DATA POINTS (INPUT).
  30. C                DF     - VECTOR OF LENGTH N. (INPUT/OUTPUT)
  31. C                           DF(I) IS THE RELATIVE STANDARD DEVIATION
  32. C                           OF THE ERROR ASSOCIATED WITH DATA POINT I.
  33. C                           EACH DF(I) MUST BE POSITIVE.  THE VALUES IN
  34. C                           DF ARE SCALED BY THE SUBROUTINE SO THAT
  35. C                           THEIR MEAN SQUARE VALUE IS 1, AND UNSCALED
  36. C                           AGAIN ON NORMAL EXIT.
  37. C                           THE MEAN SQUARE VALUE OF THE DF(I) IS RETURNED
  38. C                           IN WK(7) ON NORMAL EXIT.
  39. C                           IF THE ABSOLUTE STANDARD DEVIATIONS ARE KNOWN,
  40. C                           THESE SHOULD BE PROVIDED IN DF AND THE ERROR
  41. C                           VARIANCE PARAMETER VAR (SEE BELOW) SHOULD THEN
  42. C                           BE SET TO 1.
  43. C                           IF THE RELATIVE STANDARD DEVIATIONS ARE UNKNOWN,
  44. C                           SET EACH DF(I)=1.
  45. C                N      - NUMBER OF DATA POINTS (INPUT).
  46. C                           N MUST BE .GE. 3.
  47. C                Y,C    - SPLINE COEFFICIENTS. (OUTPUT) Y
  48. C                           IS A VECTOR OF LENGTH N. C IS
  49. C                           AN N-1 BY 3 MATRIX. THE VALUE
  50. C                           OF THE SPLINE APPROXIMATION AT T IS
  51. C                           S(T)=((C(I,3)*D+C(I,2))*D+C(I,1))*D+Y(I)
  52. C                           WHERE X(I).LE.T.LT.X(I+1) AND
  53. C                           D = T-X(I).
  54. C                IC     - ROW DIMENSION OF MATRIX C EXACTLY
  55. C                           AS SPECIFIED IN THE DIMENSION
  56. C                           STATEMENT IN THE CALLING PROGRAM. (INPUT)
  57. C                VAR    - ERROR VARIANCE. (INPUT/OUTPUT)
  58. C                           IF VAR IS NEGATIVE (I.E. UNKNOWN) THEN
  59. C                           THE SMOOTHING PARAMETER IS DETERMINED
  60. C                           BY MINIMIZING THE GENERALIZED CROSS VALIDATION
  61. C                           AND AN ESTIMATE OF THE ERROR VARIANCE IS
  62. C                           RETURNED IN VAR.
  63. C                           IF VAR IS NON-NEGATIVE (I.E. KNOWN) THEN THE
  64. C                           SMOOTHING PARAMETER IS DETERMINED TO MINIMIZE
  65. C                           AN ESTIMATE, WHICH DEPENDS ON VAR, OF THE TRUE
  66. C                           MEAN SQUARE ERROR, AND VAR IS UNCHANGED.
  67. C                           IN PARTICULAR, IF VAR IS ZERO, THEN AN
  68. C                           INTERPOLATING NATURAL CUBIC SPLINE IS CALCULATED.
  69. C                           VAR SHOULD BE SET TO 1 IF ABSOLUTE STANDARD
  70. C                           DEVIATIONS HAVE BEEN PROVIDED IN DF (SEE ABOVE).
  71. C                JOB    - JOB SELECTION PARAMETER. (INPUT)
  72. C                         JOB = 0 SHOULD BE SELECTED IF POINT STANDARD ERROR
  73. C                           ESTIMATES ARE NOT REQUIRED IN SE.
  74. C                         JOB = 1 SHOULD BE SELECTED IF POINT STANDARD ERROR
  75. C                           ESTIMATES ARE REQUIRED IN SE.
  76. C                SE     - VECTOR OF LENGTH N CONTAINING BAYESIAN STANDARD
  77. C                           ERROR ESTIMATES OF THE FITTED SPLINE VALUES IN Y.
  78. C                           SE IS NOT REFERENCED IF JOB=0. (OUTPUT)
  79. C                WK     - WORK VECTOR OF LENGTH 7*(N + 2). ON NORMAL EXIT THE
  80. C                           FIRST 7 VALUES OF WK ARE ASSIGNED AS FOLLOWS:-
  81. C
  82. C                           WK(1) = SMOOTHING PARAMETER (= RHO/(RHO + 1))
  83. C                           WK(2) = ESTIMATE OF THE NUMBER OF DEGREES OF
  84. C                                   FREEDOM OF THE RESIDUAL SUM OF SQUARES
  85. C                           WK(3) = GENERALIZED CROSS VALIDATION
  86. C                           WK(4) = MEAN SQUARE RESIDUAL
  87. C                           WK(5) = ESTIMATE OF THE TRUE MEAN SQUARE ERROR
  88. C                                   AT THE DATA POINTS
  89. C                           WK(6) = ESTIMATE OF THE ERROR VARIANCE
  90. C                           WK(7) = MEAN SQUARE VALUE OF THE DF(I)
  91. C
  92. C                           IF WK(1)=0 (RHO=0) AN INTERPOLATING NATURAL CUBIC
  93. C                           SPLINE HAS BEEN CALCULATED.
  94. C                           IF WK(1)=1 (RHO=INFINITE) A LEAST SQUARES
  95. C                           REGRESSION LINE HAS BEEN CALCULATED.
  96. C                           WK(2) IS AN ESTIMATE OF THE NUMBER OF DEGREES OF
  97. C                           FREEDOM OF THE RESIDUAL WHICH REDUCES TO THE
  98. C                           USUAL VALUE OF N-2 WHEN A LEAST SQUARES REGRESSION
  99. C                           LINE IS CALCULATED.
  100. C                           WK(3),WK(4),WK(5) ARE CALCULATED WITH THE DF(I)
  101. C                           SCALED TO HAVE MEAN SQUARE VALUE 1.  THE
  102. C                           UNSCALED VALUES OF WK(3),WK(4),WK(5) MAY BE
  103. C                           CALCULATED BY DIVIDING BY WK(7).
  104. C                           WK(6) COINCIDES WITH THE OUTPUT VALUE OF VAR IF
  105. C                           VAR IS NEGATIVE ON INPUT.  IT IS CALCULATED WITH
  106. C                           THE UNSCALED VALUES OF THE DF(I) TO FACILITATE
  107. C                           COMPARISONS WITH A PRIORI VARIANCE ESTIMATES.
  108. C
  109. C                IER    - ERROR PARAMETER. (OUTPUT)
  110. C                         TERMINAL ERROR
  111. C                           IER = 129, IC IS LESS THAN N-1.
  112. C                           IER = 130, N IS LESS THAN 3.
  113. C                           IER = 131, INPUT ABSCISSAE ARE NOT
  114. C                             ORDERED SO THAT X(I).LT.X(I+1).
  115. C                           IER = 132, DF(I) IS NOT POSITIVE FOR SOME I.
  116. C                           IER = 133, JOB IS NOT 0 OR 1.
  117. C
  118. C   PRECISION/HARDWARE  - DOUBLE
  119. C
  120. C   REQUIRED ROUTINES   - SPINT1,SPFIT1,SPCOF1,SPERR1
  121. C
  122. C   REMARKS      THE NUMBER OF ARITHMETIC OPERATIONS REQUIRED BY THE
  123. C                SUBROUTINE IS PROPORTIONAL TO N.  THE SUBROUTINE
  124. C                USES AN ALGORITHM DEVELOPED BY M.F. HUTCHINSON AND
  125. C                F.R. DE HOOG, 'SMOOTHING NOISY DATA WITH SPLINE
  126. C                FUNCTIONS', NUMER. MATH. (IN PRESS)
  127. C
  128. C-----------------------------------------------------------------------
  129. C
  130.       SUBROUTINE CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
  131. C
  132. C---SPECIFICATIONS FOR ARGUMENTS---
  133.       INTEGER N,IC,JOB,IER
  134.       DOUBLE PRECISION X(N),F(N),DF(N),Y(N),C(IC,3),SE(N),VAR,
  135.      .                 WK(0:N+1,7)
  136. C
  137. C---SPECIFICATIONS FOR LOCAL VARIABLES---
  138.       DOUBLE PRECISION DELTA,ERR,GF1,GF2,GF3,GF4,R1,R2,R3,R4,TAU,RATIO,
  139.      .                 AVH,AVDF,AVAR,ZERO,ONE,STAT(6),P,Q
  140. C
  141.       DATA RATIO/2.0D0/
  142.       DATA TAU/1.618033989D0/
  143.       DATA ZERO,ONE/0.0D0,1.0D0/
  144. C
  145. C---INITIALIZE---
  146.       IER = 133
  147.       IF (JOB.LT.0 .OR. JOB.GT.1) GO TO 140
  148.       CALL SPINT1(X,AVH,F,DF,AVDF,N,Y,C,IC,WK,WK(0,4),IER)
  149.       IF (IER.NE.0) GO TO 140
  150.       AVAR = VAR
  151.       IF (VAR.GT.ZERO) AVAR = VAR*AVDF*AVDF
  152. C
  153. C---CHECK FOR ZERO VARIANCE---
  154.       IF (VAR.NE.ZERO) GO TO 10
  155.       R1 = ZERO
  156.       GO TO 90
  157. C
  158. C---FIND LOCAL MINIMUM OF GCV OR THE EXPECTED MEAN SQUARE ERROR---
  159.    10 R1 = ONE
  160.       R2 = RATIO*R1
  161.       CALL SPFIT1(X,AVH,DF,N,R2,P,Q,GF2,AVAR,STAT,Y,C,IC,WK,WK(0,4),
  162.      .            WK(0,6),WK(0,7))
  163.    20 CALL SPFIT1(X,AVH,DF,N,R1,P,Q,GF1,AVAR,STAT,Y,C,IC,WK,WK(0,4),
  164.      .            WK(0,6),WK(0,7))
  165.       IF (GF1.GT.GF2) GO TO 30
  166. C
  167. C---EXIT IF P ZERO---
  168.       IF (P.LE.ZERO) GO TO 100
  169.       R2 = R1
  170.       GF2 = GF1
  171.       R1 = R1/RATIO
  172.       GO TO 20
  173.  
  174.    30 R3 = RATIO*R2
  175.    40 CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4),
  176.      .            WK(0,6),WK(0,7))
  177.       IF (GF3.GT.GF2) GO TO 50
  178. C
  179. C---EXIT IF Q ZERO---
  180.       IF (Q.LE.ZERO) GO TO 100
  181.       R2 = R3
  182.       GF2 = GF3
  183.       R3 = RATIO*R3
  184.       GO TO 40
  185.  
  186.    50 R2 = R3
  187.       GF2 = GF3
  188.       DELTA = (R2-R1)/TAU
  189.       R4 = R1 + DELTA
  190.       R3 = R2 - DELTA
  191.       CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4),
  192.      .            WK(0,6),WK(0,7))
  193.       CALL SPFIT1(X,AVH,DF,N,R4,P,Q,GF4,AVAR,STAT,Y,C,IC,WK,WK(0,4),
  194.      .            WK(0,6),WK(0,7))
  195. C
  196. C---GOLDEN SECTION SEARCH FOR LOCAL MINIMUM---
  197.    60 IF (GF3.GT.GF4) GO TO 70
  198.       R2 = R4
  199.       GF2 = GF4
  200.       R4 = R3
  201.       GF4 = GF3
  202.       DELTA = DELTA/TAU
  203.       R3 = R2 - DELTA
  204.       CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4),
  205.      .            WK(0,6),WK(0,7))
  206.       GO TO 80
  207.  
  208.    70 R1 = R3
  209.       GF1 = GF3
  210.       R3 = R4
  211.       GF3 = GF4
  212.       DELTA = DELTA/TAU
  213.       R4 = R1 + DELTA
  214.       CALL SPFIT1(X,AVH,DF,N,R4,P,Q,GF4,AVAR,STAT,Y,C,IC,WK,WK(0,4),
  215.      .            WK(0,6),WK(0,7))
  216.    80 ERR = (R2-R1)/ (R1+R2)
  217.       IF (ERR*ERR+ONE.GT.ONE .AND. ERR.GT.1.0D-6) GO TO 60
  218.       R1 = (R1+R2)*0.5D0
  219. C
  220. C---CALCULATE SPLINE COEFFICIENTS---
  221.    90 CALL SPFIT1(X,AVH,DF,N,R1,P,Q,GF1,AVAR,STAT,Y,C,IC,WK,WK(0,4),
  222.      .            WK(0,6),WK(0,7))
  223.   100 CALL SPCOF1(X,AVH,F,DF,N,P,Q,Y,C,IC,WK(0,6),WK(0,7))
  224. C
  225. C---OPTIONALLY CALCULATE STANDARD ERROR ESTIMATES---
  226.       IF (VAR.GE.ZERO) GO TO 110
  227.       AVAR = STAT(6)
  228.       VAR = AVAR/ (AVDF*AVDF)
  229.   110 IF (JOB.EQ.1) CALL SPERR1(X,AVH,DF,N,WK,P,AVAR,SE)
  230. C
  231. C---UNSCALE DF---
  232.       DO 120 I = 1,N
  233.          DF(I) = DF(I)*AVDF
  234.   120 CONTINUE
  235. C
  236. C--PUT STATISTICS IN WK---
  237.       DO 130 I = 0,5
  238.          WK(I,1) = STAT(I+1)
  239.   130 CONTINUE
  240.       WK(5,1) = STAT(6)/ (AVDF*AVDF)
  241.       WK(6,1) = AVDF*AVDF
  242.       GO TO 150
  243. C
  244. C---CHECK FOR ERROR CONDITION---
  245.   140 CONTINUE
  246. C     IF (IER.NE.0) CONTINUE
  247.   150 RETURN
  248.       END
  249.       SUBROUTINE SPINT1(X,AVH,Y,DY,AVDY,N,A,C,IC,R,T,IER)
  250. C
  251. C INITIALIZES THE ARRAYS C, R AND T FOR ONE DIMENSIONAL CUBIC
  252. C SMOOTHING SPLINE FITTING BY SUBROUTINE SPFIT1.  THE VALUES
  253. C DF(I) ARE SCALED SO THAT THE SUM OF THEIR SQUARES IS N
  254. C AND THE AVERAGE OF THE DIFFERENCES X(I+1) - X(I) IS CALCULATED
  255. C IN AVH IN ORDER TO AVOID UNDERFLOW AND OVERFLOW PROBLEMS IN
  256. C SPFIT1.
  257. C
  258. C SUBROUTINE SETS IER IF ELEMENTS OF X ARE NON-INCREASING,
  259. C IF N IS LESS THAN 3, IF IC IS LESS THAN N-1 OR IF DY(I) IS
  260. C NOT POSITIVE FOR SOME I.
  261. C
  262. C---SPECIFICATIONS FOR ARGUMENTS---
  263.       INTEGER N,IC,IER
  264.       DOUBLE PRECISION X(N),Y(N),DY(N),A(N),C(IC,3),R(0:N+1,3),
  265.      .                 T(0:N+1,2),AVH,AVDY
  266. C
  267. C---SPECIFICATIONS FOR LOCAL VARIABLES---
  268.       INTEGER I
  269.       DOUBLE PRECISION E,F,G,H,ZERO
  270.       DATA ZERO/0.0D0/
  271. C
  272. C---INITIALIZATION AND INPUT CHECKING---
  273.       IER = 0
  274.       IF (N.LT.3) GO TO 60
  275.       IF (IC.LT.N-1) GO TO 70
  276. C
  277. C---GET AVERAGE X SPACING IN AVH---
  278.       G = ZERO
  279.       DO 10 I = 1,N - 1
  280.          H = X(I+1) - X(I)
  281.          IF (H.LE.ZERO) GO TO 80
  282.          G = G + H
  283.    10 CONTINUE
  284.       AVH = G/ (N-1)
  285. C
  286. C---SCALE RELATIVE WEIGHTS---
  287.       G = ZERO
  288.       DO 20 I = 1,N
  289.          IF (DY(I).LE.ZERO) GO TO 90
  290.          G = G + DY(I)*DY(I)
  291.    20 CONTINUE
  292.       AVDY = DSQRT(G/N)
  293. C
  294.       DO 30 I = 1,N
  295.          DY(I) = DY(I)/AVDY
  296.    30 CONTINUE
  297. C
  298. C---INITIALIZE H,F---
  299.       H = (X(2)-X(1))/AVH
  300.       F = (Y(2)-Y(1))/H
  301. C
  302. C---CALCULATE A,T,R---
  303.       DO 40 I = 2,N - 1
  304.          G = H
  305.          H = (X(I+1)-X(I))/AVH
  306.          E = F
  307.          F = (Y(I+1)-Y(I))/H
  308.          A(I) = F - E
  309.          T(I,1) = 2.0D0* (G+H)/3.0D0
  310.          T(I,2) = H/3.0D0
  311.          R(I,3) = DY(I-1)/G
  312.          R(I,1) = DY(I+1)/H
  313.          R(I,2) = -DY(I)/G - DY(I)/H
  314.    40 CONTINUE
  315. C
  316. C---CALCULATE C = R'*R---
  317.       R(N,2) = ZERO
  318.       R(N,3) = ZERO
  319.       R(N+1,3) = ZERO
  320.       DO 50 I = 2,N - 1
  321.          C(I,1) = R(I,1)*R(I,1) + R(I,2)*R(I,2) + R(I,3)*R(I,3)
  322.          C(I,2) = R(I,1)*R(I+1,2) + R(I,2)*R(I+1,3)
  323.          C(I,3) = R(I,1)*R(I+2,3)
  324.    50 CONTINUE
  325.       RETURN
  326. C
  327. C---ERROR CONDITIONS---
  328.    60 IER = 130
  329.       RETURN
  330.  
  331.    70 IER = 129
  332.       RETURN
  333.  
  334.    80 IER = 131
  335.       RETURN
  336.  
  337.    90 IER = 132
  338.       RETURN
  339.       END
  340.       SUBROUTINE SPFIT1(X,AVH,DY,N,RHO,P,Q,FUN,VAR,STAT,A,C,IC,R,T,U,V)
  341. C
  342. C FITS A CUBIC SMOOTHING SPLINE TO DATA WITH RELATIVE
  343. C WEIGHTING DY FOR A GIVEN VALUE OF THE SMOOTHING PARAMETER
  344. C RHO USING AN ALGORITHM BASED ON THAT OF C.H. REINSCH (1967),
  345. C NUMER. MATH. 10, 177-183.
  346. C
  347. C THE TRACE OF THE INFLUENCE MATRIX IS CALCULATED USING AN
  348. C ALGORITHM DEVELOPED BY M.F.HUTCHINSON AND F.R.DE HOOG (NUMER.
  349. C MATH., IN PRESS), ENABLING THE GENERALIZED CROSS VALIDATION
  350. C AND RELATED STATISTICS TO BE CALCULATED IN ORDER N OPERATIONS.
  351. C
  352. C THE ARRAYS A, C, R AND T ARE ASSUMED TO HAVE BEEN INITIALIZED
  353. C BY THE SUBROUTINE SPINT1.  OVERFLOW AND UNDERFLOW PROBLEMS ARE
  354. C AVOIDED BY USING P=RHO/(1 + RHO) AND Q=1/(1 + RHO) INSTEAD OF
  355. C RHO AND BY SCALING THE DIFFERENCES X(I+1) - X(I) BY AVH.
  356. C
  357. C THE VALUES IN DF ARE ASSUMED TO HAVE BEEN SCALED SO THAT THE
  358. C SUM OF THEIR SQUARED VALUES IS N.  THE VALUE IN VAR, WHEN IT IS
  359. C NON-NEGATIVE, IS ASSUMED TO HAVE BEEN SCALED TO COMPENSATE FOR
  360. C THE SCALING OF THE VALUES IN DF.
  361. C
  362. C THE VALUE RETURNED IN FUN IS AN ESTIMATE OF THE TRUE MEAN SQUARE
  363. C WHEN VAR IS NON-NEGATIVE, AND IS THE GENERALIZED CROSS VALIDATION
  364. C WHEN VAR IS NEGATIVE.
  365. C
  366. C---SPECIFICATIONS FOR ARGUMENTS---
  367.       INTEGER IC,N
  368.       DOUBLE PRECISION X(N),DY(N),RHO,STAT(6),A(N),C(IC,3),R(0:N+1,3),
  369.      .                 T(0:N+1,2),U(0:N+1),V(0:N+1),FUN,VAR,AVH,P,Q
  370. C
  371. C---LOCAL VARIABLES---
  372.       INTEGER I
  373.       DOUBLE PRECISION E,F,G,H,ZERO,ONE,TWO,RHO1
  374.       DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/
  375. C
  376. C---USE P AND Q INSTEAD OF RHO TO PREVENT OVERFLOW OR UNDERFLOW---
  377.       RHO1 = ONE + RHO
  378.       P = RHO/RHO1
  379.       Q = ONE/RHO1
  380.       IF (RHO1.EQ.ONE) P = ZERO
  381.       IF (RHO1.EQ.RHO) Q = ZERO
  382. C
  383. C---RATIONAL CHOLESKY DECOMPOSITION OF P*C + Q*T---
  384.       F = ZERO
  385.       G = ZERO
  386.       H = ZERO
  387.       DO 10 I = 0,1
  388.          R(I,1) = ZERO
  389.    10 CONTINUE
  390.       DO 20 I = 2,N - 1
  391.          R(I-2,3) = G*R(I-2,1)
  392.          R(I-1,2) = F*R(I-1,1)
  393.          R(I,1) = ONE/ (P*C(I,1)+Q*T(I,1)-F*R(I-1,2)-G*R(I-2,3))
  394.          F = P*C(I,2) + Q*T(I,2) - H*R(I-1,2)
  395.          G = H
  396.          H = P*C(I,3)
  397.    20 CONTINUE
  398. C
  399. C---SOLVE FOR U---
  400.       U(0) = ZERO
  401.       U(1) = ZERO
  402.       DO 30 I = 2,N - 1
  403.          U(I) = A(I) - R(I-1,2)*U(I-1) - R(I-2,3)*U(I-2)
  404.    30 CONTINUE
  405.       U(N) = ZERO
  406.       U(N+1) = ZERO
  407.       DO 40 I = N - 1,2,-1
  408.          U(I) = R(I,1)*U(I) - R(I,2)*U(I+1) - R(I,3)*U(I+2)
  409.    40 CONTINUE
  410. C
  411. C---CALCULATE RESIDUAL VECTOR V---
  412.       E = ZERO
  413.       H = ZERO
  414.       DO 50 I = 1,N - 1
  415.          G = H
  416.          H = (U(I+1)-U(I))/ ((X(I+1)-X(I))/AVH)
  417.          V(I) = DY(I)* (H-G)
  418.          E = E + V(I)*V(I)
  419.    50 CONTINUE
  420.       V(N) = DY(N)* (-H)
  421.       E = E + V(N)*V(N)
  422. C
  423. C---CALCULATE UPPER THREE BANDS OF INVERSE MATRIX---
  424.       R(N,1) = ZERO
  425.       R(N,2) = ZERO
  426.       R(N+1,1) = ZERO
  427.       DO 60 I = N - 1,2,-1
  428.          G = R(I,2)
  429.          H = R(I,3)
  430.          R(I,2) = -G*R(I+1,1) - H*R(I+1,2)
  431.          R(I,3) = -G*R(I+1,2) - H*R(I+2,1)
  432.          R(I,1) = R(I,1) - G*R(I,2) - H*R(I,3)
  433.    60 CONTINUE
  434. C
  435. C---CALCULATE TRACE---
  436.       F = ZERO
  437.       G = ZERO
  438.       H = ZERO
  439.       DO 70 I = 2,N - 1
  440.          F = F + R(I,1)*C(I,1)
  441.          G = G + R(I,2)*C(I,2)
  442.          H = H + R(I,3)*C(I,3)
  443.    70 CONTINUE
  444.       F = F + TWO* (G+H)
  445. C
  446. C---CALCULATE STATISTICS---
  447.       STAT(1) = P
  448.       STAT(2) = F*P
  449.       STAT(3) = N*E/ (F*F)
  450.       STAT(4) = E*P*P/N
  451.       STAT(6) = E*P/F
  452.       IF (VAR.GE.ZERO) GO TO 80
  453.       STAT(5) = STAT(6) - STAT(4)
  454.       FUN = STAT(3)
  455.       GO TO 90
  456.  
  457.    80 STAT(5) = DMAX1(STAT(4)-TWO*VAR*STAT(2)/N+VAR,ZERO)
  458.       FUN = STAT(5)
  459.    90 RETURN
  460.       END
  461.       SUBROUTINE SPERR1(X,AVH,DY,N,R,P,VAR,SE)
  462. C
  463. C CALCULATES BAYESIAN ESTIMATES OF THE STANDARD ERRORS OF THE FITTED
  464. C VALUES OF A CUBIC SMOOTHING SPLINE BY CALCULATING THE DIAGONAL ELEMENTS
  465. C OF THE INFLUENCE MATRIX.
  466. C
  467. C---SPECIFICATIONS FOR ARGUMENTS---
  468.       INTEGER N
  469.       DOUBLE PRECISION X(N),DY(N),R(0:N+1,3),SE(N),AVH,P,VAR
  470. C
  471. C---SPECIFICATIONS FOR LOCAL VARIABLES---
  472.       INTEGER I
  473.       DOUBLE PRECISION F,G,H,F1,G1,H1,ZERO,ONE
  474.       DATA ZERO,ONE/0.0D0,1.0D0/
  475. C
  476. C---INITIALIZE---
  477.       H = AVH/ (X(2)-X(1))
  478.       SE(1) = ONE - P*DY(1)*DY(1)*H*H*R(2,1)
  479.       R(1,1) = ZERO
  480.       R(1,2) = ZERO
  481.       R(1,3) = ZERO
  482. C
  483. C---CALCULATE DIAGONAL ELEMENTS---
  484.       DO 10 I = 2,N - 1
  485.          F = H
  486.          H = AVH/ (X(I+1)-X(I))
  487.          G = -F - H
  488.          F1 = F*R(I-1,1) + G*R(I-1,2) + H*R(I-1,3)
  489.          G1 = F*R(I-1,2) + G*R(I,1) + H*R(I,2)
  490.          H1 = F*R(I-1,3) + G*R(I,2) + H*R(I+1,1)
  491.          SE(I) = ONE - P*DY(I)*DY(I)* (F*F1+G*G1+H*H1)
  492.    10 CONTINUE
  493.       SE(N) = ONE - P*DY(N)*DY(N)*H*H*R(N-1,1)
  494. C
  495. C---CALCULATE STANDARD ERROR ESTIMATES---
  496.       DO 20 I = 1,N
  497.          SE(I) = DSQRT(DMAX1(SE(I)*VAR,ZERO))*DY(I)
  498.    20 CONTINUE
  499.       RETURN
  500.       END
  501.       SUBROUTINE SPCOF1(X,AVH,Y,DY,N,P,Q,A,C,IC,U,V)
  502. C
  503. C CALCULATES COEFFICIENTS OF A CUBIC SMOOTHING SPLINE FROM
  504. C PARAMETERS CALCULATED BY SUBROUTINE SPFIT1.
  505. C
  506. C---SPECIFICATIONS FOR ARGUMENTS---
  507.       INTEGER IC,N
  508.       DOUBLE PRECISION X(N),Y(N),DY(N),P,Q,A(N),C(IC,3),U(0:N+1),
  509.      .                 V(0:N+1),AVH
  510. C
  511. C---SPECIFICATIONS FOR LOCAL VARIABLES---
  512.       INTEGER I
  513.       DOUBLE PRECISION H,QH
  514. C
  515. C---CALCULATE A---
  516.       QH = Q/ (AVH*AVH)
  517.       DO 10 I = 1,N
  518.          A(I) = Y(I) - P*DY(I)*V(I)
  519.          U(I) = QH*U(I)
  520.    10 CONTINUE
  521. C
  522. C---CALCULATE C---
  523.       DO 20 I = 1,N - 1
  524.          H = X(I+1) - X(I)
  525.          C(I,3) = (U(I+1)-U(I))/ (3.0D0*H)
  526.          C(I,1) = (A(I+1)-A(I))/H - (H*C(I,3)+U(I))*H
  527.          C(I,2) = U(I)
  528.    20 CONTINUE
  529.       RETURN
  530.       END
  531. C CUBGCV TEST DRIVER
  532. C ------------------
  533. C
  534. C AUTHOR          - M.F.HUTCHINSON
  535. C                   CSIRO DIVISION OF WATER AND LAND RESOURCES
  536. C                   GPO BOX 1666
  537. C                   CANBERRA ACT 2601
  538. C                   AUSTRALIA
  539. C
  540. C LATEST REVISION - 7 AUGUST 1986
  541. C
  542. C COMPUTER        - VAX/DOUBLE
  543. C
  544. C USAGE           - MAIN PROGRAM
  545. C
  546. C REQUIRED ROUTINES - CUBGCV,SPINT1,SPFIT1,SPCOF1,SPERR1,GGRAND
  547. C
  548. C REMARKS   USES SUBROUTINE CUBGCV TO FIT A CUBIC SMOOTHING SPLINE
  549. C           TO 50 DATA POINTS WHICH ARE GENERATED BY ADDING A RANDOM
  550. C           VARIABLE WITH UNIFORM DENSITY IN THE INTERVAL [-0.3,0.3]
  551. C           TO 50 POINTS SAMPLED FROM THE CURVE  Y=SIN(3*PI*X/2).
  552. C           RANDOM DEVIATES IN THE INTERVAL [0,1] ARE GENERATED BY THE
  553. C           DOUBLE PRECISION FUNCTION GGRAND (SIMILAR TO IMSL FUNCTION
  554. C           GGUBFS).  THE ABSCISSAE ARE UNEQUALLY SPACED IN [0,1].
  555. C
  556. C           POINT STANDARD ERROR ESTIMATES ARE RETURNED IN SE BY
  557. C           SETTING JOB=1.  THE ERROR VARIANCE ESTIMATE IS RETURNED
  558. C           IN VAR.  IT COMPARES FAVOURABLY WITH THE TRUE VALUE OF 0.03.
  559. C           SUMMARY STATISTICS FROM THE ARRAY WK ARE WRITTEN TO
  560. C           UNIT 6.  DATA VALUES AND FITTED VALUES WITH ESTIMATED
  561. C           STANDARD ERRORS ARE ALSO WRITTEN TO UNIT 6.
  562. C
  563.       PARAMETER (N=50, IC=49)
  564. C
  565.       INTEGER            JOB,IER
  566.       DOUBLE PRECISION   X(N),F(N),Y(N),DF(N),C(IC,3),WK(7*(N+2)),
  567.      *                   VAR,SE(N)
  568.       DOUBLE PRECISION   GGRAND,DSEED
  569. C
  570. C---INITIALIZE---
  571.       DSEED=1.2345D4
  572.       JOB=1
  573.       VAR=-1.0D0
  574. C
  575. C---CALCULATE DATA POINTS---
  576.       DO 10 I=1,N
  577.       X(I)=(I - 0.5)/N + (2.0*GGRAND(DSEED) - 1.0)/(3.0*N)
  578.       F(I)=DSIN(4.71238*X(I)) + (2.0*GGRAND(DSEED) - 1.0)*0.3
  579.       DF(I)=1.0D0
  580.   10  CONTINUE
  581. C
  582. C---FIT CUBIC SPLINE---
  583.       CALL CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
  584. C
  585. C---WRITE OUT RESULTS---
  586.       WRITE(6,20)
  587.   20  FORMAT(' CUBGCV TEST DRIVER RESULTS:')
  588.       WRITE(6,30) IER,VAR,WK(3),WK(4),WK(2)
  589.   30  FORMAT(/' IER =',I4/' VAR =',F7.4/
  590.      *        ' GENERALIZED CROSS VALIDATION =',F7.4/
  591.      *        ' MEAN SQUARE RESIDUAL         =',F7.4/
  592.      *        ' RESIDUAL DEGREES OF FREEDOM  =',F7.2)
  593.       WRITE(6,40)
  594.   40  FORMAT(/' INPUT DATA',17X,'OUTPUT RESULTS'//
  595.      *         '   I    X(I)    F(I)',6X,'    Y(I)   SE(I)',
  596.      *          '      C(I,1)      C(I,2)      C(I,3)')
  597.       DO 60 I=1,N-1
  598.       WRITE(6,50) I,X(I),F(I),Y(I),SE(I),(C(I,J),J=1,3)
  599.   50  FORMAT(I4,2F8.4,6X,2F8.4,3E12.4)
  600.   60  CONTINUE
  601.       WRITE(6,50) N,X(N),F(N),Y(N),SE(N)
  602.       STOP
  603.       END
  604.       DOUBLE PRECISION FUNCTION GGRAND(DSEED)
  605. C
  606. C DOUBLE PRECISION UNIFORM RANDOM NUMBER GENERATOR
  607. C
  608. C CONSTANTS: A = 7**5
  609. C            B = 2**31 - 1
  610. C            C = 2**31
  611. C
  612. C REFERENCE: IMSL MANUAL, CHAPTER G - GENERATION AND TESTING OF
  613. C                                     RANDOM NUMBERS
  614. C
  615. C---SPECIFICATIONS FOR ARGUMENTS---
  616.       DOUBLE PRECISION DSEED
  617. C
  618. C---SPECIFICATIONS FOR LOCAL VARIABLES---
  619.       DOUBLE PRECISION A,B,C,S
  620. C
  621.       DATA A,B,C/16807.0D0, 2147483647.0D0, 2147483648.0D0/
  622. C
  623.       S=DSEED
  624.       S=DMOD(A*S, B)
  625.       GGRAND=S/C
  626.       DSEED=S
  627.       RETURN
  628.       END
  629.  
  630. CUBGCV TEST DRIVER RESULTS:
  631.  
  632. IER =   0
  633. VAR = 0.0279
  634. GENERALIZED CROSS VALIDATION = 0.0318
  635. MEAN SQUARE RESIDUAL         = 0.0246
  636. RESIDUAL DEGREES OF FREEDOM  =  43.97
  637.  
  638. INPUT DATA                 OUTPUT RESULTS
  639.  
  640.   I    X(I)    F(I)          Y(I)   SE(I)      C(I,1)      C(I,2)      C(I,3)
  641.   1  0.0046  0.2222        0.0342  0.1004  0.3630E+01  0.0000E+00  0.2542E+02
  642.   2  0.0360 -0.1098        0.1488  0.0750  0.3705E+01  0.2391E+01 -0.9537E+01
  643.   3  0.0435 -0.0658        0.1767  0.0707  0.3740E+01  0.2175E+01 -0.4233E+02
  644.   4  0.0735  0.3906        0.2900  0.0594  0.3756E+01 -0.1642E+01 -0.2872E+02
  645.   5  0.0955  0.6054        0.3714  0.0558  0.3642E+01 -0.3535E+01  0.2911E+01
  646.   6  0.1078  0.3034        0.4155  0.0549  0.3557E+01 -0.3428E+01 -0.1225E+02
  647.   7  0.1269  0.7386        0.4822  0.0544  0.3412E+01 -0.4131E+01  0.2242E+02
  648.   8  0.1565  0.4616        0.5800  0.0543  0.3227E+01 -0.2143E+01  0.6415E+01
  649.   9  0.1679  0.4315        0.6165  0.0543  0.3180E+01 -0.1923E+01 -0.1860E+02
  650.  10  0.1869  0.5716        0.6762  0.0544  0.3087E+01 -0.2985E+01 -0.3274E+02
  651.  11  0.2149  0.6736        0.7595  0.0542  0.2843E+01 -0.5733E+01 -0.4435E+02
  652.  12  0.2356  0.7388        0.8155  0.0539  0.2549E+01 -0.8486E+01 -0.5472E+02
  653.  13  0.2557  1.1953        0.8630  0.0537  0.2139E+01 -0.1180E+02 -0.9784E+01
  654.  14  0.2674  1.0299        0.8864  0.0536  0.1860E+01 -0.1214E+02  0.9619E+01
  655.  15  0.2902  0.7981        0.9225  0.0534  0.1322E+01 -0.1149E+02 -0.7202E+01
  656.  16  0.3155  0.8973        0.9485  0.0532  0.7269E+00 -0.1203E+02 -0.1412E+02
  657.  17  0.3364  1.2695        0.9583  0.0530  0.2040E+00 -0.1292E+02  0.2796E+02
  658.  18  0.3557  0.7253        0.9577  0.0527 -0.2638E+00 -0.1130E+02 -0.3453E+01
  659.  19  0.3756  1.2127        0.9479  0.0526 -0.7176E+00 -0.1151E+02  0.3235E+02
  660.  20  0.3881  0.7304        0.9373  0.0525 -0.9889E+00 -0.1030E+02  0.4381E+01
  661.  21  0.4126  0.9810        0.9069  0.0525 -0.1486E+01 -0.9977E+01  0.1440E+02
  662.  22  0.4266  0.7117        0.8842  0.0525 -0.1756E+01 -0.9373E+01 -0.8925E+01
  663.  23  0.4566  0.7203        0.8227  0.0524 -0.2344E+01 -0.1018E+02 -0.2278E+02
  664.  24  0.4704  0.9242        0.7884  0.0524 -0.2637E+01 -0.1112E+02 -0.4419E+01
  665.  25  0.4914  0.7345        0.7281  0.0523 -0.3110E+01 -0.1140E+02 -0.3562E+01
  666.  26  0.5084  0.7378        0.6720  0.0524 -0.3500E+01 -0.1158E+02  0.5336E+01
  667.  27  0.5277  0.7441        0.6002  0.0525 -0.3941E+01 -0.1127E+02  0.2479E+02
  668.  28  0.5450  0.5612        0.5286  0.0527 -0.4310E+01 -0.9980E+01  0.2920E+02
  669.  29  0.5641  0.5049        0.4429  0.0529 -0.4659E+01 -0.8309E+01  0.3758E+02
  670.  30  0.5857  0.4725        0.3390  0.0531 -0.4964E+01 -0.5878E+01  0.5563E+02
  671.  31  0.6159  0.1380        0.1850  0.0531 -0.5167E+01 -0.8307E+00  0.4928E+02
  672.  32  0.6317  0.1412        0.1036  0.0531 -0.5157E+01  0.1499E+01  0.5437E+02
  673.  33  0.6446 -0.1110        0.0371  0.0531 -0.5091E+01  0.3614E+01  0.3434E+02
  674.  34  0.6707 -0.2605       -0.0927  0.0532 -0.4832E+01  0.6302E+01  0.1164E+02
  675.  35  0.6853 -0.1284       -0.1619  0.0533 -0.4640E+01  0.6812E+01  0.1617E+02
  676.  36  0.7064 -0.3452       -0.2564  0.0536 -0.4332E+01  0.7834E+01  0.4164E+01
  677.  37  0.7310 -0.5527       -0.3582  0.0538 -0.3939E+01  0.8141E+01 -0.2214E+02
  678.  38  0.7531 -0.3459       -0.4415  0.0540 -0.3611E+01  0.6674E+01 -0.9205E+01
  679.  39  0.7686 -0.5902       -0.4961  0.0541 -0.3410E+01  0.6245E+01 -0.2193E+02
  680.  40  0.7952 -0.7644       -0.5828  0.0541 -0.3125E+01  0.4494E+01 -0.4649E+02
  681.  41  0.8087 -0.5392       -0.6242  0.0541 -0.3029E+01  0.2614E+01 -0.3499E+02
  682.  42  0.8352 -0.4247       -0.7031  0.0539 -0.2964E+01 -0.1603E+00  0.2646E+01
  683.  43  0.8501 -0.6327       -0.7476  0.0538 -0.2967E+01 -0.4132E-01  0.1817E+02
  684.  44  0.8726 -0.9983       -0.8139  0.0538 -0.2942E+01  0.1180E+01 -0.6774E+01
  685.  45  0.8874 -0.9082       -0.8574  0.0542 -0.2911E+01  0.8778E+00 -0.1364E+02
  686.  46  0.9139 -0.8930       -0.9340  0.0566 -0.2893E+01 -0.2044E+00 -0.8094E+01
  687.  47  0.9271 -1.0233       -0.9723  0.0593 -0.2903E+01 -0.5258E+00 -0.1498E+02
  688.  48  0.9473 -0.8839       -1.0313  0.0665 -0.2942E+01 -0.1433E+01  0.4945E+01
  689.  49  0.9652 -1.0172       -1.0843  0.0766 -0.2989E+01 -0.1168E+01  0.1401E+02
  690.  50  0.9930 -1.2715       -1.1679  0.0998
  691. Documentation:
  692. C   COMPUTER            - VAX/DOUBLE
  693. C
  694. C   AUTHOR              - M.F.HUTCHINSON
  695. C                         CSIRO DIVISION OF MATHEMATICS AND STATISTICS
  696. C                         P.O. BOX 1965
  697. C                         CANBERRA, ACT 2601
  698. C                         AUSTRALIA
  699. C
  700. C   LATEST REVISION     - 15 AUGUST 1985
  701. C
  702. C   PURPOSE             - CUBIC SPLINE DATA SMOOTHER
  703. C
  704. C   USAGE               - CALL CUBGCV (X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
  705. C
  706. C   ARGUMENTS    X      - VECTOR OF LENGTH N CONTAINING THE
  707. C                           ABSCISSAE OF THE N DATA POINTS
  708. C                           (X(I),F(I)) I=1..N. (INPUT) X
  709. C                           MUST BE ORDERED SO THAT
  710. C                           X(I) .LT. X(I+1).
  711. C                F      - VECTOR OF LENGTH N CONTAINING THE
  712. C                           ORDINATES (OR FUNCTION VALUES)
  713. C                           OF THE N DATA POINTS (INPUT).
  714. C                DF     - VECTOR OF LENGTH N. (INPUT/OUTPUT)
  715. C                           DF(I) IS THE RELATIVE STANDARD DEVIATION
  716. C                           OF THE ERROR ASSOCIATED WITH DATA POINT I.
  717. C                           EACH DF(I) MUST BE POSITIVE.  THE VALUES IN
  718. C                           DF ARE SCALED BY THE SUBROUTINE SO THAT
  719. C                           THEIR MEAN SQUARE VALUE IS 1, AND UNSCALED
  720. C                           AGAIN ON NORMAL EXIT.
  721. C                           THE MEAN SQUARE VALUE OF THE DF(I) IS RETURNED
  722. C                           IN WK(7) ON NORMAL EXIT.
  723. C                           IF THE ABSOLUTE STANDARD DEVIATIONS ARE KNOWN,
  724. C                           THESE SHOULD BE PROVIDED IN DF AND THE ERROR
  725. C                           VARIANCE PARAMETER VAR (SEE BELOW) SHOULD THEN
  726. C                           BE SET TO 1.
  727. C                           IF THE RELATIVE STANDARD DEVIATIONS ARE UNKNOWN,
  728. C                           SET EACH DF(I)=1.
  729. C                N      - NUMBER OF DATA POINTS (INPUT).
  730. C                           N MUST BE .GE. 3.
  731. C                Y,C    - SPLINE COEFFICIENTS. (OUTPUT) Y
  732. C                           IS A VECTOR OF LENGTH N. C IS
  733. C                           AN N-1 BY 3 MATRIX. THE VALUE
  734. C                           OF THE SPLINE APPROXIMATION AT T IS
  735. C                           S(T)=((C(I,3)*D+C(I,2))*D+C(I,1))*D+Y(I)
  736. C                           WHERE X(I).LE.T.LT.X(I+1) AND
  737. C                           D = T-X(I).
  738. C                IC     - ROW DIMENSION OF MATRIX C EXACTLY
  739. C                           AS SPECIFIED IN THE DIMENSION
  740. C                           STATEMENT IN THE CALLING PROGRAM. (INPUT)
  741. C                VAR    - ERROR VARIANCE. (INPUT/OUTPUT)
  742. C                           IF VAR IS NEGATIVE (I.E. UNKNOWN) THEN
  743. C                           THE SMOOTHING PARAMETER IS DETERMINED
  744. C                           BY MINIMIZING THE GENERALIZED CROSS VALIDATION
  745. C                           AND AN ESTIMATE OF THE ERROR VARIANCE IS
  746. C                           RETURNED IN VAR.
  747. C                           IF VAR IS NON-NEGATIVE (I.E. KNOWN) THEN THE
  748. C                           SMOOTHING PARAMETER IS DETERMINED TO MINIMIZE
  749. C                           AN ESTIMATE, WHICH DEPENDS ON VAR, OF THE TRUE
  750. C                           MEAN SQUARE ERROR, AND VAR IS UNCHANGED.
  751. C                           IN PARTICULAR, IF VAR IS ZERO, THEN AN
  752. C                           INTERPOLATING NATURAL CUBIC SPLINE IS CALCULATED.
  753. C                           VAR SHOULD BE SET TO 1 IF ABSOLUTE STANDARD
  754. C                           DEVIATIONS HAVE BEEN PROVIDED IN DF (SEE ABOVE).
  755. C                JOB    - JOB SELECTION PARAMETER. (INPUT)
  756. C                         JOB = 0 SHOULD BE SELECTED IF POINT STANDARD ERROR
  757. C                           ESTIMATES ARE NOT REQUIRED IN SE.
  758. C                         JOB = 1 SHOULD BE SELECTED IF POINT STANDARD ERROR
  759. C                           ESTIMATES ARE REQUIRED IN SE.
  760. C                SE     - VECTOR OF LENGTH N CONTAINING BAYESIAN STANDARD
  761. C                           ERROR ESTIMATES OF THE FITTED SPLINE VALUES IN Y.
  762. C                           SE IS NOT REFERENCED IF JOB=0. (OUTPUT)
  763. C                WK     - WORK VECTOR OF LENGTH 7*(N + 2). ON NORMAL EXIT THE
  764. C                           FIRST 7 VALUES OF WK ARE ASSIGNED AS FOLLOWS:-
  765. C
  766. C                           WK(1) = SMOOTHING PARAMETER (= RHO/(RHO + 1))
  767. C                           WK(2) = ESTIMATE OF THE NUMBER OF DEGREES OF
  768. C                                   FREEDOM OF THE RESIDUAL SUM OF SQUARES
  769. C                           WK(3) = GENERALIZED CROSS VALIDATION
  770. C                           WK(4) = MEAN SQUARE RESIDUAL
  771. C                           WK(5) = ESTIMATE OF THE TRUE MEAN SQUARE ERROR
  772. C                                   AT THE DATA POINTS
  773. C                           WK(6) = ESTIMATE OF THE ERROR VARIANCE
  774. C                           WK(7) = MEAN SQUARE VALUE OF THE DF(I)
  775. C
  776. C                           IF WK(1)=0 (RHO=0) AN INTERPOLATING NATURAL CUBIC
  777. C                           SPLINE HAS BEEN CALCULATED.
  778. C                           IF WK(1)=1 (RHO=INFINITE) A LEAST SQUARES
  779. C                           REGRESSION LINE HAS BEEN CALCULATED.
  780. C                           WK(2) IS AN ESTIMATE OF THE NUMBER OF DEGREES OF
  781. C                           FREEDOM OF THE RESIDUAL WHICH REDUCES TO THE
  782. C                           USUAL VALUE OF N-2 WHEN A LEAST SQUARES REGRESSION
  783. C                           LINE IS CALCULATED.
  784. C                           WK(3),WK(4),WK(5) ARE CALCULATED WITH THE DF(I)
  785. C                           SCALED TO HAVE MEAN SQUARE VALUE 1.  THE
  786. C                           UNSCALED VALUES OF WK(3),WK(4),WK(5) MAY BE
  787. C                           CALCULATED BY DIVIDING BY WK(7).
  788. C                           WK(6) COINCIDES WITH THE OUTPUT VALUE OF VAR IF
  789. C                           VAR IS NEGATIVE ON INPUT.  IT IS CALCULATED WITH
  790. C                           THE UNSCALED VALUES OF THE DF(I) TO FACILITATE
  791. C                           COMPARISONS WITH A PRIORI VARIANCE ESTIMATES.
  792. C
  793. C                IER    - ERROR PARAMETER. (OUTPUT)
  794. C                         TERMINAL ERROR
  795. C                           IER = 129, IC IS LESS THAN N-1.
  796. C                           IER = 130, N IS LESS THAN 3.
  797. C                           IER = 131, INPUT ABSCISSAE ARE NOT
  798. C                             ORDERED SO THAT X(I).LT.X(I+1).
  799. C                           IER = 132, DF(I) IS NOT POSITIVE FOR SOME I.
  800. C                           IER = 133, JOB IS NOT 0 OR 1.
  801. C
  802. C   PRECISION/HARDWARE  - DOUBLE
  803. C
  804. C   REQUIRED ROUTINES   - SPINT1,SPFIT1,SPCOF1,SPERR1
  805. C
  806. C   REMARKS      THE NUMBER OF ARITHMETIC OPERATIONS REQUIRED BY THE
  807. C                SUBROUTINE IS PROPORTIONAL TO N.  THE SUBROUTINE
  808. C                USES AN ALGORITHM DEVELOPED BY M.F. HUTCHINSON AND
  809. C                F.R. DE HOOG, 'SMOOTHING NOISY DATA WITH SPLINE
  810. C                FUNCTIONS', NUMER. MATH. (IN PRESS)
  811. C
  812. C-----------------------------------------------------------------------
  813. C
  814.  
  815.  
  816. ALGORITHM
  817.  
  818. CUBGCV calculates a natural cubic spline curve which smoothes a given set
  819. of data points, using statistical considerations to determine the amount
  820. of smoothing required, as described in reference 2.  If the error variance
  821. is known, it should be supplied to the routine in VAR.  The degree of
  822. smoothing is then determined by minimizing an unbiased estimate of the true
  823. mean square error.  On the other hand, if the error variance is not known, VAR
  824. should be set to -1.0.  The routine then determines the degree of smoothing
  825. by minimizing the generalized cross validation.  This is asymptotically the
  826. same as minimizing the true mean square error (see reference 1).  In this
  827. case, an estimate of the error variance is returned in VAR which may be
  828. compared with any a priori approximate estimates.  In either case, an
  829. estimate of the true mean square error is returned in WK(5).  This estimate,
  830. however, depends on the error variance estimate, and should only be accepted
  831. if the error variance estimate is reckoned to be correct.
  832.  
  833. If job is set to 1, bayesian estimates of the standard error of each
  834. smoothed data value are returned in the array SE.  These also depend on
  835. the error variance estimate and should only be accepted if the error
  836. variance estimate is reckoned to be correct.  See reference 4.
  837.  
  838. The number of arithmetic operations and the amount of storage required by
  839. the routine are both proportional to N, so that very large data sets may be
  840. analysed.  The data points do not have to be equally spaced or uniformly
  841. weighted.  The residual and the spline coefficients are calculated in the
  842. manner described in reference 3, while the trace and various statistics,
  843. including the generalized cross validation, are calculated in the manner
  844. described in reference 2.
  845.  
  846. When VAR is known, any value of N greater than 2 is acceptable.  It is
  847. advisable, however, for N to be greater than about 20 when VAR is unknown.
  848. If the degree of smoothing done by CUBGCV when VAR is unknown is not
  849. satisfactory, the user should try specifying the degree of smoothing by
  850. setting VAR to a reasonable value.
  851.  
  852. References:
  853.  
  854. 1.  Craven, Peter and Wahba, Grace, "Smoothing noisy data with spline
  855.     functions", Numer. Math. 31, 377-403 (1979).
  856. 2.  Hutchinson, M.F. and de Hoog, F.R., "Smoothing noisy data with spline
  857.     functions", Numer. Math. (in press).
  858. 3.  Reinsch, C.H., "Smoothing by spline functions", Numer. Math. 10,
  859.     177-183 (1967).
  860. 4.  Wahba, Grace, "Bayesian 'confidence intervals' for the cross-validated
  861.     smoothing spline", J.R.Statist. Soc. B 45, 133-150 (1983).
  862.  
  863.  
  864. Example
  865.  
  866. A sequence of 50 data points are generated by adding a random variable with
  867. uniform density in the interval [-0.3,0.3] to the curve y=sin(3*pi*x/2).
  868. The abscissae are unequally spaced in [0,1].  Point standard error estimates
  869. are returned in SE by setting JOB to 1.  The error variance estimate is
  870. returned in VAR.  It compares favourably with the true value of 0.03.
  871. The IMSL function GGUBFS is used to generate sample values of a uniform
  872. variable on [0,1].
  873.  
  874.  
  875. INPUT:
  876.  
  877.       INTEGER          N,IC,JOB,IER
  878.       DOUBLE PRECISION X(50),F(50),Y(50),DF(50),C(49,3),WK(400),
  879.      *                 VAR,SE(50)
  880.       DOUBLE PRECISION GGUBFS,DSEED,DN
  881.       DATA DSEED/1.2345D4/
  882. C
  883.       N=50
  884.       IC=49
  885.       JOB=1
  886.       VAR=-1.0D0
  887.       DN=N
  888. C
  889.       DO 10 I=1,N
  890.       X(I)=(I - 0.5)/DN + (2.0*GGUBFS(DSEED) - 1.0)/(3.0*DN)
  891.       F(I)=DSIN(4.71238*X(I)) + (2.0*GGUBFS(DSEED) - 1.0)*0.3
  892.       DF(I)=1.0D0
  893.   10  CONTINUE
  894.       CALL CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
  895.        .
  896.        .
  897.        .
  898.       END
  899.  
  900. OUTPUT:
  901.  
  902. IER = 0
  903. VAR = 0.0279
  904. GENERALIZED CROSS VALIDATION = 0.0318
  905. MEAN SQUARE RESIDUAL = 0.0246
  906. RESIDUAL DEGREES OF FREEDOM = 43.97
  907. FOR CHECKING PURPOSES THE FOLLOWING OUTPUT IS GIVEN:
  908.  
  909. X(1)  = 0.0046    F(1)  =  0.2222     Y(1)  =  0.0342     SE(1)  = 0.1004
  910. X(21) = 0.4126    F(21) =  0.9810     Y(21) =  0.9069     SE(21) = 0.0525
  911. X(41) = 0.8087    F(41) = -0.5392     Y(41) = -0.6242     SE(41) = 0.0541
  912.  
  913.