Computes a smooth cubic spline approximation to noisy data using cross-validation to estimate the smoothing parameter.
Required Arguments
XDATA Array of length NDATA containing the data point abscissas. (Input) XDATA must be distinct.
FDATA Array of length NDATA containing the data point ordinates. (Input)
IEQUAL A flag alerting the subroutine that the data is equally spaced. (Input)
BREAK Array of length NDATA containing the breakpoints for the piecewise cubic representation. (Output)
CSCOEF Matrix of size 4 by NDATA containing the local coefficients of the cubic pieces. (Output)
Optional Arguments
NDATA Number of
data points. (Input)
NDATA must be at least
3.
Default: NDATA = size (XDATA,1).
FORTRAN 90 Interface
Generic: CALL CSSCV (XDATA, FDATA, IEQUAL, BREAK, CSCOEF [, ])
Specific: The specific interface names are S_CSSCV and D_CSSCV.
FORTRAN 77 Interface
Single: CALL CSSCV (NDATA, XDATA, FDATA, IEQUAL, BREAK, CSCOEF)
Double: The double precision name is DCSSCV.
Description
The routine CSSCV is designed to produce a C2 cubic spline approximation to a data set in which the function values are noisy. This spline is called a smoothing spline. It is a natural cubic spline with knots at all the data abscissas x = XDATA, but it does not interpolate the data (xi, fi). The smoothing spline Ss is the unique C2 function that minimizes
subject to the constraint
where σ is the smoothing parameter and N = NDATA. The reader should consult Reinsch (1967) for more information concerning smoothing splines. The IMSL subroutine CSSMH solves the above problem when the user provides the smoothing parameter σ. This routine attempts to find the optimal' smoothing parameter using the statistical technique known as cross-validation. This means that (in a very rough sense) one chooses the value of σ so that the smoothing spline (Ss) best approximates the value of the data at xi, if it is computed using all the data except the i-th; this is true for all i = 1, , N. For more information on this topic, we refer the reader to Craven and Wahba (1979).
Comments
1. Workspace may be explicitly provided, if desired, by use of C2SCV/DC2SCV. The reference is:
CALL C2SCV (NDATA, XDATA, FDATA, IEQUAL, BREAK, CSCOEF, WK, SDWK, IPVT)
The additional arguments are as follows:
WK Work array of length 7 * (NDATA + 2).
SDWK Work array of length 2 * NDATA.
IPVT Work array of length NDATA.
2. Informational error
Type Code
4 2 Points in the data point abscissas array, XDATA, must be distinct.
Example
In this example, function values are computed and are contaminated by adding a small random amount. The routine CSSCV is used to try to reproduce the original, uncontaminated data.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER NDATA
PARAMETER (NDATA=300)
!
INTEGER I, IEQUAL, NOUT
REAL BREAK(NDATA), CSCOEF(4,NDATA), ERROR, F,&
FDATA(NDATA), FLOAT, FVAL, SVAL, X,&
XDATA(NDATA), XT, RN
INTRINSIC FLOAT
!
F(X) = 1.0/(.1+(3.0*(X-1.0))**4)
!
CALL UMACH (2, NOUT)
! Set up a grid
DO 10 I=1, NDATA
XDATA(I) = 3.0*(FLOAT(I-1)/FLOAT(NDATA-1))
FDATA(I) = F(XDATA(I))
10 CONTINUE
! Introduce noise on [-.5,.5]
! Contaminate the data
CALL RNSET (1234579)
DO 20 I=1, NDATA
RN = RNUNF ()
FDATA(I) = FDATA(I) + 2.0*RN - 1.0
20 CONTINUE
!
! Set IEQUAL=1 for equally spaced data
IEQUAL = 1
! Smooth data
CALL CSSCV (XDATA, FDATA, IEQUAL, BREAK, CSCOEF)
! Print results
WRITE (NOUT,99999)
DO 30 I=1, 10
XT = 90.0*(FLOAT(I-1)/FLOAT(NDATA-1))
SVAL = CSVAL(XT,BREAK,CSCOEF)
FVAL = F(XT)
ERROR = SVAL - FVAL
WRITE (NOUT,'(4F15.4)') XT, FVAL, SVAL, ERROR
30 CONTINUE
99999 FORMAT (12X, 'X', 9X, 'Function', 7X, 'Smoothed', 10X,&
'Error')
END
Output
X
Function
Smoothed
Error
0.0000
0.0123
0.2528
0.2405
0.3010
0.0514
0.1054
0.0540
0.6020
0.4690
0.3117
-0.1572
0.9030
9.3312
8.9461
-0.3850
1.2040
4.1611
4.6847
0.5235
1.5050
0.1863
0.3819
0.1956
1.8060
0.0292
0.1168
0.0877
2.1070
0.0082
0.0658
0.0575
2.4080
0.0031
0.0395
0.0364
2.7090
0.0014 -0.2155
-0.2169
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |