稀疏逆协方差估计#

使用 GraphicalLasso 估计器从少量样本中学习协方差和稀疏精度矩阵。

要估计概率模型（例如高斯模型），估计精度矩阵（即逆协方差矩阵）与估计协方差矩阵同样重要。实际上，高斯模型是由精度矩阵参数化的。

为了处于有利的恢复条件，我们从一个具有稀疏逆协方差矩阵的模型中抽取数据。此外，我们确保数据相关性不过高（限制精度矩阵的最大系数），并且精度矩阵中没有无法恢复的小系数。另外，在观测数量较少的情况下，恢复相关矩阵比恢复协方差矩阵更容易，因此我们对时间序列进行了缩放。

在此，样本数量略大于维度数量，因此经验协方差矩阵仍然可逆。然而，由于观测值强相关，经验协方差矩阵是病态的，导致其逆矩阵——经验精度矩阵——与真实值相去甚远。

如果使用 l2 收缩（例如 Ledoit-Wolf 估计器），由于样本数量较少，我们需要进行大量收缩。因此，Ledoit-Wolf 精度与真实精度相当接近，即离对角线不远，但其非对角线结构丢失了。

l1 惩罚估计器可以恢复部分非对角线结构。它学习稀疏精度矩阵。它无法恢复精确的稀疏模式：它检测到过多的非零系数。然而，l1 估计器中最高的非零系数与真实值中的非零系数相对应。最后，l1 精度估计的系数偏向于零：由于惩罚的存在，它们都比相应的真实值小，这在图中可以看出来。

请注意，精度矩阵的颜色范围经过调整，以提高图形的可读性。经验精度矩阵的完整值范围未显示。

GraphicalLasso 中设置模型稀疏性的 alpha 参数由 GraphicalLassoCV 的内部交叉验证设置。如图 2 所示，用于计算交叉验证分数的网格在最大值附近进行迭代细化。

# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause

生成数据#

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")

Empirical covariance, Ledoit-Wolf covariance, GraphicalLassoCV covariance, True covariance, Empirical precision, Ledoit-Wolf precision, GraphicalLasso precision, True precision

# 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()

脚本总运行时间： (0 分 0.419 秒)

稀疏逆协方差估计#

生成数据#

估计协方差#

绘制结果#

本页