jgomezdans · September 18, 2015 06:30
diff --git a/HessianCalculations.ipynb b/HessianCalculations.ipynb
 {
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating the uncertainty in a calibration scenario\n",
    "\n",
    "After you minimize a cost function $J(\\vec{x})$ (where $\\vec{x}$ is the vector of your parameters), you end up having what is called the maximum *a posteriori* (MAP) estimate of the parameters, $\\vec{x}_{MAP}$. If the model is linear and the statistics are normal (Gaussian), then the *a posteriori* distribution is effectively Gaussian, with a mean given by $\\vec{x}_{MAP}$ and an uncertainty given by the inverse of the Hessian calculated at $\\vec{x}_{MAP}$, $\\mathbf{C}_{posterior}$. \n",
    "\n",
    "In the case of DALEC, the model isn't linear, but the statistics are Gaussian, so we pay a price for this convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.optimize as opt\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "def calculate_hessian ( func, x, epsilon=1e-8, args=()):\n",
    "    \"\"\"This function calculates the Hessian of ``func`` at point ``x`` using \n",
    "    finite differences. You might need to change the step size to different\n",
    "    valuews (from 1e-3 to 1e-10) if the solution is unstable and you get NaNs\n",
    "    \"\"\"\n",
    "    # First, we calculate the Jacobian using finite differences:\n",
    "    jac = opt.approx_fprime ( x, func, epsilon, *args )\n",
    "    npar = len( x )\n",
    "    hessian = np.zeros (( npar, npar))\n",
    "    xp = 1.*x\n",
    "    for i in xrange(npar):\n",
    "        xx0 = 1.*x[i]\n",
    "        xp[i] = xx0 + epsilon\n",
    "        jac1 = opt.approx_fprime ( xp, func, epsilon, *args )\n",
    "        hessian[i, :] = ( jac1 - jac)/epsilon\n",
    "        xp[i] = xx0\n",
    "    return hessian\n",
    "\n",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We now try an example function, the [Rosenbrock function](https://en.wikipedia.org/wiki/Rosenbrock_function), a function commonly used in optimisation problems. There is an implementation in ``scipy.optimize``, and it is useful because we have an analytical version of the Hessian to compare against (given by [scipy.optimize.rosen_hess](http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.optimize.rosen_hess.html). We will calculate the Hessian for the five dimensional version of the Rosenbrock function at $[0.9, 1.0, 0.2, 0.1, -1]^\\top$, and compare the result with the analytic version:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_map point is [ 0.9  1.   0.2  0.1 -1. ]\n",
      "Numerical approximation to the Hessian:\n",
      "[[  574.0216  -360.0019    -0.0003     0.0003     0.    ]\n",
      " [ -360.0019  1322.0239  -400.002      0.         0.    ]\n",
      " [   -0.0003  -400.002    210.0043   -80.0023     0.    ]\n",
      " [    0.0003     0.       -80.0023   614.0019   -40.0019]\n",
      " [    0.         0.         0.       -40.0019   200.0004]]\n",
      "Analytic Hessian is:\n",
      "[[  574.  -360.     0.     0.     0.]\n",
      " [ -360.  1322.  -400.     0.     0.]\n",
      " [    0.  -400.   210.   -80.     0.]\n",
      " [    0.     0.   -80.   614.   -40.]\n",
      " [    0.     0.     0.   -40.   200.]]\n",
      "Are they similar? If so the norm of the difference vs the true norm should be very small: 1.86707213248e-05\n"
     ]
    }
   ],
   "source": [
    "# Change the way we print numbers to the screen so it's easier to see\n",
    "np.set_printoptions (linewidth=132, precision=4, suppress=True )\n",
    "# Set the point around which we'll calculate the Hessian\n",
    "x_map = np.array([0.9, 1.0, 0.2, 0.1, -1.])\n",
    "print \"x_map point is\", x_map\n",
    "print \"Numerical approximation to the Hessian:\"\n",
    "print calculate_hessian ( opt.rosen, x_map, epsilon=1e-5)\n",
    "print \"Analytic Hessian is:\"\n",
    "print opt.rosen_hess ( x_map )\n",
    "print \"Are they similar? If so the norm of the difference vs the true norm should be very small:\", \\\n",
    "        np.linalg.norm ( calculate_hessian ( opt.rosen, x_map, epsilon=1e-5)- \\\n",
    "        opt.rosen_hess ( x_map ))/np.linalg.norm(opt.rosen_hess(x_map))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the code above shows that the Hessian works. In order to get the parameter uncertainties, you would \n",
    "\n",
    "1. Invert the hessian to obtain $\\mathbf{C}_{post}$,\n",
    "2. Extract the main diagonal elements from $\\mathbf{C}_{post}$, and take the square root to find the standard deviation for each parameters\n",
    "3. Maybe plot $\\mathbf{C}_{post}$, to see the correlations between parameters etc.\n",
    "\n",
    "**NOTE** Sometimes, the calculation of the Hessian is unstable, and you end up with a Hessian with stupid values such as ``NaN`` or somesuch. Some other times, the numerical errors result in a matrix that cannot be inverted as we do below with ``np.linalg.inv``. If that happens, try other values of ``epsilon``, as mentioned above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "STandard deviations:  [ 0.0555  0.0583  0.1367  0.0445  0.0713]\n"
     ]
    }
   ],
   "source": [
    "hess = calculate_hessian ( opt.rosen, x_map, epsilon=1e-5)\n",
    "post_covariance = np.linalg.inv ( hess )\n",
    "sigma = np.sqrt ( post_covariance.diagonal())\n",
    "print \"STandard deviations: \", sigma "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's useful to also plot the correlation matrix between the different parameters. Parameters that are strongly correlated (either positively or negatively) are parameters that compensate for each other and an indication that their values cannot be obtained from the calibration directly, but must be helped out with priors, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.      0.6591  0.5638  0.2289  0.0286]\n",
      " [ 0.6591  1.      0.8555  0.3474  0.0433]\n",
      " [ 0.5638  0.8555  1.      0.4061  0.0507]\n",
      " [ 0.2289  0.3474  0.4061  1.      0.1248]\n",
      " [ 0.0286  0.0433  0.0507  0.1248  1.    ]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdIAAAHVCAYAAABFZXqOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG99JREFUeJzt3X+M5PV93/HnzoZzMexdUwXSqLHiRvnhVA5O3WsMdwif\ny10SqE8ltt+2wLHhAur6aC2gf9jEIbB2iVzFOds4DohczzoojRO9A8Qidq70atbgbYJInEIS2T2B\nWjWKHAqO8V4wnM3O9I+dtafD7e53Zz6z8539Ph/SyPv9fme+3/cs1r339fl+vt8vSJIkSZIkSZIk\nSZIkSZIkSZIkSbU1Ne4CJEkahYhoAbcD5wGngGsy86me7TcAVwPPdFfNZuaJjR7newrUKklSHV0G\nbMvMXRHxBuBQd92K1wPvysw/G+YgrWE+LElSje0GjgFk5qPAzr7t/wz4QEQ8EhE3DnoQG6kkaava\nDiz2LC91h3tXfBqYBf4FcGFE/MtBDuLQriSpmM/d95nOtu1nFd9vu90+efjw4Y/2rZ7PzPk1PrYI\nzPQstzKz3bN8W2YuAkTEZ4F/Cnx2o7XZSCVJxWzbfha59+ri+43jR2Yyc26DH1sA9gMZEecDT3xn\nfxE7gCci4p8A32Q5lR4ZpDYbqSSpqOn6XA9yP7AvIha6ywci4nLg7Mw83D0v+hDLM3qPZ+axQQ5i\nI5UkbUmZ2QEO9q0+0bP90yyfJx2KjVSSVNT0VH0i6WZw1q4kSUMwkUqSiqrROdJNYSOVJBXl0K4k\nSarMRCpJKqppQ7smUkmShmAilSQV5TlSSZJUmYlUklRU086R2kglSUU5tCtJkiozkUqSimpaQmva\n95UkqSgTqSSpKM+RSpKkykykkqSivPxFkqQhOLQrSZIqM5FKkopq2tCuiVSSpCGYSCVJRXmOVJIk\nVWYilSQV1bRzpDZSSVJRDu1KkqTKTKSSpKKaNrRrIpUkaQgmUklSUSZSSZJUmYlUklRU02bt2kgl\nSUU5tCtJkiozkUqSimra0K6JVJKkIZhIJUlFeY60ASJiz7hrqCt/N2vz97M6fzdr8/ezdTU1ke4B\n5sdcQ13twd/NWvbg72c1e/B3s5Y9NOT3U5dzpBHRAm4HzgNOAddk5lOned9vAV/LzF8a5DiNTKSS\npNGZnir/GtBlwLbM3AXcCBzqf0NEzAKvBTqDHsRGKknaqnYDxwAy81FgZ+/GiNgF/DRwJzBwu7aR\nSpKKmp6aKv4a0HZgsWd5qTvcS0T8AHAz8G8Zooky7Ier+u1Dn+yc+7rXbMahJEnr2H7WK/npC3aN\n5N//48ePd57/heuK7/ese27jzjvv/GDf6vnMnF/tMxFxCPjjzMzu8l9l5qu6P78XuBI4CfxD4JXA\nr2Tm3RutbVMmG537uteQe6/ejENJjXLrx9867hJqa+bd7xt3CbX1xS/9xUj33xrRZKPMnNvgRxaA\n/UBGxPnAEz37+g3gNwAi4krgNYM0UWjurF1J0ohM1edC0vuBfRGx0F0+EBGXA2dn5uG+9w482chG\nKknakjKzAxzsW33iNO+7a5jj2EglSUW16pNIN4WzdiVJGoKJVJJU1NR0szJas76tJEmFmUglSUXV\naNbuprCRSpKKcrKRJEmqzEQqSSpqqtWsjNasbytJUmEmUklSUaM4Rzrw/fs2gYlUkqQhmEglSUWN\n4vKXOidSG6kkqSjvbCRJkiozkUqSihrFZKOl4nssx0QqSdIQTKSSpKKmWt4iUJIkVWQilSQV1WrY\nrF0bqSSpqKY9Rq1ZfzZIklSYiVSSVJSJVJIkVWYilSQV1bTJRs36tpIkFWYilSQV1bRzpDZSSVJR\nLe9sJEmSqlo3kUZEC7gdOA84BVyTmU+NujBJ0mTyeaQvdxmwLTN3ATcCh0ZbkiRJk6NKI90NHAPI\nzEeBnSOtSJI00VrTU8VfdVZlstF2YLFneSkiWpnZPt2bI2IPsGf40iRJoxIRc32r5jNzfgylTLwq\njXQRmOlZXrWJAnT/Q8z3rpudnb1lkOIkSaORmXOj2reXv7zcArAfyIg4H3hitCVJkiZZ0yYbVWmk\n9wP7ImKhu3xghPVIkjRR1m2kmdkBDm5CLZKkLaDuk4NK885GkqQtab37IETEW4H3Ax3gP2fmJwY5\nTrMGsiVJIzfVmir+GtCq90GIiGngw8DFwAXAtRHxDwY5iI1UkrRVrXofhMxcAl6TmSeBc4Bp4FuD\nHMRGKkkqqjXdKv4a0Gnvg7CykJntiHgL8GfAQ8A3BzmI50glSUWN6jrSAW4ise59EDLzvoi4HzgK\nvLv7vxtiI5UkTYQBbiKx6n0QImI78ACwLzO/FRHPA0uD1GUjlSQVVaMbMrzsPggRcTlwdmYejoh7\ngIcj4tvA48A9gxzERipJ2pJWuQ/CiZ7th4HDwx7HRipJKmqqVZtEuima9W0lSSrMRCpJKmqIy1Um\nko1UklRUjSYbbYpmfVtJkgozkUqSijKRSpKkykykkqSivPxFkiRVZiKVJBU1NT097hI2lY1UklSU\nk40kSVJlJlJJUlEtJxtJkqSqTKSSpKI8RypJkiozkUqSimpaIrWRSpKK8s5GkiSpMhOpJKmopg3t\nNuvbSpJUmIlUklRU0xKpjVSSVFTLRippUpy477Fxl1BbO684Ne4S1BA2UklSUV7+IkmSKjORSpKK\natpko2Z9W0mSCjORSpKKaloitZFKkopyspEkSarMRCpJKqo1PT3uEjaViVSSpCGYSCVJRTVtslGz\nvq0kSYWZSCVJRTUtkdpIJUlFNe3yFxupJGlLiogWcDtwHnAKuCYzn+rZfjlwHfAS8OfAtZnZ2ehx\nmvVngyRp5KamW8VfA7oM2JaZu4AbgUMrGyLiTODfA3sy80JgB/DmQQ5iI5UkbVW7gWMAmfkosLNn\n24vABZn5Ynf5e4AXBjmIQ7uSpKJqNNloO7DYs7wUEa3MbHeHcJ8BiIj3Amdl5vFBDmIjlSRNhIiY\n61s1n5nza3xkEZjpWW5lZrtnfy3g14AfAd46aF02UklSUaOatZuZcxv8yAKwH8iIOB94om/7nSwP\n8f78IJOMVthIJUlFTbVqc6/d+4F9EbHQXT7Qnal7NvAnwC8CDwOfjwiA2zLz9zd6EBupJGlL6qbM\ng32rT/T8XKTj20glSWXVJ5FuitpMrZIkaRKZSCVJZTXsFoHN+raSJBVmIpUkFTU13axzpDZSSVJZ\nTjaSJElVmUglSWWZSCVJUlUmUklSUaO6125dbejbRsQbIuKhURUjSdKkqZxII+J9wC8Afze6ciRJ\nE69h50g3MrT7JPAW4D+NqBZJ0lbQsEZaeWg3M+8DXhphLZIkTZzik40iYg+wp/R+JUnlRMRc36r5\nzJwvse+mTTYq3ki7/yHme9fNzs7eUvo4kqTBZebcuGvYKgZppJ3iVUiSto5RnCOtcefZUCPNzP8N\n7BpNKZIkTR5vyCBJKmsUiXSp/C5LsZFKkooayWPUatxImzW1SpKkwkykkqSyGnb5S7O+rSRJhZlI\nJUlleYtASZJUlYlUklTUVMMSqY1UklSWk40kSVJVJlJJUlFNG9o1kUqSNAQTqSSpLBOpJEmqykQq\nSSqrYbN2baSSpKJG8vSXGmvWnw2SJBVmIpUklTWSyUb1fSCpiVSSpCGYSCVJZdUkkUZEC7gdOA84\nBVyTmU/1veeVwH8FfjEz/+cglZlIJUlb1WXAtszcBdwIHOrdGBE7gYeBfwx0Bj2IjVSSVNRUq1X8\nNaDdwDGAzHwU2Nm3fRvLzXagJLrCoV1JUlkjurNRRMz1rZrPzPk1PrIdWOxZXoqIVma2ATLzv3f3\nO1RdNlJJ0kTIzLkNfmQRmOlZ/k4TLcmhXUlSWVOt8q/BLACXAkTE+cATpb5iLxOpJGmruh/YFxEL\n3eUDEXE5cHZmHi51EBupJKmswRNkUZnZAQ72rT5xmve9aZjj2EglSUV1atJIN0uzvq0kSYWZSCVJ\nZZlIJUlSVSZSaYL9zVe+Nu4S6qtT/HJBVTU1Ne4KNpWJVJKkIZhIJUllDX5v3IlkI5UkFeXlL5Ik\nqTITqSSpLBOpJEmqykQqSSrLRCpJkqoykUqSympYIrWRSpKK8vIXSZJUmYlUklSWiVSSJFVlIpUk\nleXTXyRJUlUmUklSWQ07R2ojlSQV5eUvkiSpMhOpJKmshj3Yu1nfVpKkwkykkqSyPEcqSZKqMpFK\nkspqWCK1kUqSympYI23Wt5UkqTATqSSpKG/IIEmSKjORSpLKMpFKkqSqTKSSpLIa9jxSG6kkqayG\nDe2u2Ugj4gzgU8APAa8Abs3MBzajMEmShhERLeB24DzgFHBNZj7Vs30/8CvAS8CnMvM/DnKc9f5s\neCfwTGZeBPwc8MlBDiJJao7OVKv4a0CXAdsycxdwI3BoZUM3KH4U2Ae8EfjXEXHuIAdZr7oEbu55\n70uDHESSpDHYDRwDyMxHgZ09234CeDIzv5GZ3wa+CFw0yEHWHNrNzOcBImKG5ab6y4McRJLUIPU5\nR7odWOxZXoqIVma2u9u+0bPtJLBjkIOsO9koIl4F3Af8Zmb+ToX37wH2DFKMJGlzRMRc36r5zJwf\nQymVDVDzIjDTs7zSRGG5ifZumwG+Pkhd6002+n7gQeDazHyoyg67X2q+d93s7OwtgxQnSRqNzJwb\n1b47I7r8ZYCaF4D9QEbE+cATPdu+AvxoRHwv8DzLw7ofGaSu9RLpB1iOujdHxMq50ksy88VBDiZJ\n2vo6nXFX8B33A/siYqG7fCAiLgfOzszDEfHvgP/C8hygI5n51UEOst450uuA6wbZsSRJ45SZHeBg\n3+oTPdv/APiDYY/jDRkkSUW1axRJN0NtplZJkjSJTKSSpKKalUdNpJIkDcVEKkkqqt2wSGojlSQV\n1XGykSRJqspEKkkqqmlDuyZSSZKGYCKVJBXVsEBqIpUkaRgmUklSUU07R2ojlSQV5eUvkiSpMhOp\nJKmo9rgL2GQmUkmShmAilSQV1bBTpCZSSZKGYSKVJBXl5S+SJA3By18kSVJlJlJJUlFe/iJJkioz\nkUqSimrYKVIbqSSprHbDOqlDu5IkDcFEKkkqqll51EQqSdJQTKSSpKKadmcjE6kkSUMwkUqSimrY\npF0bqSSprHbDphttWiO99eNv3axDTZQT9z027hJq7W++8rVxl1BrD/7f58ddQm3t/+svj7uEGts2\n7gK2FBOpJKmopg3tOtlIkqQhmEglSUV5+YskSarMRCpJKqrO50gj4kzgHuAc4CRwZWY+e5r3nQMs\nAK/NzG+ttU8TqSSpqDad4q+CDgKPZ+ZFwN3ATf1viIifBR4Ezq2yQxupJKlJdgPHuj8fA/ae5j1L\nwMXA16vs0KFdSVJRdRnajYirgev7Vj8NLHZ/Pgns6P9cZh7vfr7ScWykkqSJEBFzfavmM3N+tfdn\n5hHgSN8+7gVmuoszwHPD1mUjlSQV1R5RJM3MuQK7WQAuBR4DLgEeHnaHNlJJUpPcAdwVEY8Ap4Ar\nACLiBuDJzHyg572V/iKwkUqSilpqj7uC1WXmC8DbT7P+Y6dZ98NV9mkjlSQVNaqh3bry8hdJkoZg\nIpUkFbVkIpUkSVWZSCVJRXmOVJIkVWYilSQVVefLX0bBRipJKsqhXUmSVJmJVJJUlJe/SJKkykyk\nkqSi2iMIpFPld1mMiVSSpCGYSCVJRS2NIJLWuVnVuTZJ0gTy8hdJklSZiVSSVNRSswKpiVSSpGGY\nSCVJRTXtHOm6jTQipoHDwI8BHeA9mfmXoy5MkqRJUGVo981AOzMvBG4CfnW0JUmSJtlSu1P8VWfr\nNtLM/Aww2118NfD1URYkSZps7U6n+KvOKp0jzcyliDgK/DzwtpFWJEnSBKk82Sgzr4qI9wOPRsRP\nZOYLp3tfROwB9pQpT5I0ChEx17dqPjPnS+y7aZe/VJls9C7gBzPzw8ALQLv7Oq3uf4j53nWzs7O3\nDFWlJKmozJwbdw1bRZVE+nvA0Yj4AnAGcF1mnhptWZKkSVX3c5qlrdtIu0O479iEWiRJmjjekEGS\nVFS75perlGYjlSQV1bTJRt5rV5KkIZhIJUlFNW2ykYlUkqQhmEglSUUtmUglSVJVJlJJUlFe/iJJ\n0hC8/EWSJFVmIpUkFVXny18i4kzgHuAc4CRwZWY+2/eeG/jurXE/l5kfWmufJlJJUpMcBB7PzIuA\nu4GbejdGxA8DVwAXZOb5wM9ExE+utUMbqSSpqKVOp/iroN3Ase7Px4C9fdv/D/Czmbly0DNYfoTo\nqhzalSRNhI0+jDwirgau71v9NLDY/fkksKN3Y2a+BPxtREwBHwG+lJlPrlWXjVSSVNTSiC5/2ejD\nyDPzCHCkd11E3AvMdBdngOf6PxcRfw/4FPAN4Nr1jmMjlSQVNapGWsgCcCnwGHAJ8HDvxm4S/Qzw\n3zLz16rs0EYqSWqSO4C7IuIR4BTLE4tWZuo+CUwDFwFnRMQl3c/8Umb+8Wo7tJFKkoqqcyLNzBeA\nt59m/cd6Fs/cyD6dtStJ0hBMpJKkouqcSEfBRCpJ0hBMpJKkopqWSG2kkqSiRtJIazx+WuPSJEmq\nPxOpJKkoE6kkSarMRCpJKsrJRpIkDaFpjdShXUmShmAilSQVZSKVJEmVmUglSUW91LBEummNdObd\n79usQ02UnVecGncJ9dZpj7uCWtv/118edwm19d7XHxx3CbUVx4+Mu4QtxUQqSSqqaedIbaSSpKKa\n1kidbCRJ0hBMpJKkopY6JlJJklSRiVSSVJTnSCVJUmUmUklSUU1LpDZSSVJRTWukDu1KkjQEE6kk\nqaildrNu7WkilSRpCCZSSVJRniOVJEmVmUglSUU1LZHaSCVJRTXtwd4O7UqSNAQTqSSpqKYN7ZpI\nJUkagolUklSUiVSSJFVmIpUkFVXnRBoRZwL3AOcAJ4ErM/PZvvf8G+BKoAP8embmWvs0kUqSilpq\nd4q/CjoIPJ6ZFwF3Azf1boyI7wPeA1wAXAwcWm+HNlJJUpPsBo51fz4G7O3d2E2nr8vMJeAHgBfX\n26FDu5KkouoytBsRVwPX961+Gljs/nwS2NH/ucxsd4d3Pwjctt5xbKSSpIkQEXN9q+Yzc36192fm\nEeBI3z7uBWa6izPAc6t89jcj4reAP4yIR9Y6jo1UklRUZ0SJNDPnCuxmAbgUeAy4BHi4d2NE/Djw\n4cx8C/AScApYWmuHNlJJUpPcAdwVEY+w3CSvAIiIG4AnM/OBiPgfEfFHLM/a/VxmPrLWDm2kkqSi\n2jU5R3o6mfkC8PbTrP9Yz88fAj5UdZ82UklSUZ1OfRvpKFRqpBFxLvCnwMWZeWK0JUmSNDnWbaQR\ncQZwJ/D86MuRJE26UU02qqsqN2T4CMsnZ7864lokSZo4azbSiLgKeCYzH+yumhp5RZKkidZud4q/\n6my9od0DQCci9gI/xfKU4X+VmU+v9oGI2APsKVahJKm4jd7cQKtbs5Fm5htXfo6Ih4DZtZpo9zPz\nwHzvutnZ2VsGL1GSVFqhmxucVqc9qj3Xk5e/SJKK8vKXVWTmm0ZZiCRJk8hEKkkqqu6Tg0rzeaSS\nJA3BRCpJKmo0N2So79WXJlJJkoZgIpUkFdW0RGojlSQV1W7Y5S8O7UqSNAQTqSSpKJ/+IkmSKjOR\nSpKKMpFKkqTKTKSSpKKadotAG6kkqaimPf3FoV1JkoZgIpUkFdW0B3ubSCVJGoKJVJJUVNMmG5lI\nJUkagolUklRU027IYCOVJBXVtEbq0K4kSUMwkUqSivJ5pJIkqTITqSSpqKadI7WRSpKKalojdWhX\nkqQhmEglSUV5ZyNJklSZiVSSVFSdn0caEWcC9wDnACeBKzPz2dO8rwV8Fvj9zLxzrX2aSCVJTXIQ\neDwzLwLuBm5a5X23An8fWPevAhupJKmoTrtT/FXQbuBY9+djwN7+N0TE24Cl7vap9Xbo0K4kqai6\nTDaKiKuB6/tWPw0sdn8+Cezo+8xrgcuBtwG3VDnOZjXSL3zxS3/xxk06ltQg28ZdQG3F8SPjLqHO\nvjDuAgYREXN9q+Yzc36192fmEeD/+z9CRNwLzHQXZ4Dn+j72LuAfAZ8HXg18KyL+V2Y+uNpxNqWR\n7t27d89mHKeqiJjLzLlx11FH/m7W5u9ndf5u1tak30+nvTSS/Rb6/S0AlwKPAZcAD/cd4/0rP0fE\nLcBX12qi4NCuJKlZ7gDuiohHgFPAFQARcQPwZGY+sNEd2kglSUWNKpGWkJkvAG8/zfqPnWbdB6vs\n01m7kiQNoamJdH7cBdTY/LgLqLn5cRdQY/PjLqDm5sddwGapcyIdhXWvj5Ekqarjx493rjr6V8X3\ne/SqV7F3795a9iyHdiVJGkJTh3YlSSPStKFdE6kkSUMwkUqSijKRSpKkykykkqSimpZIG9VIuw9q\nvR04j+VbQ12TmU+Nt6p6iYg3AP8hM9807lrqIiLOAD4F/BDwCuDWQW4jtlVFxDRwGPgxlp/d+J7M\n/MvxVlUvEXEu8KfAxZl5Ytz1jFrTGmnThnYvA7Zl5i7gRuDQmOuplYh4H8v/IL5i3LXUzDuBZ7oP\nAv454JNjrqdu3gy0M/NClh+S/KtjrqdWun+I3Qk8P+5aNBpNa6TfeaBrZj4K7BxvObXzJPAWvFFH\nvwRu7v7cAl4aYy21k5mfAWa7i68Gvj6+amrpIyzfKP2r4y5ks3TaS8Vfdda0Rrqd7z7QFWCpO9wr\nIDPvwybxMpn5fGb+XUTMsNxUf3ncNdVNZi5FxFHgE8Bvj7mc2oiIq1gezVh5DJd/pG5BTWsii3z3\nga4Arcxsj6sYTY6IeBXLD/q9OzN/Z9z11FFmXsXyedLDEXHmmMupiwPAvoh4CPgplh/f9f1jrmnk\n2u2l4q86a9RkI5Yf6LofyIg4H3hizPVoAnT/4XsQuDYzHxp3PXUTEe8CfjAzPwy8ALS7r8bLzDeu\n/NxtprOZ+fQYS9IINK2R3s/yX4cL3eUD4yymxjrjLqBmPgDsAG6OiJVzpZdk5otjrKlOfg84GhFf\nAM4ArsvMU2OuSWNU93OapTleL0kq5vjx4513fPyx4vv93ev/uU9/kSRpK2ra0K4kacQ6S80a2jWR\nSpI0BBOpJKmopk02MpFKkjQEE6kkqaimJVIbqSSpqKY1Uod2JUkagolUklRUp92sO0SaSCVJGoKJ\nVJJUlOdIJUlSZSZSSVJRTUukNlJJUlF1fxB3aQ7tSpI0BBOpJKkon/4iSZIqM5FKkopq2mQjE6kk\nSUMwkUqSimpaIrWRSpKKqnMjjYgzgXuAc4CTwJWZ+Wzfe24Ddne3d4DLMnNxtX3aSCVJTXIQeDwz\nPxQR7wBuAq7ve8/rgZ/JzL+tskPPkUqSiuq0l4q/CtoNHOv+fAzY27sxIlrAjwKHI+KLEXFgvR2a\nSCVJRf3hr7+z+D7b7fbJiJjrWz2fmfOrfSYirublafNpYGWY9iSwo2/7K4FPAB9luUc+FBF/kpl/\nvtpxbKSSpGL27t07Ne4aVmTmEeBI77qIuBeY6S7OAM/1feybwCcy88Xu+z8PvA6wkUqSBCwAlwKP\nAZcAD/dt/3Hg0xHxemAauBA4utYObaSSpCa5A7grIh4BTgFXAETEDcCTmflARNwN/BHwbeBoZn55\nbNVKkiRJkiRJkiRJkiRJkiRJkiRJkiRpbP4fxdAX36HwQ58AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f21d552a890>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D = np.eye(5)*np.sqrt ( post_covariance.diagonal() )\n",
    "R = np.linalg.inv ( D )\n",
    "correlation_matrix =  R.dot(post_covariance.dot(R))\n",
    "plt.figure(figsize=(8,8))\n",
    "plt.imshow ( correlation_matrix,interpolation='nearest', vmin=-0.5, vmax=0.5, cmap=plt.cm.RdBu_r )\n",
    "plt.colorbar()\n",
    "print correlation_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
 }
	{
	"cells": [
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"# Calculating the uncertainty in a calibration scenario\n",
	"\n",
	"After you minimize a cost function $J(\\vec{x})$ (where $\\vec{x}$ is the vector of your parameters), you end up having what is called the maximum a posteriori (MAP) estimate of the parameters, $\\vec{x}_{MAP}$. If the model is linear and the statistics are normal (Gaussian), then the a posteriori distribution is effectively Gaussian, with a mean given by $\\vec{x}_{MAP}$ and an uncertainty given by the inverse of the Hessian calculated at $\\vec{x}_{MAP}$, $\\mathbf{C}_{posterior}$. \n",
	"\n",
	"In the case of DALEC, the model isn't linear, but the statistics are Gaussian, so we pay a price for this convenience."
	]
	},
	{
	"cell_type": "code",
	"execution_count": 1,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"import numpy as np\n",
	"import scipy.optimize as opt\n",
	"import matplotlib.pyplot as plt\n",
	"\n",
	"%matplotlib inline\n",
	"\n",
	"def calculate_hessian ( func, x, epsilon=1e-8, args=()):\n",
	" \"\"\"This function calculates the Hessian of ``func`` at point ``x`` using \n",
	" finite differences. You might need to change the step size to different\n",
	" valuews (from 1e-3 to 1e-10) if the solution is unstable and you get NaNs\n",
	" \"\"\"\n",
	" # First, we calculate the Jacobian using finite differences:\n",
	" jac = opt.approx_fprime ( x, func, epsilon, *args )\n",
	" npar = len( x )\n",
	" hessian = np.zeros (( npar, npar))\n",
	" xp = 1.*x\n",
	" for i in xrange(npar):\n",
	" xx0 = 1.*x[i]\n",
	" xp[i] = xx0 + epsilon\n",
	" jac1 = opt.approx_fprime ( xp, func, epsilon, *args )\n",
	" hessian[i, :] = ( jac1 - jac)/epsilon\n",
	" xp[i] = xx0\n",
	" return hessian\n",
	"\n",
	"\n",
	"\n",
	" "
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"\n",
	"We now try an example function, the [Rosenbrock function](https://en.wikipedia.org/wiki/Rosenbrock_function), a function commonly used in optimisation problems. There is an implementation in ``scipy.optimize``, and it is useful because we have an analytical version of the Hessian to compare against (given by [scipy.optimize.rosen_hess](http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.optimize.rosen_hess.html). We will calculate the Hessian for the five dimensional version of the Rosenbrock function at $[0.9, 1.0, 0.2, 0.1, -1]^\\top$, and compare the result with the analytic version:"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 2,
	"metadata": {
	"collapsed": false
	},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"x_map point is [ 0.9 1. 0.2 0.1 -1. ]\n",
	"Numerical approximation to the Hessian:\n",
	"[[ 574.0216 -360.0019 -0.0003 0.0003 0. ]\n",
	" [ -360.0019 1322.0239 -400.002 0. 0. ]\n",
	" [ -0.0003 -400.002 210.0043 -80.0023 0. ]\n",
	" [ 0.0003 0. -80.0023 614.0019 -40.0019]\n",
	" [ 0. 0. 0. -40.0019 200.0004]]\n",
	"Analytic Hessian is:\n",
	"[[ 574. -360. 0. 0. 0.]\n",
	" [ -360. 1322. -400. 0. 0.]\n",
	" [ 0. -400. 210. -80. 0.]\n",
	" [ 0. 0. -80. 614. -40.]\n",
	" [ 0. 0. 0. -40. 200.]]\n",
	"Are they similar? If so the norm of the difference vs the true norm should be very small: 1.86707213248e-05\n"
	]
	}
	],
	"source": [
	"# Change the way we print numbers to the screen so it's easier to see\n",
	"np.set_printoptions (linewidth=132, precision=4, suppress=True )\n",
	"# Set the point around which we'll calculate the Hessian\n",
	"x_map = np.array([0.9, 1.0, 0.2, 0.1, -1.])\n",
	"print \"x_map point is\", x_map\n",
	"print \"Numerical approximation to the Hessian:\"\n",
	"print calculate_hessian ( opt.rosen, x_map, epsilon=1e-5)\n",
	"print \"Analytic Hessian is:\"\n",
	"print opt.rosen_hess ( x_map )\n",
	"print \"Are they similar? If so the norm of the difference vs the true norm should be very small:\", \\\n",
	" np.linalg.norm ( calculate_hessian ( opt.rosen, x_map, epsilon=1e-5)- \\\n",
	" opt.rosen_hess ( x_map ))/np.linalg.norm(opt.rosen_hess(x_map))"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"So the code above shows that the Hessian works. In order to get the parameter uncertainties, you would \n",
	"\n",
	"1. Invert the hessian to obtain $\\mathbf{C}_{post}$,\n",
	"2. Extract the main diagonal elements from $\\mathbf{C}_{post}$, and take the square root to find the standard deviation for each parameters\n",
	"3. Maybe plot $\\mathbf{C}_{post}$, to see the correlations between parameters etc.\n",
	"\n",
	"NOTE Sometimes, the calculation of the Hessian is unstable, and you end up with a Hessian with stupid values such as ``NaN`` or somesuch. Some other times, the numerical errors result in a matrix that cannot be inverted as we do below with ``np.linalg.inv``. If that happens, try other values of ``epsilon``, as mentioned above."
	]
	},
	{
	"cell_type": "code",
	"execution_count": 3,
	"metadata": {
	"collapsed": false
	},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"STandard deviations: [ 0.0555 0.0583 0.1367 0.0445 0.0713]\n"
	]
	}
	],
	"source": [
	"hess = calculate_hessian ( opt.rosen, x_map, epsilon=1e-5)\n",
	"post_covariance = np.linalg.inv ( hess )\n",
	"sigma = np.sqrt ( post_covariance.diagonal())\n",
	"print \"STandard deviations: \", sigma "
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"It's useful to also plot the correlation matrix between the different parameters. Parameters that are strongly correlated (either positively or negatively) are parameters that compensate for each other and an indication that their values cannot be obtained from the calibration directly, but must be helped out with priors, etc."
	]
	},
	{
	"cell_type": "code",
	"execution_count": 4,
	"metadata": {
	"collapsed": false
	},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"[[ 1. 0.6591 0.5638 0.2289 0.0286]\n",
	" [ 0.6591 1. 0.8555 0.3474 0.0433]\n",
	" [ 0.5638 0.8555 1. 0.4061 0.0507]\n",
	" [ 0.2289 0.3474 0.4061 1. 0.1248]\n",
	" [ 0.0286 0.0433 0.0507 0.1248 1. ]]\n"
	]
	},
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAdIAAAHVCAYAAABFZXqOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG99JREFUeJzt3X+M5PV93/HnzoZzMexdUwXSqLHiRvnhVA5O3WsMdwif\ny10SqE8ltt+2wLHhAur6aC2gf9jEIbB2iVzFOds4DohczzoojRO9A8Qidq70atbgbYJInEIS2T2B\nWjWKHAqO8V4wnM3O9I+dtafD7e53Zz6z8539Ph/SyPv9fme+3/cs1r339fl+vt8vSJIkSZIkSZIk\nSZIkSZIkSZIkSbU1Ne4CJEkahYhoAbcD5wGngGsy86me7TcAVwPPdFfNZuaJjR7newrUKklSHV0G\nbMvMXRHxBuBQd92K1wPvysw/G+YgrWE+LElSje0GjgFk5qPAzr7t/wz4QEQ8EhE3DnoQG6kkaava\nDiz2LC91h3tXfBqYBf4FcGFE/MtBDuLQriSpmM/d95nOtu1nFd9vu90+efjw4Y/2rZ7PzPk1PrYI\nzPQstzKz3bN8W2YuAkTEZ4F/Cnx2o7XZSCVJxWzbfha59+ri+43jR2Yyc26DH1sA9gMZEecDT3xn\nfxE7gCci4p8A32Q5lR4ZpDYbqSSpqOn6XA9yP7AvIha6ywci4nLg7Mw83D0v+hDLM3qPZ+axQQ5i\nI5UkbUmZ2QEO9q0+0bP90yyfJx2KjVSSVNT0VH0i6WZw1q4kSUMwkUqSiqrROdJNYSOVJBXl0K4k\nSarMRCpJKqppQ7smUkmShmAilSQV5TlSSZJUmYlUklRU086R2kglSUU5tCtJkiozkUqSimpaQmva\n95UkqSgTqSSpKM+RSpKkykykkqSivPxFkqQhOLQrSZIqM5FKkopq2tCuiVSSpCGYSCVJRXmOVJIk\nVWYilSQV1bRzpDZSSVJRDu1KkqTKTKSSpKKaNrRrIpUkaQgmUklSUSZSSZJUmYlUklRU02bt2kgl\nSUU5tCtJkiozkUqSimra0K6JVJKkIZhIJUlFeY60ASJiz7hrqCt/N2vz97M6fzdr8/ezdTU1ke4B\n5sdcQ13twd/NWvbg72c1e/B3s5Y9NOT3U5dzpBHRAm4HzgNOAddk5lOned9vAV/LzF8a5DiNTKSS\npNGZnir/GtBlwLbM3AXcCBzqf0NEzAKvBTqDHsRGKknaqnYDxwAy81FgZ+/GiNgF/DRwJzBwu7aR\nSpKKmp6aKv4a0HZgsWd5qTvcS0T8AHAz8G8Zooky7Ier+u1Dn+yc+7rXbMahJEnr2H7WK/npC3aN\n5N//48ePd57/heuK7/ese27jzjvv/GDf6vnMnF/tMxFxCPjjzMzu8l9l5qu6P78XuBI4CfxD4JXA\nr2Tm3RutbVMmG537uteQe6/ejENJjXLrx9867hJqa+bd7xt3CbX1xS/9xUj33xrRZKPMnNvgRxaA\n/UBGxPnAEz37+g3gNwAi4krgNYM0UWjurF1J0ohM1edC0vuBfRGx0F0+EBGXA2dn5uG+9w482chG\nKknakjKzAxzsW33iNO+7a5jj2EglSUW16pNIN4WzdiVJGoKJVJJU1NR0szJas76tJEmFmUglSUXV\naNbuprCRSpKKcrKRJEmqzEQqSSpqqtWsjNasbytJUmEmUklSUaM4Rzrw/fs2gYlUkqQhmEglSUWN\n4vKXOidSG6kkqSjvbCRJkiozkUqSihrFZKOl4nssx0QqSdIQTKSSpKKmWt4iUJIkVWQilSQV1WrY\nrF0bqSSpqKY9Rq1ZfzZIklSYiVSSVJSJVJIkVWYilSQV1bTJRs36tpIkFWYilSQV1bRzpDZSSVJR\nLe9sJEmSqlo3kUZEC7gdOA84BVyTmU+NujBJ0mTyeaQvdxmwLTN3ATcCh0ZbkiRJk6NKI90NHAPI\nzEeBnSOtSJI00VrTU8VfdVZlstF2YLFneSkiWpnZPt2bI2IPsGf40iRJoxIRc32r5jNzfgylTLwq\njXQRmOlZXrWJAnT/Q8z3rpudnb1lkOIkSaORmXOj2reXv7zcArAfyIg4H3hitCVJkiZZ0yYbVWmk\n9wP7ImKhu3xghPVIkjRR1m2kmdkBDm5CLZKkLaDuk4NK885GkqQtab37IETEW4H3Ax3gP2fmJwY5\nTrMGsiVJIzfVmir+GtCq90GIiGngw8DFwAXAtRHxDwY5iI1UkrRVrXofhMxcAl6TmSeBc4Bp4FuD\nHMRGKkkqqjXdKv4a0Gnvg7CykJntiHgL8GfAQ8A3BzmI50glSUWN6jrSAW4ise59EDLzvoi4HzgK\nvLv7vxtiI5UkTYQBbiKx6n0QImI78ACwLzO/FRHPA0uD1GUjlSQVVaMbMrzsPggRcTlwdmYejoh7\ngIcj4tvA48A9gxzERipJ2pJWuQ/CiZ7th4HDwx7HRipJKmqqVZtEuima9W0lSSrMRCpJKmqIy1Um\nko1UklRUjSYbbYpmfVtJkgozkUqSijKRSpKkykykkqSivPxFkiRVZiKVJBU1NT097hI2lY1UklSU\nk40kSVJlJlJJUlEtJxtJkqSqTKSSpKI8RypJkiozkUqSimpaIrWRSpKK8s5GkiSpMhOpJKmopg3t\nNuvbSpJUmIlUklRU0xKpjVSSVFTLRippUpy477Fxl1BbO684Ne4S1BA2UklSUV7+IkmSKjORSpKK\natpko2Z9W0mSCjORSpKKaloitZFKkopyspEkSarMRCpJKqo1PT3uEjaViVSSpCGYSCVJRTVtslGz\nvq0kSYWZSCVJRTUtkdpIJUlFNe3yFxupJGlLiogWcDtwHnAKuCYzn+rZfjlwHfAS8OfAtZnZ2ehx\nmvVngyRp5KamW8VfA7oM2JaZu4AbgUMrGyLiTODfA3sy80JgB/DmQQ5iI5UkbVW7gWMAmfkosLNn\n24vABZn5Ynf5e4AXBjmIQ7uSpKJqNNloO7DYs7wUEa3MbHeHcJ8BiIj3Amdl5vFBDmIjlSRNhIiY\n61s1n5nza3xkEZjpWW5lZrtnfy3g14AfAd46aF02UklSUaOatZuZcxv8yAKwH8iIOB94om/7nSwP\n8f78IJOMVthIJUlFTbVqc6/d+4F9EbHQXT7Qnal7NvAnwC8CDwOfjwiA2zLz9zd6EBupJGlL6qbM\ng32rT/T8XKTj20glSWXVJ5FuitpMrZIkaRKZSCVJZTXsFoHN+raSJBVmIpUkFTU13axzpDZSSVJZ\nTjaSJElVmUglSWWZSCVJUlUmUklSUaO6125dbejbRsQbIuKhURUjSdKkqZxII+J9wC8Afze6ciRJ\nE69h50g3MrT7JPAW4D+NqBZJ0lbQsEZaeWg3M+8DXhphLZIkTZzik40iYg+wp/R+JUnlRMRc36r5\nzJwvse+mTTYq3ki7/yHme9fNzs7eUvo4kqTBZebcuGvYKgZppJ3iVUiSto5RnCOtcefZUCPNzP8N\n7BpNKZIkTR5vyCBJKmsUiXSp/C5LsZFKkooayWPUatxImzW1SpKkwkykkqSyGnb5S7O+rSRJhZlI\nJUlleYtASZJUlYlUklTUVMMSqY1UklSWk40kSVJVJlJJUlFNG9o1kUqSNAQTqSSpLBOpJEmqykQq\nSSqrYbN2baSSpKJG8vSXGmvWnw2SJBVmIpUklTWSyUb1fSCpiVSSpCGYSCVJZdUkkUZEC7gdOA84\nBVyTmU/1veeVwH8FfjEz/+cglZlIJUlb1WXAtszcBdwIHOrdGBE7gYeBfwx0Bj2IjVSSVNRUq1X8\nNaDdwDGAzHwU2Nm3fRvLzXagJLrCoV1JUlkjurNRRMz1rZrPzPk1PrIdWOxZXoqIVma2ATLzv3f3\nO1RdNlJJ0kTIzLkNfmQRmOlZ/k4TLcmhXUlSWVOt8q/BLACXAkTE+cATpb5iLxOpJGmruh/YFxEL\n3eUDEXE5cHZmHi51EBupJKmswRNkUZnZAQ72rT5xmve9aZjj2EglSUV1atJIN0uzvq0kSYWZSCVJ\nZZlIJUlSVSZSaYL9zVe+Nu4S6qtT/HJBVTU1Ne4KNpWJVJKkIZhIJUllDX5v3IlkI5UkFeXlL5Ik\nqTITqSSpLBOpJEmqykQqSSrLRCpJkqoykUqSympYIrWRSpKK8vIXSZJUmYlUklSWiVSSJFVlIpUk\nleXTXyRJUlUmUklSWQ07R2ojlSQV5eUvkiSpMhOpJKmshj3Yu1nfVpKkwkykkqSyPEcqSZKqMpFK\nkspqWCK1kUqSympYI23Wt5UkqTATqSSpKG/IIEmSKjORSpLKMpFKkqSqTKSSpLIa9jxSG6kkqayG\nDe2u2Ugj4gzgU8APAa8Abs3MBzajMEmShhERLeB24DzgFHBNZj7Vs30/8CvAS8CnMvM/DnKc9f5s\neCfwTGZeBPwc8MlBDiJJao7OVKv4a0CXAdsycxdwI3BoZUM3KH4U2Ae8EfjXEXHuIAdZr7oEbu55\n70uDHESSpDHYDRwDyMxHgZ09234CeDIzv5GZ3wa+CFw0yEHWHNrNzOcBImKG5ab6y4McRJLUIPU5\nR7odWOxZXoqIVma2u9u+0bPtJLBjkIOsO9koIl4F3Af8Zmb+ToX37wH2DFKMJGlzRMRc36r5zJwf\nQymVDVDzIjDTs7zSRGG5ifZumwG+Pkhd6002+n7gQeDazHyoyg67X2q+d93s7OwtgxQnSRqNzJwb\n1b47I7r8ZYCaF4D9QEbE+cATPdu+AvxoRHwv8DzLw7ofGaSu9RLpB1iOujdHxMq50ksy88VBDiZJ\n2vo6nXFX8B33A/siYqG7fCAiLgfOzszDEfHvgP/C8hygI5n51UEOst450uuA6wbZsSRJ45SZHeBg\n3+oTPdv/APiDYY/jDRkkSUW1axRJN0NtplZJkjSJTKSSpKKalUdNpJIkDcVEKkkqqt2wSGojlSQV\n1XGykSRJqspEKkkqqmlDuyZSSZKGYCKVJBXVsEBqIpUkaRgmUklSUU07R2ojlSQV5eUvkiSpMhOp\nJKmo9rgL2GQmUkmShmAilSQV1bBTpCZSSZKGYSKVJBXl5S+SJA3By18kSVJlJlJJUlFe/iJJkioz\nkUqSimrYKVIbqSSprHbDOqlDu5IkDcFEKkkqqll51EQqSdJQTKSSpKKadmcjE6kkSUMwkUqSimrY\npF0bqSSprHbDphttWiO99eNv3axDTZQT9z027hJq7W++8rVxl1BrD/7f58ddQm3t/+svj7uEGts2\n7gK2FBOpJKmopg3tOtlIkqQhmEglSUV5+YskSarMRCpJKqrO50gj4kzgHuAc4CRwZWY+e5r3nQMs\nAK/NzG+ttU8TqSSpqDad4q+CDgKPZ+ZFwN3ATf1viIifBR4Ezq2yQxupJKlJdgPHuj8fA/ae5j1L\nwMXA16vs0KFdSVJRdRnajYirgev7Vj8NLHZ/Pgns6P9cZh7vfr7ScWykkqSJEBFzfavmM3N+tfdn\n5hHgSN8+7gVmuoszwHPD1mUjlSQV1R5RJM3MuQK7WQAuBR4DLgEeHnaHNlJJUpPcAdwVEY8Ap4Ar\nACLiBuDJzHyg572V/iKwkUqSilpqj7uC1WXmC8DbT7P+Y6dZ98NV9mkjlSQVNaqh3bry8hdJkoZg\nIpUkFbVkIpUkSVWZSCVJRXmOVJIkVWYilSQVVefLX0bBRipJKsqhXUmSVJmJVJJUlJe/SJKkykyk\nkqSi2iMIpFPld1mMiVSSpCGYSCVJRS2NIJLWuVnVuTZJ0gTy8hdJklSZiVSSVNRSswKpiVSSpGGY\nSCVJRTXtHOm6jTQipoHDwI8BHeA9mfmXoy5MkqRJUGVo981AOzMvBG4CfnW0JUmSJtlSu1P8VWfr\nNtLM/Aww2118NfD1URYkSZps7U6n+KvOKp0jzcyliDgK/DzwtpFWJEnSBKk82Sgzr4qI9wOPRsRP\nZOYLp3tfROwB9pQpT5I0ChEx17dqPjPnS+y7aZe/VJls9C7gBzPzw8ALQLv7Oq3uf4j53nWzs7O3\nDFWlJKmozJwbdw1bRZVE+nvA0Yj4AnAGcF1mnhptWZKkSVX3c5qlrdtIu0O479iEWiRJmjjekEGS\nVFS75perlGYjlSQV1bTJRt5rV5KkIZhIJUlFNW2ykYlUkqQhmEglSUUtmUglSVJVJlJJUlFe/iJJ\n0hC8/EWSJFVmIpUkFVXny18i4kzgHuAc4CRwZWY+2/eeG/jurXE/l5kfWmufJlJJUpMcBB7PzIuA\nu4GbejdGxA8DVwAXZOb5wM9ExE+utUMbqSSpqKVOp/iroN3Ase7Px4C9fdv/D/Czmbly0DNYfoTo\nqhzalSRNhI0+jDwirgau71v9NLDY/fkksKN3Y2a+BPxtREwBHwG+lJlPrlWXjVSSVNTSiC5/2ejD\nyDPzCHCkd11E3AvMdBdngOf6PxcRfw/4FPAN4Nr1jmMjlSQVNapGWsgCcCnwGHAJ8HDvxm4S/Qzw\n3zLz16rs0EYqSWqSO4C7IuIR4BTLE4tWZuo+CUwDFwFnRMQl3c/8Umb+8Wo7tJFKkoqqcyLNzBeA\nt59m/cd6Fs/cyD6dtStJ0hBMpJKkouqcSEfBRCpJ0hBMpJKkopqWSG2kkqSiRtJIazx+WuPSJEmq\nPxOpJKkoE6kkSarMRCpJKsrJRpIkDaFpjdShXUmShmAilSQVZSKVJEmVmUglSUW91LBEummNdObd\n79usQ02UnVecGncJ9dZpj7uCWtv/118edwm19d7XHxx3CbUVx4+Mu4QtxUQqSSqqaedIbaSSpKKa\n1kidbCRJ0hBMpJKkopY6JlJJklSRiVSSVJTnSCVJUmUmUklSUU1LpDZSSVJRTWukDu1KkjQEE6kk\nqaildrNu7WkilSRpCCZSSVJRniOVJEmVmUglSUU1LZHaSCVJRTXtwd4O7UqSNAQTqSSpqKYN7ZpI\nJUkagolUklSUiVSSJFVmIpUkFVXnRBoRZwL3AOcAJ4ErM/PZvvf8G+BKoAP8embmWvs0kUqSilpq\nd4q/CjoIPJ6ZFwF3Azf1boyI7wPeA1wAXAwcWm+HNlJJUpPsBo51fz4G7O3d2E2nr8vMJeAHgBfX\n26FDu5KkouoytBsRVwPX961+Gljs/nwS2NH/ucxsd4d3Pwjctt5xbKSSpIkQEXN9q+Yzc36192fm\nEeBI3z7uBWa6izPAc6t89jcj4reAP4yIR9Y6jo1UklRUZ0SJNDPnCuxmAbgUeAy4BHi4d2NE/Djw\n4cx8C/AScApYWmuHNlJJUpPcAdwVEY+w3CSvAIiIG4AnM/OBiPgfEfFHLM/a/VxmPrLWDm2kkqSi\n2jU5R3o6mfkC8PbTrP9Yz88fAj5UdZ82UklSUZ1OfRvpKFRqpBFxLvCnwMWZeWK0JUmSNDnWbaQR\ncQZwJ/D86MuRJE26UU02qqsqN2T4CMsnZ7864lokSZo4azbSiLgKeCYzH+yumhp5RZKkidZud4q/\n6my9od0DQCci9gI/xfKU4X+VmU+v9oGI2APsKVahJKm4jd7cQKtbs5Fm5htXfo6Ih4DZtZpo9zPz\nwHzvutnZ2VsGL1GSVFqhmxucVqc9qj3Xk5e/SJKK8vKXVWTmm0ZZiCRJk8hEKkkqqu6Tg0rzeaSS\nJA3BRCpJKmo0N2So79WXJlJJkoZgIpUkFdW0RGojlSQV1W7Y5S8O7UqSNAQTqSSpKJ/+IkmSKjOR\nSpKKMpFKkqTKTKSSpKKadotAG6kkqaimPf3FoV1JkoZgIpUkFdW0B3ubSCVJGoKJVJJUVNMmG5lI\nJUkagolUklRU027IYCOVJBXVtEbq0K4kSUMwkUqSivJ5pJIkqTITqSSpqKadI7WRSpKKalojdWhX\nkqQhmEglSUV5ZyNJklSZiVSSVFSdn0caEWcC9wDnACeBKzPz2dO8rwV8Fvj9zLxzrX2aSCVJTXIQ\neDwzLwLuBm5a5X23An8fWPevAhupJKmoTrtT/FXQbuBY9+djwN7+N0TE24Cl7vap9Xbo0K4kqai6\nTDaKiKuB6/tWPw0sdn8+Cezo+8xrgcuBtwG3VDnOZjXSL3zxS3/xxk06ltQg28ZdQG3F8SPjLqHO\nvjDuAgYREXN9q+Yzc36192fmEeD/+z9CRNwLzHQXZ4Dn+j72LuAfAZ8HXg18KyL+V2Y+uNpxNqWR\n7t27d89mHKeqiJjLzLlx11FH/m7W5u9ndf5u1tak30+nvTSS/Rb6/S0AlwKPAZcAD/cd4/0rP0fE\nLcBX12qi4NCuJKlZ7gDuiohHgFPAFQARcQPwZGY+sNEd2kglSUWNKpGWkJkvAG8/zfqPnWbdB6vs\n01m7kiQNoamJdH7cBdTY/LgLqLn5cRdQY/PjLqDm5sddwGapcyIdhXWvj5Ekqarjx493rjr6V8X3\ne/SqV7F3795a9iyHdiVJGkJTh3YlSSPStKFdE6kkSUMwkUqSijKRSpKkykykkqSimpZIG9VIuw9q\nvR04j+VbQ12TmU+Nt6p6iYg3AP8hM9807lrqIiLOAD4F/BDwCuDWQW4jtlVFxDRwGPgxlp/d+J7M\n/MvxVlUvEXEu8KfAxZl5Ytz1jFrTGmnThnYvA7Zl5i7gRuDQmOuplYh4H8v/IL5i3LXUzDuBZ7oP\nAv454JNjrqdu3gy0M/NClh+S/KtjrqdWun+I3Qk8P+5aNBpNa6TfeaBrZj4K7BxvObXzJPAWvFFH\nvwRu7v7cAl4aYy21k5mfAWa7i68Gvj6+amrpIyzfKP2r4y5ks3TaS8Vfdda0Rrqd7z7QFWCpO9wr\nIDPvwybxMpn5fGb+XUTMsNxUf3ncNdVNZi5FxFHgE8Bvj7mc2oiIq1gezVh5DJd/pG5BTWsii3z3\nga4Arcxsj6sYTY6IeBXLD/q9OzN/Z9z11FFmXsXyedLDEXHmmMupiwPAvoh4CPgplh/f9f1jrmnk\n2u2l4q86a9RkI5Yf6LofyIg4H3hizPVoAnT/4XsQuDYzHxp3PXUTEe8CfjAzPwy8ALS7r8bLzDeu\n/NxtprOZ+fQYS9IINK2R3s/yX4cL3eUD4yymxjrjLqBmPgDsAG6OiJVzpZdk5otjrKlOfg84GhFf\nAM4ArsvMU2OuSWNU93OapTleL0kq5vjx4513fPyx4vv93ev/uU9/kSRpK2ra0K4kacQ6S80a2jWR\nSpI0BBOpJKmopk02MpFKkjQEE6kkqaimJVIbqSSpqKY1Uod2JUkagolUklRUp92sO0SaSCVJGoKJ\nVJJUlOdIJUlSZSZSSVJRTUukNlJJUlF1fxB3aQ7tSpI0BBOpJKkon/4iSZIqM5FKkopq2mQjE6kk\nSUMwkUqSimpaIrWRSpKKqnMjjYgzgXuAc4CTwJWZ+Wzfe24Ddne3d4DLMnNxtX3aSCVJTXIQeDwz\nPxQR7wBuAq7ve8/rgZ/JzL+tskPPkUqSiuq0l4q/CtoNHOv+fAzY27sxIlrAjwKHI+KLEXFgvR2a\nSCVJRf3hr7+z+D7b7fbJiJjrWz2fmfOrfSYirublafNpYGWY9iSwo2/7K4FPAB9luUc+FBF/kpl/\nvtpxbKSSpGL27t07Ne4aVmTmEeBI77qIuBeY6S7OAM/1feybwCcy88Xu+z8PvA6wkUqSBCwAlwKP\nAZcAD/dt/3Hg0xHxemAauBA4utYObaSSpCa5A7grIh4BTgFXAETEDcCTmflARNwN/BHwbeBoZn55\nbNVKkiRJkiRJkiRJkiRJkiRJkiRJkiRpbP4fxdAX36HwQ58AAAAASUVORK5CYII=\n",
	"text/plain": [
	"<matplotlib.figure.Figure at 0x7f21d552a890>"
	]
	},
	"metadata": {},
	"output_type": "display_data"
	}
	],
	"source": [
	"D = np.eye(5)*np.sqrt ( post_covariance.diagonal() )\n",
	"R = np.linalg.inv ( D )\n",
	"correlation_matrix = R.dot(post_covariance.dot(R))\n",
	"plt.figure(figsize=(8,8))\n",
	"plt.imshow ( correlation_matrix,interpolation='nearest', vmin=-0.5, vmax=0.5, cmap=plt.cm.RdBu_r )\n",
	"plt.colorbar()\n",
	"print correlation_matrix"
	]
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": []
	}
	],
	"metadata": {
	"kernelspec": {
	"display_name": "Python 2",
	"language": "python",
	"name": "python2"
	},
	"language_info": {
	"codemirror_mode": {
	"name": "ipython",
	"version": 2
	},
	"file_extension": ".py",
	"mimetype": "text/x-python",
	"name": "python",
	"nbconvert_exporter": "python",
	"pygments_lexer": "ipython2",
	"version": "2.7.8"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 0
	}
No results found