The default bandwidth of the SmoothKernelDistribution function

I am researching the Silverman rule of thumb for bandwidth selection. Below is a sample of my code, from which it seems that BW=63 is close to the default parameter given by Mathematica (2 images are quite similar).

data4 = RandomReal[1024, {3600, 2}];
\[ScriptCapitalD]1 = SmoothKernelDistribution[data4];
\[ScriptCapitalD]2 = SmoothKernelDistribution[data4, 63];
{ContourPlot[PDF[\[ScriptCapitalD]1, {x, y}], {x, 0, 1024}, {y, 0, 1024}],
ContourPlot[PDF[\[ScriptCapitalD]2, {x, y}], {x, 0, 1024}, {y, 0, 1024}]}


Yet, according to articles the definition for Silverman R.O.T. should be:

H1 = (4/(dim + 2))^(1/(dim + 4))*n^(-1/(dim + 4)) StandardDeviation[data4[[All, 1]]]
(* -> 74.7864 *)


When I did this analysis for the 1D case, the Mathematica results fit well to expectations.

My question is: How does Mathematica calculate BW by default? How can I see the BW results (other than trial and error until images are alike)?

• Bill thank you. i probably missed explain my self. The default Mathematica "Silverman" is diffrent from the literature "Silverman". my question is how doeas mathematica calculates Silverman R.O.T? – Doron May 18 '13 at 7:17
• I figured that out once I realized Silverman was the default. What are you using for dim and n in your example? – bill s May 18 '13 at 7:21
• n=3600 and dim=2 – Doron May 18 '13 at 7:23
• H1 = n^(-1/(dim + 4)) StandardDeviation[Data[[All, 1]]] sorry for the mistake. Still it is qiute far from the expected Silverman BW=63 – Doron May 18 '13 at 7:48
• The documentation gives their particular version of the Silverman rule within the documentation for KernelMixtureDistribution[]. – J. M.'s ennui May 19 '13 at 1:57

I'd like to add to my previous comment that the undocumented "Bandwidth" property exists for SmoothKernelDistribution (SKD) and KernelMixtureDistribution (KMD).

Note that in the multivariate case there is a bit of inconsistency in how this is currently handled for SKD vs. KMD. SKD simplifies the bandwidth matrix when the off-diagonal elements are zeros. KMD gives the full bandwidth matrix. Note that the KernelMixtureDistribution docs give a good bit of detail on what this all means.

data4 = RandomReal[1024, {3600, 2}];

SmoothKernelDistribution[data4]["Bandwidth"]

(*{64.2393, 64.7462}*)

KernelMixtureDistribution[data4]["Bandwidth"]

(*{{64.2393, 0.}, {0., 64.7462}}*)


Now as for how things are computed. I always recommend using KernelMixtureDistribution with MaxMixtureKernels set to All this gives the precise definition of a kernel density estimator and can act on symbolic data. This will be very handy for you in determining how things are computed. For example:

FullSimplify[#["Bandwidth"], DistributionParameterAssumptions[#]] &
[KernelMixtureDistribution[Array[x, 4], "Silverman", MaxMixtureKernels -> All]]

(*(1/(20 2^(
2/5)))3 √(9 x[1]^2 + 9 x[2]^2 + 9 x[3]^2 - 6 x[3] x[4] +
9 x[4]^2 - 6 x[2] (x[3] + x[4]) - 6 x[1] (x[2] + x[3] + x[4]))*)


Note that SmoothKernelDistribution is effectively a linear interpolation of the PDF for KernelMixtureDistribution so all of the properties like bandwidth are the same.

EDIT:

I'd like to add two things. First, there is some confusion that has arisen because there are in fact two rules of thumb attributed to Silverman that differ by a constant (0.9 vs about 1.06). Mathematica uses a constant 0.9 for the default "Silverman" method. This is mainly because 1.06 assumes the underlying distribution is normal. For distributions with heavier tails or multiple modes Silverman claims that 1.06 will oversmooth and suggests 0.9 is a better choice. See B. W. Silverman, "Density Estimation for Statistics and Data Analysis", Chapman & Hall/CRC, 1986, pages 45–48 for details.

Since there is no prior knowledge on what the underlying distribution might be the more robust 0.9 constant was chosen as default.

To use the alternate Silverman rule of thumb (using the optimal constant for an underlying normal distribution) use "StandardGaussian".

Note that for the Buffalo snowfall data "StandardGaussian" achieves precicely the same value computed by wolfies.

KernelMixtureDistribution[
ExampleData[{"Statistics", "BuffaloSnow"}], "StandardGaussian"]["Bandwidth"]

(* 10.9706 *)


Second, for higher dimensions the Silverman rule is given similar to what the OP posts but with an additional constant $\displaystyle\frac{9 \times 3^{1/5}}{10\times 2^{2/5}}$. The following function should give the same result as the "Silverman" method for any dimension.

bwSilverman[data_] := Block[{dim, n, c0, c1, c2, sd, iqr},
dim = If[MatrixQ[data], Length[First[data]], 1];
n = Length[data];
c0 = (9 3^(1/5))/(10 2^(2/5));
c1 = (4/(dim + 2))^(1/(dim + 4));
c2 = n^(-1/(dim + 4));
sd = StandardDeviation[data];
iqr = (Quantile[data, .75] - Quantile[data, .25])/1.34;

If[dim == 1,
c0 c1 c2 Min[sd, iqr],
MapThread[c0 c1 c2 Min[#1, #2] &, {sd, iqr}]
]

]

bwSilverman[data4]

(*{64.2393, 64.7462}*)

• As an addition to Andy's notes, you can use "Properties" to list all the possible properties of a DataDistribution[] object. – J. M.'s ennui May 19 '13 at 3:03
• Andy wrote: Silverman ... differ by a constant (0.9 vs about 1.06). Well, in the 1.06 version ... the 1.06 is not a constant ... it only happens to be 1.06 in the case of a Gaussian kernel. For any other kernel, the 1.06 changes depending on the kernel being used. Are you saying that Mathematica uses a constant 0.9 * blah for ALL kernels, when Silverman method is selected? Surely not? – wolfies May 20 '13 at 17:40
• @wolfies it differs by kernel function and is hard coded for each. For user defined kernels there is a fixed constant .9. It is obviously possible to derive the optimal value if the underlying distribution is known but this seemed like more work than it was worth at the time. If it is critical functionality I recommend filing a suggestion with WRI. – Andy Ross May 20 '13 at 17:49
• @wolfies KernelMixtureDistribution is meant to be exact if MaxMixtureKernels is set to All. – Andy Ross May 20 '13 at 18:41
• @J.M. not really. It seems to have been an observation from experience rather than a mathematical derivation. We should all keep in mind that bandwidth selection is much more of an art than a science and few are better smoothing "artists" than Silverman. – Andy Ross May 22 '13 at 17:59

The original poster ponders why, in a bivariate model, the SmoothKernelDistribution function ... using the Silverman 'rule of thumb' estimate of bandwidth ... yields results that do not match the original poster's calculations.

I had the same issue in a univariate form.

Short version It turns out that there are two versions of the Silverman rule of thumb in common usage, at least in the univariate case.

Long version Here is Parzen's (1979) yearly 'Snowfall in Buffalo' data (63 data points collected from 1910 to 1972, and measured in inches):

data = ExampleData[{"Statistics", "BuffaloSnow"}];


According to the help file for SmoothKernelDistribution, Mathematica states that, by default:

• For bandwidth selection: 'By default the "Silverman" method is used'
• The Gaussian kernel is used; i.e., $$\frac1{\sqrt{2\pi}}\exp\left(-\frac{x^2}{2}\right),\qquad -\infty<x<\infty$$

gauss = Exp[-(x^2/2)]/Sqrt[2 Pi];
domain[gauss] = {x, -Infinity, Infinity};


Under these assumptions, (see, for instance, even the simple wiki page), the Silverman rule of thumb method reduces to:

((4*StandardDeviation[data]^5)/(3*Length[data]))^(1/5)


10.9706

This is exactly the same result we get with the Bandwidth function in the mathStatica add-on to Mathematica, specifying the Silverman method and a Gaussian kernel:

Bandwidth[data, gauss, Method -> Silverman]


10.9706

[ It is worth noting here that the mathStatica answer is not a pre-prepared/caged answer for the Gaussian case ... rather, as part of its general design, mathStatica calculates from scratch the bandwidth that minimises asymptotic mean integrated square error, according to Silverman's rule of thumb, using the user's specified kernel. ]

Mathematica 9 returns something different:

SmoothKernelDistribution[data]["Bandwidth"]


9.32148

Changing from a Gaussian kernel to another kernel doesn't seem to help resolve the discrepancy.

So what is Mathematica doing? In the Gaussian case, the theoretical Silverman rule of thumb calculation is (for our data set with 63 observations):

2*(1/24)^(1/5) *
Min[StandardDeviation[data],
(Quantile[data, .75] - Quantile[data, .25])/1.34] n^(-1/5) /. n -> 63


10.9706

where:

kk = 2*(1/24)^(1/5) // N


1.05922

If one does a quick web search, you can find this constant multiplier 1.05922 given in dozens of papers on Silverman's rule-of-thumb for the Gaussian case.

By contrast, Mathematica seems to use this result (with 0.9 instead of 1.05922 etc) ... [ kindly confirmed by J.M. - see comment below ]

0.9 *
Min[StandardDeviation[data],
(Quantile[data, .75] - Quantile[data, .25])/1.34] n^(-1/5) /. n -> 63


9.32148

which is, according to this CrossValidated question, the formula given in the STATA user manual entry for the STATA function kdensity. It's a very blunt kind of 0.9; not 0.91121212... , just 0.9.

Update

Based on the R manual, it seems there are, in fact, two versions of Silverman's rule of thumb in implementation ... an approximate value 0.9 ... and the actual theoretical value 1.05922 .... just to keep everyone on their toes :)

There is a very nice theoretical derivation of the exact Silverman 1.05922 ROT result here (see page 11 and top page 12), by Hansen. I still don't know the theoretical basis for the 0.9 usage?? ... will do some library hunting.

• The version of Silverman's rule used, as I mentioned in the comments, is in the docs for KernelMixtureDistribution[]: silvermanBandwidth[data_] := Module[{n = Length[data], iqr = Quantile[data, .75] - Quantile[data, .25], sd = StandardDeviation[data]}, .9 Min[iqr/1.34, sd] n^(-1/5)]. Note also that Quantile[data, .75] - Quantile[data, .25] is quite different from InterquartileRange[data], due to different defaults of Quartile[] and Quantile[] (also discussed here). – J. M.'s ennui May 19 '13 at 19:41
• The quantile calculation is not the point here. The point is: where does the 0.9 come from? – wolfies May 19 '13 at 19:42
• I'm correcting your formula after "Mathematica seems to use this result...", that's all. I don't know why that constant was chosen. – J. M.'s ennui May 19 '13 at 19:46
• No problem - changed from InterquartileRange to Quantile[data,.75] - Quantile[data,.25]. All outputs remain the same either way. – wolfies May 19 '13 at 19:52
• @J.M. The rule of thumb used comes from page 48 of Silverman's own book "Density Estimation for Statistics and Data Analysis" and is suggested as an improvement to 1.06 when considering a broader family of underlying distributions. – Andy Ross May 19 '13 at 22:16

Finally, I think I found it--or, at least, I found a formula "close enough" to the way Mathematica calculates the BW by default:

hi = ((4*Subscript[σ, i]^5)/(3*n))^(1/5)


where σ is the standard deviation of the data. This should be performed to each component (e.g. i, j, k) of the dimensions.

The reference came from the folowing link:

http://sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore/ebooks/html/spm/spmhtmlnode15.html

I still need to figure out: is there a way to actualy see what Mathematica has calculated for the BW?

• I had classes with Prof. Härdle... his e-book is really amazing... – Rod May 18 '13 at 14:49
• can you please put a link to these books? – Doron May 18 '13 at 15:23
• You can get the bandwidth that was used using the undocumented property "Bandwidth" via SmoothKernelDistribution["Bandwidth"] – Andy Ross May 18 '13 at 15:44
• I'm not sure if his e-books are still available... but the R- and Matlab-codes use in e-books and classes are available at this link. – Rod May 18 '13 at 16:17
• Härdle's notes for the xplore package are live ... interestingly, he also uses the 1.06 result (rather than the 0.9 version) ... sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore/ebooks/html/spm/… – wolfies May 19 '13 at 23:24