Power Law vs Stable Random Number Generator

powerlaw_vs_stabrnd

{
 "metadata": {
  "name": "powerlaw_vs_stabrnd"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Power Law vs Stable Random Number Generator: which one should be used to generate a Truncated Power Law Distribution?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "From a discussion with Vincent Gauthier (author of an implementation of the Truncated Levy Walk model):\n",
      "\n",
      "Stable random number generators are usually used to generate distributions with fat tail (specially Levy distribution).\n",
      "Stable distributions in general remain consistent  (i.e. without any bias)  against the following transformation: scaling and truncation.\n",
      "\n",
      "Quoting wikipedia : \"Stable distributions owe their importance in both theory and practice to the generalization of the central limit theorem to random variables without second (and possibly first) order moments and the accompanying self-similarity of the stable family\",\n",
      "which is exactly the case of Levy distribution, where the second moment of the distribution is infinite (or ill defined in the case of experimental data that follow levy distribution or truncated levy distribution).\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# define a Stable Random Number Generator\n",
      "# Thanks to Vincent Gauthier for this Python implementation\n",
      "def stabrnd(alpha, beta, c, delta, shape):\n",
      "    \n",
      "    import numpy as N\n",
      "    import numpy.random as R\n",
      "\n",
      "    # Error traps\n",
      "    if alpha < .1 or alpha > 2 :\n",
      "        print 'Alpha must be in [.1,2] for function stabrnd.'\n",
      "        x = N.nan * N.zeros(shape)\n",
      "        return x\n",
      "\n",
      "    if N.abs(beta) > 1 :\n",
      "        print 'Beta must be in [-1,1] for function stabrnd.'\n",
      "        x = N.nan * N.zeros(shape)\n",
      "        return x\n",
      "\n",
      "    # Generate exponential w and uniform phi:\n",
      "    w = -N.log(R.rand(*shape))\n",
      "    phi = (R.rand(*shape) - 0.5) * N.pi\n",
      "\n",
      "    # Gaussian case (Box-Muller):\n",
      "    if alpha == 2:\n",
      "        x = (2*N.sqrt(w) * N.sin(phi))\n",
      "        x = delta + c*x\n",
      "        return x\n",
      "\n",
      "    # Symmetrical cases:\n",
      "    if beta == 0:\n",
      "        if alpha == 1:   # Cauchy case\n",
      "            x = N.tan(phi)\n",
      "        else:\n",
      "            x = ((N.cos((1-alpha)*phi) / w) ** (1/alpha - 1) * N.sin(alpha * phi) / N.cos(phi) ** (1/alpha))\n",
      "    # General cases:\n",
      "    else:\n",
      "        cosphi = N.cos(phi)\n",
      "        if N.abs(alpha-1) > 1.e-8:\n",
      "            zeta = beta * N.tan(N.pi*alpha/2)\n",
      "            aphi = alpha * phi\n",
      "            a1phi = (1 - alpha) * phi\n",
      "            x = ((N.sin(aphi) + zeta * N.cos(aphi)) / cosphi)  * ((N.cos(a1phi) + zeta * N.sin(a1phi)) / (w * cosphi)) ** ((1-alpha)/alpha)\n",
      "        else:\n",
      "            bphi = (N.pi/2) + beta * phi\n",
      "            x = (2/N.pi) * (bphi * N.tan(phi) - beta * N.log((N.pi/2) * w * cosphi / bphi))\n",
      "            if alpha != 1:\n",
      "                x = x + beta * N.tan(N.pi * alpha/2)\n",
      "\n",
      "    # Finale\n",
      "    x = delta + c * x\n",
      "    return x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# define a Truncated Power Law Generator\n",
      "P = lambda alpha, lower, uppper, shape: ((upper ** alpha - 1.) * rand(*shape) + 1.) ** (1./alpha)\n",
      "\n",
      "# define a Truncated Levy Generator\n",
      "def TL(alpha, lower, upper, shape):\n",
      "    b = []\n",
      "    n = reduce(lambda x,y:x*y, shape)\n",
      "    while len(b) < n:\n",
      "        b_temp = stabrnd(alpha, 0., 1.0, 0., (n-len(b),))\n",
      "        b_temp = b_temp[b_temp > lower]\n",
      "        b_temp = b_temp[b_temp < upper]\n",
      "        b = np.append(b, b_temp)\n",
      "    b.shape = shape\n",
      "    return b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plot an histogram of both distributions to see the difference:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lower = 1.\n",
      "upper = 100.\n",
      "shape = (1000000,)\n",
      "\n",
      "a = P(-1.6, lower, upper, shape)\n",
      "b = TL(1.6, lower, upper, shape)\n",
      "\n",
      "plt.figure(figsize=(10,8))\n",
      "        \n",
      "plt.subplot(211)\n",
      "plt.hist(a, bins=np.logspace(0., 2., 1000), normed=1)\n",
      "plt.loglog()\n",
      "plt.grid()\n",
      "\n",
      "plt.subplot(212)\n",
      "plt.hist(b, bins=np.logspace(0., 2., 1000), normed=1)\n",
      "plt.loglog()\n",
      "plt.grid()\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAHnCAYAAACGzjiYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3V9oHOe5x/GfHCUlIQXFEAukFSggoT8HgS4cX7ioyBfF\nmCAfGnC0KjixZJMmRjH2IcF2e+G18UVUCuVUok0rn/ikEt26tZOqxccjEqExuSgSIRjEkYpcyMIi\nXxxwMAHTxnitc5Hudne9f2Zn35mdnf1+QETvaGfmXUd69eh5n3nfpu3t7W0BAADAFztq3QEAAIBG\nQvAFAADgI4IvAAAAHxF8AQAA+IjgCwAAwEcEXwAAAD4i+AIAAPARwRcAAICPPA2+vvjiCx07dkyH\nDh3y8jYAAAB1w9Pg64UXXtClS5e8vAUAAEBdqTj4mpiYUGtrqwYGBnKOW5al3t5edXd3a2pqylgH\nAQAAwqTi4Gt8fFyWZeUcS6VSmpyclGVZWl9fVzwe18bGhrFOAgAAhEVzpScMDQ0pkUjkHFtdXVVX\nV5c6OzslSdFoVAsLC2ptbdWPfvQj3bp1S1NTUzp9+vRj12tqanLVcQAAgFrY3t6u6nwjNV9bW1vq\n6OjItCORiLa2trRz50699957un37dsHAK217e5sPQx/nzp2reR/C9B787ItX9zJ53Wqv5fZ8N+cF\n6fsoDB9h+PcM0ntgbDF7LT/HFhOMBF/VZq9isZhs2zbRlYY3PDxc6y5ULUjvwc++eHUvk9et9lpu\nz3dzXn6GHtUJ0s+lW0F6D4wtZq/lx9hi27ZisZir++QzEny1t7crmUxm2slkUpFIxMSlUaEgDS5u\nBek9MECavZafwRfMCsP/gyC9B8YWs9eqt7GladtFDi2RSGhkZERra2uSpIcPH6qnp0dLS0tqa2vT\nnj17FI/H1dfXV74DTU3G0ngAkGbbdqB+2QIIBxNxS8WZr7GxMe3du1ebm5vq6OjQ5cuX1dzcrJmZ\nGe3fv1/9/f0aHR11FHilMe0IwDQCLwAmmZx2dJX5MonMFwAvkPkC4IWaZL68QOYLAAAEGZkvAACA\nGiDzBQAA4AMyXwBQBjVfALwQmswXAABAowhE8MW0IwDTyHoBMIlpRwAAgBpg2hEAiiCbDiCoCL4A\nAAB8FIjgi5ovAKZR8wXAJGq+AAAAaoCaLwAogmw6gKAi+AIAAPAR044AAAAOhWbakYJ7AAAQZBTc\nA0AZ7O0IwAuhyXwBAAA0ikBkvvKRCQMAAEFkIvPVbKgvBd2/f1/Hjx/Xt771LQ0PD+sHP/hBwdd9\n/fXXam5u1o4dO9TU1PRYQEYwBgAAwsLTaccPP/xQr7zyin7961/rT3/6U9HXPfXUU9qx45uubG9v\n53xIygRkhbJkAFAID/EACKqKg6+JiQm1trZqYGAg57hlWert7VV3d7empqYkSVtbW+ro6JAkPfHE\nE646WCoQIxgDAAD1puLga3x8XJZl5RxLpVKanJyUZVlaX19XPB7XxsaGIpGIksmkJOnRo0dVd5as\nGACneNIRQFBVXPM1NDSkRCKRc2x1dVVdXV3q7OyUJEWjUS0sLOjEiROanJzU9evXdfDgwaLXPHLk\nSObclpYWDQ4OZgbO9NRBofb29namvW/fvoIB2PLyctHzadOmTZs2bdq0S7XTn+fHPtVw9bRjIpHQ\nyMiI1tbWJElXr17V4uKiZmdnJUnz8/NaWVnR9PR0+Q54uM4XT1ICjcu27cwgCgCmBOZpx2qn/GKx\nmIaHh40PlPn/OPnTkwRiAADACdu2c7Jh1TASfLW3t2dquyQpmUwqEomYuLRR2cEWgRgQbmS9AASV\nkWnHhw8fqqenR0tLS2pra9OePXsUj8fV19dXvgMB2F6oWOau1v0CAADBUpPthcbGxrR3715tbm6q\no6NDly9fVnNzs2ZmZrR//3719/drdHTUUeCVVuuNtfOfomRZC6D+1XJMARA+ts3G2r4hKwbUJwru\nAXghMAX31fKq4N6EQv/AbIEEBF8QxxMA9ctkwT2ZLwMIxAAAaAw1qfnyQq1rvqpVrlaMejHAf/U8\npgAIHmq+6gyZMcB/1HwB8IKJuIXgy2cEYgAA1C8K7utQqYVeC70GAADUHgX3IUZmDDCDaUcAXghN\n5gv/whZIAACEG5mvOsH0JAAAtReazFcj1Xy5lf8/moVeAQDwDzVfeAyZMSAXNV8AvBCazBeqR2YM\nAID6QOarQZAZAwCgemS+4Fj21kdpPEkJAID/2NuxwWTvQ1lqP0qg3jGmADCJvR3hqVLBF/+vUC8o\nuAfgBfZ2hC/SwRj/nwAAjY6aL/iiUL1Yoa8DAIDyPK35+uKLL3Ts2DEdOnTIy9vAJ9SLoZ5Q8wUg\nqDwNvl544QVdunTJy1ugxgjGAACojKPga2JiQq2trRoYGMg5blmWent71d3drampKU86iPpSLhgD\n/EKxPYCgchR8jY+Py7KsnGOpVEqTk5OyLEvr6+uKx+Pa2NjQ3NycTp06pTt37njSYdSXYivvE4wB\nABqVo4L7oaEhJRKJnGOrq6vq6upSZ2enJCkajWphYUFnzpzR4cOHJUlffvmlfvSjH+nWrVuamprS\n6dOnC17/yJEjmeu0tLRocHAw81drum6Ddv22l5eXc9r79u1TWn4Atry8XPP+0g5HO/15UPpDmzbt\n+mynP8+Pg6rheKmJRCKhkZERra2tSZKuXr2qxcVFzc7OSpLm5+e1srKi6enpyjrAUhP4J/aihEm2\nbWcGUQAwpaZLTZicMorFYhoeHmagbHDZ38xsDI5qMZ4AMMm27ZxsWDV2uD2xvb1dyWQy004mk4pE\nIkY6BWQX7KdRJwYACAPX044PHz5UT0+PlpaW1NbWpj179igej6uvr6+yDjDtiAqw0CucYtoRgBdM\nxC2OMl9jY2Pau3evNjc31dHRocuXL6u5uVkzMzPav3+/+vv7NTo6WnHglcbG2nCq1EKv6f+SIQMA\nmGbbbKwN5MgOtra3tx9rAwBgQmj2dqTgHtXK/0HIX+C12OsAAHDCZME9mS80jFJTkXwPhg81XwC8\nQOYLqECxH5bs7BhBGACgEDJfgAfIjAEAyglN5gsIgvwnJwm4AABecL3IqkksNYEgyV/gNX8zcJax\nqA+MKQBMYqkJoIacBF98T9ceBfcAvGAibiH4AgwpFJTxvQ0A4ULNFxAghX4YqR8DAOSj5gvwSHYm\njJox/zGmADCJmi+gjpEN8wc1XwC8QM0XUKco2geA+kTNF1CnSq22X+hzAjEACI9A1HwB+EZ6jbHs\nD4maMTeo+QIQVIHIfLG3I1BcftarUADG6vwA4C32dgRQNAPGzxMAeIeaL6CBOakbK/U6AEBtUPMF\nhEyxerFGQ80XgKDyPPO1sLCg69ev66uvvtLRo0f1ve99z+tbAvinUgHY9vY2NWIAUAO+1Xzdu3dP\nb7/9ti5dupTbAWq+AN+xDyUAuGMibnE87TgxMaHW1lYNDAzkHLcsS729veru7tbU1FTR8y9evKjJ\nyUn3PQVgTKnlLAAA3nIcfI2Pj8uyrJxjqVRKk5OTsixL6+vrisfj2tjY0NzcnE6dOqU7d+5oe3tb\np0+f1oEDBzQ4OGj8DQCoXrE1xbLbafUSpFHzBSCoHNd8DQ0NKZFI5BxbXV1VV1eXOjs7JUnRaFQL\nCws6c+aMDh8+LEn6+c9/rqWlJX311Vf629/+ph/+8IePXfvIkSOZa7S0tGhwcDCz5ld6AKVNm7b3\n7eXl5Uw7O8gqto5YrftLmzZt2l6305/nx0DVqKjmK5FIaGRkRGtra5Kkq1evanFxUbOzs5Kk+fl5\nraysaHp62nkHqPkCAq/UVkcsbQGgkdR8nS9TUw+scA8EW6GBptCTkulsGQEYgLCxbTsnG1aNHdWc\n3N7ermQymWknk0lFIpGqOwUg2ArVgGUHXaVqxvxiapAEANOqmnZ8+PChenp6tLS0pLa2Nu3Zs0fx\neFx9fX3OO8BfyUAoFZqqLFdDZpJt22TTARjn61ITY2Nj2rt3rzY3N9XR0aHLly+rublZMzMz2r9/\nv/r7+zU6OlpR4JUWi8X4KxUImeylLAoV72e/xovsGIEXAJNs21YsFjNyLTbWBuCrYgX6FO4DqAe+\nZr68ROYLaBz5C7wWKtIvVuBfSXaMMQWASWS+AIRetVsgUfMFwAs1X2rCFJaaAJCv0Mbf2e1Sa49J\n1HwBMMvkUhNkvgAETqmnIMsFXQDgpdBkvgAgW6mBrdTXik1VermkBQBUioJ7AHUtf5HX9Ed6n0oC\nLwAmUHAPAHkK1YM53YeSAA2AU0w7AsA/OVmqotDX8rdKIgAD4LVATDsCgGn5pQzFVtPPnqqUcqcx\na7UvJYBwC0Tmi6UmAPil0Ir66WNOsl5MUQKNiaUmAMAnxWrJCMKAxkTNFwB4JD/Qyg+2SmXQAKAU\nar4AhFK10wPZdWCl6sLyX1vo6wCQjcwXAFTI6cr7hTYNB4BABF8U3AMwze/xpFDAxdQkEB4U3ANA\nnShXmF9sSyQAwWQibqHmC0AoBWXLslL1YPk1Y/mZMurGgHAKxLQjADQCJ8tTZD9dSQYMCCdPpx3/\n+te/6j//8z919+5d7d+/X0ePHn28AwwwABpQfiDGumFAfTARt/hS8/Xo0SNFo1H9/ve/f7wDBF8A\n8JhiwVla/lRm/jEA3vCt5mtiYkKtra0aGBjIOW5Zlnp7e9Xd3a2pqamC5/75z3/WSy+9pGg0WlVH\nAaASQan5ciu/Viz/80ILv1IrBtQHR8HX+Pi4LMvKOZZKpTQ5OSnLsrS+vq54PK6NjQ3Nzc3p1KlT\nunPnjiRpZGREN27c0AcffGC+9wAQctnF+YUWfk0r9TUAweKo4H5oaEiJRCLn2Orqqrq6utTZ2SlJ\nikajWlhY0JkzZ3T48GFJ0s2bN/Xhhx/qH//4h/bt21f0+keOHMlcp6WlRYODg5k1etJ/vdKmTZt2\nJe3h4eFA9cdte3l5+bGvp8fT5eXlnNcXynYtLy8XfX0Q3h9t2kFvpz/Pj4Oq4bjmK5FIaGRkRGtr\na5Kkq1evanFxUbOzs5Kk+fl5raysaHp6urIOUPMFAK45rfeiLgwwo6brfJmsJ4jFYjkRJgBUqxHG\nlEoCqkL7Uhb7APA427YVi8WMXMv1Ol/t7e1KJpOZdjKZVCQSMdIpAEB5TrNd6ddW8tc6mTLAO66n\nHR8+fKienh4tLS2pra1Ne/bsUTweV19fX2UdYNoRAHxTLKhiKQvAGd+mHcfGxrR3715tbm6qo6ND\nly9fVnNzs2ZmZrR//3719/drdHS04sArjWlHAPBHsQxY9hZH+V/n6UnA7LQjG2sDCCXbtjNPLaEy\n2ZmuUou9UuSPRmQibgnE3o6xWEzD/3w0HABQW4UWdy0VSBULyAi6ECa2bRubpSPzBQBwpJIsGHtX\nIqzIfAEAfFNJEX6x1zrJoBGgIYjIfAFAGdR8Va/SYKjUE5NAWIQm8wUACJ5Kf8FUsu5YtfcC6lkg\ngi+mHQGYxngSPNUs8sqUJGqNaUcAQF1wO3XJIq8Iqpru7QgAQcbCzcHgZIHW7D0liz1JmS7WL7YH\nJftSop4QfAEAaqrQ1GJ2sJX/uvzP86/BJuEIOmq+AIQS40l9KbblUTqISn+9miUqmL5ENaj5AgA0\ntEoCKdYWg0nUfAFAEdR8hVv+NGO2YtOVhaYh82vSmK6EHwIx7QgAQKWcrK6ff6ySlfmdvB5wg2lH\nAEAouZ1uZLkLlBKaFe4puAcAmFbsF2T2tGKxqchCry+3NyXCjYJ7ACiDvR0hVZexItBCIaHJfAEA\n4IVqfkmWq/8qljUjU4ZyyHwBAFABsmmNrS6Wmrh//75efPFFXb9+3etbAQDgSv4SE6WWnHCyZVIx\n1ZyL8PA8+PrJT36i0dFRr28DADlY5wuVyA+KvAiSyq0hVskaY6xHVt8cBV8TExNqbW3VwMBAznHL\nstTb26vu7m5NTU09dt7HH3+s/v5+Pf/882Z6CwBAQJULiIrtWZn/9UquhfrkqObr008/1bPPPqtX\nX31Va2trkqRUKqWenh598sknam9v14svvqh4PK7PPvtMn3/+ud555x394he/0P3797W+vq6nn35a\nH330UcGViPkGAgDUK7frglVyHrViweHb045DQ0NKJBI5x1ZXV9XV1aXOzk5JUjQa1cLCgs6cOaPD\nhw9Lki5evChJ+uCDD/T8888XjeKPHDmSuU5LS4sGBwczj4inpw5o06ZNmzbtILbTv4ht29a+ffsy\nx+x/rgtV7Pzl5WVls21by8vLBV/v5Hq0vWmnP8+Pg6rh+GnHRCKhkZGRTObr6tWrWlxc1OzsrCRp\nfn5eKysrmp6erqwDZL4AeCD7lxRgmpuNvbNf7/R8Ml7BU9N1vkwW+rHCPQCgnlTyy9fJXpPFXktw\nFhzpzKMJroOv9vZ2JZPJTDuZTCoSiRjpFABUiz/mUA+qXYi1WL1YoWwbgsP1tOPDhw/V09OjpaUl\ntbW1ac+ePYrH4+rr66usA0w7AgBCqNqsVKGsmJ+bfpNVK8y3RVbHxsa0d+9ebW5uqqOjQ5cvX1Zz\nc7NmZma0f/9+9ff3a3R0tOLAKy0WixlL5QGAJMYU1Jzb5SDSy0ykz8++TqFr5h/LXqaimvXAWM4i\nl23bisViRq7F9kIAQomCe4RVqToxt0tcwLnQbKxNwT0A0xhPEDZOgiYTQRfTmIWZLLgn8wUAQJ1z\nG9iUO6/eAyYv1MXG2k5Q8wXANMYUNJJK67MKbSRu8rpu6syCvl8lNV8AUAY1X4BZZMG+YSJuIfgC\nACDkKinEL/c6N9cM0+95Cu4BAEBZToKFSgMKL64ZZBTcA0AZTDsCj3OTico/x3Q2y80+mbWMG0KT\n+QIAAN7IDlicBg355+Qv/FrpNUqpdp/MehSI4ItpRwCmMZ4A33CTUapkg+9iQVZYAqU0ph0BAEDN\nuJkqzH59PRfrh2adLwAwjXW+APPcTD9m709Zyflu9pYM+lphaYGYdgQAAMHnNFPldCqy1D6VpdT7\nVCfTjgAAINAqXYMs/7xKNx8vd02edgQAAKFQLBPmdjPwQucFIeETiJov9nYEYBpjChBMpeqysuvD\nnARcTveoLHdfJ9jbEQDKYJFVwD9ePpkYpKceTfWF4AsAADQkNwX/1HwBAIC6VCyL5OSpyUqvW2qa\nsxY8rfmybVtDQ0N68803dfPmTS9vBQA5qPkCgq1YXVf28UKvKVe7Veic/FoyN2uImeRp5mvHjh36\n9re/ra+//lqRSMTLWwEAgAZQ6VpglZ5X7DomgzVHNV8TExO6fv26du3apbW1tcxxy7J08uRJpVIp\nHTt2TKdPn845L50y/L//+z/9x3/8h+bn5x/vADVfAACgxpyuJebb9kLj4+OyLCvnWCqV0uTkpCzL\n0vr6uuLxuDY2NjQ3N6dTp07pzp07mTfS0tKir7/+uqqOAgCAYKn1dj4m758/HZm/lIXJezl+2jGR\nSGhkZCST+frLX/6i8+fPZ4Kyd999V5J05syZzDkfffSRFhcXde/ePR0/flzf/e53H+9AU5Nee+01\ndXZ2SvomUBscHMw8Ip6u26BNmzbtStrpz4PSH9q0afvXTgdJy8vLVV8v/XkikZAkffDBB/4tNZEf\nfF29elWLi4uanZ2VJM3Pz2tlZUXT09OVdYBpRwAesG07M4gCQD63tVw1XWrCZJoxFotpeHiYgRKA\nMYwnAEqpNICybTsnG1aNHW5PbG9vVzKZzLSTySRPNAIAAJThetrx4cOH6unp0dLSktra2rRnzx7F\n43H19fVV1gGmHQF4gGlHAF7w7WnHsbEx7d27V5ubm+ro6NDly5fV3NysmZkZ7d+/X/39/RodHa04\n8EpjY20AAFALhTbnLlRaZdtsrA0AAOA73zJfXiPzBQAAgozMFwCUQc0XAC+Q+QIAAPABmS8AAIAa\nCE3mCwAAoFEEIvhi2hGAaYwpAExi2hEAyqDgHoAXTMQtBF8AAAAOhabmi2lHAAAQZEw7AkAZTDsC\n8EJoMl8AAACNgswXAACAQ2S+AAAADGtqalJTU1PRdrUIvgCEEg/xAHBre3s7J7uV365WIIIvnnYE\nAABBxtOOAAAANUDNFwAAQJ0h+AIQSpQyAAgqT4Ov7e1t/fjHP9aJEyf0m9/8xstbAUCOW7du1boL\nAFCQp8HXH//4R21tbempp55SJBLx8lYAkOPevXu17gIAFOQo+JqYmFBra6sGBgZyjluWpd7eXnV3\nd2tqauqx8zY3N/Wd73xHP/3pT/XLX/7STI9RUhimWoL0Hvzsi1f3Mnndaq/l9vwgfU80qjD8PwjS\ne2BsMXutehtbHAVf4+Pjsiwr51gqldLk5KQsy9L6+rri8bg2NjY0NzenU6dO6c6dO4pEImppafnm\nRjsoL/NDkAYXt4L0HhggzV7LzwEykUi4uhcKC9LPpVtBeg+MLWavVW/Bl+OlJhKJhEZGRrS2tiZJ\n+stf/qLz589ngrJ3331XknTmzJnMOX//+9/11ltv6ZlnnlFfX5/efPPNxztgcMVYAAAAr1W71ESz\n2xO3trbU0dGRaUciEa2srOS85umnn9alS5dKXoc1vgAAQCNxPRdIxgoAAKByroOv9vZ2JZPJTDuZ\nTPJEIwAAQBmug6/du3fr9u3bSiQSevDgga5cuaKDBw+a7BsAAEDoOAq+xsbGtHfvXm1ubqqjo0OX\nL19Wc3OzZmZmtH//fvX392t0dFR9fX1e9xcAAKCu1XxjbQAAgEYSuMW37t+/r9dee02vv/66fvvb\n39a6OwBC4osvvtCxY8d06NChWncFQIgsLCzo9ddfVzQa1ccff+zonMBlvubm5rRz50699NJLikaj\n+t3vflfrLgEIkUOHDukPf/hDrbsBIGTu3bunt99+u+wSW5JPma9KtifKXj/siSee8KN7AOqU263P\nAKAUN2PLxYsXNTk56ej6vgRflWxPFIlEMktYPHr0yI/uAahTlYwtAOBUJWPL9va2Tp8+rQMHDmhw\ncNDR9X0JvoaGhvTcc8/lHFtdXVVXV5c6Ozv15JNPKhqNamFhQS+//LKuXbum48ePs3QFgJIqGVu+\n/PJLvfHGG7p16xbZMAAlVTK2zMzMaGlpSVevXtWvfvUrR9d3vb1QtYptT/TMM8/o/fffr1W3ANS5\nYmPLzp079d5779WwZwDqWbGxZXp6Wm+99VZF16rZ045sTwTAC4wtALxgcmypWfDF9kQAvMDYAsAL\nJseWmgVfbE8EwAuMLQC8YHJs8SX4YnsiAF5gbAHgBa/HlsAtsgoAABBmgdteCAAAIMwIvgAAAHxE\n8AUAAOAjgi8AAAAfEXwBAAD4iOALAADARwRfAAAAPiL4AgAA8BHBFwAAgI8IvgAAAHxE8AUAAOAj\ngi8AAAAfEXwBAAD4iOALAADARwRfAAAAPiL4AgAA8BHBFwAAgI8IvgAAAHzU7OXF79+/r+PHj+tb\n3/qWhoeH9YMf/MDL2wEAAASep5mvDz/8UK+88op+/etf609/+pOXtwIAAKgLFQdfExMTam1t1cDA\nQM5xy7LU29ur7u5uTU1NSZK2trbU0dEhSXriiScMdBcAAKC+VRx8jY+Py7KsnGOpVEqTk5OyLEvr\n6+uKx+Pa2NhQJBJRMpmUJD169MhMjwEAAOpYxcHX0NCQnnvuuZxjq6ur6urqUmdnp5588klFo1Et\nLCzo5Zdf1rVr13T8+HEdPHjQWKcBAADqlZGC++zpRUmKRCJaWVnRM888o/fff7/kuU1NTSa6AAAA\n4Ivt7e2qzjdScF9tAHXu3DktLy9re3ubjyo/zp07V/M+hOk9+NkXr+5l8rrVXsvt+W7OC9L3URg+\nwvDvGaT3wNhi9lp+jC3Ly8s6d+6cibDJTPDV3t6eqe2SpGQyqUgk4vj88+fPa9++fSa60vCGh4dr\n3YWqBek9+NkXr+5l8rrVXsvt+W7OSyQSru6FwoL0c+lWkN4DY4vZa/kxtgwPDysWi7m6T76m7e3t\ninNniURCIyMjWltbkyQ9fPhQPT09WlpaUltbm/bs2aN4PK6+vr7yHWhq0r/927/pf//3fyvu/Pb2\ndibr5uJtAAixI0eO6L//+79r3Q0AIWHbtmzb1vnz56uOOSoOvsbGxnTz5k3dvXtXu3bt0oULFzQ+\nPq4bN27o5MmTSqVSOnr0qM6ePeusA01NFb0JN1OcBGZA47FtO1CZDgDhUGncUvAabjJfJjU1Nenc\nuXMaHh72ZKCsJFgjSAMAAIXUNPNlmokI0lQ/yglCPwE4Q+YLgBdMxC2e7u3oVCwW8yzz5VS5f8im\npqaiARpBGQAA4ZbOfJlA5qtKBGQAADSO0GS+6lmh/wGlsmTFzgEAAI3ByDpf1YrFYsZSeUFQaHG2\nbOngLPsDgFlhGlMA1J5t27Vd58ukep92NKFc8NXo/z6AGxTcA/BCaJaaILgoLj8w498KAIDaCU3N\nVxCedgyq7P/BPHEJAEBt8LQjJBXOirHdEvANph0BeCE0mS+4Uyorlv05gRgAAMFB5ivkCk1T8u8N\nAIA7ZL5QVvobhKwYAADBwDpfDaLY2mOF1hxj3TGEAWMKAJNY5wvGZRfqU7SPMKDgHoAXWOcLnqFW\nDACAx1HzBc8U2xKp2NcBAIAzgaj5QvAVqxXL/xwICmq+AAQVwRcqlh+EpRGAAQBQXiCmHdleqD4V\nm5rMD8yYokQtMJ4AMInthRBY7D0JAAgzE3EL044wysl6YoAfqPkCEFSBmHZEuBWaniz2NQAAwo5p\nR/iq0GKuaXwfAACCLvDTjl988YWOHTumQ4cOeXkb1JHsqchSy1cAABBWngZfL7zwgi5duuTlLRAS\n1IfBNGq+AASVo+BrYmJCra2tGhgYyDluWZZ6e3vV3d2tqakpTzqIxlIoGwYAQJg4Cr7Gx8dlWVbO\nsVQqpcnJSVmWpfX1dcXjcW1sbGhubk6nTp3SnTt3POkwGkehTBjBGJxinS8AQeUo+BoaGtJzzz2X\nc2x1dVVdXV3q7OzUk08+qWg0qoWFBR0+fFg/+9nP1NbWpi+//FJvvPGGbt26RWYMrrBsBQAgbFwv\nNbG1taV2/iMYAAAVzUlEQVSOjo5MOxKJaGVlJec1O3fu1HvvvVf2WkeOHFFnZ6ckqaWlRYODg5m/\nWtN1G7RpS9Ly8nKmnR2ApQOzWvePdnDa6c+D0h/atGnXZzv9eSKRkCmOl5pIJBIaGRnR2tqaJOna\ntWuyLEuzs7OSpPn5ea2srGh6erqyDjQ16dy5cxoeHs68YaASrBuGQmzbZkwBYIxt27JtW+fPn6/6\nd43rzFd7e7uSyWSmnUwmFYlEXF0rFou57QZQdKkKArHGRuAFwKR0kuj8+fNVX2uH2xN3796t27dv\nK5FI6MGDB7py5YoOHjzo6lqxWCwnvQe4wbphAACv2LZtLFnkaNpxbGxMN2/e1N27d7Vr1y5duHBB\n4+PjunHjhk6ePKlUKqWjR4/q7NmzlXeAFe7hIVbRb1xMOwLwgom4xdG0YzweL3j8wIEDOnDgQFUd\nkL7JfFHzBS9k/4AwLQkAcCtd82UCezuiIeU/JQkAgBO+Zb68RuYLfktv7M20JADACTJfgAcIxMKF\nmi8AXjARt7h+2hEIm/wfJp6WBAB4gWlHIEuxAIwsWP1hPAFgEtOOgI+KZb/4vgWAxsO0I+CD/MVb\n08cQbCzcDCComHYEKpC9en6h4wCAcGLaEQgQ1gwDgMYRmnW+gHpWas0wAjMAQD5qvgADsuvC8qcm\nCbxqg5ovAEEViMwXNV8Im/w9JYt9DQBQH6j5AuoYGTEAqF/UfAF1JjsLRkYMABoTNV+Aj/Jrw/Jr\nxNjSyBxqvgAEFZkvoIaYggSAxkPNFxAwBGQAEFyhqfniaUfgX0qtGwYAqA2edgQaUKFaMH52irNt\nmz/oABgXmswXgNKYigSA8CDzBdQ5picBwD9kvoAGVmxJCoIxAAg21vkC6lSxNcOyg61ixfuNgHW+\nAASV55mvhYUFXb9+XV999ZWOHj2q733ve17fEmh4+Zt7pz8nCwYAtedbzde9e/f09ttv69KlS7kd\n4BcC4AumIwGgeibiFsfTjhMTE2ptbdXAwEDOccuy1Nvbq+7ubk1NTRU9/+LFi5qcnHTfUwBVyR8s\nsrczasRpSQCoFcfB1/j4uCzLyjmWSqU0OTkpy7K0vr6ueDyujY0Nzc3N6dSpU7pz5462t7d1+vRp\nHThwQIODg8bfAADnnOwtGRbUfAEIKsc1X0NDQ0okEjnHVldX1dXVpc7OTklSNBrVwsKCzpw5o8OH\nD0uSfv7zn2tpaUlfffWV/va3v+mHP/yhsc4DMCM/AGNKEgC8U1XB/dbWljo6OjLtSCSilZWVnNec\nOHFCJ06cKHmdI0eOZAK4lpYWDQ4OZlamTv/1Sps2be/by8vLSktnwZaXlzNfz1/stdb9LdUeHh4O\nVH9o06Zdn+305/kJqGpUVHCfSCQ0MjKitbU1SdK1a9dkWZZmZ2clSfPz81pZWdH09LTzDjQ16dy5\ncxr+52AJoLbKTT2SFQPQiGzblm3bOn/+vH8F94W0t7crmUxm2slkUpFIpOLrpDfWBlB7pdYPq6fA\nK/uvVgCo1vDwsGKxmJFrVRV87d69W7dv31YikdCDBw905coVHTx4sOLrxGIxBkqgTuQ/JclTkwAa\ngW3bxoIvx9OOY2Njunnzpu7evatdu3bpwoULGh8f140bN3Ty5EmlUikdPXpUZ8+erawDFPcCdY81\nxAA0ChNxSyA21qbmC6hP+QX4ABBWJmu+AhF8MXAD4RGUgn3btvmDDoBxJuIWz/d2dCJdcM9ACdS/\n7EGJ6UgAYZHOfJlA5guAb4KSFQMAt0KT+QLQGAplxUoNYtSUAQijqpaaMIWlJoDGk143rNASFdnb\nHLkNvBhTAJhUk6UmvMK0IwCpdH2YmwwYBfcAvBCapSYIvgBIuUFWfsBVKDhjWhKA30JT88XTjgCk\n3CAqP+jKD7ayg7HszwnEAHiBpx0BoAymHQF4wUTcEoiCewCoRH6RPvtKAqgngZh2BIBK5P/VWawu\nrNBrAaDWAhF8UfMFwIT8QCt7yQqK8wFUg5ovACijVM0XgRgAt6j5AoA8perBsjNh2W3qxQD4KRDT\njgBgSqF6sPzpx1J/uZIVA+A1gi8AoVdo/bBs2QEXQRcArzHtCCCUnBTG5u8hyRIWAPwQiMwXTzsC\nqIVi+0eW2tqo1NfJmgHhxdOOAOATAi0A2UKztyMABE2xIKtYtszp+QBAzReAUKp2esBJ8X25lfUJ\nvAAUQuYLABwqVReW//VKrgWgsVDzBQAGFMt+sd8kEC6BX+H+r3/9q95880298sor+q//+i8vbwUA\nNZWeZsz+yA+8mIoEIPmU+Xr06JGi0ah+//vfP94BMl8APFBqb8dayl+2AkB98S3zNTExodbWVg0M\nDOQctyxLvb296u7u1tTUVMFz//znP+ull15SNBqtqqMAEAb5mTH2lwQaj6PM16effqpnn31Wr776\nqtbW1iRJqVRKPT09+uSTT9Te3q4XX3xR8Xhcn332mT7//HO98847amtry1zj3//937WwsPB4B8h8\nAUDJ4IsxEggO39b5GhoaUiKRyDm2urqqrq4udXZ2SpKi0agWFhZ05swZHT58WJJ08+ZNffjhh/rH\nP/6hffv2Fb3+kSNHMtdpaWnR4OBgZrog/bg4bdq0aYe5vb29ndNOB2PLy8s5gdny8vJjXw9C/2nT\nDms7/Xl+HFQNxzVfiURCIyMjmczX1atXtbi4qNnZWUnS/Py8VlZWND09XVkHyHwB8IBt25lBFABM\nqekK9ybrE9jbEQDKr/1V6OlJAP6wbTsnG1YN18FXe3u7kslkpp1MJhWJRFxdKxaLue0GABRUj3/M\nOV1Rv5xiQRwLuwLupZNE58+fr/paO9yeuHv3bt2+fVuJREIPHjzQlStXdPDgQVfXisVixqJJAAiT\n/Cchi7ULLWHBOmOAObZtG0sWOar5Ghsb082bN3X37l3t2rVLFy5c0Pj4uG7cuKGTJ08qlUrp6NGj\nOnv2bOUdoOYLgAcaoebLZCaLrBjgjG81X/F4vODxAwcO6MCBA1V1QKLmCwDccDKtmH+sWN1Y/tcJ\nwoBcJmu+2NsRAEKgVNBUalV9gi2gMjV92tEkMl8AUJ1Svwycfo1ADCiOzBcAlNEINV+FZAdQ5Z56\nzMdYDJQXmswXAMCM7F8KxX5BUKAP1FYggi+mHQGYxnhiRrEgq1ABP1OYCDOmHQEAAGrARNziepFV\nAAgyFm6uXvYCroXqxLKPF3sNgMcx7QgAyOF0e6L8+rJ0AMZsBsKIaUcAgKeq2R+yVAaM8R71jqcd\nAQCecPukJIX2QHnUfAEIJWq+zHJa31Vu8+5C16FeDI0mEJkvar4AINicrB9WSKlFX/OXpig2nVNq\neyTAL9R8AQDqVrm1wQp9DQgKar4AAHWn2GKs5Yr7qSdDWFDzBSCUqPmqD8VqxLKzYOUCNGrGUG/I\nfAEAAqOS7Fa1mTAyaagVar4AAKFEcAUvhKbmi6cdAQBOFSrYLzTt6Ef2DI2Dpx0BoAzbtvmDDgRX\nMC40mS8AAEwot5ZYPrdrlgHVIPMFAKgbTgOgSgOlUovBmrwP6h+ZLwBAQ3H6S89NRiv9eToAcxJY\nEXTBDc/X+bp//75efPFFXb9+3etbAUAG63whm5P9KLPXEctvO70W4ITnma+f/OQnGh0d9fo2AAAU\nVM3UYKFzi12HKUg45SjzNTExodbWVg0MDOQctyxLvb296u7u1tTU1GPnffzxx+rv79fzzz9vprcA\n4BBPOiKtWAarknOdZLuquQ8ai6OC+08//VTPPvusXn31Va2trUmSUqmUenp69Mknn6i9vV0vvvii\n4vG4PvvsM33++ed655139Itf/EL379/X+vq6nn76aX300UePffNScA8AAOqFbwX3Q0NDSiQSOcdW\nV1fV1dWlzs5OSVI0GtXCwoLOnDmjw4cPS5IuXrwoSfrggw/0/PPPM0cOwDes84ViqnkS0ovXo/G4\nrvna2tpSR0dHph2JRLSyslLwta+99lrJax05ciQTxLW0tGhwcDAzaKaLZmnTpk2bNm0T7XRQFJTX\n0w52O/15fhKqGo7X+UokEhoZGclMO167dk2WZWl2dlaSND8/r5WVFU1PT1fWAaYdAQB1Jj+7VSzb\nVassGNk379R0na/29nYlk8lMO5lMKhKJuLoWezsCAILAadCS//ViQVh2sb6fgRBBl3m2bedkw6qx\nw+2Ju3fv1u3bt5VIJPTgwQNduXJFBw8eNNIpAKiWqUESjcXNE4v5C7IWCsxKXbPUk5SsKRZOjqYd\nx8bGdPPmTd29e1e7du3ShQsXND4+rhs3bujkyZNKpVI6evSozp49W3kHmHYE4AHbtsmmwzjT03lM\nD9YfE3FLIPZ2PHfuHNOOAIBQKRRYlQu28rNcQakhw7+mHc+fPx+O4ItvIgBAmHgZJBGA1ZaJuMV1\nzRcABBk1X/Cak/0inby2UvnF/eVQNxY8nu/t6ARPOwIA6k0l2Y9qMiWV7C9p+t74F5NPOzLtCACA\nD0plyaq9XiXXKHYe05nO1HSdL5PIfAEAwq7cL2ynC7cWul4lgVj29bMDCYKu0sh8AUAZLDWBIAt6\nlino/aul0GS+AABoJF4FNZVmz4qdx3Skt8h8AQBQY9UEOG4DrFoJSj/cCk3mi5ovAEAjcxo4FXqt\n00DA6eu8Do7qNeii5gsAyqDmC43OSRDlJtCq98xVtUKT+QIAAGYVChBK1XZVc91iGj1QK4bMFwAA\nqEi5PSjDLDSZL2q+AACoH15kv4KeJaPmCwDKoOYL8E924BT2rFhoMl8AAKB+ud13slGR+QIAAEaY\nynoFeQqSzBcAAPCEmwDI1LIW+a8xtXJ/UOyodQcAwAumCmOBRrW9vW00aKn0eumNv9OZpvypTadB\nXKXvIX1PLwUi88XTjgAAhJfpLJrba1ZzX552BAAABQV1qi1fvfQzHzVfAAAgR70EM/XSTy9Q8wUg\nlKj5AirnR72T1/fx6z1Ug8wXAACQ5F82qtjTjCb6Ue5JySDwNPNl27aGhob05ptv6ubNm17eCgBy\n8AAP4K3spxGrzTSln0os93Ri/r2c3NvpE49+Zsw8zXzt2LFD3/72t/X1118rEol4eSsAAFCEF9kf\nN9eqth/551W67lcl1/aSo8zXxMSEWltbNTAwkHPcsiz19vaqu7tbU1NTj503NDSk//mf/9G7776r\nc+fOmekxADhAzRfwL6bX7DLVD1PZJifvL0i1YI6Cr/HxcVmWlXMslUppcnJSlmVpfX1d8XhcGxsb\nmpub06lTp3Tnzp3Mm2xpadHXX39tvvcAAMAz+VOLpoMXN0GT235UEoB6Hag5mnYcGhpSIpHIOba6\nuqquri51dnZKkqLRqBYWFnTmzBkdPnxYkvTRRx9pcXFR9+7d01tvvWW04wBQCjVfQPUqXZG+0nPc\n9KFcsX6p17q9p5N7VcJ1zdfW1pY6Ojoy7UgkopWVlZzXfP/739f3v//9stc6cuRIJohraWnR4OBg\nZuBMTx3Qpk2bNm3atIPbXl5eVja/7p8OlAp93bbtoufnB4vFrp/+PJFI6LXXXtMHH3ygajle4T6R\nSGhkZERra2uSpGvXrsmyLM3OzkqS5ufntbKyounp6co6wAr3ADyQPegCqB9+Lw2Rnc1y+lRktX3b\n4fbE9vZ2JZPJTDuZTLp+ojEWi+VEmAAAoDEVq83yqg6r3BIXXtS8uQ6+du/erdu3byuRSOjBgwe6\ncuWKDh48aKRTAFAtsl6AM0F6CrCUWj2xmb7v8vKysZUbHE07jo2N6ebNm7p796527dqlCxcuaHx8\nXDdu3NDJkyeVSqV09OhRnT17tvIOMO0IAAAMKzR9aWIlfRNxi+OaL680NTXp3LlzGh4e5i9VAMZQ\n8wWgnErqy9IF/OfPnw9H8EXmC4BpBF8AvFDTgnuTKLgHYBqBF1Af6qXmzLZtxWIxI9ci8wUAAOBQ\naDJfAGAa2XQAQRWI4ItpRwAA4Ae305xMOwIAABjm5OlHE3GL670dAQAAwqTcpt2mkkVMOwIIJcYU\nANXK3npoeXmZaUcAKIV1vgB4ITQr3BN8AQAQfJWsCB9W1HwBAADfNHLQZVIgar4AwDRqvgAEVSCC\nLwruAQBAkLHOFwAAQA2wvRAAAEAZQdu8m+ALQChRygAgLb1WV1AQfAEAAOTxMltGzRcAAIBDoan5\n4mlHAAAQZDztCABlsL0QAC+EJvMFAADQKMh8AQAAOETmCwAAoM54Gnxtb2/rxz/+sU6cOKHf/OY3\nXt4KAHLwEA+AoPI0+PrjH/+ora0tPfXUU4pEIl7eCgBy3Lp1q9ZdAICCHAVfExMTam1t1cDAQM5x\ny7LU29ur7u5uTU1NPXbe5uamvvOd7+inP/2pfvnLX5rpMQA4cO/evVp3AQAKchR8jY+Py7KsnGOp\nVEqTk5OyLEvr6+uKx+Pa2NjQ3NycTp06pTt37igSiailpeWbG+2gvMwPYZhqCdJ78LMvXt3L5HWr\nvZbb84P0PdGowvD/IEjvgbHF7LXqbWxxFBENDQ3pueeeyzm2urqqrq4udXZ26sknn1Q0GtXCwoIO\nHz6sn/3sZ2pra9PLL7+sxcVFnThxgvV2fBKkwcWtIL0HBkiz1/JzgEwkEq7uhcKC9HPpVpDeA2OL\n2WvVW/DleKmJRCKhkZERra2tSZKuXr2qxcVFzc7OSpLm5+e1srKi6enpyjoQoF3GAQAAyql2qYlm\ntyeaCppY4wsAADQS14VY7e3tSiaTmXYymeSJRgAAgDJcB1+7d+/W7du3lUgk9ODBA125ckUHDx40\n2TcAAIDQcRR8jY2Nae/evdrc3FRHR4cuX76s5uZmzczMaP/+/erv79fo6Kj6+vq87i8AAEBdq/ne\njgAAAI2ExbcAAAB8FLjg6/79+3rttdf0+uuv67e//W2tuwMgJL744gsdO3ZMhw4dqnVXAITIwsKC\nXn/9dUWjUX388ceOzgnctOPc3Jx27typl156SdFoVL/73e9q3SUAIXLo0CH94Q9/qHU3AITMvXv3\n9Pbbb+vSpUtlX+tL5quSvSG3trbU0dEhSXriiSf86B6AOuV231kAKMXN2HLx4kVNTk46ur4vwVcl\ne0NGIpHM+mGPHj3yo3sA6lQlYwsAOFXJ2LK9va3Tp0/rwIEDGhwcdHR9X4KvSvaGfPnll3Xt2jUd\nP36cdcMAlFTJ2PLll1/qjTfe0K1bt8iGASipkrFlZmZGS0tLunr1qn71q185ur7r7YWqlT29KEmR\nSEQrKyt65pln9P7779eqWwDqXLGxZefOnXrvvfdq2DMA9azY2DI9Pa233nqromvV7GlHNtQG4AXG\nFgBeMDm21Cz4Ym9IAF5gbAHgBZNjS82CL/aGBOAFxhYAXjA5tvgSfLE3JAAvMLYA8ILXY0vgFlkF\nAAAIs8BtLwQAABBmBF8AAAA+IvgCAADwEcEXAACAjwi+AAAAfETwBQAA4COCLwAAAB8RfAEAAPiI\n4AsAAMBH/w+WeFK5JU4tEgAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}