2. 程式人生 > >PCA的原理及MATLAB實現

PCA的原理及MATLAB實現

阿新 • • 發佈：2019-02-01

相關文章

－－－－－－－－－－－－－－－－－－－－－－－－－－－－

關於PCA原理可以直接參考下面的文章
深入理解PCA

PCA是經常用來減少資料集的維數，同時保留資料集中對方差貢獻最大的特徵來達到簡化資料集的目的。
PCA的原理就是將原來的樣本資料投影到一個新的空間中，相當於我們在矩陣分析裡面學習的將一組矩陣對映到另外的座標系下。通過一個轉換座標，也可以理解成把一組座標轉換到另外一組座標系下，但是在新的座標系下，表示原來的原本不需要那麼多的變數，只需要原來樣本的最大的一個線性無關組的特徵值對應的空間的座標即可。
PCA演算法核心思想就是將 n 維特徵對映到 k 維上（k < n），這 k 維是全新的正交特徵。我們將這 k 維成為主元，是重新構造出來的 k 維特徵，而不是簡單地從 n 維特徵中取出其餘 n-k 維特徵。

PCA實現的步驟
（1）把原始資料中每個樣本用一個向量表示，然後把所有樣本組合起來構成一個矩陣。為了避免樣本的單位的影響，樣本集需要標準化（一般都是去均值化）。
（2）求該矩陣的協方差矩陣。
（3）求步驟2中得到的協方差矩陣的特徵值和特徵向量。
（4）將求出的特徵向量按照特徵值的大小進行組合形成一個對映矩陣，並根據指定的PCA保留的特徵個數取出對映矩陣的前n行或者前n列作為最終的對映矩陣。
（5）用步驟4的對映矩陣對原始資料進行對映，達到資料降維的目的。

function [Y,V,E,D] = newpca(X)
%其中X為輸入資料，X的每一列是一個輸入樣本。返回值Y是對X進行PCA分析後的投影矩陣。
 

%V是與X有關的協方差矩陣特徵向量的白化矩陣，E是對應的特徵向量（列）構成的矩陣，
%D是對應的特徵值構成的對角矩陣（特徵值處於對角線上）。
%返回值中的白化矩陣，特徵向量和特徵值都是按照對應特徵值大小進行排序後了的。
% do PCA on image patches
%
% INPUT variables:
% X                  matrix with image patches as columns
%
% OUTPUT variables:
% Y                  the project matrix of the input data X without whiting
 

% V                  whitening matrix
% E                  principal component transformation (orthogonal)
% D                  variances of the principal components

%去除直流成分
X = X-ones(size(X,1),1)*mean(X);

% Calculate the eigenvalues and eigenvectors of the new covariance matrix.
covarianceMatrix = X*X'/size(X,2); %求出其協方差矩陣
%E是特徵向量構成，它的每一列是特徵向量，D是特徵值構成的對角矩陣
%這些特徵值和特徵向量都沒有經過排序
[E, D] = eig(covarianceMatrix); 

% Sort the eigenvalues  and recompute matrices
% 因為sort函式是升序排列，而需要的是降序排列，所以先取負號,diag(a)是取出a的對角元素構成
% 一個列向量，這裡的dummy是降序排列後的向量，order是其排列順序
[dummy,order] = sort(diag(-D));
E = E(:,order);%將特徵向量按照特徵值大小進行降序排列，每一列是一個特徵向量
Y = E'*X;
d = diag(D); %d是一個列向量
%dsqrtinv是列向量，特徵值開根號後取倒,仍然是與特徵值有關的列向量
%其實就是求開根號後的逆矩陣
dsqrtinv = real(d.^(-0.5)); 
Dsqrtinv = diag(dsqrtinv(order));%是一個對角矩陣，矩陣中的元素時按降序排列好了的特徵值（經過取根號倒後）
D = diag(d(order));%D是一個對角矩陣，其對角元素由特徵值從大到小構成
V = Dsqrtinv*E';%特徵值矩陣乘以特徵向量矩陣
function [lowData,reconMat] = PCA(data,K)
[row , col] = size(data);
meanValue = mean(data);
%varData = var(data,1,1);
normData = data - repmat(meanValue,[row,1]);
covMat = cov(normData(:,1),normData(:,2));%求取協方差矩陣
[eigVect,eigVal] = eig(covMat);%求取特徵值和特徵向量
[sortMat, sortIX] = sort(eigVal,'descend');
[B,IX] = sort(sortMat(1,:),'descend');
len = min(K,length(IX));
eigVect(:,IX(1:1:len));
lowData = normData * eigVect(:,IX(1:1:len));
reconMat = (lowData * eigVect(:,IX(1:1:len))') + repmat(meanValue,[row,1]); % 將降維後的資料轉換到新空間
end

下面給出PCA的Matlab實現部分：Matlab自帶的函式
Matlab中已經包含了實現了的PCA演算法，可以通過princomp函式呼叫。其形式為：
[COEFF,SCORE, latent]=princomp(X);
X：為要輸入的n維原始資料。帶入這個matlab自帶函式，將會生成新的n維加工後的資料（即SCORE）。此資料與之前的n維原始資料一一對應。
COEFF為主成分分量，即變換空間中的那些基向量，是係數矩陣。通過COEFF可以知道x是怎樣轉換成SCORE的。
在n行p列的資料集X上做主成分分析。返回主成分系數。X的每行表示一個樣本的觀測值，每一列表示特徵變數。COEFF是一個p行p列的矩陣，是X矩陣所對應的協方差陣V的所有特徵向量組成的矩陣，即變換矩陣或稱投影矩陣，COEFF每列對應一個特徵值的特徵向量，列的排列順序是按特徵值的大小遞減排序。
SCORE為主成分，即X的低維表示。生成的n維加工後的資料存在SCORE裡。它是對原始資料進行的分析，進而在新的座標系下獲得的資料。他將這n維資料按貢獻率由大到小排列。（即在改變座標系的情況下，又對n維資料排序）
latent為一個包含樣本協方差矩陣的本徵值的向量。是一維列向量，每一個數據是對應SCORE裡相應維的貢獻率，因為資料有n維所以列向量有n個數據。由大到小排列（因為SCORE也是按貢獻率由大到小排列）。
princomp這個matlab自帶的函式，在降維之前就將每一個樣本減去了一個所有樣本的平均值。princomp這裡使用一行表示一個樣本，每行包括這個樣本的所有的特徵值。

其實我們要的是由函式[COEFF,SCORE, latent] = princomp(X)它所產生的pc和latent。由latent可以算出降維後的空間所能表示原空間的程度，只要這個累積的值大於95%就行了
cumsum(latent)./sum(latent)

舉例說明：

load hald; %載入matlab內部資料
[COEFF,SCORE,latent,tsquare] = princomp(ingredients); %呼叫pca分析函式
% 下面我們具體說明COEFF，SCORE，latent

%下面為計算ingredients所對應的協方差矩陣(也就是cov_ingredients矩陣)的特徵值和特徵向量，下面的矩陣V為特徵向量，D為特徵值(對比上面的latent)組成的對角線矩陣
[V,D] = eig(cov_ingredients)

V =

    0.5062    0.5673    0.6460   -0.0678
    0.4933   -0.5440    0.0200   -0.6785
    0.5156    0.4036   -0.7553    0.0290
    0.4844   -0.4684    0.1085    0.7309


D =

    0.2372         0         0         0
         0   12.4054         0         0
         0         0   67.4964         0
         0         0         0  517.7969

%說明1：對比一下矩陣V和矩陣COEFF，現在很容易明白為什麼COEFF是按列遞減順序排列的了！(V中第三列與COEFF中倒數第三列差個負號，學過線性代數的人都知道這沒問題)

%下面再驗證一下說明2
diag(cov(SCORE))

ans =

  517.7969
   67.4964
   12.4054
    0.2372
%說明2：以上結果顯示latent確實表示SCORE矩陣每列的方差，517.7969表示第一列方差

[COEFF,SCORE,latent,tsquare] = princomp(ingredients);
% 上面這個函式，COEFF矩陣是返回的轉換矩陣，也就是把樣本轉換到新的空間中的準換矩陣，這個準換矩陣式比較大的。
% SCORE是原來的樣本矩陣在新的座標系中的表示，也就是原來的樣本乘上轉換矩陣，但是,還不是直接乘，要減去一個樣本的均值。將原來的資料轉換到新的樣本空間中的演算法是這樣實現的：
x0 = bsxfun(@minus,ingredients,mean(ingredients,1));
NewSCORE1 = x0 * COEFF;

NewSCORE1 =

   36.8218   -6.8709   -4.5909    0.3967
   29.6073    4.6109   -2.2476   -0.3958
  -12.9818   -4.2049    0.9022   -1.1261
   23.7147   -6.6341    1.8547   -0.3786
   -0.5532   -4.4617   -6.0874    0.1424
  -10.8125   -3.6466    0.9130   -0.1350
  -32.5882    8.9798   -1.6063    0.0818
   22.6064   10.7259    3.2365    0.3243
   -9.2626    8.9854   -0.0169   -0.5437
   -3.2840  -14.1573    7.0465    0.3405
    9.2200   12.3861    3.4283    0.4352
  -25.5849   -2.7817   -0.3867    0.4468
  -26.9032   -2.9310   -2.4455    0.4116

NewSCORE2 = ingredients * COEFF;

NewSCORE2 =

   25.9105   -7.0189  -35.8545   48.5266
   18.6960    4.4628  -33.5112   47.7341
  -23.8931   -4.3530  -30.3613   47.0039
   12.8034   -6.7821  -29.4088   47.7514
  -11.4645   -4.6098  -37.3510   48.2723
  -21.7238   -3.7946  -30.3506   47.9950
  -43.4995    8.8318  -32.8699   48.2117
   11.6951   10.5779  -28.0270   48.4543
  -20.1739    8.8373  -31.2805   47.5862
  -14.1953  -14.3053  -24.2171   48.4705
   -1.6913   12.2380  -27.8352   48.5651
  -36.4962   -2.9297  -31.6503   48.5768
  -37.8145   -3.0790  -33.7091   48.5416

%對比NewSCORE1和NewSCORE2 與 SCORE，可以看成 SCORE與NewSCORE1是一樣。
% 下面我們就來確定最後的轉換矩陣
cumsum(latent)./sum(latent)
ans =

    0.8660
    0.9789
    0.9996
    1.0000

 SelectNum = cumsum(latent)./sum(latent)
 index = find(SelectNum >= 0.95);
 ForwardNum = index(1);
%選擇由latent可以算出降維後的空間所能表示原空間的程度，只要這個累積的值大於95%的數目ForwardNum，就選擇前ForwardNum個主成分作為變換矩陣。
%做主成分變換矩陣，也就是降維矩陣
tranMatrix = COEFF(:,1:ForwardNum)

tranMatrix =

   -0.0678   -0.6460
   -0.6785   -0.0200
    0.0290    0.7553
    0.7309   -0.1085

% Newingredients = x0 * tranMatrix
Newingredients =

   36.8218   -6.8709
   29.6073    4.6109
  -12.9818   -4.2049
   23.7147   -6.6341
   -0.5532   -4.4617
  -10.8125   -3.6466
  -32.5882    8.9798
   22.6064   10.7259
   -9.2626    8.9854
   -3.2840  -14.1573
    9.2200   12.3861
  -25.5849   -2.7817
  -26.9032   -2.9310
  % 這樣Newingredients就是ingredients經過PCA降維的最終的結果。

測試樣本降維
如果你需要對測試樣本降維，一般情況下，使用matlab自帶的方式，肯定需要對測試樣本減去一個訓練樣本均值，因為你在給訓練樣本降維的時候減去了均值，所以測試樣本也要減去均值，然後乘以coeff這個矩陣，就獲得了測試樣本降維後的資料。使用訓練樣本得到的轉換矩陣，保證訓練樣本和測試樣本轉換到相同的樣本空間中。比如說你的測試樣本是1*1000000，那麼乘上一個1000000*29的降維矩陣，就獲得了1*29的降維後的測試樣本的降維資料。

princomp(x)使用的行表示一個樣本，每行的所有的列資料都是這個樣本的特徵值。降維以後比如是30*29，那麼每一行就是降維以後的資料。每個樣本有29個特徵值。

%X為訓練集，Xtest為測試集
[COEFF,SCORE, latent] = princomp(X);
 SelectNum = cumsum(latent)./sum(latent)
 index = find(SelectNum >= 0.95);
 ForwardNum = index(1);
%選擇由latent可以算出降維後的空間所能表示原空間的程度，只要這個累積的值大於95%的數目ForwardNum，就選擇前ForwardNum個主成分作為變換矩陣。
%做主成分變換矩陣
NewtrainScores = SCORE(:,1:ForwardNum);
**NewtrainScores就是X經過PCA降維後的資料**

%
tranMatrix = COEFF(:,1:ForwardNum)
[row , col] = size(Xtest);
meanValue = mean(X);
normXtest = Xtest - repmat(meanValue,[row,1]);
% 第二種方法
%normXtest = bsxfun(@minus,Xtest,meanValue); % center the test data
NewtestScores = normXtest*tranMatrix;
%**NewtestScores 就是Xtest經過PCA降維後的資料。**

備註說明：
對測試樣本進行降維的時候，一定要減去訓練樣本的均值，使用訓練樣本得到的轉換矩陣，保證訓練樣本和測試樣本轉換到相同的樣本空間中，這樣才有意義。
很多人就選擇還是使用princomp這個函式處理測試樣本，那麼這樣測試樣本被對映到一個新的空間中，和原來的訓練樣本完全不是在一個空間，一點意義都沒有，還是要使用測試樣本減去訓練樣本的均值，然後乘上訓練樣本降維的時候獲得降維矩陣，轉換到相同的空間中。

project the test data on the reduced set of the principal 
components. Note that princomp centers data before computing the principal components. You need to center the test data as well. It is best to do that using the means obtained from the training data. Something like this would 
do:

mu = mean(Xtrain);
[coefs,scores,variances] = princomp(Xtrain,'econ');
% do your thing to select the components
Xtest = bsxfun(@minus,Xtest,mu); % center the test data
testScores = Xtest*coefs(:,1:d);

相關推薦

PCA的原理及MATLAB實現

相關文章 －－－－－－－－－－－－－－－－－－－－－－－－－－－－ 關於PCA原理可以直接參考下面的文章 深入理解PCA PCA是經常用來減少資料集的維數，同時保留資料集中對方差貢獻最大的特徵來達到簡化資料集的目的。 PCA的原理就是將

區域生長演算法原理及MATLAB實現

% Segment based on area, Region Growing; clear all; close all; clc [fileName,pathName] = uigetfile('*.*','Please select an image');%檔案筐，選擇檔案 if(fileName)

直方圖均衡原理及Matlab實現

直方圖均衡（Histogram Equalization） 預備知識（PMF、CDF） 概述 實現 Matlab程式碼 預備知識（PMF、CDF） 1.PMF 概率質量函式（Probability Mass Function）是離散隨機變數在各特

PCA 原理及 Python 實現

原文發於個人部落格 前言 說好的要做個有逼格的技術部落格，雖然這篇依然沒什麼水平，但總算走出了第一步，希望以後每天都能進步一點點吧！ 接觸機器學習也一年多了，也學了很多演算法，而PCA是資料預處理中一個很重要的演算法，當時學習的時候也在網上看了

K-means原理及Matlab實現

前言 作為勵志在機器學習界闖出一片天地的小女子，在此整理學習到的機器學習方法，並使用Matlab及Python實現。希望可以和大家互相交流和探討。 聯絡方式：[email protected] K-means原理 k-means應該是入

邏輯迴歸原理及matlab實現

個人部落格文章連結：  http://www.huqj.top/article?id=163 對於某些分類問題，自變數可能是連續的，但是因變數卻可能是離散的，例如：根據腫瘤大小判斷該腫瘤是否是良性。這種問題不適合用線性迴歸來解決，雖然可以將連續的因變數值對映到離散的分類上，但

吉布斯取樣——原理及matlab實現

MCMC(馬爾可夫鏈蒙特卡洛方法)：the Gibbs Sampler(吉布斯取樣)         在之前的部落格中，我們對比了在一個多元概率分佈p(x)中，分別採用分組（block-wise）和分成分（component-wise）實現梅特羅波利斯哈斯廷斯演算法。對

LDA演算法原理及matlab實現

注：本文非筆者原創，原文轉載自：http://blog.csdn.net/porly/article/details/8020696 1. LDA是什麼 線性判別式分析（Linear Discriminant Analysis），簡稱為LDA。也稱為Fisher線性

標記符控制的分水嶺演算法原理及matlab實現

-------------------------------------------------------------------------------------------------------------------- 附錄A  教程【3】給出的mat

深入淺出PID控制演算法（二）————PID演算法離散化和增量式PID演算法原理及Matlab實現

引言 上篇介紹了連續系統的PID演算法，但是計算機控制是一種取樣控制，他只能根據取樣時刻的偏差來計算控制量，因此計算機控制系統中，必須對公式進行離散化，具體就是用求和代替積分，用向後差分來代替微分，使模擬PID離散化為數字形式的差分方程。 準備工

BP神經網絡（原理及MATLAB實現）

gen 性能 oct temp png error alt 交叉驗證 對比 人工神經網絡概述： 人工神經元模型： 神經網絡的分類： 按照連接方式，可以分為：前向神經網絡 vs. 反饋（遞歸）神經網絡； 按照學習方式，可以分為：有導師學習神經網絡 vs. 無導

增強影象對比度演算法原理及matlab程式碼實現

clc;   close all;   clear all;    % -------------Gamma Transformations-----------------     %f = imread('Fig0316(4)(bottom_left).tif');     

用Burg法估計AR模型的引數原理詳解及matlab實現

用Burg法估計AR模型的引數。藉助如圖所示的格型預測誤差濾波器，伯格法通過求出前向預測誤差和後向預測誤差的平均功率來選取最佳的反射係數k，使誤差的平均功率取得最小值，進而通過反饋求出模型係數和噪聲方差，該演算法可直接通過分析觀測資料得到需要的模型引數，而不用取求解計算量較大

斯坦福機器學習筆記1：GDA高斯判別分析演算法的原理及matlab程式實現

ps：我本身沒有系統的學過matlab程式設計，所以有的方法，比如求均值用mean（）函式之類的方法都是用很笨的方法實現的，所以有很多需要改進的地方，另外是自學實現的程式，可能有的地方我理解錯誤，如果有錯誤請提出來，大家一起學習，本人qq553566286 首先，本文用到的

邏輯迴歸原理介紹及Matlab實現

一、邏輯迴歸基本概念 1. 什麼是邏輯迴歸 邏輯迴歸就是這樣的一個過程：面對一個迴歸或者分類問題，建立代價函式，然後通過優化方法迭代求解出最優的模型引數，然後測試驗證我們這個求解的模型的好壞。 Logistic迴歸雖然名字裡帶“迴歸”，但是它實際上是一種分類方法，主

模糊PID基本原理及matlab模擬實現（新手！新手！新手！）

有關模糊pid的相關知識就把自己從剛接觸到模擬出結果看到的大部分資料總結一下，以及一些自己的ps 以下未說明的都為轉載內容 在講解模糊PID前，我們先要了解PID控制器的原理（本文主要介紹模糊PID的運用,對PID控制器的原理不做詳細介紹）。PID控制器（比例

短波通道模型--多徑瑞利通道原理詳解及matlab實現

        瑞利衰落是一種小尺度衰落效應，它總是疊加於大尺度衰落效應上如衰減、陰影等。發射端和接收端相對運動速度的大小決定了通道衰落的快慢。相對運對導致接收訊號存在多普勒頻移，即通道衰落的快慢與多普勒頻移的大小有關。        基於頻率單調衰落的基帶等效Jakes模型

中值濾波原理及MATLAB演算法實現

中值濾波是一種非線性濾波方式，它依靠模板來實現。 對於一維中值濾波，設模板的尺寸為 M ，M=2*r+1，r為模板半徑，給定一維訊號f(i)，i = 1，2，3……N，則中值濾波輸出為： g(i) = median[ f(j-r)，f(j-r+1),…………，f(j)，f(

基於DCT變換（變換域）實現資訊隱藏(數字水印)原理詳解及matlab實現

主要就是實現了數字水印的嵌入提取和在不同攻擊如旋轉剪下噪聲等下的提取效果差異的比對1 DCT變換的原理2  DCT變換的特點    在基於DCT的變換編碼中，影象是先經分塊（8×8或16×16）後再經DCT，這種變換是區域性的，只反映了影象某一部分的資訊。當然也可以對整幅影象

雙目相機標定以及立體測距原理及OpenCV實現

ada eps 左右 flag 並且 length turn iostream 中一 作者：dcrmg 單目相機標定的目標是獲取相機的內參和外參，內參（1/dx,1/dy,Cx,Cy,f）表征了相機的內部結構參數，外參是相機的旋轉矩陣R和平移向量t。內參中dx和d