主元分析
主成分分析 ( Principal Component Analysis , PCA )或者主元分析。是一種掌握事物主要矛盾的統計分析方法,它可以從多元事物中解析出主要影響因素,揭示事物的本質,簡化複雜的問題。計算主成分的目的是將
高維數據投影到較低
維空間。給定 n 個變數的 m 個觀察值,形成一個 n ′ m 的數據
矩陣, n 通常比較大。對於一個由多個變數描述的複雜事物,人們難以認識,那么是否可以抓住事物主要方面進行重點分析呢?如果事物的主要方面剛好體現在幾個主要變數上,我們只需要將這幾個變數分離出來,進行詳細分析。但是在一般情況下,並不能直接找出這樣的關鍵變數。這時我們可以用原有變數的
線性組合來表示事物的主要方面, PCA 就是這樣一種分析方法。
PCA 主要 用於數據降維,對於一系列例子的特徵組成的多維向量,多維向量里的某些元素本身沒有區分性,比如某個元素在所有的例子中都為1,或者與1差距不大,那么這個元素本身就沒有區分性,用它做特徵來區分,貢獻會非常小。所以我們的目的是找那些變化大的元素,即
方差大的那些維,而去除掉那些變化不大的維,從而使特徵留下的都是“精品”,而且計算量也變小了。 對於一個k維的特徵來說,相當於它的每一維特徵與其他維都是正交的(相當於在多維
坐標系中,坐標軸都是垂直的),那么我們可以變化這些維的坐標系,從而使這個特徵在某些維上方差大,而在某些維上方差很小。例如,一個45度傾斜的橢圓,在第一坐標系,如果按照x,y坐標來投影,這些點的x和y的屬性很難用於區分他們,因為他們在x,y軸上坐標變化的方差都差不多,我們無法根據這個點的某個x屬性來判斷這個點是哪個,而如果將坐標軸旋轉,以橢圓
長軸為x軸,則橢圓在長軸上的分布比較長,方差大,而在短軸上的分布短,方差小,所以可以考慮只保留這些點的長軸屬性,來區分橢圓上的點,這樣,區分性比x,y軸的方法要好!
所以我們的做法就是求得一個k維特徵的投影矩陣,這個投影矩陣可以將特徵從高維降到低維。投影矩陣也可以叫做
變換矩陣。新的低維特徵必須每個維都正交,
特徵向量都是正交的。通過求樣本矩陣的協方差矩陣,然後求出協方差矩陣的特徵向量,這些特徵向量就可以構成這個投影矩陣了。特徵向量的選擇取決於協方差矩陣的
特徵值的大小。
舉例:
對於一個訓練集,100個對象模板,特徵是10維,那么它可以建立一個100*10的矩陣,作為樣本。求這個樣本的協方差矩陣,得到一個10*10的協方差矩陣,然後求出這個協方差矩陣的特徵值和特徵向量,應該有10個特徵值和特徵向量,我們根據特徵值的大小,取前四個特徵值所對應的特徵向量,構成一個10*4的矩陣,這個矩陣就是我們要求的特徵矩陣,100*10的樣本矩陣乘以這個10*4的特徵矩陣,就得到了一個100*4的新的降維之後的樣本矩陣,每個特徵的
維數下降了。
當給定一個測試的特徵集之後,比如1*10維的特徵,乘以上面得到的10*4的特徵矩陣,便可以得到一個1*4的特徵,用這個特徵去分類。
所以做PCA實際上是求得這個投影矩陣,用高維的特徵乘以這個投影矩陣,便可以將高維特徵的維數下降到指定的維數。
PCA 的目標是尋找 r ( r<n )個新變數,使它們反映事物的主要特徵,壓縮原有數據矩陣的規模。每個新變數是原有變數的
線性組合,體現原有變數的綜合效果,具有一定的實際含義。這 r 個新變數稱為“主成分”,它們可以在很大程度上反映原來 n 個變數的影響,並且這些新變數是互不相關的,也是正交的。通過
主成分分析,壓縮數據空間,將多元數據的特徵在低
維空間里直觀地表示出來。例如,將多個
時間點、多個實驗條件下的
基因表達譜數據( N 維)表示為 3維空間中的一個點,即將數據的維數從 RN 降到 R3 。
在進行基因表達數據分析時,一個重要問題是確定每個實驗數據是否是獨立的,如果每次實驗數據之間不是獨立的,則會影響基因表達數據分析結果的準確性。對於利用基因晶片所檢測到的基因表達數據,如果用 PCA 方法進行分析,可以將各個基因作為變數,也可以將實驗條件作為變數。當將基因作為變數時,通過分析確定一組“主要基因元素”,它們能夠很好地說明基因的特徵,解釋實驗現象;當將實驗條件作為變數時,通過分析確定一組“主要實驗因素”,它們能夠很好地刻畫實驗條件的特徵,解釋基因的行為。下面著重考慮以實驗條件作為變數的 PCA 分析方法。假設將數據的維數從 R N 降到 R 3 ,具體的 PCA 分析步驟如下:
(2) 第二步計算協方差矩陣S的本徵向量e1,e2,…,eN的本徵值, i = 1,2,…,N 。本徵值按大到小排序:;
(3)第三步投影數據到本徵矢
張成的空間之中,這些本徵矢相應的本徵值為。數據可以在
三維空間中展示為雲狀的點集。
對於 PCA ,確定新變數的個數 r 是一個兩難的問題。我們的目標是減小 r ,如果 r 小,則數據的
維數低,便於分析,同時也降低了噪聲,但可能丟失一些有用的信息。究竟如何確定 r 呢?這需要進一步分析每個
主元素對信息的貢獻。
貢獻率表示所定義的主成分在整個數據分析中承擔的主要意義占多大的比重,當取前 r 個主成分來代替原來全部變數時,累計貢獻率的大小反應了這種取代的可靠性,累計貢獻率越大,可靠性越大;反之,則可靠性越小。一般要求累計貢獻率達到 70% 以上。
經過 PCA 分析,一個多變數的複雜問題被簡化為低維空間的簡單問題。可以利用這種簡化方法進行作圖,形象地表示和分析複雜問題。在分析基因表達數據時,可以針對
基因作圖,也可以針對實驗條件作圖。前者稱為 Q 分析,後者稱為 R 分析。
PCA在matlab中的實現舉例
以下資料來自matlab的help,翻譯和註解部分由筆者添加:(重點部分添加了翻譯!)
函式名稱
Principal component analysis (PCA) on data
[COEFF,SCORE] = princomp(X)
[COEFF,SCORE,latent] = princomp(X)
[COEFF,SCORE,latent,tsquare] = princomp(X)
[...] = princomp(X,'econ')
函式描述
COEFF = princomp(X)performs principal components analysis (PCA) on the n-by-p data matrix X, and returns the principal component coefficients, also known as loadings. Rows of X correspond to observations, columns to variables. COEFF is a p-by-p matrix, each column containing coefficients for one principal component. The columns are in order of decreasing component variance.
在n行p列的數據集X上做主成分分析。返回主成分係數。X的每行表示一個樣本的觀測值,每一列表示特徵變數。COEFF是一個p行p列的矩陣,每一列包含一個主成分的係數,列是按主成分變數遞減順序排列。(按照這個翻譯很難理解,其實COEFF是X矩陣所對應的
協方差陣V的所有
特徵向量組成的矩陣,即
變換矩陣或稱投影矩陣,COEFF每列對應一個
特徵值的特徵向量,列的排列順序是按特徵值的大小遞減排序,後面有具體例子解釋,見
說明1)
princomp centers X by subtracting off column means, but does not rescale the columns of X. To perform principal components analysis with standardized variables, that is, based on correlations, use princomp(zscore(X)). To perform principal components analysis directly on a covariance or correlation matrix, use pcacov.
計算PCA的時候,MATLAB自動對列進行了去均值的操作,但是並不對數據進行規格化,如果要規格化的話,用princomp(zscore(X))。另外,如果直接有現成的
協方差陣,用函式pcacov來計算。
[COEFF,SCORE] = princomp(X)returns SCORE, the principal component scores; that is, the representation of X in the principal component space. Rows of SCORE correspond to observations, columns to components.
返回的SCORE是對主分的打分,也就是說原X矩陣在主成分空間的表示。SCORE每行對應樣本觀測值,每列對應一個主成份(變數),它的行和列的數目和X的行列數目相同。
[COEFF,SCORE,latent] = princomp(X)returns latent, a vector containing the eigenvalues of the covariance matrix of X.
返回的latent是一個向量,它是X所對應的協方差矩陣的特徵值向量。
[COEFF,SCORE,latent,tsquare] = princomp(X)returns tsquare, which contains Hotelling's T2 statistic for each data point.
返回的tsquare,是表示對每個樣本點Hotelling的T方統計量(我也不很清楚是什麼東東)。
The scores are the data formed by transforming the original data into the space of the principal components. The values of the vector latent are the variance of the columns of SCORE. Hotelling's T2 is a measure of the multivariate distance of each observation from the center of the data set.
所得的分(scores)表示由原數據X轉變到主成分空間所得到的數據。latent向量的值表示SCORE矩陣每列的方差(見說明2)。Hotelling的T方是用來衡量多變數間的距離,這個距離是指樣本觀測值到數據集中心的距離。
When n <= p, SCORE(:,n:p) and latent(n:p) are necessarily zero, and the columns of COEFF(:,n:p) define directions that are orthogonal to X.
[...] = princomp(X,'econ')returns only the elements of latent that are not necessarily zero, and the corresponding columns of COEFF and SCORE, that is, when n <= p, only the first n-1. This can be significantly faster when p is much larger than n.
當維數p超過樣本個數n的時候,用[...] = princomp(X,'econ')來計算,這樣會顯著提高計算速度
舉例
(下面樣本數據集為ingredients,matlab自帶)
Compute principal components for the ingredients data in the Hald data set, and the variance accounted for by each component.
load hald; 載入matlab內部數據
[pc,score,latent,tsquare] = princomp(ingredients); 調用pca分析函式
ingredients,score,pc,latent,tsquare 顯示得到的結果
ingredients =
7 26 6 60
1 29 15 52
11 56 8 20
11 31 8 47
7 52 6 33
11 55 9 22
3 71 17 6
1 31 22 44
2 54 18 22
21 47 4 26
1 40 23 34
11 66 9 12
10 68 8 12
score =
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
pc =
-0.0678 -0.6460 0.5673 0.5062
-0.6785 -0.0200 -0.5440 0.4933
0.0290 0.7553 0.4036 0.5156
0.7309 -0.1085 -0.4684 0.4844
latent =
517.7969
67.4964
12.4054
0.2372
tsquare =
5.6803
3.0758
6.0002
2.6198
3.3681
0.5668
3.4818
3.9794
2.6086
7.4818
4.1830
2.2327
2.7216
驗證
計算ingredients協方差矩陣:
cov_ingredients=cov(ingredients)
cov_ingredients =
34.6026 20.9231 -31.0513 -24.1667
20.9231 242.1410 -13.8782 -253.4167
-31.0513 -13.8782 41.0256 3.1667
-24.1667 -253.4167 3.1667 280.1667
下面為計算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和矩陣pc,易明白為什麼COEFF是按列遞減順序排列的。
下面再驗證說明2
diag(cov(score))
ans =
517.7969
67.4964
12.4054
0.2372
說明------2:以上結果顯示latent確實表示SCORE矩陣每列的方差,517.7969表示第一列方差
下面做圖表示結果:
我們要的是由函式[pc,score,latent,tsquare] = princomp(ingredients)所產生的pc和latent。由latent可以算出降維後的空間所能表示原空間的程度,只要這個累積的值大於95%。
The following command and plot show that two components account for 98% of the variance:
cumsum(latent)./sum(latent)
ans =
0.86597
0.97886
0.9996
1
%由以上ans值可以看出前兩個主成分就能表示原空間的97.886%,所以取pc中的前兩列可
%做主成分變換矩陣tranMatrix = pc(:,1:2)。則從原來的4
維空間降到2維空間。對任意一個
%原空間樣本,例如a=(7 ,26 ,6 ,60)變到低維空間的表達式為a1 = a*tranMatrix。(當然你也可
%以取pc中的前三列,由原來的4維空間變到3維空間)
biplot(pc(:,1:2),'Scores',score(:,1:2),'VarLabels',...
{'X1' 'X2' 'X3' 'X4'})