sampling_inversion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inversion Sampling\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.integrate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We wish to compare several sampling methods. We will implement each of them individually as functions. \n",
    "For reproducibility, let's fix the seed of the (pseudo-)random number generator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start with the simples Monte Carlo sampling method: inversion.\n",
    "The idea is to swap the variables $x \\leftrightarrow y$ and draw random values from the cumulative distribution function:\n",
    "$$F(x)=y \\Rightarrow F^{-1}(y)=x.$$\n",
    "Let's implement this as a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inversionSampling(x, y, n):\n",
    "    \"\"\"\n",
    "    Applies the inverse sampling method to a given function.\n",
    "    Here y = f(x) is the probability distribution function to be sampled. \n",
    "    The number of sampling points is n.\n",
    "    Note that even if f is an analytical function, we should bin the data in bins of dx.\n",
    "    \"\"\"\n",
    "    # If the points in y are a PDF, we compute the CDF and then normalise it.\n",
    "    y1 = np.cumsum(y)\n",
    "    y1 /= y1[-1]\n",
    "\n",
    "    # Normally the interpolation would be y0 = (x0, x, y), but here we invert it!\n",
    "    ysample = np.random.uniform(0, 1, n)\n",
    "    xsample = np.interp(ysample, y1, x)\n",
    "    \n",
    "    # Sorting the arrays for plotting\n",
    "    idx = np.argsort(xsample)\n",
    "    xsample = xsample[idx]\n",
    "    ysample = ysample[idx]\n",
    "\n",
    "    return xsample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try out this sampling method using a simple gaussian function. We know exactly what to expect:\n",
    "$$ f(x) = \\dfrac{1}{\\sqrt{2\\pi}} \\exp\\left( -\\dfrac{x^2}{2} \\right).$$\n",
    "Below we define this (expected) function, and the we take $n$ samples from this distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 5000\n",
    "x = np.linspace(-3, 3, 10000)\n",
    "y = np.exp(-x ** 2 / 2) / np.sqrt(2 * np.pi)\n",
    "xs = inversionSampling(x, y, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will plot the sampled points (xs, ys) and the analytic function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11276f2e0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVhV1f7H8ff3MA8CApKzaGlqDpg45FyZQ5qWaY6l3Uqrqw2aSWlqTlfNrCyra2nX/KWpjU6pWVlWTlhcZ8sUE2cRRURAYP3+4OhFBAWFszmH7+t5eJ4z7LPPZ5/sfM9ae621xRiDUkoplRub1QGUUkoVX1oklFJK5UmLhFJKqTxpkVBKKZUnLRJKKaXypEVCKaVUnrRIKMuIyA4RaWN1juxEJElEqlmdoyiIyH9EZIL9dksR2WN1JlX8uVsdQJVcxpjbrM6QkzHG3+oMjmCMWQfcanUOVfxpS0KVKCKiP4yUKgAtEsoyIhIrIm3tt8eKyCIR+VhEztq7oiLtz0WJyGc5XvuWiMyw3w4UkdkickREDonIBBFxsz83QER+EZE3ROQUMFZEbhGRH0XkjIicFJGF2fZrROSWbPv9WEROiMgBERklIrZs+/1ZRKaJSIKI7BeRjlc51hH2bGdFZI+I3G1/vLGIrBeR0/b874iIZ448T4vIn/bXjheRm+2vSbR/Zp72bduISJyIvGw/rlgR6ZtHnjYiEpfjv8ULIrLV/rksFBHvbM+/aM93WEQez/45KdemRUIVJ12AT4EgYAnwjv3xBcC9IhIAYC8ADwHz7c/PBdKBW4AGQDvg8Wz7bQLsA8KAicB4YDVQGqgIvJ1HnreBQKAa0Bp4BHg0x373AKHAVGC2iEjOnYjIrcBgoJExphTQHoi1P50BPG/fxx3A3cDTOXbRAWgINAVeBGYBfYFKQB2gd7Zty9r3VQHoD8yyv39+PGR/r6pAPWCAPX8HYCjQlqzPuHU+96dcgBYJVZz8bIxZYYzJAOYB9QGMMQeA34D77dvdBSQbYzaIyE1AR+A5Y8w5Y8xx4A2gV7b9HjbGvG2MSTfGnAcuAFWA8saYFGPMzzmD2AtRT+AlY8xZY0ws8DrwcLbNDhhjPrDnnQuUA27K5bgyAC+gtoh4GGNijTF/2Y9tizFmgz1bLPBvrvwSnmKMSTTG7AC2A6uNMfuMMWeAb8gqjNm9YoxJNcb8CCwn68s/P2YYYw4bY04BS4EI++MPAR8ZY3YYY5KBV/O5P+UCtEio4uRottvJgHe2cwjz+d8v5j78rxVRBfAAjti7bE6T9UUblm1fB3O8z4uAAJvs3Vr/yCVLKOAJHMj22AGyfqFfkdf+5QlwxYlvY8xe4DlgLHBcRD4VkfIAIlJDRJaJyFERSQQm2d87u2PZbp/P5X7290wwxpzLkbl8LseXm5yf/8X9lufyzzDn56lcmBYJ5SwWA21EpCLwAP8rEgeBVCDUGBNk/wvIMXLqsqWOjTFHjTFPGGPKA4OAd3PpXz/J/1ocF1UGDl1PeGPMfGNMC/v+DDDF/tR7wG6gujEmAHiZrAJ2vUqLiF+2+5WBwzewP4AjZHXLXVTpBvennIgWCeUUjDEngLXAR8B+Y8wu++NHyDq/8LqIBIiIzX5iN89+cxHpYS82AAlkfWln5Hi/DGARMFFESolIFbL65f+voNlF5FYRuUtEvIAUsn79X3y/UkAikCQiNYGnCrr/XLwqIp4i0hLoTFaBvRGLgEdFpJaI+AKjbzihchpaJJQzmU/WydP5OR5/hKyuoZ1kfel/Rtb5gbw0AjaKSBJZJ8ifNcbsz2W7IcA5sk56/2x/3znXkdsLmExW6+QoWV1hL9ufe4Gs7rOzwAfAwtx2UABHyfoMDgOfAE8aY3bfyA6NMd8AM4AfgL3AevtTqTeyX+UcRC86pJRrkKzZ6/9njKl4rW1v8H1qkXUC3csYk16U76Wspy0JpdQ1icgD9i6s0mSdT1mqBaJkcGiREJEO9olEe0UkKpfnB9gnLsXY/x7PbT9KKYcbBJwA/iLrfEphnDtRTsBh3U32ced/APcAccBmoLcxZme2bQYAkcaYwQ4JpZRS6qoc2ZJoDOy1TwJKI2tmbVcHvr9SSqkCcuRiZxW4fBJOHFnLGuT0oIi0IqvV8bwx5oqJOyIyEBgI4Ofn17BmzZpFEFcppVzXli1bThpjylxrO0cWidwmCOXs61oKLDDGpIrIk2QtdXDXFS8yZhZZ69cQGRlpoqOjCzurUkq5NBE5cO2tHNvdFMflMzUrkmMmqDEm3hhzcez1B2QtaqaUUsoijiwSm4HqIlLVvrRxL7ImMl0iItknQHUBdjkwn1JKqRwc1t1kjEkXkcHAKsANmGOM2SEi44BoY8wS4BkR6ULWss+nsC9VrJRSyhpOP+Naz0koVTJcuHCBuLg4UlJSrI7iVLy9valYsSIeHh6XPS4iW4wxkdd6vV7KUSnlFOLi4ihVqhTh4eHkcm0nlQtjDPHx8cTFxVG1atXr2ocuy6GUcgopKSmEhIRogSgAESEkJOSGWl9aJJRSTkMLRMHd6GemRUIppVSetEgopVQ+ubm5ERERQZ06dejRowfJycmXPX7bbbdRv359pk+fTmZmJgBr164lMDCQiIgIIiIiaNu2rZWHUGBaJJRSKp98fHyIiYlh+/bteHp68v7771/2+I4dO/j2229ZsWIFr7766qXXtWzZkpiYGGJiYlizZo1V8a+LFgmllLoOLVu2ZO/evVc8HhYWxqxZs3jnnXdw9ikGoENglVJO6LnnniMmJqZQ9xkREcGbb76Zr23T09P55ptv6NChQ67PV6tWjczMTI4fPw7AunXriIiIAKBHjx6MHDmycEI7gBYJpZTKp/Pnz1/6sm/ZsiWPPfZYnttmb0W0bNmSZcuWFXm+oqBFQinldPL7i7+wXTz3cC379u3Dzc2NsLAwdu1y7iXo9JyEUkoVohMnTvDkk08yePBgl5jXoS0JpZS6QRe7oS5cuIC7uzsPP/wwQ4cOtTpWodAioZRS+ZSUlJTr4xkZGXm+pk2bNrRp06aIEhU97W5SSimVJy0SSiml8qRFQimlVJ60SCillMqTFgmllFJ50iKhlFIqTzoEVinllMKjlhfq/mIndyrU/d2oNm3aMG3aNCIjr3kZaiBrSfJp06YV+vIf2pJQSimVJy0SSimVT+fOnaNTp07Ur1+fOnXqsHDhQsaNG0ejRo2oU6cOAwcOvLSwX5s2bXj++edp1aoVtWrVYvPmzXTr1o3q1aszatQoAGJjY6lZsyb9+/enXr16dO/e/dKFjLJbvXo1d9xxB7fffjs9evS4NKlv5cqV1KxZkxYtWvDFF18UyTFrkVBKqXxauXIl5cuX57///S/bt2+nQ4cODB48mM2bN7N9+3bOnz9/WXePp6cnP/30E08++SRdu3Zl5syZbN++nf/85z/Ex8cDsGfPHgYOHMjWrVsJCAjg3Xffvew9T548yYQJE1izZg2//fYbkZGRTJ8+nZSUFJ544gmWLl3KunXrOHr0aJEcsxYJpZTKp7p167JmzRpGjBjBunXrCAwM5IcffqBJkybUrVuX77//nh07dlzavkuXLpded9ttt1GuXDm8vLyoVq0aBw8eBKBSpUo0b94cgH79+vHzzz9f9p4bNmxg586dNG/enIiICObOncuBAwfYvXs3VatWpXr16ogI/fr1K5Jj1hPXSimVTzVq1GDLli2sWLGCl156iXbt2jFz5kyio6OpVKkSY8eOJSUl5dL2Xl5eANhstku3L95PT08HuGKl2Jz3jTHcc889LFiw4LLHY2JiHLLKrLYklFIqnw4fPoyvry/9+vXjhRde4LfffgMgNDSUpKQkPvvsswLv8++//2b9+vUALFiwgBYtWlz2fNOmTfnll18uXSo1OTmZP/74g5o1a7J//37++uuvS68tCtqSUEo5JSuGrG7bto3hw4djs9nw8PDgvffe46uvvqJu3bqEh4fTqFGjAu+zVq1azJ07l0GDBlG9enWeeuqpy54vU6YM//nPf+jduzepqakATJgwgRo1ajBr1iw6depEaGgoLVq0YPv27YVynNmJs1+oOzIy0kRHR1sdQ7mw7OPxi9tY+pJk165d1KpVy+oYhSo2NpbOnTsXyZd7drl9diKyxRhzzUkY2t2klFIqT9rdpJST0ZaN6wgPDy/yVsSN0paEUsppOHv3uBVu9DPTIqGUcgre3t7Ex8droSgAYwzx8fF4e3tf9z60u0kp5RQqVqxIXFwcJ06csDqKU/H29qZixYrX/XotEkopp+Dh4UHVqlWtjlHiaHeTUkqpPDm0JSEiHYC3ADfgQ2PM5Dy26w4sBhoZY3QShHIpeY1O0lFLqjhyWEtCRNyAmUBHoDbQW0Rq57JdKeAZYKOjsimllMqdI1sSjYG9xph9ACLyKdAV2Jlju/HAVOAFB2ZT6oadP3+ePXv2sG/fPk6ePElaWhqenp6EhYVRtWpVateujYeHh9UxlSoQRxaJCsDBbPfjgCbZNxCRBkAlY8wyEcmzSIjIQGAgQOXKlYsgqlK5y9klFBsby/z581m5ciUbN24kLS0tz9d6eXnRpEkTEt1uwffW5rgH3uSIyErdEEcWidzWtL004FlEbMAbwIBr7cgYMwuYBVlrNxVSPqXyxRhDyoH/0rbtG3z33XcANGzYkGeffZZGjRpx8803c9NNN+Hl5UVKSgrHjh3jjz/+YMuWLXz77bckbJ1Dwg8f4XNzJN/d481dd911xZLPOa/frOcolFUcWSTigErZ7lcEDme7XwqoA6y1/w9TFlgiIl305LUqLlIP7ebUdx+QdmQPu8qXZ/z48fTr14/w8PDLtsvZ4mjYsCG9e/cGoMKTszm3bQ1nY1bStm1b7rzzTqZMmXLV981ZNJRyFEcWic1AdRGpChwCegF9Lj5pjDkDhF68LyJrgRe0QKjiIOP8WRJ+mMO5bd/i5h9McPvB7Pt62mUXkskvj6CyBLXsR+AdPXmhykHGjx9PkyZN8G/QiaBWj2Dz8i2CI1Dq+jisSBhj0kVkMLCKrCGwc4wxO0RkHBBtjFniqCxKFUTK31s5uWw6GecSCGjyIIHNemHz9LmuApGduHvw+qFqePedif+6eZzdsozkvZsIve8FvCteMfBPKUs4dJ6EMWYFsCLHY6Pz2LaNIzIpBbnPUcjMzGTcuHEcWzAO9+DylH34dbzK3lLo723z8iW47SB8a7Yifvl0ji14idKtB1Cq0f0OuTylUlejy3IolYvk5GT69+/PZ599hl+duwi+5ylsnj5F+p7eFWtRbsCbnFzxJgk/zCb16F5C730Ocddhs8o6WiSUyuH48ePce++9/Pbbb0ybNo0Zx2912C96m5cfZe5/mcQNizn908ccO3eKsAdGYvP2d8j7K5WTrt2kVDbpZ0/SunVrdu7cyZIlSxg2bJjDu3xEhMA7HiKk8zBS43Zx9JMRZJxLcGgGpS7SIqGUXfqZYxybH0VcXByrVq2ic+fOlubxv+1OwnqMJf3MUY4tGKmFQllCi4RSQPrZeI4ueJnM82fxf+BVHl6eSHjUcsvnJ/iERxDWfQzpice0UChLaJFQJd6pU6c4vmg0mecTCes5Hq/yt1od6TLeletdKhTHF48lMzXZ6kiqBNEioUq0c+fO0blzZy4kHKLMAyPxKlfD6ki58q5cjzJdXyLt+H5OfDkJk3HB6kiqhNAioUqszMxM+vbty8aNGylz34v4hEdYHemqfG6OJKTjM6QciCH+mxkYk2l1JFUC6BBYVWKVbt6LxA1fU/rugfje2szqOPniX7ctGWfjOb1uHu6lywP3WR1JuThtSagS6ZNPPiFxw2L863egVEPn+qINuOMh/G67kzM/f8KSJbqajSpa2pJQJc6mTZt47LHH8KpUh+B7Bl33PAirlvMWEYLbD+ZCfBz9+vVj48aN1KpVyyHvrUoebUmoEuXUqVP06NGDsmXLUub+lxA351zywubhRZkHXsbHx4euXbty5swZqyMpF6VFQpUYxhgGDBjAkSNHWLx4MW6+gVZHuiHuAWX4/PPP2bdvH4MGDcIYvf6WKnza3aRKjOnTp7N06VLeeustGjVqBJ87ZqJcUU7I67fsDKWa92Xhwo/59nQZ4le+XWTvpUombUmoEmH9+vVERUXx4IMPMmTIEKvjFKqApt3xDm9Awnez2L59u9VxlIvRIqFc3tmzZ+nbty+VKlVi9uzZLneNBhEboZ2HIl6+PPTQQyQn64xsVXi0SCiXc3HNpYvdPEOHDuXAgQPMmzePwEDnPg+RFze/0oR2foFdu3bx0ksvWR1HuRAtEsqlhXUfzYcffoh/4240b97c6jhFyic8glIN72PGjBnc1HuS1XGUi9AioVxWRvIZ4r95G4+wqgS16Gt1HIcIat0f9+AKxK94k8TERKvjKBegRUK5JGMMp1bNJDM1idDOw5x2PkRB2Ty8Cb33eTLOxvP8889bHUe5AC0SyiUl7/6Z5D9+JahlPzzLhFsdx6G8KtQkoMmDzJkzh2XLllkdRzk5cfYJOJGRkSY6OtrqGKoYqfTspxz+8CncA8pQ9uFpiM3N6kgOZ9IvELDqFU6fPs3OnTspVaqU1ZFUMSMiW4wxkdfaTlsSyuUk/DCHzPOJhHQcUiILBIC4e/DBBx9w6NAhRo4caXUc5cS0SCiX8t1333Fu27cENOmGZ1g1q+NYqmnTpvzzn//knXfeYePGjVbHUU5Ku5uUyzh//jx169blwKlkyj36DjYPL6sjWSp2cicSExOpXbs2wcHBbNmyBQ+Pq5/At2plW+V42t2kSpxXX32Vv/76i5D2Q0p8gbgoICCAd999l23btvH6669bHUc5IS0SyiXs2LGDadOm8dhjj+FdpZ7VcYqVLl268OCDD14qokoVhBYJ5ZSyL71hjGHw4MEEBgYyefJkq6MVSzNmzMDd3Z3nnnvO6ijKyWiRUE5v4cKFrF27lokTJxIaGmp1nGKpfPnyjB07lmXLluncCVUgWiSUU8tMTWbYsGHcfvvtPPHEE1bHKdaeeeYZatWqxbPPPktKSorVcZST0CKhnNqZXz/l8OHDzJw5Eze3kjknIr88PDyYMWMG+/bt47XXXrM6jnISWiSU07pw8iCJ0V/zj3/8g6ZNm1odxym0bduW7t27M2nSJGJjY62Oo5yAFgnllIwxnFrzPjYPbz1ZXUCvv/46NpuNoUOHWh1FOQEtEsopnf9jPSkH/ktgy36UKVPG6jjFUs6LL11UuXJlRo4cyZdffsnq1astSqechUOLhIh0EJE9IrJXRKJyef5JEdkmIjEi8rOI1HZkPuUcUlNTSVg7B4/QypRqcK/VcZzSsGHDqFatGkOHDiU9Pd3qOKoYc3fUG4mIGzATuAeIAzaLyBJjzM5sm803xrxv374LMB3o4KiMyjm88847pJ8+SthD4xCb2xW/lNWVcvuMztXryb6v/sXs2bMZNGiQBamUM3BkS6IxsNcYs88YkwZ8CnTNvoExJvultPwA515YShW6EydOMH78eLyrNcSn6u1Wx3FqvjWa4VXxNl555RW9ip3KkyOLRAXgYLb7cfbHLiMi/xSRv4CpwDO57UhEBopItIhEnzhxokjCquLp1VdfJSkpidJ3PmZ1FKcnIpS+63FOnDjBpEm5XxM7r/MaquRwZJGQXB67oqVgjJlpjLkZGAGMym1HxphZxphIY0yknrQsOXbu3Mn777/PoEGD8AytbHUcl+BVrjqPPPIIb7zxBvv377c6jiqGHFkk4oBK2e5XBA5fZftPgfuLNJFyKsOHD8ff35+xY8daHcWlTJw4ETc3N1566SWro6hiyJFFYjNQXUSqiogn0AtYkn0DEame7W4n4E8H5lPF2OrVq1mxYgWjRo3SIa+FrGLFirz44ossXLiQlLhdVsdRxYzDioQxJh0YDKwCdgGLjDE7RGScfSQTwGAR2SEiMcBQoL+j8qniKz09nWHDhlG1alWGDBlidRyXNHz4cMqXL0/C9x9gTKbVcVQx4rAhsADGmBXAihyPjc52+1lH5lHFX3jUcs7GrOTU9u0sXrwYLy+9mFBR8PPzY9KkSQwYMIDkXevwq93a6kiqmNAZ16pYy0xL4czPn+BVoTYPPvig1XFc2sMPP4xHWFVOr5uHybhgdRxVTGiRUMXa2eivyTiXQFCbRxHJbYCcKiw2m43SrQeQfvooZ2NWWh1HFRNaJFSxdfLkSc5s/Byf6k3xrljL6jglgnfV2/GuUo8zvywgMzXZ6jiqGHDoOQmlCmLSpEmYCykEtXoEyH1pCVW4RISg1o9y9OPnSdz0BUEt+1kdSVlMWxKqWDpw4AAzZ87Er85dOnHOwbzKVce3ZksSN39JetIpq+Moi2mRUMXSmDFjsn7VtuhrdZQSKajVw5iMdM78Mt/qKMpiWiRUsbNt2zY+/vhjhgwZgnuATpyzgkfp8pSK6EjSf1dzIT7O6jjKQlokVLHz8ssvExAQoMtEWCywWS/Ew4vTP31sdRRloQIXCRHxs18bQqlCt27dOpYtW0ZUVBTBwcFWxynR3PyCCGjcjeQ/fiX10G6r4yiLXLNIiIhNRPqIyHIROQ7sBo7Yl894Lcd6S0pdN2MMI0aMoHz58jzzTK6rxCsHC2h0Pza/IBLWfoQxenmXkig/LYkfgJuBl4CyxphKxpgwoCWwAZgsIjpOTt2wJUuWsH79esaMGYOvr6/VcRRg8/QhqHkfUuN2sGzZMqvjKAvItX4diIiHMeaqc/Tzs01RiYyMNNHR0Va8tboB2ec8xE7uRHp6OvXr1yc9PZ0dO3bg7u5+xXaqaMRO7nTpdm6ft8lI5/Dsp6lRLoitW7fi5qa9za5ARLYYYyKvtd01WxIXv/xF5E3JY10EqwqEch0ff/wxO3fuZNKkSZcKhCoexM2doFaPsHPnTubNm2d1HOVgBTlxnQQsERE/ABFpJyK/FE0sVZKcP3+eMWPG0LhxY7p162Z1HJUL31ub06hRI0aPHk1KSorVcZQD5fsnmzFmlIj0AdaKSCpwDogqsmSqxCh/72BOx8VxoeXTuoifBfLTpSciTJkyhbvuuouZM2cybNgwByRTxUG+WxIicjfwBFnFoQzwjDFmXVEFUyVDZkoSiRsW4V21Id6V61kdR13FnXfeSfv27Zk4cSKnT5/Oc7vwqOWX/pTzK0jn70jgFWPMzyJSF1goIkONMd8XUTZVApzZ+BmZKUmUbp11EUL9YineJk+eTIMGDZg6dSqTJk2yOo5ygHy3JIwxdxljfrbf3gZ0BCYUVTDl+tLPnuRs9FL8arfB86ZqVsdR+RAREUGfPn148803OXTokNVxlAPkZzJdXiOajgB3X20bpa7mzC8LMJkZBOpy1E5l/PjxpKen8+qrr1odRTlAvibTicgQEblsvWYR8QTuEJG5QP8iSadc1oX4OJK2fkupBh3xCCprdRxVANWqVePJJ59kzpw57N6ty3W4uvwUiQ5ABrBARI6IyE4R2Q/8CfQG3jDG/KcIMyoXdPqnjxEPLwLv6Gl1FHUdRo0ahY+PDyNHjrQ6iipi1zxxbYxJAd4F3hURDyAUOG+MyXt4g1JXsWnTJpL/+JXA5r1x8wuyOo7Kp5yz5IcPH86YMWPYsGEDTZs2tTCZKkoFGQLbEVgHrAVmiYj+q1AFZowhKioKm28gAY0esDqOugFDhw4lLCyMqKgoXfzPhRVkxvW7wDCgKTALmCYivYsklXJZq1ev5ocffiDwjp7YvHQRP2fm7+/P6NGj+fHHH1m5cqXVcVQRKcg8iWPGmIvLcKwRkfXARmBB4cdSrigzM5OoqCjCw8PJjOhodRx1Ay52PZmMCrgHlWXEiBG0a9dOF/9zQQVpScSKyAT7qCaAC8DZIsikXNTChQuJiYlh/PjxiLuH1XFUIRA3D4JaPsy2bduYP1+vh+2KClIkDNANOCgiPwN7yVrHSS86pK4pLS2NUaNGUa9ePfr06WN1HFWIfGu1pEGDBrzyyiukpqZaHUcVsoIs8NcbQES8gTpAffvfhyJSzRhTqWgiKlfw4Ycfsm/fPpYvX47NppdWdyUiNqZMmUK7du147733AP3d6EoKvHC/fUhstP1PqWtKSkpi3LhxtGrVio4d9VyEK7rnnnto27YtEyZMwPfhd7F5+VkdSRUS/Umnitybb77JsWPHmDx5si4F7sImT55MfHw8iRu/sDqKKkRaJFSROnnyJFOnTqVr167ccccdVsdRRahhw4b07NmTxOivyEhKsDqOKiR6nUhVqHLOyp00aRLnzp1jU2h7XQa8BJgwYQILF3/G6V8XENLu6Suez/nvQxV/2pJQRebAgQPMnDmT/v374xla+dovUE7vlltuwb9+B5L+u4oLp3QpcVegRUIVmds6P0ZahmG1dyuroygHCmrWC3Fz5/S6/7M6iioEWiRUkUg7Ecu57d8TcHtn3APKWB1HOZCbf2kCGj1A8u51REfrIEhn59AiISIdRGSPiOwVkahcnh9qX4p8q4h8JyJVHJlPFZ7TP32MePoQcEcPq6MoCwQ07obNJ4ARI0bo4n9OzmFFQkTcgJlkXfa0NtBbRGrn2Ox3INIYUw/4DJjqqHyq8KT8vY3zezcR2LQ7bj4BVsdRFrB5+RLYrCfff/893377rdVx1A1wZEuiMbDXGLPPGJMGfAp0zb6BMeYHY0yy/e4GoKID86lCYEwmCT/Mwa1UKKUiu177BcpllYq4l/DwcKKiosjMzLQ6jrpOjiwSFYCD2e7H2R/Ly2PAN7k9ISIDRSRaRKJPnDhRiBHVjUre9RNpR/8kqNXD2Dy8rI6jLCTuHowfP57ff/+dhQsXWh1HXSdHFoncptrm2lkpIv2ASOC13J43xswyxkQaYyLLlNGTosVFSkoKCT9+jEdYNfxuu9PqOKoY6NOnD/Xr12fUqFGkpaVZHUddB0cWiTgg+yKAFYHDOTcSkbbASKCLMUaXlHQib7/9NhmJxyl95z8Q0YFzCmw2G//617/Yt28fs2bNsjqOug6OnHG9GaguIlWBQ0Av4LI1o0WkAfBvoIMx5rgDs6kbFB8fz8SJE/GpFolPeITVcZSDXW02fYcOHWjTpg3jx4/Hq1nJXGwAABb0SURBVM87ekVCJ+Own3vGmHRgMLAK2AUsMsbsEJFxItLFvtlrgD+wWERiRGSJo/KpGzNhwgTOnj1LUJtHrY6iihkRYcqUKRw/fpzEzV9ZHUcVkEPXbjLGrABW5HhsdLbbbR2ZRxWOv/76i5kzZ/LYY4+xOlintqgrNW7cmAcffJAvln5JqYiOuPmXtjqSyidd4E9dcrXF1/J6LjxqOSe+mkw6bqzwbKn/oFSe/vWvf/H5l19zet08Qjo+Y3UclU96dlHdkNRDu0je8zMBjbvh7h9sdRxVjFWvXp2AhveRtPVb0o79ZXUclU9aJNR1y8zM5NT3H+LmV5qAxg9YHUc5gcBmPbH5lOLUdx/och1OQouEum7z588n7fAeglr3x+bpY3Uc5QRs3v4EtexH6sHtfPnll1bHUfmgRUJdl6SkJEaMGIFn2er41bnL6jjKifjXb49HaBVeeOEFUlJSrI6jrkHPM6rrMnnyZA4fPkzZfq/pxDmVp9zmT4jNjdJ3P8H+haN46623GDFihAXJVH7p/92qwPbv38+0adPo168fXhVqWR1HOSGf8Ai6dOnCxIkTOXr0qNVx1FVokVAFNnz4cNzc3Jg8ebLVUZQTmzZtGikpKbzyyitWR1FXoUVCFUjKga18/vnnvPzyy1SocLVFfJW6uurVqzNkyBBmz57N77//bnUclQctEirfTGYGp76bRXh4OEOHDrU6jnIBr7zyCiEhIQwZMkSHxBZTWiRUviX9dxUXTsQybdo0fHx0yKu6cUFBQUyZMoVffvmFefPmWR1H5UKLhMqXjOQznP7pY7wq16Vbt25Wx1EuZMCAATRt2pThw4dz+vRpq+OoHHQIrMpVzqGLCWs/IjPtPMH3PIVIbtePUqpgsv8b+/Ldd4mMjGT06NHMmDHDwlQqJ21JqGtKidvJuW1rCGj0AJ6hla2Oo1xQgwYNePrpp5k5cyYxMTFWx1HZaJFQV2UyMzi1+l3cSpUhsFkvq+MoFzZ+/HhCQkJ4+umnyczMtDqOstMioa7q7JalXDgRS3Dbgdg8va2Oo1xYUFAQr732GuvXr2fu3LlWx1F2WiRUntLPnuT0z5/gc3MjfKo3tTqOKgEefvhhmjdvzosvvsipU6esjqPQIqGuIuH72ZCZQem2g/RktXIIm83GzJkzSUhIICoqyuo4Ch3dVOLldQH78/u2kLx7HYEt+uIRVNbBqVRJVr9+fYYNG8bUqVPp27cvrVu3tjpSiaYtCXWFzNRk4le9g0dIJQKbdLc6jiqBxowZQ7Vq1XjiiSd0OXGLaZFQVzi9bh4ZiScJ7vAM4u5hdRxVAvn6+jJr1iz+/PNPxo8fb3WcEk2LhLpM6qFdnN2yjFK3d8K7oi4Drqxz9913M2DAAKZOncrWrVutjlNiaZFQl5j0C8R/MwO3UqEEtXrE6jhKMW3aNIKDg3n88cfJyMiwOk6JpEVCXXJm/SIuxB8kpP0/sXn5Wh1HKUJCQnjrrbfYvHkzb7/9ttVxSiQd3aQASDsRy5kNi/Cr3QafmyOvum1eI6KUKgojfvfH5+ZGDB0eRadOnahevbrVkUoUbUm4qPCo5Zf+rsVkXCB++RvYvP0pffcTDkinVP6JCMHtByPuHvTv31+7nRxMi4TizC+fknbsL0La/xM330Cr4yh1BfdSIQTf8yTr16/n9ddftzpOiaLdTSVMzpZF6uE9nNmwGL86d+Nbo5lFqZS6Nt9arfHd8ysjXh7Jm3v8ODz7n1ZHKhG0JVGCZV5I4eTy6bj5hxDcdqDVcZS6qqxup39i8/InfvkbpKWlWR2pRNAiUYKd/nEu6acOEdLpOWxeflbHUeqa3HwDCWn/T9KO/cXEiROtjlMiaHdTCXU+NoazW5ZSquF9+FSpb3UcpfLNt8Yd+NW5i4kTJ9K+fXuaNcvqJs3elRo7uZNV8VyOFokSKOPcaeKXvY57cEWCWve3Oo5SBf6CD277JO7n9tO7d29iYmIoXbp0UcYr0bS7qYQxJpOTy98gIyWJMl1HYPPQCwkp52Pz8uXTTz/l8OHDPPHEExhjrI7ksrRIlDCJm74iZf8Wgu9+As+wqlbHUeq6NW7cmIkTJ/L555/zwQcfWB3HZTm0SIhIBxHZIyJ7ReSKK4qISCsR+U1E0kVE16guZJs2beL0T3PxrdEM/4iOVsdR6oaERy1nZnwtvMMb8OQ/h5B24oDVkVySw4qEiLgBM4GOQG2gt4jUzrHZ38AAYL6jcpUUZ86coVevXrj5BxPc8Rm90pxyCSI2QjsNRTx9OblkCplpeu2JwubIlkRjYK8xZp8xJg34FOiafQNjTKwxZiuQ6cBcLi8zM5P+/fvz999/E3rfi7h5+1sdSalC4+ZfmtDOw7hw8iDxq97W8xOFzJGjmyoAB7PdjwOaXM+ORGQgMBCgcuXKN57MxU2ePJmvv/6aN954gzeP6uJoqni7ngUkfao2IKhlP06vm8fZcrcCnQs/WAnlyCKRW//GdZV8Y8wsYBZAZGSk/my4ivP7tjDqs7H07t2bZ599ljdfWmF1JKWKRMAdPUg9soeEH2bz888P06JFC+DqRUfnU1ybI7ub4oBK2e5XBA478P1LnAunj3Jy6WvUqVOHDz74QM9DKJd28fyEe2AYPXr04MiRI1ZHcgmOLBKbgeoiUlVEPIFewBIHvn+JkpmWwokvJ4ExfPnll/j56bIbyvXZvP0p88BIEhMT6d69O6mpqVZHcnoOKxLGmHRgMLAK2AUsMsbsEJFxItIFQEQaiUgc0AP4t4jscFQ+V5KZmcnJZdO4cHw/Ife9wM0332x1JKUcxrNMOB999BG//vorjz/+uJ7IvkEOXZbDGLMCWJHjsdHZbm8mqxtK3YCXXnqJ839uoPRdT+B7cyOr4yjlcA899BB79uxh9OjRBLaEoGa9rI7ktHTtJhfz4YcfMnXqVPwb3EupyC5Wx1HKEuFRyzGmAX633cmZdf+HR+kK+NVqaXUsp6TLcriQ7777jqeeeop27doR3HaQnqhWJZqIENLhGbwq1CZ+xRukHt5jdSSnpEXCRcTExNCtWzdq1KjBokWLEJub1ZGUspy4e1Cm20jc/IM5/tmrXDh58NovUpfRIuEC/vzzT9q3b09gYCArV64kMFCvU63URW6+gYQ9NA5sNo4teoX0xBNWR3Iqek7CCVxtrf3Dhw/Trl07MjIyWL16NZUqVcr58uuawaqUK/EoXZ6beozj6Pwoji18hbJ9p+Dmqz+m8kNbEk7s1KlTtG/fnpMnT/LNN99Qs2ZNqyMpVWx53lSNsO6jyUg8zvHPxpKZmmx1JKegRcJJJSQkUK52Y7bv2o1vpygaNdKhrkpdi3elOoR2GUHa0b84vngsZ8+etTpSsadFwgklJCTQtm1b0k4eIOyBkfiER1gdSSmn4Vu9CaFdXiT18G46duyoheIatEg4mYsFYvv27VkFQifLKVVgfjVbENrlRTZs2KCF4hq0SDiRjOQzlwrEV199pQVCqRvgV7MFn376KRs2bKB9+/YkJCRYHalY0tFNTiL9zDGOLXyF+JQEvv76azp06AA//m/Uko5gUqrgunfvzqJFi+jduzetWrVi1apVlC9f3upYxYq2JJxA2olYjv7fcDKTz7BmzZqsAqGUKhTdunVjxYoVxMbG0qxZM/744w+rIxUrWiSKuZ9++oljn4wA4Ka+U2jWrJnFiZRyPXfffTdr164lOTmZFi1asHnzZqsjFRtaJIqxOXPm0LZtW2x+pSnbbxqeZcKtjqSUy3pw8VE87p/AqTQbTZq1YOHChVZHKhb0nEQxVOXFJSSs/Yizm7/Cu0oEofdH4ebtD+i5B6WKkkdwBco9Mp0TX06kV69e7Nixg7Fjx2Kzldzf01okipmEhAROfD6e8/uiKdXwPkrf9bgu1qeUA7n5BnJTz4ncfWoJ48ePZ+fOnXz00UeUKlXK6miW0CJRjGzatImHHnqI8wcPEdzuaUo1uNfqSEq5tLxa5uLuwezZs7ntttt48cUX2bZtG4sXL6ZevXoOTmg9cfZL+0VGRpro6GirY+Tpat1DFxfrM8YwY8YMhg8fTvny5Ult9Qxe5W91VESl1FWk/L0N29oZJCQkMGPGDB5//HGXuFaLiGwxxkRea7uS29FWTBw9epT777+f5557jg4dOvDbb79pgVCqGPGuXJeYmBhatmzJwIED6d27N/Hx8VbHchgtEhZatGgRderUYdWqVUyfPp2vv/6a4OBgq2MppXIICwtj5cqVTJw4kc8//5w6deqwdOlSIKu34OKfK9IiYYGMc6c58fUUevbsSdWqVfn99995/vnnXaIJq5SrstlsvPzyy5TpN534DG+6dOmCf922ZKQkWR2tSGmRcCCTmcHZ31dw+INBJP+xnvHjx7N+/Xpq1apldTSlVD553lSNcv3fIPCOnpzb8QOHPxhE0rY1GJNpdbQioSeui9jFJmjqkT84tfo90o7+iXeVegS3fQqP0CuvIqeUch5px/dxavV7pB7ahVeFWmxYNp+ICOdYul9PXBcTFxKOcGLJVI5+PJSMpHhC7xtOWM+JWiCUcgGeYdW4qe8UQu59jgsJh2nYsCGPPvoof//9t9XRCo0WiSJy7NgxhgwZwuEPn+T83o0E3tGT8o+/j1/t1nruQSkXImLDv25byj/xb55//nkWLFhAjRo1GDZsmEuMgtLupkJ24MABpk2bxocffsiFCxfwqduOwOa9cffXUUtKubrYyZ34+++/GTt2LHPnzsXX15dBgwbx6fl6uJcKuWw7q+W3u0lnXBeS7du3M23aND755BPSM8HvtjsJbtodj+AKVkdTSjlQ5cqVmTNnDittjUhcv4jXp78BNhv+t91FQONueIRUzPe+sg+rtaqwaJG4AampqXzxxRe89957rFu3Dl9fXwYPHszi1Pq4B5SxOp5SykKeoZUJve8FAlv2I3HTlyRtXU3S1tV4V6nHZ5Hn6dq1Kx4eHlbHvCbtbiqA8KjlGGO4cGI//UIPMnfuXI4fP061atV46qmnePTRRwkJCXHZSTVKqeuXcS6BpK3fcjbmGzIST1CuXDkeffRR+vbtS+3atXN9TVG2JPLb3aRFIp/279/P7QPGcm7HWi7E/427uzv33nsvTz/9NAO/T0NExwAopa7NZGYws5WN9957j5UrV5KZmUn9+vXp06cPvXr1onLlype21SJRCIqqSGRkZLBx40aWLl3KsmXL2L59OwBeFWvjV7sNvrc2x803sNDfVynl+i5+4R89epTa/cZwbuePpB3ZA4DnTTfjc3MjfG5uhGe56pd+gFpVJPSchJ0xht27d/Pjjz+ydu1avvvuO06ePIm7uzstW7bk9ddfZ/qfQbgH3mR1VKWUiyhbtiwBkV0IiOzChYQjJP/xC+f3bubM+kWc+fVTbH5B+FSJwKtSHf78swa33HKLw4fQO31LwqtcdVOu/5uX7ue32iYmJvL7778THR3Nxo0b+fHHHzl+/DgAFSpU4M4776Rz5860b9+eoKAgQK8Kp5QqHNm/p3L7Xsk4n8j5fVs4/1c0qX9vJeNcAgDlypWjdevWNGnShMjISCIiIvD397+uDCWmu+laRSI9PZ39+/fTPGoeF+IPcuFELJUyj7Jnzx4uHnvlypVp1aoVbdq0oU2bNlSrVg0R0aKglLKcMYb0U4dIObiNlL+3kXpwBxlJWZP0RIRatWrRsGFD6tSpQ+3atalVqxbh4eG4ubld8R2W/fuxRHU3ZaaeI/3McdITj/POO/s5cOAA+/fvZ/fu3fz555+kpaVd2tbNP5gad7WgT58+REZG0rBhQ8LCwixMr5RSeRMRPEIq4hFSkVIRHQFY/+ztbNmyhejoaLZs2cKaNWuYN2/epdd4e3tTo0YNTqQG4B4YhnvgTbgHhrFjRzhVqlQpUOvDoS0JEekAvAW4AR8aYybneN4L+BhoCMQDPY0xsVfdp83NkGP1RS8vL6pUqcLBzCA8Qipl/QVnfcg27ys/nLy6qLQloZRyBrGTO5GQkMDu3bvZuXMnu3btYteuXaxeH0N64gnIuHDZ9iEhIcTHxxevloSIuAEzgXuAOGCziCwxxuzMttljQIIx5hYR6QVMAXpebb82n1IENH4wq1oGlME9MAybbxCpIlxP+0ALg1LK2Vw+VPaxyx43JpOMc6fJOHOMNzpVJDY2lgMHDvD+++/na9+O7G5qDOw1xuwDEJFPga5A9iLRFRhrv/0Z8I6IiLlKc8c9oAyBTboVTWKllHJyIjbc/YNx9w+mV6//9ZoUxyJRATiY7X4c0CSvbYwx6SJyBggBTmbfSEQGAgPtd1MPTOm8/UaCyZQbeXWRCyXH8bsYVz4+Vz420OMrlq72fZbjuVvzsz9HFoncBvfmbCHkZxuMMbOAWQAiEp2ffjVnpcfnvFz52ECPz9mJSL5mITtyLYk4IPuVdioCh/PaRkTcgUDglEPSKaWUuoIji8RmoLqIVBURT6AXsCTHNkuA/vbb3YHvr3Y+QimlVNFyWHeT/RzDYGAVWUNg5xhjdojIOCDaGLMEmA3ME5G9ZLUgeuVj17OKLHTxoMfnvFz52ECPz9nl6/icfsa1UkqpoqPrWyullMqTFgmllFJ5cokiISLjRWSriMSIyGoRKW91psIkIq+JyG77MX4pIkFWZyosItJDRHaISKaIuMxwQxHpICJ7RGSviERZnacwicgcETkuIjc0P6m4EpFKIvKDiOyy/9t81upMhUVEvEVkk4j8135sr17zNa5wTkJEAowxifbbzwC1jTFPWhyr0IhIO7JGeqWLZE2HMcaMsDhWoRCRWkAm8G/gBWOMY65FW4TsS9D8QbYlaIDeOZagcVoi0gpIAj42xtSxOk9hE5FyQDljzG8iUgrYAtzvCv/9JOtiFH7GmCQR8QB+Bp41xmzI6zUu0ZK4WCDs/MhlAp4zM8asNsak2+9uIGuOiUswxuwyxuyxOkchu7QEjTEmDbi4BI1LMMb8hAvPXzLGHDHG/Ga/fRbYRdZqEE7PZEmy3/Ww/131+9IligSAiEwUkYNAX2C01XmK0D+Ab6wOoa4qtyVoXOJLpqQRkXCgAbDR2iSFR0TcRCQGOA58a4y56rE5TZEQkTUisj2Xv64AxpiRxphKwCfAYGvTFty1js++zUggnaxjdBr5OTYXk6/lZVTxJiL+wOfAczl6K5yaMSbDGBNBVo9EYxG5apeh01x0yBjTNp+bzgeWA2OKME6hu9bxiUh/oDNwt7PNQi/AfztXkZ8laFQxZu+v/xz4xBjzhdV5ioIx5rSIrAU6AHkOQnCalsTViEj1bHe7ALutylIU7BdrGgF0McYkW51HXVN+lqBRxZT95O5sYJcxZrrVeQqTiJS5ODpSRHyAtlzj+9JVRjd9Ttayt5nAAeBJY8wha1MVHvsyJV5kXa0PYIOrjN4SkQeAt4EywGkgxhjT3tpUN05E7gXe5H9L0Ey0OFKhEZEFQBuyltI+Bowxxsy2NFQhEpEWwDpgG1nfKQAvG2NWWJeqcIhIPWAuWf8ubcAiY8y4q77GFYqEUkqpouES3U1KKaWKhhYJpZRSedIioZRSKk9aJJRSSuVJi4RSSqk8aZFQSimVJy0SSiml8qRFQqlCZr8WwT322xNEZIbVmZS6Xk6zdpNSTmQMME5EwshaQbSLxXmUum4641qpIiAiPwL+QBv7NQmUckra3aRUIRORukA5IFULhHJ2WiSUKkT2S19+QtaV6M6JiNMvVqhKNi0SShUSEfEFvgCGGWN2AeOBsZaGUuoG6TkJpZRSedKWhFJKqTxpkVBKKZUnLRJKKaXypEVCKaVUnrRIKKWUypMWCaWUUnnSIqGUUipP/w8ffiHmvzGNjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plotting the histogram with the data.\n",
    "plt.plot(x, y, color = 'black', linestyle = 'solid', label = 'PDF') # original PDF\n",
    "# plt.plot(x, y1, color = 'lightgrey', linestyle = 'dashed', label = 'CDF') \n",
    "plt.hist(xs, bins = 100, density = True, label = 'sampled')\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('$p(x)$')\n",
    "plt.title('inversion sampling')\n",
    "plt.axis([-3, 3, 0, .5])\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}