hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
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 nn observations, x1,x2,,xnx1,x2,,xn, from a distribution with unknown density function, f(x)f(x), an estimate of the density function, (x)f^(x), may be required. The simplest form of density estimator is the histogram. This may be defined by:
(x) = 1/(nh) nj ,   a + (j1) 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 njnj is the number of observations falling in the interval a + (j1)ha+(j-1)h to a + jha+jh, aa is the lower bound to the histogram and b = nshb=nsh is the upper bound. The value hh is known as the window width. To produce a smoother density estimate a kernel method can be used. A kernel function, K(t)K(t), 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
(x) = 1/(nh)K((xxi)/h).
i = 1
f^(x)=1nh i= 1nK (x-xih) .
The choice of KK 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π))et2 / 2.
K(t)=12πe-t2/2.
The smoothness of the estimator depends on the window width hh. The larger the value of hh the smoother the density estimate. The value of hh can be chosen by examining plots of the smoothed density for different values of hh 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 aa to bb the following steps are required.
(i) Discretize the data to give nsns equally spaced points tltl with weights ξlξl (see Jones and Lotwick (1984)).
(ii) Compute the fft of the weights ξlξl to give YlYl.
(iii) Compute ζl = e(1/2)h2sl2Ylζl=e-12h2sl2Yl where sl = 2πl / (ba)sl=2πl/(b-a).
(iv) Find the inverse fft of ζlζl to give (x)f^(x).
To compute the kernel density estimate for further values of hh 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 > 0n>0.
The nn observations, xixi, for i = 1,2,,ni=1,2,,n.
2:     window – double scalar
hh, the window width.
Constraint: window > 0.0window>0.0.
3:     slo – double scalar
aa, 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: slo < shislo<shi.
4:     shi – double scalar
bb, 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 YlYl 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 = trueusefft=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 = falseusefft=false.
6:     fft(ns) – double array
ns, the dimension of the array, must satisfy the constraint
  • ns2ns2
  • The largest prime factor of ns must not exceed 1919, and the total number of prime factors of ns, counting repetitions, must not exceed 2020
  • .
    If usefft = trueusefft=true, fft must contain the fast Fourier transform of the weights of the discretized data, ξlξl, for l = 1,2,,nsl=1,2,,ns. Otherwise fft need not be set.

    Optional Input Parameters

    1:     n – int64int32nag_int scalar
    Default: The dimension of the array x.
    nn, the number of observations in the sample.
    Constraint: n > 0n>0.
    2:     ns – int64int32nag_int scalar
    Default: The dimension of the array fft.
    The number of points at which the estimate is calculated, nsns.
    Constraints:
    • ns2ns2;
    • The largest prime factor of ns must not exceed 1919, and the total number of prime factors of ns, counting repetitions, must not exceed 2020.

    Input Parameters Omitted from the MATLAB Interface

    None.

    Output Parameters

    1:     smooth(ns) – double array
    The nsns values of the density estimate, (tl)f^(tl), for l = 1,2,,nsl=1,2,,ns.
    2:     t(ns) – double array
    The points at which the estimate is calculated, tltl, for l = 1,2,,nsl=1,2,,ns.
    3:     fft(ns) – double array
    The fast Fourier transform of the weights of the discretized data, ξlξl, for l = 1,2,,nsl=1,2,,ns.
    4:     ifail – int64int32nag_int scalar
    ifail = 0ifail=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 = 1ifail=1
    On entry,n0n0,
    orns < 2ns<2,
    orshisloshislo,
    orwindow0.0window0.0.
      ifail = 2ifail=2
    On entry,nag_smooth_kerndens_gauss (g10ba) has been called with usefft = trueusefft=true but the function has not been called previously with usefft = falseusefft=false,
    ornag_smooth_kerndens_gauss (g10ba) has been called with usefft = trueusefft=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 = falseusefft=false.
      ifail = 3ifail=3
    On entry, at least one prime factor of ns is greater than 1919 or ns has more than 2020 prime factors.
    W ifail = 4ifail=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.

    Further Comments

    The time for computing the weights of the discretized data is of order nn, while the time for computing the fft is of order nslog(ns)nslog(ns), 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
    
    
    

    PDF version (NAG web site, 64-bit version, 64-bit version)
    Chapter Contents
    Chapter Introduction
    NAG Toolbox

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