About Silverman's bandwidth selection in SmoothKernelDistribution

Question

I asked the following question in the site Cross Validated: https://stats.stackexchange.com/questions/444320/convergence-of-kernel-density-estimate-as-the-sample-size-grows/444479#444479. It was about estimating the probability density function via a kernel-based method.

For example, let $X\sim\text{Normal}(0,1)$ and let $f_X$ be its probability density function. Let $\hat{f}_X^M$ the kernel density estimate using $M$ realizations. Let $$\epsilon=E\left[\|f_X-\hat{f}_X^M\|_\infty\right]$$ be the error, where $E$ is the expectation. The procedure for estimating $\epsilon$ is as follows:

1. Consider a geometric sequence of lengths $M$ (in log-log scale, these points will be at the same distance).

2. For each $M$, generate $M$ realizations for constructing $\hat{f}_X^M$.

3. Compute $\|f_X-\hat{f}_X^M\|_\infty$ for all lengths $M$ (the infinity norm is approximated by discretizing the domain).

4. Repeat steps 1-3 several times, say 10 times, and compute the average to estimate $\epsilon$. Compute sample quantiles to derive an estimated $90\%$ confidence interval for $\|f_X-\hat{f}_X^M\|_\infty$.

Using the built-in function SmoothKernelDistribution with no specified option for step 2, I observed that the error $\epsilon$ goes down with $M$ until a certain length $M_0$, from which $\epsilon$ stabilizes:

It seems that the bandwidth selected is not sufficiently small, which produces bias (stabilization of the error) and small variance (more narrow confidence intervals) when $M$ grows. According to the documentation, SmoothKernelDistribution employs Silverman's rule with Gaussian kernel by default.

I implemented Silverman's rule myself with bandwidth choice $1.06\times\hat{\sigma}_X\times M^{-1/5}$, see the answer from Cross Validated. In such case, the error tends to 0 as $M$ grows with no problem:

If the bandwidth selection is $0.9\times\min\{\hat{\sigma}_X,\frac{\text{IQR}}{1.34}\}\times M^{-1/5}$, the error also tends to 0 with $M$:

Notice also that, in these last two cases, the error attained is smaller and the confidence interval has a similar width for every $M$.

I cannot post the code because it is part of a work. I just want to know which is the bandwidth selected by SmoothKernelDistribution (to know if I can use this function in the future), or if you have experienced any problems with this function as the sample size grows.

Answer 1

The Silverman rule is $0.9 \min({\hat{\sigma},IQR/1.34})n^{-1/5}$ and the rule $1.06\hat{\sigma} n^{-1/5}$ is the rule that minimizes the mean integrated square error when the underlying distribution is normal.

The default for SmoothKernelDistribution is to use the Silverman rule. That bandwidth (in Mathematica code) is

0.9*Min[StandardDeviation[x], InterquartileRange[x, {{0, 0}, {1, 0}}]/1.34]*n^(-1/5) 

However, it appears that for large sample sizes that rule is not used. What it uses for sample sizes somewhat larger than 100,000 is $0.09 \hat{\sigma}$ which remains relatively constant. Here is a check on what bandwidth SmoothKernelDistribution actually uses for various sample sizes:

SeedRandom[12345]
Do[x = RandomVariate[NormalDistribution[0, 1], n];
 (* Expected bandwidth *)
 bw = 0.9*Min[StandardDeviation[x], InterquartileRange[x, {{0, 0}, {1, 0}}]/1.34]*n^(-1/5);
 (* Bandwidth for large sample sizes *)
 bwSD = 0.09 StandardDeviation[x];
 Print[{n, SmoothKernelDistribution[x][[2, 3]], 
   SmoothKernelDistribution[x, bw][[2, 3]],
   bw, bwSD}],
 {n, {50, 100, 500, 5000, 10000, 50000, 100000, 110000, 200000, 1000000, 2000000}}]

(* {50,0.377071,0.377071,0.377071,0.0824551}
{100,0.36585,0.36585,0.36585,0.0918974}
{500,0.251404,0.251404,0.251404,0.0888302}
{5000,0.163237,0.163237,0.163237,0.0896631}
{10000,0.143387,0.143387,0.143387,0.0904713}
{50000,0.103388,0.103388,0.103388,0.0900045}
{100000,0.0895679,0.0895679,0.0895679,0.0895679}
{110000,0.0903596,0.0886534,0.0886534,0.0903596}
{200000,0.0901368,0.0784686,0.0784686,0.0901368}
{1000000,0.0899489,0.0567539,0.0567539,0.0899489}
{2000000,0.0900067,0.0494389,0.0494389,0.0900067} *)

So after the sample size reaches 100,000 the bandwidth appears to level off and SmoothKernelDistribution uses $0.09*\hat{\sigma}$. For what you're doing it would seem that you'd want to set the bandwidth explicitly with

bw = 0.9*Min[StandardDeviation[x], InterquartileRange[x, {{0, 0}, {1, 0}}]/1.34]*n^(-1/5)
SmoothKernelDistribution[x, bw]

You might want to report this issue to Wolfram, Inc.

****Note that the documentation for version 12.0 states that the default option for InterquartileRange is {{0, 0}, {1, 0}} but it appears to be {{1/2, 0}, {0, 1}}. I'll open a ticket on that. (This has to do with the zillion ways one can estimate a quantile from a finite sample.)