sampling from arbitrary distributions

sampling.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:95d343db10aa274dd70b7495d62a543e9cdce9db99587d635bb8d1d4562e70c1"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "\n",
      "import numpy as np\n",
      "from numpy.random import uniform\n",
      "from scipy.integrate import romberg\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "2 methods for sampling an arbitrary PDF"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Rejection method"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def sample_rejection(pdf, xmin, xmax, ymax, n=1000):\n",
      "    \"\"\"\n",
      "    Sample a PDF using the rejection method\n",
      "    http://en.wikipedia.org/wiki/Rejection_sampling    \n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    pdf : function\n",
      "       The probability distribution function to sample from.\n",
      "    xmin, xmax : float\n",
      "       The domain of the PDF, or a finite approximation of it\n",
      "    ymax : An upper bound on pdf(x)\n",
      "    \n",
      "    Notes\n",
      "    -----\n",
      "    Generates random points in the interval [xmin, xmax], [0, ymax], and accepts\n",
      "    the random x points if the 2D point lies inside the PDF. \n",
      "    \n",
      "    This function effectively truncates the sample to [xmin, xmax]\n",
      "    \n",
      "    The function becomes slower as xmin decreases or xmax/ymax increase\n",
      "    \"\"\"\n",
      "    result = []\n",
      "    while len(result) < n:\n",
      "        x = uniform(xmin, xmax, size=1000)\n",
      "        y = uniform(0, ymax, size=1000)\n",
      "        result.extend(x[np.where(y < pdf(x))])\n",
      "    return np.array(result[:n])\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "CDF transform sampling"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def sample_cdf(cdf, n=1000):\n",
      "    \"\"\"\n",
      "    Sample a PDF using the inverse CDF transform sampling\n",
      "    see http://en.wikipedia.org/wiki/Inverse_transform_sampling\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    cdf : function\n",
      "       The cumulative distribution function to sample from\n",
      "       \n",
      "    n : int\n",
      "        number of samples to draw\n",
      "    \"\"\"\n",
      "    return np.array([invert(cdf, u) for u in uniform(0, 1, size=n)])\n",
      "\n",
      "def invert(func, y, eps=1e-3):\n",
      "    \"\"\"\n",
      "    Find x such that func(x) ~ y\n",
      "    \n",
      "    func(x) assumed to be monotonic between [0, 1]\n",
      "    \"\"\"\n",
      "    assert y >= 0 and y <= 1, \"Invalid Y value\"\n",
      "    \n",
      "    lo = 0.0\n",
      "    hi = 1.0\n",
      "    while func(hi) < y:\n",
      "        hi *= 2\n",
      "        \n",
      "    # f(lo) < y < f(hi)\n",
      "    # narrow this range with binary search\n",
      "    mid = (lo + hi) / 2.\n",
      "    while (hi - lo) > eps:\n",
      "        z = func(mid)\n",
      "        if z > y:\n",
      "            hi = mid\n",
      "        else:\n",
      "            lo = mid\n",
      "        mid = (lo + hi) / 2.\n",
      "        \n",
      "    return mid\n",
      "            \n",
      "def numerical_cdf(f):\n",
      "    \"\"\"\n",
      "    Given a PDF, return a function that numerically evaluates the CDF\n",
      "    \n",
      "    Caches previous results for efficiency\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    f : function\n",
      "       The PDF\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    cdf : function\n",
      "       The CDF\n",
      "    \"\"\"\n",
      "    cache = {0:0}    \n",
      "    def wrapper(x):\n",
      "        # lookup a previous result cdf(best) for best < x,\n",
      "        best = max(c for c in cache if c < x)\n",
      "        # compute cdf(best) + integral(pdf, best, x)\n",
      "        result = cache[best] + romberg(f, best, x, divmax=20)\n",
      "        # store for future use\n",
      "        cache[x] = result\n",
      "        return result\n",
      "    \n",
      "    return wrapper"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def viz(pdf, samples):\n",
      "    plt.hist(samples, normed=True, facecolor='#565656', edgecolor='w', bins=40, log=True)\n",
      "    x = np.linspace(0, max(samples), 5000)\n",
      "    plt.plot(x, pdf(x), 'r-', lw=2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example with a lognormal"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Rejection method (samples will be truncated to x <= 200)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import lognorm\n",
      "\n",
      "pdf = lognorm(1.2).pdf\n",
      "samples = sample_rejection(pdf, 0, 200.0, 1.0)\n",
      "viz(pdf, samples)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Users/beaumont/anaconda/lib/python2.7/site-packages/scipy/stats/_continuous_distns.py:2689: RuntimeWarning: divide by zero encountered in log\n",
        "  return -log(x)**2 / (2*s**2) + np.where(x == 0, 0, -log(s*x*sqrt(2*pi)))\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEDCAYAAAA2k7/eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VeW1x/HvSkJCwowkgZBIEBIhIoIoKIIERcU6oNfW\niK21WMcKWm71OheoFq22DjhU64AUvUrrVHtlKCopKINYkXlQMEIIEIYwyRze+8cONFKGJOcke59z\nfp/nOU/O2TnD2tmw8mbtd6/XnHOIiEj0i/M7ABERqRtK+CIiMUIJX0QkRijhi4jECCV8EZEYoYQv\nIhIjlPBFRGKEEr6ISIxIqM03N7MGwHPAbqDQOfe/tfl5IiJyZLU9wv8v4C/OuRuBS2v5s0RE5Ciq\nnfDN7BUzW2dm8w/Z3t/MlpjZV2Z2V8Xm1sCqivvlIcYqIiIhqMkIfzTQv/IGM4sHnqnYngcMNLOO\nQDGQFcJniYhImFQ7CTvnpgFlh2zuDnztnCtyzu0F3gQGAO8AV5jZc8D7oQYrIiI1F66TtpVLN+CN\n7Hs453YA14XpM0REJAThSvg17rFsZurPLCJSA845q87zw1VXX82/a/VU3C+u6oudc1F7GzZsmO8x\naP+0f7G4f9G8b87VbJwcroT/OZBjZtlmlggUUI2a/fDhwyksLAxTKCIi0auwsJDhw4fX6LU1mZb5\nBjAdyDWzVWY2yDm3DxgMTAIWAeOcc4ur+p7Dhw8nPz+/uqGIiMSc/Pz8Gif8atfwnXMDj7B9AjCh\nRlFEsWj/Rab9i2zRvH/RvG81ZTWtBYUtADM3bNgw8vPzdYBERI6hsLCQwsJCRowYgavmSdtAJHy/\nYxARiTRmVu2EH5yrX3fvhrJDr+cSEZFwCUTCH3777RS2agWpqfDSS36HIyISWKHM0glGSWfwYHjm\nGW9DQgLMnQt5eb7GJSISZJFb0pkyxfvauTPs2wd33ulvPCIiUSgYI/z4eHAOVqyATp1g+3b4/HPo\n1s3X2EREgipiR/jDy8spbNkS2rSBm27yNj72mL9BiYgEUOTX8AH69IHCQli1Ck44Afbvh6+/hrZt\nfY1PRCSIInaED1DeogW7d+9md1oa5VdeCfv3s//xx9m7d6/foYmIRIXAjPAntm/P6FNPBaBNWRmP\nTp7MnqQkyouKSG7Z0tcYRUSCJmJH+MOBmfv2HXz8bbNmzE9LI3H3bhJefdWvsEREAicqavh/PuUU\nPjjxxIPbTy0p4a5PPmF/VhZxK1Z48/NFRASI4BE+wO5DEvqcVq0oS0sjbtUqeOcdn6ISEYkegUn4\nuw5J+M6MuX37eg8ef9yHiEREoktgEv7u+Pj/2La0e3dc8+YwaxbMmOFDVCIi0SMwCX/PYWr0+5KS\n2Hfddd4DjfJFREISiIQ/HFiydethv7fv5puhXj2vjl9UVJdhiYgETp2uaVsbhgOtUlMP/82MDCgo\n8K68HTWqLsMSEQmcUNa0DUTCh/+cpfM9Q4d6X196CY7wl4CIiBxdYBL+nrijhHLqqZCfD9u2wcsv\n11lMIiLRJDAJv/wICT85Odm7c2CUP2qU1zNfRESqJTCXr5bbkS8YKygowJzjiYYNaVVUxONnn82s\nrCwAxo0bV1chiohEtMCP8A9wZozPyQHg4mXLvAVTRESkygKR8IcDqzdsOObz/tm2LdsSE8nduJFO\npaW1HpeISNBExbTM1PT0Yz5vd0IC/5ebC8APFy7UKF9EYk5UTMs8Wg2/skk5OWxLTKTjhg2cpFG+\niEiVBSbhuyom/J316vFBxSj/Rxrli4hUWSAS/r64OKhiwgeYmJPD9opRPoWFtReYiEgUCUTC33+M\nGTqHqjzKZ9gwjfJFRKogEAn/WFMyD2dCRS2fadNg/PhaiEpEJLoEIuFXd4QP3ij/3Y4dvQd33w3l\n5WGOSkQkukRswgeY1L49tGkDCxbA2LFhjkpEJLoEIuHXpKQDsC8+Hh56yHvwwAOwc2cYoxIRiS6B\nSPgP79lDaU3n1F99NZxyChQXw9NPhzcwEZGAifgrbe9ISSEtLa1mL46LY2STJgB8d//9XD9gAAUF\nBRQUFIQxQhGRYIj4K21rWsM/YG56OnNatqTB3r0MnD8/TFGJiESXQCT8ql5le0RmjOnShX1xcfT9\n5htO2LQpPIGJiESR6Ej4wJrGjRmfk0McMGjOHEwXY4mIfE/UJHyAt/PyKKtfn9yNG+n97bdheU8R\nkWgRVQl/V716vN65MwA/njcPtmwJy/uKiESDYCT8ML7XtDZtWNKiBU137fKuwBURESAoCT9MI3wA\nzHixWzevA+fzz8Mnn4TvvUVEIlggEn64FTdpwnsdOngPbrgBdu/2NyARkQAIRMIP6wi/wnsdO8KJ\nJ8KSJTByZNjfX0Qk0tRqwjeztmb2kpn99WjPq42Evzc+Hl580Xvw8MNegzURkRhWqwnfOfeNc+76\n2vyMo+rdG266CfbuhUGDvK8iIjGqSgnfzF4xs3VmNv+Q7f3NbImZfWVmd9U0iNoY4R/06KNeC+XP\nP1dpR0RiWlVH+KOB/pU3mFk88EzF9jxgoJl1NLNrzOwJM8sIb6g11LgxvPqqd//BB2H2bF/DERHx\nS5USvnNuGlB2yObuwNfOuSLn3F7gTWCAc26sc26oc67EzJqb2fNAl6P9BbC/Nkf4APn5MHSotyrW\nNdeob76IxKSEEF7bGlhV6XEx0KPyE5xzm4Cbj/VGT23dyoKKk6ppaWk1b5V8NCNHwqRJsGgR3HMP\nPPlk+D9DRKSWFBYWUlhYGNJ7hJLww3aB7G1NmrC2U6dwvd33VO6L3zYjg4cWLybhqaegf3/vJiIS\nAfLz88nPzz/4eMSIEdV+j1AS/mogq9LjLLxRfrXV6knbSr5p3py/durk9cy/5hpuPvNMypKTD/vc\ncePG1UlMIiJ1JZRpmZ8DOWaWbWaJQAHwfk3eaNSWLTVf4rCa/tahA3PT02HDBobMnInt318nnysi\nEg61vsShmb0BTAdyzWyVmQ1yzu0DBgOTgEXAOOfc4poEMaRp09qp2x+GM+PZHj0gPZ2T1q/nikWL\n6uRzRUTCIZQlDqtU0nHODTzC9gnAhBp9so+21K8Pr7/O/n79uGLRIhanprIwPd3vsEREalUgeunU\nZUnnoHPP5Z28POKA22bNotmOHXX7+SIiNVDrJZ3aNrhZszor6VT2dl4eC9LSaLprF/89YwYJ5eV1\nHoOISHWEUtIJRML3y/64OJ484ww2pKSQu3Ejg+bM8TskEZFaE4iEP2rz5rov6VTYVr8+v+/Zkz1x\ncfRbsYJzly/3JQ4RkapQSSdE3zRvzounnQbAdXPmkLNhg2+xiIgcTcSXdOrqwqujmZqdzYT27UnY\nv59fTZ8Oq1Yd+0UiIhFECb+SsV26sDA1lWa7dsHFF8O2bX6HJCISNoFI+M9u2uRbDb+y8rg4/tCz\nJyWNGsG8eVBQAPv2+R2WiMhBEV/D/8Vxx/law6/su6QkHundG447DiZMgF/+ElzY+sSJiIQk4mv4\nQbOuYUN47z1ITIRnn4VRo/wOSUQkZIFI+EGp4X9Pr14werR3f+hQeOstf+MREQlRKO2Rw+a5jRsp\nLS0NTFnnoKuvhm++gfvvhx//GJo3h3PO8TsqEYlhoSyEEogR/s0tWgQv2R9w770wZAjs2QMDBsC/\n/uV3RCISwyK+hh/Iks4BZt5yiAMHwvbtcOGF8NVXfkclIlJtgUj4gRcXB6++CuefD+vXe19LSvyO\nSkSkWgKR8AM9wj8gMRHefhu6d4eiIujXDwJw7YCISFUFIuFHjIYN4YMPoFMnWLzYS/rquyMiESIQ\ns3Se37AhmLN0DqdFC/jwQ1bn5NB6/ny+ycnhwT59+C4p6eBTtAC6iNSWiJ+lc2NaWmQk+wPS0/lN\nfj4lDRvSdvNm7ps6lZQ9e/yOSkRiQMTP0olEm5OTeTA/n7UNGtCurIx7p04lee9ev8MSETmiQJR0\ngnrStqCg4LDbD5RsNqWk8GB+PsOmTCFn0ybu/+c/Gdm7d12GKCJSZcFI+H4HEIINDRrwm759eaCw\nkPabNjGssBBKSykYMuSwz1d9X0T8EoiSTlBH+FW1vkEDhvXty+pGjWizZQv06UOznTv9DktE5HuU\n8MOkLCWFEX378m2TJrBkCcOnTKHFd9/5HZaIyEGBSPgvrlsXiAVQQrWlfn1+k58P3brRcvt2RkyZ\nQqutW/0OS0SiSMQvgHJ9y5aRNS3zKLYnJcFHH7H0uONosWMHv/n4Y9pv3Oh3WCISJTQtM2iaNOG3\nffowp2VLGu/ZwwOFhXRZs8bvqEQkxgUi4UdDDf9QuxMSeKxXLwqzs6lfXs7/fPIJfb75xu+wRCSG\nBSPh+x1ALSmPi+OPp5/Oux06EO8cv5g9Gx55RGvkiogvApHwicIR/kFmvNm5M6O7dmU/wD33wODB\nsG+f35GJSIwJRsKPARNzcnjqzDO9NsvPPQeXXAJbtvgdlojEECX8OjQzKws+/tjruDlxIpx1ltdb\nX0SkDgQi4cdURfuss2DWLOjYERYuhB49YOZMv6MSkRgQiF46L69ZQ2nTplEzF/+YTjgBpk+HK6+E\nyZMhP99bQvGqq4BjN20TkdgV8f3wr8vIiJ1kf0DTpt7qWTffDLt3e4uk33cflJf7HZmIBFgoF14F\nYoQfa743gneOC7t04Zq5c4kfORK++IIG9evzXWKifwGKSFQKxAg/ppkxITeX3559NluTkmDiREZ+\n+CFZmzf7HZmIRBkl/IBYmJ7OPf36QdeutNy+nYc++ogzVq3yOywRiSJK+AGyoUED+PRTprZpQ/3y\ncobOmMGP584lbv9+v0MTkSighB80yck82707o7t0odyMS5cu5deFhTTbscPvyEQkwgUi4Udj87SQ\nmDExN5cH+/RhU/36dNywgd9Nngz/+IffkYlIBAtEwpfDW5yWxt3nn8+89HSa7N4N/fvD/ferD4+I\n1IgSfsBtqV+fkb17M65TJ6/J3G9/C/36gfrri0g1KeFHABcXxzt5efDhh9CyJfzzn9Cli9ePR0Sk\nipTwI0nfvvDll3DuuVBaChdeCLffDrt2+R2ZiESAWk/4ZjbAzP5kZm+a2XmHe45O2lZDejpMmuQt\npJKQAKNGwemnw/z5fkcmIgFX6wnfOfc359yNwM3A4buCSfXEx8Ndd3ldNnNzYcECOO00ePJJ0Jx9\nETmCKvfSMbNXgIuAUufcyZW29weeBOKBl5xzvzvCW9wPPBNCrHKobt3giy/gV7+CF16AoUOZ+8gj\nPNe9O5uTk7/3VHXaFJHqjPBHA/0rbzCzeLwk3h/IAwaaWUczu8bMnjCzDPP8DpjgnPsybJGLp0ED\neP55eO89aNGCU9at47FJk+i5cqXWzhWR76lywnfOTQPKDtncHfjaOVfknNsLvAkMcM6Ndc4Ndc6V\nAEOAc4EfmtlN4QpcDjFgAMybx5ctW9J4zx5unzmToTNm0FgndEWkQqjtkVsDlTt8FQM9Kj/BOTcK\nGHW0NxldUsKCnTsBSEtLi73e+OHSqhUP9+7NOStWcM3cuZxRXEze+vW8fOqpfkcmIiEKZeGTA0JN\n+GGpGfwsI4Md7dqF463EjI/btWNey5bcPHs2J5eWMnTGDG91rWefhdRUvyMUkRrIz88nPz//4OMR\nI0ZU+z1CnaWzGsiq9DgLb5QvPtvQoAEP9enDi6eeyq6EBPjrX+Gkk+Ctt1TbF4lRoSb8z4EcM8s2\ns0S8aZfvV/dNRpeUUFpaGmIo8h/M+LB9e+44/3zvoq316+FHP4LLL4di/V4WiUSFhYU1XuKwygnf\nzN4ApgO5ZrbKzAY55/YBg4FJwCJgnHNucXWDGNS6ter2tWh9w4ZeW4bnnoNGjeBvf4O8PHjmGa2h\nKxJhQlnTtjqzdAY65zKcc0nOuSzn3OiK7ROccyc659o75x6uURRS++Li4JZbYPFiuOwy2LYNhgyB\nXr10la5IjAhEL53Rq1erpFNXWreGd9+Fd96BjAzvat1TT4X77oOKmVIiElx1UtKpTT9TSafuXX45\nLFoEv/iFV9YZORI6d9YiKyIBF0pJJ9RpmRLJmjTxpmr++Mdw442wcCFccAEzMzP58ymnsLFBA+DY\nbRkKCo7cIkktHUSCIxAj/Fc1S8dfPXvCF1/weufO7EpI4IziYh6fOJHLFi8mQSd1RQJFJR0JXWIi\n73fowND+/ZmelUX98nIGzp/PY5Mmee2YRSQQ6mSWTm3SZUDBsSklhafOPJMH+/RhdaNGZGzf7q2l\ne8UVsHKl3+GJSAgCkfAleBakp3Pn+efzWufOXkfOd96BDh1g2DD47ju/wxORGghEwn9V0zIDqTw+\nnr936ABLlkBBgTdt8ze/8RZdGTNGi62I+CCUGn4gZulcm5nJC6rh17oaz6bJzIQ334TBg2HoUPj8\nc/jZz7wrdZ94IvyBisgRHWii5kfzNIklvXrBrFne6D4jw0v8vXszdPp0Urdv9zs6ETkGJXypnrg4\n+OlPYdkyr56fnMwZxcU8MXEiA+fNI2XPHr8jFJEjCETCH1NcrBp+pGnQAIYPh2XLmNqmDfX27+ey\nJUsYNX48Fy1dSj3N3xepFRE/D//azEzNw49UmZk826MH9557LotSU2m0Zw8/nTuXJyZM4OyiInXj\nFAmziJ+HL5Fv+XHHMSI/n4d79+bbJk1I3bGDWz/7DLp0gQ8+0KIrIgGghC/hY8aXrVpx13nn8Uz3\n7qxPSYEFC+Dii6FPH5gxw+8IRWJaIKZlStVFQqMyFxfHtOxsZmZl8dpZZ7H17rtpPG0a9OzJ7IwM\n/tKpEyubNq1SvEfa33HjxoX0s4iEn2M4He3nKLEjEAl/THExpSkpquNHmb3x8TB0KLdNncqlS5dy\n0bJlnF5SwuklJczIzPTaM+fl+R2mSEQpLCyksLCwRq8NRElHJ22j287ERMadfDJDfvADxufksCcu\njjOLi6FTJ68187JlfocoEjF00lYiwpbkZMZ07crtP/gB/2jXDhIS4H//Fzp2hEGDYMUKv0MUiWpK\n+FLnNqWk8HK3bvDVV3D99WAGr74KJ57oLcTy7bd+hygSlZTwxT9t2sCLL8LSpV5vnv37vcc5OXDD\nDbB8ud8RikQVJXzxX7t2MHo0LF7s1fTLy+GllyA3l8EzZ9J6yxa/IxSJCkr4Ehy5ufDaa1475uuu\ng7g4eq9cye8nTWLo9Olkl5X5HaFIRAtEwlcvHfmenBx4+WX4+msmtWvHvrg4zigu5neTJ/M/06aR\ns2GD3xGK+Eb98KXW+XLhTps2vNKtG+/k5XHx0qWct3w53dasoduaNXDuuXQqL2dBWpp30jcoMYvU\nslD64Qci4YsczebkZF7r0oW/dezID5Yto//XX5Py8cc8AHzTtCnvd+jAzMxM9scF4g9WkcDS/xCJ\nGNuSkhh38sncetFF8NBDbE5Kou3mzdw+cyZPTZhA/6++ImnfPr/DFAksjfAl4uxITIT77mPwnDmc\nXVTExUuXkrF9O4PmzOGHCxcyqX17WL/e7zBFAkcJXyLW3vh4PmrXjo/btuW0khIuXbqU3I0b+eGi\nRXD88fy8dWv+LzeXdY0a+R2qSCAo4UvEc3FxzM7MZHbr1py4YQOXLl3KaSUlnL98Of2WL2d269aM\nz81lSYsWRzzBKxILVMOX6GHG0tRUHuvVCxYu5OO2bdkfF0eP1asZMWUKD3/4obcK1+7dfkcq4gsl\nfIlOeXm8cPrp3HrRRbyVl8fWpCROKCvzVuFq04YrFi6k8a5dfkcpUqdU0pGotjk5mb926sR7HTty\n1sqV/GDZMtqsW8eV69Zx+eLFfHr88YzPyeHbZs38DlWk1gUi4WsBFKlte+PjKWzblsLsbMbdcguz\nf/ITupWUkF9URH5REQtTUxmfm8u/WrXyO1SRo9ICKCJVZQZ9+/L7Xr34ZcWCLDsTEjhp/Xru/PRT\nnh4/Hh5+GNTqQwJKC6CI1MC6hg0Z07Urt1x8MWO6dGFtw4ak7tgB994LmZle585PPwXn/A5VJCyU\n8CXm7UxMZHxuLr+88EJGnn02XHIJ7NvnrcbVqxd07Qp/+hNs3+53qCIhUcIXqeDMmNuyJbz/vrfc\n4j33QGoqzJ0LN90ErVvDbbd5fftFIpASvsjhZGfDyJGwahW8/jqcdRZs3QpPPw15eXDOOfCXv8Ce\nPX5HKlJlSvgiR5OUBFdfDZ98Al9+6a25m5ICU6ZAQYE36r/jDm+ZRpGAU8IXqapTToEXXoCSEhg1\nCk4+GTZsgD/8ATp0gLPPhrFjYedOvyMVOaxAzMOX2HWkRUogwAuVNGkCQ4bA4MHw2WfewutvvgnT\npnm3226Dn/zEW4i9c+c6CSmUn2MQj4FfMQXxZxFOGuGL1JQZ9OjhLbi+Zo03+j/tNNi8GZ55xvuL\n4MD3NcNHAkAJXyQcGjXy6vuzZ8OcOXDrrd5fAp995o30W7WC66/3zgVoXr/4pFYTvpl1MLM/mtlf\nzOzntflZIoHRpYs3wi8pgTFjvLn827d7C7P37u0t0v7QQ7Bypd+RSoyp1YTvnFvinLsFuAq4oDY/\nSyRwUlLgpz/16vqLF8Pdd0NGBixfDg884E397NcPXnsNvvvO72glBlQp4ZvZK2a2zszmH7K9v5kt\nMbOvzOyuI7z2EuAD4M3QwxWJUB06eD16Vq6EiRPhqqsgMRE++giuuUYlH6kTVR3hjwb6V95gZvHA\nMxXb84CBZtbRzK4xsyfMLAPAOfd359yFwLVhjFskMsXHwwUXwBtvwNq18PzzcMYZsG3b90s+Dz4I\n337rd7QSZaqU8J1z04CyQzZ3B752zhU55/bijeAHOOfGOueGOudKzKyPmT1lZi8AU8IbukiEa9rU\na9kwY8Z/lnx+/Wuv5NOnj9fHp+zQ/34i1RdKDb81sKrS4+KKbQc55/7pnLvdOXeTc+7JED5LJLod\nruRTvz5Mner9UmjZEi6/HN56C7RSl9RQKBdeha3QOKa4mAUV85TT0tLUG19i14GSzwUXeL173n3X\nO6n78cfw3nverXFj+OEPvfbNffp4r5GoF8rCJweEkvBXA1mVHmfhjfKr7drMTHZlZ4cQikgUatwY\nrr3Wu61Z413N+9pr8MUX8Mor3q11axg4kDZlZXzbtKl3MZhEpfz8fPLz8w8+HjFiRLXfI5SE/zmQ\nY2bZQAlQAAwM4f1E5EhatYKhQ73b4sVeB8/XX4eiIvj973kUWNW4MZ8cfzzTjz+e0oYN/Y5YAqiq\n0zLfAKYDuWa2yswGOef2AYOBScAiYJxzrkaNwscUF1OqJeVEqqZjR+/CrRUrvBW5brmFrYmJZG3d\nysAFC3h6/HhGTp7MxUuWcJzm90edwsLCGi9xWKURvnPusCN359wEYEKNPrmSazMzeUF1e5HqMYOe\nPaFnT24uLeWUdes4a+VKTispoV1ZGe3Kyrhm3jw46yz679jBrKwsypKT/Y5aQnSgtFPXJZ2wGVNc\nTGlKik7WitRQeXw8X2Rk8EVGBon79tF1zRp6rlpF1zVrSJo+nUHAtV9+yeLUVGZmZTEzM5Ot9ev7\nHbbUQCgnbwOR8DXCFwmfPQkJzMrKYlZWFkl79/LnH/2I2XfcQZe1azlp/XpOWr+eQXPmsDA1FV58\nkYa7d7M9KcnvsKWKIn6ELyK1Y3e9ejBwIL9/7z2S9+zh9JISzly1is5r13JyaSnceCMvmDE/PZ1Z\nmZl83ro125T8o5YSvkiM2JmYyNTsbKZmZ9Ng925OX72aW5o3xyZPpuvatXRdu5b9//oXi1JTmdW6\nNbMzM/0OWcIsEAlfNXyRuvVdUhKFJ5zALePGcdOAAZy+ejU9Vq+m07p1dCotpVNpKT+fM8dbvP2/\n/guuuMJr9SC+C6WGH4gFUK7NzFSyF/HJtvr1+bhdOx4++2xuHDCAZ7p3Z3ZGBnvi42H6dG+R9rZt\noVs3GDlSC7b7LD8/v3anZYpIbPguMZFp2dlMy872TvhedRW88w588IF3he8XX8B990Fenjfqv+IK\nb91eXeEbEQKR8FXSEQme3fXqwZVXerddu+Af/4C334b334dFi7zbgw/CCSfApZfCgAHe6l4JgUgr\nUUslHRGpXfXre0l9zBgoLYVJk7w1fNPSvCt+n3wS+vaF9HRvla+339bC7bUklJJOIBK+iESQevXg\n/PPhhRe8dXs/+QTuvBNyc2HTJhg71uvm2aIFXHSR189/zRq/oxaU8EUkFPHx3kyeRx/1TuYuXgyP\nPOK1fNizB8aP9/r5Z2RAjx7eSd+FC7WMo0+U8EUkfDp0gLvu8pq6rVkDL70El1zilYQ++8w74dup\nk7eM469+5S3wsm+f31HHjEAkfHXLFIlC6enw8597J3k3bPAWc/nZz7xSz/Ll8Pjj3gIuqaneCl9j\nx3rPk6MKpVtmIBK+TtqKRLkGDeCyy2D0aG/x9mnTvPn9J54ImzfDuHHeyd60NOjZk8sXLaJNWZlK\nP4ehefgiEjni473pm716wWOPeaP9Dz7wboWFMGMGVwFXLVjAxuRk5rRqxRetWrEgPZ3dmvIZEv30\nRMRf7drBbbd5t+3b4cMP+ei//5uua9dy3M6d9Fuxgn4rVrAnLo5FaWnw9NPe7J8TTvA78oijhC8i\nwdGwIVx2GX964w1wjuzNmzl1zRpOLSmh3aZNdFm79t+/HDp08BL/hRd6fy2oy+cxBSLh60pbEfkP\nZhQ1a0ZRs2a8k5dHo1276Lp2LbdmZ3sXfi1Z4t3+8AdISfEu/Orf37u1b+939LVGV9qKSNTbVr8+\nU7OzvRO869d79f477/R6+ezY4Z0DGDLEm/LZrh3ceiv8/e9Rd8WvrrQVkdhSr543pfPRR2HuXFi9\n2psBVFAAzZt77R6ee85rB9G8OZxzzr+fG8Mzf5TwRSTyZWR4c/zffNPr9TNzJgwfDmeeCeXlMGWK\nd0FYly7QujUMGuT9pbBxo9+R16lA1PBFRMImPt5r49CjBwwb5vX3+fBDr+4/caLX/+fVV72bGXTv\nDuedB+edR3x5OeXx8X7vQa1RwheR6Na8+b/bPDvn9fKZONG7TZsGs2Z5t4ce4uWEBBalpjIvPZ35\nLVuyulHmgqTuAAAFt0lEQVSjqOr1r4QvIrHDzOvl06mTd6Xvd995J38nT4bJk0letIhua9bQraK7\n58bkZOanpzMvPd0rFUX45JJAJHxNyxQRXzRo4M3lv+giAG6+5BJOLi2l89q1nFxaynE7d5JfVER+\nUZFX9pk71994CW1aZiAS/rWZmbygZC8iPitLSWFqdrY3/dM5jt+yhc7r1tF57VpOufRSv8MDvGmZ\n+fn5jBgxotqvDUTCFxEJHDNWNm3KyqZN+b8TT2Tcgw/6HVHINC1TRCRGKOGLiMQIJXwRkRihhC8i\nEiOU8EVEYoQSvohIjFDCFxGJEYGYh68rbUVEqkYLoIiIxAgtgCIiIsekhC8iEiOU8EVEYoQSvohI\njFDCFxGJEUr4IiIxQglfRCRGKOGLiMQIJXwRkRhR6wnfzBqY2Wwzu6i2P0tERI6sLkb4/wOMq4PP\nERGRo6hSwjezV8xsnZnNP2R7fzNbYmZfmdldh3ndecAiYH14wo08paWlfodQq2raxClSaP8iVzTv\nW01VdYQ/GuhfeYOZxQPPVGzPAwaaWUczu8bMnjCzDKAPcAZwNXCDmVn4Qo8MSviRTfsXuaJ532qq\nSu2RnXPTzCz7kM3dga+dc0UAZvYmMMA59wgwtuI591d871pgvXPOhSFmERGpgVD64bcGVlV6XAz0\nONwTnXNjQvgcEREJA6vqoLtihP9359zJFY+vAPo7526oePwToIdzbki1AjDTqF9EpAacc9Uqk4cy\nwl8NZFV6nIU3yq+W6gYsIiI1E8q0zM+BHDPLNrNEoAB4PzxhiYhIuFV1WuYbwHQg18xWmdkg59w+\nYDAwCW/q5Tjn3OLaC1VEREJRpYTvnBvonMtwziU557Kcc6Mrtk9wzp3onGvvnHu4uh9+rHn8kc7M\nisxsnpnNMbPP/I4nVIe7HsPMmpvZZDNbZmb/MLOmfsZYU0fYt+FmVlxx/OaYWf+jvUeQmVmWmU0x\ns4VmtsDMbqvYHi3H70j7FxXH0Mzqm9ksM/vSzBaZ2cMV26t1/Kp80jbcKubxLwX64Z0PmA0MjKa/\nEszsG6Cbc26T37GEg5n1BrYDf6508v5RYINz7tGKX9rNnHN3+xlnTRxh34YB25xzj/saXBiYWUug\npXPuSzNrCPwLuAwYRHQcvyPt35VEzzFMcc7tMLME4BPgDuBSqnH8/GyednAev3NuL/AmMMDHeGpL\n1JyUds5NA8oO2XwpcGDa7Ri8/2QR5wj7BlFy/Jxza51zX1bc3w4sxptaHS3H70j7B9FzDHdU3E0E\n4vH+vVbr+PmZ8A83j7/1EZ4bqRzwoZl9bmY3+B1MLUl3zq2ruL8OSPczmFowxMzmmtnLkVruOFTF\nFOuuwCyi8PhV2r+ZFZui4hiaWZyZfYl3nKY45xZSzePnZ8KPhfn3ZznnugIXArdWlA2iVsWV1NF0\nXP8ItAW6AGuAP/gbTugqyh1vA7c757ZV/l40HL+K/XsLb/+2E0XH0Dm33znXBcgEzjazvod8/5jH\nz8+EH5Z5/EHmnFtT8XU98C5eGSvarKuon2JmrYCoaR7knCt1FYCXiPDjZ2b18JL9WOfcexWbo+b4\nVdq/1w7sX7QdQwDn3BbgA6Ab1Tx+fib8qJ7Hb2YpZtao4n4D4Hxg/tFfFZHeB66tuH8t8N5RnhtR\nKv4DHXA5EXz8KhoXvgwscs49WelbUXH8jrR/0XIMzazFgXKUmSUD5wFzqObx822WDoCZXQg8iXcC\n4uWaTO0MKjNrizeqB++K5tcjff8qrsfoA7TAqxf+Gvgb8BfgeKAIuNI5t9mvGGvqMPs2DMjHKwU4\n4Bvgpkr10ohiZr2AqcA8/v1n/z3AZ0TH8Tvc/t0LDCQKjqGZnYx3Ujau4jbWOfeYmTWnGsfP14Qv\nIiJ1R2vaiojECCV8EZEYoYQvIhIjlPBFRGKEEr6ISIxQwhcRiRFK+CIiMUIJX0QkRvw/RGS4jwZl\npkMAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a2a1850>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Inverse CDF method using analytical cdf"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cdf = lognorm(1.2).cdf\n",
      "samples = sample_cdf(cdf, n=1000)\n",
      "viz(pdf, samples)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEDCAYAAAA2k7/eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FdX5x/HPkwABgbDHEAiEH1shKqAsitsFQSkqaLUi\ntuLPpWpfRa3VVuvPBWxt3Vptpe7UilrBpa4YQSoRLC4gLqyyBiJb2BdZZDm/P84NREogyc3NzM39\nvl+veeXO5GbmYUieOfeZM+eYcw4REan+UoIOQEREqoYSvohIklDCFxFJEkr4IiJJQglfRCRJKOGL\niCQJJXwRkSShhC8ikiRqxHPnZlYXeBTYBeQ75/4Zz+OJiEjp4t3C/xHwknPuamBQnI8lIiKHUe6E\nb2Z/N7M1ZjbroO0DzGy+mS00s1uim1sAhdHXe2OMVUREYlCRFv4zwICSG8wsFRgV3d4ZGGpmnYBv\ngOwYjiUiIpWk3EnYOTcV2HjQ5p7AIudcgXNuNzAWGAz8C7jAzB4F3ow1WBERqbjKumlbsnQDvmXf\nyzm3Hbiiko4hIiIxqKyEX+Exls1M4zOLiFSAc87K8/7Kqquv4ECtnujrb8r6w8650C133XVX4DEo\nJsWUjHEpprItFVFZCX8G0N7McsysFjCEctTsR4wYQX5+fiWFIiJSfeXn5zNixIgK/WxFumW+CEwD\nOphZoZld7pzbAwwHJgBzgXHOuXll3eeIESOIRCLlDUVEJOlEIpEKJ/xy1/Cdc0NL2Z4H5FUkiOKE\nH6akH6ZYiimmslFMZRfGuBTT4eXn51e4ImIVrQVVFjNzQccgIpJozAwX0E3b2Fx6KeTkwKRJQUci\nIlJtxXXwtLIa8fzzRIDI5ZfD4sVQq1bQIYmIhFLil3RKbhgzxrf4RUSkVIlb0gHIzfVfH3002DhE\nRKqpcJR0gMiZZxIpKICPP4aFC6F9+4CjEhEJn+pR0snLg7Fj4dln4c47YeTIQOMSEQmzxC7pNGly\noHY/Zgzs2xdsPCIi1UwoEv4IYPwXX7Cpa1f2ZWVBQQFb8/LYtGkTu3fvDjo8EZHQiGVohdCUdK4+\n91w216nD0K++4rz583m/TRue6NGDG264gd69ewcao4hI2CR0SWd3aioAU3JyADixsJBae/YEGJGI\nSPUSuoS/Ij2dRY0bc9SePXRfuTLgqEREqo/QdMtcsW4dGUcfDcAHrVvTbsMGTi8oYHugkYmIhEvC\nd8v8LjWVn15wwf5t9Xbt4om33iLFOT577TV6DB4cYIQiIuGTsDX8PdFyTrFtaWl81rw5Kc7RbOLE\ngKISEaleQpHwd9f478pS8c3bjLw80PDJIiIxC0XCP7iFD/B5ZiZb0tI4aulSmDkzgKhERKqXUCT8\nvSn/Hcbe1FQ+bNXKrzz5ZBVHJCJS/YQi4d+3axdFRUX/tX1S27b+xQsvwObNVRyViEj4JPyTtqsa\nN+aGfv0O+f0nFy6kweefwyOPwPDhVRydiEg4JWwvHWelx7z6Rz/yLx57TDdvRURiEI6Ef5jvbTjt\nNMjMhLlzYerUKotJRKS6CUfCP0wL39WoAVdd5Vf+9rcqikhEpPoJRcLfd4heOt9zzTVQowa8+ioU\nFFRJTCIi1U0oEv4RK/MtW8LFF8PevfDww1URkohItROOhH+Yks5+N9/svz79NGzcGN+ARESqoVAk\n/Ie2bTtkP/zv6dIF+veHb7+Fxx+vmsBEREImln74oUj416enk5GRceQ3/vrX/utf/wo7d8Y3KBGR\nEIpEIhVO+KEYD78sJZ0hQ4aAc9zXsCE5q1cz+uSTmdiuHQDjxo2Ld4giIgkvFC38MtXwAcx4tXNn\nAM6bN4+ae/fGMSoRkeolFAl/X1kTPjC9RQuWNWhAkx076LNkSRyjEhGpXkKR8Mvcwo++95XcXADO\nmz9frXwRkTJKuIQPvpVf0LAhTXbsoK9a+SIiZZKQCd+Z8Uq0ln/+vHmwbVs8whIRqVbimvDNrI2Z\nPW1mLx/ufeWp4Reb3qIFCxs3ptHOnfCnP1U4RhGRZBHXhO+cW+qcu+qI76tAwseM57t08a/vvx9W\nrSr/PkREkkiZEr6Z/d3M1pjZrIO2DzCz+Wa20MxuqWgQFUr4wPxmzZielQXbt8Ndd1X08CIiSaGs\nLfxngAElN5hZKjAqur0zMNTMOpnZpWb2kJlllTWIipR0ir1w3HGQmgqjR8OcORXej4hIdVemhO+c\nmwocPGJZT2CRc67AObcbGAsMds4955y70Tm30swam9njQNfDfQKoaAsfYFV6uh8+ed8+uPFGzYol\nIlKKWIZWaAEUllj/BuhV8g3OuQ3AtUfa0WMbNjB79mwAMjIyyjauTkkjR8KLL8J778Err8CPf1y+\nnxcRCbn8/Hzy8/Nj2kcsCb/SmtLXNG3KpmOOqfgOmjaFP/4Rrr3Wt/J/+EOoV6+ywhMRCVwkEiES\niexfHzlyZLn3EUvCXwFkl1jPxrfyy+3x9espKioqf8u+hIsnTeL3jRvTbsUK3uzenReKe/CgwdVE\npPqIpaUfS7fMGUB7M8sxs1rAEODNiuzoZ82axZTsAVxKCqOPP559wMAFC2i5eXNM+xMRCaNYhkcu\na7fMF4FpQAczKzSzy51ze4DhwARgLjDOOTevIkHEctO2pCWNGzOpbVtqOMfVM2Zg+/ZVyn5FRKqD\nMpV0nHNDS9meB+TFGsSTa9fGXNIp9uKxx9J9xQo6rl/PDxct4p0OHWLep4hIWARV0qk0V1WkZ04p\ntteqxVPduwNw8axZHL11a6XsV0QkDBJ+xqunKrGFDzAzK4uprVpx6vLlXDNjhu+jnxKKa5uISEwS\nvoV/5dFHV1qyL/aPbt3YlJZG7tq18NhjlbpvEZGgxP2mbbzF49nYbWlpjD7hBL/y61/DvArdTxYR\nqTZCkfBHFxVRVFRU6fv9tGVL8nNyYMcOGDoUdu2q9GOIiFSl/Pz8xG7hx6OkU+yZbt2gbVv48ku4\n7ba4HENEpKqopHMYO2vWhBde8CNq/vnPMHFiHI8mIhJeoUj4cderF9x9t389bJgmSxGRpBSKhD96\nzZq41PC/55ZbIBKBNWtgyBDYvTu+xxMRiYOEr+FfkZkZtxr+fqmpfgjl5s1h6lT47W/jezwRkThI\n+Bp+lcnMhJdegho1/MTnr7wSdEQiIlUmFAm/sgZPK5NTToEHHvCvr7gCvv666o4tIhKgUCT8Z1av\njn8Nv6QbboCLLoKtW2HQINh48OyNIiLhlPA1/MurooZfkpmf9Py442DBAj8lom7iikgCSPgafpWW\ndIrVqwdvvgkZGfDvf8Mvf1n1MYiIVKFQJPzAtG4Nr78OaWnw6KMwalTQEYmIxE0ohkeuCkOGDCn1\ne6d07cp1n3zia/tt2/pJ0EVEqplQtPDjObRCWXzYujXcfrsfN//CC+HTTwOOSESk8oUi4f9j1aqq\n7aVzKHffDZddBtu3w9ln+5u5IiIhk/C9dP43K6tqe+kcihk89RQMHAjr1sFZZ2nMHREJnYTvpRMa\nNWv6J3F79YKCAhgwADZvDjoqEZFKEYqEH3QN/3vq1oW334YOHeCrr3yLXxOhi0g1EI6EH0Q//MNp\n2tSPm5+dDdOmwbnn+tq+iEgCC0XCD6XWreH99yErCz74AM47D3buDDoqEZEKU8I/nHbt/FO4GRnw\n3nu+y+Z33wUdlYhIhYTiwat/rFpFUXp64D11Sns4K7tbNx6cMQPGj/fj7rz0kn86V0SkiuXn55Of\nn1+hnw1FC/+yMHTLPIzCBg18C79RIz/+zuDBqumLSCBi6ZYZihZ+QujWjd/06MHtH3xA+oQJzMnJ\n4b5TTmFXzZoAjBs3LuAARUQOLxQt/ESxrGFDRvTpw4batcldu5b/mzKFOqrpi0iCCEXCD1U//CNY\nkZ7OiD59WHvUUXRcv547PviA+uq9IyIJIBQJP9GsqV+fEX36sLpePdpu3Mjd778PS5cGHZaIyGEp\n4VfQurp1ubNPH5Y2bEjWtm1w0knw+edBhyUiUiol/BhsrlOHkX36MCsjA9asgdNPh0mTgg5LROSQ\n4p7wzWywmT1pZmPNrH+8j1fVdtSsyR9PPRWGDvVj7gwcCC++GHRYIiL/Je4J3zn3hnPuauBaoPRp\npxLY3tRUeP55uPFGPxn6JZfA738PLpFuR4tIdVfmhG9mfzezNWY266DtA8xsvpktNLNbDrOL24Hq\nO2lsSgr8+c/wpz/5sfXvuAMuvVTj74hIaJSnhf8MMKDkBjNLxSfxAUBnYKiZdTKzS83sITPLMu8+\nIM8590WlRR5Wv/qVnxi9bl144QXo29fX90VEAlbmhO+cmwpsPGhzT2CRc67AObcbGAsMds4955y7\n0Tm3ErgOOAO40MyuqazAQ23QIPjPf/zwyh99BD17wqxZR/45EZE4irWG3wIoLLH+TXTbfs65vzrn\nujvnfu6ceyLG4yWOLl38ZOi9esHy5dC7N/zrX0FHJSJJLNaxdCrlruSYlSuZvWMHABkZGaEeSK1c\nMjNh8mS48krfc+eCC+DWW/0N3dTUoKMTkQQSyyiZxWJN+CuA7BLr2fhWfrkMy8piV9u2MYYSvNKG\nV8Y5xj34INxyC9x7L8yY4S8ATZtWbYAikrAikQiRSGT/+siRI8u9j1hLOjOA9maWY2a18N0u3yzv\nTsasXElRUVGMoYSYGdx0k38oq1kz//WEE3ziFxEph/z8/AoPj1yebpkvAtOADmZWaGaXO+f2AMOB\nCcBcYJxzbl55gxgW8vHwK8uQxx7j5716sbBxY1i+nO969eKJ7t0ZctFFQYcmIgmiSsbDd84NLWV7\nHpBXoaNHjVm5kqL69ZMi6W846ihG9OnD/37+Of2XLOGazz7jmKIi2LIF0tODDk9EQi7hZ7xKlhZ+\nsT2pqTzdvTujevZkZ40anFxYCMcfrxKPiBxRLC38UCT8ZDU1J4db+/WjoGFDWLzYd918+GENySAi\ncRGKKQ6TqaRzsFXp6dx+xhk837w5jBoFN97IZ/ffz2M9e7L1oInSyzKNYqk9hcr48yISbirpJLjd\nqanwyCM82Ls322rW5IRVq3hgwgS6rloVdGgiEjIq6VQT01u25JYzz2Re06Y02rmT306dylUzZpC2\ne3fQoYlINaCSTsisq1uXkZEI5yxYwJDZs+m/ZAnHFhUxqmfPoEMTkRBQSaeacSkpvPWDH3Bbv34s\na9CAzG3buHvyZLj9dvjuu6DDE5EAqaRTTS1v2JDb+vXjjY4dfc+de+7xI2/OnBl0aCKSgFTSCbk9\nqan8s0sXZmZlMbKgAL78kr3du/NWx4680rkzu2sc+C9ULxyR6k8lnSQwv1kzmDWL8e3bY85x3vz5\n3P/ee3Rcuzbo0ESkCqmkkyzq1mVMt27cecYZfJOeTtbWrYyYPJnLZ86ktnryiMgRKOEnoIVNmnBL\n//682qkT+8wYsGgRD06YAG+9FXRoIhJiSvgJak9qKi8deyy39e/P4kaNaLZ9OwwaxM0ffkiTb78N\nOjwRCSHdtE1wyxo25PYzzuCsRYv438WL6bFyJceuWcOrubmM79CBvSm6potUJ7ppm+T2paSQ16ED\nzJvHtOxsau/dy0+++or7Jk6kk27qilQrumkrXosW/OWkk7jntNNYVa8e2Vu2MGLyZH7+6ac0iM4Z\nLCLJSwm/GvoqM5Nfn3UWL+fm8l1KCpGCAh7Oy4P774ddu4IOT0QCooRfTe1OTeWV3FxuPusspmdl\ncdSePX4S9dxceOMNjbkvkoRCkfCr/STmAVpTvz4PnnIK95x2GnTu7CdaOe886N8fZs8OOjwRKacq\nmcQ8nnTTNv6+ysyEL76Av/4VGjWCf/8bunSB4cNh3bqgwxORMtJNWymbmjXhuutg4UL4xS/ADP72\nN2jbFv74R9i+PegIRSSOQtEPX6pYkyZ+OsVrr4Wbb4YJE+C223zyv/tuuOwySE0t9cc1jaJIYlIL\nP5kdcwy8+y689x506wYrVsCVV/pSz9tv68auSDWjhC/Qrx/MmAEvvAA5OTBnDpx7LkQi8MknQUcn\nIpVECV+8lBS45BKYPx/+/Gdo3BimTIETT4Tzz4evvgo6QhGJkRK+fF9aGtx4o+++eeutUKcOvP66\nL/MMGQLz5gUdoYhUUCgSvvrhh1DDhr7nzpIlcP31UKsWvPQSHHMMv/jkE47eujXoCEWSkvrhS/xk\nZsJf/gKLFsE110BKCqctW8ZD777LNdOn01RDMYtUKfXDl/jLzobHH4cFC5ickwNA36VL+UteHlfN\nmEGzbduCjU9Ejkj98JNMzH3o27Th8Z49eaNTJy6cM4fey5fTf8kS+i5dytTWreHrrxly552HPUYi\n9ONPhBhFykstfKmQVfXr88iJJ3LTgAFMad0agEhBAXTqxA0ffUT2pk3BBigi/0UtfInJyvR0/tar\nFy/n5jJ4/nz6FRbSO7pMz8ritc6dWdy4cdBhighq4UslKapXj6e6d4fFi3mnfXu+S02lx8qV/GHS\nJH47ZQq5RUV6clckYGrhS+XKzubZbt14vVMnzv76a85cvJiuq1fTdfVqFjdqBC+9RMq+fezTXLsi\nVS6uf3Vm9gMze8zMXjKzK+N5LAmXzbVr888uXRh+9tm8lJvL5rQ02m7cCEOG8HBeHmctXEjanj1B\nhymSVOLawnfOzQd+bmYpwFhgdDyPJ+GzLS2NV3NzebNjR04vKOBnW7Zw9KJFXPH55/x4zhwmtGvH\nhHbt2FK7dtChilR7ZWrhm9nfzWyNmc06aPsAM5tvZgvN7JZSfvZcYDw+4UuS2l2jBpPatYP58/lT\n794saNKE+t99x4Vz5/K38eO5asYMmm/ZEnSYItVaWUs6zwADSm4ws1RgVHR7Z2ComXUys0vN7CEz\nywJwzr3lnPshcFklxi2JKjWVT1u25I6+fbmzTx+mZ2VRa+9e+i9ZwsPvvgsDB/ohm/ftCzpSkWqn\nTCUd59xUM8s5aHNPYJFzrgDAzMYCg51z9wLPRbedDvwIqA1MrpyQpVow4+tmzfi6WTOytmzh7AUL\nOG3ZMmrl5UFeHnTs6MfwGTYM6tULOlqRaiGWm7YtgMIS699Et+3nnPvAOXeDc+4a59zDMRxLqrGV\n6ek81b07Pz/nHD9gW8uW8PXXfhrGli3hppv8IG4iEpNYbtpWWqfqMStXMnvHDgAyMjI0kFqS2paW\n5odkvvlmeO01P2jbf/7jx+d/6CEYNAhuuMFPzGIWdLgiVSo/P5/8/PyY9hFLwl8BZJdYz8a38stt\nWFYWu9q2jSEUqVZq1IAf/9gvn33mE//YsfDGG37p3NnPx3vppX4YZ5EkEIlEiEQi+9dHjhxZ7n3E\nkvBnAO2jtf2VwBBgaEV2NGblSorq11fLPgQON2gYBDBw2AknwJgxcP/98MQTfsTOuXPh+uvZeeON\nTGvViolt27K0xPANZYlRg6N5iXAeEiHGylDWf2csLf2ydst8EZgGdDCzQjO73Dm3BxgOTADmAuOc\ncxWaDknj4csRZWbCXXfB8uXw8svMysig9t699F26lHsnTeKeSZOILFlCLT3MJdVcLOPhl7WXziFb\n7s65PCCvQkcuQS18KbOaNeHCC/n9yy/TfMsW+i1ZQqSggHYbNtBuwwYu/fJL+OUv/WQtnToFHa1I\npYt7Cz/e1MKXiliVns5zXbty7Tnn8GiPHixs3Jh6u3f7mn/nzv7m7nPPwfbtQYcqUmni3sKPN7Xw\nJRa7a9TggzZt+KBNG3I2buS+nBx44QX44AO/XHcdDB0KV17p7wmoh48kMLXwRaIKGjWCJ5+EVav8\nTd6ePWHzZn+zt0cP6NoVHnmEurt2BR2qSIVoTluRg6Wnw9VXwyefwFdf+bp+kyb+9fXX8/hbb3H9\nRx9x7Jo1mMbplyShko5Uf8ce6x/cuvde349/9GhqTJzIyYWFnFxYSFHdunyQk8OU1q0p0jAOEnIq\n6YiURVoaXHQRTJjAdWefzcu5uRQddRQZ337Lj+fM4ZF33mHE++/DU09x1HffBR2tyCEl/E1bkaq2\nrm5dXsnN5dXOnTlmzRpOW7aMnt98Q6d16+Dqq3kiJYUZLVowtXVrvszMZK9m6JJqQAlfkpozY1Zm\nJrMyMxl9/PH0XLGCX9SvT41//3v/ZOyb09L4T6tWTMnJYamGcpAEFoqErxq+hMHOmjWZkpPDL8aN\nY/i553LKsmWctmwZLbdsYeDChQxcuJDC9HS47z645BLIzj7yTkUqmWr4IpVs/VFH8UanTtx01lnc\n2q8f77Rvz5a0NLK3bPEjerZqBaeeCo8+CmvXBh2uJBF1yxSJFzOWNm7Ms926ce2553LfKaf4G791\n6sCHH/ox+5s3hwED4B//8H3+RUJKCV+kjPampDAzKwvGjYM1a+D55+Hss/2TuxMmwOWXw9FHw49+\nBC+/DNE5HkTCIhQJf8zKlRQVFQUdhkjZ1a8PP/kJvP02rF7tn+qNROC77/zkLRddBBkZfsz+d96B\n3buDjliqifz8/MQu6aiGLwmtSRP/VO/kyVBY6Gfo6tEDtm078Cng6KPhiit88lcff4mBavgiYdGi\nBdx4I3z6KSxcCL/7nR+5c+NGeOYZn/wzMuCyy+Ctt2DnzqAjliSihC8SL+3awe23w5w5frn7bj/M\nw+bNfhavQYN88v/JT3wZSDV/iTMlfJGq0Lkz3HGHH7xt/ny45x7o1g22boV//tPf6G3WDC6+GF55\nBb79NuiIpRpSwhepah07wm23wcyZsGiRH9Ste3ef5MeN85O3N2sGF1zgJ3DZsCHoiKWaCEXCVy8d\nSVpt28Itt8D06bB0KTz4IJx4oi/v/OtfMGyYL/uccQY88oif01eSmnrpiFQHOTlw003w0Uc+sT/y\niE/0AO+/D9dfD61b+1m7fvc7mDULNJZ/0lEvHZHqJjsbhg+HSZOgqMiXdi64AOrW9aWgO++E447z\nnxB+9SuYMgX27g06agk5JXyRsGvcGH76U38zd906/7DXVVf5Ov/SpX5yl9NPh8xM39f/9df9MwAi\nB1HCF0kktWv7vvxPPeXn7f3wQ7j5Zt/SX7fO9/U//3z/MNhZZ/my0JIlQUctIaGEL5KoUlPh5JPh\ngQf8Q16zZ/vunied5IdymDjR1/3btvXdQn/zG1/62bMn6MglIKEYD19EYmQGubl+ue02P2RzXh6M\nHw/vvgvz5vnlgQegYUM/uuc551Bv1y62paUFHb1UESV8keqoWTPfpXPYMN/a//BDn/zffhu+/hrG\njoWxY3nKjAVNmjCzeXNmNm9OYYMG/uIh1VIoEr5mvBKJo5o1oU8fvzz4oH/YK5r8973/Pj9Yt44f\nrFvHJbNmsb5OHb7MzOSLzEzYtMl/GpBQ0YxXIlJ27drBDTfAe+9x1eDB/Kl3bybn5LCpdm2a7NhB\n36VL+dVHH0HTpnDKKfD738OMGbBvX9CRC7H1ww9FC19EgrGjZk0+bdmST1u2xJyj1aZNdF29mi6r\nV5O7cSP85z9+ueMOfwE480xf/z/zTD/ksySUULTwRSR4zoxljRrxRqdO3N2nD6xf70fxvOYa/4Tv\nunV+oLdhw3yf/xNOgP/7P5g6VRO8JAi18EXk0NLT4bzz/OIcLFjge/y8+y7k5/snfmfOhD/8wb+3\nb1/o188vHTro5m8IKeGLyJGZ+VE+O3b09f8dO3yf/uILwPz5/gnf11/378/OPpD8zzhD5Z+QUElH\nRMqvTh3/JO9DD/n+/QUF8PTTfjz/pk39VI/PPOMnd8nM9OP+3HSTfzZAY/0HRi18EYld69Zw5ZV+\n2bfPT/QyaZJfpkzxI3vOmuXn+61Z0z8N3L+//wTQvTvUUCqqCnFv4ZtZXTObbmZnx/tYIhICKSnQ\ntasf4+fdd/18vpMn+yeAe/b0o3pOmeJ7/px0kh/357zzYNQoPxWkhnyOm6q4rP4GGFcFxxGRMEpL\ng0jEL/fc4y8A+fm+9f/ee34coDfe8Av4p4QjkQMPizmnG8CVpEwJ38z+DpwNFDnnji2xfQDwMJAK\nPO2cu++gn+sPzAVqV1rEIpLYGjXyI3qef75fX7bMJ//Jk/2yciW8/LJfgMdr12Zus2bMychgTkYG\nq+vV0wWggsrawn8GeAQYU7zBzFKBUUA/YAUw3czeBLoDxwMPAKcDdYHOwA4ze8c5fV4TkRJK1v+d\n8y3+4uSfn0+jNWs4ubCQkwsLAVhfp87+CwBLlkCbNroAlFGZEr5zbqqZ5Ry0uSewyDlXAGBmY4HB\nzrl7geei77k9+r3LgLVK9iJyWGa+D3+HDv6BL+f41cCB5BYVkVtUROe1a2myYwenLl/OqcuX+6Gf\nW7X6fgmodeug/xWhFUsNvwVQWGL9G6DXod7onHv2cDsas3Ils3fsACAjI0Pj6oiIZ8aK9HRWpKcz\nsV07zDlabtmy/wLQc8cOP//vmDF+AX8BOO00OPVU/7Vjx2rxCSCWQdOKxZLwK621Piwri11t21bW\n7kSkmnJmFDZoQGGDBrzbvj3jXnzRd/csLgFNmeIvAM8/7xfwN4FPOeXARaBLl4TsBhqJRIhEIvvX\nR44cWe59xPKvXgFkl1jPxrfyy03DI4tIhaSk+ATepQv88pf+GYDZs33inzrVf1292o8J9Npr/mfq\n14fevQ9cAHr08FNHJohYWvqxJPwZQPtobX8lMAQYWpEdDcvK4iklexGJVUqKf6r3uONg+HB/E3jx\n4u9fAJYsgQkT/AK+22jPngcuAL17+4tCSBW39OPWwjezF/E9bpqYWSFwp3PuGTMbDkzAd8sc7Zyb\nV+4IUAtfROLEzI//364dXHGF37ZihU/+xReA2bMPrIO/aHTr5stAvXv7eYNbtAju33CQuLfwnXOH\nbLk75/KAvAoduQS18EWkyrRo4cf8ufhiv75hg58CsvgC8NlnB5a//MW/p1WrA8m/d2//CSKg+wBx\nb+GLiFRbjRvDoEF+Adi2DT7+GKZN88tHH/kbwcuX+7mAAerWhV69DlwATjwxIaaDDEXCV0lHREKj\nXr0DQzuDH/tn7lw/89e0af7rkiXw/vt+AV86ys09cAHo3ds/IxCH7qBB3bStNCrpiEhopabCscf6\n5dpr/bZH7J1bAAAHO0lEQVTVq33Lv/gi8Nln/l7A7NnwxBP+PRkZB8pAJ50Exx/vh5WOkUo6IiJV\nKTPz++MB7dzpk37xBWDaNCgq+v6kMDVq+FFETzzxwPI//1OlD4WFIuGrpCMiCa12bd+SP/lkv+4c\nLFp0oAT08ce+9T9jhl9GjfLva9rU3ws48USOXbOGRY0asaNWrcMeSiUdEZEwMYP27f1y2WV+29at\nMH06fPKJvwB8/LH/FDB+PIwfz+3APmBFejoLmzTZv3xz0DMBKumIiIRd/fp+ove+ff26c35qyGjy\nX/j887TZtInsLVvI3rKFvkuXArCjRg0YONBfGGIs/yjhi4gEwcwP7dymDQwdyu2rV1Nz715yNm2i\n/fr1tF+/nnbr15OxfbufNKYSav2hSPiq4YuIwO7U1P2lnGINduzgyT/8Yf96LDX8uM9pWxbDsrKU\n7EVEDmFznTpwzDH71yORCCNGjKjQvkKR8EVEJP6U8EVEkoRq+CIiCUQ1fBGRJKEavoiIHJESvohI\nklDCFxFJEkr4IiJJQr10REQSiHrpiIgkCfXSERGRI1LCFxFJEkr4IiJJQglfRCRJKOGLiCQJdcsU\nEUkg6pYpIpIk1C1TRESOSAlfRCRJKOGLiCQJJXwRkSShhC8ikiSU8EVEkkRcE76ZRcxsqpk9Zman\nx/NYIiJyePFu4e8DtgJpwDdxPpaIiBxGmRK+mf3dzNaY2ayDtg8ws/lmttDMbjnEj051zg0EbgVG\nVkK8VaaoqCjoEP5LRZ+uiyedp7IJY0yg/7+yCmNMFVHWFv4zwICSG8wsFRgV3d4ZGGpmnczsUjN7\nyMyynHMu+vZN+FZ+wtAfQtnoPJVNGGMC/f+VVRhjqogyjaXjnJtqZjkHbe4JLHLOFQCY2VhgsHPu\nXuC56LbzgbOAhsAjlROyiIhURCyDp7UACkusfwP0KvkG59xrwGsxHENERCqJHai6HOGNvoX/lnPu\n2Oj6BcAA59zPous/BXo5564rVwBmZQtARES+xzln5Xl/LC38FUB2ifVsKtATp7wBi4hIxcTSLXMG\n0N7McsysFjAEeLNywhIRkcpW1m6ZLwLTgA5mVmhmlzvn9gDDgQnAXGCcc25e/EIVEZFYlCnhO+eG\nOueynHNpzrls59wz0e15zrmOzrl2zrk/lvfgZejHX+XMrMDMvjKzz83s04Bi+K/nHsyssZm9Z2YL\nzGyimTUMSVwjzOyb6Pn63MwGHG4fcYgp28wmm9kcM5ttZtdHtwd2vg4TU2Dnysxqm9knZvaFmc01\nsz9Gtwd5nkqLKdDfqWgMqdFjvxVdD/zvr5S4ynWuynzTtrJF+/F/DfTD3w+YDgwN+lOCmS0FTnDO\nbQgwhlOBbcCYEjfJ7wfWOefuj14cGznnbg1BXHcBW51zf67KWErElAlkOue+MLN6wGfAecDlBHS+\nDhPTRQR7ro5yzm03sxrAh8DNwCAC/L0qJaYzCPA8ReP6FXACUN85NygMf3+lxFWuv78gB0/b34/f\nObcbGAsMDjCekgK9keycmwpsPGjzIODZ6Otn8QmkSpUSFwR4vpxzq51zX0RfbwPm4bsMB3a+DhMT\nBHuutkdf1gJS8f+Xgf5elRITBHiezKwlMBB4ukQcgf/9lRKXUY5zFWTCP1Q//halvLcqOWCSmc0w\ns58FHUwJRzvn1kRfrwGODjKYg1xnZl+a2eigPurC/q7D3YBPCMn5KhHTx9FNgZ0rM0sxsy/w52Oy\nc24OAZ+nUmKCYH+nHgJ+jR8LrFgYfp8OFZejHOcqyIQf1v73JzvnugE/BH4RLWOESnTIirCcv8eA\nNkBXYBXwpyCCiJZOXgVucM5tLfm9oM5XNKZXojFtI+Bz5Zzb55zrCrQETjOzPgd9v8rP0yFiihDg\neTKzc4Ai59znlNJyDuI8HSaucp2rIBN+pfTjr2zOuVXRr2vxTwn3DDai/dZEa8OYWXMgFIOgOOeK\nXBT+o2aVny8zq4lP9s85516Pbg70fJWI6fnimMJwrqJxbAbG42vBofi9KhFT94DPU29gUPRe3otA\nXzN7juDP06HiGlPecxVkwg9dP34zO8rM6kdf1wXOBGYd/qeqzJvAZdHXlwGvH+a9VSb6y1/sfKr4\nfJmZAaOBuc65h0t8K7DzVVpMQZ4rM2ta/HHfzOoA/YHPCfY8HTKm4sQaVaXnyTl3W7QnYhvgYuB9\n59ylBPz3V0pcw8r7OxXLk7Yxcc7tMbPifvypwOige+jg63Kv+b9XagAvOOcmVnUQ5p97OB1oamaF\nwJ3AvcBLZnYlUIDv8RF0XHcBETPriv+IuxS4porDOhn4KfCVmX0e3fZbgj1fh4rpNvyIskGdq+bA\ns2aWgm/oPeec+3c0vqDOU2kxjQn4d6qk4tJN4H9/JRgH4rrfzLpQxnMVWLdMERGpWprTVkQkSSjh\ni4gkCSV8EZEkoYQvIpIklPBFRJKEEr6ISJJQwhcRSRJK+CIiSeL/AZ+wDwOGOY3YAAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10a494e90>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Inverse CDF method using numerical integration of PDF"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cdf = numerical_cdf(pdf)\n",
      "samples = sample_cdf(cdf, n=1000)\n",
      "viz(pdf, samples)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEDCAYAAAA2k7/eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX5x/HPSUgiqxDIwhIIyiJLEJFNEQiKCiLQqi+R\n1lZxX1CrrVVbKmBdUasitWhFVKqCdQNU5AdKFIyAKCA7yBIIS9iXABISzu+PEzBSAklmJvfO5Pt+\nveblzM1k7sNFnpw895znGGstIiIS+aK8DkBERMqHEr6ISAWhhC8iUkEo4YuIVBBK+CIiFYQSvohI\nBaGELyJSQSjhi4hUEJVC+eHGmKrAS8AhIMNa+3YozyciIsUL9Qj/CuBda+0tQL8Qn0tERE6i1Anf\nGPOaMSbHGLPouOO9jDHLjTGrjDEPFB6uD2wofF4QYKwiIhKAsozwxwK9ih4wxkQDowqPtwQGGmNa\nANlASgDnEhGRICl1ErbWzgR2HXe4I/CjtXadtfYwMB7oD3wAXGmMeQmYFGiwIiJSdsG6aVu0dANu\nZN/JWnsAuCFI5xARkQAEK+GXuceyMUb9mUVEysBaa0rz/mDV1Tfyc62ewufZJf1ma23YPoYOHep5\nDIrf+zgUf/g9wjl2a8s2Tg5Wwp8HNDXGpBpjYoEBlKJmP2zYMDIyMoIUiohI5MrIyGDYsGFl+t6y\nTMt8B8gEmhljNhhjBllr84HBwFRgKTDBWruspJ85bNgw0tPTSxuKiEiFk56eXuaEX+oavrV2YDHH\npwBTyhRFGAv3H1SK31uK3zvhHHtZmbLWgoIWgDF26NChpKenV8i/ABGR0sjIyCAjI4Phw4djS3nT\n1hcJ3+sYRETCjTGm1AnfP6tfd+6EAnVfEBEJFV8k/GFXXEFGnTrQurVL/CIickKBzNLxR0mneXNY\nscIduOUWePllT2MSEfG7spR0/JHwix6IiYHVqyElpbhvERGp8MK7hn/uuXD11XD4MLz6qtfRiIhE\nHF8k/GFARuXKcPvt7sCYMZCf72VIIiK+FP41fIB77oHnnoPmzWHVKpg0Cfr29TQ2ERG/CuuSzq6Y\nGNZlZbHziisA2D9qFOvWrePAgQMeRyYiEhl8k/A//PxzHnjgAR5cupQCY4ibNo3H7rmHZctK3JJH\nREROIlj98AMyDNh18CAAeypXZn7durTftImuWVmexiUi4jdHWyuUhS9G+MOAlPj4Y69nNG4MQI+1\na0FtF0REjgmkW6YvEj7AwZiYY8/n163Lnrg4UvbupcqSJR5GJSISOfyT8Cv9XF0qiIpiZqNGANSZ\nPNmrkEREIopvEn5epV/eTsgoLOvEf/YZaKaOiEjAfJPwD0f9MpQNp5/Oj/HxRO/fDx9+6FFUIiKR\nwxcJfxiwfteu/zk+IzXVPRk7tjzDERHxrYhYaXtb377sqlz5F1+rkpfHmClTiDp0CNauhaM/AERE\nKriwXmmbb/437gOxsezq0cO9eOONco5IRCSy+CfhR504lB39+rknY8fCkSPlGJGISGTxTcIvKCbh\n72vfHho1gqwsmDGjnKMSEYkcvkn4xY3wiYqC6693z3XzVkSkzHyT8I+coIZ/zNGE//77sHt3ucQj\nIhJpfJHwHzaGrdu2Ff+G1FS48EL46Sd4++1yi0tExG8CmZbpi4Q/JCaGxMTEk7/pllvcf0eNUkM1\nEamwwr55WkF09KnfdMUVUL8+LFsG06aFPigRkQjji374R4q7YVvEgGuv5Ve1ajFw40a+v/56nura\n9djXJkyYEMrwREQigj9G+JVK9nPn8zPPJC86mnabN1N3374QRyUiEln8kfBLMMIH2BcXx6yGDQHo\ntWpVKEMSEYk4vkj4toQJH2BK06YAdF+3jip5eaEKSUQk4vgj4Z9sDv5x1tesyaLERCrn53PJ6tUh\njEpEJLKEXcIHmHTWWQD0XrmSmPz8UIQkIhJxwjLh/5CUxJpatah56JDb6FxERE7JFwl/RG4uW7du\nLfk3GMOHLVoA0HfFCjh8OESRiYj4S9ivtL2/evVTr7Q9zrf167OxenUSDxyA8eNDFJmIiL+E/Urb\n0pZ0jn7P0Vo+Tz6pXvkiIqcQtgkfYGbDhmyrUgWWLoV33w1yVCIikcUXCb+sCqKj+aBlS/di6FDQ\njB0RkWL5IuGXdYQP8GVqKpx5JqxcCePGBS8oEZEIE/YJvyAqCoYPdy+GDwetvhUROaGwT/gAXHMN\ntGzp9r199dXgBCUiEmEiI+FHR8Mjj7jnjz4K+/cHHpSISITxRcIn0IQPboOU9u1h82Z45pnAP09E\nJML4IuEHPMIH90Pj2Wfd8xEjYNOmwD9TRCSChDThG2MaG2NeNcb8N5TnOaZbNzfSP3AAhgwpl1OK\niISLkCZ8a+1aa+1Np3xfMEb4Rz31FMTEwOuvw/z5wftcEZEwV6KEb4x5zRiTY4xZdNzxXsaY5caY\nVcaYB8oaRFATfpMmMHgwWAt//KP7r4iIlHgT87HAi8CbRw8YY6KBUUBPYCPwrTFmEtAeaAc8ba0t\nUSE9GAl/wIABx55XzcvjhdhYqs+YwXNdunBvZmbAny8iEu5KNMK31s4Edh13uCPwo7V2nbX2MDAe\n6G+tHWetvddau8kYE2+MGQ20PdlvAEEd4QP7Y2MZn5YGwO8XLABteC4iUuIR/onUBzYUeZ0NdCr6\nBmvtTuC2U33Q8zt3snjxYgASExNL3Sr5RD5v3Jgea9fSZOdOtwJXUzVFJIxlZGSQkZER0GcEkvCD\nVhy/p3ZttrRuHayPA9zG6K+2a8fj06cT9fzzcP31EORziIiUl/T0dNLT04+9Hn60pUwpBDJLZyOQ\nUuR1Cm6UX2rBLukctTY+nmlnngkFBXD77bqBKyIVWiAJfx7Q1BiTaoyJBQYAk8ryQS/s2FG6LQ5L\nYUJaGiQmwqxZMHZsSM4hIlJeQr7FoTHmHSATaGaM2WCMGWStzQcGA1OBpcAEa+2ysgRxV0JCUOr2\nJ7I/Nhb+8Q/34r77tAJXRMJayLc4tNYOtNbWs9bGWWtTrLVjC49PsdY2t9Y2sdY+UaYIysNvfgOX\nXw579sBtt6m0IyIVki966Yzcvj1kJR3A9dkZPRpq1IDJk+Gdd0J3LhGREAp5SSfUBgdpKubJDLjv\nPkY3awbAvkGDuKV/fwYMGPCLBVsiIn4X8pJOyIVols7xZjRuzMKkJKrn5XHD99+rtCMiFYovEv6L\n27aFtqRzlDG80r49BytVonN2Nt2yskJ/ThGRIAr7ks6dSUkhL+kctb1qVcaecw4AN3z/PYm5ueVy\nXhGRYAj7kk6oFl4V58vUVGY3aEDl/HwGz5kD+fnlen4RES/4IuGXO2N45dxz2VG5Ms137IDHH/c6\nIhGRkPNFwh+Vk1M+Nfwi9sfF8VLHju7FI4/AN9+U6/lFRMoi7Gv4d9StW241/KIWJyUxqXlz12vn\nmmtgx45yj0FEpDRUww/A+NatoWNHWL8efv97OHLEs1hERELJFwnfSwXR0fDuu1CrFnz6KYwY4XVI\nIiIh4YuE/9KWLeVew/+FRo1g3Dj3/K9/hS+/9C4WEZGTCPsa/m316nlSwy9qwJtv8tFZZ8GRI+zs\n1UutF0TEl1TDD5IJrVuzJCGB+J9+4t7MTKILCrwOSUQkaJTwizgSFcULnTuzo3JlWmzfzg3z56vf\njohEDF8kfD/ZU7kyz3bpQl5UFD3XrOGS1au9DklEJCj8kfB9MsI/anV8PC936ADAdfPnQ4A7xYuI\n+IEvEv7ojRu9naVzArMaNWJS8+ZUshauugrWrvU6JBGR8J+lc2uDBp7P0jmRt9PSmJ+c7FbgXn45\n7N7tdUgiUsGF/Swdv7JRUYzs3BlatoSlS+GKKyAvz+uwRETKxBcJ38/zYA7ExroVuMnJMGMG3Hij\nZu6ISFjyRcL3vUaN4JNPoGpV+M9/4OGHvY5IRKTU/JHwfTZL54TatXM9d6Ki4NFHYcwYryMSESkV\nfyT8cHHZZfDSS+75rbfCxInexiMiUgq+SPgvZ2f7blpmsW69FYYMcT30BwxwdX0RkXIS9tMyb/Hp\ntMxiPfII3HEHHDoE/frB3LleRyQiFUQg0zIrBTeUCsIYePFFNy//7behd2/u69CBjaeffsK3T5gw\noZwDFBH5X74Y4YfFTdvjRUXB669Dnz6wcydDvvqKhNxcr6MSESmWPxJ+uIqJgf/+F7p1I/7gQR7+\n8kvq7N/vdVQiIifki4Tvl/bIZVK5MkyezKr4eBL37+fhjAxqK+mLiA/5IuGHvRo1eLxbN36Mjydp\n/36GZmRQ+8ABr6MSEfkFJfwgORAby2PdurG6Vi2S9u/n4RkzlPRFxFeU8IPoQGwsj3XvzupatUgu\nLO/EK+mLiE/4IuG/smFD+Cy8OoX9RZN+bi7DZ8wA7ZolIkESyMIrX8zDv7lhQ571+cKrAQMGFPu1\n4+fZH036D331FU137oSuXWHaNGjVKtRhikiES09PJz09neHDh5f6e30xwo9E+2NjebR7d5YkJMDm\nzdC9O3z3nddhiUgFpoQfQj/FxPBE165ucdaOHdCjB8yc6XVYIlJBKeGH2OFKleDDD12jtX374NJL\n4bPPvA5LRCogXyT8sF54VRIxMfDWW3DTTXDwIPTtC2+84XVUIlLB+CLhVwjR0fDKK/Dgg5CfD9df\nD489pu0SRaTcKOGXJ2PgiSdg1Cj3fMgQuP129wNARCTElPC9cOed8P77cNpp8PLLcMUVoAVaIhJi\nSvhe+fWvYfp0qFULJk+GCy+ECFl8JiL+5IuEH/E3bYvTpQt8/TU0agRz5kDHjrBokddRiUiE8kXC\nr9BatIDZs12yz8qC8893I34RkSALecI3xvQ3xrxijBlvjLk41OcLS8nJkJEB11wDubnQvz8884xm\n8IhIUIU84VtrJ1prbwFuA4pvSFPRVa7s9sf9+99dor//frjxRsjL8zoyEYkQJW6eZox5DegDbLXW\nphU53gt4HogGXrXWPlXMRwwBRgUQa1grUfO1o1M1zzoLfv97GDsWVq1y2ygmJ5dTpCISqUozwh8L\n9Cp6wBgTjUvivYCWwEBjTAtjzO+MMc8ZY+oZ5ylgirV2wYk+uMLetC3OVVe5njv16sGsWXDuufDN\nN15HJSJhrsQJ31o7E9h13OGOwI/W2nXW2sPAeKC/tXactfZea+0m4C7gIuAqY8ytwQo84p17ruuu\necEFsGmT67Y5erTq+iJSZoH2w68PbCjyOhvoVPQN1tqRwMiTfcirWVks3rEDgMTERBJ93hs/FIor\n+UQnJ/P23XfDyJFuVe7cufDSS27RVgk/A/63Z7+IhJeMjAwyMjIC+oxAE35Qhps3NWpEbt26wfio\niFMQFQUvvOCmbd58s6vr//CDW6nbqJHX4YlIOTm68clRZdkAJdCEvxFIKfI6BTfKLxXV8E9twKRJ\nNOralT9+/TVJ333H3mbN+GenTiwo/EGpEbyInEqg0zLnAU2NManGmFjctMtJpf2QMevWRcyetqGU\nVbMmD/Xsyfd161IjL4+HZs7ktwsXEn3kiNehiUg5CWRP2xInfGPMO0Am0MwYs8EYM8hamw8MBqYC\nS4EJ1tplpQ3ixtTUClm3L4v9cXGMuOAC3kpLo8AY+q1YwbAZM9wqXRGJeOnp6aHfxNxaO7CY41OA\nKWU6u5SJNYZJLVqwPCGBe2bPptmOHdC2Le1btWJe/fpehyciPuWLXjoq6ZTNyjp1eODii5lXrx7s\n3s39X3/NdfPnU6mgwOvQRCREyqWkE0oq6ZRdblwcT3fpAv/4B/lRUVy2ahWPff45Dfbs8To0EQmB\nQEo6vkj4EiBj4N57ebhHD7ZUq0bq7t08MW0avVauxGihlogU8kXCV0knOFbXrs0DF1/MF40bE3vk\nCIMWLOChr76i1sGDXocmIkGiko4c81NMDC936MAz55/P3thYzs7J4empU91CLREJeyrpyP/4tkED\n7r/0UhYkJ1M9L881ZBs0CFTbF6mwlPAj2O7KlXmia1deO+cc13vn9dehdWuYOtXr0ETEA75I+Krh\nh5AxTG3aFL7/Hjp0gOxs6NULbrpJo32RMKQavpxaixaQmQlPPgmxsTBmjBvtT9GaOZFwohq+lEyl\nSvDAAzB/vuu+mZ0Nl10GN9wAu3d7HZ2IhJgSfkXUsiV8/TWMGAFxca7lcqtWMKnUfe9EJIz4IuGr\nhu+BSpXcRunz50OnTm5Xrf794corYeNGr6MTkWKohi9l16KFG+0//zxUrQoffOCO/fOfoJ48Ir5T\nLt0yJfydcgvEK66AwYNdaWfwYBg3Dl55Bdq0KccoRSRUfDHCF59ISYGPPnKrcuvVgzlz3GbqDz4I\nBw54HZ2IBEgJX37JGDfSX7oU7rjDlXWeespN4Zw82evoRCQAvkj4umnrQ6ef7ur4X38NaWmwdi30\n6weXXw4//uh1dCIVlm7aSuicdx589527qVujBnzyiZvCOWSIyjwiHtDCKwmtmBi45x5YsQKuuw7y\n8uCxx9xsnvffB/XcFwkLSvhScsnJrgHbrFnQti2sX++6cF56KSxf7nV0InIKSvhSel26wLx5rsZf\nsyZMm+bq/PfcAzt3eh2diBRDCV/KJjrazeJZuRJuvtnN5hk5Epo0cfX+vDyvIxSR4yjhS2ASEtzi\nrPnz4aKLYNcuuPdeN41z0iTV90V8xBcJX9MyI8DZZ7vSzuTJ0Lw5rFrlevP07AkLFngdnUjE0LRM\n8Qdj3Dz9RYtceSc+Hr74Atq1cxuuqCmbSMA0LVP8JSYG7rrLjfL/8AdX7x8zBpo2hYceUu99EY8o\n4UvoxMfDc8/BkiWu7fLBg27HrTPOgKefdq9FpNwo4UvoNWsG773nmrH16OFu7P75z+74mDGQn+91\nhCIVghK+lJ+OHeHzz+Gzz9zCrexsV9tPS4MPP9SMHpEQU8KX8mWMW5n73Xfw9tuuvLN8uevQed55\nMH26Er9IiGgDFDnmlBukBFNUFAwc6Gr7r7wCf/+7K/lcfDHL6tTh3datWVpk5lZJzn+q+Mv1zyfi\nQxrhi7diY93uWqtX83ZaGvtiY2mxfTtDMzIYkpFB8+3bvY5QJGJohC/+UK0aE1u04P+aNKH3qlX0\nWbGCtK1bSfviCxYmJcHs2dC5s9dRioQ1X4zwtdJWjjoYE8MHLVtyV58+vNeyJQcqVeLsnBxX3+/T\nxzVtE6nAtNJWIs6B2Fj+27o1g/v04YMWLaBqVfj0U+jQwa3m/eYbr0MU8YRW2krE2h8Xx4SjWyz+\n+c9QpYrbdev8812zthkzNKtHpISU8CU8JCS4zdTXrYO//MVtt/jFF3DhhXDBBW70r8QvclJK+BJe\nEhLc9opZWW4qZ3w8ZGZCnz48OW0aHbOzMUr8IiekhC/hqWZNt5F6VhY88wwkJ9N4927+mJnJ01On\nckFWFlFHjngdpYivKOFLeKtWDf74R1i7ljHt2rGtShVS9u7lrjlzeOHTT+m9ciVxhw97HaWIL2ge\nvkSG007j/5o04fPGjemalUX/FSuot28f1y9YwFVLl8KQIZz+00/sOe00ryMV8YwSvkSUguhoMs44\ngy8bN+bcTZvot3w5zXfsgMceY1RUFF+lpvJxs2ZsrlHD61BFyp0SvkQkawzz6tdnXv36NN++nUeq\nVyd24kR6rlnDhWvWMK9+fSY1b86qOnW8DlWk3CjhS8RbUacOTJjAvb17c/mKFXTLyqLjxo103LiR\n5XXq8HGzZlBQ4HbmEolgIb1pa4w5yxjzL2PMu8aYG0N5LpFT2VSjBq906MCdhat3c2NiOGv7dv6U\nmem2X3zuOdizx+swRUImpAnfWrvcWns7cA1waSjPJVJSeypXZkJaGndefjmvt23LlqpV3Ure++6D\nBg3g7rvdfrwiEaZECd8Y85oxJscYs+i4472MMcuNMauMMQ8U8719gU+A8YGHKxI8P8XEMKVZM/7Q\nuzd89JHbfjE3F158EZo3h3793A5dWsglEaKkNfyxwIvAm0cPGGOigVFAT2Aj8K0xZhLQHmgHPG2t\n3WStnQxMNsZMBD4IZvBSvgLdYMSvG5DYqCjo3989Fi6EkSPhrbdg8mT3aN0a/vAH+M1voHLlMp/n\nZH9+CM4mL6EW6vP7/fPD/e+wRCN8a+1MYNdxhzsCP1pr11lrD+NG8P2tteOstfdaazcZY7obY14w\nxrwMzAhu6CIhcPbZbmP1DRtc64bkZFi82O29m5ICDz3k+vmIhKFAavj1gQ1FXmcXHjvGWvultfYe\na+2t1trnAziXSPlKSPi5dcO4cXDuubBjBzz5pNuHt29ftxm72jdIGAlkWmbQCptj1q1jceFWdomJ\nieqNL/4RGwvXXgu//a3rwf+vf8G778LHH7vHGWfA7bfDoEFQu7bX0UoEy8jIICMjI6DPCCThbwRS\nirxOwY3yS+3G1FRyk5MDCEUkxIxxPfjPPx+efRZeew1Gj4Y1a+D++91vA9dcA3fcAR07eh2tRKD0\n9HTS09OPvR4+fHipPyOQks48oKkxJtUYEwsMACYF8Hki4SExER58EFavdjd1e/eGvDx44w3o1Mnt\nyjV2LBw44HWkIr9Q0mmZ7wCZQDNjzAZjzCBrbT4wGJgKLAUmWGuXlSUI7WkrYSk62m23+Omnbt7+\nn/7k+vPPmwc33AD168Ndd8EPP3gdqUSQkO9pa60daK2tZ62Ns9amWGvHFh6fYq1tbq1tYq19okwR\noD1tJQKceSY8/TRkZ8Prr7tR/u7dMGqUm/nTqRMXrlmjVs0SsED2tPVFL50x69axNSpKSV/CX+XK\ncN117rFgAfz7325O/9y53Ar8fsECvm7YkC/OOIPVtWq5ewMipRDIzVtfbICiEb5EpLZt4Z//hE2b\n4PXXWV6nDpXz8+m5Zg2PT5/OU9OmccmqVVTJy/M6UgkjYT/CF4loVarAddcx9NNPqb9nDxeuXUu3\ndetI3b2bG+fP59offmB2gwYwc6bbkF2jfgkRX4zwRSqKjaefzri2bbm9b19e6NyZRYmJxBUU0D0r\nC7p1g2bN4NFHYf16r0OVCOSLEb5q+FLR5EdHk9mwIZkNG5KUm0uPNWv49d698OOP8Le/wcMPw0UX\nwfXXw69/7X5LEEE1fJGwllOtGuPbtHGj+s8+cwu4YmNh+nS3yrduXbjlFsjMVOdOCaiG74uELyK4\nef2XXgrvvAObN7s2Dp06wd69brZPly48N2UKv1q6lHgt6pIyUElHxI9q1YLbbnOPZcvcKt4336Te\n5s0MXLyYAYsXsygpia9SU/m2fn0OVfLFP2UpByrpiESyFi1cl87163mia1cyU1IoiIri7Jwc7poz\nh1cmTmTw7NkwZQrk53sdrYSYpmWKVASVKrGgbl0W1K1L1UOHOH/DBi5Yv56ztm+n6/r1cNllrs/P\nNde47p4dOmiKp/yCEr5IGNofF8e0Jk2Y1qQJibm5XLB+PQMOH4YVK9yOXSNHuo3Zf/tb92jSxOuQ\nxQd8UdIRkbLbWq0aH7Rs6Wr98+a57RiTklxDt2HDXOLv3Nn19dm2zetwxUO+SPjqlikSBMa4nbme\ne841cZs6FX73O6haFebMcZ0769Z1pZ8334Q9e7yOWMog5N0yQ003bUWCrFIluOQSl9hzcuDtt6FP\nH/e1KVNcc7fERPjVr2D8eMjN9TZeKTHNwxeR4lWtCgMHui0ZN2+Gl16C7t3h8GGYONF9LTERrr4a\n3n8fDh70OmIJESV8kYokIcHtwZuR4co+L7wA553nkvx//wtXXeWS/7XXut28Dh3yOmIJIiV8kYqq\nXj24+27XsiEry23g0r69K++89Rb06+du/t5wg7sfoM1bwp4vEr5u2op4rGFDt0Xjt9+62T2PPQZt\n2rgbu2PHQq9eULcut377LW03b6ZSQYHXEVdYumkrIsHTpAn85S+wcKGb6jlsGJx1FuzYwYVr1/LQ\nzJn8e9IkBs+eTYfsbGK0urdc6aatiITGWWfB0KGwdCksXsy7rVqRdfrpVDl8mK7r1/OnzExenTSJ\nP2RmwoQJsG+f1xHLSWilrYicmjHQqhXvFz6S9+2jY3Y2nbOzOXPXLs7LznYtHeLiXMfPK6909wBq\n1vQ6cilCCV9ESm1L9epMatGCSS1aUGf/fjpmZ3NdtWruBvCkSe4RE+M2cbnySujf380QEk+ppCMi\nAdletSqfNm8Os2bBxo1u4/YePaCgwG3ocvPNkJwMXbvCM8+4Xb3EE0r4IhI8devCHXfAF1/Ali1u\n45ZevdzmLrNmwf33u94+rVrBX/8Kc+fCkSNeR11h+KKkow1QRCJQQgLcdJN77N3rRvsTJ8Inn7ib\nwEuXwuOPux8S/ftD//5UKiggPzra68h9TRugiIi/1ajhWje89RZs3QrTpsHgwZCS4to9jB4NvXvz\n74kTueebb+iSlUWVvDyvo/YlbYAiIuEjNhZ69nSPkSNh/nw38p84kSoLF3L+hg2cv2ED+cawLCGB\nb+vX5/t69byOOiIo4YuId4yBdu3cY/hwBvfpQ/tNm2i/aRMttm0jbetW0rZu5Yb586F1a9fx8/LL\nXf8f7eNbarpiIuIb26pVY0qzZkxp1oyqhw7RbvNm2m/aRJucHKosWQJLlsCIEW6T9169XPLv1Qvi\n470OPSwo4YuIL+2Pi2NmaiozU1OJLijg7TvucDd8P/4YVq6Ed95xj6goOP/8n0f/rVppL99i+OKm\nrYjIyRRER8OFF8Kzz7p9e1eudDt7XXTRz1M+H3oI0tIgNRXuvBM+/VS9/Y+jhC8i4adpU7d37/Tp\nsH07vPceDBrkevmvX+82eenTB2rXhr594V//grVrvY7acyrpiEh4q1HDtW+48kq3iOu771zZ5+OP\n4fvvf34O/KN6dRYmJ7MgOZmlCQkcrmA3fivWn1ZEIltUFHTo4B7Dh8OmTa6089lnMG0a9ffupf6+\nfVy2ahV50dEsTUhgYXIy85OTwdqIr/37IuFrpa2IhES9ej+v9j18mIcvvZS2W7Zw9pYtnLlrF223\nbKHtli1cB3DGGW7GT69e7n5B9epeR39Cgay09UXCvzE1lReU7EUklGJiWJGQwIqEBCakpVHjp59o\nk5PDOZuiMAJLAAAG+0lEQVQ302bLFmqsW+dW/I4e7eb4X3DBzz8A2rTxzeg/PT2d9PR0hg8fXurv\n9UXCFxEpb3tPO41ZjRoxq1EjzJEjjP/zn13p57PPYPZst9F7RgY8+KDr93PJJVyQlcWipCT2nHaa\n1+GXiRK+iFR4tmjt/29/g127XL+foz8ANm+GN97grsL3Z51+Oj8kJbEoKYllCQnkhcnN3/CIUkSk\nPNWq5Zq9XX21u5m7aBFMm8aCZ56hxbZtNNqzh0Z79tB35UoOR0Wxok4dFiUlwbx5cM45bm2ADynh\ni4icjDGuht+mDU/MnUtMQQHNtm+nTU4OaTk5NN61i9Zbt9J661b3G0J8vLvpe/HF7tG4sdd/gmOU\n8EVESuFwdDRLkpJYkpTEO0C1Q4dovXUraTk59DxyBNatcwvB3nvPfcOZZ/6c/Hv08DJ0JXwRkUDk\nxsUxOyWF2Skp9Bw/HlavdvX/6dPdzl+rV7vH6NEQFcVjNWuyODGRxUlJLK9du1wXfynhi4gEizHQ\npIl73H475Oe7lb9HfwBkZtJk506a7NzJr5Yv53BUFCtr12ZxUhKLEhPh8GG3+XuIKOGLiIRKpUrQ\nqZN7DBkCubk8ftllx0pAjXbvptW2bbTato0BAN27Q2Zm6MIJ2SeLiMgvVavGwrp1WVi3rnt56BAt\nt20jLSeH1lu3Uq9z55CePuQJ3xhTFcgAhllrPwn1+UREwkVuXBxzGzRgboMGAEwYMSKk5yuP9sh/\nBiaUw3lERMJbiG/glijhG2NeM8bkGGMWHXe8lzFmuTFmlTHmgRN838XAUmBbcML1n61bt3odQkDC\nPf6yNpHyi3C//uEcfzjHXlYlHeGPBXoVPWCMiQZGFR5vCQw0xrQwxvzOGPOcMaYe0B3oDPwGuNkY\nn3QfCqJw/58m3ONXwvdWOMcfzrGXVYl+f7DWzjTGpB53uCPwo7V2HYAxZjzQ31r7JDCu8D1DCr92\nHbDNWmuDELOIiJRBIAWj+sCGIq+zgU4neqO19o0AziMiIkFgSjroLhzhT7bWphW+vhLoZa29ufD1\ntUAna+1dxX7IiT9Xo34RkTKw1paqTB7ICH8jkFLkdQpulF8qpQ1YRETKJpBpmfOApsaYVGNMLDAA\nmBScsEREJNhKOi3zHSATaGaM2WCMGWStzQcGA1NxUy8nWGuXhS5UEREJRIkSvrV2oLW2nrU2zlqb\nYq0dW3h8irW2ubW2ibX2idKe/FTz+P3OGLPOGPODMWa+MWau1/GcyonWUxhj4o0x04wxK40x/2eM\nqelljCdTTPzDjDHZhX8H840xvU72GV4xxqQYY2YYY5YYYxYbY+4uPB4W1/8k8YfL9T/NGDPHGLPA\nGLPUGPNE4fFwuf7FxV+q61/im7bBVjiPfwXQE3c/4FtgYDj9lmCMWQuca63d6XUsJWGM6QrkAm8W\nufk+AthurR1R+EO3lrX2QS/jLE4x8Q8F9llr/+FpcKdgjEkGkq21C4wx1YDvgF8BgwiD63+S+K8m\nDK4/gDGmirX2gDGmEjAL+BPQjzC4/lBs/BdRiutfHq0VinNsHr+19jAwHujvYTxlFTY3na21M4Fd\nxx3uBxydNvsG7h+xLxUTP4TB34G1dou1dkHh81xgGW5qc1hc/5PED2Fw/QGstQcKn8YC0bj/l8Li\n+kOx8UMprr+XCf9E8/jrF/Nev7LAdGPMPGPMzV4HU0ZJ1tqcwuc5QJKXwZTRXcaYhcaYMX79lbyo\nwinO5wBzCMPrXyT+2YWHwuL6G2OijDELcNd5hrV2CWF0/YuJH0px/b1M+JEw/76LtfYcoDdwZ2HJ\nIWwVroQOt7+XfwGNgbbAZuBZb8M5ucJyyPvAPdbafUW/Fg7XvzD+93Dx5xJG199ae8Ra2xZoAHQz\nxvQ47uu+vv4niD+dUl5/LxN+UObxe8lau7nwv9uAD3FlqnCTU1ifxRhTFwirBiPW2q22EPAqPv47\nMMbE4JL9OGvtR4WHw+b6F4n/P0fjD6frf5S1dg/wCXAuYXT9jyoSf/vSXn8vE35Yz+M3xlQxxlQv\nfF4VuARYdPLv8qVJwHWFz68DPjrJe32n8B/pUb/Gp38HxhgDjAGWWmufL/KlsLj+xcUfRte/ztFy\nhzGmMnAxMJ/wuf4njP/oD6tCp7z+ns3SATDG9Aaex92AGFOWqZ1eMcY0xo3qwa1Yfsvv8Ru3nqI7\nUAdXB3wYmAi8CzQE1gFXW2t3exXjyZwg/qFAOu7XWQusBW4tUpP1DWPMBcBXwA/8XDZ4CJhLGFz/\nYuL/CzCQ8Lj+abibslGFj3HW2qeNMfGEx/UvLv43KcX19zThi4hI+fGypCMiIuVICV9EpIJQwhcR\nqSCU8EVEKgglfBGRCkIJX0SkglDCFxGpIJTwRUQqiP8HGeGJLbl6mIoAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10ab38150>"
       ]
      }
     ],
     "prompt_number": 7
    }
   ],
   "metadata": {}
  }
 ]
}