ColCarroll/1d sampler.ipynb Secret

## 1d sampler.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sample from a pdf\n",
    "\n",
    "The goal here is to take a pdf, and generate samples from that pdf.  This is in support of a javascript library doing the same, so I've only used the standard library, instead of leaning too hard into `numpy`.  The first part defines the functions used, and the last bit is a demo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def area(x, y, x_prev, y_prev):\n",
    "    \"\"\"Area of the trapezoid with vertices ((x_prev, 0), (x_prev, y_prev), (x, y), (x, 0))\"\"\"\n",
    "    return (x - x_prev) * 0.5 * (y + y_prev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def integrate(tuples):\n",
    "    \"\"\"Trapezoidal rule?\"\"\"\n",
    "    x_prev, y_prev = tuples[0]\n",
    "    tot = 0\n",
    "    for x, y in tuples[1:]:\n",
    "        tot += area(x, y, x_prev, y_prev)\n",
    "        x_prev, y_prev = x, y\n",
    "    return tot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_cdf(tuples):\n",
    "    \"\"\"Actually turning the tuples into a list of dictionaries for future processing\"\"\"\n",
    "    tot = integrate(tuples)\n",
    "    data = sorted([(x, y/tot) for x, y in tuples])\n",
    "    x_prev, y_prev = data[0]\n",
    "    cdf = [{'x': x_prev, 'x_prev': x_prev, 'cdf': 0, 'a': 0, 'b': 0, 'c': 0}]\n",
    "    for x, y in data[1:]:\n",
    "        prev_cdf = cdf[-1]['cdf']\n",
    "        delta_cdf = area(x, y, x_prev, y_prev)\n",
    "        # computing interpolating quadratic on the fly, using the equation\n",
    "        # y(z) = a (z - x_prev)^2 + b(z - x_prev) + c; x_prev < z <= x\n",
    "        cdf.append({\n",
    "            'x': x,\n",
    "            'x_prev': x_prev,\n",
    "            'cdf': prev_cdf + delta_cdf,\n",
    "            'c': prev_cdf,\n",
    "            'b': y_prev,\n",
    "            'a': 0.5 * (y - y_prev) / (x - x_prev),            \n",
    "        })\n",
    "        x_prev, y_prev = x, y\n",
    "    \n",
    "    return cdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# not used, but might as well compute\n",
    "def get_cdf_fn(cdf):\n",
    "    def cdf_(x):\n",
    "        for data in cdf:\n",
    "            if data['x'] >= x:\n",
    "                m = x - data['x_prev']\n",
    "                return data['a'] * m * m + data['b'] * m + data['c']\n",
    "        return 1\n",
    "    return cdf_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_inv_cdf_fn(cdf):\n",
    "    \"\"\"Get a function from a probability back into feature space\"\"\"\n",
    "    def inv_cdf(y):\n",
    "        \"\"\"This could probably be a binary search instead ¯\\_(ツ)_/¯\"\"\"\n",
    "        if y < 0:\n",
    "            return float('-inf')\n",
    "        for data in cdf:\n",
    "            if data['cdf'] > y:\n",
    "                #  This gets a little hairy because the quadratic formula\n",
    "                #  Has to shift back to be centered around 0, not x_prev\n",
    "                x_prev = data['x_prev']\n",
    "                a = data['a']\n",
    "                b = data['b'] - 2 * data['a'] * x_prev\n",
    "                c = a * x_prev * x_prev - data['b'] * x_prev + data['c']\n",
    "                return (-b + (b * b - 4 * a * (c - y)) ** 0.5) / (2 * a)\n",
    "        return float('inf')\n",
    "    return inv_cdf\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sample_from_pdf(pdf, n):\n",
    "    \"\"\"Generate n samples from a list of (possibly unnormalized) x,y tuples.  Assumes all y > 0.\"\"\"\n",
    "    inv_cdf = get_inv_cdf_fn(get_cdf(pdf))\n",
    "    for x in np.random.random(n):\n",
    "        yield inv_cdf(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example usage\n",
    "\n",
    "We can generate a pdf, then sample from it, and demonstrate that the histogram of samples looks like the \n",
    "plot of the pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 20, 1000)\n",
    "y = np.maximum(np.sin(x), 0)\n",
    "pdf = np.vstack((x, y)).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmQZNdVp7+TW2VV1r713uqWumWpJcu2aAt7BMbYwpbM\njMSYTYohMOBA4cBiMGuYgFA4zB+DceCJYPDYGOwwEHiFARojj9AYgzFYtlqW3FJLlrvV6q26a8vM\nWjKzcr/zx3svK7u6sioz67182/0iKroy81W+2zdv/u65555zriil0Gg0Gk2wiLjdAI1Go9HYjxZ3\njUajCSBa3DUajSaAaHHXaDSaAKLFXaPRaAKIFneNRqMJIFrcNRqNJoBocddoNJoAosVdo9FoAkjM\nrRtPTk6qQ4cOuXV7jUaj8SVPP/30olJqarvrXBP3Q4cOcfLkSbdur9FoNL5ERC60c512y2g0Gk0A\n0eKu0Wg0AUSLu0aj0QQQLe4ajUYTQLS4azQaTQDZVtxF5FMiMi8iz7d4XUTkj0TkrIicEpE77W+m\nRqPRaDqhHcv908C9W7x+H3DU/HkY+NjOm6XRaDSanbCtuCulvgZktrjkAeAvlMGTwKiI7LGrgb2k\nXlf846mrfPLrrzC/WnS7OaFjMVfiU19/hS+dukKtro9/7DUnz2f4k399me/OrrjdFI0N2JHEtA+4\n1PT4svnc1Y0XisjDGNY9Bw8etOHW9lGvK973+Wc58Z0rAHzsX17mi+95I4cnUy63LBxcTBd458f+\ng8VcCYAfvWOWP37odYiIyy0LB3/xjfM8+venAfjw4y/xsZ/5Pn7k2C53G6XZET3dUFVKfUIpdVwp\ndXxqatvs2Z7yxacvceI7V/jVe27mS7/8A1RqdX71889S1xak4yil+K2/+Q7lao0Tj9zNb7ztZv7x\n1FU++61L2/+xZsecW8jxe196gbfeMs2//dYPc2zvML/+hWdZWC253TTNDrBD3GeAA02P95vP+YZy\ntc4f/tP3OH7DGP/9rUe4fd8Iv/ujt/LspSUePz3rdvMCz1denOfJcxl+8+2v4o79o7z3h49w16Fx\n/ugrZyhVa243L/D8/pe/S18syu//+B0cGB/gIz/1WnKlKn/29XNuN02zA+wQ9xPAz5pRM28AlpVS\n17lkvMyXn7/K/GqJ977lSMMN8M4793NwfIA//8Z5V9sWBj79H+fZN9rPQ3cZrjoR4b1vOcLsSpG/\nf/aKy60LNpcyBZ54cY6fv/sQU0N9AByZHuS+V+/hM09eZLVYcbmFmm5pJxTys8A3gFeJyGURebeI\nvEdE3mNe8hhwDjgL/CnwS4611iH+6smLHJ5M8UNH111F0Yjw068/wJPnMlxI511sXbC5mC7w9bOL\n/PTrDxCLrg/HNx2d5KapFH/99GUXWxd8Pv/UJQR48K5r98B+4e7DrJaq/N/n9crVr7QTLfOQUmqP\nUiqulNqvlPqkUurjSqmPm68rpdR7lVI3KaVerZTyVanH+dUiT13I8MBr9xKJXLt59+N37kcE/vYZ\nX3mZfMU/nDIs85/4vv3XPC8i3P+afTx1PsPV5TU3mhZ4lFL8w6kr3H1kkn2j/de8dufBUfaP9fOl\nU75ahGuaCH2G6uOn51AK3vHq66M3d48kufPgGP/83XkXWhYOnnhhjtccGGXvBnEBuP+1e1EKvvyc\nth6d4Mx8jgvpAm+/bfd1r4kI/+U1e/n3s4ssFcoutE6zU0Iv7l95cY7DkymOTg9u+vpbbpnm1OVl\n5lZ03LvdzK8WefbSEvfcMr3p64cnU9w4leJrZxZ63LJw8MQLcwDcc+vmIY/33DpNta74j5fTvWyW\nxiZCLe6VWp2nXslw95GJlvHUbzGF519f0gJjN18/swjAD7cQd4A3HZ3iyXNpHTXjAF8/s8ixPcPs\nHklu+vpr9o8y1Bfj3/Tk6ktCLe7PzSyTL9d4442TLa+5ZfcQ46kET76irRe7+dYrGYaTMY7tGW55\nzQ8enaRYqfP0+WwPWxZ8ytU6376Y5Q03TrS8JhaN8MabJvja9xZRSud7+I1Qi/s3zOXm99843vIa\nEeGuQ+N865WtKjBouuFb5zO8/tD4dRvZzdx1eBwReEqLu608N7NEqVrnrsOtxz7A3UcmmVla48qy\ndkv6jVCL+5Pn0ty8a5DJwb4tr/v+G8e5nF1jZklHbdjFwmqJcwt5Xr+NuAwl47xq1xDfvqjF3U6+\naRorrz80tuV1dx40Xn9G97/vCK241+uKZy4ubWu5ALz+kHHNyfPaercLqy/b6f/XHRzlmYtZXQrC\nRr71Soaj04NMbGPY3LJniL5YhG9fWOpRyzR2EVpxP7eYJ1eqcsf+0W2vfdXuIRKxCM/PLPegZeHg\n6QtZ+mIRbt87su21rzs4xkqxyrnFXA9aFnyUMgyb49tY7QDxaIQ79o/wzCVtufuN0Iq7JdR37N9e\nXOLRCMf2DHPqshZ3u3j+yjK37hkmEdt+CFquAW092sPl7BrLaxVu37f92Aej/0/PrOiIJZ8RWnE/\ndXmZZDzCkanN49s3csf+EU5fWdGuARuo1xWnZ1a4fV/rKJlmbpxMMZyM8exlLe52cPqKYaS0s2oC\nwy1WrtV54Yqu8+4nQivuz80scdvekWvqmWzF7ftGyJWqvKLrzOyYS9kCq6Vq2+ISiQi37hnmxata\nXOzguZllohHhVbuH2rr+NvNzevHqqpPNCgVKKb506gqzPYg+CqW41+qK52dWeHWby1JYd988p10z\nO+b5GUOk23ULANy6Z5iXZlf1yskGnp9Z4ej0IMl4tK3r94/1M9gX05OrDVxZLvLIZ57hiRfnHL9X\nKMX9fDrPWqXGbXvbcwsAHJkapC8WaSxpNd3z/JVl4lHh6K72XGIAx/YMUyjXuJApONiy4KOU4vmZ\n5Y4mVhHhlt1DWtxtwHJtbZW4ZxehFPczc8byst1lKRjZekemB/nenI7Y2Cmnr6xwdHqIvlh7liMY\nljugBWaHLKyWSOfLHRk2YPT/d2dXdabqDvleF9rTLaEUd0ugb2pzM9Xi5l1DjQ9H0z1n51Y7HtxH\ndw0SjQjf1eK+I87MG2P/5l2d9f8te4bIlapczupEvp1wZm6VfaOGm8tpQinuZ+Zz7B/rJ9VhB9+8\na4iry0VW9Ok0XZMrVbmyXORIiyqcrUjGo9w4meIFvam3I86a4t6qCmor9MrJHs4u5Lipw77vlnCK\n+9xqx5YLwM2mj/iMtt675mVTXDoVdzCWsmfmdd/vhDPzqwwlY40j9drlVeb35aVZ3f/dUq8rzs7n\nOp5YuyV04l6t1Tm3kO+qg60JQfvdu+fMDsT9pqlBLmUKFCs6maZbLHFpVeK6Fam+GHtHkpxb1KHA\n3TKztEaxUtfi7hTn0wXKtTpHu7Dc9432M5CIar/7Djg7nyMeFW4YH+j4b2+aHqSu4EJaR8x0iyHu\n3W3m3TQ9yMsL2rDpFmvV2UmU2E4InbifNTv45i46OBIRjk4PckZb7l1zdj7H4clU28ljzdw0lQLQ\nAtMl2XyZxVy5q1UTGCunl+dzOmKmSyzdODLlfKQMhFDcu42UsTg8meIVvTTtmrPzq12Ly42Txt9Z\nfntNZ5w1J8UjXVqON02lyJdrzK2U7GxWaDgzn2NqqI+RgXhP7hc6cX95Ice+0c4jZSwOTaa4srym\n/b5dUKzUuJgptF3PZyP9iSj7Rvu15d4l65Zj95Y7wDnd/11xpoebqRBCcT+/mOfwZKrrvz88mUIp\nuKgzJTvmlcU8dQVHutjvsLhperBhgWo64+x8jv64MUF2w42muOvJtXOUUpybz3W9au2G8Il7usAN\nE51v5lkcmjAmBu2a6Ryrz27cweR601SKl+fz2u/bBefTeW6YGNjyWMOt2DXcRyoR5eUFPfY7JZMv\ns1qqcsNE92O/U0Il7kuFMstrlYZAd8MhU5jOa3HvGCvKZSeT6+HJFGuVGgur2u/bKRfS+R2NfRHh\n8FRKh0N2gVUT6dAOxn6nhErcz9sgLiP9ccZTCc7r0r8dczGTZzyVYCjZ/YbSQTOEUhcQ64xaXXEp\ns7ajsQ9ww3iKS7rvO+aCqRc77f9OCJW4Wx18aAduATBmX+2W6ZwL6UJDnLvF+vuLOta9I2ZXipRr\ndQ7uUFwOTgxwOVugpksvd8SFdAER2D+mxd0Rzi8agrBTgTk0mWq8l6Z9LuxwvwOML4eIttw7pWE5\nju/MsDk4PkClpri6rAuIdcLFdIE9w8m2a+jbQajE/UI6z96RnXfw4YkUsytF1so6HLJdytU6V5fX\nuspMbSYRi7B3pJ+L2i3WERdtcEkCjc9Pr5w643w6v+NVU6eEStyNaIGd71Y3NlW1wLTN5WyBuoKD\nNvT/wfEBHYraIefTBeJRYW+XYZAWlkDplVNnXMwUdrSZ3Q1tibuI3CsiL4nIWRF5/yavHxSRr4rI\nMyJySkTeYX9Td86FdIFDkzufPS3rRwtM+1hiYMeG0g0TWtw75WImz/6xAaJdhkFa7BnpJx4V3f8d\nkCtVWcyVvWe5i0gU+ChwH3AMeEhEjm247HeBLyilXgc8CPxvuxu6U1aKFdL5si2W+wFzU0QfXNA+\nDbfADt0yAAfGB1jMlcmVqjt+r7Bgx34HQDQi7B8b0G6ZDrBrv6NT2rHc7wLOKqXOKaXKwOeABzZc\nowDr3K4R4Ip9TbQHK3xrp5upAKMDcVKJKJezeoC3y4V0gf54tOM64pthiZQOyWsPpRQX0wVbJlYw\nJldtubePXfsdndKOuO8DLjU9vmw+18wHgJ8RkcvAY8Av29I6G5kxrexuU6+bETGsF225t8/FTJ6D\n4wMd1xHfDMsC0qV/28PKjrRjvwOM1dcFvd/UNpZL0nNumTZ5CPi0Umo/8A7gL0XkuvcWkYdF5KSI\nnFxYWLDp1u0xs2SK+9jOxR1g/1i/FvcOuJRZ44BNlmMj1j2jBaYdrLF/wKaxf3B8gJVilaVC2Zb3\nCzpXltYYTsYY3kHyXje0I+4zwIGmx/vN55p5N/AFAKXUN4AkMLnxjZRSn1BKHVdKHZ+amuquxV0y\nk10jGY8wkUrY8n6GuGvLsV2uLK2x3yZxGRmIM5yM6cm1Ta6Y4r7TSBmLgzqgoCNmsmvs62HykkU7\n4v4UcFREDotIAmPD9MSGay4CbwUQkVsxxL23pvk2zCytsXe03xa3ABjJNKvFKssFfVj2diyvVVgt\nVdk7mrTtPfeO9jdES7M11iRo1+RquTZ1/7fHzNIa+2wc++2yrbgrparAI8DjwIsYUTGnReSDInK/\nedmvA78oIt8BPgv8nPJY2T6jg+0Z3AAHxo33uqSt922xRGDfqH3Wi3aLtc+VpSIDiSgj/fa4Bazv\n0cxS0Zb3Czp2a0+7tHVihVLqMYyN0ubnHm36/QXgbnubZi8z2TVu2zu8/YVtsr8pHPL2fSO2vW8Q\nWXcL2Gu5f+uVjG3vF2Rmlgq2rlpHB+L0x6Pacm+DlWKF1WLVNpdYJ4QiQ3WtXCOdL9s6e1pLXO13\n3x67N7PBsB5XilVWi9otth1Xloq2jn0RYd9YfyMCTdOaRpSejWO/XUIh7k6Iy0h/nME+vanXDjNL\naySiESZTO49xt9jb8Ptq18B2WPtNdrJ3tJ8runjYtqy7JLW4O8KMAz5fI9Zd+33bYSa7xp7RZNcn\nAG2GNVHPLOmV01aslWtk8mXbNlMt9ukN7baY0eLuLE4tjfaN9jc+PE1rrjiwoaQ39dpjxoH9DoB9\no0kWc2V9UPw2zGTNVeugfavWdgmHuC8ViEaEXTakvjezZzTJrF6absuVpaLtboGpwT7iUdF+321w\nIlIJmt1iuv+3wnCJ2btqbZdwiHt2jd3DSWJRe/+7e0b6yRYq2nrZgnK1ztyq/eIeiQh7RrRrYDuc\ns9z1nkc7OLHf0S7hEHeH4kx3DxtfmNllPcBbMbdSRCnY70D/7x1NarfYNlxZWiMi62PVLvaO6j2P\ndnDCJdku4RD37JojoUh7RowvjI4aaI214eyE9bJvdEBb7tsws+TMqnX3SJKI6D2PrShX68yvlrTl\n7hS1umJutWT7shRgj/mhacu9NVccCEO12DeaZG6lSKVWt/29g8JM1hm3QDwaYddwUk+uW3B1eQ2l\n3IlxhxCI+2KuRK2u2D3inFvmqhb3llhffmuVYyd7RvupK5hfLdn+3kHhyrIzq1YwPlMt7q1xMwwS\nQiDullVtt88RoD8RZXQgri33LbiyvMbkYMKRU9/1nsfW1OuK2eUiexwwbMBwzcyu6L5vhTUunTBs\n2iH44r7inLhb73tV+9xbcnW5yG6HBvcuLe5bkimUqdQUu4edibHeNZxkTvd9Syzt2eWQ9mxH4MV9\nzurgEWcG+N7Rfu2W2YLZ5SK7hhyaWM1JQ1uPm9NYtTo0ue4eTpIv13R9nxbMLRcZ6ouR6murPqPt\nBF7cZ5eLxCJia12TZnaPJLXluAXzqyV2OSQuYwNxErFIYwLXXMucw5ajNWno/t+c2ZWiY2O/HYIv\n7itFpof6HMsQ2zOcJJ3XadibUaoadU2ccomJCLuH9eTaioZL0nG3mN7Q3ozZlZJjY78dAi/ucw7P\nntp6ac38ivGld3KA7x7Wm3qtmFsuImKUanCCxoa27v9NmVsuuuZvhxCI++xy0VFxsSIRtN/9ehob\nSg5OrrtGknpibcHsSpHJwT7bE5gstGHTmlpdsZArsduhvb52CLy4z62UHJ0994zqiI1WOBmGarF7\nuI/Z5SIeO9XREzjtFkjGjaP79Ni/nkZ+jbbcnSFXqpIrVR3zOcK6cOkSBNezvqHnnPWyazhJqVpn\neU1HbGxkfsV5t4B2i22O05vZ7RBoce+F5ZjqizGcjGnrZRPmVor0xSK2Hcy8GTocsjWzK0XH3QLa\nLbY5ToehtkOgxb1Xs+euYT3AN2N2pcTukaRtBzNvhs5S3ZxipcZSoeK4W8Byi2muZc7h5Ml2CLS4\n92r2nB7u0/VNNqEX0QLW++vJ9Vp6ZdjsHk6ymCtR1cXbrmF2pUg0Iky4cAKTRbDFvUez5/RQshH2\np1lnbtXZSCVYFy8drXQtvTJsdo0kqStYyOnx38zsconpoT6iLpzAZBFocZ9bKTKcjNGfsL9oVTPT\nw30srJZ0xEYTShlFq5zcTAVIxCJMDia05b6BXhk22i22OXM92MzejkCL+6yDRauamR5KUq7VWSro\niA2L5bUKpWq9JwN8l85SvY65HuQYgHaLtWJ2xflV63YEWtx7NXtOmwdva7/7Ok6nvjdjhOPpvm9m\nbqVEfzzKkMNFqxrRSnpyvYa5HhmWWxFoce/V7Lku7nqAW/QiDNViejjJvLYcr8EIg3Q2UglgfCBB\nPCp6cm0iX6qyWqpqt4xTVGt1FlZLvXHLmB+i3lRdx+qLXq2cjNrlOmLDYq4H+x0AkYgwNWjsOWkM\nrBV8L/p/KwIr7ul8mbpaF14n0W6Z67FWMVNDzg/wqaE+lIJMvuz4vfzCQq7EtEN19DcyNZzUq9Ym\nrImuV/3fisCK+3oHOy8uqb4YqURUbyo1sbBaYqQ/7sjxehtpTK565dRgYbXUk4kV0Jb7Bqy+6FX/\ntyLw4t6rDt41nNQDvImFXA/FxbzPQk5PrmD4fAvlWk/7f1HHuTdY6OGqdSvaEncRuVdEXhKRsyLy\n/hbX/JSIvCAip0XkM/Y2s3Ma4t6jDLGpoT69NG1iYbXU07637qnp/difHuojnS/rLFWThVyJWEQY\ndbCmUjtsK+4iEgU+CtwHHAMeEpFjG645Cvw2cLdS6jbgfQ60tSOsjLlezZ7Tw0ntc29iYbXEZI8t\nd+2WMej12Lf2PNJ6zwMwx/6gc6e/tUs7lvtdwFml1DmlVBn4HPDAhmt+EfioUioLoJSat7eZnbOw\nWmIoGeuJzxcM62V+RWepWvTScu+LGXXFdQq8Qa9dknrldC293O/YinbEfR9wqenxZfO5Zm4GbhaR\nfxeRJ0Xk3s3eSEQeFpGTInJyYWGhuxa3SS99vmCI+1qlRq5U7dk9vUq+VCXfQ58vGAKjxcXALXHX\nbkmDeR+JezvEgKPAm4GHgD8VkdGNFymlPqGUOq6UOj41NWXTrTenl5YjGPVlwMgMDDuLPXYLgLly\n0uIOGGM/GhHGBhI9ud+0ttyvodfa04p2xH0GOND0eL/5XDOXgRNKqYpS6hXgexhi7xqLPZ49d5kx\nrdp6cScUTFvu6yyslphIJXpWkXByUIu7Ra2uSOfLvrHcnwKOishhEUkADwInNlzzdxhWOyIyieGm\nOWdjOzum134vy3LXA7z30RrWvXRlToNeuyST8SjDyZheOQHZQplaXflD3JVSVeAR4HHgReALSqnT\nIvJBEbnfvOxxIC0iLwBfBX5TKZV2qtHbsVausVqq9thy1CUILHodrQHG5Kr3PAzc2NCb1nkegHcS\nmMDwlW+LUuox4LENzz3a9LsCfs38cR3L5zvZQ8txOBmjLxbRbhmMAR4RGE/1xucL10ZsDCXdjS92\nm4XVErfsHurpPXWWqoGXxD2QGarzLnSwiDA93Kc3VDF9voO9PYVmajDZuHeYqdcViz12y4CVxBfu\nvgd3XJKtCKS4u9XB00O6gBK4Ey1g7XmEXWCW1ipUXfD5Tg/pPQ9wxyXZimCKe653RcOamRxMkM7p\nLL1eb+jB+kQedsvdLbfAlJnnkS/Xenpfr7GwWmIgESXl8CEp7RBMcV8tIT32+YLh49cFlNzZ0Bsd\niBOPSuizVN0qN7teAiLcK1evZKdCgMV9IpUgFu3tf29ysI9soRLqQyPc8vmKGIdGhD1ayaqM6Ybl\nDnrl5JUEJgiouC/mSj2NlLGYHDRWCtkQF1BaXqtQqSlXBvjUUJ+23F1yy1grhdD3vwuGTSsCKe5u\nLY0amXohHuBubihNDelY64VV42DsVKI3BfMsdGVOA+2WcRjXxN2852KIN1XdjPM1ShBon+/UUJ/j\nB2NvZLQ/TiwS7j2PUrXG8lpFu2WcQinl2tLIstzTIR7glri74Rab0odGuDb2IxEJfX0fy6jTlrtD\nrBSrlKt1V2bPCdPnHuaIGVct98EESkG2UOn5vb2Cmxt6kyHPUrUihbS4O4Sb4jLUFyMRi4TaLbOY\nK5GIRRhO9j7Od8JaOeXDKzBu+nwnBhOh73vQ4u4YbnawFY63GGLrxbIce+3zBZgw8xrCmkhWrtbJ\nFiquicvkYF9o+x68lZ0KQRR3l7JTLSYHEyyGOBTSzVAwy3IPq1vMsppdtdxz5dCWILAMy4mUFndH\ncHNDDwyBCb3l7pK4TA2GO1rJ7bE/meqjXKuzGtKyy4u5EqMDcRIxb8iqN1phIwurJeJRYaTfnbKv\nk4OJ0FqOYCWQ9bbsg8Vwf4xYREIbrWS5RNzqfyugIKyumUy+3HANeoHAiXs6V2Ii5Y7PFwyrKZMv\nU6+Hb2laqysy+bJrlqOINFwDYcQyKtxyC0yEPBR4MVdu9IEXCJy4Z/LlhgXhBhODfVTriuW18IXj\nLRXK1BWuWi8Tqb7QRmxkzL0et8a/9bmHdeWadnHVuhmBE/fFvLuz52SIY93TpriMu9j/E4MJFkJq\nuafzZZLxCAM9Lj1gMRXyDO10vuyZzVQIoLincyUmXbQcp0JcX6bh83W5/8PrFnDXJTk2EF6fe6VW\nZ6lQcdVrsJHAiXsmX+55HfdmrPoyYRzgljvEzZVTmH3u6VzZVbdAIhZhpD8eSrdYtmC6xPSGqjMU\nylUK5Zq74hJiv6Mlqm7veaxVahTK4QvHS+dLrho2EN7JdX3sa7eMI3hBXMYGEkQjElJxN07Aspbn\nbhDmLNWMB6I1JlPhPI2soT3acneGRrSAix0ciQjjqZBaL/lyY3Jzi7DW1FdKmcEEHrDcQ5ih7QWX\n5EYCJe5e6eCwnqWazrmfxLFedjlcApMrGdVQJ12O1gjr2F90OYFsMwIl7oseWRpNhjQcL50vecJy\nhPAl0njBJWndfymE5whn8iWiEWE46U5m/GYESty9MsAnQ1pfxrDc3bUcrQ3FsLkGGjkGrm+oGp9/\n2M4RTueMKL2Iiy7JjQRK3DN54/zIgUTva4k3M2nWtQ5bdby0B3y+yXiUob5Y6FwD1krFrdIPFpON\naLFwifuiB1ySGwmUuKdz7osLGNZLsVInX6653ZSeUa7WWV6ruG65g7FyC5u4pF0uPWAR1gNT0vmS\n6xPrRgIl7m6XHrCwPuQwuWYaSRwemFwnQ5ilav1/3XfLhDMU1SuGZTOBEvdMvuSJpVFjgIfI7+il\nON8wJtIs5soM9cXoi7lTV8aiYdiEcHL1wqq1mbbEXUTuFZGXROSsiLx/i+t+XESUiBy3r4nt44VQ\nPFgXuEyYxN0jYahWG8LmFnC7GqrFcDJGPCqhcosVKzXy5Zon+r+ZbcVdRKLAR4H7gGPAQyJybJPr\nhoBfAb5pdyPbQSllLo3cF5fxhriHR2C8EqkExqZeJl+mFqKa+kYYqvtjX0SMssshstzTHkie3Ix2\nLPe7gLNKqXNKqTLwOeCBTa77PeBDQNHG9rXNaqlKuVb3RBKBtTwLk1vGWoa7nUQDhuVeV0Z9+bDg\nlVUrhC9L1ZrIvDC5NtOOuO8DLjU9vmw+10BE7gQOKKX+0ca2dUQm5404X4D+RJT+eLTRpjCQyZeJ\nRYThfnfDUKHZ7xue/vfSKUATIdvQ9tKqtZkdb6iKSAT4CPDrbVz7sIicFJGTCwsLO731NXjJ5wvG\nJBMqn7uZxOFWLfFmwpalWq8rsgXvWO6TqXCFonpp1dpMO+I+AxxoerzffM5iCLgd+BcROQ+8ATix\n2aaqUuoTSqnjSqnjU1NT3bd6E7xSesAidEtTj/h8oek0rJD0//JahVpdecZynBzqC1USX6ZxApk3\n+t+iHXF/CjgqIodFJAE8CJywXlRKLSulJpVSh5RSh4AngfuVUicdaXEL3D4/ciNjA+Gy3BddPiii\nGWvPIyx5Bl5btU6kEqFK4kvny/TFIqRcOt6wFduKu1KqCjwCPA68CHxBKXVaRD4oIvc73cB28UoS\nh8VEyNwybp+A1cxIf5xoRBqJVUGnUZHQI/3fqO8TErfYYs7ITvWCS7KZtna/lFKPAY9teO7RFte+\neefN6pypYPLQAAAUgklEQVTFXJmhpPtJHBbjqUSoYq29lMQRiQhjA/HQuMW8dgqQtXrO5MvcMJFy\nuTXO48XsVAhQhmo6X/ZUbYfxQWNpGobj3tbK3kviGBtIhCZaycqn8MrKadyc5MOyck17JDN+I4ER\n94wHzo9sJkzHvVkrFK/43MGMVgqRW8Y43tAbtcTHB8KVoe2F4w03IzDi7qUkDgiX9dKIFvCIWwYM\n10AY+h6MyXVsIEEs6o2v8/hgeMS9cbyhh7THwhujwQa8lMQBzSUIgj/AvZjEMTaQCM2BEV4zbFKJ\nKIloJBQrJ+t4Qy+NfYtAiLuVxOElt0CYiod5MYljIpUgWyhTD0F9mbSHIpXAqC8zngrHnsd6NVTv\njH2LQIi7lcThpQEepqWpVw6KaGYslTDqy6xV3G6K46Rz3jsoYiwkocDrOQbeGfsWgRB3ryVxAAz1\nGaVPwxCOlzGTOAY8lMQRKreYxyx3MPM8QuCWSXuoptVGAiHuXis9AE1L0xDEui/mjFAwLyVxTIRk\nQ7taM4839JjlGJbaSlainBZ3h8h60C0ARvRIKAZ4vuy5uhpjKSMsMOiT6/JaBaW8Jy5hEfdM3nD7\nea3/ISDibi3/rPharzCRCkfxsEyhwpjn+t6y3IPtc7cE1Gv9P55KsFo0IkmCTLZQJhmPMJBwv9T1\nRgIh7pblPurBAR4G6yXrQZ9vWCz39RwDr/W/0Z6gH5iSyZc9Z1RaBELcM/kKQ30xEjFv/XfCEg7m\nRXHvi0UZ7IsFfuVk+Xy9Zrk3MrSD3v/5cmMi8xreUsMuyeRLnuzg8VSC1VKVUjW4pU/L1Tqrpaon\nrZfxVPATmbzq87Umm6D3vxcjlSyCIe6FimfFHSAbYL+vtez2Yv+PhWDPo2G5p7xRV8aicRpWCPrf\na6smi0CIezZfZtwjRZOaWV+aBtfvm/FwKJiVpRpkMvkyg33eKXVtEZY8Ay+dY7CRQIh7xqN+rzAM\ncGtPwYvWSxj2PAyfr/cMm9F+a0M7uP1fqdVZLVY9OfYhIOKeLXhzx3oiBCUIvGy5j5tumSCf5Znx\n6NiPRSOMDsQDPfYbCUwey/Gw8L24Fys1CuWaJzvYKoEb5Jru1oaZF63H8VSCUrXOWiW4G9peXbWC\nkXcS5BIE1l6aFydXCIC4Zz2awATG0jQiAbfczQHuxaWpNSaCPLl6Oc466G6xjIcNGwiAuK93sPcG\nuHGWZ8Ctl4Jxdm3cIwdFNBOGPQ8vx1mPBXxD28t1ZSAA4p71aJyvRRisFy8VbGumUXY5oAJTrBhn\n13p17Ae9/Ib1f/Pqysn34p7xaIaeRdBLEGQL3rUcG2d5BnRyXSp41yUGpuUe4A1tr5Y9sfC/uOe8\ndfL7RiYGE8GOc/eyzzfg0UrrdWW86fOdSCWo1hUrxarbTXGETL7sybInFt5sVQdkChVEYKTfmwM8\n8Ja7h32+1oEpQXXLrPt8vXNITTNB3/PIFrxX6roZ34t7Nl9mtD9ONOKdgyKaGU/1sWQeAxg0lFKe\nrq0hYm5oB9Qt43XLfSzg4p7Je7f0AARA3DMe9vmCsTRVikBGDaxVapSqdU8P8PEAb+p5tZa7RdAP\nic8WvGvYQADEPethny8E23rxuuUIZmXIAE6sYPS/l12SQa8Mmc1775CaZnwv7l7O0IOm4mEBdA1k\nPZzAZBHkPY9socxIf5yYB3MMIPiVIY2iYd6cWCEA4p4teDfOGtY3lYJ4Io2X68pYTARY3L0cqQTQ\nH4/SF4sEcuW0Vq6xVql52rD0tbgrpYylkYc7uBExEMABns17X9zHUgmW1ypUasE7y9PLOQZgbGhP\npBKBXLV69dzmZnwt7vlyjXKt7ukOHjXrzAfR7+jV8zubsVZ1QbQeMx73+YKRaxDEc2zXC+Z5t//b\nEncRuVdEXhKRsyLy/k1e/zUReUFETonIV0TkBvubej1+6GDrLM9MAE9jyhbKRASGk971O1pjI4in\nYWU9XPrBwqitFLy+twwbL/f/tuIuIlHgo8B9wDHgIRE5tuGyZ4DjSqk7gL8G/sDuhm6GH6I1wKga\nF0zL0YjzjXg0xwDWVxVByxJWSnk+DBisPY9g9T00H2/o3f5vx3K/CzirlDqnlCoDnwMeaL5AKfVV\npVTBfPgksN/eZm6O1+N8LcYHgrmp53WfL8CEmb0ZtP4vlGuUq3UfGDaJQK6aGoalh7WnHXHfB1xq\nenzZfK4V7wa+vNkLIvKwiJwUkZMLCwvtt7IFfvD5QnBLn6Zz3o7WgPVa20Hb8/CLYTORSpArVSlV\ng3VgSjZvuiQ9mmMANm+oisjPAMeBD2/2ulLqE0qp40qp41NTUzu+nx+WRhB0y927gxuaDuwIWP/7\nybCB4O15ZAplRgcSni17Au2J+wxwoOnxfvO5axCRe4DfAe5XSvXEyZbJl4lFhKG+WC9u1zVW6dOg\nkclXPC8usWiE4WQscP2f8YlhMxHQPQ8jO9Xbhk074v4UcFREDotIAngQONF8gYi8DvgTDGGft7+Z\nm2P5fEW8O3uCYV3lyzWKATrLUyll9L/H3QIAE4N9gYvYyPrA5wvrFSuDtnJN50ueN2y2FXelVBV4\nBHgceBH4glLqtIh8UETuNy/7MDAIfFFEnhWREy3ezla8fApQM5YALgVIYFaKVWp15fkBDjA2EA9c\nxIaXj5dsxtrwDZq4Z/2wam3nIqXUY8BjG557tOn3e2xuV1t4vXCPRfMA3z2SdLk19uCH7FSL8VSC\nmaWi282wlWzBcEkOJ73tkgyq5Z4plLkzNep2M7bE1xmqGY+X3LRoVMcLUMSMX3y+YFaGDJq4mGU3\nvO6SHOmPIxKsaCWj7In3XZK+FnfjFCBvb2pAME+k8YvPF4wJKBOwszy9XuraIhoxDkwJUrTSaqlK\n1QcuSd+Ke71ubOj5YYA3wsGCZLn7yS0zkKBcq5MvB2dDO+MTwwaCV3Y565McA9+K+0qxQl35wy0w\n2h+8TSW/bOhBc6x1gPrfJy5JCF6eh18MG9+Ke9onHQxGrPVIfzxw4pKIRkglom43ZVuCeNybH3y+\nFoGz3H2y3+RbcffL0sgiaGd5WvsdXt/Qg+AdddhwSXpcXCyMsr/B6HtYP1XN6y5h34q7X5ZGFmMD\nwaoM6Yda4hbWlzAoAtNwSfqo/7OFMvV6MDa0re/x+KC3+9+34u6XpZGFsTQNThKT3yxHCM6Gtt8M\nm/FUgrqC5bVgjP9MvuILl6Rvxd0SSq8vjSzGBoIVa53N+0fch/pixCISGLdYw3L0Sf8H7aBsv7gk\nfSvu2UKZ/niUfo/PnhbjqQSZQnBirf0UrSEigSre1jBsfNL/4wELBc74pKaSb8U94yPLEQz3Ubla\npxCAWOtqrc7ymn987mCdCBQQcTHr5PjFJWmNk6AclO2XVatvxd0v2akWQdrUW16roJR/LEcwz/IM\nQN+D/1ySllsmMP3vgxPIwMfi7pelkUWQslT9tpkN626xIJAtlEnGI75xSY41DJtgVOb0S+kH34q7\nX5ZGFkEqfeqXON9mglQ8LOMTcbFIxqOkEtFARIvV6oqlNe+X+wUfi3vaRxl6EKzKkOuWu3/cYmOp\nBEtrFWoBiLU2XJL+GftgJTL533JfKpR945L0pbhXanVWi1VfdLDFemVI/1svfovWABgfiKOU8eX0\nO36KVLIYT/UFIhTSTy5JX4q7nzrYYjgZJxKQutaN/vfTyilIex4+c0mCMbkGoe/9tJntT3H3UQdb\nRMy61kHY1MvkywwkoiTj/tjQA5honAgUhJWTv1ySYFjumQCEQq5XQ/W+S9KX4u6nDm4mKJt6frQc\nxxob2v72+1ZqdVZ85pIEIxwyHYADU/yUHexLcbc62LLG/MJYQBJp/OnzDcaehx9dkmC48ErVOmsV\nfyfxNQxLH6ycfCnuvrXczep4fsdPtcQtghKt5EeXJKzX1Pd7lmo2XyblE5ekL8Xdb7XcLcYCUhnS\nj5a7FWvtd3HxrWETkJr6GR+FofpS3DOFMkPJGPGov5o/njIiBvzud8zk/Ge5gzG5+t5y95HPtxmr\n7LLfAwr8ZNj4Sx1N/LihB8ZKo1ZXrBSrbjela4qVGvlyrZFx6yeCUDysUcvdZ5Nro7aSz1dOfnJJ\n+lLcMwV/VSS0aJQ+9bHALBUMt5JflqbNBGFD2xo7oz4b/+MBKR6mLXeHyeRLvungZhpnefp4aepX\nyxGMNgdBXIaSMRIxf311h/pixKPi67EPxoa2XwxLf40QEz91cDOWIPrZcvdrKB6YeQa+Fxf/WI7N\niJhJfD52y5SqNXKlqm9ckr4Ud+OgDn90cDNBiBiw2j7hQ4EZSyUolGsUfRxr7beCec2MpxK+ri/j\nN5ek78R9rVxjrVLzTQc3E4T6Jn633MHfk6ufDibfyMSgv1dOfjNsfCfujVAwH1ovqUSURDTi61h3\na4CP9vtv5bR+aIR/BcavLknw/2lYfsuvaUvcReReEXlJRM6KyPs3eb1PRD5vvv5NETlkd0MtGht6\nPpk9mzEOao772+eeLzPSHyfmsxwDWD/uze/Wox9dkmBYvOmcf2v7pH2mPdt+Q0UkCnwUuA84Bjwk\nIsc2XPZuIKuUOgL8T+BDdjfUwq9JHBZ+rwyZKfjjFJrN8Lvl7meXJBiVIVeKVSq1uttN6Qq/uSTb\nMb/uAs4qpc4ppcrA54AHNlzzAPDn5u9/DbxVRMS+Zq6znn7tjw7eiN8rQ2byJcYG/Gk5+t3nvmha\nvX50ScJ6rLtfV04z2TUSsYhv3DKxNq7ZB1xqenwZ+P5W1yilqiKyDEwAi3Y0spmsj+OswZiUnnhh\njh/5yL+63ZSOUcDZ+Rw/9tq9bjelK0b7jQNT/vifz/KZb150uzkdUygbUT637xtxuSXdYX1nf/Lj\n3yDhQ7fe7HKRm6YGiUYcsVttpx1xtw0ReRh4GODgwYNdvcfe0X7edmwXwz7c0AN46PUHfV1b5tY9\nwzzylqNuN6MrIhHhV++5mRdnV9xuSte87bZdHNsz7HYzuuINN47zzjv3+TYU9eiuQe67fY/bzWgb\n2U5oROSNwAeUUm83H/82gFLqfzRd87h5zTdEJAbMAlNqizc/fvy4OnnypA3/BY1GowkPIvK0Uur4\ndte1szZ6CjgqIodFJAE8CJzYcM0J4F3m7z8B/PNWwq7RaDQaZ9nWLWP60B8BHgeiwKeUUqdF5IPA\nSaXUCeCTwF+KyFkggzEBaDQajcYl2vK5K6UeAx7b8NyjTb8XgZ+0t2kajUaj6Rb/bVlrNBqNZlu0\nuGs0Gk0A0eKu0Wg0AUSLu0aj0QQQLe4ajUYTQLZNYnLsxiILwIUu/3wSB0ob2IBuV2fodnWOV9um\n29UZO2nXDUqpqe0uck3cd4KInGwnQ6vX6HZ1hm5X53i1bbpdndGLdmm3jEaj0QQQLe4ajUYTQPwq\n7p9wuwEt0O3qDN2uzvFq23S7OsPxdvnS567RaDSarfGr5a7RaDSaLfC0uHvpYO6mex4Qka+KyAsi\nclpEfmWTa94sIssi8qz58+hm7+VA286LyHPmPa8rli8Gf2T21ykRubMHbXpVUz88KyIrIvK+Ddf0\nrL9E5FMiMi8izzc9Ny4iT4jIGfPfsRZ/+y7zmjMi8q7NrrGxTR8Wke+an9Pfishoi7/d8jN3qG0f\nEJGZps/rHS3+dsvvrwPt+nxTm86LyLMt/taRPmulDa6NL6WUJ38wygu/DNwIJIDvAMc2XPNLwMfN\n3x8EPt+Ddu0B7jR/HwK+t0m73gx8yYU+Ow9MbvH6O4AvAwK8AfimC5/pLEacriv9BbwJuBN4vum5\nPwDeb/7+fuBDm/zdOHDO/HfM/H3MwTa9DYiZv39osza185k71LYPAL/Rxme95ffX7nZteP0PgUd7\n2WettMGt8eVly91TB3NbKKWuKqW+bf6+CryIcYasH3gA+Atl8CQwKiK9PDfsrcDLSqluk9d2jFLq\naxhnDjTTPI7+HPixTf707cATSqmMUioLPAHc61SblFL/pJSqmg+fBPbbca9OadFf7dDO99eRdpka\n8FPAZ+26X5ttaqUNrowvL4v7ZgdzbxTRaw7mBqyDuXuC6QZ6HfDNTV5+o4h8R0S+LCK39ahJCvgn\nEXlajPNqN9JOnzrJg7T+wrnRXxa7lFJXzd9ngV2bXONm3/0CxoprM7b7zJ3iEdNl9KkWbgY3++sH\ngTml1JkWrzveZxu0wZXx5WVx9zQiMgj8DfA+pdTGE5e/jeF6eA3wv4C/61GzfkApdSdwH/BeEXlT\nj+67LWIc0Xg/8MVNXnarv65DGWtkz4SQicjvAFXgr1pc4sZn/jHgJuC1wFUMF4iXeIitrXZH+2wr\nbejl+PKyuM8AB5oe7zef2/QaMQ7mHgHSTjdMROIYH95fKaX+z8bXlVIrSqmc+ftjQFxEJp1ul1Jq\nxvx3HvhbjKVxM+30qVPcB3xbKTW38QW3+quJOcs9Zf47v8k1Pe87Efk54D8D/80Uheto4zO3HaXU\nnFKqppSqA3/a4p6ujDVTB94JfL7VNU72WQttcGV8eVncPXkwt+nP+yTwolLqIy2u2W35/kXkLox+\ndnTSEZGUiAxZv2NsyD2/4bITwM+KwRuA5ablotO0tKbc6K8NNI+jdwF/v8k1jwNvE5Ex0w3xNvM5\nRxCRe4HfAu5XShVaXNPOZ+5E25r3af5ri3u28/11gnuA7yqlLm/2opN9toU2uDO+7N4xtvMHI7rj\nexi77r9jPvdBjAEPkMRY5p8FvgXc2IM2/QDGsuoU8Kz58w7gPcB7zGseAU5jRAg8CfynHrTrRvN+\n3zHvbfVXc7sE+KjZn88Bx3v0OaYwxHqk6TlX+gtjgrkKVDD8mu/G2Kf5CnAG+H/AuHntceDPmv72\nF8yxdhb4eYfbdBbDB2uNMSsqbC/w2FafeQ/66y/N8XMKQ7j2bGyb+fi676+T7TKf/7Q1rpqu7Umf\nbaENrowvnaGq0Wg0AcTLbhmNRqPRdIkWd41GowkgWtw1Go0mgGhx12g0mgCixV2j0WgCiBZ3jUaj\nCSBa3DUajSaAaHHXaDSaAPL/Ada2BWAQ3CzzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10dab6390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(pdf[:, 0], pdf[:, 1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "samples = list(sample_from_pdf(pdf, 100000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE8FJREFUeJzt3X+spNV93/H3J/xwKhsZCFuKYZ3FzqYRjhSMVuDUbkRD\nw69EWVyl7qIo3jiom6ig2mqqZu1IseuUCtraVl25ROuw8jpyDTS2yypeB68JkuU/+LFQDCzY4RqD\n2NUCa0PAllW3i7/9Y87i8eX+mLt3Zu69e94vaTTPnOc8M2eemTufPec8z7OpKiRJ/fmplW6AJGll\nGACS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTp240g1YyBlnnFEbNmxY6WZI0ppy\n//33f6eq1i1Wb1UHwIYNG9i3b99KN0OS1pQkT41SzyEgSeqUASBJnTIAJKlTBoAkdcoAkKROGQCS\n1CkDQJI6ZQBIUqcMAEnq1Ko+E3g12rD9i3OWP3nDr0+5JZK0PIv2AJL8dJJ7k3w9yf4k/76Vn5vk\nniQzSW5NcnIrf017PNPWbxh6rve38m8muWxSb0qStLhRegA/BH61qr6f5CTga0m+BPwb4GNVdUuS\nPwOuAW5q9y9U1c8l2QLcCPyLJOcBW4C3AG8AvpLk56vq5Qm8Lx3H7IVJ47FoD6AGvt8entRuBfwq\n8JetfBdwVVve3B7T1l+SJK38lqr6YVV9G5gBLhzLu5AkLdlIcwBJTgDuB34O+ATwLeDvqupIq3IA\nOLstnw08DVBVR5K8CPxMK7976GmHt5G0Bgz3vuxxrX0jHQVUVS9X1fnAOQz+1f4Lk2pQkm1J9iXZ\nd/jw4Um9jCR1b0lHAVXV3yW5C/hl4NQkJ7ZewDnAwVbtILAeOJDkROD1wHeHyo8a3mb4NXYAOwA2\nbdpUS3s7kzHfmLMkrWWjHAW0LsmpbfnvAb8GPAbcBfxWq7YVuL0t726Paev/pqqqlW9pRwmdC2wE\n7h3XG5E0XRu2f/GVm9amUXoAZwG72jzATwG3VdVfJXkUuCXJfwD+N3Bzq38z8BdJZoDnGRz5Q1Xt\nT3Ib8ChwBLjWI4AkaeUsGgBV9RDw1jnKn2COo3iq6v8A/3ye57oeuH7pzZQkjZuXgpCkTnkpCEkL\ncoz/+GUAaE3wR0gaPwNgTDxBRtJa4xyAJHXKAJCkTjkEpOPG7HkCh+KkhRkAWrWc+JUmywCQtGwe\nBLE2GQDz8F+fko53TgJLUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTnkYqKRX8TDoPtgDkKRO2QOY\nAM+KlLQW2AOQpE7ZA9Cq4tizND32ACSpU/YAJI2Vc2Brhz0ASeqUASBJnVo0AJKsT3JXkkeT7E/y\n3lb+oSQHkzzYblcObfP+JDNJvpnksqHyy1vZTJLtk3lLkqRRjDIHcAT4w6p6IMkpwP1J9rZ1H6uq\n/zJcOcl5wBbgLcAbgK8k+fm2+hPArwEHgPuS7K6qR8fxRiRJS7NoAFTVIeBQW/5ekseAsxfYZDNw\nS1X9EPh2khngwrZupqqeAEhyS6trAEjSCljSUUBJNgBvBe4B3g5cl+TdwD4GvYQXGITD3UObHeDH\ngfH0rPKLjqnVknQcmvYRVCNPAid5HfA54H1V9RJwE/Bm4HwGPYSPjKNBSbYl2Zdk3+HDh8fxlJKk\nOYwUAElOYvDj/5mq+jxAVT1bVS9X1Y+AT/LjYZ6DwPqhzc9pZfOV/4Sq2lFVm6pq07p165b6fiRJ\nIxrlKKAANwOPVdVHh8rPGqr2TuCRtrwb2JLkNUnOBTYC9wL3ARuTnJvkZAYTxbvH8zYkSUs1yhzA\n24HfAR5O8mAr+wBwdZLzgQKeBH4foKr2J7mNweTuEeDaqnoZIMl1wB3ACcDOqto/xvciaRm8DlN/\nRjkK6GtA5li1Z4Ftrgeun6N8z0LbSZKmxzOBJalTBoAkdcoAkKROeTnoIU6CSeqJPQBJ6pQBIEmd\nMgAkqVPOAei45X9NKC3MHoAkdcoAkKROGQCS1CnnACRpBa3k+UcGgKSJcSJ+dXMISJI6ZQBIUqcc\nApowu8CSVisDQCvOi/BJK8MhIEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnFg2A\nJOuT3JXk0ST7k7y3lZ+eZG+Sx9v9aa08ST6eZCbJQ0kuGHqura3+40m2Tu5tSZIWM8qZwEeAP6yq\nB5KcAtyfZC/wu8CdVXVDku3AduCPgCuAje12EXATcFGS04EPApuAas+zu6peGPebWgrPQu2Dl+SQ\nXm3RHkBVHaqqB9ry94DHgLOBzcCuVm0XcFVb3gx8ugbuBk5NchZwGbC3qp5vP/p7gcvH+m4kSSNb\n0hxAkg3AW4F7gDOr6lBb9QxwZls+G3h6aLMDrWy+8tmvsS3JviT7Dh8+vJTmSZKWYOSLwSV5HfA5\n4H1V9VKSV9ZVVSWpcTSoqnYAOwA2bdo0lueUNDeHQPs2Ug8gyUkMfvw/U1Wfb8XPtqEd2v1zrfwg\nsH5o83Na2XzlkqQVMMpRQAFuBh6rqo8OrdoNHD2SZytw+1D5u9vRQG8DXmxDRXcAlyY5rR0xdGkr\nkyStgFGGgN4O/A7wcJIHW9kHgBuA25JcAzwFvKut2wNcCcwAPwDeA1BVzyf5U+C+Vu/DVfX8WN6F\nJGnJFg2AqvoakHlWXzJH/QKunee5dgI7l9JASdJkeCawJHXKAJCkThkAktQpA0CSOjXyiWCSpPFY\nLSfgGQBT5AXJJK0mDgFJUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpzwPQilgtJ8JIPbMH\nIEmdMgAkqVMGgCR1yjkASVPhtbBWH3sAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVOLHgaaZCfw\nG8BzVfWLrexDwL8EDrdqH6iqPW3d+4FrgJeBf11Vd7Tyy4H/CpwA/HlV3TDetyJpFF6GQ0eN0gP4\nFHD5HOUfq6rz2+3oj/95wBbgLW2b/57khCQnAJ8ArgDOA65udSVJK2TRHkBVfTXJhhGfbzNwS1X9\nEPh2khngwrZupqqeAEhyS6v76JJbLEkai+XMAVyX5KEkO5Oc1srOBp4eqnOglc1XLklaIccaADcB\nbwbOBw4BHxlXg5JsS7Ivyb7Dhw8vvoG0RBu2f/GVm9SzYwqAqnq2ql6uqh8Bn+THwzwHgfVDVc9p\nZfOVz/XcO6pqU1VtWrdu3bE0T5I0gmMKgCRnDT18J/BIW94NbEnymiTnAhuBe4H7gI1Jzk1yMoOJ\n4t3H3mxJ0nKNchjoZ4GLgTOSHAA+CFyc5HyggCeB3weoqv1JbmMwuXsEuLaqXm7Pcx1wB4PDQHdW\n1f6xvxtJ0shGOQro6jmKb16g/vXA9XOU7wH2LKl1E+LYryR5JrAkdcsAkKROGQCS1CkDQJI6ZQBI\nUqcMAEnqlAEgSZ1a9DwASdLyrcbzj+wBSFKnDABJ6pQBIEmdcg5A0tQNj4c/ecOvr2BL+mYArBD/\nACStNIeAJKlTBoAkdcoAkKROGQCS1CkngTU1q/FMSKln9gAkqVMGgCR1ygCQpE4ZAJLUKQNAkjpl\nAEhSpwwASerUogGQZGeS55I8MlR2epK9SR5v96e18iT5eJKZJA8luWBom62t/uNJtk7m7UiSRjVK\nD+BTwOWzyrYDd1bVRuDO9hjgCmBju20DboJBYAAfBC4CLgQ+eDQ0JEkrY9EAqKqvAs/PKt4M7GrL\nu4Crhso/XQN3A6cmOQu4DNhbVc9X1QvAXl4dKpKkKTrWOYAzq+pQW34GOLMtnw08PVTvQCubr/xV\nkmxLsi/JvsOHDx9j8yRJi1n2JHBVFVBjaMvR59tRVZuqatO6devG9bSSpFmONQCebUM7tPvnWvlB\nYP1QvXNa2XzlkqQVcqxXA90NbAVuaPe3D5Vfl+QWBhO+L1bVoSR3AP9xaOL3UuD9x95saTx6+K85\nvQqr5rNoACT5LHAxcEaSAwyO5rkBuC3JNcBTwLta9T3AlcAM8APgPQBV9XySPwXua/U+XFWzJ5Yl\nSVO0aABU1dXzrLpkjroFXDvP8+wEdi6pdZKkienmP4SxGyxJP8lLQUhSpwwASeqUASBJnTIAJKlT\nBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqVDeXgpCkaVvtl6CxByBJnTIAJKlTBoAkdco5gFWg\nh/+VStLqYwBoolb7JJjUM4eAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnq1LICIMmTSR5O\n8mCSfa3s9CR7kzze7k9r5Uny8SQzSR5KcsE43oAk6diM40Swf1JV3xl6vB24s6puSLK9Pf4j4Apg\nY7tdBNzU7iV1zDPhV84khoA2A7va8i7gqqHyT9fA3cCpSc6awOtLkkaw3AAo4MtJ7k+yrZWdWVWH\n2vIzwJlt+Wzg6aFtD7QySdIKWO4Q0Duq6mCSvw/sTfKN4ZVVVUlqKU/YgmQbwBvf+MZlNk+SNJ9l\n9QCq6mC7fw74AnAh8OzRoZ12/1yrfhBYP7T5Oa1s9nPuqKpNVbVp3bp1y2meJGkBxxwASV6b5JSj\ny8ClwCPAbmBrq7YVuL0t7wbe3Y4Gehvw4tBQkSRpypYzBHQm8IUkR5/nf1TVXye5D7gtyTXAU8C7\nWv09wJXADPAD4D3LeG1J0jIdcwBU1RPAL81R/l3gkjnKC7j2WF9PkjRengksSZ3yfwSTjkP+T2wa\nhT0ASeqUASBJnTIAJKlTzgFIjRclU2/sAUhSpwwASerUcT0E5KFwkjQ/ewCS1CkDQJI6ZQBIUqcM\nAEnqlAEgSZ0yACSpUwaAJHXquD4PQJKmbS2df2QPQJI6ZQ9glfGCZJKmxR6AJHXKAJCkTjkEpLFb\nS5NgUs/sAUhSpwwASeqUQ0CSVg2PgpuuqfcAklye5JtJZpJsn/brS5IGphoASU4APgFcAZwHXJ3k\nvGm2QZI0MO0ewIXATFU9UVX/F7gF2DzlNkiSmP4cwNnA00OPDwAXTbkN0nHJw2+1VKtuEjjJNmBb\ne/j9JN9cxtOdAXxn+a0au5HalRun0JKftKb31ziNuO/dX0uzpHZN8fu/KvdXblxWu352lErTDoCD\nwPqhx+e0sldU1Q5gxzheLMm+qto0jucaJ9u1NLZraWzX0vTcrmnPAdwHbExybpKTgS3A7im3QZLE\nlHsAVXUkyXXAHcAJwM6q2j/NNkiSBqY+B1BVe4A9U3q5sQwlTYDtWhrbtTS2a2m6bVeqatKvIUla\nhbwWkCR1as0HwGKXlkjymiS3tvX3JNkwhTatT3JXkkeT7E/y3jnqXJzkxSQPttufTLpdQ6/9ZJKH\n2+vum2N9kny87bOHklwwhTb9w6F98WCSl5K8b1adqeyzJDuTPJfkkaGy05PsTfJ4uz9tnm23tjqP\nJ9k6hXb95yTfaJ/TF5KcOs+2C37mE2jXh5IcHPqsrpxn24ldGmaedt061KYnkzw4z7aT3F9z/j6s\nyHesqtbsjcFE8reANwEnA18HzptV518Bf9aWtwC3TqFdZwEXtOVTgL+do10XA3+1QvvtSeCMBdZf\nCXwJCPA24J4V+FyfAX52JfYZ8CvABcAjQ2X/CdjelrcDN86x3enAE+3+tLZ82oTbdSlwYlu+ca52\njfKZT6BdHwL+7Qif84J/v+Nu16z1HwH+ZAX215y/DyvxHVvrPYBRLi2xGdjVlv8SuCRJJtmoqjpU\nVQ+05e8BjzE4C3qt2Ax8ugbuBk5NctYUX/8S4FtV9dQUX/MVVfVV4PlZxcPfo13AVXNsehmwt6qe\nr6oXgL3A5ZNsV1V9uaqOtId3Mzi3Zqrm2V+jmOilYRZqV/sNeBfw2XG93qgW+H2Y+ndsrQfAXJeW\nmP1D+0qd9ofyIvAzU2kd0Iac3grcM8fqX07y9SRfSvKWabUJKODLSe7P4Mzr2UbZr5O0hfn/MFdq\nn51ZVYfa8jPAmXPUWen99nsMem5zWewzn4Tr2tDUznmGM1Zyf/1j4Nmqenye9VPZX7N+H6b+HVvr\nAbCqJXkd8DngfVX10qzVDzAY4vgl4L8B/2uKTXtHVV3A4Kqs1yb5lSm+9oIyOEHwN4H/Ocfqldxn\nr6hBX3xVHT6X5I+BI8Bn5qky7c/8JuDNwPnAIQbDLavJ1Sz8r/+J76+Ffh+m9R1b6wGw6KUlhusk\nORF4PfDdSTcsyUkMPtzPVNXnZ6+vqpeq6vtteQ9wUpIzJt2u9noH2/1zwBcYdMWHjbJfJ+UK4IGq\nenb2ipXcZ8CzR4fB2v1zc9RZkf2W5HeB3wB+u/1wvMoIn/lYVdWzVfVyVf0I+OQ8r7dS++tE4J8B\nt85XZ9L7a57fh6l/x9Z6AIxyaYndwNGZ8t8C/ma+P5JxaeOLNwOPVdVH56nzD47ORSS5kMFnMY1g\nem2SU44uM5hEfGRWtd3AuzPwNuDFoa7ppM37L7OV2mfN8PdoK3D7HHXuAC5Nclob8ri0lU1MksuB\nfwf8ZlX9YJ46o3zm427X8JzRO+d5vZW6NMw/Bb5RVQfmWjnp/bXA78P0v2OTmOWe5o3BESt/y+Bo\ngj9uZR9m8AcB8NMMhhNmgHuBN02hTe9g0H17CHiw3a4E/gD4g1bnOmA/gyMf7gb+0ZT215vaa369\nvf7RfTbctjD4j3u+BTwMbJpS217L4Af99UNlU99nDALoEPD/GIyxXsNg3uhO4HHgK8Dpre4m4M+H\ntv299l2bAd4zhXbNMBgTPvo9O3rE2xuAPQt95hNu11+0785DDH7Yzprdrvb4VX+/k2xXK//U0e/U\nUN1p7q/5fh+m/h3zTGBJ6tRaHwKSJB0jA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE79\nf/HgxUDx7uNBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10db485f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(samples, bins=100);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}