G10 Chapter Contents
G10 Chapter Introduction
NAG Library Manual

# NAG Library Routine DocumentG10ACF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

## 1  Purpose

G10ACF estimates the values of the smoothing parameter and fits a cubic smoothing spline to a set of data.

## 2  Specification

 SUBROUTINE G10ACF ( METHOD, WEIGHT, N, X, Y, WT, YHAT, C, LDC, RSS, DF, RES, H, CRIT, RHO, U, TOL, MAXCAL, WK, IFAIL)
 INTEGER N, LDC, MAXCAL, IFAIL REAL (KIND=nag_wp) X(N), Y(N), WT(*), YHAT(N), C(LDC,3), RSS, DF, RES(N), H(N), CRIT, RHO, U, TOL, WK(7*(N+2)) CHARACTER(1) METHOD, WEIGHT

## 3  Description

For a set of $n$ observations $\left({x}_{\mathit{i}},{y}_{\mathit{i}}\right)$, for $\mathit{i}=1,2,\dots ,n$, the spline provides a flexible smooth function for situations in which a simple polynomial or nonlinear regression model is not suitable.
Cubic smoothing splines arise as the unique real-valued solution function $f$, with absolutely continuous first derivative and squared-integrable second derivative, which minimizes
 $∑i=1nwiyi-fxi2+ρ ∫-∞∞f′′x2dx,$
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 fitted see Hutchinson and de Hoog (1985) and Reinsch (1967).
The fitted values, $\stackrel{^}{y}={\left({\stackrel{^}{y}}_{1},{\stackrel{^}{y}}_{2},\dots ,{\stackrel{^}{y}}_{n}\right)}^{\mathrm{T}}$, and weighted residuals, ${r}_{i}$, can be written as:
 $y^=Hy and ri=wiyi-y^i$
for a matrix $H$. The residual degrees of freedom for the spline is $\mathrm{trace}\left(I-H\right)$ and the diagonal elements of $H$ are the leverages.
The parameter $\rho$ can be estimated in a number of ways.
(i) The degrees of freedom for the spline can be specified, i.e., find $\rho$ such that $\mathrm{trace}\left(H\right)={\nu }_{0}$ for given ${\nu }_{0}$.
(ii) Minimize the cross-validation (CV), i.e., find $\rho$ such that the CV is minimized, where
 $CV=1∑i=1nwi ∑i=1n ri1-hii 2.$
(iii) Minimize the generalized cross-validation (GCV), i.e., find $\rho$ such that the GCV is minimized, where
 $GCV=n2∑i=1nwi ∑i=1nri2 ∑i=1n1-hii 2 .$
G10ACF 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 algorithm is based on Hutchinson (1986). C05AZF is used to solve for $\rho$ given ${\nu }_{0}$ and the method of E04ABF/E04ABA is used to minimize the GCV or CV.

## 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 cross-validation 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  Parameters

1:     METHOD – CHARACTER(1)Input
On entry: indicates whether the smoothing parameter is to be found by minimization of the CV or GCV functions, or by finding the smoothing parameter corresponding to a specified degrees of freedom value.
${\mathbf{METHOD}}=\text{'C'}$
Cross-validation is used.
${\mathbf{METHOD}}=\text{'D'}$
The degrees of freedom are specified.
${\mathbf{METHOD}}=\text{'G'}$
Generalized cross-validation is used.
Constraint: ${\mathbf{METHOD}}=\text{'C'}$, $\text{'D'}$ or $\text{'G'}$.
2:     WEIGHT – CHARACTER(1)Input
On entry: indicates whether user-defined weights are to be used.
${\mathbf{WEIGHT}}=\text{'W'}$
User-defined 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:     N – INTEGERInput
On entry: $n$, the number of observations.
Constraint: ${\mathbf{N}}\ge 3$.
4:     X(N) – 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 ,n-1$.
5:     Y(N) – REAL (KIND=nag_wp) arrayInput
On entry: the values ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
6:     WT($*$) – REAL (KIND=nag_wp) arrayInput
Note: the dimension of the array WT must be at least $1$ if ${\mathbf{WEIGHT}}=\text{'U'}$ and at least ${\mathbf{N}}$ if ${\mathbf{WEIGHT}}=\text{'W'}$.
On entry: if ${\mathbf{WEIGHT}}=\text{'W'}$, WT must contain the $n$ weights.
If ${\mathbf{WEIGHT}}=\text{'U'}$, 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:     YHAT(N) – REAL (KIND=nag_wp) arrayOutput
On exit: the fitted values, ${\stackrel{^}{y}}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
8:     C(LDC,$3$) – REAL (KIND=nag_wp) arrayOutput
On exit: the spline coefficients. More precisely, the value of the spline approximation at $t$ is given by $\left(\left({\mathbf{C}}\left(i,3\right)×d+{\mathbf{C}}\left(i,2\right)\right)×d+{\mathbf{C}}\left(i,1\right)\right)×d+{\stackrel{^}{y}}_{i}$, where ${x}_{i}\le t<{x}_{i+1}$ and $d=t-{x}_{i}$.
9:     LDC – INTEGERInput
On entry: the first dimension of the array C as declared in the (sub)program from which G10ACF is called.
Constraint: ${\mathbf{LDC}}\ge {\mathbf{N}}-1$.
On exit: the (weighted) residual sum of squares.
11:   DF – REAL (KIND=nag_wp)Output
On exit: the residual degrees of freedom. If ${\mathbf{METHOD}}=\text{'D'}$ this will be $n-{\mathbf{CRIT}}$ to the required accuracy.
12:   RES(N) – REAL (KIND=nag_wp) arrayOutput
On exit: the (weighted) residuals, ${r}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
13:   H(N) – REAL (KIND=nag_wp) arrayOutput
On exit: the leverages, ${h}_{\mathit{i}\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
14:   CRIT – REAL (KIND=nag_wp)Input/Output
On entry: if ${\mathbf{METHOD}}=\text{'D'}$, the required degrees of freedom for the spline.
If ${\mathbf{METHOD}}=\text{'C'}$ or $\text{'G'}$, CRIT need not be set.
Constraint: $2.0<{\mathbf{CRIT}}\le {\mathbf{N}}$.
On exit: if ${\mathbf{METHOD}}=\text{'C'}$, the value of the cross-validation, or if ${\mathbf{METHOD}}=\text{'G'}$, the value of the generalized cross-validation function, evaluated at the value of $\rho$ returned in RHO.
15:   RHO – REAL (KIND=nag_wp)Output
On exit: the smoothing parameter, $\rho$.
16:   U – REAL (KIND=nag_wp)Input
On entry: the upper bound on the smoothing parameter. If ${\mathbf{U}}\le {\mathbf{TOL}}$, ${\mathbf{U}}=1000.0$ will be used instead. See Section 8 for details on how this parameter is used.
17:   TOL – REAL (KIND=nag_wp)Input
On entry: the accuracy to which the smoothing parameter RHO is required. TOL should preferably be not much less than $\sqrt{\epsilon }$, where $\epsilon$ is the machine precision. If ${\mathbf{TOL}}<\epsilon$, ${\mathbf{TOL}}=\sqrt{\epsilon }$ will be used instead.
18:   MAXCAL – INTEGERInput
On entry: the maximum number of spline evaluations to be used in finding the value of $\rho$. If ${\mathbf{MAXCAL}}<3$, ${\mathbf{MAXCAL}}=100$ will be used instead.
19:   WK($7×\left({\mathbf{N}}+2\right)$) – REAL (KIND=nag_wp) arrayWorkspace
20:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to $0$, $-1\text{​ or ​}1$. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction 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 parameter, 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}}={\mathbf{0}}$ or $-{\mathbf{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{METHOD}}\ne \text{'C'}$, $\text{'G'}$ or $\text{'D'}$, or ${\mathbf{WEIGHT}}\ne \text{'W'}$ or $\text{'U'}$, or ${\mathbf{METHOD}}=\text{'D'}$ and ${\mathbf{CRIT}}\le 2.0$, or ${\mathbf{METHOD}}=\text{'D'}$ and ${\mathbf{CRIT}}>{\mathbf{N}}$.
${\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 ,n-1$.
${\mathbf{IFAIL}}=4$
${\mathbf{METHOD}}=\text{'D'}$ and the required value of $\rho$ for specified degrees of $\mathrm{freedom}>U$. Try a larger value of U; see Section 8.
${\mathbf{IFAIL}}=5$
${\mathbf{METHOD}}=\text{'D'}$ and the accuracy given by TOL cannot be achieved. Try increasing the value of TOL.
${\mathbf{IFAIL}}=6$
A solution to the accuracy given by TOL has not been achieved in MAXCAL iterations. Try increasing the value of TOL and/or MAXCAL.
${\mathbf{IFAIL}}=7$
${\mathbf{METHOD}}=\text{'C'}$ or $\text{'G'}$ and the optimal value of $\rho >{\mathbf{U}}$. Try a larger value of U; see Section 8.

## 7  Accuracy

When minimizing the cross-validation or generalized cross-validation, the error in the estimate of $\rho$ should be within $±3\left({\mathbf{TOL}}×{\mathbf{RHO}}+{\mathbf{TOL}}\right)$. When finding $\rho$ for a fixed number of degrees of freedom the error in the estimate of $\rho$ should be within $±2×{\mathbf{TOL}}×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{RHO}}\right)$.
Given the value of $\rho$, the accuracy of the fitted spline depends on the value of $\rho$ and the position of the $x$ values. The values of ${x}_{i}-{x}_{i-1}$ and ${w}_{i}$ are scaled and $\rho$ is transformed to avoid underflow and overflow problems.

The time to fit the spline for a given value of $\rho$ is of order $n$.
When finding the value of $\rho$ that gives the required degrees of freedom, the algorithm examines the interval $0.0$ to U. For small degrees of freedom the value of $\rho$ can be large, as in the theoretical case of two degrees of freedom when the spline reduces to a straight line and $\rho$ is infinite. If the CV or GCV is to be minimized then the algorithm searches for the minimum value in the interval $0.0$ to U. If the function is decreasing in that range then the boundary value of U will be returned. In either case, the larger the value of U the more likely is the interval to contain the required solution, but the process will be less efficient.
Regression splines with a small $\left( number of knots can be fitted by E02BAF and E02BEF.

## 9  Example

This example uses the data given by Hastie and Tibshirani (1990), which consists of the age, ${x}_{i}$, and C-peptide concentration (pmol/ml), ${y}_{i}$, from a study of the factors affecting insulin-dependent diabetes mellitus in children. The data is input, reduced to a strictly ordered set by G10ZAF and a spline with $5$ degrees of freedom is fitted by G10ACF. The fitted values and residuals are printed.

### 9.1  Program Text

Program Text (g10acfe.f90)

### 9.2  Program Data

Program Data (g10acfe.d)

### 9.3  Program Results

Program Results (g10acfe.r)