7
$\begingroup$

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:

enter image description here

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:

enter image description here

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ What do you mean by "not sufficiently small" ? If you are estimating non-normal distributions, then it usually isn't small enough but you don't mention that in your question. $\endgroup$ – JimB Jan 13 at 0:35
  • $\begingroup$ Showing your code (a minimal working example) would be helpful. $\endgroup$ – JimB Jan 13 at 0:42
  • $\begingroup$ The question has been edited with more information. $\endgroup$ – user69238 Jan 13 at 11:32
  • $\begingroup$ Thanks. I've edited my answer. Rather than dealing with the error ($\epsilon$) which you don't/can't show, it's more direct to look at the resulting bandwidth for large samples. You should report the issue to Wolfram, Inc. $\endgroup$ – JimB Jan 13 at 14:09
6
$\begingroup$

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.)

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Your answer is the key. I estimated the mean error and the confidence interval using as bandwidth $0.9\min\{\hat{\sigma}_X,\text{IQR}/1.34\}M^{-1/5}$ when $M\leq 100,000$, and $0.09\hat{\sigma}_X$ when $M>100,000$, and I obtained the corresponding plot to SmoothKernelDistribution (and also KernelMixtureDistribution) with no options: stabilization of the error and narrow confidence interval for large $M$. So the conclusion is that, for large sample size, the built-in functions fix the bandwidth as $0.09\hat{\sigma}_X$ and the estimates are not improved anymore. Thank you. $\endgroup$ – user69238 Jan 13 at 16:40
  • $\begingroup$ And I forgot to mention that the same issue occurs when the sample size is large for any other built-in option of kernel and bandwidth. $\endgroup$ – user69238 Jan 13 at 17:02
  • $\begingroup$ Thanks for accepting my answer. FYI: I've submitted a ticket to Wolfram, Inc. on this issue. $\endgroup$ – JimB Jan 13 at 17:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.