第12章：主成分分析
主成分分析（Principal Component Analysis, PCA）是最重要的降维技术之一。它通过线性变换将高维数据投影到低维空间，在保持数据主要信息的同时减少特征维度。本章将详细介绍PCA的原理、实现和应用。
12.1 什么是主成分分析？
PCA是一种无监督的降维技术，它寻找数据中方差最大的方向（主成分），并将数据投影到这些方向上，从而实现降维。
12.1.1 PCA的目标
- 降维：减少特征数量，简化数据结构
- 去相关：消除特征间的线性相关性
- 保持信息：尽可能保留原始数据的信息
- 可视化：将高维数据投影到2D或3D空间
12.1.2 PCA的应用
- 数据压缩：减少存储空间和计算复杂度
- 噪声去除：通过保留主要成分去除噪声
- 特征提取：提取最重要的特征组合
- 数据可视化：高维数据的可视化展示
- 预处理：为其他机器学习算法准备数据
12.1.3 PCA的数学原理
PCA通过以下步骤实现降维：
- 数据标准化：使各特征具有相同的量纲
- 计算协方差矩阵：衡量特征间的相关性
- 特征值分解：找到协方差矩阵的特征向量和特征值
- 选择主成分：根据特征值大小选择主要方向
- 数据投影：将原始数据投影到主成分空间
12.2 准备环境和数据
python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_iris, load_wine, load_breast_cancer, make_classification
from sklearn.decomposition import PCA, IncrementalPCA, KernelPCA
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, classification_report
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')
# 设置随机种子
np.random.seed(42)
# 设置图形样式
plt.style.use('seaborn-v0_8')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False12.3 PCA基本原理演示
12.3.1 二维数据的PCA
python
def demonstrate_pca_2d():
"""演示二维数据的PCA过程"""
# 创建相关的二维数据
np.random.seed(42)
n_samples = 200
# 生成相关数据
x1 = np.random.normal(0, 1, n_samples)
x2 = 0.8 * x1 + 0.6 * np.random.normal(0, 1, n_samples)
X_2d = np.column_stack([x1, x2])
print("PCA基本原理演示:")
print("1. 原始数据具有相关性")
print("2. PCA找到方差最大的方向作为第一主成分")
print("3. 第二主成分与第一主成分正交")
print("4. 数据投影到主成分空间实现降维")
# 标准化数据
scaler = StandardScaler()
X_2d_scaled = scaler.fit_transform(X_2d)
# 应用PCA
pca_2d = PCA(n_components=2)
X_2d_pca = pca_2d.fit_transform(X_2d_scaled)
# 获取主成分
components = pca_2d.components_
explained_variance_ratio = pca_2d.explained_variance_ratio_
print(f"\n主成分分析结果:")
print(f"第一主成分解释方差比: {explained_variance_ratio[0]:.3f}")
print(f"第二主成分解释方差比: {explained_variance_ratio[1]:.3f}")
print(f"累积解释方差比: {np.cumsum(explained_variance_ratio)}")
# 可视化
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
fig.suptitle('PCA原理演示', fontsize=16)
# 原始数据
axes[0, 0].scatter(X_2d_scaled[:, 0], X_2d_scaled[:, 1], alpha=0.6, color='blue')
axes[0, 0].set_title('标准化后的原始数据')
axes[0, 0].set_xlabel('特征 1 (标准化)')
axes[0, 0].set_ylabel('特征 2 (标准化)')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].axis('equal')
# 显示主成分方向
axes[0, 1].scatter(X_2d_scaled[:, 0], X_2d_scaled[:, 1], alpha=0.6, color='blue')
# 绘制主成分向量
mean_point = np.mean(X_2d_scaled, axis=0)
for i, (component, variance_ratio) in enumerate(zip(components, explained_variance_ratio)):
# 缩放向量以便可视化
vector = component * np.sqrt(pca_2d.explained_variance_[i]) * 2
axes[0, 1].arrow(mean_point[0], mean_point[1], vector[0], vector[1],
head_width=0.1, head_length=0.1, fc=f'C{i}', ec=f'C{i}',
linewidth=3, label=f'PC{i+1} ({variance_ratio:.2f})')
axes[0, 1].set_title('主成分方向')
axes[0, 1].set_xlabel('特征 1 (标准化)')
axes[0, 1].set_ylabel('特征 2 (标准化)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].axis('equal')
# PCA变换后的数据
axes[1, 0].scatter(X_2d_pca[:, 0], X_2d_pca[:, 1], alpha=0.6, color='red')
axes[1, 0].set_title('PCA变换后的数据')
axes[1, 0].set_xlabel('第一主成分')
axes[1, 0].set_ylabel('第二主成分')
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].axis('equal')
# 只保留第一主成分的重构
X_1d_pca = pca_2d.transform(X_2d_scaled)
X_1d_pca[:, 1] = 0 # 将第二主成分设为0
X_reconstructed = pca_2d.inverse_transform(X_1d_pca)
axes[1, 1].scatter(X_2d_scaled[:, 0], X_2d_scaled[:, 1], alpha=0.3, color='blue', label='原始数据')
axes[1, 1].scatter(X_reconstructed[:, 0], X_reconstructed[:, 1], alpha=0.6, color='red', label='重构数据')
axes[1, 1].set_title('降维重构 (仅保留第一主成分)')
axes[1, 1].set_xlabel('特征 1 (标准化)')
axes[1, 1].set_ylabel('特征 2 (标准化)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].axis('equal')
plt.tight_layout()
plt.show()
return X_2d_scaled, X_2d_pca, pca_2d
X_2d_scaled, X_2d_pca, pca_2d = demonstrate_pca_2d()
```###
12.3.2 协方差矩阵和特征值分解
```python
def analyze_covariance_and_eigenvalues():
"""分析协方差矩阵和特征值分解"""
print("协方差矩阵和特征值分解分析:")
print("=" * 40)
# 计算协方差矩阵
cov_matrix = np.cov(X_2d_scaled.T)
print("协方差矩阵:")
print(cov_matrix)
# 手动进行特征值分解
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
# 按特征值大小排序
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
print(f"\n特征值: {eigenvalues}")
print(f"特征向量:\n{eigenvectors}")
# 计算解释方差比
explained_variance_ratio_manual = eigenvalues / np.sum(eigenvalues)
print(f"\n手动计算的解释方差比: {explained_variance_ratio_manual}")
print(f"PCA计算的解释方差比: {pca_2d.explained_variance_ratio_}")
# 可视化特征值和特征向量
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# 协方差矩阵热力图
sns.heatmap(cov_matrix, annot=True, cmap='coolwarm', center=0, ax=axes[0])
axes[0].set_title('协方差矩阵')
# 特征值
axes[1].bar(range(1, len(eigenvalues) + 1), eigenvalues, color='skyblue', alpha=0.7)
axes[1].set_title('特征值')
axes[1].set_xlabel('主成分')
axes[1].set_ylabel('特征值')
axes[1].grid(True, alpha=0.3)
# 解释方差比
axes[2].bar(range(1, len(explained_variance_ratio_manual) + 1),
explained_variance_ratio_manual, color='lightgreen', alpha=0.7)
axes[2].set_title('解释方差比')
axes[2].set_xlabel('主成分')
axes[2].set_ylabel('解释方差比')
axes[2].grid(True, alpha=0.3)
# 添加数值标签
for i, (val, ratio) in enumerate(zip(eigenvalues, explained_variance_ratio_manual)):
axes[1].text(i, val + 0.01, f'{val:.3f}', ha='center')
axes[2].text(i, ratio + 0.01, f'{ratio:.3f}', ha='center')
plt.tight_layout()
plt.show()
return cov_matrix, eigenvalues, eigenvectors
cov_matrix, eigenvalues, eigenvectors = analyze_covariance_and_eigenvalues()12.4 高维数据的PCA
12.4.1 鸢尾花数据集PCA
python
def pca_iris_analysis():
"""鸢尾花数据集的PCA分析"""
# 加载鸢尾花数据集
iris = load_iris()
X_iris = iris.data
y_iris = iris.target
feature_names = iris.feature_names
target_names = iris.target_names
print("鸢尾花数据集PCA分析:")
print(f"原始数据形状: {X_iris.shape}")
print(f"特征名称: {feature_names}")
# 标准化数据
scaler = StandardScaler()
X_iris_scaled = scaler.fit_transform(X_iris)
# 应用PCA
pca_iris = PCA()
X_iris_pca = pca_iris.fit_transform(X_iris_scaled)
# 分析结果
explained_variance_ratio = pca_iris.explained_variance_ratio_
cumulative_variance_ratio = np.cumsum(explained_variance_ratio)
print(f"\n各主成分解释方差比:")
for i, ratio in enumerate(explained_variance_ratio):
print(f"PC{i+1}: {ratio:.4f}")
print(f"\n累积解释方差比:")
for i, cum_ratio in enumerate(cumulative_variance_ratio):
print(f"前{i+1}个主成分: {cum_ratio:.4f}")
# 可视化
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('鸢尾花数据集PCA分析', fontsize=16)
# 解释方差比
axes[0, 0].bar(range(1, len(explained_variance_ratio) + 1),
explained_variance_ratio, color='skyblue', alpha=0.7)
axes[0, 0].set_title('各主成分解释方差比')
axes[0, 0].set_xlabel('主成分')
axes[0, 0].set_ylabel('解释方差比')
axes[0, 0].grid(True, alpha=0.3)
# 累积解释方差比
axes[0, 1].plot(range(1, len(cumulative_variance_ratio) + 1),
cumulative_variance_ratio, 'ro-', linewidth=2, markersize=8)
axes[0, 1].axhline(y=0.95, color='red', linestyle='--', alpha=0.7, label='95%阈值')
axes[0, 1].set_title('累积解释方差比')
axes[0, 1].set_xlabel('主成分数量')
axes[0, 1].set_ylabel('累积解释方差比')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# 主成分载荷图
components_df = pd.DataFrame(
pca_iris.components_[:2].T,
columns=['PC1', 'PC2'],
index=feature_names
)
axes[0, 2].scatter(components_df['PC1'], components_df['PC2'], s=100)
for i, feature in enumerate(feature_names):
axes[0, 2].annotate(feature, (components_df.iloc[i, 0], components_df.iloc[i, 1]),
xytext=(5, 5), textcoords='offset points')
axes[0, 2].set_title('主成分载荷图')
axes[0, 2].set_xlabel('PC1')
axes[0, 2].set_ylabel('PC2')
axes[0, 2].grid(True, alpha=0.3)
# 2D投影
colors = ['red', 'blue', 'green']
for i, target_name in enumerate(target_names):
mask = y_iris == i
axes[1, 0].scatter(X_iris_pca[mask, 0], X_iris_pca[mask, 1],
c=colors[i], alpha=0.7, label=target_name)
axes[1, 0].set_title('前两个主成分的2D投影')
axes[1, 0].set_xlabel(f'PC1 ({explained_variance_ratio[0]:.2%})')
axes[1, 0].set_ylabel(f'PC2 ({explained_variance_ratio[1]:.2%})')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# 3D投影
ax_3d = fig.add_subplot(2, 3, 5, projection='3d')
for i, target_name in enumerate(target_names):
mask = y_iris == i
ax_3d.scatter(X_iris_pca[mask, 0], X_iris_pca[mask, 1], X_iris_pca[mask, 2],
c=colors[i], alpha=0.7, label=target_name)
ax_3d.set_title('前三个主成分的3D投影')
ax_3d.set_xlabel(f'PC1 ({explained_variance_ratio[0]:.2%})')
ax_3d.set_ylabel(f'PC2 ({explained_variance_ratio[1]:.2%})')
ax_3d.set_zlabel(f'PC3 ({explained_variance_ratio[2]:.2%})')
ax_3d.legend()
# 原始特征vs主成分
feature_comparison = pd.DataFrame({
'萼片长度': X_iris_scaled[:, 0],
'PC1': X_iris_pca[:, 0],
'类别': y_iris
})
for i, target_name in enumerate(target_names):
mask = y_iris == i
axes[1, 2].scatter(feature_comparison.loc[mask, '萼片长度'],
feature_comparison.loc[mask, 'PC1'],
c=colors[i], alpha=0.7, label=target_name)
axes[1, 2].set_title('原始特征 vs 第一主成分')
axes[1, 2].set_xlabel('萼片长度 (标准化)')
axes[1, 2].set_ylabel('PC1')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
return X_iris_scaled, X_iris_pca, pca_iris
X_iris_scaled, X_iris_pca, pca_iris = pca_iris_analysis()12.4.2 确定主成分数量
python
def determine_n_components():
"""确定最佳主成分数量的方法"""
# 使用葡萄酒数据集
wine = load_wine()
X_wine = wine.data
y_wine = wine.target
# 标准化
scaler = StandardScaler()
X_wine_scaled = scaler.fit_transform(X_wine)
# 计算所有主成分
pca_full = PCA()
pca_full.fit(X_wine_scaled)
explained_variance_ratio = pca_full.explained_variance_ratio_
cumulative_variance_ratio = np.cumsum(explained_variance_ratio)
print("确定主成分数量的方法:")
print("=" * 30)
# 方法1: 累积解释方差比阈值
threshold_95 = np.argmax(cumulative_variance_ratio >= 0.95) + 1
threshold_90 = np.argmax(cumulative_variance_ratio >= 0.90) + 1
print(f"方法1 - 累积解释方差比:")
print(f" 保留90%方差需要: {threshold_90} 个主成分")
print(f" 保留95%方差需要: {threshold_95} 个主成分")
# 方法2: Kaiser准则 (特征值 > 1)
eigenvalues = pca_full.explained_variance_
kaiser_components = np.sum(eigenvalues > 1)
print(f"\n方法2 - Kaiser准则 (特征值>1): {kaiser_components} 个主成分")
# 方法3: 碎石图分析
print(f"\n方法3 - 碎石图: 观察特征值的'肘部'")
# 可视化不同方法
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
fig.suptitle('确定主成分数量的方法', fontsize=16)
# 解释方差比
axes[0, 0].bar(range(1, len(explained_variance_ratio) + 1),
explained_variance_ratio, alpha=0.7, color='skyblue')
axes[0, 0].set_title('各主成分解释方差比')
axes[0, 0].set_xlabel('主成分')
axes[0, 0].set_ylabel('解释方差比')
axes[0, 0].grid(True, alpha=0.3)
# 累积解释方差比
axes[0, 1].plot(range(1, len(cumulative_variance_ratio) + 1),
cumulative_variance_ratio, 'bo-', linewidth=2, markersize=6)
axes[0, 1].axhline(y=0.90, color='red', linestyle='--', alpha=0.7, label='90%')
axes[0, 1].axhline(y=0.95, color='orange', linestyle='--', alpha=0.7, label='95%')
axes[0, 1].axvline(x=threshold_90, color='red', linestyle=':', alpha=0.7)
axes[0, 1].axvline(x=threshold_95, color='orange', linestyle=':', alpha=0.7)
axes[0, 1].set_title('累积解释方差比')
axes[0, 1].set_xlabel('主成分数量')
axes[0, 1].set_ylabel('累积解释方差比')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# 特征值 (Kaiser准则)
axes[1, 0].bar(range(1, len(eigenvalues) + 1), eigenvalues, alpha=0.7, color='lightgreen')
axes[1, 0].axhline(y=1, color='red', linestyle='--', alpha=0.7, label='Kaiser阈值')
axes[1, 0].set_title('特征值 (Kaiser准则)')
axes[1, 0].set_xlabel('主成分')
axes[1, 0].set_ylabel('特征值')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# 碎石图
axes[1, 1].plot(range(1, len(eigenvalues) + 1), eigenvalues, 'ro-', linewidth=2, markersize=8)
axes[1, 1].set_title('碎石图')
axes[1, 1].set_xlabel('主成分')
axes[1, 1].set_ylabel('特征值')
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 比较不同主成分数量的分类性能
print(f"\n不同主成分数量对分类性能的影响:")
print("主成分数\t累积方差比\t分类准确率")
print("-" * 40)
X_train, X_test, y_train, y_test = train_test_split(
X_wine_scaled, y_wine, test_size=0.2, random_state=42
)
n_components_list = [2, 3, 5, 8, 10, 13] # 13是总特征数
for n_comp in n_components_list:
# PCA降维
pca_temp = PCA(n_components=n_comp)
X_train_pca = pca_temp.fit_transform(X_train)
X_test_pca = pca_temp.transform(X_test)
# 分类
clf = LogisticRegression(random_state=42, max_iter=1000)
clf.fit(X_train_pca, y_train)
y_pred = clf.predict(X_test_pca)
accuracy = accuracy_score(y_test, y_pred)
cum_var_ratio = np.sum(pca_temp.explained_variance_ratio_)
print(f"{n_comp}\t\t{cum_var_ratio:.4f}\t\t{accuracy:.4f}")
return threshold_90, threshold_95, kaiser_components
threshold_90, threshold_95, kaiser_components = determine_n_components()12.5 PCA的变体
12.5.1 增量PCA
python
def incremental_pca_demo():
"""演示增量PCA"""
print("增量PCA (Incremental PCA):")
print("适用于大数据集，可以批量处理数据")
# 创建大数据集
X_large, y_large = make_classification(
n_samples=5000, n_features=50, n_informative=30,
n_redundant=10, random_state=42
)
# 标准化
scaler = StandardScaler()
X_large_scaled = scaler.fit_transform(X_large)
# 比较标准PCA和增量PCA
import time
# 标准PCA
start_time = time.time()
pca_standard = PCA(n_components=10)
X_pca_standard = pca_standard.fit_transform(X_large_scaled)
time_standard = time.time() - start_time
# 增量PCA
start_time = time.time()
pca_incremental = IncrementalPCA(n_components=10, batch_size=500)
X_pca_incremental = pca_incremental.fit_transform(X_large_scaled)
time_incremental = time.time() - start_time
print(f"\n性能比较:")
print(f"标准PCA时间: {time_standard:.4f}秒")
print(f"增量PCA时间: {time_incremental:.4f}秒")
# 比较结果相似性
# 由于符号可能不同，比较绝对值
correlation = np.corrcoef(np.abs(X_pca_standard[:, 0]), np.abs(X_pca_incremental[:, 0]))[0, 1]
print(f"第一主成分相关性: {correlation:.4f}")
# 可视化比较
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# 标准PCA结果
axes[0].scatter(X_pca_standard[:, 0], X_pca_standard[:, 1],
c=y_large, cmap='viridis', alpha=0.6, s=10)
axes[0].set_title('标准PCA')
axes[0].set_xlabel('PC1')
axes[0].set_ylabel('PC2')
# 增量PCA结果
axes[1].scatter(X_pca_incremental[:, 0], X_pca_incremental[:, 1],
c=y_large, cmap='viridis', alpha=0.6, s=10)
axes[1].set_title('增量PCA')
axes[1].set_xlabel('PC1')
axes[1].set_ylabel('PC2')
# 解释方差比比较
x_pos = np.arange(10)
width = 0.35
axes[2].bar(x_pos - width/2, pca_standard.explained_variance_ratio_,
width, label='标准PCA', alpha=0.7)
axes[2].bar(x_pos + width/2, pca_incremental.explained_variance_ratio_,
width, label='增量PCA', alpha=0.7)
axes[2].set_title('解释方差比比较')
axes[2].set_xlabel('主成分')
axes[2].set_ylabel('解释方差比')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
return pca_standard, pca_incremental
pca_standard, pca_incremental = incremental_pca_demo()12.5.2 核PCA
python
def kernel_pca_demo():
"""演示核PCA"""
print("核PCA (Kernel PCA):")
print("使用核技巧处理非线性数据")
# 创建非线性数据
from sklearn.datasets import make_circles, make_moons
# 同心圆数据
X_circles, y_circles = make_circles(n_samples=400, noise=0.1, factor=0.3, random_state=42)
# 月牙形数据
X_moons, y_moons = make_moons(n_samples=400, noise=0.1, random_state=42)
datasets = [
('同心圆', X_circles, y_circles),
('月牙形', X_moons, y_moons)
]
# 不同的核函数
kernels = ['linear', 'poly', 'rbf', 'sigmoid']
for dataset_name, X, y in datasets:
print(f"\n{dataset_name}数据集:")
fig, axes = plt.subplots(1, len(kernels) + 1, figsize=(20, 4))
fig.suptitle(f'{dataset_name}数据的核PCA', fontsize=16)
# 原始数据
axes[0].scatter(X[:, 0], X[:, 1], c=y, cmap='viridis', alpha=0.7)
axes[0].set_title('原始数据')
axes[0].set_xlabel('特征 1')
axes[0].set_ylabel('特征 2')
for i, kernel in enumerate(kernels):
try:
# 核PCA
if kernel == 'poly':
kpca = KernelPCA(n_components=2, kernel=kernel, degree=3, random_state=42)
elif kernel == 'rbf':
kpca = KernelPCA(n_components=2, kernel=kernel, gamma=1, random_state=42)
else:
kpca = KernelPCA(n_components=2, kernel=kernel, random_state=42)
X_kpca = kpca.fit_transform(X)
# 可视化
axes[i + 1].scatter(X_kpca[:, 0], X_kpca[:, 1], c=y, cmap='viridis', alpha=0.7)
axes[i + 1].set_title(f'{kernel.upper()}核')
axes[i + 1].set_xlabel('第一核主成分')
axes[i + 1].set_ylabel('第二核主成分')
print(f" {kernel}核: 成功")
except Exception as e:
axes[i + 1].text(0.5, 0.5, f'错误:\n{str(e)[:30]}...',
ha='center', va='center', transform=axes[i + 1].transAxes)
axes[i + 1].set_title(f'{kernel.upper()}核 (失败)')
print(f" {kernel}核: 失败 - {str(e)[:50]}")
plt.tight_layout()
plt.show()
kernel_pca_demo()
```#
# 12.6 PCA在机器学习中的应用
### 12.6.1 降维对分类性能的影响
```python
def pca_classification_performance():
"""分析PCA降维对分类性能的影响"""
# 使用乳腺癌数据集
cancer = load_breast_cancer()
X_cancer = cancer.data
y_cancer = cancer.target
print("PCA降维对分类性能的影响分析:")
print(f"原始数据形状: {X_cancer.shape}")
# 标准化
scaler = StandardScaler()
X_cancer_scaled = scaler.fit_transform(X_cancer)
# 分割数据
X_train, X_test, y_train, y_test = train_test_split(
X_cancer_scaled, y_cancer, test_size=0.2, random_state=42, stratify=y_cancer
)
# 测试不同的主成分数量
n_components_list = [2, 5, 10, 15, 20, 25, 30] # 原始30个特征
results = {
'n_components': [],
'variance_ratio': [],
'train_accuracy': [],
'test_accuracy': [],
'train_time': []
}
print("\n主成分数\t累积方差比\t训练准确率\t测试准确率\t训练时间")
print("-" * 70)
for n_comp in n_components_list:
# PCA降维
pca = PCA(n_components=n_comp)
X_train_pca = pca.fit_transform(X_train)
X_test_pca = pca.transform(X_test)
# 训练分类器
start_time = time.time()
clf = LogisticRegression(random_state=42, max_iter=1000)
clf.fit(X_train_pca, y_train)
train_time = time.time() - start_time
# 预测
y_train_pred = clf.predict(X_train_pca)
y_test_pred = clf.predict(X_test_pca)
# 计算准确率
train_acc = accuracy_score(y_train, y_train_pred)
test_acc = accuracy_score(y_test, y_test_pred)
# 累积方差比
cum_var_ratio = np.sum(pca.explained_variance_ratio_)
# 保存结果
results['n_components'].append(n_comp)
results['variance_ratio'].append(cum_var_ratio)
results['train_accuracy'].append(train_acc)
results['test_accuracy'].append(test_acc)
results['train_time'].append(train_time)
print(f"{n_comp}\t\t{cum_var_ratio:.4f}\t\t{train_acc:.4f}\t\t{test_acc:.4f}\t\t{train_time:.4f}s")
# 添加原始数据的结果
start_time = time.time()
clf_original = LogisticRegression(random_state=42, max_iter=1000)
clf_original.fit(X_train, y_train)
train_time_original = time.time() - start_time
y_train_pred_original = clf_original.predict(X_train)
y_test_pred_original = clf_original.predict(X_test)
train_acc_original = accuracy_score(y_train, y_train_pred_original)
test_acc_original = accuracy_score(y_test, y_test_pred_original)
print(f"原始数据\t1.0000\t\t{train_acc_original:.4f}\t\t{test_acc_original:.4f}\t\t{train_time_original:.4f}s")
# 可视化结果
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
fig.suptitle('PCA降维对分类性能的影响', fontsize=16)
# 准确率 vs 主成分数
axes[0, 0].plot(results['n_components'], results['train_accuracy'],
'bo-', label='训练准确率', linewidth=2, markersize=6)
axes[0, 0].plot(results['n_components'], results['test_accuracy'],
'ro-', label='测试准确率', linewidth=2, markersize=6)
axes[0, 0].axhline(y=test_acc_original, color='green', linestyle='--',
alpha=0.7, label='原始数据测试准确率')
axes[0, 0].set_title('准确率 vs 主成分数')
axes[0, 0].set_xlabel('主成分数')
axes[0, 0].set_ylabel('准确率')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# 累积方差比 vs 主成分数
axes[0, 1].plot(results['n_components'], results['variance_ratio'],
'go-', linewidth=2, markersize=6)
axes[0, 1].axhline(y=0.95, color='red', linestyle='--', alpha=0.7, label='95%阈值')
axes[0, 1].set_title('累积方差比 vs 主成分数')
axes[0, 1].set_xlabel('主成分数')
axes[0, 1].set_ylabel('累积方差比')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# 训练时间 vs 主成分数
axes[1, 0].plot(results['n_components'], results['train_time'],
'mo-', linewidth=2, markersize=6)
axes[1, 0].axhline(y=train_time_original, color='orange', linestyle='--',
alpha=0.7, label='原始数据训练时间')
axes[1, 0].set_title('训练时间 vs 主成分数')
axes[1, 0].set_xlabel('主成分数')
axes[1, 0].set_ylabel('训练时间 (秒)')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# 准确率 vs 方差比
axes[1, 1].scatter(results['variance_ratio'], results['test_accuracy'],
c=results['n_components'], cmap='viridis', s=100, alpha=0.7)
axes[1, 1].scatter(1.0, test_acc_original, c='red', s=150, marker='*',
label='原始数据')
# 添加颜色条
scatter = axes[1, 1].scatter(results['variance_ratio'], results['test_accuracy'],
c=results['n_components'], cmap='viridis', s=100, alpha=0.7)
plt.colorbar(scatter, ax=axes[1, 1], label='主成分数')
axes[1, 1].set_title('测试准确率 vs 累积方差比')
axes[1, 1].set_xlabel('累积方差比')
axes[1, 1].set_ylabel('测试准确率')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 找出最佳主成分数
best_idx = np.argmax(results['test_accuracy'])
best_n_components = results['n_components'][best_idx]
best_test_acc = results['test_accuracy'][best_idx]
best_var_ratio = results['variance_ratio'][best_idx]
print(f"\n最佳配置:")
print(f"主成分数: {best_n_components}")
print(f"测试准确率: {best_test_acc:.4f}")
print(f"累积方差比: {best_var_ratio:.4f}")
print(f"相比原始数据的性能变化: {(best_test_acc - test_acc_original):.4f}")
return results
classification_results = pca_classification_performance()12.6.2 PCA用于数据可视化
python
def pca_visualization_example():
"""使用PCA进行高维数据可视化"""
# 使用葡萄酒数据集
wine = load_wine()
X_wine = wine.data
y_wine = wine.target
feature_names = wine.feature_names
target_names = wine.target_names
print("使用PCA进行高维数据可视化:")
print(f"原始数据维度: {X_wine.shape[1]}")
print(f"类别数: {len(target_names)}")
# 标准化
scaler = StandardScaler()
X_wine_scaled = scaler.fit_transform(X_wine)
# PCA降维到2D和3D
pca_2d = PCA(n_components=2)
pca_3d = PCA(n_components=3)
X_wine_2d = pca_2d.fit_transform(X_wine_scaled)
X_wine_3d = pca_3d.fit_transform(X_wine_scaled)
print(f"\n2D PCA解释方差比: {pca_2d.explained_variance_ratio_}")
print(f"2D PCA累积方差比: {np.sum(pca_2d.explained_variance_ratio_):.4f}")
print(f"3D PCA累积方差比: {np.sum(pca_3d.explained_variance_ratio_):.4f}")
# 可视化
fig = plt.figure(figsize=(20, 15))
# 原始特征的相关性矩阵
plt.subplot(3, 3, 1)
correlation_matrix = np.corrcoef(X_wine_scaled.T)
sns.heatmap(correlation_matrix, cmap='coolwarm', center=0, square=True,
xticklabels=False, yticklabels=False)
plt.title('原始特征相关性矩阵')
# 主成分载荷图
plt.subplot(3, 3, 2)
loadings = pca_2d.components_.T
plt.scatter(loadings[:, 0], loadings[:, 1], alpha=0.7)
# 标注重要特征
for i, feature in enumerate(feature_names):
if abs(loadings[i, 0]) > 0.3 or abs(loadings[i, 1]) > 0.3:
plt.annotate(feature[:10], (loadings[i, 0], loadings[i, 1]),
xytext=(2, 2), textcoords='offset points', fontsize=8)
plt.xlabel(f'PC1 ({pca_2d.explained_variance_ratio_[0]:.2%})')
plt.ylabel(f'PC2 ({pca_2d.explained_variance_ratio_[1]:.2%})')
plt.title('主成分载荷图')
plt.grid(True, alpha=0.3)
# 2D PCA可视化
plt.subplot(3, 3, 3)
colors = ['red', 'blue', 'green']
for i, (target_name, color) in enumerate(zip(target_names, colors)):
mask = y_wine == i
plt.scatter(X_wine_2d[mask, 0], X_wine_2d[mask, 1],
c=color, alpha=0.7, label=target_name, s=50)
plt.xlabel(f'PC1 ({pca_2d.explained_variance_ratio_[0]:.2%})')
plt.ylabel(f'PC2 ({pca_2d.explained_variance_ratio_[1]:.2%})')
plt.title('2D PCA可视化')
plt.legend()
plt.grid(True, alpha=0.3)
# 3D PCA可视化
ax_3d = fig.add_subplot(3, 3, 4, projection='3d')
for i, (target_name, color) in enumerate(zip(target_names, colors)):
mask = y_wine == i
ax_3d.scatter(X_wine_3d[mask, 0], X_wine_3d[mask, 1], X_wine_3d[mask, 2],
c=color, alpha=0.7, label=target_name, s=30)
ax_3d.set_xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%})')
ax_3d.set_ylabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%})')
ax_3d.set_zlabel(f'PC3 ({pca_3d.explained_variance_ratio_[2]:.2%})')
ax_3d.set_title('3D PCA可视化')
ax_3d.legend()
# 原始特征的分布（选择几个重要特征）
important_features_idx = [0, 6, 9, 12] # 选择几个代表性特征
for idx, feature_idx in enumerate(important_features_idx):
plt.subplot(3, 3, 5 + idx)
for i, (target_name, color) in enumerate(zip(target_names, colors)):
mask = y_wine == i
plt.hist(X_wine_scaled[mask, feature_idx], alpha=0.6,
color=color, label=target_name, bins=15)
plt.xlabel(feature_names[feature_idx][:15])
plt.ylabel('频次')
plt.title(f'特征分布: {feature_names[feature_idx][:15]}')
plt.legend()
plt.grid(True, alpha=0.3)
# 解释方差比
plt.subplot(3, 3, 9)
all_pca = PCA()
all_pca.fit(X_wine_scaled)
plt.bar(range(1, len(all_pca.explained_variance_ratio_) + 1),
all_pca.explained_variance_ratio_, alpha=0.7, color='skyblue')
plt.xlabel('主成分')
plt.ylabel('解释方差比')
plt.title('所有主成分的解释方差比')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 分析主成分的含义
print(f"\n主成分分析:")
print("前两个主成分的主要贡献特征:")
# 获取载荷矩阵
loadings_df = pd.DataFrame(
pca_2d.components_.T,
columns=['PC1', 'PC2'],
index=feature_names
)
print("\nPC1的主要贡献特征:")
pc1_contributions = loadings_df['PC1'].abs().sort_values(ascending=False)
for i, (feature, contribution) in enumerate(pc1_contributions.head(5).items()):
print(f" {i+1}. {feature}: {contribution:.3f}")
print("\nPC2的主要贡献特征:")
pc2_contributions = loadings_df['PC2'].abs().sort_values(ascending=False)
for i, (feature, contribution) in enumerate(pc2_contributions.head(5).items()):
print(f" {i+1}. {feature}: {contribution:.3f}")
return X_wine_2d, X_wine_3d, pca_2d, pca_3d
X_wine_2d, X_wine_3d, pca_2d, pca_3d = pca_visualization_example()12.7 PCA的局限性和注意事项
12.7.1 PCA假设和局限性
python
def pca_limitations_demo():
"""演示PCA的局限性"""
print("PCA的局限性和注意事项:")
print("=" * 30)
# 1. 线性假设的局限性
print("1. 线性假设局限性:")
# 创建非线性数据
np.random.seed(42)
t = np.linspace(0, 4*np.pi, 300)
x1 = t * np.cos(t) + 0.5 * np.random.randn(300)
x2 = t * np.sin(t) + 0.5 * np.random.randn(300)
X_spiral = np.column_stack([x1, x2])
# 标准化
scaler = StandardScaler()
X_spiral_scaled = scaler.fit_transform(X_spiral)
# 应用PCA
pca_spiral = PCA(n_components=2)
X_spiral_pca = pca_spiral.fit_transform(X_spiral_scaled)
# 2. 方差不等于重要性
print("2. 方差不等于重要性:")
# 创建数据：一个特征有高方差但对分类无用
np.random.seed(42)
n_samples = 500
# 有用特征（低方差但有分类信息）
useful_feature = np.random.normal(0, 0.1, n_samples)
# 无用特征（高方差但无分类信息）
useless_feature = np.random.normal(0, 5, n_samples)
# 标签基于有用特征
y_synthetic = (useful_feature > 0).astype(int)
X_synthetic = np.column_stack([useful_feature, useless_feature])
# 标准化
X_synthetic_scaled = scaler.fit_transform(X_synthetic)
# PCA
pca_synthetic = PCA(n_components=2)
X_synthetic_pca = pca_synthetic.fit_transform(X_synthetic_scaled)
# 可视化
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('PCA的局限性演示', fontsize=16)
# 非线性数据
axes[0, 0].scatter(X_spiral_scaled[:, 0], X_spiral_scaled[:, 1],
c=t, cmap='viridis', alpha=0.7, s=20)
axes[0, 0].set_title('非线性螺旋数据')
axes[0, 0].set_xlabel('特征 1')
axes[0, 0].set_ylabel('特征 2')
# PCA后的螺旋数据
axes[0, 1].scatter(X_spiral_pca[:, 0], X_spiral_pca[:, 1],
c=t, cmap='viridis', alpha=0.7, s=20)
axes[0, 1].set_title('PCA后的螺旋数据')
axes[0, 1].set_xlabel('PC1')
axes[0, 1].set_ylabel('PC2')
# 螺旋数据的主成分方向
axes[0, 2].scatter(X_spiral_scaled[:, 0], X_spiral_scaled[:, 1],
c=t, cmap='viridis', alpha=0.5, s=20)
# 绘制主成分方向
mean_point = np.mean(X_spiral_scaled, axis=0)
for i, component in enumerate(pca_spiral.components_):
vector = component * np.sqrt(pca_spiral.explained_variance_[i]) * 2
axes[0, 2].arrow(mean_point[0], mean_point[1], vector[0], vector[1],
head_width=0.1, head_length=0.1, fc=f'C{i}', ec=f'C{i}',
linewidth=3, label=f'PC{i+1}')
axes[0, 2].set_title('主成分方向（螺旋数据）')
axes[0, 2].set_xlabel('特征 1')
axes[0, 2].set_ylabel('特征 2')
axes[0, 2].legend()
# 方差vs重要性问题
axes[1, 0].scatter(X_synthetic_scaled[:, 0], X_synthetic_scaled[:, 1],
c=y_synthetic, cmap='RdYlBu', alpha=0.7, s=30)
axes[1, 0].set_title('原始数据（颜色=类别）')
axes[1, 0].set_xlabel('有用特征（低方差）')
axes[1, 0].set_ylabel('无用特征（高方差）')
# PCA后的数据
axes[1, 1].scatter(X_synthetic_pca[:, 0], X_synthetic_pca[:, 1],
c=y_synthetic, cmap='RdYlBu', alpha=0.7, s=30)
axes[1, 1].set_title('PCA后的数据')
axes[1, 1].set_xlabel(f'PC1 ({pca_synthetic.explained_variance_ratio_[0]:.2%})')
axes[1, 1].set_ylabel(f'PC2 ({pca_synthetic.explained_variance_ratio_[1]:.2%})')
# 特征方差比较
original_var = np.var(X_synthetic_scaled, axis=0)
axes[1, 2].bar(['有用特征', '无用特征'], original_var,
color=['green', 'red'], alpha=0.7)
axes[1, 2].set_title('原始特征方差')
axes[1, 2].set_ylabel('方差')
# 添加数值标签
for i, var in enumerate(original_var):
axes[1, 2].text(i, var + 0.01, f'{var:.3f}', ha='center')
plt.tight_layout()
plt.show()
print(f" 螺旋数据PCA解释方差比: {pca_spiral.explained_variance_ratio_}")
print(f" 合成数据PCA解释方差比: {pca_synthetic.explained_variance_ratio_}")
print(f" 注意：PCA选择了高方差但无用的特征作为第一主成分")
# 3. 可解释性问题
print(f"\n3. 可解释性问题:")
print(f" 主成分是原始特征的线性组合，难以直接解释")
print(f" 例如：PC1 = {pca_synthetic.components_[0][0]:.3f} * 有用特征 + {pca_synthetic.components_[0][1]:.3f} * 无用特征")
return X_spiral_scaled, X_synthetic_scaled, pca_spiral, pca_synthetic
X_spiral_scaled, X_synthetic_scaled, pca_spiral, pca_synthetic = pca_limitations_demo()12.7.2 PCA最佳实践
python
def pca_best_practices():
"""PCA最佳实践指南"""
print("PCA最佳实践指南:")
print("=" * 25)
practices = {
"数据预处理": [
"总是标准化数据（除非特征已经在相同尺度）",
"检查并处理缺失值",
"考虑异常值的影响",
"确保数据质量"
],
"主成分数量选择": [
"使用累积解释方差比（通常85%-95%）",
"应用Kaiser准则（特征值>1）",
"观察碎石图的肘部",
"考虑下游任务的性能"
],
"模型验证": [
"使用交叉验证评估降维效果",
"比较降维前后的模型性能",
"检查主成分的稳定性",
"分析主成分的可解释性"
],
"应用场景": [
"数据可视化（降到2D/3D）",
"噪声去除和数据压缩",
"特征预处理（减少维度诅咒）",
"探索性数据分析"
],
"注意事项": [
"PCA假设线性关系",
"方差大不等于重要性高",
"主成分难以解释",
"对异常值敏感"
]
}
for category, items in practices.items():
print(f"\n{category}:")
for i, item in enumerate(items, 1):
print(f" {i}. {item}")
# 实践示例：完整的PCA工作流
print(f"\n完整的PCA工作流示例:")
print("-" * 30)
# 使用乳腺癌数据集
cancer = load_breast_cancer()
X, y = cancer.data, cancer.target
print(f"1. 数据概览: {X.shape}")
# 2. 数据预处理
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
print(f"2. 数据标准化完成")
# 3. 初步PCA分析
pca_full = PCA()
pca_full.fit(X_scaled)
# 4. 选择主成分数量
cumsum_var = np.cumsum(pca_full.explained_variance_ratio_)
n_components_95 = np.argmax(cumsum_var >= 0.95) + 1
print(f"3. 保留95%方差需要 {n_components_95} 个主成分")
# 5. 应用PCA
pca_final = PCA(n_components=n_components_95)
X_pca = pca_final.fit_transform(X_scaled)
print(f"4. 降维后数据形状: {X_pca.shape}")
# 6. 验证效果
X_train, X_test, y_train, y_test = train_test_split(
X_scaled, y, test_size=0.2, random_state=42
)
X_train_pca, X_test_pca, _, _ = train_test_split(
X_pca, y, test_size=0.2, random_state=42
)
# 原始数据分类
clf_original = LogisticRegression(random_state=42, max_iter=1000)
clf_original.fit(X_train, y_train)
acc_original = clf_original.score(X_test, y_test)
# PCA数据分类
clf_pca = LogisticRegression(random_state=42, max_iter=1000)
clf_pca.fit(X_train_pca, y_train)
acc_pca = clf_pca.score(X_test_pca, y_test)
print(f"5. 性能比较:")
print(f" 原始数据准确率: {acc_original:.4f}")
print(f" PCA数据准确率: {acc_pca:.4f}")
print(f" 维度减少: {X.shape[1]} → {X_pca.shape[1]} ({X_pca.shape[1]/X.shape[1]*100:.1f}%)")
# 可视化工作流结果
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# 累积解释方差比
axes[0].plot(range(1, len(cumsum_var) + 1), cumsum_var, 'bo-', linewidth=2)
axes[0].axhline(y=0.95, color='red', linestyle='--', alpha=0.7, label='95%阈值')
axes[0].axvline(x=n_components_95, color='red', linestyle=':', alpha=0.7)
axes[0].set_title('累积解释方差比')
axes[0].set_xlabel('主成分数量')
axes[0].set_ylabel('累积解释方差比')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 前两个主成分的2D可视化
axes[1].scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='RdYlBu', alpha=0.7)
axes[1].set_title('前两个主成分可视化')
axes[1].set_xlabel(f'PC1 ({pca_final.explained_variance_ratio_[0]:.2%})')
axes[1].set_ylabel(f'PC2 ({pca_final.explained_variance_ratio_[1]:.2%})')
# 性能比较
methods = ['原始数据', 'PCA数据']
accuracies = [acc_original, acc_pca]
colors = ['blue', 'orange']
bars = axes[2].bar(methods, accuracies, color=colors, alpha=0.7)
axes[2].set_title('分类性能比较')
axes[2].set_ylabel('准确率')
axes[2].set_ylim(0.9, 1.0)
# 添加数值标签
for bar, acc in zip(bars, accuracies):
height = bar.get_height()
axes[2].text(bar.get_x() + bar.get_width()/2., height + 0.005,
f'{acc:.4f}', ha='center', va='bottom')
plt.tight_layout()
plt.show()
return X_scaled, X_pca, pca_final
X_scaled, X_pca, pca_final = pca_best_practices()
```#
# 12.8 练习题
### 练习1：基础PCA
1. 使用鸢尾花数据集进行PCA分析
2. 确定保留90%和95%方差所需的主成分数量
3. 可视化前两个主成分，并分析它们的含义
### 练习2：PCA与分类
1. 使用葡萄酒数据集，比较不同主成分数量对分类性能的影响
2. 绘制准确率vs主成分数量的曲线
3. 找出最佳的主成分数量
### 练习3：PCA变体比较
1. 创建一个大数据集（样本数>5000）
2. 比较标准PCA和增量PCA的性能和结果
3. 分析两种方法的优缺点
### 练习4：非线性数据
1. 创建一个非线性数据集（如螺旋形或S形）
2. 比较线性PCA和核PCA的效果
3. 尝试不同的核函数
### 练习5：实际应用
1. 使用手写数字数据集（load_digits）
2. 应用PCA进行降维和可视化
3. 分析降维对数字识别性能的影响
## 12.9 小结
在本章中，我们深入学习了主成分分析的各个方面：
### 核心概念
- **降维技术**：在保持主要信息的同时减少特征维度
- **方差最大化**：寻找数据中方差最大的方向
- **线性变换**：通过正交变换实现降维
### 主要技术
- **标准PCA**：经典的主成分分析方法
- **增量PCA**：适用于大数据集的批量处理
- **核PCA**：使用核技巧处理非线性数据
- **主成分选择**：多种方法确定最佳主成分数量
### 实践技能
- **数据预处理**：标准化的重要性
- **可视化技术**：高维数据的2D/3D投影
- **性能评估**：降维对下游任务的影响
- **参数调优**：主成分数量的选择策略
### 关键要点
- PCA是无监督的线性降维技术
- 数据标准化对PCA结果至关重要
- 主成分数量需要在信息保留和维度减少间平衡
- PCA假设线性关系，对非线性数据效果有限
### 应用场景
**适合使用PCA的情况**：
- 数据维度很高需要降维
- 特征间存在线性相关性
- 需要数据可视化
- 要进行噪声去除
- 存储空间或计算资源有限
**不适合使用PCA的情况**：
- 数据维度本身就很低
- 特征间相关性很弱
- 需要保持特征的可解释性
- 数据存在强非线性关系
- 所有特征都很重要
### 最佳实践总结
1. **数据预处理**
- 总是标准化数据
- 处理缺失值和异常值
- 检查数据质量
2. **主成分选择**
- 使用多种方法确定主成分数量
- 考虑下游任务的性能
- 平衡信息保留和计算效率
3. **结果验证**
- 分析主成分的可解释性
- 验证降维对任务性能的影响
- 检查结果的稳定性
4. **实际应用**
- 根据具体问题选择合适的PCA变体
- 结合领域知识解释主成分
- 考虑PCA的局限性
## 12.10 下一步
现在你已经掌握了主成分分析这个重要的降维技术！在下一章[异常检测](./sklearn-anomaly-detection.md)中，我们将学习如何识别数据中的异常点和离群值，这在数据质量控制和欺诈检测等领域有重要应用。
---
**章节要点回顾**：
- ✅ 理解了PCA的数学原理和几何直觉
- ✅ 掌握了PCA的实现和参数调优
- ✅ 学会了确定最佳主成分数量的方法
- ✅ 了解了PCA的变体和适用场景
- ✅ 掌握了PCA在数据可视化中的应用
- ✅ 认识了PCA的局限性和最佳实践
- ✅ 能够在实际项目中合理使用PCA技术