如何手动计算得到与R中prcomp()$x一致的主成分得分?
问题原因
结果不一致的核心问题出在手动标准化的步骤:R的矩阵运算默认按列做向量回收,你直接对矩阵减长度等于列数的中心向量、除以长度等于列数的缩放向量时,并不会按预期逐列匹配对应统计量计算,而是会把短向量按列循环补到和矩阵等长后再做运算,得到的标准化矩阵从第二个元素开始就已经不符合中心化缩放的要求,后续乘旋转矩阵自然和pca.result$x的结果不匹配。
修正方案
直接用R内置的scale()函数完成标准化即可,该函数会自动逐列做中心化和缩放,和prcomp()内部的预处理逻辑完全一致,修正后的代码如下:
scale(as.matrix(USArrests), center = pca.result$center, scale = pca.result$scale) %*% pca.result$rotation
如果不想依赖scale()函数,也可以用sweep()显式指定按列(MARGIN=2)执行运算,结果完全等价:
# 逐列减均值 std_mat <- sweep(as.matrix(USArrests), MARGIN = 2, STATS = pca.result$center, FUN = "-") # 逐列除以标准差 std_mat <- sweep(std_mat, MARGIN = 2, STATS = pca.result$scale, FUN = "/") # 乘旋转矩阵得到PC得分 std_mat %*% pca.result$rotation
结果验证
运行以下代码可以确认手动计算结果和prcomp输出完全一致:
pca.result <- prcomp(USArrests, scale=TRUE) manual_x <- scale(as.matrix(USArrests), center = pca.result$center, scale = pca.result$scale) %*% pca.result$rotation # 计算两个结果的最大绝对差,返回0即完全一致 max(abs(manual_x - pca.result$x))
补充说明:PCA的特征向量本身符号是不固定的,同方向的特征向量整体取反也属于合法结果,但由于这里直接使用prcomp输出的固定旋转矩阵计算,不会出现符号差异,结果可以做到完全相等。
内容的提问来源于stack exchange,提问作者Fanta
相关产品推荐
相关产品推荐

