Lindo Systems
MODEL:
! LINGO model to fit a spline interpolation curve to a set of points.
Roughly speaking, find the smoothest curve that goes through
all the given (x, y) points. In any interval i, between two adjacent
given points, a cubic polynomial is derived to
interpolate between the two points.
Inputs:
Data:
1) NOBS = number of observations,
2) X(i) for i = 1, 2, ..., N,for each observation,
The X values of the set of points.
These should be in increasing order.
3) Y(i) for i = 1, 2, ..., N, .
The Y values of the set of points.
Outputs:
5) NERR = error code, = 0 if all OK, >0 if error, e.g.,
points are not in increasing order.
6) A(i) for i = 1, 2, ..., N-1,
7) B(i) for i = 1, 2, ..., N-1,
8) C(i) for i = 1, 2, ..., N-1,
9) D(i) for i = 1, 2, ..., N-1,
Details:
Given a specific XS, X(i)<= XS <= X(i+1),
the predicted(or corrected or adjusted) estimate of Y is given by:
delta = XS - X(i1),
Predict = D(i1)*delta^3 + C(i1)*delta^2 + B(i1)*delta + A(i1),
Roughly speaking, a spline, Y = F(X), is the "smoothest possible" curve that
a) passes through each of a given set of points. I.e., If X is
the measured value at one of the given points, then
F(X) = the true value at that point.
b) is continuous everywhere, i.e., F(X) has no jumps.
c) is smooth, i.e., the first derivative is continuous,-
it has no abrupt changes, F(x) has no sharp changes.
d) the second derivative is continuous everywhere, e.g.,
F(x) is not "curvy" in one interval and then abruptly
straight in the next interval ;
! Keywords: Spline fit, Interpolation;
SETS:
OBS: Y, X, A, B, C, D, DELTA, ABS2;
ENDSETS
DATA:
! Here are the calibration data points where we have
observations on Y = F(X);
! X = 1 2 3 4 5;
! Y = 1 3 2 5 4;
! Case: Twin Peaks;
X = 1 2 3 4 5;
Y = 1 3 2 3 1;
! Case: Quadratic;
! X = 1 2 3 4 5;
! Y = 1 4 9 16 25;
! Case: Step function;
! X = 2 3 4 5 6 7;
! Y = 1 1 1 2 2 2;
! Case: sin(x);
! X = 1 2 3 4 5 6 7 8 9;
! Y = 0.84147 0.90923 0.14112 -0.75680 -0.95892 -0.27942 0.65699 0.98936 0.41212;
! Case: Industrial measurement calibration data 1;
! X = 6.5 16.16378 36.20551 54.07717 71.98346;
! Y = 6.5 16.128 36.154 54.022 72.012 ;
! Case: Industrial measurement calibration data 2;
! X = 6.5 13 21.5 32.5 65;
! Y = 6.5 13 33.5 47.5 58.5;
! Internal parameters;
! Set TESTON > 0 if want code to generate output, else 0;
TESTON = 1;
! Wgt to be applied to minimizing sum of second derivatives,
in case there are alternate optima;
WGTSUM = 0.1;
ENDDATA
SUBMODEL FitSpline:
! This presumes we know the Y and want to predict X.
! Fit a cubic spline through the points.;
@FOR( OBS(i):
@FREE( A(i)); @FREE( B(i)); @FREE( C(i)); @FREE( D(i));
);
! There are 4*(NOBS-1) parameters to be estimated;
@FOR( OBS(i) | i #LT# NOBS:
! The curve in interval i must pass through X(i), Y(i);
[PASSTHRU] A(i) = Y(i); ! This fixes NOBS -1 of the paramaters;
! The curve must be continuous, i.e., in interval i must pass through X(i+1), Y(i+1),
This gives NOBS-1 more equations;
[CONT0] D(i)*DELTA(i)^3 + C(i)*DELTA(I)^2 + B(i)*DELTA(i) + A(i) = Y(i+1);
! Curve must be smooth, i.e., continuous first derivative at i+1,
(this gives NOBS-2 more equations, note B(NOBS) is free), so;
[CONT1] 3*D(i)*DELTA(i)^2 + 2*C(i)*DELTA(i) + B(i) = B(i+1);
! To avoid abrupt changes in curviness, we also want,
continuous 2nd derivative at i+1,
(this gives NOBS-2 more equations, note C(NOBS) is free) so ;
[CONT2] 6*D(i)*DELTA(i) + 2*C(i) = 2*C(i+1);
);
! One possible boundary restriction is that second derivatives be 0 at the end points...
( giving 4*(NOBS-1) linear equations in 4*(NOBS-1) variables;
C(1) = 0;
2*C( NOBS-1) + 6*D( NOBS-1)* DELTA( NOBS-1) = 0;
! Look at maximumum absolute second derivative;
@FOR( OBS(i) | i #LT# NOBS:
! The maximum in an interval occurs at one of the end points;
ABS2(i) >= 2*C(i);
ABS2(i) >= -2*C(i);
ABS2(i) >= 2*C(i) + 6*D(i)*DELTA(i);
ABS2(i) >= -(2*C(i) + 6*D(i)*DELTA(i));
MXABS2 >= ABS2(i);
);
! Minimize the maximum second derivative;
MIN = MXABS2;
ENDSUBMODEL
PROCEDURE EvalCurves:
! Input:
XS = a specific X value,
Output:
YL = corresponding Y value, using
linear interpolation,
YC = corresponding Y value, using
cubic spline interpolation;
! Find the interval containing XS;
i1 = 1;
@FOR( OBS(i) | i #GT# 1 #AND# X(i) #LT# XS:
i1 = i1 + 1;
);
! If falls above highest point, we extraplate using
the highest interval;
@IFC( i1 #GE# NOBS: i1 = NOBS-1);
! Similarly, if falls below lowest point, we extrapolate
using the lowest interval.
! We use the interval X(i1), X(i1+1);
! Get distance above low end of interval;
XDEL = XS - X(i1);
YL = (Y(i1)*(X(i1+1)- XS) + Y(i1+1)*(XS - X(i1)))/(X(i1+1)- X(i1));
YC = D(i1)*XDEL^3 + C(i1)*XDEL^2 + B(i1)*XDEL + A(i1);
ENDPROCEDURE
CALC:
@SET('TERSEO',3); ! Turn off all output;
! Check input for errors;
NERR = 0;
NOBS = @SIZE( OBS);
@FOR( OBS(i) | i #GT# 1:
! @IFC( Y(i) #LE# Y(I-1): NERR = 1);
@IFC( X(i) #LE# X(I-1): NERR = 1;); ! We expect the X(i)'s ordered;
DELTA( i-1) = X(i) - X(i-1); ! Compute length of each interval;
);
! @GEN( FitSpline);
! Fit a spline;
@SOLVE( FitSpline);
@IFC( TESTON:
@WRITE(' A simple report of the prediction coefficients', @NEWLINE(1));
@WRITE(' Interval D(i) C(i) B(i) A(i)', @NEWLINE(1));
@FOR( OBS(i) | i #LT# NOBS:
@WRITE( @FORMAT(i,'4.0F'),' ', @FORMAT( D(i),'10.6F'),' ', @FORMAT(C(i),'10.6F'),' ',
@FORMAT(B(i),'10.6F'),' ',@FORMAT(A(i),'10.5F'), @NEWLINE(1));
);
@WRITE(@NEWLINE(1));
);
XLO = @MIN( OBS(i): X(i));
XHI = @MAX( OBS(i): X(I));
! Plot some curves;
@CHARTPCURVE( 'Fitting a Line Through Points', 'x-axis', 'y-axis', EvalCurves, XS, XLO, XHI,
'Piecewise Linear', YL,
'Cubic Spline', YC );
ENDCALC
END