09_global_south.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:e75a98dd920bad796171b3657e4b1b625a5225bcd391dab3bb958501b5d056e6"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import HTML\n",
      "import folium\n",
      "import random as RD\n",
      "import pandas as pd\n",
      "import pysal as ps\n",
      "%pylab inline\n"
     ],
     "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": [
      "def inline_map(map):\n",
      "    \"\"\"\n",
      "    Embeds the HTML source of the map directly into the IPython notebook.\n",
      "    \n",
      "    This method will not work if the map depends on any files (json data). Also this uses\n",
      "    the HTML5 srcdoc attribute, which may not be supported in all browsers.\n",
      "    \"\"\"\n",
      "    map._build_map()\n",
      "    return HTML('<iframe srcdoc=\"{srcdoc}\" style=\"width: 100%; height: 510px; border: none\"></iframe>'.format(srcdoc=map.HTML.replace('\"', '&quot;')))\n",
      "\n",
      "def embed_map(map, path=\"map.html\"):\n",
      "    \"\"\"\n",
      "    Embeds a linked iframe to the map into the IPython notebook.\n",
      "    \n",
      "    Note: this method will not capture the source of the map into the notebook.\n",
      "    This method should work for all maps (as long as they use relative urls).\n",
      "    \"\"\"\n",
      "    map.create_map(path=path)\n",
      "    return HTML('<iframe src=\"files/{path}\" style=\"width: 100%; height: 510px; border: none\"></iframe>'.format(path=path))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Binding Attribute Data to the Map"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import json\n",
      "\n",
      "f = open(r'data/south.json')\n",
      "q = json.load(f)\n",
      "f.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q.keys()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "[u'type', u'features']"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "features = q['features']\n",
      "len(features)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "1412"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "feature_0 = features[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "feature_0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "{u'geometry': {u'coordinates': [[[-80.6280517578125, 40.39815902709961],\n",
        "    [-80.60203552246094, 40.480472564697266],\n",
        "    [-80.62545776367188, 40.504398345947266],\n",
        "    [-80.6336441040039, 40.53913879394531],\n",
        "    [-80.6688232421875, 40.568214416503906],\n",
        "    [-80.66793060302734, 40.58207321166992],\n",
        "    [-80.63754272460938, 40.61391830444336],\n",
        "    [-80.61175537109375, 40.619998931884766],\n",
        "    [-80.57462310791016, 40.615909576416016],\n",
        "    [-80.52220916748047, 40.63713836669922],\n",
        "    [-80.52456665039062, 40.47871780395508],\n",
        "    [-80.52377319335938, 40.4029655456543],\n",
        "    [-80.6280517578125, 40.39815902709961]]],\n",
        "  u'type': u'Polygon'},\n",
        " u'properties': {u'BLK60': 3.839454752,\n",
        "  u'BLK70': 3.2554278095,\n",
        "  u'BLK80': 2.5607402642,\n",
        "  u'BLK90': 2.5572616581,\n",
        "  u'CNTY_FIPS': u'029',\n",
        "  u'COFIPS': 29.0,\n",
        "  u'DNL60': 6.1681225056,\n",
        "  u'DNL70': 6.1714993547,\n",
        "  u'DNL80': 6.1714631077,\n",
        "  u'DNL90': 6.0508978146,\n",
        "  u'DV60': 2.2779893943,\n",
        "  u'DV70': 2.5591397849,\n",
        "  u'DV80': 5.0619350519,\n",
        "  u'DV90': 7.2636377003,\n",
        "  u'FH60': 9.9812973718,\n",
        "  u'FH70': 7.8,\n",
        "  u'FH80': 9.7857968181,\n",
        "  u'FH90': 12.604551644,\n",
        "  u'FIPS': u'54029',\n",
        "  u'FIPSNO': 54029.0,\n",
        "  u'FP59': 9.6,\n",
        "  u'FP69': 5.9,\n",
        "  u'FP79': 6.5327526442,\n",
        "  u'FP89': 10.17311807,\n",
        "  u'GI59': 0.2236450331,\n",
        "  u'GI69': 0.2953773833,\n",
        "  u'GI79': 0.3322512119,\n",
        "  u'GI89': 0.3639335641,\n",
        "  u'HC60': 0.6666666667,\n",
        "  u'HC70': 1.6666666667,\n",
        "  u'HC80': 2.6666666667,\n",
        "  u'HC90': 0.3333333333,\n",
        "  u'HR60': 1.6828642349,\n",
        "  u'HR70': 4.1929776011,\n",
        "  u'HR80': 6.5977204876,\n",
        "  u'HR90': 0.9460827444,\n",
        "  u'MA60': 28.9,\n",
        "  u'MA70': 30.0,\n",
        "  u'MA80': 31.4,\n",
        "  u'MA90': 37.7,\n",
        "  u'MFIL59': 8.8410143105,\n",
        "  u'MFIL69': 9.2471543451,\n",
        "  u'MFIL79': 10.073356901,\n",
        "  u'MFIL89': 10.327970666,\n",
        "  u'NAME': u'Hancock',\n",
        "  u'PO60': 39615,\n",
        "  u'PO70': 39749,\n",
        "  u'PO80': 40418,\n",
        "  u'PO90': 35233,\n",
        "  u'POL60': 10.586963113,\n",
        "  u'POL70': 10.590339963,\n",
        "  u'POL80': 10.607030509,\n",
        "  u'POL90': 10.469738422,\n",
        "  u'PS60': 1.218684208,\n",
        "  u'PS70': 1.1368342185,\n",
        "  u'PS80': 1.0385705291,\n",
        "  u'PS90': 0.8964534429,\n",
        "  u'RD60': -1.394676863,\n",
        "  u'RD70': -1.307438562,\n",
        "  u'RD80': -1.159302086,\n",
        "  u'RD90': -0.399028376,\n",
        "  u'SOUTH': 1,\n",
        "  u'STATE_FIPS': u'54',\n",
        "  u'STATE_NAME': u'West Virginia',\n",
        "  u'STFIPS': 54.0,\n",
        "  u'UE60': 3.1,\n",
        "  u'UE70': 2.7,\n",
        "  u'UE80': 7.0763827919,\n",
        "  u'UE90': 6.8578070515},\n",
        " u'type': u'Feature'}"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pandas as pd\n",
      "indices = []\n",
      "values = []\n",
      "for feature in features:\n",
      "    indices.append(str(feature['properties']['FIPS']))\n",
      "    values.append(feature['properties']['HR80'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df = pd.DataFrame({'HR80': values,\n",
      "                   'FIPS': indices} )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "map_osm = folium.Map(location=[33.7550, -87.3900], zoom_start=5,)\n",
      "map_osm.geo_json(geo_path='data/south.json',\n",
      "                 key_on='feature.properties.FIPS',\n",
      "                 data_out='data.json', data=df,\n",
      "                 columns=['FIPS','HR80'],\n",
      "                 fill_color='YlGn', fill_opacity=0.7,\n",
      "                 line_opacity=0.2,\n",
      "                 legend_name='Homicide Rate HR80')\n",
      "embed_map(map_osm)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe src=\"files/map.html\" style=\"width: 100%; height: 510px; border: none\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "<IPython.core.display.HTML at 0x10d011190>"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = np.array(df.HR80.tolist())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins= ps.Quantiles(y, 5).bins.tolist()\n",
      "map_osm = folium.Map(location=[33.7550, -87.3900], zoom_start=5,)\n",
      "map_osm.geo_json(geo_path='data/south.json',\n",
      "                 key_on='feature.properties.FIPS',\n",
      "                 data_out='data.json', data=df,\n",
      "                 columns=['FIPS','HR80'],\n",
      "                 fill_color='YlGn', fill_opacity=0.7,\n",
      "                 line_opacity=0.2,\n",
      "                 threshold_scale=bins[:-1],\n",
      "                 legend_name='Homicide Rate HR80 (Quantiles k=5)')\n",
      "embed_map(map_osm)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe src=\"files/map.html\" style=\"width: 100%; height: 510px; border: none\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "<IPython.core.display.HTML at 0x10d035e50>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Global Spatial Autocorrelation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w = ps.queen_from_shapefile('data/south.shp')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w.n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "1412"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w.transform = 'r'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.random.seed(12345)\n",
      "mi_hr80 = ps.Moran(y, w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80.p_sim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "0.001"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80.I"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "0.2521152652175852"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80.EI"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "-0.0007087172218284905"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80.p_norm"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "0.0"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": [
      "Create a random map"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "yr = np.random.permutation(y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df['yrandom'] = yr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "type(y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 25,
       "text": [
        "numpy.ndarray"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins= ps.Quantiles(y, 5).bins.tolist()\n",
      "map_osm = folium.Map(location=[33.7550, -87.3900], zoom_start=5,)\n",
      "map_osm.geo_json(geo_path='data/south.json',\n",
      "                 key_on='feature.properties.FIPS',\n",
      "                 data_out='data.json', data=df,\n",
      "                 columns=['FIPS','HR80'],\n",
      "                 fill_color='YlGn', fill_opacity=0.7,\n",
      "                 line_opacity=0.2,\n",
      "                 threshold_scale=bins[:-1],\n",
      "                 legend_name='Homicide Rate HR80 (Quantiles k=5)')\n",
      "embed_map(map_osm)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe src=\"files/map.html\" style=\"width: 100%; height: 510px; border: none\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "<IPython.core.display.HTML at 0x10d0541d0>"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins= ps.Quantiles(y, 5).bins.tolist()\n",
      "map_osm = folium.Map(location=[33.7550, -87.3900], zoom_start=5,)\n",
      "map_osm.geo_json(geo_path='data/south.json',\n",
      "                 key_on='feature.properties.FIPS',\n",
      "                 data_out='data.json', data=df,\n",
      "                 columns=['FIPS','yrandom'],\n",
      "                 fill_color='YlGn', fill_opacity=0.7,\n",
      "                 line_opacity=0.2,\n",
      "                 threshold_scale=bins[:-1],\n",
      "                 legend_name='Spatially Random Homicide Rate HR80 (Quantiles k=5)')\n",
      "embed_map(map_osm)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe src=\"files/map.html\" style=\"width: 100%; height: 510px; border: none\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 27,
       "text": [
        "<IPython.core.display.HTML at 0x10cfeaa50>"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80_random = ps.Moran(df['yrandom'],w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80_random.I"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "-0.0075100395407744691"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80_random.p_sim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 30,
       "text": [
        "0.34799999999999998"
       ]
      }
     ],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import gaussian_kde"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "density = gaussian_kde(mi_hr80.sim)\n",
      "xs = linspace(mi_hr80.sim.min(),mi_hr80.sim.max(),200)\n",
      "plot(xs, density(xs))\n",
      "axvline(x=mi_hr80.I, color='r')\n",
      "xlabel('I')\n",
      "ylabel('f(I)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 32,
       "text": [
        "<matplotlib.text.Text at 0x10d6ea350>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEPCAYAAACtCNj2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGVZJREFUeJzt3XuQXFdh5/FvSxrZ1gO90cuyRi9LMviFN0YGjAeXTcxW\n7JhUQdY4lQAbClIbWB7FI1CUZ1lXCghOCJUqNgGbMhQYKgWbsr3WFibJBGGEib16vzWSLFmPka3I\nwpKNLdu9f5xuaTSaR/d033vO7f5+qqbU09O378/N0L8559x7GyRJkiRJkiRJkiRJkiRJkqRkLQD+\nFdgCbAY+Vrm/G3gaWFf5uiVGOElSvuYAV1VuTwJ2ACuBu4BPxgolSRrauAyf+0jlC+AksA2YX/m+\nlOF+JUmJ6wSeIowY7gL2ARuAe4Gp0VJJknI3CXgCuL3y/esJI4UScDehGCRJCch6GqcDeBhYDXx9\nkJ93Ag8Bl/e/c8mSJeXe3t6Mo0lSy+kFljbyBGOaFGQwJcIoYCvnFsLcfrffDWwauGFvby/lcjn5\nr7vuuit6hurX1q1lpkwpM2NGmR070s1ZlNfTnOYcNCNEzzLcF7Ck0TfuLBea3wr8EbCRcOgpwOeB\nOwhHJZWBvcCHM8zQNv7qr+Czn4WxY+GLX4Qf/Sh2IklFlGUp/ILBRyKrM9xnWyqX4ZFHYO1amDkT\nFi+G3l5Y0vDfDJLaTZbTRy2vq6srdgQANm2CyZNh0aLw7/vfD3//92d/nkrOkZizuczZPEXI2Cyp\nni9QrsyPqQZ/8zewcyd885vh+x074O1vhwMHYPz4uNmkllIqhaF5okqlEjT4vu5IoQWsWwfXXnv2\n++XLYcUKePDBeJkkFZOl0AI2b4Y3vvHc+z70IfjWt+LkkVRcTh8V3KuvhnWEZ56BiRPP3v/ii7Bg\nATzxBHR2RosntRanj5S63l6YO/fcQgC46CJ43/vgXs8Xl1QHS6HgtmyBN7xh8J995CPwD/8AJ07k\nm0lScVkKBTfc+QiXXQa33gpf+lK+mSQVl6VQcHv3hvMThvKXfwnf/35YW5CkkVgKBTdSKbz+9XDP\nPfDBD8LLL+eXS1IxWQoFN1IpQFhwXrAgXB9JkobjIakFVi7DhAnhcNRJk4Z/7O7d8Ja3hLOcL7gg\nn3xSy/GQVKXsyJFwjsJIhQCwdGk4wc2znCUNx1IosP374ZJLan/8nXfCT36SXR5JxWcpFNjBgzB/\nfu2Pv/FG+Ld/S3r0KykyS6HADh2qrxQ6O6GjA3btyiySpIKzFArs4EGYN6/2x5dK4ZLaa9Zkl0lS\nsVkKBVbvSAHgyivDh/JI0mAshQKrd6QAcPnl4VLbkjQYS6HARjNSeOMbLQVJQ7MUCuzQofpHCvPm\nwUsvhRPeJGkgS6GgTp0Kb+7TptW3XakUrp66dWs2uSQVm6VQUNVRQmkUJ7QvXhyumSRJA1kKBVXv\niWv9LVpkKUganKVQUKNZT6iyFCQNxVIoqMOHw2czj4alIGkolkJBHT0Ks2ePbltLQdJQLIWCOno0\nfKraaFx8cTgk9be/bW4mScVnKRRUI6UwdizMmROmoCSpP0uhoBopBQiL1IcONS+PpNZgKRSUpSAp\nC5ZCAZXLoRRmzRr9c1gKkgZjKRTQyZMwZgxMnDj655g/31KQdD5LoYCeeaaxqSNwpCBpcJZCATW6\nngCWgqTBWQoFZClIykqWpbAA+FdgC7AZ+Fjl/unAo8BO4KfA1AwztCRLQVJWsiyF08AngDcAq4D/\nBqwEPkcohUuBf658rzo0oxSmTIHTp8OitSRVZVkKR4D1ldsngW3AfOA24P7K/fcDt2eYoSU1oxRK\npTBa8KxmSf3ltabQCVwNPA7MBvoq9/dVvlcdmlEK4BSSpPONy2Efk4AfA/8deH7Az8qVr/N0d3ef\nud3V1UVXV1c26QrIUpAE0NPTQ09PT1OfcxQf5liXDuBhYDXw9cp924EuwvTSXMJi9IoB25XL5UG7\nQsAVV8B3vwtXXdXY83z849DZGf6VVINSKVxSIFGl8Pm8Db2vZzl9VALuBbZythAAHgT+pHL7T4B/\nyjBDS2rWSGH2bDhypPHnkdQ6siyFtwJ/BLwDWFf5ugX4MnAz4ZDUGyvfq0avvQbHjsHMmY0/15w5\nloKkc2W5pvALhi6dmzLcb0s7fhwmT4bx4xt/rjlzoK9v5MdJah+e0VwwzZo6AkcKks5nKRRMM0vB\nNQVJA1kKBdPMUpg1C559Fl59tTnPJ6n4LIWCaWYpdHTAtGlh4VqSwFIonGaWAriuIOlclkLBNLsU\nXFeQ1J+lUDCOFCRlyVIomKNHwwJxs1gKkvqzFAqmGZ/P3J8nsEnqz1IoGNcUJGXJUiiQl1+G558P\nh5E2i9NHkvqzFArk2WfDhfDGNPF/NUtBUn+WQoE0e+oIXFOQdC5LoUCyKIXp0+E3vwlTU5JkKRRI\nFqUwZkw4xPXo0eY+r6RishQKJItSANcVJJ1lKRRIlqXguoIksBQKpdlnM1d5roKkKkuhQJp9NnOV\n00eSqiyFAnFNQVLWLIUCcU1BUtYshQLJqhRcU5BUZSkUxKlTUC7DxInNf26njyRVWQoFUR0llErN\nf25LQVKVpVAQWU0dAbzudXD6NLzwQjbPL6k4LIWCyLIUSiUXmyUFlkJBZFkK4GKzpMBSKIiszmau\ncl1BElgKhZHV2cxVloIksBQKI+vpI9cUJIGlUBiuKUjKg6VQEHmMFCwFSZZCQVgKkvJgKRTAa6+F\nheasjz5yTUGSpVAAzz0HkybB+PHZ7aO6plAuZ7cPSenLuhTuA/qATf3u6waeBtZVvm7JOEPhZT11\nBDBhAnR0wG9+k+1+JKUt61L4Due/6ZeBvwaurnz934wzFF4epQCuK0jKvhTWAMcHuT+Da322rqzX\nE6pcV5AUa03ho8AG4F5gaqQMhZHXSMFzFSTFKIVvAouAq4DDwD0RMhSK00eS8jIuwj6P9rv9beCh\nwR7U3d195nZXVxddXV2ZhkrZ0aOwcmX2+5k3Dw4ezH4/kpqjp6eHnp6epj5nHnP7nYQ3/ssr388l\njBAAPgH8DvC+AduUyx4becZ73hO+3vvebPfzgx/Agw/CD3+Y7X6kwiqVkj5uuxQ+mrGh9/WsRwoP\nADcAM4EDwF1AF2HqqAzsBT6ccYbCy2v6qLMT9u3Lfj+S0pV1KdwxyH33ZbzPlmMpSMqLZzQXQJ4L\nzc89By++mP2+JKXJUkjc6dPhLOPp07Pf15gxcMklsH9/9vuSlCZLIXHPPgszZoQ37DwsXOgUktTO\nLIXE5XU2c5XrClJ7sxQSl9d6QlVnJzz1VH77k5QWSyFxfX3h8hN5caQgtTdLIXExRgp79uS3P0lp\nsRQS19eXbylceins2pXf/iSlxVJI3NGj+U4fzZwZPv7z2LH89ikpHZZC4vIeKZRKsGwZ7NyZ3z4l\npcNSSFzeIwUIU0iWgtSeLIXE5b3QDJaC1M4shYSVy/lPH4GLzVI7q+UqqVOB6wifi1AG9gFrgROZ\npRIAzz8PHR0wYUK++3VNQWpfw40UrgceBH4O/BfgEkIx3AGsqfzsbRnna2t5n7hWtWxZGCm89lr+\n+5YU13AjhXcDnwKGmki4FPgI8Itmh1IQYz0B4HWvgylT4Omnw1VTJbWP4UYKn2ToQgDYWXmMMhLj\nyKOqlSth+/Y4+5YUz3AjhU8N87My8NdNzqIBYiwyV61YAdu2wTvfGWf/kuIYrhQmE978FUms6SMI\npbB1a5x9S4pnuFLoziuEBtfXF6ZxYli5En7ykzj7lhTPcGsK3cBwM9pzgf/R1DQ6R+yRwrZtcfYt\nKZ7hRgpPAD8ExgP/DzgMlIA5wJuAl4CvZR2wncU6JBVg/nw4dQqOH4dp0+JkkJS/4UYKfwi8A1hN\nOOz0VeB05fYfAjcCj2QdsJ3FHCmUSmG04BFIUnsZbqRwDTAPeC/QRRglVLkAnYOYh6TC2VK47rp4\nGSTla7hS+F/APwOLgScH/KxcuV8ZefllOHkSpk6Nl8GRgtR+hps++gawEvgOsGjAl4WQsWeegVmz\nYEzESxauXOlis9RuannL+UjmKXSemCeuVTlSkNqPl85OVOz1BIClS2H/fnjppbg5JOXHUkhUCiOF\n8eNh4ULYvTtuDkn5sRQSdfgwzJ0bO4VTSFK7sRQSZSlIisFSSNShQzBvXuwUHoEktRtLIVGOFCTF\nYCkkKpWRQrUU/GhOqT1YCgkql9MZKUydCpMnw8GDsZNIyoOlkKDnnguHg06YEDtJ4BSS1D6yLoX7\ngD5gU7/7pgOPEj7j+adAxKv7pOnw4TSmjqpcbJbaR9al8B3glgH3fY5QCpcSLrj3uYwzFM6hQ2lM\nHVU5UpDaR9alsAY4PuC+24D7K7fvB27POEPhpDZS8FPYpPYRY01hNmFKicq/ka/wk57URgpLl8Ke\nPbFTSMpD7IXmMn5gz3lSGyksWABHjoTPeJDU2ob7kJ2s9BE+5/kIMBc4OtiDuru7z9zu6uqiq6sr\nh2hpOHQIVq2KneKsjo7wmc1PPQXLlsVOI6mqp6eHnp6epj5naeSHNKwTeAi4vPL9V4FjwFcIi8xT\nOX+xuVwut+8A4vrr4e674YYbYic566ab4NOfht/93dhJpIhKpXAiUaJKpRI0+L6e9fTRA8AvgeXA\nAeADwJeBmwmHpN5Y+V79pHI2c3+LF7uuILWDrKeP7hji/psy3m9hpXQ2c3+WgtQeYi80a4DnnoNx\n42DSpNhJzmUpSO3BUkjMgQNwySWxU5xv8WLYuzd2CklZsxQSs39/uqXQ25v0GpukJrAUEpNqKUyb\nFg68OD7w/HRJLcVSSEyqpVAqua4gtQNLITH794cziFNkKUitz1JITKoLzWApSO3AUkhMqtNHYClI\n7cBSSMgrr4QT1+bPj51kcJaC1PoshYQcPgyzZoWP4kxR9bBUSa3LUkhIylNHELIdOgSnT8dOIikr\nlkJCUi+F8ePDNZkOHIidRFJWLIWEpF4K4LqC1OoshYSkfDhqlaUgtTZLISEpn7hWZSlIrc1SSMie\nPbBoUewUw7MUpNZmKSSiXA6Xpl68OHaS4VkKUmuzFBLR1wcTJ8LkybGTDM9SkFqbpZCIPXvSHyUA\nzJgRzrz2EtpSa7IUElGUUqheQttPYZNak6WQiN7eYpQCOIUktTJLIRF79sCSJbFT1MZSkFqXpZCI\nokwfgaUgtTJLIRGWgqQUWAoJePFF+I//gHnzYiepjaUgtS5LIQF798LChTB2bOwktVm4MFyn6ZVX\nYieR1GyWQgKKNHUEcMEFMHs2PP107CSSms1SSECRDketcgpJak2WQgJ27oRLL42doj6LF8Pu3bFT\nSGo2SyEBRSyF5ctDbkmtxVJIwK5dxSuFFStg+/bYKSQ1m6UQ2YsvwpEj6X/i2kCWgtSaLIXIenvD\nB+uMGxc7SX0WLw5HH730UuwkkprJUoisiFNHAB0d0NnpYrPUaiyFyHbuhGXLYqcYHaeQpNZjKURW\nxCOPqpYvhx07YqeQ1EwxS2EfsBFYB/w6Yo6oijp9BI4UpFYUsxTKQBdwNXBtxBxROX0kKSWxp49K\nkfcf1YkTcPJkca6OOlB1+qhcjp1EUrPEHin8DHgC+FDEHNHs2hVGCaWCVuP06XDRRXDwYOwkkpol\n5tHxbwUOA7OAR4HtwJrqD7u7u888sKuri66urnzT5aBaCkV2xRWwYQNcfHHsJFL76enpoaenp6nP\nmcrfqHcBJ4F7Kt+Xy20wJ9HdHT6T4O67YycZvU9/GqZNg89/PnYSKQelUtLzpaUw7dDQ+3qs6aMJ\nwOTK7YnAO4FNkbJEs20brFwZO0VjrrwSNm6MnUJSs8QqhdmEqaL1wOPAw8BPI2WJZutWuOyy2Cka\nU50+ktQaUpk+Gqjlp49eeQUmT4Zjx2DChNhpRu/ll2Hq1PDfcdFFsdNIGXP6SFnp7Q2Hoha5EADG\njw8n323eHDuJpGawFCJphamjqiuvdApJahWWQiStsMhcZSlIrcNSiKSVRgpXXQXr1sVOIakZLIVI\nWqkUrrkG1q+H06djJ5HUKEshgldfDReSa5XpoylTwseJbtkSO4mkRlkKETz1FMyYEQ5JbRXXXguP\nPx47haRGWQoRbNgQFmdbyZvfDL9u20/FkFqHpRDB+vVhcbaVOFKQWoOlEMG6dXD11bFTNNcVV8De\nvfD887GTSGqEpRBBK44UOjrClNiTT8ZOIqkRlkLOjh0Ln7i2aFHsJM23ahU89ljsFJIaYSnkbP36\n8Bf1mBZ85d/+dlizZuTHSUpXC741pe3JJ1tvPaHqbW+DX/4yXAFWUjFZCjlbuxbe8pbYKbIxc2Y4\niW39+thJJI2WpZCjcjn8Jd2qpQBhCunnP4+dQtJoWQo52rMnHKWzYEHsJNmxFKRisxRy1MpTR1XX\nXx8Wm197LXYSSaNhKeToscdavxTmzw9rC64rSMVkKeToZz+Drq7YKbL3rnfB6tWxU0gaDUshJ7t2\nwalTrXchvMFYClJxWQo5Wb06vFmWSrGTZO+GG2DjRjh+PHYSSfWyFHJSLYV2cOGFYcH50UdjJ5FU\nL0shBydOhPMTbropdpL8OIUkFZOlkIMf/xhuvBGmTo2dJD+33goPPeTnNktFYynk4HvfgzvvjJ0i\nXwsXwvLlTiFJRWMpZGzzZtixA267LXaS/N1xBzzwQOwUkuqR6rEw5XK5HDtDU3zgA7B4MXzxi7GT\n5K+vD1asgP37YfLk2GmkJiiVwkXMElUKhzc29L7uSCFDW7bAI4/Axz4WO0kcs2fDzTfDfffFTiKp\nVo4UMlIuw+/9XnhT/PjHY6eJ51e/CtNIu3fD2LGx00gNcqSg0br/fjhwAP7sz2IniWvVKpg3D77/\n/dhJJNXCkUIGtmwJ1zj6l3+Byy+PnSa+tWvhPe8Ji+7tdFiuWpAjBdVr//5w4tbXv24hVF13Hbz7\n3fD+9/tRnVLqLIUm+vd/D5fG/sxn2u+8hJF87WvwwgthxHDsWOw0koZiKTTB4cPwp38aFpa/8Q34\n8z+PnSg9F1wADz4YTmpbtgw++lF4/PGkR+JSW4pVCrcA24FdwGcjZWjYCy/A3XeHaaIZM2DnTviD\nP4idKl0XXhim1TZuDK/XH/8xLF0KX/gCHDwYO50kiFMKY4G/IxTDZcAdwMoIORqyZQt0dvawaVOY\nNvrKV2DKlNipBtfT0xM7wjkuvhi6u2H7dvjHf4STJ+Gaa+Dhh3tiR6tJaq/nUMzZPEXI2CwxSuFa\nYDewDzgN/BD4/Qg5GrJkCdx6aw8/+hEsWhQ7zfBS/YUuleBNb4K//Vvo7YUnnuiJHakmqb6eA5mz\neYqQsVlilMJ84EC/75+u3FcoF14ICxbETtE6Jk6MnUASxCkFlxYlKVExTl5bBXQT1hQA/gJ4DfhK\nv8fsBpbkG0uSCq8XWBo7RL3GEYJ3AuOB9RRwoVmS1DzvAnYQRgR/ETmLJEmSpJRMBx4FdgI/BYa6\nVNp9QB+waZTbN6rW/Qx1Ql434QirdZWvW87bsjG1nAj4jcrPNwBX17ltszSScx+wkfD6/Tq7iCNm\nXAGsBX4LfKrObZupkZz7yOe1hJFz3kn433oj8BhwRR3bNlMjOfeRzuv5+4Sc64AngRvr2DYJXwU+\nU7n9WeDLQzzuesIbxMBSqHX7RtWyn7GEqbBOoINz10nuAj6ZUbbh9lv1n4FHKrffDPyqjm1TyAmw\nl1DOWaol4yzgPwF3c+6bbWqv5VA5IZ/XEmrLeR1QPeXzFtL93RwqJ6T1evY/sPvyyuNr3faMmNc+\nug24v3L7fuD2IR63BjjewPaNqmU/I52Ql9VRXrWcCNg//+OEkc6cGreNnXN2v59nfaRcLRmfAZ6o\n/LzebVPIWZXHUYe15FwLnKjcfhy4uI5tU8hZlcrrearf7UnAs3Vse0bMUphNmBai8u/sYR6bxfbN\n3M9IJ+R9lDCsu5fmTnPVciLgUI+ZV8O2zdJITgjntvyM8Eb3oYgZs9i2Xo3uK4/XEurP+V85O1JM\n+fXsnxPSez1vB7YBq4HqBwHX9d84rqGYI3uU8FfpQF8Y8H2Zxk5qa3T7RnMOt+9vAl+q3P6fwD2E\nX6xmqPW/OfaHKTWa823AIcK0yKOEudE1TcjVX6O/f3lpdF9vBQ6T7WsJ9eV8B/BBQrZ6t21UIzkh\nvdfznypf1wPfI6wv1SXrUrh5mJ/1Ed6IjwBzgaN1Pnej2/fXaM6DQP+LXiwgtDEDHv9t4KHRx6xr\nv0M95uLKYzpq2LZZRpuzeu3UQ5V/nwH+N2E43Oz/49WSMYtt69Xovg5X/s3ytYTac14BfIswV1+d\nJk7x9RwsJ6T3elatIby/T688Lq/XsyFf5ewq+OcYfqG4k8EXmmvdvhG17Ge4E/Lm9nvcJ4AfNDFb\nLScC9l/AXcXZRbI8TyJsJOcEYHLl9kTC0R/vjJSxqptzF3BTey2rujk3Z16vJdSW8xLCXPeqUWyb\nQs7UXs8lnB1tv6ny+Fq3TcJ0wlzcwEM95wH/p9/jHiD8pfgSYV7sAyNsHyvnUCfkfZdwyNoGwrCu\n2Wsfg+33w5Wvqr+r/HwD4ZdlpMxZGG3OxYRf4vXA5oxzjpRxDuF38AThr8X9hAW9obZNLWeer2Ut\nOb8NHOPs4dq/HmHb1HKm9np+ppJjHWGk8DsjbCtJkiRJkiRJkiRJkiRJkiRJiu9k7ABSlmJeEE8q\nojyvyyPlzlKQJJ1hKUiSzrAUJElnWAqSpDMsBUnSGZaCVB+PPpIkSZIkSZIkSZIkSZIkSZIkSZIk\nSarH/wc4mhKxFdmaPgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10b1d9f90>"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "density = gaussian_kde(mi_hr80.sim)\n",
      "xs = linspace(mi_hr80.sim.min(),mi_hr80.sim.max(),200)\n",
      "plot(xs, density(xs))\n",
      "axvline(x=mi_hr80.I, color='r')\n",
      "axvline(x=mi_hr80_random.I, color='g')\n",
      "xlabel('I')\n",
      "ylabel('f(I)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 33,
       "text": [
        "<matplotlib.text.Text at 0x10d341c90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEPCAYAAACtCNj2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGRRJREFUeJzt3XuQXNVh5/FvC40APYweyHohNHohCZuX2cUiNmZMGRfe\nCgSSSrKEVOx444qpjb0xLhM/KqWRl0qZBBLsSlUeNriwy0tSKXtdwKKUcSoTy5jFgdX7zeiJnqBI\nMhIYhNT7x+mWRqN5dE/3vefc7u+nako9PX37/miP+zfnnHtvgyRJkiRJkiRJkiRJkiRJkpSs2cC/\nAhuA9cBnK/d3A68Aqypft8UIJ0nK13Tg2srt8cAWYAmwDLgvVihJ0uBGZ/jcBypfAMeBTcCsyvel\nDPcrSUpcJ7CLMGJYBuwE1gCPAhOjpZIk5W488CJwZ+X7dxNGCiXgAUIxSJISkPU0TgfwNLACeGSA\nn3cCTwFX9b1z/vz55d7e3oyjSVLL6QUWNPIEo5oUZCAlwihgI+cWwow+t+8C1vXfsLe3l3K5nPzX\nsmXLomeofm3cWIZumDKlzJYt6eYsyutpTnMOmBGiZxnqC5jf6Bt3lqXwAeB3gQ9z9vDTjwEPAmsJ\nawo3A5/LMEPb+Iu/CP/efz/86Z/GzSKpuLI8+uinDFw6KzLcZ1sql+GZZ4B74d57Yd486O2F+Q3/\nzSCp3WQ5Umh5XV1dsSMAsG4dTJgQbk+YAJ/4BPzd3539eSo5h2PO5jJn8xQhY7Oker5AuTI/phr8\n1V/B1q3wt9NLlJeV2bIFPvQh2LMHxoyJnU5qIaVSGJonqlQqQYPv644UWsCqVXDDDWe/X7QIFi+G\nJ5+Ml0lSMVkKLWD9enjve8+971Ofgm9+M04eScVlKRTcqVOweTNceeW59//Gb8BLL8HOnVFiSSoo\nS6HgenthxgwYN+7c+y++GH7nd+BRzxeXVAdLoeA2bID3vGfgn3360/D3fw/HjuWbSVJxWQoFN9T5\nCFdeCbffDl/9ar6ZJBWXpVBwO3bA3LmD//zP/gy+9z148cX8MkkqLkuh4IYrhXe/Gx5+GD75SXj7\n7fxySSomS6HghisFCAvOs2efvT6SJA3GUiiwcjkcctrZOfTjSiX4+tfD11tv5ZFMUlFZCgV24EC4\n1tH48cM/dsGCcIKbZzlLGoqlUGC7d8Pll9f++HvugR/8ILs8korPUiiwvXth1qzaH3/LLfBv/5b0\n9bwkRWYpFNi+ffWVQmcndHTAtm2ZRZJUcJZCge3dCzNn1v74UilcUnvlyuwySSo2S6HA6h0pAFxz\nTfhQHkkaiKVQYPWOFACuuipcaluSBmIpFNhIRgrvfa+lIGlwlkKB7dtX/0hh5sxwAturr2aTSVKx\nWQoFdeJEeHOfNKm+7UqlcPXUjRuzySWp2CyFgqqOEkoj+IjuefPCNZMkqT9LoaDqPXGtr7lzLQVJ\nA7MUCmok6wlVloKkwVgKBbV/f/hs5pGwFCQNxlIoqEOHYNq0kW1rKUgajKVQUIcOhU9VG4nLLguH\npP7yl83NJKn4LIWCaqQULrgApk8PU1CS1JelUFCNlAKERep9+5qXR1JrsBQKylKQlAVLoYDK5VAK\nU6eO/DksBUkDsRQK6PhxGDUKxo0b+XPMmmUpSDqfpVBAr77a2NQROFKQNDBLoYAaXU8AS0HSwCyF\nArIUJGUly1KYDfwrsAFYD3y2cv9k4FlgK/AjYGKGGVqSpSApK1mWwkngc8B7gKXAfweWAF8klMIV\nwL9UvlcdmlEKl1wCJ0+GRWtJqsqyFA4Aqyu3jwObgFnAHcDjlfsfB+7MMENLakYplEphtOBZzZL6\nymtNoRO4DngBmAYcrNx/sPK96tCMUgCnkCSdb3QO+xgPfB/4H8Dr/X5Wrnydp7u7+8ztrq4uurq6\nsklXQJaCJICenh56enqa+pwj+DDHunQATwMrgEcq920GugjTSzMIi9GL+21XLpcH7AoBV18N3/kO\nXHvtufeXlpcoL6v9dfvjP4bOzvCvpBqUSuGSAokqhc/nbeh9PcvpoxLwKLCRs4UA8CTw8crtjwM/\nzDBDS2rWSGHaNDhwoPHnkdQ6siyFDwC/C3wYWFX5ug34GnAr4ZDUWyrfq0anT8Phw3DppY0/1/Tp\nloKkc2W5pvBTBi+dj2S435Z25AhMmABjxjT+XNOnw8GDwz9OUvvwjOaCadbUEThSkHQ+S6FgmlkK\nrilI6s9SKJhmlsLUqfDaa3DqVHOeT1LxWQoF08xS6OiASZPCwrUkgaVQOM0sBXBdQdK5LIWCaXYp\nuK4gqS9LoWAcKUjKkqVQMIcOhQXiZrEUJPVlKRRMMz6fuS9PYJPUl6VQMK4pSMqSpVAgb78Nr78e\nDiNtFqePJPVlKRTIa6+FC+GNauL/apaCpL4shQJp9tQRuKYg6VyWQoFkUQqTJ8MvfhGmpiTJUiiQ\nLEph1KhwiOuhQ819XknFZCkUSBalAK4rSDrLUiiQLEvBdQVJYCkUSrPPZq7yXAVJVZZCgTT7bOYq\np48kVVkKBeKagqSsWQoF4pqCpKxZCgWSVSm4piCpylIoiBMnoFyGceOa/9xOH0mqshQKojpKKJWa\n/9yWgqQqS6Egspo6AnjXu+DkSXjjjWyeX1JxWAoFkWUplEouNksKLIWCyLIUwMVmSYGlUBBZnc1c\n5bqCJLAUCiOrs5mrLAVJYCkURtbTR64pSAJLoTBcU5CUB0uhIPIYKVgKkiyFgrAUJOXBUiiA06fD\nQnPWRx+5piDJUiiAo0dh/HgYMya7fVTXFMrl7PYhKX1Zl8JjwEFgXZ/7uoFXgFWVr9syzlB4WU8d\nAYwdCx0d8ItfZLsfSWnLuhS+zflv+mXgL4HrKl//nHGGwsujFMB1BUnZl8JK4MgA92dwrc/WlfV6\nQpXrCpJirSl8BlgDPApMjJShMPIaKXiugqQYpfA3wFzgWmA/8HCEDIXi9JGkvIyOsM9DfW5/C3hq\noAd1d3efud3V1UVXV1emoVJ26BAsWZL9fmbOhL17s9+PpObo6emhp6enqc8ZoxRmEEYIAHdx7pFJ\nZ/QthXZ36BDcfHP2+5kzB9asyX4/kpqj/x/My5cvb/g5sy6FJ4CbgUuBPcAyoIswdVQGdgB/mHGG\nwstr+qizE3buzH4/ktKVdSncPcB9j2W8z5ZjKUjKi2c0F0CeC81Hj8Kbb2a/L0lpshQSd/JkOMt4\n8uTs9zVqFFx+Oezenf2+JKXJUkjca6/BlCnhDTsPc+Y4hSS1M0shcXmdzVzluoLU3iyFxOW1nlDV\n2Qm7duW3P0lpsRQSd/BguPxEXhwpSO3NUkhcjJHC9u357U9SWiyFxB08mG8pXHEFbNuW3/4kpcVS\nSNyhQ/lOH116afj4z8OH89unpHRYConLe6RQKsHChbB1a377lJQOSyFxeY8UIEwhWQpSe7IUEpf3\nQjNYClI7sxQSVi7nP30ELjZL7ayWq6ROBG4EOgmXu94JPA8cyyyVAHj9dejogLFj892vawpS+xpq\npHAT8CTwE+C/ApcTiuFuYGXlZx/MOF9by/vEtaqFC8NI4fTp/PctKa6hRgp3AZ8HBptIuAL4NPDT\nZodSEGM9AeBd74JLLoFXXglXTZXUPoYaKdzH4IUAsLXyGGUkxpFHVUuWwObNcfYtKZ6hRgqfH+Jn\nZeAvm5xF/cRYZK5avBg2bYKPfjTO/iXFMVQpTCC8+SuSWNNHEEph48Y4+5YUz1Cl0J1XCA3s4MEw\njRPDkiXwgx/E2bekeIZaU+gGhprRngEsb2oanSP2SGHTpjj7lhTPUCOFF4F/AMYA/w/YD5SA6cD7\ngLeAh7IO2M5iHZIKMGsWnDgBR47ApElxMkjK31Ajhd8GPgysIBx2ego4Wbn928AtwDNZB2xnMUcK\npVIYLXgEktRehhopXA/MBH4L6CKMEqpcgM5BzENS4Wwp3HhjvAyS8jVUKfwt8C/APOClfj8rV+5X\nRt5+G44fh4kT42VwpCC1n6Gmj74BLAG+Dczt92UhZOzVV2HqVBgV8ZKFS5a42Cy1m1recj6deQqd\nJ+aJa1WOFKT246WzExV7PQFgwQLYvRveeituDkn5sRQSlcJIYcwYmDMHXn45bg5J+bEUErV/P8yY\nETuFU0hSu7EUEmUpSIrBUkjUvn0wc2bsFB6BJLUbSyFRjhQkxWApJCqVkUK1FPxoTqk9WAoJKpfT\nGSlMnAgTJsDevbGTSMqDpZCgo0fD4aBjx8ZOEjiFJLWPrEvhMeAgsK7PfZOBZwmf8fwjIOLVfdK0\nf38aU0dVLjZL7SPrUvg2cFu/+75IKIUrCBfc+2LGGQpn3740po6qHClI7SPrUlgJHOl33x3A45Xb\njwN3ZpyhcFIbKfgpbFL7iLGmMI0wpUTl38hX+ElPaiOFBQtg+/bYKSTlIfZCcxk/sOc8qY0UZs+G\nAwfCZzxIam1DfchOVg4SPuf5ADADODTQg7q7u8/c7urqoqurK4doadi3D5YujZ3irI6O8JnNu3bB\nwoWx00iq6unpoaenp6nPGaMUngQ+DjxY+feHAz2obym0m1TOUehr3rwwhWQpSOno/wfz8uXLG37O\nrKePngB+BiwC9gC/D3wNuJVwSOotle/VRypnM/dVLQVJrS3rkcLdg9z/kYz3W1gpnc3cl6UgtYfY\nC83q5+hRGD0axo+PneRcloLUHiyFxOzZA5dfHjvF+ebNgx07YqeQlDVLITG7d6dbCr29YXpLUuuy\nFBKTailMmgSlEhzpf366pJZiKSQm1VIolVxXkNqBpZCY3bvDGcQpshSk1mcpJCbVhWawFKR2YCkk\nJtXpI7AUpHZgKSTknXfCiWuzZsVOMjBLQWp9lkJC9u+HqVPDR3GmqHpYqqTWZSkkJOWpIwjZ9u2D\nkydjJ5GUFUshIamXwpgx4ZpMe/bETiIpK5ZCQlIvBXBdQWp1lkJCUj4ctcpSkFqbpZCQlE9cq7IU\npNZmKSRk+3aYOzd2iqFZClJrsxQSUS6HS1PPmxc7ydAsBam1WQqJOHgQxo2DCRNiJxmapSC1Nksh\nEdu3pz9KAJgyJZx57SW0pdZkKSSiKKVQvYS2n8ImtSZLIRG9vcUoBXAKSWpllkIitm+H+fNjp6iN\npSC1LkshEUWZPgJLQWpllkIiLAVJKbAUEvDmm/Af/wEzZ8ZOUhtLQWpdlkICduyAOXPgggtiJ6nN\nnDnhOk3vvBM7iaRmsxQSUKSpI4ALL4Rp0+CVV2InkdRslkICinQ4apVTSFJrshQSsHUrXHFF7BT1\nmTcPXn45dgpJzWYpJKCIpbBoUcgtqbVYCgnYtq14pbB4MWzeHDuFpGazFCJ78004cCD9T1zrz1KQ\nWpOlEFlvb/hgndGjYyepz7x54eijt96KnURSM1kKkRVx6gigowM6O11sllqNpRDZ1q2wcGHsFCPj\nFJLUeiyFyIp45FHVokWwZUvsFJKaKWYp7ATWAquAn0fMEVVRp4/AkYLUimKWQhnoAq4DboiYIyqn\njySlJPb0USny/qM6dgyOHy/O1VH7q04flcuxk0hqltgjhR8DLwKfipgjmm3bwiihVNBqnDwZLr4Y\n9u6NnURSs8Q8Ov4DwH5gKvAssBlYWf1hd3f3mQd2dXXR1dWVb7ocVEuhyK6+Gtasgcsui51Eaj89\nPT309PQ09TlT+Rt1GXAceLjyfbncBnMS3d3hMwkeeKA5z1daXqK8LN/X7QtfgEmT4MtfznW3Uhyl\nUtLzpaUw7dDQ+3qs6aOxwITK7XHAR4F1kbJEs2kTLFkSO0VjrrkG1q6NnUJSs8QqhWmEqaLVwAvA\n08CPImWJZuNGuPLK2CkaU50+ktQaYq0p7ACujbTvJLzzTrhExKJFsZM0ZvFi2LUrXNjv4otjp5HU\nqNiHpLat3t5wKOrYsbGTNGbMmHDy3fr1sZNIagZLIZJWmDqquuYap5CkVmEpRNIKi8xVloLUOiyF\nSFpppHDttbBqVewUkprBUoiklUrh+uth9Wo4eTJ2EkmNshQiOHUqXEiuVaaPLrkkfJzohg2xk0hq\nlKUQwa5dMGUKTJgw/GOL4oYb4IUXYqeQ1ChLIYI1a8LibCt5//vh5237qRhS67AUIli9OizOthJH\nClJrsBQiWLUKrrsudormuvpq2LEDXn89dhJJjbAUImjFkUJHR5gSe+ml2EkkNcJSyNnhw+ET1+bO\njZ2k+ZYuheeei51CUiMshZytXh3+oh7Vgq/8hz4EK1cO/zhJ6WrBt6a0vfRS660nVH3wg/Czn4Ur\nwEoqJkshZ88/D7/yK7FTZOPSS8NJbKtXx04iaaQshRyVy+Ev6VYtBQhTSD/5SewUkkbKUsjR9u3h\nKJ3Zs2MnyY6lIBWbpZCjVp46qrrpprDYfPp07CSSRsJSyNFzz7V+KcyaFdYWXFeQislSyNGPfwxd\nXbFTZO9jH4MVK2KnkDQSlkJOtm2DEyda70J4A7EUpOKyFHKyYkV4syyVYifJ3s03w9q1cORI7CSS\n6mUp5KRaCu3goovCgvOzz8ZOIqlelkIOjh0L5yd85COxk+THKSSpmCyFHHz/+3DLLTBxYuwk+bn9\ndnjqKT+3WSoaSyEH3/0u3HNP7BT5mjMHFi1yCkkqGkshY+vXw5YtcMcdsZPk7+674YknYqeQVA9L\nIWMPPwz33gtjxsROkr/f/E14+mk/jU0qEkshQxs2wDPPwGc/GztJHNOmwa23wmOPxU4iqVaWQkbK\nZbj/fvjSl+CSS2Kniee+++CRR+DUqdhJJNXCUsjI44/Dnj1h6qidLV0KM2fC974XO4mkWlgKGdiw\nAb7whfBGeOGFsdPE99BD8OUvw9GjsZNIGo6l0GS7d4cTtx55BK66KnaaNNx4I9x1F3ziE35Up5Q6\nS6GJ/v3fw6Wx77+//c5LGM5DD8Ebb4Qjkg4fjp1G0mAshSbYvx/+4A/gV38VvvEN+KM/ip0oPRde\nCE8+GU5qW7gQPvMZeOGFsCAvKR2xSuE2YDOwDfiTSBka9sYb8MADYZpoyhTYuhV+/ddjp0rXRReF\nabW1a8Pr9Xu/BwsWwFe+Anv3xk4nCeKUwgXAXxOK4UrgbmBJhBwN2bABOjt7WLcuTBs9+GC6h572\n9PTEjnCOyy6D7m7YvBn+6Z/g+HG4/np4+ume2NFqktrrORhzNk8RMjZLjFK4AXgZ2AmcBP4B+LUI\nORoyfz7cfnsP//iPMHdu7DRDS/UXulSC970Pvv516O2FF1/siR2pJqm+nv2Zs3mKkLFZYpTCLGBP\nn+9fqdxXKBddBLNnx07ROsaNi51AEsQpBZcWJSlRMT4ccinQTVhTAPgScBp4sM9jXgbm5xtLkgqv\nF1gQO0S9RhOCdwJjgNUUcKFZktQ8HwO2EEYEX4qcRZIkSVJKJgPPAluBHwGDfYLxY8BBYN0It29U\nrfsZ7IS8bsIRVqsqX7edt2VjajkR8BuVn68Brqtz22ZpJOdOYC3h9ft5dhGHzbgYeB74JfD5Ordt\npkZy7iSf1xKGz3kP4X/rtcBzwNV1bNtMjeTcSTqv568Rcq4CXgJuqWPbJPw5cH/l9p8AXxvkcTcR\n3iD6l0Kt2zeqlv1cQJgK6wQ6OHedZBlwX0bZhtpv1X8Bnqncfj/wf+vYNoWcADsI5ZylWjJOBf4T\n8ADnvtmm9loOlhPyeS2htpw3AtVTPm8j3d/NwXJCWq9n3wO7r6o8vtZtz4h57aM7gMcrtx8H7hzk\ncSuBIw1s36ha9jPcCXlZHeVVy4mAffO/QBjpTK9x29g5p/X5edZHytWS8VXgxcrP6902hZxVeRx1\nWEvO54FjldsvAJfVsW0KOatSeT1P9Lk9Hnitjm3PiFkK0wjTQlT+nTbEY7PYvpn7Ge6EvM8QhnWP\n0txprlpOBBzsMTNr2LZZGskJ4dyWHxPe6D4VMWMW29ar0X3l8VpC/Tn/G2dHiim/nn1zQnqv553A\nJmAFUP0g4Lr+G0c3FHN4zxL+Ku3vK/2+L9PYSW2Nbt9ozqH2/TfAVyu3/yfwMOEXqxlq/W+OcT5K\nX43m/CCwjzAt8ixhbnRlE3L11ejvX14a3dcHgP1k+1pCfTk/DHySkK3ebRvVSE5I7/X8YeXrJuC7\nhPWlumRdCrcO8bODhDfiA8AM4FCdz93o9n01mnMv0PeiF7MJbUy/x38LeGrkMeva72CPuazymI4a\ntm2WkeasXjt1X+XfV4H/TRgON/v/eLVkzGLbejW6r/2Vf7N8LaH2nFcD3yTM1VeniVN8PQfKCem9\nnlUrCe/vkyuPy+v1bMifc3YV/IsMvVDcycALzbVu34ha9jPUCXkz+jzuc8D/amK2Wk4E7LuAu5Sz\ni2R5nkTYSM6xwITK7XGEoz8+GiljVTfnLuCm9lpWdXNuzrxeS6gt5+WEue6lI9g2hZypvZ7zOTva\nfl/l8bVum4TJhLm4/od6zgT+T5/HPUH4S/EtwrzY7w+zfaycg52Q9x3CIWtrCMO6Zq99DLTfP6x8\nVf115edrCL8sw2XOwkhzziP8Eq8G1mecc7iM0wm/g8cIfy3uJizoDbZtajnzfC1ryfkt4DBnD9f+\n+TDbppYztdfz/kqOVYSRwn8eZltJkiRJkiRJkiRJkiRJkiRJUnzHYweQshTzgnhSEeV5XR4pd5aC\nJOkMS0GSdIalIEk6w1KQJJ1hKUiSzrAUpPp49JEkSZIkSZIkSZIkSZIkSZIkSZIkSVI9/j/OnO53\ngRNC1AAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10d035f90>"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}