降维与度量学习

最后发布时间 : 2023-10-10 23:08:01 浏览量 :

学习资料

使用函数`prcomp`计算主成分

# 加载 iris 数据集
data(iris)

iris数据集总共有5列，其中前四列是特征，最后一列是标签

    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1            5.1         3.5          1.4         0.2     setosa
2            4.9         3.0          1.4         0.2     setosa
3            4.7         3.2          1.3         0.2     setosa
4            4.6         3.1          1.5         0.2     setosa
5            5.0         3.6          1.4         0.2     setosa

提取前四列作为PCA的输入，并为每一个样本生成一个行名

# 提取特征变量
x <- iris[, 1:4]
rownames(x) <- paste("sample",1:nrow(x),sep = "")

        Sepal.Length Sepal.Width Petal.Length Petal.Width
sample1          5.1         3.5          1.4         0.2
sample2          4.9         3.0          1.4         0.2
sample3          4.7         3.2          1.3         0.2
sample4          4.6         3.1          1.5         0.2
sample5          5.0         3.6          1.4         0.2
sample6          5.4         3.9          1.7         0.4

使用函数prcomp执行pca

# 执行主成分分析
pca <- prcomp(x, scale= TRUE,center = TRUE)

打印pca变量，可以看到分析结果包括主成分的标准差和Rotation矩阵，其中Rotation矩阵的行是特征，列是主成分

稍后可以看到

主成分的标准差是：相关系数矩阵特征值的标准差
Rotation矩阵是：相关系数矩阵特征向量

Standard deviations (1, .., p=4):
[1] 1.7083611 0.9560494 0.3830886 0.1439265

Rotation (n x k) = (4 x 4):
                    PC1         PC2        PC3        PC4
Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

提取降维后的主成分矩阵

res1 <- as.data.frame(pca$x)

              PC1        PC2         PC3          PC4
sample1 -2.257141 -0.4784238  0.12727962  0.024087508
sample2 -2.074013  0.6718827  0.23382552  0.102662845
sample3 -2.356335  0.3407664 -0.04405390  0.028282305
sample4 -2.291707  0.5953999 -0.09098530 -0.065735340
sample5 -2.381863 -0.6446757 -0.01568565 -0.035802870
sample6 -2.068701 -1.4842053 -0.02687825  0.006586116

这个主成分矩阵是如何计算得到？

scale(x, center = TRUE, scale = TRUE) %*%  as.matrix(pca$rotation)

可以看到是中心化后的特征变量(150 * 4)乘以rotation矩阵(4 * 4)，也就是说，只要计算得到rotation矩阵就可以得到主成分矩阵

> scale(x, center = TRUE, scale = TRUE)
          Sepal.Length Sepal.Width Petal.Length   Petal.Width
sample1    -0.89767388  1.01560199  -1.33575163 -1.3110521482
sample2    -1.13920048 -0.13153881  -1.33575163 -1.3110521482

> as.matrix(pca$rotation) 
                    PC1         PC2        PC3        PC4
Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

提取主成分的方差贡献率

summary(pca)

其中第一行是相关系数矩阵特征值的标准差
第二行是相关系数矩阵特征值的比例，也就是主成分的方差贡献率，也可以通过第一行计算得到pca$sdev^2 / sum(pca$sdev^2)
第三行是累计比例：0.7296 +0.2285 =0.9581

Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

使用`相关系数`计算`Rotation`矩阵

首先计算相关系数矩阵

cor_mat <- cor(x)

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000

求解相关系数矩阵的特征值和特征向量

rs_mat <- eigen(cor_mat)

可以看到其中特征向量矩阵与Rotation矩阵是相同的

eigen() decomposition
$values
[1] 2.91849782 0.91403047 0.14675688 0.02071484

$vectors
           [,1]        [,2]       [,3]       [,4]
[1,]  0.5210659 -0.37741762  0.7195664  0.2612863
[2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
[3,]  0.5804131 -0.02449161 -0.1421264 -0.8014492
[4,]  0.5648565 -0.06694199 -0.6342727  0.5235971

计算主成分的方差贡献率

对特征值取平方根，即为函数prcomp结果的Standard deviation

sqrt(rs_mat$values)

[1] 1.7083611 0.9560494 0.3830886 0.1439265

计算特征值的比例，即为主成分的方差贡献率

rs_mat$values/sum(rs_mat$values)

[1] 0.729624454 0.228507618 0.036689219 0.005178709

使用奇异值分解(SVD)计算`Rotation`矩阵

# 对特征变量进行中心化
x_centered <- scale(x, center = TRUE, scale = TRUE)
svd_result <- svd(x_centered)

# 提取奇异值、左奇异向量和右奇异向量
singular_values <- svd_result$d
left_singular_vectors <- svd_result$u
right_singular_vectors <- svd_result$v

可以看多SVD的右矩阵与Rotation矩阵相同

> singular_values (1 * 4)
[1] 20.853205 11.670070  4.676192  1.756847

> head(left_singular_vectors) (150*4)
            [,1]        [,2]         [,3]         [,4]
[1,] -0.10823953 -0.04099580  0.027218646  0.013710648
[2,] -0.09945776  0.05757315  0.050003401  0.058435855
[3,] -0.11299630  0.02920003 -0.009420891  0.016098333
[4,] -0.10989710  0.05101939 -0.019457133 -0.037416661
[5,] -0.11422046 -0.05524180 -0.003354363 -0.020379051
[6,] -0.09920300 -0.12718049 -0.005747892  0.003748828

> right_singular_vectors (4*4)
           [,1]        [,2]       [,3]       [,4]
[1,]  0.5210659 -0.37741762  0.7195664  0.2612863
[2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
[3,]  0.5804131 -0.02449161 -0.1421264 -0.8014492
[4,]  0.5648565 -0.06694199 -0.6342727  0.5235971

计算主成分的方差贡献率

singular_values^2/sum(singular_values^2)

[1] 0.729624454 0.228507618 0.036689219 0.005178709

主成分可视化

使用函数biplot

biplot(pca)

使用ggplot

library(ggplot2)  # 加载ggplot2库

summ1 <- summary(pca)  # 对pca进行汇总统计
xlab1 <- paste0("PC1(",round(summ1$importance[2,1]*100,2),"%)")  # 设置x轴标签
ylab1 <- paste0("PC2(",round(summ1$importance[2,2]*100,2),"%)")  # 设置y轴标签

ggplot(data = res1, aes(x = PC1, y = PC2, color = iris$Species)) +  # 创建ggplot对象，设置数据和映射
  stat_ellipse(aes(fill = iris$Species), type = "norm", geom = "polygon", alpha = 0.25, color = NA) +  # 添加置信椭圆
  geom_point(size = 3.5) +  # 添加散点图
  labs(x = xlab1, y = ylab1, color = "Condition", title = "PCA Scores Plot") +  # 设置标签和标题
  guides(fill = "none") +  # 隐藏填充颜色的图例
  theme_bw() +  # 设置主题为白色背景
  scale_fill_manual(values = c("purple", "orange", "pink")) +  # 设置填充颜色的手动比例尺
  scale_colour_manual(values = c("purple", "orange", "pink")) +  # 设置颜色的手动比例尺
  theme(plot.title = element_text(hjust = 0.5, size = 15),
        axis.text = element_text(size = 11),
        axis.title = element_text(size = 13),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 13),
        plot.margin = unit(c(0.4, 0.4, 0.4, 0.4), 'cm'))  # 设置图表的主题样式

PCA Loadings

pca$rotation

因子负荷量是一个矩阵，每一列代表一个主成分，每一行代表一个原始变量。因子负荷量的元素表示了每个原始变量对应主成分的权重或相关性。

                    PC1         PC2        PC3        PC4
Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

聚类特征选择与稀疏学习

降维与度量学习

使用函数prcomp计算主成分

提取降维后的主成分矩阵

提取主成分的方差贡献率

使用相关系数计算Rotation矩阵

计算主成分的方差贡献率

使用奇异值分解(SVD)计算Rotation矩阵

计算主成分的方差贡献率

主成分可视化

PCA Loadings

使用函数`prcomp`计算主成分

使用`相关系数`计算`Rotation`矩阵

使用奇异值分解(SVD)计算`Rotation`矩阵