成大研發快訊

在統計與機率相關研究領域中，如何估計密度函數是重要且具有挑戰性的問題。密度函數的估計方法可以被應用於許多其他的研究領域，例如：區別分析 (discriminant analysis)、影像處理 (image processing)、投影跟蹤 (projection pursuit) 和凸塊獵取 (bump hunting)。在高能物理領域，數據密度函數的凸塊(bumps)可提供有關基本粒子的證據，詳情參見Good and Gaskins [Journal of American Statistical Association (1980), 75, 42-73]。

令 x 為隨機向量 x=(x1,...,xd)所組成之 n×d (d≥1) 資料矩陣，其中 x1,...,xn 是來自 d 維度密度函數 ƒ(x) 之獨立隨機向量。在本論文中，除非有特別定義， ∑ 加總範圍是由1到 n，Π 連乘範圍是由1到 d，而 ∫ 積分區域為

。令xij代表 X 第 i 行第 j 列之元素，則 ƒ(x) 的可變乘積核密度估計量 (variable product kernel estimator) 可表示為

，其中

為單變量對稱機率密度 (稱為核函數) 且 h=(h1,...,hd) 為可變帶寬 (variable bandwidth)，選取適當帶寬 h 是估計 ƒ(x) 最關鍵的一環。有別於一般固定帶寬(fixed bandwidth) 的核密度估計量，此估計量

允許每一個觀測值有不同的帶寬，較固定帶寬的估計量更具有彈性。在資料密集的區域，可選取較小的帶寬以降低估計的偏誤 (bias)；在資料稀疏的區域，可選取較大的帶寬以降低估計的變異數(variance)。

我們使用均方積分誤差 (mean integrated square error，縮寫為MISE)

來衡量

估計 ƒ(x) 的總體正確程度。令 h=hb 其中 b=(b1,...,bn) ，我們稱 b 為局部帶寬因子 (local bandwidth factor)，h 為總體平滑參數 (global smoothing parameter)。對於估計 b 和 h 的問題，一個常用的策略是先由資料選出 b ，將 b 視為固定常數代回 M(h) 可得 M(hb) ，再最小化 M(hb) 的估計

來選取。

Abramson [Annals of Statistics (1982), 10, 1217–1223] 建議使用平方根法則(square-root law)來選局部帶寬因子 b，使用平方根法則之估計量的局部偏誤收斂的速度較未使用平方根法則之估計量來的快。另外，Sain and Scott [Journal of American Statistical Association (1996), 91, 1525–1534] 和Sain [Computational Statististics & Data Analysis (2002), 39, 165–186]等人也提出不同的估計方法，並獲得類似的理論結果。

本篇論文中，我們提出一個全新的估計方法來估計 b 及 h。創造出此估計方法的靈感，既不是從降低估計偏誤的角度出發，也不是由正確的數學演算分析而來。具體而言，我們的方法 (i) 應用了群集分析 (cluster analysis) 的技巧，將資料分為數個群集 (cluster)，然後根據每個資料 xi 在樹狀圖 (dendrogram) 上的平均距離等級來估計 bi；(ii)修飾Wu and Tsai [Probability Theory & Related Fields (2004), 129, 537–558] 所提出之穩定固定帶寬選取法 (stabilized fixed-bandwidth selector) 來選取總體平滑參數 h。在我們的方法中，選取 bi 的部分我們稱之為平均群集法 (average cluster method)。

選取 bi 的兩個步驟為：(i) 應用凝聚式階層群集分析演算法 (agglomerative hierarchical clustering algorithm) 之平均連結法 (average linkage method) 將個觀測值做分群的動作，得到樹狀圖；(ii) 令 bi=xi 在樹狀圖中的平均距離等級。準確地來說，若 ni 代表 xi 所在群體被合併到更大群體之總次數， l1,...,ln 各別代表 ni 次合併發生時的距離等級，則

。

上述方法中，群集是由個別觀測值合併與他最接近觀測值或群集而形成， bi 可以視為在形成巢狀序列群集

的過程中之平均合併距離等級。明顯地，當 ni 較大 ( xi 所在群體被合併的總次數較多) bi 會較小，且多數的 lj 也會偏小，代表 xi 在這組資料中是處於資料相對密集的區域；相反地，當 ni 較小 ( xi 所在群體被合併的總次數較少) bi 會較大，多數的 lj 就會偏大，代表 xi 在這組資料中是處於資料相對稀疏的區域。這足以解釋本論文中所提出的方法為何確實可行。

本文剩下的部份，我們將把均方積分誤差 (MISE) 表示為

，即

為 h 的函數，b 被視為已知常數或是未知常數但其值固定。

表示任意函數

的 d 維度傅立葉轉換 (Fourier transform)，其中 t=(t1,...,td)。

代表資料的樣本特徵函數 (sample characteristic function)。此外，定義四個後面會用到的符號，

，

和

。我們可以證明

是

的漸近不偏(asymptotically unbiased)估計量(大樣本時 E[

] 和

僅差一個常數)，其中

=A+B-C ，

，

及

。此處，

為 d 維度之隨機矩形區域且

為Wu and Tsai (2004) 所提出之截斷頻率 (cut-off frequency)。令 h0 表示最小化

所得到之未知帶寬 (也就是最佳總體平滑參數)。我們將利用最小化

所得到之帶寬

來估計 h0。

我們執行大量的模擬研究 (simulation studies) 來說明我們所提出之可變帶寬選取法

(此處， b 代表平均群集法所選出之局部帶寬因子)確實優於其它的可變帶寬選取法。比較的對象有Abramson (1982) 之平方根準則帶寬選取法

與Sain and Scott (1996) 及Sain (2002) 提出之儲存箱法 (binned method)

。為了數值運算與直覺上的好處，模擬結果中核函數都使用常態分配且所有最小化的過程都在 h∈(0,2] 之間搜尋。

對於不同樣本數 n，皆由下列五種常態混合分配 (normal mixture densities) 中重複抽取100組資料：(#1) N(0,1) ；(#2) (強偏態分配，strongly skewed)

；(#3) (分開的雙峰分配，separated bimodal) 0.5N(-3/2,1/4) + 0.5N(3/2,1/4) ；(#4) N(0,0,1,1,0) (#5) (雙峰分配，bimodal) 0.5N(-6/5,0,1,1,0) + 0.5N(6/5,0,1,1,0) 。其中，(#1) 到 (#3) 是單變量分配而 (#4) 到 (#5) 是二變量分配。我們所考慮之樣本數在單變量分配 (#1) 到 (#3) 為 n=50,400,900，二變量分配 (#4) 到 (#5) 為 n=100,400,900。所有隨機樣本都是利用FORTRAN內建函數RAND生成，並且對於每一組樣本皆利用快速傅立葉轉換 (fast Fourier transform) 來計算樣本特徵函數。

為了比較各個選取法所選出帶寬

之間的表現，我們選擇比較

，即MISE M(

) 之樣本平均。表1及表2匯集了這些模擬結果並顯示出，無論是單變量分配或是二變量分配的狀況下，(i) 對於任意樣本數 n ，

都是最小的，也就是

表現最佳；(ii) 當樣本數 n 增加，

收斂到零的速度遠比其他兩個方法快。因此，

在許多實務上的狀況應該能有相當可靠穩定的表現。

表1. 單變量分配 #1 到 #3的模擬結果

表2. 二變量分配 #4 到 #5的模擬結果

在大部分所考慮的狀況中， M(

) 的樣本變異最大。三個單變量分配的模擬情況下， M(

) 和 M(

) 的樣本變異差不多：分配#1和#3的情況下， M(

) 樣本變異較小；分配#2的情況下， M(

) 表現較佳。在二變量分配的情況下， M(

) 在樣本數 n 為400及900時表現最佳，而在 n=100 與 M(

) 表現差不多。

在本篇論文中，雖然我們無法從理論角度證明所提出之估計量優於其他估計量，但是模擬結果顯示出此估計量是優於其他估計量，也因此我們未來的工作將包含：(i) 建立所提出之估計量的相關理論結果及 (ii) 推廣我們的方法以觧決密度函數偏導數估計時可變帶寬的選取問題和密度函數偏導數平方積分 (integrated squared density partial derivatives)的估計問題。