cub_spline_smooth
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
#include <imsl.h>
Imsl_f_ppoly *imsl_f_cub_spline_smooth (int ndata, float xdata[], float fdata[], …, 0)
The type Imsl_d_ppoly function is imsl_d_cub_spline_smooth.
Required Arguments
int ndata (Input)
Number of data points.
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
A pointer to the structure that represents the cubic spline. If a smoothed cubic spline cannot be computed, then NULL is returned. To release this space, use imsl_free.
Synopsis with Optional Arguments
#include <imsl.h>s
Imsl_f_ppoly *imsl_f_cub_spline_smooth (int ndata, float xdata[], float fdata[],
IMSL_WEIGHTS, float weights[],
IMSL_SMOOTHING_PAR, float sigma,
0)
Optional Arguments
IMSL_WEIGHTS, float weights[] (Input)
This option requires the user to provide the weights.
Default: all weights are equal to 1.
IMSL_SMOOTHING_PAR, float sigma (Input)
This option sets the smoothing parameter σ = sigma explicitly.
Description
The function imsl_f_cub_spline_smooth 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.
Consider first the situation when the optional argument IMSL_SMOOTHING_PAR 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 (xi, fi). The smoothing spline s is the unique C2 function which minimizes
subject to the constraint
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 fi is available, then wi 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,
The function imsl_f_cub_spline_smooth 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 a pointer to the structure Imsl_f_ppoly. The calling program must receive this in a pointer Imsl_f_ppoly *pp. 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 = imsl_f_cub_spline_value (x, pp, 0);
Examples
Example 1
In this example, function values are contaminated by adding a small “random” amount to the correct values. The function imsl_f_cub_spline_smooth is used to approximate the original, uncontaminated data.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define NDATA 90
/* Define function */
#define F(x) (float)(1.+ sin(x)+7.*sin(3.0*x))
int main()
{
int i;
float fdata[NDATA], xdata[NDATA], *random;
Imsl_f_ppoly *pp;
/* Generate random numbers */
imsl_random_seed_set(123457);
random = imsl_f_random_uniform(NDATA, 0);
/* Set up data */
for (i = 0; i < NDATA; i++) {
xdata[i] = 6.*(float)i /((float)(NDATA-1));
fdata[i] = F(xdata[i]) + .5*(random[i]-.5);
}
pp = imsl_f_cub_spline_smooth(NDATA, xdata, fdata, 0);
printf(" x error \n");
for(i = 0; i < 10; i++){
float x, error;
x = 6.*i/9.;
error = F(x) - imsl_f_cub_spline_value(x, pp, 0);
printf("%10.3f %10.3f\n", x, error);
}
}
Output
x Error
0.000 -0.201
0.667 0.070
1.333 -0.008
2.000 -0.058
2.667 -0.025
3.333 0.076
4.000 -0.002
4.667 -0.008
5.333 0.045
6.000 0.276
Example 2
Recall that in the first example, function values are contaminated by adding a small “random” amount to the correct values. Then, imsl_f_cub_spline_smooth is used to approximate the original, uncontaminated data. This example explicitly inputs the value of the smoothing parameter to be 5.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define NDATA 90
/* Define function */
#define F(x) (float)(1.+ sin(x)+7.*sin(3.0*x))
int main()
{
int i;
float fdata[NDATA], xdata[NDATA], *random;
Imsl_f_ppoly *pp;
/* Generate random numbers */
imsl_random_seed_set(123457);
random = imsl_f_random_uniform(NDATA, 0);
/* Set up data */
for (i = 0; i < NDATA; i++) {
xdata[i] = 6.*(float)i /((float)(NDATA-1));
fdata[i] = F(xdata[i]) + .5*(random[i]-.5);
}
pp = imsl_f_cub_spline_smooth(NDATA, xdata, fdata,
IMSL_SMOOTHING_PAR, 5.0,
0);
printf(" x error \n");
for(i = 0; i < 10; i++){
float x, error;
x = 6.*i/9.;
error = F(x) - imsl_f_cub_spline_value(x, pp, 0);
printf("%10.3f %10.3f\n", x, error);
}
}
Output
x Error
0.000 -0.593
0.667 0.230
1.333 -0.116
2.000 -0.106
2.667 0.176
3.333 -0.071
4.000 -0.171
4.667 0.196
5.333 -0.036
6.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. |
