NAG CL Interfaceg10abc (fit_​spline)

Settings help

CL Name Style:

1Purpose

g10abc fits a cubic smoothing spline for a given smoothing parameter.

2Specification

 #include
 void g10abc (Nag_SmoothFitType mode, Integer n, const double x[], const double y[], const double wt[], double rho, double yhat[], double c[], Integer pdc, double *rss, double *df, double res[], double h[], double comm[], NagError *fail)
The function may be called by the names: g10abc, nag_smooth_fit_spline or nag_smooth_spline_fit.

3Description

g10abc 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 real-valued solution function $f$, with absolutely continuous first derivative and squared-integrable second derivative, which minimizes:
 $∑i=1nwi(yi-f(xi))2+ρ ∫-∞∞(f′′(x))2dx,$
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, $\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=wi(yi-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$, ${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) cross-validation can be minimized to give $\rho$; see g10acc for further details.
g10abc 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 g10zac.
The computation is split into three phases.
1. (i)Compute matrices needed to fit spline.
2. (ii)Fit spline for a given value of $\rho$.
3. (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).

4References

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

5Arguments

1: $\mathbf{mode}$Nag_SmoothFitType Input
On entry: indicates in which mode the function is to be used.
${\mathbf{mode}}=\mathrm{Nag_SmoothFitPartial}$
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 g10abc or when the coefficients are not required.
${\mathbf{mode}}=\mathrm{Nag_SmoothFitQuick}$
Fitting only is performed. Initialization must have been performed previously by a call to g10abc with ${\mathbf{mode}}=\mathrm{Nag_SmoothFitPartial}$. This quick fit may be called repeatedly with different values of rho without re-initialization.
${\mathbf{mode}}=\mathrm{Nag_SmoothFitFull}$
Initialization and full fitting is performed and the function coefficients are calculated.
Constraint: ${\mathbf{mode}}=\mathrm{Nag_SmoothFitPartial}$, $\mathrm{Nag_SmoothFitQuick}$ or $\mathrm{Nag_SmoothFitFull}$.
2: $\mathbf{n}$Integer Input
On entry: $n$, the number of distinct observations.
Constraint: ${\mathbf{n}}\ge 3$.
3: $\mathbf{x}\left[{\mathbf{n}}\right]$const double Input
On entry: the distinct and ordered values ${x}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
Constraint: ${\mathbf{x}}\left[\mathit{i}-1\right]<{\mathbf{x}}\left[\mathit{i}\right]$, for $\mathit{i}=1,2,\dots ,n-1$.
4: $\mathbf{y}\left[{\mathbf{n}}\right]$const double Input
On entry: the values ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
5: $\mathbf{wt}\left[\mathit{dim}\right]$const double Input
Note: the dimension, dim, of the array wt must be at least ${\mathbf{n}}$.
On entry: optionally, the $n$ weights.
If weights are not provided then wt must be set to NULL, in which case unit weights are assumed.
Constraint: ${\mathbf{wt}}\left[\mathit{i}-1\right]>0.0$, for $\mathit{i}=1,2,\dots ,n$.
6: $\mathbf{rho}$double Input
On entry: $\rho$, the smoothing parameter.
Constraint: ${\mathbf{rho}}\ge 0.0$.
7: $\mathbf{yhat}\left[{\mathbf{n}}\right]$double Output
On exit: the fitted values, ${\stackrel{^}{y}}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
8: $\mathbf{c}\left[{\mathbf{pdc}}×3\right]$double Input/Output
Note: the $\left(i,j\right)$th element of the matrix $C$ is stored in ${\mathbf{c}}\left[\left(j-1\right)×{\mathbf{pdc}}+i-1\right]$.
On entry: if ${\mathbf{mode}}=\mathrm{Nag_SmoothFitQuick}$, c must be unaltered from the previous call to g10abc with ${\mathbf{mode}}=\mathrm{Nag_SmoothFitPartial}$. Otherwise c need not be set.
On exit: if ${\mathbf{mode}}=\mathrm{Nag_SmoothFitFull}$, c contains the spline coefficients. More precisely, the value of the spline at $t$ is given by $\left(\left({\mathbf{c}}\left[2×{\mathbf{pdc}}+i-1\right]×d+{\mathbf{c}}\left[1×{\mathbf{pdc}}+i-1\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}$.
If ${\mathbf{mode}}=\mathrm{Nag_SmoothFitPartial}$ or $\mathrm{Nag_SmoothFitQuick}$, c contains information that will be used in a subsequent call to g10abc with ${\mathbf{mode}}=\mathrm{Nag_SmoothFitQuick}$.
9: $\mathbf{pdc}$Integer Input
On entry: the stride separating matrix row elements in the array c.
Constraint: ${\mathbf{pdc}}\ge {\mathbf{n}}-1$.
10: $\mathbf{rss}$double * Output
On exit: the (weighted) residual sum of squares.
11: $\mathbf{df}$double * Output
On exit: the residual degrees of freedom.
12: $\mathbf{res}\left[{\mathbf{n}}\right]$double Output
On exit: the (weighted) residuals, ${r}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
13: $\mathbf{h}\left[{\mathbf{n}}\right]$double Output
On exit: the leverages, ${h}_{\mathit{i}\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
14: $\mathbf{comm}\left[9×{\mathbf{n}}+14\right]$double Communication Array
On entry: if ${\mathbf{mode}}=\mathrm{Nag_SmoothFitQuick}$, comm must be unaltered from the previous call to g10abc with ${\mathbf{mode}}=\mathrm{Nag_SmoothFitPartial}$. Otherwise comm need not be set.
On exit: if ${\mathbf{mode}}=\mathrm{Nag_SmoothFitPartial}$ or $\mathrm{Nag_SmoothFitQuick}$, comm contains information that will be used in a subsequent call to g10abc with ${\mathbf{mode}}=\mathrm{Nag_SmoothFitQuick}$.
15: $\mathbf{fail}$NagError * Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

6Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_ARRAY_SIZE
On entry, ${\mathbf{pdc}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdc}}\ge {\mathbf{n}}-1$.
On entry, argument $⟨\mathit{\text{value}}⟩$ had an illegal value.
NE_INT
On entry, ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{n}}\ge 3$.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
NE_NOT_STRICTLY_INCREASING
On entry, x is not a strictly ordered array.
NE_REAL
On entry, ${\mathbf{rho}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{rho}}\ge 0.0$.
NE_WEIGHTS_NOT_POSITIVE
On entry, at least one element of ${\mathbf{wt}}\le 0.0$.

7Accuracy

Accuracy 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.

8Parallelism and Performance

g10abc is not threaded in any implementation.

The time taken by g10abc is of order $n$.
Regression splines with a small $\left( number of knots can be fitted by e02bac and e02bec.

10Example

The data, given by Hastie and Tibshirani (1990), is 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 g10zac and a series of splines fit using a range of values for the smoothing parameter, $\rho$.

10.1Program Text

Program Text (g10abce.c)

10.2Program Data

Program Data (g10abce.d)

10.3Program Results

Program Results (g10abce.r)