在统计学中，协方差矩阵的估计是一个重要的问题。通常，会使用最大似然估计方法，例如经验协方差矩阵（EmpiricalCovariance），它在样本数量足够多时能够无偏地逼近真实（总体）协方差矩阵。然而，为了降低其方差，引入一定的偏差进行正则化也是有益的。本文将探讨在正则化协方差估计中如何设置正则化强度，即如何选择偏差-方差权衡。
首先，生成一些样本数据来模拟实际情况。使用NumPy库生成正态分布的样本数据，并使用一个随机矩阵对样本进行“着色”，以模拟协方差结构。
import numpy as np
n_features, n_samples = 40, 20
np.random.seed(42)
base_X_train = np.random.normal(size=(n_samples, n_features))
base_X_test = np.random.normal(size=(n_samples, n_features))
coloring_matrix = np.random.normal(size=(n_features, n_features))
X_train = np.dot(base_X_train, coloring_matrix)
X_test = np.dot(base_X_test, coloring_matrix)
接下来，计算测试数据上的似然度。这里使用SciPy库中的线性代数模块和scikit-learn库中的协方差模块来计算经验协方差矩阵和似然度。
from scipy import linalg
from sklearn.covariance import ShrunkCovariance, empirical_covariance, log_likelihood
shrinkages = np.logspace(-2, 0, 30)
negative_logliks = [-ShrunkCovariance(shrinkage=s).fit(X_train).score(X_test) for s in shrinkages]
real_cov = np.dot(coloring_matrix.T, coloring_matrix)
emp_cov = empirical_covariance(X_train)
loglik_real = -log_likelihood(emp_cov, linalg.inv(real_cov))
在比较不同的正则化参数设置方法时，考虑了三种方法：通过三折交叉验证似然度来设置参数，使用Ledoit和Wolf提出的公式计算渐近最优的正则化参数（最小化均方误差标准），以及Chen等人提出的OAS改进方法，它在假设数据为高斯分布的情况下，尤其是对于小样本，收敛性更好。
from sklearn.covariance import OAS, LedoitWolf
from sklearn.model_selection import GridSearchCV
tuned_parameters = [{"shrinkage": shrinkages}]
cv = GridSearchCV(ShrunkCovariance(), tuned_parameters)
cv.fit(X_train)
lw = LedoitWolf()
loglik_lw = lw.fit(X_train).score(X_test)
oa = OAS()
loglik_oa = oa.fit(X_train).score(X_test)
为了量化估计误差，绘制了不同正则化参数值下未见数据的似然度，并展示了交叉验证或LedoitWolf和OAS估计的选择。
import matplotlib.pyplot as plt
fig = plt.figure()
plt.title("正则化协方差：似然度和正则化系数")
plt.xlabel("正则化参数：正则化系数")
plt.ylabel("误差：测试数据上的负对数似然度")
plt.loglog(shrinkages, negative_logliks, label="负对数似然度")
plt.plot(plt.xlim(), 2*[loglik_real], "--r", label="真实协方差似然度")
# 绘制Ledoit-Wolf和OAS估计的垂直线
plt.vlines(lw.shrinkage_, ymin, -loglik_lw, color="magenta", linewidth=3, label="Ledoit-Wolf估计")
plt.vlines(oa.shrinkage_, ymin, -loglik_oa, color="purple", linewidth=3, label="OAS估计")
plt.vlines(cv.best_estimator_.shrinkage, ymin, -cv.best_estimator_.score(X_test), color="cyan", linewidth=3, label="交叉验证最佳估计")
plt.legend()
plt.show()
注意，最大似然估计对应于无正则化，因此表现不佳。Ledoit-Wolf估计表现非常好，因为它接近最优且计算成本不高。在这个例子中，OAS估计稍远一些。有趣的是，这两种方法都优于交叉验证，后者的计算成本要高得多。