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

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} ] 

enter image description here

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?