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.
KernelMixtureDistribution[]. $\endgroup$ – J. M.'s ennui♦ May 19 '13 at 1:57