可视化：RStudio 多元数据分析 - 主成分分析

Alkohol
568 阅读5分钟

本文已参加「新人创作礼」活动，一起开启掘金创作之路。

这次带来的是 RStudio 的多元数据分析 - 主成分分析法。

各个知识点后面都有对应的小练习哦，大家可以利用刚刚学到的知识来试着写写看！

主成分分析 Principle Component Analysis

主成分分析(PCA: Principle Component Analysis)是一种可将多个变量 / 指标 / 特征转化为少数几个新变量 / 指标 / 特征的方法. PCA由统计学家Pearson在1901年提出,后被Hotelling等人发展起来.

现在PCA已经是一种非常重要的降维方法, 它能帮助我们在处理高维度数据时寻找变量之间的关系和规律.

PCA的原始想法

高维(数值型的,定量的)数据在空间中常常会呈现出一种类似于椭球的分布(为什么? )

image.png

  • 数据椭球的长轴方向代表了数据在这个方向上分布较广,差异较大
  • 数据椭球的短轴方向代表了数据在这个方向上分布较窄,差异较小

由于数据椭球的各个轴通常不会正好沿着坐标轴方向, 那么如果沿着椭球的轴向来设置新的坐标系, 会不会更加便于我们对其进行分析, 进而便于挖掘数据潜在的模式和规律?

所以我们要去寻找椭球的各个轴向, 然后建立新坐标系?

其实重新回顾上面的想法之后, 我们发现并不需要数据是高维空间中的一个椭球这个假设, 其实我们只需要知道数据在某些方向上分布更广/差异更大,而在另一些方向上分布更窄/差异更小.

所以我们要去寻找数据在各个方向上的分布差异.

改进后的想法

假设数据是Rp空间中的一些点. 变量为 x1,x2,...,xp

  • 寻找高维空间中的一个方向 y1=(a1,a2,...,ap)T, 使得数据在这个方向上的差异最大, 也就是数据在 y1 方向上的投影的方差最大
  • 如果我们对数据在 y1 方向上携带的信息已经足够满意(椭球很扁长), 那么找到 y1 就够了
  • 我们可以用数据在 y1 方向上的投影来代替原始的p维数据. 这样就将p维数据变成了 1 维数据.
  • 如果还不满意(在垂直 y1 的空间中还有很多差异)呢?
  • 那就需要在垂直 y1 的空间中继续寻找方差最大的方向 y2=(a1,a2,...,ap)T, 这样得到的 y2 跟 y1 一定是正交的, 数据在这两个方向上的投影不会互相影响.
  • 现在我们有两个方向 y1,y2, 数据在 y1 方向上的投影方差是最大的, 剩余的方差在 y2 方向上的投影方差是最大的.
  • 如果我们对数据在 y1,y2 方向上携带的信息已经足够满意, 那么就可以用数据在 y1,y2 方向上的投影来代替原始的 p 维数据. 这样就将 p 维数据变成了 2 维数据.
  • 如果还不满意, 那么可以重复这个过程, 直到找到所有的 p 个新方向: y1,y2,...yp

数学化的正式定义

image.png

  • 求出一个 a 使得 y=aTx 的方差 Var(aTx)=aTΣa 最大. 此时的 y 称为数据的第一主成分(The First Principle Component)
  • 这里 Σ=Cov(x) 为 x 的协方差阵
  • 在跟第一主成分的方向 a 所正交的空间中, 继续寻找方差最大的投影方向, 得到第二主成分

计算方法

根据线性代数的理论, 我们可以得到如下计算各个主成分的方向 a 方法:

  • 求出 x 的协方差阵 Σ=Cov(x) 的特征值(从大到小排列) λ1, λ2,..., λp
  • 以及对应的特征向量 u1, u2,... up
  • 这些特征向量就是各个主成分的方向 a, 它们之间两两正交
  • 这些特征值就代表了原始数据在每一个主成分方向上的投影方差
# 例:
# 计算iris数据的前四列(数值型变量)的各个主成分以及所携带的方差

# data matrix
X = as.matrix(iris[,1:4])
eig = eigen(cov(X))

# show eigen values and eigen vectors
eig

# get PCs (principle components), or, the matrix of loadings
evectors = eig$vectors

# the PC1 (the first principle component)
evectors[,1]

# the PC2 (the second principle component)
evectors[,1]

# get variance of each PC
evalues = eig$values

# the proportion of variance of each PC
evalues / sum(evalues)
barplot(evalues / sum(evalues), 
        names.arg = paste("PC",1:4),
        main = "the proportion of variance of each PC",
        xlab = "Principle Components", ylab = "Proportion of Variance")

# 使用R包{stats}中的princomp函数

prc = princomp(X, scores = T)

# print the principle components information
prc
summary(prc)

# get standard deviation of each PC
prc$sdev

# get variance of each PC
prc$sdev^2

# get PCs (principle components)
prc$loadings
unclass(prc$loadings)

sum(prc$loadings[,1]^2)

# get the scores
# score is the coordinates of a observation under the new axes system

# the score on PC1
prc$scores[,1]

# the score on PC2
prc$scores[,2]

# the scores on all 4 PCs
prc$scores

# 原始变量如何用主成分表示?
# 原始变量方向的score?
# biplot
biplot(prc)

# the inverse of loadings matrix
solve(unclass(prc$loadings)) -> xs
xs # each column is a variable in terms of a combination of PCs

主成分分析的应用举例 An Application of Principle Component Analysis

当数据的变量/维数太多, 模型可能会出现过度拟合(overfitting); 当各个变量之间存在多重共线性时, 拟合结果可能会不稳定.

主成分分析有助于我们解决以上这些困难. 具体做法是:

  • 将原始变量做主成分分析
  • 用部分(或者全部的)主成分作为新变量, 替换掉原始变量
  • 基于这些主成分进行建模

数据集下载：MNIST数据集

# 回顾一个老问题: 手写数字识别

# 从数据中抽取一个子集, 将原问题(10个手写数字的识别)转化为一个二元分类问题
# 使用已学内容建立分类模型
# 评估模型的准确率

# 从数据中抽取一个子集, 将原问题(10个手写数字的识别)转化为一个二元分类问题
# 对784个维度做主成分分析并降维, 思考降到几维比较合适?
# 使用降维后的数据建立分类模型
# 评估模型的准确率
# 比较两种方法的准确率差异

# 读取MNIST数据集
train = read.csv("./data/mnist_train.csv", header=F)
dim(train)
train[1,]
Ytrain = train[,1]
table(Ytrain)
Xtrain = as.matrix(train[,-1])

# 从数据中抽取一个子集, 将原问题(10个手写数字的识别)转化为一个二元分类问题
# 作为一个演示, 为了减少计算量, 将子集限定在2000行, 其中1000行是数字0, 其余1000行是其他数字
index_0 = sample(which(Ytrain == 0), 1000)
index_other = sample(which(Ytrain != 0), 1000)
train_small = train[c(index_0, index_other),]
dim(train_small)
train_small$V1 = factor(ifelse(train_small$V1 == 0, "Digit_Zero", "Other"))


# 分割数据集
train_idx = sample(2000, 1300)
dat0_train = train_small[train_idx,]
dat0_test = train_small[-train_idx,]

# 建立模型
mod1 = glm(V1 ~ ., data = dat0_train, family = "binomial")
pred1 = round( predict(mod1, newdata = dat0_test, type = "response") )
table(pred1)
pred1 = ifelse(pred1 == 0, "Digit_Zero", "Other")
cm = table(real = dat0_test$V1, pred = pred1)
cm
sum(diag(cm)) / sum(cm) # accuracy

# ==== 利用主成分分析 降维 ====
pc = princomp(train_small[,-1]) # 注意要对dat0_train和dat0_test一起进行主成分分析
cum_pop_var = cumsum(pc$sdev^2) / sum(pc$sdev^2)
plot(cum_pop_var, type = "l")
abline(h = 0.9, col="red", lty=2)

# 要选择前多少个主成分?
which(cum_pop_var > 0.9)[1]

# 得到原始数据样本在前N个主成分上的得分, 并合并V1成为新数据集pc_train,pc_test
N = 74
pc_train = data.frame(V1=dat0_train$V1, pc$scores[train_idx,1:N])
pc_test = data.frame(V1=dat0_test$V1, pc$scores[-train_idx,1:N])

# 重新建模
mod2 = glm(V1 ~ ., data = pc_train, family = "binomial")
pred2 = round( predict(mod2, newdata = pc_test, type = "response") )
table(pred2)
pred2 = ifelse(pred2 == 0, "Digit_Zero", "Other")
cm2 = table(real = pc_test$V1, pred = pred2)
cm2
sum(diag(cm2)) / sum(cm2) # accuracy
avatar
文章
阅读
粉丝