本文目录导读：
我来详细讲解如何使用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中协方差估计的多种方法和实际应用场景,根据具体需求选择合适的方法非常重要。
