scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy\n",
    "\n",
    "\n",
    "<a id='index-1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [SciPy](#SciPy)  \n",
    "  - [SciPy versus NumPy](#SciPy-versus-NumPy)  \n",
    "  - [Statistics](#Statistics)  \n",
    "  - [Roots and Fixed Points](#Roots-and-Fixed-Points)  \n",
    "  - [Optimization](#Optimization)  \n",
    "  - [Integration](#Integration)  \n",
    "  - [Linear Algebra](#Linear-Algebra)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[SciPy](http://www.scipy.org) builds on top of NumPy to provide common tools for scientific programming, such as\n",
    "\n",
    "- [linear algebra](http://docs.scipy.org/doc/scipy/reference/linalg.html)  \n",
    "- [numerical integration](http://docs.scipy.org/doc/scipy/reference/integrate.html)  \n",
    "- [interpolation](http://docs.scipy.org/doc/scipy/reference/interpolate.html)  \n",
    "- [optimization](http://docs.scipy.org/doc/scipy/reference/optimize.html)  \n",
    "- [distributions and random number generation](http://docs.scipy.org/doc/scipy/reference/stats.html)  \n",
    "- [signal processing](http://docs.scipy.org/doc/scipy/reference/signal.html)  \n",
    "- etc., etc  \n",
    "\n",
    "\n",
    "Like NumPy, SciPy is stable, mature and widely used\n",
    "\n",
    "Many SciPy routines are thin wrappers around industry-standard Fortran libraries such as LAPACK, BLAS, etc.\n",
    "\n",
    "It’s not really necessary to “learn” SciPy as a whole\n",
    "\n",
    "A more common approach is to get some idea of what’s in the library and then look up [documentation](http://docs.scipy.org/doc/scipy/reference/index.html) as required\n",
    "\n",
    "In this lecture we aim only to highlight some useful parts of the package"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SciPy versus NumPy\n",
    "\n",
    "SciPy is a package that contains various tools that are built on top of NumPy, using its array data type and related functionality\n",
    "\n",
    "In fact, when we import SciPy we also get NumPy, as can be seen from the SciPy initialization file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Import numpy symbols to scipy name space\n",
    "import numpy as _num\n",
    "linalg = None\n",
    "from numpy import *\n",
    "from numpy.random import rand, randn\n",
    "from numpy.fft import fft, ifft\n",
    "from numpy.lib.scimath import *\n",
    "\n",
    "__all__  = []\n",
    "__all__ += _num.__all__\n",
    "__all__ += ['randn', 'rand', 'fft', 'ifft']\n",
    "\n",
    "del _num\n",
    "# Remove the linalg imported from numpy so that the scipy.linalg package can be\n",
    "# imported.\n",
    "del linalg\n",
    "__all__.remove('linalg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, it’s more common and better practice to use NumPy functionality explicitly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "a = np.identity(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is useful in SciPy is the functionality in its subpackages\n",
    "\n",
    "- `scipy.optimize`, `scipy.integrate`, `scipy.stats`, etc.  \n",
    "\n",
    "\n",
    "These subpackages and their attributes need to be imported separately"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "from scipy.optimize import brentq\n",
    "# etc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s explore some of the major subpackages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistics\n",
    "\n",
    "\n",
    "<a id='index-4'></a>\n",
    "The `scipy.stats` subpackage supplies\n",
    "\n",
    "- numerous random variable objects (densities, cumulative distributions, random sampling, etc.)  \n",
    "- some estimation procedures  \n",
    "- some statistical tests  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Variables and Distributions\n",
    "\n",
    "Recall that `numpy.random` provides functions for generating random variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.55669544, 0.50908178, 0.63340083])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.beta(5, 5, size=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This generates a draw from the distribution below when `a, b = 5, 5`\n",
    "\n",
    "\n",
    "<a id='equation-betadist2'></a>\n",
    "$$\n",
    "f(x; a, b) = \\frac{x^{(a - 1)} (1 - x)^{(b - 1)}}\n",
    "    {\\int_0^1 u^{(a - 1)} (1 - u)^{(b - 1)} du}\n",
    "    \\qquad (0 \\leq x \\leq 1) \\tag{1}\n",
    "$$\n",
    "\n",
    "Sometimes we need access to the density itself, or the cdf, the quantiles, etc.\n",
    "\n",
    "For this we can use `scipy.stats`, which provides all of this functionality as well as random number generation in a single consistent interface\n",
    "\n",
    "Here’s an example of usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/quantecon/anaconda3/lib/python3.7/site-packages/matplotlib/axes/_axes.py:6571: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAFpCAYAAABeYWb6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmcjvX+x/HX977v2ZgVY08oqWhREwmd0CJKnKPQaVEqbSralPPTeip1aFFx2kQL6lTqFJUSIjmWUNkrpOwzY2aM2b+/P2aOg4a5h3vmey/v5+NxP5qZ+3Lf77qa8Z7PdV3fy1hrEREREZEj53EdQERERCRcqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIjP1RvXqVPHNm3a1NXbi4iIiPht8eLFO6y1qRVt56xYNW3alEWLFrl6exERERG/GWM2+LOdDgWKiIiIBIiKlYiIiEiAqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIiKlYiIiEiAqFiJiIiIBIjPdQAREVeaDvvE723XP9GjCpOISLjQxEpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAJExUpEREQkQFSsRERERAKkwmJljDnKGPOVMWalMeZHY8zt5WxzjjFmlzFmadljRNXEFREREQlePj+2KQLutNYuMcYkAIuNMTOstSsO2O5ra+1FgY8oIiIiEhoqLFbW2s3A5rKPs40xK4FGwIHFSkQk4JoO+8Tvbdc/0aMKk4iIVKxS51gZY5oCbYAF5Tzd3hizzBgz3RjTKgDZREREREKKP4cCATDGxAPvAXdYa7MOeHoJcLS1NscY0x2YCrQo5zVuAG4AaNKkyWGHFhEREQlGfk2sjDFRlJaqt6y17x/4vLU2y1qbU/bxNCDKGFOnnO1estamWWvTUlNTjzC6iIiISHDx56pAA7wKrLTWjj7INvXLtsMY07bsdXcGMqiIiIhIsPPnUGAH4Erge2PM0rKv3Q80AbDWjgP6ADcZY4qAPUA/a62tgrwiInvZ4kIK03+jKP13bHEhkydn7/d8o0aNOPHEE6ldu7ajhCISafy5KnAuYCrY5nng+UCFEhE5kLWWgs1r2PPzYgq3r6dgx0aKMn4HW7J3m/7/Lv/P1qtXj1atWtGqVSs6dOhAjx49iI+Pr6bkIhJJ/D55XUSkullryf99Nbmr5rJ79TyKs7YdsIXBl9yAqNqNMVGxXHxKw73PFBcXs379elasWMHWrVvZunUrM2fOZMyYMcTGxtK9e3d2Fx5D3LFt8UTHVe+/mIiELRUrEQk6OTk5vPjii7z44ots2bBh79e98bWocdxZRDdsSXSdJvhqNcITFbv3+cnlrGNVUlLCxo0b+fHHH1m+fDkff/wx33zzDe+/X3odjvFFE9eiPUntLyU6tWmV/7uJSHgzrk6FSktLs4sWLXLy3iISnLKzs3n++ecZNWoUO3eWXv/ija9FjZYdqHF8R2IanYAxB7/mxt8FQjdt2sR7773HvU+9RP5v/1vruEbLDiR16F9uwdLioyKRzRiz2FqbVtF2mliJiHPZ2dmMGTOGUaNGkZ6eDkD79u0ZMWIEg74qPGSZOhyNGzfm9ttv5+nNx1KUtY2sBe+TvewzclfPI3f1PGocd1ZpwarbLKDvKyLhL7A/rUREKmn69OmceOKJDB8+nPT0dDp27MiMGTOYN28e3bp1C3ipOpAvsS61zruRRoNeIeH0izG+aHLXfMPm128n46vXKCnMr9L3F5HwoomViDiRnp7O0KFDmTBhAgBpaWmMHDmSzp07U7YsXqVV5r6CB/Il1KbWuYNIPPNSsr59l+wln5D1n/fJXfsttS+8DdChQBGpmCZWIlLtpk6dSqtWrZgwYQKxsbE89dRTzJ8/ny5duhx2qQoUX3wtap07iPpX/oOoOk0oyvidrW8PY/DgweTk5DjNJiLBT8VKRKpNbm4uV111Fb1792bLli107NiRZcuWcdddd+HzBdcAPabBcTS4+lmSzuoPHi/PP/88J510EosXL3YdTUSCmIqViFSLjRs30rFjR9544w1q1KjBc889x+zZsznuuONcRzso44siudNfaXD107Rp04b169fTsWNH3n77bdfRRCRIqViJSJX7+uuvSUtL47vvvqN58+YsWLCAwYMH4/GExo+g6LrNmT9/PgMHDiQvL4+//vWv3HPPPRQXF7uOJiJBJjR+qolIyBo3bhxdunRh+/btnHfeeSxcuJDWrVu7jlVpMTExvPzyyzz//PP4fD6eeuopLrroIjIyMlxHE5EgElwnNYhI2CguLmbw4MGMHTsWgKFDhzJy5MigO5eqMowx3HLLLbRq1YpLL72UTz/9lHbt2jF9+nSOOeaYw3rNyl7JqIVKRYKbJlYiEnCFhYVcddVVjB07lpiYGCZOnMioUaNCulTt65xzzmHhwoWccsoprF27lk6dOrFy5UrXsUQkCKhYiUhAFRQU0LdvX95++23i4+P5/PPPufLKK13HCrimTZsyd+5czjnnHDZv3syf/vQnli1b5jqWiDimYiUiAbNnzx569+7NBx98QHJyMl988QVnn32261hVJj4+nmnTptGtWze2b99O586dWbhwoetYIuKQipWIBMTu3bu56KKLmDZtGnXq1OGrr76iXbt2rmNVubi4OKZOnUqvXr3IyMiga9euzJ0713UsEXFExUpEjtju3bvp1q0bM2fOpH79+syaNYtTTz3VdaxqExMTwzvvvEO/fv3Izs7mggsuYPbs2a5jiYgDKlYickQKCwu59NJLmTt3Lo0bN2bOnDm0atXKdaxqFxUVxZtvvsk111xDbm4uPXv21DlXIhEoPC7RERG/VObSfn8u67fWct111zF9+nQ8cYmUXPh/nPfqGmDNEb92KPJ6vbz88svs3r2bd955h27duvHNN9/QrFkz19FEpJpoYiUih23YsGFMnDgRExVL3T4PEFWrketIznm9XiZOnEiXLl3YsmUL559/Ptu2bXMdS0SqiYqViByW0aNH8+STT+Lz+UjtdR8xDVu6jhQ0YmJi+OCDD2jTpg3r1q2jR48eZGdnu44lItVAxUpEKu2tt97izjvvBGD8+PHENT/dcaLgk5iYyLRp02jevDmLFi3iL3/5CwUFBa5jiUgVU7ESkUqZPXs2AwYMAGDUqFFcccUVbgMFsfr16/PZZ59Rt25dZsyYwaBBg7DWuo4lIlVIxUpE/LZx40b69OlDUVERQ4YMYejQoa4jBb1jjz2W6dOnExcXx+uvv86YMWNcRxKRKqRiJSJ+yc3NpVevXuzYsYPzzz+fp556ynWkkHHaaacxfvx4oPRm1DNnznScSESqioqViFTov8sqfPfddxxzzDFMnjwZr9frOlZI6du3L8OGDaO4uJjLLruM9evXu44kIlVA61iJSIVGjRrFpEmTiI+P58MPPyQlJeWIXq8y62mFk0cffZRly5Yxffp0evXqxbx581xHEpEA08RKRA7ps88+49577wVg4sSJEbmqeqB4vV7efvttWrRowbJlyxg4cKBOZhcJMypWInJQP//8M/369aOkpIQRI0bQu3dv15FCXnJyMh9++CEJCQlMmTKFrP+85zqSiASQipWIlKugoIB+/fqRmZlJz549eeCBB1xHChsnnHACb7zxBgCZsyeS/9tKx4lEJFBUrESkXMOHD2fhwoUcffTRvP7663g8+nERSJdccgl33XUX2BK2f/QUJXk5riOJSADoJ6WI/MGenxfzj3/8A6/Xy6RJk474ZHUp39///neiG7SgOGsbO6c/p/OtRMKAipWI7KcoJ50dn4wG4OGHH6Z9+/aOE4Wv6Oho6lx8DyY6jtw135Cz7FPXkUTkCGm5BRHZy9oSdn48mpLcXcQefQr/zGzNSxG6NEJ1iUppQO0LbmXHv58i48uXiWl0AtGpTV3HEpHDpImViOyVteA98jYsxVMjidoX3YnxaBHQ6lDzxD9R86RzsUUF7PjoSUoK81xHEpHDpGIlIgDk/7aKzDmlV6rV6T4EX3wtx4kiS61zb8RXqzGFOzaSMfMV13FE5DDpUKCIUFKYV3pelS0h4YxexB2T5jpS0KnMavHrn+hR6df3RMeSesk9bJ44lJyln1Lj2DO1H0RCkCZWIkLmnDcoyvidqDpNSDn7atdxIlZ03eYkd7oCgJ2fPkexlmAQCTkqViIRLu/XH8he9BEYD7V7DMX4olxHimiJZ/QmumFLinPSyfjyJddxRKSSVKxEIlhJwR52TnsGsCS1v4yY+se6jhTxjMdLne5DML5odv8wk9y1C1xHEpFKULESiWAZs16nKHMLUXWbk3RWX9dxpExU7cYklx2S3fnZGIr3ZDlOJCL+UrESiVB71i8l57tPwOOjTo87MF4dAgwmCWkXE9O4FSW7M0mfMc51HBHxk4qVSAQqyc9l5/RnAUjq0I/ous0dJ5IDGeOhdvc7MFEx5K6cw+5Vc11HEhE/qFiJRKCMWa9RnLWd6PotSDrzUtdx5CCiUhqQcs61AKR//qIOCYqEAK1jJRJh8jatIGfpp+DxlU5EtLp6wFVmzauKxLe5kN2r55G/cTkZX70G9A/Ya4tI4GliJRJBbHEh6Z8+D0BSu78QnXq040RSEWM81L7gFvBGsfv7L/jqq69cRxKRQ1CxEokgWQvep3DnRnwpDXUVYAiJqtWIpPaXATBo0CDy8nQvQZFgVWGxMsYcZYz5yhiz0hjzozHm9nK2McaY54wx64wxy40xp1VNXBE5XGvXriXzm8kA1LrgFowv2nEiqYykM/sQVfso1q5dy2OPPeY6jogchD8TqyLgTmvtCcCZwC3GmBMP2OZCoEXZ4wZgbEBTisgRsdZy4403QnEhNVt3Je7oU1xHkkoy3ihqdbsVgCeeeIIVK1Y4TiQi5amwWFlrN1trl5R9nA2sBBodsNklwERb6lsg2RjTIOBpReSwvPHGG8ycORNPXCIpna91HUcOU2zjVlx//fUUFhYyaNAgSkpKXEcSkQNU6hwrY0xToA1w4D0WGgG/7vP5Jv5YvkTEgR07djB06FAAUroMxFsjyXEiORIjR46kXr16zJ07l1dffdV1HBE5gN/LLRhj4oH3gDustQcupmLK+SO2nNe4gdJDhTRp0qQSMUXkcN1zzz3s3LmTLl26sK5VF9dx5Ai1GfkNxW2vgn8/xY2Dh/Dw8hp4ayaXu+36J3pUczoR8WtiZYyJorRUvWWtfb+cTTYBR+3zeWPg9wM3sta+ZK1Ns9ampaamHk5eEamE+fPnM378eKKjoxk7dizGlPc7kISaGiecTWyz0ynJ303G7Amu44jIPvy5KtAArwIrrbWjD7LZR8BVZVcHngnsstZuDmBOEamk4uJibr219GTnu+++m+OOO85xIgkUYwy1zr0BPD52fz+D/N9Xu44kImX8mVh1AK4EuhhjlpY9uhtjbjTG3Fi2zTTgZ2Ad8DJwc9XEFRF/vfLKKyxZsoSjjjqK++67z3UcCbCoWo1IbNsLgPQZY7ElxY4TiQj4cY6VtXYu5Z9Dte82FrglUKFE5Mjs3LmT+++/H4DRo0dTs2ZNx4mkKiS178vuH76iYMs6cpbPIOHUbq4jiUQ83StQJMSVd1+6nZ+9QE56OrFHn8KdC2O5a1Hg7l0nwcMTHUdKl4Hs+OhJMudMpEbLDnjjElzHEolouqWNSJjJ37Ku7CbLXlLOHaQT1sNcjeM7EdPkJEr2ZLFr7puu44hEPBUrkTBibQkZM8YBloTTLya6jpY1CXelJ7IPAuMh+7vpFGz92XUkkYimYiUSRnb/8BX5v6/CUzOZ5A6Xu44j1SQ6tSkJp18MtoT0GeMoPe1VRFxQsRIJEyX5uWTMHg9AyjnX4omp4TiRVKfkjpfjqZFM/m8ryF0523UckYilYiUSJnYt+BcluzOJbtiSmq06u44j1cwTU5Pks68CIGPWBEoK8x0nEolMKlYiYaBo1zay/vMBALW6XK8T1iNU/EldiarbnOLs7WQvnOo6jkhEUrESCQMZsydAcSE1TjibmEbHu44jjhiPl5QuAwHY9e27bN6sG2CIVDcVK5EQl//bKnJXzsb4okn50wDXccSxuKNPIe7YdtjCPP7v//7PdRyRiKNiJRLCrLWkz3wZgIQzeuFLqus4kQSDlM7XgsfLa6+9xtKlS13HEYkoKlYiIWzKlCkU/L4aT81kktr1cR1HgkRUrUYknHYR1lqGDh2q5RdEqpGKlUiI2rNnD/feey8AyZ2u1PIKsp+ks/qRkpLCV199xb///W/XcUQihoqVSIh65pln2LhxI1GpTYk/6VzXcSTIeOMSePDBBwG46667KCgocBtIJEKoWImEoG3btvHYY48BkNLlOozH6ziRBKObbrqJli1bsnbtWsaNG+c6jkhEULESCUGPPPIIOTk5dO/enbimp7qOI0EqKiqKJ598Eij9fyYrK8txIpHwp2IlEmLWrVvHuHHj8Hg8jBw50nUcCXIXX3wxHTt2ZMeOHXtLlohUHRUrkRAzfPhwioqKuPrqq2ndurXrOBLkjDF7C9Xo0aP5/fffHScSCW8qViIh5D//+Q/vvPMOsbGxPPzww67jSIho3749f/7zn9mzZ8/eE9pFpGqoWImECGst99xzDwC33347jRs3dpxIQsljjz2G1+vl1VdfZeXKla7jiIQtFSuREDFt2jRmz55NrVq1GDZsmOs4EmJatmzJ9ddfT0lJCffdd5/rOCJhS8VKJAQUFxfvLVPDhw8nOTnZcSIJRQ888AA1a9bkww8/ZN68ea7jiIQlFSuREDBx4kR++OEHjj76aG655RbXcSRE1a9fnzvvvBOAu+++W7e6EakCKlYiQW7Pnj2MGDECgEcffZSYmBjHiSSU3XXXXdStW5f58+czdepU13FEwo7PdQCRUNV02Cd+b7v+iR6H/T4vvvgimzZt4pRTTuHyyy8/7NcRAUhISGDEiBHceuutDB8+nJ49e+L1auV+kUDRxEokiGVlZfH4448DpVd1eTz6lpUjd/3119OsWTNWrlzJm2++6TqOSFjRT2mRIDZ69Gh27txJx44dufDCC13HkTARHR3NQw89BJSe0J6fn+84kUj4ULESCVLbt29n1KhRQOm0yhjjOJGEk8svv5xWrVqxYcMGXn75ZddxRMKGipVIkHr88cfJycnhwgsvpFOnTq7jSJjxer08+uijQOlFEbt373acSCQ8qFiJBKFff/2VF198EYC///3vjtNIuLrkkkto27YtW7du5bnnnnMdRyQsqFiJBKGHH36Y/Px8+vbtS5s2bVzHkTBljOGxxx4D4MknnyQjI8NxIpHQp2IlEmTWrFnD+PHj8Xq9utGyVLmuXbvStWtXMjMzeeqpp1zHEQl5KlYiQWbEiBEUFxdzzTXXcNxxx7mOIxHgv4ebn332WbZs2eI4jUhoU7ESCSJLly5lypQpxMTE8MADD7iOIxGiXbt29OrVi9zc3L0ntIvI4TGu7hWVlpZmFy1a5OS9RQKhMiuv+2vbew+zZ91/SEi7hFpdrw/460tkqcyK/z/++CMnnXQSUVFRrF27liZNmlRhMpHQY4xZbK1Nq2g7TaxEgkT+76vZs+4/mKgYks681HUciTCtWrWiX79+FBQUaGolcgRUrESCRObXpbcWSTj9Yrw1kx2nkUj04IMP4vF4eO211/jpp59cxxEJSSpWIkEg79cfyFv/HSY6jsS2f3YdRyLUcccdx1VXXUVxcbGuSBU5TCpWIo5Za/dOqxLTeuGNS3ScSCLZiBEj8Pl8vPnmm6xatcp1HJGQo2Il4ljehmXk//oDnth4Es+4xHUciXDNmjVj4MCBlJSU8OCDD7qOIxJyVKxEHCqdVr0BQGLbP+OJjXecSAT+9re/ERMTw5QpU1i+fLnrOCIhRcVKxKG8nxdR8PtqPHGJJJx+ses4IgA0btyYQYMGAWg9NZFK8rkOIBKp9j23KunMPnii4xwnkkh24LpsxcVnYHwxTJ06lQZXP0NMgxb7PV+ZNbJEIokmViKO7Fk7n4KtP+GtmUJ8m+6u44jsxxufQsJppeUpc+6bjtOIhA4VKxEHrC0hc+7bACSeeSmeqFjHiUT+KLHdXzDRceT9vJj833SFoIg/VKxEHMhdM5/C7evxxtcm4dRuruOIlMtbI4mE0y4CIHPe247TiIQGFSuRamZtCbvKplVJ7S/F+KIdJxI5uMS2vUunVr8sIW/TStdxRIJehcXKGPOaMWabMeaHgzx/jjFmlzFmadljROBjioSP3FXzKNyxAW9CHeJPvsB1HJFD8sYlknB6TwB2zX3LcRqR4OfPxOp1oKJjFV9ba08te+g+CCIHYUuK2TVvEgBJ7S/D+KIcJxKpWOIZvTDRNcjbsJS8TT+6jiMS1CosVtbaOUB6NWQRCXu5q+ZSuHMj3sRU4k8+z3UcEb944xJITCu9K4CmViKHFqhzrNobY5YZY6YbY1oF6DVFwootKSZz77SqL8araZWEjsQzLsHE1CRvw3Lyfi33zBARITDFaglwtLX2FGAMMPVgGxpjbjDGLDLGLNq+fXsA3lokdOSu+pqi9E14E+sSf1JX13FEKsUTG793apWpqZXIQR1xsbLWZllrc8o+ngZEGWPqHGTbl6y1adbatNTU1CN9a5GQUTqtmgxA0lmaVkloSjzjEjwxNcnf+D2zZs1yHUckKB1xsTLG1DfGmLKP25a95s4jfV2RcLJ75RyK0jfhS6pHfGtNqyQ0eWJqktC2NwAjRozAWus4kUjw8We5hUnAfKClMWaTMWagMeZGY8yNZZv0AX4wxiwDngP6WX23iexlS4rZ9c2+0yrdolNCV+LpPfHExvP1119raiVSjgp/wltr+1fw/PPA8wFLJBJmSqdVv+FLrk/NVl1cxxE5Ip6YGiSe0ZvMr9/gwQcf5JxzzqHsoIWIoJXXRarUftOq9ppWSXhIOP1iUlJSmDNnjqZWIgdQsRKpQvtPqzq7jiMSEJ6YGgwdOhSAhx56yHEakeCiYiVSRTStknA2ePBgUlJSmD17tqZWIvtQsRKpIppWSThLSkraO7V68MEH3YYRCSIqViJVQNMqiQSaWon8kYqVSBXQtEoigaZWIn+kX6NFAkzTKgkWTYd9UuXvMXjwYEaPHr13anXOOedU+XuKBDNNrEQCTNMqiSSaWonsT8VKJIBKp1VTAEhqf5mmVRIRBg8eTHJyss61EkHFSiSgcld9vfeegFplXSLFvlMrrWslkU7FSiRAbEkxmfNKz61K1LlVEmFuu+02kpOTmTVrFrNnz3YdR8QZFSuRAMldNZei9E14k+oR31rTKoksSUlJDBkyBNDUSiKbipVIAOx/JaDOrZLIdNttt5GUlMRXX33F119/7TqOiBMqViIBkLt6HoU7f8WbWFfTKolYycnJ3HHHHYCmVhK59Gu1SJnDXfPH2hJ2zdt3WhUVyFgiIeWOO+7gmWee4csvv2Tu3Ll07NjRdSSRaqWJlcgRyl01j8KdG/EmphJ/UlfXcUScSk5O5vbbbwc0tZLIpGIlcgSsLWHXN5OA/66yrmmVyB133EFiYiJffPEF8+bNcx1HpFrpUKDIEchd/Q2FOzStkshT0aFz07o7fDOZc68YzJ5fllRTKhH3NLESOUyl51aVTavOvFTTKpF9JKRdgomOI2/9d3zzzTeu44hUGxUrkcNUOq3agDchlfiTznMdRySoeOMSSDy9J6BzrSSyqFiJHIb9plXtL8X4NK0SOVDCGb0w0XF8/vnnmlpJxFCxEjkMmlaJVExTK4lEOnldpJJKrwT877pVmlaJHErCGb3IWvwRn3/+OQ2u/AcxjU445Pbrn+hRTclEqoYmViKVlLtmPoXb1+NNqKNplUgFvHEJJJRNrTLLDp+LhDMVK5FK2P/cqss0rRLxQ+IZZVcI/rKE/N9Wuo4jUqVUrEQqYc+abzWtEqkkb1wiCadfDGhqJeFPxUrET9aWkDnvbaBs3SpNq0T8llh2hWDp1GqV6zgiVUbFSsRPe9aWTaviaxN/8vmu44iEFE2tJFKoWIn4oXRapXOrRI7E/6ZWizW1krClYiXihz1rvqVw2y9l0yqdWyVyOPafWr3tOI1I1VCxEqnAfudWtb8M44t2nEgkdO1/rpWuEJTwo2IlUoHc1d+UXQmYqnOrRI5Q6dSqbF2ruZpaSfhRsRI5BN0TUCTw9k6t1n9H3qYVruOIBJSKlcgh5K6ap3sCigTYvvcQ3DX3LcdpRAJLxUrkIGxJ8f+mVWfpSkCRQEpo2xsTXYO8DcvI+/UH13FEAkbFSuQgclfPo3DnRryJdYk/6VzXcUTCijc2nsS0SwDYpSsEJYyoWImUY79pVfvLMF5Nq0QCLfGMSzAxNcnbsFxTKwkbKlYi5chdNZfCnb/iTapH/EldXccRCUuefaZWmTrXSsKEipXIAWxJ8f6rrGtaJVJlEtN6YmJqkr/xe/I2LHcdR+SIqViJHGD3itkUpW/Cl1SP+NaaVolUJU9sPIln9AIgc+6bWGsdJxI5MipWIvuwxUX/O7eqQ3+M1+c4kUj4S0y7BE9sAvmbVjBjxgzXcUSOiIqVyD52/ziToszN+Go1omarzq7jiEQET0wNEtv9GYARI0ZoaiUhTcVKpIwtLiRz3mQAkjv0x3i8jhOJRI6E0y7CUyOJBQsWMG3aNNdxRA6bipVImZzlMyjO2kZU7SbUOL6T6zgiEcUTHUdSuz6AplYS2lSsRIC8vDx2fTMFgKSOl2taJeJAfJvuNGjQgCVLljB16lTXcUQOi4qVCPCClyRkAAAfG0lEQVTSSy9RnLOTqNSm1Gh5lus4IhHJExXD/fffD5ROrUpKShwnEqk8FSuJeLm5uTz22GMAJHe6AmP0bSHiyvXXX0/jxo354YcfePfdd13HEam0Cq8lN8a8BlwEbLPWti7neQM8C3QHcoEB1tolgQ4qUlXGjh3L1q1bia5/LHHHtnMdRySitXzgC3JbXQKbXuCKm+/insVxBz00v/6JHtWcTqRi/vxq/jrQ7RDPXwi0KHvcAIw98lgi1SM7O5uRI0cCkNzxr5T+niAiLsWfdC7epHoUpW9i94rZruOIVEqFxcpaOwdIP8QmlwATbalvgWRjTINABRSpSs8++yzbt2+nffv2xDZPcx1HRADjjSK5Q38Ads17G1tc5DiRiP8CcTJJI+DXfT7fVPY1kaCWnp7OP/7xDwD+/ve/a1olEkRqtuqMr1ZjijK3kPO9VmOX0BGIYlXe30blLkBijLnBGLPIGLNo+/btAXhrkcP31FNPsWvXLs4991w6d9Yq6yLBxHi8JHe6AoBd8yZRUpjvOJGIfwJxI7RNwFH7fN4Y+L28Da21LwEvAaSlpWn1N6lyTYd9Uu7Xi3My+O2fTwPwQ4MLD7qdiLhTo+VZRNc7hoKtP5Hz3TQS2/Z2HUmkQoGYWH0EXGVKnQnsstZuDsDrilSZXfOnYIvyiWtxJjENW7qOIyLlMMbzv6nVt+9Skp/rOJFIxSosVsaYScB8oKUxZpMxZqAx5kZjzI1lm0wDfgbWAS8DN1dZWpEAKNq1jeylnwJm7w9tEQlOsc3TiGl0IiV7ssha9KHrOCIVqvBQoLW2fwXPW+CWgCUSqWKZ896GkiJqnngO0alNXccRkUMwxpD8p6vY+vYwsv7zAQmn9cAbl+g6lshBaYlpiSiFO39l9w8zweMlqePlruOIiB9ij2pNbLPTsAW5ZC14z3UckUNSsZKIkjn3bbAlxJ90HlEpDV3HERE/JXe6EoDsxR9TlHOopRVF3FKxkohRsPUncld9Dd4oks7q5zqOiFRCTIMW1DjuLGxRPru+meI6jshBqVhJxMiY9ToACaf1wJdYx20YEam0pE5XgPGQs+xTCjPKXdVHxDkVK4kIe9YvJW/9d5joGiS1v8x1HBE5DNF1mlCzdVcoKSZzzhuu44iUS8VKwp61JWTOfh2ApDP76IoikRCW3PFy8EaRu+prFi9e7DqOyB+oWEnYy101j4It6/DG1yLh9J6u44jIEfAlppJ4+sUADBs2zHEakT9SsZKwZouLyPx6IgBJHS7HEx3rOJGIHKnEMy/FE1OTL774ghkzdINmCS4qVhLWcpZ/TlHGZny1GhF/8nmu44hIAHjjEkg8sw9QOrUqKSlxnEjkf1SsJGzl5OSUrrIOJJ99FcbjdZxIRAIl4fSLadiwIUuWLOHdd991HUdkLxUrCVvPPPMMJbsziW5wHDWOO8t1HBEJIE9ULA8++CAAw4cPp6CgwG0gkTIqVhKWduzYwZNPPglAyjkDMMY4TiQigXbNNdfQsmVLfvrpJ15++WXXcUQAFSsJUw8//DDZ2dnENjud2CYnu44jIlXA5/Px2GOPAfDQQw+RlZXlOJGIipWEoTVr1jB27Fg8Hg8pna9xHUdEqlDv3r0566yz2L59OyNHjnQdR0TFSsLPsGHDKCoq4pprriE6tanrOCJShYwxjBo1CoDRo0fz66+/Ok4kkU7FSsLKnDlz+OCDD6hRowYPP/yw6zgiUg3OPPNM+vbtS15eHsOHD3cdRyKcipWEjZKSEu666y4A7rnnHho2bOg4kYhUl8cff5zo6GjeeOMNlixZ4jqORDAVKwkbU6ZMYeHChTRo0GBvwRKRyNCsWTMGDx4MwJ133om11nEiiVQqVhIW8vLyuO+++wB45JFHqFmzpuNEIlLdhg8fTq1atZg1axaffPKJ6zgSoXyuA4hUVtNhf/yBuWvBv8jcsIGo1KY8vDqVR8rZRkTCW0pKCiNGjOCOO+7g7rvv5oILLiAqKsp1LIkwmlhJyCvO3cWub94BIKXztbp1jUgEu+mmmzj22GNZtWoVr7zyius4EoFUrCTkZc59G1uQS2yz04hrdprrOCLiUHR09N71rEaMGEFmZqbjRBJpVKwkpBVsX0/O0ulgPKR0vtZ1HBEJAr179+bss89mx44dWnZFqp2KlYQsay3pX7wEtoSENt21GKiIAKWLhj777LMYYxgzZgyrVq1yHUkiiIqVhKw9a+aTv3E5ntgEkjr+1XUcEQkip556Ktdffz1FRUUMGTJEyy9ItVGxkpBkiwpI/+pVAJI7XYE3LsFxIhEJNo8++ihJSUl8+umnTJs2zXUciRAqVhKSshZOpXjXVqJSmxJ/ajfXcUQkCKWmpvLggw8CMGTIEAoKCtwGkoigYiUhpyh7B7vmly2v0PV6La8gIgd1yy23cPzxx7N27VrGjBnjOo5EABUrCTmZsydgC/OocdxZxB19ius4IhLEoqKiePrppwF4+OGH2bp1q+NEEu5UrCSkfPvtt+z+8SvwRpGs5RVExA/dunXjoosuIisri/vvv991HAlzKlYSMoqLi7nlllsASGzbm6jk+o4TiUioGD16NFFRUYwfP54FCxa4jiNhTMVKQsa4ceNYsmQJ3sRUks68zHUcEQkhLVq04M4778Ray0033URxcbHrSBKmVKwkJGzZsoXhw4cDUKvrDXiiYx0nEpFQ87e//Y0mTZrw3XffMXbsWNdxJEypWElIuPvuu9m1axfdu3cnrsWZruOISAiqWbMmzz33HADDhw9ny5YtjhNJOFKxkqA3a9Ys3nzzTWJjYxkzZgzGGNeRRCRE9ezZkx49epCVlcVdd93lOo6EIRUrCWoFBQV7T1i/7777aN68ueNEIhLKjDE899xzxMbG8tZbbzFr1izXkSTMqFhJUHvmmWdYsWIFxx57LPfcc4/rOCISBpo3b7532YWbb75ZK7JLQKlYSdDauHEjDz30EADPP/88sbE6YV1EAuPuu++mRYsWrFy5cu8CoiKBoGIlQclay2233UZubi59+vThggsucB1JRMJIbGwszz//PFC6Ivv69evdBpKwoWIlQem9997jww8/JCEhQb9NikiVOP/88+nbty+5ubkMGjQIa63rSBIGfK4DiBwoIyODW2+9FYAnnniCxo0bO04kIuGg6bBP/vC14tSL8MR+wueff07qxXcR37oLAOuf6FHd8SRMaGIlQeeuu+5i69atdOjQgRtvvNF1HBEJY96aKaR0uQ6AjC9fpnh3puNEEuo0sZKg8N/fJPesX8q2Ka+B18cvJ1xB8/unO04mIuGuZuuu7P5xFnkblpL+5cuk9rzbdSQJYZpYSdAoKcwj/bPSk0mTzupHVO2jHCcSkUhgjKFWt1sxUTHkrpxN7k8LXUeSEKZiJUFj19y3KcrcQlSdo0lq9xfXcUQkgkQl1ye54xUApH/2AllZWY4TSahSsZKgkL9lHVkLpwKG2hfehvFGuY4kIhEmIa0n0fVbUJy9Y+8CoiKV5VexMsZ0M8asNsasM8YMK+f5AcaY7caYpWWP6wIfVcJVfn4+O6c9A7aEhLSexDRs6TqSiEQg4/FS+8LB4PHywgsvMGfOHNeRJARVWKyMMV7gBeBC4ESgvzHmxHI2nWKtPbXs8UqAc0oYe+ihhyjcvh5fcn2SO13hOo6IRLDous1JatcHgKuvvprs7GzHiSTU+DOxaguss9b+bK0tACYDl1RtLIkU8+fPZ+TIkYChdo8heKLjXEcSkQiX1KEfbdq0Yf369dx5552u40iI8adYNQJ+3efzTWVfO9BfjDHLjTH/Msboci6p0O7du7n66qspKSkhsW1vYhu3ch1JRATjjWLixIlER0fz8ssvM23aNNeRJIT4U6xMOV87cN3/fwNNrbUnA18AE8p9IWNuMMYsMsYs2r59e+WSStgZNmwYa9eupVWrVjoEKCJBpXXr1jz66KMAXHfddaSnpztOJKHCn2K1Cdh3AtUY+H3fDay1O621+WWfvgycXt4LWWtfstamWWvTUlNTDyevhIkvv/yS559/Hp/PxxtvvIHxRbuOJCKyn6FDh9KxY0c2b97MLbfc4jqOhAh/Vl5fCLQwxjQDfgP6AZfvu4ExpoG1dnPZpz2BlQFNKWElMzOTAQMGAPDAAw/Qpk0bmPL7of+QiMgByrv3XyB5vV5ef/11TjnlFCZPnkyvXr3o27dvlb6nhL4KJ1bW2iLgVuAzSgvTO9baH40xDxtjepZtdpsx5kdjzDLgNmBAVQWW0Hf77bezadMm2rZty7Bhf1i9Q0QkaBxzzDGMGjUKgJtvvpnff9cvgXJofq1jZa2dZq09zlp7jLX272VfG2Gt/ajs4/usta2stadYaztba1dVZWgJXW+//TYTJ04kNjaWCRMm4PPpdpUiEtxuuOEGunXrRnp6OldeeSXFxcWuI0kQ08rrUm3WrVvHoEGDAHj22Wc5/vjjHScSEamYMYbx48dTt25dZs6cyRNPPOE6kgQxFSupFvn5+fTt25ecnBwuvfRSrr/+eteRRET8Vr9+fd544w2g9NzQuXPnOk4kwUrFSqrFfffdx5IlS2jatCkvvfQSxpS3ioeISPA6//zzuffeeykuLubyyy/XEgxSLhUrqXIff/wxTz/9ND6fj8mTJ5OcnOw6kojIYXnkkUc488wz+fXXXxk4cCDWHriso0Q6FSupUps2bdq7tMJjjz1Gu3bt3AYSETkCUVFRTJo0iaSkJKZOncoLL7zgOpIEGRUrqTJFRUX89a9/ZefOnVxwwQW655aIhIWmTZvyyiuvAHDnnXeyZMkSx4kkmKhYSZW59957mTNnDvXr12fChAl4PPrfTUTCQ58+fbjxxhspKCjgz3/+Mzt27HAdSYKE/qaTKjFp0iRGjx6Nz+fjnXfeoV69eq4jiYgE1NNPP80ZZ5zBhg0b6NevH0VFRa4jSRDQ6owScMuWLWPgwIEAPPPMM3Tq1MlxIhGRyvH3djlFaTdTd8O9fPnll9x///08+eSTVZxMgp0mVhJQO3fupHfv3uzZs4cBAwZw8803u44kIlJlfImpvPvuu/h8Pp566immTJniOpI4pmIlAVNcXEz//v355ZdfSEtLY+zYsVqvSkTC3tlnn733foLXXnsty5cvd5xIXFKxkoAZPnw4M2bMIDU1lffff5/Y2FjXkUREqsXgwYO58soryc3NpXfv3lo8NIKpWElATJw4kZEjR+L1ennnnXc46qijXEcSEak2xhj++c9/0qZNG37++WcuvfRSCgoKXMcSB4yrVWPT0tLsokWLnLy3HJ6Dncy5Z8Mytr0zAkqKqXXejSScdlE1JxMRCQ5Fu7ax+Y2hlOzOpGbrLtTuPuQPp0Ssf6KHo3RyJIwxi621aRVtp4mVHJGC7RvY/sFjUFJMwhm9VKpEJKL5kupS9y8PYKJi2P3DTHbNe9t1JKlmKlZy2Ipy0tn2rwex+bupcdxZpHS+1nUkERHnYhq0oE7Pe8F42DVvEjnff+E6klQjFSs5LCUFe9j+r4coztpOTMPjqX3RnRij/51ERABqHNuWWucOAmDnp2PYs36p40RSXfQ3oVSaLSlmx4cjKdj6E77kBqT+5f/wRMW4jiUiElQSTutBYts/Q0kx2z94jILt611HkmqgYiWVYkuK2fnJ0+z5eRGeuETqXvog3hpJrmOJiASl5HMGUOP4TtiCXLZN+T8K039zHUmqmIqV+M1aS/pnL7B7xSxMVCx1/zKCqFqNXMcSEQlaxnio02MIMU1Opnh3Blsn/40NGza4jiVVSMVK/GKtJePLl8hZ/jnGF03dPiOIaXS861giIkHP+KKp+5f/I6bRiRRnb6dLly789psmV+FKxUoqZK3lvvvuI3vxv8HrI7X3cGKbnOw6lohIyPBEx1H30geIrt+Cn3/+ma5du7J161bXsaQKqFhJhR599FFGjhwJHi+pl9xHXPPTXUcSEQk5npia1L3sYU4++WRWr17Neeedp1vfhCEVKzkoay2PPPIII0aMwOPxUOeiO6nRop3rWCIiIcsbl8CMGTM4/vjj+f777zn33HPZtm2b61gSQCpWUq6SkhKGDBmyt1SNHz+emiec7TqWiEjIq1u3Ll988QUtWrTgu+++o1OnTmzcuNF1LAkQFSv5g6KiIq699lqeffZZoqKimDJlCldddZXrWCIiYaNRo0Z8/fXXnHLKKaxZs4YOHTqwatUq17EkAFSsZD95eXn06dOHCRMmULNmTT755BP69OnjOpaISNipV68es2bNokOHDmzatIlOnTqxePFi17HkCKlYyV7Z2dl0796dDz/8kJSUFL744gvOO+8817FERMJWcnIyn3/+Od26dWPHjh107tyZWbNmuY4lR8BYa528cVpaml20aJGT95b/aTrsEwAKM7ew/b1HKNyxAW/NFOr2fYTo1KZuw4mIhKH1T/T4w9cKCgq46qqrmDJlCtHR0fzzn/9kwIAB1R9ODsoYs9ham1bRdppYCXkbv2fLxKEU7tiAr1Zj6v31SZUqEZFqFB0dzVtvvcXgwYMpKCjgmmuuYejQoRQVFbmOJpWkiVWEq93tVtJnjIOSYmKbn05qz3vwxNR0HUtEJGJlL/2U9BljoaSYCy64gMmTJ5OcnOw6VsTTxEoOqbCwkFtvvZX0z16AkmIS2/6Zun8ZoVIlIuJYwqndqNfv73jiEvnss89o164dq1evdh1L/KRiFYHWr19P586deeGFF8Dro3b3IaR0vhbj8bqOJiIiQOxRrWlw9dOcfPLJrFmzhrZt2zJ58mTXscQPOhQYYSZPnsygQYPIysqiYcOGlHQZQkyjE1zHEhGRcpQU7GHntGfIXT0PgJqtu1Dr3BvxxNT4w7blnRQvgaNDgbKf7Oxsrr76avr3709WVha9evVi+fLlKlUiIkHMEx1HnUuGUeuCWzG+GHb/MJPNr99O/u86NBisVKwiwPz582nTpg0TJ04kLi6OcePG8f7771O7dm3X0UREpALGGBJO7UaDq58hqm5zijI3s+XNu9n1zRRssa4aDDYqVmEsIyODm266iQ4dOvDTTz9xyimnsHjxYgYNGoQxxnU8ERGphKg6R9HgylEknNELbAmZX7/B5gl3kP/bStfRZB8qVmHIWstbb73F8ccfz7hx4/B6vQwbNowFCxZwwgk69CciEqqML4paXa6j7mWP4EuqR+H29Wx58252fjqG9PR01/EEFauws2rVKs4991yuuOIKtm3bRqdOnVi6dCmPP/44MTExruOJiEgAxDVrQ4OBL5DYvi94fOQs+4yWLVsyYcIESkpKXMeLaCpWYWLjxo1cd911tG7dmpkzZ1K7dm3Gjx/P7NmzadWqlet4IiISYJ6oWFLOvpKG14whpslJ7NixgwEDBnDaaafx8ccf4+qq/0inYhXitmzZwm233UaLFi149dVXARg0aBCrV69mwIABOpdKRCTMRdU5inr9HmPixIk0atSIZcuWcfHFF3PWWWcxc+ZM1/EijopViNq4cSN33303zZs3Z8yYMRQWFnL55ZezcuVKxo0bpyv+REQiiDGGK6+8knXr1vH000+TmprKt99+S9euXenatSuff/65JljVRMUqhFhrmTNnDn369KFZs2b84x//YM+ePfTq1Ytly5bx1ltv0aJFC9cxRUTEkdjYWO644w5+/vlnHn30UZKSkpg5cyYXXHABJ554Ii+++CI5OTmuY4Y1rbweAprcMYXdq+eRveRjCrf9UvpFj5caLTuSeEYvYhr8r0xVduXdpsM+CWRUEREJIsV5OeQs+YTs76ZRnLMTABNTkztuvoGBAwfqHNxK8HfldRWrIJWbm8vHH3/MpEmTmPrvj6FsEThPjWQSTr2Q+FO74Uv44+E+FSsRETmQLS4id818shf/m/zfVuz9+kknnUT//v3p168fzZo1c5gw+PlbrHzVEUb8s2XLFr744gumT5/Ohx9+yO7du8ueMcQ0OZn4k7pS8/izMb4opzlFRCS0GK+Pmid0ouYJncjfso6Lo3/kX//6F99//z3ff/89999/P+3ataNXr16cf/75nHrqqXg8OlvocGhi5dCuXbtYsGABM2bM4PPPP2f58uX7Pd+uXTv69+/PU2tSyp1OiYiIHC5bXMieX75j98rZ7Fm7AFuYt/c5T1wisU1PJa5pG2KOao0vub7fV5mH682gAzqxMsZ0A54FvMAr1tonDng+BpgInA7sBPpaa9dXNnQ4y8zMZMWKFSxcuHDvY82aNfttExcXx5/+9CfOO+88evXqRfPmzQF4WofrREQkwIw3ihrHtqXGsW0pKchjz8+LyPtlCXvWf0dx1nZyV84hd+UcADyxCUTXP5boBscR06AF0XWb402sgzGaah2owomVMcYLrAHOAzYBC4H+1toV+2xzM3CytfZGY0w/oLe1tu+hXjfcJlZ79uzh2DveojhnJ8XZOynK3k7hzt8oyviNwvTfKcnN/OMf8vqIrtuc2CYnE9usDbGNTsD4oqs/vIiISBlrLUUZv7PnlyXkrV9K/u+ry/07zPii8aU0JCqlIb7ajfEl1ceXUJsv/q83DRs2pFatWmG1lmLATl43xrQHHrTWXlD2+X0A1trH99nms7Jt5htjfMAWINUe4sWrulht3bqVX3/9de+6Hfv+01pLSUnJfo+ioiIKCwv3e+Tn55Obm7vfIzs7m4yMDNLT08nIyCAjI4MdO3ZUeI8m44vBV6sR0fWOKW37DY4jOvVojFfnS4mISPCy1lKcvZ38zWsp2LyWgi1rKNzxK8W7Mw7552JjY6lXrx61atUiJSVl7yM5OZkaNWrs94iLiyMqKuoPD4/H84fHf8vavv+MjY2ldevWVfrfIZCHAhsBv+7z+Sag3cG2sdYWGWN2AbWBHf7FDbxJkyYxZMiQans/n8+HrZGCL7423vhaeBPqEFWrEb5ajYiq1QhvQm2NTEVEJOQYY/Al1sWXWJeaLTvs/XpJ/m4K038vPTKz8zeKsrZSnL2TZjUK+O2338jKymLDhg1s2LChyjOecMIJrFixouINq4E/xaq8Od6Bkyh/tsEYcwNwQ9mnOcaY1X68v7/q4LDIFRUVQdZ2irO2u4oQ7JzuHzkk7Zvgpv0TvLRvyrHSxXuuXFneYcdA75+j/dnIn2K1CThqn88bA78fZJtNZYcCk4A/HBuz1r4EvORPsMoyxizyZ0Qnbmj/BC/tm+Cm/RO8tG+Cm6v948+xqYVAC2NMM2NMNNAP+OiAbT4Cri77uA8w81DnV4mIiIiEowonVmXnTN0KfEbpcguvWWt/NMY8DCyy1n4EvAq8YYxZR+mkql9VhhYREREJRn6tY2WtnQZMO+BrI/b5OA+4NLDRKq1KDjFKwGj/BC/tm+Cm/RO8tG+Cm5P942zldREREZFwo+v/RURERAIk5IqVMaabMWa1MWadMWZYOc/HGGOmlD2/wBjTtPpTRi4/9s9QY8wKY8xyY8yXxhi/Ll+VI1fRvtlnuz7GGGuM0dVO1cif/WOMuazs++dHY8zb1Z0xUvnxc62JMeYrY8x3ZT/burvIGYmMMa8ZY7YZY344yPPGGPNc2b5bbow5raozhVSxKru9zgvAhcCJQH9jzIkHbDYQyLDWHgs8DYys3pSRy8/98x2QZq09GfgX8GT1poxMfu4bjDEJwG3AgupNGNn82T/GmBbAfUAHa20r4I5qDxqB/Pze+RvwjrW2DaUXb71YvSkj2utAt0M8fyHQouxxAzC2qgOFVLEC2gLrrLU/W2sLgMnAJQdscwkwoezjfwFdTTjdrCi4Vbh/rLVfWWtzyz79ltJ10aTq+fO9A/AIpWU3r5znpOr4s3+uB16w1mYAWGu3VXPGSOXPvrFAYtnHSfxxrUepItbaOZSzbuY+LgEm2lLfAsnGmAZVmSnUilV5t9dpdLBtrLVFwH9vryNVz5/9s6+BwPQqTST/VeG+Mca0AY6y1n5cncEE8O975zjgOGPMPGPMt8aYQ/2WLoHjz755ELjCGLOJ0ivoB1dPNPFDZf9eOmJ+LbcQRAJ2ex2pEn7/tzfGXAGkAX+q0kTyX4fcN6b0RpZPAwOqK5Dsx5/vHR+lhzPOoXTS+7UxprW1NrOKs0U6f/ZNf+B1a+0oY0x7Std1bG2tLan6eFKBau8EoTaxqsztdTjU7XWkSvizfzDGnAsMB3paa/OrKVukq2jfJACtgVnGmPXAmcBHOoG92vj7s+1Da22htfYXYDWlRUuqlj/7ZiDwDoC1dj4QS+l96sQ9v/5eCqRQK1a6vU5wq3D/lB1u+ielpUrniFSfQ+4ba+0ua20da21Ta21TSs9/62mtXeQmbsTx52fbVKAzgDGmDqWHBn+u1pSRyZ99sxHoCmCMOYHSYrW9WlPKwXwEXFV2deCZwC5r7eaqfMOQOhSo2+sENz/3z1NAPPBu2TUFG621PZ2FjhB+7htxxM/98xlwvjFmBVAM3G2t3ekudWTwc9/cCbxsjBlC6WGmAfqFvnoYYyZReni8Ttk5bg8AUQDW2nGUnvPWHVgH5ALXVHkm7XsRERGRwAi1Q4EiIiIiQUvFSkRERCRAVKxEREREAkTFSkRERCRAVKxEREREAkTFSkRERCRAVKxEREREAkTFSkRERCRA/h8F4B8vgzRZCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import beta\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "q = beta(5, 5)      # Beta(a, b), with a = b = 5\n",
    "obs = q.rvs(2000)   # 2000 observations\n",
    "grid = np.linspace(0.01, 0.99, 100)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "ax.hist(obs, bins=40, normed=True)\n",
    "ax.plot(grid, q.pdf(grid), 'k-', linewidth=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this code we created a so-called `rv_frozen` object, via the call `q = beta(5, 5)`\n",
    "\n",
    "The “frozen” part of the notation implies that `q` represents a particular distribution with a particular set of parameters\n",
    "\n",
    "Once we’ve done so, we can then generate random numbers, evaluate the density, etc., all from this fixed distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.26656768000000003"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q.cdf(0.4)      # Cumulative distribution function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0901888000000013"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q.pdf(0.4)      # Density function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6339134834642708"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q.ppf(0.8)      # Quantile (inverse cdf) function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general syntax for creating these objects is\n",
    "\n",
    "> `identifier = scipy.stats.distribution_name(shape_parameters)`\n",
    "\n",
    "\n",
    "where `distribution_name` is one of the distribution names in [scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html)\n",
    "\n",
    "There are also two keyword arguments, `loc` and `scale`, which following our example above, are called as\n",
    "\n",
    "> `identifier = scipy.stats.distribution_name(shape_parameters, loc=c, scale=d)`\n",
    "\n",
    "\n",
    "These transform the original random variable $ X $ into $ Y = c + d X $\n",
    "\n",
    "The methods `rvs`, `pdf`, `cdf`, etc. are transformed accordingly\n",
    "\n",
    "Before finishing this section, we note that there is an alternative way of calling the methods described above\n",
    "\n",
    "For example, the previous code can be replaced by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/quantecon/anaconda3/lib/python3.7/site-packages/matplotlib/axes/_axes.py:6571: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8VFXex/HPbya9QwhIj4K6dIQooK6g2BAQdbEv6iPoYkWUl8uqq4iLiitFmsqDrviIiorSBLFLFwKIFGkCShMSCKmkTOY8fySwEQMZyGTOlN/79crLJHOZ+V4Tvpy599xzxRiDUkqp4OKwHUAppZT3abkrpVQQ0nJXSqkgpOWulFJBSMtdKaWCkJa7UkoFIS13pZQKQlruSikVhLTclVIqCIXZeuE6deqY1NRUWy+vlFIBadWqVZnGmJSqtrNW7qmpqaSnp9t6eaWUCkgi8osn2+lhGaWUCkJa7kopFYS03JVSKghpuSulVBDScldKqSCk5a6UUkFIy10ppYKQlrtSSgUhLXellApC1q5QVSpQpA79tNLv73yxp4+TKOU5HbkrpVQQ0nJXSqkgpOWulFJBSMtdKaWCkJa7UkoFIS13pZQKQjoVUilOPN1RqUClI3ellApCWu5KKRWEtNyVUioIabkrpVQQqrLcRaSxiHwjIj+JyAYRGVTJNt1EJFtEfij/eLpm4iqllPKEJ7NlXMBjxpjVIhIPrBKRL4wxG4/bbpExppf3IyqllDpVVY7cjTH7jDGryz/PBX4CGtZ0MKWUUqfvlI65i0gqcB7wfSUPdxGRtSIyX0RaeSGbUkqp0+TxRUwiEgfMAB4xxuQc9/BqoKkxJk9ErgFmAmdX8hz3AvcCNGnS5LRDK3UyJ7sgSddgV6HCo3IXkXDKin2aMebj4x+vWPbGmHkiMklE6hhjMo/bbjIwGSAtLc1UK7lSXmRcJRRn/oI7//Cx70lkDBEpqTgiYywmU+r0VFnuIiLAG8BPxpjRJ9jmDGC/McaIyAWUHe456NWkSnmRMYbVq1fz1ltvsWjRIn5dtx7cpZVsKYTVqk9kg3OJbXUpUU3bIQ6nz/Mqdao8GblfBPQD1onID+XfewJoAmCMeQ3oC9wnIi7gCHCLMUZH5srvFBcXM2XKFF5//XV+/PHHCo8IYbUbEZZYFxAA3AWHKc74BVfWXlxZe8nf8A3O+BTi2lxOwvl9rORXylNVlrsxZjFHf9tPvM0EYIK3QinlbcYYjvy8gtatB7N161YAkpOT6devH3379uXWGb/hiIj6458rLaEkcxcFP68gf92XuA7/RvbS98hd8ynJPy0krt1VlY7k9di+sk1XhVRBz5V7kIPzX6Fwx2oygHPOOYfhw4dz/fXXExERAYBjTuUnYcUZTkS9s4iodxaJXW6iaNd6Di9+l6Jd6zn0+SRyV39Kcs/BRJ7R3Id7pFTVtNxVUDuyYzWZc0fhLsjGERXH6Bf/xf333094ePgpP5eIg6gmbal3axuObFlG1jdvUJL5C7+9M4Ta3e8lrn0Pyk5RKWWflrsKSsa4yV78LtlLpwOGqKbtqdP7MQYN+mu1n1tEiDn3QqKbpXHo6ynkrZnHoc8nUfjrOpJ7PIwjIrr6O6BUNWm5q6BjSl1kzhtDwcbvACHx4ttJ7HKT12e5SFgEyVfeT1SjVhxcMIGCTYtwZe+n7o3DvPo6Sp0OXRVSBZUjR46Q8ckICjZ+h0REU/emZ0m66NYanb4Y27Ir9e8YgzOxHsX7trD/3aHs3bu3xl5PKU/oyF0Fjby8PHr16sWRn1fiiIqn7o3DiGxwrk9eOzy5EWfcPpID05+mJPNXmrbsQN1bRhCedMbvttNZNMpXtNxVUCguLuaGG27gu+++wxlXm7o3PUdESlOfZgiLr0O921/kwIfPULxvKwem/5Mz/voSzthax7Y50dIIWvrK2/SwjAp4brebu+66iy+++IKUlBTq3fqCz4v9KGd0AvVuHkFEvWa4Du/jwIfDcBcVWMmiQpuWuwpoxhgeffRR3nvvPeLi4pg/fz7hte2uSO2IjKHujcMIS6pP8f6fyfhkBMZVYjWTCj1a7iqgTZgwgVdeeYXw8HBmzpxJx44dbUcCwBlbi7o3DccRm0ThL2s5uGAiuiKH8iUtdxWwFi1axKOPPgrAW2+9Rffu3S0n+r3wWvWpd+OzSFgk+eu/JO+H+bYjqRCi5a4C0t69e7nxxhtxuVw89thj3HbbbbYjVSqiXjOSezwEwKEvJ1O4+yfLiVSo0HJXAae4uJi+ffuyf/9+unXrxosvvmg70knFtuxGfMdrwe0ic9YLlOZl2Y6kQoCWuwo4TzzxBMuWLaNRo0ZMnz6dsDD/n9Fb69K7iWzUitK8Q2TO/TfGuG1HUkFOy10FlK+++opRo0bhdDr58MMPqVu3ru1IHhFnGCl9huKISaTwlx/JXTnLdiQV5LTcVcDIysrizjvvBODpp5+mc+fOlhOdGmdcLZJ7PAxA1sKpFB/YYTmRCmZa7iogGGO477772LNnD507d+aJJ56wHem0xDTvRFz7q6HUReaclzGuYtuRVJDSclcB4f3332f69OnExsbyzjvvBMRx9hOpdekAwmo1oCTzFw4v/D/bcVSQ0nJXfi8jI4OHHiqbTjh27FiaNWtmOVH1OCKiqNN7CIiDnPRZFO3bYjuSCkJa7srvPfrooxw8eJDu3bvTv39/23G8IrL+OSScfx0YNwfnj6OkRJcnUN6l5a782oIFC3jnnXeIioritddeC6rb2CVefBthifUoydjJqFGjbMdRQUbLXfmt/Px8Bg4cCMCwYcNo3jy4bkLtCI+i9lUPAmX7t3XrVsuJVDDRcld+69lnn2Xnzp20b9/+2BoywSb6zPOIbX0ZRUVFDBw4UBcXU14TuFMOVFDbtGkTY8aMQUSYPHky4eHhXnneE90sw6Zal/Yncu8PfP3118yYMYO+ffvajqSCgI7cld8xxjBo0CBcLhf9+/fn/PPPtx2pRjljEvnXv/4FwGOPPUZBgd7cQ1WflrvyO3PmzOHzzz8nMTGR559/3nYcn7j33ntp164dv/76Ky+99JLtOCoIaLkrv1JYWMjgwYMBGD58OCkpKZYT+YbT6WT8+PEAjBw5kp07d9oNpAKeHnNXfmX06NFs376dVq1acf/999uO4zNHzwXEtOxKwcbvaHXV7aRcX7bEgt48W50OHbkrv3HgwAFeeOEFAMaNGxfQSwycrlrd/gcJj6Rgy1IKd2+wHUcFsCrLXUQai8g3IvKTiGwQkUGVbCMiMk5EtonIjyLSoWbiqmD27LPPkpeXR69evbjssstsx7EiLL4OCeffAEDWN2/q1Eh12jwZGrmAx4wxq0UkHlglIl8YYzZW2KYHcHb5Ryfg1fL/qhB0sumGJzrEsHnzZl5//XUcDofHd1byx2mN3pBwwfXk/jCf4r2bKdiyFOhlO5IKQFWO3I0x+4wxq8s/zwV+Ahoet1kf4G1TZjmQJCL1vZ5WBa0nnniC0tJS7r77blq1amU7jlWOyBiSLi67J+zh796iuFiXBVan7pSOuYtIKnAe8P1xDzUEdlX4ejd//AdAqUotXbqUjz/+mJiYGJ599lnbcfxCXNsrCavdCFfWPiZPnmw7jgpAHpe7iMQBM4BHjDE5xz9cyR/5w8FCEblXRNJFJD0jI+PUkqqgZIzh8ccfB8ou4GnQoIHlRP5BnGHU6nYXUHYuIjc3124gFXA8mo4gIuGUFfs0Y8zHlWyyG2hc4etGwN7jNzLGTAYmA6SlpemZIsWCBQtYsmQJycnJDBkypNJtgvXYelWim3cismELMvf8xLhx43jyySdtR1IBxJPZMgK8AfxkjBl9gs1mA3eUz5rpDGQbY/Z5MacKQsYYnnrqKQCGDh1KQkKC5UT+RURI+nM/AP7973+TlZVlOZEKJJ4clrkI6AdcJiI/lH9cIyIDRWRg+TbzgO3ANuB/gdC5+kSdtpkzZ7Jq1Srq168fUhcsnYqopm257LLLyM7OZvToE42tlPqjKg/LGGMWU/kx9YrbGOABb4VSwa+0tJR//vOfADz55JPExMRYTuS//vWvf3HhhRcyduxYHn744ZBZkkFVj16hqqyYPn06GzZsoGnTpgwYMMB2HL/WpUsXevbsSV5eHiNHjrQdRwUILXflcy6Xi2HDhgHw9NNPExkZaTdQABg+fDgAEydOZN8+PZ2lqqblrnxu+vTpbN26lbPOOos77rjDdpyA0KFDB66//noKCwt5+eWXbcdRAUDLXflUaWnpsRtTPPnkkyG5ONjpOnqO4tVXX+XAgQOW0yh/p3+zlE999NFHbNq0iaZNm9KvX79j3w/Vueyn4rzzzqN3797MmTOHUaNG6fF3dVI6clc+Y4yb5557DihbS8Zb90UNJUdH7xMnTiQzM9NyGuXPtNyVzxRsWcaGDRto3Lgxd955p+04Aen888+nR48e5OfnM2bMGNtxlB/TwzLKJ4wxZC99H4D8P/Xi3Ge+tJwocBx/yKoo+TJgPuPHj2fIkCHUqlXLTjDl13TkrnyicHs6JQd24IyrTVzbK2zHCWiRDVsQ1bQ9ubm5TJw40XYc5ae03JVPZC/7EID4tOuQsAjLaQJfQpcbARg7diz5+fmW0yh/pOWualzhrvUU7dmIIyqO+PZX244TFKKatKVTp04cPHiQKVOm2I6j/JCWu6px2cvLR+0deuGI1DVkvEFE+Mc//gHAyy+/rHdrUn+g5a5qVPH+7RRuX4WERxLfsbftOEGld+/etGzZkt27dzNt2jTbcZSf0XJXNeroqD2u7VU4YxItpwkuDoeDoUOHAvDiiy9SWlpqOZHyJ1ruqsaUZO2jYPMScDhJuOB623GC0i233EJqaipbtmxh1qxZtuMoP6LlrmpMzsqZYNzEtuxGWIKuQV4TwsPDefTRR4GyuzWV3VpBKS13VUNKC7LJX1d2oZKO2mvW3XffTe3atVm+fDlLliyxHUf5CS13VSNyV8/FuIqIPiuNiJRU23GCWmxsLA88UHYjtJdeeslyGuUvdPkB5XXukkJyV5ddMp/Q6QbLaYJXxWUJSvPPRcIimDNnDg0HvMaeKQNP8idVKNCRu/K6/HVf4T6SQ0T9s4ls3MZ2nJDgjE0itnV3AHJWfmI5jfIHWu7Kq0pLS4+VS8IFf0HkpPdWV16UcP51gJC34Wu9FZ/SclfeNXPmTFyHfyMs6QxizuliO05ICa/dsOz/eamLCRMm2I6jLNNyV141evRoAOLT+iAOp+U0oSf+/LKZSa+99pouKBbitNyV1yxfvpylS5fiiIwlrs3ltuOEpKhGLYhocC6HDh1i6tSptuMoi7TcldccHbXHte+BIyLacprQlVA+eh8zZgxut9tyGmWLlrvyih07djBjxgzCwsKI79jLdpyQFnNOF1JTU9m2bRtz5syxHUdZouWuvGLcuHG43W5uvfVWwuLr2I4T0sThZNCgQcB/302p0KPlrqotOzv72A0jjq5zouy6++67SUhIYOHChaSnp9uOoyzQclfVNmXKFPLy8rj00ktp37697TgKaPv8Iji37KKmrncOIXXop3+40bYKblWWu4i8KSIHRGT9CR7vJiLZIvJD+cfT3o+p/JXL5WL8+PEADB482HIaVVF8x94gDgo2LcaVm2k7jvIxT0bubwFV3fhykTGmffnH8OrHUoFi5syZ/PLLLzRv3pyePXvajqMqCEusS8w5F4K79NhaPyp0VFnuxpiFwCEfZFEBaOzYsQAMGjQIh0OP8vmb+LQ+AOT98BnukkLLaZQveetvYxcRWSsi80WklZeeU/m5lStXsmTJEhITE7nrrrtsx1GViGz4JyLqn4O7MJf8Dd/YjqN8yBvlvhpoaoxpB4wHZp5oQxG5V0TSRSQ9IyPDCy+tbDo6ar/nnnuIi4uznEZVRkRIKB+956ycpRc1hZBql7sxJscYk1f++TwgXEQqnehsjJlsjEkzxqSlpOht1wLZnj17+OCDD3A6nTz00EO246iTiDn3IpxxybgO7ebzzz+3HUf5SLXLXUTOkPJ1XUXkgvLnPFjd51X+bdKkSbhcLm644QaaNGliO446CXH+96rhV155xXIa5SueTIV8D1gGnCsiu0Wkv4gMFJGjt3rpC6wXkbXAOOAWo3fpDWpHjhzh9ddfB+CRRx6xnEZ5Iq7dVUhYJJ999hmbNm2yHUf5QJW32TPG3FrF4xMAXTw6hEybNo2DBw+SlpZGly66ZnsgcEYnENvqUvLWfsa4ceOYNGmS7UiqhuncNXVKjDHH3to/8sgjeqelABKfdi0AU6dOJSsry3IaVdO03NUp+frrr1m/fj3169fnxhtvtB1HnYKIOk244oorKCgoOLYWkApeWu7qlBwdtd93331ERERYTqNO1dHVIsePH4/L5bKcRtUkLXflsW3btjF37lwiIyP529/+ZjuOOg09evTg7LPPZteuXXzyySe246gapOWuPDZ+/HiMMdx2223UrVvXdhx1GhwOBw8//DCg0yKDnZa78khOTg7/+c9/gP++tVeB6c477yQhIYElS5awatUq23FUDdFyV1VKHfopqTcMITc3l8jGrenz3m5dGzyAxcfH079/f0BH78FMy11VybhLyV01F4CE8ul0KrA9+OCDiAjvv/8+v/32m+04qgZUeRGTCg0nG4kf2Z6O6/A+nIn1iG7eyYepVE0566yzuPbaa5k1axavvfYaw4YNsx1JeZmO3FWVctNnA5DQoSficFpOo7zl6LmTV199laKiIstplLdpuauTKs7YSeEva5HwKOLaXmk7jvKibt260bZtWw4cOMD06dNtx1FepuWuTuroqD22dXccUbpmezARkd9Ni9T1/oKLlrs6odKCbPI3fgtAQsfedsOoGnHbbbeRnJzM6tWrWbJkie04you03NUJ5a1dgHEVE3VWR8KTG9mOo2pAdHQ0AweWrd6t0yKDi86WUZUypS5yV5fNoEnoqNMfg9n999/PyJEj+WjGxzS67z+EJf7+6uOdL/a0lExVh47cVaUKtiylNO8gYbUbEXVmB9txVA1q0KBB2Qqfxk3u6rm24ygv0XJXlTo2/THtWl2zPQQcnRaZt3YB7uJCy2mUN2i5qz8o2ruZor2bcETGEtvqMttxlA906tSJiAbn4i7KJ3/D17bjKC/Qcld/kLOqbNQe1+4qHBFRltMoXzl6biUnfTbGuC2nUdWlJ1TV77hyMynYtBjEQXzHXifdVhcPCy4x516E89v/4Dq0m8Ida4g+q6PtSKoadOSufid39afgLiXm3IsIS9A120OJOMOI71D2D3pO+izLaVR1abmrY9wlheT98Bmgqz+Gqrh2VyFhkRTuWE1x5q+246hq0HJXx+Rv+AZ3YS4R9c8hosGfbMdRFjij44ltXXYS/eiMKRWYtNwVAMa4yVlZ9lY8Ia2PTn8MYUffteVv+JrSIzmW06jTpeWuACjcsQbXod0445KJOfci23GUReHJjYk6qyPGVXzsMJ0KPDpbJgidbBbLiS4lP3oCLb5jL8SpvxbBytMZTglpfSjcvorc1XMpLi4mIiKihpMpb9ORu2L9+vUU7liNhEcS1+5q23GUH4hKPY/wOk0ozTvEBx98YDuOOg1a7oqxY8cCENfmcpzR8ZbTKH8gIiScfx0Ao0eP1rXeA5CWe4jbv38/77zzDiDE6+qPqoLYlt1wxCSxZs0avvvuO9tx1CnScg9xR++fGX12J8JrN7QdR/kRCYsgvkPZOZrRo0dbTqNOVZXlLiJvisgBEVl/gsdFRMaJyDYR+VFEdH3YAHHkyBEmTZoEcOwtuFIVxZ93DZGRkcyZM4ctW7bYjqNOgScj97eAk51l6wGcXf5xL/Bq9WMpX5g2bRoZGRl07NiRyEatbMdRfsgZk8gdd9wBwJgxYyynUaeiynI3xiwEDp1kkz7A26bMciBJROp7K6CqGW63m1GjRgEwePBgvWhJndDgwYMBeOutt8jIyLCcRnnKG8fcGwK7Kny9u/x7fyAi94pIuoik6y+JXZ9++imbNm2icePG3HTTTbbjKD/WokULevXqRWFh4bHDeMr/eaPcKxvyVTpvyhgz2RiTZoxJS0lJ8cJLq9P18ssvA/DII48QHh5uOY3yd0OGDAFgwoQJFBQUWE6jPOGNct8NNK7wdSNgrxeeV9WQFStWsHDhQhITE7nnnntsx1EB4JJLLuH8888nMzOTt99+23Yc5QFvlPts4I7yWTOdgWxjzD4vPK+qIUdH7X/729+Ij9eLllTVROTY6H3UqFGUlpZaTqSq4slUyPeAZcC5IrJbRPqLyEARGVi+yTxgO7AN+F/g/hpLq6pt+/btzJgxg/DwcB5++GHbcVQAueGGGzjzzDPZtm0bs2bpzTz8XZUrRBljbq3icQM84LVEqkaNHj0at9tNv379aNhQL1pSVau42FhOsythx+vc+tATFF5/vc6y8mN6hWoIKc0/zBtvvAH89wSZUqcirs0VOKLiKd67mUWLFtmOo05Cyz2E5K6aQ2FhIb169aJ169a246gA5IiIOnbj9BdffNFyGnUyunB3iHAXFZC7ei4AQ4cOtZxGBbL4Dr3IWfEx8+fPp8HdE4ioe+YftjnRfQOU7+jIPUTkrf0Md1E+F198MRddpHdaUqfPGZNIXNurAMj+/iPLadSJaLmHAOMqIWflTEBH7co7Ei64DhxOCn5aRMnh32zHUZXQcg8B+Ru/oTTvEOF1mnLNNdfYjqOCQFhCXWJbdgXjJnflJ7bjqEpouQc54y4le3nZW+eEzn116prymoQL/gJA3o9fUJqXZTmNOp6We5Ar2LQYV9ZewhLrEdviEttxVBCJSGlK9DldMK5icnT07ne03IOYMW6yl5Xd3Dih842Iw2k5kQo2iV1uBiB3zTxKj+RYTqMq0nIPYke2fk9J5i844+sQ17q77TgqCEWe0ZyoMztiSgrJTZ9tO46qQOe5ByljDNnLpgOQ0OkvSJgu66tqRuKFN1O4YxU5q+aQcMH1OCJjf7dkQUU6/913dOQepAp3rKb4t204YpKIa3ul7TgqiEU1aklkkzaYonxyV1de6sr3tNyDkDGG7KXvA5SNpMIjLSdSwe7osfeclTNxF+nNPPyBHpYJQoU711C05ycc0QnEt+/xu8dO9HZZqeqIatqOyIYtKNrzE7mr55LYRW/daJuO3IOMMYbDi6cBkHDBDTgiYywnUqFAREi8+HYAclZ8oqN3P6DlHmQWLFhA8d7NZaP2DnrySvlOVNN2RDZqibswl9xVc2zHCXla7kHEGMMzzzwDQEKnvjgioi0nUqHk96P3j3EX5VtOFNq03IPIvHnzWLFiBY6YJOLP0zVklO9FNWlLZOPWuIvyydF571ZpuQeJiqP2xM59cUREWU6kQpGIkHR09L5yJqWFeZYThS4t9yAxY8YMVq1aRf369Ylrf7XtOCqERTVpQ1TTtpiifHJ0vXdrtNyDgMvl4qmnngLg6aefxhGuo3ZlV9IldwKQmz4HV94hy2lCk85zDwJTp05l8+bNNGvWjP79+/PiPz+3HUmFuMgG5xJ9TheObFlG9tL3Sb7yfuDk11no0gTepSP3AHfkyBGGDRsGwHPPPUd4uK4ho/xD0p/7gTjIW7uAkqx9tuOEHC33ADdp0iR2795Nu3btuPnmm23HUeqYiDpNiG11GbhLObz4HdtxQo6WewDLysri+eefB+D555/H4dAfp/IvSRffBs4wCjYupHj/z7bjhBRtgwA2YsQIDh06RNeuXenRo0fVf0ApHwtLrEv8eT0Bw6Gv38AYYztSyNByD1A///wz48aNQ0QYPXq03htV+a3Ei27FERVP0a8/cmTbCttxQoaWe4AaOnQoJSUl9OvXjw4dOtiOo9QJOaPiSLzoFgCyvn0TU+qynCg0aLkHoMWLF/PRRx8RHR3NiBEjbMdRqkrx511DWK0GuA7tIfeH+bbjhASPyl1ErhaRzSKyTUSGVvL4XSKSISI/lH8M8H5UBeB2u3nssccAGDJkCI0aNbKcSKmqiTOcWt3+B4Dsxe/qsgQ+UGW5i4gTmAj0AFoCt4pIy0o2nW6MaV/+McXLOVW5qVOnsmLFCs444wwef/xx23GU8lj02Z2JbNIGd2Eu2Yt0amRN82TkfgGwzRiz3RhTDLwP9KnZWKoyWVlZ/P3vfwfg3//+N3FxcZYTKeU5EaH25X8DcZC7Zh7F+7fbjhTUPFl+oCGwq8LXu4FOlWz3FxG5BNgCDDbG7KpkG1UNzzzzDBkZGVx88cXcfvvtess8FXAiUlKJ79CL3FWzOfTFa9S7faTO9KohnozcK/s/f/xk1TlAqjGmLfAlMLXSJxK5V0TSRSQ9IyPj1JKGuLVr1zJx4kQcDgcTJkzQvxAqYCVdfBuOmCSK9mwkf8M3tuMELU/KfTfQuMLXjYC9FTcwxhw0xhSVf/m/QMfKnsgYM9kYk2aMSUtJSTmdvCHJGMODDz6I2+3mgQceoF27drYjKXXaHFFxx06uZn37pt6xqYZ4Uu4rgbNF5EwRiQBuAX53ixURqV/hy2uBn7wXUb3xxhssXryYlJQUhg8fbjuOUtUW2/pSIhu2wJ1/mKzvKn2jr6qpynI3xriAB4EFlJX2B8aYDSIyXESuLd/sYRHZICJrgYeBu2oqcKjZt28fQ4YMAeCVV14hKSnJciKlqk/EQe2rHgCHk7w18yjcvcF2pKAjttZ6SEtLM+np6VZeO5D07duXGTNmcM011zB37tzfHWvXE6oq0B1e9A7ZS98nrHYjcvdsJSpKbzRTFRFZZYxJq2o7vULVj33yySfMmDGD2NhYXn31VT2JqoJOYpebCavdCNeh3cdWOFXeoSN3P3X48GFatWrF3r17qXX530jo2Nt2JKVqROHuDeyf9ndwOKl/51gi6p75u8f1Dk2/pyP3APfggw+yd+9eIhv8ifjzrrEdR6kaE9WoFXHnXQPuUjLnjsK4SmxHCgpa7n5o+vTpTJs2jZiYGJJ7DkYcTtuRlKpRtbr9D2FJ9SnJ2MnhxdNsxwkKWu5+Zs+ePdx3330AjBo1ivDaDS0nUqrmOSKiqdPrURAHOd/PoHDXetuRAp6Wu59IHfopTf8+l+bHrYvGAAALuklEQVR/vpasrCyizurICzt0xUcVOiIbtiCh842AIfPTMbiLCmxHCmha7n4kd+UnFO5YjSM6geQeg3R2jAo5SRfdQkS9ZpRm7+fgggl6W75q8GThMOUDhbs3kvXtWwAk9xhEWFxtu4GUskCc4dTpPYR9UwdT8NNC8hq3JnVo5YMcnUVzcjpy9wMZGRlkzhoJxk3CBTcQc3Zli24qFRrCkxuTfPVDABz6ajJFv22znCgwablb5na7+etf/0pp3kEiG7Yk6ZI7bEdSyrrYll3LpkeWusic+YLeuek0aLlb9sQTT/D555/jiE6gzrWPI049UqYUQO3LBhBRrxmu7P1kzn4J4y61HSmgaJP4WMX1YPLWfcXBeWNAHNS59nHCEupYTKaUf5GwCOpc9w9+e/tRCnesJuvrKWV3clIe0ZG7JYW7N3Dws/EA1L5iINGp7S0nUsr/hCedQcr1T4AjjNxVc8hdM892pICh5W5BSdY+Mj4eAW4X8R176/ICSp1EVOPW/z3B+sVrHNmxxnKiwKDl7mOu3EwOTH8K95Ecos7sSK3LBtiOpJTfi2vTnYROfcG4yfhkBEV7N9uO5Pe03H0oMzOTA9OfxpW9n4j655DS5++6boxSHkrqegexrS/DlBRy4MNnWLdune1Ifk3L3UdycnLo0aMHJQd/JbxOU+re+CyOyBjbsZQKGCIOknsMIvrszrgL87jiiivYunWr7Vh+S8vdBw4ePEj37t1JT08nLOkM6t78HM7oeNuxlAo44nCScu3jRDVtx/79++natSsbN260Hcsv6VTIGlBxuqMr7xAHpj9FSeavhCWdQb1bRujSAkpVg4RFkHLDUzRbM5Fvv/2WSy65hAULFtCxY0fb0fyKjtxrUMnh39g/7e+UZP5KeHIT6t02krDEerZjKRXwHBHRzJs3j549e3Lw4EEuvfRSvvvuO9ux/IqWew0p/HUdv739KK7D+4g4ozn1bnuBsPhk27GUChotnv2adS0GEPOnP5Obm0u3yy4nuccgvXF8OT0sUwNyf5jPoS9eA3cpUWd1JOXax3FExtqOpVTQObqKZFZcbXLTZ3Hos3GUZOzE5bqKsLDQrjcduXtRXl4eAwYM4NCCieAuJeGCG6j7l6e12JWqQeJwUrv7PST3eLj8StbZXHHFFezatct2NKu03L3k+++/p3379rzxxhvgDCf5msHUuvRunceulI/Etb2SereOwBGTxLfffkvbtm354IMPbMeyRsu9mvLz83nyySe56KKL+Pnnn2nbti317xxLXJvutqMpFXKiGrWiwd3j6dmzJ4cPH+bmm2/m9ttvZ9++fbaj+ZyW+2kyxjBjxgxatGjB888/T2lpKY899hgrVqwgIqWp7XhKhSxnbC3WtRpI7SvvR8Iieffdd2mY2oza3e+hpKTEdjyf0XI/RcYYvvrqK7p27Urfvn3ZtWsXHTp0YOnSpbz88stERkbajqhUyBMR4s+7hvr9JxLd/AJM8RGyvp5CmzZteOedd3C5XLYj1jgtdw+5XC5mz57NxRdfzOWXX86iRYuoVasWkyZNYsWKFXTp0sV2RKXUccKTzqDuX54mpe8zhCXVZ/PmzfTr148WLVowZcoU8vPzbUesMVruVdixYwfPPPMMqamp9OnTh6VLl1K7dm1GjBjBzp07ue+++3A69aSpUv4sptn5NBjwKm+++SbNmjVj27Zt3HPPPTRo0IAHHniAVatWYYyxHdOrxJMdEpGrgVcAJzDFGPPicY9HAm8DHYGDwM3GmJ0ne860tDSTnp5+mrFrjsvlYs2aNXz66afMnDmTtWvXHnusefPmDBw4kHvvvZf4+Hi9WEKpAGTcpeT/tJC81Z9StHfTse83bdqU6667jt69e9OlSxdiYvxzYT8RWWWMSatyu6rKXUScwBbgCmA3sBK41RizscI29wNtjTEDReQW4HpjzM0ne15/KPeSkhK2bt3K+vXrWbduHcuXL2fZsmW/e6sWHx9P7969GTBgAN26dUNEjj2m5a5UYCvO2Ene2gUUbFpMaX7Wse+Hh4fTsWNHLrzwQtq2bUvr1q1p0aKFXxS+N8u9CzDMGHNV+df/ADDGvFBhmwXl2ywTkTDgNyDFnOTJT7fcf/31Vw4cOHDsLZQxBrfbjdvtprS0lJKSEoqLiykuLqagoID8/Hzy8vI4dOgQmZmZZGZmsmfPHn755Rf27t2L2+3+w2uE1WpAVJO2xJzdmaim7ZCw8FPOqZQKHMa4ef+6OnzyySd8+eWXrF27ttJuqFu3Lk2bNqVx48akpKSQnJxMcnIycXFxxMXFERsbS0REBBEREYSHh+N0OnE6nTgcjmMDQxEhNjaWFi1anFZWT8vdk+tzGwIVL/XaDXQ60TbGGJeIZAPJQKZncT03cuRIJk2a5JXnEhHOPPNM2rRpw7cZUUTUa0ZUo5Y4Y2t55fmVUoFBxEGXLl2OTYzIyclh2bJlrFixgvXr17N+/Xq2bNnCgQMHOHDgACtXrqzW63Xq1Inly5d7I/oJeVLuUsn3jh+Re7INInIvcG/5l3kicir3yqqDl/+xMMawY8cOduzY4c2n9Tav73cACdV91/22QEb67rW+//77iod4T3W/PbqQxpNy3w00rvB1I2DvCbbZXX5YJhE4dPwTGWMmA5M9CXY8EUn35K1IsAnV/YbQ3Xfd79BSU/vtyVTIlcDZInKmiEQAtwCzj9tmNnBn+ed9ga9PdrxdKaVUzapy5F5+DP1BYAFlUyHfNMZsEJHhQLoxZjbwBvB/IrKNshH7LTUZWiml1Ml5tOCxMWYeMO+47z1d4fNC4EbvRvuD0zqcEwRCdb8hdPdd9zu01Mh+e3QRk1JKqcCiyw8opVQQ8rtyF5GrRWSziGwTkaGVPB4pItPLH/9eRFJ9n9L7PNjvR0Vko4j8KCJfiUhQrCtc1X5X2K6viBgRCYrZFJ7st4jcVP4z3yAi7/o6Y03x4He9iYh8IyJryn/fr7GR05tE5E0ROSAi60/wuIjIuPL/Jz+KSIdqv6gxxm8+KDth+zNwFhABrAVaHrfN/cBr5Z/fAky3ndtH+30pEFP++X2hst/l28UDC4HlQJrt3D76eZ8NrAFqlX9d13ZuH+77ZOC+8s9bAjtt5/bCfl8CdADWn+Dxa4D5lF0z1Bn4vrqv6W8j9wuAbcaY7caYYuB9oM9x2/QBppZ//hHQXSou+BKYqtxvY8w3xpiC8i+XU3a9QaDz5OcN8BzwElDoy3A1yJP9vgeYaIzJAjDGHPBxxpriyb4bIKH880T+eF1NwDHGLKSSa38q6AO8bcosB5JEpH51XtPfyr2ypQ4anmgbY4wLOLrUQSDzZL8r6k/Zv/KBrsr9FpHzgMbGmLm+DFbDPPl5nwOcIyJLRGR5+cqswcCTfR8G/FVEdlM2S+8h30Sz6lQ7oEoeTYX0Ia8tdRBgPN4nEfkrkAZ0rdFEvnHS/RYRBzAGuMtXgXzEk593GGWHZrpR9i5tkYi0NsYcruFsNc2Tfb8VeMsYM6p84cL/K9/3P67kFTy83mv+NnI/laUOONlSBwHGk/1GRC4HngSuNcYU+ShbTapqv+OB1sC3IrKTsmORs4PgpKqnv+ezjDElxpgdwGbKyj7QebLv/YEPAIwxy4AoytZfCWYedcCp8LdyD9WlDqrc7/LDE69TVuzBcvz1pPttjMk2xtQxxqQaY1IpO9dwrTHG/+7ycmo8+T2fSdlJdESkDmWHabb7NGXN8GTffwW6A4hIC8rKPcOnKX1vNnBH+ayZzkC2MWZftZ7R9lnkE5w13kLZGfUny783nLK/1FD2g/4Q2AasAM6yndlH+/0lsB/4ofxjtu3Mvtjv47b9liCYLePhz1uA0cBGYB1wi+3MPtz3lsASymbS/ABcaTuzF/b5PWAfUELZKL0/MBAYWOHnPbH8/8k6b/ye6xWqSikVhPztsIxSSikv0HJXSqkgpOWulFJBSMtdKaWCkJa7UkoFIS13pZQKQlruSikVhLTclVIqCP0/gLNW8YLbuvQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "obs = beta.rvs(5, 5, size=2000)\n",
    "grid = np.linspace(0.01, 0.99, 100)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(obs, bins=40, normed=True)\n",
    "ax.plot(grid, beta.pdf(grid, 5, 5), 'k-', linewidth=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Goodies in scipy.stats\n",
    "\n",
    "There are a variety statistical functions in `scipy.stats`\n",
    "\n",
    "For example, `scipy.stats.linregress` implements simple linear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.0061459977845484, -0.005573394248197847)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import linregress\n",
    "\n",
    "x = np.random.randn(200)\n",
    "y = 2 * x + 0.1 * np.random.randn(200)\n",
    "gradient, intercept, r_value, p_value, std_err = linregress(x, y)\n",
    "gradient, intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the full list, consult the [documentation](http://docs.scipy.org/doc/scipy/reference/stats.html#statistical-functions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Roots and Fixed Points\n",
    "\n",
    "A *root* of a real function $ f $ on $ [a,b] $ is an $ x \\in [a, b] $ such that $ f(x)=0 $\n",
    "\n",
    "For example, if we plot the function\n",
    "\n",
    "\n",
    "<a id='equation-root-f'></a>\n",
    "$$\n",
    "f(x) = \\sin(4 (x - 1/4)) + x + x^{20} - 1 \\tag{2}\n",
    "$$\n",
    "\n",
    "with $ x \\in [0,1] $ we get"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAHVCAYAAADPfJ4+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8VeXhx/Hvk51ABpAQVsLeGwIouEUFZ92IoCCKoHW3ddbaWrWtbbVuUZSlDHGhorgnIiQk7BVmwkxIAhlk3uf3R4I/kACBJPfc8Xm/XnmRe88x59ue5Oab5zz3OcZaKwAAALhHgNMBAAAA/AnlCwAAwI0oXwAAAG5E+QIAAHAjyhcAAIAbUb4AAADciPIFAADgRpQvAAAAN6J8AQAAuFGQ0wGOJTY21rZp08bpGAAAAMeVkpKSba2NO95+Hl2+2rRpo+TkZKdjAAAAHJcxZmtN9uOyIwAAgBtRvgAAANyI8gUAAOBGlC8AAAA3onwBAAC4EeULAADAjShfAAAAbkT5AgAAcCPKFwAAgBtRvgAAANyI8gUAAOBGlC8AAAA3onwBAAC4EeULAADAjShfAAAAbkT5AgAAcCPKFwAA8Fl5RaUqKi13OsZhKF8AAMBnvf7DZvV//EsVl1U4HeVXlC8AAOCzFm3aqy7NIxUWHOh0lF9RvgAAgE8qKi3Xssw8ndKuidNRDkP5AgAAPilla67KKizlCwAAwB0WbdqroACjpNaNnI5yGMoXAADwSYs25ahXq2g1CA1yOsphKF8AAMDnFJaUa1mG5833kihfAADAB6VszVW5y/Pme0mULwAA4IN+ne/VxrPme0mULwAA4IMWbdqr3gkxigjxrPleEuULAAD4mMKSci3P3KdT2jV2Okq1KF8AAMCnJHvwfC+J8gUAAHzMok17FRxo1N/D1vc6iPIFAAB8yqJNe9W7lWfO95IoXwAAwIcU/DrfyzMvOUqULwAA4EOSt+SowoPne0l1VL6MMW8YY/YYY1YeZbsxxjxnjEk3xiw3xvSri+MCAAAcatGmHI+e7yXV3cjXFEnDjrF9uKSOVR/jJb1cR8cFAAD41aJNe9UnIUbhIYFORzmqOilf1trvJeUcY5fLJE2zlRZJijHGNK+LYwMAAEiV871WbPfs+V6S++Z8tZSUccjjzKrnjmCMGW+MSTbGJGdlZbklHAAA8H5LvGC+l+S+8mWqec5Wt6O1dpK1NslamxQXF1fPsQAAgK9YtGmvQgID1C/Rc+d7Se4rX5mSEg553ErSDjcdGwAA+IFv1u5Rv9aePd9Lcl/5mifphqp3PZ4iaZ+1dqebjg0AAHzc5uxCrd9doAu6N3M6ynHVydKvxpiZks6SFGuMyZT0F0nBkmStfUXSfEkXSkqXVCRpbF0cFwAAQJIWrNolSTrfX8qXtfa642y3km6vi2MBAAD81mcrd6lXq2i1jAl3OspxscI9AADwarv2FSstI88rLjlKlC8AAODlvlhdecnxgu7xDiepGcoXAADwagtW7Va7uAbq0DTS6Sg1QvkCAABeK6+oVIs27fWaS44S5QsAAHixr9bsUbnLahjlCwAAoP4tWLVLzaPD1KtVtNNRaozyBQAAvNKB0gp9vyFL53eLlzHV3cnQM1G+AACAV/pufZaKy1xeNd9LonwBAAAv9fmqXYqJCNbAto2djnJCKF8AAMDrlFW49OWa3RraNV5Bgd5VZ7wrLQAAgKRFm/Zqf3G5111ylChfAADACy1YtUsRIYE6vWOs01FOGOULAAB4lZLyCs1fsUtnd26qsOBAp+OcMMoXAADwKp+t3KWcwlJdNzDR6SgnhfIFAAC8yluLtqlNkwgNbt/E6SgnhfIFAAC8xobd+Vq8JUfXDUxUQID3LKx6KMoXAADwGm/9sk0hgQG6qn8rp6OcNMoXAADwCgdKK/Tu0kwN79lMTRqGOh3npFG+AACAV/ho+Q7lF5fr+kGtnY5SK5QvAADgFd76ZZs6NG2oAW0aOR2lVihfAADA463cvk/LMvJ0/aBEGeOdE+0PonwBAACP9/bibQoLDtAVfb13ov1BlC8AAODRCkrK9WHqdl3Sq4WiI4KdjlNrlC8AAODRPkjdrsLSCl1/indPtD+I8gUAADyWtVYzFm1Vt+ZR6t0q2uk4dYLyBQAAPNanK3dp7a58jTutrddPtD+I8gUAADxSeYVL//58nTrFN9Tv+rZ0Ok6doXwBAACP9N7S7dqUVaj7zu+sQC+9j2N1KF8AAMDjFJdV6Nkv16t3QozO7xbvdJw6RfkCAAAe561ftmnHvmL96YLOPjPX6yDKFwAA8CgFJeV66Zt0DenQREM6xDodp85RvgAAgEd548fN2ltYqj9e0MXpKPWC8gUAADxGbmGpXvt+ky7oHq8+CTFOx6kXlC8AAOAxXv5uowpKy/WH8zs7HaXeUL4AAIBH2JhVoKkLt+jyvi3VMT7S6Tj1hvIFAAAcV1bh0j2z0xQeEqgHhvnmXK+DgpwOAADwLy6XVUFpufYfKFNBSblKy10qq7Aqq3CpvMKqzOVSgDEKCqj6CDQKDAhQg5BANQwLUsPQIDUICVKADy26Cem5rzZoeeY+vXx9PzWNCnM6Tr2ifAEA6oS1Vln5JdqcXaiM3APavb9Ye/YXa/f+Eu3OL1ZWfon2HyhTfkm5rK3dsYyRGoYEKaZBsGIbhiq2YajiIiv/bRYVpoTG4WrVKEItYsIUGhRYN/8DUW9StubqxW/SdWW/Vhres7nTceod5QsAcEJcLquM3CKt2blfq3fmK31PvjZnF2nr3kIVlVYctm9UWJDio8IUHxWmNm0aKDo8WFFhQYoKD1ZkWJAahgYrNChAwUEBCg4wCg4KUGCAkbWV9/WrcFmVuypHxYpKK1RQUq6C4nLll5Qrv7hMuYWlyi4o1ba9RVq6NVc5RaWHFTtjpPjIMCU2iVD7uIbq0PT/P1pEh/nc4p3eqLCkXPfOSVPz6HA9dmk3p+O4BeULAHBU1lrt2FeslK25Wro1Vyu279PanftVWFWyAozUukkDtWkSoVPaNVbb2AZq3aSBEhtHqFlUmMJD3DvqVF7h0u78EmXmFCkj94Ayc4uUmXtAW7IL9enKncorKvt13wYhgeraPErdW0Spe4todWsRpU7xkQoJYjq0Oz3+8WptyynS7PGnKjIs2Ok4bkH5AgD8ylqr9D0F+jE9W8lbcpWyNVe79hdLksKDA9WjZZSu6t9KXZtHqWvzyrLi7oJ1LEGBAWoZE66WMeEa9Jtt1lrtLSzVxj0FSs8q0Ppd+Vq1Y7/mpmRq6s9bJUkhgQHq1iJKfRNj1DexkfomxKhVo3BGyOrJF6t3a9aSDE04s70Gtm3sdBy3Mba2F97rUVJSkk1OTnY6BgD4tD35xfpxQ7Z+TM/WT+nZ2r2/RJLUMiZc/Vs3+vWjS7NIBQX63qiQy2W1ZW+hVu/crxWZ+5SakaflmXkqLnNJkmIbhmpg20Ya1LaJBrVrrE5NI5nsXwcyc4t02Qs/qWlUmD68fYhPjDgaY1KstUnH24+RLwDwM9Zabcwq0IJVu/XF6t1Ky8iTJDWKCNaQDrE6rUOshnSIVULjCIeTukdAgFG7uIZqF9dQF/dqIaly2YN1u/KVmpGnpVtztXhzjuav2CVJiokI1oA2jTWkfROd1jFW7eMaMjJ2gvbsL9b1r/+isgqXnhvRxyeK14lg5AsA/IC1Vqt27NdHy3bo89W7tTm7UJLUu1W0zusWr7M6N1W35lGM6BxDRk6RFm/O0S+b92rRphxtyymSJDWNDP21sJ7eKVZNI317mYTayiks1YhJPysz94Bm3DxI/RIbOR2pztR05IvyBQA+bHN2oeal7dCHy7ZrU1ahggONTm0fq/O6xeu8rvFqFk1ROFkZOUX6Kb3ycu3CjXuVU1gqSerWPEpndY7TmZ3i1K91IwX74KXak7W/uEwjX1ukDbsL9ObYARrcPtbpSHWK8gUAfmp/cZk+TNuhuckZWpa5T8ZIg9o21mV9Wmp4j2aKiQhxOqLPcbmsVu/cr+/WZ+m79VlK2ZqrCpdVZGiQTusYq3O6NNVZnZsqLjLU6aiOKSot1w2TFystI0+Tbuivc7rEOx2pzlG+AMCPWGu1dFuuZi7O0CfLd+pAWYW6NIvUlf1a6eLezdU8OtzpiH5lf3GZFqZn69t1Wfpm3Z5f38TQu1W0zukSr3O7NlX3FlF+M1fsQGmFbpmWrIUbs/X8df10US/fXEiV8gUAfqCwpFzvLc3UtJ+3asOeAjUICdSlfVpoxIBE9WoV7Te/3D3Zwfl236zdo6/W7tGyzDxZK8VHheqcLvEa2rWphnSIVViw5yzZUZc27M7X7W8v1YY9BfrXlb10dVKC05HqDeULAHxYRk6Rpv28RbOWZCi/uFy9WkVr5MBEXdK7hRqE8kZ2T5ZdUKJv12XpqzW79f36LBWWVigsOEBD2sfqnK5NdU6Xpj4xUmmt1TvJmXp03ko1DA3Sf6/pozM6xTkdq15RvgDAB6VszdVr32/S56t3yRijYT2a6aYhbdUvMYZRLi9UUl6hxZtz9NWaPfpq7W5l5ByQVDlp/9yqItarVYwCvexdqAUl5Xrk/RX6IG2HBrdvomdH9PGLd4FSvgDAR1hr9WN6tl78Jl2LNuUoOjxYIwclavQprdUixvtHSFDp4N0Fvlq7R1+v3fPrpP1GEcE6vWOczuocpzM6xSm2oedO2rfW6ocN2frLvFXaurdQ9wztpNvO7uB15fFkUb4AwMu5XFafr96tl75N1/LMfYqPCtUtp7fTdQMTubToB/KKSvX9hmx9u26Pvl+fpeyCyqUseraM1pAOsRrcvokGtGnsEbd3stbq67V79NzX6VqWkaeWMeH67zW9NahdE6ejuRXlCwC8lLVWX6zerf9+sV5rd+WrdZMITTizva7o11KhQc7/ooX7HVzKorKIZSs1I1dlFVYhgQHqmxijwe1jNaBNI/VOiHFrMT/4B8LzX2/Qqh371TImXLed3V5X9W/ll9+rlC8A8DIHL9n85/N1Wpa5T21jG+iuczvq4l7NffKeijh5RaXlWrw5Rws37tVP6dlavXO/rJUCjNS5WZT6JcaoX2IjdW0epfZNG9RpEdpXVKbvN1QuoXFwRK51kwjdfnYHXd63pV8vKkv5AgAvkrwlR//6bJ0Wb8lRy5hw3TW0o67o25LShRrZV1Sm1IxcLd2Wp9RtuUrblqf8knJJUmCAUesmEeocH6mO8ZFqFROuuKhQxTUMVdOoUDVpEFrtnKyi0nJtzz2gjNwiZeQcUGZukVK35Wnptly5bOU9Ls/oGKfzu8drWPdmfK+K8gUAXmFLdqH+8elafbZql5pGhuqOczromgEJfnnJBnWnwlV58/R1u/K1fvfBjwJt3Vso129+7QcYHXFja2ulknLXYc+FBgWoc7NIndkpTmd1bqo+Cd73Lsz6VtPyxYxNAHBAXlGpnvsqXdMXbVFwYIDuPa+Tbjm9nUdMnob3Cwww6hQfqU7xkYc9X1xWoaz8Eu3JL1FWfvGvn5ceWrSq+lRUWLBaNQpXQuMItWoUrriGoSxnUkfqpHwZY4ZJ+p+kQEmvW2v/8ZvtYyQ9LWl71VMvWGtfr4tjA4A3KatwadrPW/XcVxuUX1yma5ISdO95ndQ0yvfXQILzwoIDldA4QgmNI5yO4tdqXb6MMYGSXpR0nqRMSUuMMfOstat/s+tsa+3va3s8APBWizbt1aMfrtT63QU6vWOsHr6oq7o0i3I6FgA3q4uRr4GS0q21myTJGDNL0mWSflu+AMAv7d5frCc+WaN5y3aoZUy4Jo3ur/O6xXMJB/BTdVG+WkrKOORxpqRB1ex3pTHmDEnrJd1jrc2oZh8A8BnlFS5NWbhFz3yxXmUuqzvP6aCJZ3VgXhfg5+qifFX3p9tv30L5kaSZ1toSY8wESVMlnVPtFzNmvKTxkpSYmFgH8QDA/VZu36cH3luuldv36+zOcfrLJd3VJraB07EAeIC6KF+ZkhIOedxK0o5Dd7DW7j3k4WuS/nm0L2atnSRpklS51EQd5AMAtzlQWqFnv1yv13/crMYNQvTS9f00vEczLjEC+FVdlK8lkjoaY9qq8t2MIySNPHQHY0xza+3OqoeXSlpTB8cFAI/y44ZsPfT+Cm3LKdKIAQl6cHhXRUcEOx0LgIepdfmy1pYbY34vaYEql5p4w1q7yhjzN0nJ1tp5ku40xlwqqVxSjqQxtT0uAHiKgpJyPfHJGs1cvE1tYxto1vhTdIqf3VAYQM2xwj0A1MLPG/fqj3OXaXveAY0/vZ3uOa+TwoKZUA/4I1a4B4B6dKC0Qv/8bK2mLNyiNk0iNHfCqerfurHTsQB4AcoXAJygtIw83TM7TZuzCzVmcBvdP6wLy0cAqDHKFwDUUIXL6pXvNuqZL9araWSo3r5lkAa3j3U6FgAvQ/kCgBrYkXdA98xO0y+bc3Rxr+Z64vKeig7nnYwAThzlCwCO45PlO/Xge8tV4bL699W9dWW/lqzbBeCkUb4A4CiKyyr0149Wa+biberdKlr/G9GXVeoB1BrlCwCqsTm7ULe9tVRrdu7XrWe20x/O76zgwACnYwHwAZQvAPiNj5fv0APvrlBQoNGbYwbo7C5NnY4EwIdQvgCgSkl5hf7+8RpNX7RV/RJj9PzIfmoZE+50LAA+hvIFAKp8N+PEGSlalrlP489opz9ewGVGAPWD8gXA7y3cmK073k5VSblLr47urwu6N3M6EgAfRvkC4LestXr9h8166tM1ahfXUK+O7q/2cQ2djgXAx1G+APilwpJy/end5fpk+U4N79FMT1/dWw1DeUkEUP94pQHgdzJyinTLtGSt352v+4d10YQz27FoKgC3oXwB8CuLN+dowowUlVW49ObYgTqzU5zTkQD4GcoXAL8xe8k2PfLBSrVqFKHXb0xifhcAR1C+APi88gqXnpy/Vm/8tFmnd4zVC9f1U3QEN8UG4AzKFwCfll9cpt+/narv1mdp7JA2evjCrgpi/S4ADqJ8AfBZO/IO6KYpS7RhT4GeuqKnrhuY6HQkAKB8AfBNK7fv001TluhAaYWmjB2g0zsysR6AZ6B8AfA5X67erTtnpapRRIimTxykzs0inY4EAL+ifAHwKVMXbtFfP1qlHi2j9foNSWoaFeZ0JAA4DOULgE9wuaz+uWCtXv1uk87rFq//jeijiBBe4gB4Hl6ZAHi90nKX7n93ud5P3a7Rp7TWY5d2V2AAK9YD8EyULwBeraCkXBNnpOiHDdn6w/mddPvZHbhVEACPRvkC4LWy8ks0dspirdmZr39d2UvXDEhwOhIAHBflC4BX2ra3SKMm/6Ks/BK9dkN/ndMl3ulIAFAjlC8AXmftrv26YfJilVa49PYtg9Q3sZHTkQCgxihfALxKytZcjX1zsSJCgvTOraeqYzxreAHwLpQvAF7j+/VZunV6iuKjQjV93CAlNI5wOhIAnDDKFwCv8Mnynbp7dqo6NI3UtJsGKi4y1OlIAHBSKF8APN6c5Aw98O5y9UtspMljBig6PNjpSABw0ihfADzatJ+36NEPV+n0jrGaNDpJ4SGBTkcCgFqhfAHwWJO+36gn56/V0K7xevH6vgoNongB8H6ULwAex1qr575K1zNfrtfFvZrrmWv7KDgwwOlYAFAnKF8APIq1Vv/8bJ1e+W6jrurfSv+8shf3aQTgUyhfADyGtVaPf7xGb/y0WaNOSdTfLu2hAIoXAB9D+QLgEay1+utHqzVl4RaNHdJGj17cjRtkA/BJlC8AjrPW6rF5qzT1560ad1pbPXJRV4oXAJ9F+QLgKGutHv1wlaYv2qpbTm+rhy6keAHwbZQvAI5xuawenbdSMxZt061ntNMDw7tQvAD4PMoXAEdY+//Fa8KZ7XX/sM4ULwB+gYVzALjdwcn1MxZt061ntqN4AfArlC8AbmWt1ZPz12jKwi0ad1pbPTCMS40A/AvlC4DbWGv1rwXr9NoPmzVmcBve1QjAL1G+ALjNM1+s18vfbtT1gxL1l0tYxwuAf6J8AXCLF77eoOe+Tte1SQl6/LIeFC8AfovyBaDeTf5xs/79+Xpd3relnrqiJ7cMAuDXKF8A6tXMxdv0+MerNbxHMz19VS+KFwC/R/kCUG8+TNuuh95fobM6x+l/I/oqKJCXHADglRBAvViwapfunbNMg9o21iuj+iskiJcbAJAoXwDqwffrs3TH26nq1Spar984QGHBgU5HAgCPQfkCUKdStubo1ukpat+0oaaMGaiGodzFDAAORfkCUGfW7NyvsW8uUbPoME27aaCiI4KdjgQAHofyBaBObMku1OjJixUREqTp4wYqLjLU6UgA4JEoXwBqbde+Yo2a/IsqXC7NuHmgWjWKcDoSAHgsyheAWsktLNXoyb8ot7BUU8YOVIemkU5HAgCPxkxYACetqLRcY6cs0dacIk0ZO0C9E2KcjgQAHq9ORr6MMcOMMeuMMenGmAeq2R5qjJldtf0XY0ybujguAOeUVbh021tLtTwzT8+N6KvB7WOdjgQAXqHW5csYEyjpRUnDJXWTdJ0xpttvdhsnKdda20HSM5L+WdvjAnCOy2V1/9zl+nZdlp64vKeG9WjmdCQA8Bp1cdlxoKR0a+0mSTLGzJJ0maTVh+xzmaTHqj6fK+kFY4yx1tpjfeF169bprLPOOuy5a665RrfddpuKiop04YUXHvHfjBkzRmPGjFF2drauuuqqI7ZPnDhR1157rTIyMjR69Ogjtt9333265JJLtG7dOt16661HbH/kkUc0dOhQpaWl6e677z5i+5NPPqnBgwdr4cKFeuihh47Y/uyzz6pPnz768ssv9fe///2I7a+++qo6d+6sjz76SP/5z3+O2D59+nQlJCRo9uzZevnll4/YPnfuXMXGxmrKlCmaMmXKEdvnz5+viIgIvfTSS5ozZ84R27/99ltJ0r///W99/PHHh20LDw/Xp59+Kkl6/PHH9dVXXx22vUmTJnr33XclSQ8++KB+/vnnw7a3atVKM2bMkCTdfffdSktLO2x7p06dNGnSJEnS+PHjtX79+sO29+nTR88++6wkadSoUcrMzDxs+6mnnqqnnnpKknTllVdq7969h20/99xz9ec//1mSNHz4cB04cOCw7RdffLH+8Ic/SNIR33cS33uHfu89+MR/tCPvgBIaRejVReF6VXzv8b3H6x7fe3zv1VRdXHZsKSnjkMeZVc9Vu4+1tlzSPklNqvtixpjxxphkY0xyWVlZHcQDUJe+XLNbO/IOqFlUmFo2Cnc6DgB4HXOcwafjfwFjrpZ0gbX25qrHoyUNtNbeccg+q6r2yax6vLFqn73Vfc2DkpKSbHJycq3yAag77y3N1L1zlumiXs313Ii+CgwwTkcCAI9hjEmx1iYdb7+6GPnKlJRwyONWknYcbR9jTJCkaEk5dXBsAG7y/fos/Wnucg1u30T/vaY3xQsATlJdlK8lkjoaY9oaY0IkjZA07zf7zJN0Y9XnV0n6+njzvQB4jpXb92nijBR1jI/Uq6P7KzSIG2UDwMmq9YR7a225Meb3khZICpT0hrV2lTHmb5KSrbXzJE2WNN0Yk67KEa8RtT0uAPfIyCnS2ClLFBMRoiljBygyjPs1AkBt1Mkiq9ba+ZLm/+a5Rw/5vFjS1XVxLADuk1tYqhvfXKzScpdm3jJI8VFhTkcCAK/HCvcAqlVcVqFxU5coM/eA3rp5ELcNAoA6wr0dARyhwmV116xUpWbk6X/X9tGANo2djgQAPoPyBeAIf/9ktRas2q1HL+6m4T2bOx0HAHwK5QvAYd74cbPe/GmLxp3WVmOHtHU6DgD4HMoXgF99tnKXHv9ktYZ1b6aHL+zqdBwA8EmULwCSpNRtubprVqr6JMTo2RF9FMAiqgBQLyhfALRtb5Funpqs+KgwvXZDksKCWUQVAOoL5Qvwc3lFpRozZbEqrNWUsQMU2zDU6UgA4NMoX4AfKy13acKMFGXmHNCk0UlqF9fQ6UgA4PNYZBXwU9ZaPfjeCi3alKNnr+2jgW1ZywsA3IGRL8BPvfTtRr27NFN3D+2o3/Vt6XQcAPAblC/AD320bIeeXrBOl/dtqbvO7eh0HADwK5QvwM+kbM3Vfe8s04A2jfSPK3vKGJaUAAB3onwBfiQjp0jjpyWreXSYXh2dpNAglpQAAHejfAF+Ir+4TOOmLlFZhUtvjBmgxg1CnI4EAH6JdzsCfqC8wqU7ZqZqY1ahpt00UO1ZUgIAHMPIF+AHnpi/Rt+uy9LfLuuuIR1inY4DAH6N8gX4uLd+2ao3f9qisUPa6PpBrZ2OAwB+j/IF+LCf0rP16IerdFbnOD1yUTen4wAARPkCfNamrAJNnJGi9nEN9Px1fRUYwJISAOAJKF+AD9pXVKabpyYrKDBAk28coMiwYKcjAQCqUL4AH1Ne4dLvZy5VRm6RXhnVXwmNI5yOBAA4BEtNAD7myflr9cOGbP3zyp7cLBsAPBAjX4APmb1km974abNuGtJW1w5IdDoOAKAalC/ARyzenKNHPlipMzrF6aELuzgdBwBwFJQvwAdk5BRpwowUJTSO0PPX9VVQID/aAOCpeIUGvFxhSblumZas8gqXXr8hSdHhvLMRADwZE+4BL+ZyWf3hnWVavztfU8YOVDvu2QgAHo+RL8CLPf91uj5duUsPXdhVZ3SKczoOAKAGKF+Al/ps5S498+V6XdGvpcad1tbpOACAGqJ8AV5o7a79undOmnonxOjJy3vKGG4dBADegvIFeJncwlLdMi1ZDUODNGl0f4UFBzodCQBwAphwD3iR8gqXbn97qXbvL9Hs8acoPirM6UgAgBPEyBfgRZ76dK0WbtyrJ37XQ30TGzkdBwBwEihfgJd4b2mmJv+4WWMGt9HVSQlOxwEAnCTKF+AFlmfm6YH3VuiUdo318EVdnY4DAKgFyhfg4bLyS3Tr9BTFNQzViyP7KZhbBwGAV2PCPeDByipcuv2tpcotKtXcCYPVpGGo05EAALVE+QI82OMfr9biLTn634g+6tEy2uk4AIA6wPULwEPNSc7QtJ+3avwZ7XRZn5ZOxwEA1BHKF+CBlmXk6ZEPVuq0DrH60wWdnY4DAKjYHlxxAAAbw0lEQVRDlC/Aw2QXlGjCjBQ1jQzV89f1VRAT7AHApzDnC/AgByfY5xSW6t2Jg9WoQYjTkQAAdYzyBXiQJ+ev0S+bc/TMtb2ZYA8APorrGYCHeG9ppt78aYtuGtJWl/dt5XQcAEA9oXwBHmDl9n16sGoF+wcv7OJ0HABAPaJ8AQ7LLSzVhBkpatwgRC+wgj0A+DzmfAEOqnBZ3TkrVXv2l2jOhFMVywr2AODzKF+Ag/77xTr9sCFbT13RU30SYpyOAwBwA65vAA5ZsGqXXvxmo0YMSNB1AxOdjgMAcBPKF+CAjVkFum/OMvVuFa3HLu3udBwAgBtRvgA3Kygp163TUxQaFKCXR/VXWHCg05EAAG7EnC/Ajay1+tPcZdqUVaAZ4wapRUy405EAAG7GyBfgRq//sFnzV+zSn4Z10eAOsU7HAQA4gPIFuMnPG/fqH5+t1bDuzXTrGe2cjgMAcAjlC3CDXfuKdcfMpWrdJEJPX91LxhinIwEAHEL5AupZablLt72VoqLSCr06qr8iw4KdjgQAcFCtypcxprEx5gtjzIaqfxsdZb8KY0xa1ce82hwT8DZPfLJaS7fl6emreqtjfKTTcQAADqvtyNcDkr6y1naU9FXV4+ocsNb2qfq4tJbHBLzGB6nbNfXnrbr5tLa6qFdzp+MAADxAbcvXZZKmVn0+VdLvavn1AJ+xdtd+PfDecg1s21gPDO/idBwAgIeobfmKt9bulKSqf5seZb8wY0yyMWaRMeaYBc0YM75q3+SsrKxaxgOcsb+4TBOmpygqLFgvjOyroECmVwIAKh13kVVjzJeSmlWz6eETOE6itXaHMaadpK+NMSustRur29FaO0nSJElKSkqyJ3AMwCNYa/WHOcuUmXtAM8efoqaRYU5HAgB4kOOWL2vt0KNtM8bsNsY0t9buNMY0l7TnKF9jR9W/m4wx30rqK6na8gV4u1e+26TPV+/Wny/upgFtGjsdBwDgYWp7LWSepBurPr9R0oe/3cEY08gYE1r1eaykIZJW1/K4gEdauDFbTy9Yq4t7NddNQ9o4HQcA4IFqW77+Iek8Y8wGSedVPZYxJskY83rVPl0lJRtjlkn6RtI/rLWUL/icXfuKdefMVLWLa6h/XslCqgCA6tXqxtrW2r2Szq3m+WRJN1d9vlBSz9ocB/B0BxdSPVBaoVnj+6lBKPesBwBUj98QQB146tM1WrotTy+M7KsOTVlIFQBwdLz/Hailj5fv0Js/bdHYIW10ca8WTscBAHg4yhdQC+l78nX/3OXq37qRHhze1ek4AAAvQPkCTlJhSbkmzFiqsOBAvTiyn0KC+HECABwfc76Ak2Ct1QPvrdCmrALNGDdIzaJZSBUAUDP8qQ6chKkLt+ijZTt03/mdNbhDrNNxAABehPIFnKCl23L1xPw1OrdLU008s73TcQAAXobyBZyAnMJS3f7WUjWLDtN/r+mjgAAWUgUAnBjmfAE1VOGyumtWqvYWluq9iYMVHRHsdCQAgBdi5Auooee+2qAfNmTrr5d2V4+W0U7HAQB4KcoXUAPfrtuj577eoCv7tdKIAQlOxwEAeDHKF3Ac2/MO6J7ZaeocH6m//64HN8wGANQK5Qs4hsobZi9VWYXVy6P6Kzwk0OlIAAAvx4R74Bie+GS1lmXk6ZVR/dQ2toHTcQAAPoCRL+Ao5i3boak/b9XNp7XVsB7NnY4DAPARlC+gGul78vXAu8uV1LqR7h/exek4AAAfQvkCfuPgDbPDgwP1wsh+Cg7kxwQAUHeY8wUcwlqrh99foY3cMBsAUE/4kx44xIxftumDtB26d2gnDeGG2QCAekD5Aqosy8jT4x+t1pmd4nT72R2cjgMA8FGUL0BSXlGpbntrqeIiQ/XstdwwGwBQf5jzBb/nclndO2eZ9uQX650Jg9WoQYjTkQAAPoyRL/i9l7/bqK/X7tGfL+6mPgkxTscBAPg4yhf82sKN2frP5+t0Se8WGn1Ka6fjAAD8AOULfmv3/mLdOTNVbWMb6KkrenLDbACAWzDnC36pvMKlO2amqrCkQm/fcooahvKjAABwD37jwC89/fk6Ld6co2eu7a1O8ZFOxwEA+BEuO8LvfL5ql179bpOuH5Soy/u2cjoOAMDPUL7gV7btLdJ97yxTz5bR+vPF3ZyOAwDwQ5Qv+I3isgrd9naKjKSXru+nsOBApyMBAPwQc77gN/728Wqt3L5fr9+QpITGEU7HAQD4KUa+4BfeW5qpt3/ZpolntdfQbvFOxwEA+DHKF3zeul35evj9lRrUtrHuO6+T03EAAH6O8gWfVlBSrolvpahhWJCeH9lXQYF8ywMAnMWcL/gsa63uf3e5tu4t0ls3D1LTyDCnIwEAwMgXfNfUhVv0yfKd+uMFnXVKuyZOxwEAQBLlCz5q6bZcPTF/jYZ2jdf409s5HQcAgF9RvuBzcgpL9fu3lqpZdJj+c3VvBQRww2wAgOdgzhd8istldffsNGUXluq9iYMVHRHsdCQAAA7DyBd8yvNfp+v79Vn6yyXd1KNltNNxAAA4AuULPuP79Vl69qv1uqJvS40cmOh0HAAAqkX5gk/YkXdAd81KVcemDfX3y3vIGOZ5AQA8E+ULXq+03KXb316qsgqrl0f1V0QIUxkBAJ6L31Lwek/OX6PUbXl6cWQ/tY9r6HQcAACOiZEveLWPl+/QlIVbNHZIG13Uq7nTcQAAOC7KF7zWxqwC3T93ufolxujB4V2djgMAQI1QvuCVCkvKNWF6ikKDA/XCyH4KCeJbGQDgHZjzBa9jrdVD769QelaBpt80SC1iwp2OBABAjTFcAK8zY9FWfZi2Q/cO7aTTOsY6HQcAgBNC+YJXScvI098+Xq2zO8fp9rM7OB0HAIATRvmC18gpLNVtM1LUNDJMz1zbhxtmAwC8EnO+4BV+vWF2QanmTjxVMREhTkcCAOCkMPIFr/C/rzZU3jD70m7q1SrG6TgAAJw0yhc83jfr9ui5rzfoin7cMBsA4P0oX/BoGTlFumd2mjrHR+qJ3/XkhtkAAK9H+YLHKi6r0G1vLVVFhdUro/orPCTQ6UgAANQaE+7hsf760Wqt2L5Pk0b3V5vYBk7HAQCgTtRq5MsYc7UxZpUxxmWMSTrGfsOMMeuMMenGmAdqc0z4h3eSMzRz8TZNPKu9zu/ezOk4AADUmdpedlwp6QpJ3x9tB2NMoKQXJQ2X1E3SdcaYbrU8LnzYqh379MgHK3Vquya677xOTscBAKBO1eqyo7V2jaTjTYIeKCndWrupat9Zki6TtLo2x4Zv2ldUpokzliomIljPXddXQYFMSwQA+BZ3/GZrKSnjkMeZVc9Vyxgz3hiTbIxJzsrKqvdw8Bwul9U9c9K0c98BvXR9f8VFhjodCQCAOnfc8mWM+dIYs7Kaj8tqeIzqhsXs0Xa21k6y1iZZa5Pi4uJqeAj4ghe/SdfXa/fozxd3U//WjZyOAwBAvTjuZUdr7dBaHiNTUsIhj1tJ2lHLrwkf8936LP33y/X6XZ8WGn1Ka6fjAABQb9xx2XGJpI7GmLbGmBBJIyTNc8Nx4SUyc4t016xUdY6P1JNXsJAqAMC31XapicuNMZmSTpX0iTFmQdXzLYwx8yXJWlsu6feSFkhaI2mOtXZV7WLDV/y6kKqrciHViBCWngMA+LbavtvxfUnvV/P8DkkXHvJ4vqT5tTkWfNNj81ZpeeY+vXZDEgupAgD8Au/jh2NmLd6mWUsydPvZ7XVet3in4wAA4BaULzhiWUaeHv1wlU7vGKt7z+vsdBwAANyG8gW321tQookzUhQXGarnRvRVYAAT7AEA/oPZzXCr8gqX7pyVquzCUr03cbAaNQhxOhIAAG7FyBfc6t+fr9dP6Xv199/1UI+W0U7HAQDA7ShfcJvPVu7UK99t1MhBibomKeH4/wEAAD6I8gW32LA7X/fNWaY+CTH6yyXdnI4DAIBjKF+od/uLy3Tr9BSFhwTq5VH9FBoU6HQkAAAcw4R71CuXy+q+Ocu0LadIb908SM2jw52OBACAoxj5Qr168Zt0fbF6tx6+qKsGtWvidBwAABxH+UK9+WbdHv33y/W6vG9LjRncxuk4AAB4BMoX6sWW7ELdNTNVXZtF6cnLe8oYFlIFAECifKEeFJWWa8KMFAUEGL06ur/CQ5hgDwDAQZQv1Clrrf74znKt352v56/rq4TGEU5HAgDAo1C+UKde/m6jPlmxU/cP66LTO8Y5HQcAAI9D+UKd+XbdHj29YJ0u6d1C489o53QcAAA8EuULdWJLdqHunJmqzvGR+ueVTLAHAOBoKF+otcKSco2fnqyAAKPXbkhSRAhr9wIAcDSUL9SKtVZ/nLtM6XsK9MJ1/ZhgDwDAcVC+UCsvfpOu+St26YHhXXRax1in4wAA4PEoXzhpX67erf98sV6/69NCt5zOBHsAAGqC8oWTkr4nX3fPTlOPFtH6x5W9mGAPAEANUb5wwvYdKNMt01IUFhygV0f3V1gwK9gDAFBTvC0NJ6TCZXXnzFRl5hbp7VtOUYuYcKcjAQDgVShfOCFPL1in79Zn6YnLe2hAm8ZOxwEAwOtw2RE19mHadr3y3UaNHJSo6we1djoOAABeifKFGlmWkac/zV2uQW0b67FLujsdBwAAr0X5wnHt2V+s8dOTFRcZqpeu76eQIL5tAAA4Wcz5wjEVl1Xolukpyi8u17sTB6tJw1CnIwEA4NUoXzgqa60eem+FlmXk6ZVR/dW1eZTTkQAA8HpcP8JRvfbDJr2Xul33ntdJw3o0czoOAAA+gfKFan29dree+nStLurZXHec08HpOAAA+AzKF46wfne+7pyZpm7No/T01dw6CACAukT5wmFyCks1buoShYcE6vUbkxQRwrRAAADqEuULvyotd2nC9BTt2V+i125IUvNobh0EAEBdY1gDkirf2fjIByu0eEuO/jeij/okxDgdCQAAn8TIFyRJk3/crDnJmbrjnA66rE9Lp+MAAOCzKF/Q12t368n5azS8RzPdM7ST03EAAPBplC8/t2bnft3xdqq6tYjSf67prYAA3tkIAEB9onz5sT37izVuyhJFhgVr8o0DeGcjAABuwG9bP3WgtEK3TEtWblGZ3plwquKjwpyOBACAX6B8+SGXy+q+d9K0fPs+vTqqv3q0jHY6EgAAfoPLjn7ov1+s1/wVu/Tg8C46vzv3bAQAwJ0oX35mbkqmXvgmXSMGJOiW09s5HQcAAL9D+fIjCzdm68H3lmtIhyZ6/Hc9uGcjAAAOoHz5ifQ9+ZowPUVtmjTQS9f3V3Agpx4AACfwG9gPZBeUaOyUJQoJCtQbYwYoOjzY6UgAAPgtypePKy6r0M1Tk5WVX6LJNyYpoXGE05EAAPBrLDXhw1wuq3tmp2lZZp5eGdVfvblZNgAAjmPky4f947O1+nTlLj18YVddwJISAAB4BMqXj5q6cIsmfb9JN57aWuNOa+t0HAAAUIXy5YMWrNqlxz5apfO7xevRS7qzpAQAAB6E8uVjUrbm6s6ZqeqTEKP/jeirwACKFwAAnoTy5UM2Zxfq5qlL1Dw6TK/fkKTwkECnIwEAgN+gfPmI7IISjXlzsYwxmjJ2oJo0DHU6EgAAqAblywcUlZZr3NRk7d5frNdvTFKb2AZORwIAAEdB+fJyZRUu3f7WUq3IzNNzI/qqX2IjpyMBAIBjYJFVL2at1YPvrdA367L05OU9dT5reQEA4PFqNfJljLnaGLPKGOMyxiQdY78txpgVxpg0Y0xybY6J//fvz9dpbkqm7h7aUSMHJTodBwAA1EBtR75WSrpC0qs12Pdsa212LY+HKtN+3qIXv9mo6wYm6q5zOzodBwAA1FCtype1do0kFvF0s09X7NRf5q3S0K7xevwyFlEFAMCbuGvCvZX0uTEmxRgz/lg7GmPGG2OSjTHJWVlZbornPRamZ+uuWWnql9hIz1/XV0GBvGcCAABvctyRL2PMl5Kqm8n9sLX2wxoeZ4i1docxpqmkL4wxa62131e3o7V2kqRJkpSUlGRr+PX9worMfbplWrLaxjbQGzcOYBFVAAC80HHLl7V2aG0PYq3dUfXvHmPM+5IGSqq2fKF6m7IKNObNxYqJCNHUmwYqOiLY6UgAAOAk1Ps1K2NMA2NM5MHPJZ2vyon6qKFd+4o1evJiSdKMmwepWXSYw4kAAMDJqu1SE5cbYzIlnSrpE2PMgqrnWxhj5lftFi/pR2PMMkmLJX1irf2sNsf1J3lFpRo9+RftO1CmqTcNVFtWrwcAwKvV9t2O70t6v5rnd0i6sOrzTZJ61+Y4/qqwpFxjpyzR1r1FmnLTAPVoGe10JAAAUEuscO+hissqdMu0ZC3P3KcXR/bT4PaxTkcCAAB1gHUKPFBZhUu/fztVCzfu1b+u7KVhPbhtEAAAvoLy5WFcLqs/vrNMX67Zrb9d1l1X9m/ldCQAAFCHKF8exFqrR+et1AdpO/THCzrrhlPbOB0JAADUMcqXh7DW6l8L1mnGom269cx2uu2s9k5HAgAA9YDy5SGe/zpdL3+7USMHJeqBYV24XyMAAD6K8uUBXv1uo/77xXpd2a+V/n5ZD4oXAAA+jPLlsCk/bdZTn67Vxb2a619X9VJAAMULAABfRvly0MzF2/TYR6t1frd4PXNtHwVSvAAA8HmUL4e8tzRTD72/Qmd1jtPzI/sqOJBTAQCAP+A3vgM+TNuuP7yzTKe2a6JXRvVXaFCg05EAAICbUL7c7KNlO3TP7DQNaNNYr9+YpLBgihcAAP6E8uVGHy/fobtnpympTWO9OXaAIkK4tSYAAP6G8uUmnyzfqbtmpal/YiO9OYbiBQCAv6J8ucH8FTt156xU9UuM0ZtjB6hBKMULAAB/RfmqZ58s36k7Zqaqb0KM3hw7kOIFAICfo3zVow9St+uOmUvVLzFGU24aqIYULwAA/B7lq568k5yhe+akaVDbJpoyluIFAAAq0Qjqwdu/bNND76/Q6R1jNWl0ksJDWE4CAABUonzVsWk/b9GjH67S2Z3j9PKo/qzjBQAADkP5qkOTvt+oJ+ev1Xnd4vXCyL6sXA8AAI5A+aoD1lo988V6Pfd1ui7q1VzPXtuHezUCAIBqUb5qyeWyevyT1Xrzpy26NilBT17RU4EBxulYAADAQ1G+aqHCZfXAu8v1Tkqmxp3WVo9c1FXGULwAAMDRUb5OUmm5S/fMTtMnK3bq7qEddde5HSleAADguChfJ6GwpFwTZqTohw3ZeuSirrr59HZORwIAAF6C8nWCcgpLNXbKEq3IzNO/ruylawYkOB0JAAB4EcrXCdied0CjJ/+i7bkH9OroJJ3XLd7pSAAAwMtQvmpow+58jZ68WIWl5Zo+bpAGtm3sdCQAAOCFKF81kLI1RzdNSVZoUIDm3HqqujaPcjoSAADwUpSv4/hs5U7dNStNzaPDNH3cICU0jnA6EgAA8GKUr2N448fNevyT1eqTEKPXb0hSk4ahTkcCAABejvJVDZfL6on5azT5x826oHu8nr22r8JDuE8jAACoPcrXbxSXVejeOWmav2KXxgxuoz9f3I3bBQEAgDpD+TpEdkGJxk9LVmpGHounAgCAekH5qrJuV77GTV2i7IISvTSyn4b3bO50JAAA4IMoX5K+WbdHd7ydqoiQQM259VT1ahXjdCQAAOCj/Lp8WWs1deEW/e3j1erSLEqTxySpeXS407EAAIAP89vyVeGyemzeKk1ftFXndYvXs9f2UYNQv/2/AwAAuInfto0AI+07UKZbz2yn+y/oogDe0QgAANzAb8uXMUbPXNuHZSQAAIBbBTgdwEkULwAA4G5+Xb4AAADcjfIFAADgRpQvAAAAN6J8AQAAuBHlCwAAwI0oXwAAAG5E+QIAAHAjyhcAAIAbUb4AAADciPIFAADgRpQvAAAAN6J8AQAAuBHlCwAAwI0oXwAAAG5E+QIAAHAjyhcAAIAbUb4AAADcyFhrnc5wVMaYLElb6/kwsZKy6/kYODGcE8/EefE8nBPPxHnxPO46J62ttXHH28mjy5c7GGOSrbVJTufA/+OceCbOi+fhnHgmzovn8bRzwmVHAAAAN6J8AQAAuBHlS5rkdAAcgXPimTgvnodz4pk4L57Ho86J38/5AgAAcCdGvgAAANyI8gUAAOBGflO+jDHDjDHrjDHpxpgHqtkeaoyZXbX9F2NMG/en9C81OCf3GmNWG2OWG2O+Msa0diKnvzneeTlkv6uMMdYY4zFv3/ZVNTknxphrqn5eVhlj3nZ3Rn9Ug9ewRGPMN8aY1KrXsQudyOkvjDFvGGP2GGNWHmW7McY8V3W+lhtj+rk740F+Ub6MMYGSXpQ0XFI3SdcZY7r9ZrdxknKttR0kPSPpn+5N6V9qeE5SJSVZa3tJmivpX+5N6X9qeF5kjImUdKekX9yb0P/U5JwYYzpKelDSEGttd0l3uz2on6nhz8ojkuZYa/tKGiHpJfem9DtTJA07xvbhkjpWfYyX9LIbMlXLL8qXpIGS0q21m6y1pZJmSbrsN/tcJmlq1edzJZ1rjDFuzOhvjntOrLXfWGuLqh4uktTKzRn9UU1+ViTpcVWW4WJ3hvNTNTknt0h60VqbK0nW2j1uzuiPanJerKSoqs+jJe1wYz6/Y639XlLOMXa5TNI0W2mRpBhjTHP3pDucv5SvlpIyDnmcWfVctftYa8sl7ZPUxC3p/FNNzsmhxkn6tF4TQarBeTHG9JWUYK392J3B/FhNflY6SepkjPnJGLPIGHOsv/5RN2pyXh6TNMoYkylpvqQ73BMNR3Giv3fqTZATB3VAdSNYv11joyb7oO7U+P9vY8woSUmSzqzXRJCOc16MMQGqvCw/xl2BUKOflSBVXko5S5UjxD8YY3pYa/PqOZs/q8l5uU7SFGvtf4wxp0qaXnVeXPUfD9XwmN/z/jLylSkp4ZDHrXTk8O+v+xhjglQ5RHys4UvUTk3OiYwxQyU9LOlSa22Jm7L5s+Odl0hJPSR9a4zZIukUSfOYdF+vavr69aG1tsxau1nSOlWWMdSfmpyXcZLmSJK19mdJYaq8wTOcUaPfO+7gL+VriaSOxpi2xpgQVU58nPebfeZJurHq86skfW1ZgbY+HfecVF3eelWVxYs5LO5xzPNird1nrY211rax1rZR5Vy8S621yc7E9Qs1ef36QNLZkmSMiVXlZchNbk3pf2pyXrZJOleSjDFdVVm+styaEoeaJ+mGqnc9niJpn7V2pxNB/OKyo7W23Bjze0kLJAVKesNau8oY8zdJydbaeZImq3JIOF2VI14jnEvs+2p4Tp6W1FDSO1Xvfdhmrb3UsdB+oIbnBW5Uw3OyQNL5xpjVkiok/dFau9e51L6vhuflPkmvGWPuUeXlrTH8UV9/jDEzVXnpPbZqnt1fJAVLkrX2FVXOu7tQUrqkIkljnUnK7YUAAADcyl8uOwIAAHgEyhcAAIAbUb4AAADciPIFAADgRpQvAAAAN6J8AQAAuBHlCwAAwI3+D006ixA64ywPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = lambda x: np.sin(4 * (x - 1/4)) + x + x**20 - 1\n",
    "x = np.linspace(0, 1, 100)\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "plt.plot(x, f(x))\n",
    "plt.axhline(ls='--', c='k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The unique root is approximately 0.408\n",
    "\n",
    "Let’s consider some numerical techniques for finding roots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bisection\n",
    "\n",
    "\n",
    "<a id='index-6'></a>\n",
    "One of the most common algorithms for numerical root finding is *bisection*\n",
    "\n",
    "To understand the idea, recall the well known game where\n",
    "\n",
    "- Player A thinks of a secret number between 1 and 100  \n",
    "- Player B asks if it’s less than 50  \n",
    "  \n",
    "  - If yes, B asks if it’s less than 25  \n",
    "  - If no, B asks if it’s less than 75  \n",
    "  \n",
    "\n",
    "\n",
    "And so on\n",
    "\n",
    "This is bisection\n",
    "\n",
    "Here’s a fairly simplistic implementation of the algorithm in Python\n",
    "\n",
    "It works for all sufficiently well behaved increasing continuous functions with $ f(a) < 0 < f(b) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def bisect(f, a, b, tol=10e-5):\n",
    "    \"\"\"\n",
    "    Implements the bisection root finding algorithm, assuming that f is a\n",
    "    real-valued function on [a, b] satisfying f(a) < 0 < f(b).\n",
    "    \"\"\"\n",
    "    lower, upper = a, b\n",
    "\n",
    "    while upper - lower > tol:\n",
    "        middle = 0.5 * (upper + lower)\n",
    "        # === if root is between lower and middle === #\n",
    "        if f(middle) > 0:\n",
    "            lower, upper = lower, middle\n",
    "        # === if root is between middle and upper  === #\n",
    "        else:              \n",
    "            lower, upper = middle, upper\n",
    "\n",
    "    return 0.5 * (upper + lower)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact SciPy provides it’s own bisection function, which we now test using the function $ f $ defined in [(2)](#equation-root-f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4082935042806639"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import bisect\n",
    "\n",
    "bisect(f, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Newton-Raphson Method\n",
    "\n",
    "\n",
    "<a id='index-8'></a>\n",
    "Another very common root-finding algorithm is the [Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method)\n",
    "\n",
    "In SciPy this algorithm is implemented by `scipy.optimize.newton`\n",
    "\n",
    "Unlike bisection, the Newton-Raphson method uses local slope information\n",
    "\n",
    "This is a double-edged sword:\n",
    "\n",
    "- When the function is well-behaved, the Newton-Raphson method is faster than bisection  \n",
    "- When the function is less well-behaved, the Newton-Raphson might fail  \n",
    "\n",
    "\n",
    "Let’s investigate this using the same function $ f $, first looking at potential instability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4082935042793568"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import newton\n",
    "\n",
    "newton(f, 0.2)   # Start the search at initial condition x = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7001700000000282"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newton(f, 0.7)   # Start the search at x = 0.7 instead"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second initial condition leads to failure of convergence\n",
    "\n",
    "On the other hand, using IPython’s `timeit` magic, we see that `newton` can be much faster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "70.6 µs ± 333 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit bisect(f, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18.2 µs ± 680 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit newton(f, 0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hybrid Methods\n",
    "\n",
    "So far we have seen that the Newton-Raphson method is fast but not robust\n",
    "\n",
    "This bisection algorithm is robust but relatively slow\n",
    "\n",
    "This illustrates a general principle\n",
    "\n",
    "- If you have specific knowledge about your function, you might be able to exploit it to generate efficiency  \n",
    "- If not, then the algorithm choice involves a trade-off between speed of convergence and robustness  \n",
    "\n",
    "\n",
    "In practice, most default algorithms for root finding, optimization and fixed points use *hybrid* methods\n",
    "\n",
    "These methods typically combine a fast method with a robust method in the following manner:\n",
    "\n",
    "1. Attempt to use a fast method  \n",
    "1. Check diagnostics  \n",
    "1. If diagnostics are bad, then switch to a more robust algorithm  \n",
    "\n",
    "\n",
    "In `scipy.optimize`, the function `brentq` is such a hybrid method, and a good default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.40829350427936706"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "brentq(f, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18.1 µs ± 57.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit brentq(f, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the correct solution is found and the speed is almost the same as `newton`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Root Finding\n",
    "\n",
    "\n",
    "<a id='index-9'></a>\n",
    "Use `scipy.optimize.fsolve`, a wrapper for a hybrid method in MINPACK\n",
    "\n",
    "See the [documentation](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html) for details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed Points\n",
    "\n",
    "\n",
    "<a id='index-10'></a>\n",
    "SciPy has a function for finding (scalar) fixed points too"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(1.)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import fixed_point\n",
    "\n",
    "fixed_point(lambda x: x**2, 10.0)  # 10.0 is an initial guess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don’t get good results, you can always switch back to the `brentq` root finder, since\n",
    "the fixed point of a function $ f $ is the root of $ g(x) := x - f(x) $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization\n",
    "\n",
    "\n",
    "<a id='index-12'></a>\n",
    "Most numerical packages provide only functions for *minimization*\n",
    "\n",
    "Maximization can be performed by recalling that the maximizer of a function $ f $ on domain $ D $ is\n",
    "the minimizer of $ -f $ on $ D $\n",
    "\n",
    "Minimization is closely related to root finding: For smooth functions, interior optima correspond to roots of the first derivative\n",
    "\n",
    "The speed/robustness trade-off described above is present with numerical optimization too\n",
    "\n",
    "Unless you have some prior information you can exploit, it’s usually best to use hybrid methods\n",
    "\n",
    "For constrained, univariate (i.e., scalar) minimization, a good hybrid option is `fminbound`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import fminbound\n",
    "\n",
    "fminbound(lambda x: x**2, -1, 2)  # Search in [-1, 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Optimization\n",
    "\n",
    "\n",
    "<a id='index-13'></a>\n",
    "Multivariate local optimizers include `minimize`, `fmin`, `fmin_powell`, `fmin_cg`, `fmin_bfgs`, and `fmin_ncg`\n",
    "\n",
    "Constrained multivariate local optimizers include `fmin_l_bfgs_b`, `fmin_tnc`, `fmin_cobyla`\n",
    "\n",
    "See the [documentation](http://docs.scipy.org/doc/scipy/reference/optimize.html) for details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration\n",
    "\n",
    "\n",
    "<a id='index-15'></a>\n",
    "Most numerical integration methods work by computing the integral of an approximating polynomial\n",
    "\n",
    "The resulting error depends on how well the polynomial fits the integrand, which in turn depends on how “regular” the integrand is\n",
    "\n",
    "In SciPy, the relevant module for numerical integration is `scipy.integrate`\n",
    "\n",
    "A good default for univariate integration is `quad`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.33333333333333337"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "integral, error = quad(lambda x: x**2, 0, 1)\n",
    "integral"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact `quad` is an interface to a very standard numerical integration routine in the Fortran library QUADPACK\n",
    "\n",
    "It uses [Clenshaw-Curtis quadrature](https://en.wikipedia.org/wiki/Clenshaw-Curtis_quadrature),  based on expansion in terms of Chebychev polynomials\n",
    "\n",
    "There are other options for univariate integration—a useful one is `fixed_quad`, which is fast and hence works well inside `for` loops\n",
    "\n",
    "There are also functions for multivariate integration\n",
    "\n",
    "See the [documentation](http://docs.scipy.org/doc/scipy/reference/integrate.html) for more details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Algebra\n",
    "\n",
    "\n",
    "<a id='index-17'></a>\n",
    "We saw that NumPy provides a module for linear algebra called `linalg`\n",
    "\n",
    "SciPy also provides a module for linear algebra with the same name\n",
    "\n",
    "The latter is not an exact superset of the former, but overall it has more functionality\n",
    "\n",
    "We leave you to investigate the [set of available routines](http://docs.scipy.org/doc/scipy/reference/linalg.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<a id='sp-ex1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Previously we discussed the concept of [recursive function calls](https://lectures.quantecon.org/py/python_advanced_features.html#recursive-functions)\n",
    "\n",
    "Write a recursive implementation of the bisection function described above, which we repeat here for convenience"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def bisect(f, a, b, tol=10e-5):\n",
    "    \"\"\"\n",
    "    Implements the bisection root finding algorithm, assuming that f is a\n",
    "    real-valued function on [a, b] satisfying f(a) < 0 < f(b).\n",
    "    \"\"\"\n",
    "    lower, upper = a, b\n",
    "\n",
    "    while upper - lower > tol:\n",
    "        middle = 0.5 * (upper + lower)\n",
    "        # === if root is between lower and middle === #\n",
    "        if f(middle) > 0:\n",
    "            lower, upper = lower, middle\n",
    "        # === if root is between middle and upper  === #\n",
    "        else:              \n",
    "            lower, upper = middle, upper\n",
    "\n",
    "    return 0.5 * (upper + lower)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test it on the function `f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1` discussed above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Here’s a reasonable solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def bisect(f, a, b, tol=10e-5):\n",
    "    \"\"\"\n",
    "    Implements the bisection root finding algorithm, assuming that f is a\n",
    "    real-valued function on [a, b] satisfying f(a) < 0 < f(b).\n",
    "    \"\"\"\n",
    "    lower, upper = a, b\n",
    "    if upper - lower < tol:\n",
    "        return 0.5 * (upper + lower)\n",
    "    else:\n",
    "        middle = 0.5 * (upper + lower)\n",
    "        print(f'Current mid point = {middle}')\n",
    "        if f(middle) > 0:   # Implies root is between lower and middle\n",
    "            bisect(f, lower, middle)\n",
    "        else:               # Implies root is between middle and upper\n",
    "            bisect(f, middle, upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test it as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current mid point = 0.5\n",
      "Current mid point = 0.25\n",
      "Current mid point = 0.375\n",
      "Current mid point = 0.4375\n",
      "Current mid point = 0.40625\n",
      "Current mid point = 0.421875\n",
      "Current mid point = 0.4140625\n",
      "Current mid point = 0.41015625\n",
      "Current mid point = 0.408203125\n",
      "Current mid point = 0.4091796875\n",
      "Current mid point = 0.40869140625\n",
      "Current mid point = 0.408447265625\n",
      "Current mid point = 0.4083251953125\n",
      "Current mid point = 0.40826416015625\n"
     ]
    }
   ],
   "source": [
    "f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1\n",
    "bisect(f, 0, 1)"
   ]
  }
 ],
 "metadata": {
  "filename": "scipy.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  },
  "title": "SciPy"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}