naginterfaces.library.smooth.kerndens_​gauss

naginterfaces.library.smooth.kerndens_gauss(x, comm, wtype=2, window=1.0, slo=None, shi=None, ns=512)

kerndens_gauss performs kernel density estimation using a Gaussian kernel.

For full information please refer to the NAG Library document for g10bb

https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/g10/g10bbf.html

Parameters

xfloat, array-like, shape $\left(n\right)$

$x_{\text{i}}$, for $\text{i}=1,2,\dots ,n$.

On a continuation call, $x$ must be unchanged since the last call to kerndens_gauss.

commdict, communication object, modified in place

Communication structure.

If not initialized on entry then the values of $Y_l$ are to be calculated by this call to kerndens_gauss.

Otherwise, this is a continuation call and it is assumed that the values of $Y_l$ were calculated by a previous call to this function and the relevant information is stored in $comm$[‘rcomm’].

wtypeint, optional

How the window width, $h$, is to be calculated:

$wtype=1$

$h$ is supplied in $window$.

$wtype=2$

$h$ is to be calculated from the data, with

$$h=m\times \left(\frac{0.9\times min(q_{75}-q_{25},\sigma )}{n^{0.2}}\right)$$

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 stat.quantiles. This is the ‘rule-of-thumb’ suggested by Silverman (1990).

windowfloat, optional

If $wtype=1$, then $h$, the window width. Otherwise, $m$, the multiplier used in the calculation of $h$.

sloNone or float, optional

Note: if this argument is None then a default value will be used, determined as follows: $3.0$.

If $slo<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:

$$\begin{array}{c}\begin{array}{ccc}a& =& {\text{min}}_i\left\{x_i\right\}-slo\times h\\ b& =& {\text{max}}_i\left\{x_i\right\}+slo\times h\end{array}\end{array}$$

where $h$ is the window width.

For most applications $a$ should be at least three window widths below the lowest data point.

On a continuation call, a supplied $slo$ will be ignored and the appropriate value from the previous call to kerndens_gauss will be extracted from $comm$[‘rcomm’].

shiNone or float, optional

Note: if this argument is None then a default value will be used, determined as follows: $0.0$.

If $slo<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.

On a continuation call, a supplied $shi$ will be ignored and the appropriate value from the previous call to kerndens_gauss will be extracted from $comm$[‘rcomm’].

nsint, optional

$n_s$, the number of points at which the estimate is calculated.

On a continuation call, $ns$ must be unchanged since the last call to kerndens_gauss.

Returns

windowfloat

$h$, the window width actually used.

slofloat

$a$, the lower limit actually used.

shifloat

$b$, the upper limit actually used.

smoothfloat, ndarray, shape $\left(ns\right)$

$\hat{f}\left(t_{\text{l}}\right)$, for $\text{l}=1,2,\dots ,n_s$, the $n_s$ values of the density estimate.

tfloat, ndarray, shape $\left(ns\right)$

$t_{\text{l}}$, for $\text{l}=1,2,\dots ,n_s$, the points at which the estimate is calculated.

Raises

NagValueError

(errno $11$)

On entry, $n=⟨value⟩$.

Constraint: $n>0$.

(errno $12$)

On entry, $n=⟨value⟩$.

On entry at previous call, $n=⟨value⟩$.

Constraint: on a continuation call, $\text{n}$ must be unchanged since previous call.

(errno $31$)

On entry, $wtype=⟨value⟩$.

Constraint: $wtype=1$ or $2$.

(errno $41$)

On entry, $window=⟨value⟩$.

Constraint: $window>0.0$.

(errno $51$)

On entry, $slo=⟨value⟩$.

On exit from previous call, $slo=⟨value⟩$.

Constraint: on a continuation call, $slo$ must be the same value returned by the previous call.

(errno $62$)

On entry, $shi=⟨value⟩$.

On exit from previous call, $shi=⟨value⟩$.

Constraint: on a continuation call, $shi$ must be the same value returned by the previous call.

(errno $71$)

On entry, $ns=⟨value⟩$.

Constraint: $ns\ge 2$.

(errno $74$)

On entry, $ns=⟨value⟩$.

On entry at previous call, $ns=⟨value⟩$.

Constraint: on a continuation call, $ns$ must be unchanged since previous call.

(errno $111$)

$comm$[‘rcomm’] has been corrupted between calls.

Warns

NagAlgorithmicWarning

(errno $61$)

On entry, $slo=⟨value⟩$ and $shi=⟨value⟩$.

On entry, $min\left(x\right)=⟨value⟩$ and $max\left(x\right)=⟨value⟩$.

Expected values of at least $⟨value⟩$ and $⟨value⟩$ for $slo$ and $shi$.

Notes

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, $\hat{f}\left(x\right)$, may be required. The simplest form of density estimator is the histogram. This may be defined by:

$$\hat{f}\left(x\right)=\frac{1}{nh}{n}_j\text{,}a+(j-1)h<x<a+jh\text{,}j=1,2,\dots ,n_s\text{,}$$

where $n_j$ is the number of observations falling in the interval $a+(j-1)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:

$${\int }_{-\infty }^{\infty }K\left(t\right)dt=1\text{and}K\left(t\right)\ge 0\text{.}$$

The kernel density estimator is then defined as

$$\hat{f}\left(x\right)=\frac{1}{nh}\sum _{i=1}^{n}K\left(\frac{x-x_i}{h}\right)\text{.}$$

The choice of $K$ is usually not important but to ease the computational burden use can be made of the Gaussian kernel defined as

$$K\left(t\right)=\frac{1}{\sqrt{2\pi }}{e}^{-t^2/2}\text{.}$$

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.

1. Discretize the data to give $n_s$ equally spaced points $t_l$ with weights ${\xi }_l$ (see Jones and Lotwick (1984)).

2. Compute the FFT of the weights ${\xi }_l$ to give $Y_l$.

3. Compute ${\zeta }_l=e^{-\frac{1}{2}{h}^2{s}_l^2}{Y}_l$ where $s_l=2\pi l/(b-a)$.

4. Find the inverse FFT of ${\zeta }_l$ to give $\hat{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