Computes Gaussian kernel estimates of a univariate density via the fast Fourier transform over a fixed interval.

Routine DNFFT computes Gaussian kernel estimates of the density function for a random sample of (scalar‑valued) observations using a Gaussian kernel (normal density). The computations are comparatively fast because they are performed through the use of the fast Fourier transform. Routine DESKN should be used in place of DNFFT if a kernel other than the Gaussian kernel is to be used, if a irregular grid is desired, or if the approximations in DNFFT are not acceptable. Because of its speed, DNFFT will usually be preferred to DESKN.

A Gaussian kernel estimate of the density at the point y is given by:

where

is the estimated density at y, xi denotes the i‑th observation, n is the number of observations, and h is a fixed constant called the window width supplied by the user. If density estimates for several different window sizes are to be computed, then DNFFT performs a fast Fourier transform on the data only during the first call (when IFFT is zero). On subsequent calls (with IFFT set at 1), the Fourier transform of the data need not be recomputed.

If the same value of NXPT is to be used with several different input vectors X, then the computations can be made faster by the use of D2FFT. In D2FFT, it is assumed that some constants required by the Fourier transform and its inverse have already been computed via routine FFTRI (IMSL MATH/LIBRARY) in work array WFFTR. If D2FFT is called repeatedly with the same value of NXPT, WFTTR need only be computed once.

Routine DNFFT is an implementation of Applied Statistics algorithm AS 176 (Silverman 1982) using IMSL routines for the fast Fourier transforms. Modification to algorithm AS 176, as discussed in Silverman (1986, pages 61–66), gives the details of the computational method. The basic idea is to partition the support of the density into NXPT equally‑sized nonoverlapping intervals. The frequency of the observations within each interval is then computed, and the Fourier transform of the frequencies obtained. Since the kernel density estimate is the convolution of the frequencies with the Gaussian kernel (for given window size), the Fourier coefficients for the Gaussian kernel density estimates are computed as the product of the coefficients obtained for the frequencies, times the Fourier coefficients for the Gaussian kernel function. The discrete Fourier coefficients for the Gaussian kernel may be estimated from the continuous transform. The inverse transform is then used to to obtain the density estimates.

Because the fast Fourier transform is used in computing

the computations are relatively fast (providing that NXPT is a product of small primes). To maintain precision, a large number of intervals, say 256, is usually recommended. Tapia and Thompson (1978), Chapter 2, give some considerations relevant to the choice of the window size parameter WINDOW. Generally, several different window sizes should be tried in order to obtain the best value for this parameter.

In this example, the density function is estimated at 64 points using a random sample of 150 points from a standard normal distribution. The actual density for the standard normal density is also reported in the output for comparison. The random sample is generated using routines RNSET and RNNOR in Chapter 18, “Random Number Generation”.