成千累万网

主成分分析法(PCA)及其python实现

主成分分析法(PCA)及其python实现

主成分分析法(Principal Component Analysis,主成PCA)是分分一种用于把高维数据降成低维,使分析变得更加简便的析法现分析方法。比如我们的主成一个样本可以由 n n n维随机变量 ( X 1 , X 2 , . . . , X n ) (X_1,X_2,...,X_n) (X1​,X2​,...,Xn​)来刻画,运用主成分分析法,分分我们可以把这些分量用更少的析法现、这 n n n个分量的主成线性组合来表示。本文多为学习后的分分个人理解,如有错误还请指出。析法现

基本思想

我们把降维后的主成变量称为主成分(Principal Component),设其为 Z 1 ,分分 Z 2 , . . . , Z n Z_1,Z_2,...,Z_n Z1​,Z2​,...,Zn​(我们并不取这全部的 n n n个变量,否则降维就没有意义了)。析法现 Z i Z_i Zi​称为第 i i i主成分。主成每个主成分都是分分原来 n n n个变量的线性组合,即
{ Z 1 = a 11 X 1 + a 12 X 2 + . . . + a 1 n X n Z 2 = a 21 X 1 + a 22 X 2 + . . . + a 2 n X n . . . Z n = a n 1 X 1 + a n 2 X 2 + . . . + a n n X n \begin{ cases}Z_1=a_{ 11}X_1+a_{ 12}X_2+...+a_{ 1n}X_n\\ Z_2=a_{ 21}X_1+a_{ 22}X_2+...+a_{ 2n}X_n \\ ...\\Z_n=a_{ n1}X_1+a_{ n2}X_2+...+a_{ nn}X_n \\ \end{ cases} ⎩ ⎨ ⎧​Z1​=a11​X1​+a12​X2​+...+a1n​Xn​Z2​=a21​X1​+a22​X2​+...+a2n​Xn​...Zn​=an1​X1​+an2​X2​+...+ann​Xn​​
或者写成更简便的析法现线性代数形式:设 Z = ( Z 1 Z 2 ⋮ Z n ) , X = ( X 1 X 2 ⋮ X n ) , A = ( a 11 ⋯ a 1 n ⋮ ⋱ ⋮ a n 1 ⋯ a n n ) , Z=\begin{ pmatrix}Z_1\\Z_2\\ \vdots \\Z_n\end{ pmatrix},X=\begin{ pmatrix}X_1\\X_2\\ \vdots \\X_n\end{ pmatrix},A=\begin{ pmatrix}a_{ 11} & \cdots & a_{ 1n} \\ \vdots & \ddots & \vdots \\ a_{ n1} & \cdots & a_{ nn}\end{ pmatrix}, Z=⎝ ⎛​Z1​Z2​⋮Zn​​⎠ ⎞​,X=⎝ ⎛​X1​X2​⋮Xn​​⎠ ⎞​,A=⎝ ⎛​a11​⋮an1​​⋯⋱⋯​a1n​⋮ann​​⎠ ⎞​,则这个关系可被表示为 Z = A X . Z=AX. Z=AX.
为了达到降维的目的,我们需要保证

( 1 ) Z 1 , Z 2 , . . . , Z n (1)Z_1,Z_2,...,Z_n (1)Z1​,Z2​,...,Zn​是线性无关的,这要求 A A A是正交阵。
如果存在线性相关的关系( ∃ \exists ∃ 不全为 0 0 0的 k 1 , k 2 , . . . , k n k_1,k_2,...,k_n k1​,k2​,...,kn​使得 k 1 Z 1 + k 2 Z 2 + . . . + k n Z n = 0 k_1Z_1+k_2Z_2+...+k_nZ_n=0 k1​Z1​+k2​Z2​+...+kn​Zn​=0),我们得到的结果中就存在着冗余信息(某个主成分可以由其它主成分表示),这种情况应该被排除。

( 2 ) (2) (2)选出 n n n个主成分中能显著代表原本变量的 k ( k < n ) k(k < n) k(k这一条提出了一个问题:我们按照什么标准来衡量主成分的好坏关系呢?统计学认为,一组数据越分散,它的方差越大,它所包含的信息就越多。(知乎的这个问题讨论了这一点)因此,PCA选出这 n n n个主成分中方差最大的 k k k个作为新的变量。

数学推导

设我们的 m m m个样本对应的数据矩阵为 R = ( r 11 ⋯ r 1 m ⋮ ⋱ ⋮ r n 1 ⋯ r n m ) R=\begin{ pmatrix}r_{ 11} & \cdots & r_{ 1m}\\ \vdots & \ddots & \vdots \\ r_{ n1} & \cdots & r_{ nm}\end{ pmatrix} R=⎝ ⎛​r11​⋮rn1​​⋯⋱⋯​r1m​⋮rnm​​⎠ ⎞​,每列代表一个样本;我们先将其标准化,使每行的均值为 0 0 0,并消除量纲的影响,便于进一步处理:

x i j = ( r i j − r i ˉ ) / s i , x_{ ij}=(r_{ ij}-\bar{ r_i})/s_i, xij​=(rij​−ri​ˉ​)/si​,

其中 r i ˉ \bar{ r_i} ri​ˉ​和 s i s_i si​分别为第 i i i行的均值和样本标准差,记处理后的矩阵为 X = ( x i j ) n × m ; X=(x_{ ij})_{ n\times m}; X=(xij​)n×m​;对该矩阵做线性变换的结果为 F = A X = ( z i j ) n × m . F=AX=(z_{ ij})_{ n\times m}. F=AX=(zij​)n×m​.

我们用 D ( Z ) D(Z) D(Z)表示 n n n维随机变量 Z = ( Z 1 Z 2 ⋮ Z n ) Z=\begin{ pmatrix}Z_1\\Z_2\\ \vdots \\Z_n\end{ pmatrix} Z=⎝ ⎛​Z1​Z2​⋮Zn​​⎠ ⎞​的协方差矩阵,那么有

D ( Z ) = ( C o v ( Z 1 , Z 1 ) ⋯ C o v ( Z 1 , Z n ) ⋮ ⋱ ⋮ C o v ( Z n , Z 1 ) ⋯ C o v ( Z n , Z n ) ) D(Z)=\begin{ pmatrix} Cov(Z_1,Z_1) & \cdots & Cov(Z_1,Z_n)\\ \vdots & \ddots & \vdots \\ Cov(Z_n,Z_1) & \cdots & Cov(Z_n,Z_n)\end{ pmatrix} D(Z)=⎝ ⎛​Cov(Z1​,Z1​)⋮Cov(Zn​,Z1​)​⋯⋱⋯​Cov(Z1​,Zn​)⋮Cov(Zn​,Zn​)​⎠ ⎞​

由上一部分的条件 ( 1 ) (1) (1),应当有 Z i , Z j ( i ≠ j ) Z_i,Z_j(i\neq j) Zi​,Zj​(i=j)线性不相关,即 C o v ( Z i , Z j ) = 0 Cov(Z_i,Z_j)=0 Cov(Zi​,Zj​)=0.而 C o v ( Z i , Z i ) = D ( Z i ) Cov(Z_i,Z_i)=D(Z_i) Cov(Zi​,Zi​)=D(Zi​),即 Z i Z_i Zi​的方差(注意跟上文的 D ( Z ) D(Z) D(Z)的记号意义不同,上面的 D D D指的是协方差矩阵),因此 D ( Z i ) D(Z_i) D(Zi​)就是一个对角矩阵,即

D ( Z ) = ( D ( Z 1 ) D ( Z 2 ) ⋱ D ( Z n ) ) D(Z)=\begin{ pmatrix}D(Z_1) & & \\ & D(Z_2) & \\ & & \ddots & \\ & & & D(Z_n)\end{ pmatrix} D(Z)=⎝ ⎛​D(Z1​)​D(Z2​)​⋱​D(Zn​)​⎠ ⎞​

由协方差的定义, C o v ( Z i , Z j ) = E ( Z i Z j ) − E ( Z i ) E ( Z j ) ; Cov(Z_i,Z_j)=E(Z_iZ_j)-E(Z_i)E(Z_j); Cov(Zi​,Zj​)=E(Zi​Zj​)−E(Zi​)E(Zj​);

由于 X X X经我们处理,对任意的 X i X_i Xi​都有 E ( X i ) = 0 E(X_i)=0 E(Xi​)=0,故 E ( Z i ) = E ( a i 1 X 1 + . . . + a i n X n ) = 0 , E(Z_i)=E(a_{ i1}X_1+...+a_{ in}X_n)=0, E(Zi​)=E(ai1​X1​+...+ain​Xn​)=0,

C o v ( Z i , Z j ) = E ( Z i Z j ) = 1 m ∑ t = 1 m z i t z j t . Cov(Z_i,Z_j)=E(Z_iZ_j)=\frac{ 1}{ m}\sum\limits_{ t=1}^{ m}z_{ it}z_{ jt}. Cov(Zi​,Zj​)=E(Zi​Zj​)=m1​t=1∑m​zit​zjt​.那么上面的 D ( Z ) D(Z) D(Z)还可表示为

D ( Z ) = ( C o v ( Z 1 , Z 1 ) ⋯ C o v ( Z 1 , Z n ) ⋮ ⋱ ⋮ C o v ( Z n , Z 1 ) ⋯ C o v ( Z n , Z n ) ) = ( D ( Z 1 ) D ( Z 2 ) ⋱ D ( Z n ) ) = 1 m F F T . D(Z)=\begin{ pmatrix} Cov(Z_1,Z_1) & \cdots & Cov(Z_1,Z_n)\\ \vdots & \ddots & \vdots \\ Cov(Z_n,Z_1) & \cdots & Cov(Z_n,Z_n)\end{ pmatrix}=\begin{ pmatrix}D(Z_1) & & \\ & D(Z_2) & \\ & & \ddots & \\ & & & D(Z_n)\end{ pmatrix}=\frac{ 1}{ m}FF^T. D(Z)=⎝ ⎛​Cov(Z1​,Z1​)⋮Cov(Zn​,Z1​)​⋯⋱⋯​Cov(Z1​,Zn​)⋮Cov(Zn​,Zn​)​⎠ ⎞​=⎝ ⎛​D(Z1​)​D(Z2​)​⋱​D(Zn​)​⎠ ⎞​=m1​FFT.

又因为 F = A X , D ( Z ) = 1 m ( A X ) ( A X ) T = 1 m ( A X ) ( X T A T ) = A ( 1 m X X T ) A T = A D ( X ) A T F=AX,D(Z)=\frac{ 1}{ m}(AX)(AX)^T=\frac{ 1}{ m}(AX)(X^TA^T)=A(\frac{ 1}{ m}XX^T)A^T=AD(X)A^T F=AX,D(Z)=m1​(AX)(AX)T=m1​(AX)(XTAT)=A(m1​XXT)AT=AD(X)AT,这里 D ( X ) D(X) D(X)表示 n n n维随机变量 ( X 1 X 2 ⋮ X n ) \begin{ pmatrix}X_1\\X_2\\ \vdots \\X_n\end{ pmatrix} ⎝ ⎛​X1​X2​⋮Xn​​⎠ ⎞​的协方差矩阵。

由于 D ( X ) D(X) D(X)是实对称矩阵(这一点由协方差矩阵的定义即得: C o v ( X i , X j ) = C o v ( X j , X i ) Cov(X_i,X_j)=Cov(X_j,X_i) Cov(Xi​,Xj​)=Cov(Xj​,Xi​)),由实对称矩阵的性质, D ( X ) D(X) D(X)一定可以相似对角化,即存在正交矩阵 P , P, P,使得

P T D ( X ) P = ( λ 1 λ 2 ⋱ λ n ) ( 1 ) P^TD(X)P=\begin{ pmatrix}\lambda_1 & & & \\ & \lambda_2 & & \\ & & \ddots & \\ & & & \lambda_n\end{ pmatrix}\hspace{ 4.5cm}(1) PTD(X)P=⎝ ⎛​λ1​​λ2​​⋱​λn​​⎠ ⎞​(1)

其中 λ 1 . . . λ n \lambda_1...\lambda_n λ1​...λn​为 D ( X ) D(X) D(X)的 n n n个特征值。

这时候我们取出上面得到的式子

D ( Z ) = ( D ( Z 1 ) D ( Z 2 ) ⋱ D ( Z n ) ) = A D ( X ) A T ( 2 ) D(Z)=\begin{ pmatrix}D(Z_1) & & \\ & D(Z_2) & \\ & & \ddots & \\ & & & D(Z_n)\end{ pmatrix}=AD(X)A^T\hspace{ 1.5cm}(2) D(Z)=⎝ ⎛​D(Z1​)​D(Z2​)​⋱​D(Zn​)​⎠ ⎞​=AD(X)AT(2)

由基本思想部分的前提, A A A应当是正交矩阵,于是我们得到 ( D ( Z 1 ) D ( Z 2 ) ⋱ D ( Z n ) ) ∼ D ( X ) . \begin{ pmatrix}D(Z_1) & & \\ & D(Z_2) & \\ & & \ddots & \\ & & & D(Z_n)\end{ pmatrix} \sim D(X). ⎝ ⎛​D(Z1​)​D(Z2​)​⋱​D(Zn​)​⎠ ⎞​∼D(X).

由线性代数定理(若 n n n阶方阵 A A A与对角矩阵 D D D相似,则 D D D对角线上的元素即为 A A A的 n n n个特征值), λ 1 . . . λ n \lambda_1...\lambda_n λ1​...λn​即 D ( Z 1 ) . . . D ( Z n ) . D(Z_1)...D(Z_n). D(Z1​)...D(Zn​).更巧妙的是,我们可以求得变换矩阵 A = P T A=P^T A=PT,即将 D ( X ) D(X) D(X)相似对角化所需的正交阵之转置。得到了这个变换矩阵,我们就能得到 n n n个变换后的主成分了。具体地说,若

A = P T = ( a 11 ⋯ a 1 n ⋮ ⋱ ⋮ a n 1 ⋯ a n n ) , A=P^T=\begin{ pmatrix}a_{ 11} & \cdots & a_{ 1n} \\ \vdots & \ddots & \vdots \\ a_{ n1} & \cdots & a_{ nn}\end{ pmatrix}, A=PT=⎝ ⎛​a11​⋮an1​​⋯⋱⋯​a1n​⋮ann​​⎠ ⎞​,
则第 i i i个主成分 Z i = a i 1 X 1 + a i 2 X 2 + . . . + a i n X n : Z_i=a_{ i1}X_1+a_{ i2}X_2+...+a_{ in}X_n: Zi​=ai1​X1​+ai2​X2​+...+ain​Xn​: X i X_i Xi​前面的系数即 A A A的第 i i i行, P P P的第 i i i列,正好是 D ( X ) D(X) D(X)的第 i i i个特征值对应的特征向量(指的是相似对角化矩阵 P P P中标准的特征向量)。

推导结论

经过上面一大串推导,我们得到如下结论:

若 X X X是标准化处理过的数据矩阵,那么 D ( X ) D(X) D(X)的 n n n个特征值 λ 1 , λ 2 , . . . , λ n \lambda_1,\lambda_2,...,\lambda_n λ1​,λ2​,...,λn​即为线性变换后的 n n n个随机变量(即我们提到的主成分) Z 1 , Z 2 , . . . , Z n Z_1,Z_2,...,Z_n Z1​,Z2​,...,Zn​的方差 D ( Z 1 ) , D ( Z 2 ) , . . . , D ( Z n ) D(Z_1),D(Z_2),...,D(Z_n) D(Z1​),D(Z2​),...,D(Zn​);线性变换所需矩阵 A = P T A=P^T A=PT,其中 P P P为将 D ( X ) D(X) D(X)相似对角化所需的正交阵(由线性代数知识,这也是 D ( X ) D(X) D(X)的 n n n个特征向量组成的矩阵)。

这个结论使我们能够求出线性变换所需要的矩阵 A A A;此外我们可以根据特征值将 n n n个主成分排序,求得方差最大的 k k k个主成分。更具体地,求排序后的 n n n个主成分的算法如下:

1. 将原始数据矩阵 R 标准化 : x i j = ( r i j − r i ˉ ) / s i , 得到矩阵 X ; 1. 将原始数据矩阵R标准化:x_{ ij}=(r_{ ij}-\bar{ r_i})/s_i,得到矩阵X; 1.将原始数据矩阵R标准化:xij​=(rij​−ri​ˉ​)/si​,得到矩阵X;
2. 求 X 的协方差矩阵 D ( X ) ; 2.求X的协方差矩阵D(X); 2.求X的协方差矩阵D(X);
3. 求 D ( X ) 的特征值 λ 1 . . . λ n 和将其相似对角化需要的正交矩阵 P ; 3.求D(X)的特征值\lambda_1...\lambda_n和将其相似对角化需要的正交矩阵P; 3.求D(X)的特征值λ1​...λn​和将其相似对角化需要的正交矩阵P;
4. 方差第 i 大的主成分系数即第 i 大特征值对应的(单位化了的)特征向量 . 4.方差第i大的主成分系数即第i大特征值对应的(单位化了的)特征向量. 4.方差第i大的主成分系数即第i大特征值对应的(单位化了的)特征向量.

举个例子:(数据来自https://zhuanlan.zhihu.com/p/454447043)
我们要将如下的数据中5个变量(设能力,品格,担保,资本,环境为 X 1 . . . X 5 X_1...X_5 X1​...X5​)降维:
在这里插入图片描述
首先我们写出它对应的原始数据矩阵:
R = ( 66 65 57 ⋯ 62 64 64 63 58 ⋯ 63 66 ⋮ ⋮ 65 64 66 ⋯ 66 67 ) 5 × 15 R=\begin{ pmatrix}66 & 65 & 57 & \cdots & 62 & 64\\ 64 & 63 & 58 & \cdots & 63 & 66 \\ \vdots & & & & & \vdots\\ 65 & 64 & 66 & \cdots &66 & 67\end{ pmatrix}_{ 5 \times 15} R=⎝ ⎛​6664⋮65​656364​575866​⋯⋯⋯​626366​6466⋮67​⎠ ⎞​5×15​
然后将其标准化:
X = ( 0.7200823 ⋯ 0.0 − 0.06996503 ⋯ 0.62968524 ⋮ ⋮ 0.29580399 ⋯ 1.77482393 ) 5 × 15 X=\begin{ pmatrix} 0.7200823& \cdots & 0.0\\ -0.06996503 & \cdots &0.62968524 \\ \vdots & &\vdots\\ 0.29580399 & \cdots &1.77482393\end{ pmatrix}_{ 5 \times 15} X=⎝ ⎛​0.7200823−0.06996503⋮0.29580399​⋯⋯⋯​0.00.62968524⋮1.77482393​⎠ ⎞​5×15​
求出 X X X的协方差矩阵:
D ( X ) = ( 1.0 ⋯ 0.01901814 0.8816601 1.0 ⋯ 0.20695934 ⋮ ⋮ 0.01901814 ⋯ 1.0 ) 5 × 5 D(X)=\begin{ pmatrix} 1.0& \cdots & &0.01901814\\ 0.8816601 & 1.0 & \cdots &0.20695934 \\ \vdots & & &\vdots\\ 0.01901814 & & \cdots &1.0\end{ pmatrix}_{ 5 \times 5} D(X)=⎝ ⎛​1.00.8816601⋮0.01901814​⋯1.0​⋯⋯​0.019018140.20695934⋮1.0​⎠ ⎞​5×5​
求出 D ( X ) D(X) D(X)的特征值(从大到小排列)和特征向量:
λ 1 = 3.45317841 x 1 = ( 0.48198 0.51227 0.45384 0.51336 0.18914 ) T \lambda_1=3.45317841 \hspace{ 1cm}\pmb x_1=\begin{ pmatrix}0.48198 &0.51227 &0.45384& 0.51336& 0.18914\end{ pmatrix}^T λ1​=3.45317841xx1​=(0.48198​0.51227​0.45384​0.51336​0.18914​)T
λ 2 = 1.22308928 x 2 = ( 0.33297 0.13247 − 0.39212 0.20476 − 0.82213 ) T \lambda_2=1.22308928 \hspace{ 1cm}\pmb x_2=\begin{ pmatrix}0.33297 &0.13247 &-0.39212 &0.20476 &-0.82213\end{ pmatrix}^T λ2​=1.22308928xx2​=(0.33297​0.13247​−0.39212​0.20476​−0.82213​)T
λ 3 = 0.17872745 x 3 = ( 0.42459 0.1072 − 0.72892 − 0.05405 0.52344 ) T \lambda_3=0.17872745 \hspace{ 1cm}\pmb x_3=\begin{ pmatrix}0.42459 & 0.1072& -0.72892& -0.05405 & 0.52344\end{ pmatrix}^T λ3​=0.17872745xx3​=(0.42459​0.1072​−0.72892​−0.05405​0.52344​)T
λ 4 = 0.09923816 x 4 = ( − 0.39138 0.84166 − 0.11708 − 0.34902 − 0.05398 ) T \lambda_4=0.09923816 \hspace{ 1cm}\pmb x_4=\begin{ pmatrix}-0.39138 & 0.84166 &-0.11708 &-0.34902 &-0.05398\end{ pmatrix}^T λ4​=0.09923816xx4​=(−0.39138​0.84166​−0.11708​−0.34902​−0.05398​)T
λ 5 = 0.0457667 x 5 = ( − 0.56866 0.01252 − 0.3086 0.75485 0.1069 ) T \lambda_5=0.0457667 \hspace{ 1cm}\pmb x_5=\begin{ pmatrix}-0.56866 &0.01252& -0.3086 & 0.75485 & 0.1069\end{ pmatrix}^T λ5​=0.0457667xx5​=(−0.56866​0.01252​−0.3086​0.75485​0.1069​)T

我们怎么确定最终取出几个主成分呢?一般认为当取出的 k k k个主成分方差贡献比例之和达到 85 % 85\% 85%时就可以较好地代替原来的 n n n个变量了。因此我们还需要计算每个特征值所对应的方差贡献比例。由于 λ i = D ( Z i ) , \lambda_i=D(Z_i), λi​=D(Zi​),特征值 λ i \lambda_i λi​的方差占比即 λ i ∑ j = 1 n λ i . \frac{ \lambda_i}{ \sum\limits_{ j=1}^{ n}\lambda_i}. j=1∑n​λi​λi​​.按照上述方法,计算方差占比如下:

特征值 λ 1 \lambda_1 λ1​的方差贡献率 0.69064 0.69064 0.69064,累计方差贡献率 0.69064 ; 0.69064; 0.69064;
特征值 λ 2 \lambda_2 λ2​的方差贡献率 0.24462 0.24462 0.24462,累计方差贡献率 0.93525 ; 0.93525; 0.93525;(已到达 85 % 85\% 85%)
特征值 λ 3 \lambda_3 λ3​的方差贡献率 0.03575 0.03575 0.03575,累计方差贡献率 0.971 ; 0.971; 0.971;
特征值 λ 4 \lambda_4 λ4​的方差贡献率 0.01985 0.01985 0.01985,累计方差贡献率 0.99085 ; 0.99085; 0.99085;
特征值 λ 5 \lambda_5 λ5​的方差贡献率 0.00915 0.00915 0.00915,累计方差贡献率 1.0. 1.0. 1.0.

可以看到前 2 2 2个特征值对应的主成分即达到了 85 % 85\% 85%的方差贡献率,因此我们可以把原来的 5 5 5个变量“浓缩”表示为下面的 2 2 2个主成分(系数即为特征值对应的特征向量的各分量):

z 1 = 0.48198 x 1 + 0.51227 x 2 + 0.45384 x 3 + 0.51336 x 4 + 0.18914 x 5 ; z_1=0.48198x_1+0.51227x_2+0.45384x_3+0.51336x_4+0.18914x_5; z1​=0.48198x1​+0.51227x2​+0.45384x3​+0.51336x4​+0.18914x5​;
z 2 = 0.33297 x 2 + 0.13247 x 2 − 0.39212 x 3 + 0.20476 x 4 − 0.82213 x 5 . z_2=0.33297x_2+0.13247x_2-0.39212x_3+0.20476x_4-0.82213x_5. z2​=0.33297x2​+0.13247x2​−0.39212x3​+0.20476x4​−0.82213x5​.
这样,我们就成功实现了降维,以后我们分析数据的时候就可以用 z 1 z_1 z1​和 z 2 z_2 z2​来代替原来的五个指标了。

python实现

下面我们用python来实现一下上述的过程。

import numpy as npfrom numpy import linalg class PCA:    '''    dataset 形如array([样本1,样本2,...,样本m]),每个样本是一个n维的ndarray    '''    def __init__(self, dataset):    	# 这里的参数跟上文是反着来的(每行是一个样本),需要转置一下        self.dataset = np.matrix(dataset, dtype='float64').T    '''    求主成分;    threshold可选参数表示方差累计达到threshold后就不再取后面的特征向量.    '''    def principal_comps(self, threshold = 0.85):    	# 返回满足要求的特征向量        ret = []        data = []		# 标准化        for (index, line) in enumerate(self.dataset):            self.dataset[index] -= np.mean(line)            # np.std(line, ddof = 1)即样本标准差(分母为n - 1)            self.dataset[index] /= np.std(line, ddof = 1)        # 求协方差矩阵        Cov = np.cov(self.dataset)		# 求特征值和特征向量        eigs, vectors = linalg.eig(Cov)		# 第i个特征向量是第i列,为了便于观察将其转置一下        for i in range(len(eigs)):            data.append((eigs[i], vectors[:, i].T))        # 按照特征值从大到小排序        data.sort(key = lambda x: x[0], reverse = True)        sum = 0        for comp in data:            sum += comp[0] / np.sum(eigs)            ret.append(                tuple(map(                	# 保留5位小数                    lambda x: np.round(x, 5),                    # 特征向量、方差贡献率、累计方差贡献率                    (comp[1], comp[0] / np.sum(eigs), sum)                ))            )            print('特征值:', comp[0], '特征向量:', ret[-1][0], '方差贡献率:', ret[-1][1], '累计方差贡献率:', ret[-1][2])            if sum >threshold:                return ret              return ret 

测试一下刚才的例子:

p = PCA([[66, 64, 65, 65, 65], [65, 63, 63, 65, 64], [57, 58, 63, 59, 66], [67, 69, 65, 68, 64], [61, 61, 62, 62, 63], [64, 65, 63, 63, 63], [64, 63, 63, 63, 64], [63, 63, 63, 63, 63], [65, 64, 65, 66, 64], [67, 69, 69, 68, 67], [62, 63, 65, 64, 64], [68, 67, 65, 67, 65], [65, 65, 66, 65, 64], [62, 63, 64, 62, 66], [64, 66, 66, 65, 67]])lst = p.principal_comps()print(lst)

输出结果:

# 这部分是运行时输出的特征值: 3.4531784074578318 特征向量: [0.48198 0.51227 0.45384 0.51336 0.18914] 方差贡献率: 0.69064 累计方差贡献率: 0.69064特征值: 1.2230892804516 特征向量: [ 0.33297  0.13247 -0.39212  0.20476 -0.82213] 方差贡献率: 0.24462 累计方差贡献率: 0.93525# 这部分是返回的结果,为了美观稍微调整了一下格式[(array([0.48198, 0.51227, 0.45384, 0.51336, 0.18914]), 0.69064, 0.69064), (array([ 0.33297,  0.13247, -0.39212,  0.20476, -0.82213]), 0.24462, 0.93525)]

我这里设置的返回结果是一个三元组 ( 特征向量 , 方差贡献率 , 累计方差贡献率 ) (特征向量,方差贡献率,累计方差贡献率) (特征向量,方差贡献率,累计方差贡献率),如果有需要也可以自己调整一下返回的结果格式。此外,通过调整threshold可选参数可以设置累计方差贡献率到达多少时停止取主成分向量(默认为0.85)。

使用sklearn的PCA模块实现

这种经典的算法sklearn库也已经帮我们实现好了。对于上面的例子,其等价的使用sklearn库的代码如下:

import numpy as npfrom sklearn.decomposition import PCAX = np.array([[66, 64, 65, 65, 65], [65, 63, 63, 65, 64], [57, 58, 63, 59, 66], [67, 69, 65, 68, 64], [61, 61, 62, 62, 63], [64, 65, 63, 63, 63], [64, 63, 63, 63, 64], [63, 63, 63, 63, 63], [65, 64, 65, 66, 64], [67, 69, 69, 68, 67], [62, 63, 65, 64, 64], [68, 67, 65, 67, 65], [65, 65, 66, 65, 64], [62, 63, 64, 62, 66], [64, 66, 66, 65, 67]])# n_components 指明了降到几维pca = PCA(n_components = 2)# 利用数据训练模型(即上述得出特征向量的过程)pca.fit(X)# 得出原始数据的降维后的结果;也可以以新的数据作为参数,得到降维结果。print(pca.transform(X))# 打印各主成分的方差占比print(pca.explained_variance_ratio_)

运行结果:

[[-1.51394918 -0.21382815] [ 0.25137676 -1.8134245 ] [10.61577071  2.68155382] [-6.48520841 -1.16575919] [ 5.53026102 -1.52083322] [ 0.70154125 -1.8544697 ] [ 1.82460091 -1.29624147] [ 2.44281085 -1.60484093] [-1.40146605 -0.59189041] [-7.76925956  3.34817657] [ 1.8850487   0.61749314] [-5.41819247 -0.9163256 ] [-1.764172    0.155228  ] [ 3.06230672  1.51679123] [-1.96146925  2.65837042]]# 下面是方差贡献率[0.82399563 0.11748567]

我们发现这个方差贡献率(也就是特征值占比)跟我们手写的很不一样。(手写的是[0.69064, 0.24462])。这里可能会出现分歧的地方就是是否对原始数据除以样本标准差。当我把手写代码部分的

self.dataset[index] /= np.std(line, ddof = 1)

这一行注释掉后,发现运行结果与sklearn库的基本一致了:

特征值: 22.075235372070864 特征向量: [0.56177 0.58975 0.27868 0.50573 0.05644] 方差贡献率: 0.824 累计方差贡献率: 0.824特征值: 3.1474971323292125 特征向量: [ 0.37826 -0.06431 -0.61334  0.06946 -0.68686] 方差贡献率: 0.11749 累计方差贡献率: 0.94148

因此可以得出sklearn库在训练时似乎没有消除量纲,即没有对数据除以其样本标准差。当然,这仅仅是个人理解,不过与sklearn库结果基本吻合大致上印证了这个猜想。在一篇文章里我找到了关于量纲的讨论:若各个变量的单位一致,则各个属性是可以比较的,此时可以直接求协方差;但当各个变量单位不同时(如身高/体重),这时不同变量之间没有可比性,这时就应该消除量纲的影响(即除以样本标准差)。

参考资料

1.https://www.cnblogs.com/Luv-GEM/p/10765574.html
2.https://www.zhihu.com/question/274997106/answer/1055696026
3.https://zhuanlan.zhihu.com/p/454447043
4.https://www.jianshu.com/p/c21c0e2c403a

未经允许不得转载:成千累万网 » 主成分分析法(PCA)及其python实现