08_global_south.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import HTML\n",
      "import folium\n",
      "import random\n",
      "import pandas as pd\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 50
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Base layer for south"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import folium\n",
      "map_osm = folium.Map(location=[33.7550, -84.3900], zoom_start=7)\n",
      "map_osm.geo_json(geo_path = 'data/south.json')\n",
      "map_osm.create_map(path='south_base.html')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "open up a browser and load south_base.html"
     ]
    },
    {
     "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": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q.keys()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 53,
       "text": [
        "[u'type', u'features']"
       ]
      }
     ],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "features = q['features']\n",
      "len(features)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 54,
       "text": [
        "1412"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "feature_0 = features[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "feature_0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 56,
       "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": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 56
    },
    {
     "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": 57
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df = pd.DataFrame({'HR80': values,\n",
      "                   'FIPS': indices} )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "map_osm = folium.Map(location=[33.7550, -84.3900], zoom_start=7,)\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",
      "map_osm.create_map('homocide_1.html')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = np.array(df.HR80.tolist())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pysal as ps"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 61
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "max_breaks = ps.Maximum_Breaks(y, 5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 62
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins = max_breaks.bins.tolist()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "map_osm = folium.Map(location=[33.7550, -84.3900], zoom_start=7,)\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 (Maximum Breaks)')\n",
      "map_osm.create_map('homocide_MB.html')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 64
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q5 = ps.Quantiles(y, 5).bins.tolist()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 65
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "map_osm = folium.Map(location=[33.7550, -84.3900], zoom_start=7,)\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=q5[:-1],\n",
      "                 legend_name='Homicide Rate HR80 (Quintiles)')\n",
      "map_osm.create_map('homocide_q5.html')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fj= ps.Fisher_Jenks(y, 5).bins"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 67
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "map_osm = folium.Map(location=[33.7550, -84.3900], zoom_start=7,)\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=fj[:-1],\n",
      "                 legend_name='Homicide Rate HR80 (Fisher Jenks)')\n",
      "map_osm.create_map('homocide_fj.html')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Now in the browser"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Prior to starting this notebook, issue the following at a shell:\n",
      "\n",
      "```\n",
      "python -m SimpleHTTPServer 8015 &ipython notebook --pylab=inline\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Function that shows the map - call this functionto redraw the map.\n",
      "def showmap(name):\n",
      "    sint = str(random.randint(0,10000))\n",
      "    out_file = name+\"_\"+sint+\".html\"\n",
      "    map_osm.create_map(path=out_file)\n",
      "    iframe = \"<iframe src=http://localhost:8015/%s width=1000 height=500></iframe>\"%out_file\n",
      "    print iframe\n",
      "    return HTML(iframe)\n",
      "\n",
      "def resetmap():\n",
      "    return folium.Map(location=[40, -99], zoom_start=4)\n",
      "\n",
      "map_osm = resetmap()\n",
      "showmap('map_0')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<iframe src=http://localhost:8015/map_0_9024.html width=1000 height=500></iframe>\n"
       ]
      },
      {
       "html": [
        "<iframe src=http://localhost:8015/map_0_9024.html width=1000 height=500></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 69,
       "text": [
        "<IPython.core.display.HTML at 0x10c9baa10>"
       ]
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins= ps.Quantiles(y, 5).bins.tolist()\n",
      "map_osm = folium.Map(location=[33.7550, -84.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",
      "showmap('test')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<iframe src=http://localhost:8015/test_8040.html width=1000 height=500></iframe>\n"
       ]
      },
      {
       "html": [
        "<iframe src=http://localhost:8015/test_8040.html width=1000 height=500></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 71,
       "text": [
        "<IPython.core.display.HTML at 0x10c9709d0>"
       ]
      }
     ],
     "prompt_number": 71
    },
    {
     "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": 72
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w.n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 73,
       "text": [
        "1412"
       ]
      }
     ],
     "prompt_number": 73
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w.transform = 'r'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 74
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.random.seed(12345)\n",
      "mi_hr80 = ps.Moran(y, w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 75
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80.p_sim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 76,
       "text": [
        "0.001"
       ]
      }
     ],
     "prompt_number": 76
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80.I"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 77,
       "text": [
        "0.2521152652175852"
       ]
      }
     ],
     "prompt_number": 77
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80.EI"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 78,
       "text": [
        "-0.0007087172218284905"
       ]
      }
     ],
     "prompt_number": 78
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80.p_norm"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 79,
       "text": [
        "0.0"
       ]
      }
     ],
     "prompt_number": 79
    },
    {
     "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": 80
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df['yrandom'] = yr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 81
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "type(y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 82,
       "text": [
        "numpy.ndarray"
       ]
      }
     ],
     "prompt_number": 82
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins= ps.Quantiles(y, 5).bins.tolist()\n",
      "map_osm = folium.Map(location=[33.7550, -84.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",
      "showmap('test')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<iframe src=http://localhost:8015/test_4745.html width=1000 height=500></iframe>\n"
       ]
      },
      {
       "html": [
        "<iframe src=http://localhost:8015/test_4745.html width=1000 height=500></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 83,
       "text": [
        "<IPython.core.display.HTML at 0x10cb72810>"
       ]
      }
     ],
     "prompt_number": 83
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins= ps.Quantiles(y, 5).bins.tolist()\n",
      "map_osm = folium.Map(location=[33.7550, -84.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",
      "showmap('test')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<iframe src=http://localhost:8015/test_3989.html width=1000 height=500></iframe>\n"
       ]
      },
      {
       "html": [
        "<iframe src=http://localhost:8015/test_3989.html width=1000 height=500></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 84,
       "text": [
        "<IPython.core.display.HTML at 0x10c9bcd50>"
       ]
      }
     ],
     "prompt_number": 84
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80_random = ps.Moran(df['yrandom'],w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 85
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80_random.I"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 86,
       "text": [
        "-0.0075100395407744691"
       ]
      }
     ],
     "prompt_number": 86
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi_hr80_random.p_sim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 87,
       "text": [
        "0.34799999999999998"
       ]
      }
     ],
     "prompt_number": 87
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import gaussian_kde"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 88
    },
    {
     "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": 89,
       "text": [
        "<matplotlib.text.Text at 0x10cbf0990>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEPCAYAAACtCNj2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHi5JREFUeJzt3X9wVNWhB/DvBQIICRB+uEnYyJKQkADJJkIJ0MTGYmJh\nxOKzrwOvw0slOiNTn+NAFec5luVNi2GU+oR2is9Ry3t1GOhUkL5qWi2uIsZiMZHKr0AIJCS7AQIR\nAj4SyHl/XHdNIAm72XvvOXf3+5nJCMnevV+3dr855+y5VxNCCBAREQEYJDsAERGpg6VARERBLAUi\nIgpiKRARURBLgYiIglgKREQUZFopNDY24u6778b06dMxY8YMbNy4EQDg8XjgdDqRn5+P/Px8VFZW\nmhWBiIjCpJm1T8Hv98Pv9yMvLw/t7e2YOXMmdu7cie3btyMhIQErV64047RERBSBIWY9cVJSEpKS\nkgAA8fHxyM7ORlNTEwCA++WIiNRkyZrCyZMnUV1djTlz5gAANm3aBLfbjfLycrS1tVkRgYiIQmB6\nKbS3t+MHP/gBXnrpJcTHx2PFihWor69HTU0NkpOTsWrVKrMjEBFRqISJOjo6RGlpqXjxxRd7/Xl9\nfb2YMWPGTd9PT08XAPjFL37xi19hfKWnp0f8vm3aSEEIgfLyckybNg1PPPFE8Ps+ny/45x07diAn\nJ+emY+vq6iCEUP5rzZo10jMEvg4dEhg9WmDcOIGjR9XNaZfXkzmZs9eMX7+3qfpVV1cX8Xu3aQvN\ne/fuxe9+9zvk5uYiPz8fALBu3Tps3boVNTU10DQNkydPxssvv2xWhJjy/PPA6tXA4MHAs88C27bJ\nTkREdmRaKRQWFqKrq+um7y9YsMCsU8YsIYC33waqqoDx44G0NKCuDkhPl52MiOyGO5ojUFxcLDsC\nAOAf/wASEoDJk/V//vjHQPcBmCo5b4U5jcWcxrFDRqOYtnktEpqmQcFYynrxRaC2FvjNb/S/Hz0K\n3HUX0NgIDB0qNxtRVNE0fWiuKCPeOzlSiALV1cDs2d/8fepUICsL2LVLXiYisieWQhT44gtgxoye\n33vkEeCVV+TkISL74vSRzV2/rq8jnD0LjBz5zfe/+gpITQX+/nfA5ZIWjyi6cPqIVFdXByQn9ywE\nALjtNuBf/gV49VU5uYjInlgKNnfwIDB9eu8/e/RR4L/+C/jyS2szEZF9sRRsrr/9CNOmAYsWAf/x\nH9ZmIiL7YinYXH29vj+hL+vWAW+8oa8tEBHdCkvB5m5VCrffDmzYACxfDnR0WJeLiOyJpWBztyoF\nQF9wTk3Vr49ERNQffiTVxoQARozQP44aH9//Y48fB+bN03c5DxtmTT6iqMOPpJLK/H59j8KtCgEA\npkzRN7hxlzMR9YelYGMNDcAdd4T++B/9CHjzTfPyEJH9sRRsrKkJmDgx9Md/97vABx8oPfolIslY\nCjbW3BxeKbhcQFwccOyYaZGIyOZYCjbW1ASkpIT+eE3TL6m9Z495mYjI3lgKNhbuSAEA3G79pjxE\nRL1hKdhYuCMFAMjJ0S+1TUTUG5aCjQ1kpDBjBkuBiPrGUrCx5ubwRwopKcDVq/qGNyKiG7EUbOry\nZf3NPTExvOM0Tb966qFD5uQiIntjKdhUYJSgaeEfm5amXzOJiOhGLAWbCnfjWneTJ7MUiKh3LAWb\nGsh6QgBLgYj6wlKwKZ9PvzfzQLAUiKgvLAWbOnMGcDgGdixLgYj6wlKwqTNn9LuqDYTTqX8k9f/+\nz9hMRGR/LAWbiqQUBg8GkpL0KSgiou5YCjYVSSkA+iJ1c7NxeYgoOrAUbIqlQERmYCnYkBB6KUyY\nMPDnYCkQUW9YCjbU3g4MGgSMHDnw55g4kaVARDdjKdjQ2bORTR0BHCkQUe9YCjYU6XoCwFIgot6x\nFGyIpUBEZjGtFBobG3H33Xdj+vTpmDFjBjZu3AgAOH/+PEpKSpCZmYnS0lK0tbWZFSFqsRSIyCym\nlUJcXBxefPFFHDx4EJ988gl+/etf4/Dhw6ioqEBJSQlqa2sxf/58VFRUmBUhahlRCqNHA52d+qI1\nEVGAaaWQlJSEvLw8AEB8fDyys7PR1NSEXbt2oaysDABQVlaGnTt3mhUhahlRCpqmjxa4q5mIurNk\nTeHkyZOorq5GQUEBWlpa4Pj6Sm4OhwMtLS1WRIgqRpQCwCkkIrrZELNP0N7ejgcffBAvvfQSEhIS\nevxM0zRofdw6zOPxBP9cXFyM4uJiE1PaC0uBiADA6/XC6/Ua+pymlkJnZycefPBBLFu2DIsXLwag\njw78fj+SkpLg8/lwex/vbt1LgXqKdDdzQFISwIEakX3d+Avz2rVrI35O06aPhBAoLy/HtGnT8MQT\nTwS/f//992PLli0AgC1btgTLgkJn1EjB4QD8/sifh4iihyaEEGY88UcffYS77roLubm5wSmi5557\nDrNnz8YPf/hDNDQ0wOVyYfv27RgzZkzPUJoGk2LZXlcXMGwYcPkyMHRoZM/1+uvABx8Av/2tIdGI\nop+m6RcfU5QR752mTR8VFhaiq6ur15+99957Zp026l24ACQkRF4IAKePiOhm3NFsM0ZNHQF6KXD6\niIi6YynYjJGlwDUFIroRS8FmjCyFCROAc+eA69eNeT4isj+Wgs0YWQpxcUBiItDaaszzEZH9sRRs\nxshSALiuQEQ9sRRsxuhS4LoCEXXHUrAZjhSIyEwsBZsx6hIXASwFIuqOpWAzRtyfuTtuYCOi7lgK\nNsM1BSIyE0vBRjo6gEuX9I+RGoXTR0TUHUvBRs6dA8aPBwYZ+L8aS4GIumMp2IjRU0cA1xSIqCeW\ngo2YUQpjxwIXL+pTU0RELAUbMaMUBg3SP+J65oyxz0tE9sRSsBEzSgHgugIRfYOlYCNmlgLXFYgI\nYCnYitG7mQO4V4GIAlgKNmL0buYATh8RUQBLwUa4pkBEZmMp2AjXFIjIbCwFGzGrFLimQEQBLAWb\nuHwZEAIYOdL45+b0EREFsBRsIjBK0DTjn5ulQEQBLAWbMGvqCABGjQI6O4ErV8x5fiKyD5aCTZhZ\nCprGxWYi0rEUbMLMUgC42ExEOpaCTZi1mzmA6wpEBLAUbMOs3cwBLAUiAlgKtmH29BHXFIgIYCnY\nBtcUiMgKLAWbsGKkwFIgIpaCTbAUiMgKLAUb6OrSF5rN/vQR1xSIiKVgA21tQHw8MHSoeecIrCkI\nYd45iEh9ppbC8uXL4XA4kJOTE/yex+OB0+lEfn4+8vPzUVlZaWaEqGD21BEAjBgBxMUBFy+aex4i\nUpuppfDQQw/d9KavaRpWrlyJ6upqVFdX43vf+56ZEaKCFaUAcF2BiEwuhaKiIiQmJt70fcE5irCY\nvZ4QwHUFIpKyprBp0ya43W6Ul5ejra1NRgRbsWqkwL0KRGR5KaxYsQL19fWoqalBcnIyVq1aZXUE\n2+H0ERFZZYjVJ7y927vbww8/jEWLFvX6OI/HE/xzcXExiouLTU6mrjNngOxs88+TkgI0NZl/HiIy\nhtfrhdfrNfQ5LS8Fn8+H5ORkAMCOHTt6fDKpu+6lEOvOnAG+8x3zzzNpEvD55+afh4iMceMvzGvX\nro34OU0thaVLl+KDDz7AuXPnkJqairVr18Lr9aKmpgaapmHy5Ml4+eWXzYwQFayaPnK5gJMnzT8P\nEalLEwp+FEjTNH5CqZvsbOAPfwCmTTP3PM3NwJ13cl2BqE+apvQOTyPeO7mj2QasXGhuawO++sr8\ncxGRmlgKiuvs1HcZjx1r/rkGDQLuuANoaDD/XESkJpaC4s6dA8aN09+wrTBpEtcViGIZS0FxVu1m\nDuBiM1FsYykozqr1hACXCzh1yrrzEZFaWAqKa2nRLz9hFY4UiGIbS0FxMkYKJ05Ydz4iUgtLQXEt\nLdaWQmYmcOyYdecjIrWwFBR35oy100fjx+u3/2xtte6cRKQOloLirB4paBqQkQHU1lp3TiJSB0tB\ncVaPFAB9ComlQBSbWAqKs3qhGWApEMUyloLChLB++gjgYjNRLLvlpbPb2tpQVVWFkydPQtM0uFwu\nzJ07F6NHj7YiX0y7dAmIiwNGjLD2vFxTIIpdfV46e8+ePXj++edx8uRJ5OfnIyUlBUII+Hw+VFdX\nw+Vy4amnnkJhYaHxoXjpbAD6b+sLFgDHj1t73osXgeRkvZSsuuYSkS3EwKWz+xwp7NixAxs2bEBG\nRkavP6+trcXmzZtNKQXSyVhPAIBRo4DRo4HTp/WrphJR7OBNdhS2Ywfw3/+t/9Nq8+cDq1cDpaXW\nn5tIWbE8UtiwYUO/J165cmVEJ6Zbk7HIHJCVBRw+zFIgijV9lsKlS5egaZqVWegGsqaPAL0UDh2S\nc24ikqfPUvB4PBbGoN60tOj3Z5YhOxt480055yYiefr8bInH40FLS0ufB/p8PqxZs8aUUKSTPVI4\nfFjOuYlInj5HCrNmzcKSJUvQ0dGBO++8E8nJyRBCwO/347PPPsOwYcPw05/+1MqsMcfqeyl0N3Ei\ncPkycOECkJgoJwMRWa/PkcK2bdvw/vvvY8GCBSgsLMTgwYMRFxeHwsJCbNu2Dbt378bChQutzBpz\nZI4UNE0fLRw5Iuf8RCRHnyOF/fv3o7m5Gdu3b4fX6+3xMScuQFtDxsXwuguUwty58jIQkbX6LIVH\nH30U8+fPx4kTJzBz5sweP9M0DSd4ey5TdXQA7e3AmDHyMnCkQBR7brl57dFHH8XmzZutygOAm9cA\noKkJmD1b/6csb74J/Pa3wK5d8jIQKSUGNq/d8so2VhcC6WRuXAvgSIEo9vByZ4qSvZ4AAFOmAA0N\nwNWrcnMQkXVYCopSYaQwdCgwaZL1V2klInlYCory+fTLV8vGKSSi2MJSUBRLgYhkYCkoqrkZSEmR\nnUK/BhIvd0EUO1gKiuJIgYhkYCkoSpWRQqAUurpkJyEiK7AUFCSEOiOFMWOAhAS5m+iIyDosBQW1\ntekfBx0xQnYSHaeQiGKHqaWwfPlyOBwO5OTkBL93/vx5lJSUIDMzE6WlpWhrazMzgi35fGpMHQVw\nsZkodphaCg899BAqKyt7fK+iogIlJSWora3F/PnzUVFRYWYEW2puVmPqKIAjBaLYYWopFBUVIfGG\nO7Ts2rULZWVlAICysjLs3LnTzAi2pNpIgXdhI4odlq8ptLS0wPH1RX0cDke/t/yMVaqNFKZMAXil\ndKLYIHWhWdM03rCnF6qNFFJTAb9fv8cDEUW3Pm+yYxaHwwG/34+kpCT4fD7c3sdV3zweT/DPxcXF\nKC4utiagApqbgTlzZKf4Rlycfs/mU6eAjAzZaYgowOv1wuv1GvqclpfC/fffjy1btmD16tXYsmUL\nFi9e3OvjupdCrFFlj0J3aWn6FBJLgUgdN/7CvHbt2oif09Tpo6VLl2LevHk4evQoUlNT8frrr+Pp\np5/Gu+++i8zMTOzevRtPP/20mRFsSZXdzN0FSoGIotstb8cpQyzfjlMIYORI/SY78fGy03yjogJo\nbQWef152EiKJeDtOslpbGzBkiFqFAHCkQBQrWAqKaWwE7rhDdoqbpaUB9fWyUxCR2VgKimloULcU\n6uqUHjkTkQFYCopRtRQSE/Xp1AsXZCchIjOxFBSjailoGtcViGIBS0ExDQ36DmIVsRSIoh9LQTGq\nLjQDLAWiWMBSUIyq00cAS4EoFrAUFHLtmn6Ji4kTZSfpHUuBKPqxFBTi8wETJui34lRR4GOpRBS9\nWAoKUXnqCNCzNTcDnZ2ykxCRWVgKClG9FIYO1a/e2tgoOwkRmYWloBDVSwHgugJRtGMpKETlj6MG\nsBSIohtLQSEqb1wLYCkQRTeWgkJOnAAmT5adon8sBaLoxlJQhBD6panT0mQn6R9LgSi6sRQU0dKi\n33EtIUF2kv6xFIiiG0tBESdOqD9KAIBx4/Sd17yENlF0Yikowi6lELiENu/CRhSdWAqKqKuzRykA\nnEIiimYsBUWcOAGkp8tOERqWAlH0Yikowi7TRwBLgSiasRQUwVIgIhWwFBTw1VfA+fNASorsJKFh\nKRBFL5aCAurrgUmTgMGDZScJzaRJ+nWarl2TnYSIjMZSUICdpo4AYNgwwOEATp+WnYSIjMZSUICd\nPo4awCkkoujEUlBAbS2QmSk7RXjS0oDjx2WnICKjsRQUYMdSmDpVz01E0YWloIBjx+xXCllZwJEj\nslMQkdFYCpJ99RXg96t/x7UbsRSIohNLQbK6Ov3GOkOGyE4SnrQ0/dNHV6/KTkJERmIpSGbHqSMA\niIsDXC4uNhNFG5aCZLW1QEaG7BQDwykkoujDUpDMjp88Cpg6FTh6VHYKIjKStJlsl8uFUaNGYfDg\nwYiLi8O+fftkRZHq2DFg2TLZKQYmKwt4/33ZKYjISNJKQdM0eL1ejB07VlYEJdh9+ug3v5GdgoiM\nJHX6SAgh8/TSffkl0N5un6uj3igwfRTj/zMSRRVppaBpGu655x7MmjULr7zyiqwYUh07po8SNE12\nkoEZOxa47TagqUl2EiIyirTpo7179yI5ORlnz55FSUkJsrKyUFRUFPy5x+MJ/rm4uBjFxcXWhzRZ\noBTsLDcX+PxzwOmUnYQo9ni9Xni9XkOfUxMKzOGsXbsW8fHxWLVqFQB9FKFALNN5PPo9CX7+c9lJ\nBu7JJ4HERODf/112EiILaJrS86VGvHdKmT66cuUKLl26BAC4fPky/vKXvyAnJ0dGFKkOHways2Wn\niIzbDRw4IDsFERlFSim0tLSgqKgIeXl5KCgowH333YfS0lIZUaQ6dAiYNk12isgEpo+IKDooMX10\no1iYPrp2DUhIAFpbgREjZKcZuI4OYMwY/d/jtttkpyEyGaePyCx1dfpHUe1cCAAwdKi+I/uLL2Qn\nISIjsBQkiYapowC3m1NIRNGCpSBJNCwyB7AUiKIHS0GSaBop5OUB1dWyUxCREVgKkkRTKcycCdTU\nAJ2dspMQUaRYChJcv67fhyBapo9Gj9ZvJ3rwoOwkRBQploIEp04B48bpH0mNFrNnA3/7m+wURBQp\nloIEn3+uL85Gk4ICIEZviUEUVVgKEtTU6Iuz0YQjBaLowFKQoLoayM+XncJYublAfT3w9SWtiMim\nWAoSRONIIS5OnxLbv192EiKKBEvBYq2t+h3XJk+WncR4c+YAe/fKTkFEkWApWKymRv+NelAUvvJ3\n3QXs2SM7BRFFIgrfmtS2f3/0rScEFBYCH3+sXwGWiOyJpWCxqipg3jzZKcwxfry+ia2mRnYSIhoo\nloKFhNB/k47WUgD0KaQPP5SdgogGiqVgoRMn9E/ppKbKTmIelgKRvbEULBTNU0cBRUX6YnNXl+wk\nRDQQLAUL7d0b/aUwcaK+tsB1BSJ7YilY6L33gOJi2SnMt2AB8M47slMQ0UCwFCxy7Bhw+XL0XQiv\nNywFIvtiKVjknXf0N0tNk53EfN/5DnDgAHDhguwkRBQuloJFAqUQC4YP1xec331XdhIiChdLwQJf\nfqnvT7jnHtlJrMMpJCJ7YilY4A9/AL77XWDMGNlJrLNoEfDHP/K+zUR2w1KwwP/8D/CjH8lOYa1J\nk4CpUzmFRGQ3LAWTffEFcPQocP/9spNYb+lSYOtW2SmIKByaEELIDnEjTdOgYKwBeeghIC0NePZZ\n2Ums19ICZGUBDQ1AQoLsNEQG0DT9ImaKMuK9kyMFEx08CLz9NvD447KTyOFwACUlwGuvyU5CRKHi\nSMEkQgD33ae/KT7xhOw08nzyiT6NdPw4MHiw7DREEeJIgQZqyxagsRFYsUJ2ErnmzAFSUoA33pCd\nhIhCwZGCCQ4e1K9xtHs3kJMjO418VVXAP/+zvugeSx/LpSjEkQKFq6FB37j1n//JQgiYOxd44AHg\nxz/mrTqJVMdSMNCnn+qXxn7qqdjbl3ArL7wAXLmijxhaW2WnIaK+sBQM4PMBDz+sLyxv3Ag89pjs\nROoZNgzYtUvf1JaRAfzbvwF/+5vSI3GimCSlFCorK5GVlYWMjAysX79eRgRDXLkC/Pzn+jTRuHFA\nbS3wT/8kO5W6hg/Xp9UOHNBfr3/9V2DKFOCZZ4CmJtnpiAiQUArXr1/HY489hsrKShw6dAhbt27F\n4cOHrY4RsYMHAZfLi3/8Q582Wr8eGD1adqreeb1e2RF6cDoBjwc4cgT4/e+B9nZg5kzgf//XKzta\nSFR7PfvCnMaxQ0ajWF4K+/btw5QpU+ByuRAXF4clS5bgrbfesjpGxNLTgUWLvNi2DZg8WXaa/qn6\nH7SmAXfeCbz0ElBXB/z9717ZkUKi6ut5I+Y0jh0yGsXyUmhqakJqamrw706nE002nDsYPhzo9q9B\nERo5UnYCIgIklIIWC7ceIyKyK2Gxqqoqce+99wb/vm7dOlFRUdHjMenp6QIAv/jFL37xK4yv9PT0\niN+jLd/RfO3aNUydOhV//etfkZKSgtmzZ2Pr1q3Izs62MgYREfViiOUnHDIEv/rVr3Dvvffi+vXr\nKC8vZyEQESlCyWsfERGRHNJ2NJ8/fx4lJSXIzMxEaWkp2traen3c8uXL4XA4kHPDhYRCPd6qnH1t\nyPN4PHA6ncjPz0d+fj4qKysNzRfKRsDHH38cGRkZcLvdqK6uDutYFXK6XC7k5uYiPz8fs2fPlpbx\nyJEjmDt3LoYPH44NGzaEdawqOa16LUPJ+cYbb8DtdiM3Nxff/va3ceDAgZCPVSWnSq/nW2+9Bbfb\njfz8fMycORO7d+8O+dgeIl6VGKAnn3xSrF+/XgghREVFhVi9enWvj/vwww/FZ599JmbMmDGg463I\nee3aNZGeni7q6+tFR0eHcLvd4tChQ0IIITwej9iwYYMp2fo7b8Cf/vQnsWDBAiGEEJ988okoKCgI\n+VgVcgohhMvlEq2traZkCyfjmTNnxKeffiqeeeYZ8cILL4R1rAo5hbDmtQw158cffyza2tqEEEK8\n8847yv632VdOIdR6Pdvb24N/PnDgQHDROdzXU9pIYdeuXSgrKwMAlJWVYefOnb0+rqioCImJiQM+\n3oqct9qQJ0yaoQtlI2D3/AUFBWhra4Pf77d0E+FAc7a0tAR/btZrGE7GCRMmYNasWYiLiwv7WBVy\nBpj9Woaac+7cuRj99WUACgoKcPr06ZCPVSFngCqv58hum33a29sxfvz4kI/tTloptLS0wOFwAAAc\nDkePNwArjjfyPLfakLdp0ya43W6Ul5cbOs0VykbAvh7T3Nxs2SbCSHIC+t6We+65B7NmzcIrr7wi\nLaMZx4Yr0nNZ8VoC4ed89dVXsXDhwgEdKysnoN7ruXPnTmRnZ2PBggXYuHFjWMcGmPrpo5KSEvj9\n/pu+/4tf/KLH3zVNi2hTW6THR5qzv3OvWLECP/vZzwAAzz77LFatWoVXX311wFlDPW93Vvwm059I\nc3700UdISUnB2bNnUVJSgqysLBQVFRkZMeL//qwS6bn27t2L5ORkU19LILyc77//Pl577TXs3bs3\n7GMjFUlOQL3Xc/HixVi8eDH27NmDZcuW4ciRI2Gfy9RSePfdd/v8mcPhgN/vR1JSEnw+H26//faw\nnjvS443MOXHiRDQ2Ngb/3tjYCKfTCQA9Hv/www9j0aJFA84Zznn7eszp06fhdDrR2dl5y2Nl55w4\ncSIAICUlBYA+LfLAAw9g3759hv8fL5SMZhwbrkjPlZycDMDc1zKcnAcOHMAjjzyCysrK4DSxiq9n\nbzkB9V7PgKKiIly7dg3nz5+H0+kM7/U0ekEkVE8++WRwJ/Nzzz3X70JxfX19rwvNoR5vds7Ozk6R\nlpYm6uvrxdWrV3ss5DQ3Nwcf98tf/lIsXbrUsGz9nTeg+wJuVVVVcJEslGNVyHn58mVx8eJFIYS+\nkDZv3jzx5z//WUrGgDVr1vRYwFXttewrp1WvZag5T506JdLT00VVVVXYx6qQU7XX8/jx46Krq0sI\nIcT+/ftFWlpayMd2J60UWltbxfz580VGRoYoKSkRFy5cEEII0dTUJBYuXBh83JIlS0RycrIYOnSo\ncDqd4rXXXuv3eFk53377bZGZmSnS09PFunXrgt9ftmyZyMnJEbm5ueL73/++8Pv9hubr7bybN28W\nmzdvDj7mJz/5iUhPTxe5ubli//79t8xshoHmrKurE263W7jdbjF9+nRTc94qo8/nE06nU4waNUqM\nGTNGpKamikuXLvV5rGo5rXwtQ8lZXl4uxo4dK/Ly8kReXp741re+1e+xquVU7fVcv369mD59usjL\nyxOFhYVi3759/R7bF25eIyKiIN6Ok4iIglgKREQUxFIgIqIglgIREQWxFIiIKIilQEREQSwFojDE\nx8fLjkBkKpYCURisvC4PkQwsBSIiCmIpEBFREEuBiIiCWApERBTEUiAioiCWAlEY+Okjina8dDYR\nEQVxpEBEREEsBSIiCmIpEBFREEuBiIiCWApERBTEUiAioiCWAhERBbEUiIgo6P8B3r/EU9k9jioA\nAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c9321d0>"
       ]
      }
     ],
     "prompt_number": 89
    },
    {
     "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": 90,
       "text": [
        "<matplotlib.text.Text at 0x10cba3750>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEPCAYAAACtCNj2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHdhJREFUeJzt3X9wVNWhB/DvBRJQEiD8cJOwkZCQkADJZoUSsImNxcTC\nCMVnOwPt8FKIzsjU5zhQlNaxLG9aGuaJFmjn4XO0wzydFN+8grxW02JhFSMWi4lUfgVDgJBkw88o\nAUsCOe+P664JZMNu9t57zt39fmZ2Jmz27v161f1yztl7ryaEECAiIgIwSHYAIiJSB0uBiIgCWApE\nRBTAUiAiogCWAhERBbAUiIgowLRSaGpqwv3334+pU6di2rRp2LRpEwDA4/HA6XTC7XbD7Xajurra\nrAhERBQmzazzFHw+H3w+HwoKCtDR0YHp06djx44deOONN5CYmIgVK1aYsVsiIorAELPeODk5GcnJ\nyQCAhIQE5Obmorm5GQDA8+WIiNRkyZrCyZMnUVtbi1mzZgEANm/eDJfLhYqKCrS3t1sRgYiIQmB6\nKXR0dOB73/seNm7ciISEBCxfvhyNjY2oq6tDSkoKVq5caXYEIiIKlTBRZ2enKCsrEy+++GKfv29s\nbBTTpk275fnMzEwBgA8++OCDjzAemZmZEX9umzZSEEKgoqICU6ZMwVNPPRV4vrW1NfDz9u3bkZeX\nd8u2DQ0NEEIo/1izZo30DP7H4cMC8ABjxggcO6ZuTrscT+Zkzj4zfvXZpuqjoaEh4s9u0xaaa2pq\n8NprryE/Px9utxsAsG7dOlRVVaGurg6apmHixIl46aWXzIoQU/7jPwBMAJ5+GnjuOWDbNtmJiMiO\nTCuFoqIidHd33/L83LlzzdplzBICeOstAMuB5cuBjAygoQHIzJSdjIjshmc0R6CkpER2BADAP/4B\nJCbqPycmAj/6EdBzAKZKztthTmMxp3HskNEopp28FglN06BgLGW9+CJQXw9sSdYg1ggcOwbcdx/Q\n1ATEx8tORxRFNE0fmivKiM9OjhSiQG0tMHPm13+ePBnIyQF27pSXiYjsiaUQBT79FJg2rfdzjz0G\nvPyynDxEZF8sBZu7cQM4ehSYMqX38488Ahw4AJw8KSUWEdkUS8HmGhqAlBRg+PDez99xB/CDHwCv\nvCInFxHZE0vB5g4dAqZO7ft3jz8O/Nd/AZ9/bm0mIrIvloLN9Xc+wpQpwPz5wL//u7WZiMi+WAo2\n19gITJwY/Pfr1gGvvw78/e/WZSIi+2Ip2NztSuGuu4ANG4Bly4DOTutyEZE9sRRs7nalAOgLzmlp\nX10fiYioHywFGxNC/8ppenr/r9M0YONG/XHtmhXJiMiuWAo25vPp1zpKSLj9aydN0k9w41nORNQf\nloKNnT4N3H136K//4Q+BP/zBvDxEZH8sBRtrbgbGjw/99d/+NvDuu0pfz4uIJGMp2FhLS3ilkJ4O\nxMUBx4+bFomIbI6lYGPNzUBqauiv1zT9ktp795qXiYjsjaVgY+GOFADA5dJvykNE1BeWgo2FO1IA\ngLw8/VLbRER9YSnY2EBGCtOmsRSIKDiWgo21tIQ/UkhN1U9gO3fOnExEZG8sBZu6ckX/cE9KCm87\nTdOvnnr4sDm5iMjeWAo25R8laFr422Zk6NdMIiK6GUvBpsI9ca2niRNZCkTUN5aCTQ1kPcGPpUBE\nwbAUbKq1Vb8380CwFIgoGJaCTZ09CzgcA9uWpUBEwbAUbOrsWf2uagPhdOpfSf3nP43NRET2x1Kw\nqUhKYfBgIDlZn4IiIuqJpWBTkZQCoC9St7QYl4eIogNLwaZYCkRkBpaCDQmhl8K4cQN/D5YCEfWF\npWBDHR3AoEHA8OEDf4/x41kKRHQrloINnTsX2dQRwJECEfWNpWBDka4nACwFIuobS8GGWApEZBbT\nSqGpqQn3338/pk6dimnTpmHTpk0AgIsXL6K0tBTZ2dkoKytDe3u7WRGiFkuBiMxiWinExcXhxRdf\nxKFDh/Dhhx/it7/9LY4cOYLKykqUlpaivr4ec+bMQWVlpVkRopYRpTByJNDVpS9aExH5mVYKycnJ\nKCgoAAAkJCQgNzcXzc3N2LlzJ8rLywEA5eXl2LFjh1kRopYRpaBp+miBZzUTUU+WrCmcPHkStbW1\nKCwsRFtbGxxfXcnN4XCgra3NighRxYhSADiFRES3GmL2Djo6OvDII49g48aNSExM7PU7TdOgBbl1\nmMfjCfxcUlKCkpISE1PaC0uBiADA6/XC6/Ua+p6mlkJXVxceeeQRLFmyBAsXLgSgjw58Ph+Sk5PR\n2tqKu4J8uvUsBeot0rOZ/ZKTAQ7UiOzr5r8wr127NuL3NG36SAiBiooKTJkyBU899VTg+QULFmDr\n1q0AgK1btwbKgkJn1EjB4QB8vsjfh4iih2mlUFNTg9deew179uyB2+2G2+1GdXU1Vq9ejV27diE7\nOxu7d+/G6tWrzYoQlbq7gQsXgLFjI3+v5GSWAhH1Ztr0UVFREbq7u/v83TvvvGPWbqPepUtAYiIQ\nHx/5e3H6iIhuxjOabcaoqSOAIwUiuhVLwWaMLAWuKRDRzVgKNmNkKYwbB5w/D9y4Ycz7EZH9sRRs\nxshSiIsDkpL0hWsiIoClYDtGlgLAdQUi6o2lYDNGlwLXFYioJ5aCzXCkQERmYinYjFGXuPBjKRBR\nTywFmzHi/sw98QQ2IuqJpWAzXFMgIjOxFGyksxO4fFn/GqlROH1ERD2xFGzk/Hn9QniDDPy3xlIg\nop5YCjZi9NQRwDUFIuqNpWAjZpTC6NHAF1/oU1NERCwFGzGjFAYN0r/ievasse9LRPbEUrARM0oB\n4LoCEX2NpWAjZpYC1xWICGAp2IrRZzP78VwFIvJjKdiI0Wcz+3H6iIj8WAo2wjUFIjIbS8FGuKZA\nRGZjKdiIWaXANQUi8mMp2MSVK4AQwPDhxr83p4+IyI+lYBP+UYKmGf/eLAUi8mMp2IRZU0cAMGIE\n0NUFXL1qzvsTkX2wFGzCzFLQNC42E5GOpWATZpYCwMVmItKxFGzCrLOZ/biuQEQAS8E2zDqb2Y+l\nQEQAS8E2zJ4+4poCEQEsBdvgmgIRWYGlYBNWjBRYCkTEUrAJlgIRWYGlYAPd3fpCs9nfPuKaAhGx\nFGygvR1ISADi483bh39NQQjz9kFE6jO1FJYtWwaHw4G8vLzAcx6PB06nE263G263G9XV1WZGiApm\nTx0BwJ13AnFxwBdfmLsfIlKbqaWwdOnSWz70NU3DihUrUFtbi9raWnznO98xM0JUsKIUAK4rEJHJ\npVBcXIykpKRbnhecowiL2esJflxXICIpawqbN2+Gy+VCRUUF2tvbZUSwFatGCjxXgYgsL4Xly5ej\nsbERdXV1SElJwcqVK62OYDucPiIiqwyxeod39fh0e/TRRzF//vw+X+fxeAI/l5SUoKSkxORk6jp7\nFsjNNX8/qalAc7P5+yEiY3i9Xni9XkPf0/JSaG1tRUpKCgBg+/btvb6Z1FPPUoh1Z88C3/qW+fuZ\nMAH45BPz90NExrj5L8xr166N+D1NLYXFixfj3Xffxfnz55GWloa1a9fC6/Wirq4OmqZh4sSJeOml\nl8yMEBWsmj5KTwdOnjR/P0SkLlNLoaqq6pbnli1bZuYuoxJLgYiswjOabcDKheb2duDLL83fFxGp\niaWguK4u/Szj0aPN39egQcDddwOnT5u/LyJSE0tBcefPA2PG6B/YVpgwgVNIRLGMpaA4q85m9uO6\nAlFsYykozqr1BL/0dODUKev2R0RqYSkorq1Nv/yEVThSIIptLAXFyRgpnDhh3f6ISC0sBcW1tVlb\nCtnZwPHj1u2PiNTCUlDc2bPWTh+NHavf/vPCBev2SUTqYCkozuqRgqYBWVlAfb11+yQidbAUFGf1\nSAHQp5BYCkSxiaWgOKsXmgGWAlEsYykoTAjrp48ALjYTxbLbXiW1vb0d+/btw8mTJ6FpGtLT0zF7\n9myMHDnSinwx7fJlIC4OuPNOa/fLNQWi2BV0pLB3714sWLAA9913H37/+9/j9OnTOHnyJKqqqlBc\nXIwFCxbg/ffftzJrzLH6xDW/rCx9pNDdbf2+iUiuoCOF7du3Y8OGDcjKyurz9/X19diyZQuKiopM\nCxfrZKwnAMCIEcDIkcCZM/pVU4kodgQthRdeeKHfDbOzs2/7GoqMjG8e+eXmAkePshSIYk3QUtiw\nYUPQjTRNw4oVK0wJRF+Tscjsl5MDHDkClJXJ2T8RyRG0FC5fvgxN06zMQjeRNX0E6KVw+LCcfROR\nPEFLwePxWBiD+tLWpk/jyJCbC/zhD3L2TUTyBP32kcfjQVtbW9ANW1tbsWbNGlNCkU72SOHIETn7\nJiJ5go4UZsyYgUWLFqGzsxP33HMPUlJSIISAz+fDxx9/jKFDh+InP/mJlVljjqyvpALA+PHAlSvA\npUtAUpKcDERkvaAjhW3btmHPnj2YO3cuioqKMHjwYMTFxaGoqAjbtm3D7t27MW/ePCuzxhyZIwVN\n00cLR4/K2T8RyRF0pHDgwAG0tLTgjTfegNfrhRAi8DsuQFtD5ldSga9LYfZseRmIyFpBS+Hxxx/H\nnDlzcOLECUyfPr3X7zRNwwnenstUnZ1ARwcwapS8DBwpEMWeoNNHTz75JI4cOYKlS5eisbGx14OF\nYL5z54Bx44BBEi9ZmJvLxWaiWHPbj5wtW7ZYkYNuIvPENT+OFIhiDy+drSjZ6wkAMGkScPo0cO2a\n3BxEZB2WgqJUGCnExwMTJgCffSY3BxFZh6WgqNZWICVFdgpOIRHFGpaColgKRCQDS0FRLS1Aaqrs\nFPwGElGsYSkoiiMFIpKBpaAoVUYK/lLgrTmJYgNLQUFCqDNSGDUKSEwEmptlJyEiK7AUFNTern8d\n9M47ZSfRcQqJKHaYWgrLli2Dw+FAXl5e4LmLFy+itLQU2dnZKCsrQ3t7u5kRbKm1VY2pIz8uNhPF\nDlNLYenSpaiuru71XGVlJUpLS1FfX485c+agsrLSzAi21NKixtSRH0cKRLHD1FIoLi5G0k13aNm5\ncyfKy8sBAOXl5dixY4eZEWxJtZEC78JGFDssX1Noa2uD46uL+jgcjn5v+RmrVBspTJoE8MK4RLFB\n6kKzpmm8YU8fVBsppKUBPp9+jwciim5Bb7JjFofDAZ/Ph+TkZLS2tuKuIFd983g8gZ9LSkpQUlJi\nTUAFtLQAs2bJTvG1uDj9ns2nTgFZWbLTEJGf1+uF1+s19D0tL4UFCxZg69ateOaZZ7B161YsXLiw\nz9f1LIVYo8o5Cj1lZOhTSCwFInXc/BfmtWvXRvyepk4fLV68GPfeey+OHTuGtLQ0/O53v8Pq1aux\na9cuZGdnY/fu3Vi9erWZEWxJlbOZe/KXAhFFN1NHClVVVX0+/84775i5W1tT6WzmnlgKRLGBZzQr\npr0dGDIESEiQnaQ3lgJRbGApKKapCbj7btkpbpWRATQ2yk5BRGZjKSjm9Gl1S6GhQZ/eIqLoxVJQ\njKqlkJQEaBpw6ZLsJERkJpaCYlQtBU3jugJRLGApKOb0af0MYhWxFIiiH0tBMaouNAMsBaJYwFJQ\njKrTRwBLgSgWsBQUcv26fuLa+PGyk/SNpUAU/VgKCmltBcaN02/FqSL/11KJKHqxFBSi8tQRoGdr\naQG6umQnISKzsBQUonopxMfr12RqapKdhIjMwlJQiOqlAHBdgSjasRQUovLXUf1YCkTRjaWgEJVP\nXPNjKRBFN5aCQk6cACZOlJ2ifywFoujGUlCEEPqlqTMyZCfpH0uBKLqxFBTR1gYMHw4kJspO0j+W\nAlF0Yyko4sQJ9UcJADBmjH7mNS+hTRSdWAqKsEsp+C+hzbuwEUUnloIiGhrsUQoAp5CIohlLQREn\nTgCZmbJThIalQBS9WAqKsMv0EcBSIIpmLAVFsBSISAUsBQV8+SVw8SKQmio7SWhYCkTRi6WggMZG\nYMIEYPBg2UlCM2GCfp2m69dlJyEio7EUFGCnqSMAGDoUcDiAM2dkJyEio7EUFGCnr6P6cQqJKDqx\nFBRQXw9kZ8tOEZ6MDOCzz2SnICKjsRQUYMdSmDxZz01E0YWloIDjx+1XCjk5wNGjslMQkdFYCpJ9\n+SXg86l/x7WbsRSIohNLQbKGBv3GOkOGyE4SnowM/dtH167JTkJERmIpSGbHqSMAiIsD0tO52EwU\nbVgKktXXA1lZslMMDKeQiKIPS0EyO37zyG/yZODYMdkpiMhI0may09PTMWLECAwePBhxcXHYv3+/\nrChSHT8OLFkiO8XA5OQAe/bITkFERpJWCpqmwev1YvTo0bIiKMHu00f/+Z+yUxCRkaROHwkhZO5e\nus8/Bzo67HN11Jv5p49i/F8jUVSRVgqapuGBBx7AjBkz8PLLL8uKIdXx4/ooQdNkJxmY0aOBO+4A\nmptlJyEio0ibPqqpqUFKSgrOnTuH0tJS5OTkoLi4OPB7j8cT+LmkpAQlJSXWhzSZvxTsLD8f+OQT\nwOmUnYQo9ni9Xni9XkPfUxMKzOGsXbsWCQkJWLlyJQB9FKFALNN5PPo9CX7xC2PeT1urQayx9rit\nWgUkJQE/+5mluyWSQ9OUni814rNTyvTR1atXcfnyZQDAlStX8Je//AV5eXkyokh15AiQmys7RWRc\nLuDgQdkpiMgoUkqhra0NxcXFKCgoQGFhIR566CGUlZXJiCLV4cPAlCmyU0TGP31ERNFBiemjm8XC\n9NH160BiInDhAnDnnca8p4zpo85OYNQo/Z/jjjss3TWR9Th9RGZpaNC/impUIcgSH6+fkf3pp7KT\nEJERWAqSRMPUkZ/LxSkkomjBUpAkGhaZ/VgKRNGDpSBJNI0UCgqA2lrZKYjICCwFSaKpFKZPB+rq\ngK4u2UmIKFIsBQlu3NDvQxAt00cjR+q3Ez10SHYSIooUS0GCU6eAMWP0r6RGi5kzgb/9TXYKIooU\nS0GCTz7RF2ejSWEhEKO3xCCKKiwFCerq9MXZaMKRAlF0YClIUFsLuN2yUxgrPx9obAS+uqQVEdkU\nS0GCaBwpxMXpU2IHDshOQkSRYClY7MIF/Y5rEyfKTmK8WbOAmhrZKYgoEiwFi9XV6X+jHhSFR/6+\n+4C9e2WnIKJIROFHk9oOHIi+9QS/oiLggw/0K8ASkT2xFCy2bx9w772yU5hj7Fj9JLa6OtlJiGig\nWAoWEkL/m3S0lgKgTyG9957sFEQ0UCwFC504oX9LJy1NdhLzsBSI7I2lYKFonjryKy7WF5u7u2Un\nIaKBYClYqKYm+kth/Hh9bYHrCkT2xFKw0DvvACUlslOYb+5c4O23ZacgooFgKVjk+HHgypXouxBe\nX1gKRPbFUrDI22/rH5aaJjuJ+b71LeDgQeDSJdlJiChcLAWL+EshFgwbpi8479olOwkRhYulYIHP\nP9fPT3jgAdlJrMMpJCJ7YilY4H//F/j2t4FRo2Qnsc78+cD//R/v20xkNywFC/z3fwM//KHsFNaa\nMAGYPJlTSER2w1Iw2aefAseOAQsWyE5ivcWLgaoq2SmIKBwsBZNt2AAsXw7Ex8tOYr3vfx/44x95\nNzYiO2EpmOjQIeCtt4Ann5SdRA6HAygtBV59VXYSIgoVS8EkQgBPPw389KfAyJGy08izYgXw618D\nN27ITkJEoWApmGTrVqCpSZ86imWzZgGpqcDrr8tOQkShYCmY4NAhYNUq/YNw6FDZaeR7/nngZz8D\n2ttlJyGi22EpGOz0af3ErV//GsjLk51GDbNnAw8/DPzoR7xVJ5HqWAoG+ugj/dLYTz8de+cl3M7z\nzwNXr+rfSLpwQXYaIgqGpWCA1lbg0UeBhx4CNm0CnnhCdiL1DB0K7Nypn9SWlQX8278Bf/ubviBP\nROqQUgrV1dXIyclBVlYW1q9fLyOCIa5eBX7xC32aaMwYoL4e+Jd/kZ1KXcOG6dNqBw/qx+tf/xWY\nNAl49lmguVl2OiICJJTCjRs38MQTT6C6uhqHDx9GVVUVjhw5YnWMiB06BKSne/GPf+jTRuvXq/vV\nU6/XKztCL04n4PEAR48C//M/QEcHMH068Mc/emVHC4lqxzMY5jSOHTIaxfJS2L9/PyZNmoT09HTE\nxcVh0aJFePPNN62OEbHMTGD+fC+2bQMmTpSdpn+q/getacA99wAbNwINDcDf/+6VHSkkqh7PmzGn\nceyQ0SiWl0JzczPS0tICf3Y6nWi24dzBsGFAj38MitDw4bITEBEgoRS0WLj1GBGRXQmL7du3Tzz4\n4IOBP69bt05UVlb2ek1mZqYAwAcffPDBRxiPzMzMiD+jNSGs/VLg9evXMXnyZPz1r39FamoqZs6c\niaqqKuTm5loZg4iI+jDE8h0OGYLf/OY3ePDBB3Hjxg1UVFSwEIiIFGH5SIGIiNQl7YzmixcvorS0\nFNnZ2SgrK0N7kKulLVu2DA6HA3k3XUgo1O2tyhnshDyPxwOn0wm32w23243q6mpD84VyIuCTTz6J\nrKwsuFwu1NbWhrWtCjnT09ORn58Pt9uNmTNnSst49OhRzJ49G8OGDcOGDRvC2laVnFYdy1Byvv76\n63C5XMjPz8c3v/lNHDx4MORtVcmp0vF888034XK54Ha7MX36dOzevTvkbXuJeFVigFatWiXWr18v\nhBCisrJSPPPMM32+7r333hMff/yxmDZt2oC2tyLn9evXRWZmpmhsbBSdnZ3C5XKJw4cPCyGE8Hg8\nYsOGDaZk62+/fn/605/E3LlzhRBCfPjhh6KwsDDkbVXIKYQQ6enp4sKFC6ZkCyfj2bNnxUcffSSe\nffZZ8fzzz4e1rQo5hbDmWIaa84MPPhDt7e1CCCHefvttZf/bDJZTCLWOZ0dHR+DngwcPBhadwz2e\n0kYKO3fuRHl5OQCgvLwcO3bs6PN1xcXFSEpKGvD2VuS83Ql5wqQZulBOBOyZv7CwEO3t7fD5fJae\nRDjQnG1tbYHfm3UMw8k4btw4zJgxA3FxcWFvq0JOP7OPZag5Z8+ejZFfXQagsLAQZ86cCXlbFXL6\nqXI8h/c42aejowNjx44NeduepJVCW1sbHA4HAMDhcPT6ALBieyP3c7sT8jZv3gyXy4WKigpDp7lC\nOREw2GtaWlosO4kwkpyAfm7LAw88gBkzZuDll1+WltGMbcMV6b6sOJZA+DlfeeUVzJs3b0DbysoJ\nqHc8d+zYgdzcXMydOxebNm0Ka1s/U799VFpaCp/Pd8vzv/zlL3v9WdO0iE5qi3T7SHP2t+/ly5fj\n5z//OQDgueeew8qVK/HKK68MOGuo++3Jir/J9CfSnO+//z5SU1Nx7tw5lJaWIicnB8XFxUZGjPi/\nP6tEuq+amhqkpKSYeiyB8HLu2bMHr776KmpqasLeNlKR5ATUO54LFy7EwoULsXfvXixZsgRHjx4N\ne1+mlsKuXbuC/s7hcMDn8yE5ORmtra246667wnrvSLc3Muf48ePR1NQU+HNTUxOcTicA9Hr9o48+\nivnz5w84Zzj7DfaaM2fOwOl0oqur67bbys45fvx4AEBqaioAfVrk4Ycfxv79+w3/Hy+UjGZsG65I\n95WSkgLA3GMZTs6DBw/iscceQ3V1dWCaWMXj2VdOQL3j6VdcXIzr16/j4sWLcDqd4R1PoxdEQrVq\n1arAmcy/+tWv+l0obmxs7HOhOdTtzc7Z1dUlMjIyRGNjo7h27VqvhZyWlpbA61544QWxePFiw7L1\nt1+/ngu4+/btCyyShbKtCjmvXLkivvjiCyGEvpB27733ij//+c9SMvqtWbOm1wKuascyWE6rjmWo\nOU+dOiUyMzPFvn37wt5WhZyqHc/PPvtMdHd3CyGEOHDggMjIyAh5256klcKFCxfEnDlzRFZWligt\nLRWXLl0SQgjR3Nws5s2bF3jdokWLREpKioiPjxdOp1O8+uqr/W4vK+dbb70lsrOzRWZmpli3bl3g\n+SVLloi8vDyRn58vvvvd7wqfz2dovr72u2XLFrFly5bAa3784x+LzMxMkZ+fLw4cOHDbzGYYaM6G\nhgbhcrmEy+USU6dONTXn7TK2trYKp9MpRowYIUaNGiXS0tLE5cuXg26rWk4rj2UoOSsqKsTo0aNF\nQUGBKCgoEN/4xjf63Va1nKodz/Xr14upU6eKgoICUVRUJPbv39/vtsHw5DUiIgrg7TiJiCiApUBE\nRAEsBSIiCmApEBFRAEuBiIgCWApERBTAUiAKQ0JCguwIRKZiKRCFwcrr8hDJwFIgIqIAlgIREQWw\nFIiIKIClQEREASwFIiIKYCkQhYHfPqJox0tnExFRAEcKREQUwFIgIqIAlgIREQWwFIiIKIClQERE\nASwFIiIKYCkQEVEAS4GIiAL+H/blWAdYqt6dAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10cba6350>"
       ]
      }
     ],
     "prompt_number": 90
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 92
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}