用Python实战PCA异常检测：手把手教你计算T²和SPE统计量（附完整代码）

用Python实战PCA异常检测：手把手教你计算T²和SPE统计量（附完整代码）

在工业过程监控、金融风控或设备故障预警等场景中，异常检测始终是数据分析的核心挑战之一。传统单变量控制图难以捕捉高维数据中的复杂关系，而主成分分析（PCA）通过降维技术将多变量空间转换为低维特征空间，配合T²和SPE两个统计量，能有效识别数据中的异常模式。本文将用Python代码完整演示从数据标准化、PCA建模到双指标计算的闭环流程，并提供可直接复用的函数库。

1. 环境准备与数据预处理

首先导入必要的库并生成模拟数据。我们使用sklearn的make_blobs创建包含5%异常点的三维数据集：

import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import make_blobs from sklearn.preprocessing import StandardScaler # 生成含5%异常点的三维数据 X, y = make_blobs(n_samples=1000, centers=1, n_features=3, cluster_std=1.0, random_state=42) X[:50] += 5 # 添加异常点 # 数据标准化 scaler = StandardScaler() X_scaled = scaler.fit_transform(X)

标准化是PCA前的关键步骤，否则方差大的变量会主导主成分方向。通过StandardScaler将各特征缩放到均值为0、标准差为1的标准正态分布。

常见预处理错误排查：

• 未处理缺失值会导致PCA计算错误 → 用SimpleImputer填充
• 未标准化使主成分偏向高方差特征 → 必须使用StandardScaler
• 测试集单独标准化 → 应使用训练集的均值和标准差

2. PCA建模与主成分选择

使用sklearn.decomposition.PCA进行建模，并通过累积解释方差比确定保留的主成分数：

from sklearn.decomposition import PCA pca = PCA(n_components=0.95) # 保留95%方差的成分 pca.fit(X_scaled) # 可视化解释方差比 plt.plot(np.cumsum(pca.explained_variance_ratio_)) plt.xlabel('Number of Components') plt.ylabel('Cumulative Explained Variance') plt.axhline(y=0.95, color='r', linestyle='--') plt.show()

关键输出参数说明：

3. T²统计量计算与可视化

Hotelling's T²统计量衡量样本在模型空间中的变异程度，计算步骤如下：

1. 计算主成分得分：T = X_scaled @ P
2. 计算T²值：T² = np.sum(T**2 / Λ, axis=1)

完整实现代码：

def calculate_t2(X, pca): # 主成分得分 scores = pca.transform(X) # 特征值倒数对角矩阵 lambda_inv = np.diag(1 / pca.explained_variance_) # T²计算 t2 = np.array([t @ lambda_inv @ t.T for t in scores]) return t2 # 计算控制限 n, k = X_scaled.shape[0], pca.n_components_ alpha = 0.01 # 显著性水平 from scipy.stats import f ucl_t2 = (n-1)*k / (n-k) * f.ppf(1-alpha, k, n-k) # 绘制T²监控图 t2_values = calculate_t2(X_scaled, pca) plt.plot(t2_values) plt.axhline(ucl_t2, color='r', linestyle='--') plt.title('T² Monitoring Chart') plt.show()

调试技巧：

• 若所有点都超限 → 检查Λ是否取倒数
• 若图形波动异常 → 确认是否先进行标准化
• 控制限不合理 → 验证F分布自由度参数

4. SPE统计量计算与阈值确定

平方预测误差(SPE)反映未被主成分解释的变异，计算流程：

1. 重构原始数据：X_reconstructed = T @ P.T
2. 计算残差范数：SPE = np.sum((X - X_reconstructed)**2, axis=1)

Python实现：

def calculate_spe(X, pca): scores = pca.transform(X) X_recon = pca.inverse_transform(scores) spe = np.sum((X - X_recon)**2, axis=1) return spe # SPE控制限计算 theta = [sum(pca.explained_variance_[k:]**i) for i in [1,2,3]] h0 = 1 - 2*theta[0]*theta[2] / (3*theta[1]**2) ca = 3 # 严格度系数 ucl_spe = theta[0] * (ca*np.sqrt(2*theta[1]*h0**2)/theta[0] + 1 + theta[1]*h0*(h0-1)/theta[0]**2)**(1/h0) # 绘制SPE监控图 spe_values = calculate_spe(X_scaled, pca) plt.plot(spe_values) plt.axhline(ucl_spe, color='r', linestyle='--') plt.title('SPE Monitoring Chart') plt.show()

实际应用建议：

• 批处理数据时建议使用滑动窗口计算
• 动态过程可考虑自适应控制限
• 结合T²和SPE的综合指标能提升检测灵敏度

5. 双指标联合分析与案例解读

将两个统计量结合能更全面识别异常类型：

# 创建联合监控图 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,6)) ax1.plot(t2_values) ax1.axhline(ucl_t2, color='r', linestyle='--') ax1.set_title('T² Chart') ax2.plot(spe_values) ax2.axhline(ucl_spe, color='r', linestyle='--') ax2.set_title('SPE Chart') plt.tight_layout() plt.show() # 异常点标记 anomalies = np.where((t2_values > ucl_t2) | (spe_values > ucl_spe))[0] print(f"Detected anomalies at indices: {anomalies}")

异常类型诊断矩阵：

6. 性能优化与生产部署

对于实时监测场景，需优化计算效率：

# 使用numba加速 from numba import jit @jit(nopython=True) def fast_t2(scores, lambda_inv): return np.array([t @ lambda_inv @ t.T for t in scores]) # 增量PCA处理流数据 from sklearn.decomposition import IncrementalPCA ipca = IncrementalPCA(n_components=2, batch_size=100) for batch in np.array_split(X_scaled, 10): ipca.partial_fit(batch)

部署架构建议：

1. 离线阶段：训练PCA模型并保存参数
2. 在线阶段：
   • 加载模型对新数据标准化
   • 实时计算T²和SPE
   • 触发报警机制

完整代码已封装为可复用类，包含以下方法：

class PCAMonitor: def __init__(self, alpha=0.01): self.alpha = alpha def fit(self, X): # 训练流程... def transform(self, X): # 计算双指标... def plot_control_charts(self): # 绘制监控图...