主成分分析（PCA）实用指南（R＆Python）-论智
来源：Analytics Vidhya编译：Bot
编者按：Too much of anything is good for nothing! 今天我们聊聊基础知识，主要介绍一种化繁为简的神奇方法——主成分分析。注意，理解本文需要具备一定统计学基础。

有时候，华婷婷如果数据集里有太多变量，你会怎么办？以下是每一位分析师都可能遇到的几种情况：
你会发现大部分变量都是相关的；
你会失去耐心，决定在所有数据上运行模型，然后准确率果然很低，你很绝望；
你开始怀疑人生：我的目的究竟是什么；
你开始思考如何找出其中的重要变量。
试想一下，当自己身处这四种境地神箓，你会怎么做？不要逃，处理这些情况并没有想象中的那么困难，一些统计方法，如因子分析、主成分分析就有助于克服它们。什么是主成分分析？
简而言之，主成分分析就是从数据集大量可用变量中以成分形式提取出重要变量的方法。利用正交变换，它可从高维数据集中提取低维特征集，同时保留尽可能多的数据信息，由于变量减少了，数据的可视化也更加清晰明阿川兰。在处理三维及以上维度的数据时，PCA不失为一种有效的方法。
由于PCA只在对称相关或协方差矩阵上执行，所以如果要用它，我们就必须要先有一个由数字组成的矩阵，而且必须是标准化的数据。
假设我们有一个维数为300 (n) × 50 (p)的数据集，n表示观测值的数量，p表示预测变量的数量。因为p = 50，这就意味着我们会有p(p-1)/2个散点图，要用超过1000个图表来分析变量关系。单从数据上看，这样的工作量就已经令人乏味。
在这种情况下，我们要做的是从预测变量p (p << 50)中抽取一个包含大部分数据信息的子集，然后用它在较低维的空间中绘制观察图。
下图显示了利用PCA将高维数据（三维）映射到低维（二维）的结果丝迪芬妮。牢记一点，每个合成维度是p个特征的线性组合。

来源：nlpca
什么是主成分？
主成分是数据集中原始预测变量的归一化线性组合。在上图中，PC1和PC2是主成分。
假设我们的原数据有p个特征，它们在空间中就有p维：X1, X2...,Xp，其主成分则可被表示为：
Z1= Φ11X1+ Φ21X2+ Φ31X3+ .... +Φp1Xp
其中，
Z1是第一主成分；
Φp1表示loadings(Φ1, Φ2…)，它们是主成分相对于各个维度坐标的比例，这些比例的平方和为1，因为如果数值过大禁忌咒纹，它们可能会导致更大的差异。它还定义了第一主成分Z1的方向，在p维空间中，它的方向就是距离n个观察值欧式距离最小的那条直线；
X1, X2...,Xp已经过归一化处理，它们的平均值是0，标准差是1。
第一主成分是原始特征维度的线性组合，它捕捉的是数据集中的最大方差，决定了数据的主要分布规律，即最高可变性方向。捕捉到的可变性越多，主成分得到的特征信息量就越多。相比其他主成分，第一主成分的信息量是最多的。
第一主成分产生一条最接近数据的直线，即它使数据点和直线之间的欧式距离之和最小。
第二主成分也是原始特征维度的线性组合，但它捕捉的是数据集中的剩余分布规律，和Z1无关。换句话说，第一主成分和第二主成分应该是正交的，它可以被表示为：
Z2= Φ12X1+ Φ22X2+ Φ32X3+ .... + Φp2Xp

以此类推，我们可以计算第三主成分、第四主成分……它们捕捉的都是剩余特征，与之前的主成分无关且正交。理论上来说，对于n×p维数据，我们可以创建min(n-1, p)个主成分。
可以发现，因变量Y与主成分方向无关，因此这是一种无监督的方法。为什么要对特征做归一化？
特征为原始预测提供了标准的预测空间，但它也存在诸多局限。设想一下，如果一个数据集的单位包括千米、升、千克、光年……它们数据尺度太大了，方差也太大了，如果在这个原始基础上执行PCA，不同特征自带的高方差会影响主成分判断整体数据的方差水平。这就有点得不偿失了。
如下图所示，我们在经过归一化和未经归一化的数据集上执行两次PCA。这个数据集有大约40个特征，可以看到，第一主成分应该是变量ItemMRP，第二个主成分则是变量ItemWeight。由于特征自带的方差高值，我们在右图中很难看出这两个主成分的优势，但在左图中，一切就很明显了。
在R＆Python中实现PCA
所以PCA的主题围绕的是坐标变换，网上家长学校对于数学狂热分子，可能他们会觉得徒手计算乐趣多多。但为了效率和质量，我们还是乖乖用函数吧。
下文演示的数据集是大商场销售：Big Mart Sales III
记住一点，PCA只能用于数值型数据，如果数据集中包含分类变量，请务必将其转换为数字。为了算法的顺利进行，数据清理是必须的！
#directory path
> path <- ".../Data/Big_Mart_Sales"
#set working directory
> setwd(path)
#load train and test file
> train <- read.csv("train_Big.csv")
> test <- read.csv("test_Big.csv")
#add a column
> test$Item_Outlet_Sales <- 1
#combine the data set
> combi <- rbind(train, test)
#impute missing values with median
> combi$Item_Weight[is.na(combi$Item_Weight)] <- median(combi$Item_Weight, na.rm = TRUE)
#impute 0 with median
> combi$Item_Visibility <- ifelse(combi$Item_Visibility == 0, median(combi$Item_Visibility),combi$Item_Visibility)
#find mode and impute
> table(combi$Outlet_Size, combi$Outlet_Type)
> levels(combi$Outlet_Size)[1] <- "Other"
到这里，我们已经用中位数补全数据集中的缺失值，解下来就是删除因变量和标识符变量（如果有的话）。正如我们上面所说的九女仙湖 ，这是一种无监督学习技术，因此必须删除因变量。
#remove the dependent and identifier variables
> my_data <- subset(combi, select = -c(Item_Outlet_Sales, Item_Identifier,
然后是检查数据集中的可用变量（各种预测特征）：
#check available variables
> colnames(my_data)
检查是否还包含非数值型变量：
#check variable class
> str(my_data)
'data.frame': 14204 obs. of 9 variables:
$ Item_Weight : num 9.3 5.92 17.5 19.2 8.93 ...
$ Item_Fat_Content : Factor w/ 5 levels "LF","low fat",..: 3 5 3 5 3 5 5 3 5 5 ...
$ Item_Visibility : num 0.016 0.0193 0.0168 0.054 0.054 ...
$ Item_Type : Factor w/ 16 levels "Baking Goods",..: 5 15 11 7 10 1 14 14 6 6 ...
$ Item_MRP : num 249.8 48.3 141.6 182.1 53.9 ...
$ Outlet_Establishment_Year: int 1999 2009 1999 1998 1987 2009 1987 1985 2002 2007 ...
$ Outlet_Size : Factor w/ 4 levels "Other","High",..: 3 3 3 1 2 3 2 3 1 1 ...
$ Outlet_Location_Type : Factor w/ 3 levels "Tier 1","Tier 2",..: 1 3 1 3 3 3 3 3 2 2 ...
$ Outlet_Type : Factor w/ 4 levels "Grocery Store",..: 2 3 2 1 2 3 2 4 2 2 ...
很不幸，这9个变量里有6个是绝对变量龚箭 ，比如年份等，无法被归一化。所以我们得先做一些工作，用一种热门编码方法将这些分类变量转换为数字。
#load library
> library(dummies)
#create a dummy data frame
> new_my_data <- dummy.data.frame(my_data, names = c("Item_Fat_Content","Item_Type",
"Outlet_Establishment_Year","Outlet_Size",
"Outlet_Location_Type"消失的女人,"Outlet_Type"))
检查一下这是不是个整数型数据集：
#check the data set
> str(new_my_data)
把数据集分为训练集和测试集：
#divide the new data
> pca.train <- new_my_data[1:nrow(train),]
> pca.test <- new_my_data[-(1:nrow(train)),]
在R语言中，使用PCA的函数是prcomp()，默认情况下它会把变量平均值处理成0，为了让它们的标准差等于1，我们还要加上scale. = T：
#principal component analysis
> prin_comp <- prcomp(pca.train, scale. = T)
> names(prin_comp)
[1] "sdev""rotation" "center" "scale" "x"
prcomp()函数产生了5个有用的维度：
1.center和scale表示在执行PCA前，变量的平均值和标准差情况；
#outputs the mean of variables
prin_comp$center
#outputs the standard deviation of variables
prin_comp$scale
2.这个rotation表示的是我们之前提到的原坐标的扭转比例，每一个特征（维度）都有一个相应比例，它们的线性组合就是我们的主成分；
> prin_comp$rotation
> prin_comp$rotation[1:5,1:4]
PC1PC2PC3 PC4
Item_Weight 0.0054429225 -0.001285666 0.011246194 0.011887106
Item_Fat_ContentLF -0.0021983314 0.003768557 -0.009790094 -0.016789483
Item_Fat_Contentlow fat -0.0019042710 0.001866905 -0.003066415 -0.018396143
Item_Fat_ContentLow Fat 0.0027936467 -0.002234328 0.028309811 0.056822747
Item_Fat_Contentreg 0.0002936319 0.001120931 0.009033254 -0.001026615
3.为了计算主成分得分向量，我们不需要将上述比例与数据相乘。相反地，我们可以理解为矩阵x在8523×44维空间中有一个主成分向量；
> dim(prin_comp$x)
[1] 8523 44
我们绘制合成主成分。
> biplot(prin_comp, scale = 0)

scale = 0这个参数就是我们之前提到的特征归一化处理。让我们关注图像顶部、底部、左侧、右侧的极值唔该借歪 ，就目测情况看，OutletTypeSupermarket、OutletEstablishmentYear 2007、OutletLocationTypeTier1和OutletSizeother都是极具潜力的第一、第二主成分，但具体我们还是要靠计算。
4.prcomp()还可以计算每个主成分的标准差：sdev；
#compute standard deviation of each principal component
> std_dev <- prin_comp$sdev
#compute variance
> pr_var <- std_dev^2
#check variance of first 10 components
> pr_var[1:10]
[1] 4.563615 3.217702 2.744726 2.541091 2.198152 2.015320 1.932076 1.256831
[9] 1.203791 1.168101
我们的目标是找到方差最大的成分，因为方差越大圣龙的共妻 ，主成分捕捉的特征信息量越多。这可以通过简单地将方差除以总方差看占比：
#proportion of variance explained
> prop_varex <- pr_var/sum(pr_var)
> prop_varex[1:20]
[1] 0.10371853 0.07312958 0.06238014 0.05775207 0.04995800 0.04580274
[7] 0.04391081 0.02856433 0.02735888 0.02654774 0.02559876 0.02556797
[13] 0.02549516 0.02508831 0.02493932 0.02490938 0.02468313 0.02446016
[19] 0.02390367 0.02371118
第一主成分占比10.3％，第二主成分占比7.3％，第三主成分占比6.2％……那我们建模要用多少主成分呢？这用一张碎石图就能解释：
#scree plot
> plot(prop_varex, xlab = "Principal Component",
ylab = "Proportion of Variance Explained",
type = "b")

上图显示前30个主成分在总方差中占比98.4％，考虑到之前我们一共可以构建44个主成分，用了PCA后，我们可以把它精简到30个，提高了效率：
#cumulative scree plot
> plot(cumsum(prop_varex), xlab = "Principal Component",
ylab = "Cumulative Proportion of Variance Explained",
type = "b")
用PCA成分进行预测建模
当我们获得主成分后，之后要做的是在测试数据进行预测。但这里有一些注意事项：
我们不能妄图把测试集和训练集整合成一个数据集，然后一次性获得所有PCA主成分。这样做的后果是测试数据会被“泄漏”到训练集中，换句话说，测试数据集将不再保持“不可见”，最后模型的通用性会变差。
我们也不能分别在测试集和训练集上单独执行PCA大隋天帝传。因为特征空间不同，最后它们获得的主成分方向也不一样我要去旅行 。简而言之，我们新建立起的两套坐标轴不是一个坐标轴。
以下是R语言的正确做法：
#add a training set with principal components
> train.data <- data.frame(Item_Outlet_Sales = train$Item_Outlet_Sales, prin_comp$x)
#we are interested in first 30 PCAs
> train.data <- train.data[,1:31]
#run a decision tree
> install.packages("rpart")
> library(rpart)
> rpart.model <- rpart(Item_Outlet_Sales ~ .,data = train.data, method = "anova")
> rpart.model
#transform test into PCA
> test.data <- predict(prin_comp, newdata = pca.test)
> test.data <- as.data.frame(test.data)
#select the first 30 components
> test.data <- test.data[,1:30]
#make prediction on test data
> rpart.prediction <- predict(rpart.model, test.data)
#For fun, finally check your score of leaderboard
> sample <- read.csv("SampleSubmission_TmnO39y.csv")
> final.sub <- data.frame(Item_Identifier = sample$Item_Identifier, Outlet_Identifier = sample$Outlet_Identifier, Item_Outlet_Sales = rpart.prediction)
> write.csv(final.sub, "pca.csv",row.names = F)
如果直接用这个东西参加比赛，性能会很差，建议你用随机森林重新优化。
下面是Python的正确做法：
import numpy as np
from sklearn.decomposition import PCA
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import scale
%matplotlib inline
#Load data set
data = pd.read_csv('Big_Mart_PCA.csv')
#convert it to numpy arrays
X=data.values
#Scaling the values
X = scale(X)
pca = PCA(n_components=44)
pca.fit(X)
#The amount of variance that each PC explains
var= pca.explained_variance_ratio_
#Cumulative Variance explains
var1=np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)
print var1
[ 10.37 17.68 23.92 29.7 34.7 39.28 43.67 46.53 49.27
51.92 54.48 57.04 59.59 62.1 64.59 67.08 69.55 72.
74.39 76.76 79.1 81.44 83.77 86.06 88.33 90.59 92.7
94.76 96.78 98.44 100.01 100.01 100.01 100.01 100.01 100.01
100.01 100.01 100.01 100.01 100.01 100.01 100.01 100.01]
plt.plot(var1）

#Looking at above plot I'm taking 30 variables
pca = PCA(n_components=30)
pca.fit(X)
X1=pca.fit_transform(X)
print X1
要点梳理
PCA能有效解决数据集中的特征冗余；
这些特征本质上是低维的；
PCA的主成分其实是原始特征维度经归一化后的线性组合；
这些主成分旨在准确描述数据分布和尽可能多地捕捉特征信息；
第一主成分方差最高，其次是第二主成分、第三主成分……
各主成分必须不相关（正交）；
预测特征变量是分类数据时，要编码处理再归一化；
PCA适用于3维及以上数据集；
PCA只能用于数值型数据集；
PCA是一种有助于生成更好的高维数据可视化的工具
原文地址：www.analyticsvidhya.com/blog/2016/03/practical-guide-principal-component-analysis-python/