phys512/prob_distributions_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d46a72f",
   "metadata": {},
   "source": [
    "# Probability distributions Exercise 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b53886bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "seed = 5885\n",
    "rng = np.random.default_rng(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "142589c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll need this a few times, so write a function to take the x values and \n",
    "# plot a histogram comparing to the analytic function func\n",
    "def plot_distribution(x, func):\n",
    "    plt.clf()\n",
    "    plt.hist(x, density=True, bins=100, histtype = 'step')\n",
    "    xx = np.linspace(0.0,max(x),100)\n",
    "    plt.plot(xx, func(xx),':')\n",
    "    plt.yscale('log')\n",
    "    plt.xlabel('x')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84da2c1a",
   "metadata": {},
   "source": [
    "## Rejection method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b2778dd0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9dXhe15XovSXbYouZ2TIKHFtkJ2ZIYnbYkKRNA820oabtFKZJ06ZpqJ0JNk3j2OGYAwZZRsmSTJJBFjMzWGxL+v74ZXU5nfZOO7dzfe+X9zyPHknve84+ey/eay+wGh0dHTWWy3JZLstluSyX5bJc1+iyvtYTsFyWy3JZLstluSzXN/uyGCOWy3JZLstluSyX5bqml8UYsVyWy3JZLstluSzXNb0sxojlslyWy3JZLstlua7pZTFGLJflslyWy3JZLst1TS+LMWK5LJflslyWy3JZrmt6WYwRy2W5LJflslyWy3Jd02vstZ7A33ONjIyY+vp6M378eGNlZXWtp2O5LJflslyWy3JZrr/jGh0dNZcuXTL+/v7G2vpv+z/+nzBG6uvrTVBQ0LWehuWyXJbLclkuy2W5/htXTU2NCQwM/Jvf/z9hjIwfP94Yw2KcnZ2v8Wwsl+WyXJbLclkuy/X3XN3d3SYoKOjPevxvXf9PGCNyNOPs7GwxRiyX5bJclstyWa7/x67/KsTCEsBquSyX5bJclstyWa5relmMEctluSyX5bJclstyXdPLYoxYLstluSyX5bJcluuaXhZjxHJZLstluSyX5bJc1/SyGCOWy3JZLstluSyX5bqml8UYsVyWy3JZLstluSzXNb0sxojlslyWy3JZLstlua7pZTFGLJflslyWy3JZLst1TS+LMWK5LJflslyWy3JZrmt6WYwRy2W5LJflslyWy3Jd0+sfNkaOHj1qli1bZvz9/Y2VlZXZuXPnf/nM4cOHTUJCgrG1tTWRkZFm06ZN/42pWi7LZbksl+WyXJbr/4/XP2yM9Pb2mtjYWPPqq6/+XfdXVFSYm266ycydO9fk5eWZRx55xHz72982+/bt+4cna7ksl+WyXJbLclmu//9d/3CjvKVLl5qlS5f+3fe/8cYbJiwszLz44ovGGGMmTpxoMjIyzMsvv2wWL178j77+n3oNDAyY8vJyExISYhwdHa/pXCyX5bJclstyWa5v6vU/HjOSlZVlFixY8LXPFi9ebLKysv7mM4ODg6a7u/trP/8TV1xcnJk8ebLZti3HjI7q56OjxpSU8Ht01JjiYn5GRvTzkRFj9u0zZu9eYwoL//P3/9Ul7xgeZpx9+/i7pOQ/jzMyYkxaGr+vfr6w0Jh33uG5v5z3X3vmL99fXGxMURG/r57H1c8MDxuzaRO/R0aM2b+fZwQ2f229f+vz/wq+f8/9fw9M/1fzkb+vhvHVuPhr65N7/xp+r75H4Ck/hYWMJ/D8a+PJ2oeHwdeVK4q3v1zP31rf1bj+e2H1l8/Luq+e01/+/Zfjynuvptu/xOVfg7d8XlT0n+FzNU1efW9x8dfh+dee/1sw+XuuvwXrv4XzvwWXv0Zr/yuY/a1x/x78/a2x/nJ+fwsnQnN/iT+h3/9qHn85378Gw+JiYwoKkFNXrvxtfroan1ev5y/x+Lfe8b/CxV/j9b/22V9b6/8KJ/+IvP1fwfDvXc9fu/+fMY//jsz4n7isRkf/+1OwsrIyO3bsMCtXrvyb90RHR5t77rnH/PjHP/7zZ19++aW56aabTF9fn7G3t/9Pz/ziF78wTz311H/6vKuryzg7O/93p/ufrptuusl8+eWX5ic/edPMnv0dExrK51VVxiQmGrN9uzH+/sZ4eRlTXW1MRoYxP/uZMU1N/F1Swv3BwcYkJ/P5wIAxtrbGhIcbExFhTHo631tZ8RMVBaE8/7wxCxYYc+CAMSEhxjQ28uzDDxuTnW1MUpIxWVl8l5nJc6Wlxtx9N+OUlEB03d3GeHsbc++9EHBfnzH29jy7cKExW7YY8+ST+owxjFVcbMzBg8b4+vKODz4wxs/PGGtrYyZPNmbRIu595x1jAgIwury8jElNNSY/35i5c3nHtGnGNDcz/w0bjBkzhvf4+RnT0MC75BodZQ2jo3yflQUsKiqA9/jxer8wiFz+/sbU1XFvbS1zCgtjXXJ/VpYxq1cbk5NjTGioMZGRxpSV6fuys1mrvz/3JiXpZ76+xjz9tDEJCcZ0dgLDpCTgZGvL94Ije3twd/CgMUFBfNbcbIynpzGffWZMbCzzysszZt481pyayjz6+7k3OdmYrVuNcXeHPs6eNWbJEgy/1FRjjh4FFw8+aMzJk8bMnw8tJSYy5+BgaHLePGNeeAHY5+ezlt5eY1pawL+VlcI9IgKYVlUxf+uvtiJpacY4OvJMe7sxrq7MJyjImPh4Y2pqjGlrMyYuDhoXXB88aMzUqdDG6tXGvPaaMRMmAHsnJ3BZUgLstm4FLi0tfF9VBa16eMBLq1czl8JCaMLa2piJExnjwAFgm50NvfX2sq6+Pj7v62M9paXg3MqKNU2dasz588DhapoaGeH9ISG8JzKS+wcGjImJYQ3FxbwnN/fr616/nnX394ODpiZocu5cY6Kjjdm8mbU0NBhz7Bi0esstxtTXG3P99cbcdZcxL70Evz30EGsKDVX+7u1lrnZ2fLd+PTz1l0ZOVRX8sGQJ7/zud5mfwLynBz4ICzPGx4dnWlqYc2YmcG9vN+bGG43Zs4c15uXxu6aG+2V9DQ3Mo7z867CrrkYmWFlBP6OjzL+5mfeOjBhz6JAxFy8as2yZMadPg9OCAta+Y4cxa9YwflUVzw4MANtbbuHz8nJjOjqMcXMzZvFilXEODsBb1lpdbYyNDePMmfN1XOTkwMuZmeBz0SJjWluBi4MD8y4ogM/HjIHXRG5UVX19jZGRykvp6Yzb2GhMZaUxU6awmduwgeeFFq+Wf4JDKyu+lzHq61lDcDDrHx1VfDg4KB9frZtSUhQGfX3M/+23jfn+9415/31jVqxgzNBQ7hseRu+sXAmdvPeeMXfeCT0uXAiPi8z4Z17d3d3GxcXlv9Tf/1caI4ODg2ZwcPDP/3d3d5ugoKB/ujHy/e99z/z7f/yHWbXyCbNh4/OmsVGV8969/O7sBMlnzkDsVlbGfOc7MFl1tTHjxnFfSoox27bx2bRpxsyaBTFFRxtz7hzKMDgYRk1Lg+D37IFxMjIg7uRkCGL+fO4pL+c5Y1AOsbHGTJrEmKOjEF5ZGUyydi2E6u2N0ElKQrk+9pgxFy4wb09P7mlrg4FCQxGkHR0Ii61bEcb29ihEa2t2M7/5DYxWVgZj3HADgiw6GmYaGQFGVlbG/OhHzPdqxSBMWFnJGuvqjNm5E6Hf0GDM4KAxLi7G3HOPMkJJiQpme3vgOjKCwr1wgXX29AAPPz/wsWqVMa++ivF1+jR4EcHZ0IBA6+hgbVZWytjbtiF0k5KM+f3vWeukSYyxciX4WbIEQebgwHxk/YWFwNXGBnynpHC/lRUK8uOPed+YMcBocBBBsHs3vxsbUTzJySgfX19jTpxgviEhwPGhh6CHRYugpbg4Y44fR3BfuADcMjIwpEZHwamPD4J71SrwKAZYXx/Czt0dwVNUhJD29+c+b28MUD8/5hoUZMzQEPMShZySAg5TUsBpUJAxf/qTMTfdhIJZsoT1hoQonD08jPniC3ijsdGYGTOAZ2YmcHJyAq5eXmoAtLZCo2LMiLEcEsIaGhvVMCorU0Pc2prNgGwErjY4ysqgKXd3cBgQwDN9fcDTzY13Zmbyd30992dnY0gcOAA+6+vh2dxcY6ZPRz4kJGAA7dsHPA8dgi5qatgs/PKXxvz0p8Y884wxn3yCQnF1ZQ1z5wKrykp9f2Ag4/7wh/BDYSH4bWiAX44fBzaLF7O2VavgKWMwamXTkZwMv0VHG3P4MDTs68u9Z88C64EBnisowAANDWXN/f3A7swZ8Jqdzbuzs4F7Wxv3RkQY8+67jOPvz2e5udB2ays0MHcudOrpCU0tWYKsEgN6yxbeNW+eMWPH8vf+/fCsqyvvGBlhTvb2wOL4cejC19eYL7+
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGwCAYAAABhDIVPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9/ElEQVR4nO3deXxU9b3/8dfMJJOEBJJATCAhEQRFWUwUCCIqoqlI1RZX7HVBbLm/9g5WDO0tdNF6q1Jr5aYXp1L31rpEa91QKRKVRbGEaFCMIFHUGCBhzZAhZJKZ+f0RmTAkkIVJzsyZ9/PxmIecb86c8xmC5M35bha/3+9HREREJEJYjS5AREREpCsUXkRERCSiKLyIiIhIRFF4ERERkYii8CIiIiIRReFFREREIorCi4iIiESUGKMLCDWfz8e2bdvo27cvFovF6HJERESkE/x+P/v37yczMxOr9djPVkwXXrZt20Z2drbRZYiIiEg3VFVVMXjw4GOeY7rw0rdvX6Dlw/fr18/gakRERKQzXC4X2dnZgZ/jx2K68HKoq6hfv34KLyIiIhGmM0M+NGBXREREIorCi4iIiEQUhRcRERGJKAovIiIiElEUXkRERCSiKLyIiIhIRFF4ERERkYii8CIiIiIRReFFREREIkpYhpelS5cyYsQITj75ZB555BGjyxEREZEwEnbbAzQ3N1NYWMjbb79NcnIyY8eO5fLLL2fAgAFGlyYiIiJhIOyevKxbt45Ro0aRlZVFUlIS06ZNY/ny5UaXJSIiImEi5OFl1apVXHbZZWRmZmKxWHjppZfanON0OhkyZAjx8fFMmDCBdevWBb62bds2srKyAsdZWVlUV1eHukwRERGJUCEPL263m9zcXJxOZ7tfLy4uprCwkDvuuIMPPviA3Nxcpk6dSm1tbbfu19jYiMvlCnr1mIZ9UHw97P685+4hIiIixxTy8DJt2jTuuusuLr/88na/vmjRImbPns2sWbMYOXIkS5YsoU+fPjz22GMAZGZmBj1pqa6uJjMz86j3W7hwIcnJyYFXdnZ2aD/Q4Zb/Gj59lYanrmfjN/vYWF1H9b6GnrufiIiItNGrY148Hg9lZWUUFBS0FmC1UlBQwNq1awHIz89n48aNVFdXU19fzxtvvMHUqVOPes0FCxZQV1cXeFVVVfVY/dvPvI21/tHM2P4DLn3gXS5dvIaC+1cqwIiIiPSiXp1ttGvXLrxeLxkZGUHtGRkZbNq0qaWgmBjuv/9+pkyZgs/n47//+7+POdMoLi6OuLi4Hq37kN3WNH7QuICiGWcwPD2Jytp6nnnuGdzV6ZByZq/UICIiEu3Cbqo0wPe+9z2+973vGV3GUVgYnp7E6KxkYtw7+LO9iJQX74d+r0D2eKOLExERMb1eDS9paWnYbDZqamqC2mtqahg4cOBxXdvpdOJ0OvF6vcd1nS6xWNnsy2Zkoo9qXw7+6joAUhPtZKUk9F4dIiIiUaRXw4vdbmfs2LGUlJQwffp0AHw+HyUlJcyZM+e4ru1wOHA4HLhcLpKTk0NQbcf6njCYK/kN8TX72P3g+kD7ibF1PD3vcgUYERGRHhDy8FJfX09lZWXgeOvWrZSXl9O/f39ycnIoLCxk5syZjBs3jvz8fIqKinC73cyaNSvUpfS4rJQEls+bwl63J9DWuP7vjPzgt3y6cg97x90A6EmMiIhIKIU8vKxfv54pU6YEjgsLCwGYOXMmTzzxBDNmzGDnzp3cfvvt7Nixg7y8PJYtW9ZmEG+kyEpJaA0mfj8NrtUkWDy8ta6cB9YOBSAh1saKeZMVYEREREIg5OHl/PPPx+/3H/OcOXPmHHc30ZEMGfNyJIuFhOueZs+//87Fgy/lYquNytp65hZ/yF63R+FFREQkBMJytlF3GDHmpV1WK/0n3kj/Q8d+H0tii0j99GvI/AlYLMbVJiIiYgJhtzGj2fT7YikX20oZtPZOqPvG6HJEREQinmmevIQr10mXcs+/3uXivJHY3f3ArenUIiIix0PhpYelJsXzpHU6D73vhffXAJBtqeHC2E+YfdvvyErtY3CFIiIikcU04SUsBuy2IyslgRXzJgemU1u8HjJf+D6pdZ9Qs7ovfO+3xhYoIiISYUwz5sXhcFBRUUFpaanRpbSRlZLA6KxkRmclMyo7jYOnXUmtP4V9I2YYXZqIiEjEMU14iRgWC7vH/IjJjYv49EA/NlbXsbG6jpot66GDKeYiIiJiom6jSJKaaIfYROYWlwMwyrKVF+130DDsQhKufQLsGgcjIiJyNAovBjhyHMzBsir8ZbD3QDOf13rA0gRoRpKIiEh7FF4Mcvi2AtWJM/lBqZ/KL/vjeuBdAKz4iI+18ea8KQowIiIihzFNeAnX2UadkZWSwOJ5NwVt8Bi78m6qP/03rt2nkZUy1MDqREREwotpwkvYbA/QTUEbPO6vwffF3xhha+Dr7e/DMIUXERGRQzTbKBz1zeDz773EH5uuxjX0u0ZXIyIiElYUXsJU44DTeMB7OZW19WysruOTr3ZQ/8+fQv1Oo0sTERExlGm6jcwmNdFOQqwtMJ36rphHGRVTgmfHh9h/skq7U4uISNRSeAlTR06nrqlM5NOSz9g9/DZStrkC52k6tYiIRBuFlzB2+CDe1MQJXLTiD7jf8sNbLRs8jrJsZX/MAJ6ZN10BRkREooZpwkskT5XujKyUBJbPmxJ4EmM7uIchz99GfUMDdV+fCClnG1yhiIhI7zBNeIn0qdKdETSdes8eDvZJxXXAyqeNaRysrgPUjSQiIuZnmvASdfoPZfd/vMH/+9NSKl/YAmwBIC22kZfnTVOAERER01J4iWBZaf3567yrA11J7vIXGbbuN9RvssJZ040tTkREpIcovES4QFeS30/9a8UkWVz4t70HTDe6NBERkR6h8GIWFgtfTX2Clx76LaNP/C+GfTsGBjQORkREzEXhxURSkvvxd+v3aXj+k29b/NwT8yirLeP49bxCBRgRETEFhRcTOXJhu36fv0LOW29xlX8lW2uugJTTDK5QRETk+JkmvJh9nZfOCppOnTGDXbUf8uCHB7k8KdPYwkRERELENBszOhwOKioqKC0tNbqU8BFjZ8fEO3jUe0lgg8dNn21m9/p/Gl2ZiIhIt5nmyYu07/ANHm14ecp+D6daP6XO9RXJF9xmdHkiIiJdpvBickHjYHxNxK2cyP4tX1KTcT7mXIdYRETMTuElChw+DmbjlN8yZeN4fu1Jx/PtdOrY+m0kZQzRbCQREYkICi9RJjXRjju2P3OLywEYZqnmFfuvedV/Huf+9BGy0lIMrU9ERKQjCi9R5sjp1P0rnqTPux6yfdso/drF3kYLoIXtREQkfCm8RKGg6dRZc9iVPpz5/9hH1XMff3uGn4TYGFbMm6wAIyIiYUfhRUjLvZhnT2wIPI3h3T+x8aMy9tXlKbyIiEjYUXgR4LCnMa7t+DYtZnSMh9KPXmFjzIzAOepKEhGRcKDwIsH6DWL35U/zxvOPcvu7mfDumsCXEmJt6koSERHDmSa8aHuA0DlhzHe4MPsczvy2G8ni9WB9r4gry89kr9uj8CIiIoYyTXhxOBw4HA5cLhfJyVp+7XgFDepd9kvY5ORv9lPAf6GxhYmISNQzzd5G0oNOvYSmxIEsab4MLBajqxERkSin8CIdGzKJz65ZyQrf2Na23Z9D437jahIRkahlmm4j6Vn+mJYupMraeqye/Qx78QpsVgsx1z0HJ5xicHUiIhJNFF6kUw7fnXq45Rv+at+PBQt4k8g0ujgREYkqCi/SKUduK1D+zTn830uruN+b2BpevE1gizWsRhERiQ4KL9JpQTOQgM3+HCpr6wFIrF5N5ppfUXfJEtJHTDSqRBERiQIKL9Ith3cjgZ9/2u9gqPVLSp66n8zrH2RAoj1wntaFERGRUFJ4kW45shvJ2vgiVaV/4t6Nk9n32LrAeVqVV0REQk3hRbotuBspGU66n9emtm7wGLvyHuZ/PIi97okKLyIiEjIKLxJSgUDzyUvw2RKetcey1T0N0KrHIiISGlqkTnrGsAuoO+lSHvReRnPiQKOrERERE1F
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return np.exp(-x)\n",
    "\n",
    "xmax = 10\n",
    "x = xmax * rng.uniform(size = 10**5)    # Generate x values up to xmax\n",
    "y = rng.uniform(size = 10**5)\n",
    "\n",
    "# Reject\n",
    "ind = np.where(y <= f(x))\n",
    "ind2 = np.where(y > f(x))\n",
    "\n",
    "# Plot the points to show which ones are accepted and which rejected\n",
    "plt.plot(x[ind2],y[ind2], 'bo', ms=0.1)\n",
    "plt.plot(x[ind],y[ind], 'ro', ms=0.1)\n",
    "# Plot the analytic solution for the boundary (see below)\n",
    "xx = np.linspace(0.0,xmax,100)\n",
    "plt.plot(xx, f(xx), 'k')\n",
    "plt.xlabel('x')\n",
    "plt.show()\n",
    "\n",
    "plot_distribution(x[ind], f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cb90035",
   "metadata": {},
   "source": [
    "## Transformation method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4ddeb108",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGwCAYAAABhDIVPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEU0lEQVR4nO3df1yV9f3/8cc5/PYHKCIgKGlZlj+CEjXLpi6Wo7JZn5rb+kG2+dn6HirDtXRbtrbStS3H1s4nP225PmvVqDatdLmSzB9liRCVkSaFZRigKRw5Ij/Oub5/kAcQUNAD1zkXz/vtxm1c17m4rhdn7vDc+6fNMAwDERERkSBhN7sAERERke5QeBEREZGgovAiIiIiQUXhRURERIKKwouIiIgEFYUXERERCSoKLyIiIhJUQs0uwN+8Xi/79u1j4MCB2Gw2s8sRERGRLjAMg8OHD5OUlITdfuK2FcuFl3379jFixAizyxAREZFTsHfvXoYPH37CaywXXgYOHAg0//LR0dEmVyMiIiJd4XK5GDFihO/v+IlYLrwc6yqKjo5WeBEREQkyXRnyoQG7IiIiElQUXkRERCSoKLyIiIhIUFF4ERERkaCi8CIiIiJBReFFREREgorCi4iIiAQVhRcREREJKgovIiIiElQCMrysWbOGMWPGcPbZZ/OXv/zF7HJEREQkgATc9gBNTU3k5OSwYcMGYmJimDhxItdccw1DhgwxuzQREREJAAHX8rJt2zbGjRtHcnIyAwYMIDMzk1deecXsskRERCRA+D28bNq0idmzZ5OUlITNZmP16tXtrnE6nYwcOZLIyEimTJnCtm3bfK/t27eP5ORk33FycjLl5eX+LlNERESClN/Di9vtJjU1FafT2eHreXl55OTkcN9991FUVERqaiqzZs2iqqrqlJ5XX1+Py+Vq89Vj6qoh70b48uOee4aIiIickN/DS2ZmJg888ADXXHNNh68vX76c+fPnM2/ePMaOHcuKFSvo168fK1euBCApKalNS0t5eTlJSUmdPm/ZsmXExMT4vkaMGOHfX6i1V34OH75E3VM3suPzanaU11BeXddzzxMREZF2enXMS0NDA4WFhWRkZLQUYLeTkZHB1q1bAZg8eTI7duygvLyc2tpaXn75ZWbNmtXpPRcvXkxNTY3va+/evT1W/xcX3sVWYzxzv/guV/3pDa56ZAsZD29UgBEREelFvTrb6MCBA3g8HhISEtqcT0hIYOfOnc0FhYby8MMPM3PmTLxeLz/5yU9OONMoIiKCiIiIHq37mC/tcXy3fjG5cy9gdPwASqtqeebZZ3CXx8OgC3ulBhERkb4u4KZKA1x99dVcffXV3foZp9OJ0+nE4/H0UFXH2BgdP4DxyTGEuiv4n/BcBq16GKJfhBGTevjZIiIi0qvdRnFxcYSEhFBZWdnmfGVlJYmJiad1b4fDQUlJCQUFBad1n26x2dnlHcHh/il84E1hR3mNxsGIiIj0sF5teQkPD2fixInk5+czZ84cALxeL/n5+WRnZ5/WvXuv5aXFwKHD+S/uJbKymi8f3e47f0ZYDU8vvIbkQVG9VouIiEhf4ffwUltbS2lpqe+4rKyM4uJiYmNjSUlJIScnh6ysLNLT05k8eTK5ubm43W7mzZt3Ws91OBw4HA5cLhcxMTGn+2t0SfKgKF5ZOJND7gbfufrtf2ds0S/48p0jMPMHvVKHiIhIX+L38LJ9+3ZmzpzpO87JyQEgKyuLJ554grlz57J//36WLFlCRUUFaWlprFu3rt0g3mCRPCiqTQtLzaubiLI14Kr8lB3lNQAM7h+uVhgRERE/sRmGYZhdhD8da3mpqakhOjrar/feUV7DVY9sYc3t0xif3HHrTvkhN3/4/VKeb7gI71dDiqLC7KxfOEMBRkREpBPd+fsdkLONToUZY146kjy4P3fm/Jybv+pKKq10EblqHmwvg8tuA5vN1PpERESCnWXCixljXjrTuisp+uMXSQkpwLv1fUifDYN6cAVgERGRPsAy4SVQuc68iqX/eYNvpo0l3B0Nbo2DEREROR0KLz1s8IBInrTP4bG3PPDWFgBG2Cq5LOwD5t/1K5IH9zO5QhERkeBimfASKGNejpc8KIr1C6f7plPbPA0k/fNbDK75gMrNA+HqX5hboIiISJCxTHgJpDEvx2szndow+OK8/6JqaznFcVeRrOnUIiIi3WKZ8BI0bDa8F/0/rnzjLA6+WAk0b5WQFrYXZ87N6kYSERE5CYUXEyQPiuKlhbN8XUn7P9rGxa//jPrV6+GGv0G4AoyIiEhnLBNeAnXMS2dadyXtLS0HoLbey2dVDWBrBNSVJCIi0hHLhJdAHvNyMvbUuXx3fR2le2Jx/emN5nN4iQwL4dWFMxVgREREWrFMeAlmyYOieGThLW02eAzb+CDlH76N68vzSB40ysTqREREAovCS4BoMyPpcCXeT/7GmJA6PvviLThL4UVEROQYhZdANDCBj69ezQt5f2F0/6/h0nRqERERH8uEl2AbsHsy/VJSedx+HXV5xQBEUs994U8z40fLGZaUYm5xIiIiJrIZhmGYXYQ/dWdL7e7aUV7DVY9sYc3t0xif3PODgsur63zjYPq/ejej9uRRFzeeKMcW7U4tIiKW0p2/35ZpebGi1uNgdk/8AR9+8iaRU37OKAUXERHpwxRegkR97BiubljG8tDzcX81BibywPv0jxvOsOSR5hYnIiLSixRegsTg/uFEhIWx4KsxMINx8e+InxKGh6ob/kn8OZPMLVBERKSXKLwEieN3pw5zfUq/l+M4UH2YkpoYRmlGkoiI9BEKL0GkzVowyedTPvQ//PAPayj9525gNwBxYfW8sDBTAUZERCxL4SWIJcfF8n8Lr/e1xriLV3HWtnup3WmHi+aYW5yIiEgPsZtdgL84nU7Gjh3LpEl9a+xH8qAoxifHMD4pmnHlecTZXPTf96bZZYmIiPQYy4QXh8NBSUkJBQUFZpdiDpuNT2c9wYON36Ny0k/MrkZERKTHqNvIQozQSP7suYpxB+rBXgOGQdKWxRjnfJMhF37L7PJERET8QuHFQgb3DycqLMQ3nfoq+1b+FP40DR/mUTF0G4kjRptboIiIiB8ovFjI8dOpbZ7JlL5ewzO7vFxjH0qiyfWJiIj4g8KLxbSZTg3suOwBHi/ZwoSqWgBC3RUMPVzCkInXmFWiiIjIaVF4sbjWXUkheHgqfCnn2j+kpmYPMV+/y+zyREREuk3hxeLadCV5G4nYOJXDu/ewLWwyw75alRe0Mq+IiAQPhZc+oHVXUvlVS8l8eAqfr60GtgAwjC+pDotn/cLpCjAiIhLwFF76mORBUeQt/JZvUG94dSln/uv75NVfTLXrQoUXEREJeJYJL06nE6fTicfjMbuUgNdmUO++dzA8Rxlpq8CwWeafg4iIWJhl/lo5HA4cDgcul4uYmBizywkek37AHm8CC1a5+NmBOt/idoMHRKgVRkREApJlwoucuvAxGdSGbfQtbvfDkJc4O6SCi29/nKShQ8wtTkRE5DgKL9JmRlKou4Kz//E8Id5G9pashek3m12eiIhIGwovArQeBxND2Tf/xoYX/4/JZ13NCLMLExEROY7Ci7TjTr6EXzYZ5O53g82GzdNA3Lv/g2eKg6T4OLPLExGRPk7hRdo5foPHn4c+yQ9CX6Zo+xqMBRtJHtzP3AJFRKRPU3iRdo7f4LHfF2HUrS/if1yzuWrPIQ4daQS0Kq+IiJhD4UU61GYtmOTL2XfGVt74YwHrv2qNGWn7gsOhQ3hx4TcVYEREpFcpvEiXJA0d4muNsTccJuX5Rew/XI+7PAUGTTS7PBER6UMUXqTLfK0xVV/QQBNhNg8ltZE0aINHERHpRQEZXq655hpef/11LrvsMp5//nmzy5HjxZ/LgZvy+X/ONby3ag+wB4BQmggLi9AGjyIi0qPsZhfQkTvvvJO//e1vZpchJ5CUmMSjC29ize3TWHP7NDZca1A85OeMbtrtG+grIiLSEwIyvMyYMYOBAweaXYacRPKgKMYnxzA+KZpR7+UywP0Z3w553eyyRETE4rodXjZt2sT
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+N0lEQVR4nO3deXhU5cH+8XtmkpnsAyRkg0BAEJBVQSLg8rOmUuVFqS8VlwqCS6VYwbxaFgW0KuBG0UKhUlGrVVFftxpEMYDKC4oSUGgRRMCwJSRAMiF7Zs7vj5TRVHaSeZKZ7+e6cumcnJm5c4wzd57nOXNslmVZAgAACBJ20wEAAAAaEuUGAAAEFcoNAAAIKpQbAAAQVCg3AAAgqFBuAABAUKHcAACAoBJmOkCg+Xw+7d27V7GxsbLZbKbjAACAk2BZlkpLS5Wamiq7/fhjMyFXbvbu3au0tDTTMQAAwGnYtWuX2rZte9x9Qq7cxMbGSqo7OHFxcYbTAACAk+HxeJSWluZ/Hz+ekCs3R6ai4uLiKDcAADQzJ7OkhAXFAAAgqFBuAABAUKHcAACAoEK5AQAAQYVyAwAAggrlBgAABBXKDQAACCqUGwAAEFQoNwAAIKhQbgAAQFAxWm4++eQTDR06VKmpqbLZbHr77bdPeJ+VK1fqvPPOk8vlUqdOnfT88883ek4AANB8GC03ZWVl6t27t+bNm3dS++/YsUNDhgzRpZdeqg0bNmjChAm69dZb9cEHHzRyUgAA0FwYvXDmFVdcoSuuuOKk91+wYIE6dOigJ598UpLUrVs3rVq1Sn/84x81ePDgxooJAA3G5/PJ6/XKZrMpLOyHl+Dq6mp5vV45nU45HA5JktfrVUVFhex2u6Kiovz7lpWVqba2VpGRkXI6nZKk2tpalZaWym63y+12+/f1eDyqqalRdHS0IiIi/PseOnRIdrtd8fHx/n1LSkpUWVmp2NhY//PV1taqqKhINptNSUlJ/n2Li4tVUVGh2NhYxcTE+PMWFBRIklJTU+vtW1ZWptjYWC5YjIBoVmtu1qxZo8zMzHrbBg8erDVr1hzzPlVVVfJ4PPW+ADQMn8+nsrIyVVZW+rfV1NRo/fr1+vzzz+vtu379er3yyivasGGDf1t1dbUef/xxzZgxQ7W1tZIky7L0yutvaeStv9HEx+brpVXf6rOt+7T223wNveZXuuKqX+ql5Rs0f9k/9ecPNmrc9Nm64KJL9cQTT2hbgUdfbCvQis37NGDQhera81y9v3qjNu8+oPdy8zRrzp/V6eyuuueee2RZlg6UVmpbvkfn9OilpNQ2+vrrr5W9YY827z6gP8yer/iEBA25ZoS27D2kL7bt1/J/7lXHzl0UFR2ttWvXSqp7Q3/k6WcU7nTqsp9fLsuy/I/duUs32Ww2vfTme/o2v0Tb8kv098VvyOl0KmPghdpfXKbcnQdUVFqhQYMGKSoqSq+/9Q8dOFypg4er9NFHHyk2NlYZFwxQUWmFvt1XrHU7D+jyX1yhFi1a6KVX35DP56vbd+UqtWrVSv369dOhsmr5fD4dKqvWddddp4SEBL3yyiv+475p0yYlJiaqe48e2pbv0XcFHvl8Pt12221KTk7Wc8895/9vseGfW5SSkqKzzz5bh8qqZVmWJOnuu+9Wamqq5s2bJ8uydKisWvn5+WrTpo3at2/v30+S7rvvPrVt29b/h2ljOZLjx88dDM/V3DSFY2N05OZU5efn1/vLQZKSkpLk8XhUUVGhyMjIn9xn5syZevDBBwMVEWiyampqZFmW/y99j8ej5cuXq7a2VsOHD5dUVzbunP6kVq9apfDOF6jVWeepU4JLe/L3K/e5aaqosXT1/Qt1efdk7fNU69U5D+jLnHd03W+yNPi6W1XgqdC33+/Ws3ddLUl6YeVmFZdXa13eQW16e75yl76mK2+4VVeNyVLPNnEqOVyp3//+95Kk5P5D1L19ktyRTr2WnaO3X3hGFw37tZJ6XaI9noPqmhyn9956Q5J07rXjdVhRqvb5tGXbdn2+aqXcrVN06fAabSusUGu3U7m5uaqurNCWfQe1rdypLimx2vDdHn337RZt3dFLxeU1Kq6oUUlljfbs3SvPoQN6YfVOjbkqTR/9a7+qqip18MABHS71aP2uEvksSz7LkudwmSrKy7U136P+kjbu8cjr9am2pkaHK6pUXF4jSSquqJH33y/u+0urtLu4QnEup8qqvJIku0364vtidU6KVd7BCnl9dftakorLa9QiyqkCT5UkyeuzlHewQp7yGkVFOFRV45MkRTkd+v5AueIiw3XkbcRnWXJHhuv7A2VqHx+tipq65yuv9v7kd8KypOLKGrkjnfr+QPmPttc9WnF5jWIjwmWz2WSz2+WODFdxeY1aRjtV7bVkt9tVWetTcXmN3JHh2revRg6HQw6Hw7+fJDkcDoWFhclub9y/p4/k+PFzB8NzNTdN4djYrCZSO202m9566y0NGzbsmPucffbZGj16tCZPnuzftmTJEg0ZMkTl5eVHLTdVVVWqqqry3/Z4PEpLS1NJSQnDo2jWCgsLtXv3brVp00aJiYmS6v4iv/r23+tApV0tLh8rh6RwSfvfmaVD36xS28HjlH7BFaqslqoO7tLGv4xVeGSM+v/+VSVFSxsPSBUfztbu9cuVdOnNir9guGp9UnpkmT58YIRsdrsue3SpfD7prMRYLX/uUX33yds6d+gYZQy/XfnF5fJ4PPrsyVtlc4Tr9wuztWb7AbkcDm1c/r86vGW12p2fqauv/bXKqr3q286tx+/Pks3u0HXjpyomOkY2Sb49Xyl7aY669+mrbudfrE6JMQpzOPTm35+V3WbTORcPUUm1XZbPp5J9O1VdlKfLLuit5M69dKi0QodrfNr91SrtOliuCy+8WK1bxem7/WWKri3R3t156tMlXeecc44OHq7SobJqbfx6vQpLKzX80vO1Jq9CHRMitK/wkByVh3TY51SnDu3kKa9RaVWtPIV7ZZOlzPPPUXR0tLxerz7bulfFBw+qU2ornZ3eVpJ08HCVtu/Ol7e2RrHuFnI6nbJJinbaFO6r1p6SKqXEt9Du4kq1i49ShM0rn8+nSp9DjjCHbLIp1mXX/uIyySa5XC4dOlwlT5VXKTFhahntVFmN1CrGpeLyGvl8Prkj6/5ePVxtyR0ZppKKWsW67CqpqFXLaKe/XBx52T/y89tsUofWMT8pH5Zl/btohUuS/99tNtsxv/ef+wXSjzM19nMH8rmam8Y6Nh6PR263+6Tev5tVubn44ot13nnnac6cOf5tzz33nCZMmKCSkpKTep5TOTiACRUVFZLkL+vr16/Xf902SZ5au+J/cad/v/1vPKiK775Q68F3qm2fX+igpJqiXdr77FjZXdFKm7DYv29R9h9VtilHLf7fzXJn1I3SOMpLtPd/H5IjIkZnXf+ASmql9Ghp24bP5a7MU1S7nopr21WdElw6UFajlMod+vaQpaQOXTW4R4r2eapVXlYmh8OhjkluVXqlAk+F9hwqU3iYUzdltNXGvYd1qKxK6/IOKjHWpZhwh9ISolXjs6tnmzhVeW1Kdru052CF9h+u1FmtY9Q+IVpb9pdr0Fnx9dakAAhtp/L+3axeOQYMGKAlS5bU27Zs2TINGDDAUCLg9FRWVmr9+vUqKipSr169dMn8TZKkwnceVfk3nyr+irsU0+tySVLNgd3au+5D2ZxR9cqNIyZe9ugW8ko6eGRbbLzcF94oR1TdgtIjIzedMm9ReOatqnZGKS5CqqyW7HFuXThljmq8Ut7BWnVyS995pEd+e41uuLi7wsPDz/jn7NPp5Pbr1rb+7eSWsWf83ABCl9Fyc/jwYW3bts1/e8eOHdqwYYNatWqldu3aafLkydqzZ4/+9re/SZLuuOMOzZ07V7///e81ZswYLV++XK+99pqys7NN/QjACeXm5mrlypVa8r1d2yI7S5KqC7/XvkXjZHNGKW3CYv/QrT0iWpJU6yny398R17qusMTEy7Is/76tBo9T/C/uVLikWNUVnDBXlEaNvkUzr8uod3YNAIQSo+Xmyy+/1KWXXuq/nZWVJUk
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = rng.uniform(size = 10**5)\n",
    "\n",
    "x = - np.log(y)\n",
    "\n",
    "# Plot the distribution of x and compare with the exponential distribution\n",
    "plot_distribution(x, f)\n",
    "\n",
    "# Plot 1-y(x) to illustrate that y is related to the cumulative distribution function\n",
    "plt.clf()\n",
    "plt.plot(x, 1-y, '.', ms=0.1)\n",
    "xx = np.linspace(0.0,xmax,100)\n",
    "plt.plot(xx, 1-f(xx), 'k:')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('1-y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9be78685",
   "metadata": {},
   "source": [
    "## Ratio of uniforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c3a3e122",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9dXiWV7bwAd9JIApxdyUEidGSBCtECLSleNupQaftVJjznulUTsel4+20c85Up0aptzgViOAxJILF3T3EE5I83x+/b72LzhyZeT/Oy3f6Pvu6ciV5nvvesvayvfYSC5PJZDLMzdzMzdzMzdzMzdy+Ic3yek/A3MzN3MzN3MzN3MztWjazcmNu5mZu5mZu5mZu36hmVm7MzdzMzdzMzdzM7RvVzMqNuZmbuZmbuZmbuX2jmlm5MTdzMzdzMzdzM7dvVDMrN+ZmbuZmbuZmbub2jWpm5cbczM3czM3czM3cvlFt2vWewP/tNjU1ZbS0tBgzZ840LCwsrvd0zM3czM3czM3czO3vaCaTyRgYGDB8fX0NS8v/3Dbz/5xy09LSYgQEBFzvaZibuZmbuZmbuZnb/0FrbGw0/P39/9Nn/p9TbmbOnGkYBsBxdHS8zrMxN3MzN3MzN3Mzt7+n9ff3GwEBAf9bjv9n7f855UauohwdHc3KjbmZm7mZm7mZ2/+w9ve4lJgdis3N3MzN3MzN3MztG9XMyo25mZu5mZu5mZu5faOaWbkxN3MzN3MzN3Mzt29UMys35mZu5mZu5mZu5vaNamblxtzMzdzMzdzMzdy+Uc2s3JibuZmbuZmbuZnbN6qZlRtzMzdzMzdzMzdz+0Y1s3JjbuZmbuZmbuZmbt+oZlZuzM3czM3czM3czO0b1czKjbmZm7mZm7mZm7l9o9p1VW6OHz9urFmzxvD19TUsLCyMvXv3/pfvHD161IiPjzdsbGyM8PBwY/v27f/t8zQ3czM3czM3czO3/zntuio3Q0NDRkxMjPHyyy//Xc/X1tYat9xyi7FixQqjuLjY+N73vmc8+OCDxqFDh/6bZ2pu5mZu5mZu5mZu/1PadS2cuXr1amP16tV/9/OvvfaaERISYvzxj380DMMwoqKijJMnTxovvviikZ6e/t81TXMzN3MzN3MzN3P7H9T+R/nc5OXlGampqV/7LD093cjLy/sP3xkbGzP6+/u/9vPf2Uwmw6io4Mdk4qeykt/y/dX/G4ZhTE4axvbt/P73npmaMozMTH7/e+//Z3P5j541mQyjvNwwMjIY96+fk3enpnQ9U1P6ztTUf73+qSnt47+CwdWf/TUM/6Nn/rrPv2du8s7UFM+Wlf3tOP/Z8+XlXx//71nbfzb/yUn2dmLib/sXGMgzsv//Hmz+oyZwOXSItZaVKYz+ozX8PbD+r9Z4dZ9//Zn8npz8+pr/o/f/Mxj+Na7J51fj9tV4IbQ0Ocl75eX6/l/D9b+itX8ULleP91/hh7SpKcM4eNAw3n7bMK5c0blfjZd/jRtXj3E13v41Lf8jeHT1Gv6at/1HNHD1fv97+/SfjfOf4Y6s8d+D5X+Gt//VWFfD8u/Zm7/n+38EZv/I3vy1zPjr8crKlOb/I5z79+jnr3nNfyYX/h758ffi1n9HszCZrufw2iwsLIw9e/YY69at+w+fmTVrlnH//fcbP/jBD/73Z19++aVxyy23GMPDw4adnd3fvPPzn//c+MUvfvE3n1++fNlwdHS8JnO/ulVUGEZpqWH4+xvG+fOGkZhoGKOjhmFnZxiWliDFyIhh2NsbRkQECLB7t2Fs3GgYOTmGsXUrfQwP6zs1NYYxYwbI5ONjGHV1hrFihWHMnv11ZKuvN4zAQMNoaDCM4GDDsLDg+bw8/Tw52TCysw2judkwXF0Nw8PDME6cMIzHHjOMggLDSE3lvcpKfdfdnXdtbQ2jttYw5s83jMFB1js1xfNBQYbR2MjvkRG+6+w0jKQk+khM5HdwMHMeHub71FTDqK5WuNjZsb6aGsNISDCM7u6vz8nbm34Mgz4LCgwjJcUwsrIMo6rKMKKj6Sc1lf9DQw3j8GHWffiwYdx4o2G89pphxMQYRn8//fv5GUZUlO6HYRhGeDhEPjoKnOvrDaOnxzB6exnPwoK57NxpGLGxhtHezl4FBdFfayv9STOZmOONNxrGq6/yjre3YRQXG0Z6umH8+MeGkZYGXsyZw/Olpfzu7TWMVatgVosWAfsDB4DPypXALyyM9ZpMzE3GzsxkTWNjrNcwWI+LC895eLAPixaBEw0NrK+iAnhNn24YXV3gR3s770xNGcbQkGGUlBjGli30mZWlOJaYCE6vX89vX1/6z8szDE9Pw+joAC9eftkwliwxjAsXwOer6aOjg3daWhSG4eHgbmIi8K2vBzcNg7ls2AA+BAYaxpEjhjFvnmGcPMl6R0YULyMiwHlnZ+Dl728YNjbg3cKF0JqFBd/5+v7tXgre79hhGOvWGcaePczVML6O11VVOl5lJYLGZAJHHBz+ts+KCvb2arj+4Q/swfg4/f3qVyg7sbGMYzLR14kTwLu8nH0UnK6vh24/+QR6njuX8Ts7mauFBc/NmqW8ROZsYaG4W1XFmktL2ZOQEMNoawPmu3bBS2xtoTeTCVi3txuGlxd76enJ50VFhuHmxp5ERv4tTLOzDSMgQHHAMMDRkhLD2LTJMPLz2fP8fPrz8TGMvj7gZWn5dd4puNjSoryvulppJCTEMH7wA8O4807DmDnTMHJzWffQEN/Jc4J3Cxey1/feC38KCTGM99/X/318GEveEbpwd4c3zp7NZ0NDwCEmBvr394cvjoygKHz0kWEsXsxnM2bQV1UVNC60/rvf6d5s3swzwkcvXjSMd94xjOXLgYOFBXg8e7bitWGAGx4efNbRwTw/+ACcrqoCp4eH+S4kBNhUVcG3du2Ch7i7gzObNimtmkz0nZQEjvw1nv//0vr7+w0nJ6e/S35/45WbsbExY2xs7H//39/fbwQEBPy3KjdDQ4bx1VcgRG0tgkmIMzcXgRkRYRj79/O7pcUwCgsN46mnQMDMTBC9u9sw4uNhElVVMO6gIBC5tRXGIgi6bx/MoqKCPmxted9kgpm0tMCAamtBeC8vELe9HQGZmWkY27bxf3j414V8VhbvfvQRjCA7G+Z6+TLrMgzDsLZm3OFhCFMIvLoahvfeeygeFhYQtKsrTOzcOWWQMTH05+zMPJuaDONnP1OBVV/PWJ6e/C4ups+DBxHeorBs2aLCpKuLeR0/jvLw298axkMPoSiEh8MYV62CgIuLleF+/rlh3HILp+XhYZhDdbXCJCXFMJ5/nufb2vg8OhrhJMqeMANpPj6G8fOfs1/Dw+xtd7cK8jfegCF9+9u8J8qNjQ1zS01FYW5tpb8rV3gvMdEw/vxnPk9OBt729sDLzQ2cy8szDCcnYNHTg2I0NsY7K1YYhpUVMI+OZozGRsO4dIk9SE8HL9avZ95FReDn2Bh75edHn7m5KOdFRbw3fTrr6ewEzgkJ4KeNDUy4oMAwpk1DWba1BSeLigwjLk6F09QU4wQGghvr1yNgkpLYk9xchHlMDLgkNDYyAvweewzcGh3l+akpaOWGG1AYioroy9aWvaysBH99fVXZE6E4NcX/qanQwNy5X1dKbWwQujEx/D86qoccWbu1NbC+WmEXDlxTA9/YtIk1m0zA4c032euEBIRwUpJhHD3KHGXM1FT6Hx1lr+LiVCk/dAg8KCnhs6Ag3q2pYQ4pKco7SkuB9dVCVZS8nBzGmjHDMAYGGPPECcNwdGS/bWwMY+1a9tvTUxXQZcvAp+ZmBOLYGL9DQxlDBO477wC3oSGEZnMzcOnoAKZ9fYZx333AfngYntbezprc3OBjwn+Li+Ere/ZAdzExwNbVlfcSExHkS5eigN91F3RZUMAeidIsimBCgmH88peGcdtt4Ne2bSjnmzYBg8WLWUNdnSpjixYx9ksvcWCxswO2n3yi+Do6yvNz56LM/OUvPNvQwLN+fsAnKckwPvsM2NnYsP4jR8Db228HbomJzOnCBejy0iXW5e5OH7298O2CAuZlGPQlSqeTE+8dPEh
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGwCAYAAABhDIVPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBRUlEQVR4nO3deXhU9f328ffMZCUhGyEJCSAobmxJZYkgFpAojYoiPy32cYmgtNrBikFbaCvWtkptK01rp6Zat1qtEa2IoBRJEURRNoMgggajxmASEJKQIWSZmeePwCQhLEmY5Myc3K/rmsucM2fmfDKlzM13tXg8Hg8iIiIiAcJqdAEiIiIi7aHwIiIiIgFF4UVEREQCisKLiIiIBBSFFxEREQkoCi8iIiISUBReREREJKAEGV2Ar7ndbvbs2UPPnj2xWCxGlyMiIiJt4PF4OHjwIMnJyVitJ29bMV142bNnD/369TO6DBEREemA4uJi+vbte9JrTBdeevbsCTT+8lFRUQZXIyIiIm1RVVVFv379vN/jJ2O68HK0qygqKkrhRUREJMC0ZciHBuyKiIhIQFF4ERERkYCi8CIiIiIBReFFREREAorCi4iIiAQUhRcREREJKAovIiIiElAUXkRERCSgKLyIiIhIQPHL8LJs2TLOPfdczj77bP7xj38YXY6IiIj4Eb/bHqChoYHs7GxWr15NdHQ0I0aM4JprrqFXr15GlyYiIiJ+wO9aXjZs2MCQIUNISUkhMjKSzMxMVq5caXRZIiIi4id8Hl7Wrl3LlClTSE5OxmKxsGTJklbXOBwOBgwYQFhYGOnp6WzYsMH73J49e0hJSfEep6SkUFJS4usyRUREJED5PLw4nU5SU1NxOBzHfT4vL4/s7Gzuv/9+tmzZQmpqKpMnT6a8vLxD96utraWqqqrFo9PUVEDejfDt7s67h4iIiJyUz8NLZmYmv/3tb7nmmmuO+/yiRYuYNWsWM2bMYPDgweTm5tKjRw+eeuopAJKTk1u0tJSUlJCcnHzC+y1cuJDo6Gjvo1+/fr79hZpb+Uv45HV4KQs8ns67j4iIiJxQl455qaurY/PmzWRkZDQVYLWSkZHB+vXrARg9ejTbt2+npKSE6upq3nzzTSZPnnzC95w/fz6VlZXeR3Fxcef9AhPmc7jfxewe8xDb91SxvaSSkoqazrufiIiItNKls4327duHy+UiMTGxxfnExER27tzZWFBQEI888ggTJ07E7Xbz05/+9KQzjUJDQwkNDe3Uuo8q8cSR8YWdms8OAusA+G7wLv4w41ISzxzWJTWIiIh0d343VRrgqquu4qqrrmrXaxwOBw6HA5fL1UlVwQFnHTX1bnKmpzEoIZKvv9zNqP/+iJgX/gBZS6HfqE67t4iIiDTq0m6j+Ph4bDYbZWVlLc6XlZWRlJR0Wu9tt9vZsWMHGzduPK33aYtBCZEMTYlmQO+e7HL3oy56ICSp5UVERKQrdGl4CQkJYcSIEeTn53vPud1u8vPzGTNmTFeW4hMNPRK4sf7nrB6Vy/byWraXVLK9pJJvvv7c6NJERERMy+fdRtXV1RQWFnqPi4qKKCgoIC4ujv79+5OdnU1WVhYjR45k9OjR5OTk4HQ6mTFjhq9L6XSxESGEBgfz4yXFQONA4f+zruW3wU9xYPIfiB17i6H1iYiImJHPw8umTZuYOHGi9zg7OxuArKwsnnnmGaZPn87evXtZsGABpaWlpKWlsWLFilaDeNurK8a8HCslJpxVc8dzwFnnPRe37CnCv6mj6KvdlJxRCTSGnJSY8C6rS0RExMwsHo+5FiypqqoiOjqayspKoqKifPre20squfLRdSy7cxxDU6KPe03JASd//tNDvFx3Ie4jvXLhwVZWzZ2gACMiInIC7fn+9svZRoEsJTaCu7J/yc1HWmMKy6oIe3UGbCqCSXeAxWJwhSIiIoFN4aUTpMSEe1tZonYvpb9tIw3rt7Gr3wTqIxv3bVJXkoiISMeYJrwYMealLWzDpvGHVespr49g8dNFQBEA4cE2Vs0drwAjIiLSTqYJL3a7Hbvd7u0z8xcpsRH8v+xHOOCsI+vIuT1Fn7DmzTwOVF+o8CIiItJOpgkv/qx5NxINdZz16j1cFryNsi3R0PdXhtYmIiISaBReupotmIpB13BwbzEF8VeSUqLp1CIiIu1hmvDir2NeWrFYcF/4Y6549yz2Ly0DGrdKSAsuxpF9MymxPYytT0RExM9pnZd2aMs6L21VUlHjXdxu76cbGPv29dQOmEjUDf+EEAUYERHpXrTOSwBoPg6muLAEgOpaN1+V14GlHlBXkoiIyPEovPgBa+p0frCqhsIv4qj667uN53ATFmzjrbkTFWBERESaUXjxAykx4Tw695YWeyQFr3mQkk8+oOrb80mJGWhgdSIiIv5F4cVPtJhOfbAM9+f/5FxbDV998z6cpfAiIiJylNXoAnzF4XAwePBgRo0aZXQpp69nIruvWsIf669jS8R32V5SyfaSSkoqaoyuTERExHCmaXnx1xV2O6pH/1SetF5LTV4BAGHUcn/IC0y4fRF9kvsbW5yIiIiBTBNezCYlJpxVc8d7x8FEvHUvA794i5r//ADs67Q7tYiIdFsKL36s+TiYz0bcxiefv0dY+i8ZqOAiIiLdmMJLgKiNO5er6hayKGg4ziNbCoTt20ZEfF/6pAwwtjgREZEupPASIGIjQggNDmbOkTEwsVTxRujPCcZF+Q2vkHCOCQYqi4iItIHCS4A4dgxMcNWX9Hgznn0VB9lRGc1AbfAoIiLdhGnCS8BszHgaWqwFkzKckt7/5Ud/XkbhK58BnwEQH1zLa3MzFWBERMS0TBNezDZVui1S4uN4du513tYYZ8GrnLXhPqp3WuHCqcYWJyIi0klME166K29rjMdD9fI8Ii1VePa8B0w1ujQREZFOofBiFhYLX05+hiWP/4qhZ/yYs46MgQGNgxEREXNReDGRmOgo/mW9mprFHx854+GhoCd5xzKSX87NVoARERFTUHgxkWNnJEXtXkr///2Paz1rKCqbBjHnG1yhiIjI6VN4MZkWM5ISp7Ov/EMe+7CWayKTjS1MRETERxRezCwohNIx9/PkpnUMK69uPOUspffBHfQacY3BxYmIiHSMacJLd1jnpSNiI0IID7YxJ68AGy6eD3mI86yfUFn5BdGX3G10eSIiIu1mmvDSHdd5aYsW42Dc9YSsGcvBz75gQ/Bo+mhGkoiIBCDThBc5sebjYEqufJDMR0bz9fIKYB0AffiWiuAEVs0drwAjIiJ+T+Glm0mJCSdv7tXeGUkhFYWc+Z9bebF2LBVVFyi8iIiI31N46YZazEja8yEe12EGWkr5bN9hPDZt8CgiIv5N4aW7G3Ub34b0Zd7LFRS/tO3ISQ/hwUHqRhIREb+k8CLEp36PF8+o8XYl8e6f2f7RZioq0xReRETE7yi8CNCsK6nqG9w7H2VoUB0bP1rK9qDp3mvUlSQiIv5A4UVaiurDt9e8wJuLn2TBu8nw7jrvU+HBNnUliYiI4RRepJXewy5lUr9xXHCkG8niqsP6Xg7/V3ABB5x1Ci8iImIohRc5rhYzklb8HHY6+GfIOeCZZGxhIiLS7VmNLkACwHlXUB+RRG7DFLBYjK5GRES6OYUXObUBF/Hp99ewyj2CwvJqtpdU8umOAvaU7zW6MhER6YZM022kjRk7V0x0tHeDx0gO8XrIL6i3WCjLWkzimcONLk9ERLoR04QXbczYuZpv8Bh64FMSl1vYf6iB/Z5oEo0uTkREuhV1G0mbpcSEMzQlmrOHjuLLa1dya909uMJimi5w1RtWm4iIdB8KL9IhrrAYdnn6e8fAFG1YRu2fR1K+a73RpYmIiMkpvEiHxEaEeMfAXPnoOxxYdj+hVV+Q//wjlFTUGF2eiIiYmGnGvEjXaj4GBsBa+yqfvbeIX20fzytayE5ERDqRwot0WIuF7Ihme+gD1G5fR2F5NQAJm/4I50wm4fxxxhUpIiKmo/AiPtO8KynT+gGPhfyF2i2P8c2sjfTpO9Do8kRExCQUXsRnmnclWetS+fqtXbz8VQ8yLHH
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return np.exp(-x)\n",
    "\n",
    "u = rng.uniform(size = 10**5)\n",
    "v = rng.uniform(size = 10**5)\n",
    "# choose which points to keep -- we renormalize f here to drop the factor of 2\n",
    "ind = np.where(u <= np.sqrt(f(v/u)))\n",
    "ind2 = np.where(u > np.sqrt(f(v/u)))\n",
    "x = v[ind]/u[ind]\n",
    "\n",
    "# Plot the points to show which ones are accepted and which rejected\n",
    "plt.plot(u[ind2],v[ind2], 'bo', ms=0.1)\n",
    "plt.plot(u[ind],v[ind], 'ro', ms=0.1)\n",
    "# Plot the analytic solution for the boundary (see below)\n",
    "uu = np.linspace(0.001,1.0,100)\n",
    "plt.plot(uu, -2*uu*np.log(uu), 'k')\n",
    "plt.xlabel('u')\n",
    "plt.ylabel('v')\n",
    "plt.show()\n",
    "\n",
    "# Plot the distribution f(x)\n",
    "plot_distribution(x, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2d2d0bf",
   "metadata": {},
   "source": [
    "The black curve above separates the accepted and rejected points. To derive this, start with\n",
    "\n",
    "$$u \\leq \\sqrt{\\exp\\left(-{v\\over u}\\right)}$$\n",
    "\n",
    "and invert this to get $v$:\n",
    "\n",
    "$$u^2 \\leq \\exp\\left(-{v\\over u}\\right)$$\n",
    "\n",
    "$$2\\ln u \\leq -{v\\over u}$$\n",
    "\n",
    "$$2 u \\ln u \\leq -v$$\n",
    "\n",
    "$$v \\leq -2 u \\ln u$$\n",
    "\n",
    "which is plotted as the black curve in the figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "922395bd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}