PublishedSolutions/Exercise06-Solutions.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:0c8e3928c748d24e54bd070267d65728b3aee2bbe2658ebb0e6c451a39b66764"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "\n",
      "import os\n",
      "import sys\n",
      "import glob\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "import pandas as pd\n",
      "%matplotlib inline\n",
      "%precision 4\n",
      "plt.style.use('ggplot')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.optimize as opt\n",
      "import scipy.stats as stats\n",
      "import scipy.linalg as la\n",
      "from IPython.display import Image"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Notes**\n",
      "\n",
      "From this point in the course onwards, homework instructions will be briefer than before with less hand-holding and explanation. You are expected to do the necessary research to solve the given problems on your own."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercise 1 (10 points)**.\n",
      "\n",
      "A milkmaid is at point A and needs to get to point B. However, she also needs to fill a pail of water from the river en route from A to B. The equation of the river's path is shown in the figure below. \n",
      "\n",
      "- What is the minimum distance she has to travel to do this?\n",
      "- Modify the code below so that the path taken by the milkmaid is shown as a dashed line.\n",
      "\n",
      "```\n",
      "x = np.linspace(0, 10, 100)\n",
      "plt.xlim([0, 10])\n",
      "plt.plot(x, 10.0/(1+x), 'b', linewidth=5)\n",
      "plt.scatter([2,8], [8,4], c=[0,1], s=100)\n",
      "plt.text(8, 5, 'B (8,4)', fontsize=16)\n",
      "plt.text(2, 9, 'A (2,8)', fontsize=16)\n",
      "plt.text(0.5, 2, 'River', fontsize=16)\n",
      "plt.text(0.5, 0.5, 'g(x y) = y - 10/(1+x)', fontsize=16);\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f(p, A, B):\n",
      "    return la.norm(p-A) + la.norm(p-B)\n",
      "\n",
      "cons = {'type': 'eq',\n",
      "        'fun': lambda p: p[1] - 10.0/(1+p[0])}\n",
      "\n",
      "A = [2,8]\n",
      "B = [8,4]\n",
      "sol = opt.minimize(f, [0,0], constraints=cons, args=(A, B))\n",
      "sol.x\n",
      "\n",
      "print \"Minimum distance is\", f(sol.x, A, B)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Minimum distance is 9.96339829964\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(0, 10, 100)\n",
      "plt.xlim([0, 10])\n",
      "plt.plot(x, 10.0/(1+x), 'b', linewidth=5)\n",
      "plt.plot([2,sol.x[0]], [8,sol.x[1]], '--')\n",
      "plt.plot([sol.x[0],8], [sol.x[1],4], '--')\n",
      "plt.scatter([2,8], [8,4], c=[0,1], s=100)\n",
      "plt.text(8, 5, 'B (8,4)', fontsize=16)\n",
      "plt.text(2, 9, 'A (2,8)', fontsize=16)\n",
      "plt.text(0.5, 2, 'River', fontsize=16)\n",
      "plt.text(0.5, 0.5, 'g(x y) = y - 10/(1+x)', fontsize=16);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEECAYAAADXg6SsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8k9X+wPFPko50pJtdCmVjoWxZUoZQmYKogAVkiMhy\ncJWpXrxeBwgUlKnXn+IAZVdAloIMARGRIWW1QJEhlO490uT3R2xCTDdt06bf9+uVF83JM04Oyfc5\nOec85yj0er0eIYQQNkFp7QwIIYQoPRLUhRDChkhQF0IIGyJBXQghbIgEdSGEsCES1IUQwobYFfTi\nypUrOXXqFG5ubixevBiAr776it9//x07Oztq1KjBlClTcHZ2LpfMCiGEKFiBNfWePXsyd+5cs7RW\nrVqxePFiFi5cSK1atdi6dWuRThQeHl7yXNoYKQsTKQsTKQsTKQuT4pZFgUG9efPmuLi4mKUFBgai\nVBp2a9y4MbGxsWWSMVsmZWEiZWEiZWEiZWFSqkG9MPv376dt27YPcgghhBClqMRBfcuWLdjZ2fHI\nI4+UZn6EEEI8AEVhc79ER0ezYMECY0cpwIEDB9i3bx9vvvkmDg4Oee4XHh5u9rNh2LBhpZRlIYSo\nWjZs2GD8OyAggICAgHy3LXD0S15Onz7Ntm3beOutt/IN6Pmd+Pbt28U9nU3SaDQkJydbOxsVgpSF\niZSFiZSFSe3atYtVKS6wpr506VIuXLhAUlISHh4ePP3004SFhaHVanF1dQWgSZMmTJgwoUgnk6Bu\nIB9YEykLEykLEykLk9q1axdr+0KbX0qTBHUD+cCaSFmYSFmYSFmYFDeoyx2lQghhQySoCyGEDZGg\nLoQQNkSCuhBC2BAJ6kIIYUMkqAshhA2RoC6EEDZEgroQQtgQCepCCGFDJKgLIYQNkaAuhBA2RIK6\nEELYEAnqQghhQySoCyGEDZGgLoQQNkSCuhBC2BAJ6kIIYUMkqAshhA2RoC6EEDZEgroQQtgQCepC\nCGFDJKgLIYQNkaAuhBA2RIK6EELYEAnqQghhQySoCyGEDZGgLoQQNkSCuhBC2BC7gl5cuXIlp06d\nws3NjcWLFwOQkpLCkiVLiImJoVq1akyfPh0XF5dyyawQQoiCFVhT79mzJ3PnzjVLCwsLIzAwkA8/\n/JAWLVoQFhZW5JOlpSnQ6UqWUSGEEIUrMKg3b97cohb+22+/0b17dwB69OjBiRMninyyJk1qEhFR\n4I8DIYQQD6DYbeqJiYl4eHgA4O7uTmJiYpH31esVnD1rX9xTCiGEKKIH6ihVKBTF3uePPySoCyFE\nWSl2W4i7uzsJCQl4eHgQHx+Pu7t7ntuFh4cTHh5ufD5s2DAAzp9Xo9HoS5hd2+Dg4IBGo7F2NioE\nKQsTKQsTKQtzGzZsMP4dEBBAQEBAvtsWO6i3b9+eAwcOMGTIEA4ePEiHDh3y3C6/E585oyQhIRmV\nqrhnth0ajYbk5GRrZ6NCkLIwkbIwkbIw0Wg0xkpxURQY1JcuXcqFCxdISkpi8uTJDBs2jCFDhrBk\nyRJ++ukn45DG4khLU3L1qh2NG2uLtZ8QQojCKfR6fbm1heQ2wS9bFs/QoenlddoKR2ohJlIWJlIW\nJlIWJrVr1y7W9la5o1RGwAghRNmwSlCXETBCCFE2rBLUz52zlztLhRCiDFglqKekKLl2rQoPf/mH\nGTNm4Ovry1tvvVXsfVetWkVwcLDx+ZUrV3jjjTfo0aMHTZo0oW3btowbN47z588X6XharZbVq1fT\ns2dPGjduTNu2bZkwYQIXLlww227ChAnMnj272PkVQpQtq83S+McfDtY6dYWSnp7O9u3bUavVhIWF\nkZOTU+R94+LiWLZsGTNnzjSmHTp0iKNHjzJixAjWrFnDe++9R2xsLIMGDeKPP/4o9Jjz58/nvffe\no3///nzxxRf85z//4fr16zz99NP89ddfxu1mzJjBt99+y6VLl4r3hoUQZcqKQV3a1QH27NlDSkoK\nr7/+OjExMfz0009F3vfLL7/Ezc2N3r17G9MGDx7M/v37mTRpEl26dKFv376sXbsWtVrNp59+Wugx\nN2/ezOOPP86MGTPo0qULgwYNYtWqVSQkJLBv3z7jdk2bNqVdu3Z8/PHHxXvDQogyZbWgLiNgDDZu\n3Ejjxo0ZO3YsNWrUYOPGjUXed926dQwePNgszcvLy2I7jUaDv78/d+/eLfSYOp3O4k4+Nzc3AP45\n+nXIkCFs27aN1NTUIudZCFG2rBbUz52zp/xGyFdMd+7c4eeff+bxxx9HoVAwaNAgfvzxxyJNkhYR\nEcHt27d5+OGHC902Pj6eixcv0rhx40K3nTBhAlu2bGHv3r0kJydz/fp15s6dS+3atRk0aJDZtg8/\n/DAZGRkcO3as0OMKIcqH1YJ6UpKSqKiq3Vm6ZcsWcnJyePzxxwFD00lmZibbtm0rdN/Tp08D0KxZ\ns0K3ffPNN1EoFEyYMKHQbV988UUmTpzIc889R/PmzenatSuXL19m48aNxtk5czVp0gQ7OzvOnj1b\n6HGFEOXDqsvZVfV29Y0bN/LQQw/RsGFDANq0aYOfn1+RmmBiYmIA8PT0LHC7ZcuWERYWxjvvvEO9\nevUKPe7y5ctZsWIF06dPZ9OmTXz88ce4uroSEhJi0XyjUCjw8PAgOjq60OMKIcqHBHUrOXPmDBER\nEfTu3ZvExETjo3fv3vz+++9cvXr1gc/x5ZdfsmDBAmbNmsXw4cML3T42NpZFixYxZcoU/vWvf9Gp\nUycGDBjAunXriI2NZdWqVRb7lOMsE0KIIrBqUD97tuoOa8ytjX/00UfGGS0DAgL47LPPANi0aVOB\n+/v4+ACG9vK8bNq0iddff51Jkybx4osvFilP165dQ6vVEhgYaJbu4eGBn58fV65cMUvX6/UkJiZS\nvXr1Ih1fCFH2rF5Tr4p3lmZlZREWFkbbtm3ZtGmT2WPjxo0EBASwefPmAo/RqlUrgDxvKtq1axev\nvvoqISEhvPHGG0XOV61atQDDr4j7xcfHExUVRc2aNc3SL1++nOdFQAhhPeW6YKiLi47UVNN1JDFR\nyaVLdjRvXrWm4d2zZw8JCQk8++yzdOrUyeL1UaNGMWfOHI4ePUqXLl3yPEaTJk2oVasWx48fp0+f\nPsb0X375halTp/LQQw/x9NNPc/LkSeNrjo6OtGjRwvj8lVdeYdOmTdy8eROAOnXq0LdvX1atWoVC\noaBjx47Ex8ezatUqtFotzz77rFkefv31VxwdHfN8D0II6yjXmnq7dlkWacePV70mmG+++QaNRsPA\ngQPzfH3IkCGo1epCm2BCQkLYtm2bWbv20aNHycrK4ty5cwwZMoTBgwcbH88//7zZ/mlpaRZNJ8uX\nL+eVV15h586djB8/nrfeeovq1auzdetWWrZsabZtWFgYgwYNwtXVtThvXwhRhsp1PvUZM5JYtMjN\nLG3w4DRWrkworyxUCKU1V3RsbCxdu3blo48+Mpv/pajatm3LxIkTmTRpUrH3vXz5MsHBwezevbtI\nwyrzI/Nmm0hZmEhZmFTo+dQ7dsyrpu5Y5W9CKilvb29efvllQkNDi73v1atXyc7OZsyYMSU698KF\nCxkxYsQDBXQhROkr15r6lSu3ad68FtnZCrP0Y8fu4udX9ImsKjuphZhIWZhIWZhIWZhU6Jq6kxME\nBmZbpFfFdnUhhCgL5T6ksWPHTIu0X3+VoC6EEKWh3IP6ww/LCBghhCgr5R7UO3TIQqEwb8a/csWe\nmBir3gclhBA2odwjqYeHnmbNLG82kiYYIYR4cFapHksTjBBClA2rBHXpLBVCiLJhlaDeoYNlTf3c\nOXtSUhR5bC2EEKKorBLUa9fWUbeuebu6TqfgxAmprQshxIOw2pCTvKYMOHDA0Qo5KV8KhUIWlhBC\nlJlynXr3ft27Z7Jpk7NZmi0H9Zs30zh1Kppt266g1+vp18+f9u1rUq+ei7WzJoSwISUO6lu3buXw\n4cMoFAr8/PyYMmUK9vZFX56ue/dMFAo9er2pHT0y0p4bN1TUrWtb88CEh8czfPgO4uMzjGm7dl3D\n1dWeb74ZSNu2PlbMnRDClpSo+SU6Opp9+/axYMECFi9ejE6n48iRI8U6hre3Ls95YH76ybZq67du\npTFihHlAz5WSks0zz+wgKirFCjkTQtiiEgV1Z2dnVCoVmZmZ5OTkkJmZiZeXV7GP06OH5dBGW2uC\nOXMmhrg4y4CeKyUlm+PH75RjjoQQtqxEQd3V1ZVBgwYxZcoUXnjhBVxcXEq0TmXPnpbB7uefHcmy\n7EOtlBQKBT/+eN0srZ1HHGDeUbply2VAhnMKIR5cidrU79y5w/fff8+KFStwdnYmNDSUw4cP061b\nN+M24eHhhIeHG58PGzYMjUZjdpygIMO0AQkJpoCWmqrk/HkPunWr/O3q949yaapJ4vVm5/FzTuOp\nY12JyTL/RWL49VP15r9xcHCw+FxUVVIWJlIW5jZs2GD8OyAggICAgHy3LVFQv3r1Kk2bNjUWeseO\nHbl06ZJZUM/rxHlNev/II3bs2OFklrZzZw6tW9vGBPn9H/GgbfgZete4w7LIJnz9Zz20evPgPXRo\nE9LSUq2UQ+uSxRBMpCxMpCxMNBoNw4YNK/L2Jaoa1q5dm4iICLKystDr9Zw9exZfX9+SHCrPJpj9\n+9UlOlZFo79yke4H3ifLTk3Pg71Yc93fIqC7utrTsWNNK+VQCGFrSlRTr1+/PkFBQcyePRuFQoG/\nvz+9e/cuUQa6d7fsLL1wwZ47d5TUrKkr0TErDL+G2L0ZSsdoR1TDd8A/RsBoNA6sWzeA+vVdrZRB\nIYStKdc1Sm/fvp1neu/e1bhwwXyMe2hoPMOHp5dHtsrFzZupnD59j+++M9x81L9/A9q1q1Hlbz6S\nn9kmUhYmUhYmxV2j1Gp3lN6vZ88Mi6D+44/qShPU9X/dhMQ4FM3yHwHk6+uCr68Lgwb54+LiQkqK\njE0XQpS+CjHcIq/x6vv3O5KWVrGH+emTEtCtXY1u4Rz0sdFF20evR6Go2O9LCFF5VYig/vDDWXh6\nmg9hzMhQsn9/xbwRSZ+ViW7XJnTzpoJKhfLtFSi7lqxPQQghSlOFCOr29tC3r+UomO+/d8pja+vT\nrV6A/tpllLMXohzxPApXN2tnSQghgArSpg4wYEAG33xj3mn444+OpKeDUwWL7cqJr6FQOxe+oRBC\nlLMKUVMHeOSRTDw8zIcwpqUpOXCg4o1Zl4AuhKioKkxQt7eH4OC8mmCsE9T1KUnoNn+Bvore6SmE\nqJwqTFAHGDDAcgjjDz+oych/ksNSp8/ORrd3K7o3p0BGOsgqRUKISqTCtKkDdOuWiZubjqQk07Um\nJUXJoUOOBAdbDnssTXq9Hv1vR9Bv+QJq+6Gc+T6KWnXL9JxCCFHaKlRN3dER+vSxrJb/c8KvMvHn\nVfS7N6F8dhqqF9+UgC6EqJQqVFAHGDjQsglm796yb4JR1GuI8o0lKJq3KtsTCSFEGapwQT0oKBNX\nV/NRMMnJSvbsKfsOU7nTUwhR2VW4oK5W530j0rffPvgwQr1Wi27fDnTfrX3gYwkhREVU4YI6wIgR\naRZphw87cuuWqkTH0+v16E/9gm7eNPR/nEDRruuDZlEIISqkCjX6JVenTlnUq6fl+nVT9vR6BRs3\nOvHKK8Wb3VAfFYFu42eQmoLymYkoWrQt7ewKIUSFUSFr6goFPP20ZW19wwZndMVcN0P/288oOvZA\n+eZSCehCCJtXIYM6wLBh6SgU5jf+XL9ux/HjDsU6jvKpcSiDHkOhKlnTjRBCVCYVNqjXqZNDt26W\nNxytX593h2k5LuAkhBAVVoUN6kCeKx/t2KEmJcU09FCv16M/ewLdWy+iv3OzPLMnhBAVToUO6o89\nlo67u3kjenq6kq1bDXeY6m9cQ7fk3+g2fo5y6BioUcca2RRCiAqjQo5+yeXkBIMHp/Pll+bzrG9e\nk0WI9kM4dxLFwBEougWjsKvQb0UIIcpFha6pA4wcef/Ut3oUSh2RVxy4kVQd5X9XoezZXwK6EEL8\nrcIH9RYttHTqZOgwtXfPpPXbh6g5OYLXMh5j540szkenkZadU8hRhBCiaqgUVdwJE1L55RdHshPV\nnPlPN5xqpRBTO5k/2t7jp2uJOKgUvNennrWzKYQQVlcpgnpwcAZ162q5ccMOXbaK1D/dSf3TnZTm\nHix6Jynf4Yxn76SyKTwWf081/p6ONPBUU8fNAZVSJu4SQtimCt/8AqBSwbhxlsvKbdjgTFKSIt/Z\nFRt4qRnczAuNg4pfb6bw/qGbjNhwmXVn75V1loUQwioqRU0dDJN8LVqkIS3NdB1KTVXy7bfOTJyY\n9zqirg4q2tVxpV0dV2NaWnYOmdq8a/a/3kzmekKmsWbv5WQn0/EKISqVSlFTB3B31zNsmOXNSJ98\n4kpmMVa6c7ZX4emU97VM46AiJUvHtotxvLwzijFbIpm3/wZ/3JXFp4UQlUOJa+qpqamsXr2amzcN\nd3FOnjyZJk2alFrG8jJ+fApr1piPWf/rLxUbNjgzerTlBGDF1by6M82rG6Yh0Ov1xKVruRafiY+z\nfZ7bn72TiqOdknoejqjtKs31UQhhw0oc1D///HPatGnDq6++Sk5ODpnFqS6XUMOGOfTrl86uXeZr\nli5b5srw4Wk4FG+urwIpFAq8ne3xziegA1y8l84vN1O4kZhJNRd7/D0d8fdU07uhOx7qStOyJYSw\nISWqXqalpXHx4kV69eoFgEqlwtn5wVcmKopXXkm2SLt1y45Nm8rn/Pcb1tKH0H71+WZYE2Z1q0P7\n2q4kZWjJb26x6JRscnQy8ZgQouyUqDoZHR2Nm5sbK1eu5Pr16/j7+zNu3DgcHR1LO38WWrTQEhyc\nzt69lrX1p59Owz7/inWZsVMqqOfhSD0PR3rinuc2er2etw/c4F5qNg28nfFzszcOs2zkrUYpHbJC\niFKg0JdgztorV67wxhtv8N///pdGjRqxZs0anJycGD58uHGb8PBwwsPDjc+HDRtGcrJlLbskfv9d\nSY8eLhbpK1emM2qUtlTOUVZSMrX8mZTNxTtJRMakcSsxg6WDm1XZUTYODg5kZWVZOxsVgpSFiZSF\niUajYcOGDcbnAQEBBAQE5Lt9iWrq3t7eeHl50ahRIwA6depEWFiY2TZ5nbi0gnrjxtCrl4r9+9Vm\n6R98YE///vFWqa0Xx0PVNNRVa+lT39BklJJiuUTfX8lZzPnhTxr83U6f215fS2NvU7V6jUZTap+L\nyk7KwkTKwkSj0TBs2LAib1+iNnUPDw98fHy4ffs2AGfPnsXX17ckhyqx6dMt/8OjouxYu7b829bL\nQk1Xez4IrsdjjTywVyo4FJXEW/tv8O99N6ydNSFEBVai5heAqKgoPv74Y7RaLTVq1GDKlCmFdpbm\nXgRKy8iRXhw4YF5b9/LK4ciRaNzcKm6H5IPUQrJzdNirLK/Ff9xN5YfIRGONvoGnI26VYASO1MhM\npCxMpCxMateuXaztSxzUS6K0g/qZM/b071/NIn3atGTmzKm4H4iy+MDGpmVz+q9UrsVnci0+g2vx\nmajtlDzxkBeDmnmV6rlKk3x5TaQsTKQsTIob1Ct+Va4ArVplM3RoGlu2mP9C+N//XHn22TTq1Kk6\nU/J6O9vzaEMP43O9Xk90ana+wytP/ZVKTGo2/p5q/DwccMij9i+EqHwqdVAHmD07me+/dyIz09R5\nmJmpYP58DcuWJVgxZ9alUCio4Zr/3Vg5Oj3notPYfimev5KzqOXqgL+nI4OaedHIW53vfkKIiq3S\nV8/q1Mnh+ectR49s2eLM6dMVfBiMFbWv48r0LrX5aIA/a59uzMtdatGypjNqu7xH1kTEpnMnOSvf\naY6FKG3r16/H19fX+PDz86Ndu3ZMmjSJK1euFPk4q1atIjg42CwtNjaW119/nc6dO9OwYUM6d+7M\nG2+8QVxcXLHyeP36dRo2bIivry/Xr183e23ChAnMnj27WMcrDZW+pg4wdWoK69Y5ExenMkufM8ed\nHTtiUKny2VEA4KBS0tBLTUOv/GvoB64lcexGMunZOup7OOLvZeiM7VxXg4uDFLAoO5988gm1atUi\nJyeHqKgoli5dyvDhw/npp5/QaDQF7hsXF8eyZcv46KOPjGl6vZ4xY8bw559/MmPGDBo3bsylS5dY\ntGgRZ86cYfv27UXO29y5c3F3d+fePcvpvGfMmMFjjz3GuHHjaNq0adHf8AOq9DV1ADc3Pa++atmp\ncvasg8UEYKJknm9fg8+eaMTHjzdgRKAP1V3s+ONOGtn5THsgSwyK0hIQEECbNm1o3749Tz31FO+/\n/z537tzh5MmThe775Zdf4ubmRu/evY1pV65c4fTp08yaNYvRo0fTqVMnxowZw8yZMzl16hRXr14t\nUr62bt1KeHg4U6ZMyfMXbNOmTWnXrh0ff/xx0d9sKbCJoA4wcmQazZplW6R/8IGGv/6ymbdpdW5q\nO1rVdGFIc2+md62d58RlOTo9L3x3lefDInnv4E2+OXuPX24kczdFmm/Eg3N1NayPkJ1t+X3/p3Xr\n1jF48GCzNJ1OB2BRy3dzczN7vSAJCQm8/fbb/Pvf/y7w18KQIUPYtm0bqanlN323zUQ7e3uYP9+y\nYzQlRcm//533fCyibKiUCr54shFvP+pHd383cnTwQ2QC7x64Ze2siUpIq9Wi1WrJzMwkIiKC+fPn\nU61aNbp06VLgfhEREdy+fZuHH37YLL1JkyZ0796dDz/8kLNnz5KamsqpU6dYsmQJvXr1Mt4pX5B3\n332XRo0aMXTo0AK3e/jhh8nIyODYsWOFv9FSYhNt6rk6dMhm5MhU1q41b3LZudOJH39Mo3fvsp8e\nWBgoFQpqaRyopXGgq1/B295MymTRz7eNN035e6qp7+mIq7TVC6B79+5mz2vUqMEXX3yBi0vBTaun\nT58GoFmzZhavffrpp7zwwgv079/fmNa7d29Wr15daH6OHz/O5s2b2bt3b6HbNmnSBDs7O86ePWvW\nBFSWbKamnmvOnCS8vS3bc2fP9iAx0XbmTLElNVwcmNqxJs18nLiVlMWXp+/x3NZIFhyWmr2Azz77\njF27drFz504+++wzmjRpwqhRo4iMjCxwv5iYGAA8PT3N0nU6HRMnTiQ8PJwFCxawZcsW5s+fz5kz\nZ5g4cWKBTYRZWVnMmjWL559/vkg1eoVCgYeHB9HR0UV4p6XDpmrqAJ6eeubNS+Kll8z/I//6S8Ub\nb7hX6bHrFZW9SkFjbycae5umU87R6UnNyruzNfxuGr/cTDbW7H3dHbFTygXbVjVt2pR69eoBEBgY\nSFBQEB06dGDx4sWsWrWq2Mf74Ycf+Omnn1i/fj1du3YFDM0kfn5+hISE8MMPP1gMgcz16aefkpSU\nxPjx40lMTAQgPd2wzGZycjIpKSnGNv9c5d2PZHNBHWDo0HQ2bnTm8GHz+d23bHHmsccyGDgww0o5\nE0WlUirynbvGw8kOD7Udv99OYXN4LNGp2dR1d2BgUy96NZD+E1unVqvx8/Pj4sWLBW7n4+MDQHx8\nvNm8VJcvXwagVatWZtu3bt0agMjIyHyDekREBNHR0bRr187itb59+xIQEMCePXuMaXq9nsTERKpX\nr16Ed1Y6bDKoKxSwaFECvXtXIznZvIVp1iwPOnSIpkaNwnu4RcVUx82BJwO8jc8ztDquJ2TiZJ93\na2L43TRSs3Pw91Tj42xXZeeutxXp6elERUXRvHnzArfLDdrnz5+nTp06xvRatWoBhjb3Rx55xJh+\n6tQpAGrWrJnvMadOnWq2bgTATz/9xIoVK1i+fDkNGzY0e+3y5ctotVoCAwOL8M5Kh00GdQBf3xze\nfjuR6dPNm2ESEpS89poHX34Zh3y3bYPaTklTH6d8X7+Xls2Ba0lcjc9Ap8e44lSfhu74upf9al3i\nwZw7d46YmBjDfEbR0Xz++eckJSUxbty4Avdr0qQJtWrV4vjx4/Tp08eY3q9fPz744ANefvllXn75\nZRo2bEhkZCRLliyhTp069OvXz7jtK6+8wqZNm7h58yYAjRo1smhLz72TtE2bNsZmoly//vorjo6O\ndOrU6YHKoDhsNqgDPP10Onv3qi0Wqt6/X83HH7swaVL5jR0V1tPD350e/oZmmbh0LdfiDLNY5rdc\n7K2kLDydVDjby+gba8r9RfXCCy8Y07y9vWnWrBlr164lKCio0GOEhITwzTff8MYbbxjTXFxc2LZt\nG6GhoaxevZq7d+9SvXp1goODefXVV3FyMsWLtLS0IjWd5PfrLywsjEGDBlm0s5elSj31blHExip5\n9NFq3Ltn/gVVqfRs3BhLx47lv2SWTCtqUhHLYuHPtzhxMwUvZzuzYZYtazjjaFd2A8YqYllYS2mV\nRWxsLF27duWjjz7Kt528IG3btmXixIlMmjSp2PtevnyZ4OBgdu/eneewyqIq7tS7Njek8Z+8vXUs\nXGg54iUnR8HkyZ7cu2fzRSCKacYjdfhmWBNe7+5L57oaUrN0bL8YR6Y2734YndwlW2F5e3vz8ssv\nExoaWux9r169SnZ2NmPGjCnRuRcuXMiIESMeKKCXhM3X1HP95z9ufPKJ5U+grl0z+eab2HKd9Etq\nZCaVvSyyc/SM2RxBbTeH+1adUlPPwzHfjtv8VPayKE1SFiZSU8/H3LlJtG9v2dRy5Igj77/vZoUc\nCVtgr1Lw6RMNea5tdep7qLkSl8Env93lpe+vWTtrooqy6Y7S+9nbw6pVcTz2WDWLKXpXrXKlceNs\nhg9Pt1LuRGXmbK+ieXVnmlc3jYXOr0nmZmImn56MNltLtpbGAZXcPCVKSZWpqQPUrq1jxYoEFArL\nL9ysWR4cO5b/SkFCFIcyn9EQXs52DGjiiZOdkqN/JvPfAzcJ2XiZ1b/eKeccCltVZdrU77dsmSvz\n51s2uXh46Pj++3vUr1+2c4FLe6GJlAWkZuWQkpVDo1reFmVx8V464dFpxrH1Hk5V48e1fC5MqtTC\n0yU1bVoKERF2bN5svmB1QoKSUaO8CQuLwcdH7jgV5cPFQZXv6lH2KgUJGVq2nE/lWnwG9koF/p5q\n+jbxoKMWITVNAAAgAElEQVRvwav+iKqpSgZ1hQIWLkzgzz9VnDhhfkfhtWt2jBzpxaZNsWg0MlRN\nWNf9ywzq9Xpi0rRcjc/Axznv9XcjYzPQo8fP3bFMx9SLiqtKBnUAR0f4v/+LZ8AAH27cMC+Gc+cc\nGDfOi6+/jkWd/7KdQpQrhUJBNRd7qrnkv6D6hXtp7LuayK2kLGq42uPvqcbf05Fu9dwK3E/YjirZ\npn6/S5fsGDrUh4QEy1pN377prF4dj30pfxekvdBEysKkNMsiO0fPzaRMrsVncjU+gz4NPajnYTnP\nTUK6Fje1Kt+OXWuRz4VJcdvUq3xQB/jtN3tGjPAmPd0ysPfvn87KlaUb2OUDayJlYWKNsnj9h+tc\nicukvqejcZhlbqesNYdZyufCRG4+KoH27bP53//isbOzvL7t3OnElCmeFGGNWyEqnXf71ON/Qxoy\nspUPtTQOXLiXxvJf7pCez5QIouKTmvp9wsKcmDbNA73esobSr5+hxu5QCkPZpRZiImVhUtHLIlOr\nY+r2q3/X6tXGmn0NV/tSab7Ra7WkXL9Oyu3b2Nnb41ynDi6+vlT1ObLLtflFp9Mxe/ZsvLy8mD17\ndqHbV/SgDrB2rTMzZ3rk+VrPnhl88kk8zs4Pdh2s6F/e8iRlYVLRy0Kv13M3JdvYTn8tPpNr8Rmo\n7ZSsGNTggY6dfvs2p5Yv58Latei0WgDsXVxo98orNBk+HEdv70KOYLvKdZz6zp078fX1Na7RZwtG\njkxDoYAZMywD+08/qRk+3Jsvv4zF01OGO4qqRaFQUFPjQE2NA539TGPkM/JpqrmRmMmm8Fga3Fer\n1zhajsfPuHePH6dO5c6vv5qlZ6em8su775Jy6xYd5s7FzsWldN+QjSpxm3psbCynTp2iV69e5b6w\nalkLCUlj8eL4PKcT+P13B4YO9eH2bemOEAIMK0/lxc1RRcsaztxNyeKbszE8H3aFCVsj+ebsPbPt\nYv/4wyKg3+/cmjUkX71aqnm2ZSWuqX/xxReMGjXKpmrp9xsxIh2FAl57zQOdzrxN7/JlewYNqsaa\nNXG0bCk9qELkxV1tR++GHvD3sp26v5tvzOalz8nh3OefA5Dg25TkmvVxuxWJ219XuP9bd+voUTxb\ntiy/zFdiJQrqJ0+exM3NDX9/f8LDw/PcJjw83Oy1YcOGodFUrtuaJ0yAmjUzGDdOTWameWC/c0fF\n0KE+fPppBgMHaot1XAcHh0pXFmVFysKkKpSF+z+mXMpMTibt7l0AtA5qYhq15Wa7YDp//KrZdqm3\nb+Pq6lplFw3fsGGD8e+AgAACAgLy3bZEHaXr1q3j8OHDKJVKsrOzSU9Pp2PHjkybNq3A/SpDR2le\njh413GGakmL5M1Oh0DNnTjJTpqQUuZO+oneIlScpC5MqWRZ6PQenT+fyxo0Fbtb9gw9oMnJkOWWq\nYin3m4/Onz/Ptm3bbGb0S37OnbNj1Chvi7VOcw0alM7ixQm4uBRenFXyy5sPKQuTqloW0ceP893Q\nofm+rlAqGbZ/P26NG5djrioOq9x8VBV+ErVooWXHjhiaN8+7DX37dicGDPAhMrLKTqcjRIl4Nm9O\ni3Hj8n292/vv41q/fvllqJKTm4+KKSVFwZQpnuzbl/dMXy4uhoWuBw/OyPcYVbVGlhcpC5OqXBZZ\nCQncOnCAXxcsIOnPPwGo1ro1HWfNonqHDqicnKycQ+uRuV/KQU4OvPNO3gtZ5xo+PI3//jcxz+aY\nqvzl/ScpCxMpC8iMiyPj3j3sHR1x8PbGzsY7jotC5n4pByoVzJuXxIoV8Tg55X3jxfr1zvTtW42z\nZ2W6UyGKytHLC/emTanRsqUE9BKSoP4AhgxJZ8eOGPz98x7SePWqHYMG+RAa6ioTggkhyoUE9QfU\nrJmWnTvvMXBg3jdhabUKFi92Y+BAHy5ckE5UIUTZkqBeCtzc9KxeHc+iRQn5NsecO+dAv37VWLhQ\nQ0b+fahCCPFAJKiXEoUCnnkmjd27YwgIyLutJTtbwdKlGrp0ceHo0VKYw1cIIf5Bgnopa9RIy44d\n93j55WRUqrwHFkVGKnn6aR9efNGDu3flv0AIUXokopQBBweYOTOZbdtiaNw4/x7SLVuc6datOqtW\nuZCVVY4ZFELYLAnqZah162x27zbU2vNaKg8gNVXJO++407Nndb7/Xo2NzWIshChnEtTLmFptqLXv\n2XOPdu3yr45HRdkxcaIXQ4b48NtvMrZdCFEyEtTLSbNmWsLCYpg/P6HAVZN++82BwYOrMXasF+fP\nyxBIIUTxSFAvR0oljB6dxqlTKYwZk4pSmX9w/+EHNX36VGfyZE8iIiS4CyGKRoK6FXh5wXvvJbJ7\n9z26ds0scNtt25zo2bMaL7zgKTV3IUShJKhbUUCAlvXrY/nyy1iaNMl/lIxer2DHDif69KnO2LFe\nnDghY9yFEHmToG5lCgU8+mgmP/xwj0WLEqhdu+Cl8X74Qc2QIT48/rgPu3apyckpp4wKISoFCeoV\nhJ2d4Y7Un3+O5p13EqheveBoffKkAxMmePHII9X55BMXkpJsf6ESIUThZD51KyjKvNnp6Qq+/tqZ\nVatcuXs37yX07ufiouOpp9J59tlUmjUr3kLY1iRziJtIWZhIWZjIIhmVQHE+sBkZ8O23zqxc6cqt\nW0XrKO3YMZPRo9Po1y8ddd4LNFUY8uU1kbIwkbIwkaBeCZTkA5udbVgHdeVKVy5cKNrNSe7uOoYO\nTWPEiDRatKiYtXf58ppIWZhIWZhIUK8EHuQDq9fDoUOOfPqpC/v3F70a/tBD2Tz1VBpPPJFO9ep5\nTw9sDfLlNZGyMJGyMJGgXgmU1gc2MtKOzz93YcMGJ9LSitbnrVLpCQrKZMiQdPr2zcDV1bqTzciX\n10TKwkTKwsSm1yhdv349vr6+xoe/vz+PPPIIS5YsQavVWmx369YtK+a27DVqpOXddxP5/fe7vPtu\nAs2aFb5mXk6Ogp9+UvPyy560alWTF17wZPt2NWlpMnpGCFtQKW9R/OSTT6hVqxYpKSns2rWLxYsX\nk5GRwZw5cwDo06cP27dvp1q1albOafnQaPSMHZvGmDFpnDjhwLp1zmzfriYjo+BrdkaG4aamHTuc\ncHLS8eijmfTvn06vXploNDJdpBCVUaVqflm/fj2vvvoqR44coV69esb0Z555hlOnTnHx4sUHzWKJ\nZWZm4ujoWKRty+OnZVKSgu++c2L9emdOnSreHagODnoeeSSTvn0z6N07gxo1yq4NXn5mm0hZmEhZ\nmNh080t+WrRoQUpKCnFxcYBl88vo0aPp27evxX53797Fz8+PTz/91Jj2559/Mm3aNAIDA2nQoAHB\nwcHs3r3bbL/Fixfj6+vLpUuXCAkJoUmTJkyePLkM32HxubnpGT06jR07Yjh4MJoXX0wu9G7VXFlZ\nCvbvVzNzpgdt29ZkwAAflixx5exZe3QVp49VCJEHmwjqN27cwN7eHo1Gk+frTz31FOfOnSMiIsIs\nfevWrSiVSp544gkAbt26xcCBA7l48SL/+c9/WLNmDS1btuT5559n7969FscdP348Xbp0Yc2aNUyc\nOLH031gpadRIy+zZyRw/Hs2WLTE8+2wqnp5Fn1/g9GkHFi1yo1+/arRtW4Pp0z347js1cXHSDi9E\nRVMp29S1Wi1arZaUlBR2797Nrl27GDp0KPb2eY/fDg4ORqPRsGnTJmO7O8DmzZsJCgrC29sbgNDQ\nUBQKBZs2bcLDwwOAoKAgbt++zaJFiwgODjY77nPPPcf48ePL6F2WPqUSOnbMomPHLN5+O5EjRxzZ\nsUPNrl1OJCQU7fp+756KDRuc2bDBGYVCT6tW2XTrlkm3bpm0a5dV4W92EsLWVcqg3r17d7PnQ4YM\nYf78+flu7+TkxIABA9i6dasxqF+4cIELFy4wbdo043YHDhygV69eaDQas9E03bt355133iE1NRUX\nFxdjel5NOpWFvT306JFJjx6ZvP9+IkePOrJrl5o9e9RERxc+LQEYZo88fdqB06cdWLZMg1qtp0OH\nLLp0yaRLl0xatcomn+usEKKMlCiox8TEsGLFChITE1EoFDz66KP079+/tPOWr88++4xatWoRExPD\nJ598wg8//MClS5cIDAzMd58nn3ySb7/9lqNHj9KlSxc2b96MRqMxC8wxMTFs3LiRjRs3WuyvUCiI\nj483C+o1atQo3TdmJfb20L17Jt27Z/Lee4mcPGnPjz+q2btXzeXLRY/KGRkKDh925PBhQ4exs7OO\nDh0Mvww6dcqidessitiXLIQooRIFdTs7O8aMGUP9+vXJyMhg1qxZBAYG4uvrW9r5y1PTpk2No1+6\ndu1K7969eemll9i3bx8qVd61zM6dO1OnTh22bNlC586d2bp1KwMGDDAbseLl5UXHjh2ZOnVqnseo\nXr262XOFwvbalJVK6NAhmw4dspkzJ5moKBX796vZt8+RY8ccycws+ntOS1Ny8KCagwcNbTKOjnoC\nA7N4+OEs2rfPon37bPLpBhFClFCJgrqHh4exzVmtVlOnTh3i4+PLLajfz8HBgTfeeIPx48fz7bff\nMnLkyHy3HTp0KGvWrOGxxx7j7t27PPnkk2av9+jRg5MnT9K4cWPU0jgMQP36OYwfn8r48amkpSk4\ncsSBQ4ccOXBAzdWrxfv4ZGYqOHHCkRMnTBfSBg10tGmjom3bLNq2zaZ5c2myEeJBPHCbenR0NFFR\nUTRu3Lg08lMiwcHBtG7dmg8//JCnn3463+2efPJJli1bxuzZs/H19aVz585mr7/22msMHDiQJ598\nkrFjx+Lr60tiYiIXL17kxo0bLF68uKzfSoXm7KynT59M+vTJBJK4cUPF4cOO/PyzAz//7EhsbNHa\n4u939aqSq1ed2bzZGQC1Wk9AQDatW2fRqlU2rVpl06CBFqVNjNMSouw9UFDPyMggNDSUsWPHllvN\nNr8mj5kzZzJy5EjWrl2Li4tLnts1atSIVq1acfbsWbMO0lx16tRh586dhIaGMn/+fOLi4vD09KRZ\ns2ZmFwuFQmGTTS/FVbduDiEhaYSEpKHTwcWLdhw96sixYw788otjkUfU3C8jQ8HJkw6cPGm6YcrF\nRUeLFtm0aJFNy5aGfxs31mJXKbv5hShbJb6jVKvVsmDBAlq3bs2AAQMsXg8PDyc8PNz4fNiwYXKH\n2N8cHBzIysqydjbKlE4H588rOXJExdGjKo4cUREdXXrVbUdHPQ89pKNlyxxatNDRsqWOgIAc/m4V\nrJSqwueiqKQsTDQaDRs2bDA+DwgIICAgIN/tSxTU9Xo9K1aswNXVlbFjxxZ5P5ml0aAq3gKt18P1\n6ypOnHDgxAlDTfzSJTv0+tL9xVOrVg7Nmxva5ps00dKsWTaNGmkrxfj5qvi5yI+UhUm5TL178eJF\n5s2bh5+fn7EZIiQkhNatWxe4nwR1A/nAGiQlKbh0yZ2ff9Zy6pQDp07ZExdX/Hb5wiiVevz8cmja\n1NBsk/to1EiLi0vFmbhMPhcmUhYmMp96JSAfWJP7y0Kvhxs3VJw5Y8+ZMw6cPm3PuXP2JCeXXS9p\nrVo5NGqkpWFD06NBAy116uSUe+esfC5MpCxMJKhXAvKBNSmsLHQ6iIpS8ccf9vzxhwPh4XacO1c2\nNfr7qdV66tXT4u+vxd8/B39/LfXra6lfP4datcom4MvnwkTKwkSCeiUgH1iTkpSFXg9//aXk/Hl7\n4+PCBTuuXrVDpyv7UUmOjnrq1tXi55dD/fpa6tXLwc/P8NzPLwdn55J9peRzYSJlYVLcoC6DwkSl\no1BA7do6atfOpHfvTGN6ejpERhoC/OXL9ly6ZMfFi3bcvl26H/PMTAWRkfZERuZ9l5S3tyG4+/rm\nULeuFl/fHLNHRWrHF7ZHgrqwGU5O0LKlYSw7pBvTk5MVRETYcfmyHRER9kRE2BEZaceff6pKffQN\nQGysithYFadO5f26h4cOX19Du33t2jnGfxs1UuHhoaJGjRwZgy9KTJpfrEB+WppYsyzS0yEqyhDg\nIyMNzTfXrhn+TUy03i2sSqWe6tV11KqVY3zUrKmjZs0c46NWLR1OTrZb45fviIm0qVcC8oE1qYhl\noddDXJySq1dVXLtmCPRRUXZcv254npRUMeYscHPTUaNGDjVq5P6bQ/XqOqpXN6RVr2547upa+YJ/\nRfxcWIu0qQvxgBQK8PbW4e2to0OHbLPX9HpISFBw/bohyEdF2XHjhorr1w3NObdvq8qlsxYgKUlJ\nUpKSfyzoZcHZWUf16jp8fAyBvlo1HdWq5eDjo6NaNR0+PrlpuhJ38oqKQ4K6EMWgUICnpx5Pz2xa\nt862eD07G/76S8Wff6q4eVPFzZuGoH/zpopbtwxBPyenfOcNSktTEhWlJCqq8G2dnAzB38dHZ7yw\n+fjk4OVlen7/w5abgCorCepClCJ7e4xDG/OSkwN37ii5fduO27eV3Lpl93ewV3L3rgM3blDmY/AL\nkp6u5MYNJTduFG17tdoQ3L28TA9PT9O/uQ/Dcz2enoYLgcyHV3YkqAtRjlQqqFNHR506lpNV5bYj\nZ2TAnTsq/vrL8DD8rTSm3b2rJDpahVZr/ciYkaHk1i0lt24VfR9HRz0eHoZg7+Fx/0Nv/LtmTTsc\nHR3w8NDj7q7D3V2Hm5tepmAuAgnqQlQwarVhcZL69fOu7YPhTtvYWCV37xqC/d27hlkwc/+Njjb8\ne++eiqws6wf/+2VmKrh715DngjmZPVMo9Li56f8O8Drc3c0D/v1/azSG193cdGg0hjQXl6pxUZCg\nLkQlpFRi7Nxs0UKb73a5HbsxMaZgf++ekpgYw9+xsUru3TME/9hYJdnZFesCcD+9XkFioqLEw02V\nSj0ajSHgazS5Ad/0b266q6vhIuDqakp3ddUbX6vo6+zafFDfvn078+bN4+jRo+WykIder6dfv34M\nGjQo37VOK7PU1FQWL17M2bNn+eOPP0hNTWXjxo0Wq0iBYRGVhQsXsmXLFpKSkggICOD111+nY8eO\neR575MiR1K5dm4ULFwKQk5PDY489xogRI5gwYQJgGBa7cuVKzpw5w/nz58nMzOT48ePUqVOn7N70\nP+zdu5eZM2fy888/4+rqWm7nLQlTx66WwhYn0+sNM2fGxCiJiTEE+ZgYJbGxuQ9DWlyc6VGRLwL/\npNM92EUhl4ODIbi7uur/fpj/7eJi+beLy/1/m567uOjJZ1nlErPpoJ6VlcV7773HSy+9VK4rM82Y\nMYOpU6cyYsQIvL29y+W85SUuLo7169cTGBhIUFAQu3btyncVqNdee439+/fz5ptv4ufnx5o1awgJ\nCWHbtm0Wk/wnJydz7Ngx/ve//xnTvv32WxISEnj22WeNaVFRUezYsYPAwEA6derEwYMHy+aNFiA4\nOJiPPvqIZcuWMWfOnHI/f1lRKPi7SSOHhg3zb/rJpdcb7ta9P8jHxiqJjzc84uKUJCSYnuc+Klpz\nUHFlZSmIi1MRF1c6x1OrdcZg7+xsHvSdnfVs3Fi849l0UN++fTsxMTEMGzasXM/76KOP4uHhwZdf\nfsn06dPL9dxlrW7dusYVrQ4dOsSuXbvy3C48PJywsDBCQ0ON5d+pUyd69uzJokWL+Pzzz822379/\nP/b29nTr1g0w/OJZtWoVI0aMwMHBtLRd586dOX36NADr1q0rVlB/6qmnqFu3LkuWLCn6G87H6NGj\nmTdvHi+//DLOzs4PfLzKSKEANzc9bm4Ft//fT6+HtDQF8fFKEhIUxqCfmGi4ABieK0hNdSAmRkdC\ngtJYu05Ntc0G8YwMJRkZEBNTOserlKUUFhZGUFAQDRs2pHfv3uzdu5ennnqKp556ymy7devWERwc\nbPal++abb/D19WXPnj3GtJycHJ588km6du1Kampqnuc8e/Ysvr6+7N271+K1V155hfbt23P/zbmD\nBw9m7dq1D/pWi2zMmDE89thjFul//vknvr6+fP311+WWFzA0Udjb2/P4448b01QqFYMHD+bgwYNk\nZ5uP8d69ezc9evQwBvCjR48SFRXFE088Ybbdg64NW9D+hw4dwtfXl08//dQsfdq0aQQEBPDXX38Z\n0wYMGEBmZiZhYWEPlJ+qRqEAFxc9vr45tGih5ZFHshg0KINRo9KYNi2FN95IYvHiRNauzWDTplh+\n/PEeJ05Ec/nyHaKibnP27B0OH77Ljh33WLs2llWr4pg/P4E5c5KYNi2Z0aNTefzxdHr2zKBduywa\nN86mRo0cnJx01n7r5abS1dQPHTrEtGnT6Nu3L//5z3+IiYnhrbfeIjMzkwYNGhi3S0lJ4eTJk7z1\n1ltm+z/zzDMcPHiQ1157jVatWlGzZk2WLl3KyZMnCQsLw8XFJc/zBgYG0rp1a7766iuCg4ON6YmJ\niWzfvp1p06aZBYwOHTqwYsUKLl26RNOmTQt8T1pt/h1duRQKBaoCGt/GjBnDs88+y+nTp81WoFq7\ndi2urq4MHTq00HOUpsuXL+Pn52fR7NW4cWOysrKIioqi8d+NvFlZWRw4cID33nvPuN2BAwfw9vam\nYcOG5ZbnoKAgJk+ezPvvv0+XLl146KGHWL9+PWFhYfzvf/+jVq1axm1dXV156KGHOHDgACEhIeWW\nx6rM3j73Tl+Aov0yuF92NiQnK0lKUuT7b2qqkuRkxd8PJSkphn+TkxWkpBj+Le+bx4qr0gX1RYsW\n0axZM7PaVLNmzejXr59ZUD937hxarZZmzZpZHOODDz6gT58+vPTSS0yfPp2PPvqImTNnFroc37PP\nPstrr73GrVu3jB1zmzZtQqvV8swzz5ht27x5cwDOnDlTYFA/evRokZqHOnfuzMYCGtd69uxJvXr1\n+Prrr43vIzs7m/Xr1/PEE0+UexNBQkIC7u7uFukef68MnZCQYEw7cuQI6enp9O7d25h2+vRpmjRp\n8kB5+OfFUq/Xo9PpyMnJMf6q+ufFctasWRw9epSpU6eybNky3nzzTUaNGkW/fv0sjt+sWTN+/vnn\nB8qjKD/29vx9gxSU5KIAhuajjAwFKSm5j9yAb7ogpKYa0lNSFKSlmf5OTTW9lppqeK0smpQqVVDP\nycnhjz/+4MUXXzRLb9myJX5+fmZp9+7dA8DT09PiOG5ubixfvpynnnqKUaNG0alTpyKNVBk8eDBv\nv/02a9euZebMmQB8/fXX9O7dm5o1a5pt62X45BjzkZ9WrVrl2y59v/x+QeRSKBSMGjWK0NBQ5s2b\nh0ajYffu3cTExDBq1KgC9/1n8LMr53lfd+/eTefOndFoNMa0mJgYY02+JG7cuJHniJzjx4+zadMm\n4/O6dety7Ngx43M7OzuWL19O3759GTx4MPXq1bP4tZfLy8uLmNJqCBWVgkIBTk56nJz0VKsGJb04\n5NLpDH0MhgCf+1CaPQfLGFaQShXU4+LiyM7OznNESXFHmbRt25YGDRoQERHBc889V6R9HB0dGT58\nOOvXr+fVV1/lt99+IyIignnz5hXr3PdzcXEx1uoLUpS25OHDh7No0SI2b97M2LFj+eqrr2jTpo3F\nSJN/ql+/vtnzTZs20alTp0LPVxB3d3du5XGbYW4NPbfGrtfr+eGHH3jppZfMtnvQtvOaNWuaXSz1\nej2zZs2iZs2a/Otf/zKm398Jm8vf35927dpx6NAhRo4cWW4jp0TVo1RiHA6ZPxsO6l5eXtjb2xMb\nG2vxWkxMDL6+vsbn1QyXUeLj4/M8VmhoKFFRUTz00EPMmzfPoqaYn9GjR/PJJ5+wZ88edu3ahZ+f\nHz169LDYLu7v8U65+chPaTW/gKF8BgwYwNdff01QUBDHjh1j0aJFhR77n78U7m/GKqmmTZuyZ88e\nMjIyzIJiREQEDg4OxgvJiRMniI6Opm/fvmb7+/j4GMuwJOzt7WnZsqVZmouLC56enhbp//Ttt99y\n6NAhWrVqRWhoKH379s1zHHxcXFyh/79ClLdKFdRVKhWBgYF8//33ZrWts2fPcuPGDbOg3qJFC1Qq\nFRcuXLCodR4/fpxly5Yxd+5cBg0aRHBwMHPmzGH58uWF5qF+/foEBQWxcuVKLly4kO+QxQsXLgCG\n5pWClFbzS64xY8YwePBgZsyYgZubG4MHDy50n8KCXEn06dOHxYsXs337dp5++mnA0Myzbds2unfv\njr29YSm477//3thhfb9WrVrxzTfflHq+CnPlyhXefPNNxo4dy+zZs+nTpw8vvvgimzZtQvmPe8wv\nXLhAYGBguedRiIJUqqAO8OqrrxISEsJzzz1HSEgIcXFxLFmyhOrVq5t96VxdXWnbti2//PIL48aN\nM6YnJCQwbdo0unXrxqRJkwBYsGABkyZNonv37sYAVJAxY8Ywfvx4HBwcLDpIc504cYIaNWoUOvLF\nxcWlVINqu3btaNGiBcePH2f8+PFl0nSwf/9+0tLSuHjxIgDHjh0jJiYGZ2dnevXqBRguqo8//jhv\nvfUWWq0WX19fvvzyS27evMmKFSuMx9qxY0eeI3OCgoJYtWoVly9ftugw3bFjB2C4mAPs27cPLy8v\nfHx8Cm02KmhNmKysLKZOnUq9evX497//jYODA8uXL+fJJ5/kww8/NLuAp6SkcOHCBUaOHFng+YQo\nb5UuqAcFBbF8+XJCQ0N5/vnn8ff3Z968eYSGhuLm5ma2bUhICHPmzCElJcV4O/fMmTPJyspi6dKl\nxu0GDhzIM888w5tvvsnDDz9MvXr1CsxDr169UKvV9O7dO8+2fL1ez3fffWe1oW79+/fn3LlzjB49\nukyOP3fuXG7evAkY2r4XL14MWHY6hoaGsmDBAj744AOSkpJ46KGHWLt2LS1atAAgMjKSyMhIi6YX\ngG7dulGvXj3CwsKMndK5ci/GueefO3cuULQmqoLa6ufPn09kZCQ7d+40trW3a9eO6dOns2TJEoKC\ngmjXrh1g+IXh4ODAkCFDCjyfEOXNJpazu337Nt26deOll17i5ZdfNqZnZWURFBTExIkTGT9+fKmd\n7277ssYAAAZLSURBVNChQ4SEhLB+/Xq6du1q8fr+/fuZPHkyR44cwcfHx+L1sl6qa/DgwdjZ2bF5\n8+YyO0dpWL58ORs3bsz3rtCvvvqKpUuXcvToURwr2CxKAwcOpHPnzrz++uuldkxZws1EysKkuMvZ\nVbo7SjMyMpg9ezbff/89x44dY/369TzzzDM4OTlZ1IwdHBx4/fXXWbFiBRkZGQ987qioKA4dOsRb\nb71FYGBgngEdDGPpX3rppTwDelnJysrit99+Y8mSJZw8edKsNltRTZs2jd9//z3f10eOHImnpydf\nffVVOeaqcHv37uXmzZtmFQghKopKV1PPzs5m8uTJ/P7778THx+Ps7EzHjh2ZPXv2A9+sUphXXnmF\nrVu3EhAQwIcffljicdRlUQvJHZft7u7OmDFjLJosKiqpkZlIWZhIWZgUt6Ze6YK6LZAPrImUhYmU\nhYmUhUlxg3qJO0pPnz7NmjVr0Ol09OrVSzqMhBCiAihRm7pOp+P//u//mDt3LqGhoRw5csQ4GkII\nIYT1lCioR0ZGUrNmTapXr46dnR1du3blt99+K+28CSGEKKYSBfW4uDiz8dleXl4PdEu3EEKI0lHp\nhjQKIYTIX4k6Sr28vMwm1YqNjTVONZsrPDzcuOwZwLBhw4rdi2vLijJ5WFUhZWEiZWEiZWGyYcMG\n498BAQEFz7yqLwGtVqufNm2a/u7du/rs7Gz9a6+9pr9x40aB+6xfv74kp7JJUhYmUhYmUhYmUhYm\nxS2LEtXUVSoV48eP59133zUOabx/hkQhhBDWUeJx6m3atKFNmzalmRchhBAPqNw6SgtbfacqkbIw\nkbIwkbIwkbIwKW5ZlOs0AUIIIcqWDGkUQggbIkFdCCFsSLmsfCSTfxnExMSwYsUKEhMTUSgUPPro\no/Tv39/a2bIanU7H7Nmz8fLyYvbs2dbOjtWkpqayevVq4/xJkydPLvNppCuqrVu3cvjwYRQKBX5+\nfkyZMsW4nq2tW7lyJadOncLNzc24mlhKSgpLliwhJiaGatWqMX369ELXKy7zmrpM/mViZ2fHmDFj\nCA0N5d1332XPnj1VtiwAdu7cia+vb4FLzFUFn3/+OW3atGHJkiUsWrSoyg4Pjo6OZt++fSxYsIDF\nixej0+k4cuSItbNVbnr27GlcmjFXWFgYgYGBfPjhh7Ro0YKwsLBCj1PmQV0m/zLx8PCgfv36AKjV\naurUqUN8fLx1M2UlsbGxnDp1il69ehW4GLSty13AO3fBbpVKhbOzs5VzZR3Ozs6oVCoyMzPJyckh\nMzPT4k51W9a8eXOLWvhvv/1G9+7dAejRowcnTpwo9Dhl3vyS1+RfkZGRZX3aCi86OpqoqKgSr55U\n2X3xxReMGjWK9PR0a2fFqqKjo3Fzc2PlypVcv34df39/xo0bV+HWZC0Prq6uDBo0iClTpuDg4ECr\nVq0IDAy0drasKjExEQ8PDwDc3d1JTEwsdB/pKLWCjIwMQkNDGTt2LGq12trZKXcnT57Ezc0Nf3//\nKl1LB8jJyeHatWsEBwezYMEC1Gp1kX5i26I7d+7w/fffs2LFCj7++GMyMjI4fPiwtbNVYRS1mbLM\ng3pRJv+qSrRaLYsXL6Zbt248/PDD1s6OVVy6dImTJ08ydepUPvzwQ8LDw1m+fLm1s2UV3t7eeHl5\n0ahRIwA6derEtWvXrJwr67h69SpNmzZFo9GgUqno2LEjly5dsna2rMrd3Z2EhAQA4uPjcXd3L3Sf\nMg/qDRs25M6dO0RHR6PVajl69Cjt27cv69NWSHq9ntWrV1OnTh0GDBhg7exYTUhICKtWrWLFihW8\n8sorBAQEMG3aNGtnyyo8PDzw8fExrt979uzZKttRWrt2bSIiIsjKykKv11fpssjVvn17Dhw4AMDB\ngwfp0KFDofuUyx2lp06dMhvS+MQTT5T1KSukixcvMm/ePPz8/Iw/pUJCQmjdurWVc2Y958+fZ/v2\n7cyaNcvaWbGaqKgoPv74Y7RaLTVq1GDKlClVtrP0u+++4+DBgygUCvz9/Zk0aRJ2duUy8trqli5d\nyoULF0hKSsLDw4Nhw4bRoUOHYg9plGkChBDChkhHqRBC2BAJ6kIIYUMkqAshhA2RoC6EEDZEgroQ\nQtgQCepCCGFDJKgLIYQNkaAuhBA25P8BIoAZtAZRA0sAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x114266810>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercise 2 (10 points)**.\n",
      "\n",
      "You are given 5 samples of 100 rolls of a die. The die is chosen at random from one of two dice, either of which may be biased. Given the following counts, estimate the parameters for each die (i.e. the probabilites of getting a 1, 2, 3, 4, 5 or 6)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df = pd.read_csv('dice.csv', index_col=0)\n",
      "df"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>1</th>\n",
        "      <th>2</th>\n",
        "      <th>3</th>\n",
        "      <th>4</th>\n",
        "      <th>5</th>\n",
        "      <th>6</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>Sample 1</th>\n",
        "      <td>  2</td>\n",
        "      <td> 14</td>\n",
        "      <td> 10</td>\n",
        "      <td> 20</td>\n",
        "      <td> 26</td>\n",
        "      <td> 28</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>Sample 2</th>\n",
        "      <td>  1</td>\n",
        "      <td> 13</td>\n",
        "      <td> 12</td>\n",
        "      <td> 17</td>\n",
        "      <td> 21</td>\n",
        "      <td> 36</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>Sample 3</th>\n",
        "      <td> 21</td>\n",
        "      <td> 18</td>\n",
        "      <td> 16</td>\n",
        "      <td> 12</td>\n",
        "      <td> 14</td>\n",
        "      <td> 19</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>Sample 4</th>\n",
        "      <td> 15</td>\n",
        "      <td> 20</td>\n",
        "      <td> 20</td>\n",
        "      <td> 15</td>\n",
        "      <td> 12</td>\n",
        "      <td> 18</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>Sample 5</th>\n",
        "      <td>  5</td>\n",
        "      <td> 10</td>\n",
        "      <td> 11</td>\n",
        "      <td> 24</td>\n",
        "      <td> 17</td>\n",
        "      <td> 33</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "           1   2   3   4   5   6\n",
        "Sample 1   2  14  10  20  26  28\n",
        "Sample 2   1  13  12  17  21  36\n",
        "Sample 3  21  18  16  12  14  19\n",
        "Sample 4  15  20  20  15  12  18\n",
        "Sample 5   5  10  11  24  17  33"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def em(xs, thetas, max_iter=100, tol=1e-6):\n",
      "    \"\"\"Expectation-maximization for multinomial data with latent class.\"\"\"\n",
      "    \n",
      "    ll_old = -np.infty\n",
      "    for i in range(max_iter):\n",
      "        ll = np.array([np.sum(xs * np.log(theta), axis=1) \n",
      "                       for theta in thetas])    \n",
      "        lik = np.exp(ll)\n",
      "        ws = lik/lik.sum(0)\n",
      "        exps = np.array([w[:, None] * xs for w in ws])\n",
      "        thetas = np.array([expr.sum(0)/expr.sum() for expr in exps])\n",
      "        ll_new = np.sum([w*l for w, l in zip(ws, ll)])\n",
      "        if np.abs(ll_new - ll_old) < tol:\n",
      "            break\n",
      "        ll_old = ll_new\n",
      "    return i, thetas, ll_new"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xs = df.values\n",
      "\n",
      "thetas = np.random.random((2,6))\n",
      "thetas /= thetas.sum(1)[:, None]\n",
      "\n",
      "i, thetas, ll = em(xs, thetas)\n",
      "thetas"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "array([[ 0.18  ,  0.19  ,  0.18  ,  0.135 ,  0.13  ,  0.185 ],\n",
        "       [ 0.0267,  0.1233,  0.11  ,  0.2033,  0.2133,  0.3233]])"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercise 3 (10 points)**\n",
      "\n",
      "You have some data about the intervals between arrivals (in minutes) of customers to a bank. You want to know how many customers to expect from 9 AM to 5 PM. Simulate 1000 sequences of customer arrivals over an 8 hour period and esitmate the mean and variance. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arrivals = np.concatenate([np.random.exponential(3, 175), \n",
      "           np.random.normal(25, 3, 100),\n",
      "           np.random.exponential(60, 25)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.savetxt('arrival_intervals.txt', arrivals)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.loadtxt('arrival_intervals.txt')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "array([  2.0952e-01,   2.4595e-01,   6.5462e-01,   3.9017e-01,\n",
        "         2.8626e+00,   1.4996e+00,   5.5581e+00,   1.6387e+00,\n",
        "         6.4469e+00,   4.2354e-01,   1.0948e+00,   4.3440e-01,\n",
        "         8.0157e+00,   2.2653e+00,   1.4452e+00,   3.2204e-01,\n",
        "         3.7704e+00,   7.9621e+00,   7.9283e-01,   5.3954e-01,\n",
        "         2.7608e+00,   1.3300e+00,   3.7619e+00,   4.6724e-01,\n",
        "         6.8437e+00,   1.5951e-01,   9.5755e-01,   3.8534e+00,\n",
        "         1.0908e+00,   8.6424e-01,   5.2758e+00,   3.4553e+00,\n",
        "         1.0507e+00,   1.6807e+00,   8.7264e-01,   1.0299e+00,\n",
        "         3.8108e+00,   5.4847e+00,   3.7906e+00,   1.4270e-01,\n",
        "         3.9698e+00,   6.0365e+00,   5.6018e+00,   7.0588e+00,\n",
        "         9.2178e-01,   4.9560e+00,   7.6875e-01,   1.0845e+00,\n",
        "         3.3594e-01,   6.8316e+00,   7.2840e+00,   3.6679e-02,\n",
        "         7.6652e-01,   6.9307e-01,   5.9371e+00,   7.0195e+00,\n",
        "         3.4020e+00,   6.3075e-01,   3.0068e-01,   2.5030e+00,\n",
        "         5.7195e-01,   6.0101e+00,   2.2153e+00,   3.2337e-01,\n",
        "         2.6913e+00,   1.4517e+00,   1.8445e+00,   2.8519e+00,\n",
        "         1.5494e+00,   1.7431e+00,   4.8032e-01,   1.2604e+00,\n",
        "         3.6388e+00,   1.3269e+00,   1.0399e+00,   6.0045e+00,\n",
        "         1.1249e+00,   1.7349e+00,   4.6314e+00,   1.1772e+00,\n",
        "         1.7972e+00,   1.3460e+00,   7.0102e+00,   4.4128e+00,\n",
        "         1.2411e+00,   2.3825e+00,   3.8255e+00,   1.8780e+00,\n",
        "         9.1881e-01,   1.0481e+01,   3.5568e+00,   2.9665e-01,\n",
        "         2.7583e+01,   9.1318e-01,   2.1824e-01,   7.7931e-01,\n",
        "         9.6784e+00,   1.9411e+00,   4.3931e+00,   4.2959e+00,\n",
        "         1.3618e+00,   4.1625e+00,   1.9358e+00,   5.5990e+00,\n",
        "         5.1171e-01,   8.9539e+00,   1.8964e-01,   4.9843e+00,\n",
        "         1.0700e-01,   7.2117e-01,   5.6948e-01,   1.5690e+00,\n",
        "         2.3711e-01,   3.3224e+00,   7.2713e+00,   1.8412e+00,\n",
        "         1.4490e+00,   6.7471e-01,   2.4715e+00,   1.8328e+00,\n",
        "         4.0518e-01,   2.0578e+00,   4.9456e+00,   1.1607e+00,\n",
        "         1.2465e+01,   7.9455e-01,   2.3402e-01,   4.2066e-02,\n",
        "         5.6096e-01,   1.1255e+00,   1.4611e+00,   4.3470e+00,\n",
        "         3.7590e+00,   2.3997e+00,   4.9479e+00,   5.8804e-02,\n",
        "         5.6593e+00,   2.8572e+00,   8.4518e+00,   3.0810e+00,\n",
        "         5.6536e+00,   4.2691e+00,   1.7056e+00,   5.1083e-01,\n",
        "         1.4994e-01,   1.8974e-02,   3.7252e+00,   6.7738e-01,\n",
        "         1.4886e+00,   9.6532e-01,   1.4823e+00,   6.6099e+00,\n",
        "         1.2337e-01,   4.0870e+00,   1.9671e+00,   3.9297e-01,\n",
        "         2.4114e+00,   1.8343e-01,   1.1479e+01,   1.3296e+00,\n",
        "         6.2945e+00,   1.4959e+00,   3.7373e+00,   3.2958e+00,\n",
        "         1.1885e-01,   9.9590e-01,   7.8460e-01,   5.7697e-01,\n",
        "         2.7764e+00,   7.1881e-02,   1.1787e+00,   2.7527e-01,\n",
        "         6.1433e+00,   3.2017e+00,   2.7623e-01,   2.6197e+01,\n",
        "         2.4764e+01,   3.0039e+01,   2.2057e+01,   2.3315e+01,\n",
        "         1.9768e+01,   2.7899e+01,   2.4702e+01,   2.5727e+01,\n",
        "         2.8652e+01,   2.8336e+01,   2.1437e+01,   2.3371e+01,\n",
        "         2.5870e+01,   2.6159e+01,   2.1821e+01,   2.4942e+01,\n",
        "         2.0771e+01,   2.9257e+01,   2.3980e+01,   2.4899e+01,\n",
        "         3.0414e+01,   2.5214e+01,   2.0920e+01,   2.9176e+01,\n",
        "         2.6557e+01,   2.6493e+01,   2.8092e+01,   2.7286e+01,\n",
        "         2.0833e+01,   3.1187e+01,   2.3857e+01,   2.1259e+01,\n",
        "         3.1332e+01,   2.6869e+01,   2.1290e+01,   2.5214e+01,\n",
        "         2.8245e+01,   2.5613e+01,   2.4463e+01,   2.1112e+01,\n",
        "         2.2957e+01,   2.2499e+01,   2.4447e+01,   2.5072e+01,\n",
        "         2.7396e+01,   2.0077e+01,   2.2116e+01,   1.4658e+01,\n",
        "         2.4592e+01,   2.8681e+01,   2.0961e+01,   2.9223e+01,\n",
        "         2.3651e+01,   2.6711e+01,   2.4583e+01,   2.4088e+01,\n",
        "         2.7668e+01,   1.7236e+01,   2.3593e+01,   2.4987e+01,\n",
        "         2.3347e+01,   2.8019e+01,   2.6786e+01,   2.1396e+01,\n",
        "         2.5403e+01,   2.3843e+01,   1.9855e+01,   3.0139e+01,\n",
        "         2.7830e+01,   2.1973e+01,   2.8394e+01,   2.2361e+01,\n",
        "         2.5312e+01,   2.5856e+01,   2.1760e+01,   2.4431e+01,\n",
        "         2.5879e+01,   2.7804e+01,   2.3485e+01,   2.5750e+01,\n",
        "         2.4964e+01,   2.2832e+01,   2.2027e+01,   2.4498e+01,\n",
        "         2.6071e+01,   2.1524e+01,   2.7493e+01,   2.6689e+01,\n",
        "         2.5433e+01,   1.6391e+01,   2.6128e+01,   2.7973e+01,\n",
        "         2.5907e+01,   2.6497e+01,   2.8025e+01,   2.7738e+01,\n",
        "         2.7909e+01,   2.4498e+01,   2.7469e+01,   4.6460e+01,\n",
        "         9.0623e+00,   1.6338e+02,   2.1691e+01,   1.2971e+00,\n",
        "         1.7102e+01,   1.6916e+02,   8.3688e+01,   2.0345e+02,\n",
        "         3.0059e+01,   2.2643e+01,   5.9453e+01,   1.5119e+02,\n",
        "         1.3375e+01,   6.8167e+01,   3.3235e+01,   1.1708e+01,\n",
        "         1.6652e+01,   1.4009e+01,   1.5237e+01,   3.3904e+01,\n",
        "         7.5804e+00,   9.7186e+01,   4.8057e+01,   2.2868e+01])"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist(arrivals, bins=25, histtype='step');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEECAYAAADAoTRlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG6pJREFUeJzt3X9sU+f9L/D3cVxDCCc2dklDyDJQQ7dgbjK2hBalI4Tu\nO2kta8PVenob2C0pG1JS1i2jbNxsbFeiqOpQfqkjMA3uusI64WnYotoPoa1JWiVXmtOlo5i0VW7L\nvkT5hiy24ziQBByf+0eFz9KEcOIcY+Ln/ZIq2cfP4/Pxh9N3jo+PjyVVVVUQEZFQTMkugIiI7jyG\nPxGRgBj+REQCYvgTEQmI4U9EJCCGPxGRgMyzPdjS0oLu7m5kZmaivr4eANDb24sTJ05gcnISaWlp\n2LVrF/Lz8wEAbrcbra2tMJlMqKqqQlFRUeJfARERzdmse/7l5eWoq6ubsuzUqVN46qmn8LOf/QyK\nouA3v/kNAKCvrw+dnZ1oaGhAXV0djh8/jmg0etsCfD7fPMpPLeyFhr3QsBca9kIz317MGv4FBQXI\nyMiYssxms+HatWsAgKtXr2LZsmUAAK/Xi9LSUpjNZmRlZSE7Oxu9vb23LYD/mBr2QsNeaNgLDXuh\nmW8vZj3sM5Pt27fjJz/5CU6ePAlVVfHiiy8CAILBINasWRMb53A4EAgE5lUcERElxpw/8D127Biq\nqqpw9OhRPPPMMzh69Ogtx0qSNK/iiIgoMea859/b24sDBw4AAB566CEcO3YMAGC32+H3+2Pj/H4/\n7Hb7tPk+n2/K2xVFUeZcdKpiLzTshYa90LAXGkVR4HK5YvedTiecTqfu+XMO/+zsbFy8eBFr167F\nhQsXkJOTAwAoLi5Gc3Mztm7dikAggIGBgdhZQP9upgL7+/vnWkZKkmUZ4XA42WXcFdgLDXuhYS80\nOTk58/pjOGv4NzU1oaenByMjI6iuroaiKNi9ezdOnDiBGzduwGKxYPfu3QCA3NxcbNy4EbW1tbFT\nQHnYh4jo7iTdDZd05p7/J7hXo2EvNOyFhr3Q3DzqEi9+w5eISEAMfyIiATH8iYgExPAnIhIQw5+I\nSEAMfyIiATH8iYgENOdv+CaC+sF7+gZmLoO0IjexxRARCWBB7fmrnX9NdglERCnhrtjzlz7333SN\n0/0OgYiIZrWg9vyJiMgYDH8iIgEx/ImIBMTwJyISEMOfiEhADH8iIgHNeqpnS0sLuru7kZmZifr6\n+tjyP/3pTzh37hxMJhPWr1+PHTt2AADcbjdaW1thMplQVVWFoqKixFZPRERxmTX8y8vL8bWvfQ0/\n//nPY8suXLiArq4uHD58GGazGSMjIwCAvr4+dHZ2oqGhAYFAAAcPHkRzczNMJr65ICK628yazAUF\nBcjIyJiy7Ny5c9i2bRvM5k/+bmRmZgIAvF4vSktLYTabkZWVhezsbPT29iaobCIimo85f8N3YGAA\nFy9exG9/+1vcc889+OY3v4n7778fwWAQa9asiY1zOBwIBAKGFktERMaYc/hPTk7i6tWrOHToEHp7\ne9HY2DjlsNC/kyRp2jKfzwefzxe7rygKZFnWte4xyyKk6xy7EFksFt29SHXshYa90LAXU7lcrtht\np9MJp9Ope+6cw9/hcODBBx8EAOTn50OSJIyMjMBut8Pv98fG+f1+2O32afNnKjAcDutad/T6BCI6\nxy5Esizr7kWqYy807IWGvdDIsgxFUeKeP+dPY0tKSnDhwgUAQH9/PyKRCDIzM1FcXIyOjg5EIhEM\nDg5iYGAA+fn5cRdGRESJM+uef1NTE3p6ehAOh1FdXQ1FUVBeXo6jR49i7969MJvN2LNnDwAgNzcX\nGzduRG1tLdLS0rBr164ZD/sQEVHySaqqqskuor+/X9e46NnXYXq8MsHVJA/f0mrYCw17oWEvNDk5\nOfOaz5PwiYgExPAnIhIQw5+ISEAMfyIiATH8iYgExPAnIhIQw5+ISEAMfyIiATH8iYgExPAnIhIQ\nw5+ISEAMfyIiATH8iYgExPAnIhIQw5+ISEAMfyIiAc0a/i0tLfj2t7+NvXv3TnvsjTfewFNPPYXR\n0dHYMrfbjeeffx7f+9738I9//MP4aomIyBCzhn95eTnq6uqmLR8aGsL58+dx7733xpb19fWhs7MT\nDQ0NqKurw/HjxxGNRo2vmIiI5m3W8C8oKEBGRsa05a+99hp27NgxZZnX60VpaSnMZjOysrKQnZ2N\n3t5eY6slIiJDzPmYv9frhd1ux2c/+9kpy4PBIBwOR+y+w+FAIBCYf4VERGQ481wGT0xMwO1248c/\n/nFs2Wy//y5J0rRlPp8PPp8vdl9RFMiyrGv9Y5ZFSNc5diGyWCy6e5Hq2AsNe6FhL6ZyuVyx206n\nE06nU/fcOYX/lStX8K9//Qv79u0DAAQCAezfvx+HDh2C3W6H3++PjfX7/bDb7dOeY6YCw+GwrvVH\nr08gonPsQiTLsu5epDr2QsNeaNgLjSzLUBQl7vlzCv+8vDz88pe/jN1/7rnn8PLLL2Pp0qUoLi5G\nc3Mztm7dikAggIGBAeTn58ddGBERJc6s4d/U1ISenh6Ew2FUV1dDURSUl5fHHv/3wzq5ubnYuHEj\namtrkZaWhl27ds142IeIiJJPUmc7aH+H9Pf36xoXPfs6TI9XJria5OFbWg17oWEvNOyFJicnZ17z\n+Q1fIiIBMfyJiATE8CciEhDDn4hIQAx/IiIBMfyJiATE8CciEhDDn4hIQAx/IiIBMfyJiATE8Cci\nEhDDn4hIQAx/IiIBMfyJiATE8CciEhDDn4hIQLP+kldLSwu6u7uRmZmJ+vp6AMDJkyfx97//HWaz\nGffddx9qamqwZMkSAIDb7UZraytMJhOqqqpQVFSU+FdARERzNuuef3l5Oerq6qYsKyoqQn19PQ4f\nPowVK1bA7XYDAPr6+tDZ2YmGhgbU1dXh+PHjiEajiauciIjiNmv4FxQUICMjY8qywsJCmEyfTFuz\nZg38fj8AwOv1orS0FGazGVlZWcjOzkZvb2+CyiYiovmY1zH/N998E1/84hcBAMFgEA6HI/aYw+FA\nIBCYX3VERJQQsx7zn82ZM2dgNpvx8MMP33KMJEnTlvl8Pvh8vth9RVEgy7KudY5ZFiFd59iFyGKx\n6O5FqmMvNOyFhr2YyuVyxW47nU44nU7dc+MK/7a2NnR3d+PAgQOxZXa7PXYICAD8fj/sdvu0uTMV\nGA6Hda03en0CEZ1jFyJZlnX3ItWxFxr2QsNeaGRZhqIocc+f82Gfd999F2fPnsW+fftgsVhiy4uL\ni9HR0YFIJILBwUEMDAwgPz8/7sKIiChxZt3zb2pqQk9PD0ZGRlBdXY0nn3wSHo8HkUgEL774IgDg\ngQcewLe+9S3k5uZi48aNqK2tRVpaGnbt2jXjYR8iIko+SVVVNdlF9Pf36xoXPfs6TI9XJria5OFb\nWg17oWEvNOyFJicnZ17z+Q1fIiIBMfyJiATE8CciEhDDn4hIQAx/IiIBMfyJiATE8CciEhDDn4hI\nQAx/IiIBMfyJiATE8CciEhDDn4hIQAx/IiIBMfyJiATE8CciEtCsP+bS0tKC7u5uZGZmor6+HgAw\nOjqKxsZGDA0NYfny5aitrUVGRgYAwO12o7W1FSaTCVVVVSgqKkr8KyAiojmbdc+/vLwcdXV1U5Z5\nPB4UFhaiubkZ69atg8fjAQD09fWhs7MTDQ0NqKurw/HjxxGNRhNXORERxW3W8C8oKIjt1d/U1dWF\nsrIyAMDmzZvh9XoBAF6vF6WlpTCbzcjKykJ2djZ6e3sTVDYREc3HnI/5h0Ih2Gw2AIDVakUoFAIA\nBINBOByO2DiHw4FAIGBQmUREZKR5feB7ux9o5w+4ExHdnWb9wHcmVqsVw8PDsNlsCAaDsFqtAAC7\n3Q6/3x8b5/f7Ybfbp833+Xzw+Xyx+4qiQJZlXesesyxCus6xC5HFYtHdi1THXmjYCw17MZXL5Yrd\ndjqdcDqduufOOfyLi4vR1taGiooKtLe3o6SkJLa8ubkZW7duRSAQwMDAAPLz86fNn6nAcDisa93R\n6xOI6By7EMmyrLsXqY690LAXGvZCI8syFEWJe/6s4d/U1ISenh6MjIyguroaiqKgoqICjY2NaG1t\njZ3qCQC5ubnYuHEjamtrkZaWhl27dvGwDxHRXUpSVVVNdhH9/f26xkXPvg7T45UJriZ5uFejYS80\n7IWGvdDk5OTMaz6/4UtEJCCGPxGRgBj+REQCYvgTEQmI4U9EJCCGPxGRgBj+REQCYvgTEQmI4U9E\nJCCGPxGRgBj+REQCYvgTEQmI4U9EJCCGPxGRgBj+REQCYvgTEQlozj/jeJPb7cbbb78NSZKQl5eH\nmpoaTExMoLGxEUNDQ7Ff+crIyDCyXiIiMkBce/6Dg4P461//ipdffhn19fWIRqPo6OiAx+NBYWEh\nmpubsW7dOng8HqPrJSIiA8QV/kuWLEFaWhomJiYwOTmJiYkJ2O12dHV1oaysDACwefNmeL1eQ4sl\nIiJjxHXYZ+nSpfj617+OmpoaWCwWFBUVobCwEKFQCDabDQBgtVoRCoUMLZaIiIwR157/wMAA/vCH\nP+DIkSP4xS9+gfHxcbz11ltTxkiSZEiBRERkvLj2/D/66CN87nOfgyzLAIAHH3wQH374IWw2G4aH\nh2Gz2RAMBmG1WqfN9fl88Pl8sfuKosSe53bGLIuQrnPsQmSxWHT3ItWxFxr2QsNeTOVyuWK3nU4n\nnE6n7rlxhX9OTg5+//vf4/r167jnnntw/vx55OfnY/HixWhra0NFRQXa29tRUlIybe5MBYbDYV3r\njV6fQETn2IVIlmXdvUh17IWGvdCwFxpZlqEoStzz4wr/VatWYdOmTdi/fz8kScLq1avxla98BePj\n42hsbERra2vsVE8iIrr7SKqqqskuor+/X9e46NnXYXq8MsHVJA/3ajTshYa90LAXmpycnHnN5zd8\niYgExPAnIhIQw5+ISEAMfyIiATH8iYgExPAnIhIQw5+ISEAMfyIiATH8iYgExPAnIhIQw5+ISEAM\nfyIiATH8iYgEFNclnVOJGo0C18f1TzBbIJmFbxsRLXBMsathRE8egZRfcPuxwwFI+WuBL25MfF1E\nRAnE8AcgFRTBVP7Ybceplz8G/jVwByoiIkqsuMP/6tWrOHbsGPr6+gAANTU1WLFiBRobGzE0NBT7\nJa+MjAzDiiUiImPEHf6/+tWvsH79euzduxeTk5OYmJjAmTNnUFhYiCeeeAIejwcejwfbt283sl4i\nIjJAXGf7XLt2De+//z62bNkCAEhLS8OSJUvQ1dWFsrIyAMDmzZvh9XqNq5SIiAwT157/4OAgMjMz\n0dLSgn/+859YvXo1du7ciVAoBJvNBgCwWq0IhUKGFktERMaIK/wnJyfx8ccf49lnn0V+fj5effVV\neDyeKWMkSZpxrs/ng8/ni91XFAWyLOta75hlEdJ1jtUrqk7ixqJFWKTjeSczlmByNB0Wg2u4yWKx\n6O5FqmMvNOyFhr2YyuVyxW47nU44nU7dc+MKf4fDAbvdjvz8fADAQw89BLfbDZvNhuHhYdhsNgSD\nQVit1mlzZyowHA7rWm/0+gQiOsfqpY6OQp2YwHUdz6tevQaMjWHC4BpukmVZdy9SHXuhYS807IVG\nlmUoihL3/LiO+dtsNtx7773o7+8HAJw/fx6f+cxn8KUvfQltbW0AgPb2dpSUlMRdGBERJU7cZ/tU\nVVXhlVdeQSQSwX333YeamhpEo1E0NjaitbU1dqonERHdfeIO/1WrVuGll16atvzAgQPzKoiIiBKP\nF3YjIhIQw5+ISEAMfyIiATH8iYgExPAnIhIQw5+ISEAMfyIiATH8iYgExPAnIhIQw5+ISEAMfyIi\nATH8iYgExPAnIhIQw5+ISEAMfyIiATH8iYgEFPePuQBANBrF/v37YbfbsX//foyOjqKxsRFDQ0Ox\nX/LKyMgwqlYiIjLIvPb8//jHPyI3NxeSJAEAPB4PCgsL0dzcjHXr1sHj8RhSJBERGSvu8Pf7/eju\n7saWLVugqioAoKurC2VlZQCAzZs3w+v1GlNlHNT33oH6Xtft/7v4btJqJCJKlrgP+/z617/Gjh07\nMDY2FlsWCoVgs9kAAFarFaFQaP4Vxknt/r+QvvzV246TlmYCjuV3oCIiortHXOH/zjvvIDMzE6tX\nr4bP55txzM1DQZ/m8/mmzFEUBbIs61rvmGUR0vWOzVqB9MIv6Rqr12TGEkyOpsOis4a5slgsunuR\n6tgLDXuhYS+mcrlcsdtOpxNOp1P33LjC/4MPPsA777yD7u5u3LhxA2NjY3jllVdgtVoxPDwMm82G\nYDAIq9U6be5MBYbDYV3rjV6fQCQBY/VSr14DxsYwYfDz3iTLsu5epDr2QsNeaNgLjSzLUBQl7vlx\nhX9lZSUqKysBABcvXsTZs2fxne98B6dOnUJbWxsqKirQ3t6OkpKSuAsjIqLEMeQ8/5uHeCoqKvDe\ne+/hu9/9Li5cuICKigojnp6IiAw2r/P8AWDt2rVYu3YtAGDp0qU4cODAvIsiIqLE4jd8iYgExPAn\nIhIQw5+ISEAMfyIiATH8iYgExPAnIhIQw5+ISEAMfyIiATH8iYgExPAnIhIQw5+ISEAMfyIiATH8\niYgExPAnIhIQw5+ISEAMfyIiAcX1Yy5DQ0M4cuQIQqEQJEnCI488gkcffRSjo6NobGzE0NAQli9f\njtraWmRkZBhdMxERzVNc4W82m/HMM89g1apVGB8fxw9/+EMUFhaira0NhYWFeOKJJ+DxeODxeLB9\n+3ajayYionmK67CPzWbDqlWrAACLFy/GypUrEQgE0NXVhbKyMgDA5s2b4fV6DSuUiIiMM+9j/oOD\ng7h06RLWrFmDUCgEm80GALBarQiFQvMukIiIjDevH3AfHx9HfX09du7cifT09CmPSZI04xyfzwef\nzxe7rygKZFnWtb4xyyKkJ2CsXpMZSzA5mg6Lwc97k8Vi0d2LVMdeaNgLDXsxlcvlit12Op1wOp26\n58Yd/pFIBPX19di0aRM2bNgA4JO9/eHhYdhsNgSDQVit1mnzZiowHA7rWmf0+gQiCRirl3r1GjA2\nhgmDn/cmWZZ19yLVsRca9kLDXmhkWYaiKHHPj+uwj6qqOHbsGFauXInHHnsstry4uBhtbW0AgPb2\ndpSUlMRdGBERJU5ce/4ffPAB3n77beTl5eEHP/gBAKCyshIVFRVobGxEa2tr7FRPIiK6+8QV/p//\n/Odx+vTpGR87cODAvAoiIqLE4zd8iYgExPAnIhIQw5+ISEAMfyIiAc3rS1533NAg1Hc69I1V1cTW\nQkS0gC2o8Jce/QYwOalv7Jq1Ca6GiGjhWljhn52b7BKIiFICj/kTEQmI4U9EJCCGPxGRgBj+REQC\nYvgTEQmI4U9EJCCGPxGRgBj+REQCMvxLXu+++y5effVVRKNRbNmyBRUVFUavgoiI5snQPf9oNIoT\nJ06grq4ODQ0N6OjoQF9fn5GrICIiAxi659/b24vs7GxkZWUBAEpLS9HV1YXc3BS5LMPSTKjet6H2\nfWzo06r9/wnTU98GZNnQ5yW6Kfp/moB7s/QNVgHTE5WJLYiSztDwDwQCcDgcsft2ux29vb1GriKp\npGUOSP/9fxr+vNG3z+m+Cmn0b28B/9UHSIaXAdPj+v6Hj579LQCDr5qqAliRC9OGTcY+7wITbfsT\nMBKcsmzMsgjR6xPTB5tMMG39H/qe+N6sOfz7vq7vORNE/a/LUN86B6SnT3tsWi+uXYX0HxWQHMvv\nYIWpYUFd2C2VqW/9GePLHIhOjM8+zvcuTN/9KSSTsZ/VR0+1IHrOrW9wKAjTN2sMXb+qqog2/W9E\nh/0AgPFFi2/bi1Sk/r/3kVb9v6YsS5dlRMLhaWOjr/1c/79ZeER/Df/5kf7nTYShK5C+8CCkteun\nPfTpXqgXu4HJyJ2sLmVIqmrche8//PBD/O53v8OPfvQjAIDb7YYkSVM+9PX5fPD5fLH7iqIYtXoi\nIqG4XK7YbafTCafTqX+yaqBIJKLu2bNHvXLlinrjxg31hRdeUC9fvjzrnNOnTxtZwoLGXmjYCw17\noWEvNPPthaGHfdLS0vDss8/i0KFDsVM9U+bDXiKiFGL4Mf/169dj/frpx+qIiOjukfRv+M7pGFWK\nYy807IWGvdCwF5r59sLQD3yJiGhhSPqePxER3XkMfyIiASX1S16iXwTuueeeQ3p6OkwmE9LS0vDS\nSy9hdHQUjY2NGBoawvLly1FbW4uMjIxkl2q4lpYWdHd3IzMzE/X19QAw62t3u91obW2FyWRCVVUV\nioqKklm+oWbqhcvlwptvvonMzEwAwNNPPx07kSJVezE0NIQjR44gFApBkiQ88sgjePTRR4XcLm7V\nC0O3C0NOOI3D5OTknL8TkGpqamrUcDg8ZdnJkydVj8ejqqqqut1u9dSpU8koLeEuXryofvTRR+r3\nv//92LJbvfbLly+rL7zwgnrjxg31ypUr6p49e9TJycmk1J0IM/XC5XKpb7zxxrSxqdyLYDCofvzx\nx6qqqurY2Jj6/PPPq5cvXxZyu7hVL4zcLpJ22OffLwJnNptjF4ETjfqpz9u7urpQVlYGANi8eTO8\nXm8yykq4goKCae9obvXavV4vSktLYTabkZWVhezs7JS6ZtRMvQCmbxtAavfCZrNh1apVAIDFixdj\n5cqVCAQCQm4Xt+oFYNx2kbTwn+kicDdfnCgkScLBgwexf/9+/OUvfwEAhEIh2Gw2AIDVakUoFEpm\niXfUrV57MBicsq04HA4htpU///nP2LdvH44ePYqrV68CEKcXg4ODuHTpEtasWSP8dnGzFw888AAA\n47YLXtgtiQ4ePIhly5ZhZGQEBw8exMqVK6c8LkkJuHTnAnG7157qvfnqV7+Kb3zjGwCA06dP47XX\nXkN1dfWMY1OtF+Pj46ivr8fOnTuR/qkre4q2XYyPj6OhoQE7d+7E4sWLDd0ukrbnb7fb4ff7Y/f9\nfj/sdnuyykmKZcuWAQAyMzOxYcMG9Pb2wmq1Ynh4GMAnf82tVmsyS7yjbvXaRdxWrFYrJEmCJEnY\nsmVL7C18qvciEomgvr4emzZtwoYNGwCIu13c7MWXv/zlKb0wartIWvjff//9GBgYwODgICKRCDo7\nO1FcXJyscu64iYkJjI2NAfjkr/v58+eRl5eH4uJitLW1AQDa29tRUlKSxCrvrFu99uLiYnR0dCAS\niWBwcBADAwPIz89PYqWJFwxq1/T/29/+hry8PACp3QtVVXHs2DGsXLkSjz32WGy5iNvFrXph5HaR\n1G/4dnd3TznVc9u2bckq5Y4bHBzE4cOHAXzy85cPP/wwtm3bJsypnk1NTejp6cHIyAhsNhsURUFJ\nScktX/uZM2fQ2tqKtLQ07Ny5E1/4wheS/AqM8+lePPnkk7h48SIuXboESZKwfPly7N69O3bcO1V7\n8f777+OnP/0p8vLyYocsKisrkZ+fL9x2MVMvnn76aXR0dBi2XfDyDkREAuI3fImIBMTwJyISEMOf\niEhADH8iIgEx/ImIBMTwJyISEMOfiEhADH8iIgH9fwyjBzXf4BCfAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x114159b50>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.interpolate import interp1d\n",
      "\n",
      "def extrap1d(interpolator):\n",
      "    \"\"\"From StackOverflow http://bit.ly/1BjyRfk\"\"\"\n",
      "    xs = interpolator.x\n",
      "    ys = interpolator.y\n",
      "\n",
      "    def pointwise(x):\n",
      "        if x < xs[0]:\n",
      "            return ys[0]+(x-xs[0])*(ys[1]-ys[0])/(xs[1]-xs[0])\n",
      "        elif x > xs[-1]:\n",
      "            return ys[-1]+(x-xs[-1])*(ys[-1]-ys[-2])/(xs[-1]-xs[-2])\n",
      "        else:\n",
      "            return interpolator(x)\n",
      "\n",
      "    def ufunclike(xs):\n",
      "        return np.array(map(pointwise, np.array(xs)))\n",
      "\n",
      "    return ufunclike"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.interpolate import UnivariateSpline \n",
      "from statsmodels.distributions.empirical_distribution import ECDF\n",
      "\n",
      "ecdf = ECDF(x)\n",
      "inv_cdf = extrap1d(interp1d(ecdf.y, ecdf.x, \n",
      "                            bounds_error=False, assume_sorted=True))\n",
      "r = np.random.uniform(0, 1, 1000)\n",
      "ys = inv_cdf(r)\n",
      "\n",
      "def make_gen(xs, positive=False):\n",
      "    \"\"\"Make a random number generator with same distribution as xs.\"\"\"\n",
      "    ecdf = ECDF(xs)\n",
      "    inv_cdf = UnivariateSpline (ecdf.y[1:], ecdf.x[1:], ext='extrapolate')\n",
      "    def gen(n=1):\n",
      "        r = np.random.uniform(0, 1, n)\n",
      "        if positive:\n",
      "            return np.clip(inv_cdf(r), 0, None)\n",
      "        else:\n",
      "            return inv_cdf(r)\n",
      "    \n",
      "    return gen\n",
      "\n",
      "arrival_gen = make_gen(arrivals, positive=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nreps = 1000\n",
      "counts = np.zeros(nreps)\n",
      "for i, rep in enumerate(range(nreps)):\n",
      "    arrivals = []\n",
      "    while np.sum(arrivals) < 8*60:\n",
      "        arrivals.append(arrival_gen())\n",
      "    counts[i] = len(arrivals)\n",
      "    \n",
      "print counts.mean(), counts.var()          "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "35.218 8.92594398369\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercsie 4 (10 points)**\n",
      "\n",
      "The Cauchy distribution is given by \n",
      "$$\n",
      "f(x) = \\frac{1}{\\pi (1 + x^2)}, \\ \\ -\\infty \\lt x \\lt \\infty \n",
      "$$\n",
      "\n",
      "1. Integrate the tail probability $P(X > 2)$ using Monte Carlo\n",
      "    1. Sampling from the Cauchy distribution directly\n",
      "    2. Sampling from the uniform distribution using an appropriate change of variables\n",
      "2. Plot the 95% CI for the Monte Carlo estimates for n = 1 to 1000\n",
      "    1. For sampling from the Cauchy distribution using mulitple Monte Carlo sequences\n",
      "    2. For sampling from the uniform distribution using bootstrap samples of a single Monte Carlo sequence\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Substituting $y= 2/x$, we get\n",
      "\n",
      "$$\n",
      "\\int_0^1 \\frac{2}{\\pi(4 + y^2)} dy\n",
      "$$\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "1 - stats.cauchy().cdf(2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 44,
       "text": [
        "0.1476"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f(x):\n",
      "    return x > 2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 70
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 10000\n",
      "x = stats.cauchy().rvs(n) \n",
      "np.sum(f(x))/n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 71,
       "text": [
        "0.1427"
       ]
      }
     ],
     "prompt_number": 71
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def g(x):\n",
      "    return 2.0/(np.pi*(4 + x**2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 73
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = stats.uniform(0, 1).rvs(n)\n",
      "np.sum(g(x))/n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 75,
       "text": [
        "0.1476"
       ]
      }
     ],
     "prompt_number": 75
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 1000\n",
      "reps = 1000\n",
      "xp = np.arange(n)\n",
      "x = f(stats.cauchy().rvs((n, reps)))\n",
      "y = 1/np.arange(1, n+1)[:, None] * np.cumsum(x, axis=0)\n",
      "upper, lower = np.percentile(y, [2.5, 97.5], axis=1)\n",
      "\n",
      "plt.plot(xp, y, c='grey', alpha=0.02)\n",
      "plt.plot(xp ,y[:, 0], c='red', linewidth=1);\n",
      "plt.plot(xp, upper, 'b', xp, lower, 'b');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEECAYAAADAoTRlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXmQXVd95z/33fv2td/rVa1dsmS7ZWTZwjbIBmLHMCEB\n7CRoGMJUTVwUFEzAJAOZAcLiCa5UJpBimAkBBgNTWZiYFEsICcRJIIANGGHLttqSZVmLrd777fu7\n2/zx+ty+/dRabEvqhv59qlR6y7nnnnf6ve/vnN/5nd/RXNd1EQRBENYUgZVugCAIgnD5EfEXBEFY\ng4j4C4IgrEFE/AVBENYgIv6CIAhrEBF/QRCENYhxvgKf/vSnefTRR0mlUnziE59YtswXvvAFDh48\nSDgc5p3vfCdbtmy56A0VBEEQLh7nHfn/0i/9Eh/4wAfO+v4jjzzCzMwMn/rUp3jb297G5z//+Qu6\n8fj4+IW38hcc6YtFpC8Wkb5YRPpikYvVF+cV/6uuuop4PH7W9w8cOMArX/lKAK644grq9TqlUum8\nN5Y/5iLSF4tIXywifbGI9MUil038z0ehUCCXy3nPc7kchULhxVYrCIIgXEIuyoKvZIgQBEH4+eK8\nC77nI5vNks/nvef5fJ5sNntGufHx8SXTlf3797/YW//CIH2xiPTFItIXi0hfLLJ//37uv/9+7/nY\n2BhjY2PPu54XLf579+7lO9/5Dvv27ePo0aPE43EymcwZ5ZZr4OTk5Iu9/S8EyWSSarW60s1YFUhf\nLCJ9sYj0xSLr1q27KMbwvOL/yU9+ksOHD1OpVHjHO97BG9/4RmzbBuD222/nuuuu49FHH+Vd73oX\nkUiEd7zjHS+6UYIgCMKlRVvJlM4y8u8io5pFpC8Wkb5YRPpikXXr1l2UemSHryAIwhpExF8QBGEN\nIuIvCIKwBhHxFwRBWIOI+AuCIKxBRPwFQRDWICL+giAIaxARf0EQhDWIiL8gCMIaZEXF37Kslby9\nIAjCmkVG/oIgCGsQEX9BEIQ1iIi/IAjCGkTEXxAEYQ0i4i8IgrAGEfEXBEFYg4j4C4IgrEFE/AVB\nENYgIv6CIAhrEBF/QRCENYiIvyAIwhpExF8QBGENIuIvCIKwBhHxFwRBWIOI+AuCIKxBRPwFQRDW\nICL+giAIaxARf0EQhDWIiL8gCMIaRMRfEARhDSLiLwiCsAYR8RcEQViDiPgLgiCsQUT8BUEQ1iAi\n/oIgCGsQEX9BEIQ1iHG+AgcPHuRLX/oSjuNw6623cscddyx5v1Kp8L/+1/+iVCrhOA6ve93reNWr\nXnWp2isIgiBcBM4p/o7jcN999/GhD32IbDbL+9//fvbu3cv69eu9Mt/+9rfZsmULb37zm6lUKrzn\nPe/hlltuQdf1S954QRAE4YVxTrfPsWPHGB4eZnBwEMMw2LdvHwcOHFhSpq+vj0ajAUCz2SSZTIrw\nC4IgrHLOKf6FQoFcLuc9z2azFAqFJWVuu+02Tp8+zdvf/nbe97738Z/+03+6JA0VBEEQLh7n9fmf\nj6997Wts3ryZj370o0xPT/Oxj32MP/mTPyEajS4pNz4+zvj4uPd8//79hMNhYrHYi23Czz2hUIhk\nMrnSzVgVSF8sIn2xiPTFUu6//37v8djYGGNjY8+7jnOKfzabJZ/Pe8/z+TzZbHZJmaNHj3LnnXcC\neC6iyclJtm3btqTccg2s1+vYtv28G/2LRjKZpFqtrnQzVgXSF4tIXywifbFIMplk//79L7qec7p9\ntm3bxvT0NLOzs1iWxUMPPcTevXuXlFm3bh1PPPEEAKVSicnJSYaGhl50wwRBEIRLxzlH/rquc9dd\nd3Hvvfd6oZ7r16/ngQceAOD222/nzjvv5NOf/jTve9/7cByHt7zlLSQSicvSeEEQBOGFobmu667U\nzY8fP04kElmp268aZEq7iPTFItIXi0hfLLJu3bqLUo/s8BUEQViDiPgLgiCsQUT8BUEQ1iAi/oIg\nCGsQEX9BEIQ1iIi/IAjCGmRFxb9Wq63k7QVBENYsKyr+juOs5O0FQRDWLOL2EQRBWIOI+AuCIKxB\nRPwFQRDWICL+giAIaxARf0EQhDWIiL8gCMIaZEXF33VdOclLEARhBVhx8RcEQRAuP+L2EQRBWIOI\n+AuCIKxBRPwFQRDWICL+giAIa5AVFX9N01by9oIgCGsWifYRBEFYg6y420cMgCAIwuVnxcVfEARB\nuPysqPh3OnKYiyAIwkqwouL/7W9XJL2DIAjCCrDCbh9XfP6CIAgrgPj8BUEQ1iAi/oIgCGsQEX9B\nEIQ1iIi/IAjCGkTEXxAEYQ0i4i8IgrAGEfEXBEFYg6yo+DcaDYnzFwRBWAGM8xU4ePAgX/rSl3Ac\nh1tvvZU77rjjjDLj4+P83//7f7Ftm2QyyUc/+tFL0VZBEAThInFO8Xcch/vuu48PfehDZLNZ3v/+\n97N3717Wr1/vlanX69x333188IMfJJfLUalULnmjBUEQhBfHOd0+x44dY3h4mMHBQQzDYN++fRw4\ncGBJmR/+8IfceOON5HI5AFKp1KVrrSAIgnBROOfIv1AoeKIOkM1mOXbs2JIyU1NT2LbNPffcQ7PZ\n5LWvfS2veMUrLk1rBUEQhIvCeX3+58O2bU6cOMGHP/xh2u02f/AHf8AVV1zByMjIea+1LFnsFQRB\nWAnOKf7ZbJZ8Pu89z+fzZLPZJWVyuRzJZJJQKEQoFOKqq67i1KlTZ4j/+Pg44+Pj3vP9+/cTCBgk\nEglisdjF+Cw/t4RCIZLJ5Eo3Y1UgfbGI9MUi0hdLuf/++73HY2NjjI2NPe86zin+27ZtY3p6mtnZ\nWbLZLA899BB33333kjIvfelL+cIXvoDjOJimydNPP82v/dqvnVHXcg3sdCxqtdqaz+mfTCapVqsr\n3YxVgfTFItIXi0hfLJJMJtm/f/+Lruec4q/rOnfddRf33nuvF+q5fv16HnjgAQBuv/12RkdH2b17\nN+9973vRNI3bbrttSTTQuXBd2WMmCIKwEpzX579nzx727Nmz5LXbb799yfPXv/71vP71r3/eN3dd\n7XlfIwiCILx4VnTo7cgRvoIgCCvCioq/jPwFQRBWhhUVf9vu5vcRBEEQLi8rvOIawBHfjyAIwmVH\nfP6CIAhrkBX2+cvIXxAEYSVYYfEH13Ulp78gCMJlZoXdPhqO48joXxAE4TKz4gu+rutimubKNkMQ\nBGGNseI+f0EQBOHys+I+f0EQBOHys+I+f0EQBOHyI+kdBEEQ1iAy8hcEQViDyMhfEARhDbLi4q9p\nYgAEQRAuN+L2EQRBWIOsqPhPTqbodGSDlyAIwuVmRcX/oYeu4cEHTdnhKwiCcJlZ8S22tVp3p5fk\n9xEEQbh8rLj4i99fEATh8rPi4i/hnoIgCJefVSH+4vMXBEG4vKy4+IurXxAE4fKz4uLvd/vYtr2C\nLREEQVg7rALxB03TaDabcpyjIAjCZWLFxd+2A945vhLuKQiCcHlYcfG3rACO48hB7oIgCJeRFRd/\n09TI5/PecznQXRAE4dKz4uJv27rn8hHXjyAIwuVhxcX/s5+9npmZyEo3QxAEYU2x4uIPcPy44T1W\nI3/x/wuCIFw6VoX4t1oWtm0v2ekr4i8IgnDpWBXib9uabPASBEG4jKwK8bcsaLVaZyz6yuKvIAjC\npWFViL/jBAgEuk0Rd48gCMKl57zif/DgQd7znvfw7ne/m69//etnLXfs2DHe9KY38ZOf/OR5N+IT\nn7iRclkHRPwFQRAuB+cUf8dxuO+++/jABz7An/7pn/Lggw9y+vTpZcv91V/9Fddee+0LFu9Dh6Le\ntWoWIAiCIFwazqmyx44dY3h4mMHBQQzDYN++fRw4cOCMcv/4j//ITTfdRCqVesENqVbbaJoc7CII\ngnA5OKf4FwoFcrmc9zybzVIoFM4oc+DAAV796lcDvGABt6yAd63rumIIBEEQLiEv2r/ypS99iTe/\n+c1omvaikrO1291IH8uyzhB/ifoRBEG4uBjnejObzS5JupbP58lms0vKHD9+nE9+8pMAVKtVDh48\niGEY7N27d0m58fFxxsfHvef79+/n3/27f+Tb3/4VAHQ9TiDgkEgkiMViGIaxRPR/kdcBQqEQyWRy\npZuxKpC+WET6YhHpi6Xcf//93uOxsTHGxsaedx3nFP9t27YxPT3N7Ows2WyWhx56iLvvvntJmf/9\nv/+39/jTn/40119//RnCfyENLJdtKpUquq5jmqYn/kr01cziF9EIJJNJqtXqSjdjVSB9sYj0xSLS\nF4skk0n279//ous5p/jrus5dd93Fvffei+M43Hrrraxfv54HHngAgNtvv/1F3dzv1m80dGzbptVq\nkUgkvB2/ygBICKggCMLF45ziD7Bnzx727Nmz5LWzif473/nO53n7RUFvtXQ0TcOyrO47C35/teO3\nd8Qvi8KCIAgvnBX1oSwd+RsLr2m0Wq2zXqMWlWUmIAiC8MJZUfEfGztEJlMEoF5ftAS2bWPbtjfy\n70WEXxAE4cWxouIfjzd52ct+BHRH/kr0oSvw6rFy7/QaAzECgiAIL4wVD53JZEoAtFqGJ/j+EE/l\n2w8EAkt8/GebFQiCIAjnZ8XFf+fOp/nd3/1HGg0dx3EwTdMTdiX6y53tK4u9giAIL5wVF3+AwUGb\nYjHEd78b9nb5Qlfg/SkfJPGbIAjCxWFVqGg83mFmJsFHPnIDH//4Bu91x3HOGPH3Cr+kfhAEQXj+\nrArxj0YXQzt/8INBms0mlmUtOdlrufBO/6xAEARBuHBWhfiHQovn9+bzEVw3QDAYBPAWes81wlfi\nv9zagCAIgnAmq0L8e9duSyXD997ZF3b9awIgMwBBEIQLZVWIf++O3tlZnWazieu6tNttbNvGsiwv\n3bNyA8Hyrh//+4IgCMKZrArxn56eXvL89GmHZrNJu93GcZwlSd6WE3V/zL+sAwiCIJyfVSH+vX76\n+fkolrWY1M2/w1dFAJ1N3CX+XxAE4fysCvHv5c/+bCdf/GL/GUKvIoBgqcHwGwdYDAeV0b8gCMLy\nrDrxHxqqYlkBTp6MemLfbDaBxZG/EvVeA3Ah6R8kGkgQBOEC8vlfLn7nd/6Kvr4I5fIAn/rUK5ic\nDNNut3FdF13Xz1jQVWKv0kD4H8Oi+KvXlssXJAiCsFZZNeLf33+MkZERAoEEAFNTUe/YRv9mLxX3\n7z/oRdd1gCVC73f9+A2H3yhImghBENYqq079UqkKAMVihPvuG6ZQ6Gb5bLVaXpI327bPCPlUYaD+\nDV+9o3vZFyAIgtBlVYl/IBCgv9/0nn/hC1s5eDDsvQeLi7p+N44//UOv+Kvr/Cmh1WMRf0EQ1iqr\nSvwnJibQ9a4gJxIdAE6diqNpGp1OZ4kP3z/y7zUEinNlAfW7jgRBENYaq0r8NU2jXC7zmc98lb17\nTwFw4kSUYDDoiXyn0/Fy/itR9xsB/6aw8wm7iL8gCGuVVSX+amTfaJwileqmfDh5MrZE1NUCb28E\njxrJL5cGWrHce+dLGicIgvCLyKoSfwDDMNA0jVSqG9t/8mSMz31ugFJJ90I/VdoH/ylflmUtG77Z\nuzAMS9NESDoIQRDWIqtO/Kenp7Esi1e+ssgrX/ksuu7ymc+s48CBhFfGPxPoFfdAIODtC/CfCaz+\n+cVeXSfuH0EQ1hqrTvyhO/qPROb43d99lG3b6riuxoEDCb761awn1qZp4roupml6hsC2bQKBgLc3\nQL2mjACcudDrDwl9MQZAjIcgCD9PrErxDwS6h7nYts3QUB6Af/iHIT7+8U00Gouj/Ha7jWmaXhqI\n3hkBdMVerREAS2YAftHvjSB6PvSGmQqCIKx2VqX4l8tlWq0WmqYxPNz1/bdaOp1OgD/8w1H+8i/7\n6XQ63sYu5e9Xxz+q2YBKDaHEXs0I/Au/ahahjIaqU+GPHuplOdeRGABBEH4eWJXiHwwGCQaDzM/P\nc8stM7z85c967/3DPwzwzW+mPZFuNpueXx/wHnc6Hc/37xdm5QpSi8S2bXsGwDRN733Lsuh0Op64\nq8Nk/Nf3po2QtQNBEH5eWJXi3263vRj+YHCad7/7RwSDi6PvI0fi/OQnMX72M53vfCeFaS7uCrYs\nC2DJngA1M1BCr1AGAJb67NU1ymCo3EGqnN+lpNYY1HMJHRUE4eeBVZPYrZennnqKTZs20el0iEaj\nfOpT/0AgkObtb38FiYTFH/zBdnTdoVgMk0g8TSIRYN8+09sA1ul0sG2bUCgELE0P0ZveQQm8el2N\n/IPBoHedYRhLDMS5ooSWO4RGEARhNbEqxV+J6vz8PKlUCtd1GRqycJw53va2wziOzec/v8sr/5GP\nbGFoqM0NNzTIZGw0zea3f3uWQCDguXJ0XScUCmFZFqFQaMmoH/AWhV3XJRgMekZCuYr82UQBzxW0\nXPZQ/zqAqseflvr5GANJPy0IwqVgVYq/P1NnOBxmZmaGTCZDIBDgjjtOUiyG+PzndxGNWmSzbSYm\n4lQqQU6ciGHbGuGwzW/8RoFstivEKixU4Y8IUkLvjwJSG8iUfx9YYhD8awjKACiXkj+PkG3bXtSS\nqsdvRPwRRv4F6N7Mo35DIgiCcDFYleKvaLfblEolTNMkkUhgGAb1ep1QqM6v//oRdu8uEIs5/NEf\nXc/OnXUefjiLrruk0yb/5/8MMDUV47/9tylmZ10efzzLG984QzBoeG6hcDjsRQNB1ygYhrHkDAHo\nCrZaQDYMw3vPnylUhZS6rothLHar4zgEg0Ha7TbQnWH0GgI/yhCodQb/TGOlzyBYDW0QBOHisKrF\nPxwOUyqV0HUdy7IoFAoMDAwQDod5y1sOe26ZT3zi+wSDBsnkVq6/vsqGDTbvfOceAI4ciVAqGTSb\nOuGwSbkc5Td+Y5ZyOcBjj8W48846mrZ4OphyCfkXeVWkjxI/NTPwbyjTdX1JWKlf1G3b9gyCMg5K\n5HsNhf+a3pnEcnmJ1P+qvt7U1+qey80olluz6DVGyx2Gs9z5yctdKwjC6mVVi78aLUM3AicajdJo\nNAiFQl6YZzQaZXjYpVot8Du/0yAWi9Fut0kkdlGrBZmaCpNMWmzf3uR//I9NdDo63/xmhvn5ILWa\nQaPxHLfcUqHT0ZmbC3PTTXWi0egSN4wS9E6nc0Y2UbUQrEb0yhj0RgKp1/yLweoa5YZSMwd1ndq8\n5hd5dQ//+oHf4PTmLep1U/k522J1b9TS2c5H7q3zfIZEjIMgrB40dwUD0++5557nVd4wDPr7+7Ft\nm2g06rlowuGwJ3RqUdc0o7RaLj/60U5uvrlEMuny2c+uY3IyxuOPxxkaMtmwweTw4TCxmEu1qtPp\naHzxi8fYtauN6y4KuGEYnjDX63XvPkr0DcMgGAyeMSJW1/iNAeAtIqv3Y7EY9Xp9iaunN+2EX4DV\nLKTXsPhzGvmvP9+JZn6WG+k/X9H2R1Etx7nqTCaTVKvV53W/X1SkLxaRvlhk3bp1F6WeCxr5Hzx4\nkC996Us4jsOtt97KHXfcseT9H/zgB/zd3/0drusSjUZ561vfyqZNmy5KA/2oEMxGo0EkEvE2YkWj\nUVqtlic23YNfmiQSEe6442kAUqkU733v0zSbDocPZxgbqxOPB/na19J88pPb2L69zd69dd773s2s\nW9dh584q4+MJ/vt/n2TTpjamqZPP62zeHKPZbHqj/larRTAY9GYmanFX0zRarRahUGiJ2KkoIzWy\n13XdCyn17yT2j7yVkDuO49XnF3oVpup3FalZgrqn/3mvO8lvJC7GwvK5QlzPlgpjuU1yEukkCJeO\n84q/4zjcd999fOhDHyKbzfL+97+fvXv3sn79eq/M0NAQ99xzD7FYjIMHD/K5z32Oe++995I0uNFo\nLBExld+n3W5jGAb5fJ5IJOKFdubzefr6+igWi2iahmHo7NjxHPF4hmazxK23zrF3b4MNG+pEImFC\noWG++MVRpqaCNBoB3v72zYyOmuTzOnNzId761hlOnBjiwx+e5J//OcJf/MUG/uRPnqZSCXPVVRWg\nOyJXxslxHBKJxJK0Ekr0lVirqCH/iNm/Ic3vzmm324TDYYLB4Bnv9S4O+5+rjXDKaKj9EH73lH8h\nuzeE9WxG4fkKdG+4rHrNH/W03GlsEu0kCBeX84r/sWPHGB4eZnBwEIB9+/Zx4MCBJeK/Y8cO7/H2\n7dvJ5/OXoKldIVA+/6mpKSKRCPF4nHK57AlItVqlr68Px3GYmZkhlUoxPz9PMBj03DeBQIB6vU6r\n1SKRSBCNFqjVOtTrNd785jJ79hQYGKjSbCZ49NEkn/70ZnbtqvGqV1X48z8fJps1eetbN3HwYJxr\nr63z2799FTMzIV7zmgLvfvcppqbCXHddjWAwga47lEolT0xDodASd029XkfTNBqNhjcDUP96E84p\nF1S73faMnQpBVSjh7N2VrBaWlVHxj7JVKKzaj6BmWGqGok5SU+1RdZ8roZ3/HmcT7d51BFW3/7Or\nz+Rf+zjb9auVs+349hv35fpJUoUIl5Lzin+hUCCXy3nPs9ksx44dO2v5f/3Xf2XPnj0XdPOBmRnm\nhoYuqKyfTqd7vq9KA6FpGolEgna7jW3bPP3002zYsMGbGbTbbVqtFvF4nHQ6ja7rFAoFDMOg0Wh4\nQtdutwmFQmzZ0sAwDNLpMnfcEebmm2dYt667Oey226YYGenwmc9sIJXq8MlPPsvf/V2Kyck4p04Z\n7N//EhqNANdfX+HYsRgf+9hzbNpkEY9DKmVRq5XQ9QiRSABNi9BqWQQCbUqlDtFoB9MMEo06xGIx\ngCWuGyWKavSudjGraCI1elduKP+CsnIXqRmDEhv/kZfLZTTVdX3JsZn+WUR3JmWcIVy9Ytc7el8u\nIkldp9rhv265sv5y/vf9i9W9i929MyJ/O/yL27339t/jfO+fTeh7jbn/c5xtjaT3M/Z+zt6y5+J8\nxvGFbEAUfr65qNE+hw4d4rvf/S5/+Id/eMZ74+PjjI+Pe8/379/PO//8z7nnox99wfcLh8NUKpUl\nxikej1OpVACIxWJUKhUsyyKZTJJKpTyBSSQSdDodb+bgdznE43Fs2/Yyi46O2kQi3Qig3bt1Wi2N\n//pfJwiFQoTDEX7zNxtYVgXLsrn//hybNlW4//5RfvM3Z/nYxzbQbGrousbwcIdg0OX48TBXXdXi\n8OEI69Z10LQAExMGW7e2efrpCK9/fYHHH4/ze79XZMcOm1LJ4Ec/CrJpU4dm0+CaayxM02B4uEWt\nZtHXt5iKotVqeQvSSqTViN6/MK02uWmaRiQS8Ub//uMy1WxBXafq8kcYqdmBclWpk9iWizBSswuV\nIM+fMkPVads24XDYu0a1w1+3P+TWf6Sn/36991efpTeSye9687u7ev9fznio+tT7/ud+ke9971wL\n4X5UahK/0ViuXn+Zs9XlL7ecEfJHkPnr6C3b2y/+envbcTENSSgUIplMXrT6ft65//77vcdjY2OM\njY097zrOG+1z9OhRvvKVr/DBD34QgK997WtomnbGou+pU6f4+Mc/zgc/+EGGh4cv8O7aixJ/RTab\npVKpeKkgWq3u+b/JZJJ4PE6xWCSVShGLxbzsneFwmEQiQaPRWBKyWSwWCYVCpFIpLzeQWlRWO3aj\n0eiSSB4lWr0x+d3MpBrf+MZ6brihxpe/vJk9ewrMz0fodBw2bmxw5EiSG2+0eO65AIbhsGkT/NM/\nJRgYqPPTn/ZTKoVot3VGRjpMTITYvLnJs89GiUZtgkGXVivAq19d5Ikn4uzfP89LXtLmHe/YzMte\nVmLXrhqlUohGI8DAQJPDh1PcfPMMR4+m+eVfLvPkk3Fe85oG1WoHywoxOxvi5S9vMDenMTtrMDTk\n0miY5HIakcjimoFyXanRtUqcp2YlfqFebu1A1aPWOrpfhW7ZVCpFPp/3DJa/P1U5/w5s/0Y59bh3\n/cK/c7u3Lv//vWsnfvecqkcZEf+sYjnBV/jrUu3pNYpnmxX0Rricq24/ZzNUyxnGsxmMXkPqN1y9\nMy7/mlBvX/Ze38vZZoG9JBIJarXaOcss9zlfLGeTx+XcdMv9/S/U0J+N5cpdrGif84q/bdu85z3v\nWbLge/fddy/x+c/Pz3PPPffwrne9a4n///x3vzjirwiFQkSjUTqdDq1Wi4GBATqdDu12m1gshmVZ\nXqTOtm3bsG2bvsceY+OTT3LgN36DeqNBMpmkUqkQCoU8f3o4HKZer5NIJHBdl0gk4oWUtlotDMPw\nZhJK9ILBoOdeUiLmz++jRtexWIxwOLzkR6V87I2GzQ9/OMC115bIZHRqNYtUyuLIkSQbN1Y5cGA9\nV1wxzze+sZls1mZ8PMFDD2XZt6/Arl1Vvva1IXbvrjI1FSKVskkkbI4fjxKPO5w+HSEetygWg9h2\ngGYzwMhIB9PUmJvrjjhjMZtGQyeZtHjFK4rUajobNzZ53esqPP10iBtuaHLsWJjRUZNnnw1x000N\nXBeCQbCsxUyr3UVujWQyskSMVZ/4fzTxeJxqtbrEsChD6t9I56+ndx1DzRSUsKjZiV/M1cxI/b38\nBwApo+QXUf/fTrVHlfX/bf3rIWcbHcOZI+jljEMikaBSqZwhrKouv0D33ksZYf9n6jU8vQZM/e/f\nN9I7M+r9PH4RXM6Y+P+2vdFrqsxyYt17r1QqRbVaPedn7qXX2Pnb1Xu/5cr11uN/3juD7GU5g3eh\nrrfz4dfeF8MFxfk/+uijS0I977zzTh544AEAbr/9dj7zmc/w8MMP09/fD3S/eH/0R390AXe/eOKv\n/Nzr169nfn7eG/13b7P0j5pOd88D6O/v5/q/+Auu/v73Gf/VX+Wxl7+c0OQkrZ07sefnafT1EQ6H\nl4x0HcehVquRTqcpl8uk02lvPSGVStFoNIjH40D3y9loNDwRa7fbNJtNwuGwN0tRC9GAt/MXFgVL\n/fjq9brnjlKuGuXv13WdVquFrutMT4cYHjaJxaLejAS6LrJOp0MoFFoQxjCdjolhJDl+PMTwcIO+\nPpvx8RjxeIBwuEk2q3HyZJSREYtvfjPNwYM5Nmyo861vDWOaAWxbw3W7X+hk0kLXXUqlILlcB9AW\nXnOYnw/3Nf4FAAAgAElEQVRTLhvccEOZ06dD7NxZZ3o6zK5dTQIBk2g0wNVXN/n61/vodIJMTga5\n8so60ahDIAC33FJlbi7Mxo0N5ucN5uY0NM0gmbSYnQ2Sz4fYvbvBTTdVqdVckkmXZ58NccUVLSYm\nQlx5pUWrBYmEhuPo6LoJuAtRVwawKH69i9vqNVjeeCvXl8Ivkv6Fe/V3NQzDyxirZin+w4eU8dI0\njb6+Pm+Hu3rPP+PxG7heF1zvyNtvBNT3rDfhoHKj9T7vnUkp/IbOj7qHv57etZezzZrU99U/wwoE\nAiQSCW9Q0Ov6U4EQ5zJEve1Un2M543Y2V59f9M+Fv6/VZzuXG05xNmPRa8Avq/hfMjSNP/zQh3B8\n0/sXi/ryRqNR6vW697r60QUsi30PPsjBV78aKxDgP3z2szz1yley95vfJGiaRGo1TlxxBetPneLL\nv/VbmC99KYZhMDc3RzKZ9BaI1Wyi0WiQy+WoVqveTCMajdJut70ZguN0zxYolUqk02mq1Sq5XA7L\nsshkMt5+BbWw2mq1yOVy3gK0Em5d12k2m6RSKTqdjmdUms2mZxCUYVLvt9ttksmkJ/z+sE7lXzdN\nk3g87u2YrtVqxGIxqtUq8XicRqNBIpHAtm3a7Ta6HqLdtjCMKJVKg1QqiWE0KRZTHDoUZds2i/Hx\nOJs3d3jyySjbt7dZv77G4cM5Wi1ot6HTCdFu2xw9mmZwsMljj/WzbVuNm2+u0W5buK7GoUMxRkct\nDhxIU68HCQRsyuUQIyMthobanDwZZWioxeBgh6mpII8+mqNSCREIuPT1tcnnI8TjJrat0WoZxOMW\n9bpBJtMhnbaoVg0qFYPh4Ta5nEW9HmBgwKFW0+jvt6lWA/T3W1QqAXbubDI7q1MohEgmLXbvbtLp\naMTjDsPDFnNzOtPTQUZG2szMRHjmmTCbNrUYG6sDAUIhl127WmiaC9hMT+scPx5neLhBJKJRKBiU\nSiGuvLLNpk1NajWNQCBBPl+jVDKoVHRs22HTJpP5+RDhsIZlwZYtHaang4RCDpZlUCxqDA1pbNrU\notMxiERMLCtIJOJSKATo6zMJBDRcFxwHNG3xxDtYmr+pK3Y6tu2i62CaNuVykEzGJBA4c3G9+5PW\nvAHJ4s/8zBlZ76zN/51URsAvfslkkmazucSQ+t1evTNK/yxPPfcb196d8/46/e7D3lF+r7Hy1+t/\n7Dei6r3ezZ5+N6TfuCz3WfwzjG3btr04kVR/l5UW/z/+/d+ntRDZ8mIZHBxkdnaWSCSyZOSvWxY3\n/fjH/Pimm7jtn/+Zl/34xxzcvZtvvOEN/Nc//mP+7Hd/l2Cjwbs/9Sl+cPPNvPSnP+XJsTF2PvUU\nX3/jGzl95ZXdhVTLItHpUPG117/Bq9lsMlqpcOVPfsKjr30tBdP0Zgh9fX00KxXiC7ODaCxGtVaj\nv7+fZrPppY5wHMdLOKdpmreTORQKoes6+Xzei2xyHId2u00mk/GMhq7rhMNh7331T9d1T+BjsRhz\nc3Nks1lP4Ov1uif4sViMWq1GNptdss7RbDaJRCILyfVCNBbcZOVy2evzZDJJu932UnGk02kajcYZ\nPwpd1702qufKGKqFdhWxpEbeytWmRMOfWqNLgFbLxjBcHMdC12MEAhaNRpBq1SUW06nXA4RCOuPj\nQTZtcmg2beLxIEeOROnvbzMxEaG/v830dJyBAYupqTB9fW1On85gGCbxuE0i0eaZZwaZnQ2STnc4\ncSJJNtshk+nQbhtoms22bTXqdTh6NMPsbIxsts0zz6RoNIIEgw6BgMuWLVU6HZ2JiSgjI036+trM\nzEQX1oS6nymT6eC6kEqZ5HImx47FyWY71Go6uVyHZ55JkMt1qNUM+vvb2LaGrrvMzHTrCIdt2m2d\nQMAlEHAJhRxiMZtmM0CzaZBOm4yMNCmVQgwMdCgWgwwMdMjng5RKQTqdAK1WgHjcoloNEQw66LrD\nxo1NOp0A0ahDp6MRidhYlkYk4rJ+fR3HCdBs6jQaIdavbwA6jYZGtaqTTpuUyyFKJYNczqJY1Onv\nN9m4sUVfX4fJyTCVikFfn02lYqBp3UHd0FCVctkgl+u+nsvZDA83qNWC1Ovw5JMJ0mmT/n6LmZkg\npVKIoaHmQt8nMU2d0dEWwaCDbQcIh11ME4aGmrRaOrOzEWIxi3DYxjQDDA+3GBxsYFlBEoluWvla\nLczgoEmnoxEMOiQSATodjUTCod2GUMjFNDVaLR1dt4hEHCyrazxt22V21iCb7db/zDMhhoZgcNAk\nEtGxLBvL0gkEHIJBzXMv9u4L2r1790XRyxUX/0++5z2UM5lLepvhyUne/rnP8dSOHew8epT5XI5O\nKMTM8DBbn3mGT/7e7wFgmCaWYRCv1WjE42x47jn2/83f8I+/8iuMX3MNt3z/++z74Q/5xhvewOGe\n1fVQu83A3Bxv/fznKWYyaK7LN1//eo4vWGlN0/j1+++nv1ikEQ4TabX4zmtew7ObN3fvvTAzCYfD\nSxaR/VN2tQitDIpKfLdp0yZKpZIX6aSENRwOY9s2lUrFm+Irw9L7pVJCGwgEFnZIL47UVLhnLBbz\nIqSq1aq3IJlKpbx1knq9Tn9/P5VKxVtsTyQSWJZFLBaj1WoxODhIuVwmFAotObVNfX41u1LrNMFg\nkFarRSwW89rT6XS852rvg9qToNodCoVoNptks1kajW7eJ7UHpd1ue4bbb7iUcVQzt3q9TjAYpNPp\neBv31N8KFkdl/gR9aj1BudjU+oFtu5RKNokEOI5LMNg1YqlUimaz6d2/1QovGMIgltXCdcG2LS/4\nQKUz6X6ODtFoBE3TvH4Jh8PUajYQxHUd2m0bcEindUoljUIhQTrdFWTXdZmYSJNM2lQqISKRNtVq\niGi0QzIJ4XCT/n6LfD7G0JBFINDBssI8+2yMViuErjdpt8MEAm1MM4qut5mZyVCtBsjlIJFokM8n\nqVZ1crkAiUSNajVJPF4nkbBoNDRCIWi1YGIizuxs17gNDDQplSIkkx0sK0hfHzz7bJhAwMY0dTTN\nwXF0JidjaJpLJGKxc2cZx9GYmYmSSnUYGmqTz4eYnY0zMFBneLjD9HSEUilMItGhWo2QTJrMzERx\nHI2hoRau63j3nZ2NMDUVwXE0wGVmJkoiYVGpBLGsALru0G4bGIaNrrs0m0FiMZNGI4hh2AQCEArZ\n1GohkskO1WqIeNyk09FxXRgdbdJsBigUwgsGMeh9h/r6OiSTFvV6t2wgwMIgwOLo0YsT9bTi4v/n\n73gHsy8g1v/5MPbEE1z3yCNkCwXyuRxfv/NOAP7LJz7BfC7Hn73rXWe9dmh6mjd9+cscvuoqXvbj\nH/PDm29m7NAhDl91Ff/yy7/MluPHObFlCzc8/DCv+ad/AuBP3vc+RiYn+bW//3uOb9vGA7/8y9i6\nzu/96Z/ywO23c/3BgzyyezeveOghJoaGeOzGG7ni0CEe/aVfYiKRIFcqse973+MnN97IzMgIAds+\nq2tMCYIiEAgQLZdpxGK4KleQ4zA6P096Zoand+xg06lTjMzNUY1GObJrF82FWYVyA6kNZH6h672P\nQgltNpulWCx6s4YLQdXpjxRSU3J1DgIsJvhT+ZTUbERNp9VsQIW2+iO0ms2m9xlSqRT1et0LGYxE\nIhSLRc9NpxZ61dRbGR7/lFvNMpQ7z78Bze83VkKsZiv+86TVrE6VU4ZEfS5VXzAY9BZ8lY/f7/bw\nR1yp4AH1t1LGR6UYMU3T6zM1E1SL534DpXak+xfRVf3+jLfKCCujq4ywZVlEIhGvjaZpeoEQnU7H\nS28eiUS85+p7pD6TP3pMuSpTqRTFYtFro3JX+tcgXNf1XlcDAvUZ/QOG3u+N6iP1t1DfQfV51fdJ\n9adhBLFth3A4RK1mE4m4lMsG8bhDqxUmHG6jaQ6GoVGpRIjHber1INGoia63MIwQpqlhGObCdz1K\nPh8mHLYwDJdYLEA+H6FQCJNMthdmAy0cJ0osZvKf//OvXdDv63ysuPjfd9ddnN648ZLe5pXf+x4B\n2+boDTfgVCpMjY4CkCkUCLfbzIyMnPP6eLXKv/+bv6EVifDXv/VbRJtN7vzqV4nX66ybmmJqIbT1\n4Rtv5NTGjRQX9h2EWi1u+5d/4eonn+TY9u0kajX+6j/+R69ewzS56cc/5uYf/IDp4WH65+d5+oor\nCDgOyWqV/vl5JtetY3RigunhYR7at4/1p08zOzjI0R07cNUClWXx7779beqJBDueeoq+YpFOKMSj\ne/bQVyySy+dZPzHB5MgIA3NzdEIh2uEw08PDbD1+nGPbt5Pv76d/bo4nr76ap3fswFyIMb/cnC3i\n4kJQP/SV/EpfapQQK1edQs1e/OdFCM+fXv/6akPTND784Q9fnLpWWvz/8i1v4Znt2y/pbX79b/+W\nY9u38/i1177gOgKWRdCyaEci3Rcch5f/6EcYlkU1kWDfgw/yube9jY5638fw5CSv+c53+OkNN/Dk\nMpsxIq0WpmFgWBY3/ehHXP+zn/HFu+6imkyy55FHSJfLzA8M8PIHHyTSalHKZEjUavz0pS/l0euu\n47Xf+hbXHDrEobExCtks8wMDzA4Ocv3PfkY2n+e5DRs4tGsX+f5+UpUKlq7TSCS69240GBsf5+on\nn+T41q1sPnmS9adPewZhcGaGozt3cnTHDhoLUUyCIKwcH/nIRy5KPSsu/l954xuXFcQXin/kMzg9\nTbpc5lXf+x7/8Ku/ysRFCpFaMRyHcLtNOxpl3enT3PDww+w4ehTNdfmfd9990RbOY/U6Vx0+zM4j\nRzixZQvrT59m6/HjzAwNceTKKzly5ZVYwSC/+ZWvMDk6Sj6XY2Rqime2bWN2cJB1k5MUMxlGpqd5\nduNGZgYHu05LQbgM6KaJbRjQEzppmCYBx6Hj20H+80TAsnB0nY9cpPD4FRf/b7z+9Ry87rpLUv2r\nv/tdbvjBD3A1jU+89720otFLcp+zofznl5J4tUpfqcTpDRsu6X1002TriRNceeQIO556ikS9ziN7\n9lBJp9n6zDM8t3EjQ9PTbD55konRUUampji+dSsDc3NEWi1ObNnCia1bOb5lC6Vs9pK29eeVZLmM\nZRj05/NkCwVq8Ti1RIK+Uol8NssVx46RKZWYz+VIVyo0YjHmczmCpsnY+DiNeJx0qUSxr498Lkc+\nlyNZq9EKh8nnchSyWcxQiKHpafpKJdLlMuV0mni9TicUwtZ1Iq0WgzMzDE9PU8hmKWSzRJtNLMPA\nNgx2PfEElVSKSKtFpNUiUat55Up9fbRDIbacOIGraYQ6HVxNo9jXR6mvj2ImQzsSIV6rUcpkKGUy\nNGMxthw/TsB1aUYiJGs1dMsi4LqU0mmvzKZTp9Ati2StRjMSIeC66LZNKZ2mGYuRzefZ8NxzbD55\nklSlgq3rlNNpyuk0lXSaaKPB6MQEkVbrjPfq8TjJapVqIkEller+S6epxeOkqlX65+YYmJ8nUa1S\nTySoJJNUUinqC+9vf/pp+ufnsXWdajLZ/ZdKUVn4v5pM0g6FWD8xgQteGTMUIlMsEms2qSaTdEIh\nYvU6/fk8O556inS5TDscJl6vkymVSNRqlNNpMqXSRfm+rbj4f/s1r+EnL3vZRa86YFl86GMf47GX\nvIShmRm+ePfdyy5YCs8fzXEYnp5mdnCwO8LyEbBtnG4QuPdaqlRi64kTbDl+nC0nTmAZBuV0mnWT\nk5zYsoVKKtVd3xgdJWiaBE2TZzdsoJjNMjQzQ6TZJNzpUE6nGZidZXRigkoqxeS6dUyMjjI5Okrz\nAmc9sXodMxh8QWsaufl5BmdmmB0aohGLkS0UmBsYWHYkqZsm1z72GNVkkrn+fvrzeerxOLl83nut\nnkig2zYveewxrv/ZzxiZmiLguuSzWWqJRNcQzM9TymSI1+tUk0me2rmT/vl53EAAc+H9eL3Ogb17\niTWbVJJJoq0Wufl5cguRTZ1QiNyCQbF1Hc11mR4eppxOE263CTgOrUiEVKVCKZPBDIV48uqr6SsW\nyRYK6LZNKxIhWanwxDXXYFgWhmVRSyQo9vWRKZW6ZYtFUpUKE6OjVJNJDNOkGY3SVyqRKRbJlEqk\nKhUqqRSJWo1MqUSmVGJy3TrK6bT3XsBxqMfjpMtlMgtGamZoiFImg6XCftttGvE4qXLZu25qZIRn\ntm8nn8vhaBrpSoV0uUy6XMbRNKZGRpgZHibabJJaeD1dqZCqVKjH44TbbVKVildnulymHo9TTSZ5\ndtMmGrEYoU6HVKVCcuG6RizGqc2bmR4exgwGu69XqyQrFZLVKslqlVSlQrxeZ3p4uNuPC69HWy2q\niYRngOP1OrVEgno8zpNjYxT7+gjYNp1QiHImQymdRgM+cJHS5a+4+H//llv47m23XbQq1Wh7eGqK\nt3/2s/zPd7+bUl8fLERPqIgQYYVwXfrn5tj07LNMDQ+TKZXY8NxzFHI5YvU6uXye6ZERbxYxNzCA\nq2lMDw8zMDdHKZOhkMtRSSYZnZhgdGKCdZOT1ONxJkZHvX/TIyOeUEB3cf/Xv/pVNpw+jWkYlDIZ\npoeHu+V0nWStxszQEKZhsPfAASzDoJjN0j83R7Gvj2yhwPZnnuHkpk2kKhWyxSKFvj4StRr1eJzZ\nwUHmBgdphcPse/BBbF1nvr8fMxhkeHqaViSCq2m0w2GcQICBuTk01yXU6TA3MMAj11/Pwd27MYPB\nS+cic5yu0CUSZxht4RLjume4oc4VxXcufmF8/o/s2cM33/CGF1xFvFrlTf/v//FXv/Vb9M/PkymX\nOXTNNVz55JPsefxxvvymNy0p/4t4HNxaj/DQHIfc/LxnDEYnJhiYm6OQzTI1MkK+v59Xffe7PLhv\nHz996UtpxOMMzM0xPDXFyPQ02XyeuYEBMqUSuXyep3fsYGp4mOGZGWrxOJrrEmm1utcuLJQHLAvH\nMNAch0yxyODsLIOzs+TyeX62dy/NaJR8NuuF256B6xKr14m2WuQX0qIIwoVwscR/Rc2/peskLjAm\n/GxkCwXWT0zwa3//9zSHhrjue9+jHo/TVyxSyGS8U7wU1WqVDRs28Nxzz73Y5q8aVpPwq5ho/yzr\nUs+43ECA+cFB5gcHeWzhLAndshiYnWVkaoqRqSn+5k1v4mlf0sGZ4WFmhod57Bz1Prlr11nfcxZG\nzm4gQDGXo5jL8dRVV114ozWNRiJBI5FYslHMnyfoF53l9o4s1weRSMTb36HyF0UikSUb8vz959/1\nrTYt+utX+0XUOR4q7l+1pzfU81xhtL1JCdVu9rOh8hCpfQOqbf4QU7UZE/A2X6r9B9GLuG65oiP/\nYl8fzWiU//P2t7Ph1ClsXWfyeUbkXPP441wzPk66UGBwbo5HrruOHU89RX7DBp7cvJmHb7qJkZER\npqamllyndqWejcs9mj7bJqrLyfli7F9IG89X56Xq51QqhWma2LZNf38/c3Nz3g/HcRwGBgaYn59H\n0zQvOV80GvXyNqkNTGo3dKPR8DaSBYNBUqmUJzrNZtMr1263vU1doVAIwzBoNptLkrCpMyRisZi3\noUidH6GSunU6HdavX+/t2A6FQl4qDtd1qVQqnlENBoPeLmSVBwrwBFLTNJrN5pKcNYFAwNtN7v+7\nGoaBaQYIhwM4joVp6sRiAWq17mvQpNlMYNsOp08H2LBBI5lsEgiArhtYlonjRDEMqNdruG4WXc9T\nqeRIpWyCwSqNRh/hcIupqW5q8qmpCMmkzcBAh1rNwbajTE2FiMUs+vstJiZ0otEo09MG6XSTZFJj\nejpEuZwmGOyQTjcpl5M4jk2plCSVahCP16hU4oTDLuVyjGjUJR6vU60mqFbDtFohstkWyWSdTKZF\nqRQiHg8CGrVaAAhQqxnUahF03SWTaaPrGo1GmGSygaa5VCpBarUQuq6RTDZoNIIEAlCvR9A0l0DA\nWdjd29213f1emuTz3Z3CyaRJPh8mFjNJp20Khe6GsFSqw9xcN6VEpxPAsnQKhQiVik4k4vD4430X\n5TeyoiP/ajLprVxfeeQI0WaTv3se4h+JRPj1r36V6aEhDrzvfbz293+fg7t3MzUywq9+61s8+apX\nATA1NcWVV17JkSNHvGvPJfzw/EfTajegf9TSmy73bKTTae8Hn81mKRQK3uvlchmA/v5+5ufnz1qH\nuk7tao1Go2ha93hIWEx0pXYpxmIxgsEg9XqdSCSyZEetSs0QDocpFAreEZ6zs7Oe20zXdS+NtTrf\nIJ/Pk0qlKJVKRCIR+vr6mJqa8s5HmJ+f94yuruteErlwOEw2m8WyLC8lhEqprbKgql2kKk12OBwl\nGDSYnY0zNNRB0ywsq5s3p9GoeTtYm02bej1ELpfDtm2vf0zTJJPJeCM9lYohHo+jad0TygYGBqjV\nahiGQTKZXNhhm8Q020QiNonEEKbZYmDA9naidneSOjiORTDYTYrW6eicPBlhaKhNraZTq2ls3mzx\n3HNpTLPC0FATTbMol1OYZoRSSWPDBovnnuvDcUpkMvDMMxGGhkxs2+TUqTgTExlSKY1crsbsbALD\naDM5GSMScZmcjFAoRKnXw2SzNYaH27iuxsxMjMnJOJlMC9fVqFa7aQemphJkMi3m56MEgw7lchhw\nSadN8vkI0Wg3JYGmuQSDLpalYVkB1q+vUy4HaTaDZLMtKpXuIrppKt92N+dPvW4Qj5uARrOpe7l1\nMpludtVUykTTHPL5CM2mzuBgi3i8u9t1fj5MPG4TjWpks3WOHjUoFIIMDbXJ5Uxc1+Gpp0KEw92d\ntlu3WlSrGebmQoTDLpoGg4MOhmEyPz9INOoyMuIwONimUIhTLGY5ftzAMFza7W4K86GhDs2mzpYt\nFsGgTSzmMDcXp1rV2b69TamUpNkMs22bSTDYTX2ez6fRNJdYzAashTQQBsPDNYpFHddtk0w6lEoh\nUqkaluVSLIbYu7eDaTqUSgY7drSxLJdSKcTgoInjuLiuTjLZIpGokE67DA7awC+A+NcSCUYnJsBx\n0G2bjc8+e0HXqR+wSt72jTe8gYZhcPi//BdqiQTGq17FV2IxjvX1MTo6ysTEBEeOHGHXrl0cOnQI\n6B46PzMzc0H3O990XO22TKVSGIbhrSmoUZlKuZxYmOJXKhX6+/uJRqOcPn2acDjsjR7VfaLRqCeO\naqSpprJ9fX3k83lM06TRaJDNZnFdl6GhoYWRm+mNWtU9DcOgUCjQ19f94miatuR0JGUUNE1jeHiY\nqakpUqkUQ0NDzM/Pk0gkuHIhwd2GDRuwbZtiscjGjRup1+u4rsvAwADFYpG+vj5vO/727dtJpdJU\nqxXi8QEgwPr1Lo4TwnUD1GphKpUIzzzTplJJc9VVJo8/niCRaLNlSxXLCvPoozmOH88wPNzk6acz\ntFoahw/3kUqZNBrdpGUjIzVOnUoTiVhs2FClv9/iqadStFrdBGaZTJuNG2ts3Njkuefi6LqN4wSY\nmIiyeXMNCDA5GWbDhganT8cxjG4itImJOJs21envb3HsWJLZ2chCQi+LUimIrkNfX5vR0SazsxEc\nB4rFMLYNIyMtCoUQtVqQoaEWtZpBvW6weXODiYkIqZRFLGYzMxNe+Bs42DZs29ZkaipMPO4Qj5tM\nTUXYuLGbB2Zysvt4ZKRDItHhJz8ZIRx2icdNIhHQdZtNm0xuuaXJxo1FZmYMZmbCFAoa11xTY/Pm\nEq0WlEoGW7dazM/rDA/P0Ww6hEIBMhkHTXPIZGzm5sIMDJg4jo6mGaRSNoWCTjLZPUgoENAWZj4w\nPd2VkkymWz6VMqlUNDIZ6HS6qQ7AodXSCYc1AgFnwVWzmPaiNwW2P6tnIpHwvmfK0C4mCmwtSR/t\nP+TH/7r6f/E+Jrre9rLuqkSB/qNN/dk8/Wmq1fu9+NvsT+m8mBK76Svd7ElMuLQefxp59Xk07eJl\nQF5Zn79h0A6HiTUaGLZNrlAgVq8v2Um63HbrZrPbgVccPQp0/bdupcLgtm3UZmc5ceIE6Ze/nE65\nzMTEBJs2beLUqVMcOnSIHTt2cPToUWZmZhgbG/OOllSPdV3n6quv5oknnlhs54Kfcblt35lMhnK5\nTDKZJBQKMT8/TzabRdM08vk88Xh8yWhbJSPr7++nVqsxNDTEwMAAExMT5HI5zx+ZSCRIp9PeQTEq\nX4pKvpZMJj1xVwai1WphmibJZJJgMOjVFYlE0fXFox2VQWi32yQSCXRdp163SSZDnD7dT6HQYcuW\nEO12iGPHAmzfbnD0aDd52DXXZPj+9+OMjjpcddUg5XKN+fnNpNMu//IvSbZta/Jv/zZKLFbn4YeH\n6O+3OHCgj0ymQ7MZpFIxGBpq0WwGKJVCmKbOunUN5ufD3HBDgW99K044bLNxY4tvfCPGc8/FeNnL\nioyN1XCcMDfc0GTrVpNrrpml2TTYvLmbx+Xo0SChUJH1622eeMLh6NEQd91VIBSyeMlLbObmghw9\nqnPyZJzt21usW1cln4+wa9cck5PdjI5XXjnP7KyBbTfYsqXF5GSQsbEGk5M6J06E+Pf/fo6NG02G\nh03m5sK0Wg6Dgw3m5mKcOhXEsoJs29ah0XAYHW0yNRWl0wmwbVuLRCKAZdmAQSBgA9rCj97FcTSK\nxRCZTAvDWEz3m06nqdVqSw6jUXlrAC95nT+tcm/+IMfpEAr5Z7Lq3ObOQjoMa2FWqM4H0HHdAGCQ\nTruAErxuuVxuMaGdEq1IJEBfn5pdLp59kUh08wPFYq43a4tGu+1QAx1/OmZVb6/4q1w+/j0z6tjR\n3vw+/kN7YFGM/Qcq+VOAqBxGKlmeP6WzPwWzWhdQ/ehP37xcvv7eg2Z60zSrtvWey3A2g3IpvPMr\n6vN/bPduBmdn+fvXvY6XPvww1z72GF9+05s4euWVF3S92un25b/+a44uGAK/f18JPSymewa86X6j\n0SAajXr58QcGBqhWqzSbTS9zZavVIp1Oe/7VRCJBsVhccsCEclEot0skEmF2dnbJCHhkIX9QqVTy\nDiDn0YQAABWfSURBVL2pVCrecZL5fB7XdRkcHKRarWLbNplMH8ViYSE5GBSLdR59dCM33NAiGi3Q\naCTpdJIkk2X+9V+3sWXLBD/72TbCYZdm0+DQoT5e8pICP/3pIIGAy+7dRTodnccey3LNNQU0DZ54\nIsvVV5c4fDhNs6ljWQEymQ7lcpBk0qLZNOh0ul/uwUGV8rhDJGIzOdn11QYCUCwGueGGEvm8SqNr\nMzLSpq/PZGjIJJdzmJ/XuPHGJidPxpiZ0Xj5yztUq0FyuUVBU4nhVBIzXTcIBBZTOKtDVPz50v2H\ne/h/5L0nV/nPLfaXVz88tRDXm/jLf1CIf/SnHiuxUGco+E8bU+mnAU+slEuqdzbZ+5mWW5fyZ0BV\n9+lN1KYERt1XCVnvecf+8v779vaPX+T8fdorHf7Bkf+95Q4y6R1Mne3AE3VfdcKe6ic//jb3ni/Q\ni7/sckai92jQ3s9zrgNizvV5nw/n67vLdozjpeTgtddiWBZHd+7kiqNHGZ6e5vi2bXz7V35lSbl1\n69YxOTkJLHXXfPDee/nr//AfOLF16xKfvt+9s3v3bk6ePEm5XCaTyXi+6f7+fjZu3MixY8fIZDLs\n2LGDkydPeotqaiqYyWSo1WreiCKRSHh59mdnZzEMA3VilzrBq1wuk0gkPF91d/TdzfnzzDNB1q83\nicV0HnkkxIYNbUZHw/zbv6UZGWlx333bGBxsMTzc4i//civXX19gw4YGf/u3GxkaqmOaBs2mTjbb\nwnW7h4s4DgwNtXFd6HQC3HbbPIcPp9m/f4ZTp+I0m/+/vXMPiuo++/hn7zeBZWFVRC0IeAkVg4KN\nmlcb4yWXJtGmtTPGdjSZ3Ihm0oydOOk0ybRpms7U1DSvRt+pSUzsOx3emWJjM5PEpsZYrVEBYwZF\nQLyAIgisi6zs/ff+gedwWBHBKBL295lx5Cxnz++cZ9nvec7ze37PY+Cuu5r58ksblZVJLFrUxMmT\nJo4dc3LPPY2cPm3E77exYIGX8+f15OdH8XjMHDgguPPOEKGQgYYGyMkJ4vdbMRhCmEwRAgETFRWC\niRP96PWdNcsVD0oRI22pY2VyVxEiRYS0X3yl2Q10CV2st6sIktLgRhFSraAB3Rq+K8SKuNa7hO7z\nNNquXtp9r8a19umLgGibfShtHLX7xHrFsTc45TVtgxLFu1TsrP05tjmJ1u4KsU+8sZLR2zX3dGPp\n77EU8e+tn3FPjsC1zqun49wsbuR4Q6KTV3l+Pm2JiaDT4W5q4vzw4Uw5dIg3n3uu24IIbfNmJTPh\n9rIyHvrwQ/5361aqa2oAmDBhAseOHUOv1zN+/HgaGhpob28nOTmZcePGceTIEVJTU8nIyODIkSMI\nIcjKysLr9aphmVAopNaqT0tL4/z581y6dAmDwU043IbJFCE5ORmj0YLP18ZXX43C5xPMmdPGjh2j\nuHjRybx5Z9i+PZX6+mTmzTvDvn3JVFcnMW3aeYqLx5GQEOaOO5rZvduNwdApnkeP2khKCuNwRLjn\nnhb27HFy771tWK0hdu1KZM4cD3Z7iPHjQ6SkXKKqKonDhxP4yU9a8XiGYbf7cLlAp9MTDoe6PT5D\n5x+f0m4SUEMIOp2uW3ldbSln7WOw9j1a71krNLFt/bTioxDbxUh7fED18LQepyJ2PYmQVsS0X/ZY\nYVdu3lpvuKf47bXS9JT3av+PPYb2dz2NEXueVyMhIaHH8tjaa+1N5PrqlV7tGL2Jdm/jfVPPtye+\nLetzYm9CPdHT7/pjsyHh+ZdPncqJjAxyqqsxB4OUTZ3Kwk8+oXjJkm5llocNG4bdblfDNm63m6Jn\nngHgfzZtYtSoUZSWlmK1WiksLKSsrAy/38/06dOpr6+nqamJxMRE8vLyOHXqFOFwGJfLpaYCOhwO\nUlJSOH78OD6fD7d7DEZjlGDQT3X1SIzGAFu25HL6dCI//nEzc+bU8cQTBdx2WwulpSPJzr5IS4uF\ncFjP7Nle9uxJ4OJFI08+eZp9+5L5+usEfvWrKvbvTyE5OcKyZS3s2OFEr49y//3nKC1NIRwWLFjg\nw2CwoteH1dr6wWBQfcTXPnYrIRAFbZhBWydfmQBW6qYroQnld9rYZqwIaGPIyhhawdcKd6y3FevR\naz2y3ia5tI/32vi19oYS+5jf07lfbWIudhJO+bknQb/a8WLbCF5L9K52I9COqf38tKGOb4PgDQTS\nFl0MCfEvmzaN0mnTeODDD/E5HPxn5kzG1dYS1ev5bN68blk2yoevi0Zxpaby1AsvsPGxx2hJTSU/\nP5/09HT+8Y9/kJqaSnZ2Ng0NDWpMPycnh6qqKvR6PZMmTaKxsVFtW5iVlaVuX7iQRGJihF//+rs0\nN9t55JFGiotH0NRkpqPDxPvvn6a42Monn6SSne1jwYJTVFSk8corJ7lwwcXevRYWLz5PNGqgqspM\nbm7gcpN1M0ZjSG0G397erj5+WywWnE4nHo9HDUUoH4mSux0bPlEEQtvkvSchVwQ/Vqy1/VG1aOPB\n2m3oEi7thJq2Dyt0NeSI7XWq0JfQybW+5Nd6fO7N41J+fzM802vFaa8HKXhdSFt0caPE/5Zm+wid\njsbhw9WCU2GDgbL8fJa/9x57Fy6kQzMhZjQayczMZOTWrZ1ljEMhchcv5t979nDs2DFGjBjBXXfd\nxb59+6itrWXixImcOXMGn89Ha2srBQUF1NfXc+rUKex2O5MnT6a+vp4TJ07x739/hwULgqxZM466\nOgehkIGXX65h2zY3DQ02Pv74MM3NyUyY4CMv7yJr1rTT0nIJl6udRx45jcORgNstyMi4iBCdMe/8\nfAGYVQ89EukUHSWkokw0Kp66snhHWfmnjZvHCpoi9Mp8g/IkcDXRjaUnr1sr+LETYbEetoJ2orSn\n0I72uMqxrhUT14ZwevOo+xJy6e+N55twM48tkdwMbmmRdaHTETGZOO92k372LBGDgRa3m/NuN1mH\nDgGdIR5AbeE2wWZjZGMjHXY7oXCYmTNn4nA4qKioIBqNUlhYiMlk4q9/DeP35zBz5kyi0SgbNkQp\nKxtDXt4UAgErDz44hg8/zCASyWXdukksXvxdjh9P5J136nn++TP84AcB3n23gX37TuB0Gpg4sXOx\nVOeCpQAZGXpGjBihLgry+Xxqdogi3oqAm0wmEhISSEhIUJdoK1lBWnHXhmK0MXcl1GI0GjGZTKpA\nK2sItC3wrtez1T4ZaAVYyV6JXR3aU+int+P25bxi9+vPtcSOE/s+Kc4SSXduuecPUD96dKf4Xxax\nvTNnMn/HDs7eeSdCCPLy8jh8+HBndo3Xy4knn+TrUIiGo0fJycmhsLCQmpoaNm+OkJcHs2dP5qWX\nJvH++1Yefvg0zz6bwfbt4wmHoaQkwMKFDqJRG199lcq77yYzY0Ybr77awY4dIXJyvEyebMViGXY5\nM8WPyWRRe4MqS/rtdjt6vV712JXm5ko2i/YGoA3TKIKthFRMJlM3gVVS4LTxfOVGADcvbHE9DJbz\nkEgk/eeWe/4AtVlZQGehN4DmwkL8VitjP/8cna6znorVaiUQCGBqa6M9OZn0pUupqJjAK6/cxsWL\ngokTJ/Lpp9N56aX/4u23k2huHkZJSR3Hjyfxwx9m09xs5V//amT+fC+/+91osrMvsmnTOd599ySP\nP34Wh6ORH/2oTV1Y5fP5OsczmdQFVkajEYfDQXJysprSqM1/tlqtOBwOtXSCIvSxnrI2C0fx6JXQ\niVLHRfnX0zEkEonkmzIoxP9EZmb314GP772Xu//5T5Scn6lTpzI+MZHh5eWELi/K8nimceBABj/9\n6SQqKy14vQ62bWvkyy/TCIX0jBx5kS1bLvCLXzTw4INnCAQ8PPLIRSoqTvKb3zTQ0dHBbbf5mDpV\nj8vlUlfihsNh7Ha7Wh9HG3/XeupaYY4lNkyiPUZs6ESLDFdIJJKBYFCEfUJmM1t+9jOa3W41fp26\ncCFllZXkrlvH16++ihCCiZcXSpXq0kiqs+Dz6XnttTO0t7exbNlEolE9WVkBPvqolZqazpLN4XCI\nRYscPPRQhI4OGz6fD4sFUlIsmM3DCIfDapjFYrGoqzwB1ZvXroTsLb7dG1LEJRLJYGJQeP4AJ8eN\nQ2hCI5FIhMiLL2IOBkn/742dk6BtbXhnzOBg+Qwef3w6Bw8mkJISYtky+PTTU/zyl7W0tbURCPjJ\nyrKRmpqq1rwJBALYbDbS0tKw2+3qUnuj0YjVau2WOaOEesxms5pj/k2EXyKRSAYbg0b87Zd7sGqb\nfkQMBlo3b+bp/6zl/1Z6MHrbCCcmEgol8fDDbaSnd5CW1kY0GmXMmBCPP27A4XAQCAQIBoOEQiHM\nZjM2m01NuVTi+EpYR5lgVRYe9ZTqGFssSiKRSL7tDIqwTyypmrZ24YQEziTdhrfWie8PglfvP0Vr\nq55589p45RU/gYBOLYcMXaWQldo8iqArpQuU4l5KNUAldn+txUESiUQylBg0nr8QgpkzZ/L119/l\n/fdTGXa5VypA+yUj721vI+h0Mu9fr1NZacLpjKgpl263GyGE2jEJOvPxlawbJZSkLKHXdllSxu6p\nyJhEIpEMVQaN+CvodAV8+ulC/vIXpcQqBAJ6XKNNvPbFOB5b3UEgAC5Xmyro0WgUm82mVtwMBoNq\nhyYllVK7mEobw+8pFVMikUiGOrdW/GO3RWfLu4ceOsv27bPYsiUJn8+A3R5Fr++sGbN0qZ79+5uZ\nMMGI1+slEAhcUfrAbrerE7WxNcN7y5eXwi+RSOKFQRXzF0IQCOiZMSPKqlXnWbEii717Dfh8XROu\nitA7HA51dW1PhcqUnHqJRCKRXMmgCvt0hmz0mM1hxowJ8/HHHqZN8zF5soeOjo5uGTeKR2+xWLDZ\nbN3KKUgkEomkdwaN59/VtUmP2dwZl7dY4A9/CGMwdNDWZqW9vR2bzXblceRErUQikfSLa4r/oUOH\neO+994hGo8ydO5dFixZdsc8777zDoUOHsFgsFBUVkRlTruFqxHr+BoOBQECPxdJVQ14J8yQkJDBs\nWGexNaWblFJzRyKRSCT9o9ewTzQaZfPmzbz44ou88cYb7Nmzh/r6+m77lJWV0djYyJ/+9CeeeOIJ\n/vznP/d58J6yfYJBA5erOHTfV9PgRCmgpi17LL1+iUQi6Tu9in9NTQ0jR45k+PDhGI1GZs2axcGD\nB7vtc/DgQebMmQNATk7O5Y5YF/o0uFb83W735V6yeqzW3tvi9adBiEQikUiupFfxb21tJSUlRd12\nuVy0trb2uk9KSsoV+1yNnj1/PVZrV0tArbBLkZdIJJIbww2Z8L3esgi/OvQbPKeT1W2r1Up9vRW/\nvxWv14jZbCYajaopm8pCLXkTkEgkkm9Gr+LvcrloaWlRt1taWnC5XP3eB6CiooKKigp1e8mSJew7\nf8dVRp7el3MfUiQkJNzqUxg0SFt0IW3RhbRFF8XFxerPubm55Obm9vsYvYZ9srKyOHfuHE1NTYTD\nYfbu3UtBQUG3fQoKCvjiiy8AqKqqwuFw4HQ6rzhWbm4uS5YsUf9pTz7ekbboQtqiC2mLLqQtuigu\nLu6mpdcj/HANz99gMPDoo4/y29/+Vk31HD16NDt27ABg/vz5TJ06lfLyclatWoXVauXpp5++rhOR\nSCQSycBxzZh/fn4++fn53V6bP39+t+3HHnvsxp6VRCKRSG4qt6y8w/U+qgxFpC26kLboQtqiC2mL\nLm6ULXRCdjCRSCSSuEP2JpRIJJI4RIq/RCKRxCEDXtWzL4XihhLNzc2sX78er9eLTqfj7rvv5r77\n7qO9vZ0//vGPNDc343a7+fnPf47D4QCgpKSEnTt3otfrWbFiBVOmTLnFV3FjiUajrFmzBpfLxZo1\na+LWFj6fj40bN6r1soqKikhLS4tLW5SUlLB79250Oh1jx46lqKiIQCAQF7bYsGED5eXlJCYmsnbt\nWoDr+k7U1tayfv16QqEQ+fn5rFixoveBxQASiUTEypUrRWNjowiFQmL16tWirq5uIE9hwPF4POLE\niRNCCCE6OjrEs88+K+rq6sQHH3wgtm3bJoQQoqSkRGzdulUIIURdXZ1YvXq1CIVCorGxUaxcuVJE\nIpFbdfo3he3bt4s333xTvP7660IIEbe2eOutt8Rnn30mhBAiHA4Ln88Xl7ZobGwUzzzzjAgGg0II\nId544w2xc+fOuLHFkSNHRG1trXj++efV1/pz7dFoVAghxJo1a0R1dbUQQojXXntNlJeX9zrugIZ9\n+lIobqjhdDrJyMgAOstXpKen09ra2q0g3ve//30OHDgAwIEDB5g1axZGo5Hhw4czcuRIampqbtXp\n33BaWlooLy9n7ty5almQeLTFpUuXqKysZO7cuUDnmhq73R6XtlDargYCASKRCIFAAJfLFTe2mDRp\nkurVK/Tn2qurq/F4PPj9frKzswGYPXs2+/fv73XcAQ379FQo7tv8ofWXpqYmTp48SU5ODl6vV10J\nnZSUhNfrBcDj8ZCTk6O+pz+F8r4NbNmyhWXLltHR0aG+Fo+2aGpqIjExkQ0bNnDq1CkyMzNZvnx5\nXNpi2LBhPPDAAxQVFWE2m5kyZQp5eXlxaQuF/l670WjsVlanpyKcscgJ3wHC7/ezdu1ali9ffkU3\nsmsVqhsqhexKS0tJTEwkMzPzqsUA48UWkUiEEydOsGDBAn7/+99jtVrZtm1bt33ixRbnzp3jo48+\nYv369WzatAm/36+WjFGIF1v0xM26tgH1/PtaBG6oEQ6HWbt2LbNnz2b69M6idUlJSVy4cAGn04nH\n4yEpKQkY2jY6duwYpaWllJeXEwqF6Ojo4K233opLW6SkpOByudTH9DvuuIOSkhKcTmfc2aK2tpYJ\nEyaohdu+973vUVVVFZe2UOjPd0L5W9J6+n2xyYB6/n0pFDfUEEKwceNG0tPTuf/++9XXCwoK+Pzz\nzwHYtWsXhYWF6ut79uwhHA7T1NTEuXPnVIH4trN06VLefvtt1q9fz3PPPUdubi6rVq2KS1s4nU5S\nU1M5e/YsAIcPH2bMmDFMmzYt7mwxatQoqqurCQaDCCE4fPgwo0ePjktbKPT3O+F0OrHZbFRXVyOE\nYPfu3aqjeTUGfIVveXl5t1TPxYsXD+TwA05lZSUvv/wyY8eOVR/fli5dSnZ29lVTuf72t7+xc+dO\nDAYDy5cv5/bbb7+Vl3BTOHLkCNu3b+eFF17oNa1tKNvi5MmTbNq0iXA4zIgRIygqKiIajcalLf7+\n97+za9cudDodmZmZPPXUU/j9/riwxbp16zh69ChtbW04nU6WLFlCYWFhv69dSfUMBoPk5+fz6KOP\n9jquLO8gkUgkcYic8JVIJJI4RIq/RCKRxCFS/CUSiSQOkeIvkUgkcYgUf4lEIolDpPhLJBJJHCLF\nXyKRSOIQKf4SiUQSh/w/8hNj3RHS460AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1198fa250>"
       ]
      }
     ],
     "prompt_number": 76
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 1000\n",
      "reps = 1000\n",
      "\n",
      "x = g(stats.uniform(0, 1).rvs(n))\n",
      "xb = np.random.choice(x, (n, reps), replace=True)\n",
      "yb = 1/np.arange(1, n+1)[:, None] * np.cumsum(xb, axis=0)\n",
      "upper, lower = np.percentile(yb, [2.5, 97.5], axis=1)\n",
      "\n",
      "plt.plot(xp, yb, c='grey', alpha=0.02)\n",
      "plt.plot(xp ,yb[:, 0], c='red', linewidth=1);\n",
      "plt.plot(xp, upper, 'b', xp, lower, 'b');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEECAYAAAA72gP/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXuUXNV15/+5z3pXd1erW28hqTEgBBiwsJ3AGFsBnHHA\nlmN+GkbxY6KZcUJMwnitWRPLGWyyGOzFiseKY4/GWbP4mZA4ayF7xorHZpghBuMYxj8kQAa3DLgB\nCdCr1e963+fvj1vn9u2iummBkFxif9bq1V1V954691TX+d69z977aGEYhgiCIAjC66Cf6Q4IgiAI\n3YEIhiAIgrAoRDAEQRCERSGCIQiCICwKEQxBEARhUYhgCIIgCIvCfL0D9u/fzz333EMQBGzevJkt\nW7bMef3w4cPs2rWLgwcPctNNN3HDDTfEr1WrVb75zW/y6quvAnDzzTdz3nnnUalU2LlzJ2NjYwwM\nDPDZz36WXC53ii9NEARBOJUsaGEEQcDdd9/N5z//eb761a/y6KOPxpO/olAosH379jlCofjWt77F\nZZddxs6dO/nKV77CqlWrANizZw+XXHIJX/va17jooovYs2fPojo7PDy82Os665GxiJBxmEXGYhYZ\ni1lO5VgsKBgjIyMsW7aMwcFBTNPkyiuvZN++fXOOKRaLDA0NYRjGnOdrtRrPPvssmzdvBsAwDLLZ\nLAD79u3j6quvBuD9738/e/fuXVRn5Z9gFhmLCBmHWWQsZpGxmOVUjsWCLqmJiQn6+/vjx6VSiZGR\nkUU1PDo6SrFYZNeuXRw6dIh169bx+7//+6RSKaanp+nt7QWgp6eH6enpN3EJgiAIwungLVv09n2f\nl156ieuuu4677rqLdDrd0fWkadpb1QVBEAThFLKghVEqlRgfH48fj4+PUyqVFtVwf38/pVKJc889\nF4D3vOc9/MM//AMQWRVTU1P09vYyOTlJT09PxzaGh4fnmFNbt25d1Hu/HZCxiJBxmEXGYhYZi1m2\nbt3K7t2748cbN25k48aNb6itBQVjaGiIY8eOMTo6SqlU4rHHHuPWW2/teGx7DcPe3l6WLFnCkSNH\nWLFiBc8880y86L1p0yZ+/OMfs2XLFh555BGuuOKKjm12urAjR44s+uLOZgqFAuVy+Ux344wj4zCL\njMUsMhazrFix4pQJqPZ61WqfeuqpOWG1H/3oR3nwwQcBuPbaa5mammLHjh3UajV0XSedTrNz507S\n6TQHDx7kr//6r/E8j6VLl/JHf/RHZLPZNxVWK4IRIV+ICBmHWWQsZpGxmGXFihWnrK3XFYxfN0Qw\nIuQLESHjMIuMxSwyFrOcSsGQTG9BEARhUYhgCIIgCItCBEMQBEFYFCIYgiAIwqIQwRAEQRAWhQiG\nIAiCsChEMARBEIRFIYIhCIIgLAoRDEEQBGFRiGAIgiAIi0IEQxAEQVgUIhiCIAjCohDBEARBEBZF\n1wlGlxXXFQRBOGsQwRAEQRAWRdcJRhAEZ7oLgiAIb0u6TjDEwhAEQTgzdJ1guK57prsgCILwtqTr\nBKNer5/pLgiCILwt6TrBaDQa+L5/prshCILwtqPrBMN1XYIgkMVvQRCE00zXCYbvR5FSmqad6a4I\ngiC8reg6wSiXPTzPEwtDEAThNGO+3gH79+/nnnvuIQgCNm/ezJYtW+a8fvjwYXbt2sXBgwe56aab\nuOGGG+LXPvOZz5DJZNB1HcMw+PKXvwzA7t27eeihhygWiwBs27aNSy+9dFEdrlQktFYQBOFMsKBg\nBEHA3XffzW233UapVGLHjh1s2rSJVatWxccUCgW2b9/O3r17O7Zx++23k8/n5zynaRrXX389119/\n/Ul3uFzWCMNQXFKCIAinmQVdUiMjIyxbtozBwUFM0+TKK69k3759c44pFosMDQ1hGEbHNuazBt6o\nlVCpGLiuK1aGIAjCaWZBC2NiYoL+/v74calUYmRkZNGNa5rGHXfcga7rXHPNNVxzzTXxaw888AA/\n+clPWL9+PZ/85CfJ5XKLanNmhjhKaj6REgRBEE49r7uG8Wa444476OvrY2ZmhjvuuIOVK1eyYcMG\nrrvuOm688UYA7rvvPu69915uvvnm15w/PDzM8PBw/Hjr1q00GilSqRT5fB7btt/K7v9aY9s2hULh\nTHfjjCPjMIuMxSwyFnPZvXt3/PfGjRvZuHHjG2pnQcEolUqMj4/Hj8fHxymVSotuvK+vD4jcVu9+\n97sZGRlhw4YN9PT0xMds3ryZu+66q+P5nS5sehrK5TKWZZFOpxfdl7ONQqFAuVw+090448g4zCJj\nMYuMxSyFQoGtW7eekrYWXMMYGhri2LFjjI6O4nkejz32GJs2bep4bPuaQrPZjMt4NBoNnn76adas\nWQPA5ORkfNzjjz8eP78YJiZMwjCUNQxBEITTzIIWhmEYbN++nTvvvDMOq121ahUPPvggANdeey1T\nU1Ps2LGDWq2Gruvcf//97Ny5k+npab7yla8A0ZrDVVddxTvf+U4Avv3tb3Pw4EE0TWNgYIBPf/rT\ni+7wkSMZNE0KEAqCIJxutLDLbtVXrhxn9+7nGRjof0247tsJMbkjZBxmkbGYRcZilhUrVpyytrou\n07taTfPqq1G3JdtbEATh9NF1glEsVhgb0/F9nyAIZC1DEAThNNF1gpHLVThxQsPzPHzfF8EQBEE4\nTXSdYGQyFaamUmhaJBoiGIIgCKeHrhOMbLbM1JQdh9ZKmRBBEITTQ9cJRi5XZWLCxPf9uMy5CIYg\nCMJbT9cJRj5fZWoqheM4sViIYAiCILz1dJ1g5HIVJiZsyuUyzWZTypwLgiCcJrpOMAYGTnDoUBHP\niyKkGo0GIJsqCYIgvNV0nWAUixV03Wd0NI3v+zSbTVzXxff9M901QRCEs5quEwyAlStHOXAgg+u6\nuK5Lo9GI1zNEOARBEN4aulIwli8/wvPP5+LSIEokJGJKEAThraMrBWNw8BX+x/+4kKeeyqLrs5fg\neR4g6xmCIAhvBV0pGOvXvwjAl760Ad/3JSdDEAThNNCVgmHbLrfd9k0aDSMWDN/30TQtLkoolWwF\nQRBOLV0pGACZzBT1usk//mOWZrMZWxgw65qSBXBBEIRTR9cKRhD4fOITT/LDHy7HcRx8339N9ncY\nhmJpCIIgnCK6VjA8z+Oyy37F8PASZmY8Go0G9Xo9LnmedEvJuoYgCMKbp2sFA8CyKqxdO8mTT/bT\naDRwHCfeJ0MJhrI6xNIQBEF4c3S1YIRhyAc+8Dy7d6+h0WhgGAa1Wi1O6GtP5hNLQxAE4Y3T1YJR\nr9e58MKX+OUv+zl2zMTzPOr1+hxLwzCM17ioQJL8BEEQTpauFgzHcbCsOuedd5xPfWozDz+8hEYj\nqi3lOA5hGMbRU8k9wGVtQxAE4eTpasFoNpvU63VuueUB1q6d4j/+x/P4l//yn82xLnRdj8NslXAA\nczLEBUEQhNfHfL0D9u/fzz333EMQBGzevJktW7bMef3w4cPs2rWLgwcPctNNN3HDDTfEr33mM58h\nk8mg6zqGYfDlL38ZgEqlws6dOxkbG2NgYIDPfvaz5HK5k+68Km/e01PnP/2nh/jTP/0ghw/nqFSa\nZLNZGo0GlmXh+z6maWIYBkEQxL+DIEDX9djSmG9vjTAMZd8NQRDe9iwoGEEQcPfdd3PbbbdRKpXY\nsWMHmzZtYtWqVfExhUKB7du3s3fv3o5t3H777eTz+TnP7dmzh0suuYSPfOQj7Nmzhz179vB7v/d7\nb+gCHMcBIAxr/OVf/i++8IXN/J//08fWrQ66ruM4DpqmYVkWuq7HFoamafGCeFIwwjDEMAw8z0PX\n9VgokseIeAiC8HZkQb/MyMgIy5YtY3BwENM0ufLKK9m3b9+cY4rFIkNDQxiG0bGNTusE+/bt4+qr\nrwbg/e9//7xisxiCIMDzPCqVCpVKhY9/fD//7b9dwAMPpKhUKpTLZYIgoFqtUi6XqVQqNJvN2D3V\nbDbjWlSu6+J5Xrxo7rouzWYTIBYbCc8VBOHtyoIWxsTEBP39/fHjUqnEyMjIohvXNI077rgDXde5\n5ppruOaaawCYnp6mt7cXgJ6eHqanp99I32PK5TK5XA7btjn//En+3b/7J772tX/Gj340wYc/fJRL\nL62Qy6WwLIsgCLBtm0ajEVsKykXVHlUFYJpmbIVomkYQBGiaNmcNpJPLSh0n1oggCGcLr7uG8Wa4\n44476OvrY2ZmhjvuuIOVK1eyYcOGOccsNKEODw8zPDwcP966dWvH46JoKQsA13XZuHGST33qRXbu\nvJCnn17CwIDDrl2/pL9fJ5vNxq6mVCpFGIaYpolt27GVpCZ6TdPwPI8wDNF1nVQqFVsj6XT6NS6u\n9t+qPEnyGk+ViNi2TaFQeNPtdDsyDrPIWMwiYzGX3bt3x39v3LiRjRs3vqF2FhSMUqnE+Ph4/Hh8\nfJxSqbToxvv6+oDIbfXud7+bF154gQ0bNtDT08PU1BS9vb1MTk7S09PT8fyTubBqtYqu6wRBQDqd\n5qqrnuPii59lcDDNF75wBddeewW33PICv//7kTWjFsQzmQyZTCaudJsUi6Q7LQgCKpVKbGVUq9U5\nlkb74vnrCUO7wJwshUKBcrn8hs49m5BxmEXGYhYZi1kKhcK8N9sny4JrGENDQxw7dozR0VE8z+Ox\nxx5j06ZNHY9tX6tQIa8QTc5PP/00q1evBmDTpk38+Mc/BuCRRx7hiiuueLPXAUSuKeU+mpqaJJ/3\nCcOA227bx5e+9Bjf+MYQH/zgRbz0ksPMDHHpEPUTBAG1Wo1arRavY5imOWfh2zAMNE2jXq/TbDbj\nNZCk28r3/TjbXK2JqNfVsWofcpUnon46rfkkiykmX5c8EkEQTicLWhiGYbB9+3buvPPOOKx21apV\nPPjggwBce+21TE1NsWPHDmq1Grquc//997Nz506mp6f5yle+AkQT81VXXcU73/lOALZs2cLOnTt5\n+OGH47DaU4XjOPE6xPT0dOx+Wrdumo99bIAHHzyHG298L5YV8IEPnOBP/uQgy5aFpFKp12SDG4ZB\no9FA1/XYklAuKiUkyTUPZeEA8XqJcmclBSFpfSRFRPVVvd98gqDeL+nykrUSQRDearSwy25T//zP\n/3zB19XaRCaTIZ1Ok8/n4zt527bRNJ3R0V4GBx3+9m8v52c/W8pv/uYkl19eYcOGOqmUzvPPZ7nm\nmgqW5WEYBrZtx5N5exhuMqdDiYia8DutYygh6DTZt7evckeSKEEqFotUKpWOUVtJsUl+vElX29ki\nMOJ6mEXGYhYZi1lWrFhxyto66wQjyeDgIK7romkaPT09NJtNDMMgnU5jGAaGYfDLXy7nkUfW8qMf\nrQSgp8dletpC10Pe//5JLr98hq1bp0injVaiYMjPf57ive91sSwrtjiUq0oJiJq0LcuKLQI1+avw\n3KS4KEtFTejJKCslROq9NE2jWCxSLpfnWCqdxKCTaCVf7/aMdzUxLObf+GwRyfmQSXIWGYtZTqVg\nvKVRUmcaVfI8l8tRqVTwPI9CoTCneu2aNSPcfPMx/u2/DdB1m2w2RTptMTkZ8uijq/jud9fwla+s\nZeXKOpOTFoYRUi5brF9f5frrx3nf+yoMDjpkMi6mqceWjJr0LcuKJypVQTeVSgGzFoUSGiUuSmhU\niROI3GOpVGpOqROYDd9NilJSONqTEpO7EiYtIXUMzIrIQtbIQpNv+3snn0+e2/56J9da0lLqdJ6y\n7pLntV9z+/OLZb7rP9uFRxDm46wWjJmZGYA5ZUFUQp6u62QyGRzHoVqttibbBkFg02iY6LrGtde+\nyNVXj/L88xmaTZ0lS6o0GilyuYDh4Sz/+I+r+Ku/WhO/34YNZV5+OcO551Z597uneec7a6xZ41Iq\n+YBOoaBjWTZjYwEjIzqXXeYBPppGXJ5dCYRlWdi2DUSTu+d5NJtNNE0jnU7Hi+ZJkgL0euG+StQ6\nbWObDBc2TXOOi0uJSXIBvtME2u4KS/5OuuKS58636N9+XLuF1G4lzXdspz4uVA5Gnd9J+Noj6URE\nhLcDZ7VgKCqVCgCpVJT9nZxgHMeJ8zFUNJNlWXHCnu/PsHatjmlGQ+X7PsVikXXrKlx99QutCSNF\nuZxh//4VvOMdHiMjNsPDPfzwhwNMT1u4ro7nafT3u2hayOhoZGH09LhUqwbnnNMgCDQ+/vHjTEzY\nTEyYDAzUuOqqaWZmTJYu9VmzRo8tjnq9jmEY1Ov1WAw7rVskt6lVx6gFeLUY335e8nklrsr6USQX\n7zsJQzLUOMnrLdKfbvfYYsKeOx232Ii1TmIz3/ss9vzFiuFbzULXJeJ59vK2EAyFWsOwLAvHcZiZ\nmaFYLNJoNABi91Gj0SCfz1OpVOZkf6vF70qlQhiGsaUShi6pVJmrr57BsiyGhkxuuMGmXp+1GhqN\nLC++CK5rcd55M5imxdNP9/Oe99QZGTE4dCjNgw8uZe1al6NH0xw4UGL37mXYdsjEhMXy5XVAZ8UK\nh3Xr6ui6werVdS69tE4265HLBUxPm6xaZcQuLWVpmKYZl0PxPC+esHVdj4VAWRTt0VxKPBuNRiwm\nahxgdlJLBgWon6R7LClG6jzoLC7td/5Jcep0t7/QxN0+IZ+KiLKF2jgZN1unfi7U3043BEnao/xO\n9jo7WUvzuRcX0047r9eOuP5+/TmrF70Xg5o0TdOMvxypVArDMMhkMnGeRBhGobfNZjOeHDOZTDyZ\nQmStmKaJZVmxyMDspKYyxZXLR0VvZbNZgNhNpPqlLAPXzXL8eA9jYyaNRoqRkSK2bXHwoMFzzxWY\nnrbxPJ102mflyhqbNs0AGv39Dvk8HDtmMjOjk80GgI5lhUDI9DRoGug6TE2ZhGH0mu9rrF7tUiz6\ngMOGDU3WrvXp6/MZGwspFj1sm3hNJekCU8EEyYgv9cU3TTMOCFDFIJUV02nxXllFql31XsmJRy3+\nK0sqaTVpmjZnAk2KUjKooP24Tus/yd/JvxcT1rwYt9XrRa61i2KnCLt8Ph+7YdvpJAILMd/azUJi\n197PZDudxmC+9t7IlNTe106L3icbHXi2CJYsep9CgiCIE/AU1WqVTCYTu6/UWkIy2a6vr49qtRpb\nLLVabc4aie/7ZDKZuE21yKyEKSk8zWaTTCYT17dSda08z2u5yyosX15m+fLIErjoIoeenp54Ed22\nU1SrJpbl8uSTK3n66R6KxQYvvpglCAz6+hoMDHg0mxbpNNRqYBgay5c7NBqg6xYbNtRxXXCcAM8z\nOHIky4kTKWq1DA88YHP8eJqZGYts1qfRMAhDWLeuzJIlTRoNg3zeZcUKl2w2IJsNKRR8TFPHMDQK\nBYd0OuS885rk8x667uP7LpOTJpYVACl6ewNSKSP+kibL0CuRVVWElfWkxlpl+avPKXmnrSwpYI5Y\nJ8OV2ye0ZKHJMAxjd2RSzNXfKoIteZ46ttMkr/qnAh3aXXrJNaLkue3HtIsfRBay67qv6WOyjfY1\no/bovOR4JPOKkiSPTwpWu5jPN77q93xtt49LkvlCxpPnq5uH9oCJ9mufr7/Jthbi18U9eDp52wvG\nfKgsdYBarRZbCMlMb9M047vjZCl15b5xXXeOe6hcLsd1qJRwqC9Is9mMBUq5zNQ/dLPZZGZmhmaz\nGa+teJ5HuVyO9/tQX5Dlyw+xZo0VW02qjUwmE4uj+kLm8/n4bzWZKitAXZf6skXuNwPbttA0g3LZ\n4NVXM4yOWliWRqViMDWVolIxGBszqFYtgiBDrWZQq/XiuhbHjmWp1018H4JAxzQDDCPAsgI0LVrT\nWbmyThCAYYRoWkilYjEzY6FpIZmMj2EEFIsuphkwOOiydKnB0aM2tZqGYYCuh/T1eeRyPoVCSLWq\noWnR5JDLheh6CNix9ZOMLPN9n1wOfD8kk/GwLGg0DEwzIJuFahX6+jx0XcMwQnQ9Ejld9whDD9OM\n+p1K6fHflqURhgGm6WHbJqBhmjpTU1Au6yxdGmBZIdWqi2Ho6Do4Dvh+ZO1lMjqmOVvjLLkGBbOL\n/sq6Ujcu7RM5EOcKJQWzXSjVJJ6MnlOPgTluRnV+0hXZLhxKGOZzKap22y08RbKt9p9kP5XAqe+T\n4zjxTWC7uLX3IzlG7YLdyXpqd/vN5zLsJEbtblk1PqrvyXGfrwL4fJwOZ9FZLRibHn+csYEBDq5b\n96baUR+EilLSdZ1qtRq/rooXqvWAZrMZ79MBkatJubM0TYtLjxQKhfhuWH3Rk+ek02kmJiZi8VHP\nA3E5EyU2qh9qcrdtG9d144xzFYqby+Wo1+uxiADxXbv6YqtJVC2oJwszqvEwDIOlSw36+6PyKckS\nJ0n3kW3beJ5Ho9GIj/P9kDAMCAI/FtPDh3vQdZPjx/PUailyOZ9q1aS3t0lvbx3Litxxvh8yNWUQ\nBAbj43lefrlIoTBDLucAPq5rc+xYGt/XKJdTFIsOvh8COo2GRRiCaWrougmEgJrQQNdtpqej0GnP\nM9A0g1TKw3F0HEfHNEPKZRtd9wkCDd/XW781gsBoexwFOvi+jufphKGGaQY4joGmRcJVKHjU6waW\nFdBsGhhG9H/WaBjYtk8YaniexuBgA8/TgRDDCKnVzJbYqkkzxDAiy9A0l6DrkM36WFYkapYVYJrR\ne1hWgO/TEs5IpJtNg3rdoL/fJZdzqdUsNM1H18F19ZZVCMViQC7XpKfHI5cLCQIT0/SxbbCsgN5e\nHU0LsO2QSsWgtzfA80IajQDfNwAD8DCMENMETQvIZDx832BmJrJai8WQZjMa72PHTNLpkKVLHVxX\np9HQaTQ0HMeg0TAwDB/T9HHdSMTS6QaaBpWKDTTRNI++PpdcDoIgxPdDXBfCUCeVUhZtQBSpGGCa\nGmHoYxiQSpkYBti2Afit8YrGxDCIX9M0H02ba3m2h6onJ/KkOHRar1KPk2t+6vurbu7a69clRTkp\nYMpCT+5f9GY5a9cwrGaTz3/5y0z09fGNW24hPEm1Ft4YapFfiQPMLRGvwnSV9QWzd4mdQnxh7toO\nzJZnAchkMtTrdWzbxnEcbNuOv5SapsWClM1m43phatE+6bpKZtgrlEgrl5Raq/J9n2w2G19jOp2O\nv+CqTpiySCPXmU8Q6EAKy3LQNAPXbWIYOo2G2ZpcXXw/TRhq5HLNlstKJwh0Dh/uIZdz8H0PsMjl\nfILAxDBS6Lq6SdDo7U0zPu7iOAGuq+O6oGkmnqfjeQbpdECzqWMY4HnRBFOpZMnlfCzL5/jxFJ5n\nkM26hKGJ5wUYho/jpAgCqFQsqtUUlUq6ZXlFk7XvG3ieQb1u4TgGrqtTKDSpVu2WlRUJlq4Ti2kQ\naHhedP2m6VMoeK0+65imh2HAwECNRsNgaiqNaQakUj627ce/63UTXQfTDAgCjUYj+pxyOYdsFioV\njenpNK6roTx4hhFZtY4TCXwY6q0+qR/afs/+hKH+mufDEEqlSNBrNYsg0CiVov+TINDQ9bAl6NGP\nGh+AZjO6ecjlPDIZH9fV0bSQMFTCMfu/GP2bRu8XhhqaFl2LbUc3BLY9+2NZIanU7Gt79ohgvC4f\n3rOHy/bvjx+/tHYt3/voRynPUxlXEIS5JCsTLOTrP52035Wr53RdnyPc+Xyeer0eW74wN5AEiCMG\nldDncrlY7FV7MHszoVx5yYoLQWAwPV0gnQbLapBOG4yPZ9A0ME1awhi2rE2NVCpA05qEYUgU6xIJ\nbb1uYVkeYai13L60rHBaCcC0rAwN04yssSAwCQKTZlMjCGwcR8NxdDzPxPMMXNfA80y+852rT9n4\nn5UuKS0IuGz/fn561VVc9dOfArDu4EHOfeEFnrr88jPcO0HoDpIW36+DWMD8C92+789xEyfXIBdL\ne1TVYttIp6M8L9eNflpeY2DWfdUq7jAHNby2Pfec+TATs7XWim5U5y/MqROM7i4kNA/veuIJmrbN\nj665hnv+1b+Kn1967NhJt5VqNLh83z5lEwqCILxtOSstjHUvvsjed78bgENr13LPpz7F2oMHWfPy\ny/ExuueRchzqrRwIhd1s0js5iWvbTPX2suV73+OC557jhh/8gB/+zu9wdPlyeicneXFo6DXnCoIg\nnM2clYIxcOIEP3nf++LHh9at4/jSpbz3Zz/DdBxWv/oqn7z3XgD+/Pbb55z7B9/8JqXJSQAme3uZ\n6u3luzfeyMe++11+54c/BCDQNA5ceCE/+43fYM2hQ6x69VX++403EhgGpfFxqrkczZb/UxCEXz+0\nIMD0PHqnpihOT5OvVvFMkxMDA4wODICmYTsOuWqVbLWK7bqUCwUsx2Gmp4dGOo1vGJFvaLGEIalG\ng0DXcdt9VK3+eIvxTSWuIWz5pTLVKpbn4VgWWhiSbjTI1OvkE266U8FZJxjnPv88tuMwOjg45/lG\nNsvR5cv5rYceon9sLH5+8Phxeqan+Z0f/ID//rGPxWJxYskSfrlhAw9v3gyaxnPnn09xepo/+Ou/\n5hcXXcTlTz3FRcPD1DIZji1bxm/96Ef83/e+lz/++td5Yf16/u4Tnzi5f6YuQm9FP2mA5br0TUxw\n/nPPkanXGe/vZ9+mTQTm/P9ahudhOw71TGb+MQrD17zWMzVF7+QkYwMDVPP5U3U5wskSBFieh+H7\nOJa14GcNUcTiiiNHMH0fo2XZ245DcXqaqd5ePNPEcl0c28axbWzXpXdqinSjQblQINVosGRsjEy9\njm8Y+IbBeH8/xZkZShMTBLqO1+pDJZ8n3WhQy2YJNQ0N0H2fVLNJtlajZ3qaTL2OHgSEmsZ0Tw9T\nvb1Uczlsx2HZsWNkazUM3yfQdaq5HLVsFs80yVcqOLZN79QUtuMQahrVXA6AcqGAa1kEuk660cB2\nHLQwZHRwEMt1WX70KJZKqgwCKvk8rmVheh65ahUtDDE9D8e28UyTSj5PJZ8nU69juW7cn0brRrR/\nfJxUs0kln8fw/fhcuxXOX8tmqWcy1E6xF6Tro6Q03ydsrQAtOXGCz/yX/xJN2J/85GvOveCXv+Rf\n3HcfAGP9/YSaxkBCPCb6+nh11Sq+97GPzd+B1kRmOQ6f+cY32LdpEz+/9FL+8L/+V+qZDL867zze\n8fzz/OqCHaEIAAAgAElEQVS889h/6aW864knePzd72air493Pv00h1euZHTp0lMzGKeJbLXKO371\nK2rZLP3j4/zmo49SaBV0VDx7/vl4psnyo0cxfJ+Da9fy6urV9I+PYzkOlXwe0/fRg4ALh4fpnZ5m\nqqeHSj5PI51mpljk1VWrWDI+TrZa5aJf/ILRwUEObNyI6XkMjYwwODrKVG8vxZkZKvk8WhjyxKZN\nEIb0T0xgOQ6liQn0IODZDRt4ec0aXl69Gtt1cVKpjiKkyM/M0DMzw/KjRymNj1PJ5zm6fDmubVMo\nl6llsxxasyY6X9PQPY9Q1+M7PMKQvokJDN9nrHWHOi8L9GMhtCCgf3ycZUeP0j8xQb5cRgtDxpcs\noZFKseLoUfQgwGtNqHZrgirMzHB82TKOL11Kqtlk2bFjBJrGeH8/1VyOiVIJz7LIVauseuUVCuUy\nhUoFzzCY7u1l3YsvMnDiBHoQYDsOPTMz+LpO0Po5MTDAZF9fdGmaRjWfp29iAt8wKJTLLDt2jNHB\nQXzDiCdhJ5WiXCjQOzmJ5bo00mks18V2HFzLYrq3F1/XydZq1LNZxpYsoZbNovs+tuMwcOIEMz09\njPf3xxOxFgTkKxUamQz5SgVf1zF8H98wqGezlPN5Znp6cCwL17bxEgU1k2Sq1QVfV5+f2bJAAHqn\npzF8H8KQRjodCZimsfT4cRzL4tiyZTTTaRqZDJrv0zc5GYmAYVDJ5wk1DceyyLYEIlOrUahUqGWz\nuJaFZxiYvk+mViPUtPgz75mexjdNpnt65v2f+uIXv3jS/2vz0fWC8Zmvf53DK1ey53d/l42/+AW/\n+eij/P22bVQLhdeeHIZ88c//nJ+95z3873/+z8lUq1zw7LN8+H/+T56++GIueeYZ/t/t23llzZrX\nntuB/MwM1XyeUNf5jUcf5bKnnuKbf/iHFCoV/uCb3yTTaOBYFkdatVyytRqliQkOr1zJQ5s38/La\ntSd17brvc/mTT/Lk5ZcTtOWVaL7PkrEx/p/vfIdGOk0tm+V/fvjDLD96lBeGhrAdh/6xMU4MDOCm\nUuTKZa766U9ZfvQoP/7ABziRuGvXgoArf/pTLjxwAMe2WXHkCBOlEloYcmjtWl5cv55KPk8zlSLQ\nNKZ6e/HVlysMWf/iiyw9fpzzn302mhxsOzK3TZNQ03jm4os5uG4dy44epTQ5Sc/0NKlmk3MOHuT4\nsmVM9/RwYMMG1h46xPIjR9DDkJGhIUbOPRffsuidmGDp8eMEus6FBw7gGwa1XI5qLsd4fz+uaXLR\nL37BuSMjmJ5HtlZjvL+fvslJXlm9GsP3STcaZGs1TgwM4Ns2K19+mVo2y5EVKzi6fDl9k5O841e/\nwjNNyoUCpYkJ0o0GVss1ka3VmCkWAcjU64SahuH7eKaJZ5o8ffnlmI5DoGkMjI5SzecpjY8zODqK\n5bpM9fbSTKV4cf16ji1bRrH1v5StVqN2LIsVhw8zODpKpl6nkU7TPz5OtVDgyLJlTJZKVAsFtDBk\n6ZEjaGHI5MqVOJ5Htl4nW63i2TZ6GNIsFimOjdE/NoaTSjE6OIjp+/SPjpJu3XXrQUBgGLy6di21\nvj4q2Sx6o0FhZoZXL7iAowMDuLpOYNtUenrwLAs3CMg0m/QdO8aS6WkMzwPXJVevM9bfj+n7zORy\nvLxqFa5l4/s66bSeyE4HCHFdp1VlwCIIouifKK/GAGqtRFCzlTdhtpL4ouRNw9BpNtOtpMwKlqVj\nGDpgt0JqfVIpgyAIqVRA19N4XkgQpOjpaWKaUU5KrZZifLynlfvg02hkqFYj91Aq1aBeT+M4JkEA\nmhYlUfq+Qa2WIvJMRWGzoO4pVFUBjXrdbuVNqBwLDdP0gaCVx6IBUZZ/lAgYtRHlt0SJn0EQnVer\nRUmpvq8T2fhhKwQ3qp6gaVESJ4DnRQmg4+OnzhrvesH44u23M9HXx9dvvZUr/+mfyNbrPHjddfOe\nv/allzi8ciVu0lfYGoL+8XHGlyw5Jf0sjY2Rq1Y5vHIlf/JXf8Whc87hex/9KMWZGT7xt3/LkvFx\nAKaLRV5es4ae6Wl6p6b47zfeSD2TYXB0lHUvvcRUTw9jAwMUZma4+JlnWHHkCEYQsGfLFl5at44l\nJ05w6f79XPyLXwDwv377twl0nXc98QTLjh+P7wT1IMAIAgJN41fveAcrjhzh2QsuwHYcLn7mGfQw\n5Il3vYt3PP88pucx2dfHzy+9lHomwyurVzPd29vxOjuVdpjvOE3T5hR5VKXkPc+Ls8qT8f4q/j2d\nTscx86rgoyr0mEqlqNfrcVKgCrEEMHWdgZdfZnLpUpa88grNgQGW/epX1Pr6cPN5AsOgeOQIYaXC\nwd/4DcJUCsdx4qx313XjRMOgHJCt16nrGbyyS351jqWHnsPRLKo9SzCb00wsX4vT1Flz7Hku+Okj\neMUige9zeMUFMNNksn813vnLsfMhxVcO47shK/7vPvRDE7zcO8TYRJ5yvp8j9aVonkdj6RJm8v2k\nsjrVBrjFPnzLol43mZlJkctFCWOapqNp4Dgpmk2NYtFrTSBgmj4nTmTQNI1q1SSd9igUHMrlFLVa\nVKpEI0AnBEND1wPy+QBd93AcgyDQSKU8PM/Etn0yGXBdDd83qFbNVoa1x9hYGs/TyGZ9DCOkWrWo\n1w0cR6fZ1KnXTRxHjyc1lR2dzXqYZkijYbSy3QOCQEfXo2Q804yOdV29lWUdJaOlUtF5YQipVJT4\np2qcRUU0AQLCEDwvStxLp6Ps96hNjUrFjJPj8nmPFSsa6HqI70ePCwWHMNSpVk0KBY9Mxmv9L0cT\ntWVBT09IvQ6WRUsENHw/KvKpxCWfD1qfT4BhaK3KDVFyYCoVtjLOQ8JQb+VsRP2OCnxCKqUD0XuX\nStHzUZZ/2Oq/RhD4mGaYqDSgYRhRMuSWLZed3OS1AF0tGKsPHWL7t74VP973rncxunRpHCH160Jx\naopmby/NxHPrXnyRDQcO0D8+TjOVwrFt3vn00/Hr46USLw4NccXevUz09VGanOTh3/otnr3sMi56\n8kmufPhh9DBkpqeHA+96F89feCG2afLqypUYhkGjXidXreL29pIfG0PXdaazWfobDa7+7nf5+W/9\nFi+dey6+71NwXUzX5QPf/S5jS5dy9JxzeOGCC/BbJUCAOVnTyYleJUOlUqk4o1uVSalWq9TrOpal\nU6loNJspgsBg6VKXMHQxzZDJSZ1SSZVS1+PsaFXmJCp9EjI9naK/v4Hvw+HDecplG8OwgAwvvGBR\nq0VfOsfRmZpK4XnRF73RMAlDlZkbkstFSVCGEaBpIamUT6MRTWi1mkku5zM2lmZ6OtpIK5WKJqhy\nObKgohpXtCbLMM649Twdwwix7ajWVT7vMjmZolaLMntTqSjbWbUT3QlGP+m0T7HosnRpE8sKWbnS\nwTCirGjbDmk2NXwf0mlwHJ3eXo98PqDRANsO4jvX3l6TMKwzPa2i5aPM4qVLm5imQTrt0WzqzMyY\n9PUF5HIuqVSU0R1NdFHpjOlpDceJJjPH0fF9A9v2Wnf2GqmUj2lqZDJeq1yKxuCgg2Vp1GpR5ndv\nr0ehAIbhkU4H5PNg2x6mqWonRYUSy2ULXY/6lkr5eJ6JZYV4XlRrzHU1wlBrJca9tuChSrJTNxSq\nXEYul4uLUiYLUSbraqnnko+TZTcWKrCYLJ6YPFY97lT3Sv0/q/M6FVHsVFYkmWzYiTDxPU32V3Hu\nuefOe+7J0rWCsWR0lM/s2jXntZGhIf6/97yHkfPOO+Xvq2orJavaJklOqOoDTGaN2rYdl8xQmxKp\ngm2qvLlZqTDoOHj5PM1ikcCy0B0HLZVCc10yfX2zReAaDdzJSaqmiZHL4ft+vBmUqjeVrPIahmGr\nfIFLoRB9iV59tYdSKSSddqnVci3T2eaVV3JMT0clJIIAfvGLJeh6yOrVNep1g3I5xeSkDWhUq3Zc\no8iyovpA9bqBaYYcPZrFdaO7vqj8QVRH6MSJVOt5Dcvy8X29VVNIi2sq2XZAJhPVUyqXTdLpIM6W\n7e9vUiq5eJ5OX1/IypVVcrlowkilAvr6ovIPqRQsWeKhadGdne/D5GRUTNDzoru8et0mn3exLJfe\n3pCZGY3BwZAlSwIKBRfXNWg0NPr6XNLp2SqnyqAyDL1VADJ6D9d1OXAgQxiG9Pe7FAohuZzWqlM0\nuzmV+nyUOKq/1f+SOkbVC1OfZfvkkrTsenp64uSz9vpCMH+J9E4VcpPVcJMTYrJWUnKvmPYKscm+\nqqKHnWomJf9OVs9Nfq+S1mPyeXVesv/q+UKhEG+c1qnoYKdJPnmN7RV22wsqto9XJ9onbnV8p2vs\nRPvntVBF3E6vqfNPZS2proyS0nyf9//4x695vm9ykqnW4tupwrbteAJWX65UKjWnIqj6Mqs77mw2\nG9czUpVkHcehWCzGd85BEJDJZGIXinK1qLLnfa1oCFXIT9UtUlu3ptNppqamyLZKfDtOFsOAF1/M\n8fLLecbH04DBSy+VMAyo1QyOHcvhOJGbICqwp5PNRqb8kiXN2M+6enWdFSucuALsBz9YBQJ++csc\na9cGFIsN+vrqrTvgyJ/b0xOZ3UEQmcXlssX55x9j5Uof29bi6rezk21kuqv6QuWyjmH41OsB6XTk\njpiZMbBtyOWiInXT0yHptE4mMzsJ5vN5arVa/Pkkf6svWfLLFASqTpOynCL3la4bcyrXRjWgZhc9\nNS0Vfx6qbfXZJSvI2rbNFVeoG4jMnLtJ9f+R/FyV2y1pySUna7W/iDomWflVtavIZDJxgcv2SrPJ\naq6dJpdObSZFJ3nNyXPaq9omxzb5GSRf7/T7VKP2svl1ZTHX/WbH5q0Y264UjCv27mXjgQMA/O0n\nPsG6l17iqp/+lP6JCaZOQa0oVWU1WchOVZdVPnVVEE/5ulV9Gdu259w15vP5eGc+VVk2DKOS5+l0\nmlwuF99Bzsw06e3N4Hk+IyNLmJmJFvVWrvQ5eDBPtRoyNhZV6jx6tMTLL6cpFn2OHMkwMZEin3cZ\nGqqyZk2NZctcDEPjve+dJJ2GdNpj/XqX/n6HZtND0ywyGb21uBZimnp8J6dKpqvHs+b77EZHybvJ\n5NrDbAntRvw4CLS44F/7FrFAy4UQomlmPAGpyWr2S5+iWDTiSVNNwOozar+TTt4RdtrNL+lSSB6b\nnFiT+2B0uqtMbl2rJk9V0FD1XbWX3Ayq/Y68fWOo9jvmpItD9SNpqSQrA2cymTnXN584JMej/TlB\nmI/XFYz9+/dzzz33EAQBmzdvZsuWLXNeP3z4MLt27eLgwYPcdNNN3HDDDXNeD4KAz33uc5RKJT73\nuc8BsHv3bh566CGKrSiTbdu2cemlly6605lW1VCAai7Hj665hjWHDrHmlVdOKvGlHfWF831/TuGx\nZGXSTCYT3xkC8aSv9tlWoqEmnFQqxeRkmiee6CcMTV5+uUCtptPTEzA6mub4cZtazWRmxmZiIs2a\nNVUmJmwKBZ8VK2o891yBSsWiUHA555wqK1bUWbq0zoc+NMPg4GEmJ20KhYBLLnFJp6PCZGqiUxNi\ntHBrthYTM6+ZQJM/7b5cJX5qMlR3kGpCSm7nmpyokpOpsr7ixeiWlZVc3E62qSYv5QZM9jPpYrNt\nm7yK7Gq7i01OwMlrTVaoVf1IXq86rt2P3C5C87kllEXQ6dzk8UmBm89llKTdauhE+8ZKyb/nQ4RC\nOBkWFIwgCLj77ru57bbbKJVK7Nixg02bNs3xiRUKBbZv387evXs7tnH//fezatWqOYW8NE3j+uuv\n5/rrr39Tnfd1nYlSCYCfXH01W1s5FieD2sxI+f51XY9FwzAMGo0GlmVRaIXpKutAbUSkNlUaG1vC\nc8+twLKq1OsmkOJnPxug0dA5dCjPRRdVKRZd1q+vY9vRhkMXXtjkQx+awTR9Vq70KBRCXnrJZP16\nn8FBryU8R6jVvFZkxKx7KtoR0I73mwBzjgtMrakkfbMqEkldR/udtVrMhtn3UhZWpy1MVdtqElbt\nqeeToqM2dFLnJ90wyXaSk3hSjJKTYfKuvH2RUP1OWkBJlOCczETZfmwnF8tiz13ssTKRC7+OLCgY\nIyMjLFu2jMFW1vSVV17Jvn375ghGsVikWCzy5JNPvub88fFxnnrqKT760Y/ygx/8YM5rb2qtvXXu\n7n/xL+Lw2BfOPZcv/9mfdTy8U+hnoVCIJ3vbtjFNk3Q6jabN3atb3cGq59RdcqPRy0sv9XP4sM0j\nj/Szf38vV189huctIZXymZ5OceONY0xOWnz4w69QLNLabjWMF85VuGi0IB5NyuecA2BiWZn4jriv\nLxVPpOpuXW1+pAQMZu/Uk3fwaiJXrqXkHbUSFyB2RSna77yTwtPuykn61tWx7f5u9VwnN4kSjqRr\nayF/e/u57T7zpFsnaWXIJCwIb44FBWNiYoL+/v74calUYmRkZNGN/83f/A0f//jHO5YJfuCBB/jJ\nT37C+vXr+eQnP0mulWK/GNTXvrlI95OmabE4qOilTCZDLpfDsqw4lFM9ryKZLMtmfLzA1JRNtWow\nPl7kO99Zyfh4Cs+LInre855pNm2a4mtfO0ihMOtfT07eSdeKWuBWloGaJJN39DB7d5284w7DMN66\nNRmFknThtC92JiMyktEtyQidTnRyj3Ryp3SajOfbq7lTH+f7vNrbXCydXDciFIJwanjLFr2feOIJ\nisUi69atY3h4eM5r1113HTfeeCMA9913H/feey8333zza9oYHh6ec+7WrVsB0FqTj7sIwUin05RK\npXibVIgmbeWCSqVS5PN5pqZ0dD3a6evoUZPnn+/jvvtW89JLOdaubVAouBw4kOdf/+sTfPCDhzHN\nNENDNUwzEiPfL8ZbtaoQ16QLSNO0OEpGCUQ6nY4fq8VdteAMszvVtYuAuq5kEpx6j4Umx3Z/+nzH\nJBeEO7maFkNSSN7KO3zbtmN34dsdGYtZZCzmsnv37vjvjRs3snHjxjfUzoKCUSqVGG9lJEPkYiq1\n1gxej+eee44nnniCp556Ctd1qdfrfOMb3+CWW26hJxHJtHnzZu66666Obcx3Yap2TLDABKRpGrlc\nLl6wdRyHMAzjBWzTNCmXTY4dy/Dznw+ya9c7orataHvL3/zNCQoFl2996ynOPbcat6MS1MLQATQM\nIxVHSSUXdtV2oAqVe6HcUCqLuN2N02g0YmtDhUh2ivQxDCPeMKbT5Jyk03Pt4rKQFfHrTKFQeM3G\nN29XZCxmkbGYpVAoxDfbb5YFBWNoaIhjx44xOjpKqVTiscce49Zbb+14bPuktG3bNrZt2wbAgQMH\n+P73v88tt9wCwOTkJH2tfInHH3+cNYus3aRQFRnrLT9+J3K5HKlUCl3X43UItQgchlCv5/njP76K\nyck0555b5t//+5eYnLT41KdOEO2r7MeWiK73xK4sdWdfq9ValkUUfqrWJ9SkrZLxVHilWrRNlrjo\nFIGjxlIJz3wTu2r79VxBqo0kyaioTsd2g1AIgnD6WVAwDMNg+/bt3HnnnXFY7apVq3jwwQcBuPba\na5mammLHjh3UajV0Xef+++9n586dcViqIjkJffvb3+bgwYNomsbAwACf/vSnT6rTluPwDx/5CNPz\nJOmpBex8Ph+vC0QL2hb/+I+r+d//ex0jI0U+8IET3HXXL2g2G/Hkn8/nW5N1FImkXETpdJogCPB9\nP16oVtekFs2T15pctFUL00m3TjISqJ1khNF8d/7tmaqKkwmlTCaBLeY8QRDe3nRlaZCb/v7veery\ny3nuggs6HtPf3x8vaEMU4dRo9PDVr15GuWzzu797lAsvbHL++Q1se7ZmUbTZuhWvQ7RHH1lWVP9G\nuZeUm6o92SxpOSj3UfK15OswN3FMPX49xOSOkHGYRcZiFhmLWVa0qmWfCroy09t2nHkjpCzLIpPJ\nxC6hVCqF49h84QvvoVgM+PrXX2bZMq9V1ykVC4TrunGZBxU9lHQXKetCRVwlw1EVSWugXRDU705C\noRBXkCAIv850rWA4bYKhLIBCoUA2m41dTLVayPe+dyG9vU3+8i+fpbe3hzCczR9Iuo7UgnSz2YxD\nbdtzDYA5YtKpKBzMdR0ln+v0tyAIQjdwVgiGpmlks9k4CkmFox49muGrX70Sw9D5kz85SKGQj+/y\nlTsJZpPfktnDSTdTMiNZvZ96nFygTmYyJ49LIkIhCEK30r2C0cqpUGU9VNSRSrz7z/95E48+upo/\n+qNn2b69hq5raJoeu5XUng0wN5tZ1/V4Ax3gNQl17aKRzGtILmDLIrIgCGcbXSkY6UaDZksw1OK2\nWouYnMzxmc98GIBbbvk5N944haYV4rLjEE3sKrIpWWEU5kY4JbOS5xOA5PPtm7IIgiCcTXSdYOie\nh+l5sWDk8/l4Idp1Df7sz6LtWf/Df3iGj360TjrdE0/+yfIcwByXVHvtoZPJao77JlaFIAhnMV0n\nGNl6PUrYSywqqzIfDz98AZWKzX33/YjBQYPe3mXArCAo4WivsQSzu2BJoTpBEITOdJ1gZGq1OMM7\nk8nEIa4///l6vvOdC/ja155gxYoUqVQqjnpaKNWk3boQoRAEQehM1/lQ8pUK1VZlW8uy4nWGv/u7\n87j11he45JIm2Ww23pwJiBPyOu0r3P63IAiC0JmuE4wVR49ybPlygDiM9tixJVSrFr/925NYlkU2\nm41LeSQFIhkN1W1F9gRBEM40XScYpYkJxpYsAWb3d7j33svZvPkEmhbGloTaJxuIy5oncywUslAt\nCIKwOLputrQcJ46Q0nWdsbGA4eEB/s2/eRnTNMlms3P2kFYJecld2ToJhyAIgrAwXTdjpppNnFYO\nhW3bvPhiiY0bxykWo8fAnL2ok4X/uqzOoiAIwq8VXScYpYkJmq0s7UajyJe+dDVDQ7NVKZO5FUpA\n2ivJCoIgCCdP1wnGkvFxnNbWqs8/P8DAQJWPfeyVOSU81N9qjUNKdgiCILx5ui4PA8BJpfj7v/8Q\njz/+Dj70oZdYuTJaq1DbpypLIllhNll5VhAEQTh5ulIwfF3n8OEoUmpwMNo2Nbmo3R4yK0IhCILw\n5ulKwZjq7SWXc7nllp+zefMklpWf44aC2ZpQ6m8RDEEQhDdH1zn0jyxfToDBK6/0ctFFkyxbFpUB\nSRYLnG+/a0EQBOGN03UWxujgIM89dz7ptM8554SvEYo3UmVWEARBeH26zsL4/oc/zKFDQ/z2bx+i\ntzcT14dK5lwky4CIO0oQBOHU0HWCERoGTz99MR/60HhcklxlcysXVNIVJWIhCIJwaug6l5TrGjiO\nydq1kQvKNM14xz3ovHWqIAiC8OZ5XcHYv38/99xzD0EQsHnzZrZs2TLn9cOHD7Nr1y4OHjzITTfd\nxA033DDn9SAI+NznPkepVOJzn/scAJVKhZ07dzI2NsbAwACf/exnybVKlr8e1WqeYrGBYejx5kni\nfhIEQXjrWdAlFQQBd999N5///Of56le/yqOPPsqrr74655hCocD27dtfIxSK+++/n1WrVs2ZyPfs\n2cMll1zC1772NS666CL27Nmz6A5XKjn6+prx7nkqnLaTWIh4CIIgnDoWFIyRkRGWLVvG4OAgpmly\n5ZVXsm/fvjnHFItFhoaG5iTOKcbHx3nqqafYvHnznHWFffv2cfXVVwPw/ve/n7179y66w5VKnt5e\nB+gsCMmIKUEQBOHUseCsOjExQX9/f/y4VCoxMTGx6Mb/5m/+ho9//OOvmbynp6fp7e0FoKenh+np\n6UW3Wa3mKZVcTNN8TaJee4SUIAiCcOp4yxa9n3jiCYrFIuvWrWN4eHje4xZyGw0PD885d+vWrXhe\nicHBgEKhQC6XI51OxxsmvZ2EwrZtCoXCme7GGUfGYRYZi1lkLOaye/fu+O+NGzeycePGN9TOgoJR\nKpUYHx+PH4+Pj1MqlRbV8HPPPccTTzzBU089heu61Ot1vvGNb3DLLbfQ09PD1NQUvb29TE5O0tPT\n07GNThc2PZ1h9eoKtVoDwzDwPA/Lst52lkWhUKBcLr/+gWc5Mg6zyFjMImMxS6FQYOvWraekrQUF\nY2hoiGPHjjE6OkqpVOKxxx7j1ltv7XhsexmObdu2sW3bNgAOHDjA97//fW655RYANm3axI9//GO2\nbNnCI488whVXXLHoDlereZYsacSP1cK3IAiC8NayoGAYhsH27du5884747DaVatW8eCDDwJw7bXX\nMjU1xY4dO6jVaui6zv3338/OnTtJp9Nz2kq6nrZs2cLOnTt5+OGH47DaxVIu5xgcrGOapkRBCYIg\nnEa0sMsq9C1dOslf/MUIl16aIp1Ok0ql4qS9txNickfIOMwiYzGLjMUsK1asOGVtdZ3Tv1zOUiq5\nss+FIAjCaabrBMNxTPJ5F+BtFxklCIJwJum62bZYbKDrmix2C4IgnGa6TjB6eppxKXOxLgRBEE4f\nXTfj9vQ0sG07jpIS0RAEQTg9dN1s29fnxBslCYIgCKePrhMMFSElCIIgnF66TjD6+lxxRwmCIJwB\num7GLRSCOVVqBUEQhNND1826mUwAyOZIgiAIp5uuE4xsNqpkIjkYgiAIp5euEwxlYQiCIAinl64T\njHTaj/fwFgRBEE4fXSkYIhaCIAinHxEMQRAEYVF0nWCkUiIYgiAIZ4KuE4x83pMcDEEQhDNA1828\nuZxszSoIgnAm6DrBMM0FtyEXBEEQ3iK6TjCkLIggCMKZoetmXiltLgiCcGboOsEIw/BMd0EQBOFt\nSdcJhiAIgnBmeN0V5P3793PPPfcQBAGbN29my5Ytc14/fPgwu3bt4uDBg9x0003ccMMNADiOw+23\n347runiexxVXXMG2bdsA2L17Nw899BDFYhGAbdu2cemlly6uw6YpaxiCIAhngAUFIwgC7r77bm67\n7TZKpRI7duxg06ZNrFq1Kj6mUCiwfft29u7dO+dc27b54he/SCqVwvd9vvCFL/Dss89ywQUXoGka\n1y1BHCEAAA2pSURBVF9/Pddff/1Jd9gwDKlUKwiCcAZY8FZ9ZGSEZcuWMTg4iGmaXHnllezbt2/O\nMcVikaGhoY6TeCqVAsDzPIIgIJ/Px6+90bUIsS4EQRDODAtaGBMTE/T398ePS6USIyMji248CAL+\n9E//lOPHj3PdddfNsUweeOABfvKTn7B+/Xo++clPksvlFtWmCIYgCMKZ4S3NgtN1nb/4i7+gVqtx\n5513Mjw8zMaNG7nuuuu48cYbAbjvvvu49957ufnmm19z/vDwMMPDw/HjrVu3ks/nyefzb3u3lG3b\nFAqFM92NM46MwywyFrPIWMxl9+7d8d8bN25k48aNb6idBQWjVCoxPj4ePx4fH6dUKp30m2SzWS67\n7DJeeOEFNm7cSE9PT/za5s2bueuuuzqe1+nCqtUqmqa97QWjUChQLpfPdDfOODIOs8hYzCJjMUuh\nUGDr1q2npK0F/TtDQ0McO3aM0dFRPM/jscceY9OmTR2PbV+TmJmZoVqtAlHE1DPPPMO6desAmJyc\njI97/PHHWbNmzaI7LFFSgiAIZ4YFLQzDMNi+fTt33nlnHFa7atUqHnzwQQCuvfZapqam2LFjB7Va\nDV3Xuf/++9m5cyeTk5Ps2rWLIAgIw5D3ve99XHzxxQB8+9vf5uDBg2iaxsDAAJ/+9KcX3WHJ8hYE\nQTgzaGGXpU6PjIyQzWbPdDfOOGJyR8g4zCJjMYuMxSwrVqw4ZW11nW9HLAxBEIQzgwiGIAiCsCi6\nTjBkwVsQBOHM0HWz79s9nFYQBOFM0XWCIQiCIJwZRDAEQRCERSGCIQiCICwKEQxBEARhUXSdYEiU\nlCAIwplBZl9BEARhUXSdYEjiniAIwpmh6wRDEARBODOIYAiCIAiLQgRDEARBWBQiGIIgCMKiEMEQ\nBEEQFoUIhiAIgrAoRDAEQRCERSGCIQiCICwKEQxBEARhUYhgCIIgCItCBEMQBEFYFCIYgiAIwqIw\nX++A/fv3c8899xAEAZs3b2bLli1zXj98+DC7du3i4MGD3HTTTdxwww0AOI7D7bffjuu6eJ7HFVdc\nwbZt2wCoVCrs3LmTsbExBgYG+OxnP0sul3sLLk8QBEE4VSwoGEEQcPfdd3PbbbdRKpXYsWMHmzZt\nYtWqVfExhUKB7du3s3fv3jnn2rbNF7/4RVKpFL7v84UvfIFnn32WCy64gD179nDJJZfwkY98hD17\n9rBnzx5+7/d+7625QkEQBOH/b+9+YppIwzCAP6WNImvbcVAWAQmIRk3jYrNFTYhoMHDwT+JeGkOM\nIgdjkN3oxggcjCc26xpEo6AeTPx7kINojMbEuIhED4qp0dR1hWANaqCxlLoLFCj99kCc2hXLoCCF\neX6JSWc6nZnvie2br515GRMRv5JqaWlBYmIiEhISYDAYkJ2djaamprBtTCYTMjIyoNfrP3n99OnT\nAQCBQADBYBAzZ84EADQ1NWH16tUAgDVr1nxSbIiIKPpEnGF0dnYiPj5eWZZlGS0tLap3HgwGUVpa\nio6ODuTn5yszE5/PB0mSAABmsxk+n+9Lzp2IiL6hEX/D+BoxMTE4dOgQenp6UFFRAafTCYvFErZN\npD+I5HQ64XQ6lWW73Q6j0Thu5zuZTJs2jVmAOXyMWYQwi3C1tbXKY4vF8snnsFoRC4Ysy/B4PMqy\nx+OBLMujPkhcXBysVitaW1thsVhgNpvR1dUFSZLg9XphNpuHfd1wA/vnn39GffypyGg0Mgswh48x\nixBmEWI0GmG328dkXxF/w8jIyEB7ezvcbjcCgQDu378Pm8027LZCiLDl9+/fo7u7G8DQFVNPnz5F\nWloaAMBms+HOnTsAgIaGBmRlZX3lMIiIaLxFnGHo9XoUFRWhoqJCuaw2JSUFt27dAgDk5eWhq6sL\n5eXl6OnpQUxMDG7cuIGqqip4vV7U1NQgGAxCCIGcnBwsXboUALBp0yZUVVWhvr5euayWiIiim078\nf2oQ5d6+fTvRpxAVOOUewhxCmEUIswhJSkoas33xTm8iIlKFBYOIiFRhwSAiIlVYMIiISBUWDCIi\nUoUFg4iIVGHBICIiVVgwiIhIFRYMIiJShQWDiIhUYcEgIiJVWDCIiEgVFgwiIlKFBYOIiFRhwSAi\nIlVYMIiISBUWDCIiUoUFg4iIVGHBICIiVVgwiIhIFRYMIiJShQWDiIhUMYy0wePHj3HmzBkEg0Hk\n5uZi06ZNYc+/efMGNTU1cLlc2Lx5MzZu3AgAePfuHaqrq+Hz+aDT6bB27VqsW7cOAFBbW4s///wT\nJpMJAFBQUIBly5aN9diIiGgMRSwYwWAQp0+fxv79+yHLMsrLy2Gz2ZCSkqJsYzQaUVRUhIcPH4bv\n2GDAtm3bkJaWBr/fj9LSUmRmZiI5ORk6nQ4bNmzAhg0bxmdUREQ05iJ+JdXS0oLExEQkJCTAYDAg\nOzsbTU1NYduYTCZkZGRAr9eHrZckCWlpaQCA2NhYJCcno7OzU3leCDFGQyAiom8h4gyjs7MT8fHx\nyrIsy2hpaRn1QdxuN1wuFxYuXKisu3nzJu7evYv58+dj69at+O6770a9XyIi+nbG/Udvv9+Pw4cP\no7CwELGxsQCA/Px8HD9+HH/88QdmzZqFc+fOjfdpEBHRV4o4w5BlGR6PR1n2eDyQZVn1zgOBACor\nK7Fq1SosX75cWW82m5XHubm5OHjw4LCvdzqdcDqdyrLdbkdSUpLq4091RqNxok8hKjCHEGYRwixC\namtrlccWiwUWi+WL9hNxhpGRkYH29na43W4EAgHcv38fNptt2G3//5uEEAInT55EcnIy1q9fH/ac\n1+tVHj948ACpqanD7tNiscButyv/Ph601jGLIcwhhFmEMIuQ2trasM/RLy0WwAgzDL1ej6KiIlRU\nVCiX1aakpODWrVsAgLy8PHR1daG8vBw9PT2IiYnBjRs3UFVVBZfLhcbGRqSmpmLfvn0AQpfPXrx4\nES6XCzqdDnPmzMGOHTu+eABERPRtjHgfhtVqhdVqDVuXl5enPJYkCSdOnPjkdYsXL8alS5eG3WdJ\nScloz5OIiCbYpLrT+2umUlMNsxjCHEKYRQizCBnLLHSCN0QQEZEKk2qGQUREE4cFg4iIVBnxR+9o\nMFIDxKnmc40b//33X1RVVeHdu3eYM2cO9uzZo9whX1dXh/r6esTExGD79u3IzMyc4FGMrWAwiLKy\nMsiyjLKyMs1m0d3djZMnT+L169cAgOLiYsydO1eTWdTV1aGxsRE6nQ6pqakoLi5GX1+fJrKoqamB\nw+GAyWRCZWUlAHzRe6K1tRXV1dUYGBiA1WrF9u3bIx9YRLnBwUFRUlIiOjo6xMDAgNi7d69oa2ub\n6NMaV16vV7x8+VIIIURvb6/45ZdfRFtbmzh//ry4cuWKEEKIuro6ceHCBSGEEG1tbWLv3r1iYGBA\ndHR0iJKSEjE4ODhRpz8url27Jo4ePSp+//13IYTQbBbHjh0Tt2/fFkIIEQgERHd3tyaz6OjoELt2\n7RL9/f1CCCEOHz4s6uvrNZPFs2fPRGtrq/j111+VdaMZezAYFEIIUVZWJpqbm4UQQvz222/C4XBE\nPG7UfyWlpgHiVPO5xo1NTU1YvXo1AGDNmjVKh+CHDx8iOzsbBoMBCQkJSExM/KKeX9HK4/HA4XAg\nNzdXuUFUi1n09PTg+fPnyM3NBTB0n1RcXJwms4iLi4Ner0dfXx8GBwfR19cHWZY1k8WSJUs+6b83\nmrE3NzfD6/XC7/djwYIFAICcnBw8ePAg4nGj/iupsWqAOFl93LjR5/NBkiQAQ+1VfD4fgKE75z9u\n7BgfHx/WGXiyO3v2LLZs2YLe3l5lnRazcLvdMJlMqKmpwatXr5Ceno7CwkJNZjFz5kxs3LgRxcXF\nmDZtGjIzM/HDDz9oMosPRjt2g8EQ1upJluURM4n6GYaW+f1+VFZWorCwEDNmzAh7TqfTRXztSM9P\nFo8ePYLJZEJ6evpnW+JrJYvBwUG8fPkS+fn5OHjwIGJjY3HlypWwbbSSRXt7O65fv47q6mqcOnUK\nfr8fd+/eDdtGK1kMZ7zGFvUzjK9tgDhZfWjcmJOTozRuNJvN6OrqgiRJ8Hq9ShPHqZzR33//jUeP\nHsHhcGBgYAC9vb04duyYJrOIj4+HLMvKVwgrV65EXV0dJEnSXBatra1YtGiR0mBwxYoVePHihSaz\n+GA074kP/5c+nlGoySTqZxijaYA4VYjPNG602Wy4c+cOAKChoQFZWVnK+nv37iEQCMDtdqO9vV35\nUJnsCgoKcOLECVRXV2P37t2wWCz4+eefNZmFJEmYPXs23r59CwB48uQJ5s2bhx9//FFzWSQlJaG5\nuRn9/f0QQuDJkydISUnRZBYfjPY9IUkSZsyYgebmZggh0NjYGNZVfDiT4k5vh8MRdlntTz/9NNGn\nNK6eP3+OAwcOIDU1VZlaFhQUYMGCBZ+9bO7y5cuor6+HXq9HYWHhlPwb6c+ePcO1a9dQWloa8RLC\nqZyFy+XCqVOnEAgE8P3336O4uBjBYFCTWVy9ehUNDQ3Q6XRIT0/Hzp074ff7NZHFkSNH8Ndff+H9\n+/eQJAl2ux1ZWVmjHvuHy2r7+/thtVpRVFQU8biTomAQEdHEi/qvpIiIKDqwYBARkSosGEREpAoL\nBhERqcKCQUREqrBgEBGRKiwYRESkCgsGERGp8h/J73/CT+G+4wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x113cbecd0>"
       ]
      }
     ],
     "prompt_number": 78
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Exercise 5 (10 points)**.\n",
      "\n",
      "Estimate the following integral using Monte Carlo integration\n",
      "\n",
      "$$\n",
      "\\int_{-\\infty}^{\\infty} x^2 \\frac{1}{2}e^{-|x|} dx\n",
      "$$\n",
      "\n",
      "Hint: See figure.\n",
      "\n",
      "![Hint](./dexp.png)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Estimate using importance sampling"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def p(x):\n",
      "    \"\"\"Double exponential density.\"\"\"\n",
      "    return 0.5*np.exp(-np.abs(x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 10000\n",
      "x = stats.norm(0, 2).rvs(n)\n",
      "np.mean(x**2 * p(x)/stats.norm(0, 2).pdf(x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "2.0242"
       ]
      }
     ],
     "prompt_number": 14
    }
   ],
   "metadata": {}
  }
 ]
}