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 w = WEIGHT, σ = SMPAR is the smoothing parameter, and N = NDATA.

Recommended values for σ depend on the weights w. 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

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

Visual Numerics, Inc. Visual Numerics - Developers of IMSL and PV-WAVE http://www.vni.com/ PHONE: 713.784.3131 FAX:713.781.9260