編輯部擬在微信端開啟《朝花夕拾》欄目，目的是推送2013年（含）之前主站發表的優秀文章，微信端與主站的同步始於2013年年初，然而初期用戶量有限故優質文章可能被埋沒。
那年當我看到「O(h) 項的偏差就神奇地消失了」這裡時，心裡不禁大讚，厲害了我的小軒哥！
一、一個例子核密度估計應該是大家常用的一種非參數密度估計方法，從某種程度上來說它的性質比直方圖更好，可以替代直方圖來展示數據的密度分布。但是相信大家會經常遇到一個問題，那就是有些數據是嚴格大於或等於零的，在這種情況下，零附近的密度估計往往會出現不理想的情況。下面以一個指數分布的模擬數據為例（樣本量為1000），R程序代碼為：
set.seed(123)
x <- rexp(1000, 1)
plot(density(x, kernel = "epanechnikov"), ylim = c(0, 1.2))
lines(seq(0, 8, by = 0.02), dexp(seq(0, 8, by = 0.02), 1), col = "blue")
abline(v = 0, col = "red")
可以看出，理論上應該單調遞減的密度函數在0附近有明顯的「陡坡」，而且不應該有密度的小於零的區域也有著正的估計值。當樣本量增大時，這種現象也不會得到明顯好轉，下圖是將樣本量改為10000時的情形。
set.seed(123)
x <- rexp(10000, 1)
plot(density(x, kernel = "epanechnikov"), ylim = c(0, 1.2))
lines(seq(0, 8, by = 0.02), dexp(seq(0, 8, by = 0.02), 1), col = "blue")
abline(v = 0, col = "red")
我們也許從平時看的書中了解到，當樣本量趨於無窮時，核密度估計值將是收斂到真實的密度函數的，但我們可能不會特意去研究那些結論成立的條件。以上這兩個簡單的例子似乎給了我們一個直觀的感覺，那就是當真實密度函數的支撐集（函數的支撐集指的是使得的x的集合）有邊界時，核密度估計值可能會出現一些不理想的情況。下面就簡單地給出一些理論的結果。
二、理論分析在一些必要的條件下（真實的密度函數 f 二階導絕對連續，三階導平方可積），核密度估計值的偏差有表達式，其中h是帶寬，，是支撐集為[-1,1]的核函數（即在[-1,1]上有值，其餘的地方取零）。可以看出這個偏差是隨著帶寬h的減小以的速度趨於零的。
而假設密度函數以0為邊界，那麼上述表達式將不再成立，而是代之以
其中。可以看出，當時，，，此時的偏差跟之前的那個表達式沒有區別；但當時，和都是非零的，於是偏差總是存在。
也許你會提議說，將估計值除以，偏差就可以減小了吧？的確，這樣是一種改進的辦法，但也要注意到，此時h的一次項不會消除，也就是說原來的衰減速度放慢到了h，從效率上說相對於理想的情況是大打了折扣。
這時候一個巧妙的辦法是，用另外一個核函數對f也做一次估計，那麼就有
其中的和意義類似，只不過是針對的。
對以上兩個式子進行線性組合，則會有
然後把的係數移到等式左邊，項的偏差就神奇地消失了。
通過觀察核密度估計的表達式，我們可以將上面這個過程等價地認為是對用了一個新的核函數進行估計，這個新的核函數是
特別地，如果將取為，那麼將有一個簡單的形式
當時，這個新的核函數就是，而當時（也就是在邊界），它會對最初的核函數進行調整。當時，不管算出來的估計值是多少，都只需將密度的估計值取為0即可。
三、程序實現下面這段程序是對本文的第一幅圖進行「整容」，代碼及效果圖如下：
k <- function(x) 3 / 4 * (1 - x^2) * (abs(x) <= 1)
a0 <- function(u, h) {
lb <- -1
ub <- pmin(u / h, 1)
0.75 * (ub - lb) - 0.25 * (ub^3 - lb^3)
}
a1 <- function(u, h) {
lb <- -1
ub <- pmin(u / h, 1)
3 / 8 * (ub^2 - lb^2) - 3 / 16 * (ub^4 - lb^4)
}
a2 <- function(u, h) {
lb <- -1
ub <- pmin(u / h, 1)
0.25 * (ub^3 - lb^3) - 0.15 * (ub^5 - lb^5)
}
kernel.new <- function(x, u, h) {
k(x) * (a2(u, h) - a1(u, h) * x) / (a0(u, h) * a2(u, h) - a1(u, h)^2)
}
den.est <- function(u, ui, h) {
sapply(u, function(u) ifelse(u < 0, 0, mean(kernel.new((u - ui) / h, u, h)) / h))
}
set.seed(123)
dat <- rexp(1000, 1)
x <- seq(0, 8, by = 0.02)
y <- den.est(x, dat, 2 * bw.nrd0(dat))
plot(x, y, type = "l", ylim = c(0, 1.2), main = "Corrected Kernel")
lines(x, dexp(x, 1), col = "red")修正核函數（樣本量1000）
從中可以看出，邊界的偏差問題已經得到了很好的改進。
如果真實的密度函數的支撐集不是，而是某一個閉區間，那麼偏差修正的過程與上面類似，只不過是要將定義為。在編程序的時候，也只需把積分的上下限進行相應的調整即可。
四、參考文獻Jeffrey S. Simonoff, 1998. Smoothing Methods in Statistics. Springer-Verlag：http://pages.stern.nyu.edu/~jsimonof/SmoothMeth/
邱怡軒，普渡大學統計系博士，現為卡耐基梅隆大學博士後，感興趣的方向包括計算統計學、機器學習、大型數據處理等，參與翻譯了《應用預測建模》《R 語言編程藝術》《ggplot2：數據分析與圖形藝術》等經典書籍，是 showtext、RSpectra、recosystem、prettydoc 等流行 R 軟體包的作者。
編輯：向悅
統計之都：專業、人本、正直的中國統計學社區。
關注方式：掃描下圖二維碼。或查找公眾帳號，搜索 統計之都 或 CapStat 即可。
往期推送：進入統計之都會話窗口，點擊右上角小人圖標，查看歷史消息即可。