最近在学习李航老师的《统计学习方法》第16章-主成分分析（PCA），结合运用python代码实现，记录一下。

1.基础知识

1.1 基本想法

1.2 主要性质及定理

1.2.1 主成分的表示方法

1.2.2 主成分的总体性质

1.2.3 什么是因子负荷量

1.2.4 主成分的方差贡献率

1.3 规范化变量的总体主成分与样本主成分步骤

总体主成分分析是定义在样本总体上的。在实际问题中，需要在观测数据上进行主成分分析，这就是样本主成分分析。

假设对m维随机变量【既有m个特征变量】进行n次独立观测，下面的m、n 同假设。

1.4 举例理解主成分分析

2. python代码实现PCA

2.1 手动实现PCA

参考的一些优秀的文章
如何通俗理解PCA（主成分分析）算法的数学原理和代码实现
Python机器学习——主成分分析PCA
审稿人：PCA的误区就是"分类"，但Python可以画得很漂亮！

2.1.0 导入所需库包和数据集

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler 
from matplotlib.patches import Ellipse
import seaborn as sns
plt.rcParams['font.sans-serif'] = ['SimHei']#解决图表中中文显示问题
plt.rcParams['axes.unicode_minus']=False# 解决图表中负号不显示问题
import warnings
warnings.filterwarnings("ignore")

# 设置Seaborn风格
sns.set(style="whitegrid", palette="muted", font_scale=1.2)

# 加载Iris数据集
data = load_iris()
X = data.data  # 特征矩阵
y = data.target  # 标签
target_names = data.target_names

2.1.1 数据规范化处理

from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.3,stratify=y,random_state=0)

# 标准化
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)

2.1.2 计算相关矩阵

主成分分析从协方差出发和从相关矩阵出发的区别

由于这里已经把数据做了标准化，所以下面计算一下相关矩阵，

cov_mat = np.corrcoef(X_train_std.T)   # 计算相关矩阵
# cov_mat = np.cov(X_train_std.T)   # 计算协方差矩阵
cov_mat

2.1.3 求相关矩阵的特征值及其对应的单位特征向量

#计算相关矩阵的特征值和特征向量
eigen_vals,eigen_vecs = np.linalg.eig(cov_mat)
eigen_vals_sorted =eigen_vals [np.argsort(-eigen_vals)]  # 把特征值从大到小排序
print('特征值从大到小排序',eigen_vals_sorted,'\n')

eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]
# 特征对按特征值降序排序
eigen_pairs.sort(key=lambda k: k[0], reverse=True)  #特征值与其对应的单位特征向量
print('特征值与其对应的单位特征向量\n',eigen_pairs)

2.1.4 主成分的累计方差贡献率

计算出相关矩阵的特征值，按大小排序，根据特征值即可计算出主成分的累计方差贡献率。

# 累加解释方差（特征值）之和
def cum_var(eigen_vals):
    k = 0
    cum_var_exp = 0
    tot = sum(eigen_vals)
    var_exp_list = []
    cum_var_exp_list = []
    for i in sorted(eigen_vals, reverse=True):
        k = k+1
        var_exp = i / tot
        var_exp_list.append(var_exp)
        print(f'第{k}个主成分解释的方差比例:',var_exp)
        cum_var_exp = cum_var_exp + np.cumsum(var_exp)
        cum_var_exp_list.append(cum_var_exp)
        print(f'前{k}个主成分累计解释方差之和:',cum_var_exp[0])
    return var_exp_list,cum_var_exp_list
var_exp_list,cum_var_exp_list = cum_var(eigen_vals)

# 绘制解释方差
plt.bar(range(1,5), var_exp_list, alpha=0.5,
        align='center', label='individual explained variance')
plt.step(range(1,5), cum_var_exp_list, where='mid',
         label='cumulative explained variance')
plt.ylabel('Explained variance ratio')  #特征值的解释方差比累计和
plt.xlabel('Principal component index')   # 主成分顺序
plt.title('PCA-cumulative explained variance') 
plt.legend(loc='best')
plt.show()

2.1.5 主成分的分析结果

var_ma = pd.DataFrame(eigen_vecs.T, columns=data.feature_names,index=['第一主成分-y1','第二主成分-y2','第三主成分-y3','第四主成分-y4'])
var_vals =  pd.DataFrame(var_exp_list, columns=['方差贡献率'],index=['第一主成分-y1','第二主成分-y2','第三主成分-y3','第四主成分-y4'])
cum_var =  pd.DataFrame(cum_var_exp_list, columns=['累计方差贡献率'],index=['第一主成分-y1','第二主成分-y2','第三主成分-y3','第四主成分-y4'])

#单位特征向量和主成分的方差贡献率
Ma = pd.concat([var_ma,var_vals],axis=1)
Mat = pd.concat([Ma,cum_var],axis=1)
Mat

其中，x1、x2、x3、x4分别表示sepal length (cm)、sepalwidth (cm)、petal length (cm)、petal width (cm)。

2.1.6 主成分的因子负荷量

# 使用累计方差贡献率大于90%的标准确定主成分的个数
for m in range(1,len(eigen_vals_sorted)):
    if eigen_vals_sorted[:m].sum()/eigen_vals_sorted.sum() >= 0.9:
        print('特征值的比重大于90%时，','m={}'.format(m))
        break
        

a1 = np.sqrt(eigen_vals[0])*eigen_vecs[:,0]
a2 = np.sqrt(eigen_vals[1])*eigen_vecs[:,1]
A_pcm = np.vstack((a1,a2)) 
A_pcm = np.mat(A_pcm)
print('因子载荷矩阵：\n')
# 计算主成分的因子负荷量
PA_pcm = pd.DataFrame(A_pcm,columns=data.feature_names)
PA_pcm

因子负荷量的分布图

plt.scatter(PA_pcm.loc[0,:], PA_pcm.loc[1,:]) 
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.show() 

2.1.7 将数据投影到新特征上

# 从“前” 2个特征向量中构建投影矩阵 
w = np.hstack((eigen_pairs[0][1][:, np.newaxis], eigen_pairs[1][1][:, np.newaxis]))
# 前两个特征值的特征向量，行拼接为投影矩阵
print('Matrix W:\n', w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']

X_train_pca = X_train_std.dot(w)
for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train==l, 0], 
                X_train_pca[y_train==l, 1], 
                c=c, label=l, marker=m) 
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.show() 

2.2 scikit-learn实现PCA

2.2.1 StandardScaler标准化处理

# 加载Iris数据集
data = load_iris()
X = data.data  # 特征矩阵
y = data.target  # 标签
target_names = data.target_names

from sklearn.preprocessing import StandardScaler  
scaler = StandardScaler() 
X_scaled = scaler.fit_transform(X)

2.2.2 利用sklearn.decomposition的 PCA降维

# PCA降维
pca = PCA(n_components=4)  # n_components=3表示降到3维，也可以 n_components=0.9 表示需要保留多少信息确定需要降到多少维
X_pca = pca.fit_transform(X_scaled)
print('维数：',pca.n_components_)

# print('保留原数据90%的信息需要的维数：',pca.n_components_)

# 各个特征的可解释性方差
eigenvalues = pca.explained_variance_
eigenvecs = pca.components_.T
print(f'\n特征值: ',eigenvalues)
#单位特征向量
print(f'特征值对应的单位特征向量:\n ',eigenvecs)

explained_variance_ratio = pca.explained_variance_ratio_
print(f'\n各个特征的可解释性方差: ',explained_variance_ratio)
print(f'PC1解释的方差比例: {explained_variance_ratio[0]:.4f}')
print(f'PC2解释的方差比例: {explained_variance_ratio[1]:.4f}')
print(f'PC3解释的方差比例: {explained_variance_ratio[2]:.4f}')
print(f'PC4解释的方差比例: {explained_variance_ratio[3]:.4f}')

2.2.3 原始数据在新特征上的展现

# 计算置信椭圆
def plot_confidence_ellipse(ax, x, y, n_std=2.0, facecolor='none', edgecolor='black', alpha=0.3, **kwargs):
    """
    在给定的轴上绘制置信椭圆
    """
    if x.size != y.size:
        raise ValueError("x 和 y 的尺寸必须相同")
    
    cov = np.cov(x, y)  # 计算协方差矩阵
    mean_x, mean_y = np.mean(x), np.mean(y)  # 计算均值
    
    # 计算椭圆的主轴长度和角度
    lambda_, v = np.linalg.eig(cov)
    lambda_ = np.sqrt(lambda_)
    angle = np.degrees(np.arctan2(*v[:, 0][::-1]))
    
    # 绘制椭圆
    ellipse = Ellipse((mean_x, mean_y), width=lambda_[0]*n_std*2, height=lambda_[1]*n_std*2,
                      angle=angle, facecolor=facecolor, edgecolor=edgecolor, alpha=alpha, linewidth=2, **kwargs)
    ax.add_patch(ellipse)

# 创建图形
fig, ax = plt.subplots(figsize=(10, 7))

# 调色板
colors = sns.color_palette("husl", n_colors=3)

# 绘制数据点和置信椭圆
for i, target_name in enumerate(target_names):
    x_pca = X_pca[y == i, 0]
    y_pca = X_pca[y == i, 1]  #从4降维到 2维，任选两个画椭圆图
    ax.scatter(x_pca, y_pca, s=100, label=target_name, color=colors[i], edgecolor='k', linewidth=0.8, alpha=0.7)
    plot_confidence_ellipse(ax, x_pca, y_pca, n_std=2.0, edgecolor=colors[i], alpha=0.2)

# 美化坐标轴
ax.set_xlabel('PCA Component 1', fontsize=14, fontweight='bold')
ax.set_ylabel('PCA Component 2', fontsize=14, fontweight='bold')
ax.set_title('PCA of IRIS Dataset with 95% Confidence Ellipses', fontsize=16, fontweight='bold', pad=15)
# 设置图例位置在坐标系外部
ax.legend(
    title='Species', 
    loc='upper left', 
    bbox_to_anchor=(1.05, 1),  # 将图例放在右侧外部
    fontsize=12, 
    title_fontsize=13
)

# 美化网格与边框
ax.grid(True, linestyle='--', alpha=0.5)
sns.despine(trim=True, offset=10)
# 显示图形
plt.tight_layout()
plt.show()

关于因子负荷量和投影矩阵的计算参考2.1.6和2.1.7。