NAG Library Routine Document
g10abf
(fit_spline)
1
Purpose
g10abf fits a cubic smoothing spline for a given smoothing parameter.
2
Specification
Fortran Interface
Subroutine g10abf ( 
mode,
weight,
n,
x,
y,
wt,
rho,
yhat,
c,
ldc,
rss,
df,
res,
h,
comm,
ifail) 
Integer, Intent (In)  :: 
n,
ldc  Integer, Intent (Inout)  :: 
ifail  Real (Kind=nag_wp), Intent (In)  :: 
x(n),
y(n),
wt(*),
rho  Real (Kind=nag_wp), Intent (Inout)  :: 
c(ldc,3),
comm(9*n+14)  Real (Kind=nag_wp), Intent (Out)  :: 
yhat(n),
rss,
df,
res(n),
h(n)  Character (1), Intent (In)  :: 
mode,
weight 

C Header Interface
#include nagmk26.h
void 
g10abf_ (
const char *mode,
const char *weight,
const Integer *n,
const double x[],
const double y[],
const double wt[],
const double *rho,
double yhat[],
double c[],
const Integer *ldc,
double *rss,
double *df,
double res[],
double h[],
double comm[],
Integer *ifail,
const Charlen length_mode,
const Charlen length_weight) 

3
Description
g10abf fits a cubic smoothing spline to a set of $n$ observations (${x}_{i}$, ${y}_{i}$), for $i=1,2,\dots ,n$. The spline provides a flexible smooth function for situations in which a simple polynomial or nonlinear regression model is unsuitable.
Cubic smoothing splines arise as the unique realvalued solution function
$f$, with absolutely continuous first derivative and squaredintegrable second derivative, which minimizes:
where
${w}_{i}$ is the (optional) weight for the
$i$th observation and
$\rho $ is the smoothing parameter. This criterion consists of two parts: the first measures the fit of the curve, and the second the smoothness of the curve. The value of the smoothing parameter
$\rho $ weights these two aspects; larger values of
$\rho $ give a smoother fitted curve but, in general, a poorer fit. For details of how the cubic spline can be estimated see
Hutchinson and de Hoog (1985) and
Reinsch (1967).
The fitted values,
$\hat{y}={\left({\hat{y}}_{1},{\hat{y}}_{2},\dots ,{\hat{y}}_{n}\right)}^{\mathrm{T}}$, and weighted residuals,
${r}_{i}$, can be written as
for a matrix
$H$. The residual degrees of freedom for the spline is
$\mathrm{trace}\left(IH\right)$ and the diagonal elements of
$H$,
${h}_{ii}$, are the leverages.
The parameter
$\rho $ can be chosen in a number of ways. The fit can be inspected for a number of different values of
$\rho $. Alternatively the degrees of freedom for the spline, which determines the value of
$\rho $, can be specified, or the (generalized) crossvalidation can be minimized to give
$\rho $; see
g10acf for further details.
g10abf requires the
${x}_{i}$ to be strictly increasing. If two or more observations have the same
${x}_{i}$value then they should be replaced by a single observation with
${y}_{i}$ equal to the (weighted) mean of the
$y$ values and weight,
${w}_{i}$, equal to the sum of the weights. This operation can be performed by
g10zaf.
The computation is split into three phases.
(i) 
Compute matrices needed to fit spline. 
(ii) 
Fit spline for a given value of $\rho $. 
(iii) 
Compute spline coefficients. 
When fitting the spline for several different values of
$\rho $, phase
(i) need only be carried out once and then phase
(ii) repeated for different values of
$\rho $. If the spline is being fitted as part of an iterative weighted least squares procedure phases
(i) and
(ii) have to be repeated for each set of weights. In either case, phase
(iii) will often only have to be performed after the final fit has been computed.
The algorithm is based on
Hutchinson (1986).
4
References
Hastie T J and Tibshirani R J (1990) Generalized Additive Models Chapman and Hall
Hutchinson M F (1986) Algorithm 642: A fast procedure for calculating minimum crossvalidation cubic smoothing splines ACM Trans. Math. Software 12 150–153
Hutchinson M F and de Hoog F R (1985) Smoothing noisy data with spline functions Numer. Math. 47 99–106
Reinsch C H (1967) Smoothing by spline functions Numer. Math. 10 177–183
5
Arguments
 1: $\mathbf{mode}$ – Character(1)Input

On entry: indicates in which mode the routine is to be used.
 ${\mathbf{mode}}=\text{'P'}$
 Initialization and fitting is performed. This partial fit can be used in an iterative weighted least squares context where the weights are changing at each call to g10abf or when the coefficients are not required.
 ${\mathbf{mode}}=\text{'Q'}$
 Fitting only is performed. Initialization must have been performed previously by a call to g10abf with ${\mathbf{mode}}=\text{'P'}$. This quick fit may be called repeatedly with different values of rho without reinitialization.
 ${\mathbf{mode}}=\text{'F'}$
 Initialization and full fitting is performed and the function coefficients are calculated.
Constraint:
${\mathbf{mode}}=\text{'P'}$, $\text{'Q'}$ or $\text{'F'}$.
 2: $\mathbf{weight}$ – Character(1)Input

On entry: indicates whether userdefined weights are to be used.
 ${\mathbf{weight}}=\text{'W'}$
 Userdefined weights should be supplied in wt.
 ${\mathbf{weight}}=\text{'U'}$
 The data is treated as unweighted.
Constraint:
${\mathbf{weight}}=\text{'W'}$ or $\text{'U'}$.
 3: $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of distinct observations.
Constraint:
${\mathbf{n}}\ge 3$.
 4: $\mathbf{x}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the distinct and ordered values
${x}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
Constraint:
${\mathbf{x}}\left(\mathit{i}\right)<{\mathbf{x}}\left(\mathit{i}+1\right)$, for $\mathit{i}=1,2,\dots ,n1$.
 5: $\mathbf{y}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the values
${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
 6: $\mathbf{wt}\left(*\right)$ – Real (Kind=nag_wp) arrayInput

Note: the dimension of the array
wt
must be at least
${\mathbf{n}}$ if
${\mathbf{weight}}=\text{'W'}$.
On entry: if
${\mathbf{weight}}=\text{'W'}$,
wt must contain the
$n$ weights. Otherwise
wt is not referenced and unit weights are assumed.
Constraint:
if ${\mathbf{weight}}=\text{'W'}$, ${\mathbf{wt}}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,n$.
 7: $\mathbf{rho}$ – Real (Kind=nag_wp)Input

On entry: $\rho $, the smoothing parameter.
Constraint:
${\mathbf{rho}}\ge 0.0$.
 8: $\mathbf{yhat}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the fitted values,
${\hat{y}}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
 9: $\mathbf{c}\left({\mathbf{ldc}},3\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: if
${\mathbf{mode}}=\text{'Q'}$,
c must be unaltered from the previous call to
g10abf with
${\mathbf{mode}}=\text{'P'}$. Otherwise
c need not be set.
On exit: if
${\mathbf{mode}}=\text{'F'}$,
c contains the spline coefficients. More precisely, the value of the spline at
$t$ is given by
$\left(\left({\mathbf{c}}\left(i,3\right)\times d+{\mathbf{c}}\left(i,2\right)\right)\times d+{\mathbf{c}}\left(i,1\right)\right)\times d+{\hat{y}}_{i}$, where
${x}_{i}\le t<{x}_{i+1}$ and
$d=t{x}_{i}$.
If
${\mathbf{mode}}=\text{'P'}$ or
$\text{'Q'}$,
c contains information that will be used in a subsequent call to
g10abf with
${\mathbf{mode}}=\text{'Q'}$.
 10: $\mathbf{ldc}$ – IntegerInput

On entry: the first dimension of the array
c as declared in the (sub)program from which
g10abf is called.
Constraint:
${\mathbf{ldc}}\ge {\mathbf{n}}1$.

On exit: the (weighted) residual sum of squares.
 12: $\mathbf{df}$ – Real (Kind=nag_wp)Output

On exit: the residual degrees of freedom.
 13: $\mathbf{res}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the (weighted) residuals,
${r}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
 14: $\mathbf{h}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the leverages,
${h}_{\mathit{i}\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
 15: $\mathbf{comm}\left(9\times {\mathbf{n}}+14\right)$ – Real (Kind=nag_wp) arrayCommunication Array

On entry: if
${\mathbf{mode}}=\text{'Q'}$,
comm must be unaltered from the previous call to
g10abf with
${\mathbf{mode}}=\text{'P'}$. Otherwise
comm need not be set.
On exit: if
${\mathbf{mode}}=\text{'P'}$ or
$\text{'Q'}$,
comm contains information that will be used in a subsequent call to
g10abf with
${\mathbf{mode}}=\text{'Q'}$.
 16: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

On entry,  ${\mathbf{n}}<3$, 
or  ${\mathbf{ldc}}<{\mathbf{n}}1$, 
or  ${\mathbf{rho}}<0.0$, 
or  ${\mathbf{mode}}\ne \text{'Q'}$, $\text{'P'}$ or $\text{'F'}$, 
or  ${\mathbf{weight}}\ne \text{'W'}$ or $\text{'U'}$. 
 ${\mathbf{ifail}}=2$

On entry,  ${\mathbf{weight}}=\text{'W'}$ and at least one element of ${\mathbf{wt}}\le 0.0$. 
 ${\mathbf{ifail}}=3$

On entry,  ${\mathbf{x}}\left(i\right)\ge {\mathbf{x}}\left(i+1\right)$, for some $i$, $i=1,2,\dots ,n1$. 
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
Accuracy depends on the value of $\rho $ and the position of the $x$ values. The values of ${x}_{i}{x}_{i1}$ and ${w}_{i}$ are scaled and $\rho $ is transformed to avoid underflow and overflow problems.
8
Parallelism and Performance
g10abf is not threaded in any implementation.
The time taken by g10abf is of order $n$.
Regression splines with a small
$\left(<n\right)$ number of knots can be fitted by
e02baf and
e02bef.
10
Example
The data, given by
Hastie and Tibshirani (1990), is the age,
${x}_{i}$, and Cpeptide concentration (pmol/ml),
${y}_{i}$, from a study of the factors affecting insulindependent diabetes mellitus in children. The data is input, reduced to a strictly ordered set by
g10zaf and a series of splines fit using a range of values for the smoothing parameter,
$\rho $.
10.1
Program Text
Program Text (g10abfe.f90)
10.2
Program Data
Program Data (g10abfe.d)
10.3
Program Results
Program Results (g10abfe.r)