非引數估計之 kernel density estimation (核密度估計)

摘要： 在概率密度估計過程中，如果我們隊隨機變數的分佈是已知的，那麼可以直接使用引數估計的方法進行估計，如最大似然估計方法。 然而在實際情況中，隨機變數的引數是未知的,因此需要進行非引數估計.核密度估計是非引數估計的一種方法，也就是大家經常聽見的parzen 窗方法了. 本文主要介紹 非引數...

在概率密度估計過程中，如果我們隊隨機變數的分佈是已知的，那麼可以直接使用引數估計的方法進行估計，如最大似然估計方法。

然而在實際情況中，隨機變數的引數是未知的,因此需要進行非引數估計.核密度估計是非引數估計的一種方法，也就是大家經常聽見的parzen 窗方法了.

本文主要介紹 非引數估計的過程 以及 parzen窗方法估計概率密度的過程.

非引數估計

過程

如圖1所示，對於一個未知的概率密度函式 ,某一個隨機變數 落在區域 裡的概率可以表示成如下形式：

由上面的公式我們可以得到 的期望為：

同時，當 足夠大時，我們可以近似地認為 可以作為 的一個近似值。

然後，假設 足夠大， 足夠小，並且假設 是連續的，那麼我們可以得到：

其中 是區域 中的一個點， 是 的面積(體積)，結合上述4個式子，得：

因此，某一小區域內的概率密度函式就可以用上述公式表示了。

問題

我們再看一下公式（６）:

顯然我們估計的這個概率密度是一個平滑的結果,即當 選擇的越大,估計的結果和真實結果相比就越平滑;因此看起來我們需要把 設定的小一點,然而如果我們把 選擇的過小,也會出現問題:太小的 會導致這塊小區域裡面沒有一個點落在裡面,因此就會得到該點的概率密度為0;另外,假設剛好有一個點落在了這個小區域裡,那由於V過於小,我們計算得到的概率密度可能也會趨近於無窮,兩個結果對於我們來說都是沒有太大意義的．

從 實際的角度 來看，我們獲取的資料量一定是有限的，因此體積 不可能取到無窮小，我們可以總結下，使用非引數概率密度估計有以下兩方面限制，且是不可避免的：

1. 在有限資料下，使用非引數估計方法計算的概率密度一定是真實概率密度平滑後的結果.
2. 在有限資料下，體積趨於無窮小計算的概率密度沒有意義．

從 理論的角度 來看，我們希望知道如果有無限多的取樣資料，那麼上述兩個限制條件應該怎樣克服？假設我們使用下面的方法來估計點 處的概率密度: 構造一系列包含 的區域 ,其中 中包含一個樣本， 中包含 個樣本．則：

其中 表示第 次估計結果，如果要求 能夠收斂到 ,則需要滿足下面三個條件：

• 

• 

• 

parzen窗法

原理

假設 是一個 維的超立方體(hypercube),且其邊長為 , 那麼我們可以用如下公式表示 :

然後我們再定義一個窗函式（window function）:

定義了一個以圓點為中心的單位超立方體,這樣我們就可以用 來表示體積 內的樣本個數:

解釋

網上經常會見到這樣的解釋: 某一點的概率密度是其他樣本點在這一點的概率密度分佈的平均值 .還有這樣一張圖:

從公式(13)(14)我們可以看出,當 很大的時候, 就是一個 矮胖 的函式,由公式(14)將每個樣本點在點 處的貢獻取平均之後,點 處的概率密度就是一個非常平滑的結果; 當 太小的時候, 就是一個 高瘦 的函式,由公式(14)將每個樣本點在點 處的貢獻取平均之後,點 處的概率密度就是一個受噪聲影響非常大的值,因此估計的概率密度平滑性就很差,反而和真實值差的很遠.這兩點和1.2節總結的兩點缺陷正好吻合.

模擬實驗

如下我做了兩個模擬

1. 第一個是生成了均值是0，方差是2的服從高斯分佈的資料，分別使用bandwidth為0.1, 1, 5三個值進行估計
2. 第二個是生成了100個服從高斯混合分佈的資料，分別是均值為0，方差為1以及均值為5，方差為1的兩個高斯混合模型，兩者相互獨立。

這裡核函式選的是高斯核。

bandwidth = 0.1	bandwidth = 1	bandwidth = 5
0.1-1.png	1-1.png	5-1.png
bandwidth = 0.1	bandwidth = 0.5	bandwidth = 1
{ 0.1-2.png	0.5-2.png	1-2.png

可以很明顯得看到估計的概率密度是如何受到bandwidth影響的，當bandwidth選擇的太小，則估計的密度函式受到噪聲影響很大，這種結果是不能用的；當bandwidth選擇過大，則估計的概率密度又太過於平滑。總之，無論bandwidth過大還是過小，其結果都和實際情況相差的很遠，因此合理地選擇bandwidth是很重要的。

下面是高斯混合分佈密度估計的程式碼。

import matplotlib.pyplot as plt 
import numpy as np 
from sklearn.neighbors.kde import KernelDensity
from scipy.stats import norm

if __name__ == "__main__":
np.random.seed(1)

x = np.concatenate((np.random.normal(0, 1, int(0.3*100)), np.random.normal(5, 1, int(0.7*100))))[:, np.newaxis]
plot_x = np.linspace(-5, 10, 1000)[:, np.newaxis]
true_dens = 0.3*norm(0, 1).pdf(plot_x) + 0.7*norm(5, 1).pdf(plot_x)

log_dens = KernelDensity(bandwidth=1).fit(x).score_samples(plot_x)

plt.figure(),
plt.fill(plot_x, true_dens, fc='#AAAAFF', label='true_density')
plt.plot(plot_x, np.exp(log_dens), 'r', label='estimated_density')
for _ in range(x.shape[0]):
plt.plot(x[:, 0], np.zeros(x.shape[0])-0.01, 'g*') 
plt.legend()

plt.show()

參考資料

[1.] scikit-learn.org/stable/auto_examples/neighbors/plot_kde_1d.html#sphx-glr-auto-examples-neighbors-plot-kde-1d-py" target="_blank" rel="nofollow,noindex">sklearn:Simple 1D Kernel Density Estimation

[2.] Richard O. Duda, Peter E. Hart, and David G. Stork. 2000. Pattern Classification (2nd Edition). Wiley-Interscience, New York, NY, USA.

[3. ] Kernel density estimation

[4.] 邊肇祺, 張學工, 2000. 模式識別. 清華大學出版社.