模式識別：PCA主分量分析與Fisher線性判別分析 · 模式分類理論與方法

本實驗的目的是學習和掌握PCA主分量分析方法和Fisher線性判別方法。首先了解PCA主分量分析方法的基本概念，理解利用PCA 分析可以對數據集合在特征空間進行平移和旋轉。實驗的第二部分是學習和掌握Fisher線性判別方法。了解Fisher線性判別方法找的最優方向與非最優方向的差異，將高維分布的數據進行降維，并通過Fisher線性判別方法實現高維數據在一維中分類。 ##一、技術論述 **1.統計分析方法中的降維思想** 在模式識別的研究過程中，往往需要對反映事物的多個變量進行大量的觀測，收集大量數據以便進行統計分析。多變量大樣本為研究者提供了豐富的信息，但也在一定程度上增加了統計分析的運算量。而在多數情況下，變量與變量之間存在著相關性，因而又增加了問題分析的復雜性。如果分別對每個指標進行分析，分析往往是孤立的，而不是綜合的。盲目減少數據會損失很多信息，容易產生錯誤的判決結果。因此需要找到一個合理的方法，在減少需要分析的信息和保證分析質量兩者中取得平衡。前面說到由于各變量間存在一定的相關關系，這使得這種思想成為可能。可以采用特征線性組合的方法減少特征空間中的特征維數。將高維數據投影到較低維空間上。主成分分析(principal component analysis)和Fisher判別分析法是兩類有效的線性組合變換方法。 **2.主成分分析PCA** 主成分分析的主要思想是對于原先提出的所有變量，建立盡可能少的新變量，使得這些新變量是兩兩不相關的，而且這些新變量在反映課題的信息方面盡可能保持原有的信息。主分量分析的目的是通過特征的線性組合降低特征空間的維數。假設一個由n個樣本組成的數據集合，每個樣本是一個d維特征矢量：{ x1,x2,…,xn }；希望使用一個單一的模式空間的點x0來表示整個樣本集合；并希望所選取的模式空間的點x0與數據集合中各個樣本點之間的距離都盡可能地小。即： ![](https://box.kancloud.cn/2016-01-05_568b368f794d1.jpg) 使均方誤差最小的矢量x_0是均值矢量m。樣本均值的優點是表達簡單，缺點是其不能表達集合中的樣本之間的差異，同時樣本集合的0維表示可以提供的信息很少。將樣本集合的全部樣本向通過均值矢量的一條直線作投影，可以得到樣本空間中表示樣本集合的一條直線，稱之為樣本集合的1維表達： ![](https://box.kancloud.cn/2016-01-05_568b368f91d1d.jpg) 其中，e表示通過均值的直線方向上的單位矢量；a是一個實數標量，表示直線的某一點x距離均值點m的距離。假設使用該直線上的一個點(m+ak*e)來表示樣本點xk，則可以構造樣本集合的平方誤差準則函數： ![](https://box.kancloud.cn/2016-01-05_568b368f9f62c.jpg) 求解準則函數的最優解： ![](https://box.kancloud.cn/2016-01-05_568b368fb0012.jpg) 其中||e||=1，對ak求偏導數并使結果為0，可以得到： ![](https://box.kancloud.cn/2016-01-05_568b368fbf7a5.jpg) 為向量(xk-m)在直線x=m+a*e方向上的投影長度。在PCA中，定義樣本集合的散布矩陣為： ![](https://box.kancloud.cn/2016-01-05_568b368fcc01a.jpg) 由上述定義可以看出，散布矩陣是樣本協方差矩陣的(n-1)倍。最后，尋找最優的直線方向e，使得平方誤差準則函數J1值最小。將上述結論從1維推廣到高維，使用![](https://box.kancloud.cn/2016-01-05_568b368fd7953.jpg)維的特征矢量來近似原空間中的d維矢量： ![](https://box.kancloud.cn/2016-01-05_568b368fe58ea.jpg) 其中，e1,e2,…ed’是散布矩陣S的對應于d’個最大特征值的特征矢量，這些矢量之間相互正交，構成了d’維空間的基矢量。而系數ai是矢量x對應于基ei的系數，稱為主分量(Principal Component)。 **3.Fisher線性判別方法** PCA方法尋找的是用來有效表示同一類樣本共同特點的主軸方向，這對于表示同一類數據樣本的共同特征是非常有效的。但PCA不適合用于區分不同的樣本類。Fisher線性判別分析（FDA）是用于尋找最有效地對不同樣本類進行區分的方向。其主要思想是考慮將d維空間中的點投影到一條直線上。通過適當地選擇直線的方向，有可能找到能夠最大限度地區分各類樣本數據點的投影方向。設一組由n個d維樣本組成的數據集合{x1,x2,…,xn}，它們分屬于兩個不同的類ω1和ω2，其中ω1類的樣本子集為D_1，包含n1個樣本點；ω2類的樣本子集為D2，包含n2個樣本點。將每一個樣本點投影到方向為w的直線上，其中||w||=1，有： ![](https://box.kancloud.cn/2016-01-05_568b368ff3378.jpg) 得到n個投影標量y1,y2,…,yn，對應于ω1和ω2，分別屬于集合Y1和Y2，下圖給出同一組樣本點在不同方向投影后可分性不同： ![](https://box.kancloud.cn/2016-01-05_568b36900d9e4.jpg) ![](https://box.kancloud.cn/2016-01-05_568b36902853c.jpg) Fisher判別分析的目的就是要找到使得該準則函數J(.)最大化時的直線w。經過推導，J(.)最大化時的w為： ![](https://box.kancloud.cn/2016-01-05_568b369044518.jpg) ##二、實驗結果討論第一部分PCA主分量分析部分的實驗步驟主要包括以下幾個部分： 1)參考實驗Proj01-01，對二維高斯分布情況，給定均值矢量和協方差矩陣如下： ![](https://box.kancloud.cn/2016-01-05_568b3690514f6.jpg) 2)生成1000個二維樣本矢量的數據集合X。并繪出該樣本集合的二維散點圖；編寫函數[Y,mu_X,scatterMatrix,V,D] = PCA(X)，輸入：散點矩陣X，輸出：Y=D_x (x-m_x) 、mu_X：均值向量、scatterMatrix：散布矩陣、V：特征向量、D：由特征值組成的對角矩陣，計算上述樣本集合的均值向量和散布矩陣，計算散布矩陣的特征值和特征向量（可使用MATLAB中的eig函數）。在二維坐標中，以均值向量為中心，畫出每一個特征向量方向的直線。將表示特征向量方向的直線疊加在在(1)中的二維散點圖上，如下圖所示。 ![](https://box.kancloud.cn/2016-01-05_568b369061542.jpg) 3)假設數據集合X的均值向量為mx，散布矩陣的特征向量矩陣為D_x=[e1 e2]，對集合X中的每個向量x進行下面的變換，生成集合Y：Y=Dx (x-mx)。用不同的顏色分別繪出集合X和集合Y的二維散點圖，如下圖所示。經觀察Y是集合X經平移旋轉后形成的散點圖，且其分布與X正交。 ![](https://box.kancloud.cn/2016-01-05_568b369076b03.jpg) 4)在(1)中，對三維高斯分布情況，生成三維樣本集合，重復實驗(2)和(3)，輸出圖像如下圖所示。同二維散點分布的情況相似，三維散點集合Y是集合X經平移旋轉后形成的散點圖，且其分布與X正交。 ![](https://box.kancloud.cn/2016-01-05_568b36908a90e.jpg) ![](https://box.kancloud.cn/2016-01-05_568b36909ea18.jpg) **第二部分**是Fisher 線性判別分析，主要了解Fisher線性判別方法找的最優方向與非最優方向的差異，將高維分布數據降維為一維，通過Fisher線性判別方法實現高維數據在一維中分類。本實驗的數據見以下表格： ![](https://box.kancloud.cn/2016-01-05_568b3690ba75b.jpg) 1)編寫通用函數fisher(w)實現Fisher線性判別方法。 2)對表格中的類別w2和w3，計算最優方向w并畫出表示最優方向w的直線，并且標記出投影后的點在直線上的位置，如下圖所示。 ![](https://box.kancloud.cn/2016-01-05_568b3690d8e45.jpg) 3)在這個子空間中，對每種分布用一維高斯函數擬合，即把投影后的每一類數據看作滿足一維高斯分布，求出其均值和方差。并且求出分類決策面（兩個一維高斯分布的交點處），如下圖所示。 ![](https://box.kancloud.cn/2016-01-05_568b3690f1ab2.jpg) 4)計算得到的分類器的訓練誤差，如圖8所示。 ![](https://box.kancloud.cn/2016-01-05_568b369116eb7.jpg) 5)為了比較，使用非最優方向w=(1.0,2.0,-1.5)重復(4)(5)兩個步驟，在這個非最優子空間中，得到兩類數據到非最優方向的投影和訓練誤差，如下圖所示。 ![](https://box.kancloud.cn/2016-01-05_568b369127b39.jpg) ![](https://box.kancloud.cn/2016-01-05_568b369145eeb.jpg) **初步結論：**對比最佳方向與非最佳方向的分類結果，前者的效果不如后者，一個原因是因為訓練樣本只有兩類20個樣本，無法體現真實的誤差情況。另一個原因是由于樣本集中出現了距離樣本均值很大的樣本點，當使用非最優方向時，擬合后的一維高斯函數方差非常大（如下圖所示），對于快速收斂的紅色曲線，平坦的藍色曲線在更多的區域的值要高于紅色線。綜上所述，需要使用更多的樣本點才能準確驗證最優方向的分類結果確實是最優。 ![](https://box.kancloud.cn/2016-01-05_568b369156106.jpg) PCA主分量分析與Fisher線性判別分析涉及大量的數學推理，需要進一步的學習！以下是MATLAB源碼： ~~~ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 產生模式類函數 % N：生成散點個數 C：類別個數 d：散點的維數 % mu：各類散點的均值矩陣 % sigma：各類散點的協方差矩陣 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function result = MixGaussian(N, C, d, mu, sigma) color = {'r.', 'g.', 'm.', 'b.', 'k.', 'y.'}; % 用于存放不同類數據的顏色 % if nargin <= 3 & N < 0 & C < 1 & d < 1 % error('參數太少或參數錯誤'); % figure; if d == 1 for i = 1 : C for j = 1 : N/C r(j,i) = sqrt(sigma(1,i)) * randn() + mu(1,i); end end elseif d == 2 for i = 1:C r(:,:,i) = mvnrnd(mu(:,:,i),sigma(:,:,i),round(N/C)); end %title('N類二維隨機散點的分布圖'); elseif d == 3 for i = 1:C r(:,:,i) = mvnrnd(mu(:,:,i),sigma(:,:,i),round(N/C)); end else disp('維數只能設置為1，2或3'); end result = r; ~~~ ~~~ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PCA算法函數 % 輸入參數： % X：樣本點矩陣 % 輸出參數： % Y = Dx(x-mx) % mu_X：均值向量 % scatterMatrix：散布矩陣 % V：特征向量 % D：由特征值組成的對角矩陣 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Y,mu_X,scatterMatrix,V,D] = PCA(X) [N,d] = size(X); % N個d維樣本 % 計算輸入參數的樣本均值和散布矩陣 mu_X = sum(X) / N scatterMatrix = zeros(d,d); for i = 1:N scatterMatrix = scatterMatrix + (X(i,:) - mu_X)' * (X(i,:) - mu_X); end % 計算散布矩陣的特征值和特征向量 % V為特征向量，它的每一列都是特征向量 % D是特征值構成的對角矩陣 % 這些特征值和特征向量都沒有經過排序 [V,D] = eig(scatterMatrix); % 生成集合Y for k = 1:N Y(k,:) = (V * (X(k,:) - mu_X)')'; end % 畫圖線 if d == 2 plot(X(:,1),X(:,2),'r+'); grid on; hold on; title('1000個二維隨機散點的分布圖及特征向量的方向'); for j = 1:d % 畫出帶有箭頭的向量 quiver(mu_X(1), mu_X(2), 8*V(1,j), 8*V(2,j),'b'); hold on; end figure,plot(X(:,1),X(:,2),'r+'); hold on; for j = 1:d % 畫出帶有箭頭的向量 quiver(mu_X(1), mu_X(2), 8*V(1,j), 8*V(2,j),'b'); hold on; end plot(Y(:,1),Y(:,2),'b+'); grid on; title('集合 X (紅色)和 Y = Dx(x-mx) (藍色)的二維隨機散點分布圖'); else if d ==3 plot3(X(:,1),X(:,2),X(:,3),'r+'); grid on; hold on; % 畫出帶有箭頭的向量 for j = 1:d quiver3(mu_X(1),mu_X(2),mu_X(3),15*V(1,j), 15*V(2,j), 15*V(3,j)); hold on; end title('1000個三維隨機散點的分布圖及特征向量的方向'); figure,plot3(X(:,1),X(:,2),X(:,3),'r+'); hold on; % 畫出帶有箭頭的向量 for j = 1:d quiver3(mu_X(1),mu_X(2),mu_X(3),15*V(1,j), 15*V(2,j), 15*V(3,j)); hold on; end plot3(Y(:,1),Y(:,2),Y(:,3),'b+'); grid on; title('集合 X 和 Y = Dx(x-mx) 的三維隨機散點分布圖'); end end ~~~ ~~~ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fisher線性判別方法函數 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function fisher(w) [x,y,z] = size(w); for i = 1:z % wi類的樣本均值 mu(:,:,i) = sum(w(:,:,i)) ./ x; % wi類的協方差矩陣 sigma(:,:,i) = cov(w(:,:,i)); end % 總類內散布矩陣? % sw = sigma(:,:,1) + sigma(:,:,2); s = zeros(y,y,2); for q = 1:z for qt = 1:x s(:,:,q) = s(:,:,q) + (w(qt,:,i) - mu(:,:,q))'*(w(qt,:,i) - mu(:,:,q)) end end sw = s(:,:,1) + s(:,:,2); %投影向量的計算公式,計算最優方向的直線 wk = inv(sw) * (mu(:,:,1)-mu(:,:,2))';% 最優方向的特征向量 % wk = [1; 2; -1.5];% 若不是最優方向 wk = wk/sqrt(sum(wk.^2)); % 特征向量單位化 t1 = zeros(1,3,10); t2 = zeros(1,3,10); %計算將樣本投影到最佳方向上以后的新坐標 cm1=w(:,1,1)*wk(1)+w(:,2,1)*wk(2)+w(:,3,1)*wk(3); cm2=w(:,1,2)*wk(1)+w(:,2,2)*wk(2)+w(:,3,2)*wk(3); new1=[wk(1)*cm1 wk(2)*cm1 wk(3)*cm1]; new2=[wk(1)*cm2 wk(2)*cm2 wk(3)*cm2]; % 投影后的點在直線上的位置 weigh = 2; % 繪制較長的向量 wk = wk * weigh; plot3([-wk(1),0.5 * wk(1)],[-wk(2),0.5 * wk(2)],[-wk(3),0.5 * wk(3)],'g'); for r = 1:10 plot3(new1(r,1),new1(r,2),new1(r,3),'ro'); plot3([w(r,1,1) new1(r,1)],[w(r,2,1) new1(r,2)],[w(r,3,1) new1(r,3)],'r'); hold on; plot3(new2(r,1),new2(r,2),new2(r,3),'bo'); plot3([w(r,1,2) new2(r,1)],[w(r,2,2) new2(r,2)],[w(r,3,2) new2(r,3)],'b'); hold on; end title('最優方向的直線及投影后各點在直線上的位置'); legend('第一類數據到最優方向的投影',... '第二類數據到最優方向的投影',... '最優方向'); % 投影后的標量y1和y2，對應兩類數據 wk = wk / weigh; y1 = wk' * w(:,:,1)'; y2 = wk' * w(:,:,2)'; % 一維高斯擬合 [mu1, sigma1] = normfit(y1); [mu2, sigma2] = normfit(y2); x=-2:0.01:1.5; y1 = normpdf(x,mu1,sigma1); y2 = normpdf(x,mu2,sigma2); figure,plot(x,y1,'r'); hold on; plot(x,y2,'b'); hold on; % 決策面 if (mu2 > mu1) xt = mu1:0.01:mu2; else if (mu2 < mu1) xt = mu2:0.01:mu1; end end % 產生決策面 y1 = normpdf(xt,mu1,sigma1); y2 = normpdf(xt,mu2,sigma2); judge = find(abs(y1-y2) < 0.05); judge = mu2 + 0.01 * judge(1); y=-1:0.1:3; x= judge + 0.*y; % 產生一個和y長度相同的x數組 plot(x,y,'g.-'); hold on; % % 正確分類的散點個數 right1 = 0; right2 = 0; % % 正確率 rightRate1 = 0; rightRate2 = 0; for t = 1:10 P1_c1(t) = normpdf(cm1(t), mu1, sigma1); P2_c1(t) = normpdf(cm1(t), mu2, sigma2); P1_c2(t) = normpdf(cm2(t), mu1, sigma1); P2_c2(t) = normpdf(cm2(t), mu2, sigma2); if (P1_c1(t) > P2_c1(t)) plot(cm1(t),P1_c1(t),'ro'); hold on; right1 = right1 + 1; else plot(cm1(t),P1_c1(t),'b+'); end if (P2_c2(t) > P1_c2(t)) plot(cm2(t),P2_c2(t),'bo'); hold on; right2 = right2 + 1; else plot(cm2(t),P2_c2(t),'r+'); end legend('第一類數據投影的一維高斯擬合',... '第二類數據投影的一維高斯擬合',... '決策面',... '被正確分類的第一類點',... '被正確分類的第二類點'); end rightRate1 = right1 / 10; rightRate2 = right2 / 10; disp(['使用最佳方向時第一類的錯分點個數為：',num2str(10 - right1)]); disp(['使用最佳方向時第一類的分類準確率為：',num2str(rightRate1 * 100),'%']); disp(['使用最佳方向時第二類的錯分點個數為：',num2str(10 - right2)]); disp(['使用最佳方向時第二類的分類準確率為：',num2str(rightRate2 * 100),'%']); grid on; title('兩類數據投影的一維高斯擬合與決策面(最佳投影方向)'); ~~~ ~~~ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PCA 主分量分析與Fisher 線性判別分析實驗主函數 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 生成1000個二維樣本矢量的數據集合X，并繪出該樣本集合的二維散點圖。 N = 1000; % 1000個點 C = 1; % 種類 d = 3; % 維數 % mu = [5 7 ]; % sigma = [9 0.4;0.4 1]; mu = [5 7 9]; sigma = [7 1.2 0.4;1.2 3 0;0.4 0 1]; X = MixGaussian(N, C, d, mu, sigma); PCA(X); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fisher 線性判別分析部分 % w1，w2，w3三類散點 w = zeros(10,3,2); % w(:,:,1) = [ 0.42 -0.087 0.58;... % -0.2 -3.3 -3.4 ;... % 1.3 -0.32 1.7 ;... % 0.39 0.71 0.23;... % -1.6 -5.3 -0.15;... % -0.029 0.89 -4.7 ;... % -0.23 1.9 2.2 ;... % 0.27 -0.3 -0.87;... % -1.9 0.76 -2.1 ;... % 0.87 -1.0 -2.6]; w(:,:,1) = [-0.4 0.58 0.089; ... -0.31 0.27 -0.04; ... 0.38 0.055 -0.035;... -0.15 0.53 0.011;... -0.35 0.47 0.034;... 0.17 0.69 0.1; ... -0.011 0.55 -0.18; ... -0.27 0.61 0.12; ... -0.065 0.49 .0012;... -0.12 0.054 -0.063]; w(:,:,2) = [ 0.83 1.6 -0.014;... 1.1 16. 0.48; ... -0.44 -0.41 0.32; ... 0.047 -0.45 1.4; ... 0.28 0.35 3.1; ... -0.39 -0.48 0.11; ... 0.34 -0.079 0.14; ... -0.3 -0.22 2.2; ... 1.1 1.2 -0.46; ... 0.18 -0.11 -0.49]; figure,plot3(w(:,1,1),w(:,2,1),w(:,3,1),'ro'); grid on; hold on; plot3(w(:,1,2),w(:,2,2),w(:,3,2),'bo'); hold on; fisher(w); ~~~