案例代码实现
1.导入必要库
首先导入Python中用于数据处理、矩阵运算和可视化的库:
numpy:用于生成模拟数据和矩阵运算;
pandas:用于数据格式化(可选,方便查看);
matplotlib:用于绘制碎石图(选择主成分数量);
sklearn.preprocessing.StandardScaler:用于数据标准化(必选,消除量纲影响)。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
2.模拟人脸像素数据
根据案例背景,模拟5个样本(5张人脸)、100个特征(10×10像素亮度)的数据。
每个样本的特征值是递减的等差数列(模拟额头→下巴的亮度变化,比如额头亮、下巴暗);
加入少量随机噪声(模拟真实图像的细微差异),确保相邻特征相关性高(符合案例中的“冗余信息多”特点)。
def generate_face_data(n_samples=5, n_features=100, noise=0.1):
"""
模拟人脸像素数据(亮度值)
参数:
n_samples:样本数量(人脸数量)
n_features:特征数量(像素点数量,10×10=100)
noise:噪声强度(控制数据的细微差异)
返回:
data:模拟的人脸数据,形状为(n_samples, n_features)
feature_names:特征名称(特征1~特征100)
"""
# 生成基础亮度序列(递减,模拟额头→下巴的亮度变化)
base = np.linspace(200, 50, n_features) # 200(亮)→50(暗)
# 为每个样本生成带噪声的亮度值
data = []
for i in range(n_samples):
# 每个样本的基础序列加上随机噪声(噪声强度为noise×base)
sample = base + noise * base * np.random.randn(n_features)
data.append(sample)
# 转换为numpy数组(形状:n_samples×n_features)
data = np.array(data)
# 生成特征名称(特征1~特征100)
feature_names = [f"特征{i+1}" for i in range(n_features)]
return data, feature_names
# 生成模拟数据(5个样本,100个特征)
X, feature_names = generate_face_data(n_samples=5, n_features=100, noise=0.1)
# 转换为DataFrame(方便查看)
df = pd.DataFrame(X, columns=feature_names, index=[f"P{i+1}" for i in range(5)])
print("模拟人脸数据(前5个特征):")
print(df.iloc[:, :5].round(2)) # 只显示前5个特征,保留2位小数
3.数据标准化(必做步骤)
PCA对量纲敏感,必须将数据标准化(均值为0,方差为1),否则量纲大的特征会主导主成分。使用sklearn的StandardScaler实现标准化:
def standardize_data(data):
"""
数据标准化(均值为0,方差为1)
参数:
data:原始数据,形状为(n_samples, n_features)
返回:
data_std:标准化后的数据
scaler:标准化器(保存均值和标准差,可用于后续逆转换)
"""
scaler = StandardScaler()
data_std = scaler.fit_transform(data)
return data_std, scaler
# 标准化数据
X_std, scaler = standardize_data(X)
print("\n标准化后的数据(前5个特征,均值≈0,方差≈1):")
print(pd.DataFrame(X_std, columns=feature_names).iloc[:, :5].round(2))
4.计算协方差矩阵
协方差矩阵反映特征之间的相关性(标准化后,协方差矩阵等于相关矩阵)。使用numpy的cov函数计算,rowvar=False表示每列是一个特征:
# 计算协方差矩阵(形状:n_features×n_features)
cov_matrix = np.cov(X_std, rowvar=False)
print("\n协方差矩阵形状:", cov_matrix.shape) # 输出:(100, 100)
5.特征值分解(核心步骤)
对协方差矩阵进行特征值分解,得到特征值(主成分的方差)和特征向量(主成分的方向)。使用numpy的linalg.eig函数实现:
# 特征值分解(eigvals:特征值,eigvecs:特征向量)
eigvals, eigvecs = np.linalg.eig(cov_matrix)
# 注意:特征值可能为复数(因计算误差),取实部
eigvals = np.real(eigvals)
eigvecs = np.real(eigvecs)
# 按特征值从大到小排序(主成分按信息多少排列)
sorted_idx = np.argsort(eigvals)[::-1] # 逆序排列的索引
sorted_eigvals = eigvals[sorted_idx] # 排序后的特征值(从大到小)
sorted_eigvecs = eigvecs[:, sorted_idx] # 排序后的特征向量(对应特征值的顺序)
print("\n前5个特征值(从大到小):", sorted_eigvals[:5].round(2))
6.选择主成分数量(关键:累计解释方差比)
通过累计解释方差比选择主成分数量,通常要求累计解释方差比≥80%~95%(保留大部分信息)。同时,绘制碎石图(特征值折线图),找“拐点”(特征值突然下降的位置)辅助判断。
def select_components(eigvals, threshold=0.95):
"""
根据累计解释方差比选择主成分数量
参数:
eigvals:排序后的特征值(从大到小)
threshold:累计解释方差比阈值(默认0.95,即保留95%以上信息)
返回:
n_components:选择的主成分数量
cumulative_er:累计解释方差比(数组)
explained_var_ratio:单个主成分解释方差比(数组)
"""
# 计算总方差(特征值之和)
total_var = np.sum(eigvals)
# 计算每个特征值的解释方差比(单个主成分的信息占比)
explained_var_ratio = eigvals / total_var
# 计算累计解释方差比
cumulative_er = np.cumsum(explained_var_ratio)
# 找到满足累计解释方差比≥threshold的最小主成分数量
n_components = np.argmax(cumulative_er >= threshold) + 1 # +1因为索引从0开始
return n_components, cumulative_er, explained_var_ratio
# 选择主成分(累计解释方差比≥95%)
n_components, cumulative_er, explained_var_ratio = select_components(sorted_eigvals, threshold=0.95)
print(f"\n选择的主成分数量:{n_components}")
print(f"前{n_components}个主成分的累计解释方差比:{cumulative_er[n_components-1]:.2%}")
# 绘制碎石图(特征值折线图)
plt.figure(figsize=(8, 4))
plt.plot(range(1, len(sorted_eigvals)+1), sorted_eigvals, marker='o', linestyle='--')
plt.title('碎石图(Scree Plot)')
plt.xlabel('主成分编号')
plt.ylabel('特征值(方差)')
plt.axvline(x=n_components, color='r', linestyle='-', label=f'选择的主成分数量:{n_components}')
plt.legend()
plt.show()
7.投影降维(得到主成分)
将标准化后的数据投影到前n_components个特征向量上,得到降维后的主成分矩阵。
def project_data(data_std, eigvecs, n_components):
"""
将数据投影到前n_components个特征向量上(降维)
参数:
data_std:标准化后的数据(n_samples×n_features)
eigvecs:排序后的特征向量(n_features×n_features)
n_components:选择的主成分数量
返回:
principal_components:主成分矩阵(n_samples×n_components)
"""
# 取前n_components个特征向量(主成分方向)
top_eigvecs = eigvecs[:, :n_components]
# 投影计算(矩阵乘法:数据×特征向量)
principal_components = data_std @ top_eigvecs
return principal_components
# 投影降维,得到主成分矩阵
pc_matrix = project_data(X_std, sorted_eigvecs, n_components)
# 转换为主成分DataFrame(方便查看)
pc_df = pd.DataFrame(pc_matrix, columns=[f"PC{i+1}" for i in range(n_components)], index=[f"P{i+1}" for i in range(5)])
print("\n降维后的主成分矩阵(5个样本×{n_components}个主成分):".format(n_components=n_components))
print(pc_df.round(2))
8.主成分解释(结合业务)
主成分是原始特征的线性组合,权重为特征向量的分量。通过分析特征向量的权重,可以解释主成分的含义。
# 提取第一个主成分的特征向量(权重)
pc1_weights = sorted_eigvecs[:, 0] # 第一个主成分的方向(特征向量)
# 转换为DataFrame(方便查看)
pc1_weights_df = pd.DataFrame(pc1_weights, index=feature_names, columns=["权重"])
# 按权重绝对值从大到小排序(查看对PC1贡献大的特征)
pc1_weights_df = pc1_weights_df.sort_values(by="权重", ascending=False)
print("\n第一个主成分(PC1)的特征权重(前10个特征):")
print(pc1_weights_df.head(10).round(3))
print("\nPC1的含义:所有特征的权重均为正,且数值较大(比如前10个特征的权重都在0.09以上),说明PC1是**整体亮度的综合**(比如额头、眼睛、嘴巴等部位的亮度总和),能有效区分不同人的整体亮度差异(比如P3的PC1值最大,说明他的脸最亮;P4的PC1值最小,说明他的脸最暗)。")
9.完整主程序(整合所有步骤)
将上述步骤整合为一个主程序,方便运行:
if __name__ == "__main__":
# 1. 模拟数据
X, feature_names = generate_face_data(n_samples=5, n_features=100, noise=0.1)
# 转换为DataFrame(方便查看)
df = pd.DataFrame(X, columns=feature_names, index=[f"P{i+1}" for i in range(5)])
print("模拟人脸数据(前5个特征):")
print(df.iloc[:, :5].round(2))
# 2. 标准化数据
X_std, scaler = standardize_data(X)
print("\n标准化后的数据(前5个特征):")
print(pd.DataFrame(X_std, columns=feature_names).iloc[:, :5].round(2))
# 3. 计算协方差矩阵
cov_matrix = np.cov(X_std, rowvar=False)
print("\n协方差矩阵形状:", cov_matrix.shape)
# 4. 特征值分解
eigvals, eigvecs = np.linalg.eig(cov_matrix)
eigvals = np.real(eigvals)
eigvecs = np.real(eigvecs)
sorted_idx = np.argsort(eigvals)[::-1]
sorted_eigvals = eigvals[sorted_idx]
sorted_eigvecs = eigvecs[:, sorted_idx]
print("\n前5个特征值:", sorted_eigvals[:5].round(2))
# 5. 选择主成分数量
n_components, cumulative_er, explained_var_ratio = select_components(sorted_eigvals, threshold=0.95)
print(f"\n选择的主成分数量:{n_components}")
print(f"累计解释方差比:{cumulative_er[n_components-1]:.2%}")
# 6. 绘制碎石图
plt.figure(figsize=(8, 4))
plt.plot(range(1, len(sorted_eigvals)+1), sorted_eigvals, marker='o', linestyle='--')
plt.title('碎石图(Scree Plot)')
plt.xlabel('主成分编号')
plt.ylabel('特征值(方差)')
plt.axvline(x=n_components, color='r', linestyle='-', label=f'选择的主成分数量:{n_components}')
plt.legend()
plt.show()
# 7. 投影降维
pc_matrix = project_data(X_std, sorted_eigvecs, n_components)
pc_df = pd.DataFrame(pc_matrix, columns=[f"PC{i+1}" for i in range(n_components)], index=[f"P{i+1}" for i in range(5)])
print("\n降维后的主成分矩阵:")
print(pc_df.round(2))
# 8. 解释主成分
pc1_weights = sorted_eigvecs[:, 0]
pc1_weights_df = pd.DataFrame(pc1_weights, index=feature_names, columns=["权重"]).sort_values(by="权重", ascending=False)
print("\n第一个主成分(PC1)的特征权重(前10个特征):")
print(pc1_weights_df.head(10).round(3))
print("\nPC1的含义:整体亮度的综合,能区分不同人的整体亮度差异。")
10.如何使用?
运行环境:确保安装了numpy、pandas、matplotlib、sklearn库(可通过pip install numpy pandas matplotlib scikit-learn安装)。
修改参数:
若要模拟更多样本或特征,可调整generate_face_data函数的n_samples(样本数量)和n_features(特征数量)参数;
若要改变累计解释方差比阈值(比如保留80%信息),可调整select_components函数的threshold参数(如threshold=0.8)。
结果解读:
模拟数据:前5个特征是递减的,符合额头→下巴的亮度变化;
标准化后的数据:均值≈0,方差≈1,消除了量纲影响;
特征值:前几个特征值很大,后面的很小(碎石图有明显拐点);
主成分矩阵:降维后的结果,每个样本的主成分值代表其“综合特征”(比如PC1代表整体亮度);
主成分解释:通过特征向量的权重,可解释主成分的业务含义(比如PC1是整体亮度的综合)。
11.注意事项
必须标准化:否则量纲大的特征会主导主成分(比如若有一个特征是“身高(cm)”,另一个是“体重(kg)”,未标准化会导致体重的方差远大于身高,PC1会主要反映体重的变化);
主成分顺序:必须按特征值从大到小排列,否则会丢失最重要的信息;
主成分解释:需结合业务背景(比如人脸数据的主成分通常代表“整体亮度”“五官位置”等);
局限性:PCA是线性降维,无法捕捉非线性关系(比如“U”型数据);若需保留与标签相关的信息,可使用有监督降维(如LDA)。
通过这个案例,数模小白可以掌握PCA的完整流程,包括数据模拟、标准化、特征值分解、主成分选择、投影降维及结果解释,能快速将PCA应用到实际问题(如高维数据降维、可视化、去除冗余)中。