noisy_kriging_skiptest.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:2cffd02a292fe739c5419f40f87bd8396494439dc63700e54d54a42dae8bb164"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Noisy kriging observations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Y process is conditioned on heterogenous noisy observations:\n",
      "\n",
      "$$\\hat{y}_i = y(x^{(i)})+\\epsilon_i, \\epsilon_i \\sim \\mathcal{N}(0, \\tau_i^2)$$\n",
      "\n",
      "$$C^{-1} \\leftarrow (C+\\Delta)^{-1}, \\Delta=I_N diag(\\tau_1^2, \\dots, \\tau_N^2)$$\n",
      "\n",
      "<!--- # $$\\hat{\\beta}=(F^\\top C^{-1} F)^{-1} F^\\top C^{-1} y$$\n",
      "\n",
      " #$$-2log L (\\beta, \\sigma^2, \\theta;y) = nlog(2\\pi)+log(|C|)+(y-F\\hat{\\beta})^\\top C^{-1} (y-F\\hat{\\beta})$$\n",
      " -->"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import print_function\n",
      "import openturns as ot\n",
      "from openturns.viewer import View\n",
      "import numpy as np\n",
      "import math as m\n",
      "import matplotlib.pyplot as pl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dimension = 1\n",
      "\n",
      "inputs = np.array([0, 2, 4, 6, 8, 9])\n",
      "outputs = np.array([0.25, 0.28, 0.32, 0.27, 0.22, 0.19])\n",
      "eps = 0.20\n",
      "cov = np.array([eps, 0.15, eps, 0.15, eps, 0.15])\n",
      "\n",
      "t = np.linspace(0, 10, 100)\n",
      "\n",
      "noise = (outputs * cov / 2)**2\n",
      "X = ot.Sample(inputs.reshape([-1,1]))\n",
      "Y = ot.Sample(outputs.reshape([-1,1]))\n",
      "\n",
      "# create algorithm\n",
      "basis = ot.ConstantBasisFactory(dimension).build()\n",
      "#covarianceModel = ExponentialModel([0.2], [0.1])\n",
      "covarianceModel = ot.SquaredExponential([0.2], [0.1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "algo = ot.KrigingAlgorithm(X, Y, covarianceModel, basis)\n",
      "algo.setNoise(noise)\n",
      "algo.run()\n",
      "\n",
      "# perform an evaluation\n",
      "result = algo.getResult()\n",
      "\n",
      "\n",
      "def suk_f(result):\n",
      "    def inner(X):\n",
      "        muk = result.getMetaModel()(X)[0]\n",
      "        sk2 = m.sqrt(result.getConditionalCovariance(X)[0,0])\n",
      "        Y = [muk+sk2, muk-sk2]\n",
      "        #print(X, sk2, Y)\n",
      "        return Y\n",
      "    return inner\n",
      "\n",
      "\n",
      "suk = ot.PythonFunction(1, 2, suk_f(result))\n",
      "\n",
      "fig = pl.figure()\n",
      "pl.errorbar(inputs, outputs, np.sqrt(noise), fmt='b.', markersize=10, label=u'Observations')\n",
      "g = ot.Graph()\n",
      "N = 100\n",
      "xmin = X.getMin()[0]-0.5\n",
      "xmax = X.getMax()[0]+0.5\n",
      "curve = result.getMetaModel().draw(xmin, xmax, N).getDrawable(0)\n",
      "curve.setColor('red')\n",
      "g.add(curve)\n",
      "curve = suk.getMarginal(0).draw(xmin, xmax, N).getDrawable(0)\n",
      "curve.setColor('green')\n",
      "g.add(curve)\n",
      "curve = suk.getMarginal(1).draw(xmin, xmax, N).getDrawable(0)\n",
      "curve.setColor('green')\n",
      "g.add(curve)\n",
      "View(g, figure=fig)\n",
      "\n",
      "\n",
      "print(\"X=\", X)\n",
      "print(\"f(X)=\", Y)\n",
      "print(\"sk2=\",suk(X))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "X= 0 : [ 0 ]\n",
        "1 : [ 2 ]\n",
        "2 : [ 4 ]\n",
        "3 : [ 6 ]\n",
        "4 : [ 8 ]\n",
        "5 : [ 9 ]\n",
        "f(X)= 0 : [ 0.25 ]\n",
        "1 : [ 0.28 ]\n",
        "2 : [ 0.32 ]\n",
        "3 : [ 0.27 ]\n",
        "4 : [ 0.22 ]\n",
        "5 : [ 0.19 ]\n",
        "sk2=     [ y0       y1       ]\n",
        "0 : [ 0.272988 0.232692 ]\n",
        "1 : [ 0.295325 0.2627   ]\n",
        "2 : [ 0.311605 0.275542 ]\n",
        "3 : [ 0.282149 0.250623 ]\n",
        "4 : [ 0.227762 0.203353 ]\n",
        "5 : [ 0.210868 0.185984 ]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcjXX/x/HXx66UnZswiUgkI0SFqYREQ9mikFuS/BLq\nlm4y1W1JSlpIC1JZQtbcJSRJjOzGNnQ3jSVLY19n+fz+uI5pMGbOmTlnzpk5n+d5nMec5Vq+lzHv\nc53v9V1EVTHGGBMccvm7AMYYY7KOhb4xxgQRC31jjAkiFvrGGBNELPSNMSaIWOgbY0wQcSv0RaS5\niOwQkV0iMjCN5R4VkSQRqZ3itUEiEi0i20WkqTcKbYwxJmPypLeAiOQC3gfuB/YDa0VknqruuGy5\nQsBzwOoUr1UD2gPVgHLAEhG5Wa1zgDHG+IU7Z/r1gGhVjVHVeGA6EJ7Kcq8DI4HzKV4LB6araoKq\n/g5Eu7ZnjDHGD9wJ/RuA2BTP97peSyYioUA5Vf1vOuvuu3xdY4wxWSfd6p30iIgAbwNdM18cY4wx\nvuRO6O8DKqR4Xs712kXXAdWB5a4PgH8A80XkYTfWBUBErI7fGGMyQFXF0xXSvAO5gd1ACJAP2AhU\nS2P5H4BQ1+NbgQ2u9Sq6tiOprKM52dChQ/1dBJ+y48vecvLx5eRjU1V1ZWe6OZ7ynu6Zvqomikgf\nYDHONYBPVXW7iLwKrFXVhZevAohr3W0i8hWwDYgHersKaowxxg/cqtNX1W+Bqpe9NvQqy9532fMR\nwIiMFtAYY4z3WI/cLBAWFubvIviUHV/2lpOPLycfW0ZJINS2iIjV+hhjjIdExOMLuXamb4wxQcRC\n3xhjgoiFvjHGBBELfWOMCSIW+sYYE0Qs9I0xJohY6BtjTBCx0DfGmCBioW+MMUHEQt8YY4KIhb4x\nxgQRC31jjAkiFvrGGBNELPSNMSaIWOgbY0wQsdA3xpggYqFvjDFBxELfGGOCiIW+McYEEQt9Y4wJ\nIhb6xhgTRCz0jTEmiFjoG2NMELHQN8aYIGKhb4wxQcRC3xhjgoiFvjHGBBELfWOMCSIW+sYYE0Qs\n9I0xJohY6BtjTBCx0DfGmCBioW+MMUHEQt8YY4KIW6EvIs1FZIeI7BKRgam8/7SIbBaRDSKyQkRu\ncb0eIiJnRGS96z7O2wdgjDHGfaKqaS8gkgvYBdwP7AfWAh1VdUeKZQqp6inX41ZAb1V9UERCgAWq\nWjOdfWh65TDGGHMpEUFVxZN18rixTD0gWlVjXDuZDoQDyaF/MfBdCgFJKcvlSYGMCQTLlzv3i4/D\nwpzHYWF/PzYmO3In9G8AYlM834vzQXAJEekN9AfyAveleOtGEVkHnACGqOrKjBfXmKyRMtxF/v4A\nMCa7cyf03aKq44BxItIRGAJ0Aw4AFVT1qIjUBuaKyK2XfTMAICIiIvlxWFgYYXY6ZYwxl1i+fDnL\nM3kG4k6dfn0gQlWbu56/BKiqvnGV5QU4qqpFUnnvB2CAqq6/7HWr0zcBSwTsv6cJRBmp03en9c5a\noLKrJU4+oCMw/7IdV07xtCXOhV9EpITrQjAichNQGfjNkwIaY4zxnnSrd1Q1UUT6AItxPiQ+VdXt\nIvIqsFZVFwJ9RKQJcAE4CnR1rd4IeE1ELuBc3H1aVY/54kCMMcakL93qnSwphFXvmABm1TsmUPmq\nescYY0wOYaFvjDFBxELfGGOCiIW+McYEEQt9Y4wJIhb6xhgTRCz0jTEmiFjoG2NMELHQN8aYIGKh\nb4wxQcRC3xhjgoiFvjHGBBELfWOMCSIW+sYYE0Qs9I0xJohY6BtjTBCx0DfGmCBioW+MMUHEQt8Y\nY4KIhb4xaYiLc34ePerfchjjLRb6xlzFmDFQu7bzODTUeW5Mdieq6u8yICIaCOUw5qK4OCfwY2L+\nfi0kBNavh2LF/FcuY1ISEVRVPFnHzvSNSUVUFMTGXvpabCxs2+af8hjjLRb6xqSiRg0oX/7S18qX\nh+rV/VMeY7zFQt+YVBQtCn37OlU64Pzs29d53ZjszOr0jUlDXBwUL+78tMA3gSYjdfoW+sakQwTs\nv6cJRHYh1xhjTJry+LsAxvjTec5zghOc5zwJJJBIInnJSz7yUZCCXMd12LmRyUks9E2OlUACu9nN\ndrazx3WLJZYDrlsccSSRxPVcTwEKkIc85CIXCSRwnvOc4QxnOQtHCnMLJSlHOcpTnpu4iVtct6pU\nJR/5/H2oxrjN6vRNjpBIIpvZTCSRrGUtv/Iru9hFWcpSjWpUpjKVqEQFKlDGdStOcQpQAOHqVaIJ\nJJC3xDGijhwi1nXbwx52spPtbCeGGG7lVupQh7u5m8Y0pgIVsvDITTCzC7kmaCSRxAY28D3f8yM/\nsopVlKUs9alPXepyB3dQgxpcy7WZ3ldaF3JPc5pNbGIta1nJSlawgmu4huY050Ee5H7u90oZjEmN\nhb7J0Y5ylG/5loUsZDGLKUEJHuAB7uVe7uEeSlLSJ/v1pPWOomxjG/913dayliY0oR3taElL1zUC\nY7zDQt/kOH/yJ3OZyyxmEUkkjWlMK1rxIA9SnvLpb8ALMtNkM4445jGPmcxkFat4mIfpSlfu5V5y\n2QVik0kW+iZHOMpRvuZrpjKVdayjBS1oS1ua05xruCbLy+OtdvqHOMRUpjKZyZzkJM/wDE/yJMUp\nnvmNm6Dks3b6ItJcRHaIyC4RGZjK+0+LyGYR2SAiK0TklhTvDRKRaBHZLiJNPSmcCR7xxLOABbSl\nLTdyI4tYRG96c4ADTGUqj/CIXwLfm0pRiud5ng1sYCpT2cIWKlOZnvRkJzv9XTwTJNI90xeRXMAu\n4H5gP7AW6KiqO1IsU0hVT7ketwJ6q+qDInIr8CVQFygHLAFuvvy03s70g9cudvEpn/IZn1GJSnSl\nK+1pTxGK+LtoyXzZI/cQhxjnut3FXQxmMHWo45udmRzHV2f69YBoVY1R1XhgOhCecoGLge9SCEhy\nPX4YmK6qCar6OxDt2p4JYhe4wAxmEEYYDWmIovzIj/zMz/SkZ0AFvq+VohQRRPA7v9OEJrShDQ/x\nEGtY4++imRzKndC/AUg5svhe12uXEJHeIrIbGAk8d5V196W2rgkOscQymMFUoAITmMCzPEsssYxi\nFFWp6u/i+dU1XEMf+rCb3bSkJW1pSxvasA0bwN94l9d65KrqOGCciHQEhgDdvLVtE3iWL3fuFx+H\nhTmPw8L+fgxOE8YVrOA93uMHfuBxHmc5y7mFWzBXyk9+nuEZutGNcYwjjDBa0Yr/8B/KUMbfxTM5\ngDuhvw8u6WJYzvXa1cwAPkyxbsp2dVddNyIiIvlxWFgYYSmTwwSclOEu8vcHwEXnOMc0pjGWsZzj\nHM/xHJOYZO3U3VSQggxgAD3owXCGcxu30d91K0ABfxfP+Mny5ctZfvkfm4fcuZCbG9iJcyH3ABAJ\nPKaq21MsU1lVd7setwKGqGq9FBdy78Sp1vkeu5Cb46S80HmIQ4x33WpRi+d5nqY0zdZt0gNhaOU9\n7OFFXmQTm3iXd3mIh/xbIBMQfHIhV1UTgT7AYiAK58LsdhF5VURauhbrIyJbRWQ98DzQ1bXuNuAr\nYBuwCKdVj6V7DrSd7fSkJ1Wpyj72sYxlfMu3NKd5tg78QFGJSnzN14xjHP3oRzjhxBCT/orGXMY6\nZ5kMu9jq5t6FoynVci29XTdfDYfgL4Fwpp/Sec4zilGMZSxDGEIf+pCb3P4ulvED65FrskQCCcxm\nNqMZzQlOsOvp/pyZ0IWCFPR30TJGFRITnbsq5M0Luf8O0UAL/Yt2spOe9OQc55jIRKpjs7YHGwt9\n41OnOc1EJvI2b1OOcrzAC7SiFbklV0CGIn/9BdHRsGcP/PEHxMbCgQNw+LBzP34cTp2CM2ecZM+d\n2/kZH+/8LFAAihQhan8RqjcsDmXLOvcKFaByZedeqZLzIeEnSSTxMR8zmMH0pz8v8iJ5bJqMoGGh\nb3ziIAd5j/eYwAQa0YgXeZH61E9+3+9nwvHxEBUFv/4KmzbBli3OPSEBqlRxgjkkBMqXd0K7ZEko\nUQKKFIFCheCaay45swecdc+ehePHqVH+GFt/OOJ8YOzbBzExsHu3c9+7F6pWhZo1oW5duPNOqFUL\n8mXtxCoxxNCDHpzgBFOYEvT9HoKFhb7xqu1s523eZjaz6UhH+tOfylS+YrksD/1jx+Dnn+Gnn2Dl\nStiwwQn1OnWcwK1ZE2rUgNKlncJlUprHd+YMbNsGGzdCZCSsWeN8s6hfH+67D5o0ccqVy/cXsxVl\nPON5hVd4jdd4hmfSnCDGZH8W+ibTLl6cfYu3iCTSrYuzPg/9+Hgn5L/7DpYuhe3boV49aNgQ7rnH\neXz99T7bvcfHd/w4rFgBy5Y5Zf7rL2jRAsLDoXlzp9rIh3ayk8d5nJKUZBKTKE1pn+7P+I+Fvsmw\neOKZyUze4i1OcYp+9KMrXd26OOuT0D9yBBYuhAULnKC/+WZo1sw5c27QAPLn9/IOry7Tx/e//znH\nMmeO862kRQvo3BmaNoU8vql/jyeeoQzlMz5jMpN5gAcued/dHtUmsFnoG4/FEcdHfMQHfMDN3Ex/\n+tOCFh61rfda6O/fD7Nnw6xZTnVJkybw8MPO2XFp/52tevVD7eBB5xg//9y5NtC5Mzz1lHPtwQeW\nsYwudOFxHud1XicvV1509vs1GZNhFvrGbVFE8R7vMYMZPMzD9KMftaiVoW1lKjSOHHFCfupU2LoV\nWrWCtm3hgQd8Xg2Sliw5E96xAyZPhkmTnOsQvXs7H3KXX1TOpMMcpgtdOMlJpjHtihnHLPSzLwt9\nk6YEEljIQt7nfaKIohe9eJqn+Qf/yNR2PQ6Nc+dg/nyYMsW5GNuiBTz2mFN9k4XVNgHj/Hnn7P/9\n951vAs8/D08+6bQs8pIkkhjFKN7hHSYzmeY0T37PQj/7stA3qTrIQT7lUyYwgRu4gT704VEeJT/e\nCVi3QkMV1q2DTz+Fr76C0FDo2hXatPFquGV7v/wCb70FP/4I//d/8NxzTtNSL1nBCh7jMXrSkyEM\nIRe5LPSzMQt9kyyJJH7gBz7iIxazmLa0pRe9uIM7vL6vNEPj6FH44gv45BM4eRK6d3fCvnzWTGqe\nbe3cCSNGOBeAe/eGAQOgcGGvbPoAB+hAB67lWr7kS4pLMQv9bMpnc+Sa7GMf+xjO8OSLsg1pyO/8\nzsd87JPAj4tzfh49muJFVVi9Grp1g4oVYdUqGDPG6cw0eLAFvjuqVnXq+yMjnZ7ElSvDG284/QIy\nqQxlWMpSqlPdmZqx5qbMl9dkGxb6OcAZzjCd6TSnObdxG7/zO9OYxkY20oc+FMY7Z4iXGzMGatd2\nHoeGwnsjTztn9Hfc4bRKqV7dGQZh2jSno1IWdFDKcW66ybnQu2KF0+O4ShXnwyAxMVObzUteRjOa\n4QyHJU2YxjTvlNcEPKveyaYSSGAZy5jKVOYxj7rUpRvdaE1rruEan+8/Ls4J/JgYqMhv9GYc3XNN\nplDTu8j3/LNO6xsLee/75RenqufsWXjnHWjcONOblNs3UXFTG9rTnmEMsxE7sxGr3snh4olnKUvp\nRS/KUpbBDKYmNdnGNhazmE50ypLAB4jaqlT+YxlzCSeSeiSRi7q6lsh/z3da4Vjg+0aDBk7v5EGD\n4IknnFZPe/dmbpubbyeSSNawhnDCOc5x75TVBCT7ywxwJznJHObwJE9ShjK8xEvcyI2sZjWRRNKf\n/lk7d+q5czBxIg2euZ1xef6PRbQghBgGMorEChWpbqP7+p4ItG/vDEdRubIz3tBbbzmDxGVQCUqw\nmMWEEEJ96rOb3V4ssAkkVr0TYBQliigWs5hv+ZbVrKYBDWhJS1rT+oqONVnm4EEYNw4+/NCps+/X\njzFbmjD2XSEmxhnvrG9f6NfPP8ULatHRTgufw4dhwgRnpE8PXN76ajzjeZVXmc50wgjzblmNV1mT\nzWxIUXaykxWsYLnrVoACNKMZTWlKE5r4dzLxqCh4+234+mvo0MFJ9mrVkt+Oi4PixZ2fRYv6r5iX\nO8MZYly3WGL503U7zGGOcpRjHOMUpzjruiWSSBJJKEo+8pGf/BSgANdzPYUpTFGKUtp1K0c5Qgjh\nRm4khJDAGL9e1enV/MILTpXPf/7jDBnthtSa3C5hCZ3pzDCG0YMePiiw8QYL/QBytW78NZv+Sf67\n1rGOdaxhDatZzXVcRyMaEUYYjWnMTdzk3yFxVZ1Bzt56yxkDp3dveOYZZwz6VPizc89RjrLFddvK\nVnayk13s4i/+ojzlCSGE8pSnjOtWghIUoxhFKUohClHQdctN7uTxhuKJ5zznOctZTnCC4xznL/7i\nEIc4yEH2spffXbeDHOQmbuIWbuF2bieUUO7gDspS1j//IEeOOB/Mq1c7HeHcGDPiar+/XexK/oY5\nkpE213EAstAPEEkkcYADRBPNTnbS6+2dNOnvBNMFLlCb2tShDnWpy13clbV18mmJj4cZM2D0aOfx\ngAHQqVO6Y+BkVeif4xzrWMcqVrGWtaxjHYc4RA1qcBu3UYMaVKMaVahCecpnSUid5SzRRLONbWxk\nIxvYwDrWcS3Xchd3cQ/3cB/3cQu3ZO0H+YIF0KuXM47RiBFpnvWn9fv7i794hEcoRjG+4Auu5Vof\nFdhkhIW+D53nPH/wByc5ySlOcYxj/OW6HeIQBzjAfvbzB38QSyxFKMLN3EwVqjDxX1X5ZpQTTOUo\nF3gTWxw/Dh99BO++67QDf+EFZ2RLNycg8VXon+EMP/NzcrXXRjZSjWo0oAH1qEcd6lCFKgHXxFBR\noolmFatYwQqWspQEEmhKU1rRigd4IGuq7OLioE8fZ/iLKVNSret3p3ruPOd5mqeJIooFLMj0WE3G\neyz0fWgLW2hDGwpRiOu4jsIUprjrVopSydUHFVy3lE0nA3Zsk5gYGDsWPvsMHnzQObMPDfV4M946\nPkXZzGb+y39ZzGLWspZa1CLMdatP/Wx5pqkoe9jDt3zLAhawilU0pCHtaU844RTFxxdDZs50wr9P\nH6epp2sM/zFjnF+/OxfiFeU//IeJTOQbvuFWbvVtmY1bLPQDVMCF/rp1Tn39d985Y+E891ymhkbI\nzPGd5zzLWMY85vEN35CPfLSgBc1oRmMa+/cito+c5CQLWchMZrKUpdzHfXShCw/xEPnw0dy6+/Y5\nYx6dOwdffEHc9Tcmd667KCQE1q+HYsWuvpnP+ZwBDGAGM7iXe31TVuO2jIQ+qur3u1OMnCsgDi8x\nUXX+fNXGjVXLl1cdPVr12DGvbNrT4zutp3WmztQO2kELa2G9R+/RN/VN3aE7NEmTvFKm7OK4HtdP\n9VNtrI21pJbUATpAd+ku3+wsMVH1zTdVS5bUqNdmaa5czu/u4j1XLtWffkp/M8t0mZbUkvqFfuGb\nchq3ubLTs7z1dAVf3C30fej0adXx41WrVFGtXVv1yy9VL1zw6i7cOb5zek7n6Bxtr+31er1eH9AH\n9CP9SA/qQa+WJTvbrbt1oA7UUlpKm2gTXaALNFETvb+jyEhNuPEm/bxQLy3AmeTQDwlRjYtzbxNb\ndItW0Ao6XIcH3Qd1IMlI6Fv1ThbwS/XOgQPwwQfOBdoGDaB/f2jUyO2Ls5642vElkcTP/MwUpjCb\n2dSkJp3oxCM8QglSb/5pnCqvmcxkDGOS5yvuRjcK4MWZxI4fZ2fY0yRG7SA8fibxITd73LluP/tp\nQQvu4i7e472Au6AeDKxOP0BlaeivW+cMxLVwodPcsm9fn82/etHlxxdDDJOZzGd8RkEK0pWudKIT\n5Sjn03LkNIryEz8xilGsZz396U8velEIL006o8rp0eM5/a8Irpn4AYWebOfxJk5wgkd4hOu4jqlM\npSAFvVM24xYL/QDl89BPSIC5c/9uitGnD/TokfYVOS8SgXN6nnnM42M+ZgMb6EhHnuRJalM78Jqo\nZkOb2MQIRrCc5bzAC/Smt9cG17tD1rGuYjsID4dRoyDvlZOnp+UCF+hOd37jNxawgOIU90q5TPos\n9AOUz0L/yBFn/Ppx46BCBeesvk2b5CZ5WWE3u7n5zQmUfPEzbuM2nuIpWtPau1URJtlWthJBBKtY\nxWAG8xRPkRfPQvpyIqBxR+Hxx+HECWc6yzKedRhMIolBDGI+8/mWbwkhJFNlMu6xoZWDxdq1zqxU\nN9/sTKs3bx6sXAnt2mVJ4CeSyDzm0ZSmNKABAD/zM0tZSkc6WuD7UA1qMItZfMM3zGUu1anOLGah\nZPKsomhRpxdv06ZQp44zfLMHcpGLN3iDXvTibu5mM5szVx7jO55e+fXFHWu9k77Tp1UnTlStW9dp\nZvHGG6qHD3thw+47okd0pI7UEA3R+lpfp+gUPatnA6NJapD6Tr/TmlpTG2pDXa/rM7SNK35/ixap\nliqlOm6capLnLXOm63QtqSV1uS7PUHmM+7Amm4EpU4cXFaXat69q8eKqLVqoLligmpDgtbK5Y4tu\n0af0KS2iRbSLdtG1uvaS93P4ry/gJWiCTtAJWlpL61P6lB5Wz04GUv39RUer1qih2r276rlzHpdp\niS7RklpSZ+tsj9c17stI6FudfhbwuE7/zBmYNQs+/tiZTLx7d+fCbMWKPivj5RRlMYt5m7fZzGae\n4Rme5mlKUxq4+iiiYWFuDexofOAYxxjKUKYznWEMozvd3Rp07qr/P0+dcqoR9++H2bM9rudfz3pa\n0pJXeIVe9PJoXeMeu5AboNwKfVWnrn7SJGekywYN4J//hFatPG5NkRnnOc9UpvIWb5GLXAxgAB3p\nSH7yZ1kZTOZsYAO96Y0gfMzHVCft6czS/P+ZlATDhjn9Pb7+GurW9agse9hDM5rxBE/wCq9YSy4v\ns9APUGn+Ue3fD19+CZMnw/nzzplV166ZGgsnI45xjA/5kHd5l5rUZAADaEIT+yPNppJI4iM+YghD\neJZnGcSgq35wu3VSMm+e823z3XedSVo8cJCDPMiD1Ke+deLyMgv9AHXFH9XJk84f0eefQ2QkPPoo\ndOkCDRv6pMdsWvaxjzGMYRKTeIiHeIEXqEnNLC2D8Z197ONZniWaaCYzmbpceabudvXj5s1OW/5O\nneD11yGX+43/TnCC1rSmOMX5gi/sm6OX+KzJpog0F5EdIrJLRAam8n4/EYkSkY0i8r2IlE/xXqKI\nrBeRDSIy15PC5SjnzjlB37EjlCvnVOF06+aMfvjJJz4bIuFqdrGLHvTgNm4jgQQ2sIEpTLHAz2Fu\n4AbmMIfBDKYlLXmZlznP+YxtrGZN5yRlxQpncpbTp91e9XquZxGLAGhBC05wImNlMJmWbuiLSC7g\nfaAZUB14TERuuWyx9cAdqloLmA28meK906paW1VDVbW1l8qdPZw7BwsWMJmuULasMzxC48awZ4/T\nJvqxx9yex9RbNrKRDnTgbu6mHOWIJpp3eIcKVMjScpisIwiP8Rib2MR2tlOXuhlvR1+yJCxZAoUL\nwz33QGys26sWoADTmU5VqnIv93KQgxkrg8mc9Jr3APWB/6Z4/hIwMI3lawE/pXh+0o19eKkBUwA4\nflx1xgzVjh1VCxfWC3c31j68q0e37fdrsX7Wn7WFttAyWkZH62g9qSf9Wh7jH0mapJN1spbQEvqG\nvqEJmpCxJrdJSaqjRqmWLasaGelxGYbqUK2slfU3/S0DOzcX4Yt2+sCjwEcpnj8OvJvG8u8BL6d4\nfgGIBFYB4VdZx7f/Mr72+++qH3yg2ry56nXXOT8//FA/jDigISGaPGzt229nbbGSNEmX6lK9V+/V\nG/VGHa/j9ayezdpCmID0P/2fNtJGGqZhyg2xGd/Q3LmqJUqozpzp8arjdJzeoDfoRt2Y8f0HOb+H\nvuu9VUDeFK+Vcf2sCPwPqJjKej7+p/Gyc+dUly5VfeEFpwNLiRKqXbo4Z/jHj6uq6l9/aXLgpxyv\n/K+/fF+8JE3SRbpIG2gDraJV9DP9TC+od8fQN9lfgiboMB2m/FlKZ+msjG9o/XrVcuVUhw3zuAfv\nV/qV9d7NhIyEvjsDteyDSyp8y7leu4SINAEGAY1UNT5F9dEB18//ichyINQV/peIiIhIfhwWFkZY\noPXwOXLEmVz6+++dcUmqVXPmlf3kE2esktyXNkOLirqyujM2FrZtc6pCfUFRFrKQ13iNs5xlCENo\nS1trImeu4HSuyw28TOjB+3lqWCdG7lnKsLNv07SRh2MnhYbCmjVOn5Jdu5w2/fncm/axHe0oTnHa\n0Y4P+ZBHeMTjYwkmy5cvZ/nFXpEZld6nApAb2A2EAPmAjUC1y5YJdS1T6bLXiwD5XI9LADuBW1LZ\nh48/D70gJka1Z0/na6wbp+txcamf6bs7M5EnkjRJ5+gcDdVQrak1dZbO8s2MSybHOqbHtL2215pa\nU7fr9oxt5NQp1fBwZ0rOI0c8WnWdrtMyWkbH6/iM7TtI4auxd4DmrsCOBl5yvfYq0NL1+HvgAE4r\nng3AXNfrDYDNrtc2Ad2usv2s+PfJcm+/rT6t00/URJ2ts/V2vV1DNVTn6BwLe5NhSZqkH+qHWkJL\n6Jf6ZcY2kpCgOmCA6s03q+7ybK7faI3WSlpJh+pQm4LRTRkJfeuc5WNxcVC8uPOzaFHvbDOJJOYy\nl1d5lTzkYShDaUUr6z1rvGIjG2lHO+7nft7hnYwNlT1hAgwdCjNnOp0O3XSQg7SgBXWowzjGWdVk\nOqxHboDy1iQqSSQxj3m8yqvkJjcRRNCSlhb2xutOcIIe9GA3u5nFLG7iJs83snixMzHLmDHQubPb\nq53kJI/wCNdyLdOYZlMwpsFCP0BlNvQVTT6zt7A3WUVR3uM9hjGMiUzkIR7yfCNRUdCypdP7/JVX\n3O51foELPMmT/M7vzGe+TcF4FRb6ASqjoa8oC1jAUIYiCBFEWDWOyXKrWEUHOtCNbkQQ4XmVy59/\nwsMPQ9WqTmu3/O6Nu5NEEgMZyEIW2hSMV2GhH6A8DX1FWcQihjKUBBKIIIJwwi3sjd8c5CAd6EBB\nCvIlX1KMYp5t4MwZeOIJOHwY5sxxLnS5aSxjeZM3WcACQglNd/lgmuvBQj9AuRv6ivId3zGUoZzh\nDBFE0IajBf5mAAASqUlEQVQ2bk2EYYyvJZDAQAYyhzl8zdfUopZnG0hKgpdegrlz4ZtvnDme3TSL\nWfSmN1/wBU1p6vZ63rqeFqgyEvoeNfXx1Z0c2mTzovQOL0mTdLEu1vpaX6tpNf1Kv7KmlyZgTdNp\nmWvW+dFHqqVLqy73rBfuT/qTltJSOlEnur1ODo8Wa7IZqK52tqEoP/ADQxnKYQ4TQQTtaGfN1EzA\n28xm2tCG1rTmDd4gj1ud+1NYssQZl//NN51Jg9y0k520oAWd6MRrvJZulaed6aeyTiCEbTCG/o/8\nyFCGso99DGUoj/GYhb3JVuKIoxOdiCeeGcygBCU828D27U7Lno4dPZqU5RCHeJiHqUQlJjIxzQlZ\nLPSvZJXFWWwlK7mf++lOd7rRje1s53Eet8A32U4xivEN31CXutSjHpvY5NkGqlWD1avhxx+hQwfn\nYq8bSlGKZSzjPOe5n/s5zOEMlD54WehnkV/4haY05QmeoDOd2cEOutHN86/FxgSQ3ORmJCMZwQia\n0IQZzPBsAyVLwtKlULCgM8HQgQNurXYN1/AVX9GIRtSnPtvZnoHSBycLfR9bwxpY9CAd6Ug72rGT\nnXSnO3nJ6++iGeM1HejAEpbwkuuWSKL7K+fPD599Bq1bw513wvr1bq2Wi1wMZziv8AqNacx3fJfB\n0gcXq9P3kUgiiSCCrWwlttfLnP+wO/lwb7hZY7KrIxyhPe3JT36mMpWieDjg1OzZ0KsXfPghPPqo\n26utZCXtac+/+Bd96Zt8gdfq9K9kZ/peFkkkD/EQj/IorWhFNNEwoZcFvgkKJSjBd3xHVapyJ3ey\njW2ebeDRR+G776BfP+firpuJfQ/38Au/MJnJ/JN/co5zGSh9cLDQ95KUYf8QD7Gb3TzDM2m2LDAm\nJ8pLXt7hHV7mZcIIYz7zPdtA7drOpCzffOO07HHzAm8IIaxkJSc5SSMasZe9GSh9zmehn0mrWc2D\nPEhb2iaHfW96W9iboNeNbixkIc/yLK/yKkkkub9ymTLOGAr58ztDM18+Dd1VFKIQX/EVj/Io9agH\njX7MWOFzMAv9DPqZn2lKUzrQgXDCiSbawt6Yy9SjHmtZy2IW8wiPcIIT7q9coIBzgfexx5wLvCtW\nuLWaIAxkIJOYBDM6MIpRKDm4Yt9DFvoeUJTlLOd+7ucJnqA97Ykmml70srA35ir+wT/4gR8oTWnu\n5E52stP9lUXghRdg0iRo1w7GjXO7nr8ZzaBeJLOZTRvacIxjGTyCnMVC302/8RuNaERPevIET7CT\nnfSgh12gNcYN+cjHBCbQj340pCELWejZBpo1g59/dkL/n/+Ec25eqI2twE/8RAghhBLKKlZ5Xvgc\nxppsuukUp1jAAtrT3q3es8E0vKsxnviFX2hHO57iKYYwxLNRZE+dgu7d4bff4OuvoUKFNBdP2WRz\nHvPoSU+e4zle4qUc0Qvext4xxmQLBzhAO9pRlKJ8zucUoYj7K6vC22/D6NEwZQo88ECqi6U2P/Ve\n9vI4j5NEEp/xGRWp6IWj8R9rp2+MyRbKUIZlLONGbqQuddnCFvdXFoEBA2DaNGeEztdfd8bqT2HM\nGKflJ0BoqPMcoBzlWMYywgmnHvX4hE+C7iKvnekbY/zqcz6nP/0Zy1g60cmzlffvd9ryX3stfP45\nlChBXJwT+DExfy8WEuKM7lAsxYRfW9lKV7pSjGJMYELGJn/3MzvTN8ZkO0/wBEtZyiu8wv/xf1zg\ngvsrly3rDNhWs6ZzSr9yJVFRVzbrj42FbZd1Dq5BDdawhmY0ox71GM1o4onP/AEFOAt9Y4zf1aQm\nv/IrscTSiEb8wR/ur5w3L7zxBkyYAG3bUmfxcELKXTrgW/nyUL36lavmIQ8v8AJrWMNSlnI7t7OE\nJZk8msBmoW+MCQhFKMIc5vAoj1KXuvyX/3q2gRYt4NdfKfjTYpbne4B6N+wDnKqdvn3/vpibmkpU\nYhGLGMEInuIpWtM6xw7XbKFvjAkYgvAiLzKLWfSkJ4MYRAIJ7m+gXDlYupQKXe/jlwu1CWcuGzY4\n47e5s+9wwtnGNu7mbhrTmB708OxbRzZgoW+MCTgNach61rOBDYQR5tngablzw+DB5Jo3l7cYQNEX\n/gknT7q9ekEK8iIvsotdlKIUoYTyJE+ygx0ZOJLAY6FvjAlIJSnJIhbRkpbUoY7no3U2aEAtNjpz\n79aqBT/95NHqRSjCcIazm91UohKNaEQLWrCABZ5NEhNgrMmmMSbgrWIVnehEK1rxJm9SgAJurZfc\nI3f+fGdylvbtYfhwuOYaj8twlrPMYAbjGc8BDtCJTnSkI7dze/KkLVnNmmwaY3Kku7iLjWzkT/6k\nHvXYylbPNvDww7BlCxw+DLff7kzG7qGCFKQb3VjDGhawAEFoQxuqUpW+9OUbvuEUpzzeblazM31j\nTLahKJOZzL/4F0MYQh/6pDl2T6rTJc6bB336wIMPwqhRUMSDISBSKc961rPYdYskkpu4ibrUpSY1\nqeS6laEMhSmcalmTSOI85ylIQY/3b2PvGGOCwm5205nOFKYwk5jEDdyQ6nJXnSP3+HEYNMj5ABg9\n2unVK5mvornABbawhbWsJYoo9rhuf/InpzlNYQqTj3zJ4X/KdWtAA1ay0uP9WegbY4JGAgmMYATv\n8R5jGUtHOl5Rt57uxOi//AK9ezuN+D/4AKpV81l544nnGMdIIIEkklCU67iOQhTK8IifFvrGmKCz\njnV0oQu3cAvjGEdpSie/l27oAyQkOOP0v/YaPP44DB2adk+uAGIXco0xQecO7mAd67iZm7md25nO\ndM9GzsyTB557zhmc59w5uOUW56w/PmeOw2Nn+saYHGMNa+hOdypTmfGM5wYp6+7sin/btAkGDoQ9\ne2DYMGjb1mnrH4B8Vr0jIs2Bd3C+GXyqqm9c9n4/oAcQDxwGuqtqrOu9rsC/AQWGqeqUVLZvoW+M\n8YrznGc4wxnHOI70HEbiRz08m53roqVLnfBPTISICKfZZyoXe/05S55PQl9EcgG7gPuB/cBaoKOq\n7kixTGNgjaqeE5FeQJiqdhSRosCvQG1AgHVAbVU9ftk+LPSNMV61hS3U/KUndzfIzQQmUJ1UhtlM\nj6rTwiciwhne4d//htatr3rm79Y1BC/yVZ1+PSBaVWNUNR6YDoSnXEBVf1TVizMVr4bk9lPNgMWq\nelxVjwGLgeaeFNAYYzLiNm6De1bSmc6EEcYgBnGa055tRMQJ+fXrYfBgGDkSbr0VPv3U/cnZA4w7\noX8DkHJKgr1wlUaxjn9C8piol6+7L511jTHGe5Jy8wzPsJnNxBDDrdzKXOZ6PkVirlzQpg2sWQPj\nx8OsWc6YzYMHw14PBoMLAHm8uTEReRy4A2js6boRERHJj8PCwgjzdWWYMSZolKEMU5nKMpbxLM8y\ngQm8wztUpapnGxKBe+917jt3wvvvw223QcOG0KMHuWmBl2P1EsuXL2f5xQsIGeROnX59IEJVm7ue\nvwRoKhdzmwBjgUaq+pfrtY449fu9XM8/BH5Q1RmXrWt1+sYYr0utjv0CF3iP9xjBCJ7kSQYzmMIU\nzvhOTp2CmTPhk0/4c9Ue/vFcB+jcGerW9Uov37T46kJubmAnzoXcA0Ak8Jiqbk+xTCgwE2imqntS\nvJ7yQm4u1+M7XPX7KfdhoW+M8Qp3W9P8yZ8MYhDf8i2v8Rrd6Z7hnrEXVZFd7IqYBl9+6bTzb93a\nud99t9MfwMt83WRzLH832RwpIq8Ca1V1oYh8D9TA+VAQIEZVW7vW7cbfTTb/Y002jTGBZB3reJ7n\nOclJRjOaJjTJ8LaSv1mowtatMHeuc//tN6dK6IEHoHFjpwOYF9r+2zAMxhiTAYoym9m8xEtUpSqj\nGJWhJp5XbbJ58CAsWQLff+9M5nLsGNx1F9SpA6GhULu2M9Wjx/uz0DfGmAy7wAU+4ANGMIJwwnmV\nVylLWbfXd7ud/v79sGqV0xR0/Xqn6mfhQo/La6FvjDFecIxjjGQkH/MxT/M0/+JfFCHtcffj4qB4\ncednVo3XZgOuGWOMFxShCCMZyUY2cohDVKEKb/ImZzmb6vJjxjg1NODU1owZk4WF9ZCd6RtjTDq2\ns53BDGYNaxjMYLrTnXzkA5wz+9q1ISbm7+VDQpxam2LFfFsuO9M3xhgfqEY1ZjObOa5bNaoxhSkk\nkkhUFMTGXrp8bKwzUnMgstA3xhg31aUu3/EdE5nIx3xMDWqwO3Qm5SokXbJc+fJQPQPju2UFC31j\njPFQYxqzghWMZSzjC71JfGQtSvacAyghIdC3b+BOvmV1+sYYkwmK8g3f8O+EV9i8Rfn85qF0LhR+\nxXy9vmBNNo0xxk8UJdfDCwidH0ESSQxhCG1ok7EJXNxkoW+MMX4kAkmqLGQhr/EaZznLv/k37Wmf\n6XF9Ut+fhb4xxvhNyh65irKYxbzO6xziEAMZyBM8kdzU0zv7s9A3xhi/SW0YBkX5iZ8YznC2spV+\n9KMnPbmO67ywP2unb4wxAUUQGtGIb/mW+cznV36lIhUZxCD2sz/Ly2Ohb4wxWaQ2tZnGNNayllOc\nogY16EIXNrAhy8pg1TvGGOMlbo+y6XKUo3zMx2xmM1/wRQb2Z3X6xhjjN56Gfub3Z3X6xhhj0mBn\n+sYYkwnuzsnrC1a9Y4wxQcSqd4wxxqTJQt8YY4KIhb4xxgQRC31jjAkiFvrGGBNELPSNMSaIWOgb\nY0wQsdA3xpggYqFvjDFBxELfGGOCiIW+McYEEQt9Y4wJIhb6xhgTRCz0jTEmiFjoG2NMELHQN8aY\nIOJW6ItIcxHZISK7RGRgKu83FJF1IhIvIo9c9l6iiKwXkQ0iMtdbBTfGGOO5dENfRHIB7wPNgOrA\nYyJyy2WLxQBdgS9T2cRpVa2tqqGq2jqzBc6Oll+cSy2HsuPL3nLy8eXkY8sod8706wHRqhqjqvHA\ndCA85QKq+oeqbgVSm/PQo6m8cqKc/h/Pji97y8nHl5OPLaPcCf0bgNgUz/e6XnNXfhGJFJFVIhKe\n/uLGGGN8JU8W7CNEVQ+ISEVgmYhsVtX/ZcF+jTHGXEZUU6uRSbGASH0gQlWbu56/BKiqvpHKspOA\nBar69VW2ler7IpJ2IYwxxqRKVT2qQnfnTH8tUFlEQoADQEfgsTSWTy6AiBQBzqjqBREpAdwFXPFh\n4WmhjTHGZEy6dfqqmgj0ARYDUcB0Vd0uIq+KSEsAEakjIrFAW+BDEdniWr0a8KuIbACWAiNUdYcv\nDsQYY0z60q3eMcYYk3METI9cERklIttFZKOIzBaR6/1dJm9Ir2NbdiYi5URkmYhEicgWEXnO32Xy\nNhHJ5epcON/fZfE2ESksIjNdf3dRInKnv8vkTSLST0S2ishmEflSRPL5u0yZISKfishBEdmc4rWi\nIrJYRHaKyHciUji97QRM6ONUH1VX1VpANDDIz+XJNDc7tmVnCUB/Va0ONACezWHHB9AX2ObvQvjI\nWGCRqlYDbge2+7k8XiMiZYH/A2qrak2c65cd/VuqTJuEkyUpvQQsUdWqwDLcyM2ACX1VXaKqSa6n\nq4Fy/iyPl6TbsS07U9U/VXWj6/EpnNDwpA9HQBORckAL4BN/l8XbXN+kG6rqJABVTVDVE34ulrfl\nBq4VkTzANcB+P5cnU1R1JXD0spfDgc9cjz8D0h31IGBC/zLdgf/6uxBekNmObdmGiNwI1ALW+Lck\nXjUGeJHUe5pndxWBIyIyyVV99ZGIFPR3obxFVfcDbwF/APuAY6q6xL+l8olSqnoQnJMwoFR6K2Rp\n6IvI9676tYv3La6frVIs828gXlWnZmXZTMaJSCFgFtDXdcaf7YnIQ8BB1zcZIecNJ5IHqA18oKq1\ngTM4VQU5gqu5eDgQApQFColIJ/+WKkuke4KSFT1yk6nqA2m9LyLdcL5O35clBfK9fUCFFM/LuV7L\nMVxfnWcBn6vqPH+Xx4vuBh4WkRZAQeA6EZmiql38XC5v2QvEquqvruezgJzU0KAJ8JuqxgGIyNc4\n/YRy2snkQREpraoHReQfwKH0VgiY6h0RaY7zVfphVT3v7/J4SXLHNlfLgY5ATmsFMhHYpqpj/V0Q\nb1LVl1W1gqrehPN7W5aDAh9XlUCsiFRxvXQ/OeuC9R9AfREpICKCc3w54UL15d865wPdXI+7Aume\neGXpmX463gPyAd87vyNWq2pv/xYpc1Q1UUQudmzLBXyqqjnhPx4AInI30BnY4uqAp8DLqvqtf0tm\n3PQc8KWI5AV+A570c3m8RlUjRWQWsAGId/38yL+lyhwRmQqEAcVF5A9gKDASmCki3XGGuG+f7nas\nc5YxxgSPgKneMcYY43sW+sYYE0Qs9I0xJohY6BtjTBCx0DfGmCBioW+MMUHEQt8YY4KIhb4xxgSR\n/wekUNG+O7BAPAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7fb9f3824470>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}