07_global_spatial_autocorrelation.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:efc5d6f4fd8e0584edbd84d495213777be994d0de54be518c524f4d159618e51"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "%pylab inline"
     ],
     "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": [
      "y = np.arange(25)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y.shape = (5,5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "matshow(y,cmap=cm.gray_r)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "<matplotlib.colorbar.Colorbar instance at 0x108491950>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAPkAAADyCAYAAACPiGNTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE19JREFUeJzt3XmMHvV9x/H37hrH8cVhXK937bJ/1AegGoONsXCoN1wy\nkEIiRQFUVI4KUATliMoVItVthZqStiCuCAJEEMIRQbhNOVJKTSG+Yi84YGwOW9hrG4OJ8e5isbaf\n/vGd3X32YeZ55vnN7M48v+fzkkb7HPPM/Nbrz/Ob+c3Md0BERERERERERERERERERCSxQhXTzoza\nCEBDlisXqWGFhoZ48SkUCpBh1kZktWKRWldlyDOjkIs4ihvyrCnkIo4aGxtjzbdv374hbkl5CrmI\nI/XkIp5TyEU8p5CLeK5WQh5v5GDoLQLWARuA6zJuy/3AduDtjNsBMBV4FfgjsBa4ItvmMApYBqwB\n3gH+Ndvm9GsCVgPPDudKGxoaYk0hov6uPwPeBTqA3wIHRqx6I/AW9jsvT+83GjpNwPtAG3AA9h/o\n8AzbcwJwNPkIeTMwO3g8FniPbP9tAEYHP0cAvwe+lWFb+vwI+DXwzDCuszBmzJhYE3bWW7Gov+sp\nDHS8Pw2mMB8Bh8RtaB568nlYyDcCvcCjwFkZtmcp8HmG6y+2DfvSA+jCvuVbsmsOAD3Bz5HYF3Sm\np2wCU4DTgXsZ5rPKGhsbY00hov6uLwP7g9eXYb9blNi/ax5C3gp8XPR8c/CaDNaGbWEsy7gdjdh/\n0O3YJuc72TaHW4BrGAjHsEmwuV6sjfC/60XAkojPFIBXgJXAxZVWkIeQZ3vOX20YCzwOXIl982dp\nP7apOQX4K6A9w7Z8B/gE2zcd9lGwFEIe9Xe9EfgKeDjicwuwL4bTgMuwXcxIeQj5Fmwgos9UrDcX\ncwDwBPAQ8FTGbSm2C3gemJthG44HzsT2UR8BTgQeHK6VR4V679697Nmzp3+KEPV3vQDb/fibMqve\nGvzcATyJ7fLm2gjgA2yzZSTZD7yBtSUPA28N2H/aW7JuSOBQ4KDg8TeB/wVOyq45gyxkeEfXCwcf\nfHCsia9vrUb9XRdhI+6HllnvaGBc8HgM8H/AqeUamoeefC9wOfAitn/3GDYQkZVHgDeA6dhYwYUZ\ntmUBcB7wbWyTdDX2HyErk4H/xr6Il2Gh+l2G7Sk1rLt+CTbXw/6upwG3Y5vwLwev3RXM34JtNYGN\nzC9l4G/wHPBS2XYm+SVF6lhhwoQJsWb87LPPQNeTi9SeuFehZU0hF3FUK6e1KuQijhRyEc/VSsgT\n71QsXLiwmqqVmjTldpo7d26BKqR0xtuQS6MFhTQK1S1evJjFixcnXk5PT0/lmWK46aabuPHGGxMv\np7u7O4XWwM0338y1116beDlptOfWW2/lqquuSrwcSOfvdeedd3LZZZclXs6RRx4J8TNRaGmJdxlB\nZ2dnNctNnTbXRRzloZeOQyEXcaRDaFVqb2/PugmDnHBC2XP+h92CBQuybkK/+fPnZ92EQY499thM\n1lsrPXlu9snTktY+eVrS2idPS97ak6e/V7X75IcddlisGTdt2lTNclOXm55cpNbUSk+ukIs4UshF\nPKeQi3iuVkbXa6OVIjk0BCWZD8GuJV+PXSN+UNiHqbKEuUIu4ihByHuBq4EjgflYnbbDgeuxkE/H\ninFcH/LZJuAOLOhHAOdSoZKSQi7iKEHIw0oyt2L16h4IXn8A+G7IZ6suYR4n5Hm6u4lIbgxBSeZJ\nWKlrgp+TQuavuoR5pYG3vk2Dk7Gqqiuwu1RkWYNNJBdSGF0fi1VsvRLYXfJe39Vxpao+86xSyIs3\nDWBg00Ahl7oXFfLu7u44Zxb2lWT+FQMlmbdjhRq3YUUzPwn5XNUlzCttruvuJiIRom6LNG7cOJqb\nm/unEA3AfVh14luLXn8GOD94fD7hdfZXAtMYKGF+NhXuAVcp5Pk5KV0kZ1IuybwIu8HhKdghtBMZ\nuOFhcUnmqkuYV9pcj7VpUFzsob29PXdXlImEWb58OStWrHD+fIJ98teJ7mBPDnmtEzij6PkLwRRL\npVaOwG6relKwouXYcbnibw5dhVZG3q76ylt78vT3qvYqtKOOOirWjB0dHdUsN3WVevLiTYMmbD9C\ng24i+HXuelWbBiL1wqeQi0gIhVzEc7VyFZpCLuJIPbmI5xRyEc8p5CKeU8hFPKeQi3hOIRfxnA6h\niXhOPbmI5xRyEc8p5CKeSxDy+7Hrwz8B/jJ47VFgRvD4IOBPWIHHUhuBL4B9WLXWeZVWlkrIk1x4\nnzZdL11e3v598taeaiQI+S+B24EHi147p+jxv2MhD1MA2oGdcVemnlzEUYKQL8VqtIUuFvgBVhoq\nctXVrKw2jgGI5FBUIcfSqUonYFVbP4h4vwC8ghV0vDjOAtWTiziK6sl37tzJ559/7rrYc4GHy7y/\nANgKTMRuqbQO2zKIpJCLOIoK+YQJE5gwYUL/8w8//DDuIkcA3wOOKTPP1uDnDuBJbOCtbMi1uS7i\nKKXbJBU7Gauh2Bnx/mhgXPB4DHAq8HalhSrkIo4ShPwR4A3s7qUfAxcGr58dvFesuOZ6M9Zrr8Hu\nnfYcdovjsrS5LuIowej6uRGvXxjyWnHN9Q+B2dWuTCEXcaQLVEQ8p9NaRTynkIt4TiEX8ZxCLuK5\nWgl5nOHB+7FzaSsedBepJ0NwMsyQiBPyX2I3SBeRIkN0gUrq4myul7ssTqRu5aGXjkP75CKO6irk\n99xzT//jOXPmMGfOnDQWKzKk1q1bx7p165w/X1chv+SSS9JYjMiwmjlzJjNnzux//vTTT1f1+boK\nuUg9qpWQxxn6i7osTqSu+XQI7VzsmtZvAFOxQ2oidS/BIbSwc08WA5uB1cEUddh6EVbyaQNwXax2\nxvt1RKRUgp487NyTAvCfWK31o4H/CvlcE3BH8NkjsA748ErtVMhFHCUI+VIgrNJjpW37ecD72A0W\nerEbMpxVqZ0KuYijIdgn/3ugA7gPu4tKqVZsXKzP5uC1sjS6LuIoKsDbtm1j27Zt1S7u58A/B4//\nBfgP4O9K5ilUu1BQyEWcRYV88uTJTJ48uf95R0dHnMV9UvT4XuDZkHm2YIPffaZivXlZ2lwXcZTy\n5vrkosffI/yqz5XANOxakpFYdddnKi1YPbmIowRXmD0CLAQOxfax/xG7ieFsbJP8I+DSYN4W4BdY\nxda9wOXAi9hI+31YnfayFHIRRymXZL4/Yt7ikswALwRTbAq5iKM8nM0Wh0Iu4kghF/GcQi7iOYVc\nxHN1FfIVK1aksZhUdHd3Z92EQfLWnp6enqybMEhXV1fWTXCWhyKNcagnF3FUVz25SD1SyEU8p5CL\neE4hF/GcQi7iOY2ui3hOPbmI52ol5LWxvSGSQwmKRoSVZP4Zdm14B/Bb4MCI1W4E3sLKNi+P006F\nXMRRyiWZXwKOBI4C1gM3RKy2gBWYOBqr3lqRQi7iKOWSzC8D+4PHy4Ap5VZdTTsVchFHQ3ibpIuA\nJRHvFYBXsHpvF8dZmAbeRBwN0SG0G4GvgIcj3l8AbAUmYr3/OmzLIFKckE8FHgT+DPsWuQe4LV57\nRfwV1Utv2rSJTZs2uSzyAuB04KQy82wNfu4AnsT2yxOHvBe4GlgDjAVWYd8gFatEivgsKuRtbW20\ntbX1P1+6tGwG+ywCrsGquO6JmGc0VqV1NzAGOBX4p0oLjrO9sQ0LOEAXFu6WGJ8T8VqCffK+24HP\nwEoyXwTcjnWiL2OHx+4K5m0Bng8eN2O99hpscO45bFS+rGr3yduwoftlVX5OxDsZlGT+EKvNXpVq\nQj4WeBy4EuvR+z3//PP9j6dNm8b06dOrbYfIsNuyZQudnZ3On6+VM97ihvwA4AngIeCp0jfPOOOM\nr31AJO9aW1tpbR24KejKlSur+rxPIW/AbsfyDnDr0DZHpHbUylVocVq5ADgP+DY2ILCar5+SJ1J3\nhvBkmFTF6clfR2fGiXxNHgIch854E3GkkIt4TiEX8ZxCLuI5hVzEc7VyCE0hF3GknlzEcwq5iOcU\nchHP1UrIa2PkQCSHUi7JfAh2Lfl67BrxgyJWuwgr+bQBuC5OOxVyEUcpl2S+Hgv5dOB3wfNSTcAd\nwWePwK5LP7xSOxVyEUeNjY2xphBhJZnPBB4IHj8AfDfkc/OA97EbLPQCjwJnVWpnKvvky5fHupHD\nsOjp6cm6CYN0dXVVnmkY5e3fJ2/tqUbK++STsE14gp+TQuZpxcpF9dkMHFdpwRp4E3E0hANvhWAK\ne71qCrmIo6iQr1+/ng0bNlS7uO1YocZtwGTgk5B5tmAl0vtMxXrzshRyEUdRIZ8xYwYzZszof75k\nSdTNUAZ5Bjgf+Lfg59fKrGF3TZmGFVTtBM4mvCjkIBp4E3GUYknmC4GfAqdgh9BODJ7D4JLMe4HL\ngRexcmyPEeP+B+rJRRwluEAlqvc9OeS14pLMAC8EU2wKuYijWjnjTSEXcaSQi3hOIRfxnEIu4jmF\nXMRzCrmI51TjTcRz6slFPOdTyEcBrwHfAEYCTwM3DGWjRGqBTyHfg93RtCeY/3XgW8FPkbrlU8jB\nAg7WkzcBO4emOSK1o1ZCHnd4sBFYg13z+ip2BYxIXfPp/uQA+4HZwIHYZW7twP/0vblmzZr+GZub\nm2lubk6tgSJD5YsvvmD37t3On/f1ENou7NrWuRSFfPbs2Sk2SWR4jB8/nvHjx/c/7+zsrOrzeeil\n44jzVXQoAzWgv4ld2L56yFokUiMSbK7PwDLUN+0CriiZpz14vW+en7i2M05PPhkrEdsYTL/C6kKL\n1LUEPfl7wNHB40asdtuTIfO9hpVqTiROyN8Gjkm6IhHfpLS5fjLwAYNLLfevIo0V1MbIgUgOpTS6\nfg7wcMjrBeB4oANYgt0xxYlOaxVxFBXgtWvXsnbt2jiLGAn8NeH3NPsDVnK5BzgNq9463aWdCrmI\no6hDaLNmzWLWrFn9zx977LGoRZwGrAJ2hLxXfGzvBeAu7KaIVZ+IppCLOEphn/xcrDxzmEnYDRYK\n2D3QGnA801QhF3GUMORjsEG3i4teuzT4eTfwfeCHWK31Hmzf3YlCLuIoYci7sXNQit1d9PjOYEpM\nIRdxVCtnvCnkIo4UchHPKeQinvP1KjQRCagnF/FcXYV8+fLlaSwmFT09PZVnGkZ5a09XV1fWTRjk\nyy+/zLoJzuoq5CL1SCEX8ZxCLuI5ja6LeE49uYjnFHIRzynkIp5TyEU8lzDkG4EvgH1AL1YYotRt\nWPWYHuACHEuhK+QijhKGvIDVVo+q9nI68BfANOA44OfAfJcVKeQijlI4hFbuW+JM7H4HAMuwG5xM\nwu5HWJXaONAnkkMJSzIXgFeAlQwuAdWnlcG12DcDU1zaqZ5cxFFUgFetWsWqVasqfXwBsBWYCLwM\nrAOWlq6i5Hmh+lYq5CLOokI+d+5c5s6d2//83nvvDZtta/BzB3aLpHkMDvkWrO56nynBa1WLu7ne\nhI3sPeuyEhEfJdhcHw2MCx6PAU7FbkdW7Bngb4PH84E/4bA/DvF78iuBd4oaJlL3EoyuT2LgBocj\ngF8DLzG4JPMSbIT9fayy64WuK4sT8inBym4CfuS6IhHfJAj5R8DskNfvLnl+uesKisUJ+S3ANcD4\nSjOK1JNauQqtUiu/g92qZTUp3UZVxBcp3dV0yFXqyY/HDsqfDozCevMHGRgQAGDHjoH7tY0ePZox\nY8ak20qRHMpDgOOoFPIfBxPAQuAfKAk4wMSJE1Nulkj++RLyUk4H40V85GPIXwsmEcHPkItIEYVc\nxHO1cghNIRdxpJ5cxHMKuYjnFHIRzynkIp5TyEU8Vyshr41jACI51NjYGGsKMRV4FfgjsBa4ImSe\ndmAXdnHYauAnru3MTU/e3d2dqwtb9uzZw6hRo7JuRr+vvvqKkSNHZt0MAPbt20dTU1PWzchcgp68\nF7gaWAOMBVZhdd7eLZnvNewCsURy05P39PRk3YRB9uzZk3UTBunt7c26Cf3279+fdRNyIcGlptuw\ngAN0YeFuCVtFGu3MTchFak1K15O3AUdjtdWLFbBLvTuwUlBHuLYzN5vrIrUmhYG3scDjWA3FrpL3\n/oDtu/dgt0p6CpjuspI0Qv7au+++uzCF5fDpp5+msZjU7Nq1K+smDJKnXZq9e/dm3YShUNVVllEh\nf+ONN3jzzTcrffwA4AngISzApXYXPX4BuAs4hOjbKkW3s9oPiAgAhc7OzlgztrS0wOCsNWC3QPoM\nG4ALMwkrvVbAarL/Btu0r5o210UcJbgKbQFwHvAWA3cq/THw58Hju4HvAz8E9mKb7Oe4rkw9uYib\nwvbt8e51MGnSJMgwa+rJRRzVyhlvCrmII4VcxHMKuYjnFHIRz6nGm4jn1JOLeE4hF/GcQi7iOYVc\nxHMKuYjnFHIRz+kQmojn1JOLeK5WQl4b2xsiOZSwxtsiYB2wAbguYp7bgvc7sDpwIjKMCr29vbEm\nrLpLsSbgfazSywFY5dbDS+Y5HSvgCHAc8HvXhqonF3GUoCefh4V8I1aD/VHgrJJ5zsRKRIFVcj0I\nKwlVNYVcxFGCkLcCHxc93xy8VmmeKS7t1MCbiKMEh9BKN9+jlH5DxP3cIAq5iKMqRtd3lzzfgtVU\n7zMV66nLzTMleE1EasAI4ANs4G0klQfe5pNg4E1EsnEa8B42AHdD8NqlwdTnjuD9DuCYYW2diIiI\niIiIiIiIiIiIiIiIiEge/T/qN7t+EnVvJAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1083f0590>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "array([[ 0,  1,  2,  3,  4],\n",
        "       [ 5,  6,  7,  8,  9],\n",
        "       [10, 11, 12, 13, 14],\n",
        "       [15, 16, 17, 18, 19],\n",
        "       [20, 21, 22, 23, 24]])"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = np.zeros((5,5))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Perfect negative autocorrelation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = np.zeros((25,))\n",
      "ids = range(0,25,2)\n",
      "y[ids] = 1\n",
      "y.shape = (5,5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "matshow(y,cmap=cm.gray)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<matplotlib.image.AxesImage at 0x10878c890>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAPYAAAD7CAYAAABZjGkWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACC1JREFUeJzt3d+LXGcdx/H3ZjchlSC96EWtCeRGwQuhhSIisXusCrXU\nH3fihULvpUGxSK4yV5ZexT+gVmpb4kULYi3qhWRsRQz+SOiPVLGiUJWm0JZCCYol24szSWPiOmcy\n5zxnns+8XzDsJDvsPvkm7z1zZrPPAUmSJEmSJEmSJEk9uwv4I/Bn4Nsjr6WLh4HzwPNjL2QBh4BT\nwIvAC8B94y5nrv3AaeAscA54YNzlLGQTOAM8NfZCxrQJvAwcBvbS/kV+ZMwFdfBJ4DbqCvtm4NbZ\n/QPAn1j9Ob9v9nYL+A1wZMS1LOKbwOPAj0t9wj2lPtECPkYb9t+A/wA/BL445oI6eBZ4c+xFLOhV\n2i+aAG8DLwG3jLecTi7M3u6jPQC8MeJaujoI3A08BGyU+qSrGPYHgVeu+PXfZ7+n4RymfcZxeuR1\nzLOH9ovRedrTiHPjLqeTE8D9wMWSn3QVw94ZewFr5gDwBHCU9si9yi7Snj4cBO4AmlFXM989wGu0\n59fFjtawmmH/g/aFnUsO0R611b+9wJPAY8CPRl7LIt4CngZuH3shc3wC+ALwV+AkcCfwg1FXNKIt\n4C+0Tw/3UceLZ9Cut6YXzzZo/5GdGHshHd0E3Di7fwPwDPDp8ZazsG3W/FVxgM/Rvkr7MnBs5LV0\ncRL4J/Bv2tcH7h13OZ0coX1qe5b2qeIZ2m8zrqqPAn+gXe9ztOetNdmm4KvikiRJkiRJkiSps+3t\n7R3a/y3mzZu3wrdZf4PYGcrx48cH+bhj/2Vcz20oQ814Z8c5X1JgxtdYxf9SKmlJhi0FWumwm6YZ\newnxnPHwxphxHz9KNnuqX4+NjaI/QdeL2mYMzrmE2YyvGfRKH7ElXR/DlgIZthTIsKVAhi0FMmwp\nkGFLgQxbCmTYUiDDlgIZthTIsKVAXcKu7VrV0tqb9+M3m7RX5PgM7TW1fgt8hfaSq5f4010F1DZj\ncM4lXO9Pd9V4rWpp7c0L22tVSxXamvP+Ts9LJpPJ5ftN07grhzSQ6XTKdDqd+7h5J0EfBya8dxXG\nY7RXaHzwisd4jl1AbTMG51zC9Z5j/w74EO9dq/rLeClQaeXNeyr+DvB14Oe0r5B/j/9+RVzSCnIz\nw0rUNmNwziW4maG0RgxbCmTYUiDDlgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVA\nhi0FMmwpkGFLgXrZQaWHj1FUbbtkgLuRlFLjnHEHFWk9GLYUyLClQIYtBTJsKZBhS4EMWwpk2FIg\nw5YCGbYUyLClQIYtBTJsKZBhS4EMWwpk2FIgw5YCdQn7YeA88PzAa5HUky5hfx+4a+iFSOpPl7Cf\nBd4ceiGS+uM5thTIsKVAhi0FMmwpUJewTwK/Bj4MvALcO+iKJC3NK4FUosYrVDjnYrwSiLQODFsK\nZNhSIMOWAhm2FMiwpUCGLQUybCmQYUuBDFsKZNhSIMOWAhm2FMiwpUCGLQUybCmQYUuBDFsKtNXH\nB6ltC5wat7+pbcbgnEvYbcYesaVAhi0FMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVAhi0FMmwpkGFL\ngQxbCmTYUiDDlgIZthTIsKVAXcI+BJwCXgReAO4bdEWSltZl75qbZ7ezwAHg98CXgJdm799J2U5m\nldU2Y3DOJcxmfM2guxyxX6WNGuBt2qBv6W1lknq36Dn2YeA24HT/S5HUl0V2KT0APAEcpT1yXzaZ\nTC7fb5qGpml6WJqkq02nU6bT6dzHdT0J2gv8BPgp8N2r3uc5dgG1zRiccwm7nWN3mfwG8AjwOvCN\n//F+wy6gthmDcy5hmbCPAM8AzwGX/tTHgJ/N7ht2AbXNGJxzCcuEPY9hF1DbjME5l7DMt7skVcaw\npUCGLQUybCmQYUuBDFsKZNhSIMOWAhm2FMiwpUCGLQUybCmQYUuBDFsKZNhSIMOWAhm2FKiXHVR6\n+BhF1bZLBrgbSSk1zhl3UJHWg2FLgQxbCmTYUiDDlgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUiDD\nlgIZthTIsKVAhi0FMmwpUJew9wOngbPAOeCBQVckaWlbHR7zL+BTwIXZ438FHJm9lbSCuj4VvzB7\nuw/YBN4YZjmS+tA17D20T8XPA6don5JLWlFdw74I3AocBO4AmqEWJGl5i74q/hbwNHD7AGuR1JMu\nYd8E3Di7fwPwWeDMYCuStLQur4p/AHiE9ovAHuBR4BdDLkrScrwSSCVqvEKFcy7GK4FI68CwpUCG\nLQUybCmQYUuBDFsKZNhSIMOWAhm2FMiwpUCGLQUybCmQYUuBDFsKZNhSIMOWAhm2FMiwpUBd9jyb\nq7YtcGrc/qa2GYNzLmG3GXvElgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVAhi0F\nMmwpkGFLgQxbCmTYUiDDlgJ1DXsTOAM8NeBaJPWka9hHgXNAXfvGSGuqS9gHgbuBh4D6NrGS1lCX\nsE8A9wMXB16LpJ7M26X0HuA12vPrZrcHTSaTy/ebpqFpdn2opCVMp1Om0+ncx817av0d4KvAO8B+\n4P3Ak8DXrnjMTsqWrausthmDcy5hNuNrBr3I5LeBbwGfv+r3DbuA2mYMzrmE3cJe9PvYdf2ppTXV\nx5dUj9gF1DZjcM4l9HXEllQBw5YCGbYUyLClQIYtBTJsKZBhS4EMWwpk2FIgw5YCGbYUyLClQIYt\nBTJsKdBKh91lCxgtxxkPb4wZG/aac8bDM2xJvTBsKVAfe9dMaTc6lFTeL/k/W4NLkiRJkiRJkpTv\nXaM2uz4qpEDvAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10855e1d0>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y.shape = (25,)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Join Counts"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pysal as ps\n",
      "w = ps.lat2W(5,5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$bb_i = \\sum_j y_i w_{i,j} y_j = y_i \\sum_j w_{i,j} y_j$ "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bbi = y * ps.lag_spatial(w,y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bbi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$ww_i = \\sum_j (1-y_i) w_{i,j} (1-y_j) = (1-y_i) \\sum_j w_{i,j} (1-y_j)$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "yc = 1 - y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "wwi = yc * ps.lag_spatial(w,yc)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "wwi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$bw_i = \\sum_{j} w_{i,j} - bb_i - ww_i$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$2 bw = \\sum_i bw_i = \\sum_i ( \\sum_j w_{i,j} - bb_i - ww_i ) = S_0 - bb - ww$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w.s0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "80.0"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bw = (w.s0 - bbi.sum() - wwi.sum())/2."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bw"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "40.0"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc = ps.Join_Counts(y, w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc.bw"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "40.0"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc.p_sim_bb"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "1.0"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc.p_sim_bw"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "0.001"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import gaussian_kde"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "density = gaussian_kde(jc.sim_bw)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xs = linspace(0,50,200)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(xs, density(xs))\n",
      "axvline(x=jc.bw, color='r')\n",
      "xlabel('bw')\n",
      "ylabel('f(bw)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "<matplotlib.text.Text at 0x109d60a90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUlPWd7/F3sQoiiyAgy6QJIluQpRVbkVguKBJEJ3MU\niYm5mUl0TtSgyU0ImSQ2OXMmiXE019EoGnNjEhNQkomQCxKjlNso0EIjyCI0NDZog6JgRBTo7vvH\n96mmKHqp7q6nfs/yeZ1Tp2t5qvqbJ1if/q0PiIiIiIiIiIiIiIiIiIiIiIiIiIiISI6mApuBrcCc\nBl4fAbwMfAx8q4HX2wNrgSV+FSgiIk3r4ONntwfuAy4FdgOrgcXApoxj9gG3Alc38hmzgY3AKf6V\nKSIiTWnn42dPBLYBlcARYAFwVdYx7wBl3uvZBgHTgF8CCd+qFBGRJvkZFAOBqozHu7zncnUP8G2g\nNp9FiYhIy/gZFHVteO90YC82PqHWhIiIQ36OUewGBmc8Hoy1KnJxPjAD63o6CegO/Aa4IfOgoUOH\n1lVUVLS9UhGReKkAzsj1YD9bFGXAMKAI6ATMxAazG5LdavgeFixDgOuAZ8kKCYCKigrq6up0q6vj\njjvucF5DUG46FzoXDZ4LcF5DUG7A0JZ8mfvZojgK3AIsx2ZAPYLNeLrJe30+0B+bDdUdG4uYDYwC\nPsz6rLZ0Y4mISBv4GRQAy7xbpvkZ96s5vnuqIc95NxERccDPricpoGQy6bqEwNC5OEbn4pik6wJC\nLOwziuq8/jYRkaYlEqDvCwASiQS04PtfLQoREWmSgkJERJrk92C2SMEdOACpFGzYANOnw9ixrisS\nCTe1KCRSnnsORo+GX/wC9uyBK66AWbOgVhvBiLSagkIiY8kSmDkTfvlLWL4c7r0XduyAqiq4+27X\n1YmEl2Y9SSRUV8P48bBoEUyadPxrO3fCOedYeIwf76Y+CQDNeqrX0llPCgqJhCuvtLGIf//3hl+/\n/34bt3jiiYKWJUGioKinoJDYWbMGrr4aKiqgY8eGjzlwAD71KXjjDejbt7D1SUAoKOppHYXEzvz5\ncOONjYcEQI8eFia//W3h6hKJCrUoJNQ++MBaChs3wumnN33siy/C175mxybC/i9fWk4tinpqUUis\n/P73cPHFzYcE2CD30aOwdq3/dYlEiYJCQm3RIrjhhCuVNCyRsHUVTz/tb00iUaOgkND65BNYuRJa\nskHqlCkKCpGWUlBIaK1aBSNG2EB1rpJJC5ePPvKtLJHIUVBIaK1YARdd1LL3nHIKjBsHL7zgT00i\nUaSgkNBqTVCAup9EWkpBIaH08cewejVccEHL3ztlCjzzTP5rEokqBYWE0sqV8JnPWFdSSxUX2wrt\nDz/Mf10iUaSgkFBas8Y2+muNTp1sX6jVq/Nbk0hUKSgklNavhzFjWv/+886DV17JXz0iUaagkFBa\nv966nlqrpERBIZKrQgTFVGAzsBWY08DrI4CXgY+Bb2U8PxhYAbwObAC+4W+ZEha1tbBpU36CQlv/\niDTP76BoD9yHhcUoYBYwMuuYfcCtwF1Zzx8BbgdGAyXAzQ28V2Jo+3bo0we6d2/9ZwwaBB06QGVl\n3soSiSy/g2IisA2oxL74FwBXZR3zDlDmvZ6pGij37n8IbAIG+FWohEdbu53A9n1S95NIbvwOioFA\nVcbjXd5zLVUEjAdW5qEmCbkNG9o2kJ1WUmLTbEWkaR18/vx89AB3AxYBs7GWxXFKS0vr7yeTSZIt\n2SFOQmn9ersIUVsVF8PixW3/HJGgS6VSpFKpVr/f78u3lACl2BgFwFygFvhpA8fegQXBf2Y81xH4\nC7AM+HkD79GFi2Jo5EhYuBDOOqttn7N/PwwebD/bt89PbRJgunBRvaBduKgMGIZ1HXUCZgKN/Q2X\nXXQCeATYSMMhITF0+DDs2AHDh7f9s3r2hH79bJW2iDTO76A4CtwCLMe+8Bdig9I3eTeA/tg4xu3A\n94E3se6mScAXgYuAtd5tKhJrO3fCgAHQuXN+Pq+4GF59NT+fJRJVfo9RgHUbLct6bn7G/WpszUS2\nF9GCQMmyfTt8+tP5+7x0UHzxi/n7TJGo0RexhEpFBQwdmr/PKy6GsrL8fZ5IFCkoJFTy3aKYMAHK\ny6GmJn+fKRI1CgoJlXy3KHr1gr59NaAt0hQFhYRKvoMCNKAt0hwFhYRGXV3+u55AQSHSHAWFhMbe\nvXDSSdCjR34/V0Eh0jQFhYSGH60JsKAoL7fty0XkRAoKCQ0/xifABrRPO00D2iKNUVBIaGzf7k9Q\ngLqfRJqioJDQqKjwp+sJFBQiTVFQSGjs3AlFRf58toJCpHEKCgmNqiq7hKkfJkyAtWs1oC3SEAWF\nhEJtLeze7V9QnHqqXYd761Z/Pl8kzBQUEgrvvAPdukHXrv79Dm0QKNIwBYWEQlWVXY3OTxqnEGmY\ngkJCYdcuBYWIKwoKCYVCtCg0oC3SMAWFhEIhgqJ3b7tpQFvkeAoKCYVCBAWo+0mkIQoKCQU/11Bk\nUlCInEhBIaFQqBbF2WcrKESyJVwX0EZ1dXV1rmsQn9XUQJcu8Pe/Q+fO/v6ufftsP6n334d2+jMq\nWhIJu/qVkEgkoAXf//pPQQKvutpWTvsdEmCD2aeeCtu2+f+7RMLC76CYCmwGtgJzGnh9BPAy8DHw\nrRa+V2KiUN1OaRqnEDmen0HRHrgP+8IfBcwCRmYdsw+4FbirFe+VmCjEYrtMCgqR4/kZFBOBbUAl\ncARYAFyVdcw7QJn3ekvfKzGhFoWIW34GxUCgKuPxLu85v98rEVOoqbFpxcWwZo1WaIukdfDxs9sy\nvSDn95aWltbfTyaTJJPJNvxaCaKqKpg4sXC/r3dvu472tm1w5pmF+70ifkmlUqRSqVa/38+g2A1k\ndhgMxloGeX1vZlBINBW66wmOdT8pKCQKsv+InjdvXove72fXUxkwDCgCOgEzgcWNHJs9n7cl75WI\ncxkUIuJvUBwFbgGWAxuBhcAm4CbvBtAfG4u4Hfg+8CbQrYn3SswcOWIXLRowoLC/V0EhcoxWZkug\n7dwJkybZFNlCevddGDpUK7QjRSuz62lltkSKi24nsOtn9+wJFRWF/90iQaOgkEBzFRSg7ieRNAWF\nBNquXYVdQ5FJO8mKGAWFBJpaFCLuKSgk0FwHxZo1Gv8UUVBIoLkMij59oEcPDWiLKCgk0FwGBaj7\nSQQUFBJgn3xi6xj69XNXg4JCREEhAbZrl63Ibt/eXQ3FxVBW5u73iwSBgkICy3W3E2hAWwQUFBJg\nQQiK006D7t01oC3xpqCQwApCUIDGKUQUFBJYCgqRYFBQSGAFJSgmTIDyctdViLijoJDACkpQjB1r\nQaEBbYkrBYUEVlCCYsAAqK2F6mrXlYi4oaCQQDp4EA4dsm00XEskYNw4WLfOdSUibigoJJCqqmx7\n8URArsGY7n4SiSMFhQTSrl3B6HZKU4tC4kxBIYEUlPGJNLUoJM4UFBJI6a6noBgxAnbutHETkbhR\nUEggBa1F0akTDB8OGza4rkSk8BQUEkhBCwpQ95PEl99BMRXYDGwF5jRyzL3e6+uA8RnPzwVeB9YD\nvwc6+1emBE0Qg0ID2hJXfgZFe+A+LCxGAbOAkVnHTAPOAIYBNwIPeM8XAV8DJgBjvM+6zsdaJWCC\nGBRqUUhc+RkUE4FtQCVwBFgAXJV1zAzgUe/+SqAn0A/4wHtPV6CD93O3j7VKgBw4YCuhe/Z0Xcnx\nxo6F116z2kTixM+gGAhUZTze5T2XyzHvAf8JvAm8BewH/uZbpRIo6dZEUBbbpZ16qoVXZaXrSkQK\nq4OPn53rFmoNfR0MBW7DuqAOAE8A1wOPZR9YWlpafz+ZTJJMJltWpQROELud0tLdT5/+tOtKRHKX\nSqVIpVKtfr+fQbEbyPzPfTDWYmjqmEHec0ngf4B93vN/As6nmaCQaAh6UKxbB5//vOtKRHKX/Uf0\nvHnzWvR+P7ueyrBB6iKgEzATWJx1zGLgBu9+CdbFtAfY4j3ugrU4LgU2+lirBEiQg2LcOA1oS/z4\nGRRHgVuA5diX/EJgE3CTdwNYCmzHBr3nA1/3ni8HfoOFzWvecw/5WKsESJCDIt2iEImTgA0Xtlhd\nna4mEzmXXAJz5sBll7mu5ES1tdCjh4VZ0GZlSTMSCV19ypOwmSI5f/9rZbYETpBbFO3awZgxalVI\nvDQXFBOAn2FrHPYA1d79n3H8KmqRvKirC94W49nU/SRx09Ssp6XA+9iA8wPYeoYEcDq2mO5/Ywvk\nPudzjRIj+/ZB587QrZvrSho3bhysXOm6CpHCaSoovoK1IrJt924LgL5+FCXxFeRup7SxY+EhTa2Q\nGGmq6ykdEl/Fprk2ZG9+y5G4C0NQjBkDmzbBkSOuKxEpjFwGs/8Bm7q6A1shfSswzs+iJL7CEBQn\nn2w1btniuhKRwsglKH4IXIztAPsi8B3gVT+LkvjasQOGDHFdRfM0oC1xkktQ/ABYBvwV2xL8Wxy/\n7YZI3lRWhiMotEJb4iSXoPg80BvbvfVPwJPYDCiRvNuxA4qKXFfRPLUoJE5yCYrx2F5Lq4ApwAas\nC0ok78LU9VReroW+Eg+57B47BpgMfBY4G9sB9nk/i5J4OnAADh+GPn1cV9K8gQNtO4/qajj9dNfV\niPgrl6D4MfACdm3r1diV50TyrrLSup2CdsGihiQSx7qfFBQSdbl0PU0Hfg4cBEZgW4aL5F1Yup3S\nNKAtcZFLUHwO2wb8XuA+oAKY5mdREk9hCwoNaEtc5BIUdwMXARd6tyRwj481SUylu57CQkEhcZFL\nUHyAtSjStnvPieRV2FoUI0dazYcOua5ExF9NDWb/k/ezDNtJ9nHv8TXecyJ5Fbag6NQJhg+HDRvg\nnHNcVyPin6aC4kogPUt8L9btBPAOcJKfRUn81NWFr+sJYPx4ePVVBYVEW1NB8b8KVYTIvn3QoUP4\nLi9aUgKvvAL/+q+uKxHxT1NjFKVAvyZePx2Yl9dqJLa2boUzznBdRculg0IkyppqUZRhFyfqBKwB\n3saucNcfu0TqJ8Bdfhco8bBtGwxr7KonATZ6tF269f33oVcv19WI+KOpFsVMbFrsMmxvpxpsVfaL\n3msXY4PcIm22dWs4g6JDBzj7bFi1ynUlIv5pqkVRDAwArsXWTmRurKCt0CSvtm6FaSFdxpnufrr8\ncteViPijqRbFg8AzwHDsQkVlWbdcTAU2A1uBOY0cc6/3+jpsp9q0nsAiYBOwESjJ8XdKCIW1RQEa\np5Doy2X7tQeB1szpaA9swbYo341tKDgL++JPmwbc4v08F/g/HAuER4HngF9hLZ+TgQNZv6OuTvs8\nh15dnfXvV1RA796uq2m56moYNQrefRfa5bKEVdxIJLQvvCdhO2/mvP1mLv+sWzvxbyK2orsSG9tY\nAFyVdcwMLBAAVmKtiH5AD2xr8195rx3lxJCQiEh/wYYxJAD694cePaxVJBJFfv79MxCoyni8y3uu\nuWMGAUOwhX3/F5tx9TDQ1bdKxamwTo3NpO4nibJcrkfRWrm28bKbP3VYXROwbqnV2Dbn3wV+mP3m\n0tLS+vvJZJJkMtnySsWpMI9PpKWD4stfdl2JyIlSqRSpVKrV7/czKHYDgzMeD8ZaDE0dM8h7LuEd\nu9p7fhEWFCfIDAoJp7CuochUUgK//rXrKkQalv1H9Lx5LVsr7WfXUxkwDCjCFu3NBBZnHbMYuMG7\nXwLsB/YA1ViX1Jnea5cCr/tYqzgUhRbFuHHwxhtw8KDrSkTyz88WxVGs62g5NgPqEWzG003e6/Ox\nBXvTsEHvg8BXMt5/K/AYFjIVWa9JhLzxRviDonNnOOss2yDws591XY1IfoXg6sRN0vTYkKupgVNO\ngT177GeY3X67XT/7O99xXYk0SNNj6/kxPVbENzt3Qp8+4Q8JsHGKl192XYVI/ikoxKmNG22xWhRM\nmgQvvaQ/WiV6FBTi1MaNtgNrFAwaBF272piLSJQoKMSpKLUoAC64AF580XUVIvmloBCnohYUkyfD\nCy+4rkIkvxQU4kxdHWzaBCNHuq4kf9SikChSUIgzVVXQvXv4rpPdlJEj7Wp3b7/tuhKR/FFQiDOv\nvx6tbiewXXAnTVKrQqJFQSHORDEoQN1PEj0KCnGmvNz2SIoaBYVEjYJCnIlqUBQXw5Yt8Pe/u65E\nJD8UFOLEoUOwfXt0Fttl6tzZwkLbeUhUKCjEifXrYfhw6NTJdSX+UPeTRImCQpyIardTmoJCokRB\nIU6sXQvjx7uuwj/nnw+rV8Phw64rEWk7BYU4EfWg6NEDhg61/50iYaegkIKrqbExirFjXVfiL3U/\nSVQoKKTgNm+2K8F17+66En9pg0CJCgWFFNwrr9jV4KIuvZWHLmQkYaegkIKLS1AMGmSXeN2yxXUl\nIm2joJCCe+UVOO8811UUhrqfJAoUFFJQH3wAO3bAWWe5rqQwNKAtUaCgkIJatcqmxXbs6LqSwlBQ\nSBT4HRRTgc3AVmBOI8fc672+DsieWd8eWAss8atAKaw4dTsBjBgB+/fDW2+5rkSk9fwMivbAfVhY\njAJmAdkXvZwGnAEMA24EHsh6fTawEdC8kYiIy0B2Wrt2alVI+PkZFBOBbUAlcARYAFyVdcwM4FHv\n/kqgJ9DPezwIC5JfAgkf65QCqamBl16y7S3iREEhYednUAwEqjIe7/Key/WYe4BvA7V+FSiF9dpr\n0L+/3eLkggs080nCrYOPn51rd1F2ayEBTAf2YuMTyabeXFpaWn8/mUySTDZ5uDi0YgXE8f+e4mLY\nuhUOHLA9oEQKLZVKkUqlWv1+P7t0SoBSbIwCYC7WOvhpxjEPAimsWwps4DsJfAP4EnAUOAnoDvwR\nuCHrd9TVadlraMyYAddfDzNnuq6k8JJJmDsXLr/cdSUxlkhombwnkUhAC77//ex6KsMGqYuATsBM\nYHHWMYs59uVfAuwHqoHvAYOBIcB1wLOcGBISIjU11v0SxxYFqPtJws3PoDgK3AIsx2YuLQQ2ATd5\nN4ClwHZs0Hs+8PVGPkt/BoRcebltBNivX/PHRtHkyRrQlvAK+2widT2FxF132TWyf/EL15W48cEH\nMGAAvPdedC//GnjqeqoXpK4nkXp//StMmeK6Cne6d4dhw+DVV11XItJyCgrx3cGD8PLLcMklritx\nS+spJKwUFOK7556zKaJRv1BRczROIWGloBDfPfUUTJ3a/HFRl76QUa2WkErIKCjEd8uWKSgABg60\nBXebN7uuRKRlFBTiq23b4MMPYexY15UEg7qfJIwUFOKrJUvgc5+zmYmihXcSTgoK8dXixXBV9p7B\nMaaZTxJGYf87TwvuAuy992DIEKiuhi5dXFcTDHV1cNpptlJ90CDX1cSMFtzV04I7CYylS+GiixQS\nmRIJa1W89JLrSkRyp6AQ3zz5pLqdGqJxCgkbBYX44tAhePppmD7ddSXBM3kyPP+86ypEcqegEF8s\nWwZnn2398XK84mJ480145x3XlYjkRkEhvnjiCbj2WtdVBFOHDtaqWLHCdSUiuVFQSN4dOmQtin/8\nR9eVBNfFF8Ozz7quQiQ3CgrJO3U7NU9BIWGioJC8e/xxdTs1Z8wYeP99qKpyXYlI8xQUklcffWS7\nxarbqWnt2tkaE7UqJAwUFJJXy5bBOeeo2ykXU6bA8uWuqxBpnoJC8uqJJ+Caa1xXEQ5XXGGXiK2p\ncV2JSNMUFJI36nZqmUGDYMAAWLXKdSUiTVNQSN6o26nlpk2zPbFEgkxBIXmj2U4tN22aBaxIkBUi\nKKYCm4GtwJxGjrnXe30dMN57bjCwAngd2AB8w98ypS3U7dQ6550HFRXw9tuuKxFpnN9B0R64DwuL\nUcAsYGTWMdOAM4BhwI3AA97zR4DbgdFACXBzA++VgFi6FM49F/r0cV1JuHTsaFcA/POfXVci0ji/\ng2IisA2oxL74FwDZG0/PAB717q8EegL9gGqg3Hv+Q2ATMMDfcqW1Hn9cs51a65pr7PyJBJXfQTEQ\nyFx7ust7rrljsq/9VYR1Sa3Mc32SBx9+aNM81e3UOpdfble8q652XYlIwzr4/Pm5Xncw+5J8me/r\nBiwCZmMti+OUlpbW308mkySTyRYVKG33pz/Zbqjqdmqdk06y63b88Y9w882uq5EoSqVSpFKpVr/f\n72tmlwCl2BgFwFygFvhpxjEPAimsWwps4PtCYA/QEfgLsAz4eQOfr2tmB8Bll8FXv6oZT22xZAnc\neaeufOcrXTO7XtCumV2GDVIXAZ2AmcDirGMWAzd490uA/VhIJIBHgI00HBISAG+9BWVlcOWVrisJ\nt8svh61bYfNm15WInMjvoDgK3AIsx77wF2KD0jd5N4ClwHZs0Hs+8HXv+UnAF4GLgLXeLd0ykYD4\nwx9sbKJLF9eVhFunTvDP/wwPPeS6EpET+d315Dd1PTk2fjzcfbfthCpts2OHrWyvqlLw+kJdT/WC\n1vUkEbZhA7z7Llx4oetKomHIEJg4UVNlJXgUFNJqv/sdXH+9XVtB8uO22+AnP9GOshIs+k9cWqW2\nFh57DL70JdeVRMuUKXDqqbBgQfPHihSKgkJaZcUKWzcxerTrSqIlkYB58+BHP4KjR11XI2IUFNIq\nDz4IN97ouopouuQSu07F/PmuKxExmvUkLfbWW/CZz8DOnXDKKa6riabXX4dkEtavh/79XVcTEZr1\nVE+znsR3Dz8M112nkPDT6NG22v2221xXIqIWhbTQ4cM2jfOpp2DMGNfVRNtHH8GECXDHHTBrlutq\nIkAtinpqUYivfvMb63ZSSPiva1db+T57Nmzf7roaiTO1KCRnR4/C8OHw61/bbrFSGPfcAwsX2oaB\nHTu6ribE1KKopxaF+GbhQhg0SCFRaLNn29qKjB31RQpKLQrJyccf2wDrww/DxRe7riZ+9u61fbUe\nfRQuvdR1NSGlFkU9tSjEFz/7GYwbp5BwpW9fGx+64QabnixSSGpRSLMqK6G4GF59FYqKXFcTbz/6\nETzzjN06+H19yqhRi6KeWhSSV4cP29TM735XIREE//Zv0Lkz/PCHriuROFGLQpr0zW/aldeefFK7\nxAbF3r3WwnvoIbjiCtfVhIhaFPVa2qJQUEij7rzTBq9XrrRZNxIcL7wA11wDq1fD4MGuqwkJBUU9\ndT1Jm9XU2GrgRx6BVEohEUSTJ1tr79prrXtQxE9qUchxKivt2s01NXZNhNNPd12RNKa2Fq6+2naa\nfeAB+4NZmqAWRT21KKRV3noLbr/d+r6nTIFnn1VIBF27dnaVwZdegv/6L9fVSJQpKGJu0yb4l3+x\nxXQ1Nba99dy50L6968okF927w5IldvnUhQtdVyNRpZnYMbV5M/zgB/D883DLLbBtG/Tu7boqaY2i\nItvNd8oU6NIFZsxwXZFEjd8tiqnAZmArMKeRY+71Xl8HjG/he6WFqqrsOgeTJ8PZZ8OOHRYYColw\nO+ss+Mtf7KqDjzziuhqJGj+Doj1wH/aFPwqYBYzMOmYacAYwDLgReKAF75UMqVSqydf37rUxiLFj\nbTuIN96AOXNsK+uoae5cRNU551gL8cc/hltvtetZxPVcNCTluoAQ8zMoJgLbgErgCLAAuCrrmBnA\no979lUBPoH+O75UMjX0hVFfD978PI0faGMTGjfAf/wG9ehW2vkKK85fjmWfa2or337e9ue6/P0Vt\nreuqgiHluoAQ83OMYiBQlfF4F3BuDscMBAbk8F5pwMcfw5YtNhNm6VL7ee21sGYNfOpTrquTQujV\ny2ZDPfWUTVQ480z7N3DppbYDbc+eTU+lPXwY9u+3sNm//8T7Bw/a7+jd226nnWat1L594eSTC/e/\nUwrHz6DIdcJym2Z/T5+e9Usb+K1Rf+6jj6CiwqZIHjwIQ4daN8QXvmBXSNO1reNp6lQbj7rySli0\nyBZRrlsHn3xiiyh79bILIR0+bLdDhywIjhyx13r2tFv2/S5dbKyrvBzefddue/fCnj0WQH37Qrdu\n0KmT7UuV/tnYTLrGQitfz6cVY+dCgqUEeCrj8VxOHJR+ELgu4/FmoF+O7wXrnqrTTTfddNOtRbdt\nBEQHoAIoAjoB5TQ8mL3Uu18CvNKC94qISARcAWzB0muu99xN3i3tPu/1dcCEZt4rIiIiIiKSH3Fe\nkPcrYA+wPuO5U4GngTeAv2JTjeNgMLACeB3YAHzDez6O5+MkbJp5ObAR+LH3fBzPBdh6rLXAEu9x\nXM8D2FKD17Dzscp7LvLnoz3WJVUEdCR+YxiTsVXsmUFxJ/Ad7/4c4CeFLsqR/sA47343rLtyJPE9\nH+kllB2wMb8LiO+5+CbwGLDYexzX8wCwAwuGTJE/H+dx/Kyo73q3OCni+KBIzxgD+/LcXOiCAuLP\nwKXofHQFVgOjiee5GAT8DbiIYy2KOJ6HtB1A9kY9OZ+PsO4e29hCvTjrh3VH4f3s18SxUVWEtbRW\nEt/z0Q5rYe/hWJdcHM/FPcC3gcx16XE8D2l1WHCWAV/znsv5fIR199g61wUEXHqudJx0A/4IzAb+\nnvVanM5HLdYV1wNYjv1FnSkO52I6sBfrj082ckwczkOmScDbwGnYuER266HJ8xHWFsVubBAzbTDW\nqoizPVjzEeB07D+UuOiIhcRvsa4niPf5ADgA/D9sQXLczsX52D5yO4A/ABdj/zbidh4yve39fAf4\nb2w/vZzPR1iDogzbcbYIW5A3k2MDVnG1GPiyd//LHPvCjLoE8Ag2y+fnGc/H8Xz04djMlS7AFOyv\n6ridi+9hfzwOwXZ+eBb4EvE7D2ldgfRGPicDl2Hjm7E4H3FekPcH4C3gMDZW8xVsRsPfiPBUt0Zc\ngHW3lGNfimuxqdNxPB9jgDXYuXgN66OHeJ6LtAs59kdkXM/DEOzfRDk2hTz9fRnX8yEiIiIiIiIi\nIiIiIiIiIiIiIiIiIiIirVfE8ZszikROWFdmi4hIgSgoRNquA/A7bBuRJ4DPYntPAVwFfOQdcxJ2\nLXiRUFFQiLTdcOB+YBTwAXAOxy6mNBnrmpoInItdTEgkVMK6zbhIkFQBL3v3f4dtdb4NGIGFxt1Y\nK6M98ILhkIJwAAAAfElEQVSLAkXaQi0KkbbL3Mc/gW1S+DwwDTgCPIO1LC5AQSEhpKAQabt/AEq8\n+1/AwuBF4Dbgf4B3sctQnoldcU4kVBQUIm1Th213fzM2mN0DeABYBfTFWhYA69A0WhERERERERER\nEREREREREREREREREREREREREZHw+v8em8Gl3l7PLgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1007e7110>"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "density = gaussian_kde(jc.sim_bb)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(xs, density(xs))\n",
      "axvline(x=jc.bb, color='r')\n",
      "xlabel('bb')\n",
      "ylabel('f(bb)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 28,
       "text": [
        "<matplotlib.text.Text at 0x108581d50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG0RJREFUeJzt3XuQXOWdn/FnJI0QQleEkEASErqALiBuayGDFAYbZIHB\ngFMxsdfOOsYOqZiNd+NKiLNlI5e3yk7WTlgvXoPXbNbluExsByglgLGJGS8mGCGhiy0kodF1dL8B\nEhdJI+j88Z72tEajnu6ePrfu51PVNd2nz+nz4wD97fd9z3kPSJIkSZIkSZIkSZIkSZIkSarQYmA9\nsBG4t5f3/xhYDawBngfmlry3NVq+ElgWa5WSpFQMBDqAKUArsAqY1WOd9wMjo+eLgd+WvLcFODve\nEiVJfRkQ42fPIwTFVqALeAS4rcc6LwBvRM9fBCb2eL8lxvokSRWIMygmAJ0lr3dEy07nLuDJktcF\n4BlgOfC5ulcnSarIoBg/u1DFutcDnwGuLVl2LbAbGAv8kjDW8VzdqpMkVSTOoNgJTCp5PYnQquhp\nLvB3hDGK10qW747+7gceI3RlnRQU06CwqV7VSlLz2ARMr3TlOLuelgMzCIPZg4E7gaU91rkAeBT4\nJGE8o2goMDx6fhawCPhdzx1sAgqFgo9Cgfvuuy/1GrLy8Fh4LDwW5R/AtGq+zONsUZwA7gGeJpwB\n9TCwDrg7ev8h4CvAaOC70bIuQsthPCFAijX+CPhFjLVKkk4jzqAAeCp6lHqo5Plno0dPm4HL4ypK\nklS5OLuelKC2tra0S8gMj0U3j0U3j0Xt8n6dQiHqb5MkVailpQWq+P63RdEA3n0Xbrop/JWkejMo\nGsDu3fDzn8PWrWlXIqkRGRQNoDO6/n3DhnTrkNSYDIoGYFBIipNB0QA6O2HoUINCUjwMigawYwcs\nXAivvpp2JZIakUHRADo74YYbbFFIiodB0QA6O+H974fXXoMjR9KuRlKjMSgaQGcnTJ4MM2bY/SSp\n/gyKnDt+HA4cgPPOg4svtvtJUv0ZFDm3axeMHw8DBxoUkuJhUORcZydMjO40fvHFsH59uvVIajwG\nRc51dsKk6D6Cl10Gq1eH5+++C9dcA/ffD11d6dUnKf8MipwrDYqZM8PrN98MLYvt22HpUrjrrnRr\nlJRvBkXOlQZFayvMng1r1sCLL8J118HXv253lKT+MShybteu7jEKgCuugJUrQ1BcfTWMGgWvv55e\nfZLyz6DIuddfD2FQdMUV8PLLsGxZCIrRo8OFeJJUK4Mi5w4fhhEjul9feSU8/3y48O6yy7pbFN4I\nUFKtDIqc6xkUl14KGzfCnDkwZAgMHhweb72VXo2S8s2gyLkjR04OiqFDw9lPV1/dvczuJ0n9YVDk\nXM8WBcDNN8OiRd2vHdCW1B+D0i5AtTtxAo4dC62IUn/1Vye/tkUhqT9sUeTYkSMwfDi0tJRfzxaF\npP4wKHKst26n3tiikNQfBkWOVRMUtigk1cqgyLFKg2LUKFsUkmpnUOTY4cNhjKIvdj1J6g+DIseq\naVHY9SSpVgZFjvW82O50bFFI6g+DIsdsUUhKgkGRY54eKykJBkWOeXqspCQYFDnm6bGSkmBQ5Fil\nQTFsGBw9Cl1d8dckqfHEHRSLgfXARuDeXt7/Y2A1sAZ4HphbxbZNr9KgaGlxQFtS7eIMioHAA4Qv\n/NnAx4FZPdbZDPwTQkB8DfheFds2vUqDAhzQllS7OINiHtABbAW6gEeA23qs8wLwRvT8RWBiFds2\nvWqDwhaFpFrEGRQTgM6S1zuiZadzF/Bkjds2pUqn8AAHtCXVLs4bFxWqWPd64DPAtdVuu2TJkj88\nb2tro62trYrd5lulV2aDLQqpmbW3t9Pe3l7z9nEGxU5gUsnrSYSWQU9zgb8jjEcUf/NWuu1JQdFM\nCoXuGxdVwhaF1Lx6/oj+6le/WtX2cXY9LQdmAFOAwcCdwNIe61wAPAp8kjAmUc22Te3tt+GMM2BQ\nhVE/ciS88Ubf60lST3G2KE4A9wBPE85iehhYB9wdvf8Q8BVgNPDdaFkXYSD7dNsqUs1ANhgUkmoX\nZ1AAPBU9Sj1U8vyz0aPSbRWpJSh27YqvHkmNyyuzc8oWhaSkGBQ5ZVBISopBkVMGhaSkGBQ59cYb\nBoWkZBgUObVvH5x7buXrGxSSamVQ5NSePTBuXOXrGxSSamVQ5NTevTB+fOXrjxgRxjUK1UysIkkY\nFLm1d291LYpBg2DIEHjzzfhqktSYDIqc2rOnuhYFhO6nw4fjqUdS4zIocqraFgU4TiGpNgZFDnV1\nhSnDx4ypbjuDQlItDIoc2r8fzjkHBg6sbjuDQlItDIocqqXbCQwKSbUxKHKo2msoigwKSbUwKHKo\n2msoigwKSbUwKHLIFoWkJBkUOWSLQlKSDIoccjBbUpIMihyy60lSkgyKHLLrSVKSDIocskUhKUkG\nRc4cPx4m9qt2+g4wKCTVxqDImd//HqZPhwE1/JsbMcKgkFQ9gyJnnn8eFi6sbdtii8KbF0mqhkGR\nM7/5DSxYUNu2gwdDayu8/XZ9a5LU2AyKHCkUQlBce23tnzFmDBw8WL+aJDU+gyJHtm2D996DqVNr\n/4zJk2Hr1rqVJKkJGBQ5Uux2ammp/TMuvBC2bKlfTZIan0GRI/0ZnygyKCRVy6DIkTVr4Mor+/cZ\nBoWkahkUObJjB0ya1L/PMCgkVcugyIl33w1Td5x/fv8+x6CQVC2DIif27YPRo8O1EP0xcWL4rGPH\n6lOXpMZnUOTEjh3hS76/Bg2CCRNg+/b+f5ak5mBQ5ES9ggJC99PmzfX5LEmNL+6gWAysBzYC9/by\n/kzgBeAo8MUe720F1gArgWXxlZgP9Q4KxykkVWpQjJ89EHgAuAHYCbwELAXWlaxzEPhT4PZeti8A\nbcChGGvMjZ07Q5dRPRgUkqoRZ4tiHtBBaBl0AY8At/VYZz+wPHq/N/24Brmx2KKQlJY4g2IC0Fny\neke0rFIF4BlCkHyujnXlUj2DYto06Oioz2dJanxxdj31964H1wK7gbHALwljHc/1t6i8qmfX08yZ\nsGFDmGCwlhsgSWoucQbFTqD0OuJJhFZFpXZHf/cDjxG6sk4JiiVLlvzheVtbG21tbVWWmX2FQmhR\n1CsoRo6Es88Os8j2ZyZaSfnQ3t5Oe3t7zdvHOQYwCNgAfBDYRThz6eOcPJhdtAQ4Anwrej2UMBh+\nBDgL+AXw1ehvqUKhCW7XdvBguP3pa6/V7zMXL4bPfx5uvbV+nykpH1rCFNQVf//H2aI4AdwDPE34\n0n+YEBJ3R+8/BIwnnA01AngP+AIwGzgXeLSkxh9xakg0jXp2OxXNmQNr1xoUkvoWZ1AAPBU9Sj1U\n8nwPJ3dPFb0JXB5XUXlTz4Hsojlz4Nln6/uZkhqTQ5k5UM/xiaJii0KS+mJQ5MD+/TBuXH0/c/Zs\nWL8+zEorSeUYFDlw4ACMGVPfzxw+HMaO9cI7SX0zKHLg4EE455z6f+6cOfDKK/X/XEmNxaDIgQMH\n4gmKWbPChXeSVI5BkQMHD9a/6wnCuMfevfX/XEmNxaDIgbhaFGPHhoFySSrHoMgBg0JSmqoJiiHA\nGXEVot51dcHbb4f5merNoJBUiXJBMQD4KPBTwgR/W4Bt0fOfAXfg/SJid/BgmMCvJYYjbVBIqkS5\noGgHrgK+CUwFziPMzTQ1WvY+4Ncx19f04up2AoNCUmXKzfV0I3Csl+XHgN9GD7uiYhbXNRQAw4aF\nK7PffhuGDo1nH5Lyr1xQlIbEVcACwgyvzwMv97KOYhDHVdlFLS3drYrJk+PZh6T8q2Qw+yvAPwBn\nE+4299+BL8dYk0rE2fUEcO65dj9JKq+SacY/CcwFjkavvw6sBr4WV1HqFmfXEzhOIalvlbQodgJn\nlrweQnW3NFU/xNn1BAaFpL6Va1H8TfT3DWAt3XeYu5FwW1Ml4MABuOyy+D7foJDUl3JBsQIoAMuB\nx6PnEE6bbfwbVWeEXU+S0lYuKP6h5PkZwExCQKwHjsdYk0ok0fXU0RHf50vKv0oGsz8MPAhsjl5P\nBe4GnoyrKHWL+6wnWxSS+lJJUPxX4Hqg+LtzGiEkDIoE2PUkKW2VnPV0mO6QgNCyOBxPOSrV1QVv\nvRXPhIBFBoWkvpRrUfzT6O9yQuvhJ9HrfxYtU8wOHIhvQsAig0JSX8oFxa10n920D7guer6fcC2F\nYrZtW/xTa4wcCe+8A8eOwRnO3CWpF+WC4tNJFaHebd4MU6fGu4+WljAGsn8/TJwY774k5VO5MYol\nwLgy758HfLWu1egkW7bEHxRg95Ok8sq1KJYDjwCDCbPF7ibcqGg8cCVh5thvxl1gM9u8GebPj38/\nBoWkcsq1KO4knBb7FPAb4F2gK3p+J/ABPEU2Vkl0PYFBIam8ci2Kq4DzgY8BbZx821On8EjAli1w\n4YXx78egkFROuaB4EPi/hCuxV/R4rxAtV0yOH4fdu2HSpPj3ZVBIKqdc19O3gVmEGxVd2ONhSMRs\n+3aYMAFaW+Pfl0EhqZxKrsz+17FXoVMk1e0EBoWk8ioJCqUgqYFsMCgklWdQZFSSLQrvmy2pHIMi\no2xRSMoKgyKj1q6Fiy5KZl+jR8ORI2G2WknqKe6gWEy4I95G4N5e3p8JvAAcBb5Y5bYNa/t22Lcv\n3ntllxowIMxSe/BgMvuTlC9xBsVA4AHCF/5s4OOE021LHQT+lFOnAqlk24b11FOweDEMHJjcPu1+\nknQ6cQbFPMINj7YSpv54BLitxzr7CXNK9ez0qGTbhvXkk3Dzzcnu06CQdDpxBsUEoLPk9Y5oWdzb\n5tqxY9DeDosWJbtfg0LS6cQZFP2ZD6pp55J67jm45BIYMybZ/RoUkk6n3FxP/bUTKJ2paBKhZVDX\nbZcsWfKH521tbbS1tVVTY+Y89xxcf33y+zUopMbV3t5Oe3t7zdvHeDdmBgEbgA8Cu4BlhEHpdb2s\nuwQ4Anyrym0LhUJjNT4+8Qm46Sb41KeS3e93vhNOyf3bv012v5KS19LSAlV8/8fZojgB3AM8TTiL\n6WHCF/3d0fsPEW6C9BIwAngP+ALhLKc3T7Ntw+vogOnTk9+vLQpJpxNniyIJDdWiKBTC9QwbN4b7\nWCfp2WdhyRL49a+T3a+k5FXbovDK7Aw5dCiERdID2QDnnQe7diW/X0nZZ1BkSEcHzJgBLSm086ZO\nhc7OcHquJJUyKDJk48Z0xicABg+GyZNh06Z09i8puwyKDElrILto5kxY1xSnDEiqhkGRIcWup7TM\nmgXr16e3f0nZZFBkSJpdT2CLQlLvDIoMyULXky0KST0ZFBlx6BCcOBEufEvLzJmwYUM4RVeSigyK\njNi2LZx1lMapsUWjRsGwYbCj0hm5JDUFgyIjdu+G889Puwq7nySdyqDIiN27w9XRaZs1ywFtSScz\nKDIiK0Exdy6sWpV2FZKyxKDIiKwExbx5sGxZ2lVIyhKDIiOyEhSXXgpbtsCRI2lXIikrDIqMyEpQ\ntLbCZZfBihVpVyIpKwyKjMjKWU9g95OkkxkUGVAowJ492WhRgEEh6WQGRQYcOgRnnglDhqRdSWBQ\nSCplUGRAVsYniqZNg7feCnVJkkGRAVkLipYWuOQSr9CWFBgUGZC1oACYMgW2bk27CklZYFBkgEEh\nKcsMigwwKCRlmUGRAQaFpCwzKDLAoJCUZQZFBuzYARMmpF3FySZODBcBdnWlXYmktBkUKevqgp07\nw93tsqS1FcaP9253kgyK1G3fHuZ4Gjw47UpOZfeTJDAoUrdpE0ydmnYVvTMoJIFBkbpNm8KUGVlk\nUEgCgyJ1BoWkrDMoUmbXk6SsMyhStnlzdlsUkycbFJKgJe0C+qlQKBTSrqFmhQIMHx5Ojx05Mu1q\nTnX8OAwbBu+8AwMHpl2NpHppaWmBKr7/bVGkaN++cLOiLIYEhFN2x471vhRSszMoUpTl8YmiCy4I\n13pIal5xB8ViYD2wEbj3NOt8O3p/NXBFyfKtwBpgJdCQN+bM8vhEkUEhaVCMnz0QeAC4AdgJvAQs\nBdaVrHMzMB2YAVwNfBeYH71XANqAQzHWmKqNGw0KSdkXZ4tiHtBBaBl0AY8At/VY5yPAD6LnLwKj\ngHEl7+d9sL2s55+Hq69Ou4ryJk0yKKRmF2dQTAA6S17viJZVuk4BeAZYDnwuphpTc/w4vPgiLFiQ\ndiXl2aKQFGfXU6XnrZ6u1bAA2AWMBX5JGOt4rudKS5Ys+cPztrY22traqqkxNStWhG6n0aPTrqQ8\ng0LKv/b2dtrb22vePs6unfnAEsKANsCXgPeA/1yyzoNAO6FbCkIYXAfs7fFZ9wFvAt/qsTy311F8\n4xvhtNO//uu0KynvwAG46CI41LAjRVLzydJ1FMsJg9RTgMHAnYTB7FJLgX8RPZ8PvE4IiaHA8Gj5\nWcAi4Hcx1pq4f/xHuO66tKvo25gxcPQoHDmSdiWS0hJnUJwA7gGeBl4B/ifhjKe7owfAk8BmwqD3\nQ8C/iZaPJ3QzrSIMcv8f4Bcx1pqoEyfCQPbChWlX0reWltD91NnZ97qSGlOcYxQAT0WPUg/1eH1P\nL9ttBi6PpaIMWLMm3Pp07Ni0K6lMcZxi9uy0K5GUBq/MTsGyZTB/ft/rZYUD2lJzMyhS8NJL8L73\npV1F5S64ALZtS7sKSWkxKFKQt6CYNQvWrk27CklpMSgS9tZbYTLAuXPTrqRyl18Oq1alXYWktBgU\nCXv5ZZgzJ0zhnRfTpoXrKLyWQmpOBkXC8tbtBDBgQGgBrV6ddiWS0mBQJCyPQQFwxRV2P0nNyqBI\n2IoV+QyKyy+HlSvTrkJSGgyKBB09Gq5wvuiitCupni0KqXkZFAnauBEuvBBaW9OupHqzZ4f6jx5N\nuxJJSTMoErRuHcycmXYVtRkyJLSEHNCWmo9BkaD16/MbFACLFsETT6RdhaSkGRQJWr8+XOWcVx/9\nKDz6aNpVSEqaQZGgPHc9Qbi/96FD8OqraVciKUkGRULeey98weY5KAYMgDvugMceS7sSSUkyKBKy\nfXu4P/bw4X2vm2V33GH3k9RsDIqE5H18oui666CjwzveSc3EoEhI3s94KmpthVtvhccfT7sSSUkx\nKBKydm1jtCjAs5+kZmNQJGTlyjANRiO48cYwXfr+/WlXIikJBkUCurrglVfydbOics48Ez70IVi6\nNO1KJCXBoEjAunUweTKcdVbaldSP3U9S8zAoEtBI3U5FN98Mzz0Hhw+nXYmkuBkUCVi1KtzPoZGM\nGAELF8KTT6ZdiaS4GRQJaMQWBdj9JDWLlrQL6KdCoVBIu4ayCoVwRXZHB5xzTtrV1Nf+/TB9OuzZ\nEwa4JeVDS0sLVPH9b4siZlu2hGk7Gi0kAMaODbd1/elP065EUpwMipg98wxce23aVcTnvvvgy1+G\nd95JuxJJcTEoYvbII/Cxj6VdRXwWLoQ/+iO4//60K5EUF8coYrRnT5i2Y9euxu7D7+iA+fPh6afh\nqqvSrkZSXxyjyJCf/QxuuaWxQwLCgPb3vx/+WTdsSLsaSfU2KO0CGtmPfwxf+lLaVSTj9tvh9ddh\n8WJ44QUYPz7tiiTVi0ERk6eegr17YdGitCtJzqc/He5Tccst0N4Ow4alXZGkenCMIgbvvAOXXAIP\nPAA33ZR2NckqFOCzn4V9+8ItUwf5U0TKHMcoMuAv/zJM2dFsIQHQ0gIPPgjHjsE994R7hUvKt7iD\nYjGwHtgI3Huadb4dvb8aKJ3oopJtM+eJJ+AHPwitiWbV2hoG8tetg+uvh02b0q5IUn/EGRQDgQcI\nX/izgY8DPe/xdjMwHZgB/Cvgu1Vsmznr1sFnPgM/+Qmcd16y+25vb092h30YMQJ+9Sv4yEdg3rww\nfvHKK8nsO2vHIk0ei24ei9rFGRTzgA5gK9AFPALc1mOdjwA/iJ6/CIwCxle4baasXAkf+AB885tw\nzTXJ7z+L/xMMHAhf/GK4zmLGjNC6uP12WL063v1m8VikxWPRzWNRuziHGicAnSWvdwBXV7DOBOD8\nCraNTVcX7NwJGzeGL7nXXgt97WPHwoQJ4SZEkyeH6yO2bQvXEHz/+/C974UZVXWy0aPhL/4C/vzP\nwzH60IdgwYLQ+vrgB+GMM9KuUFI5cQZFpacj9evMqw9/OJxpU/qAU5dVsvzYsRAQBw7AuHHhV/D0\n6d0T+i1fDo8/HsJh27YQKOPHh+6VFStCeOj0hg6FP/szuOsu+OEP4WtfC8F6wQUwalQI3iFDwhhH\nf7z6avh3lVctdTwXccOG8N+mPBZZNR/4ecnrL3HqoPSDwD8veb0eGFfhthC6pwo+fPjw4aOqRwcZ\nMQjYBEwBBgOr6H0wu3iPtPnAb6vYVpLUAG4CNhDSqziZxd3Ro+iB6P3VwJV9bCtJkiRJ9ZHLC/Lq\n5O+BvcDvSpadDfwSeBX4BeFU40Y3CXgWWAv8Hvi30fJmPBZDCKeYrwJeAb4eLW/GY1E0EFgJ/O/o\ndbMei63AGsKxWBYta4pjMZDQJTUFaKX5xjAWEq5iLw2K/wL8h+j5vcA3ki4qBeOBy6PnwwhdlbNo\nzmMBMDT6O4gw3reA5j0WAP8O+BGwNHrdrMdiCyEYSjXFsXg/J58V9R+jRzOZwslBUTxjDMIX6Pqk\nC8qAx4Eb8FgMBV4C5tC8x2Ii8AxwPd0timY9FluAMT2WVXUs8jop4Oku1Gtm4wjdUUR/x5VZtxFN\nIbSyXqR5j8UAQut6L91dcs16LP4b8O+B0mkpm/VYFAihuRz4XLSsqmOR10mgC2kXkHHFc6WbxTDg\nfwFfAI70eK+ZjsV7hK64kcDThF/TpZrlWNwC7CP0ybedZp1mORYA1wK7gbGEcYmerYc+j0VeWxQ7\nCQOZRZMIrYpmtpfQhAQ4j/A/SjNoJYTEDwldT9C8x6LoDeAJ4Cqa81hcQ5hHbgvwY+ADhP8+mvFY\nQAgJgP3AY4S59Ko6FnkNiuWEGWenEC7Iu5PuAatmtRT4k+j5n9D9pdnIWoCHCWf53F+yvBmPxTl0\nn7lyJnAj4Rd1Mx6L/0T48XghYeaHXwGfojmPxVBgePT8LGARYWyzaY5FM1+Q92NgF3CcMFbzLwln\nNTxDg5/u1sMCQnfLKsKX4krCadPNeCwuBV4mHIs1hP55aM5jUeo6un9ENuOxuJDw38Qqwinkxe/K\nZjwWkiRJkiRJkiRJkiRJkiRJkqT6mMLJEzQWbeXUWTul3MnrldlSHhQIV49LuWZQSPUxCPgfhOlE\nfkqYRgPCnP9rCLPaTkunNKl/DAqpPi4GvgPMBg4Dn4+Wvw7MJdwb/v7eN5UkNbopwLaS19cTJlnb\nEr0HYZbbA4lWJdWJLQqpPkrn82/h5Bvm9LaOlBsGhVQfFwDzo+efAH4TPb+z5O//S7ooSVI2TAbW\nEW6OUzqYvYVw0/rVhMHsqWkVKEmSJEmSJEmSJEmSJEmSJEmSJEmSpH76/6VmcNEOpdiHAAAAAElF\nTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10856ca50>"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc.bb"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "0.0"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w.s0/2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 30,
       "text": [
        "40.0"
       ]
      }
     ],
     "prompt_number": 30
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Random pattern"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import bernoulli"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = bernoulli.rvs(0.5, size=25)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 33,
       "text": [
        "array([1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0,\n",
        "       0, 1])"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc_random = ps.Join_Counts(y,w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y.shape = (5,5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "matshow(y,cmap=cm.gray_r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 36,
       "text": [
        "<matplotlib.image.AxesImage at 0x10a0b66d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAPYAAAD7CAYAAABZjGkWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAB5VJREFUeJzt3cuLnfUdx/H3ZCYhlixcuPCSQDYtuBAUpAuJ5mhbsOJt\nJy7di0FRirvZuYx/gI1oW9KFhdJU2i5KjlrE4CXBS1RUFGyLEaoIEhQl08XvJNrE9PkdznM5z+e8\nXxDmZM5h5pvM857nMsPvAUmSJEmSJEmSJEktuxV4B3gP+NXAs9Q4BJwC3hh6kDnsAY4CbwFvAg8M\nO06jncAx4ARwEnhs2HHmsg4cB44MPciQ1oH3gb3AdsoX8uohB6pwI3Ad4wr7cuDa2eNdwLss///z\nj2ZvN4CXgH0DzjKPh4DfAX/q6xNu6+sTzeGnlLA/Ar4Bfg/cNeRAFV4APh96iDl9QvmmCfAl8DZw\n5XDjVDk9e7uDsgP4bMBZau0GbgOeANb6+qTLGPZVwMff+/s/Z+9Td/ZSjjiODTxHk22Ub0anKKcR\nJ4cdp8pB4BHgTJ+fdBnD3hp6gBWzC3gGOEDZcy+zM5TTh93ATcBk0Gma3Q58Sjm/7m1vDcsZ9r8o\nF3bO2kPZa6t924E/AL8F/jjwLPP4AngWuH7oQRrcANwJfAgcBm4Bnh50ogFtAB9QDg93MI6LZ1Dm\nHdPFszXKRnZw6EEqXQZcOnt8CfA88LPhxpnbflb8qjjALylXad8HHh14lhqHgX8DX1OuD9w37DhV\n9lEObU9QDhWPU37MuKyuAV6jzPs65bx1TPbT41VxSZIkSZIkSZJWRxu/DeNviknDuqDjZfzNM0kL\nMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVAhi0Fqgl7\nbPeqllZe09JI65Q7cvycck+tl4F7KbdcPculkaRhzb000hjvVS2tvKawvVe1NEJNYXuYLY1QU9je\nq1oaoaawXwF+zHf3qr4HbwUqLb2Nhue/Be4H/ka5Qv5r/veKuKQl5J1ApPHzTiDSKjBsKZBhS4EM\nWwpk2FIgw5YCGbYUyLClQIYtBTJsKZBhS4EMWwpk2FIgw5YCGbYUyLClQIYtBWpaGqnK1paLqHRt\nba2NxW7UZGzb8sW2C/fYUiDDlgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVAhi0F\nMmwpkGFLgQxbCmTYUqCasA8Bp4A3Op5FUktqwn4SuLXrQSS1pybsF4DPux5EUns8x5YCGbYUqJXl\nhzc3N889nkwmTCaTNj6spPNMp1Om02nj62oXq94LHAGu+YHntsa2FvMYua54P8a2Lc+2iws2jppD\n8cPAi8BPgI+B+1qdTFLr2tgNuMfugXvsfoxtW15kjy1pZAxbCmTYUiDDlgIZthTIsKVAhi0FMmwp\nkGFLgQxbCmTYUiDDlgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVAhi0FMmwpkGFL\ngQxbCmTYUiDDlgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUiDDlgLVhL0HOAq8BbwJPNDpRJIWtlbx\nmstnf04Au4BXgbuBt2fPb21tbXUznc5ZW6v5UmlRY9uWZ9vFBRtHzR77E0rUAF9Sgr6ytckktW7e\nc+y9wHXAsfZHkdSWjTleuwt4BjhA2XOfs7m5ee7xZDJhMpm0MJqk802nU6bTaePrak/ctgN/Bv4C\nPH7ec55j98Bz7H6MbVu+2Dl2zdayBjwF/Ad48AeeN+weGHY/xrYtLxL2PuB54HXg7L/6UeCvs8eG\n3QPD7sfYtuVFwm5i2D0w7H6MbVte5MddkkbGsKVAhi0FMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVA\nhi0FMmwpkGFLgQxbCmTYUiDDlgIZthSolRVUWvgYajC2lT3AVV965Aoq0iowbCmQYUuBDFsKZNhS\nIMOWAhm2FMiwpUCGLQUybCmQYUuBDFsKZNhSIMOWAhm2FMiwpUCGLQWqCXsncAw4AZwEHut0IkkL\n26h4zVfAzcDp2ev/AeybvZW0hGoPxU/P3u4A1oHPuhlHUhtqw95GORQ/BRylHJJLWlK1YZ8BrgV2\nAzcBk64GkrS4ea+KfwE8C1zfwSySWlIT9mXApbPHlwC/AI53NpGkhdVcFb8CeIryTWAb8Bvg710O\nJWkx3glkJLwTiP4P7wQirQLDlgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUiDDlgIZthTIsKVAhi0F\nMmwpkGFLgQxbCmTYUiDDlgKt5NJIY1xmSP0Y6XJOLo0krQLDlgIZthTIsKVAhi0FMmwpkGFLgQxb\nCmTYUiDDlgIZthTIsKVAhi0FMmwpkGFLgQxbCmTYUqDasNeB48CRDmeR1JLasA8AJxnhMkjSKqoJ\nezdwG/AE7ayRJqljNWEfBB4BznQ8i6SWNIV9O/Ap5fzavbU0Ek1h3wDcCXwIHAZuAZ7ueihJi5ln\nL7wfeBi447z3j+6CmuuK62JWdV1xi5BGwDuBSN+zqntsSSNg2FIgw5YCGbYUyLClQIYtBTJsKZBh\nS4EMWwpk2FIgw5YCGbYUyLClQIYtBTLsFk2n06FHmJszZzLsFo1xg3PmTIYtBTJsKVAb68BMKQsd\nSurfc8Bk6CEkSZIkSZIkSRrOfwEhkw5fME28MwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10878cd50>"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc_random.p_sim_bb"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 37,
       "text": [
        "0.38200000000000001"
       ]
      }
     ],
     "prompt_number": 37
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "density = gaussian_kde(jc_random.sim_bb)\n",
      "xs = linspace(0,30,200)\n",
      "plot(xs, density(xs))\n",
      "axvline(x=jc_random.bb, color='r')\n",
      "xlabel('bb')\n",
      "ylabel('f(bb)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 38,
       "text": [
        "<matplotlib.text.Text at 0x10a0f6b10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHoZJREFUeJzt3XmUVPWd9/F3080mIIsiQgs0zW6QAApiRNO4BTWOS2Yg\n6LhiQsaoOY8zCeOMxiZmzjyjjxlNzCE8E4PEaOCZKAkTFxSlXQiLrdCEsDUNjQ2IAqIsEQS7nj++\nVXZRdFdXddet310+r3PqdNWte6u+xaX7U7/f797fBREREREREREREREREREREREREREREcnQJGAD\nUA3MaOT5G4AqYA2wFBiZ9FxtfPkqYKWnVYqIiBOFwGagBGgLrAaGp6xzHtA1fn8SsDzpua1AD29L\nFBGR5rTx8LXHYUFRCxwF5gFXp6yzDPgkfn8FcEbK8wUe1iciIhnwMiiKgbqkx9vjy5oyDXgh6XEM\nWAxUAt/KeXUiIpKRIg9fO5bFuhOB24Dzk5adD7wP9ARewcY63sxZdSIikhEvg2IH0DfpcV+sVZFq\nJPBf2BjFvqTl78d/7gYWYF1ZxwXFwIEDYzU1NbmqV0QkKmqAQa6LAAuhGmwwux2ND2b3w8Yxxqcs\nPwnoEr/fCTsi6rJG3iMWZg888IDrEjwV+s8X4v+fod93If98ZNfj42mL4hhwJ7AIOwLqCWA9MD3+\n/Gzgh0B3YFZ82VGs5XA68FxSjU8DL3tYq4iINMHLoAB4MX5LNjvp/u3xW6otwCivihIRkcx5edST\ntFJZWZnrEjwV+s/nugAPhX7fhfzzZSvo5ynEu9tEfKigAPT/U3yooKAAsvj7rxaFiIikpaAQEZG0\nFBQiIpKWgkJERNJSUIiISFoKChERSUtBISIiaSkoREQkLQWFiIikpaAQEZG0FBQiIpKWgkIkA/v3\nw403wq5drisRyT8FhUgzDh2CK6+EBQtg8WLX1Yjkn4JCIq+uLv3z//Zv0KcPzJwJf/pTfmoS8RMF\nhUTa3LnQvz9UVze9ztKlcPvtMGGCgkKiSUEhkbVwIdx7L1xwAVRUNL7OsWPw7rswdiyMHg2bN9t4\nhUiUKCgksubMgYcfhptuajoo1q2D4mLo1g3atYMxY2DFiryWKeKcgkIiq7YWhg6FsjILisYuRrdy\nJYwb1/D4K19R95NEj4JCImvrVhgwAEpLoU0b61ZKlRoU559vYxYiUaKgkEjatw/q66FHD7u09cSJ\njXc/rVwJ557b8HjMGKiqyluZIr6goJBIqq2FkhILCWjofkp26JAdDTVyZMOy3r1tMPvQofzUKeIH\nCgqJpES3U0Jj4xRvvw1nnQXt2zcsa9PGAqa2Nj91iviBgkIiKTUoBgyAwsLjz6d4/XULkFSlpbBl\ni+cliviGgkIiKdH1lFBQcGL305IlCgoRUFBIRKW2KOD4oDh8GCor7SinVAoKiRoFhURSbW3TQRGL\nwfLlMGIEdOly4rYKCokaBYVETix2YtcTWHAUFdk4RUWFHTLbGAWFRE2R6wJE8m33bjuS6eSTj19e\nUAAXXwzTptl1Jx5/vPHtBwywrqtYrOHwWpEwU1BI5DTWmkh47DF46SVYswYuvLDxdTp3ti6pXbvs\nvAqRsFNQSOTU1UHfvo0/d/LJMHmy3dIZMMC6nxQUEgUao5DI2bMHTjutda9RWmrdTyJRoKCQyNm9\nG049tXWvUVoKNTW5qUfE7xQUEjm7d0PPnq17jSFDYOPG3NQj4ncKComcPXtaHxTDh8OGDbmpR8Tv\nvA6KScAGoBqY0cjzNwBVwBpgKZA0T2ez24q0SC66noYOtRZFfX1uahLxMy+DohB4HPuDfyYwFRie\nss4W4EIsIB4E/m8W24q0SC5aFCefDN272xFUImHnZVCMAzYDtcBRYB5wdco6y4BP4vdXAGdksa1I\ni+SiRQEwbBisX9/61xHxOy+DohhI/r61Pb6sKdOAF1q4rUhGYrHcDGaDxikkOrw84a6RS9U3aSJw\nG5CYqzPjbcvLy7+4X1ZWRllj80KLxB06ZBcfOumk1r/WsGF2BreI31VUVFDR2LV+M+RlUOwAks9/\n7Yu1DFKNBP4LG4/Yl+W2xwWFSHNy1ZoAa1HMn5+b1xLxUuqX6JkzZ2a1vZddT5XAYKAEaAdMARam\nrNMPeA74e2xMIpttRbKWi4HshGHD1PUk0eBlUBwD7gQWAeuA+cB6YHr8BvBDoDswC1gFrGxmW5FW\nydVANtg8T0eOwN69xy9fvBj278/Ne4j4gdeTAr4YvyWbnXT/9vgt021FWiWXLYqCAhg1Ct54A669\n1pbNmQO33QaPPAL35OZtRJzTmdkSKblsUQDceSc8/LAdTbVgAdx3H/z85/DMM7l7DxHXFBQSKbkc\nzAa47jp7zV/9Cr79bfjDH2D6dNixI3fvIeKagkIiJZddTwCFhfCDH8Dtt8OPfwznnGPLpkzJ3XuI\nuKYLF0mk5LrrCeCmm+yKd8nhcP31wGO5fR8RV9SikEjJdYsC7Prb3/zm8dfPPucc+/nXv+b2vURc\nUFBIpHjRomhMm/hvVm2t9+8l4jUFhURKrgezm6Or4EkYKCgkMo4ehQMHbHrwfNmyJX/vJeIVBYVE\nxt690KNHQ7dQPigoJAwUFBIZXgxkN0dBIWGgoJDIyNdAdjIFhYSBgkIiw0WLYutWm95DJMgUFBIZ\nLloUXbrArl35fU+RXFNQSGTk+9BYgNJSHSIrwaegkMhw0fU0cKDGKST4FBQSGS66nkpLFRQSfAoK\niQwXLYoBA2xAWyTIFBQSGS5aFMXFujaFBJ+CQiLDxWC2gkLCQEEhkRCLWddTvlsUffrAzp35fU+R\nXFNQSCTs32/XjejQIb/v262bTUZ48GB+31cklxQUEgkuBrLBLmakVoUEnYJCIsHFQHaCxikk6BQU\nEgmuWhSgFoUEn4JCIsHFEU8JalFI0CkoJBJcdj2pRSFBp6CQSHDZ9aQWhQSdgkIiwXWLQkEhQaag\nkEjYtQtOP93NexcXq+tJgk1BIZFQVwd9+7p579694f33ob7ezfuLtJaCQiJh+3Y44ww3792xI3Tu\nDHv3unl/kdZSUEjoHTgAn30G3bu7q0HjFBJkCgoJve3brdupoMBdDTrySYJMQSGh57LbKaFfP3jv\nPbc1iLSUgkJCzw9B0b8/bNvmtgaRlvI6KCYBG4BqYEYjzw8DlgGHgX9Mea4WWAOsAlZ6V6KEXaLr\nyaWSEqitdVuDSEsVefjahcDjwCXADuBtYCGwPmmdvcBdwDWNbB8DyoCPPKxRIqCuDsaMcVuDWhQS\nZF62KMYBm7GWwVFgHnB1yjq7gcr4841xOPwoYaGuJ5HW8TIoioG6pMfb48syFQMWY0HyrRzWJRHj\nh66nPn3sPIojR9zWIdISXgZFrJXbnw+MBi4Hvgtc0OqKJJLq6ty3KAoL7RBZHfkkQeTlGMUOIPl7\nXF+sVZGp9+M/dwMLsK6sN1NXKi8v/+J+WVkZZWVlWZYpYXbwoH2L79HDdSUN3U+DB7uuRKKmoqKC\nioqKFm/v5RhAEbARuBjYiR25NJXjB7MTyoEDwCPxxydhg+EHgE7Ay8DM+M9ksVistQ0XCbONG+Gq\nq2DTJgdvXlAASf8/b70VJkyAadMc1CKSpMDOPs3477+XLYpjwJ3AIuyP/hNYSEyPPz8bOB07Gupk\noB74HnAmcBrwXFKNT3NiSIg0yw/dTgka0Jag8jIoAF6M35LNTrq/i+O7pxIOAqO8Kkqio6bGzmHw\ng/79YckS11WIZE9nZkuorV4No3zylaOkRC0KCSYFhYTa6tUwerTrKoy6niSogn5CmwazpUmffw5d\nu9qsrV27OiggZTD76FE4+WTYtw86dHBQj0hctoPZalFIaFVXQ69ejkKiEW3bwplnWitHJEgUFBJa\nfhqfSBg3DlZqiksJGAWFhNaqVf4Zn0hQUEgQKSgktPzaolixwnUVItnRYLaEUixm4xPvvuvwhLuU\nwWywAfbu3e3aFH6YVkSiycvB7A5A+2wLEnGhtrZhIj4/KSyEs8+Gt992XYlI5tIFRRvgOuC/sQn+\ntgLb4vd/B1xL8FskElKvvQZlZfal3m80TiFBky4oKoCzgf8DlAK9sbmZSuPLxgKve1yfSIssWQIX\nXeS6isZNmACLFrmuQiRz6b5vtQeau8xKJut4SWMUcoJYzMYlXn8dBg1yWEgjYxQAx47BwIHwu9/B\n2LEO6pLIy+UYRXIAnI3N7HoXMKaJdUR8oboa2rSxP8Z+VFQE3/se/Od/uq5EJDOZDGb/EHgS6AH0\nBOYA93tYk0irLFkCEyf6c3wiYdo0636qq2t+XRHXMvlV2gSMBA7HH3cEqoAhXhWVBXU9yQmmToXL\nLrMLBTnVRNdTwj33WOvioYfyWJMI3hweuwMLh4QOZHdJU5G8eucdO7LI7+6+G554Ag4ccF2JSHrp\nguJn8dsnwF+w7qcngbXxZSK+c/CgzRY7dKjrSppXUgIXXwxz5riuRCS9dE2PW4BY0nqp9+d6V1bG\n1PUkx1m2DO66CyorXVdCs11PAMuXw/XXw+bNNgAvkg+5vGb2k0n32wPDsIDYAHzWgtpEPFdVBV/+\nsusqMjd+vOXJhg02BbmIH2XyHeZKYDPwU6wrqga4wsuiRFqqqsp/EwE2Z/x4TRQo/pZJUPwEmAh8\nNX4rA3QEuPjS6tXBalGABcXy5a6rEGlaJkGxH2tRJGyJLxPxlfp6WLsWRo50XUl2zj1XLQrxt3Rj\nFN+I/6wEXgD+X/zx38WXifhKTQ2ceip06+a6kuyMGmVnkx88CJ07u65G5ETpguIqGo50+hDrdgLY\njZ1LIeIrQRvITmjXzlpBlZU2462I36QLilvyVYRILtTUwODBrqtomcSAtoJC/CjdGEU50CvN872B\nmTmtRqQV6uqgb1/XVbSMxinEz9K1KCqBeUA74F3gfewEjdOxGWSPYNelEPGFujo70zmIvvQl2LjR\ndRUijUsXFFOww2LvA6qBEmzM4i3gP9B8T+IzQW5RlJbCli125JbO0Ba/SRcUZwN9gMnYuRPJp3tr\n3gzxnbo6u2BREHXqBD16wPbt0K+f62pEjpcuKH4BvIpd+vSdlOdi8eUivnD4MOzfD6ed5rqSlhs8\n2OZ8UlCI36Rr5P4UGI5dqGhAyk0hIb6yfTsUFwe722bQIDufQsRvMvm1+o7nVYi0UpDHJxISLQoR\nvwnw9y+RBmEIikGDFBTiTwoKCYWwBIW6nsSPFBQSCmEJisQhsiJ+oqCQUAjyobEJnTrZhIY7d7qu\nROR4XgfFJOyKeNXAjEaeHwYsAw4D/5jltiJfCEOLAjROIf7kZVAUAo9jf/DPBKZih9sm2wvcxYlT\ngWSyrcgXtm8PT1BonEL8xsugGIdd8KgWOIrNG3V1yjq7sTmljrZgWxEADh2yE+5OOcV1Ja1XWgpb\nt7quQuR4XgZFMVCX9Hh7fJnX20rEJMYnCgqaX9fvBgxQUIj/eBkUrZkPSnNJScbCMj4BCgrxp3Rz\nPbXWDiD517cvmc84m/G25eXlX9wvKyujTFd+iZywBcWWLa6rkLCpqKigoqKixdt72VgvAjYCFwM7\ngZXYoPT6RtYtBw4Aj2S5bSwWU+Mj6mbOhGPH4MEHXVeSoqAAsvz/GYvZYbIffqjrZ4t3CqyfNuO/\n/152PR0D7gQWAeuA+dgf+unxG9hFkOqA/4Vd9+I9oHOabUVOEIZzKBIKCqB/f6itdV2JSAMvu54A\nXozfks1Our+L47uYmttW5AR1dXDdda6ryJ3EOMWIEa4rETE6M1sCLyznUCToEFnxGwWFBF6YBrNB\nRz6J/ygoJNA++cQGgLt2dV1J7igoxG8UFBJoidZEGE62S1BQiN8oKCTQwtbtBA1BoSO/xS8UFBJo\nYTo0NqFbNygqgj17XFciYhQUEmhhbFEADBkCmza5rkLEKCgk0MIaFMOGwYYNrqsQMQoKCbSwnUOR\nMHw4rNdcBOITCgoJNLUoRLynoJDAisUUFCL5oKCQwProI2jfPpyzrJaWWrfa4cOuKxFRUEiAhfHQ\n2IS2be18is2bXVcioqCQAAtrt1PCsGEa0BZ/UFBIYEUhKDROIX6goJDAUlCI5IeCQgIrrOdQJIwY\nAVVVrqsQUVBIgIW9RfHlL8N778Heva4rkahTUEhghT0oiopg/HhYutR1JRJ1CgoJpPp62LEjvIfH\nJlxwAbz5pusqJOoUFBJIu3fbiXYdO7quxFsXXghvvOG6Cok6BYUEUti7nRLGjYO1a+HgQdeVSJQp\nKCSQohIUHTvC6NGwfLnrSiTKFBQSSFEJCoCLLoIXX3RdhUSZgkICafNmGDTIdRX5ccMN8MwzcOyY\n60okqhQUEkjV1TB4sOsq8mPoUCgpgUWLXFciUaWgkEDatCk6QQFw880wd67rKiSqClwX0EqxWCzm\nugbJs88+gy5d4MABaNfOdTVpFBTY1ZVyYN8+m3Z861bo3j0nLykRVlBQAFn8/VeLQgJn61Y70c7X\nIZFj3bvD174G8+a5rkSiSEEhgVNdDUOGuK4i/9T9JK4oKCRwojSQneyyy2DbNti40XUlEjUKCgmc\nqA1kJxQV2aGyalVIvikoJHCi2qIAuPFGjVNI/ikoJHCiOkYBMHIkfPop1NS4rkSiREEhgfLpp/DB\nB9Cvn+tK3CgogEsugcWLXVciUaKgkEBZt866nYqKXFfizqWXwiuvuK5CosTroJgEbACqgRlNrPPT\n+PNVwOik5bXAGmAVsNK7EiVIli+3q75F2SWXwGuvweefu65EosLLoCgEHsfC4kxgKjA8ZZ0rgEHA\nYODbwKyk52JAGRYe4zysUwJk2TIFRZ8+0Ls3vPuu60okKrwMinHAZqxlcBSYB1ydss7fAImD/VYA\n3YBeSc8HfYoRyTG1KMyll8Krr7quQqLCy6AoBuqSHm+PL8t0nRiwGKgEvuVRjRIgu3fDnj0wPLVd\nGkHnnQcrVriuQqLCyyHBTGdDa6rVMAHYCfQEXsHGOk64zHx5efkX98vKyigrK8umRgmQFSvs0qBt\ndAgG48bBPfe4rkKCoqKigoqKihZv72XXznigHBujALgXqAf+I2mdXwAVWLcUWBh8Ffgg5bUeAA4C\nj6Qs1+yxEfKv/wqFhfCjH7muJEM5nD02VSwGPXtCVRUUp7bTRZrhp9ljK7FB6hKgHTAFWJiyzkLg\npvj98cDHWEicBHSJL+8EXAb82cNaJQA0kN2goMBaFW+/7boSiQIvg+IYcCewCFgHzAfWA9PjN4AX\ngC3YoPds4I748tOxbqbV2CD3H4GXPaxVfO7jj6GyEi64wHUl/jF2rIJC8iPoRxWp6ykinn4a5s+H\nhaltUj/zsOsJ4Pnn4dFHdfKdZM9PXU8iOfPcc3Ddda6r8JexY62VVV/vuhIJOwWF+N6hQza30VVX\nua7EX047Dbp2hc2bXVciYaegEN9btMgGbk85xXUl/qNxCskHBYX43m9/C3/7t66r8KexY2GlZkIT\njykoxNc++ghefhmmTHFdiT/pEFnJBwWF+Nq8eTBpEnTr5roSfzr7bDvp7uhR15VImCkoxNfmzoVb\nbnFdhX916QIlJbB2retKJMwUFOJbGzZAXZ3NlCpN0ziFeE1BIb41dy7ccEO0r2aXCY1TiNcUFOJL\nn38OTz0FN9/suhL/GztWU46LtxQU4kuvvQa9esGIEa4r8b9Ro2DbNjtCTMQLCgrxpblz1ZrIVNu2\nNqvumydcrUUkNxQU4jv79sEf/whTp7quJDi++lV4/XXXVUhYKSjEd556Cq64wi7MI5lRUIiXNM24\n+EosBl/6EsyaZX/8As3jacaTHTlic2Ht2GETBYqko2nGJdDeesv+tl54oetKgqV9eztM9q23XFci\nYaSgEF95+GG44w77Mi7ZKSuDV191XYWEUdB/HdX1FCJvvWUn2G3cCB06uK4mB/LY9QQ2jcfll9uh\nsm30FVDSUNeTBFIsBv/8zzBzZkhCwoERI6B7d1i61HUlEjYKCvGF2bNh/3648UbXlQTb1Kl2/Q6R\nXFLXkzj30ktw6612wtigQa6ryaE8dz0BbN1qg9o7d9qJeCKNUdeTBMqvf22tiGefDVlIODJgAAwb\nZv+eIrmiFoU4cfgw3H23nST27LMhndPJQYsCYPFiO3LsL39Rq0IapxaF+N6OHfCVr8Ann0BlZUhD\nwqFLLoF+/WDOHNeVSFioRSF5tW+fnUw3eTLcd1/Iz5dw1KIAuz7F1VfbZVI1FYqkyrZFEfRfUwVF\ngBw9at92x4yBn/wk5CEBToMC7HDjNWtsgkWdVyHJ1PUkvnX//dC5MzzySARCwgcefBA+/hgeesh1\nJRJ0usik5MUrr8BvfgOrVunbbb60bQvz58N558HQoXDtta4rkqAK+vc6dT0FwAcfwOjRFhQXXeS6\nmjxy3PWUUFlpU3v8z//YBY5E1PUkvlJfDzfdBLfdFrGQ8JFzzoEnn2wY3BbJloJCPHX//XDoEJSX\nu64k2q68En72M5g0CZYvd12NBI2CQjzz2GN2Mt2CBVCk0TDnJk+GX/4SrroKnnjCWnsimdAYheTc\nZ5/BjBnw+99DRQX07++6Ikd8MkaR6s9/httvt6D4p3+Cr38dOnVyXZXkk8YoxJlYDF54AcaOhZoa\neOedCIeEj511FixbBt//vp29XVwM118Pf/iDXVJVJJVaFJITb7wB//Ivdub1j38M11yjcyX82qJI\ntXu3dRHOm2cn6E2bZifrnXKK68rEK35rUUwCNgDVwIwm1vlp/PkqYHSW24pjlZXwta/BLbfA9On2\nh+baaxUSQdKzJ3znO9ZNWFUFBw/CkCHwox/BgQOuqxM/8DIoCoHHsT/4ZwJTgeEp61wBDAIGA98G\nZmWxbehVVFS4LqFJ69bBN75hh1xecw1s2GDThRcWZv4afv58uVDhuoAW6NsXZs2CFSvskrSDB8Oj\nj9psv8lCv+9C/vmy5WVQjAM2A7XAUWAecHXKOn8DzI3fXwF0A07PcNvQ8+N/1rVr7byIsjI747e6\nGv7hH6Bdu+xfy4+fL5cqXBfQCoMGwdNPw8svw5IlNo5xxx3w/PM262/o913IP1+2vDxosRioS3q8\nHTg3g3WKgT4ZbCt5cOQIrF9v3RILFsCmTfDd79ox+V27uq5OvDZypA1yb9sGzzxjkzlOnmwtx6VL\nbWqQkhJrifTrZz97986uZSn+52VQZDqK16re7K9/vYk3b+bdXT6f6bY1NfbLmO/aPv3UBqU/+sj6\nqAcOtNbDPffYCVvt26d/fQmf/v3h3nvtVl/f8H9h40YLkRUr4L33oK4O9uyBLl3si0TXrtCxowVH\nYaHN85V836/zfm3aZONv4r3xwEtJj+/lxEHpXwDfTHq8AeiV4bZg3VMx3XTTTTfdsrptxieKgBqg\nBGgHrKbxwewX4vfHA8uz2FZERELgcmAjll73xpdNj98SHo8/XwWMaWZbERERERGR3Aj7CXm1wBpg\nFbDSbSmt9ivgA+DPSct6AK8Am4CXsUOjg6qxz1eOHa23Kn6blP+ycqYvsAT4C7AWuDu+PCz7sKnP\nV07w92EH7NSD1cA64N/jy8Oy79IqxLqkSoC2hHMMYyu2M8PgAuys++Q/pA8BP4jfnwH873wXlUON\nfb4HgHvclJNzpwOj4vc7Y13CwwnPPmzq84VlH54U/1mEjQNPIMt959OD05oVlRPywjIRxpvAvpRl\nySdbzgWuyWtFudXY54Pw7L9d2JcxgIPAeux8p7Dsw6Y+H4RjH/41/rMd9iV7H1nuu6AGRVMn6oVJ\nDFgMVALfclyLF3ph3TXEf/ZyWItX7sIO0niC8DTtS7DW0wrCuQ9LsM+XOAIzDPuwDRaEH9DQxZbV\nvgtqUMRcF5AH52P/YS8Hvot1b4RV4tjuMJkFDMC6NN4HHnFbTk50Bp4FvgekThcYhn3YGfgd9vkO\nEp59WI99hjOAC4GJKc83u++CGhQ7sAGohL5YqyJM3o//3A0swLrbwuQDrG8YoDfwocNavPAhDb+A\nvyT4+68tFhJPAb+PLwvTPkx8vt/Q8PnCtg8/AZ4HzibLfRfUoKjEZpwtwfrdpgALXRaUYycBXeL3\nOwGXcfxAaRgsBG6O37+Zhl/OsOiddP9agr3/CrCul3XAo0nLw7IPm/p8YdiHp9LQZdYRuBQ7giss\n+65ZYT4hbwDWp7gaO1wv6J/vt8BO4DNsbOlW7IiuxYTj8LzUz3cb8Gvs8OYq7JcwyP33E7Dui9Uc\nf6hoWPZhY5/vcsKxD88C3sU+2xrg+/HlYdl3IiIiIiIiIiIiIiIiIiIiIiIiIiIiIpJ/JTR+MlYt\n4ZkBWCIsqGdmiwRBjHDMPioRp6AQyY0ibJ6gdcB/Y9MlgM35vwabbXWgm9JEWkdBIZIbQ4GfA2cC\n+7EZfwE+BkZi14Z/tPFNRUQk7EqAbUmPJ2JzA22NPwc2O+mevFYlkiNqUYjkRvJ8/gXYJHPp1hEJ\nDAWFSG70A8bH718PvBW/PyXp55/yXZSIiPhDf+w6y09x/GD2Vuyi9VXYYHapqwJFRERERERERERE\nREREREREREREREREREREREREpJX+P8MDg9ie8XF7AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1087b1e10>"
       ]
      }
     ],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc_random.p_sim_bw"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 39,
       "text": [
        "0.88900000000000001"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc_random.ww"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 40,
       "text": [
        "8.0"
       ]
      }
     ],
     "prompt_number": 40
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc_random.bw"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 41,
       "text": [
        "17.0"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jc_random.bb"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 42,
       "text": [
        "15.0"
       ]
      }
     ],
     "prompt_number": 42
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Continuous Variable"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = np.arange(w.n)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "yc = y.copy()\n",
      "yc.shape = (5,5)\n",
      "matshow(yc,cmap=cm.gray_r)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 44,
       "text": [
        "<matplotlib.colorbar.Colorbar instance at 0x10a417d40>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAPkAAADyCAYAAACPiGNTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE19JREFUeJzt3XmMHvV9x/H37hrH8cVhXK937bJ/1AegGoONsXCoN1wy\nkEIiRQFUVI4KUATliMoVItVthZqStiCuCAJEEMIRQbhNOVJKTSG+Yi84YGwOW9hrG4OJ8e5isbaf\n/vGd3X32YeZ55vnN7M48v+fzkkb7HPPM/Nbrz/Ob+c3Md0BERERERERERERERERERCSxQhXTzoza\nCEBDlisXqWGFhoZ48SkUCpBh1kZktWKRWldlyDOjkIs4ihvyrCnkIo4aGxtjzbdv374hbkl5CrmI\nI/XkIp5TyEU8p5CLeK5WQh5v5GDoLQLWARuA6zJuy/3AduDtjNsBMBV4FfgjsBa4ItvmMApYBqwB\n3gH+Ndvm9GsCVgPPDudKGxoaYk0hov6uPwPeBTqA3wIHRqx6I/AW9jsvT+83GjpNwPtAG3AA9h/o\n8AzbcwJwNPkIeTMwO3g8FniPbP9tAEYHP0cAvwe+lWFb+vwI+DXwzDCuszBmzJhYE3bWW7Gov+sp\nDHS8Pw2mMB8Bh8RtaB568nlYyDcCvcCjwFkZtmcp8HmG6y+2DfvSA+jCvuVbsmsOAD3Bz5HYF3Sm\np2wCU4DTgXsZ5rPKGhsbY00hov6uLwP7g9eXYb9blNi/ax5C3gp8XPR8c/CaDNaGbWEsy7gdjdh/\n0O3YJuc72TaHW4BrGAjHsEmwuV6sjfC/60XAkojPFIBXgJXAxZVWkIeQZ3vOX20YCzwOXIl982dp\nP7apOQX4K6A9w7Z8B/gE2zcd9lGwFEIe9Xe9EfgKeDjicwuwL4bTgMuwXcxIeQj5Fmwgos9UrDcX\ncwDwBPAQ8FTGbSm2C3gemJthG44HzsT2UR8BTgQeHK6VR4V679697Nmzp3+KEPV3vQDb/fibMqve\nGvzcATyJ7fLm2gjgA2yzZSTZD7yBtSUPA28N2H/aW7JuSOBQ4KDg8TeB/wVOyq45gyxkeEfXCwcf\nfHCsia9vrUb9XRdhI+6HllnvaGBc8HgM8H/AqeUamoeefC9wOfAitn/3GDYQkZVHgDeA6dhYwYUZ\ntmUBcB7wbWyTdDX2HyErk4H/xr6Il2Gh+l2G7Sk1rLt+CTbXw/6upwG3Y5vwLwev3RXM34JtNYGN\nzC9l4G/wHPBS2XYm+SVF6lhhwoQJsWb87LPPQNeTi9SeuFehZU0hF3FUK6e1KuQijhRyEc/VSsgT\n71QsXLiwmqqVmjTldpo7d26BKqR0xtuQS6MFhTQK1S1evJjFixcnXk5PT0/lmWK46aabuPHGGxMv\np7u7O4XWwM0338y1116beDlptOfWW2/lqquuSrwcSOfvdeedd3LZZZclXs6RRx4J8TNRaGmJdxlB\nZ2dnNctNnTbXRRzloZeOQyEXcaRDaFVqb2/PugmDnHBC2XP+h92CBQuybkK/+fPnZ92EQY499thM\n1lsrPXlu9snTktY+eVrS2idPS97ak6e/V7X75IcddlisGTdt2lTNclOXm55cpNbUSk+ukIs4UshF\nPKeQi3iuVkbXa6OVIjk0BCWZD8GuJV+PXSN+UNiHqbKEuUIu4ihByHuBq4EjgflYnbbDgeuxkE/H\ninFcH/LZJuAOLOhHAOdSoZKSQi7iKEHIw0oyt2L16h4IXn8A+G7IZ6suYR4n5Hm6u4lIbgxBSeZJ\nWKlrgp+TQuavuoR5pYG3vk2Dk7Gqqiuwu1RkWYNNJBdSGF0fi1VsvRLYXfJe39Vxpao+86xSyIs3\nDWBg00Ahl7oXFfLu7u44Zxb2lWT+FQMlmbdjhRq3YUUzPwn5XNUlzCttruvuJiIRom6LNG7cOJqb\nm/unEA3AfVh14luLXn8GOD94fD7hdfZXAtMYKGF+NhXuAVcp5Pk5KV0kZ1IuybwIu8HhKdghtBMZ\nuOFhcUnmqkuYV9pcj7VpUFzsob29PXdXlImEWb58OStWrHD+fIJ98teJ7mBPDnmtEzij6PkLwRRL\npVaOwG6relKwouXYcbnibw5dhVZG3q76ylt78vT3qvYqtKOOOirWjB0dHdUsN3WVevLiTYMmbD9C\ng24i+HXuelWbBiL1wqeQi0gIhVzEc7VyFZpCLuJIPbmI5xRyEc8p5CKeU8hFPKeQi3hOIRfxnA6h\niXhOPbmI5xRyEc8p5CKeSxDy+7Hrwz8B/jJ47VFgRvD4IOBPWIHHUhuBL4B9WLXWeZVWlkrIk1x4\nnzZdL11e3v598taeaiQI+S+B24EHi147p+jxv2MhD1MA2oGdcVemnlzEUYKQL8VqtIUuFvgBVhoq\nctXVrKw2jgGI5FBUIcfSqUonYFVbP4h4vwC8ghV0vDjOAtWTiziK6sl37tzJ559/7rrYc4GHy7y/\nANgKTMRuqbQO2zKIpJCLOIoK+YQJE5gwYUL/8w8//DDuIkcA3wOOKTPP1uDnDuBJbOCtbMi1uS7i\nKKXbJBU7Gauh2Bnx/mhgXPB4DHAq8HalhSrkIo4ShPwR4A3s7qUfAxcGr58dvFesuOZ6M9Zrr8Hu\nnfYcdovjsrS5LuIowej6uRGvXxjyWnHN9Q+B2dWuTCEXcaQLVEQ8p9NaRTynkIt4TiEX8ZxCLuK5\nWgl5nOHB+7FzaSsedBepJ0NwMsyQiBPyX2I3SBeRIkN0gUrq4myul7ssTqRu5aGXjkP75CKO6irk\n99xzT//jOXPmMGfOnDQWKzKk1q1bx7p165w/X1chv+SSS9JYjMiwmjlzJjNnzux//vTTT1f1+boK\nuUg9qpWQxxn6i7osTqSu+XQI7VzsmtZvAFOxQ2oidS/BIbSwc08WA5uB1cEUddh6EVbyaQNwXax2\nxvt1RKRUgp487NyTAvCfWK31o4H/CvlcE3BH8NkjsA748ErtVMhFHCUI+VIgrNJjpW37ecD72A0W\nerEbMpxVqZ0KuYijIdgn/3ugA7gPu4tKqVZsXKzP5uC1sjS6LuIoKsDbtm1j27Zt1S7u58A/B4//\nBfgP4O9K5ilUu1BQyEWcRYV88uTJTJ48uf95R0dHnMV9UvT4XuDZkHm2YIPffaZivXlZ2lwXcZTy\n5vrkosffI/yqz5XANOxakpFYdddnKi1YPbmIowRXmD0CLAQOxfax/xG7ieFsbJP8I+DSYN4W4BdY\nxda9wOXAi9hI+31YnfayFHIRRymXZL4/Yt7ikswALwRTbAq5iKM8nM0Wh0Iu4kghF/GcQi7iOYVc\nxHN1FfIVK1aksZhUdHd3Z92EQfLWnp6enqybMEhXV1fWTXCWhyKNcagnF3FUVz25SD1SyEU8p5CL\neE4hF/GcQi7iOY2ui3hOPbmI52ol5LWxvSGSQwmKRoSVZP4Zdm14B/Bb4MCI1W4E3sLKNi+P006F\nXMRRyiWZXwKOBI4C1gM3RKy2gBWYOBqr3lqRQi7iKOWSzC8D+4PHy4Ap5VZdTTsVchFHQ3ibpIuA\nJRHvFYBXsHpvF8dZmAbeRBwN0SG0G4GvgIcj3l8AbAUmYr3/OmzLIFKckE8FHgT+DPsWuQe4LV57\nRfwV1Utv2rSJTZs2uSzyAuB04KQy82wNfu4AnsT2yxOHvBe4GlgDjAVWYd8gFatEivgsKuRtbW20\ntbX1P1+6tGwG+ywCrsGquO6JmGc0VqV1NzAGOBX4p0oLjrO9sQ0LOEAXFu6WGJ8T8VqCffK+24HP\nwEoyXwTcjnWiL2OHx+4K5m0Bng8eN2O99hpscO45bFS+rGr3yduwoftlVX5OxDsZlGT+EKvNXpVq\nQj4WeBy4EuvR+z3//PP9j6dNm8b06dOrbYfIsNuyZQudnZ3On6+VM97ihvwA4AngIeCp0jfPOOOM\nr31AJO9aW1tpbR24KejKlSur+rxPIW/AbsfyDnDr0DZHpHbUylVocVq5ADgP+DY2ILCar5+SJ1J3\nhvBkmFTF6clfR2fGiXxNHgIch854E3GkkIt4TiEX8ZxCLuI5hVzEc7VyCE0hF3GknlzEcwq5iOcU\nchHP1UrIa2PkQCSHUi7JfAh2Lfl67BrxgyJWuwgr+bQBuC5OOxVyEUcpl2S+Hgv5dOB3wfNSTcAd\nwWePwK5LP7xSOxVyEUeNjY2xphBhJZnPBB4IHj8AfDfkc/OA97EbLPQCjwJnVWpnKvvky5fHupHD\nsOjp6cm6CYN0dXVVnmkY5e3fJ2/tqUbK++STsE14gp+TQuZpxcpF9dkMHFdpwRp4E3E0hANvhWAK\ne71qCrmIo6iQr1+/ng0bNlS7uO1YocZtwGTgk5B5tmAl0vtMxXrzshRyEUdRIZ8xYwYzZszof75k\nSdTNUAZ5Bjgf+Lfg59fKrGF3TZmGFVTtBM4mvCjkIBp4E3GUYknmC4GfAqdgh9BODJ7D4JLMe4HL\ngRexcmyPEeP+B+rJRRwluEAlqvc9OeS14pLMAC8EU2wKuYijWjnjTSEXcaSQi3hOIRfxnEIu4jmF\nXMRzCrmI51TjTcRz6slFPOdTyEcBrwHfAEYCTwM3DGWjRGqBTyHfg93RtCeY/3XgW8FPkbrlU8jB\nAg7WkzcBO4emOSK1o1ZCHnd4sBFYg13z+ip2BYxIXfPp/uQA+4HZwIHYZW7twP/0vblmzZr+GZub\nm2lubk6tgSJD5YsvvmD37t3On/f1ENou7NrWuRSFfPbs2Sk2SWR4jB8/nvHjx/c/7+zsrOrzeeil\n44jzVXQoAzWgv4ld2L56yFokUiMSbK7PwDLUN+0CriiZpz14vW+en7i2M05PPhkrEdsYTL/C6kKL\n1LUEPfl7wNHB40asdtuTIfO9hpVqTiROyN8Gjkm6IhHfpLS5fjLwAYNLLfevIo0V1MbIgUgOpTS6\nfg7wcMjrBeB4oANYgt0xxYlOaxVxFBXgtWvXsnbt2jiLGAn8NeH3NPsDVnK5BzgNq9463aWdCrmI\no6hDaLNmzWLWrFn9zx977LGoRZwGrAJ2hLxXfGzvBeAu7KaIVZ+IppCLOEphn/xcrDxzmEnYDRYK\n2D3QGnA801QhF3GUMORjsEG3i4teuzT4eTfwfeCHWK31Hmzf3YlCLuIoYci7sXNQit1d9PjOYEpM\nIRdxVCtnvCnkIo4UchHPKeQinvP1KjQRCagnF/FcXYV8+fLlaSwmFT09PZVnGkZ5a09XV1fWTRjk\nyy+/zLoJzuoq5CL1SCEX8ZxCLuI5ja6LeE49uYjnFHIRzynkIp5TyEU8lzDkG4EvgH1AL1YYotRt\nWPWYHuACHEuhK+QijhKGvIDVVo+q9nI68BfANOA44OfAfJcVKeQijlI4hFbuW+JM7H4HAMuwG5xM\nwu5HWJXaONAnkkMJSzIXgFeAlQwuAdWnlcG12DcDU1zaqZ5cxFFUgFetWsWqVasqfXwBsBWYCLwM\nrAOWlq6i5Hmh+lYq5CLOokI+d+5c5s6d2//83nvvDZtta/BzB3aLpHkMDvkWrO56nynBa1WLu7ne\nhI3sPeuyEhEfJdhcHw2MCx6PAU7FbkdW7Bngb4PH84E/4bA/DvF78iuBd4oaJlL3EoyuT2LgBocj\ngF8DLzG4JPMSbIT9fayy64WuK4sT8inBym4CfuS6IhHfJAj5R8DskNfvLnl+uesKisUJ+S3ANcD4\nSjOK1JNauQqtUiu/g92qZTUp3UZVxBcp3dV0yFXqyY/HDsqfDozCevMHGRgQAGDHjoH7tY0ePZox\nY8ak20qRHMpDgOOoFPIfBxPAQuAfKAk4wMSJE1Nulkj++RLyUk4H40V85GPIXwsmEcHPkItIEYVc\nxHO1cghNIRdxpJ5cxHMKuYjnFHIRzynkIp5TyEU8Vyshr41jACI51NjYGGsKMRV4FfgjsBa4ImSe\ndmAXdnHYauAnru3MTU/e3d2dqwtb9uzZw6hRo7JuRr+vvvqKkSNHZt0MAPbt20dTU1PWzchcgp68\nF7gaWAOMBVZhdd7eLZnvNewCsURy05P39PRk3YRB9uzZk3UTBunt7c26Cf3279+fdRNyIcGlptuw\ngAN0YeFuCVtFGu3MTchFak1K15O3AUdjtdWLFbBLvTuwUlBHuLYzN5vrIrUmhYG3scDjWA3FrpL3\n/oDtu/dgt0p6CpjuspI0Qv7au+++uzCF5fDpp5+msZjU7Nq1K+smDJKnXZq9e/dm3YShUNVVllEh\nf+ONN3jzzTcrffwA4AngISzApXYXPX4BuAs4hOjbKkW3s9oPiAgAhc7OzlgztrS0wOCsNWC3QPoM\nG4ALMwkrvVbAarL/Btu0r5o210UcJbgKbQFwHvAWA3cq/THw58Hju4HvAz8E9mKb7Oe4rkw9uYib\nwvbt8e51MGnSJMgwa+rJRRzVyhlvCrmII4VcxHMKuYjnFHIRz6nGm4jn1JOLeE4hF/GcQi7iOYVc\nxHMKuYjnFHIRz+kQmojn1JOLeK5WQl4b2xsiOZSwxtsiYB2wAbguYp7bgvc7sDpwIjKMCr29vbEm\nrLpLsSbgfazSywFY5dbDS+Y5HSvgCHAc8HvXhqonF3GUoCefh4V8I1aD/VHgrJJ5zsRKRIFVcj0I\nKwlVNYVcxFGCkLcCHxc93xy8VmmeKS7t1MCbiKMEh9BKN9+jlH5DxP3cIAq5iKMqRtd3lzzfgtVU\n7zMV66nLzTMleE1EasAI4ANs4G0klQfe5pNg4E1EsnEa8B42AHdD8NqlwdTnjuD9DuCYYW2diIiI\niIiIiIiIiIiIiIiIiEge/T/qN7t+EnVvJAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x108427e10>"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Moran's I\n",
      "\n",
      "$$I = \\frac{n}{S_0} \\frac{\\sum_i \\sum_j z_i w_{i,j} z_j}{\\sum_i z_iz_i}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi = ps.Moran(y,w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi.I"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "0.84000000000000019"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi.EI"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 47,
       "text": [
        "-0.041666666666666664"
       ]
      }
     ],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi.p_norm"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 48,
       "text": [
        "4.3976446928439827e-09"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mi.p_sim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 49,
       "text": [
        "0.001"
       ]
      }
     ],
     "prompt_number": 49
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "mi.sim"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "density = gaussian_kde(mi.sim)\n",
      "xs = linspace(mi.sim.min(),mi.sim.max(),200)\n",
      "plot(xs, density(xs))\n",
      "axvline(x=mi.I, color='r')\n",
      "xlabel('I')\n",
      "ylabel('f(I)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 50,
       "text": [
        "<matplotlib.text.Text at 0x10a613210>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEPCAYAAACtCNj2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XnclXP+x/HX3aJFJTSKVsVYmkqiTXRCNVKyhCw14ke2\n4RfCMHQ38zP8ZoxtMMWvbKFpLFGEajqptEy0L6MilFRStFii8/vjc+7T7XQv577vc13f67rO+/l4\nnMd9lss57273OZ/zXS8QEREREREREREREREREREREREREREJrOrAXGAhsBy4t5jjHgFWAYuAtv5E\nExERF2omf1YB5gBd0h7vBbyZvN4heYyIiDhSyePn35X8uR9QGfgq7fGzgGeS1+cCdYH6HmcSEZFi\neF0UKmHdRxuBaVg3UmENgc8K3V4HNPI4k4iIFMProrAHOA77oD8FiBVxTF7a7YTHmUREpBhVfHqd\nr4E3gBOAeKH71wONC91ulLzvZ1q0aJFYs2aNl/lERKJoDXBEWf4DL1sK9bAxAoAaQHdgQdoxrwMD\nk9c7AtuwrqafWbNmDYlEIvCXYcOGOc8QlZxhyKicylnsBZxnTCQSAC3K+sHtZUvhUGwQuVLy8hww\nFRicfHwkNvOoF7Aa2AkM8jCPiIiUwsuisAQ4voj7R6bdvt7DDCIiUgZeDzTnlFgs5jpCRsKQMwwZ\nQTmzTTndS5/5E1SJZP+YiEjw5eVBAD6z8vLyoIyf82opiIhIioqCiIikqCiIiEiKioKIiKSoKIiI\nSIqKgoiIpKgoiIhIioqCiIikqCiIiEiKioKIiKSoKIiISIqKgoiIpKgoiIhIioqCiIikqCiIiEiK\nioKIiKSoKIiISIqKgoiIpKgoiIhIioqCiIikqCiIiEiKioKIiKSoKIiISIqKgoiIpKgoiIhIioqC\niIikeFkUGgPTgGXAUuCGIo6JAV8DC5KX33uYR0RESlHFw+feDQwBFgK1gPeBycCKtOOmA2d5mEN8\nNHs2LFoExx8P7du7TiMiZeVlS+ELrCAA7MCKwWFFHJfnYQbxye7dMHAgXHopzJsH/fvDeefB9u2u\nk4lIWfg1ptAMaAvMTbs/AXQGFgFvAsf6lEey7OabYdMmWLIERo+G5cvh4IPh9NNh507X6UQkU358\nS68FxIH/AcanPVYb+AnYBZwBPAz8sojnSCQSCQ8jSkVMnw4DBsDixVC37t77Ewm4/HLYsQPGjYM8\ntQklV+Tl2RvAeYw8KOPnvNdv06rARGAS8FAGx38MtAO+Srs/MWzYsNSNWCxGLBbLUkSpiEQC2rWD\nW2+1LqN0338PsRj06gV33eV7PBE3HBWFeDxOPB5P3R4+fDgEqCjkAc8AW7AB56LUBzZh3UjtgXFY\nV1M6tRQC6u23YehQG1wuriWwYQOceCKMGgU9e/qbT8QJtRSK1AV4F1iMfegD3AE0SV4fCVwHXAP8\niHUh3QTMKeK5VBQCqmdPuOgiuOyyko+bPh0uvBDmz4dGjXyJJuKOioLnVBQCaO1aawGsWwfVqpV+\n/J/+BJMmwb/+BVWreh5PxJ0QFwWtaJZyGzMGLrggs4IAcPvtUKsW3Hmnt7lEpPxUFKRcEgl47jmb\ndZSpSpXsv3nlFXjsMe+yiUj5ebmiWSJs6VKbWdShQ9n+u3r1YPJkOPlkOOAAW+wmIsGhoiDl8uqr\ncPbZ5Vt7cPjhNmupZ09YtQruvhsqV85+RhEpO3UfSbmMHw/nnFP+/75lS5uJFI9Djx62AlpE3FNR\nkDL74gv4+GM46aSKPU+DBjBlCvTubQvcrr4a1q/PSkQRKScVBSmzKVOgWzeokoXOx6pVYcgQWLkS\nateG1q1tMdyWLRV/bhEpOxUFKbN33rEun2w66CD4y19sQ70dO6x76dVXs/saIlI6LV6TMkkkoGFD\nmDEDWrTw7nXmzLEV0FddBXfcoc30JGRCvHhNs4+kTJYtg+rVvS0IAB07wty5tvX2jz9Cof0QRcRD\nKgpSJl50HRWnQQOYOhU6d7ZprAMH+vO6IrlMYwpSJpMnQ/fu/r1e/frw+ut2Ep8PP/TvdUVyVVh6\najWmEAA//mgDwmvX2k8/PfooPPsszJqlzfQkBEI8pqCWgmRs0SJo0sT/ggBw3XV2es977/X/tUVy\niYqCZGzmzIovWCuvvDx44gl4+GFrqYiIN1QUJGOzZkGXLu5ev3FjuPFGO/WniHhDYwqSkYL1CTNn\nQvPm7nLs2gXHHGNbcJ9yirscIiXSmIJEXUGXzeGHO41BzZpw331w222BeM+JRI6KgmRk1iwbTwjC\nyuILLoBvvrHtt0Uku1QUJCMzZ7odTyiscmVb4TxsmFoLItmmoiAZKWgpBEW/frBzJ0ya5DqJSLQE\noDMgIxpodmjrVmjaFL76KjvbZWfLP/9pO6vOnRuMbi2RFA00S5TNng0nnhisggBw3nnw3Xfw5puu\nk4hEh4qClCpI4wmFVaoEd90F99wTiC9lIpGgoiClCtp4QmHnngubNllGEam4sPTEakzBkR9+sL2O\nPv8c6tRxnaZoI0bAG2/AhAmuk4gkaUxBouqDD+DII4NbEAB+8xuYPx+WLnWdRCT8VBSkREEdTyis\nRg347W9tJpKIVIyKgpQoyOMJhV1zjXUfffqp6yQi4eZlUWgMTAOWAUuBG4o57hFgFbAIaOthHimj\nRML9zqiZOvBAGDQIHnzQdRKRcPNyoLlB8rIQqAW8D5wNrCh0TC/g+uTPDsDDQMcinksDzQ58+KGd\nevOTT1wnycy6ddC6Naxe7eZEQCIpGmgu0hdYQQDYgRWDw9KOOQt4Jnl9LlAXqO9hJimDMIwnFNao\nEZx1Fvz9766TiISXX2MKzbCuoblp9zcEPit0ex3QyKdMUoqwjCcUdtNN8PjjNpVWRMrOj40LagEv\nATdiLYZ06U2bIttc+fn5qeuxWIxYLJaddFKsmTPhhuJGggKqdWs4+mjbF+mSS1ynEfFXPB4nHo9X\n6Dm8XrxWFZgITAIeKuLxEUAcGJu8vRLoCmxMO05jCj7bvNnWJ2zZYltVh8mECTB8OPz739ooTxzR\nmEKR8oBRwHKKLggArwMDk9c7AtvYtyCIA7NmQceO4SsIAGeeCV9/ra0vRMrDy+6jk4BLgcXAguR9\ndwBNktdHAm9iM49WAzuBQR7mkTIIy1TUolSqBDfeCA89FN5/g4grYWlcq/vIZ506wZ/+BN26uU5S\nPjt2QLNmtv1Fs2au00jOCXH3kYqC7GPnTqhf38YVatRwnab8hg6FPXvgr391nURyToiLgra5kH28\n9x60bRvuggBw/fXw9NOwfbvrJCLhoaIg+5g+HU45xXWKimvaFE49FZ55pvRjRcSoKMg+pk+Hrl1d\np8iOa6+FkSMD0ZIXCQUVBfmZb7+FBQugc2fXSbIjFrPVze+95zqJSDioKMjPzJkDv/oV1KrlOkl2\n5OXB4MHWWhCR0qkoyM9EqeuowG9+A6+/bquzRaRkKgryM1EsCgcfDL17w5gxrpOIBJ/WKUjK99/b\nB+jnnwf7nMzlMWUK3HYbvP++6ySSE7ROQaJg3jzbYTRqBQFsZfbGjbB0qeskIsGmoiAp774bva6j\nApUrw4AB8OyzrpOIBJuKgqTE49FYtFacgQNtXOHHH10nEQkuFQUBbH3CnDnh3QAvE8ccY6fsnDrV\ndRKR4FJREABmzIA2baI5nlBY//4wbpzrFCLBpaIgALzzDvTo4TqF9/r1g9deg927XScRCSYVBQFy\npyg0aQJHHAHTprlOIhJMKgrChg2wbh2ccILrJP44/3z45z9dpxAJJhUFYcoU22K6ipcnZw2Qfv1g\n/Hh1IYkURUVBePtt6N7ddQr/NG0KzZvblh4i8nMqCjnup5/grbegVy/XSfzVty9MmOA6hUjwqCjk\nuNmzbe5+48auk/ird28rCgHYnkYkUFQUctzEidCnj+sU/mvVylpJK1a4TiISLCoKOW7CBPvWnGvy\n8va2FkRkLxWFHPbRR/Dll3Diia6TuNGnj7WURGQvFYUcNnEinHkmVMrRv4JYDBYv1hnZRArL0Y8D\nAdvuIRfHEwpUr24bAL71luskIsGhopCjNm+G+fOhZ0/XSdzq0QMmT3adQiQ4MlnDWhfoBDQDEsBa\nYDbwtWepxHOvvWYfiDVruk7iVo8ecM89NjU1LywnpxXxUEkthZOB14F3gf5AE6wwXATMSD7WpZTn\nHw1sBJYU83gMKy4LkpffZxZbKurll227h1zXogVUqwbLl7tOIhIMJbUUzgFuBlYV8/gvgauBmSU8\nx1PA34CSToI4HTirhMcly7ZuhVmzdF4BsNZB9+62S2zLlq7TiLhXUkvhJoovCAAfJo8pyQxgaynH\nqNHuswkTbAO82rVdJwmGHj2sKIhIyS2Fm0t4LAE8kIXXTwCdgUXAeuAWQA15j40bZ2cgE3PqqXDZ\nZfDddzYjSSSXlVQUamMf2l76AGgM7ALOAMZj3VL7yM/PT12PxWLEYjGPo0XT5s0wcyaMHes6SXAc\neKB1Hb33nhUIkbCKx+PE4/EKPYcfXTfNgAlAqwyO/RhoB3yVdn8ioZ3LsuKxx+zD7/nnXScJljvu\ngMqV4Y9/dJ1EIiEvLxC7LebZlLoyfc6XNKaQD9Qv4fFDgeFlebEi1Gdv4PbJ6+kFQbJozBi49FLX\nKYInFoMKfsESiYSSuo/mA2OB/bBung3Yh3YD4Hjge+D+Up7/RaArUA/4DBgGVE0+NhLoB1wD/Ih1\nIamn20OrV9t+R7l0Qp1MnXQSLFgAu3Zp7YbktpKaFc8BA7C1A6vYu3jtE2AWsM7rcIWo+ygLhg+H\nr76Chx92nSSYunSB/Hw4/XTXSST0Qtx9VFJLoR1wGHABtsis8BO7/9dKmSQS1nX0wguukwRXLAbT\npqkoSG4rqSiMAKYCzYH30x5LJO+XkJg3z3ZDPeEE10mCq1s3uPtu1ylE3MqkWTECW7nskrqPKui3\nv4VDDoG77nKdJLh27bLf0RdfQK1artNIqIW4+yiTXVJdFwSpoN274R//gEsucZ0k2GrWhOOPty1A\nRHKVts7OAe+8A0ccAc3V4Veqbt00NVVym4pCDtDahMwVDDaL5KqwbEanMYVy+uYbaNwY1qyBevVc\npwm+776z39OGDdowUCog4mMKEmKvvgpdu6ogZKp6dZuhNbOkDeFFIkxFIeLUdVR2sRhMn+46hYgb\nKgoR9vnndh7mPn1cJwkX7YMkuUxFIcLGjoVzzoEaNVwnCZcOHWDpUtixw3USEf+pKESYuo7Kp0YN\naNdO6xUkN6koRNSyZbBpkw0yS9l17apxBclNKgoR9fzzcNFFduIYKTuNK0iu0jqFCEok4PDDYfx4\nOO4412nCqWAfpI0bYf/9XaeR0NE6BQmS2bNtH582bVwnCa+aNaFtWzt1qUguUVGIoBdegIsvti8r\nUn7qQpJcpKIQMbt3w7hxNp4gFaPBZslFKgoRM3Wq7YbaooXrJOHXqRMsXAg7d7pOIuIfFYWIKeg6\nkorbf38bqJ8923USEf+oKETIrl0wYQJceKHrJNGhLiTJNSoKETJxIrRvD/Xru04SHRpsllyjohAh\n6jrKvs6dYcECa4WJ5AIVhYjYutXOGHbOOa6TRMv++0Pr1jBnjuskIv5QUYiIl1+GHj2gTh3XSaJH\nXUiSS1QUIkJdR97RYLPkkrCsedXeRyVYvx5atbKT6lSv7jpN9OzYAQ0awObNOjeFZEh7H4lLBSfT\nUUHwRq1ato+UztssucDrojAa2AgsKeGYR4BVwCKgrcd5IkldR97r0QMmT3adQsR7XheFp4Bfl/B4\nL+AI4EjgKuDvHueJnJUrrdsoFnOdJNq6d4d33nGdQsR7XheFGcDWEh4/C3gmeX0uUBfQ0qsyePFF\n6N9fJ9PxWvv28Mkndn4FkShzPabQEPis0O11QCNHWUInkVDXkV+qVIFu3WDKFNdJRLxVxXUA9h0Z\nL3LIPj8/P3U9FosRU38J8+fbzxNOcJsjV/ToYV1Il1ziOolI0eLxOPEKLqrxY0pqM2AC0KqIx0YA\ncWBs8vZKoCs2OF2YpqQWYcgQW6w2fLjrJLlhzRro0sXGcHQCIymRpqSW2+vAwOT1jsA29i0IUoSf\nfrKpqDqZjn9atLDTdC5d6jqJiHe87j56EfvmXw8bOxgGVE0+NhJ4E5uBtBrYCQzyOE9kxONw2GFw\n9NGuk+SWgqmprYpq94pEQFgaweo+SnPFFdCyJdx0k+skueXVV2HkSHjrLddJJNBC3H2kohBC331n\nrYQlS6BhQ9dpcsu2bdC4MWzapC0vpAQhLgquxxSkHCZNsm0XVBD8V7cutG1r25SLRJGKQgiNGQOX\nXuo6Re7q08dOeyoSReo+Cplt26BpU1tdW7eu6zS5acUKG3D+9FNNTZViqPtI/PLSS7YPjwqCO0cf\nDdWqwcKFrpOIZJ+KQsio68i9vDx1IUl0qSiEyKef2sKpM85wnURUFCSqVBRC5MUX4bzzrOtC3Dr5\nZNv2Yt0610lEsktFISQSCXjuOXUdBUXVqtZaeOUV10lEsktFISQWL7ZzBZ90kuskUuD8823gXyRK\nVBRCYswYO29CJf0fC4zu3W1V+YYNrpOIZI8+YkLgp59sPEFdR8FSrRqceabthyQSFSoKIRCPwy9+\nAcce6zqJpOvXT11IEi0qCiEwejQM0qbigdSzJyxYoC4kiQ4VhYDbuhXeeEOngAyqGjXg7LOte08k\nClQUAm7sWPs2evDBrpNIcQYMsIkAIlGgohBwo0bB5Ze7TiEl6drVzq+wbJnrJCIVp6IQYIsW2YfN\n6ae7TiIlqVzZuvfUWpAoUFEIsNGj4bLL7ENHgu3SS60o/PST6yQiFaOiEFDffw8vvGBFQYKvVSto\n0AAmT3adRKRiVBQC6pVXoHVraN7cdRLJ1JVXwhNPuE4hUjFhOW9Uzp15rUsXGDLEdkWVcNi+HZo0\ngeXL4dBDXacRp3TmNcmmhQvtdJt9+7pOImVRu7atcH76addJRMpPRSGAHnsMBg+GKlVcJ5GyuvJK\nePJJDThLeKkoBMzWrbaXzpVXuk4i5XHiiXDQQTBpkuskIuWjohAwTz9tp9usX991EimPvDy48UZ4\n+GHXSUTKRwPNAbJnDxx1lBUGnUwnvL7/Hpo1gylToGVL12nECQ00SzZMmgS1akHnzq6TSEVUqwZX\nXw1/+5vrJCJl53VR+DWwElgF3FbE4zHga2BB8vJ7j/ME2l/+ArfcYl8yJNyuvhrGjbNtSkTCxMui\nUBl4FCsMxwIXAccUcdx0oG3y8j8e5gm0efPg44/hggtcJ5FsqF8f+veHhx5ynUSkbLwsCu2B1cBa\nYDcwFihq5r2+F2OthCFDoGpV10kkW4YOhZEjYds210lEMudlUWgIfFbo9rrkfYUlgM7AIuBNrEWR\nc1avtlNu/td/uU4i2XT44dC7Nzz6qOskIpnzsihkMvT+AdAYaAP8DRjvYZ7AeuABuOoqG2SWaLn9\ndnjkEdi503USkcx4uWZ2PfaBX6Ax1loobHuh65OAx4GDgK/Snyw/Pz91PRaLEYvFshTTrU2b7Oxq\ny5e7TiJeOOYYOOUU2yhvyBDXaSTq4vE48Xi8Qs/hZX9+FeA/wGnA58A8bLB5RaFj6gObsFZFe2Ac\n0KyI54rsOoVbb4UdO+Dxx10nEa988AH06QOrVkHNmq7TiC9CvE7B60HeM4CHsJlIo4B7gcHJx0YC\n1wHXAD8Cu4CbgDlFPE8ki8KmTXD00bB4MTRq5DqNeOmCC6BtW/jd71wnEV+oKHgukkVh6FD49lsN\nROaCVaugUyf4z3/g4INdpxHPqSh4LnJFYeNG629WKyF3XHst1KgBf/2r6yTiORUFz0WuKNxyi+2R\no60QcscXX9heSB98AE2buk4jnlJR8FykisKGDfbhsGQJNExfuSGRdvfdsHYtPPus6yTiKRUFz0Wq\nKFx1FdSpA/ff7zqJ+G37dptc8NJLNsYgEaWi4LnIFIVly6BbNxtwPPBA12nEhRdesG1N5s+HypVd\npxFPhLgoaOtsn912m61yVUHIXRddBAccACNGuE4isi+1FHw0bRpccQWsWGF77kvuWroUTj3Vfh5y\niOs0knUhbimoKPhkzx5o395mHfXv7zqNBMHNN8OWLXamPYmYEBcFdR/5ZNQo2G8/nS9B9srPt91x\nJ092nURkL7UUfLB5s01BnTwZ2rRxnUaC5O23YfBgm55cu7brNJI1IW4pqCj4YNAgqFsXHnzQdRIJ\nokGDYP/9td1JpKgoeC60RWHGDJttsmKFvglK0bZuhV/9yrZQP/lk12kkK0JcFDSm4KEffrD9bh58\nUAVBinfggbZ1+qBBtrhNxCW1FDz0+9/bhnevvWZfHERKcsUV9nPUKLc5JAtC3FIIy0dV6IrC3LnQ\nty8sXAgNGrhOI2GwYwccdxz8+c9w7rmu00iFhLgoqPvIAzt2wMCBtgOqCoJkqlYtGDMGrrkG1q93\nnUZylVoKWZZIwIABULUqPPWU6zQSRn/4w971C9obKaTUUpACTz4JixbBY4+5TiJhdeedUKmSbbMt\n4je1FLJo5kzrC54xA446ynUaCbPNm6FdO/ty0aeP6zRSZmopyKpV0K+f9QmrIEhF/eIX8I9/2Iyk\njz5ynUZyiYpCFqxfD716wR//CD16uE4jUdGpk3Uh9e0L33zjOo3kCnUfVdD69XbSnCuusHMliGRT\nIgHXXWethYkToUoV14kkI+o+yk0rV0LXrnDllSoI4o28PHjkEft59dWB+JyRiFNRKKdJk+CUU2ym\nyNChrtNIlFWpAuPG2alchwxRYRBvqfuojHbutNNpvvKKvVFPOsl1IskV27ZZV+Wpp9o5nivpK11w\nqfso+vbssROut2plu1ouWaKCIP6qWxemTIH33oPLLoPdu10nkihSS6EU27fbNNNHH7WdTu+7D2Ix\nJ1FEANi1Cy680P42X3wRDj3UdSLZh1oK0ZFIwIcfwujRcM450LAhTJ1qg32zZ6sgiHs1a8L48fa3\n2K4dvPmm60QSJV63FH4NPARUBv4P+N8ijnkEOAPYBVwGLCjiGM9aCrt2wfz51iQvuNSqZXPEzzgD\neveGgw7y5KVFKmzqVDudZ5s2Ns7QvLnrRAKopVCMysCjWGE4FrgIOCbtmF7AEcCRwFXA3z3MkzJz\nps3i6NDBVo4OHQqbNtnOposWwdq11iwfOLBsBSEej3sVOavCkDMMGcF9ztNOg6VLrSi0bw8XXwyz\nZtkYWGGuc2ZKOd3zsii0B1YDa4HdwFigb9oxZwHPJK/PBeoC9T3MBMCGDbal9f33w5df2rkPHnjA\ntqlo2LD8zxuWP5Qw5AxDRghGzurVbeXzRx/B8cdby6FZM/v53HN2To+33467jpmRIPw+MxGWnOXh\n5frIhsBnhW6vAzpkcEwjYKOHuTj/fC+fXcSNOnXgllvssmyZdS1NmGDdSitXwhNPwCGH2ISJOnXs\nUrMm1KgBd90FTZq4/hdIEHhZFDLtUEvv73LfEScSci1b2uWGG+z2sGG2XcaWLbaPUsFl1y749lsr\nDiLg7UBzRyAfG1MA+B2wh58PNo8A4ljXEsBKoCv7thRWAy08yikiElVrsHHbQKiCBWoG7AcspOiB\n5oIJdR2BOX6FExER/50B/Af7pv+75H2Dk5cCjyYfXwQc72s6EREREREJl4OAycCHwDvYVNWi1AVe\nAlYAy7EuKL9kmhFszcYCYIIPudJlkrMxMA1YBiwFbvAtnY05rQRWAcVtQP5I8vFFQFufcqUrLecl\nWL7FwCygtX/RfiaT3yfAicCPwLl+hCpCJjlj2PtmKTb26LfSMtYD3sK6xpdii29dGI2Nwy4p4Zgg\nvIcq5M/ArcnrtwH3FXPcM8DlyetVgAM8zlVYphkBbgKeB173OlQRMsnZADgueb0W1uWXPv7jhcpY\n12EzoCqljzt1wM24UyY5O7H37+/XBDdnwXH/AiYC5/kVLu31S8tZF/uS0ih5u55f4ZIyyZgP3Ju8\nXg/YgrczOotzMvZBX1xRCMJ7qMJWsncRW4Pk7XQHAC7PXptJRrA/6ilAN9y0FDLNWdh44DTPEu3V\nCfumVeD25KWwEcCFhW4X/vf4JZOchR2IrbnxW6Y5/xu4FngKN0Uhk5zXAn/wLdG+Msk4GHgseb05\n1hp3pRnFF4UyvYeCuiFeffZOS91I0f+Aw4HN2B/2B8CTgJ+zrTPJCPAgMBSbjutCpjkLNMO+dcz1\nMFOBohYvpq8pL26Bo58yyVnYFez9ZuanTH+ffdm7pYyLdUGZ5DwS6/qcBswHBvgTLSWTjE8CLYHP\nsW6ZG/2JVmZleg+5POPrZOyba7o7024nKPoPtwo2W+l64N/Yxnu3A3cHKGNvYBPWLxrLYq50Fc1Z\noBY2RnMjsCM70UoUlgWOZXm9bliXpouzbWSSs+B9ksB+ry62z88kZ1Xs/X0a9mVvNtbtscrDXIVl\nkvEOrFsphq2jmgy0AbZ7F6vcMn4PuSwK3Ut4bCP2IfcFcCj2wZpuXfLy7+Ttlyi5SV8eFc3YGdvf\nqRdQHagDPAsMzG7MCucEexO+DIzBuo/8sB4b5C7QmH27XdKPaZS8z0+Z5AQbXH4SG1PY6kOudJnk\nbMfexaL1sGnju/F3vCuTnJ8BXwLfJi/vYh+4fhWFTDJ2Bu5JXl8DfAwchbVsgiQI76EK+zN7R/tv\np/hB3HeBXyav51P01txeyTRjga64GVPIJGceVqwe9CtUUlgWOGaSswk2MOnnDLh0meQs7CnczD7K\nJOfR2FhcZaylsATbbdkvmWR8ABiWvF4fKxquNtpvRmYDzaFdJHwQ9geRPo3yMOCNQse1wVoKi4BX\n8Hf2UaYZC3TFzeyjTHJ2wcY8FmJdXQvYuz2J18KywLG0nP+HzT4p+P3N8ztgUia/zwKuigJklvMW\nbAbSEvydJl2gtIz1sC96i7CMF/sdMOlFbFzjB6yFdTnBfA+JiIiIiIiIiIiIiIiIiIiIiIiIiIiI\nSOn82A5ExDdB3RBPJCxcbCgn4hkVBRERSVFREBGRFBUFERFJUVEQEZEUFQUREUlRURCpGM0+EhER\nERERERG73x3NAAAAHElEQVQRERERERERERERERERERERERERESnJ/wMoDXKsUv/RfwAAAABJRU5E\nrkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a06da50>"
       ]
      }
     ],
     "prompt_number": 50
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Moran Scatter Plot"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "w.transform = \"R\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "wy = ps.lag_spatial(w, y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scatter(y,wy)\n",
      "vlines(y.mean(),wy.min(),wy.max())\n",
      "hlines(wy.mean(),y.min(),y.max())\n",
      "xlabel('y')\n",
      "ylabel('wy')\n",
      "title('Moran Scatter Plot')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 53,
       "text": [
        "<matplotlib.text.Text at 0x10a7a2cd0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEZCAYAAABxbJkKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FPX9x/HXZrM5NuFGCRUwCsjRcomIUI4I1hODFkVt\nFc9CWw/qAR78NNii1VqvWq0tVsAWLFiwCoqCYhCwgIBFLAiCXJVDMMEkLCQhu78/vhuzhCQkYXdn\nJvN+Ph77YGd2d/ad3WU+M9/vd2ZARERERERERERERERERERERERERERERBwgCJxudQixvwSrA0iD\ntQ0oBlpUmv8JZgXVLt6BqpAEPAnsBAqBrcDTUVjuNmBIxHQm5m+O5v+38mUWUpH93nos5wZgSdRS\nieOoCEishIAvgWsi5nUDUsOP1Yf3RENVcj9wJtAHaARkAaujsNwQ4KliflXzaiOxhseaYLJfAzwE\nnF/P9xARiaqtwARgZcS83wMPcPSeQBPgFeBrzBb0BCpWljcAy4CngP3ArzFNHIvC0/uAv4eXUW4b\ncDewFjgA/ANIribjXGBsDX9DW2BOONt+4Lnw/PY1ZPgbUAYEMFvo44DtHL3V3jf83JuA9UAe8A5H\n7x0FgV8CXwBbqsiWybF7FyuBuyJeX94cVN1n3AU4BBwJ58qr4bMQEamTrcBQ4HOgM2YrfidmRRdZ\nBF4BXgfSgFOBjZiVI5giUArcilnZpWBWwEMBH9ASWMzRTThbgeVABtAMs5IdU03GCZgV9C8weymR\nW+peTCF5ErP3kgz8MPxYbTJENgedyrEr7OGYFXyn8PwJmIJXLgi8CzSl6iKWGX6ON5z7h8BB4NyI\n15cXgZo+4+tRc5CIxEB5EZgAPApciFmpeakoAl5Mv0HniNeNBj4I378Bs5KuyWXAmkrv+5OI6ceB\nP1Xz2gTM1vZS4DDwFTAq/Fg/zJZzbZpMq8pwvD6B+VSsiMuzHMTsfRB+flYN71m+zHzMFvx64LaI\nx8uLQG0+YxUBF6uprVHkRIUwzSNLgNMwW6SRW9stMVvTkSv6HcApEdM7Ky2zFfAsMADTFp7Asc0Y\neyLuHwK+V02+IPBC+JYM3Ay8jGlWaUtFM05ltclwPKeGl/FkpfmnUPE3V/7bq9KimozlavMZi4up\nY1hibQemg/giTPt6pP2Y5p7MiHntgP9FTFfuRH4U0+b+A0xb93XU/DuubSd0MaYY5GPayndQsbdS\n2fEyVH7PqjLswGyRN4u4pWGasuqavSbH+4yj8R7iYCoCEg83Y5pHDlWaXwbMAh4B0jFbx3diOlqr\nk45pNinAbM2OO8571zQiZywwGNPmn4hpH0/HDGNdCewGHgP8mP6I/rXMsBfTb1BuH2ZrPXLei5hO\n8q7h6SbAlcf5W+rjeJ/xXqANZm9BRCRqKreLl0vErJjKO4abYpqMvsZsHf8fFSvu64EPK72+K7AK\nM5plDWY0zI4a3jcH0wxVlZ+Fl3UAswewHLg44vG2mA7V8lFAz9QyQzam+SWfitE6D4f/xnzg7PC8\na4FPgW/Dr38pYhll1HywV2b4OdVtyEW+vqbP2AfMA74JPy4SNW0xnU//BT4D7gjPn4jZFf0kfLvQ\ninAiIhJbGUDP8P10zLC0Lpgts7uqe5GIiMRPLEcH7aFilEYRsIGKEQn1PXJSREQcKBPTRpqO2RPY\nhjkQ56+Y9koREWmg0jGdaJeFp0/G7Al4gEmYQiAiIhaIdbNM+ciD+VSMrIiUiTl/S7fIme3btw9t\n2VLV6VJERKQGW4AOdXlBLI8T8GC28tdzdAFoHXH/cmBd5Rdu2bKFUCjk2FtOTo7lGZTf+hzK77yb\nk7OHQiE4+liUWollx/APqRgH/Ul43gOYU972xBypuJXqT+4lIiIxFssisJSq9zTmx/A9RUSkDnTa\niBjIysqyOsIJUX5rKb91nJy9vuw6Xj8Ubt8SEZFa8ng8UMf1uvYERERcTEVARMTFVARERFxMRUBE\nxMVUBEREXExFQETExVQERERcTEVARMTFVARERFxMRUBExMVUBEREXExFQETExVQERERcTEVARMTF\nVARERFxMRUBExMVUBEREXExFQETExVQERERcTEVARMTFVARERFxMRUBExMVUBEREXCzR6gAiIvUR\nDAZZsGABX3/9Nf369aNjx45WR3Ikj9UBqhEKhUJWZxARmyorK2PYsJEsXboF+D7B4AJmzpzCsGHD\nrI5mKY/HA3Vcr2tPQEQc580332Tp0h0UFX0M+IBljBo1kry8r6yO5jjqExCpwcSJE62OIFXYvXs3\nZWVnYgoAQB++/XYvwWDwmOfu2LGDKVOmMHPmTA4dOhTXnE6g5iCRGng8HvRbtJ81a9YwcOAlBALv\nA13wen9Nt27v8cknS4563scff8yQIZcQCv0Ij2cv3/teHqtWLaZRo0bWBI+x+jQHaU9ARBznzDPP\n5PnnHycl5Ry83lQ6d36bN9+ccczzRo++m6KiJzl4cDpFRQvZvr0Tzz//ggWJ7UtFQEQc6YYbRnHw\n4AEKCvL47LMVtG3b9pjn7N27F+gdnvJQXNybnTv3xDWn3akIiIhjJSQk4Pf7q3383HMHkZz8CHAY\n2InfP5nzzhsUt3xOoCIgIg3Wn//8NFlZh/F6G5OU1JkHHriZyy+/3OpYtqKOYZEaqGPYGmvWrGHx\n4sWcdNJJjBw5kqSkpBNa3pEjR/B6veUdpw1WfTqG7fqJqAiILagIxN/MmbO46aY7OHLkSny+z+ja\nFZYuffeEC4Eb2G10UFvgA+C/wGfAHeH5zYGFwCZgAdA0hhlExGHGjBlLIDCPkpLnOHjwfdavDzJ7\n9myrYzVYsSwCpcCdwPeBc4BbgS7AfZgicAbwfnhaRIRQKERR0TeY1QZAAmVlXdm/f7+VsRq0WBaB\nPcB/wveLgA3AKUA2MC08fxpwWQwziIiDeDwe+vcfis93L2a18W8SEuYwePBgq6M1WPEaHZQJ9AJW\nAK2AveH5e8PTIiIAzJnzN/r330Ji4kk0bz6Sv//9Rbp37251rAYrHieQSwdmA2OBwkqPhcK3Y0Se\nsyUrK4usrKzYpBMRW2nZsiW5uW9ZHcMRcnNzyc3NPaFlxHp0kA+YB8wHngnP+xzIwjQXtcZ0Hneu\n9DqNDhJb0OggcRK7jQ7yAH8F1lNRAADeBK4P378e+FcMM4iISA1iuScwAPgQ+JSKJp/7gZXALKAd\nsA0YCRyo9FrtCYgtaE8guvbv38+UKVMoLDxIdvYwzjrrLKsjNSg6WEwkylQEomffvn1069aX/PzB\nlJZ+j9TUl3jttSlcfPHFVkdrMOzWHCQi8p0XX/wzeXnnUVIyhVDoEQKBqYwd+6DVsVxPRUBE4iI/\nv4DS0lMj5pxKYWGBZXnEUBEQkbgYPvwSUlOfx3QVbiE19S5GjMi2OpbrqU9ApAbqE4iuf/xjJuPG\nPUwgUMQVV1zOc889oRPDRZE6hkWiTEUANm7cyEsvTSUUCjFq1E909K6NqQiIRJnbi8C6devo128I\ngcDPCIV8+P0vsGjRPPr27Wt1NKmCioBIlLm9CIwceQP//Gc3QqG7w3P+wo9+9C4LFujUznakIaIi\nElWFhQFCochzPGZQWHjQsjwSfSoCIlKt668fgd+fAywFluP338cNN1xhdSyJonicRVREHOrqq6+i\noKCQxx67jWAwyB13/JzRo2+2OpZEkfoERGrQkPsEDh8+zJNPPsO6dV/Qp083xo69jcTE+m0XBoNB\nXnjhRZYuXU3nzpmMG3cXaWlpUU4sx6OOYZEoa6hFoKysjMGDL2b16hQOH74Ev38WQ4e25I03Xi1f\nkdTJqFGjmT17PYHAKFJS3qdz569YufIDfD5fDNJLdVQERKKsoRaB1atXk5V1DUVF6zGtwodJSWnH\nxo2raNeuXZ2WlZeXR0bGqZSW7gIaAUHS03szd+7TuhhUnGl0kIjUSnFxMQkJjajoFkwmISGV4uLi\nOi+rpKQErzcZ8IfnJJCQ0KRey5L4UxEQcaFevXrRtGmAxMQHgY9JShrLaadlcPrpp9d5Wa1ataJH\njx4kJ48BVuL1/paUlO3069cv6rkl+lQERFwoNTWVjz56jwsu2MTpp48hO7uAxYvfxuv11nlZHo+H\nd9+dw4gRHtq3/zlDhnzMihUf0Lhx4xgkl2hTn4BIDRpqn4A0TOoTEBGROlEREBFxMRUBkQZozpzX\nueiikYwYMYrVq1dbHUdsTH0CIjVwYp/A3/42nZ///AECgUlAPmlpv+Gjj97XdQBcQAeLiUSZE4tA\n16792LDh18CPwnMmMWbMPl588VkrY0kcqGNYRMJFK3Kop49gMGhVHLE5FQGRBubOO2/B7x8DvA78\nFb//94wefb3VscSmdCppkQZm9OhbSE5OZvLkl0hNTSYn53XOOussq2OJTalPQKQGTuwTEPdSn4CI\niNSJioCIiIupCIg4yJYtW7j00mvo2XMw48b9n07XLCdMfQIiNbBTn8D+/fvp1KknBw7cRjB4Nqmp\nT3LJJS147bVXrI4mNqGDxUSizE5FYMaMGYwZ8xpFRa+H5xzE621OIFBIUlKSpdnEHtQxLNKAmXP9\nH46YY5qC6nNNYJFyKgIiDnHRRRfRtOmX+Hxjgen4/cP42c9+oYu5ywmx6yaEmoPEFuzUHASwb98+\ncnIeZfv23Zx//gBuv/2XJCRoW04M9QmIRFm8ikBpaSkPPPAw//rXO7Ro0Yxnn51E3759Y/6+0rDY\nsU/gZWAvsC5i3kTgf8An4duFMc4gYnu33no3L7ywnM2b/8iKFdcxdOgwNm3aZHUscYFYF4EpHLuS\nDwFPAb3Ct3dinEHE9mbMmEEgMA04BxhFSclVzJ071+pY4gKxLgJLgPwq5tu1GUrEEj5fMnDgu2mv\n9wDJycnWBRLXsKpH6XZgLfBXoKlFGURs48EH78Xvvwx4gcTEX9G48VKuvvpqq2OJC8RjizwTmAt0\nC0+fDOwL3/8N0Bq4udJr1DEsthDP0UGzZ89mzpx3aNWqOePH30lGRkZc3lcaDruODsrk6CJQm8dC\nOTk5301kZWWRlZUVm3QiNbDbEFGRSLm5ueTm5n43/fDDD4MDikBrYHf4/p1AH+AnlV6jPQGxBRUB\ncZL67AnE+spirwKDgZbATiAHyAJ6YkYJbQXGxDiDiIhUw66jdLQnILagPQFxEjseLCYiIjamIiAi\n4mIqAiIiLqYiICLiYioCIiIupiIgEmMzZ86ideuONGp0MtdcczOBQMDqSCLf0RBRkRqc6BDRf//7\n35x33o8JBP4JZJKSMpYf/7gZ06dPjl5IkTA7Hiwm4mrvvPMuhw7dBPwQgMOHn+Ttt/tbG0okgpqD\nRGKoWbOmJCdviZizhUaNdOJcsQ81B4nU4ESbg7799lt69uzP3r3fp6TkNJKTp/Lqq5PJzs6OYkoR\nw65nEa0PFQGxhWicNqKgoIBXXnmFb7/9lgsvvJDevXtHKZ3I0VQERKKspiKwc+dOLrvsWtauXc5J\nJ7Vl+vS/MGTIkDgnFKmgIiASZdUVgVAoRJcuZ7F582WUld0NLCEt7VrWr19Nu3bt4h9UhNidQO4p\n4Pv1CSTSUOXn57N16xeUlf0f4AcuwOsdyIoVK6yOJlIntSkCG4C/ACuBnwNNYppIxAHS09MJhcqA\nbeE5JQSDm2jRooWFqUTqrjZFYDJmkPMozFXC1gEzgHNjF0vE3pKSknjiid/h9w8iKekO0tIGMnBg\nZ10GVRyntm1HXuBS4EagDTALGAAEgKtikEt9AmILxxsdtGzZMlauXEmbNm0YMWIECQk69EasE6uO\n4acxBWAR8BKmWajcRqBTXd6wllQExBZ0ZTFxklgVgd3Ab4GFmP6BSE2BA3V5w1pSERBbUBEQJ4nV\n6KCfYlb2f8BcGH428KvwY7EoACK2UFpaanUEkZirTRFYBDwKPIjpJO4D/CKWoUSsFAwGue22u/H7\nGwEwfPg1HD582OJUIrFRmyLwPrAU0wG8ETiL2PQDiNjCn/70F6ZMWcqRI7sAWLjwEOPHP2RxKpHY\nqE0R+BQoBX4AdA//mxrLUCJWeu+9ZQQCvwCaA3Do0F28//4ya0OJ2EAj4HZgO1Ac4/cKOVlOTk4I\n0E033Wx0y8nJsXrVEHPhv7VOatOLfDswEOiN6RheEr4tquub1UH47xGJv2+++YbevQfyzTenUFT0\nHs2bn8KKFbl06NDB6mgiNYrVENFxwIfAGkyzUDyoCIilioqKeOedd7jyyivZt28fLVu2tDqSyHHp\nLKIiUabjBMRJYnWcgIiINFAqAuI6xcXF7Nixg+LiWI9vELE/FQFxlYULF9KyZRu6dOlHy5ansGDB\nAqsjiVhKfQLiGvn5+bRtewYHD84GBgEfkpY2gp07N9GsWbMqX6M+AXES9QmI1GDz5s14vW0xBQBg\nEF5vW7744gsrY4lYSkVAXKNNmzaUlGyj4mpg2ygp2Ubbtm2tCyViMRUBcY3WrVvz+OOTSE3tS+PG\nF5Oa2pfHH59E69atrY4mYhn1CUiDEQwG+fLLL/H5fLRr1668ffQYGzduZNOmTZxxxhl06lTzuRDV\nJyBOooPFxLUOHDjAkCGXsnHjNoLBEs49dyBvvPEqPp/vhJarIiBOYseO4ZeBvZiL05drjrlK2SZg\nAeaCNSInZOzY+/jvf7sQCGzn8OEd5OYe5IknnrI6lojtxboITAEurDTvPkwROANzrYL7YpxBXGDV\nqk8pKbkO85NO5tChq1mx4lOrY4nYXqyLwBIgv9K8bGBa+P404LIYZxAX6Nq1Iz7fXMyZdMtISXmb\nH/ygo9WxRGwvHn0CmcBcoFt4Oh8oPzLHA+RFTJdTn4DUyZ49e+jf/zz2708iFArQuXMGixe/jd/v\nP6Hlqk9AnKQ+fQKJsYlSa/W6CIJIZRkZGaxfv4o1a9bg8/no1asXiYlW/7xF7M+K/yV7gQxgD9Aa\n+LqqJ02cOPG7+1lZWWRlZcUhmjhZSkoK/fv3tzqGSNzk5uaSm5t7Qsuwojnod8A3wOOYTuGmHNs5\nrOYgsQU1B4mT2PE4gVeBwUBLzB7AQ8AbwCygHeb4/ZHAgUqvUxEQW1ARECexYxGoLxUBsQUVAXES\nOx4sJiIiNqYiICLiYioCIiIupiIgIuJiKgLiCJ999hnz58/nq6++sjqKSIOiIiC2d889Ezj77PO5\n5pqnOOOMnsybN8/qSCINhoaIiq19/PHHZGVdQSDwCeYs5Mvx+y+moGAfXq835u+vIaLiJBoiKg3O\nl19+idfbB1MAAM6hrCxEfn7lk9OKSH2oCIitde/enSNHlmCuQQQwi8aNG9OiRQsrY4k0GCoCYmtd\nunTh2Wd/S3LyWaSltaNFi7uZP39OtdcPFpG6sev/JPUJuMSePXv46KOPaNy4Meeee2617fxFRUXs\n27ePU045haSkpLjlU5+AOInOHSSOsmrVKoYMuQSPpy/B4A569sxg0aK5J3xx+GhSERAnUcewOMp1\n1/2SwsKnKSh4k6Ki1axZU8bUqVOtjiXiKioCYpldu3YCg8JTXgKB/uzYsdPKSCKuoyIglunTpy8+\n37NAENhNWto/6Nv3bKtjibiKioBYZvr0P9O160ckJTXD52vPPfeMYtiwYVbHEnEVdQyLpUKhEHl5\neaSlpZGSkmJ1nGOoY1icRKODRKJMRUCcRKODRESkThKtDiANU2lpKW+88Qb5+fkMGjSITp06WR1J\nRKqg5iCJupKSEgYMuIANG0oJBjsC85gz5+9ccMEFVkerMzUHiZPUpzlIewISdTNmzGD9ei8HD76P\naXF8j5tu+iVffbXpeC8VkThTn4BE3d69eyku7knFz+tM8vL2WhlJRKqhIiBRN3DgQJKSXgXWA6X4\nfBPp12+w1bFEpAoqAhJ1/fv35w9/mITf/0MSEvz06bORWbNetjqWiFRBHcMSM6FQiLKyMhITndv1\npI5hcRIdLCYSZSoC4iQ6WExEROrEufvpUmuhUIh58+axbt06OnbsyBVXXHFCl2f88MMPWbZsGRkZ\nGfz0pz+N65W+RCS61BzkAnfeeR+TJ7/J4cOXkpLyPtnZPZg+/aV6FYLnn3+R8eMfpbj4alJSVtO9\newIffjjf0e3+NVFzkDiJ+gTkGLt37+a007pSXLwFaA4E8Ps7s3z5W3Tr1q1OywoGg6SmNqak5D9A\nByBIeno/pk+fQHZ2dgzSW09FQJxEfQJyjAMHDuDztcQUAAA/Pl8b8vPz67yskpISyspKgdPCcxKA\nDvValojYg4pAA9e+fXsaN4aEhKeBb4CpJCRsp0ePHnVeVkpKCmee2Z/ExHHAfmA+weACBg4cGOXU\nIhIvKgINXFJSEosXz6d79zmkpJxO587Pk5s7nyZNmtRreW+9NYsBA74gNbUDbdrcxdy5Mzn99NOj\nnFpE4kV9AvKdXbt2MXXqNA4dOsyIEZfTs2dPqyNZTn0C4iTqGJZ627lzJz179qOw8FKOHGlGaupk\n3nrrNbKysqyOZikVAXESpxWBbUABUAaUAmdHPKYiEGd33XUvf/hDGWVlvw/PmUWvXi+wZk2ulbEs\npyIgTuK06wmEgCwgz8IMEpaXV0BZWZeIOadSUFBgWR4RiQ+rO4bt2hzlOlddlY3f/yTwEfA5fv89\nXHXVcKtjiUiMWbkS/hL4FtMc9GdgcsRjag6ywJQp03jwwccoKSnmuuuu4ne/m4TX67U6lqXUHCRO\n4rQ+gdbAbuAkYCFwO7Ak/FgoJyfnuydmZWW5voOyOmvXrmXq1OkkJHi45ZYb6NKly/FfJLWmIiB2\nlpubS25u7nfTDz/8MDioCETKAYqAJ8PT2hOoheXLlzN06KUEArfi8ZTh97/IsmXv1etAMKmaioA4\niZNOG+EHGoXvpwHnA+ssyuJYDz30BIHAo8BEQqHfcPDg/Uya9LTVsUTEQawaHdQKeD0iw3RggUVZ\nHKuoKACcHDGnVXieiEjtWFUEtgI6HPUE3XTTlaxdex+BwEnAEfz+h7jxxsesjiUiDtIwTwLvEjff\nfCOBwCGeeeaXeDwe7r9/AiNHXml1LBFxELt0DFemjmGxBXUMi5M4qWNYRERsQEVARMTFVARsqrCw\nkPHjJzB8+E954omnOHLkiNWRRKQBUp+ADZWUlHDmmQPZvLkzxcVD8funMWxYW2bOnGp1NNdRn4A4\nidNOG1ETVxeBDz74gOHD76GwcBXmKwqQlJTBrl1badGihdXxXEVFQJxEHcMNRGlpKR6Pn4rvMgmP\nx6cmIRGJOhUBG+rfvz/p6Xvweh8CPiQ5+Ub69DmLk08++bivFRGpCxUBG0pPT2fFig8YNmwLXbve\nx7XXNmX+/Nnlu3oiIlFj17WKI/sEtm/fzoQJk9i1ax+XXjqUsWNvJSFBddbJ1CcgTuK0y0s2KF9/\n/TW9ew8gP/8mgsELWLHiCXbs+Iqnn9a5fETEvrQnECWTJ0/mV79aRCDwanjObpKTz+DQoQI14ziY\n9gTESTQ6yEJmRRF5KUZvlSuPUCjESy+9zKBBlzJs2NWsWbMmbhlFRCpTEYiS7OxskpIWkZDwW2Ae\nfv8Ibrll9DF7Ac8880fGjn2CJUtu5K23BjBo0AVs2LDBmtAi4np2badwXHMQwObNm7n33l+ze/c+\nhg0bwr333nXMhdrbtv0+//vfFOBsADyeBxg/3sNjjz1iQWI5HjUHiZOoY9hiHTp0YPbsV2p8jvmS\nIlcqQTwefQ12lZOTY3UEkZjSnkCcPffcC9x337MEAr8BdpGW9girVy+lU6dOVkcTEYfTuYMcYtq0\nv/HKK3No0iSdnJx76NGjh9WRRKQBUBEQEXExDREVEZE6UREQEXExFQERERdTERARcTEVARERF1MR\nEBFxMRUBEREXUxEQEXExFQERERdTERARcTEVARERF1MREBFxMRUBEREXUxEQEXExFQERERezqghc\nCHwOfAHca1EGERHXs6IIeIE/YgpBV+AaoIsFOWImNzfX6ggnRPmtpfzWcXL2+rKiCJwNbAa2AaXA\nP4DhFuSIGaf/kJTfWspvHSdnry8risApwM6I6f+F54mISJxZUQR08WAREZuw4kLz5wATMX0CAPcD\nQeDxiOdsBtrHN5aIiONtATpYHeJ4EjFBM4Ek4D80sI5hERGp2UXARswW//0WZxERERERETuZiBk5\n9En4dmGNz7YPpx8Itw34FPOZr7Q2ynG9DOwF1kXMaw4sBDYBC4CmFuSqraryT8Q5v/u2wAfAf4HP\ngDvC853yHVSXfyLO+A5SgBWYJvX1wG/D853y+R9XDnCX1SHqyItp4soEfDizv2Mr5kfkBAOBXhy9\nEv0dMD58/17gsXiHqoOq8jvpd58B9AzfT8c08XbBOd9Bdfmd9B34w/8mAsuBAdTx87f7uYOsGL10\nIhrKgXBO+dyXAPmV5mUD08L3pwGXxTVR3VSVH5zz+e/BbOgAFAEbMMf8OOU7qC4/OOc7CIT/TcJs\nhOZTx8/f7kXgdmAt8FecsUvTEA6ECwHvAauAn1mcpT5aYZpYCP/bysIs9eW03z2Yvd9emOYJJ34H\nmZj8y8PTTvkOEjCFbC8VTVuO+vwXYnaFK9+ygZMx1dgDTMJ8GXY3ApgcMX0t8JxFWeqrdfjfkzA/\nroEWZqmNTI5uTqm8ZZ0Xvyj1ksnR+Z34u08HVlOxxem07yAds9FTnt+J30ETTAE7F+d9/rWSydH/\nUezqHOCdiOn7cWbncLkc4G6rQxxHJkf/Nj7HtPWCKWifxztQHWVS/W+7psfswge8C/wqYp6TvoOq\n8kfKxP7fQbkHgXuo4+dv5+ag1hH3L8cZX8QqoCMVB8JdBbxpZaA68gONwvfTgPNxxuce6U3g+vD9\n64F/WZilPpz0u/dgtpLXA89EzHfKd1Bdfqd8By2paKpKBX6EGc3klM//uF7BDFVci/kjbN2uFcHJ\nB8KdhmkC+g9myJzd878K7AJKMH0xN2JGNr2HM4bHVc5/E8763Q/AnPLlPxw9nNIp30FV+S/COd9B\nN2ANJv+nwLjwfKd8/iIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiEi1eqwOIOMTDQA/MWTIBHgF+\nEDEtIiKaobASAAAAj0lEQVQN2KmYM2WCOefWZqCZdXFEoiPR6gAiDrEd+AZzJaoMzDlbqrogjIij\nqDlIpPaKMde6GABMwewNiIiIS/ioOEOsUy4/KFIjNQeJ1F4psAjTDBSyOIuIiMRZAuac8+2tDiIS\nLXa+spiInXQFvsBcrGOLxVlEREREREREREREREREREREREREREREqvP/6K0OImHctrQAAAAASUVO\nRK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x109ebbd90>"
       ]
      }
     ],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}