gistfile1.json

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from scipy import optimize\n",
      "from astropy.modeling.fitting import Fitter\n",
      "\n",
      "# TODO: This is currently needed to avoid an error in the contructor\n",
      "from astropy.modeling.fitting import constraintsdef\n",
      "constraintsdef['CashFitter'] = ['fixed']\n",
      "\n",
      "def cash(D, M):\n",
      "    \"\"\"Cash Poisson likelihood statistic.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    D : array-like\n",
      "        \"data\", i.e. observed counts per bin\n",
      "    M : array-like\n",
      "        \"model\", i.e. expected counts per bin\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    cash : array\n",
      "        Cash statistic value per bin\n",
      "    \"\"\"\n",
      "    D = np.asanyarray(D, dtype=np.float64)\n",
      "    M = np.asanyarray(M, dtype=np.float64)\n",
      "\n",
      "    stat = 2 * (M - D * log(M))\n",
      "    stat = np.where(M > 0, stat, 0)\n",
      "\n",
      "    return stat\n",
      "\n",
      "\n",
      "class CashFitter(Fitter):\n",
      "    \"\"\"Cash Poisson likelihood fitter.\n",
      "    \n",
      "    Calls the `scipy.optimize.minimize` optimization function.\n",
      "    \"\"\"\n",
      "    \n",
      "    def __init__(self, model):\n",
      "        Fitter.__init__(self, model)\n",
      "                        \n",
      "    def errorfunc(self, fitpars, *args):\n",
      "        \"\"\"The Cash Poisson likelihood fit statistic.\n",
      "        \n",
      "        TODO: give formula and reference.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        fitpars : TODO\n",
      "            TODO\n",
      "        *args : (y, x)\n",
      "            Tuple with y counts at coordinate x\n",
      "        \n",
      "        Returns\n",
      "        -------\n",
      "        stat : float\n",
      "            Cash Poisson likelihood fit statistic\n",
      "        \"\"\"\n",
      "        self.fitpars = fitpars\n",
      "        D = args[0]\n",
      "        M = self.model(*args[1:])\n",
      "        stat = cash(D, M)\n",
      "        return stat.sum()\n",
      "    \n",
      "    def __call__(self, x, y):\n",
      "        \"\"\"Execute the likelihood minimization.\n",
      "        \n",
      "        TODO: document\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        x : array-like\n",
      "            x-coordinate\n",
      "        y : array-like\n",
      "            Observed number of counts at ``x``\n",
      "        \"\"\"\n",
      "        result = optimize.minimize(self.errorfunc,\n",
      "                                   x0=self.model.parameters[:],\n",
      "                                   args=(y, x))\n",
      "        self.fitpars = result.x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from astropy.modeling.models import Gaussian1DModel\n",
      "model = Gaussian1DModel(amplitude=10, mean=2, stddev=3)\n",
      "print('True parameters: ', model.parameters)\n",
      "x = np.arange(-10, 20, 0.1)\n",
      "y = np.random.poisson(model(x))\n",
      "\n",
      "model = Gaussian1DModel(amplitude=7, mean=1, stddev=2)\n",
      "print('Parameters before fit: ', model.parameters)\n",
      "fitter = CashFitter(model)\n",
      "fitter(x, y)\n",
      "print('Parameters after fit: ', model.parameters)\n",
      "y_model = model(x)\n",
      "\n",
      "plt.plot(x, y, 'o')\n",
      "plt.plot(x, y_model, 'r-');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "('True parameters: ', [10.0, 2.0, 3.0])\n",
        "('Parameters before fit: ', [7.0, 1.0, 2.0])\n",
        "('Parameters after fit: ', [10.574787489518142, 1.8945598462935331, 2.9540072831120914])"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "-c:27: RuntimeWarning: divide by zero encountered in log\n",
        "-c:27: RuntimeWarning: invalid value encountered in multiply\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEACAYAAABF+UbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VPX9//FnSAKJhL2QCQQhhlUIgSKIVHQsZVFkKfJF\nUEQJ/Fwqi2hdcPkSggjWUisB3FoUa7+uVaEEEZROQATRAoILgiFBwMQlEHYSknx+fwyZzCR3JnPv\nLHdm8n6cMyfJnXs/9z3b+9zcufd1o5RSCiGEEGGpgdkFCCGEME6auBBChDFp4kIIEcakiQshRBiT\nJi6EEGFMmrgQQoQxj008IyODxMRE0tLSHNO2b99O//796dOnD/369eOzzz4LeJFCCCG0eWziU6ZM\nYd26dS7THnjgAebPn8/OnTvJysrigQceCGiBQggh3PPYxAcNGkSLFi1cpiUlJXH8+HEASkpKaNeu\nXeCqE0II4VFUXWdsFhQUMHLkSPbs2QPAwYMHufLKK4mKiqKyspKtW7fSvn37oBQrhBDCle4vNqdO\nncqSJUv4/vvvefrpp8nIyAhEXUIIIbyh6pCfn6969uzp+LtJkyaO3ysrK1XTpk01l0tNTVWA3OQm\nN7nJTcctNTW1rrbsQveWeKdOncjNzQVg48aNdOnSRXO+vLw8lFIRe5s7d67pNYTD4xs69BHN9+qw\nYY9GxOMLxVskP7b68Pjy8vJ09eQYT3dOnDiR3NxcfvnlF9q3b09WVhYvvPACd999N6WlpcTHx/PC\nCy/oWqGoX2bOHEpe3iPk5S1wTEtNfZgZM4abWJUQkcNjE3/ttdc0p3/66acBKUZEnhEjrgIgO/sx\nzp2LJi6ughkzhjumCyF847GJC/esVqvZJQSUPx/fiBFXhVzTjuTXL5IfG0T+49OrzkMMDQ8cFUWA\nhhZCiIilt3dKdooQQoQxaeJCCBHGpIkLIUQYkyYuhBBhTJq4EEKEMWniQggRxqSJCyFEGJMmLoQQ\nYUyauBBChDFp4kIIEcakiQshRBiTJi6EEGFMmrgQQoQxaeJCCBHGPDbxjIwMEhMTSUtLc5menZ1N\n9+7d6dmzJw8++GBACxRCCOGex4tCTJkyhRkzZjB58mTHtP/85z+sXr2a3bt3Exsby88//xzwIkV4\nyMnZxJIl6yktjaFRo3JmzhwacheDECLSeGzigwYNoqCgwGXas88+y5w5c4iNjQWgdevWAStOhI+c\nnE3MmvWBy7U08/IeAZBGLkQA6d4nvn//fjZt2sSAAQOwWq18/vnngahLhJklS9a7NHCAvLwFZGdv\nMKkiIeoH3dfYLC8v59ixY2zbto3PPvuM8ePHc+DAAc15MzMzHb9brVa5Nl4EKy3VfiudOxcd5EqE\nCC82mw2bzWZ4ed1NPDk5mbFjxwLQr18/GjRoQHFxMa1atao1r3MTF5GtUaNyzelxcRVBrkSI8FJz\nA3fevHm6lte9O2XMmDFs3LgRgH379lFWVqbZwEX9MnPmUFJTH3GZlpr6MDNmDDGpIiHqB49b4hMn\nTiQ3N5fi4mLat29PVlYWGRkZZGRkkJaWRsOGDXnllVeCVasIYVVfXmZnP8a5c9HExVUwY8Zw+VJT\niACLUkqpgAwcFUWAhhZCiIilt3fKGZtCCBHGpIkLIUQYkyYuhBBhTJq4EEKEMWniQggRxnSf7COE\nr9wFZUmAlhD6SRMXQeUuKOuzz77k1VePSICWEDrJceIiqIYNe5T16x+vNb1VqxspLn5DY/7HWLdu\nfjBKEyIkyHHiIqS5C8oqL4/XnC4BWkJ4Jk1cBJW7oKyYmLOa0yVASwjPpImLoHIXlDV9+tUSoCWE\nAbJPXARdTs4msrM3OAVlDXEcnaI1XYj6RG/vlCYuhBAhRL7YFEKIekSauBBChDFp4kIIEcY8NvGM\njAwSExNJS0urdd/ixYtp0KABR48eDVhxQgghPPN42v2UKVOYMWMGkydPdpl+6NAhNmzYQIcOHQJa\nnNDPH/kjWmMAkmsiRAjy2MQHDRpEQUFBren33nsvf/rTnxg9enSg6hIGuMslAe/zR7TG2L17KtCM\noqK/GB5XCBEYuveJr1q1iuTkZHr16hWIeoQPlixZ79J8AfLyFpCdvcGnMYqKklwauJFxhRCBoSvF\n8MyZMzzxxBNs2FD94fV0PGNmZqbjd6vVitVq1V2g8J67XBI9+SPaY/g+rhBCm81mw2azGV5eVxPP\ny8ujoKCA9PR0AA4fPkzfvn3Zvn07bdq0qTW/cxMXgecul0RP/oj2GL6PK4TQVnMDd968ebqW17U7\nJS0tjR9//JH8/Hzy8/NJTk5mx44dmg1cBJ+7XBI9+SNaY1gsP2Cx3OvTuEKIwPC4JT5x4kRyc3Mp\nLi6mffv2ZGVlMWXKFMf9UVFRAS9QeK/qS8bs7Mec8keG6/ryUXuM23weVwgRGJKdIoQQIUSyU4QQ\noh6RJi6EEGFMmrgQQoQxaeJCCBHGpIkLIUQY03WyjxAQOgFZ/gj7EiLcSRMXuoRKQJY/wr6EiASy\nO0XoEioBWf4I+xIiEkgTF7qESkCWP8K+hIgE0sSFLqESkOWPsC8hIoE0caFLqARk+SPsS4hIINkp\nQrecnE1kZ29wCsOyN86a04JxdEqw1ylEoOntndLEhRAihEgAlhBC1CPSxIUQIoxJExdCiDBWZxPP\nyMggMTGRtLQ0x7T777+f7t27k56eztixYzl+/HhAixRCCKGtzi82N2/eTEJCApMnT2bPnj0AbNiw\ngcGDB9OgQQMeeughABYtWuQ6sHyxGVLMyBnxxzolH0XUN3p7Z53ZKYMGDaKgoMBl2pAh1cfiXn75\n5fzrX//yvkIRdGbkjPhjnZKPIkTdfN4nvmLFCq677jp/1CICxIycEX+sU/JRhKibTymGCxYsoGHD\nhtx0002a92dmZjp+t1qtWK1WX1YnDDIjZ8Qf65R8FFEf2Gw2bDab4eUNN/GXX36ZtWvX8tFHH7md\nx7mJC/OYkTPij3VKPoqoD2pu4M6bN0/X8oZ2p6xbt46nnnqKVatWERcXZ2QIEURm5Iz4Y52SjyJE\n3eo8OmXixInk5ubyyy+/kJiYyLx581i4cCFlZWW0bNkSgCuuuILly5e7DixHp4QUM3JG/LFOyUcR\n9Y1kpwghRBiT7BQhhKhHpIkLIUQYkyYuhBBhTJq4EEKEMZ9O9hGRxUhOibfLaM0HeL0+oxkqkr0i\nIp4KkAAOLQJgzZpclZr6sALluKWmPqzWrMn1eRmt+SyWDGWxzPZqfUZq82U5Icykt3dKExdKKaWG\nDn3EpdlV3YYNe9TnZbTn8359RmrzZTkhzKS3d8o+cQEYyynxdhnt+bxfn9EMFcleEfWBNHEBGMsp\n8XYZ7fm8X5/RDBXJXhH1gTRxARjLKfF2Ga35LJYfsFju9Wp9RjNUJHtF1Ady2r1wMJJT4u0yWvMB\nXq/PaIaKZK+IcCPZKUIIEcb8fnk2IQzbtw/efhu2bYMjR6BBA7j4YvjNb+B//gfatze7QiHCnjRx\n4X9798KDD9qb9/jxcOut9uZdWQkHDsBHH0Hv3jBiBDzxBCQnm12xEGFLdqcI/1EKnn0W5s6FOXPg\nrrsgPl573uPH4c9/ts+/eLG90Qsh/LtPPCMjg5ycHNq0acOePXsAOHr0KDfeeCMHDx6kY8eOvPnm\nmzRv3tznQkSYO3sW7rgDvvgC3nkHUlO9W27PHvvW+ogR8NRTEBUV2DqFCHF+zROfMmUK69atc5m2\naNEihgwZwr59+xg8eDCLFi0yVqmIHOfOwejR9p9bt3rfwAHS0mDLFvtt6lQo1z62Wwihrc7dKQUF\nBYwcOdKxJd6tWzdyc3NJTEykqKgIq9XK3r17aw8sW+JhSU9gVE7OJpb/9X3u3/Z/lJQpMjuPIbFd\nU2NBVqdPw9ixYLHAyy973CJ3N4aEXYlIoLt31nVefn5+vurZs6fj7+bNmzt+r6ysdPnbl/P/hfn0\nBEatWZOrUi+Zo15hiHqXLiqGMt+DrE6dUuqyy5TKytJd49y5yyTsSkQEvb3TpyaulFItWrTwSyHC\nfHoCo4YOfUTN5K9qB4kqntP+C7IqLFSqQwel/vlPXTW2ajVewq5ERNDbO3UfYli1G8VisVBYWEib\nNm3czpuZmen43Wq1YrVa9a5OBJGewKhuPx7mYV5kABM4y0VeLePV+BYLrF4Nv/0t9O0LXbt6NUZ5\nufZRMBJ2JUKdzWbDZrMZXl53Ex81ahQrV67kwQcfZOXKlYwZM8btvM5NXIQ+rwOjjh7l0b3vcAtv\nUsDH3i2jZ/xevSArCyZOtH9R2qhRnWPExJz1ug4hQknNDdx58+bpWt7j0SkTJ05k4MCBfPvtt7Rv\n356XXnqJhx56iA0bNtClSxc2btzIQw89ZKhwEXq8DoyaOZNTQ4fwXepmYCjgXciUrkCqu+6ynyD0\niOv87saYPv1qCbsS9ZKc7CNc1BkY9e678MADsGsXObb/kp29gcOHf6Ko6DhJSRbatWvivyCrX36x\nb5W/8w4MGFDnGBJ2JSKBBGCJwDl2DHr0gLfesuefBMNrr8HChfDf/0JsbHDWKYSJpImLwJk+3Z5/\nsnx58NapFAwfDoMH2/8DECLCSRMXgbFrFwwbBt98Ay1bBnfdBw5A//72rfEOHYK7biGCzK+n3QsB\n2LeG774bHn88+A0c4JJL7P8FzJkT/HULEeKkiYu6vf46lJZCRoZ5Ndx/P+Tm2uNthRAOsjtFeFZW\nBt26wYoVYLU68kmOHPmZoqISkpKSaNs2wW1OiZE8k5ycTTz22CsUFJxCqUakpDRm/vwJJG9YS+zL\nr/KH9Kk0iqvwKRtFclZEqPJ7dopRARxaBNOSJUpde61Syjm3JFdB3TklerJYnJexWDJqjd+8+STV\nNnGW2kFvNY43fcpGMVKXEMGit3dKExfuHT+uVGKiUrt2KaWcc0u8y1jRk8XiuozWcvZpw3hffUV3\n1YByw9koRuoSIlj09k7ZJy7ce/ppGDIE0tMB59wS7zJW9GSxuC6jtZx92gcM4ygtmchrdY7leR36\n6hIiVMk1NoW2EycgO9ueXXJBdW6JdxkoXmeleLFM9TqjeJTH+RvTeIMbDWWjGKlLiFAlW+JC29Kl\n9pNsOnd2TKrOLfEuL0VXVorTMhZLYa3xmzcvwGK5F4BcrBTQkT+2HmMoG8VIXUKEKjk6RdR26pT9\n2OzcXOje3eWuqnwSb/NSjOSZVB+dchpoSEpKAllZNwI4xup7Np+s7z+k8aGDEKP/H0rJWRGhSs7Y\nFL576in4/HN44w2zK/Fs0CD4wx/skbVCRAhp4sI3Z87Yt8I3bLBfxDiUrV1rP4tz1y6P1+QUIpzI\naffCNy+8AAMHhn4DB7j2WnvzXrvW7EqEMI1siYtqZWX2rfBVq+yXRgsHb7wBS5bAxx/L1riICEHb\nEl+4cCE9evQgLS2Nm266idLSUqNDiVDxxhv2U+zDpYEDjBsHP/0EmzebXYkQpjDUxAsKCnjxxRfZ\nsWMHe/bsoaKigtdff93ftYlgUgr+/Ge47z6zK9EnOtqeM75wodmVCGEKQyf7NG3alNjYWM6cOUN0\ndDRnzpyhXbt2/q6tXgh0EJPW+ECtgKnlN3RiQHm5/dhwqg/z27evkLNnS4F44uOb0KVLc+bPn+D3\nGrUCrzwFalUFcCUkxFB2spRdJzbz6MDbGf3IJAC/PKcSkiXCgtHz+59//nmVkJCgWrdurSZNmuTz\n+f/1UaCDmLTGt1gyVPPmk2oFTG1s2EHtmvWgYzl7CFXVreYYs/1ao1bgldY6agdwVQdxPUqWepGp\nymLJUBbLbJ+fUwnJEmbR2zsNddrvvvtOde/eXf3yyy/q/PnzasyYMerVV1/1qZD6KNBBTNrj1w6Y\n6sEedYQkdf2Qh5yWe0Rz3sDUqDdQq3YQ16/4SR2luWrNPX6pV0KyhFn09k5Du1M+//xzBg4cSKtW\nrQAYO3Ysn3zyCTfffLPLfJmZmY7frVYrVqvV0H8LkSrQQUza49eedi9/YSnTOVlW5rEuZ4GtUXsd\ntQO4qpf9hda8zTjuZBfzvRjLaF0SkiX8zWazYbPZDC9vqIl369aN+fPnc/bsWeLi4vjwww/p379/\nrfmcm7ioLdBBTNrju05LpIgxvEcnvqN/3NMe63IW2Bq111E7gMt12b9yDx9yOU9SShmNPI5ltC4J\nyRL+VnMDd968ebqWN3R0Snp6OpMnT+ayyy6jV69eANx+++1GhqrXAh3EpDW+xfIDzZsfpCpg6nZe\n4E3G08gy37He6hCqqlvNMWb7tUatwCutddQO4HIN4vqaHnzbsBl3Nh3lspyR51RCskS4kJN9TBbo\nICat8cF+dMrh/JPsOv5v7uk2iluf+oPLemsfnXIR8fEJdO3agqysGwN0dIpr4JWnQK2qAK6EhGhO\nnapwBHFlDWxKt5eeZ3yXCZwrjfHpOZWQLGEGyU4R3nvzTVi+HHzYHxdylIIePWDZMrjmGrOrEUI3\nyU4R3lu2DKZPN7sK/4qKghkz7HnoQtQDsiVeX+3eDdddB/n5EBtrdjX+dfIkdOgAX3wB7dubXY0Q\nusiWuPDOsmVwxx2R18ABmjSBSZPg+efNrkSIgJMt8fqopARSUuCbb8BiMbuawNi7F6xWOHgQGjWq\nc3YhQoXe3ikXSg4Dfs/wePll+66USG3gYE9jTEuDt9+Gm292eQ5PnDgMNKSsrIKiohKSkpJo2zaB\nK65oy9atP0hWiggr0sRDXE7OJmbN+oC8vAWOaXl59uOXDTWYykr7rpSVK/1VYui6+27405/Iad7e\n6TncBHwADLvw83mKi+HLLzexceP/UV7+nGNxn55nIYJE9omHuCVL1rs0cIC8vAVkZ28wNuD69fZ9\nxldc4YfqQtz118ORI6x9fKXTc7geWOD0E8d05wYOPj7PQgSJNPEQ5/cMj6rDCuvDVXBiYuCuuxhR\nsMN5Yo2fNae7kqwUEeqkiYc4v2Z4HDgA27bBhAk+VhVGpk5l0C97aUnxhQnauSu1/7aTrBQR6qSJ\nhzi/Zng8+yzcdhtcdJF/igsHrVtz/Kqrub/l+AsTtHNXYCgxMXe6LCpZKSIcyCGGYcAvGR5nzthP\ngPn0U/vFkOuT7ds5M2oM49KncKY0lhMnDhMV1YjS0nKKio47clcGDEhi27ZCyUoRppLsFKFtxQp4\n5x1Ys8bsSszRvz/87//av+wUIoTJGZuiNqUgOzvyclL0uPtuyVMREUmaeH2wdSucOgVDh5pdiXlu\nvBF27ID9+82uRAi/kiZeHyxbBn/4AzSoxy93XBxkZNi/3BUighjeJ15SUsK0adP46quviIqKYsWK\nFQwYMKB6YNknHhqKiqB7d/vhhS1amF2NuQoKoG9f+P57aNzY7GqE0BS07JRZs2Zx3XXX8fbbb1Ne\nXs7p06eNDiV0qsoBOXLkZ4qKSkhIiOHUqXJHBsjMmfbdJkuWrGf8vi1YElLgkz2MGHGV/3NY6qjR\nzBwSzRoGDYJ//hMuXE5QT53VVyA6hVKNSElpzPz5E+QIFmEuZUBJSYlKSUnxOI/BoUUd1qzJVamp\nDyvIVeD8UzluFkuGslhmqxjK1CHaqTS+UKmpD6u5c5ddWLZ63tTUh9WaNbkBqjGw6zFSw7b5f1aq\nVy+lKit11blmTa6yWDI0nuvZQX1cIvLp7Z2GOu3OnTtV//791W233ab69Omjpk2bpk6fPu1TIcI7\nQ4c+cqGB1PzpfLNPu4G31CaudExv1Wq8xrxKDRv2aIBqDOx6jNQwfOgjSnXtqtSmTbrqtM9r/uMS\nkU9v7zS0O6W8vJwdO3awdOlS+vXrxz333MOiRYvIyspymS8zM9Pxu9VqxWq1Gv2HQVxQnaXiLgOk\netp0lpLNDMfU8vJ4zTH9nQ/i97wXP9ZwtjTGcbhhaWl3zXm06nQ3nrv5hfCWzWbD5sN1bg018eTk\nZJKTk+nXrx8A48aNY9GiRbXmc27iwj+qs1TcZYDYp/VkD53Zz7v83jE1Juas5pj+zgfxa95LIGqY\nPBnmziUpTTtPXatOd+O5m18Ib9XcwJ03b56u5Q0dc2axWGjfvj379u0D4MMPP6RHjx5GhhI6VWep\nuMsAAYvlB+6Pv5nnuJNy7JdfS019mOnTr/ZfDotXNQZ2PYZraNYMJkxgXrsTXtc5c+ZQLJZCaj/X\nsyVfRZjK8CGGX3zxBdOmTaOsrIzU1FReeuklmjVrVj2wHGIYMFVZKocP/0RR0XESEqI5darCkQEy\ne8oAfjt1ApP7TaNQNXPJAfFLDouOGs3MIfFYw1dfwZAhrF3+D5Y8Z/OqzuqjU04DDUlJSSAr60Y5\nOkX4lWSnCHj6afj8c/uhdMK9a66BO++0n80pRIiQ7JT6rurya/U5J8Vb06dLnooIe9LEI80HH9j3\n+TqdPSvcGD3afhbnrl1mVyKEYdLEI83SpfXn8mu+iomBO+6w/+ciRJiSfeKR5LvvYOBAOHgQ4rWP\nCRc1/PgjdOsm2TIiZMg+8frs2WftSX3SwL2XmAgjRsBLL5ldiRCGGA7AEv7jl7Co06dh5Ur7USk6\n1lsV6FRaepbo6Gg6dLjYEaJl9NA5vaFSzvNecUVbtm79wWVZwKfnp656tvQeSMqjj3LzquM0jKvQ\nrEEOIxQhy4+n/LsI4NARxW9hUS+8oNTo0brWWx3oVDtEy2hgld5QKdd5c1VMzB2aYV5Ga6urnjVr\nclXqJXPU5/xaDWetZg3BDu8S9Zve3ilN3GR+CYuqrFSqRw+lNmzQuV5PIVrGgp30h0rVDu6qe5r3\ntdVVT9X9t7FCrWW4hFwJ0+ntnbJP3GR+CYtavx6io2HwYJ3r9RSiZSzYSc/jqT2v+zAvb8YzUk/V\n/a8xkT7spDvHfFqfEMEmTdxkfgmL+stf4N57dR1WaF+vpxAtY8FOeh5P7Xm1w7y8Hc9IPVX3lxLH\ncv7AbD7xaX1CBJs0cZP5HBa1Z4/9NnGi7vVWBzrVDtEyGlil5/HUnncoMTF3usxjsfyAxXKv4drq\nqsf5/me5i3F8R9voyYbXJ0SwyXHiIcCnsKiMDOjUCR5+2NB6qwKdysrO0KBBDB06tKdduyY+BVbp\neTw15x0wIIlt2wpdlgV8CtOqqx7n++/P+zeNL0lkUXxfU8O7RP0lAVj1SVERXHop7N8PrVqZXU1k\n2LsXrr7afjq+HG8vTCAn+9Qny5bZd6NIA/efbt2gf3/4xz/MrkQIr8iWeLg6cwY6doQtW6BzZ7Or\niSw2G9x1lz1zvIFs54jgki3x+uKVV+w5KdLA/e/qq+27Ut5/3+xKhKiTT028oqKCPn36MHLkSH/V\nI7xRWWm/8MO999Y9r9AvKgruu89+6KYQIc6n7JRnnnmGSy+9lJMnT/qrnnpLK8ekZct4iouPARcR\nE9OYlJTGzJ8/gRFlxfbM8EGD3I6llUdy5MjPFBWVkJSU5HM+SqjxJoOl6vJ0XuWwjB8PDz3E5mde\nZNbKTykoOIVSjapfgwvL+CX3RghfGD019NChQ2rw4MFq48aN6vrrr/f51NH6TDvHJFdB1TSnLJHE\ne9SxTl2Vevddt2Np55H4Lx8l1HiTwZKa+rCaO3eZrpyaL//fDLWmUYfar4FltlqzJtd/uTdCONHb\nOw132nHjxqkdO3Yom80mTdxH2jkmztOqb0P4QOUntFaqosLDWFrZI5GbCeJdBotSrVqN1/UcjBz8\ngCqkserOV5rL+CX3Roga9PZOQ7tT1qxZQ5s2bejTpw82m83tfJmZmY7frVYrVqvVyOoinmu+h+c8\nk0dYwKvJV/Kom6Mm3OeR+C8fJdR4l8EC5eXax327ew5OlMfzDJczh4VM5h9eLVPXfULUZLPZPPbR\nuhhq4p988gmrV69m7dq1nDt3jhMnTjB58mReeeUVl/mcm7hwzzXfw32eyW/4mGQO80n7K70cy/N4\nEBmZIN5lsEBMzFnN6e6eg0aNyllOH/J4mRQOkM8lLssoN4eBRcJzKoKn5gbuvHnz9A3g66a/7E7x\nnbf7xNcyXP2x6WCP+1z17ROfExH7b73bJz7HzT5x989B1euSxW/Uc9zutE/8Hg/7xCPjORXm0ds7\nfT7ZJzc3l8WLF7N69WqX6XKyjz5aOSYtW8ZRXFwCxNM36hxvnLOx4833uHbM7+ocSyuP5PDhnygq\nOk5SksXnfJRQ400GS9XRKXpyWHJyNrF4zgv868u3GNR0JI1SE8nKutHl6BRfcl2EqEmyUyLV6NHw\n29/CrFlmV1I/3XcflJVBdrbZlYgIJ008En36KYwbZw+6ioszu5r66aefoHt32LkTLr7Y7GpEBJPT\n7iPRo4/CY49JAzdTmzZw552QlWV2JUK4kC3xUGezwdSp9ojU2Fizq6nfjh2DLl3soWNduphdjYhQ\nsiUeSZSyb4VnZkoDDwUtWsDs2TB3rtmVCOEgTTyUrVsHR4/CTTeZXYmoMnMm/Oc/sHu32ZUIAcju\nFL/xexBSZSVcdpn9smvjxrkdv2r6kSM/c/DgIWqFZYXQ4W7BCIsyug5dy/31r7BxI9Q4rNb5tfAm\naEzCs4QW3b3TD8emawrg0CEnIEFIL72k1IABSlVWuh2/+uQVN2FZF4KaQkEwwqKMrkP3cmfPKpWS\notTGjRpjeBc0JuFZwh29vVOauB/4PQjp5Eml2rZVautWj+NXBzpph2WFUhhTMMKijK7D0HJvvaVU\nerpS5eU1xvBuLAnPEu7o7Z2yT9wPagcw2RkOQnrySbjmGhgwwOP41YFOMYR6wJXfnyM/rsPQcjfc\nAE2bwooVNcbwbqxgPB+ifvDpohDCrnYAk52hIKTvv4fly2HXrjrHrw500r7fcA0B4NfnyM/rMLRc\nVJT96krXXw/jxzuN4d1YwXg+RP0gW+J+MHPmUFJTH3GZlpr6MDNmDNE/2EMPwfTp0L59neNPn371\nhelDgULAdR6LZbaxGgLAr8+Rn9dhuLa+feHaa2HBAqcxhlLzddAaKxjPh6gf5OgUP/FLENLmzTBx\nInz7LTRu7NX4VdMPH/6JgwcPA/HExjYmJSXBJagpFAQjLMroOgzXVlgIPXvC1q3k7C9yvBbeBI1J\neJbQItmkiwA5AAAQ6klEQVQp4aq0FHr3hscft+9vFeFj8WLIyYGPPrLvZhHCB3LGZrhatAi6doWx\nY82uROg1axacOAEvv2x2JaIeki3xUPDNN/Yr1+/aBcnJZlcjjNi5E4YNgz17IDHR7GpEGJPdKeGm\nshKsVhg/3v6FpghfDzwAhw7Ba6+ZXYkIY0HbnXLo0CGuueYaevToQc+ePVmyZInRoeq3F1+0X2zg\nrrvMrkT4KjPTnv2ek2N2JaIeMbwlXlRURFFREb179+bUqVP07duX9957j+7du9sHrgdb4tWXVDuF\nUo1ISWnMqFE92br1B0pLYzhx4jDQkKZN22hnY3z3HVxxhT1utkePWmPXzNUAaq0v1PJRIoWnXBN3\nr82SJetp9+0uHj+0kevaXkdeyRlqZtlUzVe17BVXtHW8X7TeI97mq0gOS+QwLTtl9OjR6sMPP3T8\n7cehQ5LrxY2rL9AbHX2743ePGRrnz9uzUf76V82xa+ZqWCwZqnnzSSGdjxIpPOWauHttLJbZjtd8\nMf+j3qaDgjku8zVvPunCfJ4u6Fz9HvE2X0VyWCKL3t7pl06bn5+vLr74YnXy5EnDhYQbe/ZFzfyL\nR9z8rpGNkZWl1O9+p1RFhZuxay4b+vkokcJTron716b6ZyMeVLtprW5jhYf3R93vEW/zVSSHJbLo\n7Z0+n3Z/6tQpxo0bxzPPPENCQoLLfZmZmY7frVYrVqvV19WFDO3sixg3v1c7dy7avvtk2TL4/HNo\nUPtribrH1hhT+I3+XBPX3JRS4riZG/iIB/iUy/mGS2vMV3M57fV4W4fksIQ3m82GzWYzvLxPTfz8\n+fPccMMNTJo0iTFjxtS637mJRxrt7ItyN79XS4o6DjffDK+84vZwwrrHdiV5G/7lKddEae6rrJmb\nUs4eWvEgT/I24+jPdk6TQO3X0HN+irf5KpLDEt5qbuDOmzdP1/KGj05RSjF16lQuvfRS7rnnHqPD\nhK2ZM4disdTMKxlKdPQdjt9rZmh0veRBlhbbYNo0GDrU49g1czUslh9o3vxgrTFDKR8lUnjKNXH3\n2lgs91L9mtuzbF4ij20M4AVuBxTNmxdcmK/KUGJi7tRcT111eFuviHyGj075+OOPueqqq+jVqxdR\nF041XrhwIcOHD7cPXK+OTjkNNCQlJYGRI3uwbVsh585Fc+LEYaKiGtGkSWvi4ip4LuZLOlacgzVr\nINrzv7pauRpArfWFWj5KpPCUa+LutXHOTUlIiKa4uIR4FcuGc5uxJXal0wuLHPNVLTtgQJLj/aKV\nn+JtvorksEQOOdknVC1bZo+Y/eQTaNbM7GpEMB0+bM+GX7oUNHY7CuFMmngoysmBqVPtDfySS8yu\nRpjhs8/guuvg/fft104Vwg0JwAo1H38Mt90G770nDbw+69cP/vY3GDkS9u41uxoRQeTKPoH0xRf2\nWNl//tNxqTVRj40eDceO2YOyNm+Giy82uyIRAaSJB8qOHfZ/n5ct83gkiqhnbrsNjh+3h5599BGk\npJhdkQhzpuwTD9WcB+cslNLSs0RHR9OyZTzFxceoKwPD5TF89pn92ovPPQe//32t8fftK7ywXGO6\ndGku+Sdhqup9fOTIzxQVlZCUlETbtglu80+c57ujfC8Zxbu4aMtm6NzZMZ9zNk7Lludo1qyVZvaO\n0c+QZLGEPtOyU2pyN3So5jy4ZqHkOv2smY+ilYHh9BjWrlXqV79SavVqN+PXHk/yT8JP9fvYc0aO\np/nmtL5OnW3ZSqnt2zWyeNyPa/QzJFks4UFvWw56Ew/VnAfXLBTnn56yMlxvT196nVKJiUpt2eJh\n/NB8/EKf6vext/kn2vPN7TNeqV/9Ss3tPU55m71j9DMkWSzhQW8TD/o+8VDNeXCtK6bGz5pcpzek\nlKeZzfX5n8KuT6BLlzrGr83sxy/0qX49vc0/0Z7P1rQ7rP0jd185GMXPzKcSRYM6x9V7n2stnpcL\n1c+o0Bb0QwxDNefBXlftDAztfIvqaR3JJ5ersVDEPQMzNBu46/ih+fiFPtXvY2/zTzzM168fM6+Y\nwmAOkMMIfsXPHuc3+hmSLJbIFPQmHqo5D65ZKK4ZGDXzSpo3L8CSOJsprGA7/XmT8TxwSVemzr7e\ni/Frjyf5J+Gn+n1cOyNHO//E83yT7r+BmxMH8wXH2U0vRnGR2/mNfoYkiyUymXZ0SijmPDhnoZSV\nnaFBgxhatoyjuLgEiCc2tjEpKQn8+c7f0OO5Zzix/yALuo6lsFUrrx6D1tEpXbu2kPyTMFX1Pq7K\nS0lKstCuXRO3+SfezPfYY6+Q+F0By09/ylfxLVieYuVMq0u8ym/x9ugUyWIJbXLafSCVlMBTT8Hz\nz8OcOTBzJsTGml2ViERnzsCTT9rzVmbPhj/+EeLizK5KBIGcdh8IJ0/CggX243mLimDnTrjvPmng\nInAuugjmzbNfOGTHDvt3LUuXwtmzZlcmQow0cU8OHrRvcaemwtdfw5Yt8Pe/Q/v2Zlcm6ouUFHjn\nHXjzTdiwwZ6/86c/QXGx2ZWJECFNvKZz5+xhVWPGwK9/bf97yxZ7/ombI0+ECLgBA2DVKvjgA/jy\nS/uGxY03wvr1UCFHjdRnhpv4unXr6NatG507d+bJJ5/0Z03Bd+yYfWtn0iRISoJnnoERI+D77+Hp\npx2nRQthul697Jf2KyiAq6+2/6fYtq39alH//rfsbqmHDH2xWVFRQdeuXfnwww9p164d/fr147XX\nXqN79+7VA4fyF5uFhfZ8k23b7CFE33wDAwfCqFEwdixYLHUOYbPZIurCzzXJ4wsjBw7Yt9JXr4bP\nPsOWmop11Ci46iro2xdatjS7Qr+KqNdOg97eaeiMze3bt9OpUyc6duwIwIQJE1i1apVLEweIjR1N\ndHQFDRqUoVQDysrKqayspEGDxsTGKsf08vJGmvM5q7mM1jzO88VEldJGnadD2Uk6q9N04zTdOEMv\nThNPBZ/Tgv9GNyM3uhnbogZwakM5bHif+Ac+JjGxolbwEOA4PPDs2VIqK7+ncePLXAKsnIOODh48\nBFxEZeV5lDpLdHTTWgFazmFHVdND5TCuSP+g6H18vgZCOS9/4sRhjh8/ybFjUY7XftSonmzd+oPm\n/S1bngPg2LEoSkvPOt5PVe+t6s9QY1rGXkPr7/Yw6cnXuGLBMnqpExwnht005gsS2BvVkkMx8fzQ\nsBllrVrwy9ESyssra32eqj5HDRs2dITAac3nPK+3n03fP+f5NGjQM8jrDPzjhHji45vUGq8uhpr4\nkSNHaO/05V5ycjKffvpprfnKy++jvHxljakWKiuHUVq60vE3DNOYD6Kp5CLKiaclF1UOJL70DeKJ\npjHlNKeCFpTTklLHrQ3QrjKa5NICLJylhIbspyl7acZekrDRjC9pTj6dgeFQsRIqcKpjAadPb+LA\ngQ+ABY46du+eyrlzZZSUNLwwnwW4ktOnM9m5E6ZNu5c77viSV189Ql7eMGAl0Aeo+t3iGG/nTpg0\n6RYASkouBv7mmD5t2r387W+ETCMXdjk5m5g16wPy8qrfE3l59pNhvD02u3r5TdjfE52ofk9sYvfu\nf1JR8bzG/ZsoKal6Dzm/n6p+h+rP0AcUlg6jkFk8wK8v3KfoyCnSOUo6ZVyrYul4/ltSzp+l2enj\nHKIxPxLPT8TxE/H8TBw/0Yajlb05Wfotp0oHcvLkf2jERZwihpPEco5oztMAiKL259kdb+bzdqzz\nVFbeH+R1BvpxVvUfgNfqWIcrQ0286sLIdcllEtEooql0+plINC86TWtONMuc/lbEUc5FnKcBijPE\ncpZmnGElZ4nhLDGcIZYS4jhKPEeJ4yjN+IY4NjOFw7zPEaz8QAJlbh/e48CjQFKNaQDrcW7gAEVF\nSbh6HMh0uv8vLF16I8XFbziN+3iN36uVlHSosc7qcbKzH5MmHmKWLFnv0sAB8vIWeP1auS6/ntrv\nifUXGrjW/c5/a723cPp7wYWfTXB+bxdcuK1yzPcR8CjxnKcDx2nDadpwhtYXfl5KKi15hwRSaMLT\nJFBGE8ocP+MoJ5ZKztOAMuIo4y8Xfo+m7EKDLyOaCqKoJApFFJUko3jJ6e8oFNT4+xIq+T+X+7V2\nKrxKFJOYUufzrugKvFHHXF3czuO6bm/Gqp5P4alHdgXe1Jh2Pc8wiw+9WIsLIylbW7duVcOGDXP8\n/cQTT6hFixbVSOJKrXoN5CY3uclNbl7fUnX1Y0NfbJaXl9O1a1c++ugj2rZtS//+/Wt9sSmEECLw\nDO1OiYmJYenSpQwbNoyKigqmTp0qDVwIIUwQsOwUIYQQgef3MzbfeustevToQXR0NDt27HC5b+HC\nhXTu3Jlu3bqxfv16f6866DIzM0lOTqZPnz706dOHdevWmV2SzyLqJC4NHTt2pFevXvTp04f+/fub\nXY7PMjIySExMJC0tzTHt6NGjDBkyhC5dujB06FBKSkpMrNA3Wo8vUj53hw4d4pprrqFHjx707NmT\nJUuWAAZePyNfbHryzTffqG+//VZZrVb13//+1zH9q6++Uunp6aqsrEzl5+er1NRUVVFR4e/VB1Vm\nZqZavHix2WX4TXl5uUpNTVX5+fmqrKxMpaenq6+//trssvyqY8eOqri42Owy/GbTpk1qx44dqmfP\nno5p999/v3ryySeVUkotWrRIPfjgg2aV5zOtxxcpn7vCwkK1c+dOpZRSJ0+eVF26dFFff/217tfP\n71vi3bp1o4tGxsiqVauYOHEisbGxdOzYkU6dOrF9+3Z/rz7oVATtjXI+iSs2NtZxElekiaTXbNCg\nQbRo0cJl2urVq7n11lsBuPXWW3nvvffMKM0vtB4fRMZraLFY6N27NwAJCQl0796dI0eO6H79ghaA\n9cMPP5CcnOz4Ozk5mSNHjgRr9QGTnZ1Neno6U6dODet/W0H7JK5IeI2cRUVF8bvf/Y7LLruMF198\n0exyAuLHH38kMTERgMTERH788UeTK/K/SPrcARQUFLBz504uv/xy3a+foSY+ZMgQ0tLSat3+/e9/\n6xrH25OGzOTusa5evZq77rqL/Px8du3aRVJSEvfdd5/Z5fokHF4PX23ZsoWdO3fy/vvvs2zZMjZv\n3mx2SQEVFRUVca9rpH3uTp06xQ033MAzzzxDkyaup9178/oZOsRww4YNupdp164dhw4dcvx9+PBh\n2rVrZ2T1QeXtY502bRojR44McDWBVfM1OnTokMt/T5EgKcl+JmPr1q35/e9/z/bt2xk0aJDJVflX\nYmIiRUVFWCwWCgsLadOmjdkl+ZXz4wn3z9358+e54YYbuOWWWxgzZgyg//UL6O4U5/1Wo0aN4vXX\nX6esrIz8/Hz2798f9kcHFBYWOn5/9913Xb5BD0eXXXYZ+/fvp6CggLKyMt544w1GjRpldll+c+bM\nGU6ePAnA6dOnWb9+fdi/ZlpGjRrFypX2XI6VK1c6mkOkiJTPnVKKqVOncumll3LPPfc4put+/fz9\njes777yjkpOTVVxcnEpMTFTDhw933LdgwQKVmpqqunbtqtatW+fvVQfdLbfcotLS0lSvXr3U6NGj\nVVFRkdkl+Wzt2rWqS5cuKjU1VT3xxBNml+NXBw4cUOnp6So9PV316NEjIh7fhAkTVFJSkoqNjVXJ\nyclqxYoVqri4WA0ePFh17txZDRkyRB07dszsMg2r+fj+/ve/R8znbvPmzSoqKkqlp6er3r17q969\ne6v3339f9+snJ/sIIUQYk8uzCSFEGJMmLoQQYUyauBBChDFp4kIIEcakiQshRBiTJi6EEGFMmrgQ\nQoQxaeJCCBHG/j85YfjoMJTsyAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c084bd0>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}