在统计学中，估计协方差矩阵是描述数据集中变量之间关系的常见方法。通常，使用最大似然估计器，比如经验协方差矩阵，来完成这项任务。经验协方差矩阵是无偏的，也就是说，当样本数量足够多时，它能够收敛到真实的总体协方差。然而，在某些情况下，为了降低其方差，可能会考虑对其进行正则化处理，尽管这会引入一定的偏差。本文将探讨收缩协方差估计器中使用的简单正则化方法，并重点讨论如何设置正则化的程度，即如何选择偏差与方差之间的权衡。
首先，需要生成一些样本数据。在下面的代码中，使用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库中的linalg模块和sklearn库中的covariance模块来计算测试数据上的似然度。将尝试一系列可能的收缩系数值，并计算每个值对应的负对数似然度。
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("Regularized covariance: likelihood and shrinkage coefficient")
plt.xlabel("Regularization parameter: shrinkage coefficient")
plt.ylabel("Error: negative log-likelihood on test data")
plt.loglog(shrinkages, negative_logliks, label="Negative log-likelihood")
plt.plot(plt.xlim(), 2*[loglik_real], "--r", label="Real covariance likelihood")
plt.vlines(lw.shrinkage_, ymin, -loglik_lw, color="magenta", linewidth=3, label="Ledoit-Wolf estimate")
plt.vlines(oa.shrinkage_, ymin, -loglik_oa, color="purple", linewidth=3, label="OAS estimate")
plt.vlines(cv.best_estimator_.shrinkage, ymin, -cv.best_estimator_.score(X_test), color="cyan", linewidth=3, label="Cross-validation best estimate")
plt.legend()
plt.show()
需要注意的是，最大似然估计对应于无收缩，因此表现不佳。Ledoit-Wolf估计表现非常好，因为它接近最优且计算成本不高。在这个例子中，OAS估计稍远一些。有趣的是，这两种方法都优于交叉验证，后者的计算成本显著更高。