核密度估計，或Parzen窗，是非參數估計概率密度的一種。比如機器學習中還有K近鄰法也是非參估計的一種，不過K近鄰通常是用來判別樣本類別的，就是把樣本空間每個點劃分為與其最接近的K個訓練抽樣中，占比最高的類別。

直方圖

　　首先從直方圖切入。對于隨機變量$X$的一組抽樣，即使$X$的值是連續的，我們也可以劃分出若干寬度相同的區間，統計這組樣本在各個區間的頻率，并畫出直方圖。下圖是均值為0，方差為2.5的正態分布。從分布中分別抽樣了100000和10000個樣本：

　　這里的直方圖離散地取了21個相互無交集的區間：$[x-0.5,x+0.5), x=-10,-9,...,10$，單邊間隔$h=0.5$。$h>0$在核函數估計中通常稱作帶寬，或窗口。每個長條的面積就是樣本在這個區間內的頻率。如果用頻率當做概率，則面積除以區間寬度后的高，就是擬合出的在這個區間內的平均概率密度。因為這里取的區間寬度是1，所以高與面積在數值上相同，使得長條的頂端正好與密度函數曲線相契合。如果將區間中的$x$取成任意值，就可以擬合出實數域內的概率密度（其中$N_x$為樣本$x_iin [x-h,x+h),i=1,...,N$的樣本數）：

$displaystylehat{f}(x)=frac{N_x}{N}cdotfrac{1}{2h}$

　　這就已經是核函數估計的一種了。顯然，抽樣越多，這個平均概率密度能擬合得越好，正如藍條中上方幾乎都與曲線契合，而橙色則稂莠不齊。另外，如果抽樣數$N o infty$，對$h$取極限$h o 0$，擬合出的概率密度應該會更接近真實概率密度。但是，由于抽樣的數量總是有限的，無限小的$h$將導致只有在抽樣點處，才有頻率$1/N$，而其它地方頻率全為0，所以$h$不能無限小。相反，$h$太大的話又不能有效地將抽樣量用起來。所以這兩者之間應該有一個最優的$h$，能充分利用抽樣來擬合概率密度曲線。容易推理出，$h$應該和抽樣量$N$有關，而且應該與$N$成反比。

核函數估計

　　為了便于拓展，將擬合概率密度的式子進行變換：

$displaystylehat{f}(x)=frac{N_x}{2hN} = frac{1}{hN}sumlimits_{i=1}^{N}egin{cases}1/2& x-hle x_i < x+h\ 0& else end{cases}$

$displaystyle= frac{1}{hN}sumlimits_{i=1}^{N}egin{cases} 1/2,& -1le displaystylefrac{x_i-x}{h} < 1\ 0,& else end{cases}$

$displaystyle= frac{1}{hN}sumlimits_{i=1}^{N}displaystyle K(frac{x_i-x}{h}),;; where ; K(x) =egin{cases} 1/2,& -1le x < 1\ 0,& else end{cases}$

　　得到的$K(x)$就是uniform核函數（也又叫方形窗口函數），這是最簡單最常用的核函數。形象地理解上式求和部分，就是樣本出現在$x$鄰域內部的加權頻數（因為除以了2，所以所謂“加權”）。核函數有很多，常見的還有高斯核函數（高斯窗口函數），即：

$displaystyle K(x) = frac{1}{sqrt{2pi}}e^{-x^2/2}, -infty<x<infty$

　　各種核函數如下圖所示：

核函數的條件

　　并不是所有函數都能作為核函數的，因為$hat{f}(x)$是概率密度，則它的積分應該為1，即：

$displaystyleintlimits_{R}hat{f}(x) dx = intlimits_{R}frac{1}{hN}sumlimits_{i=1}^{N} K(frac{x_i-x}{h})dx=frac{1}{hN}sumlimits_{i=1}^{N}int_{-infty}^{infty} K(frac{x_i-x}{h})dx$

　　令$displaystyle t = frac{x_i-x}{h}$

$displaystyle=frac{1}{N}sumlimits_{i=1}^{N}int_{infty}^{-infty} -K(t)dt$

$displaystyle=frac{1}{N}sumlimits_{i=1}^{N}int_{-infty}^{infty} K(t)dt=1$

　　因積分部分為定值，所以可得$K(x)$需要的條件是：

$displaystyleint_{-infty}^{infty} K(x)dx=1$

　　通常$K(x)$是偶函數，而且不能小于0，否則就不符合實際了。

帶寬選擇與核函數優劣

　　正如前面提到的，帶寬$h$的大小關系到擬合的精度。對于方形核函數，$N o infty$時，$h$通常取收斂速度小于$1/N$的值即可，如$h=1/sqrt{N}$。對于高斯核，有證明指出$displaystyle h=left ( frac{4 hat{sigma}^5 }{3N} ight )^{frac{1}{5}}$時，有較優的擬合效果（$hat{sigma}^2$是樣本方差）。具體的帶寬選擇還有更深入的算法，具體問題還是要具體分析，就先不細究了。使用高斯核時，待擬合的概率密度應該近似于高斯分布那樣連續平滑的分布，如果是像均勻分布那樣有明顯分塊的分布，擬合的效果會很差。我認為原因應該是它將離得很遠的樣本也用于擬合，導致本該突兀的地方都被均勻化了。

　　Epanechnikov在均方誤差的意義下擬合效果是最好的。這也很符合直覺，越接近$x$的樣本的權重本應該越高，而且超出帶寬的樣本權重直接為0也是符合常理的，它融合了均勻核與高斯核的優點。

多維情況

　　對于多維情況，假設隨機變量$X$為$m$維（即$m$維向量），則擬合概率密度是$m$維的聯合概率密度：

$displaystyle hat{f}(x)= frac{1}{h^mN}sumlimits_{i=1}^{N}displaystyle K(frac{x_i-x}{h})$

　　其中的$K(x)$也變成了$m$維的標準聯合概率密度。另外，既然$displaystylefrac{1}{N}sumlimits_{i=1}^{N} K(frac{x_i-x}{h})$代表的是概率，$m$維的概率密度自然是概率除以$h^m$而不是$h$。

實驗擬合情況

　　分別取帶寬$h=0.1,0.3,0.7,1.4$時，使用三種核函數對分布$displaystyle p_X(x) = frac{-x+5}{50},xin [-5,5]$進行擬合：

　　抽樣數量$N=100000$，可以看出隨著$h$增大，偏差增大，而$h$太小時，方差變大了。可以發現高斯核的擬合從來都是光滑的（方差比較小），這樣看起來似乎在$h$取得很小時，高斯核是比較好的核函數。而當$h$因為抽樣較少而不得不取大時，另外兩個核函數則更能勾勒出待擬合函數的輪廓。

　　以下是實驗代碼：

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 from scipy.special import comb,perm
 4 
 5 sample_num = 100000
 6 #獲取要擬合的分布抽樣并排序 Y = 5-10*(1-X)**0.5 
 7 ran = np.random.rand(sample_num)
 8 ran = 5-10*(1-ran)**0.5 
 9 ran = np.sort(ran) 
10 
11 #高斯核
12 def ker_gass(x0):
13     return (1/(2*np.pi)**0.5)*np.e**-(x0**2/2)
14 #Epanechnikov核
15 def ker_Epanechnikov(x0):
16     return 3/4*(1-x0**2)
17 
18 #擬合概率密度函數
19 def fitting_proba_density(X,h,way): 
20     if way == 1:    #使用均勻核
21         i_X = 0
22         begin_ran = 0
23         end_ran = 0
24         sum0 = np.zeros(len(X)+1)
25         while i_X < len(X):  
26             while begin_ran < sample_num:
27                 if X[i_X] - h > ran[begin_ran]:
28                     sum0[i_X] -= 0.5
29                     begin_ran+=1 
30                 else:
31                     break
32             while end_ran < sample_num:
33                 if X[i_X] + h >= ran[end_ran]:
34                     sum0[i_X]+=0.5
35                     end_ran+=1
36                 else:
37                     break
38             sum0[i_X+1] = sum0[i_X] 
39             i_X+=1
40         return sum0[0:-1]/h/sample_num 
41     elif way == 2:    #使用高斯核 
42         sum0 = np.zeros(len(X))
43         for i in range(sample_num):
44             sum0 += ker_gass((ran[i]-X)/h)
45         return sum0/h/sample_num
46     else:    #使用Epanechnikov核 
47         i_X = 0
48         begin_ran = 0
49         end_ran = 0
50         sum0 = np.zeros(len(X))
51         while i_X < len(X):  
52             while begin_ran < sample_num:
53                 if X[i_X] - h > ran[begin_ran]: 
54                     begin_ran+=1 
55                 else:
56                     break
57             while end_ran < sample_num:
58                 if X[i_X] + h >= ran[end_ran]: 
59                     end_ran+=1
60                 else:
61                     break
62             i = begin_ran
63             while i < end_ran:
64                 sum0[i_X] += ker_Epanechnikov((ran[i]-X[i_X])/h)
65                 i+=1
66             i_X+=1
67         return sum0/h/sample_num 
68 
69 #畫出擬合概率密度
70 def paint_(a):
71     X = np.linspace(-10,10,500)
72     j=0
73     for h in a: 
74         j+=1
75         ax = plt.subplot(2,2,j)
76         ax.set_title('h='+ str(h))#設置子圖
77 
78         X0 = np.linspace(-5,5,10)
79         Y0 = (-X0+5)/50
80         plt.plot(X0,Y0,label = 'Probability density')#分布密度函數 
81 
82         Y = fitting_proba_density(X,h,1)#均勻核
83         ax.plot(X,Y,label = 'Uniform kernel') 
84         Y = fitting_proba_density(X,h,2)#高斯核
85         ax.plot(X,Y,label = 'Gassian kernel')
86         Y = fitting_proba_density(X,h,3)#Epanechnikov核
87         ax.plot(X,Y,label = 'Epanechnikov kernel')
88         ax.legend()
89 paint_([0.1,0.3,0.7,1.4])
90  
91 #圖像參數
92 plt.xlim(-10,10)
93 plt.show()

總結

以上是生活随笔為你收集整理的非参数估计——核密度估计（Parzen窗）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。