PCA（主成分分析法）的Python代码实现（numpy，sklearn）

PCA（主成分分析法）的Python代码实现（numpy，sklearn）

语言描述

PCA设法将原来众多具有一定相关性的属性（比如p个属性），重新组合成一组相互无关的综合属性来代替原属性。通常数学上的处理就是将原来p个属性做线性组合，作为新的综合属性。

PCA 中的线性变换等价于坐标变换，变换的目的是使

    n
   
  
  
   n
  
 
n 个样本点在新坐标轴 

 
  
   
    
     y
    
    
     1
    
   
  
  
   y_1
  
 
y1​ 上的离散程度（方差）最大，这样变量 

 
  
   
    
     y
    
    
     1
    
   
  
  
   y_1
  
 
y1​ 就代表了原始数据的绝大部分信息，即使忽略 

 
  
   
    
     y
    
    
     2
    
   
  
  
   y_2
  
 
y2​ 也无损大局，从而把两个指标压缩成一个指标。从几何上看，找主成分的问题就是找出 

 
  
   
    N
   
  
  
   N
  
 
N 维空间中椭球体的主轴问题。从数学上也可以证明，它们分别是相关矩阵的 

 
  
   
    k
   
  
  
   k
  
 
k 个较大的特征值所对应的特征向量。

算法描述

输入：

    n
   
  
  
   n
  
 
n 个 

 
  
   
    m
   
  
  
   m
  
 
m 维 

 
  
   
    (
   
   
    n
   
   
    ×
   
   
    m
   
   
    )
   
  
  
   (n \times m)
  
 
(n×m) 样本集 

 
  
   
    X
   
   
    =
   
   
    (
   
   
    
     x
    
    
     1
    
   
   
    ,
   
   
    
     x
    
    
     2
    
   
   
    ,
   
   
    ⋯
    
   
    ,
   
   
    
     x
    
    
     m
    
   
   
    )
   
  
  
   X = (x_1,x_2,\cdots,x_m)
  
 
X=(x1​,x2​,⋯,xm​)，低维空间维数 

 
  
   
    k
   
  
  
   k
  
 
k 。

主成分的计算步骤如下：

1. 对所有样本进行中心化： x i ← x i − 1 n ∑ i = 1 n x i x_i \leftarrow x_i-\frac{1}{n}\sum_{i=1}^{n}x_i xi​←xi​−n1​∑i=1n​xi​
2. 计算样本的协方差矩阵 X T X X^TX XTX ( m × m ) (m \times m) (m×m)
3. 计算特征值与特征向量- 解特征方程 ∣ λ E − X T X ∣ = 0 |\lambda E-X^TX|=0 ∣λE−XTX∣=0，常用雅可比 ( Jacobi ) 法求出特征值，并使其按大小顺序排列，即 λ 1 ≥ λ 2 ≥ ⋯ ≥ λ m ≥ 0 \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_m \geq 0 λ1​≥λ2​≥⋯≥λm​≥0 。- 分别求出对应于特征值 λ i \lambda_i λi​ 的特征向量 e i ( i = 1 , 2 , ⋯   , m ) e_i(i=1,2,\cdots,m) ei​(i=1,2,⋯,m)，要求 ∣ ∣ e i ∣ ∣ = 1 ||e_i|| = 1 ∣∣ei​∣∣=1，即 ∑ j = 1 m e i j 2 = 1 \sum_{j=1}^{m}e_{ij}^2=1 ∑j=1m​eij2​=1，其中 e i j e_{ij} eij​ 表示向量 e i e_i ei​ 的第 j j j 个分量。- 计算主成分贡献率及累计贡献率。 - 贡献率的公式： f i = λ i ∑ i = 1 m λ i f_i=\frac{\lambda_i}{\sum_{i=1}^{m}\lambda_i} fi​=∑i=1m​λi​λi​​- 累计贡献率： α k = ∑ i = 1 k f i \alpha_k=\sum_{i=1}^{k}f_i αk​=∑i=1k​fi​- 一般取累计贡献率达 85% ~ 95% 的特征值 λ 1 , λ 2 , ⋯   , λ l \lambda_1,\lambda_2,\cdots,\lambda_l λ1​,λ2​,⋯,λl​ 所对应的第一、第二、第 l ( l ≤ m ) l(l \leq m) l(l≤m)个主成分。
4. 计算主成分值前 k k k 个主成分值 z = ( X e 1 , X e 2 , ⋯   , X e k ) = ( z 1 , z 2 , ⋯   , z k ) z = (Xe_1,Xe_2,\cdots,Xe_k)=(z_1,z_2,\cdots,z_k) z=(Xe1​,Xe2​,⋯,Xek​)=(z1​,z2​,⋯,zk​) ( n × k 即 n 个 k 维 ) (n \times k\ 即\ n\ 个\ k\ 维) (n×k 即 n 个 k 维)

与通过保留原属性集的一个子集来减少属性集的大小不同，PCA通过创建一个能替换的、较小的变量集“组合”属性的基本要素。原数据可以投影到较小的集合中。PCA常常能够揭示先前未被察觉的联系，并允许解释不寻常的结果。

示例

学号语文数学物理化学英语历史18465617279812647777765570365676349576747480697563745847470807482

1 使用numpy降维

>>>import numpy as np

# 输入待降维数据 (5 * 6) 矩阵，6个维度，5个样本值>>> A = np.array([[84,65,61,72,79,81],[64,77,77,76,55,70],[65,67,63,49,57,67],[74,80,69,75,63,74],[84,74,70,80,74,82]])>>>print(A)[[846561727981][647777765570][656763495767][748069756374][847470807482]]# 对每一个属性的样本求均值>>> MEAN = np.mean(A, axis=0)# 沿轴0调用mean函数>>>print(MEAN)[74.272.668.70.465.674.8]# 去中心化>>> X = np.subtract(A, MEAN)>>>print(X)[[9.8-7.6-7.1.613.46.2][-10.24.49.5.6-10.6-4.8][-9.2-5.6-5.-21.4-8.6-7.8][-0.27.41.4.6-2.6-0.8][9.81.42.9.68.47.2]]>>>print(X.T)#矩阵的转置[[9.8-10.2-9.2-0.29.8][-7.64.4-5.67.41.4][-7.9.-5.1.2.][1.65.6-21.44.69.6][13.4-10.6-8.6-2.68.4][6.2-4.8-7.8-0.87.2]]# 计算协方差矩阵>>> COV = np.dot(X.T, X)>>>print(COV)[[380.8-55.6-95.248.6401.4252.2][-55.6165.2131.179.8-107.8-20.4][-95.131.160.170.-132.-34.][248.6179.8170.605.2214.8215.4][401.4-107.8-132.214.8443.2263.6][252.2-20.4-34.215.4263.6174.8]]# 计算特征值和特征向量>>> W, V = np.linalg.eig(COV)>>>print(W)# 特征值[1.22517276e+036.54041238e+023.95721181e+011.04138814e+011.50877843e-145.51899893e-14]>>>print(V)# 特征向量[[-0.532642530.20279107-0.344338060.39437042-0.61869481-0.55543331][0.00876193-0.46059524-0.815970780.021852320.258425160.34848844][0.04593605-0.473283850.378770770.70892582-0.031448860.21014772][-0.51955599-0.642385940.24891406-0.45230979-0.15412561-0.22434743][-0.551319360.327754780.09651389-0.130445260.294467280.67491022][-0.374451030.051452020.02970770.346148120.662554490.14160509]]# 计算主成分贡献率以及累计贡献率>>> sum_lambda = np.sum(W)# 特征值的和>>>print(sum_lambda)1929.1999999999994>>>f = np.divide(W, sum_lambda)# 每个特征值的贡献率（特征值 / 总和）>>>print(f)[6.35067780e-013.39021998e-012.05121906e-025.39803100e-037.82074656e-182.86077075e-17]>>> f[0]+f[1]# 前两大的特征值的累计贡献率0.974089778403108>>> f[0]+f[1]+f[2]# 前三大的特征值的累计贡献率0.9946019690025047# 0.97 > 0.85，只需要选取前两大特征值即可，以从6维降到2维# 前两大特征值对应的特征向量为：>>> e1 = V.T[0]>>>print(e1)[-0.532642530.008761930.04593605-0.51955599-0.55131936-0.37445103]>>> e2 = V.T[1]>>>print(e2)[0.20279107-0.46059524-0.47328385-0.642385940.327754780.05145202]# 计算主成分值（已去中心化）>>> z1 = np.dot(X, e1)>>>print(z1)[-16.1486052810.6167674323.40212697-0.43966353-17.43062559]>>> z2 = np.dot(X, e2)>>>print(z2)[12.48396235-15.6731742813.607117-7.77054621-2.64735885]# 输出降维后的结果（已去中心化）>>> RES = np.array([z1,z2])>>>print(RES)[[-16.1486052810.6167674323.40212697-0.43966353-17.43062559][12.48396235-15.6731742813.607117-7.77054621-2.64735885]]>>>print(RES.T)[[-16.1486052812.48396235][10.61676743-15.67317428][23.4021269713.607117][-0.43966353-7.77054621][-17.43062559-2.64735885]]

2 直接使用sklearn中的PCA进行降维

>>>import numpy as np
>>>from sklearn.decomposition import PCA

# 输入待降维数据 (5 * 6) 矩阵，6个维度，5个样本值>>> A = np.array([[84,65,61,72,79,81],[64,77,77,76,55,70],[65,67,63,49,57,67],[74,80,69,75,63,74],[84,74,70,80,74,82]])>>>print(A)[[846561727981][647777765570][656763495767][748069756374][847470807482]]# 直接使用PCA进行降维>>> pca = PCA(n_components=2)#降到 2 维>>> pca.fit(A)
PCA(n_components=2)>>> pca.transform(A)# 降维后的结果
array([[-16.14860528,-12.48396235],[10.61676743,15.67317428],[23.40212697,-13.607117],[-0.43966353,7.77054621],[-17.43062559,2.64735885]])>>> pca.explained_variance_ratio_ # 降维后的各主成分的方差值占总方差值的比例，即方差贡献率
array([0.63506778,0.339022])>>> pca.explained_variance_ # 降维后的各主成分的方差值
array([306.29319053,163.51030959])