PCA主成分分析之三維演示(Matlab)
PCA主成分分析之三維演示(Matlab)
寫這個的主要原因是實驗課上的要求，原本需要寫一個演示 PCA 原理的 demo ，按照實驗指導書上來說，在二維上演示就好了，但是為了折騰和無聊裝逼，我寫了這個程式，我覺得更能直觀的看出 PCA 的原理。可以完整的看到資料從三維降到二維，再到一維的整個過程。別看是簡簡單單的感覺，要寫出來，還是需要對原理有比較深刻的認識的，程式碼沒有參照網上，僅靠個人推導和聯想，多多指教。
目錄
原理簡介
主成分分析(principal component analysis PCA),通過正交變換將一組可能存在相關性的變數轉換成一組線性不相關的變數，轉換之後的這組變數叫做主成分。常用於提取資料主要特徵變數，常用於高維資料的降維。在用統計分析的方法研究多變數的時候，變數的個數太多，增加了複雜性，由於變數之間是有一定的相關性的，也就是說不同的變數所代表的全域性資訊是有重疊，有冗餘的。其演算法流程如下圖所示：
結果展示
如圖1所示繪製隨機的3維點，如圖2所示，選取特徵值最大的兩個特徵向量進行投影，繪製其投影平面，如圖3所示，繪製出3維到2維的投影垂線以示意這個過程的原理。如圖4所示為最後投影得到的二維點在三維座標中的位置，可以看到，這些點都位於一個平面上，如圖5所示為二維點到一維點的投影線，所選的投影基為特徵值最大的特徵向量，如圖6所示，繪製了二維到一維的投影垂線，如圖7所示為二維降為一維的點。整個過程演示了一個三維特逐漸降為一維的過程，特徵有所損失，但是仍能被區分。
圖1 隨機三維點
|
圖2 三維到二維投影面
|
 |
 圖4 三維到二維的投影點
|
 圖5 二維到一維的投影線
|
 圖6 二維到一維的投影過程
|
 圖7 二維到一維的投影點
|
原始碼
function PCAdemo
close all;clear;clc;
h_figure = figure;
mu = [0 0 0];
%協方差矩陣，對角為方差值0.3，0.35
dim = 3;
var = [0.2 0 0; 0 0.5 0; 0 0 0.8];
samNum = 100;
data = mvnrnd(mu,var,samNum);
h_plotMain = plot3(gca,data(:,1),data(:,2),data(:,3),'o',...
'MarkerSize',5,'MarkerFaceColor','k');
xlim([-2 2]);ylim([-2 2]);zlim([-2 2]);
grid on;axis square;hold on;
pause();
covdata = cov(data);%求協方差矩陣
[eigVector, eigValue] = eig(covdata);%求協方差矩陣的特徵值和特徵向量
eigValue = diag(eigValue)';
[eigValue_sort, IX] = sort(eigValue, 2, 'descend');
eigVector_sort = [];
for i = 1:dim
tempIX = IX(i);
eigVector_sort = [eigVector_sort eigVector(:, tempIX)];
end
base_3to2 = eigVector_sort(:,1:2);%三維降2維
base_3to1 = eigVector_sort(:,1);%三維降1維
base_3to2 = base_3to2';
base_3to1 = base_3to1';
%%%%%繪製3維變成2維度的投影平面%%%%%
A = base_3to2(1, :);%新基i帽(相對原基)
B = base_3to2(2, :);%新基j帽(相對原基)
C = [0 0 0];
syms x y z;
D = [ones(4, 1), [[x, y, z]; A; B; C]];%由空間解析幾何的內容知道D的行列式等於零就是平面方程。
detd = det(D);
z = solve(detd,z);
ezmesh(z,[-2, 2, -2, 2]);
pause();
%%%%%%%%%%%%
data_proj3_2 = base_3to2*data';%投影到新基座標(相對新基)
data_proj3_1 = base_3to1*data';%投影到新基座標(相對新基)
data_proj3_2 = data_proj3_2';
data_proj3_1 = data_proj3_1';
data_proj_respect_to_orienbasis3_2 = [];%新基座標用原基表示
data_proj_respect_to_orienbasis3_1 = [];%新基座標用原基表示
for i = 1:samNum
data_proj_respect_to_orienbasis3_2 = [data_proj_respect_to_orienbasis3_2;...
data_proj3_2(i,1)*base_3to2(1,:)+data_proj3_2(i,2)*base_3to2(2,:)];
data_proj_respect_to_orienbasis3_1 = [data_proj_respect_to_orienbasis3_1;...
data_proj3_1(i)*base_3to1(1,:)];
end
%繪製3維降為2維後的資料
h_plot3_2 = plot3(gca,data_proj_respect_to_orienbasis3_2(:,1),...
data_proj_respect_to_orienbasis3_2(:,2),data_proj_respect_to_orienbasis3_2(:,3),...
'o','MarkerSize',5,'MarkerEdgeColor','g','MarkerFaceColor','g');
%繪製3維降為2的投影虛線
hl3_2 = [];
hl3_1 = [];
for i = 1:samNum
dd = [data_proj_respect_to_orienbasis3_2(i,:);data(i,:)];
hl3_2(i) = plot3(gca,dd(:,1),dd(:,2),dd(:,3),'-.','markersize',10);
end
pause();
delete(h_plotMain);delete(hl3_2);
pause();
%繪製3維降為1的基準線
h_zhixian = plot3(gca,data_proj_respect_to_orienbasis3_1(:,1),...
data_proj_respect_to_orienbasis3_1(:,2),data_proj_respect_to_orienbasis3_1(:,3),...
'color','r','linewidth',2);
pause();
%繪製3維降為1的資料
h_plot3_1 = plot3(gca,data_proj_respect_to_orienbasis3_1(:,1),...
data_proj_respect_to_orienbasis3_1(:,2),data_proj_respect_to_orienbasis3_1(:,3),...
'o','MarkerSize',5,'MarkerEdgeColor','r','MarkerFaceColor','r');
%繪製3維降為1的投影虛線
for i = 1:samNum
dd=[data_proj_respect_to_orienbasis3_1(i,:);data_proj_respect_to_orienbasis3_2(i,:)];
hl3_1(i)=plot3(gca,dd(:,1),dd(:,2),dd(:,3),'-.','markersize',10);
end
pause();
delete(h_plot3_2);
delete(hl3_1);
delete(h_zhixian);
end 可以直接複製以上程式碼，或者可以下載連結：https://download.csdn.net/download/qq_33826564/10574325
GitHub: https://github.com/swq123459/PCA-3D-demo
體會
1.使用PCA對資料進行降維的思路：將輸入的樣本看成一個隨機向量，即
有M個樣本，說明發生了M次 事件， 有N維，每一維為單獨的隨機變數 ，若要將降維，降維後的需滿足其內部各維之間的協方差為0，已保證降維後的資料無冗餘，每一位表示的資訊不正交不重疊；同時各維的方差需要最大，因為方差大代表著資料集中性越差，說明越方便被辨別。隨機向量協方差矩陣中，各維均值位於對角線，協方差位於其他非對角線位置，我們需要對角線值最大，其他位置值為0，這不免讓我們想到了矩陣的對角化。對角化完成的功能真是將矩陣除對角線為其他元素變成0.矩陣的對角化表示式為
其中A為待對角化矩陣，T為A的特徵向量組成的矩陣，這裡我們將A看成是的協方差矩陣，那麼T就是這個協方差矩陣特徵向量組成的矩陣。對角化的過程已經完成了協方差0的任務，接下來使各維的方差最大，即B的對角線元素最大。因為特徵值大的特徵向量是矩陣變換的主要方向，所以其佔的變換成分應該比較大，對B的對角線元素影響也比較大，所以這裡才有根據特徵值來排序特徵向量來組成矩陣T的方法。按照成分的貢獻率，可以執行降維，維度貢獻率不能少於 。
在計算協方差矩陣的時候，A為樣本資料集矩陣，N為樣本維數，需要注意的是樣本資料集A必須去均值，即使樣本各維均值為0，這樣使用該公式計算出來的才是協方差矩陣。
在執行資料降維的時候，從二維資料降至一維的變換為A，得到的這個一維的資料是以A的基為基的座標再乘以這個基的範數，若這個變換本身就同其基一致的話，得到的一維資料為其投影長度。Matlab中計算特徵向量的函式已經完成了特徵向量單位化的操作，所以，將資料同其內積得到的便是投影長度。實驗中涉及將這個投影長度在原來的維度中顯示的問題，只要將這個投影長度乘以原來投影的基向量就行了。
Matlab中如果要在上寫文中索引當前axes，進行繪圖， 只需要進行一次hold on,就可以防止plot控制代碼的自動回收。