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 (g10ba)

## Purpose

nag_smooth_kerndens_gauss (g10ba) performs kernel density estimation using a Gaussian kernel.

## Syntax

[smooth, t, fft, ifail] = g10ba(x, window, slo, shi, usefft, fft, 'n', n, 'ns', ns)
[smooth, t, fft, ifail] = nag_smooth_kerndens_gauss(x, window, slo, shi, usefft, fft, 'n', n, 'ns', ns)

## Description

Given a sample of n$n$ observations, x1,x2,,xn${x}_{1},{x}_{2},\dots ,{x}_{n}$, from a distribution with unknown density function, f(x)$f\left(x\right)$, an estimate of the density function, (x)$\stackrel{^}{f}\left(x\right)$, may be required. 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 + j h ,   j = 1,2, … ,ns , $f^ (x) = 1nh nj , a + (j-1) h < x < a + j h , j=1,2,…,ns ,$
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 to the histogram and b = nsh$b={n}_{s}h$ is the upper bound. The value h$h$ is known as the window width. To produce a smoother density estimate a kernel method can be used. A kernel function, K(t)$K\left(t\right)$, satisfies the conditions:
 ∞ ∫ K(t)dt = 1  and  K(t) ≥ 0. − ∞
$∫-∞∞K(t)dt=1 and K(t)≥0.$
The kernel density estimator is then 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$ is usually not important but to ease the computational burden use can be made of the Gaussian kernel defined as
 K(t) = 1/(sqrt(2π))e − t2 / 2. $K(t)=12πe-t2/2.$
The smoothness of the estimator depends on the window width h$h$. The larger the value of h$h$ the smoother the density estimate. The value of h$h$ can be chosen by examining plots of the smoothed density for different values of h$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$a$ to b$b$ the following steps are required.
 (i) Discretize the data to give ns${n}_{s}$ equally spaced points tl${t}_{l}$ with weights ξl${\xi }_{l}$ (see Jones and Lotwick (1984)). (ii) Compute the fft of the weights ξl${\xi }_{l}$ to give Yl${Y}_{l}$. (iii) Compute ζl = e − (1/2)h2sl2Yl${\zeta }_{l}={e}^{-\frac{1}{2}{h}^{2}{s}_{l}^{2}}{Y}_{l}$ where sl = 2πl / (b − a)${s}_{l}=2\pi l/\left(b-a\right)$. (iv) Find the inverse fft of ζl${\zeta }_{l}$ to give f̂(x)$\stackrel{^}{f}\left(x\right)$.
To compute the kernel density estimate for further values of h$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 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:     x(n) – double array
n, the dimension of the array, must satisfy the constraint n > 0${\mathbf{n}}>0$.
The n$n$ observations, xi${x}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
2:     window – double scalar
h$h$, the window width.
Constraint: window > 0.0${\mathbf{window}}>0.0$.
3:     slo – double scalar
a$a$, the lower limit of the interval on which the estimate is calculated. For most applications slo should be at least three window widths below the lowest data point.
Constraint: ${\mathbf{slo}}<{\mathbf{shi}}$.
4:     shi – double scalar
b$b$, the upper limit of the interval on which the estimate is calculated. For most applications shi should be at least three window widths above the highest data point.
5:     usefft – logical scalar
Must be set to false if the values of Yl${Y}_{l}$ are to be calculated by nag_smooth_kerndens_gauss (g10ba) and to true if they have been computed by a previous call to nag_smooth_kerndens_gauss (g10ba) and are provided in fft. If usefft = true${\mathbf{usefft}}=\mathbf{true}$ then the arguments n, slo, shi, ns and fft must remain unchanged from the previous call to nag_smooth_kerndens_gauss (g10ba) with usefft = false${\mathbf{usefft}}=\mathbf{false}$.
6:     fft(ns) – double array
ns, the dimension of the array, must satisfy the constraint
• ns2${\mathbf{ns}}\ge 2$
• The largest prime factor of ns must not exceed 19$19$, and the total number of prime factors of ns, counting repetitions, must not exceed 20$20$
• .
If usefft = true${\mathbf{usefft}}=\mathbf{true}$, fft must contain the fast Fourier transform of the weights of the discretized data, ξl${\xi }_{\mathit{l}}$, for l = 1,2,,ns$\mathit{l}=1,2,\dots ,{n}_{s}$. Otherwise fft need not be set.

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the array x.
n$n$, the number of observations in the sample.
Constraint: n > 0${\mathbf{n}}>0$.
2:     ns – int64int32nag_int scalar
Default: The dimension of the array fft.
The number of points at which the estimate is calculated, ns${n}_{s}$.
Constraints:
• ns2${\mathbf{ns}}\ge 2$;
• The largest prime factor of ns must not exceed 19$19$, and the total number of prime factors of ns, counting repetitions, must not exceed 20$20$.

None.

### Output Parameters

1:     smooth(ns) – double array
The ns${n}_{s}$ values of the density estimate, (tl)$\stackrel{^}{f}\left({t}_{\mathit{l}}\right)$, for l = 1,2,,ns$\mathit{l}=1,2,\dots ,{n}_{s}$.
2:     t(ns) – double array
The points at which the estimate is calculated, tl${t}_{\mathit{l}}$, for l = 1,2,,ns$\mathit{l}=1,2,\dots ,{n}_{s}$.
3:     fft(ns) – double array
The fast Fourier transform of the weights of the discretized data, ξl${\xi }_{\mathit{l}}$, for l = 1,2,,ns$\mathit{l}=1,2,\dots ,{n}_{s}$.
4:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and 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.

ifail = 1${\mathbf{ifail}}=1$
 On entry, n ≤ 0${\mathbf{n}}\le 0$, or ns < 2${\mathbf{ns}}<2$, or ${\mathbf{shi}}\le {\mathbf{slo}}$, or window ≤ 0.0${\mathbf{window}}\le 0.0$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, nag_smooth_kerndens_gauss (g10ba) has been called with usefft = true${\mathbf{usefft}}=\mathbf{true}$ but the function has not been called previously with usefft = false${\mathbf{usefft}}=\mathbf{false}$, or nag_smooth_kerndens_gauss (g10ba) has been called with usefft = true${\mathbf{usefft}}=\mathbf{true}$ but some of the arguments n, slo, shi, ns have been changed since the previous call to nag_smooth_kerndens_gauss (g10ba) with usefft = false${\mathbf{usefft}}=\mathbf{false}$.
ifail = 3${\mathbf{ifail}}=3$
On entry, at least one prime factor of ns is greater than 19$19$ or ns has more than 20$20$ prime factors.
W ifail = 4${\mathbf{ifail}}=4$
On entry, the interval given by slo to shi does not extend beyond three window widths at either extreme of the dataset. This may distort the density estimate in some cases.

## 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$n$, while the time for computing the fft is of order nslog(ns)${n}_{s}\mathrm{log}\left({n}_{s}\right)$, as is the time for computing the inverse of the fft.

## Example

function nag_smooth_kerndens_gauss_example
% Initialize the seed
seed = [int64(6698)];
% genid and subid identify the base generator
genid = int64(1);
subid =  int64(1);

window = 0.1;
slo = -4;
shi = 4;
usefft = false;
fft = zeros(100,1);

% Initialize the generator to a repeatable sequence
[state, ifail] = nag_rand_init_repeat(genid, subid, seed);
[state, x, ifail] = nag_rand_dist_normal(int64(1000), 0, 1, state);
[smooth, t, fftOut, ifail] = nag_smooth_kerndens_gauss(x, window, slo, shi, usefft, fft);
[ifail, y, isort] = nag_stat_plot_scatter_2var(t, smooth, int64(40), int64(20))

.+....+....+....+....+....+....+....+....+.
0.5000+                    +                    +
.                   1.                    .
.                   2.                    .
.                    1                    .
.                  1 12                   .
0.3750+                    + 11                 +
.                  1 .11                  .
.                21  .  1                 .
.                 2  .                    .
.               1    .  1                 .
0.2500+                    +   1                +
.              11    .   12               .
.             111    .    1               .
.                    .     1              .
.             1      .     121            .
0.1250+            11      +      1             +
.           11       .       12           .
.           2        .        1           .
.         12         .         2          .
.       122          .          3231      .
0.0000+13232322..+....+....+....+....+...1323231+
.+....+....+....+....+....+....+....+....+.
-4.000    -2.000     0.000     2.000     4.000
-3.000    -1.000     1.000     3.000

ifail =

0

y =

0.4646
0.4550
0.4459
0.4259
0.4123
0.4022
0.3982
0.3970
0.3742
0.3708
0.3587
0.3544
0.3455
0.3253
0.3250
0.3194
0.3135
0.3053
0.3021
0.2782
0.2770
0.2617
0.2303
0.2286
0.2165
0.2164
0.2135
0.2107
0.2060
0.2048
0.1958
0.1821
0.1590
0.1558
0.1452
0.1405
0.1387
0.1318
0.1198
0.1147
0.1062
0.1017
0.0974
0.0945
0.0899
0.0856
0.0843
0.0672
0.0608
0.0521
0.0388
0.0386
0.0385
0.0364
0.0355
0.0294
0.0282
0.0266
0.0242
0.0219
0.0216
0.0215
0.0210
0.0173
0.0161
0.0155
0.0150
0.0110
0.0104
0.0092
0.0072
0.0072
0.0067
0.0056
0.0044
0.0044
0.0042
0.0039
0.0034
0.0030
0.0030
0.0019
0.0016
0.0015
0.0015
0.0004
0.0004
0.0002
0.0002
0.0002
0.0002
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0
0
0

isort =

48
49
47
50
52
51
53
46
56
57
45
54
55
40
44
41
58
43
42
59
39
60
38
61
62
63
35
37
64
36
34
65
69
33
66
70
68
32
67
31
71
30
29
73
72
74
28
27
75
26
76
23
25
24
22
77
78
81
21
82
80
79
84
83
19
20
85
87
18
86
14
15
88
93
7
13
6
16
17
92
94
8
5
89
12
91
95
9
4
90
11
96
10
98
3
100
2
1
97
99

function g10ba_example
% Initialize the seed
seed = [int64(6698)];
% genid and subid identify the base generator
genid = int64(1);
subid =  int64(1);

window = 0.1;
slo = -4;
shi = 4;
usefft = false;
fft = zeros(100,1);

% Initialize the generator to a repeatable sequence
[state, ifail] = g05kf(genid, subid, seed);
[state, x, ifail] = g05sk(int64(1000), 0, 1, state);
[smooth, t, fftOut, ifail] = g10ba(x, window, slo, shi, usefft, fft);
[ifail, y, isort] = g01ag(t, smooth, int64(40), int64(20))

.+....+....+....+....+....+....+....+....+.
0.5000+                    +                    +
.                   1.                    .
.                   2.                    .
.                    1                    .
.                  1 12                   .
0.3750+                    + 11                 +
.                  1 .11                  .
.                21  .  1                 .
.                 2  .                    .
.               1    .  1                 .
0.2500+                    +   1                +
.              11    .   12               .
.             111    .    1               .
.                    .     1              .
.             1      .     121            .
0.1250+            11      +      1             +
.           11       .       12           .
.           2        .        1           .
.         12         .         2          .
.       122          .          3231      .
0.0000+13232322..+....+....+....+....+...1323231+
.+....+....+....+....+....+....+....+....+.
-4.000    -2.000     0.000     2.000     4.000
-3.000    -1.000     1.000     3.000

ifail =

0

y =

0.4646
0.4550
0.4459
0.4259
0.4123
0.4022
0.3982
0.3970
0.3742
0.3708
0.3587
0.3544
0.3455
0.3253
0.3250
0.3194
0.3135
0.3053
0.3021
0.2782
0.2770
0.2617
0.2303
0.2286
0.2165
0.2164
0.2135
0.2107
0.2060
0.2048
0.1958
0.1821
0.1590
0.1558
0.1452
0.1405
0.1387
0.1318
0.1198
0.1147
0.1062
0.1017
0.0974
0.0945
0.0899
0.0856
0.0843
0.0672
0.0608
0.0521
0.0388
0.0386
0.0385
0.0364
0.0355
0.0294
0.0282
0.0266
0.0242
0.0219
0.0216
0.0215
0.0210
0.0173
0.0161
0.0155
0.0150
0.0110
0.0104
0.0092
0.0072
0.0072
0.0067
0.0056
0.0044
0.0044
0.0042
0.0039
0.0034
0.0030
0.0030
0.0019
0.0016
0.0015
0.0015
0.0004
0.0004
0.0002
0.0002
0.0002
0.0002
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0
0
0

isort =

48
49
47
50
52
51
53
46
56
57
45
54
55
40
44
41
58
43
42
59
39
60
38
61
62
63
35
37
64
36
34
65
69
33
66
70
68
32
67
31
71
30
29
73
72
74
28
27
75
26
76
23
25
24
22
77
78
81
21
82
80
79
84
83
19
20
85
87
18
86
14
15
88
93
7
13
6
16
17
92
94
8
5
89
12
91
95
9
4
90
11
96
10
98
3
100
2
1
97
99