兩種概率密度導數估計方法的比較研究*

王文武， 孔德茹

(曲阜師范大學統計與數據科學學院,273165,山東省曲阜市)

0 引 言

在探索樣本總體信息時，概率分布函數及密度函數為統計研究的基本問題. 相應地，其估計理論已經得到較好的完善，但密度函數的導數卻鮮少被關注. 雖然密度導數沒有受到足夠的重視，但是其作用不容小覷. 例如，一階導函數為密度函數的斜率，反映了密度函數的波動幅度；二階導函數為密度函數的曲率，反映了密度函數的彎曲程度以及光滑程度. 密度函數的導數不僅直觀反映了概率密度的圖像走勢，而且被廣泛應用于經濟學與統計學等領域：在經濟學中，常需要GDP(國民生產總值)進行建模并做出導數估計，利益的最值問題也需要導數或者偏導數的計算[1]；在統計學中，概率密度函數的導數也發揮了重要作用，例如，導數作為重要的一環，出現在最速下降及牛頓等算法中.

因此，密度函數導數需要得到關注并且估計方法也應被深入研究. 從密度導數本身出發，H?rdle等提出了平均導數估計并給出了實例應用[2]；次年，又提出了一種廣義交叉驗證，通過利用差商的核平滑來估計一階導數[3]. 2009年，Hansen在《Lecture Notes on Nonparametrics》書中總結了密度函數導數的簡單估計方法——根據導數定義對密度函數的導數進行核估計[4]. 2013年，Chacón等提出了多變量數據下的核密度導數估計帶寬選擇器，優化了核方法在導數估計中的應用[5]. 近年來，回歸函數的導數估計受到了重視，這為密度函數導數估計提供了新的思路.Charnigo等在2011年提出通過GCp標準選擇合適的調整參數從而提高導數估計的效果，并且該方法可用于任何非參數回歸方法[6]. 2013年，Brabanter等闡述了經驗一階導數在局部多項式回歸框架中的使用，并推導了偏差和方差[7]. 2016年，Wang等將局部加權最小二乘方法以及差分方法應用到導數估計中，得到了更小的均方誤差[8]. 2019年，Cattaneo等采用局部多項式回歸估計各階導數，給出了較為全面的軟件包，并應用到斷點回歸的臨界點檢驗中[9，10].

在估計概率密度導數時，盡管核方法和局部多項式回歸方法經常被使用，但在應用中并不清楚哪種方法更有優勢. 故本文從理論和模擬兩個角度進行比較研究. 首先，對這兩種非參數估計方法的理論知識進行了梳理總結并做出了合理的比較. 其次，選擇均勻分布、正態分布、卡方分布，在不同樣本量下對概率密度導數進行數值模擬并根據經驗積分均方誤差比較估計效果. 為了更充分地對比估計性質，每種方法均通過兩個軟件包實現. 通過觀察兩種方法在4個R包下均方誤差值的變化以及估計結果圖，最終證實局部多項式回歸方法相較于核方法具有更加優良的性質.

1 兩種估計方法

在總體信息模糊的情況下，常采用非參數方法來達到目的. 常見的密度函數估計方法有直方圖法、核方法、局部多項式回歸法、樣條估計法以及最近鄰估計法等，其中核方法和局部多項式回歸法最為常用[11]. 本部分將簡單介紹核方法和局部多項式回歸方法的一般理論.

1.1 核方法

首先采用核方法估計概率密度函數，即核密度估計. 核密度估計通過核函數來估計未知的概率密度，其基本思想是按照其余觀測值距估計點的遠近對該估計點賦予不同的權重：距估計點較近的觀測值包含該估計點較多的信息，所以賦予更大的權重；距估計點較遠的觀測值包含信息較少，故賦予較小的權重甚至0權重[10].

其中k(·)為核函數.通過核密度估計的最終形式可以看出核密度估計結果實際上是對協變量x附近的觀測值進行加權平均，核函數為權重函數. 另外，核函數需要滿足非負性、對稱性、正則性等條件. 在有關核的應用中，均勻核函數、Bieweight核函數以及Gaussian核函數較為常見.

在核密度估計的基礎上由導數定義，對估計得到的密度函數求r階導[4]，即

1.2 局部多項式回歸方法

局部多項式回歸方法作為非參數統計的一種重要估計方法，應用范圍極其廣泛. 其核心思想是在局部區域用某一多項式近似目標函數，從而得到該函數及其各階導數的估計[12]. 估計過程中用到的多項式通常是采用泰勒展開得到的. 因此，當估計同一函數時，采用不同階數的多項式將得到不同的估計效果[13].

在介紹局部多項式回歸方法估計概率密度導數之前，作出兩個一般假定[9].

假定1假定x1,…,xn是來自分布函數為F(x)的一組隨機觀測值，概率密度函數為f(x)，其中，協變量x的定義域為[xlower,xupper](xupper≠xlower)，分布函數F(x)至少p+1階可導并且F(r+1)(x)=f(r)(x),r≤p.

假定2為了優化估計效果的表達式，本文要求核函數k(·)除滿足一般性質外，還需要滿足該條件：無論核函數為哪種形式，均只在[-1,1]上有意義，其余位置點為0.

注對于假定1的條件，幾乎所有數據都可以滿足.由于距離x越近的觀測點包含的信息越多，給出的權重也應越大[14]，故假定2把核函數的取值范圍限制在了[-1,1]上，即在[x-h,x+h]上核函數非零.對于高斯核函數，雖不能滿足該條件，但在[-1,1]之外的區域，核函數取值接近于0.因此可認為高斯核函數滿足假定2中的限制.此外，當帶寬接近于0時，假定2將核函數限制在估計點附近，從而可優化邊界估計效果.

目標函數F(x)為x的分布函數，yi為F(x)的樣本值，則有局部最小二乘估計

1.3 兩種方法的估計效果

在衡量估計效果時，常采用積分均方誤差MISE(Mean Integrated Squared Error),

從上式可以看出均方誤差既包含了方差也包含了偏差，能夠很好地度量估計性質：積分均方誤差越小說明估計性質越好.

通過變量代換和r次分部積分可以得到核方法下導數估計的偏差和方差[3]，從而積分均方誤差為

通過使積分均方誤差達到最小，得到密度函數各階導數的最優帶寬為

h=Cr,v(k,f)n-1/(1+2r+2v)，Cr,v(k,f)=R(f(r+v))-1/(1+2r+2v)Ar,v(k)，

通過上述結果可以看出，密度函數導數核估計使用的最優帶寬為O(n-1/(1+2r+2v)),而密度函數使用的最優帶寬為O(n-1/(1+2v)).這是因為光滑估計的最優值取決于估計對象和分析目的，故估計的均方誤差會發生變化，最優帶寬也會隨之發生變化[16].

到得草原中來，綠色王國，花草天堂，但在漫無涯際之中，此番我的觀光焦點，不全在草原山崗河流等大景致，反而破天荒地“微觀化”了。

引理1假定密度函數f(x)至少r+v+1階可導，核函數滿足k(s)(∞)=0,k(s)(-∞)=0，s=0,1,2,…,r,h→0，當n→∞時，nh→∞，則核方法的漸近偏差與漸近方差分別為

引理2滿足假定1和假定2，2≤r≤p，h→0，當n→∞時，nh→∞且nh2p+1=O(1)，則局部多項式回歸方法的漸近偏差與漸近方差分別為

因為滿足林德伯格條件[17]，所以由林德伯格-萊維中心極限定理可得：核方法和局部多項式回歸方法下估計的概率密度導數服從漸近正態分布，關于更加細致的證明詳見參考文獻[4，9].

核方法估計密度導數的漸近偏差與漸近方差的證明過程詳見文獻[4]，局部多項式回歸方法估計性質的證明在Cattaneo和Jansson (2019)中有詳盡的敘述[9]，這里僅對結果進行分析. 從式子上看，核方法與局部多項式回歸方法下的估計偏差與方差并不能合理比較，因此在第2節中，通過數值模擬結合實際比較二者的估計性質.

2 數值模擬

本節通過數值模擬的方式，以密度函數的一階導數估計為例，分別采用核方法與局部多項式回歸方法估計，對比同一分布下2種方法的經驗均方誤差，并且直觀呈現密度導數估計結果[13].

2.1 均勻分布

從圖1可以看出，在樣本量一致的情況下，無論使用哪個包估計，局部多項式回歸方法的估計積分均方誤差始終小于核估計的誤差，故對導數估計，局部多項式回歸方法優于核方法. 橫向比較趨勢，局部多項式回歸方法隨著樣本量增加，均方誤差存在大幅度或小幅度的減少；而核方法的兩個估計結果卻越來越差. 在不同樣本量下，lpdensity包的經驗均方誤差值最小，且隨著樣本量增加，該估計的均方誤差越來越小，即大樣本量下，會有更小的均方誤差. 此外，通過圖2，可以看出在4個軟件包中，lpdensity的估計結果最接近于真實值0. 雖然局部多項式回歸方法在邊界處的估計仍存在進步空間，但相較于核方法在-1和1附近的估計結果，局部多項式估計的表現還是明顯更優的. 局部多項式對內點的估計效果也優于核估計.

圖1 均勻分布導數估計的均方誤差

圖2 均勻分布導數估計結果

2.2 正態分布

采用核估計方法和局部多項式估計方法，選擇均值為0，方差為1的正態分布進行100次數值模擬，給出樣本量為2 000時，呈現兩種方法中最優的估計結果見圖3和圖4.

圖3 正態分布導數估計的均方誤差

橫向比較圖3，當樣本量逐漸增大，4個軟件包下的估計均方誤差均存在減小的趨勢，即局部多項式回歸方法和核方法的估計效果越來越好. ks的估計均方誤差隨著樣本量增加越來越接近0，在圖4中的估計結果也十分優秀. 可見對正態分布來說，合適的核方法可以得出較好的估計. 每種樣本量下，局部多項式估計的兩個包的MISE雖然不是最優，但也不是最差. 另外，由于KernSmooth包locpoly函數的帶寬需要提前設定，這可能導致了經驗均方誤差有小幅增加. 通過圖4可以看出在正態分布的密度導數估計中，ks包中的kdde函數在邊界附近和內點的擬合均明顯優于lpdensity的估計. 雖然lpdensity的表現并不最優，但得到的估計曲線較為光滑，趨勢也接近真實概率密度導函數.

圖4 正態分布導數估計結果

2.3 卡方分布

以卡方分布作為偏態分布的代表，對χ2(2)進行100次數值模擬，兩種方法下的經驗均方誤差結果和樣本量為2 000的估計結果與真實值，見圖5和圖6.

圖5 卡方分布導數估計的均方誤差

從以上結果可以清晰地看出偏態分布下，局部多項式估計依然優于核估計. 圖中的數值趨勢與均勻分布類似：無論在哪個包下，局部多項式回歸方法的均方誤差均小于核方法的均方誤差；隨著樣本量增加，MISElp表現越來越好，估計結果越來越接近真實值，而MISEk卻越來越大. 通過圖6，可以直觀看出核估計在邊界處的估計結果與真實值的偏差較大，但局部多項式回歸方法在0附近明顯有更好的表現效果. 在內點，核估計出現了欠光滑，但局部多項式估計幾乎完全擬合真實密度導數，估計性質十分優秀.

圖6 卡方分布導數估計結果

2.4 小結

對比3個分布下的估計結果圖，lpdensity對卡方分布的密度導數估計最優，其次是均勻分布，正態分布. 雖然在正態分布下，ks的估計在邊界點和內點更勝一籌，但局部多項式估計方法在這3種分布下一直表現良好，說明該方法不僅有邊界自適應性，還具有穩健性. 核方法在均勻分布和卡方分布的估計均方誤差較大，這說明核方法不具有穩健性. 當分布未知時，核方法并不是最優選擇. 正態分布下，ks包估計效果優于lpdensity. 但在其他兩個分布下，ks的估計效果并不好，而且與事實相悖：隨著樣本量增加，均方誤差卻存在或大或小的增幅. 這可能與軟件包的內置函數有關，不同R包的最優適應范圍不同，為了更好地適應某類數據需要做出適當調整. 另外，ks包可以用于1至6維數據，在多維數據中也許會具有一定的穩健性.

通過觀察數值模擬得到的3個均方誤差圖以及估計結果對比圖，可以得到以下兩個結論.

第一，由于正態分布的取值范圍是全體實數，局部多項式回歸方法在邊界點擬合效果并沒有完美體現. 但另外兩個分布存在明確的邊界，并且相對于核估計，局部多項式估計更接近真實值. 即在一般情況下，局部多項式估計具有邊界自適應性.

第二，雖然局部多項式回歸方法在正態分布密度導數下的估計性質不是最優，但這并不能否定該方法在另外兩個分布下的良好表現. 即在大多分布下，相較于核方法，局部多項式估計方法更穩健、高效.

3 討 論

本文基于核方法和局部多項式回歸方法，得到概率密度函數的導數估計，并通過數值模擬對比兩種方法下的估計性質，最終得到結論：樣本量一定時，相較于核方法，無論分布為均勻、正態還是偏態，局部多項式回歸方法的估計性質均更加優良. 局部多項式回歸方法無需任何復雜的數據預處理就可以減輕邊界偏倚并且具有邊界自適應性、估計性質穩健等優點. 此外，該方法的理論體系愈見成熟，各階導數估計偏差和估計方差的計算和證明也都有著系統的闡述. 為了進一步對比核方法和局部多項式回歸方法，在第二部分以均勻分布、正態分布和卡方分布為例，在不同樣本量下采用兩種方法估計概率密度函數的一階導數，分析對比了均方誤差和各個點上的估計結果. 從數值與圖形兩個方面，驗證了局部多項式回歸方法不僅在內點的估計性質優越，在邊界上的估計效果也十分良好，整體估計效果明顯優于核估計. 在估計密度函數的導數時，該方法會更加穩健、高效、簡單. 此外，局部多項式估計不僅可以應用在操縱檢驗中[10]，還可以應用到極小化算法中[18]，這會使檢驗結果與算法應用更加精準有效.

在導數估計中，本文僅討論了一元局部多項式回歸模型和核方法的估計性質，未比較兩種方法在多元情形下的效果，這一比較會更加復雜. 另外，因為帶寬選擇對最終估計結果具有顯著影響，所以這將會是另一個有意義的研究課題. 在統計推斷方面，關于局部多項式回歸的置信帶也是一個不錯的探索方向.