使用Numpy实现PCA

点击上方“小白学视觉”,选择加"星标"或“置顶

重磅干货,第一时间送达

作者 | Guillermina Sutter Schneider 

编译 | VK 

来源 | Towards Data Science

PCA是一种常用于处理多重共线性的特征提取方法。在这种情况下,PCA的最大优点是,在应用它之后,每个“新”变量将彼此独立。

本节基于mattbrems的这篇文章:https://towardsdatascience.com/a-one-stop-shop-for-principal-component-analysis-5582fb7e0b。我将一步一步地解释如何只使用numpy(以及使用一点pandas来操作数据帧)来进行PCA。

第1步

首先清理数据集非常重要。因为这不是本文的目标,你可以去仓库看看我做了什么:https://github.com/glosophy/WindPowerForecasting/blob/main/windPowerPCA.ipya

第2步

将数据分为因变量(Y)和特征或自变量(X)。

import pandas as pd

features = ['WindSpeed', 'RotorRPM', 'ReactivePower', 'GeneratorWinding1Temperature', 
       'GeneratorWinding2Temperature', 'GeneratorRPM', 'GearboxBearingTemperature', 'GearboxOilTemperature']

# 将数据分离为Y和X
y = data['ActivePower']
X = data[features]

第3步

取自变量X的矩阵,对于每一列,从每个特征中减去该列的平均值。(这样可以确保每列的平均值为零。)也可以通过除以每列的标准差来标准化X。此步骤的目的是标准化特征,使其平均值等于零,标准偏差等于1。

# 减去均值
X = X - X.mean()

# 标准化
Z = X / X.std()

第4步

这是一个健全性检查步骤。让我们确保平均值和标准差分别为0和1。

# 检验均值= 0和标准差= 1
print('MEAN:')
print(Z.mean())
print('---'*15)
print('STD:')
print(Z.std())

第5步

取矩阵Z,转置它,然后将转置后的矩阵乘以Z,这就是Z的协方差矩阵。

import numpy as np

Z = np.dot(Z.T, Z)

第6步

计算的特征值数组和一个特征矩阵,特征矩阵的列是与特征值对应的归一化特征向量。

在这一步中,重要的是要确保特征值及其特征向量按降序排序(从大到小)。对特征值进行排序,然后对特征向量进行相应排序。

eigenvalues, eigenvectors = np.linalg.eig(Z)

第7步

把特征向量的矩阵赋给P,把对角矩阵赋给D,特征值在对角线上,其它地方的值都为零。D对角线上的特征值将与P中相应的列相关联。

D = np.diag(eigenvalues)
P = eigenvectors

第8步

计算Z* = ZP。这个新的矩阵,Z*,是X的中心或标准化版本,但现在每个观测值都是原始变量的组合,其中权重由特征向量确定。

关于这个新矩阵Z*的一个重要的事情是,因为P中的特征向量彼此独立,所以Z*中的列也是相互独立的!

Z_new = np.dot(Z, P)

第9步

有了Z*你就可以决定保留多少特征与删除多少特征。这通常是通过一个scree图来完成的。

scree图显示了每个主成分从数据中捕捉到的变化量。y轴代表变化量(有关scree图以及如何解释它们的更多信息,请参阅这篇文章:https://bioturing.medium.com/how-to-read-pca-biplots-and-screeplot-186246aae063#:~:text=A%20scree%20plot%20shows%20how,the%20principal%20components%20to%20keep.&text=Proportion%20of%20variance%20plot%3A%20the,least%2080%25%20of%20the%20variance.)。

左边的曲线图显示了从每个主成分中得出的变化量,以及为模型添加或考虑另一个主成分而产生的累积变化量。

选择多少个主成分,归根结底是一个经验法则:所选的主成分应该能够描述至少80-85%的方差。在本例中,仅第一个主成分就解释了大约82%的方差。加上第二个主成分,这个数字几乎达到90%。

#1. 计算每个特征解释的方差的比例
sum_eigenvalues = np.sum(eigenvalues)

prop_var = [i/sum_eigenvalues for i in eigenvalues]

#2. 计算累积方差
cum_var = [np.sum(prop_var[:i+1]) for i in range(len(prop_var))]


# 导入plt
import matplotlib.pyplot as plt

x_labels = ['PC{}'.format(i+1) for i in range(len(prop_var))]

plt.plot(x_labels, prop_var, marker='o', markersize=6, color='skyblue', linewidth=2, label='Proportion of variance')
plt.plot(x_labels, cum_var, marker='o', color='orange', linewidth=2, label="Cumulative variance")
plt.legend()
plt.title('Scree plot')
plt.xlabel('Principal components')
plt.ylabel('Proportion of variance')
plt.show()

如果你想知道如何在数据上下文中解释,我发现这篇文章特别容易理解:https://blogs.sas.com/content/iml/2019/11/04/interpret-graph-principal-component.html。他们使用iris数据集,并使用scree、profile和pattern图给出了许多示例。

结论

PCA有时很棘手。但只要有适当的资源和耐心,它就可以变得令人愉快。最重要的是,我发现回到基础知识对于全面掌握PCA是很有用的。

希望你觉得这篇文章有帮助!你可以在这里看到完整的Python脚本:https://github.com/glosophy/WindPowerForecast/blob/main/windPowerPCA.ipyn

下载1:OpenCV-Contrib扩展模块中文版教程

在「小白学视觉」公众号后台回复:扩展模块中文教程即可下载全网第一份OpenCV扩展模块教程中文版,涵盖扩展模块安装、SFM算法、立体视觉、目标跟踪、生物视觉、超分辨率处理等二十多章内容。

下载2:Python视觉实战项目52讲

在「小白学视觉」公众号后台回复:Python视觉实战项目即可下载包括图像分割、口罩检测、车道线检测、车辆计数、添加眼线、车牌识别、字符识别、情绪检测、文本内容提取、面部识别等31个视觉实战项目,助力快速学校计算机视觉。

下载3:OpenCV实战项目20讲

在「小白学视觉」公众号后台回复:OpenCV实战项目20讲即可下载含有20个基于OpenCV实现20个实战项目,实现OpenCV学习进阶。

交流群

欢迎加入公众号读者群一起和同行交流,目前有SLAM、三维视觉、传感器、自动驾驶、计算摄影、检测、分割、识别、医学影像、GAN、算法竞赛等微信群(以后会逐渐细分),请扫描下面微信号加群,备注:”昵称+学校/公司+研究方向“,例如:”张三 + 上海交大 + 视觉SLAM“。请按照格式备注,否则不予通过。添加成功后会根据研究方向邀请进入相关微信群。请勿在群内发送广告,否则会请出群,谢谢理解~

  • 2
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
好的,首先我们先来了解一下PCA降维的原理。PCA全称为Principal Component Analysis,即主成分分析,它是一种常用的数据降维方法,可以将高维数据降到低维空间中。PCA的核心思想是通过线性变换将原来的高维数据转换为一组新的低维度的数据,这组新的数据被称为主成分。通常我们只选择前几个主成分进行保留,而抛弃其余的主成分,从而达到降维的目的。 下面我们就来实现一下PCA降维,并对手写数字进行分类。我们将使用sklearn库中自带的手写数字数据集进行操作。 ``` import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import load_digits from sklearn.model_selection import train_test_split from sklearn.decomposition import PCA from sklearn.neighbors import KNeighborsClassifier # 加载手写数字数据集 digits = load_digits() X = digits.data y = digits.target # 将数据集分为训练集和测试集 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42) # 对数据进行归一化处理 mean = np.mean(X_train, axis=0) std = np.std(X_train, axis=0) X_train = (X_train - mean) / std X_test = (X_test - mean) / std # 使用PCA进行降维 pca = PCA(n_components=2) X_train_pca = pca.fit_transform(X_train) X_test_pca = pca.transform(X_test) # 使用KNN分类器进行分类 knn = KNeighborsClassifier(n_neighbors=3) knn.fit(X_train_pca, y_train) # 计算准确率 accuracy = knn.score(X_test_pca, y_test) print("Accuracy:", accuracy) # 可视化降维后的数据 plt.figure(figsize=(8, 6)) for i in range(10): plt.scatter(X_train_pca[y_train==i, 0], X_train_pca[y_train==i, 1], label=str(i), alpha=0.8) plt.xlabel("Principal Component 1") plt.ylabel("Principal Component 2") plt.legend() plt.show() ``` 在上述代码中,我们首先加载了手写数字数据集,并将其分为训练集和测试集。然后对数据进行了归一化处理,接着使用PCA将数据降到2维空间中。最后使用KNN分类器进行分类,并计算了准确率。同时,我们也将降维后的数据可视化出来,可以看到不同数字之间在2维空间中的分布情况。 下面我们来运行一下这段代码,看看最终结果如何。 ``` Accuracy: 0.562962962962963 ``` 可以看到,准确率并不是很高,这是因为我们只选择了2个主成分进行保留,而且KNN分类器对于高维度的数据表现并不是很好。如果我们选择更多的主成分进行保留,或者使用其他的分类器,可能会得到更好的结果。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值