主成分分析（pca）演算法的實現步驟及程式碼
阿新 • • 發佈：2019-01-08
%%%%%%%%%%%%開啟一個30行8列資料的txt檔案%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%第一步：輸入樣本矩陣%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename='src.txt';
fid=fopen(filename,'r');
vector=fscanf(fid,'%g',[30,8]);%輸入一個矩陣，此矩陣有8個維度30個樣本
%現對其進行pca降維
%%
%第二步：計算樣本中每一維的均值，然後計算觀察值與均值之間的偏差，再計算協方差矩陣
s=sum(vector,1);%計算每一列的均值
[nSmp,nFea] = size(vector);
vector=vector-repmat(mean(vector),nSmp,1);%把每一維的均值複製成nSmp*nFea矩陣，然後計算偏差
W=zeros(nFea,nFea);%建立一個nFea*nFea的零矩陣
W=W+vector'*vector;
W=W/(nSmp-1);%根據協方差公式計算協方差，得到協方差矩陣W
%%
%第三步：計算協方差矩陣的特徵值和特徵向量矩陣
fprintf(1,'Calculating generalized eigenvectors and eigenvalues...\n');
[eigvectors, eigvalues] = eig(W);%eigvectors為特徵向量組成的矩陣，eigvalues特徵值組成的對角矩陣
fprintf(1,'Sorting eigenvectors according to eigenvalues...\n');
d1=diag(eigvalues);%返回對角矩陣上的值
[d2 index]=sort(d1); %以升序排序，d2為排列後的值，index為索引值
cols=size(eigvectors,2);% 特徵向量矩陣的列數
for i=1:cols
vsort(:,i) = eigvectors(:, index(cols-i+1) ); % vsort 是一個M*col(注:col一般等於M)階矩陣，儲存的是按降序排列的特徵向量,每一列構成一個特徵向量
dsort(i) = d1( index(cols-i+1) ); % dsort 儲存的是按降序排列的特徵值，是一維行向量
end %完成降序排列
%%
%第四步：計算總能量，並選取貢獻率最大的特徵值
dsum = sum(dsort);%對所有的特徵值求和
dsum_extract = 0;%求前幾個特徵值之和,當前幾個特徵值之和大於90%時，可以認為這幾個特徵值可以表徵當前矩陣
p = 0;
while( dsum_extract/dsum < 0.9)
p = p + 1;
dsum_extract = sum(dsort(1:p));
end
%當p=5時，貢獻率為0.9095，所以就可以選取5個特徵值代表當前矩陣
%%
%第五步：計算前p個特徵值所對應的的特徵向量組成的矩陣，計算降維後的樣本矩陣
vsort = vsort(:,1:p);%提取前p個特徵生成8*p矩陣
fprintf(1,'Feature extraction and calculating newData...\n');
newvector=vector*vsort;%生成nSmp行p列的矩陣，達到降維的效果
%%
%第一步：輸入樣本矩陣%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename='src.txt';
fid=fopen(filename,'r');
vector=fscanf(fid,'%g',[30,8]);%輸入一個矩陣，此矩陣有8個維度30個樣本
%現對其進行pca降維
%%
%第二步：計算樣本中每一維的均值，然後計算觀察值與均值之間的偏差，再計算協方差矩陣
s=sum(vector,1);%計算每一列的均值
[nSmp,nFea] = size(vector);
vector=vector-repmat(mean(vector),nSmp,1);%把每一維的均值複製成nSmp*nFea矩陣，然後計算偏差
W=zeros(nFea,nFea);%建立一個nFea*nFea的零矩陣
W=W+vector'*vector;
W=W/(nSmp-1);%根據協方差公式計算協方差，得到協方差矩陣W
%%
%第三步：計算協方差矩陣的特徵值和特徵向量矩陣
fprintf(1,'Calculating generalized eigenvectors and eigenvalues...\n');
[eigvectors, eigvalues] = eig(W);%eigvectors為特徵向量組成的矩陣，eigvalues特徵值組成的對角矩陣
fprintf(1,'Sorting eigenvectors according to eigenvalues...\n');
d1=diag(eigvalues);%返回對角矩陣上的值
[d2 index]=sort(d1); %以升序排序，d2為排列後的值，index為索引值
cols=size(eigvectors,2);% 特徵向量矩陣的列數
for i=1:cols
vsort(:,i) = eigvectors(:, index(cols-i+1) ); % vsort 是一個M*col(注:col一般等於M)階矩陣，儲存的是按降序排列的特徵向量,每一列構成一個特徵向量
dsort(i) = d1( index(cols-i+1) ); % dsort 儲存的是按降序排列的特徵值，是一維行向量
end %完成降序排列
%%
%第四步：計算總能量，並選取貢獻率最大的特徵值
dsum = sum(dsort);%對所有的特徵值求和
dsum_extract = 0;%求前幾個特徵值之和,當前幾個特徵值之和大於90%時，可以認為這幾個特徵值可以表徵當前矩陣
p = 0;
while( dsum_extract/dsum < 0.9)
p = p + 1;
dsum_extract = sum(dsort(1:p));
end
%當p=5時，貢獻率為0.9095，所以就可以選取5個特徵值代表當前矩陣
%%
%第五步：計算前p個特徵值所對應的的特徵向量組成的矩陣，計算降維後的樣本矩陣
vsort = vsort(:,1:p);%提取前p個特徵生成8*p矩陣
fprintf(1,'Feature extraction and calculating newData...\n');
newvector=vector*vsort;%生成nSmp行p列的矩陣，達到降維的效果