GWASCook RR语言与统计分析

PCA计算原理(超简单,R代码)

2019-08-27  本文已影响0人  佛系分析师

今天介绍如下三点内容:

1. 使用R语言自带函数prcomp计算PCA 并计算每个PC所能解释的方差。

1.1 先随机生成一个矩阵用于测试

#随机生成10*5的矩阵。
set.seed(2019)
mx <- matrix(rnorm(10*5), nrow=10, ncol=5)
mx
捕获.PNG

1.2 prcomp 计算PCA

#R语言自带函数计算PCA
pca <- prcomp(mx)
str(pca)
捕获.PNG
其中pca$x是计算出来的主成分(PC), pca$rotation是特征分解得到的特征向量, pca$sdev是每个主成分的标准差。因此,可以利用这个标准差的结果,我们就可以计算每个PCA所能解释的方差。
我们可以写个代码验证下 pca$sdev是不是每个主成分的标准差
pca$sdev
apply(pca$x,2,sd)
捕获.PNG
可以看到 pca$sdev确实是每个主成分的标准差

1.3 计算每个PCA所能解释的方差比例

#解释方差比例
pcvar <- apply(pca$x,2,var)
pcvar/sum(pcvar)

#利用标准差的结果计算
pca$sdev^2/sum(pca$sdev^2)
捕获.PNG

可以看到,结果是完全一致的。

2. 自己编写R代码计算PCA

2.1 PCA计算原理

只需要三步就可以自己编码实现PCA的计算:

  1. 计算原矩阵X的方差协方差矩阵,记作Cov(X)
  2. Cov(X)进行特征分解
  3. 用原矩阵(至少需要进行中心化)乘以对应的特征向量得到PC。


    捕获.PNG

也许语言不好理解,我们用代码来理解

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)
捕获.PNG
因此,我们可以得出结论:特征值占总特征值的比例,等于每个PC的方差占总方差的比例
上一篇 下一篇

猜你喜欢

热点阅读