Parzen窗估计和K最近邻估计

非参数概率密度估计

概率密度估计即从样本中估计总体的概率密度分布。进行概率密度估计的方法有参数估计与非参数估计。参数估计即在已知总体的分布函数形式的情况下对分布函数参数进行估计。如假如已知一总体呈正态分布,在取得样本的情况下,对总体的概率密度分布的估计就变成对总体的$\mu$ 和 $\sigma^2$ 进行估计。而在未知总体的分布函数形式的情况下,我们可以通过非参数概率密度估计方法对概率密度进行估计。这篇文章会讨论非参数估计的两种方法:Parzen窗估计和K最近邻估计,并对其中的概率密度估计和预测分类过程进行Octave语言实现。

非参数概率密度估计的思想非常直接,简单来说即是将观测到事件的概率值来估计事件发生的概率再推测基本事件的概率值(从样本满足独立同分布,计算二项分布概率值再求最佳估计值得来)。假设事件发生频率 $P(A) = k / n$, 事件中元素个数(二维空间下的面积)乘以每个元素发生的概率(假设每个都相等)为 $V\times p(x)$,假设两者相等即有 $p(x) = \frac{k/n}{V}$。

当然,假如直接使用 $T = \frac{k/n}{V}$ 来估计概率密度的话是十分粗糙的,但是假如不断调整V(使其趋近0)或k的值使得 $p(x)$ 趋近一个稳定值的话那么便可以用 $T$ 来估计 $p(x)$。根据调整方法的不同可以分为 Parzen 窗口法和 K最近邻估计法。

Parzen窗估计

Parzen窗估计的核心思想是按一定规律不断缩小$V$使得 $\frac{k/n}{V}$ 趋近于真实的$p(x)$。Parzen窗泛化后可以使用其他核函数来取代事件的计数频率值。

首先来形式化以上语言。对于有$d$个纬度多维空间$R_d$,假设子空间$V_n = h_n^d$,其中n是迭代计数,$h$为$V$在坐标系中每个纬度投影的长度。假设希望求在空间中$x$坐标下的概率密度,则包含$x$且大小为$V_n$的空间中的向量$y$要满足$|x_i - y_i|/h_n\leq 1/2, i=1,2,..,d$(既$y$要位于$x$周围的一个超立方体内)。按此规则可以表达为示性函数:

$$\varphi (\mu) = \begin{cases} 1 \ \ \ \ & if |\mu_j|\leq\frac{1}{2} \ \ j=1,...,d \\ 0 \ \ \ \ & otherwise \end{cases}$$

其中$\mu_j = |x_i - y_i|/h_n$, 则对于每个样本$x_i$,是否落入$x$附近就表达为$\varphi (\frac{x-x_i}{h_n})$是否为$1$,这样$P_n(A) = k/n = \frac{\sum\varphi (\frac{x-x_i}{h_n}) }{n}$,令其等于$V_n \times p_n(x)$,即有$p_n = \frac{\sum\varphi (\frac{x-x_i}{h_n}) }{n}/V_n$。

上式中$\varphi (\mu)$可以视为核函数,这样其便能使用其他核函数进行替代,例如使用高斯核函数。

Parzen窗函数实现

原始的Parzen窗口函数即超立方体窗口是作用在每个样本点上的示性函数,对其进行推广转变为高斯窗口的形式为$\varphi_{\sigma} (\frac{x-x_i}{h_n}) = norm( (\frac{x-x_i}{h_n}), \mu, \sum)$,其中$norm$为多元高斯概率分布函数,其中$x$的维度为$d$,且$\mu=(0,…,0)^t, \sum=I_d$ 。

根据窗口函数可以计算出样本中事件发生的次数$k = \sum\varphi (\frac{x-x_i}{h_n}) $。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 超立方体窗口
# data 为 n*d, point 为1*d
function [z] = w_hypercube(data, point, h_n)
u = (data - point) ./ h_n;
z = sum(all((abs(u) <= 0.5), 2));
end
# 高斯窗口
pkg load statistics;
function [z] = w_normal(data, point, h_n)
[n, d] = size(data);
u = (data - point) ./ h_n;
z = sum(mvnpdf(u, zeros(1, d), eye(d)));
end

有了$k$值,就可计算出概率密度,即有概率密度估计方程 $p_n = \frac{k}{n}/V_n = \frac{\sum\varphi (\frac{x-x_i}{h_n}) }{n}/V_n$。

1
2
3
4
5
6
7
8
9
10
# Parzen 函数
#data is [n * d], dimension X number
#point is [1 * d]
#h_n is a scale
function [pd] = parzen (data, point, h_n, f)
[n_data, d_data] = size(data);
pd = f(data, point, h_n) / (n_data * (h_n^d_data));
end
## parzen(XX, [3, 4], 1, @w_hypercube)

两类实验数据

使用Parzn窗估计概率密度需要基于样本数据进行。下面先生成1000个二维样本。这些样本可以分成两类,每个类别分别服从独立高斯分布。我们可以先用这些数据来讨论Parzen窗进行概率估计的效果(使用所有数据进行),再基于此数据集讨论使用Parzen窗估计进行分类的效果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 生成二维高斯分布数据函数
function [x] = create_norm(mu1, sigma1, mu2, sigma2, n)
x = [normrnd(mu1, sigma1, 1, n); normrnd(mu2, sigma2, 1, n)]';
end
# X1为第一类,X2为第二类
# X1为服从正态分布且中心点在(5,4),各分量sigma为(1,0.5)的数据,样本数量为500
# 两类数据真实分布有重叠,如下图
[mu11, s11, mu12, s12, n1] = deal(5, 1, 4, 0.5, 500);
[mu21, s21, mu22, s22, n2] = deal(8, 0.5, 4, 1, 500);
X1 = create_norm(mu11, s11, mu12, s12, n1);
X2 = create_norm(mu21, s21, mu22, s22, n2);
# 合并数据
XX = [X1;X2];
# 画出数据分布
plot(X1(:,1), X1(:,2), 'r+', X2(:,1), X2(:,2), 'bo', ...
mu11, mu12, '+k', mu21, mu22, '+k');

Parzen_数据分布

Parzen窗概率密度估计测试

下面先将两类数据合为总数据集,查看Parzen概率密度估计的效果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# parzen函数测试,输入数据集,窗口宽度及窗口函数,输出概率密度分布平面图
function [z] = test_parzen(XX, h, f)
#x,y 为构造的平面坐标
x = 0:0.1:10;
y = x;
[X,Y] = meshgrid(x);
n = size(x, 2);
Z = zeros(n, n);
for i=1:n
for j=1:n
#对平面中的每个点,计算其概率密度值
Z(i,j) = parzen(XX, [x(i), y(j)], h, f);
end
end
#Z_scale = (Z-min(Z(:))) ./ (max(Z(:)-min(Z(:))));
# 画平面图
surf(X,Y,Z);
end
# 窗函数可用w_normal(高斯), w_hypercube(超立方体窗口)
# 宽度h取 0.1 ~ 2 进行测试
f = @w_hypercube;
h = 0.1;
test_parzen(XX, h, f);

超立方体窗口(样本量为1000):

超立方体窗口

高斯窗口(样本量为1000):

parzen_norm

从图中可看出Parzen窗估计能有效得还原出原始数据的分布。窗口函数和窗口宽度的选择都会影响给出的概率密度估计值。高斯窗比立方体窗平滑,并可以看出随着窗口宽度的减小(从2降到0.5),估计的概率值峰值增大(高斯窗的峰值从0.02增大到0.18,立体窗的峰值从0.1增大到0.5),且分布的广度变小。

下面看看不同的样本容量对概率密度估计的影响。

1
2
3
4
5
6
7
8
# 窗函数可用w_normal(高斯)
# 宽度h取 0.5 进行测试
# 样本容量取 100, 200, 500, 1000进行测试
f = @w_normal;
h = 0.5;
n = 100;
sub_XX = [X1(1:(n/2),:);X2(1:(n/2),:)];
test_parzen(sub_XX, h, f);

parzen_n

可以看出当样本量较少时(上图n=200的情况)概率密度的估计分布与原分布差异较大。增大样本量能提高概率密度估计分布与真实分布的显示度。另一方面,从图中可以看出改变样本量对概率分布峰值和分布广度的影响较小,图中不同的样本量下概率值峰值都在0.1左右。

基于Parzen窗估计的分类

由于Parzen窗估计能估计出概率密度值,有了概率值就能用来进行分类预测。实际应用Parzen窗估计时可以分别在每个类别下进行概率密度估计,即得到每个类别的$P(x|w_i)$分布。一般地,在进行分类的时候,只要能根据观测值$x$计算出$P(w_i|x)$便可以在各种类别中进行选择最佳类别(类别为$w_1, w_2, …, w_c$)。根据贝叶斯公式$P(w_i|X) \propto P(w_i) \times P(X|w_i)$,我们可以通过Parzen窗口估计出不同类别下的$P(x|w_i)$值,再结合先验概率$P(w_i)$就可知道后验$P(w_i|x)$概率值,从而对样本$x$进行分类。下面根据数据的类别标签将数据分为两个数据集,再分别在两个类别下进行概率密度估计,进而可以得到分类界面。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
function [z] = train_parzen(X1, X2, h, f)
x = 0:0.1:10;
y = x;
[X,Y] = meshgrid(x);
n_1 = size(X1, 1);
n_2 = size(X2, 1);
n = size(x, 2);
Z = zeros(n, n);
for i=1:n
for j=1:n
Z(i, j) = parzen(X1, [x(i), y(j)], h, f)*n_1 > ...
parzen(X2, [x(i), y(j)], h, f)*n_2;
end
end
#Z_scale = (Z-min(Z(:))) ./ (max(Z(:)-min(Z(:))));
plot(X(Z(:)==0), Y(Z(:)==0), 'bo', ...
X(Z(:)==1), Y(Z(:)==1), 'r.', ...
5, 4, '+k', 8, 4, '+k');
end
#plot(X(Z(:)==0), Y(Z(:)==0),'.r','markersize', 20);
#hold on
#plot(X(Z(:)==1), Y(Z(:)==1),'^b', 'markersize', 7);
#f = @w_normal;
f = @w_hypercube;
h = 0.5;
res = train_parzen(X1, X2, h, f);

在窗口宽度为0.5时,在不同的窗口函数下分类界面如下:

Parzen_分类界面

从图中可看出不同的窗口函数下分类界面存在明显差异。直接用Parzen窗得到的分类效果在泛化时效果不好,因为很多边缘区域并不能很好的分类(若考虑拒识的情况则会有大量空白),这是因为Parzen窗的到的概率值完全取决于其周边数据的分布,若周围样本稀疏,Parzen窗函数并不能得到一个合适的概率值。进一步可以看出在两个类别分布类似且有重叠时用超立方函数将两类分开较好。

实际测试

下面由原分布生成数据进行测试,比较在不同的窗口函数和窗口宽度下分类的正确率。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# 生成测试数据
Y1 = create_norm(mu11, s11, mu12, s12, 200);
Y2 = create_norm(mu21, s21, mu22, s22, 200);
Y = [Y1;Y2];
class_Y = [ones(200, 1); zeros(200, 1)];
# 预测类别
function [z] = predict_parzen(X1, X2, x, h, f)
n = size(x, 1);
P = zeros(n, 1);
n_1 = size(X1, 1);
n_2 = size(X2, 1);
for i=1:n
point = [x(i, 1), x(i, 2)];
# 1 if predict to class X1
P(i) = (parzen(X1, point, h, f)*n_1 > ...
parzen(X2, point, h, f)*n_2);
end
#Z_scale = (Z-min(Z(:))) ./ (max(Z(:)-min(Z(:))));
z = P;
end
function [z] = pre_rate(Y, pre_Y)
z = sum(Y == pre_Y) / size(Y, 1);
end
#f = @w_hypercube;
f = @w_normal;
h = 0.5;
pre_Y = predict_parzen(X1, X2, Y, h, f);
pre_rate(class_Y, pre_Y)

测试结果:

窗口函数 窗口宽度 预测率/识别率
normal 1 0.987
超立方体 1 0.985
normal 0.5 0.987
超立方体 0.5 0.972

从结果来看使用立方体窗口的预测效果与高斯窗口的情况差别并不大。两种窗口对样本的分类效果都较好$\geq 95\%$,这点单从分类界面是看不出来的。

K最近邻估计

Parzen窗口法和K最近邻(The k-nearest neighbor,,KNN)的目的都是计算$p_n = \frac{k}{n}/V_n$,不同的是Parzen窗口选择在迭代的过程不断固定$V_n$来计算$p_n$,而K最近邻估计选择不断固定一个$k$值来计算$p_n$。K最近邻的计算过程简单来说是对于固定的$k$值,在待测空间点附近找到一个区域$V_n$使得落入该区域的样本点刚好有$k$个再计算出$P(x)$值。

在使用概率密度来进行分类时,对于一个固定的$k$,找到一个$V_n$(将所有类别的样本合为总体)之后下一步的计算时找到最大的$P(w_i|x) \propto P(w_i) \times P(x|w_i)$,对于恒定$P(w_i)$ 的情况,即找到最大的$P(x|w_i) = \frac{k_i}{n}/V_n$,又由于此时$n, V_n$(总样本量,固定k值下找到的空间)对于所有$i$(每个类别)都是固定的,因此只需找到类别$i$使得其落入在空间$V_n$次数最多。于是,K最近邻的分类就变成在待测空间找到相邻的$k$个样本点,这些样本点中类别出现次数最多的类别$i$就会有最大的$P(w_i|x)$。

由于K最近邻计算中有一个找$V_n$的过程,即找在待测空间点附近每个维度距离小于$h_n/2$的区域,在这个过程中距离的度量除了使用欧几里得距离外还可以使用其他距离度量,如曼哈顿距离等。

在二维空间下,对每个样本点都赋不同的类别标签,在$k=1$下画得的分类界面图即是Voronoi图。下图为在欧几里得距离下的Voronoi图(来自wikipedia)。

Euclidean_Voronoi_diagram

KNN概率密度估计

KNN的实现过程有一个固定$k$找邻域的过程,这是耗时的过程,但是已经有了一些优化的方法。下面用最直接的思想找领域再计算概率密度值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
# 定义距离度量, 欧几里得距离,
#data is [n * d], dimension X number
#point is [1 * d]
function [d] = dist(data, point)
d = sqrt(sum((data - point).^2, 2));
end
# KNN pdf 函数, k > 0
function [pd] = KNN_pdf (data, point, k, dis_f)
#compute all data distance to point
#sqrt((x1-y1)^2+(x1-y2)^2)
all_dist = dis_f(data, point);
k_dist = sort(all_dist)(k);
[n_data, d_data] = size(data);
pd = k / (n_data * ((k_dist * 2)^d_data));
end
# KNN函数测试,输入数据集,k值及距离度量函数,输出概率密度分布平面图
function [z] = test_KNN(XX, k, dis_f)
#x,y 为构造的平面坐标
x = 0:0.1:10;
y = x;
[X,Y] = meshgrid(x);
n = size(x, 2);
Z = zeros(n, n);
for i=1:n
for j=1:n
#对平面中的每个点,计算其概率密度值
Z(i,j) = KNN_pdf(XX, [x(i), y(j)], k, dis_f);
end
end
#Z_scale = (Z-min(Z(:))) ./ (max(Z(:)-min(Z(:))));
# 画平面图
surf(X,Y,Z);
end
#
f = @dist;
k = 100;
test_KNN(XX, k, f);

使用Parzen窗口中使用的数据集($n=1000$)进行测试,可以得到以下概率密度分布图。

KNN_pdf

可以看出随着$k$的增大,概率密度峰值减小,且双峰分布较为明显。下面看看使用KNN分类的情况。

KNN分类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
# KNN 分类器,data(n*d)为数据特征,Y(n*1)为数据类别, c为类别数量
# 要求类别编号从1开始
# x 为待分类样本,dis_f为距离度量函数
# 输出predict_Y(n*1)
function [z] = predict_KNN(data, Y, c, x, k, dis_f)
n = size(x, 1);
P = zeros(n, 1);
for i = 1:n
#对每个待分类样本计算距离
point = [x(i, 1), x(i, 2)];
all_dist = dis_f(data, point);
#找出最近k个样本点所属的类别
# sort all distance
[a, b] = sort(all_dist);
#get first k-th point class
first_k = zeros(k, 1);
for j = 1:k
first_k(j) = Y(b(j));
end
#这k个类别总出现次数最多的即为预测类别
#count and get max hit class
count = zeros(c,1);
for j = 1:c
count(j) = sum(first_k==j);
end
# predict class
[a, b] = sort(count);
P(i) = b(c);
end
z = P;
end
function [z] = pre_rate(Y, pre_Y)
z = sum(Y == pre_Y) / size(Y, 1);
end
data = XX;
data_true = [ones(500, 1); zeros(500, 1)] + 1;
c = 2;
f = @dist;
y = Y;
y_true = class_Y + 1;
k = 100;
#predict_KNN(data, data_true, c, [8, 4], k, f)
y_pred = predict_KNN(data, data_true, c, y, k, f);
pre_rate(y_true, y_pred)

对于1000个训练样本,200个测试样本的分类预测结果如下:

k 识别率
5 0.970
20 0.980
100 0.995

可以看出KNN的分类效果较好,在k只有5的情况下都可以获得较高的识别率($97\%$),随着k值增大预测识别率有上升。

总结

Parzen窗口法和KNN都能从样本中估计出总体的概率密度分布。Parzen窗口法通过调整$V_n$的范围来逼近真实的概率密度值,Parzen窗口法的实践依赖于窗口核函数和窗口宽度的选择。而KNN通过找出最近$k$样本划分得到的空间进行概率密度估计。KNN分类的实现可以绕过概率密度估计而直接求后验概率值。相比与Parzen窗口法,KNN不需要设定窗口宽度值,只需设定好距离度量和k值就可以进行分类,并给出较好的分类结果。

参考

  1. http://www.byclb.com/TR/Tutorials/neural_networks/ch11_1.htm
  2. http://sebastianraschka.com/Articles/2014_kernel_density_est.html