ACADEMIA

The E02 Chapter - Fitting a set of points with a cubic spline

A common task in many areas of science is to find a function which approximates a set of data points. The data may contain random perturbations - such as measurement errors - which it is desirable to smooth out. In our previous article, we looked at how to use the NAG Toolbox to compute a bicubic spline approximation to a two-dimensional set of data values given on a grid in the x-y plane; the result was displayed as a surface, f(x,y). Here, we investigate the one-dimensional case, where we seek to find a cubic spline curve f(x) which approximates a set of data values given on the x-axis. The NAG routine which does this is e02be, while the supplementary routine e02bb can be used to evaluate the spline at any point, given the results from e02be.

As in the two-dimensional example, a single parameter is specified to control the trade-off between closeness and smoothness of fit. A small value for this smoothing parameter, named s in e02be, will give a closer fit to the data, but the exact meaning of “small” in this context will vary according to the particular data set. In the extremes, a value of s = 0 will produce an interpolating spline, while a very large value will produce the weighted least-squares cubic polynomial approximation to the points.

% Here's the (x,y) data to be fitted, followed by the 
% weights for each point (set these to 1 if all points
% are equally important).  
x = [0.5; 1; 1.5; 2; 2.5; 3; 4; 4.5; 5; 5.5; 6; 7; 7.5];
y = [-0.372; 0.431; 1.69; 2.11; 3.1; 4.23; 4.35;
     4.81; 4.61; 4.79; 5.23; 6.35; 7.19];
w = [2; 1.5; 1; 3; 1; 0.5; 1; 2; 2.5; 1; 3; 1; 2];

% Calculate work array lengths, according to the e02be documentation.
m = length(x);
nest = m+4;
lwrk = 4*m + 16*nest + 41;

% Now initialize some parameters required by the NAG routine.
start = 'C';
n = int32(0);
lamda = zeros(nest,1);
wrk = zeros(lwrk, 1);
iwrk = zeros(nest, 1, 'int32');

% Set the smoothness parameter and call the NAG routine to 
% compute the cubic spline approximation to the data.
s = 10000.0;
[nOut, lamdaOut, c, fp, wrkOut, iwrkOut, ifail] = ...
    e02be(start, x, y, w, s, n, lamda, wrk, iwrk);

% Set up the points at which the spline is to be evaluated.
px = [x(1) : 0.01 : x(m)];
pf = zeros(size(px));

% Now call the NAG routine e02bb to calculate the spline, and plot it.
for i = 1 : length(px)
  [pf(i), ifail] = e02bb(lamdaOut, c, px(i), 'ncap7', nOut);
end
plot(px, pf, 'Color', [.8 0 0], 'Linewidth', 2);

Figure 3 shows the data points and the spline for three values of s. It can be seen that, as s is decreased, the fit becomes less smooth, and closer to the points.

Figure 3: Fitting a cubic spline through data points.

The MATLAB script for this demo is available as the file NAGToolboxDemos/Curve_and_surface_fitting/e02be_demo.m, distributed in this archive.