我将为您详细介绍如何使用Statsmodels进行主成分分析（PCA）。
Statsmodels PCA 基础用法
安装和导入库
# 安装所需库 # pip install statsmodels pandas numpy matplotlib import pandas as pd import numpy as np import matplotlib.pyplot as plt from statsmodels.multivariate.pca import PCA import seaborn as sns # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False
创建示例数据
# 创建示例数据集
np.random.seed(42)
n_samples = 100
# 生成具有相关性的数据
X1 = np.random.normal(0, 1, n_samples)
X2 = X1 * 0.7 + np.random.normal(0, 0.5, n_samples)
X3 = X1 * 0.5 + X2 * 0.3 + np.random.normal(0, 0.3, n_samples)
X4 = np.random.normal(0, 2, n_samples)
X5 = X4 * 0.6 + np.random.normal(0, 0.4, n_samples)
# 创建DataFrame
data = pd.DataFrame({
'变量1': X1,
'变量2': X2,
'变量3': X3,
'变量4': X4,
'变量5': X5
})
print("数据集前5行:")
print(data.head())
print("\n数据集形状:", data.shape)
执行PCA分析
# 方法1：标准化数据后进行PCA
pca_standardized = PCA(data, standardize=True, ncomp=5)
# 方法2：不标准化数据（适用于相同量纲的数据）
# pca_raw = PCA(data, standardize=False, ncomp=5)
print("=== PCA 分析结果 ===")
print(f"解释方差比例:\n{pca_standardized.rsquare}")
print(f"\n累计解释方差比例:\n{np.cumsum(pca_standardized.rsquare)}")
主成分载荷和得分
# 载荷矩阵（主成分与原始变量的关系）
loadings = pca_standardized.loadings
print("主成分载荷矩阵:")
print(pd.DataFrame(
loadings,
index=data.columns,
columns=[f'PC{i+1}' for i in range(loadings.shape[1])]
).round(3))
# 主成分得分
scores = pca_standardized.scores
print("\n前5个样本的主成分得分:")
print(pd.DataFrame(
scores[:5, :],
columns=[f'PC{i+1}' for i in range(scores.shape[1])]
).round(3))
可视化分析结果
# 设置图形样式
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 图1：碎石图（解释方差）
explained_var = pca_standardized.rsquare
axes[0, 0].bar(range(1, len(explained_var)+1), explained_var, alpha=0.7)
axes[0, 0].plot(range(1, len(explained_var)+1), explained_var, 'ro-', linewidth=2)
axes[0, 0].set_xlabel('主成分')
axes[0, 0].set_ylabel('解释方差比例')
axes[0, 0].set_title('碎石图 - 各主成分解释方差')
axes[0, 0].grid(True, alpha=0.3)
# 图2：累计解释方差
cum_var = np.cumsum(explained_var)
axes[0, 1].plot(range(1, len(cum_var)+1), cum_var, 'bo-', linewidth=2)
axes[0, 1].axhline(y=0.8, color='r', linestyle='--', alpha=0.7, label='80% 阈值')
axes[0, 1].axhline(y=0.9, color='g', linestyle='--', alpha=0.7, label='90% 阈值')
axes[0, 1].set_xlabel('主成分数量')
axes[0, 1].set_ylabel('累计解释方差比例')
axes[0, 1].set_title('累计解释方差图')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# 图3：载荷图
loadings_df = pd.DataFrame(loadings[:, :2],
index=data.columns,
columns=['PC1', 'PC2'])
axes[1, 0].scatter(loadings_df['PC1'], loadings_df['PC2'], alpha=0.7)
for i, label in enumerate(loadings_df.index):
axes[1, 0].annotate(label,
(loadings_df['PC1'].iloc[i], loadings_df['PC2'].iloc[i]),
fontsize=9)
axes[1, 0].axhline(y=0, color='gray', linestyle='-', alpha=0.3)
axes[1, 0].axvline(x=0, color='gray', linestyle='-', alpha=0.3)
axes[1, 0].set_xlabel('第一主成分 (PC1)')
axes[1, 0].set_ylabel('第二主成分 (PC2)')
axes[1, 0].set_title('载荷图 - 变量在主成分上的投影')
axes[1, 0].grid(True, alpha=0.3)
# 图4：得分图
axes[1, 1].scatter(scores[:, 0], scores[:, 1], alpha=0.6)
axes[1, 1].set_xlabel('第一主成分得分')
axes[1, 1].set_ylabel('第二主成分得分')
axes[1, 1].set_title('得分图 - 样本在PC1和PC2上的投影')
axes[1, 1].axhline(y=0, color='gray', linestyle='-', alpha=0.3)
axes[1, 1].axvline(x=0, color='gray', linestyle='-', alpha=0.3)
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
实际案例：鸢尾花数据集
from sklearn.datasets import load_iris
# 加载鸢尾花数据集
iris = load_iris()
iris_data = pd.DataFrame(iris.data, columns=iris.feature_names)
iris_target = iris.target
print("鸢尾花数据集示例:")
print(iris_data.head())
# 执行PCA
pca_iris = PCA(iris_data, standardize=True, ncomp=4)
# 可视化结果
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# 碎石图
explained_var_iris = pca_iris.rsquare
axes[0].bar(range(1, len(explained_var_iris)+1), explained_var_iris, alpha=0.7)
axes[0].plot(range(1, len(explained_var_iris)+1), explained_var_iris, 'ro-')
axes[0].set_xlabel('主成分')
axes[0].set_ylabel('解释方差比例')
axes[0].set_title('鸢尾花数据集 - PCA碎石图')
axes[0].grid(True, alpha=0.3)
# 得分图（按类别着色）
scores_iris = pca_iris.scores
scatter = axes[1].scatter(scores_iris[:, 0], scores_iris[:, 1],
c=iris_target, cmap='viridis', alpha=0.7)
axes[1].set_xlabel('第一主成分')
axes[1].set_ylabel('第二主成分')
axes[1].set_title('鸢尾花数据集 - PCA得分图')
plt.colorbar(scatter, ax=axes[1])
plt.tight_layout()
plt.show()
# 显示载荷矩阵
print("\n鸢尾花数据集的载荷矩阵:")
loadings_iris = pca_iris.loadings
loadings_iris_df = pd.DataFrame(
loadings_iris[:, :3],
index=iris.feature_names,
columns=[f'PC{i+1}' for i in range(3)]
)
print(loadings_iris_df.round(3))
选择主成分数量的方法
def select_components(pca_obj, threshold=0.8):
"""
根据累计解释方差比例选择主成分数量
"""
cum_var = np.cumsum(pca_obj.rsquare)
n_components = np.argmax(cum_var >= threshold) + 1
return n_components, cum_var
# 选择主成分
n_comp, cum_var = select_components(pca_standardized)
print(f"解释80%方差需要的主成分数量: {n_comp}")
print(f"各主成分累计解释方差:\n{cum_var}")
# 使用选定的主成分数重新运行PCA
pca_optimized = PCA(data, standardize=True, ncomp=n_comp)
print(f"\n优化后的主成分数量: {pca_optimized.ncomp}")
print(f"解释方差比例: {pca_optimized.rsquare}")
完整函数封装
def complete_pca_analysis(data, standardize=True, threshold=0.8):
"""
完整的PCA分析函数
参数:
data: DataFrame, 输入数据
standardize: bool, 是否标准化
threshold: float, 解释方差阈值
返回:
PCA对象和主要结果
"""
# 执行PCA
pca = PCA(data, standardize=standardize)
# 选择主成分数量
cum_var = np.cumsum(pca.rsquare)
n_optimal = np.argmax(cum_var >= threshold) + 1
# 使用优化后的主成分数重新运行
pca_optimal = PCA(data, standardize=standardize, ncomp=n_optimal)
results = {
'pca_object': pca_optimal,
'optimal_components': n_optimal,
'explained_variance': pca_optimal.rsquare,
'cumulative_variance': cum_var,
'loadings': pca_optimal.loadings,
'scores': pca_optimal.scores,
'sdev': pca_optimal.sdev
}
return results
# 使用示例
results = complete_pca_analysis(data)
print(f"最优主成分数量: {results['optimal_components']}")
print(f"解释方差: {results['explained_variance'].round(4)}")
- 数据标准化：使用
standardize=True处理不同量纲的变量 - 主成分选择：通常选择累计解释方差达到80%-90%的主成分数量
- 解释结果：
- 载荷矩阵显示变量与主成分的关系
- 得分矩阵显示样本在主成分上的投影
- 可视化：碎石图、载荷图、得分图是常用的诊断工具
Statsmodels的PCA实现提供了详细的统计信息,适合需要深入分析的应用场景。
