PCA 降维分析

数据来源:

从这篇博客中借用下，百度云网页端居然可以直接下载文件了，稀奇

数据已经分享在百度云：客户年消费数据
密码：lehv
该份数据中包含客户id和客户6种商品的年消费额，共有440个样本

分析过程：

用python和R做个简单的PCA分析，顺便比较下结果是不是一致

numpy计算结果：

import numpy as np
import pandas as pd
data=pd.read_excel("客户年消费数据.xlsx")
data.head()

看一眼数据前几行：

#缺失检查
print(data.isnull().sum())

features = data[['Fresh', 'Milk', 'Grocery', 'Frozen', 'Detergents_Paper', 'Delicatessen']]
# 计算每一列的平均值及标准差
meandata = np.mean(features, axis=0) 
stddata = np.std(features, axis=0)
# 均值归一化
features = (features - meandata) / stddata    

# 求协方差矩阵
cov = np.cov(features.transpose()) # np.cov(features.T)

# 求解特征值和特征向量
eigVals, eigVectors = np.linalg.eig(cov) 
print('前两个主成分包含的信息百分比：{:.2%}'.format(np.sum(eigVals[:2])/np.sum(eigVals)))

# 选择前两个特征向量数据用于作图
pca_mat = eigVectors[:,:2]
pca_data = np.dot(features , pca_mat)
pca_data = pd.DataFrame(pca_data, columns=['pca1', 'pca2'])

结果：前两个主成分包含的信息百分比：72.46%

这里有些不严谨，特征值（PC解释的方差）eigVals为： [2.65099857 1.70646229 0.74175057 0.56501435 0.06311455 0.28632709]，第6个特征值比第五个特征值大，即方差更大，所以这里还需要按照解释的方差大小排序，特征向量也一样排序

# 计算特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(cov)

# 将特征向量按照特征值排序
sorted_indices = np.argsort(eigenvalues)[::-1]
sorted_eigenvalues = eigenvalues[sorted_indices]
sorted_eigenvectors = eigenvectors[:, sorted_indices]

此时方差解释率正常（从大到小）：[0.44082893 0.283764 0.12334413 0.09395504 0.04761272 0.01049519]

scikit-learn计算结果：

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import pandas as pd

# 读取数据
data=pd.read_excel("客户年消费数据.xlsx")
features = data[['Fresh', 'Milk', 'Grocery', 'Frozen', 'Detergents_Paper', 'Delicatessen']]

# 数据标准化
scaler = StandardScaler()
features_std = scaler.fit_transform(features)
# PCA计算
pca = PCA()
pca.fit(features_std)

# 获取主成分
principal_components = pca.transform(features_std)

# 创建一个新的数据框来存储主成分
principal_df = pd.DataFrame(data=principal_components, columns=['PC1', 'PC2', 'PC3', 'PC4','PC5','PC6'])

# 输出主成分的方差解释比例
explained_variance_ratio = pca.explained_variance_ratio_
print("主成分的方差解释比例:", explained_variance_ratio)

# 输出特征向量的解释方差
explained_variance = pca.explained_variance_
print("特征向量的解释方差:", explained_variance)

# 输出特征向量的累计方差解释比例
cumulative_explained_variance_ratio = np.cumsum(explained_variance_ratio)
print("特征向量的累计方差解释比例:", cumulative_explained_variance_ratio)

结果：
可以看到和numpy计算结果一样的，前两个PC解释率都是72.46%



碎石图和PCA结果图：

# 输出特征值（解释的方差）
print(sorted_eigenvalues)

fig,axx=plt.subplots(1,2,figsize=(8,3),gridspec_kw={"hspace": 0.6, "wspace": 0.4})
#ax.scatter(range(0,6),sorted_eigenvalues)
# 支持中文
plt.rc("font",family='KaiTi')
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

ax=axx[0]
ax.plot(range(1,7),sorted_eigenvalues ,color='r',marker='o',linewidth=1.5,alpha=0.7,linestyle='--')
ax.set_title('PCA碎石图')
ax.set_xlabel('PC主成分')
ax.set_ylabel('特征值')
ax.grid()
bx=axx[1]
bx.scatter(principal_df.iloc[:,0],principal_df.iloc[:,1],marker='+',alpha=0.5)
bx.set_xlabel('PC1({:.2%})'.format(explained_variance_ratio[0]))
bx.set_ylabel('PC2({:.2%})'.format(explained_variance_ratio[1]))
bx.set_title('PCA分析')
bx.grid()

R语言计算结果：

library(openxlsx)
data<-read.xlsx("客户年消费数据.xlsx")

# 去掉ID列 
df<-data[,-1]

# 数据归一化
df_scale<-scale(df)
# 计算协方差矩阵
#df_cov<-cov(df_scale)

# 计算特征值和特征向量   # 数据标准化后均值为0，方差为1，此时相关矩阵和协方差矩阵等价
value<-eigen(cor(df_scale))$value
#vec<-value<-eigen(cor(df_scale))$vectors

# 主成分方差解释比例
ratio<-value/sum(value)
#ratio结果： 0.44082893 0.28376400 0.12334413 0.09395504 0.04761272 0.01049519

# 累计解释率
cumsum(ratio)
#累计解释率：0.4408289 0.7245929 0.8479371 0.9418921 0.9895048 1.0000000

#或者使用函数 prcomp: 直接传入归一化的数据
sdev<-prcomp(df_scale)$sdev  # 和下面这句等价
#sdev<-prcomp(df, center = T, scale = T)$sdev

# 先计算特征值，即标准差的平方，然后计算各个方差解释率
ra<-sdev**2/sum(sdev**2)
# 累计解释率：
cumsum(ra)
# 结果：0.4408289 0.7245929 0.8479371 0.9418921 0.9895048 1.0000000

输出结果和python结果一致

ref: https://zhuanlan.zhihu.com/p/496823203