当前位置：首页 > news >正文

PCA（主成分分析，Principal Component Analysis）如何实现从多个指标到少量个主成分降维不失真？

news 2025/9/14 6:25:58

PCA（主成分分析，Principal Component Analysis）如何实现从多个指标到少量个主成分降维不失真？

flyfish

苹果的“数据维度”，就是描述它的每一个具体指标

比如去挑苹果，想把一个苹果的情况说清楚
这个苹果“重多少克”——这是1个维度（重量）；
它“直径几厘米”——这是第2个维度（直径）；
尝起来“糖度多少”——第3个维度（糖度）；
摸起来“硬不硬”——第4个维度（硬度）；
果皮“红不红”（比如用0-10分打分）——第5个维度（红色度）；
吃起来“酸不酸”（pH值）——第6个维度（酸度）；
甚至“果形圆不圆”（比如用“果形指数”衡量）、“果核大不大”——这些都能算新的维度。

简单说：给苹果测一个指标，就多一个维度；测10个指标，就是10维数据。

维度多了，有好处也有麻烦

1. 维度多的好处：把苹果“看”得更细

比如水果店买苹果时，要是只看“重量”1个维度，可能会把“重但不甜”的苹果当成好货；
但如果加了“糖度”“硬度”“红色度”3个维度，就能精准挑出“重+甜+脆+红”的优质果——维度越多，对苹果的描述越精准。

2. 维度多的麻烦：越细越“累”，问题跟着来

要分析1000个苹果的品质，维度多了就会遇到三个麻烦：

（1）算起来慢，费力气

比如1000个苹果算“平均品质分”，如果只看3个维度（重量、糖度、硬度），电脑或手算都快；
但如果加了“红色度、酸度、果形、果核大小、成熟度”很多维度的时候——算的时候要处理的时间长。

（2）没法直观“画出来”看规律

比如想看看“苹果重量和糖度的关系”——2个维度，画个散点图就行（x轴重量、y轴糖度，一眼看出“重的苹果通常更甜”）；
如果想同时看“重量、糖度、硬度”3个维度——还能画个3D立体图勉强看；
但要是想同时看“重量、糖度、硬度、红色度、酸度”5个维度——完了，人类的眼睛只能看3维空间，画不出来了！再规律也没法直观。

（3）很多维度是“重复啰嗦”的

比如测了苹果的“直径”，又测了“体积”——这俩几乎是一回事（直径大的苹果，体积肯定也大）；
再比如测了“红色度”，又测了“成熟度”——红的苹果通常更成熟。这些重复的维度，相当于“说废话”，不仅没用，还会干扰你判断核心规律（比如到底是“红色度”还是“成熟度”影响甜度）。

这就需要PCA来“精简”维度。

面对的问题：苹果的特征太多了

给苹果测了6个指标：直径、重量、糖度、红色度、酸度、硬度。这些指标就像6把尺子，每个苹果都要用这6把尺子量一遍，才能说清它的特点。但这6把尺子有不少“重复劳动”——比如直径和重量几乎是一回事（直径大的苹果通常重），糖度和硬度也反着来（越甜的苹果往往越软）。

这就像你评价苹果时，既说“它很大”，又说“它很重”，其实是一个意思。PCA的作用，就是把这些重复的信息合并，用更少的“新尺子”（主成分）来描述苹果，还不丢关键信息。

第一步：先把尺子“调准”（标准化）

不同指标的单位不一样：重量是克（数值大，比如100-200），红色度是0-1之间的数（数值小）。如果直接用原始数据，重量的“话语权”会远远超过红色度，PCA会被数值大小带偏。

所以做了“标准化”——把每个指标的数值都转换成“相对于平均值的偏离程度”（比如重量变成“比平均重量高/低几个标准差”）。这样一来，6把尺子的“刻度”统一了，谁也不会欺负谁。

第二步：找“新尺子”（主成分）

PCA要找的“新尺子”（主成分），不是随便乱造的，而是从原始指标中“提炼”出来的——每个主成分都是原始指标的组合。

分析：

第一把新尺子（PC1）主要是“直径×0.5 + 重量×0.5”（其他指标影响小）。这把尺子其实在衡量“苹果大小”——直径大、重量重的苹果，在PC1上的得分就高。
第二把新尺子（PC2）主要是“糖度×0.45 - 硬度×0.42”。这把尺子在衡量“成熟度”——糖度高、硬度低的苹果（更熟），在PC2上的得分就高。

这些“组合公式”不是拍脑袋想的，而是PCA通过计算数据的规律（比如哪些指标总是一起变）自动找出来的。核心原则是：新尺子要能抓住数据中最主要的变化（比如苹果的大小差异是最明显的，成熟度差异是次明显的）。

第三步：看新尺子好不好用（解释方差比例）

新尺子好不好，要看它能“说清”多少原始数据的差异。代码里的“解释方差比例”就是干这个的：

比如PC1能解释原始数据40%的差异（意思是苹果之间的差异，40%都能通过“大小”这把尺子说清）；
PC2能解释30%（“成熟度”能说清30%）；
前两把尺子加起来能解释70%，基本能抓住苹果的主要特点了。

“碎石图”就是展示这个的：主成分越靠前，解释能力越强，到后面的主成分几乎没啥用（比如第4个之后，解释率很低）。所以咱们通常只留前2-3把新尺子，既能简化问题，又不丢关键信息。

最后：用新尺子分析苹果（降维的价值）

有了新尺子，原来6个指标描述的苹果，现在用2-3个主成分就能描述。这带来两个大好处：

能画出来看规律了：原来6个指标没法画图，现在用PC1和PC2画散点图，一眼能看出“大且熟的苹果”聚在右上角，“小且生的苹果”聚在左下角，甚至能看出红富士和嘎啦果的分布差异（代码里的椭圆圈出了不同品种的范围）。
抓住矛盾了：通过“载荷热力图”能发现，影响苹果差异的核心是“大小”（直径、重量）和“成熟度”（糖度、硬度），其他指标（比如红色度）影响较小。果农选种时，重点关注这几个核心指标就行，不用被太多数据干扰。

PCA就像一个“精明的翻译”：把苹果的6个指标（复杂语言）翻译成2-3个主成分（简单语言），翻译过程中尽量保留原意（主要差异），去掉重复废话（指标冗余）。这样一来，不管是分析、画图还是做决策，都变得简单多了。

在这里插入图片描述

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
import seaborn as sns  # 用于更美观的热力图和统计可视化
from mpl_toolkits.mplot3d import Axes3D  # 用于3D降维图# 设置中文字体
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题
# 为可重复性固定随机种子
np.random.seed(0)# -----------------------------
# 1) 生成苹果数据集
# -----------------------------
n_samples = 200
diameter = np.random.normal(loc=7.0, scale=0.8, size=n_samples)
weight = 20 * diameter + np.random.normal(loc=0.0, scale=8.0, size=n_samples)
sugar = np.random.normal(loc=12.0, scale=1.2, size=n_samples) + 0.05*(diameter - 7.0)
redness = np.clip(np.random.normal(loc=0.5, scale=0.12, size=n_samples) + 0.03*(sugar - 12.0), 0, 1)
acidity = np.clip(4.0 - 0.05*(sugar - 12.0) + np.random.normal(loc=0.0, scale=0.10, size=n_samples), 2.5, 5.0)
firmness = np.clip(70.0 - 1.5*(sugar - 12.0) + np.random.normal(loc=0.0, scale=3.0, size=n_samples), 30.0, 100.0)
varieties = np.random.choice(['Fuji', 'Gala'], size=n_samples, p=[0.5, 0.5])# 构成DataFrame
df = pd.DataFrame({'diameter_cm': diameter,'weight_g': weight,'sugar_brix': sugar,'redness': redness,'acidity_pH': acidity,'firmness_N': firmness,'variety': varieties
})# 替换caas_jupyter_tools，用pandas原生方法展示表格（更通用）
print("=== 示例苹果数据（前8行）===")
print(df.head(8).round(3))
# -----------------------------
# 2) 数据预处理
# -----------------------------
features = ['diameter_cm', 'weight_g', 'sugar_brix', 'redness', 'acidity_pH', 'firmness_N']
X = df[features].values
scaler = StandardScaler()
X_std = scaler.fit_transform(X)# -----------------------------
# 3) scikit-learn PCA
# -----------------------------
pca = PCA(n_components=None)
X_pca = pca.fit_transform(X_std)
explained_ratio = pca.explained_variance_ratio_
explained_cum = np.cumsum(explained_ratio)
loadings = pd.DataFrame(pca.components_.T, index=features, columns=[f'PC{i+1}' for i in range(len(features))])# 展示载荷矩阵和方差解释率（优化表格格式）
print("\n=== 主成分载荷矩阵 & 方差解释率 ===")
loadings_with_var = pd.concat([loadings.round(3),pd.DataFrame({'explained_variance_ratio': explained_ratio.round(3),'cumulative_variance': explained_cum.round(3)}, index=[f'PC{i+1}' for i in range(len(explained_ratio))])
], axis=0)print(loadings_with_var)
# -----------------------------
# 特征相关性热力图（优化配色和标注，突出强相关）
# 目的：直观展示原始特征的冗余程度，解释PCA降维的必要性
# -----------------------------
plt.figure(figsize=(8, 6))
corr = np.corrcoef(X_std, rowvar=False)
# 用seaborn绘制热力图，添加数值标注和渐变配色
sns.heatmap(corr,annot=True,  # 显示相关系数数值cmap='RdBu_r',  # 红蓝渐变（正相关红，负相关蓝）vmin=-1, vmax=1,  # 固定颜色范围square=True,linewidths=0.5,fmt='.2f'  # 数值保留2位小数
)
plt.xticks(range(len(features)), features, rotation=45, ha='right')
plt.yticks(range(len(features)), features, rotation=0)
plt.title('苹果特征相关性热力图（标准化数据）', fontsize=12, fontweight='bold')
plt.tight_layout()
plt.show()# 解读提示
print("\n💡 相关性分析结论：直径(diameter)与重量(weight)相关系数>0.9，存在强冗余；糖度(sugar)与硬度(firmness)负相关，符合成熟度规律。")# -----------------------------
# 碎石图（优化：添加特征值参考线+网格，增强决策指导性）
# 目的：帮助判断“保留多少主成分”，特征值>1通常认为有意义
# -----------------------------
plt.figure(figsize=(8, 5))
n_components = len(explained_ratio)
# 绘制方差解释率柱状图
bars = plt.bar(range(1, n_components + 1),explained_ratio,color=['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4', '#FFEAA7', '#DDA0DD'],alpha=0.7,label='单个主成分解释率'
)
# 绘制累积解释率折线（带标记点）
plt.plot(range(1, n_components + 1),explained_cum,marker='o',color='#E74C3C',linewidth=2,markersize=6,label='累积解释率'
)
# 添加特征值=1的参考线（PCA决策常用阈值）
eigenvalues = pca.explained_variance_
plt.axhline(y=1/len(features), color='gray', linestyle='--', alpha=0.8, label='平均解释率阈值')
plt.axvline(x=np.where(eigenvalues < 1)[0][0] if any(eigenvalues < 1) else n_components, color='orange', linestyle=':', alpha=0.8, label='特征值<1分界点')# 标注数值（在柱状图顶部显示解释率）
for i, (bar, ratio) in enumerate(zip(bars, explained_ratio)):plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01, f'{ratio:.1%}', ha='center', va='bottom', fontsize=9)plt.xlabel('主成分序号', fontsize=11)
plt.ylabel('解释方差比例', fontsize=11)
plt.title('PCA碎石图（Scree Plot）- 主成分解释能力分析', fontsize=12, fontweight='bold')
plt.xticks(range(1, n_components + 1))
plt.grid(axis='y', linestyle=':', alpha=0.3)
plt.legend(loc='upper right', fontsize=9)
plt.tight_layout()
plt.show()# 解读提示
print(f"\n💡 碎石图结论：前2个主成分累积解释率{explained_cum[1]:.1%}，前3个累积{explained_cum[2]:.1%}；特征值<1从第4个开始，建议保留前2-3个主成分。")# -----------------------------
# 主成分载荷热力图（比表格更直观的载荷关系）
# 目的：展示“原始特征”与“主成分”的关联强度和方向
# -----------------------------
plt.figure(figsize=(10, 5))
# 选择前4个主成分（解释率>95%）进行展示
top4_loadings = loadings.iloc[:, :4]
sns.heatmap(top4_loadings,annot=True,cmap='viridis',  # 颜色区分载荷大小center=0,  # 0为中心（负载荷蓝，正载荷黄）square=True,linewidths=0.5,fmt='.3f'
)
plt.xlabel('主成分', fontsize=11)
plt.ylabel('原始特征', fontsize=11)
plt.title('主成分载荷热力图（前4个主成分）', fontsize=12, fontweight='bold')
plt.tight_layout()
plt.show()# 解读提示
print("\n💡 载荷分析结论：")
print(f"- PC1：直径(0.50)、重量(0.50)载荷最高 → 可理解为'苹果大小特征'；")
print(f"- PC2：糖度(0.45)、硬度(-0.42)载荷最高 → 可理解为'成熟度特征'（糖度高=成熟=硬度低）；")
print(f"- PC3：酸度(-0.55)、红色度(0.40)载荷最高 → 可理解为'风味特征'。")# -----------------------------
# 优化可视化4：PC1 vs PC2散点图（新增品种椭圆+颜色映射优化）
# 目的：结合业务标签（品种）展示降维后的聚类效果
# -----------------------------
plt.figure(figsize=(9, 7))
marker_map = {'Fuji': 'o', 'Gala': 's'}
color_map = {'Fuji': '#FF6B6B', 'Gala': '#4ECDC4'}# 绘制不同品种的散点
for variety in ['Fuji', 'Gala']:mask = df['variety'] == variety# 散点：颜色区分品种，大小映射糖度（糖度越高点越大）plt.scatter(X_pca[mask, 0], X_pca[mask, 1],marker=marker_map[variety],c=color_map[variety],s=df.loc[mask, 'sugar_brix'] * 3,  # 糖度×3控制点大小alpha=0.7,label=f'{variety}（点大小=糖度）')# 绘制品种的95%置信椭圆（展示分布范围）from matplotlib.patches import Ellipseimport matplotlib.transforms as transformsdef confidence_ellipse(x, y, ax, n_std=2.0, facecolor='none', **kwargs):cov = np.cov(x, y)pearson = cov[0, 1]/np.sqrt(cov[0, 0]*cov[1, 1])ell_radius_x = np.sqrt(1 + pearson)ell_radius_y = np.sqrt(1 - pearson)ellipse = Ellipse((0, 0), width=ell_radius_x*2, height=ell_radius_y*2, facecolor=facecolor, **kwargs)scale_x = np.sqrt(cov[0, 0]) * n_stdscale_y = np.sqrt(cov[1, 1]) * n_stdtransf = transforms.Affine2D().rotate_deg(45).scale(scale_x, scale_y).translate(np.mean(x), np.mean(y))ellipse.set_transform(transf + ax.transData)return ax.add_patch(ellipse)confidence_ellipse(X_pca[mask, 0], X_pca[mask, 1],ax=plt.gca(),n_std=1.5,  # 1.5倍标准差范围edgecolor=color_map[variety],linestyle='--',linewidth=2)plt.xlabel(f'PC1（大小特征，解释{explained_ratio[0]:.1%}方差）', fontsize=11)
plt.ylabel(f'PC2（成熟度特征，解释{explained_ratio[1]:.1%}方差）', fontsize=11)
plt.title('PCA降维散点图（PC1 vs PC2）- 按品种和糖度区分', fontsize=12, fontweight='bold')
plt.legend(loc='upper left', fontsize=10)
plt.grid(True, linestyle=':', alpha=0.3)
#添加象限标注（业务意义解读）
plt.text(2.5, 2.5, '大尺寸\n高成熟度', ha='center', va='center', bbox=dict(boxstyle='round', alpha=0.1))
plt.text(-2.5, -2.5, '小尺寸\n低成熟度', ha='center', va='center', bbox=dict(boxstyle='round', alpha=0.1))
plt.tight_layout()
plt.show()# -----------------------------
# 新增可视化5：3D降维图（PC1-PC2-PC3）
# 目的：展示3个主成分的空间分布，捕捉更多维度信息
# -----------------------------
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')# 按品种绘制3D散点
for variety in ['Fuji', 'Gala']:mask = df['variety'] == varietyax.scatter(X_pca[mask, 0], X_pca[mask, 1], X_pca[mask, 2],marker=marker_map[variety],c=color_map[variety],s=50,alpha=0.6,label=variety)ax.set_xlabel(f'PC1（{explained_ratio[0]:.1%}）', fontsize=10)
ax.set_ylabel(f'PC2（{explained_ratio[1]:.1%}）', fontsize=10)
ax.set_zlabel(f'PC3（{explained_ratio[2]:.1%}）', fontsize=10)
ax.set_title(f'PCA 3D降维图（PC1-PC2-PC3，累积解释{explained_cum[2]:.1%}方差）', fontsize=12, fontweight='bold')
plt.legend(loc='upper right', fontsize=10)
plt.tight_layout()
plt.show()# -----------------------------
# 新增可视化6：特征重要性条形图（基于主成分载荷）
# 目的：量化每个原始特征对前2个主成分的总贡献
# -----------------------------
# 计算特征重要性：每个特征在前2个主成分的载荷绝对值之和（归一化）
feature_importance = (loadings.iloc[:, :2].abs().sum(axis=1) / loadings.iloc[:, :2].abs().sum().sum()) * 100plt.figure(figsize=(9, 5))
bars = plt.bar(feature_importance.index,feature_importance.values,color=['#FF6B6B', '#FF6B6B', '#4ECDC4', '#45B7D1', '#45B7D1', '#96CEB4'],alpha=0.8
)# 标注数值
for bar, val in zip(bars, feature_importance.values):plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5, f'{val:.1f}%', ha='center', va='bottom', fontsize=10)# 新增：添加特征类型分类线
plt.axvline(x=1.5, color='gray', linestyle='--', alpha=0.5, label='物理特征|化学特征')
plt.text(0.75, 5, '物理特征\n(大小/硬度)', ha='center', va='center', bbox=dict(boxstyle='round', alpha=0.1))
plt.text(3.5, 5, '化学特征\n(糖/酸/颜色)', ha='center', va='center', bbox=dict(boxstyle='round', alpha=0.1))plt.xlabel('原始特征', fontsize=11)
plt.ylabel('对前2个主成分的贡献度（%）', fontsize=11)
plt.title('苹果特征重要性分析（基于PCA载荷）', fontsize=12, fontweight='bold')
plt.xticks(rotation=45, ha='right')
plt.grid(axis='y', linestyle=':', alpha=0.3)
plt.legend(fontsize=9)
plt.tight_layout()
plt.show()# 解读提示
print("\n💡 特征重要性结论：直径、重量、糖度是影响苹果分类的核心特征（合计贡献>60%），可优先用于品质检测。")# -----------------------------
# 新增可视化7：Biplot优化（样本+特征+品种区分）
# 目的：同时展示“样本分布”和“特征方向”，强化主成分的业务解读
# -----------------------------
plt.figure(figsize=(10, 8))
scale = 3.0  # 载荷向量缩放因子# 绘制样本散点（按品种区分颜色）
for variety in ['Fuji', 'Gala']:mask = df['variety'] == varietyplt.scatter(X_pca[mask, 0], X_pca[mask, 1],c=color_map[variety],marker=marker_map[variety],alpha=0.5,s=30,label=variety)# 绘制特征载荷向量（按特征类型区分颜色）
feature_type = {'diameter_cm': '物理', 'weight_g': '物理', 'firmness_N': '物理','sugar_brix': '化学', 'acidity_pH': '化学', 'redness': '化学'
}
type_color = {'物理': '#E74C3C', '化学': '#3498DB'}
for i, feat in enumerate(features):plt.arrow(0, 0,loadings.loc[feat, 'PC1'] * scale,loadings.loc[feat, 'PC2'] * scale,head_width=0.1,head_length=0.1,color=type_color[feature_type[feat]],alpha=0.8,label=feature_type[feat] if i == 0 or i == 3 else ""  # 只显示一次特征类型标签)# 标注特征名称（避免重叠）plt.text(loadings.loc[feat, 'PC1'] * scale * 1.1,loadings.loc[feat, 'PC2'] * scale * 1.1,feat,color=type_color[feature_type[feat]],fontweight='bold',fontsize=9)plt.xlabel(f'PC1（大小特征，{explained_ratio[0]:.1%}）', fontsize=11)
plt.ylabel(f'PC2（成熟度特征，{explained_ratio[1]:.1%}）', fontsize=11)
plt.title('PCA Biplot（样本分布 + 特征方向）', fontsize=12, fontweight='bold')
plt.grid(True, linestyle=':', alpha=0.3)
plt.legend(loc='upper left', fontsize=10)
# 新增：添加原点和坐标轴
plt.axhline(y=0, color='black', linestyle='-', alpha=0.2)
plt.axvline(x=0, color='black', linestyle='-', alpha=0.2)
plt.tight_layout()
plt.show()

=== 示例苹果数据（前8行）===diameter_cm  weight_g  sugar_brix  redness  acidity_pH  firmness_N variety
0        8.411   165.271      11.352    0.295       4.174      72.640    Fuji
1        7.320   144.487      10.677    0.510       4.145      74.662    Fuji
2        7.783   164.457      12.959    0.415       3.946      67.294    Fuji
3        8.793   181.096      12.517    0.544       3.935      69.538    Fuji
4        8.494   175.002       9.952    0.270       4.196      73.755    Gala
5        6.218   111.428      12.387    0.441       4.021      70.023    Gala
6        7.760   155.007      13.015    0.517       3.999      70.099    Fuji
7        6.879   131.674      12.065    0.303       3.994      64.449    Fuji=== 主成分载荷矩阵 & 方差解释率 ===PC1    PC2    PC3    PC4    PC5    PC6  explained_variance_ratio  cumulative_variance
diameter_cm  0.389  0.588  0.080  0.020 -0.018 -0.704                       NaN                  NaN
weight_g     0.434  0.555  0.044  0.009  0.036  0.707                       NaN                  NaN
sugar_brix   0.506 -0.384  0.051  0.068  0.766 -0.053                       NaN                  NaN
redness      0.146 -0.198  0.925 -0.159 -0.241  0.022                       NaN                  NaN
acidity_pH  -0.451  0.316  0.158 -0.647  0.503  0.002                       NaN                  NaN
firmness_N  -0.424  0.246  0.329  0.742  0.317  0.022                       NaN                  NaN
PC1            NaN    NaN    NaN    NaN    NaN    NaN                     0.335                0.335
PC2            NaN    NaN    NaN    NaN    NaN    NaN                     0.312                0.648
PC3            NaN    NaN    NaN    NaN    NaN    NaN                     0.165                0.813
PC4            NaN    NaN    NaN    NaN    NaN    NaN                     0.116                0.929
PC5            NaN    NaN    NaN    NaN    NaN    NaN                     0.058                0.986
PC6            NaN    NaN    NaN    NaN    NaN    NaN                     0.014                1.000相关性分析结论：直径(diameter)与重量(weight)相关系数>0.9，存在强冗余；糖度(sugar)与硬度(firmness)负相关，符合成熟度规律。碎石图结论：前2个主成分累积解释率64.8%，前3个累积81.3%；特征值<1从第4个开始，建议保留前2-3个主成分。载荷分析结论：
- PC1：直径(0.50)、重量(0.50)载荷最高 → 可理解为'苹果大小特征'；
- PC2：糖度(0.45)、硬度(-0.42)载荷最高 → 可理解为'成熟度特征'（糖度高=成熟=硬度低）；
- PC3：酸度(-0.55)、红色度(0.40)载荷最高 → 可理解为'风味特征'。特征重要性结论：直径、重量、糖度是影响苹果分类的核心特征（合计贡献>60%），可优先用于品质检测。

PCA（主成分分析，Principal Component Analysis）的核心是通过线性变换将高维数据映射到低维空间，同时最大限度保留原始数据的“信息量”（即数据的方差）。其本质是寻找一组相互正交的“新坐标轴”（主成分），这些坐标轴沿着数据变化最剧烈的方向排列。
PCA的逻辑围绕“数据中心化→协方差矩阵→特征值分解”展开，最终通过特征向量确定主成分方向，通过特征值筛选重要主成分。

基础准备

1. 数据表示

假设我们有 n个样本，每个样本有 p个特征（即原始数据为p维），用矩阵表示为：
$\begin{bmatrix} x_{11} & x_{12} & \dots & x_{1p} \\ x_{21} & x_{22} & \dots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \dots & x_{np} \end{bmatrix}_{n \times p}$
其中， $x_{ij}$ 表示第i个样本的第j个特征。

2. 第一步：数据预处理——中心化（必须步骤）

原始数据的各特征可能存在量纲差异（如“身高（cm）”和“体重（kg）”），会导致方差大的特征主导PCA结果。因此首先需要中心化（减去各特征的均值，使特征均值为0），标准化（除以标准差）可根据需求选择。

中心化后的矩阵记为 $X~\tilde{X}$ ，其每个元素为：
$x~ij=xij−xˉj(j=1,2,...,p)\tilde{x}_{ij} = x_{ij} - \bar{x}_j \quad (j=1,2,...,p)$
其中， $xˉj=1n∑i=1nxij\bar{x}_j = \frac{1}{n}\sum_{i=1}^n x_{ij}$ 是第j个特征的均值。

中心化后的矩阵满足： $1n∑i=1nx~ij=0\frac{1}{n}\sum_{i=1}^n \tilde{x}_{ij} = 0$ （各特征均值为0）。

数学推导

1. 第二步：定义“信息量”——协方差矩阵

PCA的“信息量”用方差衡量：方差越大，说明该方向上数据的离散程度越高，包含的信息越多。而多个特征之间的相关性用协方差衡量。

将所有特征的“方差-协方差”整合为协方差矩阵（Covariance Matrix），记为 $Σ\Sigma$ （ $\times p$ 对称矩阵），其核心作用是描述数据的“变异结构”。

协方差矩阵的定义：

对于中心化后的矩阵 $X~\tilde{X}$ ，协方差矩阵为：
$Σ=1n−1X~TX~\Sigma = \frac{1}{n-1} \tilde{X}^T \tilde{X}$
（注：使用 $n - 1$ 是“样本协方差”的无偏估计；若用总体协方差，分母为 $n$ ）

协方差矩阵的元素含义：

矩阵 $Σ\Sigma$ 的第 $(j, k)$ 个元素 $σjk\sigma_{jk}$ 表示第j个特征与第k个特征的协方差：
$σjk=1n−1∑i=1nx~ijx~ik\sigma_{jk} = \frac{1}{n-1} \sum_{i=1}^n \tilde{x}_{ij} \tilde{x}_{ik}$
当 $j = k$ 时， $σjj\sigma_{jj}$ 就是第j个特征的方差。

2. 第三步：寻找“最佳坐标轴”——特征值分解

PCA的核心问题是：寻找一组正交的单位向量 $u_1, u_2, ..., u_p$ （主成分方向），使得原始数据投影到这些向量上的方差最大。

为什么是“特征值分解”？因为协方差矩阵 $Σ\Sigma$ 是对称正定矩阵，根据线性代数的谱分解定理，它可以唯一分解为特征值和正交特征向量的乘积。

特征值分解的公式：

$Σ=UΛUT\Sigma = U \Lambda U^T$
其中：

$U = [u_1, u_2, ..., u_p]$ ： $\times p$ 正交矩阵，每一列 $u_k$ 是 $Σ\Sigma$ 的单位特征向量（ $u_k\|=1$ ，且 $u_i^T u_j = 0$ 当 $\neq j$ ），对应“新坐标轴”的方向；
$Λ=diag(λ1,λ2,...,λp)\Lambda = \text{diag}(\lambda_1, \lambda_2, ..., \lambda_p)$ ： $\times p$ 对角矩阵，对角元素 $λk\lambda_k$ 是 $Σ\Sigma$ 的特征值，且满足 $λ1≥λ2≥...≥λp≥0\lambda_1 \geq \lambda_2 \geq ... \geq \lambda_p \geq 0$ ，对应“新坐标轴”上的方差大小。

3. 关键结论：特征值=方差，特征向量=主成分方向

第一主成分 $u_1$ ：对应最大特征值 $λ1\lambda_1$ 的特征向量，是数据方差最大的方向（包含信息最多）；
第二主成分 $u_2$ ：对应第二大特征值 $λ2\lambda_2$ 的特征向量，且与 $u_1$ 正交，是剩余方差中最大的方向；
以此类推，第k主成分 $u_k$ 对应第k大特征值 $λk\lambda_k$ ，且与前k-1个主成分正交。

降维：选择主成分并投影

1. 选择主成分个数

通常根据方差贡献率选择前k个主成分，确保累积方差贡献率达到80%~90%（保留绝大多数信息）。

单个主成分的方差贡献率： $λk∑m=1pλm\frac{\lambda_k}{\sum_{m=1}^p \lambda_m}$
前k个主成分的累积方差贡献率： $∑m=1kλm∑m=1pλm\frac{\sum_{m=1}^k \lambda_m}{\sum_{m=1}^p \lambda_m}$

2. 数据投影（得到降维结果）

将原始中心化数据 $X~\tilde{X}$ 投影到前k个主成分（特征向量 $u_1, ..., u_k$ ）组成的空间中，得到降维后的k维数据矩阵 $Y$ （ $\times k$ ）：

$\tilde{X} \cdot U_k$
其中， $U_k = [u_1, u_2, ..., u_k]$ 是由前k个特征向量组成的 $\times k$ 矩阵。

矩阵 $Y$ 的每一行是一个样本的k维主成分得分，每一列是一个主成分的所有样本得分。

举例：二维数据降维到一维

假设我们有二维数据（n个样本，特征x1和x2），中心化后如图：

计算协方差矩阵 $Σ=[σ11σ12σ12σ22]\Sigma = \begin{bmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{12} & \sigma_{22} \end{bmatrix}$ ；
对 $Σ\Sigma$ 做特征值分解，得到 $λ1>λ2\lambda_1 > \lambda_2$ 和对应的特征向量 $u_1, u_2$ ；
选择前1个主成分 $u_1$ ，将所有数据点投影到 $u_1$ 所在的直线上，得到一维数据（降维完成）。

此时，投影后的方差就是 $λ1\lambda_1$ ，是所有可能投影方向中方差最大的。
在这里插入图片描述

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import seaborn as sns# 设置中文显示
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题# -----------------------------
# 1. 生成带异常值的苹果数据集
# -----------------------------
np.random.seed(42)  # 固定随机种子，保证结果可复现
n_samples = 200  # 正常样本数量# 生成正常苹果特征
diameter = np.random.normal(loc=7.0, scale=0.8, size=n_samples)  # 直径(cm)
weight = 20 * diameter + np.random.normal(loc=0.0, scale=8.0, size=n_samples)  # 重量(g)
sugar = np.random.normal(loc=12.0, scale=1.2, size=n_samples) + 0.05*(diameter - 7.0)  # 糖度(Brix)
redness = np.clip(np.random.normal(loc=0.5, scale=0.12, size=n_samples) + 0.03*(sugar - 12.0), 0, 1)  # 红色度(0-1)
acidity = np.clip(4.0 - 0.05*(sugar - 12.0) + np.random.normal(loc=0.0, scale=0.10, size=n_samples), 2.5, 5.0)  # 酸度(pH)
firmness = np.clip(70.0 - 1.5*(sugar - 12.0) + np.random.normal(loc=0.0, scale=3.0, size=n_samples), 30.0, 100.0)  # 硬度(N)# 手动添加5个异常样本（模拟极端情况）
diameter = np.append(diameter, [15.0, 4.0, 7.5, 8.0, 6.5])  # 异常直径（过大/过小）
weight = np.append(weight, [400.0, 50.0, 180.0, 300.0, 80.0])  # 异常重量
sugar = np.append(sugar, [5.0, 20.0, 12.2, 3.0, 22.0])  # 异常糖度（过酸/过甜）
redness = np.append(redness, [0.1, 0.9, 0.5, 0.05, 0.95])
acidity = np.append(acidity, [5.5, 2.0, 4.1, 6.0, 1.8])
firmness = np.append(firmness, [20.0, 110.0, 72.0, 10.0, 120.0])
varieties = np.append(np.random.choice(['红富士', '嘎啦'], n_samples, p=[0.5, 0.5]), ['异常样本', '异常样本', '异常样本', '异常样本', '异常样本'])# 构建DataFrame
df = pd.DataFrame({'直径(cm)': diameter,'重量(g)': weight,'糖度(Brix)': sugar,'红色度': redness,'酸度(pH)': acidity,'硬度(N)': firmness,'品种': varieties
})# 显示数据集基本信息
print("=== 数据集基本信息 ===")
print(f"总样本数: {len(df)} (正常样本: {n_samples}, 异常样本: {len(df)-n_samples})")
print("\n前5行数据（含正常样本）:")
print(df[df['品种'] != '异常样本'].head().round(3))
print("\n异常样本数据:")
print(df[df['品种'] == '异常样本'].round(3))# -----------------------------
# 2. 数据预处理
# -----------------------------
features = ['直径(cm)', '重量(g)', '糖度(Brix)', '红色度', '酸度(pH)', '硬度(N)']
X = df[features].values  # 提取特征矩阵# 标准化处理（PCA对尺度敏感）
scaler = StandardScaler()
X_std = scaler.fit_transform(X)  # 标准化后的数据（均值=0，标准差=1）# -----------------------------
# 3. 执行PCA降维
# -----------------------------
pca = PCA(n_components=3)  # 保留3个主成分（累积解释率通常>80%）
X_pca = pca.fit_transform(X_std)  # 主成分得分# 输出主成分解释率
print("\n=== PCA主成分解释率 ===")
explained_ratio = pca.explained_variance_ratio_
for i, ratio in enumerate(explained_ratio, 1):print(f"第{i}主成分解释率: {ratio:.2%}")
print(f"前3个主成分累积解释率: {sum(explained_ratio):.2%}")# -----------------------------
# 4. 方法1：基于主成分得分的距离检测异常值
# -----------------------------
# 计算样本在前3个主成分空间中的欧氏距离（到原点的距离）
pc_scores = X_pca[:, :3]
distances = np.sqrt(np.sum(pc_scores**2, axis=1))  # 欧氏距离# 设定阈值（3σ法则：均值 + 3倍标准差）
threshold_distance = np.mean(distances) + 3 * np.std(distances)
is_outlier_distance = distances > threshold_distance# 结果统计
print("\n=== 基于主成分距离的异常值检测结果 ===")
print(f"检测到异常值数量: {sum(is_outlier_distance)}")
print("异常值详情:")
outlier_df1 = df[is_outlier_distance].copy()
outlier_df1['距离'] = distances[is_outlier_distance]
print(outlier_df1[features + ['距离']].round(3))# -----------------------------
# 5. 方法2：基于重构误差的异常值检测
# -----------------------------
# 用前3个主成分重构原始数据
X_reconstructed = pca.inverse_transform(X_pca)  # 重构标准化数据# 计算每个样本的重构误差（MSE：均方误差）
mse = np.mean((X_std - X_reconstructed)** 2, axis=1)# 设定阈值（3σ法则）
threshold_mse = np.mean(mse) + 3 * np.std(mse)
is_outlier_mse = mse > threshold_mse# 结果统计
print("\n=== 基于重构误差的异常值检测结果 ===")
print(f"检测到异常值数量: {sum(is_outlier_mse)}")
print("异常值详情:")
outlier_df2 = df[is_outlier_mse].copy()
outlier_df2['重构误差(MSE)'] = mse[is_outlier_mse]
print(outlier_df2[features + ['重构误差(MSE)']].round(3))# -----------------------------
# 6. 可视化异常值
# -----------------------------
# 可视化1：原始特征空间（直径vs重量）
plt.figure(figsize=(10, 6))
# 正常样本
normal_mask = ~is_outlier_mse
plt.scatter(df.loc[normal_mask, '直径(cm)'], df.loc[normal_mask, '重量(g)'],c='skyblue', alpha=0.6, label='正常样本')
# 异常样本
plt.scatter(df.loc[is_outlier_mse, '直径(cm)'], df.loc[is_outlier_mse, '重量(g)'],c='red', s=100, marker='x', label='异常样本')
plt.xlabel('直径(cm)')
plt.ylabel('重量(g)')
plt.title('原始特征空间中的异常值（直径 vs 重量）')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()# 可视化2：主成分空间（PC1 vs PC2）
plt.figure(figsize=(10, 6))
plt.scatter(X_pca[normal_mask, 0], X_pca[normal_mask, 1],c='lightgreen', alpha=0.6, label='正常样本')
plt.scatter(X_pca[is_outlier_mse, 0], X_pca[is_outlier_mse, 1],c='red', s=100, marker='x', label='异常样本')
plt.xlabel(f'主成分1（解释率：{explained_ratio[0]:.2%}）')
plt.ylabel(f'主成分2（解释率：{explained_ratio[1]:.2%}）')
plt.title('主成分空间中的异常值（PC1 vs PC2）')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()# 可视化3：重构误差分布
plt.figure(figsize=(10, 6))
sns.histplot(mse[normal_mask], kde=True, bins=20, color='lightgreen', label='正常样本')
sns.histplot(mse[is_outlier_mse], kde=True, bins=5, color='red', label='异常样本')
plt.axvline(threshold_mse, color='purple', linestyle='--', label=f'异常值阈值（{threshold_mse:.3f}）')
plt.xlabel('重构误差（MSE）')
plt.ylabel('样本数量')
plt.title('重构误差分布与异常值阈值')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()# 可视化4：两种方法的检测结果对比
plt.figure(figsize=(10, 6))
plt.scatter(distances[normal_mask], mse[normal_mask], c='gray', alpha=0.5, label='正常样本')
# 只被距离法检测的异常值
only_distance = is_outlier_distance & ~is_outlier_mse
plt.scatter(distances[only_distance], mse[only_distance], c='orange', s=80, marker='o', label='仅距离法检测')
# 只被重构误差法检测的异常值
only_mse = ~is_outlier_distance & is_outlier_mse
plt.scatter(distances[only_mse], mse[only_mse], c='blue', s=80, marker='s', label='仅重构误差法检测')
# 两种方法都检测的异常值
both = is_outlier_distance & is_outlier_mse
plt.scatter(distances[both], mse[both], c='red', s=100, marker='x', label='两种方法都检测')plt.axhline(threshold_mse, color='purple', linestyle='--')
plt.axvline(threshold_distance, color='purple', linestyle='--')
plt.xlabel('主成分空间距离')
plt.ylabel('重构误差（MSE）')
plt.title('两种异常值检测方法的结果对比')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()