{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Sparse inverse covariance estimation\n\n\nUsing the GraphicalLasso estimator to learn a covariance and sparse precision\nfrom a small number of samples.\n\nTo estimate a probabilistic model (e.g. a Gaussian model), estimating the\nprecision matrix, that is the inverse covariance matrix, is as important\nas estimating the covariance matrix. Indeed a Gaussian model is\nparametrized by the precision matrix.\n\nTo be in favorable recovery conditions, we sample the data from a model\nwith a sparse inverse covariance matrix. In addition, we ensure that the\ndata is not too much correlated (limiting the largest coefficient of the\nprecision matrix) and that there a no small coefficients in the\nprecision matrix that cannot be recovered. In addition, with a small\nnumber of observations, it is easier to recover a correlation matrix\nrather than a covariance, thus we scale the time series.\n\nHere, the number of samples is slightly larger than the number of\ndimensions, thus the empirical covariance is still invertible. However,\nas the observations are strongly correlated, the empirical covariance\nmatrix is ill-conditioned and as a result its inverse --the empirical\nprecision matrix-- is very far from the ground truth.\n\nIf we use l2 shrinkage, as with the Ledoit-Wolf estimator, as the number\nof samples is small, we need to shrink a lot. As a result, the\nLedoit-Wolf precision is fairly close to the ground truth precision, that\nis not far from being diagonal, but the off-diagonal structure is lost.\n\nThe l1-penalized estimator can recover part of this off-diagonal\nstructure. It learns a sparse precision. It is not able to\nrecover the exact sparsity pattern: it detects too many non-zero\ncoefficients. However, the highest non-zero coefficients of the l1\nestimated correspond to the non-zero coefficients in the ground truth.\nFinally, the coefficients of the l1 precision estimate are biased toward\nzero: because of the penalty, they are all smaller than the corresponding\nground truth value, as can be seen on the figure.\n\nNote that, the color range of the precision matrices is tweaked to\nimprove readability of the figure. The full range of values of the\nempirical precision is not displayed.\n\nThe alpha parameter of the GraphicalLasso setting the sparsity of the model is\nset by internal cross-validation in the GraphicalLassoCV. As can be\nseen on figure 2, the grid to compute the cross-validation score is\niteratively refined in the neighborhood of the maximum.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(__doc__)\n# author: Gael Varoquaux <gael.varoquaux@inria.fr>\n# License: BSD 3 clause\n# Copyright: INRIA\n\nimport numpy as np\nfrom scipy import linalg\nfrom sklearn.datasets import make_sparse_spd_matrix\nfrom sklearn.covariance import GraphicalLassoCV, ledoit_wolf\nimport matplotlib.pyplot as plt\n\n# #############################################################################\n# Generate the data\nn_samples = 60\nn_features = 20\n\nprng = np.random.RandomState(1)\nprec = make_sparse_spd_matrix(n_features, alpha=.98,\n                              smallest_coef=.4,\n                              largest_coef=.7,\n                              random_state=prng)\ncov = linalg.inv(prec)\nd = np.sqrt(np.diag(cov))\ncov /= d\ncov /= d[:, np.newaxis]\nprec *= d\nprec *= d[:, np.newaxis]\nX = prng.multivariate_normal(np.zeros(n_features), cov, size=n_samples)\nX -= X.mean(axis=0)\nX /= X.std(axis=0)\n\n# #############################################################################\n# Estimate the covariance\nemp_cov = np.dot(X.T, X) / n_samples\n\nmodel = GraphicalLassoCV(cv=5)\nmodel.fit(X)\ncov_ = model.covariance_\nprec_ = model.precision_\n\nlw_cov_, _ = ledoit_wolf(X)\nlw_prec_ = linalg.inv(lw_cov_)\n\n# #############################################################################\n# Plot the results\nplt.figure(figsize=(10, 6))\nplt.subplots_adjust(left=0.02, right=0.98)\n\n# plot the covariances\ncovs = [('Empirical', emp_cov), ('Ledoit-Wolf', lw_cov_),\n        ('GraphicalLassoCV', cov_), ('True', cov)]\nvmax = cov_.max()\nfor i, (name, this_cov) in enumerate(covs):\n    plt.subplot(2, 4, i + 1)\n    plt.imshow(this_cov, interpolation='nearest', vmin=-vmax, vmax=vmax,\n               cmap=plt.cm.RdBu_r)\n    plt.xticks(())\n    plt.yticks(())\n    plt.title('%s covariance' % name)\n\n\n# plot the precisions\nprecs = [('Empirical', linalg.inv(emp_cov)), ('Ledoit-Wolf', lw_prec_),\n         ('GraphicalLasso', prec_), ('True', prec)]\nvmax = .9 * prec_.max()\nfor i, (name, this_prec) in enumerate(precs):\n    ax = plt.subplot(2, 4, i + 5)\n    plt.imshow(np.ma.masked_equal(this_prec, 0),\n               interpolation='nearest', vmin=-vmax, vmax=vmax,\n               cmap=plt.cm.RdBu_r)\n    plt.xticks(())\n    plt.yticks(())\n    plt.title('%s precision' % name)\n    if hasattr(ax, 'set_facecolor'):\n        ax.set_facecolor('.7')\n    else:\n        ax.set_axis_bgcolor('.7')\n\n# plot the model selection metric\nplt.figure(figsize=(4, 3))\nplt.axes([.2, .15, .75, .7])\nplt.plot(model.cv_alphas_, np.mean(model.grid_scores_, axis=1), 'o-')\nplt.axvline(model.alpha_, color='.5')\nplt.title('Model selection')\nplt.ylabel('Cross-validation score')\nplt.xlabel('alpha')\n\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}