IntroCovReg.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:00d94ccc98e502a03283196c0ef131edf4751266f832f051b6c1fb8e890cc715"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import sklearn as skl\n",
      "import scipy.linalg as lin\n",
      "# % run ../utils/stats\n",
      "from __future__ import print_function\n",
      "import matplotlib.pylab as plt\n",
      "from sklearn.datasets import make_sparse_spd_matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1,
     "trusted": true
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The problem\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Some utility functions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# construct data with a certain covariance matrix\n",
      "\n",
      "def cov2corr(C):\n",
      "    \"\"\"\n",
      "    no check that it is symetric positive\n",
      "    \"\"\"\n",
      "    C = np.asarray(C)\n",
      "    D = 1./np.sqrt(np.diag(C))\n",
      "    C = C*D[np.newaxis]\n",
      "    C = (C.T*D[np.newaxis])\n",
      "    return C \n",
      "\n",
      "\n",
      "def frob_norm(C, scaling=True):\n",
      "    \"\"\"\n",
      "    \"\"\"\n",
      "    return (C**2).sum()/np.min(np.asarray(C.shape))\n",
      "\n",
      "def Cex_vs_Cth(N,T,compare='max'):\n",
      "    \"\"\"\n",
      "    \"\"\"\n",
      "    x = np.random.normal(0,1,size=(N,T))\n",
      "    A = np.random.normal(0,1,size=(N,N))\n",
      "    y = A.dot(x)\n",
      "    Cexp = np.cov(y)\n",
      "    Ctheo = A.dot(A.T)\n",
      "    # np.max(np.abs(Cexp-Ctheo)))\n",
      "    if compare=='max':\n",
      "        return(np.mean(np.abs(cov2corr(Cexp)-cov2corr(Ctheo))))\n",
      "\n",
      "# check \n",
      "C = [[1.0, 1.0,  8.1],\n",
      "     [1.0, 16.0, 18.0],\n",
      "     [8.1, 18.0, 81.0]];\n",
      "\n",
      "np.set_printoptions(precision=3)\n",
      "print(cov2corr(C))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 1.    0.25  0.9 ]\n",
        " [ 0.25  1.    0.5 ]\n",
        " [ 0.9   0.5   1.  ]]\n"
       ]
      }
     ],
     "prompt_number": 2,
     "trusted": true
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Two cases: we can estimate a non degenerate covariance - or we cannot"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 10\n",
      "T = 20\n",
      "x = np.random.normal(0,1,size=(N,T))\n",
      "A = np.random.normal(0,1,size=(N,N))\n",
      "y = A.dot(x)\n",
      "Cexp = np.cov(y)\n",
      "print(\"np.linalg.matrix_rank(Cexp) {}\".format(np.linalg.matrix_rank(Cexp)))\n",
      "Ctheo = A.dot(A.T)\n",
      "print(np.max(np.abs(Cexp-Ctheo)))\n",
      "print(np.max(np.abs(cov2corr(Cexp)-cov2corr(Ctheo))))\n",
      "print(frob_norm(cov2corr(Cexp)-cov2corr(Ctheo)))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "np.linalg.matrix_rank(Cexp) 10\n",
        "8.90320197008\n",
        "0.464413039299\n",
        "0.350669182079\n"
       ]
      }
     ],
     "prompt_number": 3,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 20\n",
      "T = np.arange(10,500,10)\n",
      "Rep = 1\n",
      "res = [ np.mean(np.asarray([Cex_vs_Cth(N,t) for i in range(Rep)])) \n",
      "                                                   for t in T]\n",
      "line = plt.plot(res)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD9CAYAAABdoNd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1YVFXiB/DvKOOa7y8h4gw5CaOgwsiGkKtuZBLmJqb9\nWim11lyX+mVZT2+221PqbiqVlWW/NFfNVVfJtoQtxFdmK0hHDURFEhJ0GFHJRBFfYOD8/jg7QQTD\nDAzMMPf7eZ77xMy9586ZW33vmXPPPVclhAAREXm3Du6uABERtT6GPRGRAjDsiYgUgGFPRKQADHsi\nIgVg2BMRKUCTYZ+WljYhODg4T6/X5ycmJr5Yf/2mTZumGwyGw2FhYTmjR4/OyMnJCbOt0+l0RWFh\nYTnh4eFZkZGRJldXnoiIHCSEaHSxWq0dAwMDCwoLC3WVlZVqg8GQnZubG1J3m8zMzFFlZWU9hRDY\nvn37hKioqH22dTqdrvDChQt97H0GFy5cuHBp/cVuy95kMkUGBQUV6HS6IrVaXRUfH78lOTl5ct1t\nRo0a9U3Pnj0vAUBUVNT+4uJibb2Ticr1pygiInKG3bC3WCyagIAAs+21Vqsttlgsmsa2X7NmzeyJ\nEyem2l6rVCoxfvz43REREQdXr149xzVVJiIiZ/nYW6lSqRyeSyE9Pf3OtWvXPpqRkTHa9l5GRsZo\nf3//ktLSUt+YmJhdwcHBeWPHjv2quZ9BRES1nOk5sduy12g0FrPZHGB7bTabA7RabXH97XJycsLm\nzJmzOiUlJa53794Xbe/7+/uXAICvr2/plClTPjOZTJGNVJiLEHj11VfdXgdPWXgseCx4LOwvzrIb\n9hEREQfz8/P1RUVFusrKyk5JSUnT4uLiUupuc/r06VumTp366caNG2cEBQUV2N6/evVql/Ly8u4A\nUFFR0XXnzp13h4aGHnG6hkRE1GJ2u3F8fHysK1asmBsbG7ujurq64+zZs9eEhIQcX7VqVQIAJCQk\nrFq0aNErFy9e7P34449/AABqtbrKZDJFnj17tv/UqVM/BQCr1eozffr0TXfffffO1v9KRERUn6o5\nPwdcWgGVSri7Dp7CaDQiOjra3dXwCDwWtXgsavFY1FKpVBBO9Nkz7ImI2iFnw57TJRARKQDDnohI\nARj2REQKwLAnIlIAhj0RkQIw7ImIFIBhT0SkAAx7IiIFYNgTESmAR4T9tWvurgERkXfziLAvLHR3\nDYiIvJtHhH1BQdPbEBFR8zHsiYgUwCPCPj/f3TUgIvJuHhH2bNkTEbUuhj0RkQJ4RNiXlAA3bri7\nFkRE3ssjwj4gACgqcnctiIi8l0eEfVAQu3KIiFqTx4Q9R+QQEbUejwl7tuyJiFoPw56ISAEY9kRE\nCqASQri3AiqVuHFDoEcPoLwcUKvdWh0ionZBpVJBCKFydHuPaNl36gQMGACcOuXumhAReSePCHuA\nI3KIiFqTR4U9++2JiFoHw56ISAEY9kRECuAxYa/XM+yJiFqLRwy9FELg+nWgVy/gyhXAx8etVSIi\n8njtcuglAHTuDPj5AadPu7smRETex2PCHmC/PRFRa2HYExEpAMOeiEgBmgz7tLS0CcHBwXl6vT4/\nMTHxxfrrN23aNN1gMBwOCwvLGT16dEZOTk6Yo2XrY9gTEbUSIUSji9Vq7RgYGFhQWFioq6ysVBsM\nhuzc3NyQuttkZmaOKisr6ymEwPbt2ydERUXtc7Tsf0cCCZsjR4QICRFERNSE/2an3Qyvu9ht2ZtM\npsigoKACnU5XpFarq+Lj47ckJydPrrvNqFGjvunZs+clAIiKitpfXFysdbRsfYMGAYWFQHV1i85f\nRERUj90R7RaLRRMQEGC2vdZqtcX79++Pamz7NWvWzJ44cWKqs2UXLFjw09/dukWjuDgaAwc6/iWI\niLyd0WiE0Whsdnm7Ya9SqRy+4yo9Pf3OtWvXPpqRkTHa2bJ1w95olP32DHsiolrR0dGIjo7+6fXC\nhQudKm+3G0ej0VjMZnOA7bXZbA7QarXF9bfLyckJmzNnzuqUlJS43r17X3SmbH28SEtE5Hp2wz4i\nIuJgfn6+vqioSFdZWdkpKSlpWlxcXErdbU6fPn3L1KlTP924ceOMoKCgAmfKNoRhT0Tkena7cXx8\nfKwrVqyYGxsbu6O6urrj7Nmz14SEhBxftWpVAgAkJCSsWrRo0SsXL17s/fjjj38AAGq1uspkMkU2\nVrapCun1wP79rvlyREQkecxEaDaHDwMzZgBHjrixUkREHs7ZidA8LuyvXAH69ZP/7OBR9/cSEXmO\ndjvrpU23bkDPnsCZM+6uCRGR9/C4sAd4kZaIyNUY9kRECuCRYc9HFBIRuZZHhj1b9kREruWxYZ+f\n7+5aEBF5D48begkAly4BAwbI4ZcqhwcWEREpR7sfegnIoZdduwJnz7q7JkRE3sEjwx5gvz0RkSt5\nbNgPGQJkZLi7FkRE3sEj++wB4LvvgDFjgPR0YPhwN1SMiMiDeUWfPSBb9m+8AcTHA1evurs2RETt\nm8e27AFACGD6dKBHD2DlyjauGBGRB/Oalj0gh12uXAns2gX861/urg0RUfvl0S17G5MJmDQJOHAA\nuOWWNqoYEZEH86qWvU1kJPDss8BDDwFWq7trQ0TU/rSLsAeA554DunQB/vpXd9eEiKj9aRfdODZn\nzwLh4cCWLcAdd7RyxYiIPJhXduPY9O8PrFsnn1F77Zq7a0NE1H60q7AHgAkTgD59gNxcd9eEiKj9\naHdhD8gbrr77zt21ICJqP9pt2J844e5aEBG1H+027NmyJyJyHMOeiEgB2tXQS5vLl+WTrMrL+SQr\nIlImrx56adOjB9C9O2CxuLsmRETtQ7sMe4BdOUREzmDYExEpAMOeiEgBGPZERArQbsN+8GCGPRGR\no9rl0EtAzmvfrRtw8SJw002tUDEiIg+miKGXAODjA9x6K1BQ4O6aEBF5vnYb9gD77YmIHNVk2Kel\npU0IDg7O0+v1+YmJiS/WX5+Xlxc8atSobzp37nx92bJlz9Zdp9PpisLCwnLCw8OzIiMjTa6sOMCw\nJyJylI+9ldXV1R3nzp27Yvfu3eM1Go1l5MiRB+Li4lJCQkKO27bp27fvhffee+/Jbdu23Ve/vEql\nEkajMbpPnz4/tkblhwwBjMbW2DMRkXex27I3mUyRQUFBBTqdrkitVlfFx8dvSU5Onlx3G19f39KI\niIiDarW6qqF9OHMBwVls2RMROcZu2FssFk1AQIDZ9lqr1RZbLBaNoztXqVRi/PjxuyMiIg6uXr16\nTksq2hBb2Lt5QBERkcez242jUqlaFKMZGRmj/f39S0pLS31jYmJ2BQcH540dO/ar+tstWLDgp7+j\no6MRHR3t0P5vvhno2BE4fx7w82tJTYmIPJvRaISxBf3WdsNeo9FYzGZzgO212WwO0Gq1xY7u3N/f\nvwSQXT1Tpkz5zGQyRTYV9s6yPbWKYU9E3qx+Q3jhwoVOlbfbjRMREXEwPz9fX1RUpKusrOyUlJQ0\nLS4uLqWhbev3zV+9erVLeXl5dwCoqKjounPnzrtDQ0OPOFU7B7DfnoioaXZb9j4+PtYVK1bMjY2N\n3VFdXd1x9uzZa0JCQo6vWrUqAQASEhJWnT17tv/IkSMPXL58uUeHDh1qli9fPi83N3fo+fPn+02d\nOvVTALBarT7Tp0/fdPfdd+909Rdg2BMRNa3dTpdg89lnwLp1QEqDvzeIiLyTYqZLsGHLnoioae2+\nZX/jBtCzp3werVrtwooREXkwxbXsf/UrQKsFTp50d02IiDxXuw97gHPbExE1xSvCnv32RET2MeyJ\niBSAYU9EpAAMeyIiBfCKsPf3B65fl8+jJSKiX/KKsFepOCKHiMgerwh7gF05RET2MOyJiBSAYU9E\npABeFfYnTri7FkREnqndT4RmU1EB+PrKCdE6dnRBxYiIPJjiJkKz6dpVPpP29Gl314SIyPN4TdgD\n7LcnImoMw56ISAEY9kRECuB1YX/8uLtrQUTkebxmNA4gR+JotcD338uLtURE3kqxo3EAoHt3YOJE\nYOtWd9eEiMizeFXYA8D06cCmTe6uBRGRZ/GqbhwAqKoCBgwA9u8HBg1y2W6JiDyKortxAECtBqZN\nA/75T3fXhIjIc3hd2APAjBnAxo2Am3+0EBF5DK8M+6gowGoFDh1yd02IiDyDV4a9SlXbuiciIi+8\nQGuTnw+MHQsUFwM+Pi7fPRGRWyn+Aq2NXg/odMDu3e6uCRGR+3lt2APsyiEisvHabhwAKC2VLfzi\nYqBbt1b5CCIit2A3Th2+vsCYMcC2be6uCRGRe3l12AOcPoGICPDybhwAuHoV0GiAvDzAz6/VPoaI\nqE2xG6eeLl2AuDhgyxZ314SIyH2aDPu0tLQJwcHBeXq9Pj8xMfHF+uvz8vKCR40a9U3nzp2vL1u2\n7FlnyrYVjsohIqWz241TXV3dcciQId/t3r17vEajsYwcOfLA5s2bHwwJCfnpeVClpaW+p06dGrht\n27b7evfuffHZZ59d5mhZoPW7cWRdgIAAYO9eIDi4VT+KiKhNuLQbx2QyRQYFBRXodLoitVpdFR8f\nvyU5OXly3W18fX1LIyIiDqrV6ipny7aVjh2BP/4R+O1vgcmTgcWL5c1Wly65ozZERG3PbthbLBZN\nQECA2fZaq9UWWywWjSM7bknZ1rBwIXDggOzS+fFHYNEieeE2JASYPRsoK3NXzYiIWp/dWWNUKlWz\n+1ecKbtgwYKf/o6OjkZ0dHRzP9ZOfYCBA+XywAPyPasVOHYMePllYOVKYP58l38sEZFLGI1GGI3G\nZpe3G/YajcZiNpsDbK/NZnOAVqstdmTHzpStG/ZtyccHMBiABQuAKVOA557jpGlE5JnqN4QXLlzo\nVHm73TgREREH8/Pz9UVFRbrKyspOSUlJ0+Li4lIa2rb+hQJnyrrbbbfJC7gpHlk7IiIXEELYXVJT\nU+8ZPHjwd4GBgQWLFy9+SQiBlStXJqxcuTJBCIGSkpL+Wq3W3KNHj0u9evW6GBAQcLq8vLxbY2Xr\nL7IK7rd5sxDR0a7Z12OPCfHNN67ZFxFRQ/6bnU1muG3x+jtoHVVVJadE3r4dCAtr/n6uXAH69gUe\nfBD46CNX1Y6I6Od4B20zqdVAQgKwYkXL9pOeLsfyJyfLqRqIiDwBw76OhARg61Y5NLO5UlOBhx+W\nz8FNTnZd3YiIWoJhX4efH3DvvcDatc0rL4TsBpo4EZg5E9iwwbX1IyJqLvbZ12MyAdOmAQUF8s5b\nZ+TmyqAvLKydbfO77zjbJhG5HvvsWygyEujXD/jiC+fL2lr1KhXQtaucbTMpyfV1JCJyFsO+AU8+\nCbz3nvPlUlOBe+6pfT1jBrtyiMgzsBunATduyGGYe/YAQ4c6Vqa8XHbblJTIVj0gp2MICKgdoUNE\n5CrsxnGBX/0K+NOfnBuGuWcPcPvttUEPyKkXHnqIc+kTkfsx7BuRkCCfbuXoNMipqbK/vr4ZM+Qz\ncGtqXFs/IiJnMOwbMWAAEBsLrFvX9La2IZd1++ttRoyQrf3MTNfXkYjIUQx7O556Cnj3XdmHb8/R\no0CnTsDgwb9cp1LxQi0RuR/D3o5Ro4DQUOCtt+xvZ2vVqxq5VPLQQ8AnnzR90iAiai0M+ya8/Tbw\n5puA2dz4No3119vccoucXK05Y/eJiFyBYd+EQYOAuXPlg00acukS8O23QFMP15oxg6NyiMh9GPYO\nePFFYP9+YO/eX67bvRsYPRro0sX+Pu6/Xw7PbMkka0REzcWwd0CXLrI756mn5Lz3dTU2Cqe+Xr3k\n6J6tW1unjkRE9jDsHXTfffIO2fffr32v7iyXjuCoHCJyF4a9g1QqOQzztdeAs2flezk5stUfFOTY\nPiZMACwW2fVDRNSWGPZOGDIEePRRYP58+bqpUTj1deoELF8OPPGEc8MwT5wASkudqysRUV0Meye9\n/LJsmWdmOt5fX1dcnLz5qqmx+zZnzgB33CHH+7O/n4iai7NeNsPmzbI75/Rp4Nw54KabnCt/8iQw\ncqQcsjlwYOPbVVcDMTEy7GNjgVmzgOHD5XWDfv3sf8b168CpU8Ctt8pfFETkXTjrZRuIjwf69gXG\njnU+6AE5dn/ePOCZZ+xv99pr8iLwyy/LGTWzsmTZsDDg449/ub3VCuzYAfzhD4C/vzxB9OwJDBsG\nPPAAsGCBLHfsGCdmI1Iatuyb6cwZ4MqVhufDccT167KV/t57DXcF/ec/8qRy6JCclK2ufftqW/kr\nVshHKP7zn3JKhoED5fQMv/+9LHf9unw04rFj8rGJx44BBw/KkUFLljSv7kTkfs627Bn2brR9u3wq\n1tGjQOfOte+XlgK//jWwerUcwdOQ69eBV1+V4//1euDBB+XJwZGRQadPy9k4T5wAbr7ZNd+FiNoW\nw76dmTpVBu8rr8jXNTXApEmy6+X115suf/ky0L1745OwNWbOHKB/f+Cvf3W+zkTkfgz7dub0aSA8\nHDhwQPbHL1smu2O+/BJQq1vvc0+elA9XLyiQd/cSUfvCsG+HliwBMjJk6/7eewGTST4Dt7U9/LC8\nd+Avf2n9zyIi12LYt0OVlXKEzfnzwN//Lrt22kJeHvDb38pWfrdubfOZROQaDPt2KjMT+OorOcNm\nW5o2TY75b2wKZyLyTAx7ckpOjhyPf/Jk8+4ZICL34E1V5JSwMHmhds0a1+7XagW+/hpYtKh24jgi\nch+27AkHDsiHqxQUtGxqhR9/BNLS5OMXd+wAtFqgTx95sXntWpdV1+2KioB164CFC91dE1IytuzJ\naSNHAkOHAuvXN6/89u3AmDEy1LdskXP5ZGfL5bPP5Oyg2dkurbJbLV0qf7GYTO6uCZHj2LInALLL\n5eGH5V21Pj6Ol/vsM+Cxx4BVq+TdvnXvBLZ5/3253a5dzt/85WnOnQOCg+VTy7KygJQUd9eIlIot\ne2qWMWOAW26RM3o66tNPZdBv3y6f5NVQ0APAn/4EFBfLFn57t2KFnJZi/nw5x5A3/WIh78aWPf1k\nzx75YJXs7MaD2+aTT4C5c2XQh4c3ve9//1sOK83Jce6Xgye5ckVOGZ2ZKecjeust+fcnn7i7ZqRE\nLm/Zp6WlTQgODs7T6/X5iYmJDY4Cf+qpp97V6/X5BoPhcFZW1k//6+t0uqKwsLCc8PDwrMjISPZw\nerhx42T/fVCQnI3z+vWGt9u6VQb9jh2OBT0g7wzu31/eNNZerV0rr0fo9fJ1QoK8N+LYMffWi8gh\nQohGF6vV2jEwMLCgsLBQV1lZqTYYDNm5ubkhdbf54osvJt5zzz2pQgjs27cvKioqap9tnU6nK7xw\n4UIfe58hq0Ce5OBBIeLihBgwQIh33hHi6tXadVu2CNG/vxCHDzu/32+/lWUvXXJdXdtKVZUQAwcK\nsW/fz99fskSIhx5yS5VI4f6bnXYzvO5it2VvMpkig4KCCnQ6XZFara6Kj4/fkpycPLnuNikpKXGP\nPPLIegCIioraX1ZW1uvcuXN+dU4m7fySnPLcdhuQnAx8/jlgNAKBgXIq5Y8+Ap5+Gti5U47Pd1Z4\nuLyBa+lSV9fYPqtVdsG0xNat8ppGVNTP33/iCXk88vNbtn+i1ma399RisWgCAgLMttdarbZ4//79\nUU1tY7FYNH5+fudUKpUYP3787o4dO1YnJCSsmjNnzuqGPmfBggU//R0dHY3o6Ohmfh1ypfBwOYom\nO1sONTx0SI6oGT68+ft87TV5onjsMRmerc1olCNnTp+WfeyzZjk/IkgI4I035DGor3t3+UyCxYvl\n2Hui1mI0GmE0Gptd3m7Yq1Qqh66cNtZ6//rrr8cMGDDgTGlpqW9MTMyu4ODgvLFjx35Vf7u6YU+e\nZ8QIOfLGFTQa2Rr+85+BjRtds8+GmM1yvp99+4A335TDJWfOBLZtkw+F8fNreh82e/bI6xcTJza8\n/skn5XWOwkJ5AddZJpMcsWQ0crppalz9hvBCJ+/qs9uNo9FoLGazOcD22mw2B2i12mJ72xQXF2s1\nGo0FAAYMGHAGAHx9fUunTJnymclkinSqduSVXngB2LtXDl10tevXgb/9TZ6ggoOB48fl83dDQ2Wo\nhoYCBoNzJ6833gCefx7o0Mj/Lb17y18qzemeqqwEHn1U/v3SS86XJ3KYvQ79qqoqn0GDBn1fWFio\nu3HjRqemLtB+8803t9su0FZUVHS5fPlydyEErly50vU3v/lNxo4dO+6u/xngBVpF+vBDIYYPF+KF\nF4R49VUhli4VYvly+f6GDUKcOOH8PlNShLj1ViGmThXi5MnGt8vMFCIoSIiZM4W4eNH+PrOz5YXq\n69ftb1daKkTv3kKYzc7V+dVXhbj3XiF+/FEIf38hMjKcK98Wqqvlv5+6F+rJ/eDkBdomN0hNTb1n\n8ODB3wUGBhYsXrz4JSEEVq5cmbBy5coE2zZPPPHEisDAwIKwsLDDhw4d+rUQAt9///0gg8GQbTAY\nsocNG3bUVvYXFWDYK5LVKsTf/y5Hs7zyihDPPy/E3LlCzJ4txLRpctTOqVOO7+/jj4XQaoXYudOx\n7a9cEeJ//1eIgAAh1q5tfITQ9OlCJCY6ts/nnhPiyScd21YIIXJyhLj55toTRFKSEMOGCXHjhuP7\naMzRoy3fh83q1TIp3n3XdfuklnN52Lf2wrCnhrz5phC//rUQFRVNb5uVJUPz22+d/5w9e4SYPFmI\nnj2FiI8X4vPPhaislOuKioTo00eIsjLH9lVSIlv3JSVNb1tVJcTIkUKsWlX7Xk2NEBMnCvHaa85/\nj7o2bJD/Z7/+esv2I4QQFy4I0a+fPDEPGCDEtWst3ye5hrNhzztoySMJAcyYIUfObNjQ+Aia8+fl\nFM2JifJBLM114QLw8cfys77/Xk6JcP68nLnzjTcc38+zz8rRSx9/DPTt2/h2b74pp4/Ys+fn362o\nCIiIAL75pvbmLWecPCmHh65bB8ybJ0cizZvn/H5s5s4FqquBDz4AJk8Gxo+XF6TJ/Zy9g5Yte/JY\nFRVChIcLsWxZw+tv3BBi7Fgh/vxn135uQYEQCxYIcfvtzvfBV1UJ8eyz8tpBVlbD2+TnC9G3r/yc\nhrz5phDjxsmWvjMqK4WIihLi7bfl66IieSPYBx84tx+b7GzZqv/hB/n60CHnWvcbNwqxfn3zPpua\nBnbjkDc5dUr23zfUF5+QIMSkSfICoqfZvFl2LW3c+PP3q6uFuOMOId56q/GyVVXyJOdsUP7lL0JM\nmPDz41FQIK9lrFnj3L5qaoQYM0aIlSt//n5cnGN994cPy+8/eLC8FmPrGiPXYdiT1zEaZQvz++9r\n3/u//xNi6FDPnnrh8GEhAgOFmDevNuxWrpStb6vVftmDB4Xw85OjfBxhNMrRPGfP/nJdXp5skW/Y\n4HjdN26U10zq19OR1v21a3Kk1UcfydFOEycKER0txPnzjn8+NY1hT15pxQohQkOFKC+vDf/8fHfX\nqmk//ijEPffI1vzBg7K16+hImaefFuLhh5ve7sIFOapo+/bGtzl2TP5CSkpqen+XL8tAb2wYaFOt\n+3nzhHjggdpuKKtViJdekl1KjXVt2Vy7Zn/YrLOuXZMX4V3166+yUoitW10zYqqlGPbklWpqhHj0\nURmcjXXreCqrVYiXXxaiQwchFi1yvFx5uRC33CKHlTbWf19TI+8rePrppvd3+LD8tfDRR/bD7/nn\nhXjkkcbXHzzYeOs+LU2eeC5c+OW6pCR5stuy5efvnzolf/HExQnRo4e8nnHbbUK8/37T90HY8/XX\nQgwZIoRGI8T48UJYLM3flxDyl2VkpPx+o0a1fH8txbAnr3X9urxw2V7He2dny/54Z6SnCxESIoRO\nJ1vHR478fP2HHwphMDR905fNt9/KIB0yRA6nrF/u+HEZyE0NH5006Zf/Hs6flyeBPXsaL5eVJb/L\n00/LG+qGD5fhPn26EJs2yYvBVqs8aTzwgBwSO326EHv3Ot46Ly+X9zv4+wvxySfymC9YIE90n37q\n2D7qs12DeecdWb+//U1+1//8p3n7cwWGPZGXqamRJ4oXXpCtytBQIRYvlr9ubr5ZiNxc5/e3d68Q\nsbEysBIT5b0ENTVCxMTYv3hsc/CgbDHbWvc1NfJ+heeea7rs+fNCzJolf+1kZtq/flFaKgM2NFSI\nQYOEePFFeS9EYy3+HTtkd9Ef/vDLXxeZmXIff/yjPCE44soV+YtSr5fXK+pKS5Pdie+84/zIKVdg\n2BN5sepqIb78UojHHpNh6+wom/qys2XLuU8fIX7/e3nR29GRM5MmCfHee/LvDz8UYsQIx39hOKum\nRgiTSU4vMW6cEN26CREWJsQTT8huoePHZcDrdDLwG3P5stxOr5f7s+fwYSGCg+W0GpcvN7zNyZPy\nez/0kDwxtCVnw543VRERTp2Sz9e9/37g9tsdK3PokLzRKjVVPuXsyy+BoUNbt542VVXAt9/KJ4V9\n9RVw4ADwP/8jp5ru1q3p8ranrY0bB3TtKh+V2bGj/KePD3DjBpCUBCxbBjz8sP19Xb0KPP64fAD9\n2rXAwIFy9lK12jXftTHO3lTFsCeiZouLk1MzL14sw7M9sVjk8xmsVrlUV9f+XVMD3Hef43cxCyHv\nMn77beDiRaCsTD7HuVev2qVfPznFt0Yj78y2/a3RyBOOsxj2RNRmjhyRzwdYvtz5h8J4MyGAigoZ\n+mVl8gRw7pw8wVgsQHFx7d8WC1BQIEPfGQx7IqJ2xBZ/zp4snQ17u0+qIiKi1tVWv4jsPqmKiIi8\nA8OeiEgBGPZERArAsCciUgCGPRGRAjDsiYgUgGFPRKQADHsiIgVg2BMRKQDDnohIARj2REQKwLAn\nIlIAhj0RkQIw7ImIFIBhT0SkAAx7IiIFYNgTESkAw56ISAEY9kRECsCwJyJSAIY9EZECNBn2aWlp\nE4KDg/P0en1+YmLiiw1t89RTT72r1+vzDQbD4aysrHBnylIto9Ho7ip4DB6LWjwWtXgsms9u2FdX\nV3ecO3fuirS0tAm5ublDN2/e/ODx48dD6m6Tmpo6saCgICg/P1//4Ycf/unxxx//wNGy9HP8D7kW\nj0UtHovuxTRYAAAFc0lEQVRaPBbNZzfsTSZTZFBQUIFOpytSq9VV8fHxW5KTkyfX3SYlJSXukUce\nWQ8AUVFR+8vKynqdPXu2vyNliYiobdgNe4vFogkICDDbXmu12mKLxaJxZJszZ84MaKosERG1DR97\nK1UqlXBkJ0IIVUsqoVK1qLhXWbhwobur4DF4LGrxWNTisWgeu2Gv0WgsZrM5wPbabDYHaLXaYnvb\nFBcXa7VabXFVVZW6qbJAy08URETUNLvdOBEREQfz8/P1RUVFusrKyk5JSUnT4uLiUupuExcXl/KP\nf/zjYQDYt2/f7b169Srz8/M750hZIiJqG3Zb9j4+PtYVK1bMjY2N3VFdXd1x9uzZa0JCQo6vWrUq\nAQASEhJWTZw4MTU1NXViUFBQQdeuXSvWrVs3y17ZtvhSRERUjxDCbcv27dsnDBkyJC8oKCh/6dKl\nL7qzLm29zJo1a22/fv3ODR8+/IjtvQsXLvQZP378Lr1efyImJmbnxYsXe7m7nq29nD59OiA6Ojp9\n6NChx4YNG3Z0+fLlTyn1WFy7dq1zZGTkfoPBkB0SEpI7f/78JUo9FrbFarV2HDFiRNa99977byUf\ni4EDBxaFhobmjBgxImvkyJGm5hwLt91Bq/Rx+LNmzVqXlpY2oe57S5cunR8TE7PrxIkTg++66649\nS5cune+u+rUVtVpd9fbbbz9z7NixYfv27bv9/ffff+L48eMhSjwWnTt3vp6enn5ndnb2iJycnLD0\n9PQ7v/766zFKPBY2y5cvnzd06NBc22ARpR4LlUoljEZjdFZWVrjJZIoEmnEs3HWmyszMHBUbG5tm\ne71kyZL5S5Ysme/uM2hbLoWFhbq6LfshQ4bknT171k8IgZKSkv5DhgzJc3cd23qZPHnytl27do1X\n+rGoqKjoEhERceDo0aPDlHoszGaz9q677tq9d+/eO20te6UeC51OV/jDDz/0rfues8fCbS17R8bw\nK825c+f8/Pz8zgGAn5/fuXPnzvm5u05tqaioSJeVlRUeFRW1X6nHoqampsOIESOy/fz8zt15553p\nw4YNO6bUY/HMM8+8/cYbbzzfoUOHGtt7Sj0WKpVKjB8/fndERMTB1atXzwGcPxZ2L9C2JkfH8CuV\nSqUSSjpGV65c6Xb//ff/a/ny5fO6d+9eXnedko5Fhw4darKzs0dcunSpZ2xs7I709PQ7665XyrH4\n/PPP7+3Xr9/58PDwLKPRGN3QNko5FgCQkZEx2t/fv6S0tNQ3JiZmV3BwcF7d9Y4cC7e17B0Zw680\nfn5+586ePdsfAEpKSvz79et33t11agtVVVXq+++//18zZ87ccN99920DlHssbHr27Hnpd7/73ReH\nDh26TYnHIjMz8zcpKSlxt956a+GDDz64ee/eveNmzpy5QYnHAgD8/f1LAMDX17d0ypQpn5lMpkhn\nj4Xbwp7j8H8pLi4uZf369Y8AwPr16x+xBZ83E0KoZs+evWbo0KG5Tz/99Du295V4LH744Yeby8rK\negHAtWvXbtq1a1dMeHh4lhKPxeLFi/9sNpsDCgsLb92yZUv8uHHj9m7YsGGmEo/F1atXu5SXl3cH\ngIqKiq47d+68OzQ09IjTx8KdFx1SU1PvGTx48HeBgYEFixcvfsndF0HacomPj9/s7+9/Rq1WV2q1\nWvPatWtnXbhwoc9dd921W0nDyr766qsxKpWqxmAwZI8YMSJrxIgRWdu3b5+gxGORk5MTGh4e/q3B\nYMgODQ3Nef31158XQg6xU9qxqLsYjcY7Jk2alKLUY3Hy5MlbDQZDtsFgyB42bNhRW1Y6eyxUQiii\ny4uISNH4pCoiIgVg2BMRKQDDnohIARj2REQKwLAnIlIAhj0RkQL8P3pAaPj7U01rAAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x34e5390>"
       ]
      }
     ],
     "prompt_number": 4,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4,
     "trusted": true
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "An (explicit) graphical model example"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# let's make a small number of nodes gaussians, depending on each others\n",
      "N = 7\n",
      "A = np.zeros((N,N))\n",
      "A[0,:] = [1., .8, .3, 0, 0 ,0, 0] \n",
      "A[1,:] = [0, 1., .75, -.3, 0, 0 ,0] \n",
      "A[2,:] = [0, 0., 1., -.3, .4, 0 ,0] \n",
      "A[3,:] = [0, 0, 0, 2, .6, .4, 0] \n",
      "A[4,:] = [0, 0, 0, 0, 1.5, -.2, 0] \n",
      "A[5,:] = [0, 0, 0, 0, .2, 1.2, -.2] \n",
      "A[6,:] = [0, 0, 0, 0, 0, .2, .5] "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T = 10\n",
      "N = A.shape[0]\n",
      "x = np.random.normal(0,1,size=(N,T))\n",
      "y = A.dot(x)\n",
      "print(y.shape)\n",
      "np.set_printoptions(precision=2)\n",
      "\n",
      "print(\"theoretical corr:\\n\",cov2corr((A.dot(A.T))))\n",
      "print(\"experimental corr:\\n\",np.corrcoef(y),\"\\n\")\n",
      "print(\"theoretical cov:\\n\",((A.dot(A.T))))\n",
      "print(\"experimental cov:\\n\",np.cov(y))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(7, 10)\n",
        "theoretical corr:\n",
        " [[ 1.    0.61  0.2   0.    0.    0.    0.  ]\n",
        " [ 0.61  1.    0.58 -0.22  0.    0.    0.  ]\n",
        " [ 0.2   0.58  1.   -0.15  0.35  0.06  0.  ]\n",
        " [ 0.   -0.22 -0.15  1.    0.25  0.23  0.07]\n",
        " [ 0.    0.    0.35  0.25  1.    0.03 -0.05]\n",
        " [ 0.    0.    0.06  0.23  0.03  1.    0.21]\n",
        " [ 0.    0.    0.    0.07 -0.05  0.21  1.  ]]\n",
        "experimental corr:\n",
        " [[ 1.    0.82  0.62 -0.17  0.15 -0.56  0.07]\n",
        " [ 0.82  1.    0.8  -0.47 -0.03 -0.58  0.13]\n",
        " [ 0.62  0.8   1.   -0.59  0.39 -0.14  0.5 ]\n",
        " [-0.17 -0.47 -0.59  1.   -0.09  0.28 -0.25]\n",
        " [ 0.15 -0.03  0.39 -0.09  1.    0.24  0.08]\n",
        " [-0.56 -0.58 -0.14  0.28  0.24  1.    0.33]\n",
        " [ 0.07  0.13  0.5  -0.25  0.08  0.33  1.  ]] \n",
        "\n",
        "theoretical cov:\n",
        " [[ 1.73  1.02  0.3   0.    0.    0.    0.  ]\n",
        " [ 1.02  1.65  0.84 -0.6   0.    0.    0.  ]\n",
        " [ 0.3   0.84  1.25 -0.36  0.6   0.08  0.  ]\n",
        " [ 0.   -0.6  -0.36  4.52  0.82  0.6   0.08]\n",
        " [ 0.    0.    0.6   0.82  2.29  0.06 -0.04]\n",
        " [ 0.    0.    0.08  0.6   0.06  1.52  0.14]\n",
        " [ 0.    0.    0.    0.08 -0.04  0.14  0.29]]\n",
        "experimental cov:\n",
        " [[ 3.3   3.1   1.47 -0.73  0.35 -1.49  0.05]\n",
        " [ 3.1   4.29  2.16 -2.35 -0.09 -1.75  0.11]\n",
        " [ 1.47  2.16  1.72 -1.85  0.66 -0.27  0.29]\n",
        " [-0.73 -2.35 -1.85  5.71 -0.27  0.96 -0.26]\n",
        " [ 0.35 -0.09  0.66 -0.27  1.64  0.45  0.04]\n",
        " [-1.49 -1.75 -0.27  0.96  0.45  2.11  0.21]\n",
        " [ 0.05  0.11  0.29 -0.26  0.04  0.21  0.19]]\n"
       ]
      }
     ],
     "prompt_number": 6,
     "trusted": true
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Conclusion: estimation can be real bad with a small number of points.\n",
      "See how bad it gets if you have less data than regions"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Graph Lasso Estimation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "skl.__file__"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "'/home/jb/.local/lib/python2.7/site-packages/sklearn/__init__.pyc'"
       ]
      }
     ],
     "prompt_number": 7,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.covariance import GraphLassoCV, GraphLasso, LedoitWolf"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8,
     "trusted": true
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####Scikit-learn : cite:\n",
      "\n",
      "Pedregosa, Fabian, Ga\u00ebl Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, and Vincent Dubourg. \u201cScikit-Learn: Machine Learning in Python.\u201d The Journal of Machine Learning Research 12 (2011): 2825\u201330."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "GL = GraphLassoCV(alphas=10, mode='cd')\n",
      "GL.fit(y.T)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "GraphLassoCV(alphas=10, assume_centered=False, cv=None, max_iter=100,\n",
        "       mode='cd', n_jobs=1, n_refinements=4, tol=0.0001, verbose=False)"
       ]
      }
     ],
     "prompt_number": 9,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C_GL = GL.covariance_\n",
      "Prec_GL = GL.precision_"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10,
     "trusted": true
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Compare correlations"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.set_printoptions(precision=3)\n",
      "Ctheo = A.dot(A.T)\n",
      "print(\"theoretical corr:\\n\",cov2corr(Ctheo))\n",
      "print(\"experimental corr:\\n\",np.corrcoef(y))\n",
      "print(\"GL corr:\\n\",cov2corr(C_GL))\n",
      "print(\" corr_theo - corr_exp : {:2.3f}\".format(frob_norm(cov2corr(Ctheo)-np.corrcoef(y))))\n",
      "print(\" corr_theo - corr_GL : {:2.3f}\".format(frob_norm(cov2corr(Ctheo)-cov2corr(C_GL))))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "theoretical corr:\n",
        " [[ 1.     0.606  0.204  0.     0.     0.     0.   ]\n",
        " [ 0.606  1.     0.584 -0.22   0.     0.     0.   ]\n",
        " [ 0.204  0.584  1.    -0.151  0.355  0.058  0.   ]\n",
        " [ 0.    -0.22  -0.151  1.     0.255  0.229  0.07 ]\n",
        " [ 0.     0.     0.355  0.255  1.     0.032 -0.049]\n",
        " [ 0.     0.     0.058  0.229  0.032  1.     0.211]\n",
        " [ 0.     0.     0.     0.07  -0.049  0.211  1.   ]]\n",
        "experimental corr:\n",
        " [[ 1.     0.825  0.617 -0.168  0.152 -0.565  0.069]\n",
        " [ 0.825  1.     0.796 -0.474 -0.034 -0.582  0.127]\n",
        " [ 0.617  0.796  1.    -0.591  0.392 -0.141  0.502]\n",
        " [-0.168 -0.474 -0.591  1.    -0.088  0.275 -0.25 ]\n",
        " [ 0.152 -0.034  0.392 -0.088  1.     0.245  0.075]\n",
        " [-0.565 -0.582 -0.141  0.275  0.245  1.     0.325]\n",
        " [ 0.069  0.127  0.502 -0.25   0.075  0.325  1.   ]]\n",
        "GL corr:\n",
        " [[ 1.     0.727  0.469 -0.253  0.051 -0.425  0.   ]\n",
        " [ 0.727  1.     0.66  -0.4    0.086 -0.459  0.   ]\n",
        " [ 0.469  0.66   1.    -0.473  0.172 -0.291  0.   ]\n",
        " [-0.253 -0.4   -0.473  1.    -0.078  0.173  0.   ]\n",
        " [ 0.051  0.086  0.172 -0.078  1.     0.046  0.   ]\n",
        " [-0.425 -0.459 -0.291  0.173  0.046  1.     0.   ]\n",
        " [ 0.     0.     0.     0.     0.     0.     1.   ]]\n",
        " corr_theo - corr_exp : 0.526\n",
        " corr_theo - corr_GL : 0.289\n"
       ]
      }
     ],
     "prompt_number": 11,
     "trusted": true
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Compare Covariances"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"theoretical cov:\\n\",Ctheo)\n",
      "print(\"experimental cov:\\n\",np.cov(y))\n",
      "print(\"GL cov:\\n\",C_GL)\n",
      "print(\" Ctheo - Cexp : {:2.3f}\".format(frob_norm(Ctheo-np.cov(y))))\n",
      "print(\" Ctheo - C_GL : {:2.3f}\".format(frob_norm(Ctheo-C_GL)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "theoretical cov:\n",
        " [[ 1.73   1.025  0.3    0.     0.     0.     0.   ]\n",
        " [ 1.025  1.653  0.84  -0.6    0.     0.     0.   ]\n",
        " [ 0.3    0.84   1.25  -0.36   0.6    0.08   0.   ]\n",
        " [ 0.    -0.6   -0.36   4.52   0.82   0.6    0.08 ]\n",
        " [ 0.     0.     0.6    0.82   2.29   0.06  -0.04 ]\n",
        " [ 0.     0.     0.08   0.6    0.06   1.52   0.14 ]\n",
        " [ 0.     0.     0.     0.08  -0.04   0.14   0.29 ]]\n",
        "experimental cov:\n",
        " [[ 3.299  3.104  1.47  -0.729  0.352 -1.49   0.054]\n",
        " [ 3.104  4.294  2.163 -2.35  -0.089 -1.751  0.114]\n",
        " [ 1.47   2.163  1.72  -1.851  0.658 -0.268  0.287]\n",
        " [-0.729 -2.35  -1.851  5.711 -0.269  0.955 -0.26 ]\n",
        " [ 0.352 -0.089  0.658 -0.269  1.635  0.454  0.042]\n",
        " [-1.49  -1.751 -0.268  0.955  0.454  2.109  0.205]\n",
        " [ 0.054  0.114  0.287 -0.26   0.042  0.205  0.189]]\n",
        "GL cov:\n",
        " [[ 2.969  2.461  1.006 -0.988  0.106 -1.008  0.   ]\n",
        " [ 2.461  3.864  1.614 -1.782  0.204 -1.243  0.   ]\n",
        " [ 1.006  1.614  1.548 -1.334  0.26  -0.499  0.   ]\n",
        " [-0.988 -1.782 -1.334  5.14  -0.214  0.54   0.   ]\n",
        " [ 0.106  0.204  0.26  -0.214  1.472  0.077  0.   ]\n",
        " [-1.008 -1.243 -0.499  0.54   0.077  1.898  0.   ]\n",
        " [ 0.     0.     0.     0.     0.     0.     0.17 ]]\n",
        " Ctheo - Cexp : 7.549\n",
        " Ctheo - C_GL : 4.147\n"
       ]
      }
     ],
     "prompt_number": 12,
     "trusted": true
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Ledoit-Wolf estimator"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(1 - shrinkage)*cov\n",
      "        + shrinkage*mu*np.identity(n_features)\n",
      "\n",
      "with \n",
      "\n",
      "* mu from \n",
      "     mu = trace(cov) / n_features\n",
      "* shrinkage from\n",
      "    \u201cA Well-Conditioned Estimator for Large-Dimensional Covariance Matrices\u201d, Ledoit and Wolf, Journal of Multivariate Analysis, Volume 88, Issue 2, February 2004, pages 365-411."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "LW = LedoitWolf(store_precision=True, assume_centered=False, block_size=1000)\n",
      "LW.fit(y.T)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "LedoitWolf(assume_centered=False, block_size=1000, store_precision=True)"
       ]
      }
     ],
     "prompt_number": 13,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C_LW = LW.covariance_"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"theoretical cov:\\n\",Ctheo)\n",
      "print(\"experimental cov:\\n\",np.cov(y))\n",
      "print(\"LW cov:\\n\",C_LW)\n",
      "print(\" Ctheo - Cexp : {:2.3f}\".format(frob_norm(Ctheo-np.cov(y))))\n",
      "print(\" Ctheo - C_GL : {:2.3f}\".format(frob_norm(Ctheo-C_GL)))\n",
      "print(\" Ctheo - C_LW : {:2.3f}\".format(frob_norm(Ctheo-C_LW)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "theoretical cov:\n",
        " [[ 1.73   1.025  0.3    0.     0.     0.     0.   ]\n",
        " [ 1.025  1.653  0.84  -0.6    0.     0.     0.   ]\n",
        " [ 0.3    0.84   1.25  -0.36   0.6    0.08   0.   ]\n",
        " [ 0.    -0.6   -0.36   4.52   0.82   0.6    0.08 ]\n",
        " [ 0.     0.     0.6    0.82   2.29   0.06  -0.04 ]\n",
        " [ 0.     0.     0.08   0.6    0.06   1.52   0.14 ]\n",
        " [ 0.     0.     0.     0.08  -0.04   0.14   0.29 ]]\n",
        "experimental cov:\n",
        " [[ 3.299  3.104  1.47  -0.729  0.352 -1.49   0.054]\n",
        " [ 3.104  4.294  2.163 -2.35  -0.089 -1.751  0.114]\n",
        " [ 1.47   2.163  1.72  -1.851  0.658 -0.268  0.287]\n",
        " [-0.729 -2.35  -1.851  5.711 -0.269  0.955 -0.26 ]\n",
        " [ 0.352 -0.089  0.658 -0.269  1.635  0.454  0.042]\n",
        " [-1.49  -1.751 -0.268  0.955  0.454  2.109  0.205]\n",
        " [ 0.054  0.114  0.287 -0.26   0.042  0.205  0.189]]\n",
        "LW cov:\n",
        " [[ 2.73   1.54   0.729 -0.361  0.175 -0.739  0.027]\n",
        " [ 1.54   3.224  1.073 -1.165 -0.044 -0.868  0.057]\n",
        " [ 0.729  1.073  1.947 -0.918  0.326 -0.133  0.142]\n",
        " [-0.361 -1.165 -0.918  3.927 -0.133  0.474 -0.129]\n",
        " [ 0.175 -0.044  0.326 -0.133  1.905  0.225  0.021]\n",
        " [-0.739 -0.868 -0.133  0.474  0.225  2.14   0.102]\n",
        " [ 0.027  0.057  0.142 -0.129  0.021  0.102  1.188]]\n",
        " Ctheo - Cexp : 7.549\n",
        " Ctheo - C_GL : 4.147\n",
        " Ctheo - C_LW : 1.876\n"
       ]
      }
     ],
     "prompt_number": 15,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15,
     "trusted": true
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Trying out other things : not relevant today"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Simulate a sparse precision matrix.  Conclusion: if N>=50, it becomes hard to fin tril matrices that are not rank deficient"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def iround(x):\n",
      "    \"\"\"iround(number) -> integer\n",
      "    Round a number to the nearest integer.\"\"\"\n",
      "    y = round(x) - .5\n",
      "    return int(y) + (y > 0)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# create a covariance matrix of N nodes with T time values (repetitions) and sparsity of \n",
      "# precision matrix is L% \n",
      "\n",
      "\n",
      "def sim_sqrtP(N, L, scalef=None):\n",
      "    \"\"\"\n",
      "    N : create a NxN matrix\n",
      "    L : put L% zeros in the sqrtP matrix\n",
      "    scalef : scaling factor of the entries of P\n",
      "    \"\"\"\n",
      "    \n",
      "    if scalef is None:\n",
      "        scalef = N \n",
      "    \n",
      "    def arr_entries(n, scalef=scalef):\n",
      "        # returns n values drawn from N(0,1) and multiplied by f\n",
      "        return np.abs(np.random.randn(n))*scalef    \n",
      "    \n",
      "    # get the (N,N) strict lower triang matrix indices from a (N-1,N-1) lower tri matrix  \n",
      "    zeros_arr = np.zeros((N-1,N-1))\n",
      "    tri_indices = np.tril_indices_from(zeros_arr)\n",
      "    n_tri = tri_indices[0].shape[0]\n",
      "\n",
      "    # precision matrix\n",
      "    P = np.zeros((N,N))\n",
      "               \n",
      "    n_zeros = iround(n_tri*L)\n",
      "    # generate the strict low tri values \n",
      "    z = arr_entries(n_tri)\n",
      "    z[:n_zeros] = 0\n",
      "    z = np.random.permutation(z)\n",
      "    # fill P : push the indices by one towards the bottom of P\n",
      "    P[(tri_indices[0]+1, tri_indices[1])] = z\n",
      "    # fill the diagonal - sq\n",
      "    P[(np.arange(N),np.arange(N))] = arr_entries(N, scalef=scalef)\n",
      "    \n",
      "    return P\n",
      "\n",
      "#    rkP = np.linalg.matrix_rank(P)\n",
      "#    if rkP < N:\n",
      "#        print(\"rank deficient {}\".format(rkP))\n",
      "#        return None\n",
      "    \n",
      "\n",
      "#    try:\n",
      "#        assert int((N)*(N-1)/2) == n_tri\n",
      "#    except:\n",
      "#        print(\"Aie : {} differ from {}\".format(int((N)*(N-1)/2), n_tri))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 6\n",
      "Z = np.zeros((N,N))\n",
      "sqrtP = sim_sqrtP(7, .5, 1.)\n",
      "np.set_printoptions(precision=1)\n",
      "print('{}\\n'.format(sqrtP))\n",
      "Prec = sqrtP.dot(sqrtP.T)\n",
      "print(Prec,\"\\n\")\n",
      "print((Prec==0.).sum())\n",
      "Prec2 = sqrtP.T.dot(sqrtP)\n",
      "print(Prec2,\"\\n\")\n",
      "print((Prec2==0.).sum())\n",
      "\n",
      "#print(lin.inv(Prec))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.2  0.   0.   0.   0.   0.   0. ]\n",
        " [ 0.1  0.6  0.   0.   0.   0.   0. ]\n",
        " [ 0.   1.6  0.1  0.   0.   0.   0. ]\n",
        " [ 0.1  1.3  0.   0.   0.   0.   0. ]\n",
        " [ 0.   0.   0.5  0.   1.4  0.   0. ]\n",
        " [ 0.1  0.4  1.4  0.   0.   0.9  0. ]\n",
        " [ 0.4  0.   0.   0.1  0.   0.   0.5]]\n",
        "\n",
        "[[ 0.1  0.   0.   0.   0.   0.   0.1]\n",
        " [ 0.   0.4  1.   0.8  0.   0.2  0. ]\n",
        " [ 0.   1.   2.5  2.   0.   0.7  0. ]\n",
        " [ 0.   0.8  2.   1.6  0.   0.5  0.1]\n",
        " [ 0.   0.   0.   0.   2.1  0.7  0. ]\n",
        " [ 0.   0.2  0.7  0.5  0.7  2.9  0. ]\n",
        " [ 0.1  0.   0.   0.1  0.   0.   0.5]] \n",
        "\n",
        "12\n",
        "[[ 0.3  0.2  0.1  0.1  0.   0.1  0.2]\n",
        " [ 0.2  4.6  0.6  0.   0.   0.3  0. ]\n",
        " [ 0.1  0.6  2.3  0.   0.7  1.2  0. ]\n",
        " [ 0.1  0.   0.   0.   0.   0.   0.1]\n",
        " [ 0.   0.   0.7  0.   1.8  0.   0. ]\n",
        " [ 0.1  0.3  1.2  0.   0.   0.8  0. ]\n",
        " [ 0.2  0.   0.   0.1  0.   0.   0.2]] \n",
        "\n",
        "20\n"
       ]
      }
     ],
     "prompt_number": 45,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 200\n",
      "M = sim_sqrtP(N, .7, scalef=1.)\n",
      "rk = np.linalg.matrix_rank(M)\n",
      "print(rk)\n",
      "(u, s, vt) = lin.svd(M, full_matrices=False)\n",
      "plt.plot(s[N-5:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "197\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 49,
       "text": [
        "[<matplotlib.lines.Line2D at 0x41d7a90>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEECAYAAADj+mWwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH1dJREFUeJzt3XlYVPX+B/D3DOCCyOKCGpAgqDCIgBtlomNS4oKZYj/R\njB8u2a1cyN+t2709N7zdJ8s2USu9ll6zrnYzM/c1x10UBTdQCSEB9wUUQVlmfn+cxhCBGWY758y8\nX89zHhEO57yfb/nhy3c5R6HT6UBERPKgFDsAEREZj0WbiEhGWLSJiGSERZuISEZYtImIZIRFm4hI\nRswu2hMnTlzarl27K2FhYSctESg2NnaLl5fXrbi4uPU1P5+XlxcQFRWV1rlz55yxY8euqqysdLHE\n/YiI5MTsop2UlLRsy5YtsZYIAwBvvvnm3BUrVkyo/fm33nrrw1mzZn2Sk5PT2cvL69bXX389yVL3\nJCKSC7OLdnR09F4vL69bNT+Xm5sbOGTIkM29evVK79+//56zZ892NfZ6Tz/99C9ubm6lNT+n0+kU\nu3btGhgfH78aABITE5evXbt2pLnZiYjkxtkaF3355Zf/tXjx4qlBQUG/pqWlRb366qtf7Ny5c5Cp\n17tx40ZrT0/PYqVSqQUAHx+foqKiIh/LJSYikgeLF+3S0lK3gwcPPjlmzJgf9J+rqKhoAgBr1qwZ\n9e67786u/T2+vr6FmzdvHmLpLERE9sbiRVur1So9PT2LMzIyImt/bdSoUWtGjRq1xtA1FArFQw9E\nad269Y3i4mJPrVarVCqV2sLCQl8fH58iS+YmIpIDiy/5c3d3vx0QEJC3evXqeEAYjz5x4kT3xlxD\np9Mpav5doVDoBg4cuOuHH34YAwDLly9PHDly5FrLpSYikgmdTtfgMW/evBndunU7GRoaemrevHkz\nan997NixKzt06HDRxcWlwtfXt2Dp0qVJeXl5/rGxsZvDw8MzVSrV6ffee+8dQ/fRH/369dvbtm3b\nq82bNy/z9fUt2LZt2zM6nQ7nz58P6NOnT1pQUFDOCy+88H1FRYWLsdfkwYMHD3s5FDpd/Y9mPXXq\nVLeEhISVR44c6e3i4lIZGxu7ZdGiRa8EBgbm2vDnChER/a7B4ZEzZ84ER0VFpTVr1uyek5NT9YAB\nA3avWbNmlK3CERHRwxos2t26dTu1d+/e6Js3b7YqKytz3bhx47DCwkJfW4UjIqKHNbh6JDg4+Mxb\nb7314bPPPrutRYsWdyMjIzP0a6X1aq/0ICIi4+hqLbowhsHVIxMnTlyanp7ea/fu3QM8PT2Lu3bt\neraOG0v+ePfdd0XPwJzMyZzMqD9MZXCd9tWrV729vb2vXrhw4fGffvrp+bS0tCiT70ZERGYxWLTj\n4+NX37hxo7WLi0vlF1988aq7u/ttWwQjIqJHGSzae/bs6W+LINamVqvFjmAU5rQs5rQsOeSUQ0Zz\nNLhO26gLKBQ6c69BRORoFAoFdNaYiCQiIulg0SYikhEWbSIiGWHRJiKSERZtIiIZYdEmIpIRFm0i\nIhlh0SYikhEWbSIiGWHRJiKSERZtIiIZYdEmIpIRFm0iIhmxSNG+c8cSVyEiIkMsUrSXLbPEVYiI\nyBCDRXvOnDlvh4aGng4LCzs5bty4/9y/f79p7XNSU4HqausEJCKiPzRYtPPz8/2XLFky5dixYz1O\nnjwZVl1d7bRq1aqxtc9r0wZYv956IYmISNBg0XZ3d7/t4uJSWVZW5lpVVeVcVlbm6uPjU1T7vORk\n4LPPrBeSiIgEDb4jslWrVjdnzZr1yeOPP36hefPm5YMHD94aExOzo/Z5WVkpyMwEpk4FEhLUdv+O\nNiKixtJoNNBoNGZfp8F3RObm5gbGxcWt37t3b7SHh0fJmDFjfoiPj189fvz47x5c4Pd3RH70EXDi\nBLBihdmZiIjsnlXeEZment6rb9++B1q3bn3D2dm5atSoUWsOHDjQt65zp0wBNm4Eih4ZPCEiIktp\nsGgHBwefOXTo0BPl5eXNdTqdYseOHTEqlSqrrnM9PYHx44HPP7dOUCIiMjA8AgBz5859c/ny5YlK\npVLbo0ePY1999dVkFxeXygcX+H14BAB+/RV48kngt98AV1frBicikjNTh0cMFm0jbqyreY3nngOG\nDAFeecWsyxIR2TXJFG2NRijYWVmAkk82ISKqk1UmIk0xYADQvDmwZYulr0xERBYv2goF8MYb3GxD\nRGQNFh8eAYCKCsDfX+htd+9u1uWJiOySZIZHAKBJE+C114B586xxdSIix2WVnjYAXL8OdO4MnDkD\ntGtn1i2IiOyOpHragPDkvxdeAL780lp3ICJyPFbraQNAdjYwcCCQnw80a2bWbYiI7IrketoAEBIC\nREYC//mPNe9CROQ4rL79Rb/8z8wOPRERwQZFOyZG+HPHI0/hJiKixrJ60VYogJkzudmGiMgSrDoR\nqXfvHtCxo/BckpAQs25HRGQXJDkRqdesmfAQqdRUW9yNiMh+2aSnDQBXrgDBwcIzt1u3NuuWRESy\nJ+meNiDsihw5Eli82FZ3JCKyPwaL9tmzZ7tGRkZm6A8PD4+S+fPnTzflZsnJwuvIKipM+W4iImrU\n8IhWq1X6+PgUHT58uI+fn18BYPzwiF5MDJCYCEyY0PiwRET2wibDIzt27IgJDAzM1RdsUyQnc7MN\nEZGpnBtz8qpVq8aOGzfukU3pKSkpDz5Wq9VQq9X1XmPIEGDWLGDPHuEtN0REjkCj0UCj0Zh9HaOH\nRyoqKpr4+PgUZWVlqdq2bXvtwQUaOTwCCE/+27oVWLu2cWGJiOyF1YdHNm/ePKRnz55HaxZsU730\nErB/v7D8j4iIjGd00V65cmVCQkLCSkvctEULYPJkYP58S1yNiMhxGDU8cvfu3RYdO3b8LS8vL6Bl\ny5Z3HrqACcMjAFBUBISFAefPA56ejf52IiJZM3V4xGY7Iuvy4otAeDjw5z+bFYGISHZkWbSPHgWe\nf17obTs3ah0LEZG8SX4be1169gT8/YEffxQzBRGRfIhatIE/NtsQEZFhohftESOAa9eAgwfFTkJE\nJH2iF20nJ2D6dPa2iYiMIepEpN6dO8LY9tGjwp9ERPZOlhORei1bAklJwIIFYichIpI2SfS0AeC3\n34AePYD8fKGIExHZM1n3tAHhxb+DBgFLl4qdhIhIuiTT0waEFSTjxwM5OcIEJRGRvZJ9TxsAnnwS\n8PYG1q0TOwkRkTRJqmgD3GxDRNQQyRXt0aOFScn0dLGTEBFJj+SKtrMzMG0ae9tERHWR1ESkXnEx\n0KkTcPIk4ONj0UsTEUmCXUxE6nl6Cs/aXrhQ7CRERNJisGgXFxd7xsfHrw4JCclWqVRZhw4desIW\nwWbMAL76Crh71xZ3IyKSB4NFe8aMGalDhw7dlJ2dHXLixInuISEh2bYIFhgIPPUU8M03trgbEZE8\nNDimXVJS4hEZGZlx/vz5TvVewApj2np79gBTpgDZ2YBSkgM5RESmMXVMu8GXfOXl5QW0bdv2WlJS\n0rLjx4+H9+zZ82hqauoMV1fXsprnpaSkPPhYrVZDrVY3NkedoqMBNzdg82Zg2DCLXJKISBQajQYa\njcbs6zTY005PT+/15JNPHjxw4EDf3r17H5k5c+Y8d3f32//4xz/+/uACVuxpA8C33wLLlgE7d1rt\nFkRENmeV1SO+vr6Fvr6+hb179z4CAPHx8auPHTvWw9SQpnjhBeDMGeD4cVvelYhImhos2u3bt7/s\n5+dXcO7cuS4AsGPHjpjQ0NDTtokmaNIEeO01YN48W96ViEiaDG6uOX78ePjkyZO/qqioaBIYGJi7\nbNmyJA8Pj5IHF7Dy8AgA3LgBBAUJE5Lt21v1VkRENmHq8Igkd0TW5ZVXgHbtgNmzrX4rIiKrs/ui\nfeYMMGCA8DCpZs2sfjsiIquyq23sdQkOBnr1Ar77TuwkRETikU3RBv541rYNOvZERJIkq6I9aJCw\nM3L7drGTEBGJQ1ZFW6EAZs7ks7aJyHHJZiJS7949wN8f2LULCAmx2W2JiCzK7ici9Zo1E5b/cbMN\nETki2fW0AeDKFWE1SU4O0KaNTW9NRGQRDtPTBoRNNqNGAYsWiZ2EiMi2ZNnTBoT3Rw4eDOTlAU2b\n2vz2RERmcaieNgCEhQGhocD334udhIjIdmRbtAFutiEixyProh0bC5SXA7t3i52EiMg2ZF20lUpu\ntiEixyLbiUi9sjJhs83+/UDnzqLFICJqFKs+mtXf3z/f3d39tpOTU7WLi0vl4cOH+9S4sahFGwD+\n9jegpARYuFDUGERERrNq0Q4ICMg7evRoz1atWt2s48aiF+2LF4WVJOfPA15eokYhIjKK1Zf8mXJx\nW3nsMWD4cGDJErGTEBFZl1FFW6FQ6GJiYnb06tUrfcmSJVOsHcoUycnAggVAZaXYSYiIrMfZmJP2\n79//VIcOHS5du3at7TPPPLM9ODj4THR09F7911NSUh6cq1aroVarLR7UkB49gE6dgB9/BMaOtfnt\niYgapNFooNFozL5Oo1ePzJ49+103N7fSWbNmfQJIY0xbb+1aYM4c4NAh4dnbRERSZbUx7bKyMtc7\nd+60BIC7d++22LZt27NhYWEnTQlpbXFxwI0bwMGDYichIrIOg8MjV65caff888//BABVVVXO48eP\n/+7ZZ5/dZv1ojefkBMyYAXz6KdC3r9hpiIgsT/aba2orLRU22xw5AgQEiJ2GiKhuDveUv/q4uQFJ\nScJKEiIie2N3PW0AuHABiIwUnrXt7i52GiKiR7GnXcPjjwMxMcDSpWInISKyLLvsaQPCsr9x44T3\nSDo5iZ2GiOhh7GnX8sQTQPv2wM8/i52EiMhy7LZoA8LW9k8/FTsFEZHl2HXRfv55oLBQWP5HRGQP\n7LpoOzsD06bxzTZEZD/sdiJSr6RE2GRz4gTg6yt2GiIiASci6+HhAUyYwLfaEJF9sPueNgDk5gJR\nUcBvvwEtWoidhoiIPe0GBQYC/fsDy5eLnYSIyDwO0dMGgL17gUmTgDNnAKVD/KgiIiljT9uAfv2E\n55Bs3Ch2EiIi0zlM0VYohM02XP5HRHLmMEUbAMaMAc6dAzIzxU5CRGQao4p2dXW1U2RkZEZcXNx6\naweypiZNgNdeA+bNEzsJEZFpjCraqampM1QqVZZCoZD+jKMBU6cKD5G6fFnsJEREjWewaBcWFvpu\n2rRp6OTJk78yZaZTalq1AhISgC++EDsJEVHjGSzaycnJn3300Ud/ViqVWlsEsoUZM4BFi4DycrGT\nEBE1ToNvY9+wYcNwb2/vq5GRkRkajUZd33kpKSkPPlar1VCr6z1VErp2Bfr0Ab79FpgyRew0ROQI\nNBoNNBqN2ddpcHPNX//61/dXrFgxwdnZuerevXvNbt++7T569Ogfv/nmm5ceXEAmm2tq27kTmD4d\nOHVKWA5IRGRLpm6uMXpH5O7duwd8/PHH/7d+/fq4WjeWZdHW6YCICGDuXGDwYLHTEJGjscmOSHtY\nPaKnUAAzZ3KzDRHJi8M8e6Qu9+4B/v7AL78AKpXYaYjIkfDZIyZo1gx49VVutiEi+XDonjYAXL0q\nrCY5dw5o21bsNETkKNjTNpG3NzB6tLBum4hI6hy+pw0Iy/6eeQbIzweaNhU7DRE5Ava0zdCtGxAW\nBqxaJXYSIqKGsWj/Tv+sbZn/0kBEdo5F+3eDBwP37wMW2GVKRGQ1LNq/Uyr5Zhsikj5ORNZQXg50\n7Ajs2wd06SJ2GiKyZ5yItIDmzYGXXwZSU8VOQkRUN/a0a7l0SdjSnpsrvDCBiMga2NO2kA4dgLg4\nYMkSsZMQET2KPe06ZGQAI0YA588DLi5ipyEie8SetgVFRgJBQcDq1WInISJ6GIt2PZKTgU8/5WYb\nIpIWFu16DB8OFBcD+/eLnYSI6A8Gi/a9e/eaRUVFpUVERGSqVKqst99+e44tgolNqRTe2s7NNkQk\nJUZNRJaVlbm6urqWVVVVOffr12/fxx9//H/9+vXbB9jnRKReaanwZpsjR4CAALHTEJE9sepEpKur\naxkAVFRUNKmurnZq1arVzcbeSI7c3ICJE4H588VOQkQkcDbmJK1Wq+zRo8ex3NzcwD/96U9fqlSq\nrJpfT0lJefCxWq2GWq22aEgxTZsGhIcDs2cD7u5ipyEiudJoNNBY4Il0jVqnXVJS4jF48OCtH3zw\nwV/UarUGsO/hEb2EBKBPH2FFCRGRJdhknbaHh0fJsGHDNqanp/dq7I3kLDlZGCKpqhI7CRE5OoNF\n+/r1622Ki4s9AaC8vLz59u3bn4mMjMywfjTp6NMHeOwxYO1asZMQkaMzOKZ96dKlDomJicu1Wq1S\nq9UqJ0yYsGLQoEE7bRFOSvTP2o6PFzsJETkyPnvESFVVwtb2//5X6HkTEZmDzx6xMmdnYPp0brYh\nInGxp90IJSXCJpvjxwE/P7HTEJGcsadtAx4eQGIisHCh2EmIyFGxp91IeXlA795Afr6wY5KIyBTs\nadtIQAAwYADw73+LnYSIHBF72ibYtw9ISgLOnhWeBkhE1FjsadvQU08Bnp7Ahg1iJyEiR8OibQKF\n4o/NNkREtsSibaIxY4CcHCAzU+wkRORIWLRN5OIiPLaVvW0isiVORJrh1i0gMBA4fRro0EHsNEQk\nJ5yIFIGXl/Cs7c8/FzsJETkK9rTNdO4c0K8f8NtvQPPmYqchIrlgT1skXboAUVHAihViJyEiR8Ci\nbQHJycC8eYAD/8JBRDbCom0BAwcKq0m2bhU7CRHZO4NFu6CgwG/gwIG7QkNDT3fr1u3U/Pnzp9si\nmJwoFMAbb3D5HxFZn8GJyMuXL7e/fPly+4iIiMzS0lK3nj17Hl27du3IkJCQbIATkXr37wP+/sD2\n7UC3bmKnISKps9pEZPv27S9HRERkAoCbm1tpSEhI9sWLFx8zJaQ9a9oUePVVYWybiMhaDL7Yt6b8\n/Hz/jIyMyKioqLSan09JSXnwsVqthlqttkg4uXnlFWE1yfvvA97eYqchIinRaDTQaDRmX8foddql\npaVuarVa88477/xz5MiRax9cgMMjD5kyRXgV2d//LnYSIpIyU4dHjCralZWVLsOHD98wZMiQzTNn\nznxoAIBF+2GnTwMxMcKbbZo2FTsNEUmV1ca0dTqdYtKkSV+rVKqs2gWbHhUaCnTvDqxcKXYSIrJH\nBnva+/bt69e/f/893bt3P6FQKHQAMGfOnLdjY2O3AOxp12XrVuDNN4XHtioa/XOUiByBVYdHDNyY\nRbsWnU5Y9rdgAfD002KnISIp4rNHJEShAGbOBD79VOwkRGRv2NO2kvJyoGNHYO9eoGtXsdMQkdSw\npy0xzZsDU6cCqaliJyEie8KethVdugSoVEBuLtCqldhpiEhK2NOWoA4dgBEjgH/9S+wkRGQv2NO2\nssxMYPhwIC9PeHwrERHAnrZkRUQIzyP54QexkxCRPWDRtoHkZGH5H38hISJzsWjbwLBhwO3bwL59\nYichIrlj0bYBpRKYMYNvtiEi83Ei0kZKS4U32xw+DHTqJHYaIhIbJyIlzs0NmDQJmD9f7CREJGfs\nadtQYaHw2Na8PMDDQ+w0RCQm9rRlwNcXiI0Fvv5a7CREJFfsadvYkSNAfLywtd25UW/oJCJ7wp62\nTPTuLbxD8qefxE5CRHJksGhPnDhxabt27a6EhYWdtEUgR5CczOV/RGQag0U7KSlp2ZYtW2JtEcZR\njBwpPAEwLU3sJEQkNwaLdnR09F4vL69btgjjKJycgOnT2dsmosazyFRYSkrKg4/VajXUarUlLmvX\nJk0C/vlP4MIF4PHHxU5DRNam0Wig0WjMvo5Rq0fy8/P94+Li1p88eTLskQtw9YjJ3nhDWEEyd67Y\nSYjI1rh6RIamTxfWbJeWip2EiOSCRVtE/v7AwIHAsmViJyEiuTBYtBMSElb27dv3wLlz57r4+fkV\nLFu2LMkWwRxFcrLw8t/qarGTEJEccEekyHQ6ICoK+NvfgOeeEzsNEdkKx7RlSqHgZhsiMh6LtgTo\nn0WSkSF2EiKSOhZtCXBxAaZNY2+biAzjmLZE3LolvNHm9GngscfETkNE1mbqmDaLtoS8/jpw9Srw\nv/8LqFTCTkklfxciskss2nbg8mXgww+F3nZWFlBcDAQHCwW85hEQIDy/hIjki0XbDpWUANnZQgGv\neVy9CnTp8kcRDwkR/gwKEsbHiUj6WLQdSGkpcPbso8W8oAAIDHy0Z96lC9C0qdipiagmFm1CeTlw\n7tyjxTwvD+jY8dFi3rUr4Ooqdmoix8SiTfWqqABych4t5r/+KqxUqV3MQ0IANzexUxPZNxZtarSq\nKmFTT+1ifvYs0LZt3cXc01Ps1ET2gUWbLKa6GsjPf7SYZ2cDHh6PFnOVCmjdWuzURPLCok1Wp9UK\nk521i3lWFtCsWd3F3NtbeL4KET2MRZtEo9MBFy8+WshPnxYKdl3F/LHHWMzJsbFok+TodMKa8rp6\n5vfu1V3M/fy4C5Qcg9Uezbply5bY4ODgM507d8758MMP3zItnvgs8UJNW7CnnAoF0K6d8Hae114D\nPv8c2LULuHJFmACdO1d4lnhBgfCwrKeeAtzdgd69gcRE4IMPgHXrhFUupr4kwp7aUwrkkFMOGc3R\nYNGurq52ev311xdu2bIlNisrS7Vy5cqE7OzsEFuFsyS5/Id0lJxt2gDR0cDLLwPz5gHbtgGFhUBR\nEbBgATBgAHD9OrBoERATA7RsCUREAOPGCW+xX7MGOHMGqKy0bk5bYU7LkUNGczg39MXDhw/3CQoK\n+tXf3z8fAMaOHbvq559/fi4kJCTbJunI4Xh4AE88IRw1lZYKRVo/vPLvfwt/FhXVvQu0c2fuAiX7\n1GDRLioq8vHz8yvQ/93X17cwLS0tyvqxiB7m5gb06iUcNZWXP7ylf9Uq4c/8fOHFyWVlwPbtD39P\nXROgtT9nzDmW/L68PGDfPmnlrOucnBzg6NFHPy8l585JP6NZdDpdvcfq1atHT548eYn+7ytWrHjx\n9ddfX1DzHAA6Hjx48ODR+KOh+lvf0WBP28fHp6igoMBP//eCggI/X1/fwprnmDL7SUREpmlwIrJX\nr17pOTk5nfPz8/0rKiqafP/99/8zYsSIdbYKR0RED2uwp+3s7Fy1cOHC1wcPHry1urraadKkSV9z\nEpKISETGjqNs3rw5tmvXrmeCgoJyPvjgg7fqOmfatGnzg4KCcrp373782LFjkaaM15h7GMq5a9cu\ntbu7e0lERERGRERExnvvvfeOrTMmJSUt9fb2vtKtW7eT9Z0jhbY0lFMKbanT6XDhwgU/tVq9S6VS\nnQ4NDT2Vmpo6XYptakxOsdu0vLy8WZ8+fdLCw8MzQ0JCsv7yl7/MkWJbGpNT7LaseVRVVTlFRERk\nDB8+fL257Wn0DQMDA3/Ny8vzr6iocAkPD8/MysoKqXnOxo0bhw4ZMmSTTqfDoUOHoqKiog6J0TCG\ncu7atUsdFxe3Tqz/eDqdDnv27Ik+duxYZH3FUAptaUxOKbSlTqfDpUuX2mdkZETodDrcuXPHrUuX\nLmel+P+nMTml0KZ379511el0qKysdI6Kijq0d+/eflJrS2NySqEt9ccnn3zyxrhx476rK09j29Oo\nDcM112u7uLhU6tdr1zxn3bp1IxITE5cDQFRUVFpxcbHnlStX2lnjtwNzcgLiT55GR0fv9fLyulXf\n16XQloDhnID4bQkA7du3vxwREZEJAG5ubqUhISHZFy9efOid9lJoU2NyAuK3qauraxkAVFRUNKmu\nrnZq1arVzZpfl0JbGpMTEL8tAaCwsNB306ZNQydPnvxVXXka255GFe261msXFRX5GDqnsLDQ15jr\nW4oxORUKhe7AgQN9w8PDjw8dOnRTVlaWypYZjSGFtjSGFNsyPz/fPyMjIzIqKiqt5uel1qb15ZRC\nm2q1WmVERERmu3btrgwcOHCXSqXKqvl1qbSloZxSaEsASE5O/uyjjz76s1Kp1Nb19ca2p1FFW6FQ\n6Iw5r/ZPEWO/z1KMuV+PHj2OFRQU+B0/fjx82rRpC0aOHLnWFtkaS+y2NIbU2rK0tNQtPj5+dWpq\n6gw3N7fS2l+XSps2lFMKbapUKrWZmZkRhYWFvnv27Omv0WjUtc+RQlsayimFttywYcNwb2/vq5GR\nkRkN9fob055GFW1j1mvXPqewsNDXx8enyJjrW4oxOVu2bHlH/2vVkCFDNldWVrrcvHmzlS1zGiKF\ntjSGlNqysrLSZfTo0T+++OKL39b1j1MqbWoop5Ta1MPDo2TYsGEb09PTH9qHKpW21KsvpxTa8sCB\nA33XrVs3IiAgIC8hIWHlL7/88vRLL730Tc1zGt2exgyiV1ZWOnfq1Ck3Ly/P//79+00MTUQePHjw\nCTEmJ4zJefny5XZarVah0+mQlpbWp2PHjvliTEzk5eX5GzMRKVZbGpNTKm2p1WoVEyZM+GbmzJmf\n1XeOFNrUmJxit+m1a9fa3Lp1y1On06GsrKx5dHT0nh07dgySWlsak1Pstqx9aDSaAXWtHmlseza4\nTluvvvXaixcvngoAU6dOXTx06NBNmzZtGhoUFPRrixYt7i5btizJxB9OJjMm5+rVq+O//PLLPzk7\nO1e5urqWrVq1aqytcyYkJKzcvXv3gOvXr7fx8/MrmD179ruVlZUu+oxSaEtjckqhLQFg//79T337\n7bcvdu/e/URkZGQGALz//vt/vXDhwuP6rFJoU2Nyit2mly5d6pCYmLhcq9UqtVqtcsKECSsGDRq0\nU2r/1o3JKXZb1kU/7GFOe5r9EgQiIrIdviOEiEhGWLSJiGSERZuISEZYtImIZIRFm4hIRli0iYhk\n5P8B4TMxBlEO3moAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x39d7c50>"
       ]
      }
     ],
     "prompt_number": 49,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "print(s[N-5:])\n",
      "plt.plot(u[1:15,N-1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[  5.9e-04   2.6e-05   2.8e-07   8.7e-11   3.7e-15]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "[<matplotlib.lines.Line2D at 0x3ff6690>]"
       ]
      }
     ],
     "prompt_number": 19,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nulld = u[:,-1]\n",
      "print(nulld[:10])\n",
      "print(M.T.dot(nulld).T)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[  9.8e-01  -2.0e-01  -1.7e-05   3.9e-02  -3.0e-02   6.8e-03  -5.8e-05\n",
        "   3.8e-04  -1.4e-03   7.5e-04]\n",
        "[  5.3e-15  -1.2e-15   4.4e-15  -1.4e-15   4.7e-16  -2.3e-15  -1.1e-15\n",
        "  -2.5e-15  -1.7e-16  -1.1e-15   1.4e-16   3.4e-16  -1.5e-15   4.0e-16\n",
        "   2.2e-15   7.7e-16  -1.2e-16  -7.9e-16   1.9e-15   9.8e-16   4.0e-16\n",
        "  -3.4e-16  -3.9e-15  -1.4e-16  -1.6e-15   4.2e-17   2.9e-16  -2.7e-15\n",
        "  -7.0e-16  -1.7e-15  -2.2e-15  -2.7e-16   8.1e-16   1.1e-15  -8.9e-16\n",
        "  -9.7e-17   2.2e-15   7.2e-16  -1.8e-15  -1.8e-16  -1.9e-15   1.7e-16\n",
        "  -1.4e-16   6.1e-16   3.2e-16   8.6e-16  -1.0e-15  -1.2e-15  -2.9e-16\n",
        "  -2.2e-16  -4.8e-16  -1.4e-15  -9.3e-16   1.2e-16   1.2e-15   9.2e-17\n",
        "   1.6e-16   9.4e-16  -9.8e-16  -5.1e-16  -4.8e-16   1.3e-15  -2.9e-16\n",
        "  -6.0e-16  -5.4e-16  -5.8e-17  -1.0e-15   4.0e-16  -2.1e-16  -3.9e-16\n",
        "   4.2e-16   4.7e-16   3.9e-16  -3.6e-16   6.6e-16   5.5e-16   1.1e-16\n",
        "   7.6e-16   4.6e-16   2.0e-16   7.0e-16   4.0e-16   4.6e-16  -1.7e-17\n",
        "  -1.7e-16  -1.0e-16  -3.1e-16  -1.0e-16  -3.6e-16   4.9e-16  -2.4e-16\n",
        "  -1.7e-16  -9.1e-19  -6.4e-17  -1.0e-16   6.1e-17   1.2e-16  -6.1e-18\n",
        "  -7.6e-17  -9.8e-17]\n"
       ]
      }
     ],
     "prompt_number": 20,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "prng = np.random.RandomState(1)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 170,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Generate the data\n",
      "n_samples = 1000\n",
      "n_features = 500\n",
      "\n",
      "prec = make_sparse_spd_matrix(n_features, alpha=.5,\n",
      "                              smallest_coef=.4,\n",
      "                              largest_coef=.7,\n",
      "                              random_state=prng)\n",
      "np.linalg.matrix_rank(prec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 208,
       "text": [
        "493"
       ]
      }
     ],
     "prompt_number": 208,
     "trusted": true
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "trusted": true
    }
   ],
   "metadata": {}
  }
 ]
}