王文武， 孔德茹

(曲阜師范大學(xué)統(tǒng)計與數(shù)據(jù)科學(xué)學(xué)院,273165,山東省曲阜市)

0 引 言

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

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

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

1 兩種估計方法

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

1.1 核方法

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

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

在核密度估計的基礎(chǔ)上由導(dǎo)數(shù)定義，對估計得到的密度函數(shù)求r階導(dǎo)[4]，即

1.2 局部多項式回歸方法

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

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

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

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

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

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

1.3 兩種方法的估計效果

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

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

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

通過使積分均方誤差達到最小，得到密度函數(shù)各階導(dǎo)數(shù)的最優(yōu)帶寬為

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)，

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

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

引理1假定密度函數(shù)f(x)至少r+v+1階可導(dǎo)，核函數(shù)滿足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]，所以由林德伯格-萊維中心極限定理可得：核方法和局部多項式回歸方法下估計的概率密度導(dǎo)數(shù)服從漸近正態(tài)分布，關(guān)于更加細致的證明詳見參考文獻[4，9].

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

2 數(shù)值模擬

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

2.1 均勻分布

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

圖1 均勻分布導(dǎo)數(shù)估計的均方誤差

圖2 均勻分布導(dǎo)數(shù)估計結(jié)果

2.2 正態(tài)分布

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

圖3 正態(tài)分布導(dǎo)數(shù)估計的均方誤差

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

圖4 正態(tài)分布導(dǎo)數(shù)估計結(jié)果

2.3 卡方分布

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

圖5 卡方分布導(dǎo)數(shù)估計的均方誤差

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

圖6 卡方分布導(dǎo)數(shù)估計結(jié)果

2.4 小結(jié)

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

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

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

第二，雖然局部多項式回歸方法在正態(tài)分布密度導(dǎo)數(shù)下的估計性質(zhì)不是最優(yōu)，但這并不能否定該方法在另外兩個分布下的良好表現(xiàn). 即在大多分布下，相較于核方法，局部多項式估計方法更穩(wěn)健、高效.

3 討 論

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

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