Python案例如何用Scikit-learn做协方差估计

wen python案例 2026-07-16 251

本文目录导读：

基础协方差估计
稳健协方差估计
收缩协方差估计
实际案例：金融数据协方差分析
异常检测应用
性能评估与选择
关键要点

我来详细讲解如何使用Scikit-learn进行协方差估计,包含多个实际案例。

基础协方差估计

import numpy as np
import pandas as pd
from sklearn.covariance import EmpiricalCovariance
import matplotlib.pyplot as plt
import seaborn as sns
# 生成示例数据
np.random.seed(42)
n_samples = 100
n_features = 3
# 创建具有相关性的数据
mean = [0, 0, 0]
cov = [[1, 0.8, 0.3],
       [0.8, 1, 0.2],
       [0.3, 0.2, 1]]
X = np.random.multivariate_normal(mean, cov, n_samples)
df = pd.DataFrame(X, columns=['特征1', '特征2', '特征3'])
# 基本协方差估计
cov_estimator = EmpiricalCovariance().fit(X)
cov_matrix = cov_estimator.covariance_
print("协方差矩阵：")
print(np.round(cov_matrix, 3))

稳健协方差估计

from sklearn.covariance import MinCovDet, EllipticEnvelope
# 添加异常值
X_outliers = np.copy(X)
X_outliers[:10] = X_outliers[:10] + 5  # 添加10个异常值
# 比较不同方法
methods = {
    '经验协方差': EmpiricalCovariance(),
    'MCD稳健估计': MinCovDet(support_fraction=0.8),
}
# 可视化对比
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for idx, (name, method) in enumerate(methods.items()):
    method.fit(X_outliers)
    cov = method.covariance_
    # 计算马氏距离
    if name == '经验协方差':
        mahal = method.mahalanobis(X_outliers)
    else:
        mahal = method.mahalanobis(X_outliers)
    axes[idx].scatter(X_outliers[:, 0], X_outliers[:, 1], 
                     c=mahal, cmap='viridis', s=50)
    axes[idx].set_title(f'{name}\n协方差矩阵迹: {np.trace(cov):.2f}')
    axes[idx].set_xlabel('特征1')
    axes[idx].set_ylabel('特征2')
plt.tight_layout()
plt.show()

收缩协方差估计

from sklearn.covariance import ShrunkCovariance, LedoitWolf, OAS
# 生成高维数据
np.random.seed(42)
n_samples_small = 50
n_features_big = 100
X_high_dim = np.random.randn(n_samples_small, n_features_big)
# 比较不同收缩方法
shrinkage_methods = {
    '经验协方差': EmpiricalCovariance(),
    '收缩协方差(α=0.1)': ShrunkCovariance(shrinkage=0.1),
    'Ledoit-Wolf': LedoitWolf(),
    'Oracle近似(OAS)': OAS(),
}
print("不同收缩方法的协方差矩阵迹对比：")
for name, method in shrinkage_methods.items():
    method.fit(X_high_dim)
    cov = method.covariance_
    # 计算特征值分布
    eigenvalues = np.linalg.eigvalsh(cov)
    condition_number = eigenvalues[-1] / eigenvalues[0]
    print(f"{name:30s}: 迹={np.trace(cov):.2f}, "
          f"条件数={condition_number:.2f}, "
          f"最小特征值={eigenvalues[0]:.4f}")

实际案例：金融数据协方差分析

import yfinance as yf
from datetime import datetime, timedelta
# 下载股票数据
tickers = ['AAPL', 'GOOGL', 'MSFT', 'AMZN', 'META']
end_date = datetime.now()
start_date = end_date - timedelta(days=365)
# 获取股票价格数据
stock_data = yf.download(tickers, start=start_date, end=end_date)['Close']
# 计算日收益率
returns = stock_data.pct_change().dropna()
print("收益率数据形状:", returns.shape)
# 协方差估计
estimators = {
    '经验协方差': EmpiricalCovariance(),
    'Ledoit-Wolf': LedoitWolf(),
    'MCD稳健估计': MinCovDet(support_fraction=0.75),
}
# 可视化协方差矩阵
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 原始收益率相关性矩阵
sns.heatmap(returns.corr(), annot=True, fmt='.2f', 
            cmap='coolwarm', ax=axes[0, 0])
axes[0, 0].set_title('收益率相关系数矩阵')
# 不同方法的协方差估计
for ax, (name, estimator) in zip(axes.flatten()[1:], estimators.items()):
    estimator.fit(returns.values)
    cov_mat = estimator.covariance_
    # 转换为相关系数矩阵便于比较
    std = np.sqrt(np.diag(cov_mat))
    corr_mat = cov_mat / np.outer(std, std)
    im = ax.imshow(corr_mat, cmap='coolwarm', vmin=-1, vmax=1)
    ax.set_xticks(range(len(tickers)))
    ax.set_yticks(range(len(tickers)))
    ax.set_xticklabels(tickers, rotation=45)
    ax.set_yticklabels(tickers)
    ax.set_title(f'{name}')
    # 添加数值标签
    for i in range(len(tickers)):
        for j in range(len(tickers)):
            ax.text(j, i, f'{corr_mat[i,j]:.2f}', 
                   ha='center', va='center')
plt.tight_layout()
plt.show()
# 计算投资组合风险
def portfolio_risk(weights, cov_matrix):
    """计算投资组合风险（方差）"""
    return np.dot(weights.T, np.dot(cov_matrix, weights))
# 等权投资组合
n_assets = len(tickers)
equal_weights = np.ones(n_assets) / n_assets
print("\n等权投资组合风险对比：")
for name, estimator in estimators.items():
    estimator.fit(returns.values)
    risk = portfolio_risk(equal_weights, estimator.covariance_)
    print(f"{name:20s}: 年化波动率 = {np.sqrt(252 * risk):.4f}")

异常检测应用

from sklearn.covariance import EllipticEnvelope
# 生成混合数据
np.random.seed(42)
X_normal = np.random.randn(200, 2) * 2
X_outliers = np.random.uniform(low=-6, high=6, size=(20, 2))
X_mixed = np.vstack([X_normal, X_outliers])
# 应用协方差估计进行异常检测
cov_estimator = EllipticEnvelope(contamination=0.1, random_state=42)
labels = cov_estimator.fit_predict(X_mixed)
# 可视化异常检测结果
plt.figure(figsize=(10, 6))
colors = np.array(['blue', 'red'])
plt.scatter(X_mixed[:, 0], X_mixed[:, 1], 
           c=colors[(labels + 1) // 2], 
           s=50, alpha=0.7)
# 绘制协方差椭圆
from matplotlib.patches import Ellipse
import matplotlib.transforms as transforms
def plot_ellipse(ax, mean, cov, color, alpha=0.3):
    """绘制协方差椭圆"""
    eigenvalues, eigenvectors = np.linalg.eigh(cov)
    angle = np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0])
    angle = np.degrees(angle)
    width, height = 2 * np.sqrt(eigenvalues)
    ellipse = Ellipse(xy=mean, width=width, height=height, 
                     angle=angle, facecolor=color, 
                     alpha=alpha, edgecolor='black')
    ax.add_patch(ellipse)
plot_ellipse(plt.gca(), 
             cov_estimator.location_, 
             cov_estimator.covariance_,
             'green', 0.2)
'基于协方差估计的异常检测')
plt.xlabel('特征1')
plt.ylabel('特征2')
plt.legend(['正常样本', '异常样本', '协方差椭圆'])
plt.grid(True, alpha=0.3)
plt.show()

性能评估与选择

from sklearn.model_selection import cross_val_score
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
# 生成受污染的数据
np.random.seed(42)
X_clean = np.random.randn(300, 10)
X_contaminated = X_clean.copy()
# 添加异常值
outlier_indices = np.random.choice(300, 30, replace=False)
X_contaminated[outlier_indices] = X_contaminated[outlier_indices] + 8
# 分割数据
X_train, X_test = train_test_split(X_contaminated, test_size=0.3, random_state=42)
# 评估不同估计器
estimators_to_test = {
    'Empirical': EmpiricalCovariance(),
    'Shrunk(0.1)': ShrunkCovariance(shrinkage=0.1),
    'Shrunk(0.5)': ShrunkCovariance(shrinkage=0.5),
    'Ledoit-Wolf': LedoitWolf(),
    'OAS': OAS(),
    'MCD': MinCovDet(support_fraction=0.75)
}
results = {}
for name, estimator in estimators_to_test.items():
    try:
        estimator.fit(X_train)
        # 计算协方差矩阵估计误差
        true_cov = np.cov(X_clean, rowvar=False)
        estimated_cov = estimator.covariance_
        # Frobenius范数误差
        error = np.linalg.norm(estimated_cov - true_cov, 'fro')
        results[name] = {
            'error': error,
            'condition_number': np.linalg.cond(estimated_cov)
        }
    except Exception as e:
        print(f"{name}: 估算失败 - {e}")
# 显示结果
print("\n协方差估计性能对比：")
results_df = pd.DataFrame(results).T
results_df = results_df.round(3)
print(results_df)
# 可视化结果
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# 误差对比
colors = ['skyblue', 'lightgreen', 'lightgreen', 'lightcoral', 'lightcoral', 'gold']
ax1.bar(results_df.index, results_df['error'], color=colors[:len(results_df)])
ax1.set_title('Frobenius范数估计误差')
ax1.set_xticklabels(results_df.index, rotation=45)
ax1.set_ylabel('误差')
# 条件数对比
ax2.bar(results_df.index, np.log10(results_df['condition_number']), color=colors[:len(results_df)])
ax2.set_title('条件数(log10)')
ax2.set_xticklabels(results_df.index, rotation=45)
ax2.set_ylabel('log10(条件数)')
ax2.axhline(y=np.median(np.log10(results_df['condition_number'])), 
           color='red', linestyle='--', label='中位数')
ax2.legend()
plt.tight_layout()
plt.show()

关键要点

选择建议：
- 数据干净：使用经验协方差
- 高维数据：使用Ledoit-Wolf或OAS
- 存在异常值：使用MCD稳健估计
- 需要正则化：使用收缩协方差
实际应用：
- 资产组合优化
- 异常检测
- 特征选择
- 降维预处理
注意事项：
- 小样本高维数据要谨慎
- 异常值会严重影响估计
- 不同的应用场景选择不同的方法

这些案例展示了Scikit-learn中协方差估计的多种方法和实际应用场景,根据具体需求选择合适的方法非常重要。

本文地址： https://www.bmjup.com/post/22501.html