【數學建模】day10-主成分分析

0.

關於主成分分析的詳細理解以及理論推導，這篇blog中講的很清楚。

主成分分析是一種經常使用手段。這應該與因子分析等區別開來，重點在於理解主成分分析的做用以及什麼狀況下使用主成分分析，本文重點講解如何使用PCA。

 

1.

主成分分析是一種降維方法。

實際上這個降維是這樣作的：原始變量有m維，PCA主成分變量有t維（t<m），那麼就至關於把這m維分別往t維上投影。

 

例如咱們要作迴歸分析，若是自變量衆多，彼此之間又具備複雜的相關性，那麼咱們考慮對自變量個數進行「減小」。而這個減小不可以丟失有效信息，由此可使用主成分分析。

 

主成分分析的主要思想是，對原始衆多自變量進行【線性組合】（這實際上是向PCA方向投影的過程），從而獲得新的PCA變量(個數一般比原始變量少)。對於每個新變量，其線性組合的係數向量叫該主成分的方向（重點在於如何獲得這個PCA方向）；不一樣主成分的方向是正交的，從而保證新的變量彼此不相關，消除了糅合性。

具體來講，是：

原始變量Xi(i = 1,2,…p)，主成分變量Z(z = 1,2…p)，則：

其中，要求知足：

 

注意：

　　1）主成分分析的結果受量綱的影響，因爲各變量的單位可能不同，若是各自改 變量綱，結果會不同，這是主成分分析的大問題，迴歸分析是不存在這種狀況的， 因此實際中能夠先把各變量的數據標準化，而後使用協方差矩陣或相關係數矩陣進行分析。（相關係數就是標準化變量的協方差） 

　　2）使方差達到大的主成分分析不用轉軸（因爲統計軟件常把主成分分析和因子 分析放在一塊兒，後者每每須要轉軸，使用時應注意）。 

　　3）主成分的保留。用相關係數矩陣求主成分時，Kaiser主張將特徵值小於1的主成 分予以放棄（這也是SPSS軟件的默認值）。 

　　4）在實際研究中，因爲主成分的目的是爲了降維，減小變量的個數，故通常選取 少許的主成分（不超過5或6個），只要它們能解釋變異的70％～80％（稱累積貢獻率） 就好了。 

 

2.

如何選取主成分？

或者說：重點在於找到主成分方向，由於有了PCA方向，把原始變量向PCA方向上投影就得新的PCA變量。

明確目的在於線性組合原始變量獲得：

（z是向量）

在係數平方和爲1下，使得z的方差最大。（這樣就使得新的變量z1，z2…差別最大，就表明咱們抓住了原始變量的大多數信息）；同時，保證各個主成分變量的係數矩陣c兩兩正交（這表明咱們使得新的變量彼此不糅合，不相關）。

 

簡單推導：

設原始變量數據：

（X稱做設計陣）每列表明一個變量指標，每行是一組數據。

 

首先將X標準化。

 

設新變量z = c1x1 + c2x2 +… + cpxp；

係數向量C = (c1,c2,…cp).T

 

係數向量C就是主成分的方向。

 

目的在於使得原始變量X的每一行，在C的方向投影后，在C的方向上達到方差最大。也就是：

　　Z = CX （Z是X向C投影）

　　max E（(Z-E(Z))‘* (Z-E(Z))）

這個求max能夠通過變化，轉換^{【具體看推導blog，僅僅使用PCA能夠略過】}爲求max E((X - E(x))' * (X - E(X)))

因爲X標準化了（歸一化過），E(X) = 0;

從而變爲 max E（X.T*X）

而X.T * X 是一個二次型的。

那麼，最大值就在：這個矩陣X.T * X 的最大特徵值λ對應的特徵向量c = [c1,c2...cp]處取得。

c還要通過單位化處理。

也就是說，這就獲得了第一主成分的方向c；把X向c上投影獲得了第一主成分Z1.

而這個半正定矩陣X.T *X的不一樣特徵向量是正交的，從而把第二特徵向量做爲第二主成分的方向。

以此類推。

 

這裏，X標準化後，X.T*X / (n-1)其實是原始變量x1,x2…xp的相關係數矩陣。（協方差矩陣）

求解時，其實是計算這個相關係數矩陣。

 

3.

PCA步驟：

　　1. 列出設計陣X，將X標準化處理。

　　2. 求相關係數矩陣R（標準化後，這也就是協方差矩陣）。這裏要知道，這個R能夠用：

X.T*X / (n-1)求得。也可使用MATLAB求相關係數矩陣命令。

　　3. 求R的特徵值（從大到小排列），以及對應的單位標準正交特徵向量。

 

　　3. 從特徵值大到小選擇PCA變量：先選取最大的特徵值以及其特徵向量，從而構成第一個主成分變量，這個特徵向量就是PCA方向。計算累計貢獻率：已選取的特徵值之和佔全部特徵值之和的比重。

　　4. 重複步驟3，直到累計貢獻率達標或者選取了足夠的PCA變量。

　　5. 單純考慮累積貢獻率有時是不夠的，還須要考慮選擇的主成分對原始變量的貢獻 值，咱們用相關係數的平方和來表示.若是選取的主成分爲 z1,z2,…zr ，則它們對原 變量 xi 的貢獻值爲:

　　　　　　　　　　pi = ∑ (r(zj,xi))^2；（注1）

　　　　　（即每個主成分變量與xi的相關係數平方的和。）

　　6. 進而咱們能夠用主成分變量對問題作出其餘分析（如迴歸分析等）。

 

注1：兩個向量的相關係數就是兩向量夾角的餘弦，展開來講就是：

 

4.

一個例子：

 

5.

MATLAB實現

函數使用：

1. 求矩陣X列向量間的相關係數矩陣： r = corrcoef(X) param: 　　X：這個矩陣X的每一列看作一個向量 return： 　　r： 相關係數矩陣。返回每兩列之間的相關係數組成的矩陣，對稱陣，對角爲方差，一個向量的方差就是對向量的每一個數求var。 　　注意使用時沒必要對設計陣作標準化處理。由於咱們實際要求的是標準化變量的協方差矩陣，而就是原設計陣的相關係數矩陣。

2. matlab作PCA分析 [vec1,lamda,rate] = pcacov(r) param： 　　r：原始數據的相關係數矩陣 return： 　　vec1: r的特徵向量 　　lamda: 對應的特徵值 　　rate: 各個主成分的貢獻率

3. 累計求和函數cumsum y = cumsum(x,axis) 對x矩陣累計求和 param: 　　axis：軸，axis = 1則對列向量累計求和 　　axis = 2則對行向量累計求和 return: 　　y: 累和矩陣

4. 標準化化處理函數 y = zccore(x) param: 　　x： 矩陣或者向量，標準化處理，如果矩陣，是對列向量標準化處理 return: 　　y： 處理後的矩陣或者向量

5. 主成分分析作線性迴歸最小二乘估計函數 [c,s,t] = princomp(x) param: 　　x是設計陣 return: 　　c:  對主成分變量作多元線性迴歸分析，迴歸方程的係數 　　s:  這個是作主成分分析獲得的特徵向量矩陣，每一列是一個特徵向量（單位化） 　　t：相應的特徵值

 

 

 

應用主成分分析+迴歸分析案例:

分別作主成分迴歸以及原始變量直接回歸。

數據保存：

data.txt

7 26 6 60 78.5 
 1 29 15 52 74.3 
 11 56 8 20 104.3 
 11 31 8 47 87.6  
 7 52 6 33 95.9 
 11 55 9 22 109.2 
 3 71 17 6 102.7  
 1 31 22 44 72.5  
 2 54 18 22 93.1 
 21 47 4 26 115.9 
 1 40 23 34 83.8 
 11 66 9 12 113.3  
 10 68 8 12 109.4

matlab求解以下：

（注：也能夠直接用princomp函數，更簡單，下面使用標準的pcacov函數）

 1 clc,clear  2  load data.txt %導入數據  3 [m,n] = size(data);  4  x0 = data(:,[1:n-1]);  5  y0 = data(:,n);  6  hg1 = [ones(m,1),x0] \ y0; %普通多元線性迴歸係數.列向量  7 hg1 = hg1' %顯示
 8 fprintf('y = % f',hg1(1));  9  for i = 2:n 10      if hg1(i) > 0
11          fprintf('+ % f*x% d',hg1(i),i-1); 12      else
13          fprintf('% f*x% d',hg1(i),i-1); 14  end 15  end 16 
17 fprintf('\n'); 18  r = corrcoef(x0) %相關係數矩陣 19 xd = zscore(x0); 20  yd = zscore(y0); %標準化處理 21 [vec1,lamda,rate] = pcacov(r) %PCA 22  f = repmat(sign(sum(vec1)),size(vec1,1),1); 23 %產生與vec1同維數的元素爲+1/-1的矩陣 24 vec2 = vec1.*f %修改特徵值正負號，使得特徵值的全部份量和爲+
25 contr = cumsum(rate) %計算累計貢獻率 26 df = xd * vec2; %計算全部主成分的得分 27 num = input('請輸入主成分個數：') 28 hg21 = df(:,[1:num]) \ yd %主成分變量回歸係數 29 hg22 = vec2(:,1:num) * hg21 %標準化變量的迴歸係數 30 hg23 = [mean(y0) - std(y0)*mean(x0)./std(x0)*hg22,std(y0)*hg22'./std(x0)]
31  % 轉換，求原始變量的迴歸方程係數 32 
33 fprintf('y = % f',hg23(1)); 34  for i = 2:n 35      if hg23(i) > 0
36          fprintf('+ % f*x% d',hg23(i),i-1); 37      else
38          fprintf('% f*x% d',hg23(i),i-1); 39  end 40  end 41  fprintf('\n') 42  %下面計算兩種迴歸分析的剩餘標準差 43 rmse1=sqrt(sum((x0*hg1(2:end)'+hg1(1)-y0).^2)/(m-n)) 
44  rmse2=sqrt(sum((x0*hg23(2:end)'+hg23(1)-y0).^2)/(m-num)) 

結果以下：

相關係數矩陣：

PCA分析獲得的：特徵值、特徵向量、貢獻率等結果：

作代碼所示處理後的特徵向量：

累計貢獻率：

能夠看到，最後一個變量幾乎沒有貢獻，使用前三個PCA成分獲得:（其實兩個也能夠）

其中，hg21是主成分變量回歸係數、hg22是標準化X的迴歸係數(沒有常數項)，hg23是原始X的迴歸係數（第一個是常數項）

從而，PCA後的結果爲：

而在原始變量上直接作迴歸分析的結果爲：

 

代碼中還計算了兩個剩餘標準差：

可見，PCA後再回歸的剩餘標準差要小，效果更好，這也證實了使用PCA是有效的。