非參數估計:核密度估計KDE
核密度估計Kernel Density Estimation(KDE)概述
密度估計的問題
由給定樣本集合求解隨機變量的分布密度函數問題是概率統計學的基本問題之一。解決這一問題的方法包括參數估計和非參數估計。
參數估計
參數估計又可分為參數回歸分析和參數判別分析。在參數回歸分析中,人們假定數據分布符合某種特定的性態,如線性、可化線性或指數性態等,然后在目標函數族中尋找特定的解,即確定回歸模型中的未知參數。在參數判別分析中,人們需要假定作為判別依據的、隨機取值的數據樣本在各個可能的類別中都服從特定的分布。經驗和理論說明,參數模型的這種基本假定與實際的物理模型之間常常存在較大的差距,這些方法並非總能取得令人滿意的結果。
[參數估計:最大似然估計MLE][參數估計:文本分析的參數估計方法]
非參數估計方法
由於上述缺陷,Rosenblatt和Parzen提出了非參數估計方法,即核密度估計方法。由於核密度估計方法不利用有關數據分布的先驗知識,對數據分布不附加任何假定,是一種從數據樣本本身出發研究數據分布特征的方法,因而,在統計學理論和應用領域均受到高度的重視。
核密度估計(kernel density estimation)是在概率論中用來估計未知的密度函數,屬於非參數檢驗方法之一,由Rosenblatt (1955)和Emanuel Parzen(1962)提出,又名Parzen窗(Parzen window)。Ruppert和Cline基於數據集密度函數聚類算法提出修訂的核密度估計方法。
核密度估計在估計邊界區域的時候會出現邊界效應。
[https://zh. wikipedia.org/zh-hans/核密度估計]因此,一句話概括,核密度估計Kernel Density Estimation(KDE)是在概率論中用來估計未知的密度函數,屬於非參數檢驗方法之一。
在密度函數估計中有一種方法是被廣泛應用的——直方圖。如下圖中的第一和第二幅圖(名為Histogram和Histogram, bins shifted)。直方圖的特點是簡單易懂,但缺點在於以下三個方面:密度函數是不平滑的;密度函數受子區間(即每個直方體)寬度影響很大,同樣的原始數據如果取不同的子區間范圍,那么展示的結果可能是完全不同的。如下圖中的前兩個圖,第二個圖只是在第一個圖的基礎上,划分區間增加了0.75,但展現出的密度函數卻看起來差異很大;直方圖最多只能展示2維數據,如果維度更多則無法有效展示。
核密度估計有多種內核,圖3(Tophat Kernl Density)為不平滑內核,圖4(Gaussian Kernel Density,bandwidth=0.75)為平滑內核。在很多情況下,平滑內核(如高斯核密度估計,Gaussian Kernel Density)使用場景較多。
雖然采用不同的核函數都可以獲得一致性的結論(整體趨勢和密度分布規律性基本一致),但核密度函數也不是完美的。除了核算法的選擇外,帶寬(bandwidth)也會影響密度估計,過大或過小的帶寬值都會影響估計結果。如上圖中的最后三個圖,名為Gaussian Kernel Density,bandwidth=0.75、Gaussian Kernel Density,bandwidth=0.25、Gaussian Kernel Density,bandwidth=0.55.
核密度估計的應用場景
股票、金融等風險預測:在單變量核密度估計的基礎上,可以建立風險價值的預測模型。通過對核密度估計變異系數的加權處理,可以建立不同的風險價值的預測模型。
密度估計中應用較多的算法是高斯混合模型以及基於近鄰的核密度估計。高斯混合核密度估計模型更多會在聚類場景中應用。
[核密度估計Kernel Density Estimation(KDE)]
核密度分析可用於測量建築密度、獲取犯罪情況報告,以及發現對城鎮或野生動物棲息地造成影響的道路或公共設施管線。可使用 population 字段根據要素的重要程度賦予某些要素比其他要素更大的權重,該字段還允許使用一個點表示多個觀察對象。例如,一個地址可以表示一棟六單元的公寓,或者在確定總體犯罪率時可賦予某些罪行比其他罪行更大的權重。對於線要素,分車道高速公路可能比狹窄的土路產生更大的影響,高壓線要比標准電線桿產生更大的影響。[ArcGIS中的介紹]
熱力圖大家一定聽說過,其實熱力圖就是核密度估計。總而言之,核密度就是用來估計密度的,如果你有一系列空間點數據,那么核密度估計往往是比較好的可視化方法
核密度估計
所謂核密度估計,就是采用平滑的峰值函數(“核”)來擬合觀察到的數據點,從而對真實的概率分布曲線進行模擬。
核密度估計(Kernel density estimation),是一種用於估計概率密度函數的非參數方法,為獨立同分布F的n個樣本點,設其概率密度函數為f,核密度估計為以下:
K(.)為核函數(非負、積分為1,符合概率密度性質,並且均值為0)。有很多種核函數,uniform,triangular, biweight, triweight, Epanechnikov,normal等。
h>0為一個平滑參數,稱作帶寬(bandwidth),也看到有人叫窗口。
Kh(x) = 1/h K(x/h). 為縮放核函數(scaled Kernel)。
核密度函數的原理比較簡單,在我們知道某一事物的概率分布的情況下,如果某一個數在觀察中出現了,我們可以認為這個數的概率密度很大,和這個數比較近的數的概率密度也會比較大,而那些離這個數遠的數的概率密度會比較小。
基於這種想法,針對觀察中的第一個數,我們可以用K去擬合我們想象中的那個遠小近大概率密度。對每一個觀察數擬合出的多個概率密度分布函數,取平均。如果某些數是比較重要的,則可以取加權平均。需要說明的一點是,核密度的估計並不是找到真正的分布函數。
Note: 核密度估計其實就是通過核函數(如高斯)將每個數據點的數據+帶寬當作核函數的參數,得到N個核函數,再線性疊加就形成了核密度的估計函數,歸一化后就是核密度概率密度函數了。
以下面3個數據點的一維數據集為例:5, 10, 15
KDE核函數k(.)
理論上,所有平滑的峰值函數均可作為KDE的核函數來使用,只要對歸一化后的KDE而言(描繪在圖上的是數據點出現的概率值),該函數曲線下方的面積和等於1即可。
只有一個數據點時,單個波峰下方的面積為1,存在多個數據點時,所有波峰下方的面積之和為1。概而言之,函數曲線需囊括所有可能出現的數據值的情況。
常用的核函數有:矩形、Epanechnikov曲線、高斯曲線等。這些函數存在共同的特點:在數據點處為波峰;曲線下方面積為1。
單個數據點(只有一個數據時)所對應的這些核函數
[概率論:高斯/正態分布 ]
sklearn中實現的核函數
sklearn核函數形式
-
Gaussian kernel (
kernel = 'gaussian'
) -
Tophat kernel (
kernel = 'tophat'
)if
-
Epanechnikov kernel (
kernel = 'epanechnikov'
) -
Exponential kernel (
kernel = 'exponential'
) -
Linear kernel (
kernel = 'linear'
)if
-
Cosine kernel (
kernel = 'cosine'
)if
wekipedia上各種核函數的圖形
均勻核函數 k(x)=1/2,-1≤x≤1 加入帶寬h后: kh(x)=1/(2h),-h≤x≤h
三角核函數 k(x)=1-|x|,-1≤x≤1 加入帶寬h后: kh(x)=(h-|x|)/h^2,-h≤x≤h
伽馬核函數 kxi(x)=[x^(α-1)exp{-xα/xi}]/[(xi/α)^α.Γ(α)]
[https://zh.wikipedia.org/zh-hans/%E6%A0%B8%E5%AF%86%E5%BA%A6%E4%BC%B0%E8%AE%A1]
不同內核的比較
Epanechnikov 內核在均方誤差意義下是最優的,效率損失也很小。
由於高斯內核方便的數學性質,也經常使用 K(x)= ϕ(x),ϕ(x)為標准正態概率密度函數。
對於多個數據點的KDE曲線:由於相鄰波峰之間會發生波形合成,因此最終所形成的曲線形狀與選擇的核函數關系並不密切。考慮到函數在波形合成計算上的易用性,一般使用高斯曲線(正態分布曲線)作為KDE的核函數。
KDE算法:索引樹
lz發現sklearn算法實現中有一個參數是算法項,如algorithm='auto',想了一下是為了加速。
KDE的概率密度函數公式得到后
有了上述公式之后,只需遍歷輸出圖像的每一個點,計算其核密度估計值即可。
但是稍微想一下就發現這個程序太冗余了,如果有很多點(n很大),並且輸出圖像很大,那么每一個像素都需要進行n個累積的加法運算,並且大部分都是+0(因為一般來說,一個點附近的點不會很多,遠遠小於n,其余大部分點與這個像素的距離都大於r),這樣就造成了冗余計算。
解決方案當然也非常簡單,就是建立一個索引,然后在計算某個像素的核密度估計值時利用索引搜索出附近的點,然后累積這些點的核函數即可。
如Dotspatial自帶了多種空間索引,有R樹,R*樹,KD樹等;sklearn自帶了kd tree, ball tree等等。
如果只需找出附近的點,對索引要求不高,任意一個索引都能使用。
[ 空間點雲核密度估計算法的實現-以Dotspatial為基礎GIS庫]KDE帶寬h
如何選定核函數的“方差”呢?這其實是由帶寬h來決定,不同的帶寬下的核函數估計結果差異很大。
帶寬反映了KDE曲線整體的平坦程度,也即觀察到的數據點在KDE曲線形成過程中所占的比重。帶寬越大,觀察到的數據點在最終形成的曲線形狀中所占比重越小,KDE整體曲線就越平坦;帶寬越小,觀察到的數據點在最終形成的曲線形狀中所占比重越大,KDE整體曲線就越陡峭。
還是以上面3個數據點的一維數據集為例,如果增加帶寬,那么生成的KDE曲線就會變平坦:
如果進一步增加帶寬,那么KDE曲線在變平坦的同時,還會發生波形合成:
相反,如果減少帶寬,那么KDE曲線就會變得更加陡峭:
從數學上來說,對於數據點Xi,如果帶寬為h,那么在Xi處所形成的曲線函數為(其中K為核函數):
在上面的函數中,K函數內部的h分母用於調整KDE曲線的寬幅,而K函數外部的h分母則用於保證曲線下方的面積符合KDE的規則(KDE曲線下方面積和為1)。
帶寬的選擇
帶寬的選擇很大程度上取決於主觀判斷:如果認為真實的概率分布曲線是比較平坦的,那么就選擇較大的帶寬;相反,如果認為真實的概率分布曲線是比較陡峭的,那么就選擇較小的帶寬。
帶寬計算好像也有相應的方法,如R語言中計算帶寬時,默認采用”nrd0″方法。
如何選擇h?顯然是選擇可以使誤差最小的。下面用平均積分平方誤差(mean intergrated squared error)的大小來衡量h的優劣。
在weak assumptions下,MISE (h) =AMISE(h) + o(1/(nh) + h4) ,其中AMISE為漸進的MISE。而AMISE有,
,
其中,
,
為了使MISE(h)最小,則轉化為求極點問題,
當核函數確定之后,h公式里的R、m、f''都可以確定下來,有(hAMISE ~ n−1/5),AMISE(h) = O(n−4/5)。
如果帶寬不是固定的,其變化取決於估計的位置(balloon estimator)或樣本點(逐點估計pointwise estimator),由此可以產產生一個非常強大的方法稱為自適應或可變帶寬核密度估計。
[核密度估計(Kernel density estimation) ]在選擇合適的核函數及帶寬后,KDE可以模擬真實的概率分布曲線,並得到平滑而漂亮的結果。以近200個點的CPU使用率為例,使用KDE繪制的結果為:

[一維數據可視化:核密度估計(Kernel Density Estimates)]
核密度估計的實現
Python中KDE的實現:sklearn
[sklearn.neighbors.
KernelDensity
(bandwidth=1.0, algorithm='auto', kernel='gaussian', metric='euclidean', atol=0, rtol=0, breadth_first=True, leaf_size=40, metric_params=None)
from sklearn.neighbors import kde import numpy as np X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]]) kde = kde.KernelDensity(kernel='gaussian', bandwidth=0.2).fit(X) print(kde.score_samples(X)) print(np.exp(kde.score_samples(X)))[-0.41075698 -0.41075698 -0.41076071 -0.41075698 -0.41075698 -0.41076071]
[ 0.66314807 0.66314807 0.6631456 0.66314807 0.66314807 0.6631456 ]
score_samples(X)
-
Evaluate the density model on the data.
Parameters: X : array_like, shape (n_samples, n_features)
kde.score_samples(X)返回的是點x對應概率的log值,要使用exp求指數還原。
Note: 還原后的所有點的概率和范圍是[0, 無窮大],只是說一維數據線下面的面積或者二維數據面下面的體積和為1。
[Density Estimation¶]
[sklearn.neighbors.KernelDensity¶]
spark中KDE的實現
MLlib中,僅僅支持以高斯核做核密度估計。
[核密度估計]
R中KDE的實現
在R語言中,KDE的繪制是通過density()函數來實現的 — 通過density()函數計算得到KDE模型,然后再使用plot()函數對KDE曲線進行繪制:
x <- c(5, 10, 15)
plot(density(x))
出於兼容性上的考慮,R語言中density()函數在計算帶寬時,默認采用”nrd0″方法。不過,根據R語言的幫助文檔,帶寬參數bw應該顯式聲明為其它更合適的方法,比如”SJ”:
plot(density(x, bw="SJ"))
對於調整帶寬,除了修改bw參數,還可以通過設定adjust參數來進行擴大或縮小:
plot(density(x, bw="SJ", adjust=1.5))
在上面的例子中,最終使用的帶寬將是采用”SJ”方法計算得到的帶寬的1.5倍。adjust參數的默認值為1,也即既不擴大、也不縮小。
至於核函數,density()默認采用高斯曲線。可以通過設定kernel參數來更改核函數。比如:plot(density(x, bw="SJ", kernel="epanechnikov"))
density()函數接受以下7個核函數選項:
gaussian。高斯曲線,默認選項。在數據點處模擬正態分布。
epanechnikov。Epanechnikov曲線。
rectangular。矩形核函數。
triangular。三角形核函數。
biweight。
cosine。余弦曲線。
optcosine。
from: http://blog.csdn.net/pipisorry/article/details/53635895
ref: [有邊界區間上的核密度估計]