今天介绍如下三点内容:
- PCA是如何计算的
- 如何计算每个PC所能解释的方差
- 特征值与PC解释方差比例的关系
1. 使用R语言自带函数prcomp
计算PCA 并计算每个PC所能解释的方差。
1.1 先随机生成一个矩阵用于测试
#随机生成10*5的矩阵。
set.seed(2019)
mx <- matrix(rnorm(10*5), nrow=10, ncol=5)
mx
1.2 prcomp
计算PCA
#R语言自带函数计算PCA
pca <- prcomp(mx)
str(pca)
其中
pca$x
是计算出来的主成分(PC), pca$rotation
是特征分解得到的特征向量, pca$sdev
是每个主成分的标准差。因此,可以利用这个标准差的结果,我们就可以计算每个PCA所能解释的方差。我们可以写个代码验证下
pca$sdev
是不是每个主成分的标准差
pca$sdev
apply(pca$x,2,sd)
可以看到
pca$sdev
确实是每个主成分的标准差
1.3 计算每个PCA所能解释的方差比例
#解释方差比例
pcvar <- apply(pca$x,2,var)
pcvar/sum(pcvar)
#利用标准差的结果计算
pca$sdev^2/sum(pca$sdev^2)
可以看到,结果是完全一致的。
2. 自己编写R代码计算PCA
2.1 PCA计算原理
只需要三步就可以自己编码实现PCA的计算:
- 计算原矩阵
X
的方差协方差矩阵,记作Cov(X)
- 对
Cov(X)
进行特征分解 -
用原矩阵(至少需要进行中心化)乘以对应的特征向量得到PC。
也许语言不好理解,我们用代码来理解
2.2 自编代码计算PCA
#中心化
mx_n <- scale(mx, scale=FALSE)
#计算协方差矩阵,并进行特征分解
e_mx <- eigen(crossprod(mx_n))
#计算PC
pc <- mx_n %*% e_mx$vectors
# 看下和 R语言自带函数的结果是否一致
pc
pca$x
ok,完全一致。
2.3 自编代码计算每个PC解释的方差比例
特征值本身就反应了每个PC解释方差的比例关系,因此利用特征即可计算出每个PC解释方差的比例
e_mx$values/sum(e_mx$values)
# 和之前的比较一下
pcvar/sum(pcvar)
因此,我们可以得出结论:特征值占总特征值的比例,等于每个PC的方差占总方差的比例