2D kernel density estimation (SmoothKernelDistribution) with bandwidth estimation: what are the values that Mathematica chooses?

Mathematica has built in bandwidth estimation including the rules Scott, SheatherJones and Silverman (the default one); they work in both 1D and multiple dimensions. Most of the statistical documentation that I could find of these bandwidth rules are for 1D data. Their implementation for 2D or higher dimensions seems not, as far as I know, so robust.

I could not find a Mathematica documentation on how exactly these rules are implemented in any dimensions. For the Silverman case, there is a nice question about it that raises very important subtleties: About Silverman's bandwidth selection in SmoothKernelDistribution .

For 2D data, my first guess was that Mathematica uses the same 1D algorithm, but for each of the axis, thus yielding a diagonal bandwidth matrix. Hence, I extended the code provided in the previous link to 2D as follows:

Clear[data, silvermanBandwidth];
silvermanBandwidth[data_] := silvermanBandwidth[data] = Block[
{m, n},
m = MapThread[Min @ {#1, #2} &,
{
StandardDeviation @ data,
InterquartileRange[data, {{0, 0}, {1, 0}}]/1.349
}
];
n = Length @ data;
0.9 m/n^(1/5)
];

(In the statistical literature I could find different conventions for rounding the real numbers that appear in the above code, I do not know precisely which version Mathematica picks; anyway the problem below is larger than these small rounding changes).

The approach above (and a few variations I tried) is quite close to what Mathematica does in 2D, but it is not identical. Here is an example:

data = RandomReal[1, {100, 2}];
silvermanWMDist = SmoothKernelDistribution @ data;
silvermanMyDist = SmoothKernelDistribution[data, silvermanBandwidth @ data, "Gaussian"];
ContourPlot[PDF[silvermanWMDist, {x, y}],
{x, -0.1, 1.1},
{y, -0.1, 1.1}
]
ContourPlot[PDF[silvermanMyDist, {x, y}],
{x, -0.1, 1.1},
{y, -0.1, 1.1}
] My questions are: how Silverman's rule is implemented in Mathematica for 2D data? Is there a way to print out Mathematica's derived bandwidth matrix, either for Silverman or any other rule?

• I just learned that in the middle of the documentation for KernelMixtureDistribution reference.wolfram.com/language/ref/… there are definitions for some of the bandwidth rules, including the Silverman one. However, I can only see the 1D expressions.
– Davi
Sep 5 '21 at 23:53
• You can see the resulting bandwidths with SilvermanWMDist[[2,3]] and SilvermanMyDist[[2,3]]. The default does not appear to apply the Silverman rule that I know. It's sure looks closer to using n^(1/5.735) rather than n^(1/5). Your use of SilvermanMyDist matches what I learned for the Silverman rule (other than having a few more decimal places for the 0.9 constant). Asking Wolfram, Inc. directly might be the only way to know exactly what is being used for the defaults. However, I'd recommend calculating your own bandwidth that depends on the shape of the distribution.
– JimB
Sep 6 '21 at 4:37
• Thanks @JimB, this is an important part of the question. Very nice to know that it is easy to discover the derived bandwidth values. It is also good to know that my approach to 2D Silverman is reasonable. I will contact Wolfram on this. Perhaps there is some nice reason about mathematica's implementation that I am missing. Would you like to answer the question?
– Davi
Sep 6 '21 at 22:35
• Hopefully someone else knows the answer. Because the differences/ratios between the default and what we think is the Silverman approach are so consistent, I might stumble across what the default actually does. I'll spend a little time later today to see if any nice rule becomes apparent to me.
– JimB
Sep 6 '21 at 23:00
• Have you already seen this question? Sep 16 '21 at 13:40

Just to make it explicit for this question based on the link provided by @J.M.'storpor...

The two bandwidths chosen using what Mathematica calls the "Silverman" rule for 2D data is as follows:

n = 100;
SeedRandom;
xy = RandomVariate[MultinormalDistribution[{0, 0}, {{1, 1/2}, {1/2, 1}}], n];
skd = SmoothKernelDistribution[xy, "Silverman"];

bw = (9 3^(1/5))/(10 2^(2/5) n^(1/6)) Min[StandardDeviation[xy[[All, #]]],
(Quantile[xy[[All, #]], 0.75] - Quantile[xy[[All, #]], 0.25])/1.34] & /@ {1, 2}
(* {0.382494, 0.331636} *)

Besides the multiplier of 0.9 (or almost equivalently as 9/10) there is another multiplier of $$\sqrt{\frac{3}{4}}=\frac{\sqrt{3}}{2^{2/5}}\approx 0.944088$$. But that confuses me more because the reciprocal of that number is $$1.05922$$ which is either close or exactly the original multiplier (which is usually rounded to 1.06). Coincidence?

As a check one can obtain the resulting bandwidth with

skd[[2, 3]]
(* {0.382494, 0.331636} *)
• Thanks. This answers my two questions. I don't know from where these exact factors come from. If someone knows... Meanwhile, I received the answer from Wolfram, and they pointed out the same StackExchange post that you link. I asked them to include such definitions in future documentations, and they took note of it.
– Davi
Sep 26 '21 at 20:00
• I think I understand these coefficients now. In Silverman's book, from eqs. 4.14-4.15 and Table 4.1, one finds the exact coefficients, except for the factor $(9 * 3^{1/5})/(10 * 2^{2/5})$. The answer without the latter factor is useful when the underlying distribution is known to be Normal (which leads to the 1.06 coefficient in 1D problems, and which is not the coefficient used by Mathematica). The exact factor was introduced such that, when d=1, one recovers the 0.9 factor (actually 9/10), which is used when the distribution is unknown.
– Davi
Sep 26 '21 at 21:01