\documentstyle[openbib]{article}

\topmargin  -20pt

\footheight 0pt

\headheight 0pt

\headsep 0pt

\oddsidemargin 0pt

\evensidemargin 0pt

\parindent 0pt

\parskip 1.5ex plus 0.5ex minus 0.5ex

\textwidth 450pt

\textheight 57\baselineskip

\begin{document}



\section*{Procedures for fitting Cubic Smoothing Splines.}

\begin{centering}

{\bf David Baird, MAF Technology, PO Box 24, Lincoln, New Zealand.} \\

\end{centering}

\vspace{2ex}



In fitting a curve to observations of an independent variable, $y$, at design

points, $x$, one has two main purposes. Firstly exploring and representing the

relationship between the design and response variable; secondly allowing

predictions to be made for future observations of the response variable.

Typically a parametric curve is used to fit the relationship but this

forces the model for the relationship into a rigidly defined class of models.

Non-parametric models for the curve allow much more flexibility in the fitted

model. A cubic spline curve (Reinsch, 1967; Rice and Rosenblatt, 1983)

is a popular non-parametric curve which is very flexible but still has useful

smoothing properties. A good review of spline smoothing can be found in

Silverman (1985).



Cubic splines are piecewise cubic polynomials which can be fitted to a series

of $(x_{i},y_{i})$ data points. A full cubic spline with {\it knots} at $x_{i}$

(where the piecewise portions join) exactly interpolates the data points.

The piecewise portions are defined so that at the knots the function and its

first two derivatives are continuous.

A cubic spline with knots at $x_{i}$, $i=1,\ldots,n$ is defined

(letting $x_{0}=-\infty$, $x_{n+1}=\infty$) by:

\[ f(x) = a_{i}+b_{i}x+c_{i}x^{2}+d_{i}x^{3}

          \; \mbox{ if } x_{i} \leq x < x_{i+1} \]



subject to the following restrictions:

\[ \begin{array}{rcl}

   a_{i-1} + b_{i-1}x_{i} + c_{i-1}x_{i}^{2} + d_{i-1}x_{i}^{3} & = &

   a_{i} + b_{i}x_{i} + c_{i}x_{i}^{2} + d_{i}x_{i}^{3} \\

   b_{i-1} + 2c_{i-1}x_{i} + 3d_{i-1}x_{i}^{2} & = &

   b_{i} + 2c_{i}x_{i} + 3d_{i}x_{i}^{2} \\

   2c_{i-1} + 6d_{i-1}x_{i} & = &

   2c_{i} + 6d_{i}x_{i} \\

   c_{0} = d_{0} = c_{n} = d_{n} & = & 0

   \end{array} \]



The first three restrictions are to ensure the function and its first and

second derivatives respectively are continuous at the knots.

The final restriction above means that cubic spline is a linear function

outside the range of the knots.

However the cubic spline has a discontinuous third derivative:

\[ f^{\prime\prime\prime}(x) = d_{i} \; \mbox{ if } x_{i} \leq x < x_{i+1} \]



The cubic spline has the property that it is the interpolating function which

minimizes the integrated squared second derivative

($\int (f^{\prime\prime})^{2}$).

The integrated squared second derivative acts as a roughness penalty,

with smooth or slowly varying curves giving small values (for example a

straight line gives a value of 0).



However when dealing with noisy data it does not make sense to interpolate

the data points, since the resulting variance of any predicted values will

be very large. Rather it makes more sense to use a more slowly varying curve

which does not interpolate the data points but damps out the noise component

while minimizing the residual sum of squares $\sum_{i}\{y_{i}-f(x_{i})\}^{2}$.

By increasing the smoothness of the curve, the variance of predicted

values is reduced. However at the same time this can introduce bias in the

predicted values where the true underlying relationship is changing rapidly.

Fitting a cubic spline which minimizes the residual sums of squares

with a restriction on the smoothness of the resulting cubic spline can be

achieved by using a modified sum of squares:

\[S(p) = p\sum_{i}\{y_{i}-f(x_{i})\}^{2} + (1-p)\int f^{\prime\prime}(x)^{2}\]



The {\it smoothing parameter}, $p$, controls the trade off between

minimizing the residual error and minimizing local variation.

Equivalently, $p$ controls the trade off between reducing variance and

reducing possible bias in the predicted values.

A $p$ value of 1 gives the function which minimizes the residual error

(the interpolating cubic spline).

A $p$ value of 0 gives the function which minimizes the roughness penalty

(the linear regression of $y$ on $x$).

The curve which minimizes $S(p)$ is known as a {\it cubic smoothing spline}.

An added sophistication is to allow the contribution to $S(p)$ of the residual errors

at $x$ to be weighted using a vector of weights $w$. For example, differential

weights could be used to reflect differing estimated variances of the $y$

values, if this changes with $x$.



The following GAUSS procedure, CSSPLINE, fits the weighted cubic smoothing

spline to a set of points for a particular value of the smoothing parameter

$p$. The algorithm is based on one by de Boor (1978) but is not fully efficient as

it does not take advantage of the banded nature of the design matrices,

{\tt q} and {\tt r}. The procedure returns a matrix of coefficients for the

cubic smoothing spline with the columns of the matrix containing the knots,

the fitted function value and its first three derivatives evaluated at the

knots. The fitted curve can be evaluated at any set of $x$ values using the

procedure FSPLINE, which takes the spline coefficients returned by CSSPLINE

and the matrix of $x$ values and returns a corresponding matrix of fitted

$y$ values.



\pagebreak[1]

\begin{verbatim}

/* CSSPLINE - Fit Cubic Smoothing Spline

**

** Usage:   c = csspline(x,y,w,p) ;

**

** Input:   x - N x 1 vector of x values of the set of points;

**              x must be sorted into ascending order and must not contain

**              duplicate values.

**          y - N x 1 vector of y values of the set of points.

**          w - N x 1 vector of weights for the set of points;

**              for equal weighting set w to a scalar for efficiency.

**          p - Scalar giving the smoothing parameter. This defines the

**              relative weighting of the residual sum of squares and the

**              integrated squared second derivative (a roughness penalty)

**              of the curve; p = 0 gives an interpolating cubic spline and

**              p = 1 gives a linear regression fit (0 <= p <= 1).

**

** Output:  c - N x 5 matrix of coefficients of cubic smoothing spline;

**              c[.,1] = knots,

**              c[.,2..5] = function value and first 3 derivatives at knots.

**

** Reference: C. de Boor 1978, "A Practical Guide to Splines", Springer-Verlag.

*/

proc csspline(x,y,w,p) ;

     local n,n1,n2,dx,dxi,r,q,qw,t,u,c0,c1,c2,c3 ;

     n   = rows(x) ;

     w   = w./meanc(w) ;

     t   = 1 - p ;

     n1  = n - 1 ;

     n2  = n - 2 ;

     dx  = x[2:n] - x[1:n1] ;

     dxi = 1./dx ;

     @ Create R matrix which fits linear regression to y @

     r = rotater((0|dx[2:n2])~2.*(dx[1:n2]+dx[2:n1])~(dx[2:n2]|0)~

                 zeros(n2,n-5),seqa(-1,1,n2))./6 ;

     @ Create Q matrix which fits cubic spline to y @

     q = rotater(dxi[1:n2]~(-dxi[1:n2]-dxi[2:n1])~dxi[2:n1]~zeros(n2,n-3),

                 seqa(0,1,n2))' ;

     @ Solve for U where ((1-p)Q'W^(-1)Q + pR)U = Q'y @

     qw = q./sqrt(w) ;

     u   = solpd(q'y,t.*moment(qw,0) + p.*r) ;

     @ Find derivatives, c0-c3, of smoothing spline at x @

     c0  = y - ((t./w).*q)*u ;                                    @ f(x) @

     c2  = 0 | (p.*u) | 0 ;                                       @ f''(x) @

     c3  = (c2[2:n] - c2[1:n1])./dx | 0 ;                         @ f'''(x) @

     c1  = (c0[2:n] - c0[1:n1])./dx - dx.*(c2[1:n1]./2 + dx.*c3[1:n1]./6) ;

     c1  = c1 | (c1[n1] + dx[n1].*(c2[n1] + dx[n1].*c3[n1]./2)) ; @ f'(x) @

     retp(x~c0~c1~c2~c3) ;

endp ;

\end{verbatim}



\pagebreak[1]

\begin{verbatim}

/* FSPLINE - Calculates fitted values for a cubic spline.

**

** Usage:  y = fspline(c,x) ;

**

** Input:  c - N x 5 matrix of spline coefficients output by CSSPLINE.

**         x - N x 1 vector of x values for which the spline function is to

**             be evaluated at.

**

** Output: y - N x 1 vector of fitted values corresponding to the x values.

*/

proc fspline(c,x) ;

     local i,z ;

     @ i locates which interval each value of x is in @

     i = 1 + sumc(x' .> c[2:rows(c),1]) ;

     @ Calculate distance to x from start of interval that x falls in @

     z = x - c[i,1] ;

     @ Return cubic spline for x @

     retp(c[i,2] + z.*(c[i,3] + (z.>=0).*z.*(c[i,4]./2 + z.*c[i,5]./6))) ;

endp ;

\end{verbatim}



The smoothing parameter $p$ is a free variable in the above approach and

may be chosen by the user of this procedure as seen fit. Varying this

parameter between 0 and 1 will slowly allow more variation in the fitted

curve until it exactly interpolates all the data points. One approach to

choosing the value of $p$ to be used, is to manually change $p$ until the

resulting curve is as smooth as possible whilst still showing no systematic

bias in following the data points. Other procedures for automatically

choosing an optimum value of $p$ have been defined. The most popular class of

these methods is based on {\it cross validation} (CV). The basic principle of

cross validation is to leave out the data points one at a time and to choose

the value of $p$ for which missing data points are best estimated by the

cubic smoothing spline using the remainder of the data points.

Let $f_{p}^{-i}$ be the cubic smoothing spline calculated from all data pairs

except $(x_{i},y_{i})$, using the value of $p$ for the smoothing parameter.

The cross validation score XVSC$(p)$ is defined by:

\[ \mbox{XVSC}(p) = n^{-1}\sum_{i}\{y_{i} - f_{p}^{-i}(x_{i})\}^{2} \]



This is identical to the PRESS criterion for regression model selection.

The value of $p$ which minimizes XVSC$(p)$ is then the cross validation choice

for the smoothing parameter.

An alternative to this criterion is to use a related criterion called,

{\it generalized cross validation} (GCV), which reweights the contribution of

deletion residuals to the total score (Craven and Wahba, 1979). The smoothing

parameters chosen by CV and GCV are in general very similar.

An approximation of the GCV criterion (AGCV) has been developed by Silverman

(1984), which is claimed to have superior speed and properties to GCV.



Procedures written by the author which find the CV, GCV and AGCV choices of

the smoothing parameter, along with several other procedures dealing with

cubic smoothing splines have been put together in the GAUSS procedure file

SPLINE.ARC.



The following example and the accompanying graphs illustrate the use of the

spline procedures.



\pagebreak[2]

\begin{verbatim}

#include spline.arc ;                  /* Include spline procedures */

call splset ;                          /* Reset spline globals */

rndseed 123357 ;                       /* Initial random number generator */

n  = 25 ;                              /* Number of data points */

x  = pi.*seqa(0.05,.025,n) ;           /* x data points */

y  = sin(x)^2 + rndn(n,1)./10 ;        /* y data points */

w  = 1 ;                               /* Use equal weighting for points */

xp = pi.*seqa(0,.005,5*n+30) ;         /* Points to plot fitted values at */

y1 = fspline(csspline(x,y,w,1),xp) ;   /* Interpolating cubic spline (p=1)*/

y7 = fspline(csspline(x,y,w,.75),xp) ; /* Cubic smoothing spline (p=.75) */

y0 = fspline(csspline(x,y,w,0),xp) ;   /* Linear Regression (p=0) */

yc = fspline(cvspline(x,y,w),xp) ;     /* Cross validation smoothing spline */



graphset ;                             /* Default graphics options */

_ptek = "splines.tkf" ;                /* Name of tek file (Fig 1) */

title("Cubic smoothing splines") ;     /* Graph title */

xlabel("x") ;                          /* Label x axis */

ylabel("y") ;                          /* Label y axis */

_psym = x~y~ones(n,1).*(8~0~14~1~0) ;  /* Plot data points */

let _plegctl = 1 5 ;                   /* Legend for various curves */

_plegstr = "p = 1\000p = .75" $+       /* Legend labels */

           "\000p = 0\000p = min XVSC" ;

let _plwidth = 3 3 3 3 ;               /* Set width of lines */

let _pltype = 6 5 3 1 ;                /* Set up line patterns */

let _pcolor = 9 10 11 12 ;             /* Set up line colors */

xy(xp,y1~y7~y0~yc) ;                   /* Plot data and 4 splines */



w  = 1./seqa(1,.1,25) ;                /* Use unequal weighting for points */

p  = min_gcv(x,y,w) ;                  /* Smoothing parameter for min. GCV */

yg = fspline(csspline(x,y,w,p),x) ;    /* Calculate fitted values */

l  = spl_clim(x,y,w,p,.95) ;           /* Calculate 95% confidence limits */



_ptek = "gspline.tkf" ;                /* Name of graphics file (Fig 2) */

title("GCV spline and 95% confidence limits") ; /* Graph title */

yn = sin(x)^2 ;                        /* Original curve with no error */

_plegstr = "Sin[2] (no error)" $+      /* Legend labels */

           "\000Fitted curve" $+

           "\000Lower CL\000Upper CL" ;

let _plwidth = 4 4 2 2 ;               /* Set width of lines */

let _pltype =  5 6 3 3 ;               /* Set line patterns */

xy(x,yn~yg~l) ;                        /* Plot spline and confidence limits */

end ;

\end{verbatim}



\pagebreak[1]

\bibliographystyle{plain}

\begin{thebibliography}{}



\bibitem{}

Craven, P. and Wahba, G. (1979). Smoothing noisy data with spline functions.

{\it Numerische Mathematik}, {\bf 31}, 377---403.



\bibitem{}

de Boor, C. (1978). {\it A Practical Guide to Splines.}

Springer--Verlag, New York.



\bibitem{}

Reinsch, C. (1967). Smoothing by spline functions.

{\it Numerische Mathematik}, {\bf 10}, 177---183.



\bibitem{}

Rice, J. and Rosenblatt, M. (1983). Smoothing splines: regression, derivatives

and deconvolution. {\it Annals of Statistics}, {\bf 11}, 141---156.



\bibitem{}

Silverman B.W. (1984). A fast and efficient cross-validation method for

smoothing parameter choice in spline regression.

{\it Journal of the American Statistical Association}, {\bf 79}, 584---589.



\bibitem{}

Silverman B.W. (1985). Some aspects of the spline smoothing approach to

non-parametric regression curve fitting.

{\it Journal of the Royal Statistical Society} B, {\bf 47}, 1---52.



\end{thebibliography}

\end{document}