PCA（主成分分析）
2024-05-25 16:37:16
PCA（主成分分析）
1.概述
- 通过计算，计算出主要的“特征值”，达到减少特征值个数的目的。（这里的特征值并不是原来的特征值，而是经过投影后的新特征值）
- 对数据的要求：PCA成立的假设前提是数据应该满足正交分布并且是高信噪比的（低通滤波特性），而这两个特点都正是高斯分布的特点，所以，PCA输入需要数据满足高斯分布。或者说，高斯分布是PCA适用的充分条件，但不是必要条件
- 优点：理论依据充分
- 缺点：由于降维过程是纯数学分析，可能丢失掉数据的实际代表含义
- 作用：分类、对指标进行建模时候的客观赋权（如最小风险贝叶斯中的风险函数）、构造新特征（因为得到了原来数据的投影）、异常数据监测（因为计算出了源数据的方差，方差大的是好数据，方差小的是噪声数据）
2.原理
采用了yale人脸数据集作为手动计算案例
步骤
1.计算均值向量u
2.计算中心化矩阵C
3.计算协方差矩阵Sigma
4.计算变换矩阵W
5.计算投影Y
6.计算测试集投影Z
7.遍历搜索匹配，预测Z的标签
ps.步骤6和步骤7是预测步骤，降维结果是Y
2.1 手动计算见：pca主成分分析.pdf
3.手动实现PCA算法
- 采用了Iris数据集和DryBean数据集
- 也进行了均衡标签/不均衡标签的划分
import pandas as pd
import numpy as np
import struct
import os
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.decomposition import PCA
#引入数据集
iris = datasets.load_iris()
x= iris.data
y= iris.target
#拆分训练集和测试集
split = 0.7
#trianset : testset = 7:3np.random.seed(0) #固定随机结果
train_indices = np.random.choice(len(x),round(len(x) * split),replace=False)
test_indices = np.array(list(set(range(len(x))) - set(train_indices)))
np.random.seed(0) #固定随机结果
train_indices = np.random.choice(len(x),round(len(x) * split),replace=False)
test_indices = np.array(list(set(range(len(x))) - set(train_indices)))train_indices = sorted(train_indices)
test_indices =sorted(test_indices)
train_x = x[train_indices]
test_x = x[test_indices]
train_y = y[train_indices]
test_y = y[test_indices]
# print(train_indices)
# print(test_indices)
#PCA核心计算过程
def MyPca(X,k):U = np.mean(X,axis=0)C = X - Un = len(X)
# print(n)Cov = np.dot(C.T,C) / neigval, eigvec = np.linalg.eig(Cov)indexes = np.argsort(eigval)[-k:]W = eigvec[:, indexes]Y = np.dot(C, W).T
# print(pd.DataFrame(W).shape)
# print(pd.DataFrame(Y).shape)return Y,W
#使用了矩阵乘法计算欧几里得距离，缩短计算时间
#计算原理见附件：矩阵间的欧氏距离.pdf
def test(test_x,Y,test_y,train_y,W):n = len(test_x)Z = test_xU = np.mean(Z, axis=0)C = Z - UZZ = np.dot(C, W).T # k * n#求出ZZ的每一行与Y的每一行的欧式距离#ZZ的第一行表示第一个待测样本#对ZZ的第一行进行广播，与Y的每一行求欧式距离
# print(ZZ.T.shape,Y.shape)dists = np.sqrt(-2 * np.dot(ZZ.T, Y) + np.repeat([np.sum(np.square(ZZ.T), axis=1)],repeats=Y.shape[1],axis=0).T + np.repeat(np.transpose([np.sum(np.square(Y), axis=0)]),repeats=ZZ.shape[1],axis=1).T)# print(dists.shape)# print(np.argmin(dists, axis=1))pre = np.argmin(dists, axis=1)right = 0test_y = pd.DataFrame(test_y)train_y = pd.DataFrame(train_y)for i in range(len(pre)):
# print(i,pre[i],len(test_y),len(train_y))
# print(test_y[i],train_y.iloc[1212])#ub_train_y.iloc[1212][0]if test_y.iloc[i][0] == train_y.iloc[pre[i]][0]:right=right+1return right/n
#抽取不均衡标签的数据集
from numpy import *
i_train_x = train_x
i_train_y = train_y
i_test_x = test_x
i_test_y = test_y
ui_train_x = concatenate((concatenate((x[:45],x[50:95]),axis=0),x[100:115]),axis=0) #axis=0表示竖向拼接
ui_train_y = concatenate((concatenate((y[:45],y[50:95]),axis=0),y[100:115]),axis=0)
ui_test_x = concatenate((concatenate((x[:5],x[50:55]),axis=0),x[100:135]),axis=0)
ui_test_y = concatenate((concatenate((y[:5],y[50:55]),axis=0),y[100:135]),axis=0)
#比较均衡标签和不均衡标签的效果
klist = []
i_score = []
ui_score = []
for k in range(1,5): klist.append(k)Y,W = MyPca(i_train_x,k)ret = test(i_test_x,Y,i_test_y,i_train_y,W)i_score.append(ret)Y,W = MyPca(ui_train_x,k)ret = test(ui_test_x,Y,ui_test_y,ui_train_y,W)ui_score.append(ret)plt.plot(klist, i_score, marker = 'o', label = 'banlanced Iris')
plt.plot(klist, ui_score, marker = '*', label = 'unbanlanced Iris')
plt.legend() #让图例生效
plt.xlabel('k-value')
plt.ylabel('accuracy-value')
plt.title(u'Iris map')
plt.show()
3.1 接下来验证自我实现的结果与调包的结果一样
- 取k=3，画出3D图来展示
- 发现计算结果仅有方向上的差别
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from matplotlib import cm
klist = []
i_score = []
ui_score = []
for k in range(3,4): klist.append(k)Y,W = MyPca(i_train_x,k)fig = plt.figure()ax = fig.add_subplot(projection='3d')ax.scatter3D(Y[0,: ], Y[1,: ] ,Y[2,:] * -1) # 三个数组对应三个维度（三个数组中的数一一对应）plt.xlabel('X')plt.ylabel('Y', rotation=38) # y 轴名称旋转 38 度ax.set_zlabel('Z')pca = PCA(n_components=k)X_new = pca.fit_transform(i_train_x)fig = plt.figure()ax = fig.add_subplot(projection='3d')ax.scatter3D(X_new[:, 2], X_new[:, 1] * -1,X_new[:,0] * -1) # 三个数组对应三个维度（三个数组中的数一一对应）# print(Y[0,:].flatten().shape,Y[1,:].flatten().shape,Y[2,:].flatten().shape)# surf = ax.plot_trisurf(Y[0,:].tolist(),Y[1,:].tolist(),Y[2,:].tolist(),cmap=cm.coolwarm,linewidth=0, antialiased=False)plt.xlabel('X')plt.ylabel('Y', rotation=38) # y 轴名称旋转 38 度ax.set_zlabel('Z') # 因为 plt 不能设置 z 轴坐标轴名称，所以这里只能用 ax 轴来设置（当然，x 轴和 y 轴的坐标轴名称也可以用 ax 设置）plt.savefig('3D-2.jpg', bbox_inches='tight', dpi=2400) # 保存图片，如果不设置 bbox_inches='tight'，保存的图片有可能显示不全plt.show()
3.2 用DryBean数据集试一试，看看效果
import operator
from sklearn.preprocessing import StandardScaler # 均值归一化
from sklearn.metrics import confusion_matrix # 生成混淆矩阵
from sklearn.metrics import classification_report # 分类报告
def openfile(filename):"""打开数据集，进行数据处理:param filename:文件名:return:特征集数据、标签集数据""" # 打开excelsheet = pd.read_excel(filename,sheet_name='Dry_Beans_Dataset')data = sheet.iloc[:,:16].valuestarget = sheet['Class'].values
# print(data.shape)
# print(target.shape) return data, target, sheet.columns
def split_data_set(data_set, target_set, rate=0.7):"""说明：分割数据集，默认数据集的30%是测试集:param data_set: 数据集:param target_set: 标签集:param rate: 测试集所占的比率:return: 返回训练集数据、训练集标签、测试集数据、测试集标签"""# 计算训练集的数据个数train_size = len(data_set)# 随机获得数据的下标train_index = sorted(np.random.choice(train_size,round(train_size * rate), replace=False))test_index = sorted(np.array(list(set(range(train_size)) - set(train_index)))) #不用排序也行，强迫症，为了上面保持一致就排序了# 分割数据集（X表示数据，y表示标签）x_train = data_set.iloc[train_index,:] #因为这里的data_set和target_set变成DataFrame，而不是ndarray了，所以要用iloc访问x_test = data_set.iloc[test_index,:]y_train = target_set.iloc[train_index,:]y_test = target_set.iloc[test_index,:]return x_train, y_train, x_test, y_test
filename = r'D:\jjq\code\jupyterWorkSpace\datasets\DryBeanDataset\Dry_Bean_Dataset.xlsx'
o_bean_dataset = openfile(filename)
#每个类别的种子抽取400条数据，这个是每个类别的起始索引
step = 400
start_index = [0,1322,1844,3474,7020,8948,10975] #一共7类
# bean_dataset_x = pd.DataFrame(columns=o_bean_dataset[2])
# bean_dataset_y =pd.DataFrame(columns=o_bean_dataset[2])
bean_dataset_x = pd.DataFrame(columns=range(16))
bean_dataset_y =pd.DataFrame(columns=range(1))
bean_dataset_x.drop(bean_dataset_x.index,inplace=True)
bean_dataset_y.drop(bean_dataset_y.index,inplace=True)
for i in range(7):bean_dataset_x = pd.concat((bean_dataset_x, pd.DataFrame(o_bean_dataset[0][start_index[i]:(step+start_index[i])])),axis=0)bean_dataset_y = pd.concat((bean_dataset_y, pd.DataFrame(o_bean_dataset[1][start_index[i]:(step+start_index[i])])),axis=0)
# bean_dataset_y.to_excel("./123.xlsx")
#按照均衡和不均衡的方式，划分训练集和测试集#均衡
b_train_x, b_train_y, b_test_x, b_test_y = split_data_set(bean_dataset_x,bean_dataset_y)
# print(b_train_x.shape,b_train_y.shape)
# print(b_test_x.shape,b_test_y.shape)
#不均衡
steps_train = [318,318,318,318,318,212,159]
steps_test = [68,91,136,136,136,136,136]
now = 0
#初始化不均衡数组
ub_train_x = pd.DataFrame(columns=range(16))
ub_test_x = pd.DataFrame(columns=range(16))
ub_train_y = pd.DataFrame(columns=range(1))
ub_test_y = pd.DataFrame(columns=range(1))
#保证添加数据之前数组为空
ub_train_x.drop(ub_train_x.index,inplace=True)
ub_test_x.drop(ub_test_x.index,inplace=True)
ub_train_y.drop(ub_train_y.index,inplace=True)
ub_test_y.drop(ub_test_y.index,inplace=True)#开始添加数据
for i in range(7):ub_train_x = pd.concat((ub_train_x, bean_dataset_x[now:(now+steps_train[i])]),axis=0) ub_train_y = pd.concat((ub_train_y, bean_dataset_y[now:(now+steps_train[i])]),axis=0)now = now+steps_train[i]ub_test_x = pd.concat((ub_test_x, bean_dataset_x[now:(now+steps_test[i])]),axis=0)ub_test_y = pd.concat((ub_test_y, bean_dataset_y[now:(now+steps_test[i])]),axis=0)now = now+steps_test[i]
b_score = []
ub_score = []
klist = []
for k in range(1,10): klist.append(k)Y,W = MyPca(b_train_x,k)ret = test(b_test_x,Y,b_test_y,b_train_y,W)b_score.append(ret)Y,W = MyPca(ub_train_x,k)ret = test(ub_test_x,Y,ub_test_y,ub_train_y,W)ub_score.append(ret)plt.plot(klist, b_score, marker = 'o', label = 'banlanced DryBean')
plt.plot(klist, ub_score,marker = '*', label = 'unbalanced DryBean')
plt.legend() #让图例生效
plt.xlabel('k-value')
plt.ylabel('accuracy-value')
plt.title(u'DryBean map')
plt.show()
#这里发现，Iris数据集和DryBean数据集都是均衡标签的数据效果更好
#接下来寻找原因
4.寻找效果不同的原因
4.1首先查看两个数据集的分布
- 发现分布一样
from sklearn.preprocessing import StandardScaler # 均值归一化
scaler = StandardScaler()
b_train_x = pd.DataFrame(scaler.fit_transform(b_train_x))
plt.figure(figsize=(20,30))
p = 1
k=14
for i in range(k):a = plt.subplot(k,2,p)a.plot(range(b_train_x.shape[0]),b_train_x.iloc[:,i].tolist(),label="balance "+str(i),color='red')p = p+1plt.legend() #让图例生效b = plt.subplot(k,2,p)b.plot(range(ub_train_x.shape[0]),ub_train_x.iloc[:,i].tolist(),label="unbalance "+str(i))p = p+1plt.legend() #让图例生效plt.xlabel('index')
plt.ylabel('value')
plt.title(u'DryBean map')
plt.show()
# plt.subplots_adjust(hspace=100, wspace=10)
plt.figure(figsize=(15,40))
p = 1
for i in range(14):a = plt.subplot(14,2,p)a.set_title(b_train_x.columns.values.tolist()[i])a.hist(x = b_train_x.iloc[:,i].tolist(), # 指定绘图数据bins = 50, # 指定直方图中条块的个数color='red')p=p+1a = plt.subplot(14,2,p)a.set_title(ub_train_x.columns.values.tolist()[i])a.hist(x = ub_train_x.iloc[:,i].tolist(), # 指定绘图数据bins = 50, # 指定直方图中条块的个数)p=p+1
plt.show()
plt.figure(figsize=(15,5))
b = plt.subplot(1,2,1)
b.hist(x = b_test_y.iloc[:,:].values[:,0].tolist(), # 指定绘图数据bins = 50, # 指定直方图中条块的个数color="red")
c = plt.subplot(1,2,2)
c.hist(x = ub_test_y.iloc[:,:].values[:,0].tolist(), # 指定绘图数据bins = 50, # 指定直方图中条块的个数)
plt.show()
4.2根据理论，pca需要符合高斯分布的数据
- 使用函数，计算两个数据集的“正态分数”，发现效果更好的数据确实更接近正态分布一些（强行解释）
#DryBean数据集
from scipy.stats import normaltest
k2, pp = normaltest(b_train_x)
k3,ppp = normaltest(ub_train_x)
plt.plot(range(len(pp)),pp,label="balance")
plt.plot(range(len(ppp)),ppp,label="ubbalance")
plt.legend()
plt.show()
# print(k2-k3,pp-ppp)
#统计量p越接近0就越表明数据和标准正态分布拟合的越好
#Iris数据集
from scipy.stats import normaltest
k2, pp = normaltest(i_train_x)
k3,ppp = normaltest(ui_train_x)
plt.plot(range(len(pp)),pp,label="balance")
plt.plot(range(len(ppp)),ppp,label="ubbalance")
plt.legend()
plt.show()
5.sklearn中PCA
- 原型：sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False)
5.1属性
- components_ ：返回具有最大方差的成分。
- explained_variance_ratio_：返回 所保留的n个成分各自的方差百分比。
- n_components_：返回所保留的成分个数n。
- mean_：均值
- noise_variance_：噪声比
5.2 方法
- fit(X,y=None)：fit(X)，表示用数据X来训练PCA模型。函数返回值：调用fit方法的对象本身。比如pca.fit(X)，表示用X对pca这个对象进行训练
- fit_transform(X)：用X来训练PCA模型，同时返回降维后的数据
- inverse_transform()：将降维后的数据转换成原始数据，X=pca.inverse_transform(newX)
- transform(X)：将数据X转换成降维后的数据。当模型训练好后，对于新输入的数据，都可以用transform方法来降维
PCA（主成分分析）相关推荐
- R语言PCA主成分分析（Principle Component Analysis）实战2
R语言PCA主成分分析(Principle Component Analysis)实战2 目录 R语言PCA主成分分析(Principle Component Analysis)实战2 #案例分析
- R语言PCA主成分分析（Principle Component Analysis）与线性回归结合实战
R语言PCA主成分分析(Principle Component Analysis)与线性回归结合实战 目录 R语言PCA主成分分析(Principle Component Analysis)与线性回归 ...
- R语言PCA主成分分析（Principle Component Analysis）实战1
R语言PCA主成分分析(Principle Component Analysis)实战1 目录 R语言PCA主成分分析(Principle Component Analysis)实战1 #案例分析
- 【数学与算法】PCA主成分分析(降维)的通俗理解
1.PCA降维 PCA主成分分析简单的理解,就是把某物的很多个能直接获取到的特征,经过变换得到很多个新特征,这些新特征对该物体来说,有的影响很大,有的影响很小,只需要使用这些影响大的新特征,舍弃很多影 ...
- 统计学习方法第十六章作业：PCA主成分分析算法 代码实现
PCA主成分分析 import numpy as np class PCA:def __init__(self,x,R=None):self.x = np.array(x)self.dim = sel ...
- PCA主成分分析 特征降维 opencv实现
最近对PCA主成分分析做了一定的了解,对PCA基础和简单的代码做了小小的总结 有很多博客都做了详细的介绍,这里也参考了这些大神的成果: http://blog.sina.com.cn/s/blog_7 ...
- pca 主成分分析_超越普通PCA：非线性主成分分析
pca 主成分分析 TL;DR: PCA cannot handle categorical variables because it makes linear assumptions about t ...
- pca 主成分分析_六分钟的主成分分析（PCA）的直观说明。
pca 主成分分析 Principle Component Analysis (PCA) is arguably a very difficult-to-understand topic for be ...
- canoco5主成分分析步骤_R语言 PCA主成分分析
微信公众号:生信小知识 关注可了解更多的教程及生信知识.问题或建议,请公众号留言; R语言 PCA主成分分析 前言统计学背景知识协方差相关系数函数总结实例讲解1.载入原始数据2.作主成分分析3.结果解 ...
- 无监督学习 | PCA 主成分分析之客户分类
文章目录 1. 开始 2. 数据探索 2.2 特征相关性 2.3 可视化特征分布 3. 数据预处理 3.1 特征缩放 3.2 异常值检测 4. 数据转换 4.1 主成分分析(PCA) 4.2 降维 4 ...
最新文章
- pycharm重点插件
- 【模板】树链剖分 P3384
- OSS正式支持IPv6公测
- 简单几步搞定ISA ×××
- 正常情况下ffmpeg生成moov是在mdat写完成之后写入
- matlab 7.9.0 帮助翻译--zeros函数
- 阿里云ACM：云原生配置管理利器，让云上的Spring Cloud应用配置管理舞动起来
- windows下 Mysql 错误1067 Can't open and lock privilege tables: Table 'mysql.user' doesn't exist
- 【教程】如何在C＃中创建PDF417条码
- 为什么c语言输出到文件慢,【图片】今天写几个性能测试，为什么C语言跑得这么慢呢？？【c语言吧】_百度贴吧...
- php symlink,php函数symlink详解
- 转： qemu-kvm内存管理
- logo设计灵感的创意网站
- Java基础8顺序语句判断语句
- Error at hooking API “LoadStringA“ Dump first 32 bytes:
- AI：人工智能技术层企业简介(更新中)
- 以太网UDP数据协议
- python练习, 打鱼晒网问题
- 如何在rhel4上禁用不需要的相关服务
- Nginx 配置 SSL 证书 + HTTPS 站点小记