核密度估计Kernel Density Estimation(KDE)-代码详细解释
?
在介紹核密度評估Kernel Density Estimation(KDE)之前,先介紹下密度估計的問題。由給定樣本集合求解隨機變量的分布密度函數問題是概率統計學的基本問題之一。解決這一問題的方法包括參數估計和非參數估計。(對于估計概率密度,如果確定數據服從的分布類型,可以使用參數擬合,否則只能使用非參數擬合)
參數估計又可分為參數回歸分析和參數判別分析。在參數回歸分析中,人們假定數據分布符合某種特定的性態,如線性、可化線性或指數性態等,然后在目標函數族中尋找特定的解,即確定回歸模型中的未知參數。在參數判別分析中,人們需要假定作為判別依據的、隨機取值的數據樣本在各個可能的類別中都服從特定的分布。經驗和理論說明,參數模型的這種基本假定與實際的物理模型之間常常存在較大的差距,這些方法并非總能取得令人滿意的結果。
?
由于上述缺陷,Rosenblatt和Parzen提出了非參數估計方法,即核密度估計方法。由于核密度估計方法不利用有關數據分布的先驗知識,對數據分布不附加任何假定,是一種從數據樣本本身出發研究數據分布特征的方法,因而,在統計學理論和應用領域均受到高度的重視。
因此,一句話概括,核密度估計Kernel Density Estimation(KDE)是在概率論中用來估計未知的密度函數,屬于非參數檢驗方法之一。
在密度函數估計中有一種方法是被廣泛應用的——直方圖。如下圖中的第一和第二幅圖(名為Histogram和Histogram, bins shifted)。直方圖的特點是簡單易懂,但缺點在于以下三個方面:
核密度估計有多種內核,圖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.
上圖為使用Python的sklearn實現,算法為KernelDensity。代碼如下:
#coding:utf-8 import numpy as np import matplotlib.pyplot as plt from sklearn.neighbors import KernelDensity np.random.seed(1) N = 20 X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),np.random.normal(5, 1, int(0.7*N)) ))[:, np.newaxis] #np.newaxis是None的意思 #前半部分的平均值是0,方差是1 #int(0.3 * N)指的是輸出多少數量符合要求的數據 #---------------------------以上是創建數據集-----------------------------------------------------------------X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]#創建等差數列,在-5和10之間取1000個數 bins = np.linspace(-5, 10, 10)#這個的作用是,在相鄰兩個邊界時間的數據對應的y值都一樣大 print("bins=",bins) fig, ax = plt.subplots(2, 2, sharex=True, sharey=True) fig.subplots_adjust(hspace=0.05, wspace=0.05) # 直方圖 1 'Histogram' print("---------------------------------") ax[0, 0].hist(X[:, 0], bins=bins, fc='#AAAAFF', normed=True)#fc指的應該是顏色的編碼 #這里的ax[0,0]的意思是畫在第幾副圖上ax[0, 0].text(-3.5, 0.31, 'Histogram')#-3.5, 0.31的意思是每張圖的logo要畫在什么地方 # 直方圖 2 'Histogram, bins shifted' ax[0, 1].hist(X[:, 0], bins=bins + 0.75, fc='#AAAAFF', normed=True)#histogram的縮寫 ax[0, 1].text(-3.5, 0.31, 'Histogram, bins shifted')#每個子圖內畫標簽#-----------------------------------------------------------------------------------# 核密度估計 1 'tophat KDE' kde = KernelDensity(kernel='tophat', bandwidth=0.75).fit(X)#什么是帶寬 log_dens = kde.score_samples(X_plot) #所以這里有兩組數據,X和X_plot,其實是在利用X_plot對X進行采樣。 #所以想要復用這段代碼的時候,改X即可,X_plot不用修改ax[1, 0].fill(X_plot[:, 0], np.exp(log_dens), fc='#AAAAFF')#fill就是用來畫概率密度的 ax[1, 0].text(-3.5, 0.31, 'Tophat Kernel Density')#設置標題的位置 # 核密度估計 2 'Gaussian KDE' kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(X) log_dens = kde.score_samples(X_plot)#返回的是點x對應概率的log值,要使用exp求指數還原。 ax[1, 1].fill(X_plot[:, 0], np.exp(log_dens), fc='#AAAAFF')#fill就是用來畫概率密度的 #所以上面一句代碼就非常清晰了,X_plot[:, 0]是具體數據,np.exp(log_dens)指的是該數據對應的概率 ax[1, 1].text(-3.5, 0.31, 'Gaussian Kernel Density')#設置標題的位置print("ax.ravel()=",ax.ravel()) print("X.shape[0]=",X.shape[0]) print("X=",X)#這個是為了在每個子圖的下面畫一些沒用的標記,不看也罷 for axi in ax.ravel():axi.plot(X[:, 0], np.zeros(X.shape[0])-0.01, '+k')axi.set_xlim(-4, 9)#設定上下限axi.set_ylim(-0.02, 0.34)##畫圖過程是兩行兩列,這里是遍歷第1列,每個位置的左側畫一個“xNormalized Density” for axi in ax[:, 0]:print("axi=",axi)axi.set_ylabel('Normalized Density')##畫圖過程是兩行兩列,這里是遍歷第2行,每個位置畫一個“x” for axi in ax[1, :]:axi.set_xlabel('x') plt.show()#ravel函數的作用如下: # >>> x = np.array([[1, 2, 3], [4, 5, 6]]) # >>> print(np.ravel(x)) # [1 2 3 4 5 6]?
KernelDensity算法包括kd_tree和ball_tree,默認自動選擇。核模型包括gaussian、tophat、epanechnikov、exponential、linear、cosine,默認是gaussian模型。可調整的參數如下:
?
核密度估計的應用場景:
- 股票、金融等風險預測;
總結
以上是生活随笔為你收集整理的核密度估计Kernel Density Estimation(KDE)-代码详细解释的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 先验概率与后验概率、贝叶斯区别与联系
- 下一篇: 根据数据集获取概率密度图像和概率分布图像