SJ bandwidth estimation

Introduction¶

Kernel Density Estimation (KDE) is a powerful non-parametric technique used to estimate the probability density function of a random variable without relying on rigid structural assumptions (such as the bell curve of a Normal distribution). The behavior of a KDE model depends primarily on two components: the kernel function-typically a standard Gaussian distribution that assigns local probability weights based on distance, and the bandwidth parameter $h$.

The bandwidth $h$ acts as a crucial smoothing parameter that dictates sensitivity of the model. If $h$ is too small, the estimator overfits the data, creating a noisy, wiggling curve. Conversely, if $h$ is too large, the estimator undersmooths (underfits) the data, washing away genuine structural details into a flat, uninformative curve Scott, 2012.

While numerous methodologies have been developed to optimize $h$, the plug-in approach by Sheather & Jones (1991) remains a benchmark standard, reliably minimizing the Asymptotic Mean Integrated Squared Error (AMISE) across a wide array of complex, non-parametric distributions Eidous et al., 2010.

Workflow¶

Initially, the input data vector $X$ is partitioned into discrete bins to calculate a bin-step density (detailed in Bin-step density). This binning strategy dramatically reduces the computational cost of evaluating the double summations required by the pilot bandwidth estimators when performing the secondary bandwidth $\alpha$ estimation detailed in the secondary bandwidth ($\alpha$) estimator.

From that, the roughness functionals are estimated, allowing the optimal primary bandwidth $h$ to be found (as described in the primary bandwidth $h$) by solving the derivative of the AMISE equation, the core metric evaluating the accuracy of estimation. The underlying mathematical theories are detailed in the boxes below.

Bin-step density¶

First, the data is partitioned into $N$ discrete bins (the default is 1000) to determine the individual bin frequencies ($b_k$, where $0 \le k \le N-1$). These frequencies are used to efficiently compute the total number of data point pairs at each specific bin-step.

Equation (10) measures total density of unique pairs ($c_i$) between 2 bins with the $i$-bin distance.

For the baseline case where the distance between bins is zero ($i = 0$), the value $c_0$ represents the pairs residing within the exact same bin, described in Equation (11).

Secondary bandwidth ($\alpha$) estimator¶

Equation (12) is transformed from equation (9), where both kernel function $K$ and $L$ follows Probability Density Function (PDF) of Normal Distribution. Additionally, $R(f'')$ and $R(f''')$ are estimated to $\hat{S}_D(a)$ and $\hat{T}_D(b)$ perspectively using the method of Jones & Sheather (1991) with the heuristic pilot bandwidth $a$ and $b$. The practical computation of $\hat{S}_D(a)$ and $\hat{T}_D(b)$ details in The first pilot bandwidth and The second pilot bandwidth correspondently.

The first pilot bandwidth¶

By Equation (8), assuming both the kernel function $L$ and the true density function $f$ follow normal distributions, the heuristic bandwidth $a$ for bias optimization is obtained as Equation (13) with $\hat{\lambda}$ is the estimated interquartile.

The pilot roughness estimator $\hat{S}_D(a)$ of Equation (12) is estimated following Equation (14), which originates from Equation (5). Practically, due to using bin-step density ($c_i$) with the bin width $d$ calculated while dividing bins, it reduces complexity from $O(k^2)$ ($k$ is the data size) to $O(N)$ comparing to the theoretical formula.

In Equation (14), $\phi^{\text{iv}}(t)$ is the 4-th derivative of the Normal Distribution PDF. The expansion using the Hermite polynomials is shown by Equation (15).

The second pilot bandwidth¶

Simmilarly, $R(f''')$ is estimated and canceled out bias by using specific bandwidth $b$. Subsequently, assuming both kernel functions follow normal scale model, the heuristic bandwidth $b$ is computed following Equation (16).

Likewise to the first pilot bandwidth, the component $\hat{T}_D(b)$ estimated for $R(f''')$ in Equation (12) is computed by Equation (17). Note that $\phi^{\text{vi}}(t)$ is the 6-th derivative of the PDF of Normal Distribution.

Primary bandwidth ($h$) computation¶

Ultimately, the optimal bandwidth $h$ is found by Equation (18), which originates from Equation (2) with $K$ following normal scale model and the substitution of $R(f'')$ by $\hat{S}_D( \hat{\alpha}(h) )$. Noticeably, data size $n$ is replaced by the number of bins $N$.

Finally, Equation (18) is solved using Brent’s method Brent, 2013 to find optimized bandwidth $h$.

Summary¶

Overall, the Sheather-Jones method optimizes the bandwidth parameter ($h$) by minimizing the AMISE. The primary challenge of this approach lies in estimating the unknown curvature of the true underlying density function $f$. This is achieved by utilizing a pilot estimation step (a secondary KDE) with an auxiliary bandwidth, which effectively cancels the asymptotic biases of the roughness estimator.

In practice, the implementations utilize a linear binning technique to drastically reduce the computational complexity. Furthermore, while Sheather & Jones (1991) originally outlined a Newton-Raphson framework to solve for the final bandwidth, the practical implementation employs Brent’s method. This numerical root-finding approach guarantees convergence and computational robustness when solving Equation (18).

However, the primary theoretical drawback of this method is its circularity: it relies on the assumption that the real unknown density function $f$ can be approximated well enough by an initial normal-scaled kernel estimator.