NAG Library Function Document

nag_smooth_spline_fit (g10abc)

void	nag_smooth_spline_fit (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)

3 Description

nag_smooth_spline_fit (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:

\sum_{i = 1}^{n} w_{i} {(y_{i} - f (x_{i}))}^{2} + ρ \int_{- \infty}^{\infty} {(f^{''} (x))}^{2} d x,

where

w_{i}

is the (optional) weight for the

i

th observation and

ρ

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

ρ

weights these two aspects; larger values of

ρ

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} = {({\hat{y}}_{1}, {\hat{y}}_{2}, \dots, {\hat{y}}_{n})}^{T}

, and weighted residuals,

r_{i}

, can be written as

\hat{y} = H y and r_{i} = \sqrt{w_{i}} (y_{i} - {\hat{y}}_{i})

for a matrix

H

. The residual degrees of freedom for the spline is

trace (I - H)

and the diagonal elements of

H

h_{i i}

, are the leverages.

The parameter

ρ

can be chosen in a number of ways. The fit can be inspected for a number of different values of

ρ

. Alternatively the degrees of freedom for the spline, which determines the value of

ρ

, can be specified, or the (generalized) cross-validation can be minimized to give

ρ

; see nag_smooth_spline_estim (g10acc) for further details.

nag_smooth_spline_fit (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 nag_order_data (g10zac).

The computation is split into three phases.

(i)	Compute matrices needed to fit spline.
(ii)	Fit spline for a given value of $ρ$ .
(iii)	Compute spline coefficients.

When fitting the spline for several different values of

ρ

, phase (i) need only be carried out once and then phase (ii) repeated for different values of

ρ

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

1: $mode$ – Nag_SmoothFitTypeInput

On entry: indicates in which mode the function is to be used.

$mode = 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 nag_smooth_spline_fit (g10abc) or when the coefficients are not required.
$mode = Nag_SmoothFitQuick$: Fitting only is performed. Initialization must have been performed previously by a call to nag_smooth_spline_fit (g10abc) with $mode = Nag_SmoothFitPartial$ . This quick fit may be called repeatedly with different values of rho without re-initialization.
$mode = Nag_SmoothFitFull$: Initialization and full fitting is performed and the function coefficients are calculated.

Constraint:

mode = Nag_SmoothFitPartial

Nag_SmoothFitQuick

Nag_SmoothFitFull

2: $n$ – IntegerInput

On entry:

n

, the number of distinct observations.

Constraint:

n \geq 3

3: $x [n]$ – const doubleInput

On entry: the distinct and ordered values

x_{i}

, for

i = 1, 2, \dots, n

Constraint:

x [i - 1] < x [i]

, for

i = 1, 2, \dots, n - 1

4: $y [n]$ – const doubleInput

On entry: the values

y_{i}

, for

i = 1, 2, \dots, n

5: $wt [n]$ – const doubleInput

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:

wt [i - 1] > 0.0

, for

i = 1, 2, \dots, n

6: $rho$ – doubleInput

On entry:

ρ

, the smoothing parameter.

Constraint:

rho \geq 0.0

7: $yhat [n]$ – doubleOutput

On exit: the fitted values,

{\hat{y}}_{i}

, for

i = 1, 2, \dots, n

8: $c [pdc \times 3]$ – doubleInput/Output

Note: the

(i, j)

th element of the matrix

C

is stored in

c [(j - 1) \times pdc + i - 1]

On entry: if

mode = Nag_SmoothFitQuick

, c must be unaltered from the previous call to nag_smooth_spline_fit (g10abc) with

mode = Nag_SmoothFitPartial

. Otherwise c need not be set.

On exit: if

mode = Nag_SmoothFitFull

, c contains the spline coefficients. More precisely, the value of the spline at

t

is given by

((c [2 \times pdc + i - 1] \times d + c [1 \times pdc + i - 1]) \times d + c [i - 1]) \times d + {\hat{y}}_{i}

, where

x_{i} \leq t < x_{i + 1}

and

d = t - x_{i}

mode = Nag_SmoothFitPartial

Nag_SmoothFitQuick

, c contains information that will be used in a subsequent call to nag_smooth_spline_fit (g10abc) with

mode = Nag_SmoothFitQuick

9: $pdc$ – IntegerInput

On entry: the stride separating matrix row elements in the array c.

Constraint:

pdc \geq n - 1

10: $rss$ – double *Output

On exit: the (weighted) residual sum of squares.

11: $df$ – double *Output

On exit: the residual degrees of freedom.

12: $res [n]$ – doubleOutput

On exit: the (weighted) residuals,

r_{i}

, for

i = 1, 2, \dots, n

13: $h [n]$ – doubleOutput

On exit: the leverages,

h_{i i}

, for

i = 1, 2, \dots, n

14: $comm [9 \times n + 14]$ – doubleCommunication Array

On entry: if

mode = Nag_SmoothFitQuick

, comm must be unaltered from the previous call to nag_smooth_spline_fit (g10abc) with

mode = Nag_SmoothFitPartial

. Otherwise comm need not be set.

On exit: if

mode = Nag_SmoothFitPartial

Nag_SmoothFitQuick

, comm contains information that will be used in a subsequent call to nag_smooth_spline_fit (g10abc) with

mode = Nag_SmoothFitQuick

15: $fail$ – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
See Section 3.2.1.2 in the Essential Introduction for further information.
NE_ARRAY_SIZE: On entry, $pdc = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pdc \geq n - 1$ .
NE_BAD_PARAM: On entry, argument $〈value〉$ had an illegal value.
NE_INT: On entry, $n = 〈value〉$ .
Constraint: $n \geq 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.

An unexpected error has been triggered by this function. Please contact NAG.
See Section 3.6.6 in the Essential Introduction for further information.
NE_NO_LICENCE: Your licence key may have expired or may not have been installed correctly.
See Section 3.6.5 in the Essential Introduction for further information.
NE_NOT_STRICTLY_INCREASING: On entry, x is not a strictly ordered array.
NE_REAL: On entry, $rho = 〈value〉$ .
Constraint: $rho \geq 0.0$ .
NE_WEIGHTS_NOT_POSITIVE: On entry, at least one element of $wt \leq 0.0$ .

7 Accuracy

Accuracy depends on the value of

ρ

and the position of the

x

values. The values of

x_{i} - x_{i - 1}

and

w_{i}

are scaled and

ρ

is transformed to avoid underflow and overflow problems.

8 Parallelism and Performance

Not applicable.

9 Further Comments

The time taken by nag_smooth_spline_fit (g10abc) is of order

n

Regression splines with a small

(< n)

number of knots can be fitted by nag_1d_spline_fit_knots (e02bac) and nag_1d_spline_fit (e02bec).

10 Example

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 nag_order_data (g10zac) and a series of splines fit using a range of values for the smoothing parameter,

ρ

NAG Library Function Documentnag_smooth_spline_fit (g10abc)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Function Document

nag_smooth_spline_fit (g10abc)