Python案例如何用Statsmodels做主成分

wen python案例 2026-07-16 179

我将为您详细介绍如何使用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)}")