使用 matplotlib 的 PCA 基本示例

Posted 2023-02-16

【中文标题】使用 matplotlib 的 PCA 基本示例

【英文标题】：Basic example for PCA with matplotlib

【发布时间】：2013-08-20 09:46:15

【问题描述】：

我尝试使用matplotlib.mlab.PCA 进行简单的主成分分析，但使用该类的属性我无法为我的问题找到一个干净的解决方案。这是一个例子：

获取一些 2D 虚拟数据并启动 PCA：

from matplotlib.mlab import PCA
import numpy as np

N     = 1000
xTrue = np.linspace(0,1000,N)
yTrue = 3*xTrue

xData = xTrue + np.random.normal(0, 100, N)
yData = yTrue + np.random.normal(0, 100, N)
xData = np.reshape(xData, (N, 1))
yData = np.reshape(yData, (N, 1))
data  = np.hstack((xData, yData))
test2PCA = PCA(data)

现在，我只想将主成分作为原始坐标中的向量，并将它们作为箭头绘制到我的数据上。

什么是到达那里的快速而干净的方式？

谢谢，泰拉克斯

【参考方案1】：

我认为mlab.PCA 类不适合您想做的事情。特别是，PCA 类在找到特征向量之前重新调整数据：

a = self.center(a)
U, s, Vh = np.linalg.svd(a, full_matrices=False)

center 方法除以sigma：

def center(self, x):
    'center the data using the mean and sigma from training set a'
    return (x - self.mu)/self.sigma

这会产生特征向量pca.Wt，如下所示：

[[-0.70710678 -0.70710678]
 [-0.70710678  0.70710678]]

它们是垂直的，但与原始数据的主轴没有直接关系。它们是按摩数据的主轴。

也许直接编写您想要的代码可能更容易（不使用mlab.PCA 类）：

import numpy as np
import matplotlib.pyplot as plt

N = 1000
xTrue = np.linspace(0, 1000, N)
yTrue = 3 * xTrue
xData = xTrue + np.random.normal(0, 100, N)
yData = yTrue + np.random.normal(0, 100, N)
xData = np.reshape(xData, (N, 1))
yData = np.reshape(yData, (N, 1))
data = np.hstack((xData, yData))

mu = data.mean(axis=0)
data = data - mu
# data = (data - mu)/data.std(axis=0)  # Uncommenting this reproduces mlab.PCA results
eigenvectors, eigenvalues, V = np.linalg.svd(data.T, full_matrices=False)
projected_data = np.dot(data, eigenvectors)
sigma = projected_data.std(axis=0).mean()
print(eigenvectors)

fig, ax = plt.subplots()
ax.scatter(xData, yData)
for axis in eigenvectors:
    start, end = mu, mu + sigma * axis
    ax.annotate(
        '', xy=end, xycoords='data',
        xytext=start, textcoords='data',
        arrowprops=dict(facecolor='red', width=2.0))
ax.set_aspect('equal')
plt.show()

【讨论】：

常数1.618是什么意思？它来自哪里？

@joaquin：大约是golden ratio。当然，您可以选择任何您喜欢的常量，但它是often looks good。

@unutbu：这两个向量不是正交的，这里肯定有问题。

谢谢，这很有帮助。我想知道为什么 pca.Wt 有这些奇怪的值。我仍然对 pca 类不能真正用于如此简单的 pca 任务感到惊讶。 matplotlib.org/api/mlab_api.html#matplotlib.mlab.PCA 的文档恕我直言无法使用。这里有一个更好的：clear.rice.edu/comp130/12spring/pca/pca_docs.shtml，但我仍然无法真正理解发生了什么。

两个向量长度相同是巧合吗？因为它们也在我的情节中。特征值说明它们应该有多长，对吧？

【参考方案2】：

请注意，matplotlib.mlab.PCA 是 removed in 3.1。

以下是三种可选的 PCA 实现，一种基于 lastmatplotlib.mlab.PCAimplementation，一种基于unutbu's answer，一种基于doug's answer to another question。

前两种使用奇异值分解（svd）得到特征值和特征向量，后者使用协方差矩阵（cov）的方法。

here 可以找到对 svd 和 cov 方法之间关系的精彩解释。

为了便于比较，对实现进行了简化和重构。

def pca_svd(data):
    """ based on matplotlib.mlab.PCA with standardize=False """
    data -= data.mean(axis=0)
    __, singular_values, eigenvectors_transposed = numpy.linalg.svd(
        data, full_matrices=False)
    eigenvalues = singular_values ** 2 / (data.shape[0] - 1)
    eigenvectors = eigenvectors_transposed.T
    transformed_data = numpy.dot(data, eigenvectors)
    return transformed_data, eigenvalues, eigenvectors


def pca_svd_transposed(data):
    """ based on unutbu's answer """
    data -= data.mean(axis=0)
    eigenvectors, singular_values, __ = numpy.linalg.svd(
        data.T, full_matrices=False)  # note data transposed
    eigenvalues = singular_values ** 2 / (data.shape[0] - 1)
    transformed_data = numpy.dot(data, eigenvectors)
    return transformed_data, eigenvalues, eigenvectors
    
    
def pca_cov(data):
    """ based on doug's answer """
    data -= data.mean(axis=0)
    covariance_matrix = numpy.cov(data, rowvar=False)
    eigenvalues, eigenvectors = scipy.linalg.eigh(covariance_matrix)
    decreasing_order = numpy.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[decreasing_order]
    eigenvectors = eigenvectors[:, decreasing_order]
    transformed_data = numpy.dot(data, eigenvectors)
    return transformed_data, eigenvalues, eigenvectors

eigenvalues表示数据沿主轴的方差，即transformed_data的方差。

使用timeit 计时会在我的系统上显示以下内容：

array shape:  (15000, 4)
iterations:  1000
pca_svd_transposed: 4.32 s (average 4.32 ms)
pca_svd:            1.87 s (average 1.87 ms)
pca_cov:            1.41 s (average 1.41 ms)

请注意，转置输入数组的svd 对于这种数组形状来说相对较慢。