カーネル密度推定のbw.nrd0

• help(density)でdensity()のヘルプが出る
• 平滑化の幅のオプションが色々あるが、help(bw.nrd)でその色々のヘルプ記事が出る
• その中で「最適」とは言いにくいかもしれないが、歴史的経緯から良く用いられているnrd0オプションでは…などの説明がある

bw.nrd0 implements a rule-of-thumb for choosing the bandwidth of a Gaussian kernel density estimator. It defaults to 0.9 times the minimum of the standard deviation and the interquartile range divided by 1.34 times the sample size to the negative one-fifth power (= Silverman's ‘rule of thumb’, Silverman (1986, page 48, eqn (3.31)) unless the quartiles coincide when a positive result will be guaranteed.

• とあるが、0.9をかける相手、1.34で割る相手…など、文章で書かれるとわかりにくい→ので、式で
• もちろん、それぞれの関数を表示させてもよい。特に、bw.ucv(),bw.bcv(),bw.SJ()などは関数が長いことからもわかるように、単純な計算式とは違う(適切解の探索)

#nrd0
x <- precip
ssdd <- sd(x)
qq <- quantile(x,c(0.25,0.75))
min(ssdd,(qq[2]-qq[1]))*0.9 /(1.34*length(x)^(1/5))
bw.nrd0(x)
# nrd
min(ssdd,(qq[2]-qq[1]))*1.06 /(1.34*length(x)^(1/5))
bw.nrd(x)

plot(density(precip, n = 1000))
rug(precip)
lines(density(precip, bw = "nrd"), col = 2)
lines(density(precip, bw = "ucv"), col = 3)
lines(density(precip, bw = "bcv"), col = 4)
lines(density(precip, bw = "SJ-ste"), col = 5)
lines(density(precip, bw = "SJ-dpi"), col = 6)
legend(55, 0.035,
       legend = c("nrd0", "nrd", "ucv", "bcv", "SJ-ste", "SJ-dpi"),
       col = 1:6, lty = 1)