從一個連續型分布的群體隨機抽取 n 個觀測值 X1,...,Xn,我們想推論群體分布的模樣。在參數化模型,最常做的就是推論其未知參數;但群體模型若未參數化,一種方法是推論其分位數如中位數 (median)、四分位數 (quartile)、十分位數 (decile)、百分位數 (percentile) 等,另一種方法是用經驗分布函數 (empirical distribution function)推論群體的分布函數,再者是考慮機率密度函數 (probability density function)。

但 n 個觀測值本身是離散的,如何變成連續型的機率密度?直方圖 (histogram) 是類比機率密度的一個方法,單變量的情形是定 k 個切割點把資料分成 k+1 組,計算落入各組的觀測值個數 fj。此法在不等組距分組時必須以組距 (class interval) 除觀測值個數 fj 才能正確展現分布的模樣;若是等組距,則無需調整,fj 即能表現出機率密度的樣子。不過,直方圖只是圖示分布的模樣,若要將直方圖看成機率密度的估計,還要對直方圖高度單位適當設定,使直方圖上沿之下 x 軸之上的區域面積總和為 1,這是機率密度函數 f(x) 的基本特性:∫ f(x) dx = 1。

如果群體分布是很偏斜的 (skewed),或是長尾巴的 (long-tailed),做數值分組及畫直方圖常需要採不等組距。事實上我們也可採行另一辦法:把一個觀測值,以該觀測值為中心,左右各擴展半個組距,也就是觀測值 Xi 決定了一個均勻分布 (uniform distribution)

fi(x) = 1/h,   if  Xi-h/2 ≦ x ≦ Xi+h/2;   = 0  其他處

當然,必要時 h 也可以在不同 Xi 處取不同值。

想像樣本資料中有 fj 個落在第 j 組的組中點 (class mid-value),則第 j 組範圍的機率密度就是

(fj/h)/n = (1/n) Σ_{Xi in interval j} 1/h

如果樣本觀測值都是落在各組組中點,則

Σ_i fi(x)/n = fj/(nh), 當 x 在第 j 分組;  = 0  當 x 落在所有資料分組之外

這是直方圖的數學式,也是機率密度函數的一個估計。

與直方圖不同的是:我們實際抽取的觀測值 Xi 不會都恰好落在資料分組後各組的組中點,因此 Σ fi(x) 和直方圖有些差異。為了更清楚看出其差異,我們設 k+1 個等距分組的組中點是 Mj, j = 0,1,...,k, 由於是等距分組,因此 Mj = M0 + jh, h 是共同組距。令

K(x) = 1, -1/2 ≦ x ≦ 1/2;  = 0  elsewhere

則直方圖的機率密度式是

f`(x) = (1/n) Σ_i K((x-M(Xi))/h)/h

其中

M(x) = Mj  if  Mj-h/2 < x ≦ Mj+h/2

所以 M(Xi) 就是 Xi 所在分組的組中點;而新方法是

f*(x) = (1/n) Σ_i K((x-Xi)/h)/h

也就是:直方圖是以各組組中點為中心,或說是把 Xi 挪移到所在分組的組中點;而新方法直接以 Xi 為中心。我們稱上列 f*(x) 為群體機率密度函數 f(x) 的核估計,或說此法為核密度估計 (kernel density estimation, KDE),式中 K(x) 為核函數 (kernel function),h (或 h/2) 稱為帶寬 (band width)。

不只一種核函數,理論上任何連績型分布的機率密度函數都可以當作核函數,不過實際上常被使用的核函數不多,可參考 https://en.wikipedia.org/wiki/Kernel_(statistics)#In_non-parametric_statistics。如果資料沒有範圍限制,核函數可以進一步限制為對稱分布的機率密度函數如均勻、三角形、標準常態等對稱型分布;如果是有限制的,例如非負資料的核密度估計,如 gamma, log-normal 分布可能也合適,或者以 Xi 為左端點而非中心點的均勻分布;若資料限制在有限區間,無論何種核函數都會遭遇邊界問題:靠近端點的 Xi 可能使核密度估計超出範圍,需要對邊界做特殊處理。

回憶二項分布等離散型分布要以常態或其他連續型分布做近似計算時,有所謂「連續性校正」, 其觀念就是將 [X = k] 的單點機率變成 [k-1/2 < X ≦ k+1/2] 的連續型機率。核密度估計其實就是把 Xi 觀測值依 K((x-Xi)/h)/h 這個核函數擴展至其覆蓋範圍。因為 Xi 是自連續型 X 群體抽出的,沒有像二項分布那樣的自然擴展範圍;而群體分布 f(x) 未知,因此我們只能任意性地選取一個核函數 K 及一個帶寬參數 h 把 Xi 這一單點擴展成連續型的機率密度。如此一來,如果 h 小至各個以 Xi 為中心的核函數相互幾乎無重疊,則 f*(x) 幾乎就像原本分立在諸 Xi 的條柱;如果 Xi 如直方圖一般分組,而諸 Xi 分散在各組的組中點,f*(x) 就像直方圖一般,只是各條柱不是水平頂而是曲線型頂(除非核函數是均勻密度);所以一般 f*(x) 就是把直方圖變成平滑曲線圖。

核密度估計的評估採用平均總合平方誤差 (mean integrated squared error, MISE)

MISE(h) = E[∫(f*(x)-f(x))^2 dx]

式中 f(x) 是想要估計的目標:未知的群體分布 p.d.f.,而 f*(x) 就是上述 f(x) 的 KDE。這判準主要是用於帶寬 h 的選擇,因為 h 愈小  f*(x) 愈「純淨」, 基本上由觀測值決定了 f*(x) 的形狀;h 愈大則愈平滑,受核函數 K(x) 形狀影響愈大。學者得出

MISE(h) = R(K)/(nh) + (h^4/4)(m2(K))^2R(f") + o((nh)^(-1)+h^4)

前兩項之和,也就是小 o 之外的部分是漸近的平均總合平方誤差 (asymptotic MISE, AMISE)。期望 AMISE 最小,得出 h 應取

h*= {R(K)/[(m2(K))^2R(f")n]}^(1/5)

其中

R(g) = ∫ (g(x))^2 dx, g = K, f";   m2(K) = ∫ x^2 K(x) dx

帶寬採用 h* 則 AMISE 等級,比參數化方法的 O(1/n) 等級均方誤差低,也就是收斂速度較慢。但 f 既是未知,其二階導數 f" 更是不得而知,因此上述計算最適帶寬 h* 的方法沒辦法直接套用。有許多基於資料的自動化決定 h 的方法被提出;另外,也有許多調適型、變動的 h(h 非常數,而是隨著 Xi 而變)的方法被提出。

核密度估計也可以用來估計多變量聯合機率密度,例如樣本是 (Xi, Yi), i = 1, ..., n,取標準雙變量核函數

K(x) = [1/(2π)]e^(-x'x/2),   x' = (x,y)

而估計的聯合機率密度函數是

f*(x,y) = 1/(nhk) Σ K((x - Xi)/(hk)),  Xi' = (Xi, Yi)

arrow
arrow
    全站熱搜
    創作者介紹
    創作者 等死的老賊 的頭像
    等死的老賊

    劉應興的部落格

    等死的老賊 發表在 痞客邦 留言(0) 人氣()