Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox Chapter IntroductionG10 — Smoothing in Statistics

## Scope of the Chapter

This chapter is concerned with methods for smoothing data. Included are methods for density estimation, smoothing time series data, and statistical applications of splines. These methods may also be viewed as nonparametric modelling.

## Background to the Problems

### Smoothing Methods

Many of the methods used in statistics involve fitting a model, the form of which is determined by a small number of parameters, for example, a distribution model like the gamma distribution, a linear regression model or an autoregression model in time series. In these cases the fitting involves the estimation of the small number of parameters from the data. In modelling data with these models there are two important stages in addition to the estimation of the parameters; these are the identification of a suitable model, for example, the selection of a gamma distribution rather than a Weibull distribution, and the checking to see if the fitted model adequately fits the data. While these parametric models can be fairly flexible, they will not adequately fit all datasets, especially if the number of parameters is to be kept small.
Alternative models based on smoothing can be used. These models will not be written explicitly in terms of parameters. They are sufficiently flexible for a much wider range of situations than parametric models. The main requirement for such a model to be suitable is that the underlying models would be expected to be smooth, so excluding those situations where, for example, a step function would be expected.
These smoothing methods can be used in a variety of ways, for example:
1. producing smoothed plots to aid understanding;
2. identifying of a suitable parametric model from the shape of the smoothed data;
3. eliminating complex effects that are not of direct interest so that attention can be focused on the effects of interest.
Several smoothing techniques make use of a smoothing parameter which can be either chosen by you or estimated from the data. The smoothing parameter balances the two criterion of smoothness of the fitted model and the closeness of the fit of the model to the data. Generally, the larger the smoothing parameter is, the smoother the fitted model will be, but for small values of the smoothing parameter the model will closely follow the data, and for large values the fit will be poorer.
The smoothing parameter can be either chosen using previous experience of a suitable value for such data, or estimated from the data. The estimation can be either formal, using a criterion such as the cross-validation, or informal by trying different values and examining the result by means of suitable graphs.
Smoothing methods can be used in three important areas of of statistics: regression modelling, distribution modelling and time series modelling.

### Smoothing Splines and Regression Models

For a set of n$n$ observations (yi,xi${y}_{i},{x}_{i}$), i = 1,2,,n$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 function, f$f$, with continuous first derivative which minimizes
 n ∞ ∑ wi(yi − f(xi))2 + ρ ∫ (f ′ ′ (x))2dx, i = 1 − ∞
$∑i=1nwi (yi-f(xi)) 2+ρ∫-∞∞(f′′(x))2dx,$
where wi${w}_{i}$ is the (optional) weight for the i$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.
Splines are linear smoothers since the fitted values, = (1,2,,n)T $\stackrel{^}{y}={\left({\stackrel{^}{y}}_{1},{\stackrel{^}{y}}_{2},\dots ,{\stackrel{^}{y}}_{n}\right)}^{\mathrm{T}}$, can be written as a linear function of the observed values y = (y1,y2,,yn)T $y={\left({y}_{1},{y}_{2},\dots ,{y}_{n}\right)}^{\mathrm{T}}$, that is,
 ŷ = Hy $y^=Hy$
for a matrix H$H$. The degrees of freedom for the spline is trace(H)$\mathrm{trace}\left(H\right)$ giving residual degrees of freedom
 n trace(I − H) = ∑ (1 − hii). i = 1
$trace(I-H)=∑i=1n(1-hii).$
The diagonal elements of H$H$, hii${h}_{ii}$, are the leverages.
The parameter ρ$\rho$ can be estimated in a number of ways.
1. The degrees of freedom for the spline can be specified, i.e., find ρ$\rho$ such that trace(H) = ν0$\mathrm{trace}\left(H\right)={\nu }_{0}$ for given ν0${\nu }_{0}$.
2. Minimize the cross-validation (CV), i.e., find ρ$\rho$ such that the CV is minimized, where
 n CV = 1/n ∑ ((ri)/(1 − hii))2. i = 1
$CV=1n∑i=1n (ri1-hii ) 2.$
3. Minimize generalized cross-validation (GCV), i.e., find ρ$\rho$ such that the GCV is minimized, where
GCV = n (( ∑ i = 1nri2)/(
 (n ) ∑ (1 − hii)i = 1 2
)) .
$GCV=n (∑i=1nri2 (∑i=1n (1-hii) ) 2 ) .$

### Density Estimation

The object of density estimation is to produce from a set of observations a smooth nonparametric estimate of the unknown density function from which the observations were drawn. That is, given a sample of n$n$ observations, x1${x}_{1}$, x2,,xn${x}_{2},\dots ,{x}_{n}$, from a distribution with unknown density function, f(x)$f\left(x\right)$, find an estimate of the density function, (x)$\stackrel{^}{f}\left(x\right)$. The simplest form of density estimator is the histogram; this may be defined by
 f̂(x) = 1/(nh)nj;  a + (j − 1)h < x < a + jh;  j = 1,2, … ,ns, $f^(x)=1nh nj; a+(j-1)h
where nj${n}_{j}$ is the number of observations falling in the interval a + (j1)h$a+\left(j-1\right)h$ to a + jh$a+jh$, a$a$ is the lower bound of the histogram and b = nsh$b={n}_{s}h$ is the upper bound. The value h$h$ is known as the window width. A simple development of this estimator would be the running histogram estimator
 f̂(x) = 1/(2nh)nx;  a ≤ x ≤ b, $f^(x)=12nh nx; a≤x≤b,$
where nx${n}_{x}$ is the number of observations falling in the interval [xh : x + h]$\left[x-h:x+h\right]$. This estimator can be written as
 n f̂(x) = 1/(nh) ∑ w((x − xi)/h) i = 1
$f^(x)=1nh ∑i=1nw (x-xih)$
for a function w$w$ where
 w(x) = (1/2) if − 1 < x < 1 = 0 otherwise.
$w(x) =12 if-1
The function w$w$ can be considered as a kernel function. To produce a smoother density estimate, the kernel function, K(t)$K\left(t\right)$, which satisfies the following conditions can be used:
 ∞ ∫ K(t)dt = 1and ​K(t) ≥ 0.0. − ∞
$∫-∞∞K(t)dt=1 and ​K(t)≥0.0.$
The kernel density estimator is therefore defined as
 n f̂(x) = 1/(nh) ∑ K((x − xi)/h). i = 1
$f^(x)=1nh ∑i= 1nK (x-xih) .$
The choice of K( · )$K\left(·\right)$ is usually not important, but to ease computational burden use can be made of Gaussian kernel defined as
 K(t) = 1/(sqrt(2π))e − t2 / 2. $K(t)=12πe-t2/2.$
The smoothness of the estimator, (x)$\stackrel{^}{f}\left(x\right)$, depends on the window width, h$h$. In general, the larger the value h$h$ is, the smoother the resulting density estimate is. There is, however, the problem of oversmoothing when the value of h$h$ is too large and essential features of the distribution function are removed. For example, if the distribution was bimodal, a large value of h$h$ may result in a unimodal estimate. The value of h$h$ has to be chosen such that the essential shape of the distribution is retained while effects due only to the observed sample are smoothed out. The choice of h$h$ can be aided by looking at plots of the density estimate for different values of h$h$, or by using cross-validation methods; see Silverman (1990).
Silverman (1990) shows how the Gaussian kernel density estimator can be computed using a fast Fourier transform (FFT).

### Smoothers for Time Series

If the data consists of a sequence of n$n$ observations recorded at equally spaced intervals, usually a time series, several robust smoothers are available. The fitted curve is intended to be robust to any outlying observations in the sequence, hence the techniques employed primarily make use of medians rather than means. These ideas come from the field of exploratory data analysis (EDA); see Tukey (1977) and Velleman and Hoaglin (1981). The smoothers are based on the use of running medians to summarise overlapping segments; these provide a simple but flexible curve.
In EDA terminology, the fitted curve and the residuals are called the smooth and the rough respectively, so that
 Data = Smooth + Rough. $Data=Smooth+Rough.$
Using the notation of Tukey, one of the smoothers commonly used is 4253H,twice. This consists of a running median of 4$4$, then 2$2$, then 5$5$, then 3$3$. This is then followed by what is known as hanning. Hanning is a running weighted mean, the weights being 1 / 4$1/4$, 1 / 2$1/2$ and 1 / 4$1/4$. The result of this smoothing is then ‘reroughed’. This involves computing residuals from the computed fit, applying the same smoother to the residuals and adding the result to the smooth of the first pass.

## Recommendations on Choice and Use of Available Functions

The following functions fit smoothing splines:
• nag_smooth_fit_spline (g10ab) computes a cubic smoothing spline for a given value of the smoothing parameter. The results returned include the values of leverages and the coefficients of the cubic spline. Options allow only parts of the computation to be performed when the function is used to estimate the value of the smoothing parameter or as when it is part of an iterative procedure such as that used in fitting generalized additive models; see Hastie and Tibshirani (1990).
• nag_smooth_fit_spline_parest (g10ac) estimates the value of the smoothing parameter using one of three criteria and fits the cubic smoothing spline using that value.
nag_smooth_fit_spline (g10ab) and nag_smooth_fit_spline_parest (g10ac) require the xi${x}_{i}$ to be strictly increasing. If two or more observations have the same xi${x}_{i}$-value then they should be replaced by a single observation with yi${y}_{i}$ equal to the (weighted) mean of the y$y$ values and weight, wi${w}_{i}$, equal to the sum of the weights. This operation can be performed by nag_smooth_data_order (g10za).
The following function produces an estimate of the density function:
The following function produces a smoothed estimate for a time series:
The following service function is also available:
• nag_smooth_data_order (g10za) orders and weights the (x,y)$\left(x,y\right)$ input data to produce a dataset strictly monotonic in x$x$.

## Functionality Index

 Compute smoothed data sequence,
 running median smoothers nag_smooth_data_runningmedian (g10ca)
 Fit cubic smoothing spline,
 smoothing parameter estimated nag_smooth_fit_spline_parest (g10ac)
 smoothing parameter given nag_smooth_fit_spline (g10ab)
 Kernel density estimation,
 Gaussian kernel nag_smooth_kerndens_gauss (g10ba)
 Reorder data to give ordered distinct observations nag_smooth_data_order (g10za)

## References

Hastie T J and Tibshirani R J (1990) Generalized Additive Models Chapman and Hall
Silverman B W (1990) Density Estimation Chapman and Hall
Tukey J W (1977) Exploratory Data Analysis Addison–Wesley
Velleman P F and Hoaglin D C (1981) Applications, Basics, and Computing of Exploratory Data Analysis Duxbury Press, Boston, MA

Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013