使用Numpy實(shí)現(xiàn)PCA

重磅干貨，第一時(shí)間送達(dá)

作者 | Guillermina Sutter Schneider 

編譯 | VK 

來源 | Towards Data Science

PCA是一種常用于處理多重共線性的特征提取方法。在這種情況下，PCA的最大優(yōu)點(diǎn)是，在應(yīng)用它之后，每個(gè)“新”變量將彼此獨(dú)立。

本節(jié)基于mattbrems的這篇文章:https://towardsdatascience.com/a-one-stop-shop-for-principal-component-analysis-5582fb7e0b。我將一步一步地解釋如何只使用numpy（以及使用一點(diǎn)pandas來操作數(shù)據(jù)幀）來進(jìn)行PCA。

首先清理數(shù)據(jù)集非常重要。因?yàn)檫@不是本文的目標(biāo)，你可以去倉(cāng)庫(kù)看看我做了什么：https://github.com/glosophy/WindPowerForecasting/blob/main/windPowerPCA.ipya

將數(shù)據(jù)分為因變量（Y）和特征或自變量（X）。

import pandas as pd

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

# 將數(shù)據(jù)分離為Y和X
y = data['ActivePower']
X = data[features]

取自變量X的矩陣，對(duì)于每一列，從每個(gè)特征中減去該列的平均值。（這樣可以確保每列的平均值為零。）也可以通過除以每列的標(biāo)準(zhǔn)差來標(biāo)準(zhǔn)化X。此步驟的目的是標(biāo)準(zhǔn)化特征，使其平均值等于零，標(biāo)準(zhǔn)偏差等于1。

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

# 標(biāo)準(zhǔn)化
Z = X / X.std()

這是一個(gè)健全性檢查步驟。讓我們確保平均值和標(biāo)準(zhǔn)差分別為0和1。

# 檢驗(yàn)均值= 0和標(biāo)準(zhǔn)差= 1
print('MEAN:')
print(Z.mean())
print('---'*15)
print('STD:')
print(Z.std())

取矩陣Z，轉(zhuǎn)置它，然后將轉(zhuǎn)置后的矩陣乘以Z，這就是Z的協(xié)方差矩陣。

import numpy as np

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

計(jì)算的特征值數(shù)組和一個(gè)特征矩陣，特征矩陣的列是與特征值對(duì)應(yīng)的歸一化特征向量。

在這一步中，重要的是要確保特征值及其特征向量按降序排序（從大到小）。對(duì)特征值進(jìn)行排序，然后對(duì)特征向量進(jìn)行相應(yīng)排序。

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

把特征向量的矩陣賦給P，把對(duì)角矩陣賦給D，特征值在對(duì)角線上，其它地方的值都為零。D對(duì)角線上的特征值將與P中相應(yīng)的列相關(guān)聯(lián)。

D = np.diag(eigenvalues)
P = eigenvectors

計(jì)算Z* = ZP。這個(gè)新的矩陣，Z*，是X的中心或標(biāo)準(zhǔn)化版本，但現(xiàn)在每個(gè)觀測(cè)值都是原始變量的組合，其中權(quán)重由特征向量確定。

關(guān)于這個(gè)新矩陣Z*的一個(gè)重要的事情是，因?yàn)镻中的特征向量彼此獨(dú)立，所以Z*中的列也是相互獨(dú)立的！

Z_new = np.dot(Z, P)

有了Z*你就可以決定保留多少特征與刪除多少特征。這通常是通過一個(gè)scree圖來完成的。

scree圖顯示了每個(gè)主成分從數(shù)據(jù)中捕捉到的變化量。y軸代表變化量（有關(guān)scree圖以及如何解釋它們的更多信息，請(qǐng)參閱這篇文章：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.）。

左邊的曲線圖顯示了從每個(gè)主成分中得出的變化量，以及為模型添加或考慮另一個(gè)主成分而產(chǎn)生的累積變化量。

選擇多少個(gè)主成分，歸根結(jié)底是一個(gè)經(jīng)驗(yàn)法則：所選的主成分應(yīng)該能夠描述至少80-85%的方差。在本例中，僅第一個(gè)主成分就解釋了大約82%的方差。加上第二個(gè)主成分，這個(gè)數(shù)字幾乎達(dá)到90%。

#1. 計(jì)算每個(gè)特征解釋的方差的比例
sum_eigenvalues = np.sum(eigenvalues)

prop_var = [i/sum_eigenvalues for i in eigenvalues]

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


# 導(dǎo)入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()

如果你想知道如何在數(shù)據(jù)上下文中解釋，我發(fā)現(xiàn)這篇文章特別容易理解：https://blogs.sas.com/content/iml/2019/11/04/interpret-graph-principal-component.html。他們使用iris數(shù)據(jù)集，并使用scree、profile和pattern圖給出了許多示例。

PCA有時(shí)很棘手。但只要有適當(dāng)?shù)馁Y源和耐心，它就可以變得令人愉快。最重要的是，我發(fā)現(xiàn)回到基礎(chǔ)知識(shí)對(duì)于全面掌握PCA是很有用的。

希望你覺得這篇文章有幫助！你可以在這里看到完整的Python腳本：https://github.com/glosophy/WindPowerForecast/blob/main/windPowerPCA.ipyn