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

本文轉載自查看原文 2019-11-03 18:46 671 機器學習

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

張王李劉趙孫楊

0.1922018.11.22 22:17:06字數 1,642閱讀 8,195

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

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

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

非參數估計

過程

圖1

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

$P = \int_{R} P(x^{'}) dx^{'} \tag{1}$

如果 $R$ 足夠窄，我們可以用 $P$ 來表示 $p(x)$ 的一個平均后的結果。假設我們現在有 $n$ 個樣本，
且他們服從獨立同分布，那么 $n$ 個樣本中的 $k$ 個落在區域 $R$ 中的概率可以表示成下面的公式：

$P_k = C_n^k P^k (1-P)^{n-k}　　\tag{2}$

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

$E(k) = n \cdot P　\tag{3}$

同時，當 $n$ 足夠大時，我們可以近似地認為 $\frac{k}{n}$ 可以作為 $P$ 的一個近似值。

然后，假設 $n$ 足夠大， $R$ 足夠小，並且假設 $p(x)$ 是連續的，那么我們可以得到：

$\int_{R} P(x^{'}) dx^{'} \approx p(x) \cdot V 　\tag{4}$

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

$p(x)V = \int_{R} p(x^{'}) dx^{'} = P = \frac{E(k)}{n} \tag{5}$
$\therefore p(x) = \frac{E(k)}{n \cdot V} \approx \frac{k}{n \cdot V}　\tag{6}$

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

問題

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

$p(x) \approx \frac{k/n}{V} \approx \frac{P}{V} = \frac{\int_{R} p(x^{'})d{x^{'}}}{\int_{R}dx{'}}　　\tag{7}$

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

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

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

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

$p_n(x) = \frac{k_n / n}{V_n} \tag{8}$

其中 $ｐ_n(x)$ 表示第 $n$ 次估計結果，如果要求 $p_n(x)$ 能夠收斂到 $p(x)$ ,則需要滿足下面三個條件：

$\lim_{n \rightarrow \infty} V_n = 0$
$\lim_{n \rightarrow \infty} k_n = \infty$
$\lim_{n \rightarrow \infty} k_n / n = 0$

parzen窗法

原理

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

$V_n = h_n^d \tag{9}$

然后我們再定義一個窗函數（window function）:

$\varphi (u) = \begin{cases} 1 \quad |u_j| \leq 1/2 \\ 0 \quad otherwise \end{cases} \quad \tag{10}$

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

$k_n = \sum_{i=1}^{n} \varphi(\frac{x - x_i}{h_n}) \tag{11}$

好了,有了 $k_n$ 和 $V_n$ ,直接把他們帶入公式(6),我們可以得到parzen窗法的計算公式:
$p_n(x) = \frac{1}{n} \sum_{i=1}^{n} \frac{1}{V_n} \varphi(\frac{x - x_i}{h_n}) \tag{12}$

我們發現這個 $\varphi$ 不僅可以是上述的單位超立方體的形式,只要其滿足如下兩個約束就可以,因此也就出現了各種各樣更能表現樣本屬性的窗函數,比如用的非常多的高斯窗.

$\varphi(x) \geq 0$

$\int \varphi(u)du = 1$

解釋

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

圖2

上面一句話可以這樣解釋:
我們定義核函數:
$\delta(x) = \frac{1}{V_n} \varphi(\frac{x}{h_n}) \tag{13}$

那么某一點 $x$ 的概率密度可以用如下函數來表示:
$p_n(x) = \frac{1}{n} \sum_{i=1}^{n} \delta_{n}(x - x_i) \tag{14}$

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

仿真實驗

如下我做了兩個仿真

第一個是生成了均值是0，方差是2的服從高斯分布的數據，分別使用bandwidth為0.1, 1, 5三個值進行估計
第二個是生成了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.] 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

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 【686】核密度估計，kernel density estimation，KDE說明核密度估計 Kernel Density Estimation (KDE) MATLAB 非參數估計——核密度估計（Parzen窗）非參數估計：核密度估計KDE 概率密度估計筆記——非參數估計概率密度估計--參數估計與非參數估計非參數密度估計 R語言與非參數統計（核密度估計）核密度估計(parzen窗密度估計) 參數估計&非參數估計