Mathematica has builtin bin 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 bin-width 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 bin 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 bin matrix, either for Silverman or any other rule?