模式識別七--非引數估計法之Parzen窗估計和k
文章轉自：http://www.kancloud.cn/digest/prandmethod/102849
本實驗的目的是學習Parzen窗估計和k最近鄰估計方法。在之前的模式識別研究中，我們假設概率密度函式的引數形式已知，即判別函式J(.)的引數是已知的。本節使用非引數化的方法來處理任意形式的概率分佈而不必事先考慮概率密度的引數形式。在模式識別中有躲在令人感興趣的非引數化方法，Parzen窗估計和k最近鄰估計就是兩種經典的估計法。
參考書籍：《模式分類》
作者：RichardO.Duda，PeterE.Hart，DavidG.Stork
一、基本原理
1.非引數化概率密度的估計
對於未知概率密度函式的估計方法，其核心思想是：一個向量x落在區域R中的概率可表示為：
其中，P是概率密度函式p(x)的平滑版本，因此可以通過計算P來估計概率密度函式p(x)，假設n個樣本x1,x2,…,xn，是根據概率密度函式p(x)獨立同分布的抽取得到，這樣，有k個樣本落在區域R中的概率服從以下分佈：
其中k的期望值為：
k的分佈在均值附近有著非常顯著的波峰，因此若樣本個數n足夠大時，使用k/n作為概率P的一個估計將非常準確。假設p(x)是連續的，且區域R足夠小，則有：
如下圖所示，以上公式產生一個特定值的相對概率，當n趨近於無窮大時，曲線的形狀逼近一個δ函式，該函式即是真實的概率。公式中的V是區域R所包含的體積。綜上所述，可以得到關於概率密度函式p(x)的估計為：
在實際中，為了估計x處的概率密度函式，需要構造包含點x的區域R1,R2,…,Rn。第一個區域使用1個樣本，第二個區域使用2個樣本，以此類推。記Vn為Rn的體積。kn為落在區間Rn中的樣本個數，而pn (x)表示為對p(x)的第n次估計：
欲滿足pn(x)收斂：pn(x)→p(x)，需要滿足以下三個條件：
有兩種經常採用的獲得這種區域序列的途徑，如下圖所示。其中“Parzen窗方法”就是根據某一個確定的體積函式，比如Vn=1/√n來逐漸收縮一個給定的初始區間。這就要求隨機變數kn和kn/n能夠保證pn (x)能收斂到p(x)。第二種“k-近鄰法”則是先確定kn為n的某個函式，如kn=√n。這樣，體積需要逐漸生長，直到最後能包含進x的kn個相鄰點。
2.Parzen窗估計法
已知測試樣本資料x1,x2,…,xn，在不利用有關資料分佈的先驗知識，對資料分佈不附加任何假定的前提下，假設R是以x為中心的超立方體，h為這個超立方體的邊長，對於二維情況，方形中有面積V=h^2，在三維情況中立方體體積V=h^3，如下圖所示。
根據以下公式，表示x是否落入超立方體區域中：
估計它的概率分佈：
其中n為樣本數量，h為選擇的窗的長度，φ(.)為核函式，通常採用矩形窗和高斯窗。
3.k最近鄰估計
在Parzen演算法中，窗函式的選擇往往是個需要權衡的問題，k-最近鄰演算法提供了一種解決方法，是一種非常經典的非引數估計法。基本思路是：已知訓練樣本資料x1,x2,…,xn而估計p(x)，以點x為中心，不斷擴大體積Vn，直到區域內包含k個樣本點為止，其中k是關於n的某一個特定函式，這些樣本被稱為點x的k個最近鄰點。
當涉及到鄰點時，通常需要計算觀測點間的距離或其他的相似性度量，這些度量能夠根據自變數得出。這裡我們選用最常見的距離度量方法：歐幾里德距離。
最簡單的情況是當k=1的情況，這時我們發現觀測點就是最近的（最近鄰）。一個顯著的事實是：這是簡單的、直觀的、有力的分類方法，尤其當我們的訓練集中觀測點的數目n很大的時候。可以證明，k最近鄰估計的誤分概率不高於當知道每個類的精確概率密度函式時誤分概率的兩倍。
二、實驗基本步驟
第一部分，對錶格中的資料，進行Parzen 窗估計和設計分類器，本實驗的窗函式為一個球形的高斯函式，如下：
1) 編寫程式，使用Parzen 窗估計方法對一個任意的測試樣本點x 進行分類。對分類器的訓練則使用表格 3中的三維資料。同時，令h =1，分類樣本點為(0.5,1.0,0.0)，(0.31,1.51,-0.50)，(-0.3,0.44,-0.1)進行實驗。
2) 可以改變h的值，不同的h將導致不同的概率密度曲線，如下圖所示。
h=0.1時：
h=0.5時：
h=1時：
第二部分的實驗目的是學習和掌握非引數估計：k-近鄰概率密度估計方法。對前面表格中的資料進行k-近鄰概率密度估計方法和設計分類器。
編寫程式，對錶格中的3個類別的三維特徵，使用k-近鄰概率密度估計方法。並且對下列點處的概率密度進行估計：(-0.41,0.82,0.88)，(0.14,0.72, 4.1) ，(-0.81,0.61,-0.38)。
實驗程式碼如下：
% Parzen窗演算法
% w：c類訓練樣本
% x：測試樣本
% h：引數
% 輸出p：測試樣本x落在每個類的概率
function p = Parzen(w,x,h)
[xt,yt,zt] = size(w);
p = zeros(1,zt);
for i = 1:zt
hn = h;
for j = 1:xt
hn = hn / sqrt(j);
p(i) = p(i) + exp(-(x - w(j,:,i))*(x - w(j,:,i))'/ (2 * power(hn,2))) / (hn * sqrt(2*3.14));
end
p(i) = p(i) / xt;
end% k-最近鄰演算法
% w：c類訓練樣本
% x：測試樣本
% k：引數
function p = kNearestNeighbor(w,k,x)
% w = [w(:,:,1);w(:,:,2);w(:,:,3)];
[xt,yt,zt] = size(w);
wt = [];%zeros(xt*zt, yt);
if nargin==2
p = zeros(1,zt);
for i = 1:xt
for j = 1:xt
dist(j,i) = norm(wt(i,:) - wt(j,:));
end
t(:,i) = sort(dist(:,i));
m(:,i) = find(dist(:,i) <= t(k+1,i)); % 找到k個最近鄰的編號
end
end
if nargin==3
for q = 1:zt
wt = [wt; w(:,:,q)];
[xt,yt] = size(wt);
end
for i = 1:xt
dist(i) = norm(x - wt(i,:));
end
t = sort(dist); % 歐氏距離排序
[a,b] = size(t);
m = find(dist <= t(k+1)); % 找到k個最近鄰的編號
num1 = length(find(m>0 & m<11));
num2 = length(find(m>10 & m<21));
num3 = length(find(m>20 & m<31));
if yt == 3
plot3(w(:,1,1),w(:,2,1),w(:,3,1), 'r.');
hold on;
grid on;
plot3(w(:,1,2),w(:,2,2),w(:,3,2), 'g.');
plot3(w(:,1,3),w(:,2,3),w(:,3,3), 'b.');
if (num1 > num2) || (num1 > num3)
plot3(x(1,1),x(1,2),x(1,3), 'ro');
disp(['點：[',num2str(x),']屬於第一類']);
elseif (num2 > num1) || (num2 > num3)
plot3(x(1,1),x(1,2),x(1,3), 'go');
disp(['點：[',num2str(x),']屬於第二類']);
elseif (num3 > num1) || (num3 > num2)
plot3(x(1,1),x(1,2),x(1,3), 'bo');
disp(['點：[',num2str(x),']屬於第三類']);
else
disp('無法分類');
end
end
if yt == 2
plot(w(:,1,1),w(:,2,1), 'r.');
hold on;
grid on;
plot(w(:,1,2),w(:,2,2), 'g.');
plot(w(:,1,3),w(:,2,3), 'b.');
if (num1 > num2) || (num1 > num3)
plot(x(1,1),x(1,2), 'ro');
disp(['點：[',num2str(x),']屬於第一類']);
elseif (num2 > num1) || (num2 > num3)
plot(x(1,1),x(1,2), 'go');
disp(['點：[',num2str(x),']屬於第二類']);
elseif (num3 > num1) || (num3 > num2)
plot(x(1,1),x(1,2), 'bo');
disp(['點：[',num2str(x),']屬於第三類']);
else
disp('無法分類');
end
end
end
title('k-最近鄰分類器');
legend('第一類資料',...
'第二類資料',...
'第三類資料',...
'測試樣本點');clear;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parzen窗估計和k最近鄰估計
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w1(:,:,1) = [ 0.28 1.31 -6.2;...
0.07 0.58 -0.78;...
1.54 2.01 -1.63;...
-0.44 1.18 -4.32;...
-0.81 0.21 5.73;...
1.52 3.16 2.77;...
2.20 2.42 -0.19;...
0.91 1.94 6.21;...
0.65 1.93 4.38;...
-0.26 0.82 -0.96];
w1(:,:,2) = [0.011 1.03 -0.21;...
1.27 1.28 0.08;...
0.13 3.12 0.16;...
-0.21 1.23 -0.11;...
-2.18 1.39 -0.19;...
0.34 1.96 -0.16;...
-1.38 0.94 0.45;...
-0.12 0.82 0.17;...
-1.44 2.31 0.14;...
0.26 1.94 0.08];
w1(:,:,3) = [ 1.36 2.17 0.14;...
1.41 1.45 -0.38;...
1.22 0.99 0.69;...
2.46 2.19 1.31;...
0.68 0.79 0.87;...
2.51 3.22 1.35;...
0.60 2.44 0.92;...
0.64 0.13 0.97;...
0.85 0.58 0.99;...
0.66 0.51 0.88];
x(1,:) = [0.5 1 0];
x(2,:) = [0.31 1.51 -0.5];
x(3,:) = [-0.3 0.44 -0.1];
% 驗證h的二維資料
w2(:,:,1) = [ 0.28 1.31 ;...
0.07 0.58 ;...
1.54 2.01 ;...
-0.44 1.18 ;...
-0.81 0.21 ;...
1.52 3.16 ;...
2.20 2.42 ;...
0.91 1.94 ;...
0.65 1.93 ;...
-0.26 0.82 ];
w2(:,:,2) = [0.011 1.03 ;...
1.27 1.28 ;...
0.13 3.12 ;...
-0.21 1.23 ;...
-2.18 1.39 ;...
0.34 1.96 ;...
-1.38 0.94 ;...
-0.12 0.82 ;...
-1.44 2.31 ;...
0.26 1.94 ];
w2(:,:,3) = [1.36 2.17 ;...
1.41 1.45 ;...
1.22 0.99 ;...
2.46 2.19 ;...
0.68 0.79 ;...
2.51 3.22 ;...
0.60 2.44 ;...
0.64 0.13 ;...
0.85 0.58 ;...
0.66 0.51 ];
y(1,:) = [0.5 1];
y(2,:) = [0.31 1.51];
y(3,:) = [-0.3 0.44];
h = .1; % 重要引數
p = Parzen(w1,x(1,:),h);
num = find(p == max(p));
disp(['點：[',num2str(x(1,:)),']落在三個類別的概率分別為：',num2str(p)]);
disp(['點：[',num2str(x(1,:)),']落在第',num2str(num),'類']);
% 給定三類二維樣本，畫出二維正態概率密度曲面圖驗證h的作用
num =1; % 第num類的二維正態概率密度曲面圖，取值為1，2，3
draw(w2,h,num);
str1='當h=';str2=num2str(h);str3='時的二維正態概率密度曲面';
SC = [str1,str2,str3];
title(SC);
% k近鄰演算法設計的分類器
% x1和y1為測試樣本
x1 = [-0.41,0.82,0.88];
x2 = [0.14,0.72, 4.1];
x3 = [-0.81,0.61,-0.38];
y(1,:) = [0.5 1];
y(2,:) = [0.31 1.51];
y(3,:) = [-0.3 0.44];
w = w1;
%w = w1(:,1,3);
k = 5;
kNearestNeighbor(w,k,x1);
kNearestNeighbor(w,k,x2);
kNearestNeighbor(w,k,x3);