在概率模型的估计过程中，如高斯模型，精度矩阵（即协方差矩阵的逆）的估计与协方差矩阵的估计同样重要。实际上，高斯模型就是由精度矩阵参数化的。为了获得有利的恢复条件，从具有稀疏逆协方差矩阵的模型中采样数据。此外，确保数据之间的相关性不会过高（限制精度矩阵中的最大系数），并且精度矩阵中没有太小而无法恢复的系数。此外，由于观测值数量较少，恢复相关性矩阵比恢复协方差矩阵更为容易，因此对时间序列进行了缩放。
在本例中，样本数量略多于维度数量，因此经验协方差仍然是可逆的。然而，由于观测值高度相关，经验协方差矩阵的条件数很差，因此其逆——经验精度矩阵——与真实值相差甚远。如果使用Ledoit-Wolf估计器这样的l2收缩方法，由于样本数量较少，需要进行大量的收缩。结果，Ledoit-Wolf精度与真实精度相当接近，接近对角线，但失去了非对角线结构。l1惩罚估计器可以恢复部分非对角线结构。它学习了一个稀疏的精度。它无法恢复确切的稀疏模式：它检测到太多的非零系数。然而，l1估计的最高非零系数对应于真实值中的非零系数。最后，l1精度估计的系数偏向于零：由于惩罚，它们都小于相应的真实值，如图中所示。
请注意，为了提高图形的可读性，精度矩阵的颜色范围已经调整。经验精度的全部值范围没有显示。GraphicalLasso中设置模型稀疏性的alpha参数是通过GraphicalLassoCV中的内部交叉验证设置的。如图2所示，计算交叉验证分数的网格在最大值附近进行了迭代细化。
import numpy as np
from scipy import linalg
from sklearn.datasets import make_sparse_spd_matrix
n_samples = 60
n_features = 20
prng = np.random.RandomState(1)
prec = make_sparse_spd_matrix(n_features, alpha=0.98, smallest_coef=0.4, largest_coef=0.7, random_state=prng)
cov = linalg.inv(prec)
d = np.sqrt(np.diag(cov))
cov /= d
cov /= d[:, np.newaxis]
prec *= d
prec *= d[:, np.newaxis]
X = prng.multivariate_normal(np.zeros(n_features), cov, size=n_samples)
X -= X.mean(axis=0)
X /= X.std(axis=0)
from sklearn.covariance import GraphicalLassoCV, ledoit_wolf
emp_cov = np.dot(X.T, X) / n_samples
model = GraphicalLassoCV()
model.fit(X)
cov_ = model.covariance_
prec_ = model.precision_
lw_cov_, _ = ledoit_wolf(X)
lw_prec_ = linalg.inv(lw_cov_)
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.subplots_adjust(left=0.02, right=0.98)
# plot the covariances
covs = [
("Empirical", emp_cov),
("Ledoit-Wolf", lw_cov_),
("GraphicalLassoCV", cov_),
("True", cov),
]
vmax = cov_.max()
for i, (name, this_cov) in enumerate(covs):
plt.subplot(2, 4, i + 1)
plt.imshow(this_cov, interpolation="nearest", vmin=-vmax, vmax=vmax, cmap=plt.cm.RdBu_r)
plt.xticks(())
plt.yticks(())
plt.title("%s covariance" % name)
# plot the precisions
precs = [
("Empirical", linalg.inv(emp_cov)),
("Ledoit-Wolf", lw_prec_),
("GraphicalLasso", prec_),
("True", prec),
]
vmax = 0.9 * prec_.max()
for i, (name, this_prec) in enumerate(precs):
ax = plt.subplot(2, 4, i + 5)
plt.imshow(np.ma.masked_equal(this_prec, 0), interpolation="nearest", vmin=-vmax, vmax=vmax, cmap=plt.cm.RdBu_r)
plt.xticks(())
plt.yticks(())
plt.title("%s precision" % name)
if hasattr(ax, "set_facecolor"):
ax.set_facecolor(".7")
else:
ax.set_axis_bgcolor(".7")
# plot the model selection metric
plt.figure(figsize=(4, 3))
plt.axes([0.2, 0.15, 0.75, 0.7])
plt.plot(model.cv_results_["alphas"], model.cv_results_["mean_test_score"], "o-")
plt.axvline(model.alpha_, color=".5")
plt.title("Model selection")
plt.ylabel("Cross-validation score")
plt.xlabel("alpha")
plt.show()