CSSMH
Computes a smooth cubic spline approximation to noisy data.
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)
SMPAR — A nonnegative number which controls the smoothing. (Input)
The spline function S returned is such that
the sum from I = 1 to NDATA of ((S(XDATA(I)) ‑ FDATA(I)) / WEIGHT(I))**2
is less than or equal to SMPAR. It is recommended that SMPAR lie in the confidence interval of this sum, i.e., NDATA − SQRT(2 * NDATA).LE. SMPAR.LE. NDATA + SQRT(2 * NDATA).
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 2.
Default: NDATA = size (XDATA,1).
WEIGHT — Array of length NDATA containing estimates of the standard deviations of FDATA. (Input)
All elements of WEIGHT must be positive.
Default: WEIGHT = 1.0.
FORTRAN 90 Interface
Generic: CALL CSSMH (XDATA, FDATA, SMPAR, BREAK, CSCOEF [, …])
Specific: The specific interface names are S_CSSMH and D_CSSMH.
FORTRAN 77 Interface
Single: CALL CSSMH (NDATA, XDATA, FDATA, WEIGHT, SMPAR, BREAK, CSCOEF)
Double: The double precision name is DCSSMH.
Description
The routine CSSMH 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 S is the unique C2 function which minimizes
subject to the constraint
where wi = WEIGHT(I), σ = SMPAR is the smoothing parameter, and N = NDATA.
Recommended values for σ depend on the weights wi. If an estimate for the standard deviation of the error in the value fi is available, then wi should be set to this value and the smoothing parameter σ should be chosen in the confidence interval corresponding to the left side of the above inequality. That is,
The routine CSSMH is based on an algorithm of Reinsch (1967). This algorithm is also discussed in de Boor (1978, pages 235− 243).
Comments
1. Workspace may be explicitly provided, if desired, by use of C2SMH/DC2SMH. The reference is:
CALL C2SMH (NDATA, XDATA, FDATA, WEIGHT, SMPAR, BREAK, CSCOEF, WK, IWK)
The additional arguments are as follows:
WK — Work array of length 8 * NDATA + 5.
IWK — Work array of length NDATA.
2. Informational errors
Type | Code | Description |
|---|
3 | 1 | The maximum number of iterations has been reached. The best approximation is returned. |
4 | 3 | All weights must be greater than zero. |
3. The cubic spline can be evaluated using
CSVAL; its derivative can be evaluated using
CSDER.
Example
In this example, function values are contaminated by adding a small “random” amount to the correct values. The routine CSSMH is used to approximate the original, uncontaminated data.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER NDATA
PARAMETER (NDATA=300)
!
INTEGER I, NOUT
REAL BREAK(NDATA), CSCOEF(4,NDATA), ERROR, F,&
FDATA(NDATA), FLOAT, FVAL, SDEV, SMPAR, SQRT,&
SVAL, WEIGHT(NDATA), X, XDATA(NDATA), XT, RN
INTRINSIC FLOAT, SQRT
!
F(X) = 1.0/(.1+(3.0*(X-1.0))**4)
! 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
! Set the random number seed
CALL RNSET (1234579)
! Contaminate the data
DO 20 I=1, NDATA
RN = RNUNF()
FDATA(I) = FDATA(I) + 2.0*RN - 1.0
20 CONTINUE
! Set the WEIGHT vector
SDEV = 1.0/SQRT(3.0)
CALL SSET (NDATA, SDEV, WEIGHT, 1)
SMPAR = NDATA
! Smooth the data
CALL CSSMH (XDATA, FDATA, SMPAR, BREAK, CSCOEF, WEIGHT=WEIGHT)
! Get output unit number
CALL UMACH (2, NOUT)
! Write heading
WRITE (NOUT,99999)
! Print 10 values of the function.
DO 30 I=1, 10
XT = 90.0*(FLOAT(I-1)/FLOAT(NDATA-1))
! Evaluate the spline
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.1118 0.0995
0.3010 0.0514 0.0646 0.0131
0.6020 0.4690 0.2972 -0.1718
0.9030 9.3312 8.7022 -0.6289
1.2040 4.1611 4.7887 0.6276
1.5050 0.1863 0.2718 0.0856
1.8060 0.0292 0.1408 0.1116
2.1070 0.0082 0.0826 0.0743
2.4080 0.0031 0.0076 0.0045
2.7090 0.0014 -0.1789 -0.1803
