{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      ":0: FutureWarning: IPython widgets are experimental and may change in the future.\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "import statsmodels.api as sm\n",
    "import george\n",
    "#import scipy.optimize as op\n",
    "import math\n",
    "import emcee\n",
    "import matplotlib.pyplot as pl\n",
    "from __future__ import division, print_function\n",
    "from george import kernels\n",
    "import triangle\n",
    "from lmfit import Parameters, minimize, fit_report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xval = np.array([262.5, 283.0, 305.0, 327.5, 351.5, 377.0, 403.5, 431.0, 459.5, 489.0, 519.5, 550.5, 582.5,\n",
    "                   616.0, 650.5, 686.5, 724.0, 762.5, 802.0, 843.0, 885.5, 929.5, 975.5, 1023.5, 1073.0, 1124.0, \n",
    "                   1176.5, 1231.0, 1287.5, 1346.0, 1406.5, 1469.0, 1534.0, 1601.0, 1670.5, 1742.5, 1816.5, 1893.0, \n",
    "                   1972.0, 2053.5, 2138.0, 2225.0, 2314.5, 2407.0, 2502.5, 2600.5, 2701.5, 2806.5, 2915.0, 3027.0, \n",
    "                   3143.0, 3263.0, 3387.0, 3515.5, 3648.5, 3785.5, 3927.0, 4074.0, 4226.0, 4383.0, 4546.0]) #mjj in GeV\n",
    "\n",
    "#half the bin width\n",
    "xerrminus = np.array([9.5, 11.0, 11.0, 11.5, 12.5, 13.0, 13.5, 14.0, 14.5, 15.0, 15.5, 15.5, 16.5, 17.0, 17.5, \n",
    "                        18.5, 19.0, 19.5, 20.0, 21.0, 21.5, 22.5, 23.5, 24.5, 25.0, 26.0, 26.5, 28.0, 28.5, 30.0, \n",
    "                        30.5, 32.0, 33.0, 34.0, 35.5, 36.5, 37.5, 39.0, 40.0, 41.5, 43.0, 44.0, 45.5, 47.0, 48.5,\n",
    "                        49.5, 51.5, 53.5, 55.0, 57.0, 59.0, 61.0, 63.0, 65.5, 67.5, 69.5, 72.0, 75.0, 77.0, 80.0, 83.0])\n",
    "\n",
    "#half the bin width\n",
    "xerrplus = np.array([9.5, 11.0, 11.0, 11.5, 12.5, 13.0, 13.5, 14.0, 14.5, 15.0, 15.5, 15.5, 16.5, 17.0, 17.5, 18.5, \n",
    "                       19.0, 19.5, 20.0, 21.0, 21.5, 22.5, 23.5, 24.5, 25.0, 26.0, 26.5, 28.0, 28.5, 30.0, 30.5, 32.0, \n",
    "                       33.0, 34.0, 35.5, 36.5, 37.5, 39.0, 40.0, 41.5, 43.0, 44.0, 45.5, 47.0, 48.5, 49.5, 51.5, 53.5,\n",
    "                       55.0, 57.0, 59.0, 61.0, 63.0, 65.5, 67.5, 69.5, 72.0, 75.0, 77.0, 80.0, 83.0])\n",
    "\n",
    "yval = np.array([4.917897773E8, 3.799410724E8, 2.571684927E8, 1.841859503E8, 1.35935548E8, 9.41937297E7, \n",
    "                   6.79501892E7, 4.90487119E7, 3.49414712E7, 2.5617166E7, 1.85988426E7, 1.31836779E7, 1.0136639E7, \n",
    "                   7453854.2, 5553447.2, 4264932.3, 3163415.4, 2362702.7, 1775332.4, 1366279.0, 1026756.5, 788737.4,\n",
    "                   603941.0, 461687.4, 346323.9, 264838.0, 198505.7, 153246.1, 115508.6, 89024.0, 66542.8, 51185.0, \n",
    "                   38806.7, 29264.0, 21907.5, 16796.6, 12379.0, 9150.3, 6763.4, 5155.4, 3874.5, 2786.4, 1945.9, \n",
    "                   1482.2, 1085.5, 795.8, 538.7, 379.5, 275.1, 174.8, 136.3, 87.5, 59.7, 37.7, 17.0, 20.2, 13.0, \n",
    "                   6.7, 4.0, 3.0, 1.0])\n",
    "\n",
    "yerr = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, \n",
    "                        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\n",
    "                        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, \n",
    "                        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])\n",
    "xerr = xerrplus+xerrminus\n",
    "newy = yval/(xerrplus+xerrminus) #scale data by bin width since curve-fit can't have the model take in an argument\n",
    "#that isn't a fit parameter\n",
    "#yerrplus = array(\"f\", [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, \n",
    "                   #    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, \n",
    "                   #    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, \n",
    "                   #    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])\n",
    "#numpoints = 61\n",
    "#p8757_d1x1y1 = ROOT.TGraphAsymmErrors(numpoints, xval, yval, xerrminus, xerrplus, yerrminus, yerrplus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10fba7650>]"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYlWX9x/H3BwHFFVdKwSUERSxCXNBQ0SwFE9zF3EDF\n",
       "3dTUckkfntzKtMwWNDWX3PqlllSWZoqVqeFG7nu5L4k7uADf3x/3AYeZYebMcOY8M+d8XtfFpeec\n",
       "+5zz9bmYz9zez70oIjAzs9rRregCzMysshzsZmY1xsFuZlZjHOxmZjXGwW5mVmMc7GZmNabVYJd0\n",
       "lKSHJD0s6aiFtDlf0lOSpksaWvkyzcysXC0Gu6T1gQOBjYAhwNck9W/UZjSwdkQMAA4CJndQrWZm\n",
       "VobWeuzrAvdExIcRMQe4A9i5UZsxwOUAEXEP0FtSn4pXamZmZWkt2B8GNpe0gqQlge2Bvo3arAa8\n",
       "0ODxi820MTOzKune0osR8bik7wO3AB8ADwBzm2mqxm+tTHlmZtZWLQY7QET8EvglgKQzgecbNXkJ\n",
       "6Nfgcd/ScwuQ5LA3M2uHiGjceW5Rq8EuaZWIeF3S6sBOwCaNmkwBjgCulTQceDsiXqtEcbVK0qSI\n",
       "mFR0HZ2Br8WnfC0+5WvxqfZ0ilsNduA6SSsCnwCHRcS7kg4GiIgLI+ImSaMlPU0arpnQ1iLMzKxy\n",
       "yhmK2aKZ5y5s9PiIShZlZmbt55WnxZhadAGdyNSiC+hEphZdQCcytegCujJV66ANSeExdjOztmlP\n",
       "drrHbmZWYxzsZmY1xsFuZlZjHOxmZjXGwW5mVmMc7GZmNcbBbmZWYxzsZmY1xsFuZlZjHOxmZjXG\n",
       "wW5mVmMc7GZmNcbBbmZWYxzsZmY1xsFuZlZjHOxmZjXGwW5mVmMc7GZmNcbBbmZWYxzsZmY1xsFu\n",
       "ZlZjWg12SSdKekTSQ5KulrR4o9dHSnpH0gOlP9/puHLNzOqHcvVsz/u6t/ih0prARGBQRHwk6dfA\n",
       "OODyRk3viIgx7SnAzMwWpFxLAAcD327P+1sMduBd4BNgSUlzgCWBl5qroz1fbmZmn1Ku7sB+QAY8\n",
       "CGwHTG/r57QY7BExQ9K5wPPALODmiLi1cTNgM0nTSaF/XEQ82tZCzMzqlXJ1A3YDvgu8DOwRWdwF\n",
       "oElt7zcrIhb+otQf+D2wOfAO8Bvguoi4qkGbZYA5ETFT0ijgxxExsJnPCiBv8NTUiJja5orNzGqE\n",
       "cgkYBZxBGh05iZzZBCMbNMsiok3p3lqw7wF8JSIOLD3eBxgeEYe38J7ngGERMaPR89HW4szMapVy\n",
       "bQGcCSwPnAzcGFnTQG5PdrY2xv44cIqkXsCHwDbAvxp9aR/g9YgISRuTflnMaPpRZmamXMNIPfSB\n",
       "pLH0qyOLOZX8jtbG2KdLugK4F5gL3A9cJOng0usXArsCh0qaDcwkzZoxM7MGlGsQcBqwKXA6cElk\n",
       "8XGHfFdLQzEV/SIPxZhZHVKuNUk98+2Bc4CfRhYzy3qvWBb0TqWHYszMrB2U6zPAScBewM+AAZHF\n",
       "O2W9VwwDDgd2as93O9jNzCpIuZYHjictMPoVMCiyeL3V94klSFMeDwc+C0wmjcO3+t7GHOxmZhWg\n",
       "XEsB3wC+CfwOGBpZPN/q+8SapF8C+5MWJZ0J/DGCOaXX28zBbma2CJRrceAg0rDLHcCIyOKJFt8j\n",
       "ugFfBQ4DNgOuADaP4MlK1ORgNzNrh9Ly/31IN0YfAUZFFg+2+B6xPDABOBR4nzT2Pi6Csm6mlsvB\n",
       "bmbWBqXl/zuTpi6+DuwVWdzZ4nvEUNLY+S7AH4F9gbsj6JBpiQ52M7MylJb/b0taXARwNHBLc6tF\n",
       "ASQWJ63zORzoC1wArBPR9puhbeVgNzNrhXKNIN3UXBn4DnBDC4Hej3Qz9EDgYeBs4A8RzK5SuQ52\n",
       "M7OFUa6hpFWig4FJwJWRRZOAlhAwEjgC2Aq4EhgZweNVK7YBB7uZWSPKNZC0he6WpKGXnSOLj5q0\n",
       "E8uQbqDO2xjxp8D4CN6rVq3NcbCbmZUo1+rAqcBY4EfAAZHFB03aiUGkMP86cBuppz61o26GtpWD\n",
       "3czqnnKtQpqHvg/pJufAyOKtBdqI7sAOpBAfDFwEfCGCF6tcbqsc7GZWt5SrN3AsaaHQlcB6kcVr\n",
       "C7QRq5BuhB4CvEAabrk+gg7ZmbESHOxmVneUa0ngSOA40ilxG0QW/12gjdiY1DvfAbgeGBvBA9Wu\n",
       "tT0c7GZWN5SrJzCRdGLRP4DNI4v5M1dKc893J4X+yqSVocdE8GYB5babg93Map5yLUbaPjcHHgO+\n",
       "FlncP/910Zc01DKRtBHXacBN8zbi6moc7GZWs0qrRXcizUV/E9gvsvgbzJ97PoLUO98GuArYsqi5\n",
       "55XkYDezmlMK9G1Iq0W7k26Q/jmyCIlepGmKRwJLkIZbDozg3aLqrTQfjWdmNUW5NiUF+meBU4Dr\n",
       "I4u5EmuQZr/sD9wN/AS4NYK5hRVbhvZkp4PdzGqCcn2BtEr0C6RVo5czKeYAW5N65yOAy4GfR/BM\n",
       "YYW2UXuy00MxZtalKdcA0k3RrUk99V2ZFN1Jc8+PKDX7CbBXBE1WkdYiB7uZdUnK1Y801LIzafn/\n",
       "QUyKVUjhvh/wd9JRdbd3lqX+1eJgN7MuRblWBk4ExgO/4KOlB3LWe8OAq0nHzF0CbBjBfworsmCt\n",
       "BrukE4G9gbnAQ8CEiAV3OZN0PjAKmAmMj4gusTrLzLoO5VqOdFD0EcA1PLHDRlwz5aukhUZzgPPp\n",
       "gGPmuqIWg13SmqQJ+4Mi4iNJvwbGkW5AzGszGlg7IgZI2gSYDAzvsIrNrK4oVy/STorfAv7E30/Y\n",
       "ib+eNRb4F+nw6MPpRDsrdgat9djfBT4BlpQ0B1gSeKlRmzGUgj4i7pHUW1KfiAU30jEzawvl6gEc\n",
       "AJxCcDd//tEJ3HP0WNK+Lb8ENojgvy1+SJ1qMdgjYoakc4HngVnAzRFxa6Nmq5F2PJvnRdL5fg52\n",
       "M2uz0vL/PYGcOd2f448/v5z7J44FBpKGW/b0cEvLWhuK6U86sHVN4B3gN5L2ioirGjdt9HghZwFq\n",
       "UoOHUyNialuKNbPaVVotOgY4nU96fcQffzaNBydsA7xPJzvIoiNJGkk6Zq/dWhuK2RD4Z0S8WfrC\n",
       "G0h3nRsG+0tAvwaP+9J0uAaAiJjU7krNrGYp15cJzuTjZVbgz+e9wQPj14FudwAbR/Bs0fVVU6nD\n",
       "O3XeY0lZWz+jtWB/HDhFUi/gQ9LeC/9q1GYK6bfptZKGA297fN3MyqFcmzC321l8vNxgbjlnFg9M\n",
       "+JBY7ErgigjeL7q+rqq1Mfbpkq4A7iVNd7wfuEjSwaXXL4yImySNlvQ08AEwoaOLNrOuTbnW5+Ol\n",
       "zmHucpvy1zPE/Qf+kzmLnwfc0tn3bukKvFeMmVWNcvXn3VV/Qs8PtuLvJwX3HnQFH/U+rxa2yu0o\n",
       "3ivGzDolfXP11Zm5wkUsv9zWTN/3faYddjrv9vt5BG+1/m5rKwe7mXUYfeWE/nzmwUtY7b3NeXrb\n",
       "l5ly8QG8suHVEcwuurZa5qEYM6s4rTV1GEMun8w6Uzbk+RFP8ugu34jp+95SdF1dkfdjN7PCSHRj\n",
       "mZe2Z+Ofnc0GFw9kRv+HeXPgQfHby+8purauzGPsZlZ1EkvS4/3xbHLxKWz+vRX5aJmHmL3EhnHx\n",
       "Xd4MsCAOdjNrF4m+aPbhDLn6cL58YjcW++Qxes0YGz94tfFaF6syB7uZtYnEhjD3GNaZMobRR37I\n",
       "km8+Q49Zx0YWtxVdmyUeYzezVkksBnwN+Cb9b16HMRM/ZJmXP6DbnJOAKZFVKUjqkMfYzayiJJYi\n",
       "nVR0NGtM/Yid95nDsi9+gDgVuDaymFNshdYc99jNrAmJVYEjgQPp+88H2GPXHiz9ygDEacAvI4tP\n",
       "Ci6xbni6o5ktEokhpOPndmDVaVPYe9QyLPnm5sDZwM8ii1nFVlh/HOxm1mYSArYFjgXW4zMPXsb+\n",
       "I/rQ84OdgJ8AP4ws3i20yDrmMXYzK5vE4sBepB76HFZ++AIO2eABFvvkEOAyYN3I4o0ia7T2cbCb\n",
       "1RmJFYBDSOco/JsVHz+BI9b7IorTgBuAIZHFi4UWaYvEwW5WJyTWAo4B9gamsOLjO3DkoBHAxcBt\n",
       "wKaRxVNF1miV4WA3q3FpQRHHkU5Au5jez36Ro/t/Bfgt8CCwbWQxvcgarbIc7GY1qMEN0W8BawPn\n",
       "0fu5gzn6c9sBfwFeBvaILO4qsEzrIA52sxoi0QMYBxxfeupsVn7k1xy+/ldIByTPJo2t3+rVorXL\n",
       "wW5WAySWAQ4kjaE/Teqp38wkjSCNn68AfAf4nQO99jnYzbowiT6kFaKHALcDu0QwTbmGATcB6wKT\n",
       "gCu9/L9+ONjNuiCJtUkLisYB1wLDI3hauQYp5zfAl4DTgbGRxccFlmoFcLCbdSESw4BvA1sDFwDr\n",
       "RvCacq2hnEuB7YFzgP0ii5kFlmoF8pYCZp1caYbL1sAJpKGVHwIXR/CecvUBTiatIJ0MnBNZvF1Y\n",
       "sVZxHbJXjKR1SP+rN8/ngFMi4vwGbUYCNwLPlp66PiJOX9TizOpZaQ/0nUg99KWB7wNXR/Cxci1P\n",
       "mpt+CHAFcFZk8XphxVqH6fBNwCR1A14CNo6IFxo8PxL4ZkSMqWRxZvVIoiewD2lmy1vA94ApEcxV\n",
       "rqWAb5D2d7kR+G5k8XxhxVqHq8YmYNsAzzQM9Ybf38bPMrMGJJYGDiKF9sPAwcAdEYRyLa6cg4CT\n",
       "gL8BIyKLJ4qr1jqztgb7OODqZp4PYDNJ00k9+uMi4tFFLc6sHkisSJqyeBhpEdGYCO4HUK7uytmb\n",
       "NGXxEWB0ZPFAQaVaF1F2sEvqCexAGu9r7H6gX0TMlDQK+B0wsJnPmNTg4dSImNqmas1qSOmUom8C\n",
       "+5N2VRwRwZMAytUN2Bk4DXgD2Duy+EdRtVr1lIa2Ry7SZ5Q7xi5pLHBoRGxXRtvngGERMaPBcx5j\n",
       "NwMk+pPGz3cDLgfOjeBFAOUS8FXgzFLzk4GbvVq0fnX0GPuewDUL+eI+wOsREZI2Jv3CmNFcW7N6\n",
       "JTGYNEa+LfBzYGAE/5v/eq4vkQJ9FdLy/xsc6NYeZQW7pKVIN04nNnjuYICIuBDYFThU0mxgJmks\n",
       "3syYv6joZGAz4EfAoRHMP2pOub4InAEM5tPl/7MLKNVqhBcomXUQiRGkQP886TDoiyOYvxpUuQYC\n",
       "3wW2JPXUfxFZfFRErdZ5+cxTs4KVVol+mTSUsjppDvqOEcwPbOVaHTgVGEtaRXpAZPFBAeVajXKw\n",
       "m1VAKdBHkwJ9edLQyjURzB9SUa5VgBOBfUn7vAyMLN4qoFyrcQ52s0Ug0Q3YkRTo3Uk7Kl4fwfwt\n",
       "cpVrOdLy/8OAq4DBkcWrBZRrdcLBbtYOpX1cdiMF+iwgB34fwdz5bXItSTqt6DjgD8AGkcV/CyjX\n",
       "6oyD3awNJLoDXydNW3yTFNo3RzB/FoJy9SSdZnQycCewRWTxeAHlWp1ysJuVoXSW6L6kQH+BNKxy\n",
       "e6NAX4wU+jnwBDAmsrivgHKtzjnYzVpQ2mlxPOmm5zPAhAj+tkCbtFp0R9L4+lvA+Mjib5gVxMFu\n",
       "1gyJxYEJpEB/HNg7gjsXaJMCfRvSDJgewPHAn7xa1IrmYDdrQGIJ4ADSaUUPA+MiuKtJu1ybkhYV\n",
       "rQqcAlwXWcxt3M6sCA52M+YH+oGkQJ8O7BLBv5q0y/UFUg99CGn5/xVe/m+djYPd6lop0CeSAv1+\n",
       "YKcIpjVpl2sA6abo1sBZwG6RxYfVrNWsXA52q0uNAv0+0uEWTWawKFdf0vL/nUkbeB0UWbxfzVrN\n",
       "2srBbnWl0ZDL/Sw80FcutZkA/IK0/N9bUVuX4GC3ulCa5XIAaZbLdNLGXPc2aZdrWeBY0orRa0nL\n",
       "/1+pZq1mi8rBbjWtwTz075BmuSzspmgv4HDSyUZ/AjaKLJ6tYqlmFeNgt5pUWim6D2kq4pPA7hHc\n",
       "3aRdrh6knvwpwN3AVpHFI9Ws1azSHOxWU0qbc+1FuuH5X2CfCJocAl1a/j+OdNDFs8COkUWT2TBm\n",
       "XZGD3WpCafvc3Uhzy/8HTIzg9ibt0mrRMaTl/+8DB0YWTdqZdWUOduvSSgdcjCX1vGcBRwF/abg5\n",
       "1/y2ubYmrRbtRbqJ+kcv/7da5GC3LqkU6NuSet6LkXZd/ONCAn0T0mrRNUlj6b/28n+rZQ5263Ik\n",
       "RpICfXkgA25oeMDF/Ha51i+1G0bq0V8WWXxSxVLNCuFgty5DYjhwGrAWaSz9moZH0M1vl6t/6fWv\n",
       "At8H9owsZlWvUrNiOdit05MYQgr0L5b+eVkETXreyrUaab76bsBPgMMji3erWatZZ9BisEtah7T6\n",
       "bp7PAadExPmN2p0PjAJmAuMj4oFKF2r1R2IgaeOtrUgbb+0eQZONt5RrRdLy/wOAS4B1Ios3q1mr\n",
       "WWfSYrBHxBPAUABJ3YCXgN82bCNpNLB2RAyQtAkwGRjeMeVaPZBYnTQPfSxp462JETTZeEu5lgGO\n",
       "Ab4B/Ab4fGTxUjVrNeuM2jIUsw3wTES80Oj5McDlABFxj6TekvpExGuVKtLqg8TKpNkt+wIXAAMj\n",
       "eKtJu1xLAIeSeum3AsMji6erWatZZ9aWYB8HXN3M86uRDved50WgL+Bgt7JILAt8EzgSuApYL6Lp\n",
       "35/S8v/xpN78fcA2kcVDVSzVrEsoK9gl9QR2AL69sCaNHje76EPSpAYPp0bE1HK+32pTaQvdQ0iL\n",
       "hW4GNozguSbtcnUD9iBNWXwe2DWyuKeatZpVi6SRwMhF+Yxye+yjgPsi4o1mXnsJ6Nfgcd/Sc01E\n",
       "xKQ2VWc1qcF+LqeRttDdJoImPe/S8v/tSYuLPgQOjSxurWatZtVW6vBOnfdYUtbWzyg32PcErlnI\n",
       "a1Mo7V0taTjwtsfXrTml1aKjgO8B7wF7NbdBF4ByjSQt/18WOBmY4uX/ZuVRtPKzImkp0i55a0XE\n",
       "e6XnDgaIiAtLj38KbAd8AEyIiPub+ZyIiMZDNlYnJDYhLRbqQ7rpOWUhy/83IvXQ+5PG0q+NLJos\n",
       "QjKrF+3JzlaDvVIc7PVJYm1Sz3sz0mrQyyKY3aRdrvVIy/83IQ3RXOLl/2bty06vPLUOIbESacOt\n",
       "vYAfAuMjmNmkXa552wOMAs4G9vLyf7NF42C3ipLoBRxNOjf0WtLUxdebtMv1WdLy/3HAT4EBkcU7\n",
       "1azVrFY52K0iSgdd7EMaRvkXsFkETzZpl2sF0rmiE0kL29aNrNnZVmbWTg52W2QS2wA/IB10sWcE\n",
       "dzZpk2tp0iEYxwDXA0MiixerWqhZnXCwW7tJfJ40Lj6ANNPl+sYzXZRrceBg0iKk24BNI4unql2r\n",
       "WT1xsFubSaxKWgW6A2lq4tgIPl6gTa7upD1fMtIipO0ii+nVrtWsHjnYrWwSSwPHkfZ0uRhYJ4K3\n",
       "F2iTlv/vShprfxkYF1ncVe1azeqZg91aVdoCYAKpl347MCyC/yzQJi3/347Ug59DWo18q1eLmlWf\n",
       "g91aJLEtcA4wgzTkMq1Jm1ybkxYhrUiawvhbB7pZcRzs1iyJwaRA70+annhjMzdGNyD10NcljaVf\n",
       "5eX/ZsVzsNsCJPqQjqPbmRTak5u5MbouaVhmBGkbgLGRxceNP8vMiuFgN2D+3uhHk26O/gpYN4IZ\n",
       "C7TJtQapZ/41Um9+fGTRZJsAMyuWg73OlbbS3YU0H/3fwKYRLDDPXLn6kI6s2xv4OTAwsni78WeZ\n",
       "WefgYK9jEhuQDovuDRwYwW0LvJ6rN3A86ZSjK4BBkUWTfV/MrHNxsNchic+Qxs+3J+15fkkE8296\n",
       "KtdSpLnqxwI3AkMji+eLqNXM2s7BXkckFieNox8PXEpaYDR/R8XS8v+JpGGXvwEjIosniqjVzNrP\n",
       "wV4HSuPoY4BzgUdoNI5eWv6/N2lf9EeB7SOLBwoo1cwqwMFe4yTWJ42jrwocGsFf5r+Wlv/vTFr+\n",
       "/wawd2TR7BmkZtZ1ONhrlERv0nz0PUlzzi+YdyRdafn/tqRxdkjDM7d4tahZbXCw15jSgRfjSaE9\n",
       "hXSC0f/mv55rBGn5/yqk5f/XO9DNaouDvYZIbEQ6Zi6AHSK4d/5ruYaSVokOJo2lXxlZNDlU2sy6\n",
       "Pgd7DZBYETiLtD/6icAVEcwFUK6BpKGYLUk99Z0ji4+KqtXMOp6DvQsrDbtMIAX2/wGD5u2Prlyr\n",
       "k+aojwV+CBwQWXxQVK1mVj2tBruk3qRDFQaT/hd//4i4u8HrI0mLWJ4tPXV9RJxe+VKtIYkhwGRg\n",
       "MWBUBPcDKNcqpF77vsAFpOX/bxVWqJlVXTk99h8DN0XErpK6A0s10+aOiBhT2dKsORLLkma77EW6\n",
       "+XlxBHOVaznSBl6HAVcBgyOLV4ur1MyK0mKwS1oO2Dwi9gOIiNnw6UrFhk07oDZroMFmXecBtwCD\n",
       "I3hDuZZUzpGkUP8DMCyy+E9xlZpZ0Vrrsa8FvCHpUmAIcB9wVMQCW7UGsJmk6cBLwHER8WiHVFun\n",
       "JPqTZrv0A8ZF8A/l6qmcw0i99juBLSKLx4qs08w6B0ULU5glbQjcBWwWEdMknQe8GxGnNmizDDAn\n",
       "ImZKGgX8OCIGNvNZQRpCmGdqREyt0H9HTSrt7XI8aQHR2cCPmKS5wNdJ1/JJ4OTI4r7iqjSzSird\n",
       "txzZ4KksIto0KtJasH8GuCsi1io9HgGcEBFfa+E9zwHDImJGo+ejrcXVM4ktSTc/nwKOZJKeB3Yk\n",
       "zUV/GzgpsrijwBLNrArak50tDsVExKuSXpA0MCKeBLYhbSLV8Ev7AK9HREjamPTLYkZzn2etk1gJ\n",
       "+AHpWh/JJN0IfBn4DdCT1IP/k1eLmtnClDMr5kjgKkk9gWeA/SUdDBARFwK7AodKmg3MBMZ1VLG1\n",
       "rHRzdDzwPeBqYD0maTDwV2A14BTgushibmFFmlmX0OJQTEW/yEMxCyWxDnAhsDRwEJM0mzTk8kXS\n",
       "qtHLvPzfrD5VfCjGOpY0f2jlGOA0vr3Cn+j1VkYaejkL2D2y+LDIGs2s6+lWdAH1SmIT0vTRzRiZ\n",
       "bc8krUevt/4JPAYMiCx+7FA3s/Zwj73KJJYmDbPswarTMiZuvA7iJtK2DQMj841nM1s0DvYqktgW\n",
       "uJAlX7+TI9f9Fb3eOgu4Flg/snil4PLMrEY42KugtK3uj+jxwRbstf0trHnHjsCfgI0ii2dbebuZ\n",
       "WZs42DtQaQrj7nT7+Dy+fNIjbPbDHihWBLaKLB5p7f1mZu3hYO8gEqvR7ZPJDL10KNsdNZseH84F\n",
       "dowsphVdm5nVNgd7hUkIzRnPoN+dx+gjPmTp155DcWJkcXvRtZlZfXCwV5BEPwZddz1bnjaYFZ96\n",
       "iR6zjgN+7+X/ZlZNXnlaARJiw8lnsP6vj2OVh96jx6yj6DHrai//N7NF1Z7sdLAvIu2y1zYs8/JV\n",
       "rPzo8ryzxlmsNu30yOKTousys9rgLQWqSKcs3p8Za19O/zc246nRN/Px0rvF1b9/v+i6zMwc7G2k\n",
       "XKvy0dKnw+J78eT27/DE2C3j+S/9vei6zMzmcbCXSblWBL7NnB6H8MD+4u6jL+HttY6NYFbRtZmZ\n",
       "NeRgb4VyLQMcQ/ANnh71Kn+Y/C7vrLFPBJ6+aGadkoN9IZRrCeBQ4ATe+8w0Lr/tLf43aDpweARv\n",
       "F1yemdlCeVZMI8rVg3SS0akE9/F/1z3FY7uMB46O4KpCizOzuuNZMYtAuboBu5NOLHqBp79yGFfe\n",
       "chSwKbBxBM8VWqCZWZnqPtiVS8D2wBnAh8ChTIqlgIuAnwFnReBj6cysy6jrYFeuLYEzgeWAk/nh\n",
       "8zfzbr+zgTHAzhH8s9ACzczaoS6DXbmGkQJ9beBU4FomxdrA3cDTwNAI3iqwRDOzdqurYFeuQcBp\n",
       "pHHz04BLmBSzSTdLfwCcDPwiAm/aZWZdVl0Eu3KtCUwCRgNnA/tEFrNK549eCgwFRkbwcFE1mplV\n",
       "SqvBLqk36aDlwUAA+0fE3Y3anA+MAmYC4yPigQ6otc2U67OkXviewE+BAZHFOwASnwd+A9wJbBTB\n",
       "zMIKNTOroHJ67D8GboqIXSV1B5Zq+KKk0cDaETFA0ibAZGB45Ustn3KtAHwLmAhcBqwbWbwx/3Ux\n",
       "gdRzPzaCKwop0sysg7QY7JKWAzaPiP0AImI28E6jZmOAy0uv3yOpt6Q+EfFaRxTcEuVaGji69OcG\n",
       "YEhk8eL818VSpCmMGwNbRvBotWs0M+torfXY1wLekHQpMAS4DzgqIhoOW6wGvNDg8YtAX6BqwV5a\n",
       "/n8wcCJwG7BpZPHUAm3EIOA64F7S0MsH1arPzKyaWgv27sAGwBERMU3SecAJpCmCDTVe7trsrBJJ\n",
       "kxo8nBoRU8svtZnPy9Ud2A/IgOnAtpHF9Kbfy+6knvoJwC8968XMOitJI4GRi/IZrQX7i8CLETGt\n",
       "9Pg6Ujg29BLQr8HjvqXnmoiISe2osYnS8v9dSVMWXwHGRRZNFhNJ9CCNpY8Fto3g/kp8v5lZRyl1\n",
       "eKfOeywpa+tntBjsEfGqpBckDYyIJ4FtgEcaNZsCHAFcK2k48HZHja+Xlv+PIi3/n1363lubOyxa\n",
       "4rPA/wHvAhtGMKMjajIz62xa3d1R0hDSdMeewDPA/sAeABFxYanNT4HtgA+ACRHRpGe8qLs7Ktfm\n",
       "pNWiKwDfAX7XXKCn72IL4BrgAuCMCHyotJl1STV5mLVybUDqoa9LGku/KrKY0/x3IOBI0tz1fSO4\n",
       "eRFKNjMrXE1t26tc65LG0L9ECvaxkcXHC20vepF66F8ENo3g2aoUambWyXS6YFeuNUg98x2Ac4Dx\n",
       "kUWLUxMl+gG/BZ4CNvNURjOrZ50m2JWrD3ASsDep5z0gsmj1CDqJLYFrgXOBcz2V0czqXeHBrly9\n",
       "geOBQ4BfAetF1vqsmtJ4+uGkG6n7RPCXDi3UzKyLKCzYlWsp0o3OY0lTJodGFs+X9V6xOPBzYCPS\n",
       "0IvH083MSqoe7MrVEziINOzyD2DzyOLxst8vViXtA/MCKdTf75BCzcy6qKoGu3LtR9oX/TFg+8ja\n",
       "tr2vxHDS6tfJwJkeTzcza6raPfYDgX0ji7+39Y0S+wPfA/aP4A8Vr8zMrEZUO9i3WNhq0YUp7fdy\n",
       "LrAtsEX9mTwqAAAFaElEQVQEZQ/bmJnVo0698lRiRdIpRx8CX4+g1emPZma1pD3Z2a2jillUpaPr\n",
       "pgH/AnZwqJuZlafweezNkdiFtEjpGxFcU3Q9ZmZdSacK9tKio1NJO0huF8F9BZdkZtbldJpgLy06\n",
       "uhgYCGwSwasFl2Rm1iV1ijF2iZWAW4FewFYOdTOz9is82CUGAncBdwK7RzCzlbeYmVkLCg320s6M\n",
       "fwe+H8EJPunIzGzRFbcJmNiHtPBozwj+WlQdZma1pvqbgKWZL5OAfYGRETxa7RrMzGpZdTcBSzNf\n",
       "LgEGAMMjaHXfdTMza5tqj7H/BViCNPPFoW5m1gGqHex34ZkvZmYdqqxNwCT9B3gXmAN8EhEbN3p9\n",
       "JHAjzD/J6PqIOL1RmzZvZGNmVu86chOwAEZGxNDGod7AHaXXhzYOdVtQ6Reh4WvRkK/Fp3wtFk1b\n",
       "hmJa+43h3nj5RhZdQCcysugCOpGRRRfQiYwsuoCurC099lsl3Stp4kJe30zSdEk3SVqvciWamVlb\n",
       "lDvd8UsR8YqklYG/SHo8YoHj7e4H+kXETEmjgN+RNvMyM7Mqa/MJSpIy4P2IOLeFNs8BwyJiRoPn\n",
       "fPC0mVk7tPXmaas9dklLAotFxHuSlgK+CuSN2vQBXo+IkLQx6RfGjIZtPCPGzKw6yhmK6QP8VtK8\n",
       "9ldFxC2SDgaIiAuBXYFDJc0GZgLjOqheMzNrRdUOszYzs+qo+MpTSdtJelzSU5K+vZA255deny5p\n",
       "aKVr6CxauxaS9ipdg39LulPSF4qos6OV83ei1G4jSbMl7VzN+qqpzJ+PkZIekPSwpKlVLrFqyvj5\n",
       "WEnSnyU9WLoW4wsosyok/VLSa5IeaqFN+bkZERX7AywGPA2sCfQAHgQGNWozGrip9O+bAHdXsobO\n",
       "8qfMa7EpsFzp37erxWtRznVo0O424A/ALkXXXeDfid7AI0Df0uOViq67wGsxCThr3nUA3gS6F117\n",
       "B12PzYGhwEMLeb1NuVnpHvvGwNMR8Z+I+AS4FhjbqM0Y4HKAiLgH6F26+VprWr0WEXFXRLxTengP\n",
       "0LfKNVZDOX8nAI4ErgPeqGZxVVbOtfg6aUuOFwEi4n9VrrFayrkWrwDLlv59WeDNiJhdxRqrJtL0\n",
       "8bdaaNKm3Kx0sK8GvNDg8Yul51prU4uBVs61aOgA4KYOragYrV4HSauRfqgnl56q1Rs/5fydGACs\n",
       "IOn20oLAfapWXXWVcy0uAgZLehmYDhxVpdo6ozblZqX3Yy/3B7Lx1Mda/EEu+79J0lbA/sCXOq6c\n",
       "wpRzHc4DToiIUJp+VatTY8u5Fj2ADYAvA0sCd0m6OyKe6tDKqq+ca3ES8GBEjJTUn7Q4ckhEvNfB\n",
       "tXVWZedmpYP9JaBfg8f9SL9ZWmrTt/RcrSnnWlC6YXoRsF1EtPS/Yl1VOddhGHBtaUrtSsAoSZ9E\n",
       "xJTqlFg15VyLF4D/RcQsYJakvwFDgFoL9nKuxWbAGQAR8Uxp4eM6wL1VqbBzaVNuVnoo5l5ggKQ1\n",
       "JfUE9gAa/3BOIR2Lh6ThwNsRUYuHbrR6LSStDtwA7B0RTxdQYzW0eh0i4nMRsVZErEUaZz+0BkMd\n",
       "yvv5uBEYIWmx0uLATaAmj48s51o8DmwD8xdBrsOnW4PXmzblZkV77BExW9IRwM2ku96XRMRjDRcz\n",
       "RcRNkkZLehr4AJhQyRo6i3KuBXAqsDwwudRbbbLXfVdX5nWoC2X+fDwu6c/Av4G5wEURUXPBXubf\n",
       "izOBSyVNJ3VCvxWNVrTXCknXAFsCK0l6AchIw3Ltyk0vUDIzqzHVPhrPzMw6mIPdzKzGONjNzGqM\n",
       "g93MrMY42M3MaoyD3cysxjjYzcxqjIPdzKzG/D9Gf9er9HmuLwAAAABJRU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fba7350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs = np.linspace(0,1,len(xval))\n",
    "def func(x, p1, p2, p3, p4):\n",
    "    return p1 + p2*x + p3*np.exp(p4*x)\n",
    "fitParams, fitCov = curve_fit(func, xs, np.log(xval))\n",
    "pl.plot(xs, np.log(xval))\n",
    "pl.plot(xs, func(xs, fitParams[0], fitParams[1], fitParams[2], fitParams[3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "first fit [0.55765009832017132, 8.6753672190427995, -5.4038968232544793, -0.042639640230075333]\n",
      "second fit [[Variables]]\n",
      "    p0:   0.55765045 +/- 0.221877 (39.79%) (init= 0.04905465)\n",
      "    p1:   8.67536775 +/- 0.346239 (3.99%) (init= 7.308738)\n",
      "    p2:  -5.40389635 +/- 0.280458 (5.19%) (init=-6.004465)\n",
      "    p3:  -0.04263955 +/- 0.051870 (121.65%) (init= 0.1347386)\n",
      "[[Correlations]] (unreported correlations are <  0.100)\n",
      "    C(p2, p3)                    =  0.993 \n",
      "    C(p0, p2)                    =  0.993 \n",
      "    C(p0, p1)                    =  0.988 \n",
      "    C(p0, p3)                    =  0.974 \n",
      "    C(p1, p2)                    =  0.967 \n",
      "    C(p1, p3)                    =  0.937 \n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEDCAYAAAA4FgP0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAGr5JREFUeJzt3X2QZXV95/H3d3qGcQZsJtNsYnyIULUaoYotla24tSuR\n",
       "3YiDceNo1MSJ61ImYUWLtO5gdDNYWbZWRjQ6sg2lqVUyBZQB1jJksFgyTBKxqFiJ+ECJD+BSpRZR\n",
       "osAA7TAEZnq++8c5d/rOndvT9/Y9p+/Deb+qpug+9/a9vz5F3V//nr6fyEwkSc21ZtgNkCQNlx2B\n",
       "JDWcHYEkNZwdgSQ1nB2BJDWcHYEkNZwdgSQ1nB2BJDVcrR1BRJwVETdHxCcj4k11vpckaWXqHhFc\n",
       "AFydme8G/nPN7yVJWoG+O4KI+LOI+ElE3Ntx/YKIuC8i/l9EfKC8fAPw1oj4KDBTQXslSRWLfmsN\n",
       "RcS5wAHg+sw8u7w2BdwPvBr4EXA3sC0zv9v2+Ocz8w0Vtl2SVIG1/f5AZt4VEad3XP4V4IHM/AFA\n",
       "RNwEbI2Ig8AO4GTgowO1VJJUi747giU8D3iw7ft/BF6RmT8E3lnRe0iSalBVR7DiWtYRYR1sSVqB\n",
       "zIwqXqeqjuBHwAvavn8BxaigJ1X9MuMuIi7PzMuH3Y5R4L1Y5L1Y5L1YVOUf0VVtH/0q8KKIOD0i\n",
       "TgJ+G7i1oteWJNVoJdtHbwS+DLw4Ih6MiHdk5mHgEmAv8B3g5taOIUnSaFvJrqFtS1y/Hbh94BY1\n",
       "253DbsAIuXPYDRghdw67ASPkzmE3YBL1fY6g8gZEpGsEktSfKj87LTonSQ1nRyBJDWdHIEkNZ0cg\n",
       "SQ1nRyBJDVfVyeKuIuL5wBzwGPC9zPxIne8nSepf3SOCsynKT/8e8LKa30uStAJ1B9N8GfgvEfE3\n",
       "wF9V0F5JUsVqDaYBtgBfKzMMPpeZb+nyeh4ok6Q+DfVAWWbeRTHn3+5oME1mHgJuArYCfwu8JyI+\n",
       "BXx/0MZKkqpXdzDNN4E3V/QekqQaDD2YBiAi7gR+UP67MzPvHLhFkjRBIuI84Dzg9PJfZUYlmOa8\n",
       "itohSROp/AP5ztb3BtNIkipjMI0kNZx5BJI0hswjkCRVxo5Akhpu7DuCiNhxaqx95NRY+0hE7Bh2\n",
       "eyRp3NRafbRuEbFjGq6YYwGAWbgiIsjMnUNumiSNjbFeLD411j4yx8LMheX31wGzTD36RB4+rboW\n",
       "StLoqXKxeKxHBC17gf8N/Bg4xJENQ26OJI2VWtcIIuKVEfGpiPh0RPxd1a8/z8KudwP/CXg9cDGw\n",
       "jtzoWoEk9W5VpoYiYivw85n56S6PDTS8mY6pJ6/myEanhyQ1yVDPEfQZTNPyO8CfD9LQJdtDPNX6\n",
       "ei/wp8AaFjY7KpCk3qxkamg3cEH7hTKY5pry+lnAtog4s3zsl4AnMvPJAdva1TwLu2aB91FMEV0M\n",
       "XAUxXewgsjOQpGXUHUwD8LvAnw3UyhO3Z+c8XHYt5MeAC8t/c8A0U9vrel9JmhRVLRZ3C6Z5HkBm\n",
       "Xp6Zf1/R+3SVmTuPMLUfiumhN1FMEbmDSJKWNzHBNOUOois2Ah8rr82WO4g8YCZp3BlM09tr7JyO\n",
       "qcs+Vu4g2kuxWHE/XB4RX8vMvYO+hyQNi8E0PWrtINpLsU5wMfAJWDcNeyJiy2q3R5LGwUQF07R2\n",
       "EF0OfIRjFo7XbwKnhySpi76nhjJz2xLXbwduH7hFA8jMnRHBffAh4JiDFgusPTsitjhFJEnHGuui\n",
       "c0u/5ilfm+aZl89xiHuBawkO80sc4OGn4eBWOwNJ467Kz84J7QhiC2zc8yxYfxJPMVduapplHfOc\n",
       "9PXMA+dU+X6StNqMqlxG8Rf/wa3reObQHMlzKFauz+IQGzj4kmG3T5JGyUR2BFB0BlMcvvdeigVj\n",
       "q5NKUncTOTXU9tpbNsFfXUXRGYDVSSVNBqeGepSZexdY87PW91YnlaTjTXRHAPAzjlxpdVJJWlqt\n",
       "U0MRERR7+p8NfDUzr+/ynNqmhtreY8cm+NBVEE4RSZoE4zQ19AaKKqTP0Eftoaq1VycFp4gkqV3d\n",
       "CWUvBv4uM98HvKuC9q6YATaS1F3dCWX/CDxePu3IAO0cWLcAm+dQNHYTfMjOQFJT1Z1Q9hfAloiY\n",
       "o6186rB0Bti0KpQ6MpDUZFXlEXRLKHtFZj4F/P5yP1xFME2vyimiK85isUJpy2wRbWmVUkkjZxyC\n",
       "aQbaelRFME0f73VchdLOxWMTzSSNmjqDaUYioWy1tTqDWbjiXoqtpGW8ZcwWU0TYGUhqiolKKOvH\n",
       "iRaPf66ItzTRTFIjTFRCWb+WWjw23lJSk0x00bke33/HdLl4fDHHFqd7L3z9sUyzCySNnHE6WTzy\n",
       "WlNE97UteLcWj4GXuaVU0qRr/IhgsR1FvOVFHGpfPGYWmIfLXDyWNEqMqqynHVtg455NPL3+KhaO\n",
       "ThG9D7gW8ghT++dZ2GWHIGkUODVUg1a85QJxqHVtL8VawVUQcyzMePpY0iRyRNChtXg8R7FO0LmA\n",
       "bOlqSaPAEUGNWovHs0w92m0B2dLVkiZN3cE05wH/E/gWcFNmfqnLc0ZqRNCuNTq4CFxAljRSxmlE\n",
       "cAT4GbCeES45sRRLV0tqgrqDae7KzF8H/hvwPypo76qzdLWkSVdrME0uzjs9TjEqGEutdLPLWSxd\n",
       "fSEwB0wXpaslaWz1XX00M++KiNM7Lh8NpgGIiJuArRHxEmALsAm4eqCWDlG30tUA91IsHp8aax/x\n",
       "jIGkcVV3MM2VwC3L/fBqBtOsVHvpaig6gU8DcxCwMGP5akl1MphmRCx2BlPb17CweQ6idcbgXuBa\n",
       "+NCpsXa7owNJVaszmKaqXUNjFUwziMzc+UQePq21gAyeQJY03hobTDOo1gLydRSLyG4vlTSuGh1M\n",
       "M4huJ5DdXippHFlrqAInCrexNpGkOozTyeJG6BZuA8duL3VkIGlUVbVrqPHcXippXDk1VLFimqjY\n",
       "XnpV2/ZSA24kVcmpoRHm9lJJ48aOoCZLbS+9ELiIYnupaweSRoFrBDXpPIVMWaOofXTg2oGkUVD7\n",
       "GkFEnExxLPryzLyty+MTtUbQjfGXkqo2bmsE7wduXoX3GVlLxV+CW0wlDV+twTQRcT7FSeOHq2nu\n",
       "+GotIj8OH2ytHbyPYoupi8iShqnvqaGIOBc4AFyfmWeX16aA+4FXUxSguxvYBrwNOJmi/M5TwBuz\n",
       "4w2bMDXUyS2mkgZV5WdnrcE0mfnB8vsLgYc7O4GmKj/gd54aax+BhRlwEVnS8NQaTNP6JjOvO9EP\n",
       "j0MwTR3KLaZXQLGI3NpiCuYbSDqWwTQTyi2mknplMM0E67aIfDnHHkCbAzYQl50aax9xd5GkqlU1\n",
       "IjgaTAP8mCKYZltFr90IS40OoJgmWkdunGNhI4AjBElVWsmuoRuBVwEzwE+BP87M3RHxWuAqYAq4\n",
       "NjM/3OPrNW7X0HLaD6ABvJfixnoITVLLsHcNdf1LPzNvB24fuEU6ZnQAkOQGOLKx9Xj7ITQXkiUN\n",
       "yjLUY6B9hNCWcwDALDAPl9kZSM1S5WenHcGYWOoQ2nXAu4iD61jzFBRbUu0UpMk3brWGVIFuOQdw\n",
       "zELyjGUqJK2EI4Ix08tCsiMEafINdbFYw9XLQrJbTSX1wxHBmHOrqdRMYzMiiIiXAO+hOHOwNzOv\n",
       "rfP9mmi5EQLAIY5sKArcOVUk6XirMiKIiDXATZn5W10ec0RQoc4Rwrspenu3m0qTZai7hvoJpimv\n",
       "/wZwG3DT4M3VctrT0GaZenSKNQfnsG6RpKWtZPvobuCC9gtlMM015fWzgG0RcSZAZn4hM1/L4rS1\n",
       "atbaavpEHj4tiKfaH3O7qaROtQbTRMTPA78JPAv44kAt1Yq0Zx4AXEsxKmjvld9VjBC2t57vtJHU\n",
       "LLUG02Tml4AvVfQeWgG3m0pazkgE0zQ1oWy1tKIxoVhMXm6EUHYadgTSCBmHhLKBgmmamlA2DG43\n",
       "lcZTnQllK9o+Wq4RfCEzzy6/XwvcD/waRTDNV4BtmfndHl7L7aND5HZTaTwN9UBZezBNRDzIYjDN\n",
       "JRRxu61gmmU7AQ1f5whhitwwx5GNnVNFEcF0+RxHCdJkscSEjnFqrH1kjoWZ44vY5UZHCdLoMI9A\n",
       "temcKpoFkjUHr24bJVjhVBq+sak1pPHTOVU0z8KuaWI74JZTaUI5ItCyrHAqjR4TyrSqOusXLbDm\n",
       "YLfnRcQOaxhJ48cRgfrWbR1hHu6Yhte4oCytDtcINFTd1xGmts+xcNwJZbedSqOv1hFBRGwFXgdM\n",
       "U5wt2NflOY4IJoDbTqXVNXbbRyNiE/CxzPz9Lo/ZEUyAXreduqgsVWOsgmlKH6TIK9CE6lxQnofL\n",
       "OrMQWlxUlkZL3yOCiDgXOABc31ZraIqi1tCrKQrQ3Q1sA+4DrgTuyMy/WeL1HBFMKBeVpfoMdUSQ\n",
       "mXcBj3VcPhpMk5mHKGIptwKXUBSie3NEvHPQxmq8dBslTDN1Tmd0Zmsx2ZGCNBx1B9P8AXB1Re+h\n",
       "MdSehQDQSkLrtDh6WAA8rSytJoNptKo6ozNncfup1AuDaTQxup1ByMyd3UYKhziywVGCVDCYRhPP\n",
       "7adSf4a9ffRG4MvAiyPiwYh4R2YeplgY3gt8B7jZYBr1o5/tp+DCslQlaw1pZC2x/fQygG7XnTJS\n",
       "k4zdyeITNsCOQCdQdAbHrSccV86iGEUsvG0zXAqwHz6emXuH1W6pbhadU2N0bj9dyiFy3TTs2QXr\n",
       "AWbhVyNiq52BtDw7Ao2dbltQ13HkwBxMt20/Xf/eogOxI5CWYTCNxk63hWWYWt/5vCNMvRBcWJaW\n",
       "4xqBJkLEKV+b5pmXz3EIgFnWMc9JX4cnP+/CsiaRi8VSh4jYAhv3nMIZ6wEO8P2n4eDWaaY+221h\n",
       "2bMIGndjk1kcEWdExGci4nN1vo9ULAof3HqAh/Yd4KF9cLCnhWKnjaTVC6b5XGa+ZYnHHBGoNkud\n",
       "RShLXSz52HBaK/Vu2CeLVxJMIw1Ft4Xl1gd9UejOktjSSqaGdgMXtF8og2muKa+fBWyLiDMHb540\n",
       "uMzc+UQePu2JPHxaL3/tt5XEnpljYWa6KHZnZ6CJVWswTURsjog/BV7qKEGjpjyPwHW0FpHbS2J3\n",
       "HylIk6iqxeJuwTTPy8z9mXlxZr4oMz9S0XtJlTjRtNGJOG2kSWMwjRqtWwmLpcJzwCQ1DY/BNNIq\n",
       "Wio8B1oLzMcnqdFDPSRpEHUG01Q1NfRV4EURcXpEnAT8NnBrRa8trbp+F5jBKSONr77PEZTBNK8C\n",
       "ZoCfAn+cmbsj4rXAVcAUcG1mfrjH1/McgcaGGQkaFZaYkIaon4wES1moLmNTYkKaRE4badKYRyBV\n",
       "wJ1GGmdODUkV6TZlBNDLtNFSPystxahKaQT1GqvZyRGDhs2OQKrZiaaNwLMJGj47AqlmJzqg1osi\n",
       "dGfzpcV3+z/eS86C1A/XCKQhWyYzoWvymp2BxuYcQUScDHwSeJqihtCfd3mOHYEab6nF4qWymDMP\n",
       "nDPE5moEjFNH8HZgf2beFhE3ZeZbuzzHjkBagjuOtJRxSihrL0+9MGBbpcZZw8IPT3TNEB1VYSW1\n",
       "hs4FDgDXZ+bZ5bUp4H7g1RSVSO8GtgHnAI+VI4IbM3Nbl9dzRCAtISK2TMOeOVgPMAtPz8PRNQJL\n",
       "WzTXUEcE/SSUAX8BvCkiPonVSKW+Zebeedi6HfZth33tnUCvImJLxMwdxb/YUldbNb6q2j7aLaHs\n",
       "FZl5EPjd5X7YYBppaeUHf9cP/+XOKCzuOvrFctfRP/9qRLjraAyNQzDNQCvOBtNIK7P8GYWTd07z\n",
       "zPo5vg3ALOvWz3PyTpboWDS66gymGYmEMkkrd6LSFtP88wuPPbV8iFmOvHC12qbxYEKZNMGW23UE\n",
       "lsjWCkYE7QllEfEgiwlll1AMN1sJZd+ttqmS+vU47JiFPRy76+joh70F7wSWmJAmXkRs2QyXAuyH\n",
       "Y2oVuf10fFmGWlLPTrTrqFeeXp5sdgRSgy23/RScPmoCp4akhlvur32nj0aT4fWSKpOZO5/Iw6c9\n",
       "kYdPG+SvfE8wjy+nhiSdUI/TR1tg+hbYtaF81isj4o2eYB4PTg1JWtZy00cRM3fArvOhfQJp+77M\n",
       "R1+z2m1tirGZGoqIMyLiMxHxuTrfR1K9qpw+mom4YybC6aMRsiojgoj4XGa+ZYnHHBFIY66XSM3l\n",
       "SmqrP6s+IugzjEZSA01zkGv4NtfwbaY5eNzjm2DnHKy/kGICaQ7Wb1qiRpJWV69TQ7uBC9ovlGE0\n",
       "15TXzwK2RcSZEfH2iPhERDy32qZKGlWb4dLOD/nWaeaWI0wdV+yu2zWnj1ZfT7uGMvOuiDi94/LR\n",
       "MBqAiLiJYph3JXBDeW0zRY//0oj4QGZ+pKJ2Sxoz8zzrh7M8MwOHAJhlHfOc1FkAb8s03LILNhTP\n",
       "wd1Hq2CQ7aNdw2jan5CZ+4GLl3shg2mk8bYfPj4Lr2TxA/ypefj4sc96csc8G/dcwovb1hGePKba\n",
       "6Wa4dBdsWCybzYbtxcii8R3BqAbTVLbKbDCNNN4yc29EvLH80Ga+o7hd23O2HuChcsro4HHPOQIz\n",
       "na/d7VoTjWowjWE0ko7qpbjdcs95nJOZ5Rk6po+qbKa6GOQcgWE0kiq2/tF5ZrmE87mE85lnFlj/\n",
       "aLdnGqhTnV63j94IfBl4cUQ8GBHvyMzDQCuM5jvAzYbRSBrM/o/Dp586wNs4wNuATz9VXDtWW0XU\n",
       "mTkWZqaLiqh2BitkiQlJI6XYMrq5XEfYf9w6AlgRFQymkTTBqgjSUX8sQy1p7JQVUbmO1mjg+Iqo\n",
       "La4lLM+pIUljqZf4zMW1hELRYXDZJKSrVfnZaUcgaWJN8lrC2JShlqRx0tQ6Ry4WS5pYvaSrtTS5\n",
       "zlHtU0MRsRV4HTANXJuZ+zoed2pIUm16WUsAmIm4Yxecf0zGGux7NHMkU9bGavtoZu4B9kTEJuBj\n",
       "wL5lfkSSKlN+8C+7ONzkOkc9rxFUEE7zQYr8AkkaOUWdo3VtW1LX8TgnD7tZq6KfEcFu4Grg+taF\n",
       "tnCaV1MUobs7Im4F/jXwcuBPgIeAK4HbM/OeitotSRVb/+g8F3MJ3wTgAP8K2N21ztGk6XlEkJl3\n",
       "AY91XD4aTpOZh4BWOM0NmflfM/PHwB8Avwa8OSLeWVXDJalavdU5aomILREzdxT/xnuH0aBrBL2E\n",
       "08wBc5yAwTSShq2VqQDbyzpH813rHEGrHtL0LbBrQ3FltvYdRqMaTAMVhdMYTCNpFPRe52jzpfCO\n",
       "DafwWQAOcNEG2F1rktqoBtOA4TSSGunpmWnmmDsaoHMn85w0tjuMBu0IjobTAD+mCKfZNuBrStJI\n",
       "28STXAUsZisf4r1lpzCO+tk+ajiNJAFr4LjdRN2ujQuLzklSn1rlKOYWy1E8NQ+rWo7C6qOSNGQR\n",
       "sWUzXAqwH5bcYVTj+9sRSFKTWYZaklQZOwJJajg7AklqODsCSWq4WjuCiHhJRHwqIv5PRPxene8l\n",
       "SVqZVdk1FBFrgJsy87e6POauIUnq06rvGhoklCYifgO4jaJEtSRpxPQ0IoiIc4EDwPWZeXZ5bQq4\n",
       "n7ZQGoo6Q0dDaco8gtZr7MnMrV1e2xGBJPVp1TOLM/OusrBcu6OhNGWjWqE0VwI3lNdeBfwm8Czg\n",
       "i1U0WJJUrUGqj/YSSvMl4EsDvIckqWaDdASVrTKbUCZJJzaqCWWVhdKYUCZJJ1ZnQtkg5wiOhtJE\n",
       "xEkUoTS3VtMsSdJq6XX7qKE0kjShLEMtSWPIMtSSpMrYEUhSw9kRSFLD2RFIUsPZEUhSw9kRSFLD\n",
       "1d4RRMTJEXF3RLyu7veSJPVvNUYE7wduXoX3kSStQK3BNBFxPsWp44erae5kK4tKCe9FO+/FIu9F\n",
       "PXodEewGLmi/UAbTXFNePwvYFhFnRsTbI+ITEfFc4FXAvwF+B7goIjxBfGLnDbsBI+S8YTdghJw3\n",
       "7AaMkPOG3YBJVGswDfDB8rELgYdz2PUsJEnHqTWYpiUzrxvgfSRJNRqVYBpHCqWI+O/DbsOo8F4s\n",
       "8l4s8l5Ub+jBNFYelaThMphGkhrOYBpJarieOoLM3JaZz83M9Zn5gszcXV6/PTN/OTP/ZWZ+uN83\n",
       "X+4cwiTodgYjIjZHxL6I+F5E3BERm9oe+6PyftwXEa9pu35ORNxbPva/Vvv3GFREvCAivhgR346I\n",
       "b0XEbHm9iffiWRHxDxFxT0R8JyI+XF5v3L1oiYipiPhGRHyh/L6R9yIifhAR3yzvxVfKa/Xfi8wc\n",
       "yj9gCngAOB1YB9wDnDms9tT4e54LvAy4t+3aR4H3l19/ALiy/Pqs8j6sK+/LAyymyH0F+JXy6/8L\n",
       "XDDs363P+/Ac4KXl16cA9wNnNvFelO3eWP53LfD3wCubei/Ktm8HPgvcWn7fyHsBfB/Y3HGt9nsx\n",
       "zKJzR88hZOYh4CZg6xDbU4vMvAt4rOPy64HWltrrgDeUX28FbszMQ1mcz3gAeEVE/CLw7Mz8Svm8\n",
       "69t+Zixk5j9l5j3l1weA71JsQW7cvQDIzIPllydR/FH0GA29FxHxfODXgc8Arc0jjbwXpc4NNLXf\n",
       "i2F2BN3OITxvSG1Zbb+QmT8pv/4J8Avl18/l2J1XrXvSef1HjPG9Kg8nvgz4Bxp6LyJiTUTcQ/E7\n",
       "fzEzv01D7wXwCeAPgSNt15p6LxL464j4akRcVF6r/V4Msn10UJ4dADIzm3SOIiJOAT4PvCczf9Ze\n",
       "daRJ9yIzjwAvjYhTgb0R8e87Hm/EvYiI/wj8NDO/EUvUEWrKvSj9u8x8KCL+BbAvIu5rf7CuezHM\n",
       "EUEl5xDG1E8i4jkA5TDup+X1znvyfIp78qPy6/brP1qFdlYqItZRdAI3ZOZflpcbeS9aMvMJ4Dbg\n",
       "HJp5L/4t8PqI+D5wI/AfIuIGmnkvyMyHyv8+DNxCMYVe+70YZkfQ5HMItwIXll9fCPxl2/W3RsRJ\n",
       "EXEG8CLgK5n5T8B8RLwiij+h3972M2OhbPe1wHcy86q2h5p4L05r7fyIiA3A+cA3aOC9yMwdWexE\n",
       "PAN4K/C3mfl2GngvImJjRDy7/Ppk4DXAvazGvRjyCvlrKXaPPAD80TDbUuPveCPwY+AZijWRdwCb\n",
       "gb8GvgfcAWxqe/6O8n7cB2xpu35O+T/FA8DcsH+vFdyHV1LMAd9D8aH3DYrKtU28F2cDXy/vxTeB\n",
       "PyyvN+5edNyXV7G4a6hx9wI4o/x/4h7gW63PxNW4F62tRpKkhjKzWJIazo5AkhrOjkCSGs6OQJIa\n",
       "zo5AkhrOjkCSGs6OQJIazo5Akhru/wP1yNc/0NqdJgAAAABJRU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bca37d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#all fitting models. First has x bin width included, second is standard for plotting purposes, \n",
    "#third is log (gives better fit). Residual is for lmfit\n",
    "\n",
    "def model1_CDF2009(t, params): \n",
    "    global xerr\n",
    "    p0, p1, p2, p3 = params\n",
    "    sqrts = 8000.\n",
    "    return (p0 * (1.-t/sqrts)**p1 * (t/sqrts)**(p2+p3*np.log(t/sqrts)))*(xerr)\n",
    "\n",
    "def model2_CDF2009(t, p0, p1, p2, p3): \n",
    "    #p0, p1, p2, p3 = params\n",
    "    sqrts = 8000.\n",
    "    return (p0 * (1.-t/sqrts)**p1 * (t/sqrts)**(p2+p3*np.log(t/sqrts)))\n",
    "\n",
    "def model3_CDF2009(t, p0, p1, p2, p3):\n",
    "    #p0, p1, p2, p3 = params\n",
    "    sqrts = 8000.\n",
    "    return np.log(p0 * (1.-t/sqrts)**p1 * (t/sqrts)**(p2+p3*np.log(t/sqrts)))\n",
    "\n",
    "def residual_CDF2009(params, t,data=None): #gaussian distribution\n",
    "    #p0, p1, p2, p3 = params\n",
    "    vals = params.valuesdict()\n",
    "    p0 =  vals['p0']\n",
    "    p1 =  vals['p1']\n",
    "    p2 = vals['p2']\n",
    "    p3 = vals['p3']\n",
    "    sqrts = 8000.\n",
    "    return (np.log(p0 * (1.-t/sqrts)**p1 * (t/sqrts)**(p2+p3*np.log(t/sqrts)))-np.log(data))\n",
    "\n",
    "\n",
    "\n",
    "#first fit version using curve_fit\n",
    "fitParams, fitCovariances = curve_fit(model3_CDF2009, xval, np.log(newy), p0=(0.0490546530927,7.3087377889,-6.00446454843,0.134738604245), maxfev=100000)#, diag=(1./xval.mean(),1./yval.mean()))\n",
    "print(\"first fit\", list(fitParams))\n",
    "fitps = list(fitParams)\n",
    "\n",
    "#second fit version using lmfit\n",
    "fit_params = Parameters()\n",
    "fit_params.add('p0', value=0.0490546530927)\n",
    "fit_params.add('p1', value=7.3087377889)\n",
    "fit_params.add('p2', value=-6.00446454843)\n",
    "fit_params.add('p3', value=0.1347386042)\n",
    "\n",
    "out = minimize(residual_CDF2009, fit_params, args=(xval,), kws={'data':newy})\n",
    "\n",
    "print(\"second fit\", fit_report(fit_params))\n",
    "\n",
    "#compare first fit results with data\n",
    "pl.scatter(xval, model2_CDF2009(xval,  fitps[0], fitps[1], fitps[2], fitps[3])) \n",
    "pl.scatter(xval, newy, c='r')\n",
    "pl.ylim([10e-5, 10e8])\n",
    "pl.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEDCAYAAAA1CHOzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+QHPV55/H3s7vSspJYw84lIca6iCtDgjiqOLkKfDlx\n",
       "1h0WK58TFOLksGL7VDh24LC8TiQ7nLHLyBUHTAot6zXBxEISwkUEoYgCnK3sKrHlwnbZGLCKNSBA\n",
       "PnAJ4UCshQxCsiy0z/3RPTu9o5nd+dE90z3zeVVNeadntqfVpua73+d5vt/H3B0REZGuVl+AiIik\n",
       "gwYEEREBNCCIiEhIA4KIiAAaEEREJKQBQUREAA0IIiIS0oAgIiJAAgOCmZ1lZneY2X2RYwvNbLuZ\n",
       "fdXM/ijuzxQRkcbFPiC4+/Pu/pGSw78P/J27/wlwWdyfKSIijatqQDCzrWb2splNlBxfZWb7zOw5\n",
       "M7t2llOcCRwIfz5R57WKiEiCqp0hbANWRQ+YWTdwa3h8KbDGzM6t8PsvAotr/EwREWmiqr6c3f1h\n",
       "4NWSwxcC+939BXc/DtwDrDazATO7HbggMmv4e+B9ZnYb8GBM1y4iIjHqaeB3o2EgCGYBF7n7JHB1\n",
       "9I3ufgT4cAOfJSIiCWtkQIh132wz0z7cIiJ1cHeL4zyNDAgHKeYFCH9+sZGLiesflXVmttHdN7b6\n",
       "OtJA96JI96JI96Iozj+mG0nwPgqcbWZLzGw+cAUN5gfMbKOZrWjkHCIincDMVpjZxjjPWW3Z6Q7g\n",
       "e8A5ZnbAzK509zeBdcAY8BRwr7s/3cjFuPtGd9/TyDlERDqBu++Je5ZUVcjI3ddUOL4L2BXXxYSj\n",
       "3R4NCuxp9QWkyJ5WX0CK7Gn1BaTInlZfQKuF0ZQVsZ4zLT2VzcyVQxARqU2c351aJCYiIoAGBBER\n",
       "CaVqQGikysjMBs1y48HDBmO+NBGRVEmiyqgtcgjBANC/E0b7giNDRyF/ubuPxXmNIiJpE2cOoZGF\n",
       "aSkysAGG+2Bt4UAfrN9AUBIrIiJVaJMBoWAM+CrwEnAs1+KLERHJlDYKGS14ABb0ws3h0aFjkF+t\n",
       "sJGItDOVnZYIvvTnPxkMBmvDx2hvEEoSEZFqNGVAMLOlZnavmd1mZu+b5X0N7GXUdajOyxMRyZzM\n",
       "VhmZ2XrgEXf/jpk94O6ry7ynoWmPKo1EpBPFGTKqe0Aws63Ae4FX3P38yPFVwAjQDdzh7jeZ2a8A\n",
       "1wNHgN929+VlztfwPyoYFAY2BAnleQSzhslNGhREpF2lZUC4GDgM3FUYEMI+y88A7ybol/BDYE1h\n",
       "F9Tw9fvd/ffKnC+Wf5RmCiLSSVKxDsHdHzazJSWHp/ssA5hZoc/yEeA6YCHwV/V+ZnW0JkFEpB5x\n",
       "r0Oo1Gf5p8BVMX+WiIjEKO4BoaEMdUnGvM6+CJObYGg5TPTBd4F9U/BaHecREUmfJPogTJ+7kSqj\n",
       "MGT0UCSH8E5go7uvCp9/Gphy95uqOFd8iRGz66D/L2A0LKtVHkFE2lMqcggVTPdZJtg/4gqgbLe1\n",
       "cuLrmDawAoa7lEcQkXaVxEyh7oVpSfRZbqSnspkN5szGc2bj2sdIRNpdEj2VU7WXEfB56pghmNlg\n",
       "P+wchT6AITiWZwFwWy9MAFumgL3w2nUKG4lIO4jMEK5v+TqEuDUSB8uZjQ/DyjMo7nX6BOw/wsI8\n",
       "dF+gXIKItCttblfGBEHG4DLgaqAH/kPwymhXZMO7Pm14JyJSXlv0Q5iETVvgkhHoWls83DXEL34j\n",
       "37rLEhHJlFQNCPVWGbn72Olme4Fl0eNdnPgpDC0A+iK5hJyZDSpsJCJZlkSVUVvkEMLfn04sTwBb\n",
       "YArY+xrcDwvfp1yCiLSjVGxuF7e4djtdCDd0wwWjYX5kCI7mWfg0/PWy4rqE7cD63e6HLm30ukVE\n",
       "WklJ5QrcfawXDo2GuYS1wCj09fOL32j1tYmIpF2qcghxGqNYgvomJ14Ncgna30hEpJK2ChmF5xlc\n",
       "AA8sgN6bw2PBQjV2QP//Uh5BRNpJ24aMGuupHHD3sfnw5M0EIaMzgKXQezp8AD6qNQki0haS6Knc\n",
       "lAHBzN5mZn9vZlvM7NpK72tkL6OoLjgEQdhoLcFCtVtgXj+jaH87EWkHSexl1KwcwvkErTPvDruo\n",
       "JWoSNg3B8qXQdxORPU85zjo2cJjdaE2CiMhMjex2utXMXjaziZLjq8xsn5k9F5kNfA/4EzP7Z+Af\n",
       "G7jeqrj7WB4ufzacKcy4Pp5+HTZPwUgXjCyD/p1BH2YRkc5Wd1LZzC4GDgN3RRrkdAPPAO8GDgI/\n",
       "JOiHMAg8FvZhvs/d/7DM+WJLjETOedJitRN0vfE6f3YqFFLOWpMgItmViqSyuz8MvFpy+EJgv7u/\n",
       "4O7HgXuA1cA3gU+Y2VeA5+v9zDqucSwPl38MHt8MUyPQ9WWmTi3mEsaA2wGWaZYgIp0u7hzCmcCB\n",
       "yPMXgYvc/QngD+b65Xh6Ks/k7mM5sw3DMza+O846ruIwbxDOFHIwtNPMVIYqIqmWZE/luAeEhhY1\n",
       "xJ0xLxVdrHaCA2/C1h612RSRLAn/UN5TeG5m18d17rjLTg8CiyPPFxPMEqoSxzqEciZh0zVw7IMU\n",
       "+yXMY6o72AFVRCR7kliH0NBKZTNbAjwUSSr3ECSVLyH4Q/wRYE01fZWTSCpHnW722AgsW0swBdgI\n",
       "7ANeo7A2TW02RSR7UpFUNrMdBOWk55jZATO70t3fBNYRfOc+BdxbzWDQDOUWq40A/Wxy+IrKUEWk\n",
       "46VqLyPg88SUTC5z/sF+2LkU+q5m5kbY6ziPw/w4ckRlqCKSbpHk8vUtnyEkIa6tKyqc+6TFaoWi\n",
       "0x5eQLlkEcmSJLauSNUMIckcQuRzBvth50ehbzvF5WlDzCPPZcA/K5cgIpnRth3TSDBkVPJZgwNw\n",
       "9zDkoqGjj9PN62wJj2h7bBFJL4WM4vucMeDxwvNC6KibUwg2zNb22CKSbgoZxft5FUJHfeT5DPB/\n",
       "gWcPweQHNEsQkbRq25BRMweE8DPLho7WMY/DbA6PKHQkIumVinUI7aBS6KiH+Sh0JCKdpqNnCOHn\n",
       "KnQkIpnVtiEjmlRlVOazFToSkUzJbJWRmS03s6+Y2WYz+26l9zWryqjM5yp0JCKZktmeyu7+HeA7\n",
       "ZraaYMO71Cn0YZ6YETp6gyEuJ8/O1l6ciEgTNKuncsEfAX9b72cmqbC1xTY4dDPBnOAMYClHOY33\n",
       "AVcdh6mcWW5cm9+JSDtqJGS0DVgVPRD2VL41PL4UWGNm54av/Xvg39z9jQY+M1HR0NHMXVHfoJ9j\n",
       "8+CPl8HwSu2IKiLtqFk9lQE+DGyt9/OaJQwdHd0I3EQwKKwFRoFFPIHyCSLSrprSUxmqa4+ZRE/l\n",
       "Wrn7mJld/izcDeQKxyeAHh5jEcs4zFGAZWY2qKojEWkm9VRusnBQ+MAQ7AT6JoDNwCiTwGS4M+pQ\n",
       "DjbvNDOVoopI06incgsUkszrYfc2ODRKNHx0PAwfKXQkIq2RRE/luAeER4GzzWyJmc0HrgAerPaX\n",
       "W7UOoRJ3Hzvkfill1yh8H/gihKGj1lyhiHSqJNYhdExP5UYUEs2fBD5IofLodfr5CXBlTlVHItIO\n",
       "tHVFlUq3txgDNgL7GOA1rgS2qQ+ziDRNEltXpGpAaMXmdrXImY0Pw8rCZhY3hceH6CHPKXmY/wOY\n",
       "3KQks4g0S9tubpf2AaGwM+pS6LuaYFCAwkZ453GYT6FN8ESkmdQPoUUKlUfPwqHo8WCNws9YxN3A\n",
       "R1V5JCKZlKoZAinOIUQVZgqjM9YoBIboIs+CPBz+n5oliEhSlENIkTDJvMHhnbfAqTMTzQt5DT8G\n",
       "R1ZrUBCRJClklAKFNQoG34dym+Ed74WFN7T0IkVEaqAZQoNmSzR/nJ7jr9O/R5VHIpIUzRBSpJBo\n",
       "fgZejx6fALrpnbeIX18JCx7QwjURSTvNEGISzhQeGIXekxPN88gz/3H3w+9o4SWKSBvK3DoEMzPg\n",
       "C8CpwKPufleZ92SmyqiSSKJ5xS0wrxA++iSwBXyK7sk8J4bdXbkFEWlIZquMzOxygkY5Pwe+4e7f\n",
       "LPOeTM8Qok43e2wElhUqjz5IoUczDAF5+IwGBRGJQypyCDX2VD4H+K67fxL43w1cbya8BtcNwbHt\n",
       "BGWohR7Nhc5r/XSvb+HliYiU1ayeyi8Cr4Vvm2rgMzMhTDSvXg+7n4Hj0dcmgC5ODLzFen5uZte1\n",
       "6BJFRE5Sd8c0d3/YzJaUHJ7uqQxgZoWeyl8CvmxmFxPp9NPOwjLTMTO7bgj+EoLBIEw2G5zIDcFf\n",
       "mhkKH4lIGjSlp7K7HwU+EvNnZYK732BmDNG9vosTA6NghWTzBLAVNubMVkyC1iqISEulqqdySTu4\n",
       "zFYblQpnADe8xXp+DidyECSbtwO3wDxg5RD8VzPTVhciMqtIdVHs4h4QGuqpHHc7uLTJc2K4ED66\n",
       "nWKyGWACerfCfTmz72u2ICKVhH8o7yk8N7Pr4zp3qnoqm9nGcPRrS+5+Qx4+M0T3oX2R2VRktnDq\n",
       "MKzsB7XkFJFZmdmKkqhK4+esdx1C2FP5XUAOeAX4nLtvM7P3ACNAN7DF3W+s8nxtsw6hGmaLHuvn\n",
       "l8tGOc7tBJviRfdBWg+7D7mrJaeIzCrO785GqozWVDi+C9hVzznD0a5tcgeze+O6PAseWMc5vT38\n",
       "BPjF9CsTgMM7c2bjCh+JSDlJ5BK0l1ELBWGhgQ3wy4v6OdY/ynFO3geJo3lQS04RKStzexlVoxMH\n",
       "hAKz3DhcuXIRT9DDY4wwSXQfpG1By87HNVsQkVKp2LpC4jS5CTYfPcwHeJNTp48Wks3DkFOyWUSS\n",
       "lqoZAhnf7bQR5cJH5ZLNH4PHe4MZA5oxiHSuzO52Wo1ODhlFzRU+2gxTo+HMTvkFEUlFlZEkZXIT\n",
       "bF5+mNE+6GeI+6df2QKMQtfa4pv71sMGguiSiEhDNCCkjLuPBf0j1m8AluXZkFvHE+GrzwI/nX5v\n",
       "uO/4MpWnikgcFDJKsSB8NLySSNCon00+CqbyVBEBlZ12jCDR3L8TRvuCI9ccA7oWcda8Hl5ghDeU\n",
       "cBbpcG2bQ+islcpzKwkfAfNzMLLsMGtZxKXA7un3TgDdcMFwMeG83Mw0YxBpU5ldqRxe+F8APwbu\n",
       "cfdvl3mPZghzmBlCGqOf32GUNwH4U5gaiSSctR+SSGfI4gxhCngd6KWG7bCl1OQmGFoO9MEEeeaz\n",
       "jt8EwHia0u6kxyCXMxsHhZBEZG6N7Ha6FXgv8Iq7nx85voribqd3uPtNVhjCzH4VGHb3D5Y5n2YI\n",
       "VSguYGMZDOfKJZwBroHjPcGahV5Q0lmkXaVl64ptwKroATPrBm4Njy8F1pjZuV4cdV4j/IKS+rj7\n",
       "mPuhS4HHS1/Ls8DXcR7rOI9uurpHoXctwZAxCn0L4Yac2XjObFxbYIhIqUa2v37YzJaUHL4Q2O/u\n",
       "LwCY2T3AajP7LWAQOA34cr2fKVHR8BHAlim4retwOGNYxH/sgien362ks4jMJe4cwpnAgcjzF4GL\n",
       "3P2LwM65frldeyon4eQKJHLAssLrh1nFEE86YQhpS7jlhVY5i2RblnoqN1Sy1O49leMW/nU/BoXc\n",
       "wtBOpmcMXzmWZ0HXOs6aB0o6i7SLJHsqxz0gHAQWR54vpoaqIq1DqN9saxYCn+waYtMUYcjoGjjW\n",
       "A+cNF5POCiGJZEjq1iGEOYSHClVGZtYDPANcArwEPAKscfenqziXqoxidPK2F9uBq/b38+bpAMaJ\n",
       "V78Eb9dKZ5FsS8U6BDPbAbwLyJnZAeBz7r7NzNYRhDG6gS3VDAaShNKk8zXHoHdxnr/pBTiVPz4d\n",
       "Tky/W0lnEUnVXkZ0cIOcJETWLABTORhZVrJuYbq3QrmVzpoxiKRXEg1yUrWXkZLK8ZqZdM6Nz3z1\n",
       "fPL0/r+hSAgJeHvhVc0YRNKtkFxOc1JZUmv2EBJcs2iII8cIk8wqUxXpPKkaEFRllJxKVUjFEBK9\n",
       "ea7aX2nGUGBmgwPBwKAwkkgLpa7KKE6qMmqucs13YPMUjIbbmVxzrJ8jlOyF9IV++OxoOMvQ/kgi\n",
       "racGOdKwk5vv/OkUjHTNLFP92OMDvDGdVB6ADcOwUolnkfRIRdmpZNtcW18UTDIw/VMpJZ5F2otm\n",
       "CAJUatfZA4yGSeeho5CfETJSUx6R1tMMQWJXRdK5D9avyDN5eVhtBBVmFUo8i2RTqgYEVRm11uzr\n",
       "FooiYaT7h+BcZiaZ9/TDzuHiMYWRRBKQ6SojM1tIsEPfRnf/epnXFTJKkZNDSEHICPo/W3psIPyP\n",
       "slLieT3sLrxWeJ8GCJF4ZDVk9OfAvU38PGnAySGk/KZgG4zhvtIw0qGggxsAObMNpec6BjnNGkTS\n",
       "r5HN7WrpqbwSeAo4pcHrlSaKhpAAzHInfdnDsVwxvDS5Cdg0BNMroofg6DxgBPq06lkk3RqZIWwj\n",
       "aId5V+FApKfyuwl6I/zQzB4k2BV1IUGf5aNm9g1PS3mT1KDc9hc958FwoRJpOeQvz8N04jkfCRWV\n",
       "UvJZJF2a0lPZ3T8bPl8L/KsGg2yqshJpgwchpMjMwiidNSj5LJI+TempXHji7ttn+2X1VE6/aiqR\n",
       "Zm67PbmpMJCUzhqGFUYSqZl6KktKlYaQho5Cfk9QnTRcOBb9y3/6y75c8hkURhKZi3oqSyrVUIl0\n",
       "0l/+k2WSzwojiVQviZlC3APCo8DZYW7hJeAKYE21v6wZQvZUV4nUeBhJMweRmVLVIEc9laW8ZMJI\n",
       "mjmIJC9Vm9uhnsptoXQ2EIaRVs7cWnv97uJrhfdBP+ws7bdQafWzNtGTTqaeypIJ1S5oK501lFvD\n",
       "4O5jlWYOwbkVSpLOlETIKFUzBO1l1J7K74vU9fTMNQzBrMEj22BEf7/czAHKzyg0KEgnyepeRnNS\n",
       "lVF7qlyNdLJyyedyCehw5jCutQzSqbJQZdQQhYza18lhJKOW5HPp789GYSTpBAoZSVupNvlcLoxU\n",
       "+P0yIaMZXd0URpJ217YhI+ks1a5hCF6rLpSkLTFE6qcBQVKk7BqGTcWk9NyhJFUkidRPISNJlfIr\n",
       "mnPj1YaSVJEknUYhI2lbtSSPobotMaqpSNLsQaRJA4KZ/RbwCSAHjLn7lgrv24jKTuUktYeSqHFQ\n",
       "0dYYkjWJbIPt7k17AF3A31V4zZt5LXpk6wEMwsB48GAwODYwDnc6ePi402FgvNLv98ORO8HvBO+H\n",
       "I4XzDMD4ncWT+J3gA1D2PHrokbZHnN+dXfUOJGa21cxeNrOJkuOrzGyfmT1nZtdGjv8u8HXgnno/\n",
       "UzqXu4+5H7o0eMz+l7uZDZrlxoOHDRZ+P9wWY/d62F1t/sDMBnNm4zmz6XOJtKu6k8pmdjFwGLjL\n",
       "3c8Pj3UDzxDpqQys8ciOp2b2gLuvLnM+dyWVpQblt8TIfwH6P1tybNYvfyWiJctSkVT2Gnoqm9mv\n",
       "Ar8PnAJ8q97PFInyBhr0lDuPtsaQTteUnsru/m3g2zF/lgje4OK2SuephiqTpN2kqqdyWGVUsMdV\n",
       "bSQ1q6siqfyZyrf53ASqTJLWSaS6qKDB7PYSYCLy/J3AP0aefxq4ttmZcj06+0GDFUml5xqA8QGY\n",
       "Ppe7KpP0SM8jzu/OVPVU1joEiYPXGf4pF06q91wKJUnSkpgpNFJlNN1TGXiFYk/l9wAjFHsq31jl\n",
       "+dxVZSQJqVCRdLkHCeWKr812PlUmSRrE+d2Zqr2MUE9lSVClWUAteyWVnq90JhBWJqn/syROPZVF\n",
       "GhBnKKmR8ymcJHFwNcgRiV+lkFHwczyhpEJoSuEkiVvbhow0IEirNLrtdum5ys0CFE6SJKRipXIS\n",
       "VGUkrdLqyqToORVSkmqkqsoobpohSNo0qzKpmNxWSElqp5CRSJM0ozKp8JpCSlKPtg0ZiaRNGiqT\n",
       "RJpFMwSROsRZmRQ951whJeUXpJRCRiIpEGdlUvSc5b70lV+QSjIXMjKz1cB7gX6C7Sx2V3jfRlRl\n",
       "JBkRdzhptnMOwIa5ejNoBtFZkqgyasqA4O4PAA+Y2WnAzUDZAUErlSX7ym+/DdEwU/VbcFdL23F3\n",
       "niRWKjey5epW4GUi21+Hx1cB+4DnKNn6mmAwuKDC+bzea9FDjzQ9KLP9dnC8vi24C+fshyN3hltt\n",
       "98MRtB23Hp6e7a+3AV8G7iocCHsq30qkp7KZPRgOEF8Edrn73gY+UyT1vMGFaZUWu5Vr8xnTJYsA\n",
       "TeqpTDBAXAL0m9nb3f1v6v1ckeyqHE6CuUNKsw00s3V3i55fOQaZTbN6Kn+cYDYh0rEKf+XD+nAG\n",
       "kC/5Uh7YEAwG06njvvC9c35xzzWDUI5BqqGeyiJNlEQ4qZpzV1OlJNmQZE/luAeEg8DiyPPFBLOE\n",
       "qriqjKSjtaZCqUAhpWwI/1DeU3ieiiqjMLu9hEiVEcEA85Pw+HxgL3ButZlyYCOwotVZez30aNWD\n",
       "BCqUCuedrUpprtf1SN+DYJawkTRUGUV7KpvZAYo9ldcRTEMLPZWfrvacrhmCdDhPaK8jnyPHoJBS\n",
       "9ngC6xAaqTJaU+H4LmBXPefUSmWRSmavUILGcgySPeqHINLBZvvCr6c/Q+m5tVdSNmlzOxGZodFN\n",
       "9YJzzJ1UVuI5fTK3uZ2ItF6jISWtZWh/qRoQlEMQqVdjq6CrocRzumR2t9NqqcpIpD6e4CpoSadU\n",
       "VRmJSLokXUVUzX5Jkm1KKot0gGqqkObKMRTeE0dSWcnp+KjKSERqlmTZaq3XoRLX+KjKSERqNntI\n",
       "Kb4cw1x//Ss5nV6pGhBUZSSSbnOFlVSa2jyJ7HrapE2YzgLuAO6b5T2xbdCkhx561PYABqH/SLBp\n",
       "3p0e/Dxzc7tq3lNNK0+q3EgPGByA8QEYL/e6HtP3yeM6V1NmCO7+PPARM7uvGZ8nIrXxOctWIa6w\n",
       "ks+x0R5optEqjex2uhV4L/CKu58fOb4KGCHY7fQOd7+p4asUkcR5DGWr1ZamzvVZyjO0RiMzhG0E\n",
       "bTHvKhwws27gVoIeygeBH5rZg17DFtgiklZz77ga/vX/hSG61wPkOTGsv+qzo6GyUzNbAjxUmCGY\n",
       "2X8Grnf3VeHz/xO+9avADcAlVJg1qOxUJP2qSSrHUb6q0tTqpbns9EzgQOT5i8BF7j4JXD3XL6un\n",
       "ski6zR1Wal6eoVNlqadyQ6vcXHsZiUio2pxGp6169gR7Ksc9IBwEFkeeLyaYJVRF6xBEsm7uPANU\n",
       "t01GNTq5Gil1HdPK5BB6gGcIcgUvAY8Aa6pJKiuHINIempVnAMiZjQ/DyhltgWD3IfeqGwNlXSpy\n",
       "CGa2A3gXkDOzA8Dn3H2bma0jmOZ1A1tUYSTSWZqVZ5D41T0guPuaCsd3AbvqOadCRiJSi1q25G63\n",
       "XEPqQkZxUshIpDNUGzKqNs9QbS/odi1j1fbXIpJpzcwzQHvnGlKRQ0iCQkYinSHteYYshJfUU1lE\n",
       "pA615hqyUMrqCfRUVshIRFKnlpBRnLkGyF54qW1DRiIiUO123NGBY7gwcFT8az6O3VzbnWYIIpJZ\n",
       "ZrlxGF5ZzDVsB9bvdj9U91/zWatI0gxBRCQhnbyxXqoGBFUZiUhtqts7CWrbPykL4aXMLkwzs4XA\n",
       "bcAxgi/8vy3zHoWMRKRm1XzRx72uIU0ytzDNzD4ETLr7183sHnd/f5n3aEAQkUQkkWtIizi/O7sa\n",
       "uIitZvaymU2UHF9lZvvM7DkzuzY8HG2cc6LezxQRkeTUPSAQ9FReFT0Q6am8ClgKrDGzcwl6IhT6\n",
       "JDTymSIidZjcFISJthM8ho4Gx8ozs0Gz3HjwsMHZzmxmgzmz8ZzZnO9Nu2b1VB4lGCh+ATzs7jvK\n",
       "nEshIxFJTC0L2GpZFNfqEtU0l51W6ql8BPhwzJ8lIlK16iuHqt9HaQA2DEPf2uKhvrBcNZPJ6lT1\n",
       "VA7LTgtUfioiUiKJctOCVPVUDmkgEJEWqn5tQy2b5sUlsqndCtK0DkE9lUWkHdWyiM3Mruunez1A\n",
       "nhPD7n5Dc65y+vNbvw4h2lMZeIViT+X3ACMUeyrfWOX5NCCISKakYcFbKgaEuJmZA59HISMRyYhW\n",
       "LniLhIyuT2uVUUPUIEdEpDpJNMhJ1YAgIpIt1SegsyBVISPlEEQka2pJQCf0+cohiIh0siRyCKka\n",
       "EDRDEBGpTSp2OxURkfaiAUFERAANCCIiEtKAICIiQJMGBDM7y8zuMLP75njfxjBzLiIiszCzFSU7\n",
       "RDd+zmZWGZnZfe7+hxVeU5WRiEiNWlZlVGMfZamTZklFuhdFuhdFuhfJqDVkVHUfZTP7kJndYmZv\n",
       "jedSO8qKVl9Aiqxo9QWkyIpWX0CKrGj1BbSjmgYEd38YeLXk8IXAfnd/wd2PA/cAq939a+7+Z+7+\n",
       "kpkNmNntwAWaQYiIpFMcm9uV7aMcfYO7TwJXx/BZIiKSkDgGhNiy0uF+RgLEuaVt1uleFOleFOle\n",
       "xC+OASGOPsqowkhEpLXiWIfwKHC2mS0xs/nAFcCDMZxXRESaqNay0x3A94BzzOyAmV3p7m8C64Ax\n",
       "4CngXnd/Ov5LFRGRJNVaZbTG3d/q7r3uvtjdt4XHd7n7b7r72939xlovohPWMZRbwxFWX+02s2fN\n",
       "bNzMTou89unwfuwzs0sjx99hZhPha19q9r+jUWa22My+ZWZPmtmPzWwoPN6J9+IUM/uBme01s6fM\n",
       "7MbweMfdiwIz6zazH5nZQ+HzjrwXZvaCmT0R3otHwmPJ3wt3b+kD6Ab2A0uAecBe4NxWX1cC/86L\n",
       "gf8ETESO/RXw5+HP1wJfDH9eGt6HeeF92U9xVfkjwIXhz98AVrX631bjfTgDuCD8eRHwDHBuJ96L\n",
       "8LoXhP/bA3wfWN6p9yK89vXA3cCD4fOOvBfA88BAybHE70UaNrcru46hxdcUOy+/huMyYHv483bg\n",
       "98KfVwM73P24u79A8H/wRWb268Cp7v5I+L67Ir+TCe7+L+6+N/z5MPA0Qelyx90LAHc/Ev44n+CP\n",
       "o1fp0HthZm8D/gdwB1AoMunIexEqLbRJ/F6kYUAot47hzBZdS7P9mru/HP78MvBr4c9vZWalVuGe\n",
       "lB4/SIbvlZktIZg1/YAOvRdm1mVmewn+zd9y9yfp0HsB3AJ8CpiKHOvUe+HAP5nZo2b20fBY4vci\n",
       "jrLTRmk62NqfAAABwUlEQVTtAeDu3knrMMxsEXA/8Al3f92s+MdQJ90Ld58iWMH/FmDMzP5byesd\n",
       "cS/M7HeAV9z9R5X2KeqUexH6L+7+MzP7FWC3me2LvpjUvUjDDCGWdQwZ9bKZnQEQTu9eCY+X3pO3\n",
       "EdyTg+HP0eMHm3CdsTKzeQSDwdfc/R/Cwx15Lwrc/d+ArwPvoDPvxW8Dl5nZ88AO4L+b2dfozHuB\n",
       "u/8s/N9/BXYShNYTvxdpGBA6eR3Dg8Da8Oe1wD9Ejr/fzOab2VnA2cAj7v4vQN7MLrLgT+oPRX4n\n",
       "E8Lr3gI85e4jkZc68V78u0KliJn1ASuBH9GB98Ldr/OgcvEs4P3AN939Q3TgvTCzBWZ2avjzQuBS\n",
       "YIJm3ItWZ9PD7Pd7CKpN9gOfbvX1JPRv3AG8BPySIGdyJTAA/BPwLDAOnBZ5/3Xh/dgHDEaOvyP8\n",
       "j2M/MNrqf1cd92E5QYx4L8GX348IdsrtxHtxPvB4eC+eAD4VHu+4e1FyX95Fscqo4+4FcFb438Re\n",
       "4MeF78Rm3IumNsgREZH0SkPISEREUkADgoiIABoQREQkpAFBREQADQgiIhLSgCAiIoAGBBERCWlA\n",
       "EBERAP4/t7U6YCYFRq0AAAAASUVORK5CYII=\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1066c4490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#fit parameters give to me. Fits ok when taking bin width into account\n",
    "\n",
    "params = (0.0490546530927,7.3087377889,-6.00446454843,0.134738604245)\n",
    "pl.scatter(xval, model1_CDF2009(xval, params)) \n",
    "#pl.ylim([1, 10e9])\n",
    "pl.scatter(xval, yval, c='r')\n",
    "pl.yscale('log')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def model_ind(params, t): \n",
    "    global xerr\n",
    "    p0, p1, p2, p3 = params\n",
    "    sqrts = 8000.\n",
    "    return (p0 * (1.-t/sqrts)**p1 * (t/sqrts)**(p2+p3*np.log(t/sqrts)))*xerr\n",
    "\n",
    "def lnlike_ind(p, t, y, yerr):\n",
    "    return -0.5 * np.sum(((y - model_ind(p, t))/np.sqrt(model_ind(p,t))) ** 2) \n",
    "    # return p1*math.log(1-t)+(p2+p3*math.log(t)\n",
    "\n",
    "def lnprior_ind(p):\n",
    "    p0, p1, p2, p3 = p\n",
    "    if (0 < p0 < 5 and  0 < p1 < 10 and -10 < p2 < 0 and\n",
    "            -1 < p3 < 1):\n",
    "        return 0.0\n",
    "    return -np.inf\n",
    "\n",
    "\n",
    "def lnprob_ind(p, x, y, yerr): \n",
    "    lp = lnprior_ind(p) \n",
    "    return lp + lnlike_ind(p, x, y, yerr) if np.isfinite(lp) else -np.inf\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fit_ind(initial, data, nwalkers=32):\n",
    "    ndim = len(initial)\n",
    "    p0 = [np.array(initial) + 1e-8 * np.random.randn(ndim)\n",
    "          for i in xrange(nwalkers)]\n",
    "    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob_ind, args=data)\n",
    "    print(\"Running burn-in\")\n",
    "    p0, lnp, _ = sampler.run_mcmc(p0, 500)\n",
    "    sampler.reset()\n",
    "\n",
    "    print(\"Running second burn-in\")\n",
    "    p = p0[np.argmax(lnp)]\n",
    "    p0 = [p + 1e-8 * np.random.randn(ndim) for i in xrange(nwalkers)]\n",
    "    p0, _, _ = sampler.run_mcmc(p0, 500)\n",
    "    sampler.reset()\n",
    "\n",
    "    print(\"Running production\")\n",
    "    p0, _, _ = sampler.run_mcmc(p0, 1000)\n",
    "    return sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def model_gp(params, t): \n",
    "    global xerr\n",
    "    _,_, p0, p1, p2, p3 = params\n",
    "    sqrts = 8000.\n",
    "    return (p0 * (1.-t/sqrts)**p1 * (t/sqrts)**(p2+p3*np.log(t/sqrts)))*xerr\n",
    "\n",
    "\n",
    "def lnlike_gp(p, t, y, yerr):\n",
    "    a, tau = np.exp(p[:2])\n",
    "    gp = george.GP(a * kernels.Matern32Kernel(tau))\n",
    "    gp.compute(t, yerr)\n",
    "    return gp.lnlikelihood(y - model_gp(p, t))\n",
    "\n",
    "def lnprior_gp(p):\n",
    "    lna, lntau ,p0, p1, p2, p3 = p\n",
    "    if (-5 < lna < 5 and  -5 < lntau < 5 and 0 < p0 < 5 and  0 < p1 < 10 and -10 < p2 < 0 and\n",
    "            -1 < p3 < 1):\n",
    "        return 0.0\n",
    "    return -np.inf\n",
    "\n",
    "\n",
    "def lnprob_gp(p, x, y, yerr): \n",
    "    lp = lnprior_gp(p) \n",
    "    return lp + lnlike_gp(p, x, y, yerr) if np.isfinite(lp) else -np.inf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fit_gp(initial, data, nwalkers=32):\n",
    "    ndim = len(initial)\n",
    "    p0 = [np.array(initial) + 1e-8 * np.random.randn(ndim)\n",
    "          for i in xrange(nwalkers)]\n",
    "    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob_gp, args=data)\n",
    "    print(\"Running burn-in\")\n",
    "    p0, lnp, _ = sampler.run_mcmc(p0, 500)\n",
    "    sampler.reset()\n",
    "\n",
    "    print(\"Running second burn-in\")\n",
    "    p = p0[np.argmax(lnp)]\n",
    "    p0 = [p + 1e-8 * np.random.randn(ndim) for i in xrange(nwalkers)]\n",
    "    p0, _, _ = sampler.run_mcmc(p0, 500)\n",
    "    sampler.reset()\n",
    "\n",
    "    print(\"Running production\")\n",
    "    p0, _, _ = sampler.run_mcmc(p0, 1000)\n",
    "    return sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting ind\n",
      "Running burn-in\n",
      "Running second burn-in\n",
      "Running production\n",
      "Making plots\n",
      "32000 61 (32000, 4)\n",
      "(24, 61)\n",
      "(61, 61)\n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAZIAAAETCAYAAAAcboCHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmclWXdx/HP79xnZhBQwB0Ql0zNLStTXAPFONDoY+YK\n",
       "aWZqZmmLT+U6SJ1MLVNbnkxT7NEURU0zjzWYioL6iIoIKosmKvui7DAzZ/k9f9xn4jQOMjBnzvp9\n",
       "v169mnPPmZuLW17nO9d1/a7rMndHRERkS0WK3QARESlvChIREekUBYmIiHSKgkRERDpFQSIiIp2i\n",
       "IBERkU5RkIiISKcoSEREpFNKOkjMbD8zu9/Mfm9mJxe7PSIi8lElHSTAMOC37v5t4GvFboyIiHxU\n",
       "wYPEzMaY2WIzm97m+jAzm2lmb5nZpdnLdwNnmNkvgO0K3VYREdk0K/ReW2Z2NLAGuMvdD8xeC4BZ\n",
       "wHHAfOAlYIS7z8j5/kPu/uWCNlZERDYpWug/0N0nmtnubS4fCrzt7u8CmNl9wIlmtg64AugB/KKA\n",
       "zRQRkQ4qeJBsRH9gbs7recBAd38PuKA4TRIRkY4olSDZ4vE1M9M++CIiW8DdLR/3KZUgmQ8MyHk9\n",
       "gLBX0iH5ehjlzsxGu/voYrejFOhZbKBnsYGexQb5/CW8VMp/Xwb2MrPdzawWOB14tKM/bGajzWxw\n",
       "VzVORKRSmNlgMxudz3sWo/x3LPA8sLeZzTWzc9w9BVwENAJvAve3Vmx1hLuPdvcJXdJgEZEK4u4T\n",
       "8t0rK0bV1oiNXP878PcCN6fSTCh2A0rIhGI3oIRMKHYDSsiEYjegEhV8HUm+Zcf5fgJM2NJeSa89\n",
       "D3tg1Tsv7gCsA0a6+4o8NlFEpGRkpwEGA1fna365IoKksw/j2FGPJNPJ5uiz150BMM7dT89P60RE\n",
       "SlM+Pjtblcpke1GtWjhnWU3dVmy/zxEvoXUrIiKbpSKCpLNVW6/cfsm+mUwmdcCpP/wMcK+Z9c5f\n",
       "60RESkdFVG11hc5Wbbn7ioVvPD8vCIIagp2HA7fmr3UiIqWjK6q2KiJI8mHmX345w905+pJfrkfD\n",
       "WyIiHaYg2WBkKtncVNt9m1pVbYmIdJyCJMvdV3TbqscRgO05/FuzzOxxzZWIiGyagiRHY0P9q5lM\n",
       "JrPr52N7A5orERHpAAVJG/OnPfN2EAQqBRYR6aCKCJJ8bto4+683H57JZFL7n/LfB6FSYBGpMFVX\n",
       "/mtmu5jZX8zsjpxz3D8in5s2uvuKFfNnfRCNRmvr+n1Gw1siUlGqsfz3QMKz2s8FPluoP/TVMZdN\n",
       "yWQyHP6NUUk0vCUi8rGKsY38GDNbbGbT21wfZmYzzeytnN7H88A3zexJ4B8FbObIea8+OSUIguAT\n",
       "wy54URVcIiIbV4weyZ3AsNwLZhYAv8te3w8YYWb7AucAV7n7EKC+UA109xWz//bbg9PpdHq3Q4ap\n",
       "gktE5GMUPEjcfSKwvM3lQ4G33f1dd08C9wEnAk8B3zOzW4A5hW0pvP3Mfa9HIhE+efxFc9EQl4hI\n",
       "u0pljqQ/MDfn9Tygv7tPc/dT3P1Cd/9xoRs1b+K4Y9PpdHLAZ4f0RRVcIiLtKvgJiRvRqUNR2pSy\n",
       "bfEBV225+4q+Bw+bvv8J3/7cp0767+EzH/7VrYTnyYuIlJWcA63yrlSCZD4wIOf1AMJeSYfku5Qt\n",
       "16IpjYv3GX4+/Q48Oj3z4V9peEtEylL2F+wJra/N7Op83btUhrZeBvYys93NrJbwt/5HO/rD+VyQ\n",
       "2I6RU+659sGImQ396WOFrBwTEcm7iliQaGZjCct69zazuWZ2jrungIuARuBN4H53n9HRe+ZzQWI7\n",
       "916xas7Lp7YkW1pSqdTnVQosIuWsKxYk6sz2jv85E44d9cigljUrmHTj13Wuu4iUNZ3Z3kYXD221\n",
       "Wrdm6fvUbd0n8/lv3py3sUURkULqiqEt9Ug6/uf0Bm49dtQjJ2cymfSEn33lSWCkDsESkXKkHkkR\n",
       "uPsKdz99wevPzQ2CoHa3IedotbuICAqSzTbzL7+ckU6n2PPIE9NotbuISGUESYHmSFqNfOXPP3/Q\n",
       "zOzI74+ZqyouESknmiNpR6HmSNo66pI/re62zbY9n7z1Ilg0V1VcIlJWNEdSAp676ZyJmUyGwef9\n",
       "WmeWiEhVU5BsuZFzXnj0xSAIgoPPu2GWhrhEpFopSLaQu6+Y88SYw5LN65O9+n1yR4I+quISkaqk\n",
       "IOmkZ6874ymAwZffoSEuEalKCpLOGzl3yj+nBEEQfE5DXCJShSoiSApc/vsfWo/lTbasT/bu98kd\n",
       "0bG8IlLCqq7818yOAr5KeG7Kfu5+ZDvvKUr5bzvtePzYUY8MT6fTLc9cc/JO2jpFREpZ1ZT/uvsk\n",
       "d78QeAz4U5Gbsykjl7716qyaaDQ69KePXVPsxoiIFEoxziMZY2aLzWx6m+vDzGymmb1lZpe2+bGR\n",
       "wL2Fa+Xmc/cV0+4d/alU6ALr21dzJSJSFYrRI7kTGJZ7wcwC4HfZ6/sBI8xs3+z3dgVWuvvaQjd0\n",
       "S0y45uQXgOCY82/RXImIVIWCB4m7TwSWt7l8KPC2u7/r7kngPuDE7Pe+AYwpYBM7a91zf7wEM0sP\n",
       "uvLBL6uKS0QqXbTYDcjqD8zNeT0PGAjhMbqb+uE2FQgTuurY3Q4a2bzwnVsXTHtmYP+DBu928Ddv\n",
       "HP7KbZfcSngOvYhIUWQrWwd3xb1LJUg6VTqW7/OHOyNbrXW6mT3eZ9d9d+vdd8/MET/889Rit0tE\n",
       "qlv2F+wJra/NLG8nvZZK1dZ8YEDO6wGEvZIOKeY6ko8x8oXffHNcOp1Od+veM64hLhEpBRWzjsTM\n",
       "dgf+5u4HZl9HgVnAEGABMBkY4e4zOnCvklhHsjFmNuHYUY8Mcneejp+k7eZFpCSU9ToSMxsLPA/s\n",
       "bWZzzewcd08BFwGNwJvA/R0JkTKx7rVxNxCJRDLHXPXQscVujIhIvpX0yvaOMDMHfkLxJ9nblR3O\n",
       "uvULl43bpqaubujKhXOWv3zb9ycDI7X6XUQKLWfS/ep89UgqIkhKeWgr19E/vGttXc/e3edPm8jM\n",
       "h2/QMJeIFE1ZD21Vs0m/OvuZlnVr6P/pozNH/fjexcVuj4hIPihICmvkxF9+dVwqlUrVbdXjO30P\n",
       "+dIrquQSkXJXEUFSouW/H+HuK9z99GeuOfmFTCYT2e9LF3yOnXfVNioiUjAVU/6bT+U0R9LKzB4H\n",
       "hh876pG0u/vT8ZOeA9ahCXgRKRDNkZS/kcC41R8sOM7MbPBVfxmEDsQSkTKlHkmRDTj6tCn7DDnz\n",
       "s6lUsuWZa07RgVgiUhDqkVSQeZMeOHb+9GenRaM10S9cOnahtlIRkXKjICkyd18x46EbDlr2zrTF\n",
       "Nd26dzvie3doiEtEyoqCpES8dnfD1A/ee5Otem+fGXzlgwcWuz0iIh1VEUFSLuW/mzBy6p8uH5dO\n",
       "JmcE0Zp9DvvO7z/QMJeI5FvVlf+amQE/A7YGXnb3u9p5T1lPtrfnqB/cubpbr+16rlk2nxf/59va\n",
       "SkVE8q6aJtu/THh6YgubcT5JuXvu5m9MXLHwHXpu3z9zzFUPHVXs9oiIfJxibCM/xswWm9n0NteH\n",
       "mdlMM3vLzC7NXt4beM7dfwhcWOi2FtHIV277wbhkS/OkSBDdedBl963XMJeIlKpi9EjuBIblXjCz\n",
       "APhd9vp+wAgz25ewF9K6riJTyEYW07+3Uvn5qYOWvjN1UU237t2+8ON7VM0lIiWp4EHi7hOB5W0u\n",
       "Hwq87e7vunsSuA84EfgLEDOz35Bz1nA1mXb31a/NnzaRmq16Zr5wxbj/MrMJ6p2ISCkplTmS/sDc\n",
       "nNfzgP7uvt7dz3P377r7LUVqW7GNnPnwDeNamtaPq6mpqx105QPaTkVESkq02A3I6lTpWJtStpI8\n",
       "KXFLZbdMOR1g7/rvfHrA54fud0zDw6kP333jgiI3TUTKSM7JiHlXKkEyHxiQ83oAm1Gl5e6j892g\n",
       "UvTW478/ctX7M/66/1e+d+R2exywbOtPfP6pNXNeOU37c4nIpmR/wZ7Q+trMrs7XvUslSF4G9jKz\n",
       "3YEFhL+Bj+joD2d7JBXVE2lPNjAGWV3dhGMuvX/QoWc1fHHu1CefBj5b7LaJSHnoip5JwRckmtlY\n",
       "YBCwHbAEGOXud5rZcOBmIADucPdrO3i/iluQuCmt55kMuvKhlppoNLpq8fvzJv/h4jnoTBMR6aB8\n",
       "fnaW9Mr2jqjSIOlNONl+wdCfPvZGOp3euXn18sikG78OoJXwIrJJCpIcZubAT6iCoa2NOeSCmxb1\n",
       "6rvnTul0OjnhZ1/ZUT0SEdmYnKGtqxUkWdXYI2nLzHrvOfy8f+5+6AmfcXcm/fHqSS0LX9Mwl4hs\n",
       "lHokORQkG8TiiT2TyeRMM4u+mbidRa88pmEuEWmXgiSHhrb+k5k9Pviqh4YHQTSTbG567JlrTzux\n",
       "2G0SkdKhoa12qEfyn1on4o+56i9DIkHQZ/3KZWufv/ncSWiYS0RyqEeSQ0GycUdc/IeV3bfrt006\n",
       "lWTCNaeOc89omEtEAAXJf1CQbJyZPT5g0Mjhew06LW0GMxr/t3HB/z3cA603Eal6CpIcCpKNax3m\n",
       "+tx5N/xo+10/9a90Oh1ZNHNy5I37rwGtNxGpagqSHAqSjht0+X3ra+q6d2tZt2r9xF+e1U89EpHq\n",
       "VU1H7UoePXvdiL6rFs2ZX9t9m7pjGh5eorNNRCQfKiJIzGx0tqRNPoa7r5j8h+/tEgTBD8wscuyo\n",
       "RwbtOvhMnW0iUkXMbHCbozc6f89SHtrKhkMceB24z92faec9GtraAlZT8/jgy+4bHgQ1mUw6veKp\n",
       "+EnbFbtNIlI41TS0lQFWA3Vsxvkk0gGp1MgJPztlXKqleUJNbW3vIVf/NdV/4Imva6hLRDZXMbaR\n",
       "HwPUA0vc/cCc68PYsI387e5+vbVGptmOwI3ufmY791OPpJNi8cQuyWRyTiQSiS6fN5spd/xIFV0i\n",
       "Fa6sq7bM7GhgDXBXa5CYWQDMAo4jPC3xJWCEu8/Ifr8WuMfdT23nfgqSPDCzx4/8/h3Du/XaPpNO\n",
       "JVMTrjnlBbTeRKRilXWQAGRPQvxbTpAcTrjvy7Ds68uyb50FxIDewO/d/dl27qUgyYPWNSdH//i+\n",
       "lbXdun3DIXj9sT+wZMo/1DsRqUD5/OwslaN2+wNzc17PAwa6+3XAw8VpUnXJ9jpOB6jbse8uR37r\n",
       "luEHnnBhJv2lbx5b5KaJSIkrlSDpVLeoTSmbdgHupJali0Y+HT/p1kGXj9utbqvuhwy5+q+pd19M\n",
       "/Oudxj/+Cw11iZSlrjir/d/3LpGhrcOA0TlDW5cDGXe/vgP30tBWF4rFEzslk8l5kUgk2rxmBZN+\n",
       "dbaGukQqQCXOkUQJ50OGAAuAyeRMtm/iXjqPpIuZ2eMHnf3z4dvvvn8G90wQjQ5pbKj/yHyViJS+\n",
       "ijiPxMzGAoOA7YAlwCh3v9PMhrOh/PcOd7+2g/dTj6SLtU7ED7z41rt77dD/kXQ6bSsW/OvDV/54\n",
       "yUtoqEukLJV9jySfFCSFN+jSe9fXbNWzWyaTYeoDv3jswxnPn1DsNonI5lGQ5NDQVuGZ2eM7H3L8\n",
       "8P2Gn5eOmFkkCN4fP+r4J4C90doTkZJWEUNb+aYeSeG1DnUBFwz96WMvAHunUimf9dS9wfxJ40Bn\n",
       "nYiUPPVIcihIii8WT3wqmUxOj0Qi0VRLU/PsxjF9F7zyj+XFbpeIbJyCJIeCpDSYWe+Dz71heu9d\n",
       "9uoXBIEvmjVl5rR7Ri1DQ10iJUlBkkNzJKUlFk/0ABan0+mt0qmWyDO//AbeslpDXSIlQnMk7VCP\n",
       "pDTtf/oV7+z8qYF7mFkmnUy+9PQ1Jx9W7DaJyAbqkeRQkJSmcEK+9tZjG+4bEq2p7ZNKpTJTx/1q\n",
       "yoczJ36AhrpEik5BkkNBUvpi8cS5qVTqVjMLmteuZNKN3xrn6bUa6hIpIgVJDgVJeTCzxw/55s3D\n",
       "t+m7RwYgCIJ/NjbUx4rdLpFqpSDJoSApD61rTw447bIf9DvgqBlAz1Qq5bOfvHvmvOceeh8Nd4kU\n",
       "lIIkh6q2ylMsnjg7lUrdYWZBqmkdz/32ggeTa1d+5ARMEckvVW21Qz2S8mVmjx90Znz49nt+OmNA\n",
       "JAgmjR91/Cy01YpIl6uqHomZ9QAmEJ5Xkmjn+wqSMtU63LXP8Rd+f7dD698EtkmlUj7tL78Jlr3+\n",
       "FGirFZEuU21B8hNgNTBDQVLZYvHE11Op1O1mFqRamlqWz5m2z2tjf/ZusdslUonKOkjMbAxQDyxp\n",
       "Pdgqe30YG84jud3drzezLwLbAt2AZQqSymdmvQ/51q9f32an3fsGQQDw6vhRJ0wB13CXSB6Ve5Ac\n",
       "DawB7so5ITEgPCHxOGA+8BIwAvgq0APYD1gPnORtGqwgqUyxeKI78B6wbSqV8neefjB4d+KfQcNd\n",
       "InlR0CAxsz8BS4HngBfcfXGn/9CPHrV7OGEFQeuZ7ZcBuPt12ddnA0vd/fF27qUgqWCxeOKEVCr5\n",
       "sFkkyGRSqeY1q7763I1fH1fsdomUu4L3SMxsX+Cw7P8OBsYBN7h7Zov+0I8GySlAzN3Pz74+Exjo\n",
       "7hd34F6t5b+tVAZcYcys9/6nXfXCzvseslcQBAasfv6WSxJrFs7uj4a7RDokp+y3VeHKf83ssOz7\n",
       "Xsi+PhV4DfiCu9++RX/oR4PkZGDYlgaJeiTVIRZPGGHPeGA6nWb9yqWR5397IaSTGu4S2Uz5/OyM\n",
       "duA9xwFJM/s+4W9/7wPLgE4PceWYDwzIeT0AmNfRHzaz0agnUvEaG+odOCIWT3RrWbty2Va9dugx\n",
       "5KoHM+lUy7bFbptIuWinZ9L5e3agR3IA0N3dJ+dcOw+Y6+6NW/SHfrRHEiWcbB8CLAAmAyPcfUYH\n",
       "7qUeSRUys97b7HHQ3Yec9ZPBQTTaHci8O/nvr8x+7H+a0HCXyCaVe9XWWGAQsB2wBBjl7nea2XA2\n",
       "lP/e4e7XdvB+CpIqF4snLgJuTKfTQSbZEnnl3p+x6t3XNNwl8jHKOkjyTXttSavPfO3n87b/xAH9\n",
       "gYx7ZnU0WnNoY0P97GK3S6SUaK+tdqhHIq3MrLdFu902+LJ7P1lTW3tQ9vKiZ3594RPNH8wdDLyJ\n",
       "hrxEAPVI/oOCRNoTiyd6AVOA3dPpNGs+WBCZfNt/4y3rNOQlgoLkPyhI5OPE4ondm9aunlHTrXs3\n",
       "IJNJp1+rqa0dOH7U8f+DdhmWKqYgyaE5EtkUM+vde5+BYz93+uVHRqPRHgDL3pv5wZQxP9qB8N+/\n",
       "eilSNTRH0g71SGRzxOKJ/wLGptPpbuCRJbNeWTzj4V99Ktm0Vj0SqSrqkeRQkMiW+MKP7/lhzVY9\n",
       "rotYxCJBkAZ+0dhQf1Wx2yVSKAqSHAoS6YxYPPEr4LtABGie/czY59598p69UIWXVDgFSQ4FiXRW\n",
       "LJ4IgAeBE9LptKWTTZGp425gxdsvae5EKpaCJIeCRPIlFk90X7nk/fd7btd/OyCD+5ogGj2tsaF+\n",
       "i7YCEillCpIcChLJJzPrXdtrhzuP+u5te0draj6Vvbxq6oM3Pbtk2pO9ULmwVAgFSQ6V/0pXicUT\n",
       "/YEngH3S6TSp9WsjU+//OSvfe11DXlK2qq7818w+BXyPcIPHRne/o533qEciXSoWT+yxbtWy1+p6\n",
       "9NkayLj78mg0ekpjQ/2EYrdNZEtVXY/EzCLAfe5+WjvfU5BIlzOz3lv33+fPh5x7/YHRaHSX7OUP\n",
       "X3vwpucXT3vyIFTlJWWmrIPEzMYA9cCS1vNIsteHsWEb+dvd/frs9ROAbwN/dPe/tHM/BYkUVCye\n",
       "2B/4O9A/nU6TXLcmMm3cdax4b7qGvKRslHuQHA2sAe7KOdgqIDzY6jjC0xJfos3BVmb2V3c/sZ37\n",
       "KUikKGLxxOfWr17xbG33rXsQVnktD6LRsxob6v9e7LaJbEpZBwm0e0Li4YQTP8Oyry/LvvUF4CtA\n",
       "N2CGu9/czr0UJFI0Zta7zycPvu+zIxsOjEajO2cvLwfOGT/q+BOAoWjYS0pQJQbJKUDM3c/Pvj4T\n",
       "GOjuF3fgXgoSKQmxeOJQwoWN/QFamtampj90Y+0Hs14EbQwpJSafn53RfNwkDzqVZmY2OuelyoCl\n",
       "KBob6icDu8biicOAcUFNt/6fOeMKUi3rWtLN6/9a7PZJdcsp+83/vUukR3IYMDpnaOtywjLL6ztw\n",
       "L/VIpCQdcfEtQ2u26v1gTfcePYIgAPgQOH/8qOO/hIa8pMgqcWgrSjjZPgRYAEymzWT7x9xLCxKl\n",
       "pMXiiUHAvcDOAM3rViffeOjGumVvvQQa8pICq4gFiWY2FhhEuMhwCTDK3e80s+FsKP+9w92v7eD9\n",
       "1CORshCLJ2LAn9Lp9I5AJNWyvoVU8kd1W/f5XWNDfabY7ZPqUvY9knxSj0TKzdGX3DkyqOv+x2jd\n",
       "Vt2yQ17rgF8A1zY21KeK2zqpdBXRI8k39UikXMXiiaHArcAAwIDmRTMm/2v6g9dv7clmzZ9Il1KP\n",
       "JIeCRMpdLJ44AriNcHPIiGcykaWzXmJ24+2PNK1YfFKx2yeVSUGSQ0EilSIWT3x61eL3nu6xfb9t\n",
       "IZIBMkEQ/B/wjcaG+reK3T6pLAqSHJojkUpiZr179t3zz4ee+4tu0dq6LxAWnzgw69Wx18xaOuOF\n",
       "z6GyYekEzZG0Qz0SqVSxeKIv8CvgJKA2nU7TtGpZZGbij3ww+8Vx7hmVDcsWU48kh4JEKl0sntgW\n",
       "iKdams63oKYGyIAvC4LopcBdKh2WLaEgyaEgkWrxicEj+tX12amx/0HH7hYEQXfCSq+VwHVP/OSk\n",
       "vTyd/CIa9pIOUpDkUJBItYnFE92AS4AfAlsDkWTz+vR7L/6tZs6E+/FUs1bLyyYpSHIoSKRaxeKJ\n",
       "WuBC4Mp0Or0dEMmkUyncH6ip6/a9xob6pUVuopQwBUkOBYlUu1g8UdO0ZsV3MLumplvP2iAIIkAL\n",
       "MBH4bmND/cwiN1FKkIIkh8p/RUKxeCIATgHiwO6EpcNpYPrUB2+Yv2TahG0It2PRHEoVq7ryXzM7\n",
       "kfB8920IN3J8op33qEcikiMWTxgQA64BDgCi6XSKphVLI2899WeWvjHpgUw6dVpxWynFVnU9EjPr\n",
       "Ddzg7ue18z0FichGxOKJo4CfpJMtg4gEAZABlgVBcDNwU2NDfVNxWyjFUtZBYmZjCHsZS1rPI8le\n",
       "H8aGbeRvzz3UysxuAP7s7lPbuZ+CRGQTBl5w4xGpVMvYPgP22z4IgjrC0uE1wN3AqPGjjr8OHbZV\n",
       "Vco9SI4m/Ad8V87BVgHhwVbHAfOBl4ARwEzgOmC8uz+5kfspSEQ6KBZP7A5cAZwG9AAiQPPSt6es\n",
       "mvHYLTs1fbgQdNhWVSjrIIF2T0g8nHDip/Wo3cuyb10LnE0YLFPd/dZ27qUgEdlMsXhiB+C/gXOB\n",
       "3ul0OpLJpCNrFr+/Ilpbe8bWO+46vrGhvvTHvWWLVWKQnALE3P387OszgYHufnEH7qUgEdlCsXhi\n",
       "a+A76XTqItz7mkU8EgRp4C3C6q+HdNhWZarEIDkZGLalQUJY/ttKZcAimym7Wv7rwPcJS4ejhLsO\n",
       "LwJuevKaM/ZPN68ZguZQylZO2W+r8i7/bSdIDgNG5wxtXQ5kcifcP+Ze6pGI5EksnqgBTiYc9joA\n",
       "qAG8ef3qzPxX/ln77qRxpNat1hxKBajEHkmUcLJ9CLAAmAyMcPcZHbiXFiSK5Fl2Lcowwv28Dk2n\n",
       "092BSDrVnIxY8Gi0tu6Kxob62cVtpWyJiliQaGZjgUHAdsASYJS732lmw9lQ/nuHu1/bwfupRyLS\n",
       "hWLxxJGplubL3TNfjERrI0EQGOEWLC8Ao4AXtJV9+Sn7Hkk+KUhECiMWTxwI/AA4AehFuBYlBcyc\n",
       "Of7OZXOf/+tenklpDqVMKEhyaGhLpLCya1G+R7gWZXsgkk4lI00rlkbef/kfLJr6xEMta1edUtRG\n",
       "ykZVxNBWvqlHIlIcsXhiR+DbwNdSyZZdzSIBRgZsWRAEDwDXNzbUzy1yM2Uj1CPJoSARKa5YPLFN\n",
       "y7pVF6Wa1/2wbpsd6rJbsGQIdxp+GrgWeFnzKKVFQZJDQSJSGmLxxFaEWxudAxxIuAULhBPzr81I\n",
       "3Lp63st/38fTmkcpBQqSHJojESktsXgiCgwHvgUMJDwGglSyJbJu+eJg3uS/sWjaMw8l16/RPEoR\n",
       "aI6kHeqRiJSm7FqUIwmPAz42nWzZHotEHU+bRRYEQXAX8NvGhvrFxW1pdVKPJIeCRKT0xeKJA1LN\n",
       "Td9PNq39am3P3h4EQQ3hFiwfAgngF8BsbRRZOAqSHAoSkfKRLR2+EDgJ6At0IwyU1cAk4HrgxcaG\n",
       "+nSx2lgtFCQ5FCQi5ScWT+wEXEC4r9duQPfst5qAqcBNQKKxob6lvZ83s9HuProATa1YCpIcChKR\n",
       "8hWLJ3oTVnmdDuxNWOnVugXLW8AfgHsbG+pXt/6Mmd0GnA/8HVV/bTEFSQ5VbYmUv1g80RM4A/hq\n",
       "snn9ITV1W0UJT29MAfOAPwO3NTbULzKzCYT79YFOc9xsqtpqh3okIpUj2q377dvtdci5ux1+wpJt\n",
       "+u3lQRDdhnAj1zSwFEg8f8t391uz8J1BhCenDlWPZMtUTY/EzPYArgR6ufupG3mPgkSkQvy7txGJ\n",
       "0u+zx0084MSLPgAOBrYlezZKOp1eufjN53fcZqfdhvbccdenNjYxr3mUj1c1QdLKzB5QkIhUPjN7\n",
       "nHAx40vA0KE/fWwVcAzh2fKHAjsBtelkMhLU1KwinJj/LfCPxob6pjb30mfDxyjrIDGzMUA9sKT1\n",
       "YKvs9WFsOI/k9tzTERUkItXBzHoDy4E+uUNW2cWNhwLfBI4gLB3eirB0uIlwYv6PwP2NDfUrNSG/\n",
       "aeUeJEcDa4C7ck5IDAhPSDwOmE/428i/T0hUkIhUj00NScXiiU8TBsrRQD9g6+y3Wggn5h949qbz\n",
       "v9i0fOHh2euakG9HWQcJtHvU7uGEFQStZ7Zfln3rbcDPCY/g/Y9eSs69FCQiVSgWT3yScMhrCLAr\n",
       "OYdtrV/9AUvfnNxjwavjX29es/wLTSuXLi9mW0tRJQbJKUDM3c/Pvj4TGOjuF3fgXq3lv61UBixS\n",
       "RWLxRH/nusJ8AAALX0lEQVTCtSjDgN2BPpl0OsikknUWrVkYBMELwP8Akza2wLEa5JT9tirv8t92\n",
       "guRkYNiWBol6JCISiyd2AL4KHA98MtXSvFO0ti5CeDbKKuBNwlGOvzU21K9p/blqre6qxB7JYcDo\n",
       "nKGty4FMe0NZ7dxLCxJF5N9i8UQv4FTgK8BehMcBt07MrwXeBe4BxmYXOFbVL6MVsyCxnSCJEk62\n",
       "DwEWAJPJmWzfxL2q6h+BiHRMLJ7oTtg7GUG4/cqOQE/CHkoLsHDh9Geb/vX0PZ9dt2x+1VV3lXWP\n",
       "xMzGEm5vsB2wBBjl7nea2XA2lP/e4e7XdvB+ChIR2ahYPFFLWBH6NWAfwtLhrQFLNq2vaV63Ivrh\n",
       "v6axePqzjYd845r6atl5uKyDJN80tCUiHRGLJwLCkuGvAQcAfVMtzTsBNWQ8ZTU1C4MgeB64BXih\n",
       "7cR8pcylVMzQVj6pRyIimyO7uPFg4OvpVPLzTas+HNit1/ZrgyAICDeJXEG4Yv5W4MnGhvr1UHmf\n",
       "NeqR5Ki0/7giUjixeGLfOZMeenOPo06eSLi4cQegljBQVgEzgDHAY+NHHb+yo5815dB7UZDkUJCI\n",
       "SGeYmQ/96WP7AN8ADiecQ9mB8PTGJLB2zdJ5vPfiozt/MHPyP5tWLTt1U5Py5fC5pCDJoTkSEemM\n",
       "3N5DLJ7YE/g6cBThBpE7A1slm5uCVMu6muZVy1g5b/b03QYefwLw/sbOmC/lINEcSTtK+T+YiJSn\n",
       "7NnyZxNOzu/Ysn713hbU1JHxZKSmZkUQRBcBjcCfgJltK73K4XNJPZIc5fAfTETKUyye2A04K51K\n",
       "Dlm75P3BPXbabUkQRLcmPGhrPfAB8CzhPMqrjQ31LeWy87CCJIeCRES6WiyeGDD7iT+9v/cXvz6e\n",
       "8JCtfoSbRLZuY78ceAX405M/O/XKdMv6o7M/WrI7DytIcihIRKQQspPyuwCnE+7CsR1hoGybfUsT\n",
       "sGr5e2/0eO//Ht3xw3+99kqqac1x6pGUAQWJiBRC7mdNLJ7YiXA/rxhhkPQl3NMrkk6nW5pWLevT\n",
       "ree2E4KamruBRxsb6pcVreEboSDJoSARkUJo77Mmu+PwSYSnvvYhrPTaKdm8rldNXfcP2LBJ5IPA\n",
       "Q40N9QsK2+qNU5DkUPmviBTCxy0yjMUT2wJfBk4Aeq9c9N7gXjvvthKIAs2Ep8LOAx4D7gfmtJYO\n",
       "F3rxYtWV/5pZD+D3hP8hJrj7ve28Rz0SESkJ2S3sT1j0xvN377z/ERMI51F2IVwt30zYQ1kEjAfu\n",
       "BWaNH3V8qhifYVXTIzGzs4AP3T1hZve5+xntvEdBIiIlJdqthw+5ctyZhMNevQkDZVegjg2Bsgx4\n",
       "+qlrzrgk2bS6rIMkko+bbA4zG2Nmi81sepvrw8xsppm9ZWaXZi/3B+Zmv66KrZ1FpPylm9fR2FB/\n",
       "D+GZ8ncCbwFTCPfuyhDOp+wGnFrXZ4eitTNfCh4khA91WO4FMwuA32Wv7weMMLN9CccUB2TfVoy2\n",
       "iohsiZ8ANDbUr8wGynmEixbnEO4s/AbhL8d9anv2KVoj86VUTkg8nHDip/Wo3cuyb/0NYcA0ARPd\n",
       "fWw799LQloiUhVg8sQ3wJcIhrz4rF72335S7Gvon13xY8FXw+fzsjObjJnmQO4QFYU9koLuvI9yR\n",
       "82OZ2eicl6reEpGS1NhQvwq4LxZPPA4Mb1m55Ndh4SnDCc8/6bJV8DnVWvm/d4n0SE4Ghrn7+dnX\n",
       "ZxIGycUduJd6JCJSlqLduv8j3bw+BrwEDC3XHkmpzDvMZ8NcCNmv5xWpLSIiBZFuXt9aiVrQEMm3\n",
       "UgmSl4G9zGx3M6sl7N492tEfNrPR2W6biEjZaA2PAvdEBreZDuj8PQs9tGVmY4FBhHXVS4BR7n6n\n",
       "mQ0HbgYC4A53v7aD99PQloiUrWJ9hlXNgsSO0BYpIlLOCh0kVbdFSkeoRyIi5awSeiSlMkciIiJl\n",
       "SkEiIiKdUioLEjslW4GgORIRkU3oioWJmiMRESkizZGIiEjVU5CIiEinKEhERKRTFCQiItIpChIR\n",
       "EekUlf+KiFSRqiv/NbM9gCuBXu5+6kbeo/JfESlbKv/tYu4+x93PK3Y7RERk4woSJGY2xswWm9n0\n",
       "NteHmdlMM3vLzC4tRFsqmc5k2UDPYgM9iw30LLpGoXokdwLDci+YWQD8Lnt9P2CEme1rZmeZ2U1m\n",
       "1q9Abaskg4vdgBIyuNgNKCGDi92AEjK42A2oRAUJEnefCCxvc/lQ4G13f9fdk8B9wInufre7/8Dd\n",
       "F5jZtmb2B+Az6rGIiJSmYlZt9Qfm5ryeBwzMfYO7fwh8q5CNEhGRzVPMIMlbuVj2lEQBzOzqYreh\n",
       "VOhZbKBnsUEpPoty/wwrZpDMBwbkvB5A2CvZLCr9FREprmKW/74M7GVmu5tZLXA68GgR2yMiIlug\n",
       "UOW/Y4Hngb3NbK6ZnePuKeAioBF4E7jf3WcUoj0iIpI/haraGuHu/dy9zt0HuPud2et/d/d93P2T\n",
       "7n7t5t630tehtLf+JlvJ9oSZzTaz8WbWO+d7l2efxUwzG5pz/WAzm5793q8L/ffIBzMbYGZPm9kb\n",
       "Zva6mX03e73qnoeZdTOzF81sqpm9aWbXZq9X3bNoZWaBmb1qZn/Lvq7KZ2Fm75rZtOyzmJy91vXP\n",
       "wt3L8n9AALwN7A7UAFOBfYvdrjz/HY8GPgtMz7n2C+DH2a8vBa7Lfr1f9hnUZJ/J22zYAmcycGj2\n",
       "68eBYcX+u23Bs9gZ+Ez2657ALGDfKn4e3bP/HwX+DziqWp9Ftu2XAPcAj2ZfV+WzAOYA27a51uXP\n",
       "oqS3SNmEdtehFLlNeeXtr7/5L+B/s1//L/Dl7NcnAmPdPenu7xL+oxhoZn2Brd19cvZ9d+X8TNlw\n",
       "90XuPjX79RpgBmEJebU+j3XZL2sJf6laTpU+CzPbBfgScDvQWnxTlc8iq20BUpc/i3IOkvbWofQv\n",
       "UlsKaSd3X5z9ejGwU/brfvxn1Vvr82h7fT5l/pzMbHfCntqLVOnzMLOImU0l/Ds/7e5vUKXPArgJ\n",
       "+BGQyblWrc/CgX+a2ctmdn72Wpc/i3LeRr6s667zwd293OvPN5eZ9QQeAr7n7qvNNvzyVU3Pw90z\n",
       "hDs+9AIazeyYNt+vimdhZscDS9z91Y3to1UtzyLrSHdfaGY7AE+Y2czcb3bVsyjnHkle1qGUocVm\n",
       "tjNAtgu6JHu97fPYhfB5zM9+nXt9fgHamXdmVkMYIne7+yPZy1X7PADcfSWQAA6mOp/FEcB/mdkc\n",
       "YCxwrJndTXU+C9x9Yfb/lwIPE04BdPmzKOcgqdZ1KI8CZ2e/Pht4JOf6GWZWa+E5LnsBk919EbDK\n",
       "zAZa+Ov7WTk/Uzaybb8DeNPdb875VtU9DzPbvrXyxsy2Ar4IvEoVPgt3v8LDStA9gDOAp9z9LKrw\n",
       "WZhZdzPbOvt1D2AoMJ1CPItiVxl0skJhOGH1ztvA5cVuTxf8/cYCC4AWwvmgc4BtgX8Cs4HxQO+c\n",
       "91+RfRYzgVjO9YOz/6DeBn5T7L/XFj6LowjHwKcSfmi+SrhzdNU9D+BAYEr2WUwDfpS9XnXPos1z\n",
       "GcSGqq2qexbAHtl/E1OB11s/EwvxLEr6hEQRESl95Ty0JSIiJUBBIiIinaIgERGRTlGQiIhIpyhI\n",
       "RESkUxQkIiLSKQoSERHpFAWJiIh0ioJEpIuY2b5mdkWx2yHS1RQkIl3nGMKtXEQqmoJEpAuY2XDg\n",
       "XGCX1p1XRSqV9toS6SJm9jd3P6HY7RDpauqRiHSBbC9kUbHbIVIIChKRrnEIMNnMDjGz7sVujEhX\n",
       "UpCIdI0FhOdc93T3dcVujEhX0hyJiIh0inokIiLSKQoSERHpFAWJiIh0ioJEREQ6RUEiIiKdoiAR\n",
       "EZFOUZCIiEinKEhERKRT/h9I4WIC2UI1fAAAAABJRU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10db47450>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEGCAYAAACAd+UpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAF0dJREFUeJzt3X+wZGV95/H3d+6duXMv42W4gwVBWGGzakiVWY0pg5FU\n",
       "ruvyI1aWBMmW2QoJEZwkld0dszO6RKytZbd2YtySASaWZRlR0Vjmx27hwmbNoIZZNanIGpmAIv7Y\n",
       "wlpEwcCAIzDC/PjuH+dp76Hn3vnRfXq65+H9qurqc55zTvfzdPf9nHOf8ysyE0lSnVaNuwKSpNEx\n",
       "5CWpYoa8JFXMkJekihnyklQxQ16SKjbykI+ID0TEwxFxz1HM+8KI+HRE/H1E3BERLxh1/SSpZsdj\n",
       "S/6DwMVHOe+7gA9l5j8F/jPwjpHVSpKeA0Ye8pn5WeCxdllE/GhEfCIivhARn4mIl5RJ5wJ/VYZ3\n",
       "Ar846vpJUs3G1Sf/PuDfZuZPAW8F3lPK/x64rAxfCjwvIk4ZQ/0kqQrTx/sNI2Id8CrgzyOiV7ym\n",
       "PL8FeHdE/AbwGeBB4MDxrqMk1eK4hzzNfw+PZ+bL+ydk5ncoW/JlZXBZZu45zvWTpGoM3V0TEd+M\n",
       "iLsj4q6IuPNI85fQvj8ifrksHxHxE2V4Q0T06vQ24KZh6ydJz2Vd9MknsJiZL8/MV/ZPjIiPAX8D\n",
       "vCQiHoiINwK/ClwVEbuALwGXlNlfA9wXEV8Fng9s7aB+kvScFcNeajgi7gd+KjMf7aZKkqSudLUl\n",
       "/6lyOOTGDl5PktSRLna8vjozvxMRzwc+GRH3lWPjJUljNnTIlyNiyMx/iIhbgFcCPwz5iPDWU5I0\n",
       "gMyMI891eEOFfETMAVOZ+f2IOAm4EPhP/fN1UdFJFRHXZua1467HqNTcvprbBrbvRNfVBvKwW/Kn\n",
       "AbeUk5qmgY9m5u1D10qS1ImhQj4z7wde1lFdJEkd83ryw9s57gqM2M5xV2CEdo67AiO2c9wVGLGd\n",
       "467AiWDo4+SP+AYRWXOfvCSNQlfZ6Za8JFXMkJekihnyklQxQ16SKmbIS1LFDPkBRcRFERtubx5x\n",
       "0bjrI0nL8RDKATShPn8LbJ9tSjbthT2XZuaO8dZMUi26ys5x3P6vAgtbYNssXNErmIXNWwBDXtJE\n",
       "sbtGkirmlvxAdl8Hm84H2t011421SpK0DPvkB9T0yy9sacZ2X2d/vKQudZWdhrwkTSCvXSNJOiJD\n",
       "XpIqZshLUsUMeUmqmCEvSRUz5CWpYoa8JFXMkJekihnyklQxQ16SKmbIS1LFDHlJqpghL0kVM+Ql\n",
       "qWKdhHxETEXEXRFxWxevJ0nqRldb8m8G7gVGe3F6SdIxGTrkI+JM4HXA+wFvDiJJE6SLLfnrgbcC\n",
       "Bzt4LUlSh4YK+Yj4BeC7mXkXbsVL0sSZHnL5nwEuiYjXAWuB+Yj4cGb+enumiLi2NbozM3cO+b6S\n",
       "VJWIWAQWO3/drm7kHRE/B7wlM/9FX7k38pakYzSpN/L26BpJmiCdbcmv+AZuyUvSMZvULXlJ0gQx\n",
       "5CWpYoa8JFXMkJekihnyklQxQ16SKmbIS1LFDHlJqpghL0kVM+QlqWKGvCRVzJCXpIoZ8pJUMUNe\n",
       "kipmyEtSxQx5SaqYIS9JFTPkJalihrwkVcyQl6SKGfKSVDFDXpIqZshLUsUMeUmqmCEvSRUz5CWp\n",
       "Yoa8JFXMkJekihnyklSxoUI+ItZGxOcjYldE3BsR7+iqYpKk4U0Ps3Bm/iAiXpOZT0XENPC5iDg/\n",
       "Mz/XUf0kSUMYursmM58qg2uAKWD3sK8pSerG0CEfEasiYhfwMHBHZt47fLUkSV0YqrsGIDMPAi+L\n",
       "iJOBHRGxmJk72/NExLWt0Z390yXpuS4iFoHFzl83M7t7sYj/AOzNzHe1yjIzo7M3kaTngK6yc9ij\n",
       "a06NiPVleBa4ALhr2EpJkroxbHfNjwA3R8QqmhXGRzLz08NXS5LUhU67a5Z9A7trJOmYTUR3jSRp\n",
       "shnyklQxQ16SKmbIS1LFDHlJqpghL0kVM+QlqWKGvCRVzJCXpIoZ8pJUMUNekipmyEtSxQx5SaqY\n",
       "IS9JFTPkJalihrwkVcyQl6SKGfKSVDFDXpIqZshLUsUMeUmqmCEvSRUz5CWpYoa8JFXMkJekihny\n",
       "klQxQ16SKmbIS1LFhgr5iDgrIu6IiC9HxJciYlNXFZMkDS8yc/CFI04HTs/MXRGxDvg74Jcy8yut\n",
       "eTIzY/iqStJzR1fZOdSWfGY+lJm7yvATwFeAM4atlCSpG531yUfE2cDLgc939ZqSpOFMd/Eipavm\n",
       "vwFvLlv0/dOvbY3uzMydXbyvJNUiIhaBxc5fd5g+eYCIWA38T+ATmXnDMtPtk5ekY9RVdg674zWA\n",
       "m4FHM/PfrTCPIS9Jx2hSQv584DPA3UDvhd6WmX/ZmseQl6RjNBEhf1RvYMhL0jGbiEMoJUmTzZCX\n",
       "pIoZ8pJUMUNekipmyEtSxQx5SaqYIS9JFTPkJalihrwkVcyQl6SKGfKSVDFDXpIqZshLUsUMeUmq\n",
       "mCE/oIi4KGLD7c0jLhp3fSRpOV5PfgBNqM/fAttnm5JNe2HPpZm5Y7w1k1SLrrKzkxt5P/csbIFt\n",
       "s3BFr2AW3vxu4EVjrJQkHcLums7EP4mIa8ZdC0lqs7tmAKW75i9heym5Grgc+OCjmY+eOsaqSaqE\n",
       "3TVjlJk7Ip63B947D2cANwMPjbtaknQIu2sG9sQ74V7gEpqA3wTs3jbeOknSs9ldM4SmD35hczO2\n",
       "e1tm/v54aySpFl1lpyEvSROoq+y0u2ZAngwl6UTglvwAPBlK0qh5dM1YLXcy1OYtgCEvaaLYXSNJ\n",
       "FXNLfiC7r4NN5wO97hpgz4PjrJEkLWfoLfmI+EBEPBwR93RRoRNB0/e+50/hd4H3AhuB+d/wsgaS\n",
       "Js3QO14j4meBJ4APZ+ZLl5le3Y5XgIgNj8C2DUv98jcDm72sgaROTMwhlJn5WeCxYV9HktQ9d7wO\n",
       "bPc2+C3grPL4LbysgaRJc1x2vEbEta3RnZm583i87+jNAP+lDG8Cnh5jXSSdyCJiEVjs/HW7OBkq\n",
       "Is4GbrNP3j55Sd2YmD55SdLk6uIQyo8BfwO8OCIeiIg3Dl+tE8HubU0Xzc3l4aWGJU0er10zoOb6\n",
       "NSe/G1a9EPbthSfe6aWGJXXFa9eM0dIFym7snfG6H/i7sVZKkpZhn/xAFrY0V6C8guaxfRZOuc0z\n",
       "XiVNGkO+My9ZDfNbDXpJk8Q++QEcej35q1m6mbeHUUoanrf/G7Mm6E+5DU5bDacBG4BzgA8a8pKG\n",
       "5o7XMcvMHRFxBxy4EH6vlG4C9rgDVtLEsE9+KCedDz8O3AqcDmwHFl4x3jpJ0hK35AfU7GCdn4Pf\n",
       "LiVXAJePs0qSdAhDfmALm2Ebrfu80txE5HHPepU0Meyu6VQ+5VmvkiaJIT+w5a5ds5+IDY94rLyk\n",
       "SeEhlENownxhMxyYhX1z8J4yZROw5+1u1UsalJcaHrPmOPmFReCLcGB/E/A/vMwBTfhL0ni543UA\n",
       "S2e8butdoAy4Z6x1kqTlGPIDWdjSBHz7yJp/Dfx1Gb4beMqjbCSNnd01Q9sBvAtYDbya5rh5152S\n",
       "JoM7Xgew1F2zcbY5suZdZYoXKpPUDa9dM0blujWXwgdug+tXP7vb5n3AJeOqmiQ9i901A8rMHfDM\n",
       "PngvcBlNtw3Atyn3e/VCZZLGzpAfUHOM/FS5ds0lNNet+R3gcWAjMH+hJ0VJGjf75AcUseER2LZh\n",
       "qaum1zd/T2vcfnlJg/FkqIm0btwVkKRnMeQHtvu25qqTrwLeQtNVcy9wFnBuGd99yLHyEXFNc30b\n",
       "r3EjafQM+QGUQyjfADfQ9Mm/G1hDcwORNwCPLLfMNRHrnoT5rU03zxs3wPqtESc9E7H+SUNf0igY\n",
       "8gNZ2NLcxPsKmjtCzbIU+H9MU/4TwMlvB4iIHU24v3Suua7N6WW+q4C51fCmOXhxCf3Ysdw7StIg\n",
       "DPmhbAV+HfgxmuC+AngnS5c3iLkmtNdf2IT7GaX8fWW++8syf0yzgrgKWH9hs2XvVr2k4RnyA9l9\n",
       "HWzc3xxN849pLmdwBUvHyn8NuI8mtE+5oFkJAPwmzVmx32691l/TBH5v6/4Gmi379VsNe0nDMuQH\n",
       "NjsNrwW+BdwEHAB+Ffg3ZXgB+CiQ0awErqa53MHlNIdZbgLOoVkhwNLWvWEvqTtDHycfERfTJNIU\n",
       "8P7MfGff9OqOk4/YcDtMXQB7gJNoQn0fzc7XZ4AZmmm/C/xRWWojzVb7fcDjtwP/u7nm/JPzMLO6\n",
       "2Wn728CtNCdXtbt/oDlaZ/VTsJ/mahRTe2H3Nm9MItWpq+wcKuQjYgr4KvDPgQeB/wP8q8z8StcV\n",
       "nSQRU4/AwQ1wCvAUMAckTcD3gj5owjjK+GqagKavfJpmJdFe7mzgZJoVQ68b6PIy3D7p6r1l+R8A\n",
       "a8vrTJe6/KAMrynl+2lWPntL2dRhhp8qr9erY68tq1rTVgFPH2hee3aqmac3jfJ+a4B9Cfv2N+89\n",
       "NQ3P7IUnPwcnnw+rZsv41maZ3o1WXHlJkxLyrwL+Y2ZeXMZ/DyAz/6Drik6SiFUHm26YBZpAfWyF\n",
       "OdfSHHmTNJc7WN03fU3r+SBNSM7Q/FfQmzf7pq/hUP3L95Z5ukzvle0t9XmmVf50ea8o0+fK9Kny\n",
       "mqtYCv7esmtKHXvB/wTNfzS9svbKZl/rNYJm5dNbmewDngTWl+HeyuTp8hpzNCuLfeW9D5T3X1vq\n",
       "1h7u1bU33GtT7z+r/hXXTHnez9JKq7dC662gVrG0Yn66zDdd6tN7rGu1ZW+Zd748f4+l76v3O3iy\n",
       "zN8enynzTLfaNV3acKD1HbVX4FNluQPl0fvOYWkF3/+Z99qzv/V6vd/MVHmfXpt6dVvX+tx738HT\n",
       "5bs5WMZX9ZX9oPWavd9NsvRfaO/9n6T53bS/vwM8+7fX+y31/lPeV6b3NoyeZun3tY+l33QAaxKe\n",
       "Ll/gzHQz/WBrvNnAmNQNikk54/UFwAOt8W+VssplNH/I24DrWfqjbpuh+RFeT9ObNU/zw+zphdBs\n",
       "mefGMryX5ge7qjz6p+8v03uP5abPAlfS/HHMteowQ7MfoVd+ZXm+sTX9n5XlN9KE2vayfG/ZOZod\n",
       "yrNluevL8GvL85vKcjeU6b0AvbEst6b1mmtadV0L/GEpX1vadj3Nzure57i9DF9V6remvM/2sszG\n",
       "8l5rW22aK22a6Su7kqWgbLfjYKuON7C0Epntqx80Abi9VTZNE1rXl8c6llaKve9gttSxPb6v9Zn1\n",
       "2tibZ3up46rWZzpT3qv9mUyX6b329T7zG8r0ja3hVa22zJRpvTas7qvbla336A3PlOHe97uxr6z9\n",
       "mjPl++r9Pt9Unq9qPbe/v/Zvb2Nrmfbvrlfe+/2uas3TK5sCrozmEOXp1c3w2miGe+OzczC3tfr9\n",
       "XZk58IPm8ot/1Bq/HPjDvnlymPeYxAeckvChhCyPDyUs9D3OXGGeU8vjvPJYbp7zDjP91CNM773v\n",
       "61eY/qOt8tevML3raeeV4ZWmHesyrz9M+UptPtp5l/vehilbOEzb2t/pkeYZdvz1HUw73PAg8x7N\n",
       "tKN5Pu8IZf3D/eMLj4w7U1bI106yc9jryT9Icx5/z1k0W/PPEhHXtkZ3ZubOId9XkqoSEYvAYucv\n",
       "POSaZhr4vzR7CtcAu4Bz69+S5wDMly2BD2UzTN9jJg+dZ1XC6vKYS1jX2orrzTNVps2XR//0mda0\n",
       "5aavLWVbVph+Wat8ywrTu5y20PoctuSzP5OF0tbllplbYZn5UtZf3m7zQl/5ZYeZt7987ijL+r/f\n",
       "5cp633l/2Vzf+Mwyn2f/MocbP2WZz3W5z6z/sz2WaSsNtz/Lo513uWX6px3L8/wRyvqH2+NzCVwz\n",
       "7kwZ5ZZ8F4dQ/jxLh1DelJnv6JueWdmOV4CIOACnlH0a7nh1x6s7Xt3x2q2JOLrmqN6g0pCH3o1D\n",
       "5rc2O46gOcFpz9ub4fmtcCbN2a3rWPpBP16en6T5B+hRmnDo7Vz7fnmeykn/EUoaHe/xOgEy8/cj\n",
       "Athcju/e88Pju5vyhzbD3nnI1c3WTdAEf28L9v/RhPne/ZD7mhOcnvAYcUmdcUu+Y83W/aEn9TTl\n",
       "J70d1szCwb3wPbfQJa3I7poJtFL3jWEu6VgZ8hNo+fu+ep9XScduUs54lSRNMHe8dmr3Nti0dWl8\n",
       "E7DnkPu8StLxYndNx1ba8SpJx8LumglkwEuaNHbXdGDp8Mj5uebKlACbtkYEBr2kcbK7ZkhLh032\n",
       "7uzkkTWShmd3zcRY2NwcF3/GuCsiSYewu6Yzv8nSVjx4ZI2kSWB3zZCefZbrPcBNQD7lZQskDcMz\n",
       "XieIR9VI6pohL0kVc8erJOmIDHlJqpghL0kVM+QlqWKGvCRVzJCXpIoZ8pJUMUNekipmyEtSxQx5\n",
       "SaqYIS9JFTPkJalihrwkVWzgkI+IfxkRX46IAxHxk11WSpLUjWG25O8BLgU+01FdTkgRsTjuOoxS\n",
       "ze2ruW1g+9QYOOQz877M/FqXlTlBLY67AiO2OO4KjNDiuCswYovjrsCILY67AicC++QlqWKHvZF3\n",
       "RHwSOH2ZSddk5m2jqZIkqStD3/4vIu4AtmTmF1eYPtr7C0pSpbq4/d9ht+SPwYoV8f6ukjQ+wxxC\n",
       "eWlEPACcB/xFRHyiu2pJkrowdHeNJGlyDXV0zeFOiIqIt0XE1yPivoi4sFX+ioi4p0y7sVU+ExF/\n",
       "Wsr/NiJeOEzdRi0iLi5t+3pEXD3u+hytiPhARDwcEfe0yhYi4pMR8bWIuD0i1remHdP3OE4RcVZE\n",
       "3FF+k1+KiE2lvJb2rY2Iz0fEroi4NyLeUcqraF9PRExFxF0RcVsZr6Z9EfHNiLi7tO/OUjba9mXm\n",
       "wA/gx4AXA3cAP9kq/3FgF7AaOBv4Bkv/NdwJvLIM/y/g4jL8O8B7yvAbgD8Zpm6jfABTpU1nlzbu\n",
       "As4dd72Osu4/C7wcuKdV9l+Bf1+Grwb+YNDvccxtOx14WRleB3wVOLeW9pW6zJXnaeBvgfNral+p\n",
       "z2bgo8CtNf0+S13uBxb6ykbavqG25HPlE6J+EfhYZu7LzG+Wyv10RPwI8LzMvLPM92Hgl8rwJcDN\n",
       "Zfi/A68dpm4j9krgG5n5zczcB/wJTZsnXmZ+Fnisr7j92d/M0ncyyPc4Npn5UGbuKsNPAF8BXkAl\n",
       "7QPIzKfK4BqajY3HqKh9EXEm8Drg/Swd0FFN+4r+g1FG2r5RnQx1BvCt1vi3aP7Y+ssfLOWU5wcA\n",
       "MnM/8L2IWBhR/Yb1w7oWvfadqE7LzIfL8MPAaWV4kO9xIkTE2TT/sXyeitoXEasiYhdNO+7IzC9T\n",
       "UfuA64G3AgdbZTW1L4FPRcQXImJjKRtp+454CKUnRC2r2r3VmZkn+rkNEbGO5r/BN2fm9yOWNpxO\n",
       "9PZl5kHgZRFxMrAjIl7TN/2EbV9E/ALw3cy8K1a4Ls2J3L7i1Zn5nYh4PvDJiLivPXEU7TtiyGfm\n",
       "BQO87oPAWa3xM2nWPA+W4f7y3jL/CPh2REwDJ2fm7gHe+3job99ZPHvNeqJ5OCJOz8yHyr+C3y3l\n",
       "x/I9PnhcanoEEbGaJuA/kpkfL8XVtK8nM78XEX8BvIJ62vczwCUR8TpgLTAfER+hnvaRmd8pz/8Q\n",
       "EbfQdP2OtH1ddte0+5luBX4lItZExDnAi4A7M/MhYE9E/HQ0m1e/BvyP1jJXlOFfBj7dYd269gXg\n",
       "RRFxdkSsodlRfOuY6zSM9md/BfDxVvnRfo8f73/R463U5Sbg3sy8oTWplvad2jvyIiJmgQuAu6ik\n",
       "fZl5TWaelZnnAL8C/FVm/hqVtC8i5iLieWX4JOBCmqv5jrZ9Q+4pvpSmb3ov8BDwida0a2h2FNwH\n",
       "XNQqf0Vp2DeA7a3yGeDPgK/THDVw9qj2cHfxAH6e5uiNbwBvG3d9jqHeHwO+DTxTvrs3AgvAp4Cv\n",
       "AbcD6wf9HsfctvNp+nJ30YTfXcDFFbXvpcAXS/vuBt5ayqtoX19bf46lo2uqaB9wTvnudgFf6uXG\n",
       "qNvnyVCSVDEvNSxJFTPkJalihrwkVcyQl6SKGfKSVDFDXpIqZshLUsUMeUmq2P8Hjt6qsF7zGuwA\n",
       "AAAASUVORK5CYII=\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10dc00950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEHCAYAAABCwJb2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAER1JREFUeJzt3VGoHNd5wPH/J7npgxu1tqEvioJeImqBCiFgPTQmAqdY\n",
       "bR/kEodETolxmiakuG5NQ11CiWRIaWta3MgmKQ2RsUNjieA2rSGJHlIuyPaD7UKJwLJrQRwkGdwE\n",
       "O0nT5CG2vz7s3t7V+t7dmTuzOztz/j+44Dkzd/b4eO+Xk+98cyYyE0lSGXZ03QFJ0vIY9CWpIAZ9\n",
       "SSqIQV+SCmLQl6SCGPQlqSAGfUkqiEFfkgrSetCPiP0RcToivhARH2j7/pKk7VvETP8w8EBm/iHw\n",
       "0QXcX5K0TZWCfkScjIhXIuLcVPvhiHg+Il6MiHvGzV8BPhwR9wHXtdxfSVIDUWXvnYi4EfgJ8Ehm\n",
       "Hhi37QReAN4PXAaeAY5m5vmJ849l5i0L6rskqaarqlyUmWcjYu9U8w3Ahcx8CSAiTgFHIuKnwGeA\n",
       "q4H7WuupJKmxSkF/C7uBixPHl4CDmfk94JPzfjki3N5TkrYhM2O7v9tkIbdx0M7MaOsHuLfN62ed\n",
       "3+zcvLbp87PODW0s6hw7Fo6FYzF3LBppMtO/DOyZON7DaLZfWUQcB9Yyc61BP9bVvce862ed3+zc\n",
       "vLbp8/M+v4m69553/azzm52bbqt73Ka69553/azzm52bbqt73Ka69553/azzm52bbqt73Ka69553\n",
       "/azzm52bbpt7HBGHgENz+jFXZMWXqIxz+o/nxkLuVYwWcm8CXgaeZmIht8L9so3/1RqCiDiemce7\n",
       "7scqcCw2OBYbHIsNTWNn1ZLNR4GngH0RcTEi7sjM14E7gTPAc8DpqgFfb7HWdQdWyFrXHVgha113\n",
       "YIWsdd2Boag802/9g0cLuffSXnpHkgZrIr1zrMlMv9Ogb3pHkupZSnpHkjQMTap3Gmu5ekeSBmvp\n",
       "1TttM70jSfWZ3pEkVWbQl6SCmNOXpB4wpy9JBTKnL0mqzKAvSQUx6EtSQVzIlaQecCFXkgrkQq4k\n",
       "qbJO0zuSNEQR8QZcM55Uv/ZmZu7stkcbDPqS1KJRwN+1A+4ft9y1IyLeWJXA70KuJLXqmnHAv32i\n",
       "7e7GqfS2FnI7Dfq+81KSqhlPjtci4liT+5jekaRWvfYm3DUxs78L+PGbnXVniiWbktSyRS7kNo2d\n",
       "zvQlqWWrsmi7Gev0JakgBn1JKoglm5LUA+69I0kd6uqpWxdyJWnJVv2p21kM+pJU22Keul2GXnRS\n",
       "ktQOZ/qSVNtqP3U7iwu5krQNLuRKUkH6sGi7GXP6klQQZ/qSxGq/7apNPpErqXh9qLv3iVxJaknE\n",
       "tXll3f3DwN1kvrpyMapp7DSnL0kFMacvST2uu6/L9I4k0Z+FXOv0JakFqxrk22bQlzRofZnBL4tB\n",
       "X9Jg9aEUc9kM+pIGrL9bIC9K0f/yklQaZ/qSBqycUsyqLNmUNGhDW8hduZLNiHgHcAJ4DfivzPyb\n",
       "tj9Dkqrqe5Bv2yJy+geAxzLz94F3L+D+kqRtqhT0I+JkRLwSEeem2g9HxPMR8WJE3DNufgr4RER8\n",
       "G/hWy/2VJDVQKacfETcCPwEeycwD47adwAvA+4HLwDPAUeBm4D8y82xEfC0zP7jFPc3pS2pkaPn6\n",
       "Kpayy2ZmnmWUo590A3AhM1/KzJ8Dp4AjwL8DfxwRXwS+u92OSdIsVz54dT+wa8eoTbM0WcjdDVyc\n",
       "OL4EHMzM7wC3VrnB+CUq63yZiqQaynjwqq2Xp6xrEvQb13pm5vGm95CkIRtPhtfWjyPiWJP7NQn6\n",
       "l4E9E8d7GM32K/N1iZK2r6wHr5b+usSI2As8PrGQexWjhdybgJeBp4GjmXm+4v1cyJXUiAu59VUt\n",
       "2XyUUSnmvoi4GBF3ZObrwJ3AGeA54HTVgC9JbcjMnZmvxuhn+AG/DZ1uwwDci+kdSZprIr1zrMlM\n",
       "3713JK2cEtM2Va3c3juS1IQvPlmsToO+1TuS3qqM+vu62qre6TToW6cvSdWs1+t3WacvSQtQVv39\n",
       "srmQK2nluJC7tV4v5JrTl7QZg/xbLf2J3LY505ek+no905dUHlM33TLoS1oaa/C7Z05f0hJZg79d\n",
       "1ulLUkGs05fUQ9bgd83qHUlL5UJuM1bvSOoVg3y3XECRpIJYvSOpFaZtFssnciWtjI36+xPjltEC\n",
       "rYG/feb0Ja0A6+/7wv8oklQQZ/qSWmD9fV+Y05fUChdyl8OcvqSVYJDvB0s2JakHLNmUtBSmbVaL\n",
       "6R1JC+P+98Nj0Jc0g/X3Q+N/PEkqiDN9STNYfz80LuRKmsmF3NXiQq6khTLID4tBXyqYs/jyGPSl\n",
       "QlmOWSafyJWKZTlmn7T1RG6nQT8zj3f5+ZLUF+PJ8VpEHGtyH9M7UrEsxyyRJZtSwVzI7R9LNiVt\n",
       "m0G+PC7aSFJBnOlLA2XqRpsx6EsDZA2+tmLQlwbJGnxtzi+BJBXEmb40SNbga3PW6UsD5ULuMK1c\n",
       "nX5EvBf4yPje+zPzN9r+DEnzGeS1mYXN9CPiCPCrmfmlLc4705e2wRl82ZrGzkoLuRFxMiJeiYhz\n",
       "U+2HI+L5iHgxIu6Z+rXbgK9ut2OS3urKUsz7gV07Rm1SNVWrdx4CDk82RMRO4MFx+37gaERcPz73\n",
       "TuBHmfm/LfZVEtfsgBOMSjFvZ/TP11iFp8oq5fQz82xE7J1qvgG4kJkvAUTEKeAIcB74GHBy3n3H\n",
       "++mvc199SZrS1j7665os5O4GLk4cXwIOQvV98t1PX6rLUszSrO+jv37c5X763dR6SgXLzJ2jHP76\n",
       "07U/diFXtTQJ+peBPRPHexjN9ivzdYlSfQb5MrWV5qlcsjnO6T+emQfGx1cBLwA3AS8DTwNHM/N8\n",
       "xftZsiltwpJMzbKsks1HgaeAfRFxMSLuyMzXgTuBM8BzwOmqAV/S5izJ1KJVrd45ukX7N4FvbvfD\n",
       "Te9I09wdU5trK73T6YZrVu9IUjXrVTxdVu9Iap0lmVqsToO+6R3pSpZkaitLr95pm9U7klTfym2t\n",
       "LGk2SzLVJYO+tES+sFxdM6cvLZUlmdoeSzYlqSCWbEq9ZEmmumX1jrRkLuSqiV5X75jTV4kM8toO\n",
       "6/SlFeQsXovW65m+NCSWY6oPDPpSayzH1OrzCylJBXGmL7XGckytPqt3pJa4Q6YWyeodSSqQ1TvS\n",
       "ElmSqb4z6EsVWZKpITDoS5VZkqn+8wsrSQVxpi9VZkmm+s+STakiSzLVJUs2JalAlmxKLbEcUyUw\n",
       "6EtYjqlyGPQlwHJMlcIvtSQVxJm+BFiOqVJYvSONuZCrPrB6R2qJQV4l8OEsFcOZvPrMh7OkGjZK\n",
       "Mk+MW0Y5ewO/+sb0jlSJJZkSWLIpSUVxpq9CWJIpgTl9FcSFXA2BOX2pIoO8ZNDXgDiTl+Yz6GsQ\n",
       "3CVTqsagr4GwJFOqwj8KSSqIM30NhCWZUhWtl2xGRACfA94OPJuZj2xxnSWbapULuSrBKpZs3gLs\n",
       "Bn4AXFrA/aVNGeSl+Srl9CPiZES8EhHnptoPR8TzEfFiRNwzbt4HPJmZnwY+1XJ/VbCIeCPi2hz9\n",
       "xBtd90fqo6oLuQ8BhycbImIn8OC4fT9wNCKuZzS7/+H4MnOqasWVJZn3A7t2GPil+iqldzLzbETs\n",
       "nWq+AbiQmS8BRMQp4AjweeCBiLgRWGuroyqdJZlSG5rk9HcDFyeOLwEHM/NnwMer3GD8EpV1vkxF\n",
       "kqa09fKUdU2CfuOyn8w83vQeKoUlmSrTeDK8tn4cEcea3K9J0L8M7Jk43oPVOlqQzNw5yuGvp3R8\n",
       "65W0HU2C/rPAu8a5/peBDwFH69zAd+RqXZUae4O8SrbUd+RGxKPA+4DrgP8GPpuZD0XEbwF/D+wE\n",
       "vpyZf1X5g304S2O+v1aqrmns9CUq6lzEtXllZc7DwN1kvur3Q5qyik/kVmZ6R5KqWWp6ZxGc6Wud\n",
       "6R2pOtM7GgQ3S5OqMb2jQTDIS7OZ3lHvOJuXmuv1TF/l8B220mow6GtJ3DBNWgXm9CWpB8zpq1cs\n",
       "y5TaYcmmesOFXKk5F3K1MuYFdYO81D1z+mqF1TnSYpnT10px0zRpOZrGTkvmJKkg5vTVEl9nKPWB\n",
       "6R21xuocafF6Xb3jQu6wGOSlxXEhV0vnTF7qXq9n+uoPSzKlYTDoqyI3TJOGwD9aSSqIM31VZEmm\n",
       "NAQu5KoyF3Kl7vV6IdeSzX4xyEvdsWRTrXIWL/VDr2f6Wg2WY0rlMOgLyzGlcviHLUkFcaYvLMeU\n",
       "yuFCrgAXcqW+cCFXrTDIS2Uw6BfCmbwkMOgXwZJMSet8IrcIlmRKfdfWE7mdBv3MPN7l50tSX4wn\n",
       "x2sRcazJfUzvFMGSTEkjlmwWwoVcaRgs2RQwP6gb5CWBQX8QrM6RVJVBfxCszpFUjYFBkgriTH8Q\n",
       "rM6RVI3VOwNhdY5UBqt3BFidI6kac/o9ERHPRrw9I67NiF0/ioibu+6TpP5pPehHxKGIOBsRX4yI\n",
       "97V9/xJFxLOw4z3wC8BO4Nd3wa5vGPgl1bWImf6bwP8AvwhcWsD9B2e8kdIMO94Dv8SoLPNvgXPA\n",
       "gR1w7Z8uvHNLNn8syuFYbHAs2lMp6EfEyYh4JSLOTbUfjojnI+LFiLhn3Hw2M38b+HPg3pb7O1SH\n",
       "Zp/+ZeAEozr828f//MKi+9SVQ113YIUc6roDK+RQ1x0Yiqoz/YeAw5MNEbETeHDcvh84GhHX50Y5\n",
       "0A8ZzfaXou5MYN71s85vdm5e2/T55jOXN4BX/65q/2ZZ9FjUPW6TY7H9ezsW1a/v01hUCvqZeRZ4\n",
       "bar5BuBCZr6UmT8HTgFHIuJ3I+IfgEeAB9rqaAWHWr5+1vnNzs1rmz4/7/MnvPbmqPb+4fHPXYzL\n",
       "Ms/U6N8s866fdX6zc9NtdY/bVPfe866fdX6zc9NtdY/bVPfe866fdX6zc9NtdY/bVPfe866fdX6z\n",
       "c9NtdY+3rXKdfkTsBR7PzAPj41uBmzPzD8bHvwcczMw/qni/bh4QkKSe66pOv1HQ9sEsSVq+JtU7\n",
       "l4E9E8d7sFpHklZak6D/LPCuiNgbEW8DPgT8WzvdkiQtQtWSzUeBp4B9EXExIu7IzNeBO4EzwHPA\n",
       "6cw8v7iuSpKa6mzDNUnS8rn3jiQVZKWCfkRcHRHPRMTvdN2XLrl/0UiM/GVEnIiIj3bdn65FxHvH\n",
       "34kvRcSTXfenKxHxjoj454j48sROAEWKiP0RcToivhARH6jyO6u2tfKfAae77sQKcP+ikVuA3cAP\n",
       "KHscAMjMJ4AnIuII8HTX/enQAeCxzPyniDjVdWc6dhh4IDOfiIh/BR6b9wsLnenX2bMnIn6T0YLw\n",
       "9xfZp664f9FIzXHYBzyZmZ8GPrX0zi5BzfFYdxvw1eX1cvFqjsNTwCci4tvAt5be2QWrORZfAT4c\n",
       "EfcB11X6gMxc2A9wI/Bu4NxE207gArCX0V7B/wlcD3yO0TaSZ4CvM15kHspPnbGYOP824Gtd973D\n",
       "78RHgA+Orznddd+7Ho/xuXcC/9h1vzv+XvwJcOP4mkH9fWznOzFx/utV7r/Q9E5mnh1v3zDp//fs\n",
       "ARj/37MjmfkX4+Pbge/n+N9kKOqMRUT8GnAz8Cssd/+ihaszDsDngQci4kZgbXm9XJ6a43Ee+Bhw\n",
       "coldXIqa4/AN4LMRcRvw3SV2cylqxoqfAp8Brgbuq3L/LnL6u4GLE8eXgIPrB5n58NJ71J1NxyIz\n",
       "/xr4l2661ImtxuFnwMe76VKntvwbyczjXXSoI1t9L74D3NpNlzqz1Vh8D/hknRt1Ub0zqBl8Q47F\n",
       "iONwJcdjxHHY0NpYdBH03bNng2Mx4jhcyfEYcRw2tDYWXQR99+zZ4FiMOA5XcjxGHIcNrY3Foks2\n",
       "3bNnzLEYcRyu5HiMOA4bFj0W7r0jSQVZqW0YJEmLZdCXpIIY9CWpIAZ9SSqIQV+SCmLQl6SCGPQl\n",
       "qSAGfUkqyP8BBYeBLQNQr94AAAAASUVORK5CYII=\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10dc70710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "if __name__ == \"__main__\":\n",
    "    np.random.seed(1234)\n",
    "    truth = [fitps[0], fitps[1], fitps[2], fitps[3]]\n",
    "    data = (xval, yval, np.sqrt(yval))\n",
    "    # Fit assuming GP.\n",
    "    print(\"Fitting ind\")\n",
    "    sampler = fit_ind(truth, data)\n",
    "    # Plot the samples in data space.\n",
    "    print(\"Making plots\")\n",
    "    samples = sampler.flatchain\n",
    "    #x = np.linspace(250, 4600, len(xval))\n",
    "    #xbin = x[1]-x[0]\n",
    "    #xwidth = np.ones(len(xval))*xbin ##need to add binwidth in?\n",
    "    pl.figure()\n",
    "    pl.errorbar(xval, yval, yerr=np.sqrt(yval), fmt=\".k\", capsize=0)\n",
    "    modMat=np.empty([24, len(xval)])\n",
    "    i=0\n",
    "    #for s in samples:\n",
    "    #    cov[i,:]=model_ind(s,xval)\n",
    "    #    i+=1\n",
    "    print(len(samples), len(xval), samples.shape)\n",
    "    for s in samples[np.random.randint(len(samples), size=24)]:\n",
    "        #randint(32000, 24) returns array of length 24 filled with random ints from 0 to 32000\n",
    "        #samples(randint) retuns an array of length 24 filled with position of the parameter in sample definited by the random int\n",
    "        modMat[i,:]=model_ind(s,xval)\n",
    "        pl.plot(xval, model_ind(s, xval), color=\"#4682b4\", alpha=0.3)\n",
    "        i+=1\n",
    "    pl.ylabel(r\"$y$\")\n",
    "    pl.xlabel(r\"$t$\")\n",
    "    pl.yscale('log')\n",
    "    print(modMat.shape)\n",
    "    covMat = np.cov(modMat.T)\n",
    "    print(covMat.shape)\n",
    "\n",
    "    ys=[]\n",
    "    xs=[]\n",
    "    for i in range(covMat.shape[0]): #row\n",
    "        for j in range(covMat.shape[1]): #column\n",
    "                ys.append(covMat[i][j])\n",
    "                xs.append(abs(xval[i] - xval[j]))\n",
    "    pl.figure()\n",
    "    pl.scatter(np.array(xs), np.array(ys))\n",
    "    \n",
    "    pl.figure()\n",
    "    for i in range(modMat.T.shape[0]):\n",
    "        if i == (modMat.T.shape[0]-1):\n",
    "            pl.scatter(modMat.T[i], modMat.T[0])\n",
    "            continue\n",
    "        pl.scatter(modMat.T[i], modMat.T[i+1])\n",
    "    pl.yscale('log')\n",
    "    pl.xscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting GP\n",
      "Running burn-in\n",
      "Running second burn-in\n",
      "Running production\n",
      "Making plots\n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEDCAYAAAA4FgP0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJztnXeUW9W5t589mhnXccXYYJvYFIMpobcAwZRYEJFAGgFD\n",
       "eiHhS28kuZHxRTedNG4KIYUEEiDJDaFEEJlmQuiEjgFjArjghvvMuEzZ3x/vPkiWNfYUaSSNfs9a\n",
       "s5C29jnaOnjtd7/dee8RQghRu9SVewFCCCHKiwSBEELUOBIEQghR40gQCCFEjSNBIIQQNY4EgRBC\n",
       "1DgSBEIIUeNIEAghRI1TUkHgnNvfOfcn59zPnXPvKuV3CSGE6B2l1ghOA/7Xe38h8P4Sf5cQQohe\n",
       "0GNB4Jz7rXNuhXPuqbzx05xzzznnXnDOXRSGrwbOcc59DxhbhPUKIYQoMq6ntYaccycAzcBV3vuD\n",
       "wlgMeB44FVgKPAyc671/Nufzv3rvzyri2oUQQhSB+p5e4L2/xzk3JW/4KGCh9/5lAOfcdcCZzrlW\n",
       "4OvAMOB7fVqpEEKIktBjQdAFE4HFOe+XAEd7718BLijSdwghhCgBxRIEva5l7ZxTHWwhhOgF3ntX\n",
       "jPsUSxAsBSbnvJ+MaQXdolg/ptpxzs3x3s8p9zoqAT2LLHoWWfQsshTzEF2s8NFHgH2cc1Occ43A\n",
       "e4Gbunuxc26Oc25Gb788nkpLkAghagLn3Azn3Jxi3rM34aPXAvcB05xzi51zH/LetwOfAjLAfOBP\n",
       "UcRQd/Dez/Hez+vpWgDiqXQL0Nqba4UQotrw3s8rtlbUY0HgvT/Xe7+7936Q936y9/7KMH6r935f\n",
       "7/3e3vtv9+SefdQIbgIa46n0ob28vpKYV+4FVBDzyr2ACmJeuRdQQcwr9wLKTSk0gh7nERQb55zv\n",
       "q4/glItvbO9o39p+97fOvhOY5b1fV6TlCSFERVKMvTNiQBSda1396sb6hkGDgNOBX5Z7PUIIUU0M\n",
       "CEHwwM8uvN97z/Ff+N0mYLRzblS51ySEENVCRQiCvkYNAbPaWjZsahw+agjwFqQVCCEGKPIR7Pg+\n",
       "t5xy8Y2nb1ixaP3Dl396ivwEQoiBjHwEhZnV3t62dcT4ycMlBIQQovsMGEHgvV83aNDgY8C5A875\n",
       "xivOuVvkKxBCiJ0zYAQBQCaZeKyjo6Nj/LQj9kARREII0S0GlCAAePmBm56pq6tjt2Pf8SyqfCqE\n",
       "EDul1D2LJznnrnfO/Sana1mheX2NGnqdl2//3UmdnZ3t+57yvr2Ba2QeEkIMJCqi1lAPOQjrTPYR\n",
       "oMsSEH2pNVTgXuuWPHHXK7FYrGHstGNlHhJCDCgqotZQD3sW3wd83Dl3B/CPIqy3W7xw02ULOjs7\n",
       "OejsL7cBE+U4FkKIrumNRnAlcFruQOhJ/NMwvj9wrnNuOvAh4Bve+1OARB/X2hNmrVzw0PxYLBYb\n",
       "O+3Y45DjWAghuqQ31UfvAdbmDb/es9h73wZcB5wJ3Al81jn3C+Clvi62B2tc9/R13zqgs7Oz86Cz\n",
       "vwzwMHIcCyFEQYrlIyjUs3ii9/5J7/27vfef9N5/pUjf1W3at2y+KhaLdU44LB5DjmMhhChI2XsW\n",
       "g0UN5bydVyzH8T+/e85HTp59w/v3e+vHD1v+aAbMPPTeYtxbCCH6kxBZOaMU966UnsVzirSO7Vj1\n",
       "wr8Xjd/3yD132fdND7/2/H0yDwkhqpJwQJ4XvXfOXVysew+InsU74unr/ufwzs7O9gPf86VDgBsU\n",
       "QSSEqGYqIo+g0noWd+Pe6xoaGv4Qi8Vig3Y/5EQUQSSEqGJKkUdQEWWogf+miL6BQpw8+4Y2oP6u\n",
       "1DseBmaqQqkQohrJ8RVcPKDKUJdSI4jY0tL8xbq6us4jPvajXZGJSAhRpVREZnG1cu8Pzr8sFott\n",
       "GLHb1EmATERCCBGoGUEAkEkmRgO8+aJrQElmQggBVIggKGXUUD7tW7fe1TB4WOcbTjq/GZmIhBBV\n",
       "hnoWF4l4Kt3W1tbGXal3RHkUf/beK9FMCFE1qGdx3zmjrq6u7vCP/xBkIhJC1Dg1KQgyyUSms6P9\n",
       "tVET9uxgwp7zkYlICFHD1KRpKCKeSre3tbV13pV6R0MYkolICFEVyDRUPD5ZV1cXO/YzV4BMREKI\n",
       "GqUiBEF/Rg3lkkkmftXetvWpoaPHd057+6c9MhEJISqcqosacs4dD5yHVTnd33t/XIE5ZTMNRcRT\n",
       "6U3t7e0Nd15yViwMyUQkhKhoqsY05L3/l/f+k8Dfgd+V8rv6yFCAGd+4HmQiEkLUGKVuXh8xC7im\n",
       "LwstJZlkwm9pbb44Fot1Hve5K/dEJiIhRA1R6ub1OOf2ANZ771v6uNaS8q/vn/fNWCz24uCRY0ZP\n",
       "OPrtqkUkhKgZSt28HuDDwG/7tMp+IpNMTOvo6OjYP/5haBqxCpgozUAIMdDplbPYOTcFuNl7f1B4\n",
       "/24g7r3/WHh/PnC09/7T3bhX1I8goqR9CXaGc27UybNveM1771WCQghRKRToWVxx/Qj6FHoU+hFE\n",
       "f/OKtKbermVdfX39rLq6uroTvnQ1yHkshKgAoj4E0V8x710sQdCn5vXlyiPoikwy8ee2LZvnNg4b\n",
       "oUY2QoiKomLyCAqYhuqB54FTgFeBh4Bzu9O3uBLyCLoinkqv7OjoGPv4td+pW73gfpCJSAhRIZQ1\n",
       "j6AUzesrlUwysWtnZ2fnwedcBGYi2uScmyftQAgxkKiIonP0Q/P63nLI+785auzUA19z4O+45Kx7\n",
       "sTaXIO1ACFEGStG8viIEQaWahiLiqXQcSG/ZsqXj7m++qxHTDmZ679eVeWlCiBqlakpMDBQyyUQG\n",
       "+GZDfX39CRdd04KZv+RAFkIMCCpCI6CCTUO5xFPpW4GZG1YsWv/Azy4cHYZlIhJC9BsyDVUA8VR6\n",
       "QUdHx16LHr2j7oWbL3sY0w6mAK3ALJmLhBD9QTH3TgmCXnDy7L+trauLjdjS2vy9f33/vGORA1kI\n",
       "0c/IR1Bm7rzkHaNjsdjmoU0jv7LHjHNHhGFlIAshqhJpBH0gnkpv6ujoaJj/j9/cu+zBm15CJiIh\n",
       "RD8hjaBCyCQTQ2KxmD8o8fHjdp1+/FTMRKTy1UKIqqIiBEGl1RrqCZlkogHwB77nS8cxYTIoA1kI\n",
       "UUIqptZQt2/unAP+B2gCHvHeX1VgTtWahiLiqbTraG/f6sFtXP3qqQ//7MI5yIEshCgh1WQaOguY\n",
       "CGylB9VIq41MMuFj9fWN9fX1fvT4PW4fOvmQ1vCRHMhCiIqn1D2LpwH3eu+/BHyyCOutWDLJhAca\n",
       "AX/sh/975q4HzPgnMBP4nsxEQohKpsemIefcCUAzcFVOGeoYVob6VKw3wcPAucBhwFbv/V+cc38q\n",
       "ZCIZCKahfOKpdBvggFlzZ59xITITCSGKTFlNQz3sWXw9EHfOXQbM6+Naq4bIgQxcM+mEs3PzDORE\n",
       "FkJUHPU7n9ItJgKLc94vwXoWbwI+urOL8zzgFV9zqDtkkomGeCq9Zd+Tzzt4+LjJTz53/Q9mAjeQ\n",
       "1Q5+CUg7EEJ0iwI9i4t37yJ1KHsXcFpvm9cPNNNQLvFUeiMwFLh57uwzGrE8A5WxFkL0iUqMGhpQ\n",
       "PYuLSSaZaAJWA2878et/2Q34MypjLYToJRWTR1ArPYuLSTyVXgDsBSyZO/uMl5ADWQjRB9SzuArJ\n",
       "JBPTgPuBSSd+/S/HhOGHgQucc1fIiSyEKBcVUXSOKmlMUwziqfQVHR0dH+ns7Oxc8/Iz4x6/6r/W\n",
       "OefmIQ1BCNEN1JhmgBBPpc8Grglv3zR39hlzyDqR1ehGCLFTKtFZLHpAJpn4MzAOyzW47/gvX/MU\n",
       "5kSeiQkBVTEVQvQb0gjKTDyVbgaGAPMyycQpzrlbUIipEGInDDiNYCCHj+6MTDIxHFgBzIin0i8B\n",
       "s8hqB6pTJITYhooJHy3qAmpcI4iIp9L3AscArSH3ADmRhRBdoeb1A5R4Kv0z4BNAJ9A4d/YZaXLM\n",
       "RMD3sIquciQLUeMMONOQMDLJxP8DzsEql7a96Qu/v5lgJgqb/jTkSBZCFBlpBBVIPJUeh5XtiAEL\n",
       "MsnEdAA5koUQETIN1QjxVPp5YG+gA9hz7uwzmjFN4AJkJhKipqka01Dwbt/jnPuFc+7EHcyr2aih\n",
       "HZFJJvYFPob9f3pl5iV/v9p7/16ZiYSoXaouasg592bgq8By4Jve+xcLzJFG0A3iqfQaYCTQkkkm\n",
       "RuSbiZCGIERNUVbTkHPut0ACWBlVHw3jpwE/xuzav/bef9dFK3VuV+CH3vvzC9xPgqCbxFPpq7A8\n",
       "A1rWrfnUvT98/0nABd571SsSosYot2noSuC0vAXFgJ+G8f2Bc51z031WyqwDBvVloQIyycT7sW5w\n",
       "ncNGjfnZzEv+vm/Oyb81/FctMYUQPaJY/QiOxSrhnRbefzVMfR6IA6OAn3vv/1ngXtIIekE8lX4R\n",
       "q0vUBjTNnX3GMLKO5NyWmNIOhBiAlFsjKEShnsUTvfd/895/wnt/TiEhIHpPJpnYC7gI6zu9aeYl\n",
       "f/95jiM5Vzu4oFxrFEJUB8VqXt8nj/NAbF7fH2SSiUuBS+Op9HrgvfFU+gzMoTyLoB0E/8EVyJEs\n",
       "RFVTDc3rjwHm5JiGvgZ0eu+/2417yTRUBOKp9B+ByAR0SSaZuCT6TI5kIQYeZU8oK3bPYmqoQ1kp\n",
       "iafSg4D1QAPQAozKJBOdCjUVYuBQER3KQs/iE4GxwEpgtvf+Sufc6WTDR3/jvf92N+8njaDIxFPp\n",
       "PwPvxGoW3T139hnvZFtT0TykIQhR1ZRdIygm0ghKQzyVdsBqzGfQCZyWSSbuANUsEqKaqQiNoNhI\n",
       "Iygt8VT6w5g2UAe8nEkm9gq5Ba9rCAByKAtRXQw4jUCCoPTEU+nngH2wCK9PZpKJX+V+LnORENVF\n",
       "JeYRiAonk0zsBxyNmYkuj6fSK/KmbJN74Jy7QtnJQtQGFaERIB9Bv5LTFtMDl2eSiU/lm4ukIQhR\n",
       "mchHIIpGPJUegVWFHQRsBsZlkolIK5BDWYgKRz4CUTRCn+QLsFDTRzLJxNEABTQEOZOFqCAkCETR\n",
       "iafSy4FxmA/h45lk4srcz2UqEqKykCAQJSGeSp8J/B8WRLARGJtJJjqgsKlIWoIQ5UOCQJSUeCp9\n",
       "G3ASZi66K5NMnNpF7sE8pCUIURYGnCBAUUMVR8hMXgOMwKKLfpBJJi7KnSMtQYj+pyqjhpxzw4B5\n",
       "WHXSdIHPpRFUMPFU+nTgb1ghu3bgg5lk4lrY3qEcxuYhLUGIklNVGoFz7r8xe/OzEgTVSzyVng0k\n",
       "Mf/BJmBqJplYlT9PYadC9A/V1Lz+LcAYYDDwmgRB9RNPpf8KnIn5DxZjAuH1f0SqYyRE/1BuQXAC\n",
       "0AxcldOPIIb1IzgVWIqdBs8FzgOGYQ3tNwHv8HlfKEFQneT0TPZAOpNMnNnVXJmLhCg+ZTcNdbd5\n",
       "vff+O+H9B4BV3vtbCtxLgqBKiafSo4BFmLDvBL6ZSSbm5M+TuUiI4lOJguDdQNx7/7Hw/nzgaO/9\n",
       "p7txryhqKELRQ1VGPJX+EGYOigEdwOczycTPos9lLhKi7xToWVzeqKECguBdwGm9FQTSCAYG8VQ6\n",
       "CczGHMrtwKxMMvHXQnNlLhKib1RiGeqlwOSc95OBJd292Dk3J0g7UcVkkolUJploAH6B/dv6czyV\n",
       "3hRPpeMFpm9T9rq/1ihEteOcm+Gcm1PUexZJI+hT83ppBAOTeCp9NXAOJhQ2AzMzycS9sENz0Uxg\n",
       "PjIXCbFDyh01pOb1okfEU+m/AGeRzUF4SyaZuD9/nsxFQnSfsjuLi4lKTNQO8VT678BpWA7CZuD0\n",
       "TDLxz+jzrqKL5FgWIktVlpjY6QKkEdQc8VT6ZmzDjwTCzEwycW8hcxFIUxCiEANOI5AgqE3iqfRN\n",
       "wFsxgdACTM8kE0vz56m4nRDbM+AEATIN1TTxVPpO4M2YQFiLla3YGH2u4nZCZJFpSAxYQtnrJ4Hp\n",
       "mEBYBeybSSbWF5qvbGVR6ww4jUCCQETEU+k6YCGwByYQ1gFvzDcZ7cCfIJORqAkkCMSAJwiEJ4H9\n",
       "MIHQDByRSSZe2NF1MhmJWkGCQNQMQSDcBxyBCYSNwKGZZOKlQvNlMhK1ggSBqDniqXQMEwiHkxUI\n",
       "B2eSiVdy58lkJGqFAScIUNSQ6CZBIDwEHEw2D+GCTDLxhx1dJ5ORGChUXdSQc24/4LNYOYqM9/43\n",
       "BeZIIxA9Jp5K1wN3A0dhpSs6gWuBD+R2TIuQyUgMNKpOI3DO1QHXee/PLvCZBIHoE/FU+lLg00A9\n",
       "1jHtTiC+sxaaYVwmI1GVlLvoXLd7FofxtwEXAr/y3l9f4H4SBKIoxFPpU4AbgSGYQLgbK1/R0dU1\n",
       "MhmJaqXcgqDbPYtzy1A752703m/X11aCQBSbeCp9AnArJhDA/j2elEkmNuXPlclIVCtlNw11t2cx\n",
       "cD/wTmAw8Kz3/scF7iVBIEpCPJWeAdwMDA1DrwFvziQTz0dzujIZhc/UH0FULJUoCNSzWFQs8VR6\n",
       "EvAgMCEMbQW+lNtXuRAyG4lKQj2LhSgCITntdrIF7jqAHwFf7UmkkRzMohKoRI3gGGBOjmnoa0Bn\n",
       "5DDeyb2URyD6nXgq/VPg41hwQyfwF+DDmWRiczRH/RFEJVIxeQTqWSwGCvFUOgn8F9CARRotAd6R\n",
       "SSYe6+oaOZhFJVDuqCH1LBYDjngqfR7wU2BEGNoKXAV8JpNMbMmduyMHsxD9RdlNQ8VEpiFRScRT\n",
       "6cFYQtqRWMayB+4CzsokEy07u16RRqLUVIxpqJhIIxCVSjyVngNcBDSGoWeBkzPJxMqurpH/QPQX\n",
       "A04jkCAQlUronHYe8DNgeBheAZyRSSYezZ8v/4HoLwacIECmIVHhxFPpBuAk4A+YfwwsfPRrwK8i\n",
       "P4IS1ESpkWlIiDITT6UHAYcBVwNTsHyENuAW4ENd9VgGmY1EcRlwGoEEgag2gsnoKOAy4FAsWg5g\n",
       "OfBF4PZMMvFa7jUyG4liIkEgRIUQGuXsDZyD9d4YGT7aCvwe+K9MMrEaZDYSxUWCQIgKJZ5KTwb+\n",
       "RDb8tBN4ArgUuD+/tWaEzEaip0gQCFHhhA5qc4AvAIOwfIRXgSTwt0wysSF3vsxGoqdIEAhRRcRT\n",
       "6bcCfySbtdwCXAL8PJNMtIKylUXPqRpB4Jw7E+tmNgIrO3FbgTkKHxU1QTyVngjcA+xBNtooA6SA\n",
       "RzPJRGdX18qHICKqNnw0nHYu9d5/tMBn0ghETRFPpcdjLTUPw6KNPLAGuAb4O/DvTDKxNvca+RBE\n",
       "PuUuOtejnsXhs0uBP3jvHy9wPwkCUZPEU+kDsUijs4FhmJawFfMT/AWYC7ycSSa27MyHII2h9ii3\n",
       "IOh2z2LgOeA7wFzv/R1d3E+CQNQ08VR6HBAHvo6FosawpjkvYJrD7Qvv/OPC/8y79vt04UOQxlB7\n",
       "lN1H0IOexS3ABzDB8Lj3/pcF7iVBIAQQT6VHApGW8HayPRLWY72Xf4uZjbargqqoo9qjEgVBn3oW\n",
       "SxAIsS2hreaXga9gSWqRc/lJ4CasB/N8YB3Q+dg1qaZVzz30v+AVdVQjFHPvrC/GTbBTS69xzs3J\n",
       "eavoIVHzhAii78ZT6eswn9wngOnA4ZiTeS3wKNaD+flDZyXbgF8AGwvdTz6E6qdA8/ri3btIGkGf\n",
       "ehZLIxBi5wRfwreA92AlsSPn8gtAGlgEPJZJJh7Iv1Y+hIFHJZqG+tSzGOURCNEtQrG744C3YtFG\n",
       "kzBfQicWgno71mLzQWBTTnls+RAGCBWRR6CexUJUBvFUuhHYFfgc8CFgFGam3QTcBzyARe69uHrh\n",
       "Y+1P3/CTZH3j4I80r1q8Jv9eMh1VH2XXCIqJNAIhikM8ld4fuBWYiJmNNmDO5TTwImZOagIeyCQT\n",
       "j+ReK9NR9VARGkGxkUYgRPEIWsInsGJ3EzANvQ3TEH4IvIKd/NOZZGJBdJ1MR9VHMffOumLcRAhR\n",
       "GWSSiagPwiVYMtojmCCYAfwNy0cYApwaT6UPgNdDVWcBfyYIgSBQRI0gjUCIAUg8lZ4EvAnYM/x3\n",
       "V2BfzDTksGTPxzENYBUW7HF7JpnYGE+lRwAnYCak1dE95UeoLOQjEELslHDiH4NF8u2JaQVNwMnh\n",
       "fQPQDjyN1TYagwmHwcAyLCfhUUyDGPrU9T/8x5r/PHmg7+ygLlZ/45u/eOUHgdaghUTfORzYNZNM\n",
       "/KdffmQNIh+BEKLXhJP+mZgQOBgYj9U2GoNFG60AFgCtmPAYgpXJXgQcueqFf39h8UO3TOzY3PzM\n",
       "folPfK5pwtQtYc5DmWRiXfiOU8L9HwSeyiQT5d1gBjADTiOQIBCifwh5CAcBRwCTMRPRUsy5vB/W\n",
       "TW0NJgDmY+UtngFOaNvc2rr0sdvOb17xyhuX/juzPNxrcrjuQazPwoxwv1cxbaNgn4V4Kj2K7bWJ\n",
       "YZjG0hj+BmF+zPmZZKKjyI+i6pEgEEL0iRBqei6wGliMla8Yj5XDPim87sAqDS/CwlCXA+MwwZAG\n",
       "dgvvn8MqD88P1zcANwCHYhv6w5lkoj3nuxvCd8zPJBNLcsaPwzSTFixjegsWCrswk0wsK8FjqGoU\n",
       "NSSE6BOZZGI+8H/AAVgHwYVYQtrdWMe0X2Eb/2isttH5WEXUg4CzsOTRjwFHAm8DXsKEwB7Y5j0E\n",
       "8y+0AMfmRSHtu2DulR8ssCyHCYcnMsnEs8HP8DKWPS1KiASBEDVKJpl4Aju5z8A2+RVY2YqjsKoB\n",
       "/8DKVazATucHAIdgzuR1mKDYLXzWghXC2+/VJ+56O7Z5NwFLsOKW08H8FE/f8JPvblj6wrn3/vTC\n",
       "q0P3wh2xDBgbtAjiqfTQeCp94k6uET2kpILAOTfVOfdr59xfdjJvTvCECyH6l1vCXwb4N3Atph2s\n",
       "w7SD1VgJi9OBJGYqOhJ4F9ZMZwTmJ9gLOPmVB9NHdrRtOfU/d//p+uZVi08A9scS246Mp9IzgLc2\n",
       "r3h5XNNue9EwZPgMYLseJbkEk9IqTOCAaRyD8+fl5z3EU+kRIYS2X+mP/Avn3Iy8is19pqSCwHv/\n",
       "UqE+xQXmzVHoqBD9TyaZ8Jlk4k+ZZOL3WBXTVVgI6SLMR7AUuB4YCtwG/Brrq/wcMAUrfnc+Vuvo\n",
       "TUNGjpvSvHIxLauW7P/A5Z/7IKYlRJrBSqChafzUurF7H8bw8Xu+CFzQ1dqcc3PiqfQYYDMwPZ5K\n",
       "N2ECZ5t9K55KT8SEUi6TsbyJbebFU+m98samxVPp0Tt/UtsTGglFr4fFU+mjML9ISfHez/Pezynm\n",
       "PXssCJxzv3XOrXDOPZU3fppz7jnn3AvOuYuKt0QhRD/xOGbmGYU5hPcCbglRPwuwjX8JFiH0c8xP\n",
       "cDOWuXwkcPLIydOamnadyug937js+M9ecQfWcW03zM9wMPD3PU9876X1Q4a17HPq+/+el5SW7/i8\n",
       "GHgzFrm0PyZ09gMODtnQEYMK/JZxQGPuZg3sggmUyTljh2AFNF8nnkrvGk+lYznvXTyVPiL3tB+i\n",
       "no6Pp9KD46n0dOB4THt6uMBaKp7eaARXAqflDoSexT8N4/sD5zrnpvd9eUKI/iKYYR7BIndmAktz\n",
       "EsOWYOXlF2Kb3QFYKYvbsQ37+8Ci+sFNbrdDT2K3g08ePXjE2DimVUTZyY3A4UNGj99/1MRp9zUM\n",
       "Hjoud8PNJWQxs+Tfme+vfvHxzcAdWDjq9Zjg6fIUH0+lB2Hmo0WYWSpiMBbZND2eSo+Np9JvwPar\n",
       "qXm3OBAr0BexKybM9si5/wws+/rzWMRUY7jPwV2tq5LpsSDw3t+DqXu5HAUs9N6/7L1vA64DznTO\n",
       "jXHOXQ4cIi1BiMonk0w0A/dilUvvzBn3mWRiEfAaFm66FYsgehkz/bwCXBKLxa5oWbXkhVgs1oz5\n",
       "FS7H6hgdgmkPTZjzuAnbbE8Kp+t8pgHUxRr2XvnsfZdhmspUzHm8Fgtv7YpdwjpfJetbAItkWoVF\n",
       "M52E9XV4HhgXT6VzuzXm+yAmY1FRU+Kp9BSsSuuo8Iwuxw7B/wjPbQFVSLFaVU7E/nFELMF6Fq/B\n",
       "KiEKIaqETDKxPJ5K35xJJgq1vXwNM7v8E9tM78c25U2YtnDziAlvWIZtlCOBYzCz0GHY6ft2TJM4\n",
       "FfMZrAKOiqfSr+R9TytAbNDQ1yYcfNKXscPmemwzX4ed0ueHuUOwU/zN4f24sM4NmHloOBbVNBjz\n",
       "N3jMnDQMq866HNgdWBTMP7lmoUZMsDyEaQBDgXuwkNn1mWQiOhS3BR/Gflil16qiWM7iPvcszvmb\n",
       "UaQ1CSF6SRdCAMzMswsWUTQSeCNZG/tt2Gl9MbbhL8XMxT/B/A9vAC7EnM27YwfR9rx7RZvwLOCZ\n",
       "XfY+9I7Rk6c/j23aG7Cw1BZsgx8S5k4A3hhPpQ8J73fB9rbTMT/HOdjGfR7wScyxfDfW1W1muGaP\n",
       "8N/onhGTMEFxFOYbWUNWY2jJm+uwZLqSEEULRX/FvHexBMFSTH2KmIxpBd0iRA3NUfSQEJVNJplo\n",
       "xez+DvMVTMfMQ4swIRAdChdgm/wbsZP8d4HZWJhqK6ZFnAvchEUcrcQ2++EA3vt13vsDYw2DHNnT\n",
       "fwt28ifMHw/wfObKj2BC4i05IaM+rOf/sAinx8K698Q2+wMwU5PDtJch4URfyCz0GuaXeAzTJPbE\n",
       "TGP5gmAkeZFKxSSKFor+innvYgmCR4B9nHNTnHONwHux/8HdQpqAEFXFa9ip+17Mxn430JlXU2gj\n",
       "Zio6Htv4J2MC4weYMHgWC1cdHcYeBN4DnJBz0gfTNqYB/wH+hZ3Qx2CCYFfn3BVtm5rPXbXg4UTr\n",
       "6lcXYglxG8N1e4e1DsdMVA4zK0WVVncP3zES02L2IEcjCBFH9ZgQaA4F9F7ChNEQimda7xEVkUcQ\n",
       "ehbfB0xzzi12zn3Ie98OfIpsoao/dadxfYQ0ASGqikgQLMWEQTtmZolwWOJZa/hsApZ4tgXbzF/E\n",
       "bO33YSagizDH7kGY8/X+eCo9K55Kj8dMLZOAZaHC6eOY43kVMNbV1U0bPHIMbZuaJ/776os/gG36\n",
       "UzAT0whsjxuEOac3YL6MwzFzz5PhN+yNCYKoNEbEHmF8ODAtnkrvg2kRx4Xf9IY+PMNeUxF5BN77\n",
       "c733u3uo3YLkAAAS6klEQVTvB3nvJ3vvrwzjt3rv9/Xe793dxvUR0giEqCpWk/ULbMJOxtv4CUMo\n",
       "6lZsg98jfL4M25wfwYTEbtgmuxb4NFbO4lFsg/1NmHdQeB/5Dh4DxmWSiTZg/aARY9sbh47Ed3Ys\n",
       "P/Adn/8WtsGDOW3rME1kCLBl1YJHjsDMSe1hHc+G9TVgG3sz2bpGDtMYFmNC5NBwzQZMmC0q8Fwa\n",
       "Me2ipFSERlAKpBEIUT1kkonN2CbfFIaGAq0Fyis0Y5rCVsy8sxQzBa3BNtQoP2BxJpn4F3bafx9m\n",
       "evl5uMdIzMF7SzyVfie2kUdCYeURH/zWLzo72hl/wPFzR79h/wmYtrEW2/z3wJzBDwMjtm5qPgbb\n",
       "3F/GIouexwTUIsxnsJrsKX88FhW0CRNWWzFfQx0WKjuM7RPghpCXnFYKKkIjEEIIzDw0FjtNN2Cm\n",
       "mqa8OR3Y6XsFpjXshW3Aa8M42GY6LCoqB5BJJpZi5qIfY/b+dVj/hN9hCWWHx1PpccDKoWN2m9rR\n",
       "tpX6xsFt4ftbsQ1+a/jbH9vcmwYNGwXm4N0rrCHKam7GhNQEzHcAphlEIfFN4V4ufN6KaR4jevTE\n",
       "KhgJAiFEb3gN2xSHYxvqRrbfGCNH61rMDNOC9U/eiJ3EN2Ib7hq2P0kPwQTMCuxEn8CqGgzFtIv5\n",
       "wB+BI12sHkxL6AzXLCIbxvlP4Kmlj95+dsOwJta8/PRbW9Ys2w0zaW0Ic1rJaihDgV0W3nH1x4Hl\n",
       "8VR6MCbQOsMaN4TfvBA4tKvM6GqjIgSBfARCVB2rMTNPE1lBkK8RRD2Rl4e/F7FN/whs427B/ASF\n",
       "BMFRWKTQJmwTHpZJJj4brr0fC08/Fjhzj6PeuinM3w0z2YAJkPYwr3FL89pR7a0t4NyEF2773f8D\n",
       "WnK6nkXmn/nYqX+P1f958uzweVP4bWAO8hbM6b0W0xIOD53aYHtTUUmQj0AIURGEFpObsRj/tdhJ\n",
       "OV8QDGJbJ/IYLFLnFbLawwZsH9olmhRO2YdjEUkeM9scHT7ehAmIjwBfAubVNw5ej5l7vgZ8NMyd\n",
       "iG3YU4AD6+rr29s2byTWMHjjtLd86BdYbgM592wK63gSOKB965bos+GY6Yiwxq2Y4POYAItaf/Yb\n",
       "8hEIISqJ1VikzDpss3zdNBQKs8UIZR7C8GhgTShk144lby0J1w0l6wQ+BFgbfAVglVB3D+GkkWBZ\n",
       "iG32dwG/wCJ5NmBlL24E3olt2k8BZ+1+8Mm/bGvdwLBxk+8dOmbCGLYXBMPDepatXfTsEbvsfQhD\n",
       "x+x2x7rFz+2OaQQxTLC5vGv/DYwMoaVV23JXgkAI0VvWAFsyycSWoCG0ky0JPRrTFDbmjI0hW7Ay\n",
       "2tCXkbW9R7V/jgAeyPmeRZj2cUjOWAPZk/xizGfxZaycRUuY+0WsOmhj47CRsbqGwdQ3NLZiQifK\n",
       "e/BkI5tGA02ta15t8J0dDBm168kr5t93abjfMEzwjSLrZHYhTPYhLEJp1+4/uspCgkAI0VtWYk7g\n",
       "iA3YJgvZTX8j5viNhf9uyL1B2EhXYo7loWRzBl7ImbYY2/Snk92zdsf8ALFwfRsmDJ7C8hH+F4tk\n",
       "ejfWY3lW04QpAPtgZSCmh2J0EVE28piGwU24uhib1q96Yurx7/orZvrZFxMGI9lWIyCTTGzBMqP3\n",
       "zPn9VYUEgRCiL7TmvN5Ito7/GExjiDSCoZi5p1CByiWYkGjCcgieCbkKEeux0/gwLMrIYaGenZi2\n",
       "UIdpFlMw4TAPywO4A4s0Wgm8YdCIcYTxvbGM5q9jDud9MEEzEfBj9nzjbZ3tbRx2/sXfaBw28unw\n",
       "OxrD75jGtr0KgNfLd88v9Fk1UBGCQFFDQlQl+Zt6dPqvY9toosHYJp7fxyRiFbapR6fphXmfd2AO\n",
       "ZrCT+SDMlDMOK1MxIdx7Eqah3B/WsQ9WBfWTwLWd7VtWhmuj3ILISX00duo/AHi5vnFwI8CwXSa5\n",
       "MKcVK6L3PLbZH4kJhDqw/Svn96/s4jcWjaqLGnLODXPO/d45d4VzblZX8xQ1JMSAIDINDcWKtHVg\n",
       "TuRII1hT6KKgJSzFonyiZK18XsL2q+Fk7fTDsI25A9MGFgE+mGqWhnvdg5327x4yYuz9mON5AWZS\n",
       "OhFzVO+Hlcs+EQtBHTV0zATCesaQDRsdgQmpuzEfxfh4Kn30sHGTLw7tL5vI6WVQKqoxauidwJ+9\n",
       "9x8H3l7i7xJClJdmbMMfTtj0gw+gkx1rBGDmoQbg1Zz4fjAzkMe0hqjsc5TN/FqoORRpC8+RzSNY\n",
       "DWwMJpvBWFmJKJv4mfD+xjDPY+aig7Fy2Xt72xqnYyahKAFtFFZ2Ymu4x2vAy50d7WD+hb3Ycee0\n",
       "iqXUzetzO5d1IIQYsIQNPDpF557+N2HRRe0FL7Rr12Mb/asFPo5CNtdip/KoQ1nU82QNlqy2kawg\n",
       "WEfWTDME0xCawvr+gZmPxmBO6WfD3/Jw3e4TDjoB4INY+eqxmLAYDqwP5akPw0xTRwwZtSuYsGkm\n",
       "m3NQVZS6ef0Ssg1rKsIfIYQoKvmx81FLyNzT/2a6t0H+h2xoZj4bwt9IbDPfje0rgEZhnq8TehHX\n",
       "YU7k3bHNfyvwp3CvVixxbRBW7fQy4HHf3r4M27veBnwDOAOLCjoJK5OxAnNQ37p+8XNggii/UU3V\n",
       "UNLm9ViBqHc5535ODxrVCCGqllZga17Uzwr65kTtDH6E5djJfDSWs5AvNLYTBJhQ2hS0lagIXgtm\n",
       "rroLMwm9CRNod2BC4Zkho3b5BNZ68xVM8OyF9V8/D3gzcD6w39yLz7y8o20Lc2ef8fv1r764kbzw\n",
       "2Gqh1M3rW4EP7+ziPA/4PDmOhaha1rH9QTHagItx79VkT/X5fZULCYIhmEYC2UY5LZiZ5wGsXWYL\n",
       "1ucgt17SZuAJ4FZMk/gMJjwOxfa7GDD14HO+OumF235P6+qlp79w++/qj3h/6tEi/M6ChMjKGaW4\n",
       "d7EEQZ+a1xfbAy6EKBvrMJNJd+hpSYZIyAzCTunbmJsyycTWeCpNTl8ET9AIwvtNmCBoxgTGi2Sj\n",
       "lJaEsUgQDA7Xv4hpCu/BymAfgwmE64DPDh4x7rBJh8+kY+vmxZOOStyKZRiXhHBAnhe9d85dXKx7\n",
       "V0TzeuURCDHgyd/0PT0/QK7DNvNNmPkptz1mdK98rWAwWY2A3Dkh+mc9ZrpqwQ7GUULYUExobCbb\n",
       "6/gqrKjds8AHgOFDx0z47ZCxk5hw8MlfHDRs5DDMj1BSKjmPoE/N65VHIITYAR4gk0y0YCaZjWQL\n",
       "2eWTLwiGkNUIXM6caMNfjkUqtYe/RrLlMDqxyKB9CUlvoWPZTzBfwvCGIcMHNw4fxbAxE76GBcoM\n",
       "7tMv7QYVkUdQiub1QoiapDfVOtdj+1ZXNX0KCYJCGsHQ0EdgM9lKqEsxraABiygagvkD7iPHxxEE\n",
       "0i2YCezGRQ/cDCacppAtZldVqHm9EKJcdFcQ5JqQ1mH71qYQz59PIdPQptwJIYJoK9ue3qOS2IMx\n",
       "H8R0TOgsDFpAPhuAFZlk4o4VT94F8Fmy4a0lpZJNQ31CpiEhRDeJQkZfxHIJ8okEQSRkCmkE0bzI\n",
       "PNSARTatwExDEzBhcCvdN/U8gSWp9SlwpjtUhGlICCFKRFebaK7msAbbtF9ix4Iguq4uOIXz7587\n",
       "bxCWa7AVC08dgWkDS+H1Jjs7JOQ5PBqurzoqQhDINCTEgKHUJ+KtwPJQkqIunkpv0x4zZ9OPYfb+\n",
       "XLNOrkCJQkjBtIBIa3gJ0wJaQ/G63NLaO2MzXWdGFw2ZhoQQlU7JTSM5LMOSvfJpwU759RQ2C0Vz\n",
       "hmPrjbqdgfkJYliUEJjAiIRNf/62LpFpSAhRSXR3Y+xrL9+uInGWYeahQvWOIkGQ7+jNDSHdxjQE\n",
       "rxe/20w2KS5XI8j9nooQCsVCgkAIUSl0JTAKjmeSibXYZj8k76MWzNyzI42gFTMBObY1DYH5BqLa\n",
       "SLkawYBFgkAIUc0swxrH5LIjjQCAkJW8GTMLDcqbl1sXqZkqbT/ZEyQIhBClpq9mlB1dvwzL/s2l\n",
       "Oz6CaF4D22sErxNyCKJS1n1da8VSEYJAUUNCVC19tf/3iUwysYbsZh4RCYIYXWgEgaiNZgM7FhjN\n",
       "FC7QWei3l1wQVF3UkHNuqnPu1865v+xonqKGhBA7YGeb62vk2PFzQkgLnfTzcwmGYDWEdtRBcSPF\n",
       "q9TcZ6ouash7/5L3/qOl/A4hRFXQW83BdePa1Wxvx9+CNbLfWmB+RBQ5tCNtAEwjaNjJnKqmW4Kg\n",
       "h32KhRC1Tf6+UuhEX0wTSgfbC4st7LxPejO21p01zenKNJRPWc1kfaG7GkG3+xQ7597nnPuRc65Q\n",
       "oocQQvSF7lb33AK0FRh/XQAFR3Anlq28IzayvUZQlU7hruiWIOhJn2Lv/dXe+8977191zo1xzl0O\n",
       "HCKNQQjRj2zChMHOaGPngqAV2yur9sS/M/riIyjUp3hi7gTv/Rrv/Se89/t477/bh+8SQgw8elOG\n",
       "urtsoHsF4DazE4ERCspBhURZloK+eMKLphqpeb0QNUl3BUEpT+LLsIqmvfm+ftUQKrV5fZ/6FBdA\n",
       "AkAIUYid+QUK9T8upbbRFSUVDFHz+lIIhL6oOn3qU5yL8giEEGWiJ4KgIhzEZcsjUJ9iIYQYuHTL\n",
       "NOS9P7eL8Vuxdm59IvgIZBoSQvSG/rLVV0TUUClMQxWRNl1sNUcI0S8U21RSaKPt7VhPNu38ud35\n",
       "XWUzE+X4Ci4u1j0HbDiUEKIs9CWLuKI34G7gqOz1dYkEgRCimJRiI6zKzbWaqAjTkHwEQgjRPeQj\n",
       "EEJUIz1qQVnBVMR65SMQQlQDpd4wi33/sjSYqSQkCIQQpaaUrSr7el1Xc3Z2bW+vq0gkCIQQosaR\n",
       "IBBCiBqnpM5i59yZQAIYAfzGe39bF/PmoKghIaqRinCg1hJVFzXkvb8RuNE5Nwq4FCgoCBQ1JIQo\n",
       "Mt3NSC40XqhPcnds//1iYSlb1FARehZ/A2trKXZAkPQCPYtcBuiz6IlT9fUy1Hek3n1CGdcxYClp\n",
       "z2JnfBe41Xv/eFFXPjCZUe4FVBAzyr2ACmJGuRdQKXR2bD2x3GsYiHS3+ug9zrkpecOv9ywGcM5F\n",
       "PYu/A1wdxj4DnAKMcM7t7b3/ZZHWLYQQokj0xUdQqGfx0bkTvPeXAZf14TuEEAOX3jqa5aAuMs77\n",
       "7pnIgkZws/f+oPD+XcBp3vuPhffnA0d77z/dowU4JxudEEL0Au99UYRi2XsWF+uHCCGE6B0V0bNY\n",
       "CCFE+VDPYiGEqHG6JQi89+d673f33g/y3k/23l8Zxm/13u/rvd/be//tnn55D/IQqpJC+RfOuTHO\n",
       "uduccwucc3NDsl302dfCs3jOOTczZ/xw59xT4bOf9PfvKAbOucnOubucc884554OEWU1+Tycc4Od\n",
       "cw865x53zs13zn07jNfcs4hwzsWcc485524O72vyWTjnXnbOPRmexUNhrPTPwntflj8gBiwEpgAN\n",
       "wOPA9HKtp0S/8QTgUOCpnLHvAV8Jry8CvhNe7x+eQUN4JgvJOvMfAo4Kr2/BnPRl/309fBYTgEPC\n",
       "6+HA88D0Gn4eQ8N/64EHgONr9VmEtX8B+CNwU3hfk88CeAkYkzdW8mdRzqJzr+cheO/bgOuAM8u4\n",
       "nqLjvb8HWJs3/Hbg9+H174GzwuszgWu9923ecjMWAkc753YDmrz3D4V5V+VcUzV475f7kFTovW8G\n",
       "nsVCkGv1ebSGl43YoWgtNfosnHOTgLcCvyYbGlqTzyKQH0BT8mdRTkFQKA9hYpnW0p+M996vCK9X\n",
       "AOPD693ZNuoqeh7540up8ucUQpEPBR6kRp+Hc67OOfc49pvv8t4/Q40+C+BHwJfJKSVB7T4LD9zu\n",
       "nHvEOfexMFbyZ1HOVpU1nz/gvfe1lkfhnBsO/BX4rPd+o3PZw08tPQ/vfSdwiHNuJJBxzp2U93lN\n",
       "PAvn3BnASu/9Y66Lmkq18iwCx3nvlznnxgG3Oeeey/2wVM+inBpBUfIQqpAVzrkJAEGFWxnG85/H\n",
       "JOx5LA2vc8eX9sM6i45zrgETAld7728IwzX7PAC89+uBNHA4tfks3gS83Tn3EnAtcLJz7mpq81ng\n",
       "vV8W/rsK+BtmQi/5syinIKjVPISbgA+E1x8AbsgZP8c51+icmwrsAzzkvV8ObHDOHe3s+Py+nGuq\n",
       "hrD23wDzvfc/zvmo5p6Hc26XKPLDOTcEeAvwGDX4LLz3X/cWiTgVOAe403v/PmrwWTjnhjrnmsLr\n",
       "YcBM4Cn641mU2UN+OhY9shD4WjnXUqLfdy3wKrAV84d8CBgD3A4sAOYCo3Lmfz08i+eAeM744eEf\n",
       "xELgsnL/rl4+i+MxG/Dj2Kb3GFa5tuaeB3AQ8Gh4Fk8CXw7jNfcs8p7LiWSjhmruWQBTw7+Jx4Gn\n",
       "oz2xP55Ft2sNCSGEGJioZ7EQQtQ4EgRCCFHjSBAIIUSNI0EghBA1jgSBEELUOBIEQghR40gQCCFE\n",
       "jSNBIIQQNc7/B1hNVCYPhO3oAAAAAElFTkSuQmCC\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bc271d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    np.random.seed(1235)\n",
    "    data = (xval, yval, np.sqrt(yval))\n",
    "    print(\"Fitting GP\")\n",
    "    truth = [1.,1.,fitps[0], fitps[1], fitps[2], fitps[3]]\n",
    "    sampler = fit_gp(truth, data) \n",
    "    print(\"Making plots\")\n",
    "    samples = sampler.flatchain\n",
    "    pl.figure()\n",
    "    pl.errorbar(xval, yval, yerr=np.sqrt(yval), fmt=\".k\", capsize=0)\n",
    "    for s in samples[np.random.randint(len(samples), size=24)]:\n",
    "        a, tau = np.exp(s[:2])\n",
    "        gp = george.GP(a * kernels.Matern32Kernel(tau))\n",
    "        gp.compute(xval, np.sqrt(yval))\n",
    "        m = gp.sample_conditional(yval - model_gp(s,xval),x) + model_ind(s[2:],x)\n",
    "        pl.plot(x, m, color=\"#4682b4\", alpha=0.3)\n",
    "    pl.yscale('log')\n",
    "    #pl.savefig('CDF2009.png')"
   ]
  },
  {
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}