【数学建模】常用算法-主成分分析PCA的Python实现

1前言

本文主要讲解主成分分析析法（PCA）的python实现，后续会跟进实例分析

2 原理-代码实现

2.1 实现步骤

主成分分析PCA是一种应用广泛的和降维方法，对其实现做以下归纳

2.2 代码实现

导入包

import numpy as np

• 定义计算协方差矩阵函数
  X为输入的数据，m为样本数据的条数，也就是X的行数。
  对X进行标准化，方法为：减去均值除以方差，这部分的原理不懂的可以百度一下。
  标准化之后的数据就是均值为0，方差为1的标准正态分布。

# 计算协方差矩阵
def calc_cov(X):
    m = X.shape[0] # 样本的数量，行数
    # 数据标准化
    X = (X - np.mean(X, axis=0)) / np.var(X, axis=0) # 标准化之后均值为0，方差为1
    return 1 / m * np.matmul(X.T, X) # matmul为两个矩阵的乘积

• 定义PCA的流程
  首先计算输入数据X的协方差，然后计算其特征值记为：eigenvalues，计算其特征向量记为：eigenvectors
  计算特征值和特征向量用的是 np.linalg.eig()函数，使用起来十分方便
  然后接下来就是计算出矩阵P，用Y=XP计算出降维后的数据Y

def pca(X, n_components):
    # 计算协方差矩阵
    cov_matrix = calc_cov(X)
    # 计算协方差矩阵的特征值和对应特征向量
    eigenvalues, eigenvectors = np.linalg.eig(cov_matrix) # eigenvalues特征值，eigenvectors特征向量
    # 对特征值排序
    idx = eigenvalues.argsort()[::-1]
    # 取最大的前n_component组
    eigenvectors = eigenvectors[:, idx]
    eigenvectors = eigenvectors[:, :n_components]
    # Y=XP转换
    return np.matmul(X, eigenvectors)

2.3鸢尾花数据集例子

导入数据

from sklearn import datasets
import matplotlib.pyplot as plt
# 导入鸢尾花数据集
iris = datasets.load_iris()
X = iris.data
y = iris.target

查看数据的形状，其结果为(150, 4)

X.shape
# (150, 4)

计算协方差矩阵

cov_matrix = calc_cov(X) # 计算特征值
cov_matrix

可以看到协方差矩阵为4*4的矩阵，然后我们计算该矩阵的特征值和特征向量

eigenvalues, eigenvectors = np.linalg.eig(cov_matrix) # eigenvalues特征值，eigenvectors特征向量

然后计算我们需要的P，这里我们保留3个主成分

idx = eigenvalues.argsort()[::-1]
# 取最大的前n_component组
eigenvectors = eigenvectors[:, idx]
eigenvectors = eigenvectors[:, :3]

得到了一个4行3列的矩阵

然后利用P求得降维后的数据

# Y=PX转换
np.matmul(X, eigenvectors)

降维后的数据为（150, 4）*（4, 3）=（150, 3）
也就是150条，3列的数据，数据由原来的4维降低到了3维

3 基于Sklearn的实现

# 导入sklearn降维模块
from sklearn import decomposition
# 创建pca模型实例，主成分个数为3个
pca = decomposition.PCA(n_components=3) # 写我们需要几个主成分
# 模型拟合
pca.fit(X)
# 拟合模型并将模型应用于数据X
X_trans = pca.transform(X)

# 颜色列表
colors = ['navy', 'turquoise', 'darkorange']
# 绘制不同类别
for c, i, target_name in zip(colors, [0,1,2], iris.target_names):
    plt.scatter(X_trans[y == i, 0], X_trans[y == i, 1], 
            color=c, lw=2, label=target_name)
# 添加图例
plt.legend()
plt.show()