cubSplineSmooth¶

Computes a smooth cubic spline approximation to noisy data by using cross-validation to estimate the smoothing parameter or by directly choosing the smoothing parameter.

Synopsis¶

cubSplineSmooth (xdata, fdata)

Required Arguments¶

float xdata[] (Input): Array with ndata components containing the abscissas of the problem.
float fdata[] (Input): Array with ndata components containing the ordinates of the problem.

Return Value¶

The structure that represents the cubic spline. If a smoothed cubic spline cannot be computed, then None is returned.

Optional Arguments¶

weights, float[] (Input)

This option requires the user to provide the weights.

Default: all weights are equal to 1.

smoothingPar, float (Input)

This option sets the smoothing parameter σ = sigma explicitly.

Description¶

The function cubSplineSmooth is designed to produce a \(C^2\) cubic spline approximation to a data set in which the function values are noisy. This spline is called a smoothing spline.

Consider first the situation when the optional argument smoothingPar is selected. Then, a natural cubic spline with knots at all the data abscissas x = xdata is computed, but it does not interpolate the data \((x_i, f_i)\). The smoothing spline s is the unique \(C^2\) function which minimizes

\[\int_a^b s''(x)^2dx\]

subject to the constraint

\[\sum_{i=0}^{n-1} \left|\left(s\left(x_i\right) - f_i\right)w_i\right|^2 \leq \sigma\]

where w = weights, σ = sigma 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 \(f_i\) is available, then \(w_i\) should be set to the inverse of this value; and the smoothing parameter σ should be chosen in the confidence interval corresponding to the left side of the above inequality. That is,

\[n - \sqrt{2n} \leq \sigma \leq n + \sqrt{2n}\]

The function cubSplineSmooth is based on an algorithm of Reinsch (1967). This algorithm is also discussed in de Boor (1978, pp. 235−243).

The default for this function chooses the smoothing parameter σ by a statistical technique called cross-validation. For more information on this topic, refer to Craven and Wahba (1979).

The return value for this function is the structure Imsl_d_ppoly. This structure contains all the information to determine the spline (stored as a piecewise polynomial) that is computed by this procedure. For example, the following code sequence evaluates this spline at x and returns the value in y.

y = cubSplineValue (x, pp)

Examples¶

Example 1¶

In this example, function values are contaminated by adding a small “random” amount to the correct values. The function cubSplineSmooth is used to approximate the original, uncontaminated data.

from __future__ import print_function
from numpy import *
from pyimsl.math.cubSplineSmooth import cubSplineSmooth
from pyimsl.math.cubSplineValue import cubSplineValue
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform

# Define function


def F(x):
    return 1. + sin(x) + 7. * sin(3.0 * x)


# Generate random numbers
ndata = 90
randomSeedSet(123457)
random = randomUniform(ndata)

# Set up data
xdata = empty(ndata)
fdata = empty(ndata)
for i in range(0, ndata):
    xdata[i] = 6. * i / (ndata - 1)
    fdata[i] = F(xdata[i]) + .5 * (random[i] - .5)
pp = cubSplineSmooth(xdata, fdata)
print("        x        Error  ")
for i in range(0, 10):
    x = 6. * i / 9.
    error = F(x) - cubSplineValue(x, pp)
    print("%10.3f  %10.3f" % (x, error))

Output¶

        x        Error  
000      -0.201
667       0.070
333      -0.008
000      -0.058
667      -0.025
333       0.076
000      -0.002
667      -0.008
333       0.046
000       0.275

Example 2¶

Recall that in the first example, function values are contaminated by adding a small “random” amount to the correct values. Then, cubSplineSmooth is used to approximate the original, uncontaminated data. This example explicitly inputs the value of the smoothing parameter to be 5.

from __future__ import print_function
from numpy import *
from pyimsl.math.cubSplineSmooth import cubSplineSmooth
from pyimsl.math.cubSplineValue import cubSplineValue
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform

# Define function


def F(x):
    return 1. + sin(x) + 7. * sin(3.0 * x)


# Generate random numbers
ndata = 90
randomSeedSet(123457)
random = randomUniform(ndata)

# Set up data
xdata = empty(ndata)
fdata = empty(ndata)
for i in range(0, ndata):
    xdata[i] = 6. * i / (ndata - 1)
    fdata[i] = F(xdata[i]) + .5 * (random[i] - .5)
pp = cubSplineSmooth(xdata, fdata,
                     smoothingPar=5.0)
print("      x          Error  ")
for i in range(0, 10):
    x = 6. * i / 9.
    error = F(x) - cubSplineValue(x, pp)
    print("%10.3f  %10.3f" % (x, error))

Output¶

      x          Error  
000      -0.593
667       0.230
333      -0.116
000      -0.106
667       0.176
333      -0.071
000      -0.171
667       0.196
333      -0.036
000       0.971

Warning Errors¶

IMSL_MAX_ITERATIONS_REACHED The maximum number of iterations has been reached. The best approximation is returned.

Fatal Errors¶

`IMSL_DUPLICATE_XDATA_VALUES`	The xdata values must be distinct.
`IMSL_NEGATIVE_WEIGHTS`	All weights must be greater than or equal to zero.