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: nag_smooth_kerndens_gauss (g10bb)

Purpose

nag_smooth_kerndens_gauss2 (g10bb) performs kernel density estimation using a Gaussian kernel.

Syntax

[window, slo, shi, smooth, t, rcomm, ifail] = g10bb(x, fcall, rcomm, 'n', n, 'wtype', wtype, 'window', window, 'slo', slo, 'shi', shi, 'ns', ns)
[window, slo, shi, smooth, t, rcomm, ifail] = nag_smooth_kerndens_gauss(x, fcall, rcomm, 'n', n, 'wtype', wtype, 'window', window, 'slo', slo, 'shi', shi, 'ns', ns)

Description

Given a sample of $n$ observations, ${x}_{1},{x}_{2},\dots ,{x}_{n}$, from a distribution with unknown density function, $f\left(x\right)$, an estimate of the density function, $\stackrel{^}{f}\left(x\right)$, may be required. The simplest form of density estimator is the histogram. This may be defined by:
 $f^ x = 1nh nj , a + j-1 h < x < a + j h , j=1,2,…,ns ,$
where ${n}_{j}$ is the number of observations falling in the interval $a+\left(j-1\right)h$ to $a+jh$, $a$ is the lower bound to the histogram, $b={n}_{s}h$ is the upper bound and ${n}_{s}$ is the total number of intervals. The value $h$ is known as the window width. To produce a smoother density estimate a kernel method can be used. A kernel function, $K\left(t\right)$, satisfies the conditions:
 $∫-∞∞Ktdt=1 and Kt≥0.$
The kernel density estimator is then defined as
 $f^x=1nh ∑i= 1nK x-xih .$
The choice of $K$ is usually not important but to ease the computational burden use can be made of the Gaussian kernel defined as
 $Kt=12πe-t2/2.$
The smoothness of the estimator depends on the window width $h$. The larger the value of $h$ the smoother the density estimate. The value of $h$ can be chosen by examining plots of the smoothed density for different values of $h$ or by using cross-validation methods (see Silverman (1990)).
Silverman (1982) and Silverman (1990) show how the Gaussian kernel density estimator can be computed using a fast Fourier transform (FFT). In order to compute the kernel density estimate over the range $a$ to $b$ the following steps are required.
 (i) Discretize the data to give ${n}_{s}$ equally spaced points ${t}_{l}$ with weights ${\xi }_{l}$ (see Jones and Lotwick (1984)). (ii) Compute the FFT of the weights ${\xi }_{l}$ to give ${Y}_{l}$. (iii) Compute ${\zeta }_{l}={e}^{-\frac{1}{2}{h}^{2}{s}_{l}^{2}}{Y}_{l}$ where ${s}_{l}=2\pi l/\left(b-a\right)$. (iv) Find the inverse FFT of ${\zeta }_{l}$ to give $\stackrel{^}{f}\left(x\right)$.
To compute the kernel density estimate for further values of $h$ only steps (iii) and (iv) need be repeated.

References

Jones M C and Lotwick H W (1984) Remark AS R50. A remark on algorithm AS 176. Kernel density estimation using the Fast Fourier Transform Appl. Statist. 33 120–122
Silverman B W (1982) Algorithm AS 176. Kernel density estimation using the fast Fourier transform Appl. Statist. 31 93–99
Silverman B W (1990) Density Estimation Chapman and Hall

Parameters

Compulsory Input Parameters

1:     $\mathrm{x}\left({\mathbf{n}}\right)$ – double array
${x}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
If ${\mathbf{fcall}}=0$, x must be unchanged since the last call to nag_smooth_kerndens_gauss2 (g10bb).
2:     $\mathrm{fcall}$int64int32nag_int scalar
If ${\mathbf{fcall}}=1$ then the values of ${Y}_{l}$ are to be calculated by this call to nag_smooth_kerndens_gauss2 (g10bb), otherwise it is assumed that the values of ${Y}_{l}$ were calculated by a previous call to this routine and the relevant information is stored in rcomm.
Constraint: ${\mathbf{fcall}}=0$ or $1$.
3:     $\mathrm{rcomm}\left({\mathbf{ns}}+20\right)$ – double array
Communication array, used to store information between calls to nag_smooth_kerndens_gauss2 (g10bb).
If ${\mathbf{fcall}}=0$, rcomm must be unchanged since the last call to nag_smooth_kerndens_gauss2 (g10bb).

Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the dimension of the array x.
$n$, the number of observations in the sample.
If ${\mathbf{fcall}}=0$, n must be unchanged since the last call to nag_smooth_kerndens_gauss2 (g10bb).
Constraint: ${\mathbf{n}}>0$.
2:     $\mathrm{wtype}$int64int32nag_int scalar
Suggested value: ${\mathbf{wtype}}=2$ and ${\mathbf{window}}=1.0$.
Default: $2$
How the window width, $h$, is to be calculated:
${\mathbf{wtype}}=1$
$h$ is supplied in window.
${\mathbf{wtype}}=2$
$h$ is to be calculated from the data, with
 $h=m× 0.9× minq75-q25,σ n0.2$
where ${q}_{75}-{q}_{25}$ is the inter-quartile range and $\sigma$ the standard deviation of the sample, $x$, and $m$ is a multipler supplied in window. The $25%$ and $75%$ quartiles, ${q}_{25}$ and ${q}_{75}$, are calculated using nag_stat_quantiles (g01am). This is the "rule-of-thumb" suggested by Silverman (1990).
Constraint: ${\mathbf{wtype}}=1$ or $2$.
3:     $\mathrm{window}$ – double scalar
Suggested value: ${\mathbf{window}}=1.0$ and ${\mathbf{wtype}}=2$.
Default: $1.0$
If ${\mathbf{wtype}}=1$, then $h$, the window width. Otherwise, $m$, the multiplier used in the calculation of $h$.
Constraint: ${\mathbf{window}}>0.0$.
4:     $\mathrm{slo}$ – double scalar
Suggested value: ${\mathbf{slo}}=3.0$ and ${\mathbf{shi}}=0.0$ which would cause $a$ and $b$ to be set $3$ window widths below and above the lowest and highest data points respectively.
Default: $3.0$
If ${\mathbf{slo}}<{\mathbf{shi}}$ then $a$, the lower limit of the interval on which the estimate is calculated. Otherwise, $a$ and $b$, the lower and upper limits of the interval, are calculated as follows:
 $a = minixi-slo×h b = maxixi+slo×h$
where $h$ is the window width.
For most applications $a$ should be at least three window widths below the lowest data point.
If ${\mathbf{fcall}}=0$, slo must be unchanged since the last call to nag_smooth_kerndens_gauss2 (g10bb).
5:     $\mathrm{shi}$ – double scalar
Default: $0.0$
If ${\mathbf{slo}}<{\mathbf{shi}}$ then $b$, the upper limit of the interval on which the estimate is calculated. Otherwise a value for $b$ is calculated from the data as stated in the description of slo and the value supplied in shi is not used.
For most applications $b$ should be at least three window widths above the highest data point.
If ${\mathbf{fcall}}=0$, shi must be unchanged since the last call to nag_smooth_kerndens_gauss2 (g10bb).
6:     $\mathrm{ns}$int64int32nag_int scalar
Default: $512$
${n}_{s}$, the number of points at which the estimate is calculated.
If ${\mathbf{fcall}}=0$, ns must be unchanged since the last call to nag_smooth_kerndens_gauss2 (g10bb).
Constraints:
• ${\mathbf{ns}}\ge 2$;
• The largest prime factor of ns must not exceed $19$, and the total number of prime factors of ns, counting repetitions, must not exceed $20$.

Output Parameters

1:     $\mathrm{window}$ – double scalar
Suggested value: ${\mathbf{window}}=1.0$ and ${\mathbf{wtype}}=2$.
Default: $1.0$
$h$, the window width actually used.
2:     $\mathrm{slo}$ – double scalar
Suggested value: ${\mathbf{slo}}=3.0$ and ${\mathbf{shi}}=0.0$ which would cause $a$ and $b$ to be set $3$ window widths below and above the lowest and highest data points respectively.
Default: $3.0$
$a$, the lower limit actually used.
3:     $\mathrm{shi}$ – double scalar
Default: $0.0$
$b$, the upper limit actually used.
4:     $\mathrm{smooth}\left({\mathbf{ns}}\right)$ – double array
$\stackrel{^}{f}\left({t}_{\mathit{l}}\right)$, for $\mathit{l}=1,2,\dots ,{n}_{s}$, the ${n}_{s}$ values of the density estimate.
5:     $\mathrm{t}\left({\mathbf{ns}}\right)$ – double array
${t}_{\mathit{l}}$, for $\mathit{l}=1,2,\dots ,{n}_{s}$, the points at which the estimate is calculated.
6:     $\mathrm{rcomm}\left({\mathbf{ns}}+20\right)$ – double array
The last ns elements of rcomm contain the fast Fourier transform of the weights of the discretized data, that is ${\mathbf{rcomm}}\left(\mathit{l}+20\right)={Y}_{\mathit{l}}$, for $\mathit{l}=1,2,\dots ,{n}_{s}$.
7:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Note: nag_smooth_kerndens_gauss2 (g10bb) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{ifail}}=11$
Constraint: ${\mathbf{n}}>0$.
${\mathbf{ifail}}=12$
Constraint: if ${\mathbf{fcall}}=0$, n must be unchanged since previous call.
${\mathbf{ifail}}=31$
Constraint: ${\mathbf{wtype}}=1$ or $2$.
${\mathbf{ifail}}=41$
Constraint: ${\mathbf{window}}>0.0$.
${\mathbf{ifail}}=51$
Constraint: if ${\mathbf{fcall}}=0$, slo must be unchanged since previous call.
W  ${\mathbf{ifail}}=61$
slo is not at least three window widths below the lowest data point or shi is not at least three window widths above the highest data point. All output values have been returned.
${\mathbf{ifail}}=62$
Constraint: if ${\mathbf{fcall}}=0$, shi must be unchanged since previous call.
${\mathbf{ifail}}=71$
Constraint: ${\mathbf{ns}}\ge 2$.
${\mathbf{ifail}}=72$
Constraint: largest prime factor of ns must not exceed $19$.
${\mathbf{ifail}}=73$
Constraint: total number of prime factors of ns must not exceed $20$.
${\mathbf{ifail}}=74$
Constraint: if ${\mathbf{fcall}}=0$, ns must be unchanged since previous call.
${\mathbf{ifail}}=101$
Constraint: ${\mathbf{fcall}}=0$ or $1$.
${\mathbf{ifail}}=111$
rcomm has been corrupted between calls.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

Accuracy

See Jones and Lotwick (1984) for a discussion of the accuracy of this method.

The time for computing the weights of the discretized data is of order $n$, while the time for computing the FFT is of order ${n}_{s}\mathrm{log}\left({n}_{s}\right)$, as is the time for computing the inverse of the FFT.

Example

Data is read from a file and the density estimated. The first $20$ values are then printed.
```function g10bb_example

fprintf('g10bb example results\n\n');

% kernel density estimation from 100 values
x  = [ 0.114 -0.232 -0.570  1.853 -0.994 ...
-0.374 -1.028  0.509  0.881 -0.453 ...
0.588 -0.625 -1.622 -0.567  0.421 ...
-0.475  0.054  0.817  1.015  0.608 ...
-1.353 -0.912 -1.136  1.067  0.121 ...
-0.075 -0.745  1.217 -1.058 -0.894 ...
1.026 -0.967 -1.065  0.513  0.969 ...
0.582 -0.985  0.097  0.416 -0.514 ...
0.898 -0.154  0.617 -0.436 -1.212 ...
-1.571  0.210 -1.101  1.018 -1.702 ...
-2.230 -0.648 -0.350  0.446 -2.667 ...
0.094 -0.380 -2.852 -0.888 -1.481 ...
-0.359 -0.554  1.531  0.052 -1.715 ...
1.255 -0.540  0.362 -0.654 -0.272 ...
-1.810  0.269 -1.918  0.001  1.240 ...
-0.368 -0.647 -2.282  0.498  0.001 ...
-3.059 -1.171  0.566  0.948  0.925 ...
0.825  0.130  0.930  0.523  0.443 ...
-0.649  0.554 -2.823  0.158 -1.180 ...
0.610  0.877  0.791 -0.078  1.412];

% Calculate window width from data.
wtype = int64(2);

% First Call
fcall = int64(1);
ns    = 512;
rcomm = zeros(ns+20,1);

% Perform kernel density estimation
[window, slo, shi, smooth, t, rcomm, ifail] = ...
g10bb( ...
x, fcall, rcomm, 'wtype',wtype);

% Display the results
fprintf('Window Width Used = %11.4e\n', window);
fprintf('Interval = (%11.4e,%11.4e)\n\n', slo, shi);
fprintf('First %2d output values:\n\n',20);
fprintf('    Time point      Density estimate\n');
fprintf('    ----------      ----------------\n');
fprintf(' %13.4f     %13.4e\n', [t(1:20), smooth(1:20)]')

fig1 = figure;
plot(t,smooth);
title('Gaussian Kernel Density Estimation');
xlabel('t');
ylabel('Density Estimate');
wind_leg = sprintf('window = %7.4f',window);
legend(wind_leg);
legend('boxoff');

```
```g10bb example results

Window Width Used =  3.7638e-01
Interval = (-4.1882e+00, 2.9822e+00)

First 20 output values:

Time point      Density estimate
----------      ----------------
-4.1811        3.8281e-06
-4.1671        4.0305e-06
-4.1531        4.4233e-06
-4.1391        5.0212e-06
-4.1251        5.8461e-06
-4.1111        6.9279e-06
-4.0971        8.3048e-06
-4.0831        1.0025e-05
-4.0691        1.2145e-05
-4.0551        1.4736e-05
-4.0411        1.7881e-05
-4.0271        2.1677e-05
-4.0131        2.6239e-05
-3.9991        3.1700e-05
-3.9851        3.8214e-05
-3.9711        4.5960e-05
-3.9571        5.5141e-05
-3.9431        6.5990e-05
-3.9291        7.8775e-05
-3.9151        9.3796e-05
```
This plot shows the estimated density function for the example data for several window widths.