主成分分析全
主成分分析
一、介绍(introduction)
多元统计处理的是多变量的问题,由于变量个数多,并且彼此之间存在着一定的相关性,因而使得所观测到的数据在一定程度上反映的信息有所重叠,这种信息的重叠有时甚至会掩盖事物的真正特征与内在规律,而且当变量较多时,在高维空间中研究样本的规律比较复杂,增加了分析问题的复杂性,自然希望能够用较少的综合变量来替代原来较多的变量,而这几个综合变量又能够尽可能多地反映原来变量的信息,并且彼此之间互不相关。
主成分分析是利用降维的思想,在损失很少信息的前提下,把多个指标转化为几个综合指标的多元统计方法。通常把转化成的综合指标称为主成分,其中每个主成分都是原始变量的线性组合,且各个主成分之间互不相关,使得主成分比原始变量具有某些更优越的性能。
在主成分分析中,首先对给定给数据进行规范化,使得数据每一变量的平均值为0,方差为1,之后对数据进行正交变换,原来由线性相关变量表示的数据,通过正交变换变成若干个线性无关的新变量表示的数据。新变量是可能的正交变换中变量的方差和最大的,这里方差表示信息量的大小,根据新变量方差的大小将新变量称为第一主成分,第二主成分···等,这就是主成分的基本思想。
二、总体主成分
设$\boldsymbol{x}=(x_1,x_2,···,x_p)^{T}$,为一个$p$维随机变量,并且假定二阶矩存在,记$\boldsymbol{u}=E(x)=(u_1,u_2,···,u_p)$, $\boldsymbol{\Sigma}=cov(\boldsymbol{x,x})=E[(\boldsymbol{x-u})(\boldsymbol{x-u})^T]$,考虑如下线性变换:
$$\enclose{box}{ \begin{cases} y_1=\alpha_{11}x_1+\alpha_{21}x_2+···+\alpha_{p1}x_p=\boldsymbol{\alpha_1^Tx} \\ y_2=\alpha_{12}x_1+\alpha_{22}x_2+···+\alpha_{p2}x_p=\boldsymbol{\alpha_2^Tx} \\ ··· \\ y_p=\alpha_{1p}x_1+\alpha_{2p}x_2+···+\alpha_{pp}x_p=\boldsymbol{\alpha_p^Tx} \end{cases} \tag{1.1} } $$
可以得出:
$$\enclose{box}{ \begin{cases} E(y_i)=\boldsymbol{\alpha_{i}^Tu} \\ var(y_i)=\boldsymbol{\alpha_{i}^T{\Sigma}\alpha_i} \\ cov(y_i,y_j)=\boldsymbol{\alpha_{i}^T\Sigma\alpha_{j}}, i=1,2,···,m; j=1,2,···,m \\ \end{cases} \tag{1.2} } $$
定义1:总体主成分的定义,给定一个如式子(1.1)所示的线性变换,如果满足下列 条件:
-
(1) 系数向量$\boldsymbol{\alpha_i^T}$是单位向量,即$\boldsymbol{\alpha_i^T\alpha}=1,i=1,2,···,m$
-
(2) 新变量$y_i与y_j$互不相关,也就是$cov(y_i,y_j)=0(i{\neq}j)$
-
(3) 变量$y_1$是$\boldsymbol{x}$的所有线性变换中方差最大的,$y_2$是与$y_1$不相关的$\boldsymbol{x}$的所有线性变换中方差最大的,···
这时分别称$y_1,y_2,···,y_p$为$\boldsymbol{x}$的第一主成分、第二主成分,···,第$p$主成分。
$$\bbox[#9ff, 5px]{ (1)给出了一组标准正交基,(2)(3)给出了求主成分的方法。 } $$
主成分的主要性质
1.主成分与协方差矩阵的特征值和特征向量的关系:从主成分思想出发求解主成分
设$\boldsymbol{x}$是$p$维随机变量,$\boldsymbol{\Sigma}$是$\boldsymbol{x}$的协方差矩阵,$\Sigma$的特征值分别是:
$\lambda_1\geq{\lambda_2}\geq{···}\geq{\lambda_p}\geq{0}$,特征值对应的单位特征向量分别是$\boldsymbol{\alpha_1,\alpha_2,···,\alpha_p}$,则$\boldsymbol{x}$的第$k$主成分是:
$\boldsymbol{y_k}=\boldsymbol{\alpha_k^Tx}=\alpha_{1k}x_1+\alpha_{2k}x_2+···+\alpha_{pk}x_p$
第$k$主成分的方差是:
$$var(y_k)=E(\boldsymbol{(\alpha_k^T(x-u))^2}) \\ =E(\boldsymbol{\alpha_k^T(x-u)\alpha_k^T(x-u)}) \\ = E(\boldsymbol{\alpha_k^T(x-u)(x-u)^T\alpha_k)} \\ = \boldsymbol{\alpha_k^T\Sigma\alpha_k}=\lambda_k\\ Cov(y_i,y_j)=\boldsymbol{\alpha_i^T \Sigma\alpha_j}=0,(i\neq j) $$
证明1:$\bbox[#9ff, 5px]{特征值和特征向量作为已知的条件,可以直接使用。}$
设$\boldsymbol{A}$是$p$阶对称阵,它的特征值是$\lambda_1\geq\lambda_2\geq···\geq\lambda_p$,相应的一组正交特征向量是$\boldsymbol{\alpha_1,\alpha_2,···,\alpha_p}$,则有:
$$\max_{x\neq 0}\frac{\boldsymbol{x^TAx}}{\boldsymbol{x^Tx}}=\lambda_1 \quad (当x=\boldsymbol{\alpha_1})\\ \max_{x\neq 0}\frac{\boldsymbol{x^TAx}}{\boldsymbol{x^Tx}}=\lambda_p\quad(当x=\boldsymbol{\alpha_p})\\ $$ $$\max_{x\bot \alpha_1,···,\alpha_k}\frac{\boldsymbol{x^TAx}}{\boldsymbol{x^Tx}}=\lambda_{k+1} \\\quad(当x=\boldsymbol{\alpha_{k+1}} ,(k=1,2,···,p-1时达到) $$
证明过程参见$\bbox[skyblue,5px]{单位球面上点的二次型的极大化(p62)}$
将这里的$\boldsymbol{A}$换成$\boldsymbol{\Sigma}$可得:
$$\max_{x\neq 0}\frac{\boldsymbol{x^T\Sigma x}}{\boldsymbol{x^Tx}}=\lambda_1 \quad (当x=\alpha_1)\\ $$
由于$\boldsymbol{\alpha_1}$是单位特征向量,$\boldsymbol{\alpha_1^T\alpha_1}=1$,所以:
$$\max_{x\neq 0}\frac{\boldsymbol{x^T\Sigma x}}{\boldsymbol{x^Tx}}=\lambda_1 = \frac{\boldsymbol{\alpha^T\Sigma \alpha}}{\boldsymbol{\alpha^T\alpha}}=\boldsymbol{\alpha^T\Sigma\alpha}=var(y_1) $$
由1.4
式可得:
$$\max_{x\bot \alpha_1,···,\alpha_k}\frac{\boldsymbol{x^TAx}}{\boldsymbol{x^Tx}}=\lambda_{k+1}\quad k=1,2,···,p-1 $$
于是取$\boldsymbol{x=\alpha_{k+1}}$使得$\boldsymbol{\alpha_{k+1}^T\alpha_i}=0$,$k=1,2···,p-1$。故有:
$$\frac{\boldsymbol{\alpha_{k+1}^T\Sigma\alpha_{k+1}}}{\boldsymbol{\alpha_{k+1}^T\alpha_{k+1}}}=\boldsymbol{\alpha_{k+1}^T\Sigma\alpha_{k+1}}=var(y_{k+1}) $$
又由于$\boldsymbol{\alpha_{k+1}}$为$\boldsymbol{\Sigma}$的特征值$\lambda_{k+1}$对应的特征向量,所以:
$$\boldsymbol{\alpha_{k+1}^T(\Sigma\alpha_{k+1})}=\lambda_{k+1}\boldsymbol{\alpha_{k+1}^T\alpha_{k+1}}=\lambda_{k+1} $$
故:
$$var(y_{k+1})=\lambda_{k+1} $$
当$\Sigma$的特征值均不相等时,对应的特征值相互正交,而如果出现相同的特征值时,可以在取特征向量时取正交的,也可以利用施密特正交变换进行变换,总之能够保证各特征向量之间两两正交$\boldsymbol{\alpha_i^T\alpha_k}=0\quad i\neq k$。因此:
$$Cov(y_i,y_k)=\boldsymbol{\alpha_i^T\Sigma\alpha_k}=\lambda_k\boldsymbol{\alpha_i^T\alpha_k}=0 $$
证明2:$\bbox[#9ff, 5px]{\boldsymbol{特征值和特征向量是需要证明的内容,不能直接使用。}}$
$\color{blue}$求主成分的方法是拉格朗日乘子法:首先求第一主成分$\boldsymbol{\alpha_1^Tx}$,根据主成分的定义可知$\alpha_1$ 是在$\alpha_1^{T}\alpha_1=1$
的条件下,$\boldsymbol{x}$的所有线性变换中使方差$var(\boldsymbol{\alpha_1^Tx}=\alpha_1^T\Sigma\alpha_1) $达到最大的,由此可以写出一个优化问题:
$$\enclose{box}{ \color{blue} \max_{\alpha_1}\quad\boldsymbol{\alpha_1^T\Sigma\alpha_1}\\ \color{blue}\qquad s.t.\quad\boldsymbol{\alpha_1^T\alpha_1=1}\\ } $$
于是引入参数$\lambda$写出拉格朗日函数:
$$\color{blue} \boldsymbol{\alpha_1^T\Sigma\alpha_1-\lambda(\alpha_1^T\alpha_1-1)}\\ $$
对其求偏导令其为0,得:(这里涉及到对向量的导数,参考了求导公式,如下表)
$$\color{blue} \boldsymbol{(\Sigma+\Sigma^T)\alpha_1-2\lambda\alpha_1=0}\\ $$
由于$\Sigma$是实对称矩阵,因此有$\Sigma^T=\Sigma$,化简得到:
$$\color{blue}\boldsymbol{\Sigma\alpha_1-\lambda\alpha_1=0}\\ $$
根据特征值与特征向量的定义知:$\lambda$是$\Sigma$的特征值,$\alpha_1$是对应的特征向量,则第一主成分对应的方差为:$\lambda_1\\$ 设$\lambda_1$是$\Sigma$的最大特征值,则其对应特征向量为$\alpha_1\\$
$$\color{blue}var(\boldsymbol{\alpha_1^Tx})=\boldsymbol{\alpha_1^T\Sigma\alpha_1=\lambda_1} $$ $$\color{blue}根据向量的求导公式可以得出: \frac{\partial(\alpha_1,\alpha_2,···,\alpha_p)}{\partial(\alpha_1,\alpha_2,···,\alpha_p)}=\begin{pmatrix} \frac{\partial\alpha_1}{\partial\alpha_1} & \frac{\partial\alpha_1}{\partial\alpha_2} & \cdots & \frac{\partial\alpha_1}{\partial\alpha_p} \\ \frac{\partial\alpha_2}{\partial\alpha_1} & \frac{\partial\alpha_2}{\partial\alpha_2} & \cdots &\frac{\partial\alpha_2}{\partial\alpha_p} \\ \vdots & \vdots & \vdots & \vdots \\ \frac{\partial\alpha_p}{\partial\alpha_1}& \frac{\partial\alpha_p}{\partial\alpha_2} & \cdots & \frac{\partial\alpha_p}{\partial\alpha_p} \\ \end{pmatrix}\\ 可以得到:\frac{\partial(\alpha_1,\alpha_2,···,\alpha_p)^T}{\partial(\alpha_1,\alpha_2,···,\alpha_p)}=E(单位阵)\\ 以及 : \frac{\partial(\alpha_1,\alpha_2,···,\alpha_p)}{\partial(\alpha_1,\alpha_2,···,\alpha_p)}=E(单位阵),由此可以得到下表的求导公式: $$
image-20221010212829463
继续求第二主成分:$y_2=\boldsymbol{\alpha_2^Tx}$,第二主成分满足的条件是,在$\boldsymbol{\alpha_2^T\alpha_2=1}$,且 使得$\boldsymbol{\alpha_2^Tx}与\boldsymbol{\alpha_1^Tx}$不相关的$\boldsymbol{x}$的所有线性变换中,使得方差$var(\alpha_2^Tx)=\boldsymbol{\alpha_2^T\Sigma\alpha_2}$最大,可以得到优化问题:
$$\enclose{box}{ \color{blue}\max_{\alpha_2}\quad\alpha_2^T\Sigma\alpha_2\\ \color{blue}s.t. \begin{cases}\quad \left. \begin{array}{l} \alpha_1^T\Sigma\alpha_2=0;\\ \color{blue}\alpha_2^T\Sigma\alpha_1=0;\\ \color{blue}\quad\alpha_2^T\alpha_2=1; \end{array} \right\} \Longleftrightarrow \begin{cases} \alpha_1^T\alpha_2=0;\\ \alpha_2^T\alpha_1=0;\\ \alpha_2^T\alpha_2=1; \end{cases} \end{cases}\\ } $$
解得$\lambda_2$是$\Sigma$的特征值,$\alpha_2$是对应的单位特征向量,可见,第二主成分对应的方差等于第二大特征值:
$$\color{blue}var(\alpha_2^Tx)=\alpha_2^T\Sigma\alpha_2=\lambda_2 $$
根据不同特征值对应的特征向量相互正交,如果有共同特征值时,也可以将特征值选为正交的,所以有:
$$Cov(y_i,y_j)=\alpha_i^T\Sigma\alpha_j=\alpha_i^T\lambda_j\alpha_j =\lambda_j\alpha_i^T\alpha_j=0, (i\neq j) $$
性质1:
$$\enclose{box}{ 主成分Y的协方差阵为对角阵\Lambda } $$
由于$y_i$与$y_j(i\neq j)$,所以相关系数为0,即协方差为0,故其协方差阵为对角阵。
性质2:
$$\enclose{box}{ 记\Sigma=(\sigma_{ij})_{p\times p},有\sum_{i=1}^{p}\lambda_i=\sum_{i=1}^{p}\sigma_{ii} } $$
根据矩阵的迹有:
$$\sum_{i=1}^p\sigma_{ii}=tr(\Sigma) $$
又根据实对称矩阵的对角化:
$$\Sigma=A\Lambda A^T $$
其中$A$是由$\Sigma$的单位特征向量构成的矩阵,$\Lambda$是由$\Sigma$的特征构成的对角阵,特征值和特征向量的相互对应。则:
$$tr(\Sigma)=tr(A\Lambda A^T)=tr(\Lambda A^TA)=tr(\Lambda)=\sum_{i=1}^p\lambda_i=\sum_{i=1}^pvar(y_i) $$
定义2:称
$$\enclose{box}{ r_k=\frac{\lambda_k}{\lambda_1+\lambda_2+···+\lambda_p}(k=1,2,···,p) } $$
为第$k$个主成分$y_k$的方差贡献率。
$$\frac{\sum_{i=1}^{m}\lambda_i}{\sum_{i=1}^{p}\lambda_i} $$
称为主成分$y_1,y_2,···,y_m$的累计贡献率。方差贡献率越大表示所含的信息量越多,对变量的解释能力越强。主成分分析的目的之一是为了降维,减少变量的数量,所以会选择$m(m<p)$个主成分,主成分选择的个数m就是根据累计贡献率的大小选择的,通常所取的m能够使得累计贡献率的大小在85%以上。< p="">
定义3:第$k$个主成分$y_k$与原始变量$x_i$的相关系数$\rho(y_i,x_i)$称为因子载荷量。
因子载荷量是主成分解释的重要依据,因子载荷量的绝对值大小刻画了主成分的主要意义及成因。
性质3:
$$\enclose{box}{ \rho(y_k,x_i)=\frac{\sqrt{\lambda_k}\alpha_{ik}}{\sqrt{\sigma_{ii}}} } $$
证明:不妨设单位向量$e_i^T=(0,···,0,1,0,···,0)$,$\alpha_k^T=(\alpha_{1k},···,\alpha_{ik},··,\alpha_{pk})$,则可得出:
$$x_i=e_i^T\boldsymbol{x}\qquad and\qquad y_i=\alpha_k^T\boldsymbol{x} $$
因此有:
$$\rho(y_k,x_i)&=\frac{cov(y_k,x_i)}{\sqrt{var(y_k)var(x_i)}}\\ &=\frac{cov(\alpha_k^T\boldsymbol{x},e_i^T\boldsymbol{x})}{\sqrt{var(y_k)var(x_i)}}\\ &=\frac{\alpha_k^T\boldsymbol{\Sigma}e_i}{\sqrt{\lambda_k\sigma_{ii}}}\\ &=\frac{e_i^T\Sigma\alpha_k}{\sqrt{\lambda_k\sigma_{ii}}}\\ &=\frac{\lambda_k e_i^T\alpha_k}{\sqrt{\lambda_k\sigma_{ii}}}\\ &=\frac{\lambda_k\alpha_{ki}}{\sqrt{\lambda_k\sigma_{ii}}}\\ &=\frac{\sqrt{\lambda_k}\alpha_{ki}}{\sqrt{\sigma_{ii}}} $$
由该性质可以得出因子载荷量$\rho(y_k,x_i)$与系数$\alpha_{ki}$成正比,与$x_i$的标准差成反比。在解释主成分的成因或第$i$个变量对第$k$个主成分的重要性时,应当根据因子载荷量。
性质4:
$$\enclose{box}{ \sum_{i=1}^{p}\rho^2(y_k,x_i)\sigma_{ii}=\lambda_k } $$
证明:
$$\sum_{i=1}^{p}\rho^2(y_k,x_i)\sigma_{ii}=\sum_{i=1}^{p}\lambda_k\alpha_{ki}^2\\ =\lambda_k\sum_{i=1}^{p}\alpha_{ki}^2\\ =\lambda_k\\ (因为\alpha_k是\boldsymbol{x}的第k个单位特征向量,因此模长为1) $$
性质5:
$$\color{blue} \enclose{box}{ \sum_{i=1}^{p}\rho^2(y_k,x_i)=1 } $$
证明1:
由于任意$y_k,k=1,2,···,p$线性无关,而且$y_k$均由$x_i,i=1,2,···,p$线性表示,所以$x_i$也可由$y_1,y_2,···,y_p$线性表出,易知$x_i=\sum_{j=1}^{p}\alpha_{ji}y_i$,由此可得:
$$\sum_{k=1}^{p}\rho^2(y_k,x_i)=\sum_{i=1}^{p}\frac{cov(y_k,\sum_{j=1}^{p}\alpha_{ji}y_i)^2} {var(y_k)var(\sum_{j=1}^{p}\alpha_{ji}y_i)}\\ =\sum_{k=1}^{p}\frac{\lambda_{k}^2 \alpha_{ki}^2} {\lambda_k \sum_{j=1}^{p}\alpha_{ji}^2\lambda_j}\\ =\frac{\sum_{k=1}^{p}\lambda_{k} \alpha_{ki}^2}{\sum_{j=1}^{p}\alpha_{ji}^2\lambda_j}\\ =1 $$
注意:
$$A= \begin{pmatrix} {\boldsymbol{\alpha_1}}\\ {\boldsymbol{\alpha_2}}\\ \vdots\\ {\boldsymbol{\alpha_p}} \end{pmatrix}= \begin{pmatrix} {\alpha_{11}}&{\alpha_{12}}&\cdots&{\alpha_{1p}}\\ {\alpha_{21}}&{\alpha_{22}}&\cdots&{\alpha_{2p}}\\ \vdots&\vdots&\vdots&\vdots\\ {\alpha_{p1}}&{\alpha_{p2}}&\cdots&{\alpha_{pp}} \end{pmatrix}; \boldsymbol{x}=\begin{pmatrix} {x_1}\\ {x_2}\\ \vdots\\ {x_p} \end{pmatrix}\\ A\boldsymbol{x}=\begin{pmatrix} {y_1}\\ {y_2}\\ \vdots\\ {y_p} \end{pmatrix} =\boldsymbol{y}\\ 由于特征向量为单位特征向量,\\且不同特征值对应的特征向量相互正交,有:\\ AA^T=E,则可以得出\boldsymbol{x}=A^T\boldsymbol{y} $$
证明2:
回顾复相关系数的定义:$y$和$\boldsymbol{x}$的线性函数$\boldsymbol{\alpha_i^Tx}$间的最大相关系数为$y$和$\boldsymbol{x}$之间的复相关系数(多重相关系数),它度量了一个变量$y$和一组变量$x_1,x_2,···,x_p$之间的相关程度,其中$\boldsymbol{\alpha_i^Tx}$表示的是最大限度的提取的关于$y$的信息。
因为$y_k$是$x_1,x_2,···,x_p$的线性组合,所以$x_i$也可以精确的表示为$y_1,y_2,···,y_p$的线性组合,也就是说它们之间存在确定的线性关系,所以复相关系数为1。
总体判定系数:表示$y$的方差能够由$x_1,x_2,···,x_p$联合解释的比例。
$$\sum_{i=1}^{p}\rho^2(y_k,x_i)=\rho^2((y_1,y_2,···,y_p),x_i),根据决定系数的定义:\\ =\frac{SSR}{SST}=\frac{\sum_{i=1}^p(\hat{x_i}-\bar{x})^2}{\sum_{i=1}^p(\hat{x_i}-\bar{x})^2+\sum_{i=1}^p(x_i-\hat{x_i})^2}\\由于x_i可由y_1,y_2,···,y_p线性表出,所以x_i-\hat{x_i}=0\\所以可得\frac{SSR}{SST}=1,也就是\sum_{i=1}^{p}\rho^2(y_k,x_i)=1,因此:\\ x_i的变差可以由y_1,y_2,···,y_p完全解释。 $$
从多元正态随机变量导出主成分
两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。
设$\boldsymbol{X}$的分布是$N_p(\boldsymbol{u,\Sigma})$,在以$u$为中心的椭球上,可得$X的密度$:其轴为$\pm c\sqrt{\lambda_i}\alpha_i$
$$\left(x-u\right)^T\Sigma^{-1}\left(x-u\right)=c^2 $$
$(\lambda_i,\alpha_i)$是$\boldsymbol{\Sigma}$的特征值-标准特征向量对,为了方便研究,令$u=0$,根据谱分解定理得:
$$\underset{(p\times p)}{\Sigma^{-1}}=\sum_{i=1}^p\frac{1}{\lambda_i}\underset{(p\times 1)}{\alpha_i}\underset{(1\times p)}{\alpha_i^T} $$ $$c^2=x^T\sum_{i=1}^p\frac{1}{\lambda_i}\alpha_i\alpha_i^Tx=\sum_{i=1}^p\frac{1}{\lambda_i}(x^T\alpha_i)^2\\ =\sum_{i=1}^p\frac{1}{\lambda_i}(\alpha_i^Tx)^2=\sum_{i=1}^{p}\frac{1}{\lambda_i}\alpha_i^Txx^T\alpha_i $$
1.令$\alpha_i^Tx=y_i$,得到:
$$c^2=\sum_{i=1}^p\frac{1}{\lambda_i}y_i^2 $$
由于实对称矩阵的特征值均是正的,这是以$y_i$为坐标轴的标准椭球的表达式,根据椭球的轴都在特征向量的方向上,知主成分均在特征向量方向上,主成分构成了一组正交基。
2.这里$\alpha_i^Tx$也可以看作是一个坐标变换,也就是将$x$变换到特征向量方向上,平方和表示到中心$u$的欧式距离的平方。
$\bbox[skyblue,5px]{如果不懂坐标和基变换的话,说一下}$
坐标变换****和基变换
在$p$维空间中$V$,线性无关的$p$个向量可以作为$V$的一个基,基一般用一维行向量表示。坐标的是指关于基的坐标,$V$中同一个向量关于不同的基的坐标是不一样的。坐标一般用一维列向量表示。
$$基变换\\ 设\alpha_1,···,\alpha_p \quad and \quad \beta_1,···,\beta_p是V的两个基,\\将\beta_1,···,\beta_p表示成\alpha_1,···,\alpha_p的线性组合\\(\beta_1,···,\beta_p)=(\alpha_1,···,\alpha_p)A,其中\\ A=\begin{pmatrix} \alpha_{11}&\cdots&\alpha_{1p}\\ \alpha_{21}&\cdots&\alpha_{2p}\\ \vdots&\vdots&\vdots\\ \alpha_{p1}&\cdots&\alpha_{pp} \end{pmatrix} 也就是\\ \boldsymbol{\beta=\alpha A}这就是基变换 $$
这里$A$就是$\alpha_1,···,\alpha_p$到$\beta_1,···,\beta_p$的过度矩阵。
$$坐标变换\\ 设\gamma=(\alpha_1,···,\alpha_p) \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_p \end{pmatrix} =(\beta_1,···,\beta_p) \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_p \end{pmatrix}\\ =(\alpha_1,···,\alpha_p)A \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_p \end{pmatrix} $$
由同一组基对应的坐标是一样的,有:
$$\begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_p \end{pmatrix} =A \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_p \end{pmatrix} $$
这就是坐标变换。可以观察到基变换都是将基右乘一个过度矩阵
坐标变换都是左乘一个过度矩阵。
3.其中$x_1\alpha_i$表示向量的内积:
$$x_p\alpha_i=|x_p||\alpha_i|\cos\theta,\quad\theta=<x_p,\alpha_i data-tool="mdnice编辑器"> $$
由于$|\alpha_i|=1$,所以:
$$x_p\alpha_i=|x_p|\cos\theta,\quad\theta=<x_p,\alpha_i> $$
如下图1,是二维的情形,可以直观地得到$x_p\alpha_i$是$x_p$在特征向量$\alpha_i$方向上的投影的模长,我们需要的主成分是投影的向量,需要再乘以标准特征向量。这里投影的平方也就是在特征向量上的投影点到原点(中心点$u$)的欧式距离,再除以特征值,根据主成分的方差与特征值相等,可知这是在对投影在特征向量方向上的数据的方差进行归一化,消除方差(也就是量纲)对数据的影响。综上:$c^2$就是$X$在各特征向量方向上的投影到中心点$u$的欧式距离特征值的和。
内积示意图.drawio
根据样本数据从相关矩阵
出发求解主成分
由于不同维度的数据之间存在不同的量纲,为了避免数据量纲不同对数据分析的影响,我们先对数据进行标准化
$$Z_i=\frac{x_i-\bar{x_i}}{\sqrt{s_{ii}}},i=1,2,···,p\\ \bar{x_i}和s_{ii}分别表示x_i的均值和方差:\\ \bar{x_i}=\frac{1}{n}\sum_{j=1}^nx_{ij},i=1,2,···,p\\ s_{ii}=\frac{1}{n-1}\sum_{j=1}^n(x_{ij}-\bar{x_i})^2,i=1,2,···,p\\ 由此可得:E(Z_i)=0,var(Z_i)=1,所以:\\ cov(Z)=S=\begin{pmatrix} 1&s_{12}\rho_{12}&\cdots&s_{1p}\rho_{1p}\\ s_{12}\rho_{12}&1&\cdots&s_{2p}\rho_{2p}\\ \vdots&\vdots&\vdots&\vdots\\ s_{10}\rho_{1p}&s_{2p}\rho_{2p}&\cdots&1 \end{pmatrix}\\= \begin{pmatrix} 1&\rho_{12}&\cdots&\rho_{1p}\\ \rho_{12}&1&\cdots&\rho_{2p}\\ \vdots&\vdots&\vdots&\vdots\\ \rho_{1p}&\rho_{2p}&\cdots&1 \end{pmatrix}=R $$
即此时样本协方差阵等于相关系数矩阵。因此在提取主成分时,两者是等价的。
有特殊结构的协方差矩阵的主成分
$$\Sigma=\begin{pmatrix} \sigma_{11}&0&\cdots&0\\ 0&\sigma_{22}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\sigma_{pp} \end{pmatrix} $$
可见变量之间的相关系数为0,标准特征向量构成的矩阵为单位阵,因此所得到的主成分也就是原变量,并没有发生改变。
$$\Sigma=\begin{pmatrix} \sigma^2&\rho\sigma^2&\cdots&\rho\sigma^2\\ \rho\sigma^2&\sigma^2&\cdots&\rho\sigma^2\\ \vdots&\vdots&\ddots&\vdots\\ \rho\sigma^2&\rho\sigma^2&\cdots&\sigma^2 \end{pmatrix}\Rightarrow\rho=\begin{pmatrix} 1&\rho&\cdots&\rho\\ \rho&1&\cdots&\rho\\ \vdots&\vdots&\ddots&\vdots\\ \rho&\rho&\cdots&1 \end{pmatrix}\\ =\begin{pmatrix} 1-\rho&0&\cdots&0\\ 0&1-\rho&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&1-\rho \end{pmatrix}+\begin{pmatrix} \rho&\rho&\cdots&\rho\\ \rho&\rho&\cdots&\rho\\ \vdots&\vdots&\ddots&\cdots\\ \rho&\rho&\cdots&\rho \end{pmatrix} $$
可以直接得出最大特征值为$1-\rho+p\rho=1+(p-1)\rho$,其余$(p-1)$个特征值都为$(p-1)\rho$
1.相关矩阵的特征值分解求解主成分(传统方法)
(1)对样本数据进行标准化,得到一个标准化后的数据矩阵用$\boldsymbol{X}$表示。
(2)根据标准化矩阵$\boldsymbol{X}$计算出样本相关系数矩阵$\boldsymbol{R}$。
(3)求样本相关系数矩阵$\boldsymbol{R}$的$k$个特征值和对应的$k$个单位特征向量。,根据方差贡献率确定主成分的个数。
(4)求$k$个样本主成分:将单位特征向量作为系数进行线性变换。
(5)计算$k$个主成分与原始变量的相关系数,以及$k$个主成分对原始变量的贡献率。
(6)计算主成分的值
$\enclose{box}{下面用鸢尾花数据集iris给出一个简单的例子}$
数据集中的变量有花萼长(sepal.Length),花萼宽(Sepal.width),花瓣长(Petal.Length),花瓣宽(Petal.width),以及种类(Species)知乎文章[1]
导入数据并对数据进行标准化
data(iris)
dt <- as.matrix(scale(iris[,1:4]))
head(dt)
image-20221109145756570
计算相关系数
r=cor(dt)
r
image-20221109145927912
计算特征值和特征向量,贡献率以及累计贡献率
rs=eigen(r)
# 提取特征值,也就是方差
val=rs$values
# 变换成标准差
std_val=sqrt(val)
# 计算方差贡献率和累计贡献率
proportion_of_variance=val/sum(val) # 方差贡献率
cumulate_proportion=cumsum(proportion_of_variance) # 累计贡献率
proportion_of_variance
cumulate_proportion
image-20221109150007695
确定主成分的个数k
plot(val/sum(val),type="b",
cex=1,
cex.lab=2,
cex.axis=2,
lty=2,
lwd=2,
xlab = "主成分编号",
ylab="特征值(主成分方差)")
image-20221109150031863
求出主成分以及各主成分对原始变量的贡献率
#提取结果中的特征向量(也称为Loadings,载荷矩阵);
U<-as.matrix(rs$vectors)
#进行矩阵乘法,获得PC score;
PC <-dt %*% U # 利用线性变换计算出主成分
colnames(PC) <- c("PC1","PC2","PC3","PC4")
head(PC) # 主成分
cor(PC,dt) # 各主成分对原始变量的贡献率
image-20221109150122311
绘制主成分的散点图
library(tidyverse)
library(ggplot2)
#将iris数据集的第5列数据合并进来;
df<-data.frame(PC,iris$Species)
#提取主成分的方差贡献率,生成坐标轴标题;
xlab<-paste0("PC1(",round(proportion_of_variance[1]*100,2),"%)")
ylab<-paste0("PC2(",round(proportion_of_variance[2]*100,2),"%)")
#绘制散点图并添加置信椭圆;
df %>%
ggplot(aes(x=PC1,y=PC2,color=iris.Species)) +
stat_ellipse(aes(fill=iris.Species),type ="norm", geom ="polygon",alpha=0.2,color=NA)+
geom_point()+labs(x=xlab,y=ylab,color="")+
guides(fill=F)+
theme_classic()+
theme(panel.background = element_rect(color="black",size=0))
image-20221109150144012
画一下3D散点图
library(scatterplot3d)
color = c(rep('purple',50),rep('orange',50),rep('blue',50))
scatterplot3d(df[,1:3],color=color,
pch = 16,angle=30,
box=T,type="p",
lty.hide=2,lty.grid = 2)
legend("topright",c('Setosa','Versicolor','Virginica'),
fill=c('purple','orange','blue'),box.col=NA)
image-20221109150215989
2.矩阵的奇异值分解求解主成分
对于$m\times n$实矩阵$A$,设其秩为$r$,$0<k<r$,则可以将矩阵进行截断奇异值分解:< p=""> $$\enclose{box}{ \color{blue}A\approx U_{m\times k}\Sigma_kV_{n\times k}^T } $$
其中$\Sigma_k$是$k$阶对角矩阵,取自完全奇异值分解的矩阵$\Sigma$的前$k$个对角元素,。$U_{m\times k},V_{n\times k}$分别取自完全奇异值分解的矩阵$U,V$的前$k$列。
定义$n\times m$矩阵:$X$是经过标准化后的数据矩阵。
$$\color{blue}X^{’}=\frac{1}{\sqrt{n-1}}X^T $$
所以:
$$\color{blue}X^{'T}X^{'}=\frac{1}{n-1}XX^T $$
有定义知道,这就是矩阵$X$的协方差阵。所以求主成分就转化为求$X^{'T}X^{'}$的特征值和单位特征向量,根据奇异值分解:
$$\color{blue}X^{'}=U\Sigma V^T $$
矩阵$V$的前$k$列构成$k$个样本主成分:
$$\color{blue}Y=XV $$
利用iris数据实现下:结果和特征值分解是一样的
data(iris)
dt <- as.matrix(scale(iris[,1:4]))
sigma=svd(dt)$d # 从奇异值分解中取出特征值的平方根
U=svd(dt)$u
V=svd(dt)$v
cumulate_proportion=cumsum(sigma^2)/sum(sigma^2) # 平方后得到特征值,然后计算累计贡献率
cumulate_proportion
image-20221110205002522
#求出主成分
dt_=U%*%diag(sigma)%*%t(V) # 将奇异值分解还原元数据发现是一致的
PC_svd=dt%*%V # 利用奇异值分解求解主成分
colnames(PC_svd)=c("PC1","PC2","PC3","PC4")
head(PC_svd)
cor(PC_svd,dt) # 各主成分在原始变量中的贡献率
image-20221110205113740
奇异值分解求主成分的优势
1.在主成分分析中,主成分就是奇异值分解求得的右奇异值向量,求出奇异值分解也就等于求出了主成分。SVD不仅可以让数据在列向量上压缩,还可以在行向量上压缩(相当于把一些信息重合的样本抛弃了)
2.当数据维度较高时,计算协方差的特征值会特别慢,计算量比较大,而截断SVD可以在损失一定精度的条件按下不用计算协方差阵的特征向量即可快速得到相应的矩阵分解,求得主成分。
主成分分析的应用
1.数据降维:在损失信息较小的情况下,降低数据存储空间,提高机器学习算法的速度。
2.数据可视化:将高维的数据投射到2维或3维,便于观察不同主成分之间数据的差异程度。
3.特征空间转化:可以用来将原空间N个维度投射到另一个特征空间,维度仍是N,但可以帮助分类器从另一个角度切割数据。
4.异常检测:思想就是先压缩数据,然后再乘以特征向量的转置,投射回原来的维度,然后比较原数据与投射回去的数据的分布距离。原理就是特征值大的向量包含了主要信息,特征值比较小的代表了噪声,那么原数据与投射后的数据的距离就是由于噪声引起的,说明数据可能存在异常数据。
主成分的图形展示
library(FactoMineR)
library(factoextra)
library(packrat)
library(patchwork)
library(qqplotr)
data8_4=read.table("E8-4.DAT")
names(data8_4)=c("long","wide","heigh")
model3=PCA(data8_4_log,graph=TRUE)
PCA_scores=as.data.frame(model3$ind$coord)
names(PCA_scores)=c("PCA1","PCA2","PCA3")
p1=ggplot(PCA_scores, aes(sample = PCA1)) +
stat_qq_band() +
stat_qq_line() +
stat_qq_point() +
theme_bw()+
labs(x = "PCA_1_Q", y = "PCA_1")
p2=ggplot(PCA_scores, aes(sample = PCA2)) +
stat_qq_band() +
stat_qq_line()+
stat_qq_point()+
annotate("text", x = 0.5, y = 0.8, label = "outler", size = 4)+
theme_bw(base_size = 10)+
labs(x = "PCA_2_Q", y = "PCA_2")
p3=ggplot(data = PCA_scores,aes(x=PCA2,y=PCA1))+
geom_point()+
theme_bw()+annotate("text",x=0.75,y=-2.8,label="outler")
p1+p2+p3
其实总的来看,书上的那个离群点算不上太偏,可以接受。
大样本推断
的大样本的性质
样本协方差阵的特征值$[\hat{\lambda_1,···,\hat{\lambda_p}}]$和特征向量$\boldsymbol{\hat{\alpha_i},···,\hat{\alpha_p}}$的大样本理论:
1.$\enclose{box}{令\boldsymbol{\Lambda}是\boldsymbol{\Sigma}的特征值\lambda_1,···,\lambda_p的对角矩阵,则\sqrt{n}({\hat{\lambda}-\lambda})近似服从分布N_p(\boldsymbol{0},2\boldsymbol{\Lambda}^2).}$
证明:
设总体$\boldsymbol{X}\backsim N_p(\boldsymbol{u},\boldsymbol{\Sigma})$,$x_1,···,x_n$是它的样本,$\boldsymbol{\Sigma}>0,n>p$,由统计学的知识知道$\boldsymbol{\Sigma}$的极大似然估计为:
$$\boldsymbol{\hat{\Sigma}}=\boldsymbol{S}=\frac{\boldsymbol{\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x})^T}}{n}\qquad 是样本协方差阵 $$
设$\boldsymbol{S}$的特征根分别为$\hat{\lambda_1}\geq\cdots\geq\hat{\lambda_p}\geq0$,对应的标准正交特征向量分别为$\boldsymbol{\hat{\alpha_1},\cdots,\hat{\alpha_p}}$,根据极大似然估计的不变性,知$\hat{\lambda_1},\cdots,\hat{\lambda_p}$和$\boldsymbol{\hat{\alpha_i},\cdots,\hat{\alpha_p}}$分别是$\lambda_i,\cdots,\lambda_p$和$\boldsymbol{\alpha_i,\cdots,\alpha_p}$的极大似然估计。根据极大似然估计的渐进正态性知道,它们有渐进正态分布。
样本$x_1,\cdots,x_n$的联合密度为:
$$\frac{1}{(2\pi)^{\frac{np}{2}|\boldsymbol{\Sigma}|^{\frac{n}{2}}}}exp\lbrace -\frac{1}{2}tr(\boldsymbol{\Sigma^{-1}}(\boldsymbol{V}+n(\boldsymbol{\bar{x}-u})(\boldsymbol{\bar{x}-u})^T))\rbrace $$
其中,$\boldsymbol{V}=\sum_{i=1}^p(\boldsymbol{x_i-\bar{x}})(\boldsymbol{x_i-\bar{x}})^T\backsim W_p(n-1,\boldsymbol{\Sigma})$,对于$(\hat{\lambda_1},\cdots,\hat{\lambda_p})$的渐进正态分布的计算来说,其关键就是计算样本$x_1,\cdots,\_n$的Fisher信息量。$\boldsymbol{u},\lambda_1,\cdots,\lambda_p与\boldsymbol{\alpha_1,\cdots,\alpha_p}$的似然函数为:
$$L(\lambda_1,···,\lambda_p,\alpha_1,···,\alpha_p)\frac{1}{|\boldsymbol{\Sigma}|^{\frac{n}{2}}}exp\lbrace -\frac{1}{2}tr(\boldsymbol{\Sigma^{-1}}\boldsymbol{V})\rbrace\tag{*} $$
$\boldsymbol{T=(\alpha_1,···,\alpha_p)}$是正交矩阵且$\boldsymbol{T^T\Sigma T=\Lambda,\Lambda=diag(\lambda_1,···,\lambda_p)}$,从而由$(*)式可知$:
$$L(\lambda_1,···,\lambda_p,\alpha_1,···,\alpha_p)\\ =(\lambda_1···\lambda_p)^{-\frac{n}{2}}exp\lbrace{-\frac{1}{2}(\frac{\boldsymbol{\alpha_1^TV\alpha_1}}{\lambda_1}+···+\frac{\boldsymbol{\alpha_p^TV\alpha_p}}{\lambda_p})\rbrace} $$
由于协方差阵$\boldsymbol{\Sigma}$的特征根$\lambda_1>···>\lambda_p>0$与这些特征根对应的正交特征向量$\boldsymbol{\alpha_1,···,\alpha_p}$没有关系,所以可以得出:对任意的$i \neq j$,都有:
$$\frac{\partial^2 \ln{L(\lambda_1,···,\lambda_p,\alpha_1,···,\alpha_p)}}{\partial\lambda_i\partial\lambda_j}=0 $$
$\boldsymbol{\hat{\alpha_1},···,\hat{\alpha_p}}$并不渐进相互独立。
得到$\lambda_i$的似然函数为
$$L(\lambda_i)=(\lambda_i)^{-\frac{n}{2}}exp\lbrace{-\frac{1}{2}\frac{\boldsymbol{\alpha_i^TV\alpha_i}}{\lambda_i}\rbrace}\qquad i=1,···,p \tag{**} $$
由于$\boldsymbol{\alpha_i^T\Sigma\alpha_i}=\lambda_i,\boldsymbol{V\backsim W_p(n-1,\Sigma)}$,所以$\boldsymbol{\alpha_i^TV\alpha_i}\backsim\lambda_i\chi^2(n-1)$,由$(**)$式得$\hat{\lambda_i}$的Fisher信息为
$$-E(\frac{\partial^2\ln{L(\lambda_i)}}{\partial\lambda_i^2})=\frac{n-2}{2\lambda_i^2} $$
考虑到$\hat{\lambda-1},···,\hat{\lambda_p}$渐进相互独立,从而有$(\hat{\lambda_1,···,\lambda_p})$的渐进正态分布为
$$\require{AMDcd} \begin{CD} \rm{\sqrt{n-2} \begin{pmatrix} \hat{\lambda_i}-\lambda_1\\ \vdots\\ \hat{\lambda_p}-\lambda_p \end{pmatrix}} @>{\large\text{L}}>> \rm{N_p(\boldsymbol{0},diag(2\lambda_1^2,···,2\lambda_p^2))} \end{CD}\tag{\#} $$
在$\boldsymbol{\Sigma}$的特征根全不相等的时候,计算了$(\hat{\lambda_1},···,\hat{\lambda_p})$的渐进正态分布,至于在$\boldsymbol{\Sigma}$的特征根有重根时,关于$(\hat{\lambda_1,···,\hat{\lambda_p}})$的渐进正态分布的计算可以参考。。。
由$(\#)$式得$\lambda_i$的$1-\beta$的渐进置信区间为
$$\frac{\hat{\lambda_j}}{1+\sqrt{\frac{2}{n-2}}U_{1-\frac{\beta}{2}}}\leq\lambda_i\leq\frac{\lambda_j}{1-\sqrt{\frac{2}{n-2}}U_{1-\frac{\beta}{2}}} $$
由此可见,样本容量$n$必须足够大,以使得$\sqrt{\frac{2}{n-2}}U_{1-\frac{\beta}{2}}<1$。在大样本条件下$n$和$n-2$是等价的。
可以得到$\lambda_i$的一个大样本置信区间:
$$\enclose{box}{ \frac{\hat{\lambda_i}}{(1+z(\frac{\alpha}{2})\sqrt{\frac{2}{n}})}\leq\lambda_i\leq\ \frac{\hat{\lambda_i}}{(1-z(\frac{\alpha}{2})\sqrt{\frac{2}{n}})}} $$
其中$z(\frac{\alpha}{2})$是标准正态分布的上$100(\alpha/2)$百分位数
Anderson T.W.An Introduction to Multivariate Statistical Analysis(3rd).New York
2.令
$$\enclose{box}{\boldsymbol{E}_i=\lambda_i\sum_{k=1(k\neq i)}^p\frac{\lambda_k}{(\lambda_k-\lambda_i)}\boldsymbol{\alpha_k\alpha_k^T}} $$
$\enclose{box}{则\sqrt{n}(\boldsymbol{\hat{\alpha_i}-\alpha_i})近似服从分布N_p(\boldsymbol{0,E_i})。}$
3.$\enclose{box}{各个\hat{\lambda_i}与对应的\boldsymbol{\hat{\alpha_i}}的诸元独立分布。}$
用主成分监控质量(其实就是绘制主成分的置信椭圆)
检查一组已知测量值的稳定性
这里数据标准化和不标准化的结果是不一样的
library(FactoMineR)
library(factoextra)
library(patchwork)
data5_8=read.table("T5-8.DAT")
names(data5_8)=c("x1","x2","x3","x4","x5")
model4=PCA(data5_8,scale.unit=FALSE,graph=TRUE)
model5=PCA(data5_8,scale.unit=TRUE,graph=TRUE)
summary(model4)
summary(model5)
p3=fviz_eig(model4,addlabels=TRUE,ylim=c(0,100),main="数据不标准化")
p4=fviz_eig(model5,addlabels=TRUE,ylim=c(0,100),main="数据标准化")
p3+p4
image-20221027182045338
image-20221027182335565
p3=fviz_eig(model4,addlabels=TRUE,ylim=c(0,100),main="数据不标准化")
p4=fviz_eig(model5,addlabels=TRUE,ylim=c(0,100),main="数据标准化")
p3+p4
image-20221027182958709
可见数据标准化和不标准化差距挺大,所以说,在数据的量纲一样的时候,建议不对数据做标准化,有些主成分分析的包是默认对数据做标准化的。
PC4=as.data.frame(model4$ind$coord)
names(PC4)=c("PCA1","PCA2","PCA3","PCA4","PCA5")
PC5=as.data.frame(model5$ind$coord)
names(PC5)=c("PCA1","PCA2","PCA3","PCA4","PCA5")
p5=ggplot(data = PC4,aes(x=PCA1,y= PCA2))+
stat_ellipse(type="norm",geom="polygon",alpha=0.2,color="#4169E1")+
geom_point()+
#labs(x=xlab,y=ylab,color="")+guides(fill=F)+
theme_bw()+ggtitle("不标准化")
p6=ggplot(data = PC5,aes(x=PCA1,y= PCA2))+
stat_ellipse(type="norm",geom="polygon",alpha=0.2,color="#4169E1")+
geom_point()+
#labs(x=xlab,y=ylab,color="")+guides(fill=F)+
theme_bw()+
ggtitle("标准化后")
p5+p6
T方图绘制
PC4_2=PC4*PC4
PC5_2=PC5*PC5
eng_value4=model4$eig[,1]
eng_value5=model5$eig[,1]
PC4_2$SUM=PC4_2$PCA3/eng_value4[3]+PC4_2$PCA4/eng_value4[4]+PC4_2$PCA5/eng_value4[5]
PC5_2$SUM=PC5_2$PCA3/eng_value5[3]+PC5_2$PCA4/eng_value5[4]+PC5_2$PCA5/eng_value5[5]
p7=ggplot(PC4_2, aes(x=c(1:16), y=SUM))+
geom_point()+
geom_line(aes(group=""), color="green")+
labs(x="Time",y="T^2",title="T方图绘制(数据未标准化)")+
theme(plot.title = element_text(hjust = 0.5)) +
theme_bw()+
geom_hline(aes(yintercept=8.75))+
geom_hline(aes(yintercept=0.1),linetype="dashed")
p8=ggplot(PC5_2, aes(x=c(1:16), y=SUM))+
geom_point()+
geom_line(aes(group=""), color="green")+
labs(x="Time",y="T^2",title="T方图绘制(数据标准化了)")+
theme(plot.title = element_text(hjust = 0.5)) +
theme_bw()+
geom_hline(aes(yintercept=8.75))+
geom_hline(aes(yintercept=0.1),linetype="dashed")
p7+p8
image-20221027193058510
控制未来值
去除异常值之后:
data2=data5_8[-11,]
model6=PCA(data2,scale.unit=FALSE,graph=TRUE)
model7=PCA(data2,scale.unit=TRUE,graph=TRUE)
print("下面是不标准化原始数据的结果")
model4$var$cor
print("#########################################")
print("下面是去除异常值不标准化的结果")
model6$var$cor
print("##########################################")
print("下面是原始数据标准化的结果")
model5$var$cor
print("#########################################")
print("下面是去除异常值标准化的结果")
model7$var$cor
image-20221110193814224
可以明显地看到,去除异常值后,主成分与原始数据的相关系数有明显的变化,也就是各个变量对各个主成分的贡献度发生了明显的变化,所以说,在进行主成分分析之前,处理异常值是有必要的。
p5=fviz_eig(model6,addlabels=TRUE,ylim=c(0,100),main="去除异常后数据不标准化")
p6=fviz_eig(model7,addlabels=TRUE,ylim=c(0,100),main="去除异常数据标准化")
p5+p6
image-20221110194043592
p3+p4+p5+p6
image-20221110194939405
p7=ggplot(data = PC3,aes(x=PCA1,y= PCA2))+
stat_ellipse(type="norm",geom="polygon",levels=0.99,alpha=0.01,color="#4169E1")+
geom_point()+
theme_bw()+
ggtitle("去除异常不标准化")
p8=ggplot(data = PC4,aes(x=PCA1,y= PCA2))+
stat_ellipse(type="norm",geom="polygon",levels=0.99,alpha=0.01,color="#4169E1")+
geom_point()+
theme_bw()+
ggtitle("去除异常标准化")
p7+p8
image-20221110194104258
主成分的通径分析与决策分析
第$i$个主成分为:
$$\enclose{box}{ \color{blue}y_i=\alpha_{i1}x_1+\alpha_{i2}x_2+···+\alpha_{ip}x_p=\boldsymbol{\alpha_i^TX} } $$
其中$\boldsymbol{X}\backsim N_m(0,\Sigma),\boldsymbol{\alpha_1}=(\alpha_{i1},···,\alpha_{ip})^T$为$\boldsymbol{\Sigma}$的第$i$个特征值$\lambda_i$所对应的单位特征向量,$var(y_i)=\lambda_i$,下面研究$y_i$中各$x_i$对$y_i$的综合决定能力:
$y_i=\boldsymbol{\alpha_i^TX}$满足$\boldsymbol{\Sigma\alpha_i}=\lambda_i\boldsymbol{\alpha_i}$,展开来写就是:
$$\color{blue} \begin{pmatrix} \sigma_{11}^2&\sigma_{12}^2&\cdots&\sigma_{1p}^2\\ \sigma^2{21}&\sigma_{22}^2&\cdots&\sigma_{2p}^2\\ \vdots&\vdots&\ddots&\vdots\\ \sigma_{p1}^2&\sigma_{p2}^2&\cdots&\sigma_{pp}^2 \end{pmatrix} \begin{pmatrix} \alpha_{i1}\\ \alpha_{i2}\\ \vdots\\ \alpha_{ip} \end{pmatrix} =\begin{pmatrix} \lambda_i\alpha_{i1}\\ \lambda_1\alpha_{i2}\\ \vdots\\ \lambda_i\alpha_{ip} \end{pmatrix} \tag{*} $$
第$i$个方程为:
$$\sigma_{i1}\alpha_{j1}+\sigma_{i2}\alpha_{j2}+···+\sigma_{ip}\alpha_{jp}=\lambda_j\alpha_{ji} $$
两边同除以$\sigma_k\sqrt{\lambda_i}$,并令$b_{ik}=\sigma_k\alpha_{ik}/\sqrt{\lambda_i}$得:
$$\rho_{i1} b_{j1}+\rho_{i2}b_{j2}+···+\rho_{ip}b_{jp}=\frac{\sqrt{\lambda_j}\alpha_{ji}}{\sigma_i} $$
这里利用了相关系数的定义式.
由于主成分$Y=(y_1,···,y_p)=(\alpha_1^TX,···,\alpha_p X)=A^TX$
得到$X=AY,X_i=\alpha_{1i}y_1+\alpha_{2i}y_2+···+\alpha_{pi}y_p$因此有:
$$cov(x_i,y_j)=cov(\alpha_{1i}y_1+\alpha_{2i}y_2+···+\alpha_{pi}y_p,y_j)\\ =\alpha_{ji}\lambda_j\\ \rho(x_i,y_j)=\frac{cov(x_i,y_j)}{\sigma_i\sqrt{\lambda_j}}=\frac{\sqrt{\lambda_j}\alpha_{ji}}{\sigma_i} $$
于是式子(*)可以化为:
$$\color{blue} \begin{pmatrix} 1&\rho_{12}&\cdots&\rho_{1p}\\ \rho_{21}&1&\cdots&\rho_{2p}\\ \vdots&\vdots&\ddots&\vdots\\ \rho_{p1}&\rho_{p2}&\cdots&\rho_{pp} \end{pmatrix} \begin{pmatrix} b_{j1}\\ b_{j2}\\ \vdots\\ b_{jp} \end{pmatrix} =\begin{pmatrix} \rho_{x_1}y_j\\ \rho_{x_2}y_j\\ \vdots\\ \rho_{x_p}y_j \end{pmatrix} \tag{*2} $$
式子(*2)是以$y_j$为因变量的关于$x_1,x_2,···,x_p$的标准化多元线性回归方程的最小二乘方程组,由它1可以进行通径分析,分析可以分为两个部分:
-
对每一个$x_i$,它对$y_j$的总影响为$\rho_{x_i}y_j$,其中$x_i$的直接影响为$b_{ji}$;$x_i$通过$x_k$对$y_j$的间接影响为$\rho_{jk}b_{jk},k\neq i$
-
$x_i$对$y_j$的直接决定系数$R_i^2=b_{ji}^2$,相关路$x_i\leftrightarrow x_k$对$y_j$的相关决策系数为$R_{ik}=2b_{ji}\rho_{ik}b_{jk}$,显然各$x_i$对$y_j$的决定是完全的,因为:
- $$R_{xy_j}^2=\sum_{i=1}^p b_{ji}^2+2\sum_{k\neq i}b_{ji}\rho_{ik}b_{jk}=\sum_{i=1}^p b_{ji}\rho_{x_i}y_j\\ =\sum_{i=1}^p\frac{\sigma_i\alpha_{ji}}{\sqrt{\lambda_j}}(\frac{\sqrt{\lambda_j}\alpha_{ji}}{\sigma_i})=\sum_{i=1}^p\alpha_{ji}^2=1 $$
根据袁志发等(2001)提出的通径分析中的决策系数概念,各$x_i$对主成分$y_j$的决策系数为:
$$R(i)=R_i^2+\sum_{k\neq i}R_{ik}=2b_{ji}\rho_{x_i}y_j-b_{ji}^2\\ =\frac{2\sigma_i\alpha_{ji}}{\sqrt{\lambda_j}}·\frac{\sqrt{\lambda_j}\alpha_{ji}}{\sigma_i}-\frac{\sigma_i^2\alpha_{ji}^2}{\lambda_j}=\alpha_{ji}^2(2-\frac{\sigma_i^2}{\lambda_j}) $$
$R(i)$表示$x_i$对$y_j$的综合决定能力,根据$R(i)$的排序可以确定各$x_i$对$y_j$的综合决定大小和方向。
参考文献
【1】高慧璇.应用多元统计分析[M].北京:北京大学出版社,2005:265
【2】何晓群.多元统计分析(第五版)[M].北京:中国人民大学出版社,2019:106
【3】李航.统计学习方法(第二版)[M].北京:清华大学出版社,2019:297
【4】袁志发.多元统计分析(第二版)[M].北京:科学出版社,2009.
【5】王学民.应用多元统计分析(第五版) [M].上海:上海财经大学出版社,2017
参考资料
[1]
知乎: https://zhuanlan.zhihu.com/p/354086571
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)