phys512/prob_distributions_solutions_part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d46a72f",
   "metadata": {},
   "source": [
    "# Probability distributions Part 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "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, y_log=True):\n",
    "    plt.clf()\n",
    "    plt.hist(x, density=True, bins=100, histtype = 'step')\n",
    "    xx = np.linspace(min(x),max(x),100)\n",
    "    plt.plot(xx, func(xx),':')\n",
    "    if y_log:\n",
    "        plt.yscale('log')\n",
    "    plt.xlabel('x')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1389685e",
   "metadata": {},
   "source": [
    "### Power law \n",
    "\n",
    "$f(x) = C x^{-\\alpha}$ for $a<x<b$, where $C$ is the normalization factor,\n",
    "\n",
    "$$\\int_a^b C x^{-\\alpha} dx = 1\\Rightarrow C = (1-\\alpha)\\left[b^{1-\\alpha}-a^{1-\\alpha}\\right]^{-1}$$\n",
    "\n",
    "$${dy\\over dx} = C x^{-\\alpha} \\Rightarrow y = {C\\over 1-\\alpha} \\left[x^{1-\\alpha}\\right]^x_a  = {x^{1-\\alpha} - a^{1-\\alpha}\\over b^{1-\\alpha} - a^{1-\\alpha}}$$\n",
    "\n",
    "$$\\Rightarrow x =  \\left[b^{1-\\alpha}y + a^{1-\\alpha} (1-y)\\right]     ^{1/(1-\\alpha)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "64f1b187",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGwCAYAAACOzu5xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7CklEQVR4nO3deXiU5b3G8XtmMkmG7AGyDbvgAihhL4gsEqUoKFgRURGh2vZ0qGCwFrVqPUVxqcqRRtG2ltaKoraKoqCCCIIoCIaKCBhEwEASlpAhIWSZmfNHyISQBJKQ5J3l+7muXGaeefPObxwht89q8ng8HgEAAPgBs9EFAAAA1BfBBQAA+A2CCwAA8BsEFwAA4DcILgAAwG8QXAAAgN8guAAAAL8RYnQBTc3tdmv//v2KioqSyWQyuhwAAFAPHo9Hx44dU0pKiszmuvtVAi647N+/X+3btze6DAAA0Aj79u1Tu3bt6nw+4IJLVFSUpIo3Hh0dbXA1AACgPpxOp9q3b+/9PV6XgAsulcND0dHRBBcAAPzM2aZ5MDkXAAD4DYILAADwGwQXAADgNwguAADAbxBcAACA3yC4AAAAv0FwAQAAfoPgAgAA/AbBBQAA+A2CCwAA8BsEFwAA4DcILgAAwG8QXAAAgN8guDREQbZ04L9GVwEAQNAiuNTX129Kz/SQlt1jdCUAAAQtgkt9dbxUMpkrvspOGF0NAABBKcToAvxGdLI0a7sUmWB0JQAABC16XBqC0AIAgKEILo1RXioV5xtdBQAAQYfg0lCb/yn9qZu0+gmjKwEAIOgQXBog+2ix9pyIkE4cVXHWGm3NLlD20WKjywIAIGgwObeeso8WK+2p1SorC9EA8336/Mfucs9fK5vVohWzhskeazO6RAAAAh7BpZ7yi0pVXObSvIn91DVhuCQpK69QMxdnKr+olOACAEALILg0UNeESPW0x5zS4jGsFgAAgg1zXBpry2s6762rdI35M6MrAQAgaPhkcFm6dKkuuOACdevWTX/961+NLqd2R76X7dBWXWshuAAA0FJ8bqiovLxc6enpWrVqlWJiYtS3b1+NHz9erVu3Nrq06npN0v7yKN29so1eNroWAACChM/1uGzYsEE9evSQ3W5XZGSkRo8erQ8//NDosmqK76wj3W9VvqKNrgQAgKDR5MFlzZo1Gjt2rFJSUmQymfT222/XuCYjI0OdOnVSeHi4Bg4cqA0bNnif279/v+x2u/ex3W5XdnZ2U5cJAAD8UJMHl6KiIvXq1UsZGRm1Pr948WKlp6froYce0ubNm9WrVy+NGjVKeXl5jXq9kpISOZ3Oal8taZD5G7Vf8StpD3NdAABobk0eXEaPHq05c+Zo/PjxtT7/9NNP64477tDUqVPVvXt3LViwQK1atdJLL70kSUpJSanWw5Kdna2UlJQ6X2/u3LmKiYnxfrVv375p39BZjDV/ppjd70ubFrbo6wIAEIxadI5LaWmpNm3apLS0tKoCzGalpaVp/fr1kqQBAwZo69atys7OVmFhoZYtW6ZRo0bVec97771XBQUF3q99+/Y1+/s41SLXSB3ucZs0+M4WfV0AAIJRi64qOnTokFwulxITE6u1JyYmavv27RUFhYToqaee0ogRI+R2u3XPPfeccUVRWFiYwsLCmrXuM9nq6aIDg29V66SYs18MAADOic8th5aka665Rtdcc43RZdRbVl6h9/u4iFC2/wcAoJm0aHBp06aNLBaLcnNzq7Xn5uYqKSmpJUtpEnERobJZLZq5OFMdTTm6wfKJPjJdqoxZUwgvAAA0gxad4xIaGqq+fftq5cqV3ja3262VK1dq0KBBLVlKk7DH2rRi1jAt/c0QLen+iRwh7+hnno+UX1RqdGkAAASkJu9xKSwsVFZWlvfx7t27lZmZqfj4eHXo0EHp6emaMmWK+vXrpwEDBmjevHkqKirS1KlTm7qUFmGPtVX0rlx6u46VOLXq+1T1NrooAAACVJMHly+//FIjRozwPk5PT5ckTZkyRQsXLtTEiRN18OBBPfjgg8rJyVFqaqqWL19eY8JuQ2VkZCgjI0Mul+uc7tNoXYZpT1iqPp6/VunGVAAAQMBr8uAyfPhweTyeM14zffp0TZ8+vUlf1+FwyOFwyOl0KiaGFT4AAAQinzuryN+Fq0Rx21+Vcr8xuhQAAAIOwaWJ/SHkH7J/+jtp/XNGlwIAQMAhuDSx113DVRrVXkq+xOhSAAAIOD65AZ0/2+zppp03rFHP9vFGlwIAQMAhuDQ5k7IOFUvmAm8Lu+kCANA0CC5N6NSddCWP+pt2qERWfRdyvlbMGkZ4AQDgHAVMcDF8HxdV7aSbX1SqNlsWKGnDo8pte6kG7jtP+UWlBBcAAM5RwEzOdTgc2rZtmzZu3GhoHfZYm3raY5Q0aKIUGilr6w4yy21oTQAABIqA6XHxOXGdpLt3av/Bcrkz1xpdDQAAASFgelx8UmiE0RUAABBQCC4tIFmH1erAF0aXAQCA32OoqJlFZK/VmrCZ8nySIvXOlCz8KwcAoLHocWlmxxP7qUARKo1qJx0/ZHQ5AAD4Nf73v5l5QsJ1ZckTenDAUHV12iRnARvSAQDQSASXZhYXEapia9zJTekq2KwWNqQDAKARAia4+MIGdLU5dVM6U/kJFXy7Sjd9EsWGdAAANELABBeHwyGHwyGn06mYmBijy6nGHmuT3VYuzb9MKsxVF9OfjC4JAAC/xOTclhIWJdn7qTQiRSmmw0ZXAwCAXwqYHhe/MPb/tDPfpLUZ7OkCAEBj0OPSkiLbSmar0VUAAOC3CC6G8Chi/3rphNPoQgAA8CsEFwM8GfKCOr83Ufryb0aXAgCAXyG4GOBzd3e5LWFSaZHRpQAA4FeYnGuAJe7BumLkjWrXsYuUXSBJ7KYLAEA9EFxaWFxEqKzWMP3qnf2S9nvb2U0XAICzC5jg4qs7557u1J10JSmk8IDy9n2nW1eI3XQBADiLgJnj4nA4tG3bNm3cuNHoUs7KHmtTT3uMepZ8pQsXD9FPMu9ViMqNLgsAAJ8XMMHFL7UbIIXHqCzKrngdM7oaAAB8HsHFSKGtpP/5TLvHvK48xRldDQAAPo/gYrSoRKMrAADAbxBcfESoyhS743XJxVwXAADqEjCrivyax6N/hz6kdmt+0I8mk46eP4F9XQAAqAXBxQfERYZpkS5VG49Tj3+0W+9+sJZ9XQAAqAXBxQfYY2266c5HdKToIf0yJFwj8wo1c3Em+7oAAHAagouPsLeJk72N0VUAAODbmJzrg1od+EKjzL6/kR4AAC2N4OJrdixTl6UTNMf6kkxlx42uBgAAnxIwwSUjI0Pdu3dX//79jS7l3HRN04nYrvrA1U9mV6nR1QAA4FMCJrj401lFZ2Sxatd1y/X78p/LFR5rdDUAAPiUgAkugcRjCTW6BAAAfBLBxYeFHM+TlqZLR/caXQoAAD6B5dA+LHbFLCl3tfKP5it7xDx20wUABD16XHxQXESobFaL7th7hTa7u+oX3/TUmPlrlfbUamUfLTa6PAAADEOPiw+yx9q0YtYw5RcNkjyT9ZDJpCx20wUAgODiq+yxNgIKAACnYajIH7jKFbd9kf5ufVzyuI2uBgAAwxBc/EGJU0lfPKoRli2K/n6p0dUAAGAYgos/aBWv3P736I9lt8jZebTR1QAAYBiCi5840v1W/c11lWS2Gl0KAACGIbj4I49HKik0ugoAAFocwcXP5Oz6rwpfHK2CV3+urdkF7OsCAAgqLIf2E5Wb0j2+bJuWhX6hcln0y/lv6Yg1SStmDWPpNAAgKARMcMnIyFBGRoZcLpfRpTSLqk3pSpWz06Ki5EH6bXEsm9IBAIJKwAQXh8Mhh8Mhp9OpmJgYo8tpFt5N6ey3S5LKsgsMrggAgJbFHBc/18l0QJYTR4wuAwCAFkFw8WNx376iD0J/p8QNjxldCgAALYLg4sdK4s5XmKlc1uO5kqvM6HIAAGh2BBc/djypv8aUzNGeUQslCxvTAQACX8BMzg1WWz1dlHWwSDKZJFUsm2aFEQAgUBFc/Fjl3i4zF2cqVGX6uWWZ3jMP06uzxhN
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Power law by transformation\n",
    "def f(x, a, C):\n",
    "    return C * x**(-a)\n",
    "\n",
    "alpha = 2.2\n",
    "a = 1\n",
    "b = 10\n",
    "n = 1-alpha\n",
    "C = n/(b**n - a**n)\n",
    "\n",
    "y = rng.uniform(size = 10**5)\n",
    "x = (b**n*y + a**n*(1-y))**(1/n)\n",
    "\n",
    "plt.hist(x, density=True, bins=100, histtype = 'step')\n",
    "\n",
    "# plot the analytic solution, the multiplicative factors account\n",
    "# for our x-axis being in log_10(x)\n",
    "xx = np.linspace(a,b,100)\n",
    "plt.plot(xx, f(xx,alpha,C),':')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$x$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ca403ec",
   "metadata": {},
   "source": [
    "### Lorentzian (Cauchy distribution) by transformation method and ratio of uniforms\n",
    "\n",
    "$$f(x) = {1\\over \\pi}{1\\over 1+x^2}$$\n",
    "\n",
    "for $-\\infty<x<\\infty$\n",
    "\n",
    "$${dy\\over dx} = {1\\over \\pi}{1\\over 1+x^2}\\Rightarrow y = {1\\over \\pi}\\left[\\arctan(x)\\right]^x_{-\\infty} = {1\\over \\pi}\\arctan(x) + {1\\over 2}$$\n",
    "\n",
    "$$\\Rightarrow x = \\tan\\left(\\pi\\left(y-{1\\over 2}\\right)\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "eef56316",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples between +-100 = 99363\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGwCAYAAABhDIVPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABqo0lEQVR4nO3dd3gU1foH8O9sTU8IIQ1Cr6GXEHqRKCL2xlWvInqxBcvF6xUsWH4XwYY1inoF9NoQRUQRBCK9hhJaaIGE9N7b1vP7I2STJT1sMrub7+d58pA9c2bmXSY7+86ZM+dIQggBIiIiIgehkDsAIiIiouZg8kJEREQOhckLERERORQmL0RERORQmLwQERGRQ2HyQkRERA6FyQsRERE5FJXcAdia2WxGWloaPD09IUmS3OEQERFREwghUFxcjODgYCgUDbetOF3ykpaWhpCQELnDICIiohZITk5Gly5dGqzjNMlLVFQUoqKiYDQaAVS+eS8vL5mjIiIioqYoKipCSEgIPD09G60rOdv0AEVFRfD29kZhYSGTFyIiIgfRnO9vdtglIiIih8LkhYiIiBwKkxciIiJyKExeiIiIyKEweSEiIiKHwuSFiIiIHIrTJC9RUVEIDQ1FWFiY3KEQERFRK+I4L0RERCQ7jvNCRERETovJCxERETkUJi9ERETkUJi8EBERkUNh8kJEjsFsAkwGuaMgIjvA5IWI7J8QKPgkAqmf3FyZxBBRu8bkhYjsXnHSMZiyz+NiVhFMZrPc4RCRzJi8EJHdu3T+OMqhRQU0+P1kltzhEJHMVHIHYCtRUVGIioqCycQmZSJnkxNyPW7UdQQgsKCwQu5wiEhmTtPyEhkZibi4OMTExMgdChG1gi5SFj5QR2HaqYVyh0JEMnOalhcicm4SBG5R7oUxRwsIAUiS3CERkUyYvBCR3Qv+81EsUhXiv8YZKPEdhGfMJkDJ0xdRe+U0t42IyHl1zj+Ia5WHcTLodnxZGMbEhaid4xmAiOzeig5Po5OqFKF9+mFXHp82ImrvmLwQkd074DYZPm5qjDKVYoCIB0pHAO4d5Q6LiGTC20ZE5DAizryMb8wLgAt/yR0KEcmIyQsR2bWK4nyUXNiPjvo0lGr8kS585Q6JiGTG20ZEZNcyT+/BOu0iFOb2wZ7rfsMT3x7BKm0YpsgdGBHJhi0vRGTnBFKEHyTvLpg2wB8AkFeqlzkmIpITkxcismulXSZjgu5DJExfJXcoRGQnmLwQkcOQUg7iY/UHGBj3ntyhEJGM2OeFiByGVJaLG5UHkJdTKncoRCQjJi9EZNdKd3yA5eod8Eoqh7l/GBYZZmNKjxG4Ru7AiEg2TnPbKCoqCqGhoQgLC5M7FCKyIXXaIVyvjEE3ZQ6UPiH42jQd7yT2kjssIpKR0yQvkZGRiIuLQ0xMjNyhEJEN7elwK77wmgdFj0lQKRW4dVgwzELIHRYRychpkhcick4X3Idji9uNQEAoAKC7Igt9TBcAIx+XJmqv2OeFiBzKE6fvh0bogOKpQIfucodDRDJg8kJEdq1L+RmYjErAZASUKhRp/AF9GfwMFXKHRkQyYfJCRPbLbMYziY9BAQGUTQI8A/BR6A84mJiPjf795Y6OiGTCPi9EZL8MZShQ+aNMcgNcfSzFxRUG+WIiItkxeSEi+6X1wPXSJ5gd8DOg0gIA3LQqpOSX40J2iczBEZFcmLwQkd3KK9Ujq1iHvgGelrKH3fcgSv0+9MfWyhgZEcmJyQsR2S1xeTyXKf38LWXanFOYqTwIj/w4ucIiIpkxeSEiu6VK3ovl6vfQ/fTnlrLyXjdgkWE2crtcK2NkRCQnPm1ERHZLmReP65UxyM6tvm2k6zIOX5t0uM5vqIyREZGc2PJCRHbL2GUMXjLMQUrPu+UOhYjsCJMXIrJbug598I3pWuQE15hD2qhDVykTroUX5AuMiGTF5IWI7NZ3B5IAAD383C1l2swj2Kn9J/pse0SusIhIZkxeiMhuuZYkYoxXHnr7SJYyrw7+KBVa5Ot5+iJqr9hhl4js1vSEt/GY/iBwWg0M/RsAwKXLYDzeZyPK9CZ8JXN8RCQPXroQkd0yKtQogRvg2kHuUIjIjjB5ISK79VO/d3Gj+3dA3+lyh0JEdoS3jYjI4dyesxzeunQg7yPAt4fc4RBRG7PLlpfbbrsNHTp0wJ133il3KEQko7MZxTCaRa3yIaV7Mbp8F1CUJkNURCQ3u0xenn76aXz99ddyh0FEctIV456LC/GKOQowm6wWxXadjUWG2Ug0+8kUHBHJyS6TlylTpsDT07PxikTkvMpyMV15CFMNuwCF0mqR36S5+No0HXr3YJmCIyI5NTt52blzJ2666SYEBwdDkiSsW7euVp2oqCh0794dLi4uCA8Px8GDB20RKxG1Jy7eeMn4EGL7PS13JERkZ5rdYbe0tBRDhw7FQw89hNtvv73W8tWrV2P+/PlYvnw5wsPD8f7772P69Ok4e/Ys/P0rp7UfNmwYjEZjrXU3b96M4ODmXUnpdDrodDrL66Kioma+IyKyS64d8L35WgzoPhCjrlikMJSiq5QJZWkmALbSErU3zU5eZsyYgRkzZtS7fNmyZZg7dy7mzJkDAFi+fDk2bNiAFStWYMGCBQCA2NjYlkVbhyVLluC1116z2faIyP4FHv8EO7VRyD86B+j5vtzhEFEbs2mfF71ej8OHDyMiIqJ6BwoFIiIisG/fPlvuymLhwoUoLCy0/CQnJ7fKfoiojZVkoTvSoTYU11pk0nihVGhlCIqI7IFNx3nJycmByWRCQECAVXlAQADOnDnT5O1ERETg2LFjKC0tRZcuXbBmzRqMHTu2zrparRZaLU9iRE7n8FeI1vwH8WduByastFqUOegRTNg9GJsnTALH3iVqf+xykLqtW7fKHQIRyUwIM4qFK3Rq79oLpcqJGmOTC9A3gH1eiNobm9428vPzg1KpRGZmplV5ZmYmAgMDbbmrWqKiohAaGoqwsLBW3Q8RtY2fPO7FEN2XqJj8cq1lw0J8AAAxCXltHBUR2QObJi8ajQYjR45EdHS0pcxsNiM6Orre2z62EhkZibi4OMTExLTqfoiobRRVGKFWShjZvWOtZcqSDHzR4Wvckfa2DJERkdyafduopKQE8fHxltcJCQmIjY2Fr68vunbtivnz52P27NkYNWoURo8ejffffx+lpaWWp4+IiJpKrazn+sqkw7Xlm6CvYH83ovao2cnLoUOHMHXqVMvr+fPnAwBmz56NVatWYdasWcjOzsaiRYuQkZGBYcOGYdOmTbU68RIRNWTsmSX4PykLyOsL+Pa0XujeCd+73w8XT1/cJoSlDwwRtQ+SEKL2rGcOrKioCN7e3igsLISXl5fc4RBRC5W/0QOu+jzgsT1A4KBay2//ZA/OZhTjxKvToVAweSFydM35/rbLuY1agh12iZzLp4p78Y7xbsC7c53LR3X3RanehKIKQxtHRkRyc5rkhR12iZzLd8apONf3UcC17pFcwoJU6CplArqSNo6MiORml+O8EBGplRIGBNXfdBy262Hs1B5DSVJHoMMtbRgZEcmNyQsR2R9dCbqaU+FiqH8AOr3GB6VCC8lQ3oaBEZE9cJrbRuzzQuREUmKw2vAk7jr5WL1Vjo6LwkDdShhCa89uT0TOzWmSF/Z5IXIiRh2K4YZylU+9VYSCDcdE7ZXTJC9E5DxMfaZjcMV/8dOgTxqtW6Y3tUFERGRPmLwQkd3ZfCoDABDWw6/eOsNUiXhD9V9cXPtqG0VFRPaC7a5EZHeM5sqxM4d19am3TiDyca/qL6TnZrRRVERkL9jyQkR2p9e5L/Gm6nMokvfXX8l/AL5U/w2Hg/7WdoERkV1wmuSFTxsROQ+35O2YpdoORWFy/ZV8e2CVehZOd5reZnERkX1wmuSFTxsROY/vzNdhmWkW1CEj5A6FiOyQ0yQvROQ8DriMR/bweVD692uwnocohU95EmDUtVFkRGQP2GGXiBzWN+Xz0DG2AAjfBQQNkTscImojbHkhIvti1KOLMQluxsJGq5apvFEiXAA9J2ckak/Y8kJE9iU/AVEFj6Os2BN
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Lorentzian by transformation method\n",
    "\n",
    "import scipy.integrate\n",
    "\n",
    "def f(x):\n",
    "    return 1 / (np.pi * (1+x**2))\n",
    "\n",
    "y = rng.uniform(size = 10**5)\n",
    "x = np.tan(np.pi*(y-0.5))\n",
    "\n",
    "# trim the x values to focus on the center of the distribution\n",
    "xtrim = 100.0\n",
    "x = x[np.where(np.logical_and(x<xtrim, x>-xtrim))]\n",
    "# and also calculate the new normalization between these limits\n",
    "norm, err = scipy.integrate.quad(f, -xtrim, xtrim)\n",
    "print('Number of samples between +-%g = %d' % (xtrim, len(x)))\n",
    "\n",
    "plt.hist(x, density=True, bins=1000, histtype = 'step')\n",
    "xx = np.linspace(min(x),max(x),1000)\n",
    "plt.plot(xx, f(xx)/norm,':')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$x$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0f61e05d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGwCAYAAAC5ACFFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9d3SV15W4fSSBKgj13iVUaJJoamADKmBsmnEv4JLE4CSTxDVxYidOYqfZTjKJMYkruBeqG01gA2pUiabeG+pCXQLpfn88s7NxJvHkN589TrLuWUtL0r3ve8o+u519drGxWCwWY23WZm3WZm3WZm3WZm1/s9l+1ROwNmuzNmuzNmuzNmv7Z25WZcnarM3arM3arM3arO1zmlVZsjZrszZrszZrszZr+5xmVZaszdqszdqszdqszdo+p1mVJWuzNmuzNmuzNmuzts9pVmXJ2qzN2qzN2qzN2qztc5pVWbI2a7M2a7M2a7M2a/ucNu6rnsC/QxsbGzNNTU1m4sSJxsbG5quejrVZm7VZm7VZm7X9A81isZje3l4TEBBgbG3/vv3Iqix9Aa2pqckEBwd/1dOwNmuzNmuzNmuztv9Fq6+vN0FBQX/3e6uy9AW0iRMnGmMAtqur61c8G2uzNmuzNmuzNmv7R1pPT48JDg7+ixz/e82qLH0BTa7eXF1drcqStVmbtVmbtVnbv1j7n1xorA7e1mZt1mZt1mZt1mZtn9OsypK1WZu1WZu1WZu1WdvnNKuyZG3WZm3WZm3WZm3W9jnNqixZm7VZm7VZm7VZm7V9TrMqS9ZmbdZmbdZmbdZmbZ/TrMqStVmbtVmbtVmbtVnb5zSrsmRt1mZt1mZt1mZt1vY5zaosWZu1WZu1WZu1WZu1fU6zKkvWZm3WZm3WZm3WZm2f06zKkrVZm7VZm7VZm7VZ2+e0fyll6eDBg2bZsmUmICDA2NjYmO3bt/+P73zyySdm5syZxsHBwURFRZlXXnnlvz3z7LPPmrCwMOPo6GiSkpLMkSNHvvjJW5u1WZu1WZu1Wdu/ZPuXUpb6+/tNfHy8efbZZ/+h56urq83VV19tFi5caAoLC813v/td87Wvfc3s3r37L8+8/fbb5r777jM//vGPzYkTJ0x8fLxZvHixaW1t/bKWYW3WZm3WZm3WZm3/Qs3GYrFYvupJ/G+ajY2N2bZtm1m5cuXffebhhx82H374oTlz5sxfPrvppptMd3e32bVrlzHGmKSkJDNnzhzzxz/+0RhjzNjYmAkODjbf/va3zfe///1/aC49PT1m0qRJ5sKFC9ZCutZmbdZmbdZmbf8i7R+V3/9SlqX/15aXl2cyMjI+89nixYtNXl6eMcaYkZERc/z48c88Y2trazIyMv7yzN9qw8PDpqen5zM/X3azWIwpL+f3/+aZy7/7e3//9bNjY8aUlfEzOmrM3r3GXLrE77Gxz747Nqaf/yPrGBv77+P+dR9/q0+LRef0t+Ysn/31Gi9/R/4vLdWfy7/7W/D4vDEvh5N8PzZmzJ499H3597Jugefo6OfD/++t6fJnS0sZS+D01+v763n9T/v01+v9vHf+eq6Xr+uv1/w/wVfmKvD667nLsyUlxrz8MmN83pr+en+Ki4156SV+Xz7W36KZy/scHWW8kpJ/bJ8+Dz5/a5zPw5PL9/7y7/4an/8R3vD3cPjvwfqvm8y1pORvz+Pvwe/zxpLn/h5f+bx1/DVeSd8lJZ+lh9FRY155hd9/ixf8PfoTXL58Tv8TTf71c5+3j59H1/9Tu5zuL1367L78rbH+mt/8vfldvn+7d/Nz+R5evq9/C07/05w/j6/8PXly6ZIxv/wlv7/K9m9tWYqOjjZ33nmn+cEPfvCXzz766CNz9dVXm4GBAdPV1WUCAwNNbm6uSUlJ+cszDz30kPn0009NQUHB3+z3Jz/5iXn88cf/2+dfhmVpbMyYfftAoJQUY/LzjQkN1e/y841Zs8YYOzsQzd/fmOZmY6Ki+N8YYyZPNqaigu+amoyprTUmOZnnLv/bxsaYyEjG8/Y2prXVmKEh+ujooJ+XXjImM5PP580zJiDAmMZGY7ZtYx5nz/K9xcKYUVH0e/k6UlOZd0qKzrWiwpjKSmO6u41xdzcmK8uYTZuMWbzYmDNnjMnIYD01NczNGGPa2/ncxkaZZGws89y82ZiVK5lXQIAxPj6stauL8QcGmMP583wXGGiMkxPzGx7WfvbtY+5DQ8B9wgSd79gY/Zw8aYynpzEhIXw/eTKE7uJiTFsbsBseNiYpyZiiImMSEoz5+GNjwsNZz7e+xTwmT6bPTZuMiY/nXW9vY3btMmbmTMYvLaW/lBRjCgr4XVrK+BUV7EF2Nv/n5wOb5GTmYmNjTHW1MdOmMcbKlcDExoY9M4Z5NzSwB83NwCQw0JgZM2Cca9cCo+xsYxYt4reXlzGFhcYkJvL7qqtYn4eHMcHBxpw6Zcy11zLf0FD243Icramhj/x8xnJ0BN4NDYqTWVnMb/duY15/3ZhbbwVma9eiyDg5Md7ixcp0q6qAu6Mje/zii+xXcDDwDAlhvgkJwEfoJDIS/Fm8mLnv3Qt88vJ4NygIvMvOZn9aWowZHAS31q5l3IgIvh8bM8bXl70U+k1PZ5zRUea0fLkx+/fzva8v77z+ujE338xeh4UBo6Eh5tzWxjqamlhbWhowzc0FRiEhxtTXM47MZfNm9jQ4GBwMCAD/bGyM6e9nb5KSFE+Ebi9vu3YZ8+mnxtjbM+e9e9kXW1tjnJ0/+97eveDZrl3Qm42NKrpNTeCkszP0mZrKHt51F/sWFAQszp9XWouK4v3Nm8HF+nr2PDUVXJK+LRb+nzZN6eHXvzZm9WrmXlwMPIODjZkyhXn5+bG3YWGf5ZXvvQcO7N7NfAYH4YmX72N2NnDLz+f90VFjPvkEGF1/PX14eDDXzEzgUlVlzJw5xpw+rbDx99c+y8qYj/B1Y1TRqKlhHBsb1uLlZcyhQ8bMnw+upqQAV2Og9cFB/i4pATd37wYWW7ZAk++9B69wcABPhCa3bYNWx8ZY4/33w+ctFuUFKSnG9PUZs2GDMY8/Dq9cswY+fjnfv5xXDA6Cn05OxrzwgjHf+54x585BT/v2qSzKzVXef+yYMVOnIlsuE+VfWPtHLUtWZel/oSwNDw+b4eHhv/zf09NjgoODvxRlac8emLCDA4rEtGnG/P73IHpREQpLbS1EdcstxrzxhjG33WbMgQMQQV0d7wsS19SAkCK4/PwgHH9/Zere3rzX3o6wsVgYv6EBQjx40Jg77zRm4kT6yMmBOHbsMObuu2FkFgv9OTvDTKurUSTa2iAUUfgsFghjxgxjXnsNITB7Ngxh7lxjnnvOmBUrmG98PATY2grhzpgBw1mzBiY6eTLEe/IkBJmdDQPu6OBHBMjAAGOXl/NZXx+MMzAQRu7uDoMzxpjOTpiSr68K7YoK4CaM0GIBNtXVxvziFzDLgAAEQEYGDOLoUYTdlVfCDC5dYt2rVxvT2wvjsbVlf5ydWVd4uDEffWTM+PH0d+aMMePGwSAvXKDfiAgYb1MTOLFxI0y5qIh1NDUB77Q01tfQAJNascKYd98Fzh0dMNvz59nn06fp/9IlY2JigM358/Tf0sI6XVzAg8WLjXn/feYUEAC8goLY34EB8Peee4z52c+A36JFnxUAJSU8v3evMdHRjDN+PPNsb2d+AwPG3Hgje1Vezhxzc1nr668b4+oKDUyZwnzKyui3vt6YkREU3+ZmmO7AAHC55x5w3MuL52Jj+T85GUHl6cmcEhKMmT7dmMce4/e2bcYsXYoQiIlBQI0fD4wDAtjja69FeAQGMtfgYOjjpZeMueEGFKOsLOCWlGTMM88Yc8cd4FJEhDFbt6LMvf8+a/L3R6mLj4eWAwLYg7o6lMApU8Cdvj7wp7vbGDc38OOqq6AHe3v2e/Zs9s4Y9sPRkfUPD/N/eDjvOTnRpwg9i4XTvRwmHB3BmbIyBH5uLvDJzjbmwQd5Z9MmYyZNAhcTE5nDrl3QwMSJzMHJyZg//1kV7JUroW97ez1UrV6teOfsjNLj4UE/9fX0ExLCnp89C0y
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of samples between +-100 = 24693\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGwCAYAAABhDIVPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhhklEQVR4nO3dd3wUZf4H8M9sTU9IIRAghI4BJNICqBSNInIglpPzPEFU/HkGy0VR1BPEhhW50/XAAtgLp6InCEIEEaUGQgstEEggpJNs6paZ5/dHyIYlCSRhk9ndfN6vFy8yz87OfDaTnf3uM8/MSEIIASIiIiIPoVE7ABEREVFTsHghIiIij8LihYiIiDwKixciIiLyKCxeiIiIyKOweCEiIiKPwuKFiIiIPIpO7QCupigKsrOzERgYCEmS1I5DREREjSCEQGlpKaKioqDRXLhvxeuKl+zsbHTp0kXtGERERNQMWVlZ6Ny58wXn8ZrixWQywWQywW63A6h+8UFBQSqnIiIiosYwm83o0qULAgMDLzqv5G23BzCbzQgODkZJSQmLFyIiIg/RlM9vDtglIiIij8LihYiIiDwKixciIiLyKCxeiIiIyKOweCEiIiKPwuKFiIiIPIrXFC8mkwmxsbEYOnSo2lGIiIioBfE6L0RERKQ6XueFiIiIvBaLFyIiIvIoLF6IiIjIo7B4ISIiIo/C4oWIiIg8CosXIvIItvyjqLTY1Y5BRG6AxQsRub+SU9CbBqHs5e44eKpA7TREpDIWL0Tk9koydkIREvJECDakl6gdh4hU5jXFC6+wS+S9dvnGo7flI9xhfUbtKETkBrymeElMTERaWhq2b9+udhQiagEhKMfnhpdx1+/XA951YXAiaiKd2gGIiBrDDD/01xwHrACqigHfdionIiK1sHghIrfXde/bSNIdw2zbfShpNwD/MQSoHYmIVOQ1h42IyHuFH/0WD+h+hD20NzaVdgS0erUjEZGKWLwQkdtbYZiAxfYJ6NWnH0otdhSUWdSOREQqYvFCRG5vTeCt2Bv7OCb3kHCLZiN2pWxWOxIRqYjFCxF5jPC9H2CBYREi835XOwoRqYjFCxG5N9mGQLkYGmGH6Hg5NsoDUOkTqXYqIlIRzzYiIveWux//OX07irXhkKccxtSVvbGg00C1UxGRitjzQkRu7VRuLgDAbgxSOQkRuQsWL0Tk1koiR6BH1Sc4fev3akchIjfB4oWI3J4MLRRDEKRTKfjZMAtX/XGP2pGISEUc80JEnkOjRW/NKRQX8zovRG0ZixcicmuBGT/had33CDgpwzDoGjxqfB4WQwj+o3YwIlKN1xw2MplMiI2NxdChQ9WOQkQu5H9yI+7XrYR/3g7AGIDg2GtxXNdd7VhEpCKvKV4SExORlpaG7du3qx2FiFyovMtoLLZPQHnkMLWjEJGb4GEjInJrOVEJmG/3x/BOVwIAepi3ItCSBVT0A/xCVU5HRGpg8UJEbu3LbVkAgKgQXwDAn7LeQKg1Gyi8EfBjbwxRW+Q1h42IyDv52M0YGh2EiEAjACDT/3Kk6K4AdEaVkxGRWtjzQkRu7amjdyJALgHytgDtL8N3Mc9i2/Ez+KkjbxFA1Fax54WI3JcQ8JXLqn/2CVY3CxG5Dfa8EJHbsisCfS3LcHUXA5YGdHC0Hzhthl1WoNPy+xdRW8R3PhG5rZJKG+xCiyv69AA01bur28VqrDXMwqnv56mcjojUwp4XInJ7l3WsvaN0J18bQjSnkFN+SsVERKQm9rwQkdvSFKXjad1n6HT0S0dbRe/JuMP6DLL6PahiMiJSE4sXInJb2sLDuF+3ElHHv3W0ycFdsVnphz0V7VRMRkRqYvFCRG5LCemKxfYJyOkywdEWGeQDAPj1cL5asYhIZRzzQkRuS47oh/n2O9G99xD0PdtmkMvxRLdjkCsOAuAVdonaIhYvRORZSnPw4Ol/olzyB/CY2mmISAUsXojIfVnLoYPduc0vDMeNfWCWAnG5EIAkqZONiFTDMS9E5LYC1j6GdJ+p6HL4o9pGv1C82uU/eLP9fBYuRG0UixcicltSVQkAQNYHqJyEiNwJixcicls7Rr6LuKrFyOkyXu0oRORGOOaFiNzWr+lnUIxAXHlZtFP7tNz56FJ1GDixGOg6QqV0RKQWt+x5ufnmm9GuXTvcdtttakchIpXFhPnBR691agu35aCT7QRQlqtSKiJSk1sWL4888gg+/vhjtWMQkZqEwDUn/o2/WZcD1gqnh74Lm4E7rM+gLIq9LkRtkVsWL2PGjEFgYKDaMYhITZZSxOd8jvtsn9V5qN/w67FZ6Ydsq58KwYhIbU0uXjZu3IiJEyciKioKkiRhxYoVdeYxmUyIiYmBj48P4uPjsW3bNldkJaI2ZrF9Atb5jQf0vk7tNbcIIKK2qcnFS3l5OQYOHAiTyVTv41999RWSkpIwd+5c7Ny5EwMHDsS4ceOQl5fnmCcuLg79+/ev8y87O7vJL8BiscBsNjv9IyIv4BOE15S/IXf0q3Wu52IoO4lrNSnwyd2pUjgiUlOTzzYaP348xo9v+LTFBQsWYMaMGZg+fToAYNGiRVi5ciWWLFmC2bNnAwBSU1Obl7Ye8+fPx7x581y2PCJyfyEn1uBDw5so2HkIuHy02nGIqJW5dMyL1WpFSkoKEhISaleg0SAhIQGbN2925aocnnrqKZSUlDj+ZWVltch6iKh1VVRWQFJs9T4WEtUTqUp37C4Lad1QROQWXHqdl4KCAsiyjMjISKf2yMhIHDx4sNHLSUhIwO7du1FeXo7OnTtj+fLlGDGi/rMKjEYjjEbjJeUmIvez64d3ke7zAsqOTATiP3V6LDBuMl7Y3B7dw/1xrUr5iEg9bnmRunXr1qkdgYhUprdVj18L8OetAYjImUsPG4WHh0Or1SI31/nCUbm5uejQoYMrV1WHyWRCbGwshg4d2qLrIaLWsa/rVIxQPgCuf0ntKETkZlxavBgMBgwePBjJycmONkVRkJyc3OBhH1dJTExEWloatm/f3qLrIaLWISQtShAI+IfVfbC8EC8VPYak9LsBIVo9GxGpq8mHjcrKypCenu6YzsjIQGpqKkJDQxEdHY2kpCRMmzYNQ4YMwbBhw7Bw4UKUl5c7zj4iIrpkOiP62tKqf7ZVAAZ/dfMQUatqcvGyY8cOjB071jGdlJQEAJg2bRqWLVuGKVOmID8/H3PmzEFOTg7i4uKwevXqOoN4iYgupO/J5fg/HAEKegLhPZ0fNPjj9eBnEBgSjge0BnUCEpFqmly8jBkzBuIi3bQzZ87EzJkzmx2qOUwmE0wmE2RZbtX1ElHL6HPqW1ypOQAU/blu8SJJ2OJzFboH+ANavToBiUg1bnlvo+bgmBci75IReT3+K8YCod3qfVwIgRWppyArHPNC1NZ4TfFCRN7l6bxr8YR1BhDeq97Hb+pYgtFiB8pz0ut9nIi8l1te54WI6EyFDQmXNTxW7ob8DzHNsA4VGR2AqPoLHCLyTux5ISL3o8jw1dhxWcegBmcpC+yOVKU7hLHheYjIO3lN8cKL1BF5j5Ks/fjNOgX37JjU4DxHBiRhsvVF2GJvacVkROQOvKZ44YBdIu+xP/0YAMDX1/ei81rsSkvHISI34zXFCxF5j8KwIRhY9R5sf/2uwXn6dqg+XPT51szWikVEboIDdonI/UgalCAAIrhzg7PEWA7hB7/nYTzQEbhuRetlIyLVsXghIs8kFFyuHERJeYnaSYiolfGwERG5nQ6nfsZM7XfQZKc0PFNYTzxrmIWVPZ9rtVxE5B68pnjh2UZE3iMqew0e1y+H9uS2hmfyDcEG3UicCoprtVxE5B68pnjh2UZE3iM/YgS+tI+BEjlA7ShE5Ia8pnghIu+RFXMbZtvvh9z1qgvO17HsAA79+jWK8rJbKRkRuQMWL0TksZ7Hf/CB4U3s2Pqr2lGIqBXxbCMici9CQCNbGzXrARGNKsUAReKujKgt4TueiNyLpRQTfhiIUUZfaOR0AAENzvoPWyIAYF5ov1YKR0TugMULEbmXigIAgF5S4OPXcOFCRG2X1xQvJpMJJpMJsiyrHYW
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Lorentzian by ratio of uniforms\n",
    "def f(x):\n",
    "    return 1 / (np.pi * (1+x**2))\n",
    "\n",
    "u = rng.uniform(size = 10**5)\n",
    "v = -1.0 + 2.0*rng.uniform(size = 10**5)\n",
    "# choose which points to keep\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",
    "# trim the x values to focus on the center of the distribution\n",
    "xtrim = 100.0\n",
    "x = x[np.where(np.logical_and(x<xtrim, x>-xtrim))]\n",
    "# and also calculate the new normalization between these limits\n",
    "norm, err = scipy.integrate.quad(f, -xtrim, xtrim)\n",
    "print('Number of samples between +-%g = %d' % (xtrim, len(x)))\n",
    "\n",
    "# Plot the distribution f(x)\n",
    "plt.hist(x, density=True, bins=1000, histtype = 'step')\n",
    "xx = np.linspace(min(x),max(x),1000)\n",
    "plt.plot(xx, f(xx)/norm,':')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$x$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5d2b12a",
   "metadata": {},
   "source": [
    "### Gaussian by ratio of uniforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b7be41ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9dXid15XovY9ssUnMjEZJJpFZMsRJzEkDhjRp00CnEGo6nU460zZpmhSnoUnSmMIxhmxLMsQWGGWZxMwMFtvS+f74zcpyeqfTmXubm6+3Zz+PHknnvO+GtRfttRdYrFar1diardmardmardmarf0dNruvegK2Zmu2Zmu2Zmu2ZmtfVbMpQrZma7Zma7Zma7b2d9tsipCt2Zqt2Zqt2Zqt/d02myJka7Zma7Zma7Zma3+3zaYI2Zqt2Zqt2Zqt2drfbbMpQrZma7Zma7Zma7b2d9tsipCt2Zqt2Zqt2Zqt/d22sV/1BP7/3kZHR01DQ4MZP368sVgsX/V0bM3WbM3WbM3WbO2/0axWq7l69arx9/c3dnZ/3u5jU4T+QmtoaDBBQUFf9TRszdZszdZszdZs7X+j1dbWmsDAwD/7vU0R+gtt/PjxxhgAOWHChK94NrZma7Zma7Zma7b232k9PT0mKCjoczn+55pNEfoLTa7DJkyYYFOEbM3WbM3WbM3W/sbaX3JrsTlL25qt2Zqt2Zqt2drfbbMpQrZma7Zma7Zma7b2d9tsipCt2Zqt2Zqt2Zqt/d02myJka7Zma7Zma7Zma3+3zaYI2Zqt2Zqt2Zqt2drfbbMpQrZma7Zma7Zma7b2d9tsipCt2Zqt2Zqt2Zqt/d02myJka7Zma7Zma7Zma3+3zaYI2Zqt2Zqt2Zqt2drfbbMpQrZma7Zma7Zma7b2d9v+5hShF154wYSGhhonJyeTmJhoTp069V8+39XVZR5++GHj5+dnHB0dTXR0tPnkk0/+L83W1mzN1mzN1mzN1v7/3P6mao29++675pFHHjEvv/yySUxMNL/97W/N8uXLTXFxsfH29v5fnh8eHjZLly413t7e5oMPPjABAQGmurraTJo06f/+5G3N1mzN1mzN1mzt/3fNYrVarV/1JP67LTEx0cyZM8f84Q9/MMYYMzo6aoKCgsw//MM/mCeffPJ/ef7ll182zz33nCkqKjL29vb/W2P29PSYiRMnmu7ublvRVVuzNVuzNVuztb+R9t+V338zV2PDw8Pm7NmzJj09/fPP7OzsTHp6usnNzf1P39m/f79JTk42Dz/8sPHx8THTpk0zTz/9tBkZGfmz4wwNDZmenp4v/HxZbXTUmIwMfv9ps1qNKS3lu5ISfm7822rlR/4fGfnzfUl/N7574xh/2tf/rmo8OmrMwYPGfPqpMQcOGFNYyP/FxTq/kRHGKC7+82ONjBizdasx16//13P+n6zhzz03OvrF/v6zMf5cP3/6+Z+u60/H/NP9/HNr+p98diPsDx1i/P8OrOSd/wpn/qsxb3z3f4JH/zt49ufg+Odg8Zf6+q/w78bv/3QvS0qMKSpSnP5ze/Pn/v5zcLhxLX+6vht/3zh+UdF/jmf/GdxGR5nvoUNf3Ov/7L3/al+lr6IiY15/HTovKvpfYSl9jIz89/boxnH+q3f/K1w8eJAfee9P8eO/yxv+q3H+p5//6Xd/iltFRezJ8LAxzzxjzJUrCss/3f+Rkf+Vvv8SrP4nPO7P8af/Hbj9Z+PdyNP/HB/8KtvfjEWooaHBBAQEmJycHJOcnPz550888YQ5duyYOXny5P/yTmxsrKmqqjJ33323eeihh0xZWZl56KGHzHe+8x3z1FNP/afj/OQnPzH/8i//8r98/mVYhDIyjJk+3ZiLF41ZulQRSZq/vzE5OcYMDhrj52fMs88ac+edxoSHGzNuHM8XFRkTFGRMQYExK1bQV1qaMVlZxixZYkxZmTHV1caEhBjT329MXp4xYWH8VFUZ4+1tjKurMRaLMX199FlQYMzmzcbY/YmaPDpqTGamMcHBfBcZCSLn5BiTksI4V64Y09JCv5WVxoSGGuPsbIyTE/Pbto21fvSRMXFxxkyZYkx09BfHeOIJY775TWNefdWYu+9mbq6uxkRFAR9fX2Nyc+l7ZMSYY8eMWbjQmLo6Yzw9ed7Z2ZiaGmBQUUG/fX3GnDtnTGurMcuXM96hQ8YsW2ZMW5vuwdatxri5MWc7O2BnDOszhrU2NbH+sjL93denzzg58dvLi75DQ9nDXbtY941zDA42JiDAmMZG+iotZc5DQ8CmtpY9zsszZv16Y5qbgYXVquNnZBjT0cG4U6cqrPz86Dcqir3q72dci8WY7Gz25NIl1i7txn5LS1lXSws4ExXFM9u2AbdDh4xJTqbflhb69fJSJv3448aMGaN9l5SgIBtjzOTJrO9GBlpdDawEJ8rK2OPiYtYWGspY588Di7w8PouKYuzr1415/nljVq82JjZW+wgPhyasVub30UfG3HILeFVdrThtjDG9veDAkiX07erK53194N3goDEODsa4uEAnZWXQgOyNKAuOjsxJaEzWZLVCR4KrQi+CD0lJjOPtzXcpKfzv6WnMH//IPjg7GzNnjjEvvmhMYiLrDgsD1gK/0lLdl4EB3u/tZUzB6YEB1hEVxToqK42ZMQM+IjzDxUV5jb29Me+8Q1/9/eCyvz84tngx427bBn0dPGjMunXG7N5tzKZNxpSXM6bwjdxc4FdRAU3n5bHXHR3GfPABaxWcCgsz5rnn2NMpU4yJiVGcOnTImMuX6WNw0JiEBOhzeJh11dWxvwMDPO/sDO+zWlnjjXQVGcneJCYyP4vFmPR0cKOkBBw4f571ZGYa09DA/gQG8rfVynoaG43ZsoU19/bCU1NSgNnJk9BlQ4Mx7u7Mf/Vq+ktLY04tLeDpuXPGeHhAo3197Pm4cexxVRV7deAA8xeeKHvZ2ws+OTp+ERZLlhhz+DB4Jnzs0CHg09EB/ITn3ihjXFzob/Fi3h8dNcbHR3lzZiZ95uTAu2JidJ+2bjUmNdWYV14x5qmnlI6CgliP8JW/dvvvWoT+n1aEoqOjzeDgoKmsrDRj/oMb//rXvzbPPfecaWxs/E/HGRoaMkNDQ5//39PTY4KCgr4URUgUi5AQRaSBAZCjtZVnrFaQ7ac/VaGYng6RlZaieBQWQkh5eTClhgYYRnExyDY0hNBpbIQJdHSA6H5+/J48WQXl7t30fekSzCcykjlkZTHfceOMuXABRnHokDHjxxvT3c0cFixgrM5O5j48DCPYsAEiamyEMb74IoKkpwfmsGwZ8zAG5tnWxjyefJK5hobyvRCLCBFjjNmxA2WyqcmYhx+GgQYEMNeuLhj78uUwi8ZG4NvVBbH7+yNcS0uNWbWKtVVWwlRbWmA8ycnA0WJhrU5OwGTpUj2hOzjAKIODVXG1WHQvLBZjFi2CAa5eDTNYtcqYvXt5p6qKvUhNpZ/+fhiFmxvvODkBqzvvhMGJ8K2oYKyODuZZUsKzAs8/xa9Dh1hfWBg4VV3NXkl/VVUKa19f9iApib09csSYjRthgGvWsLYXXgB3HB15d3SUebS1GXP2rDGzZhlTX2/MPffoCbGqCrg2NaGchIezj/39ykAdHcFxmccLL7AGR0fWGBPDd52dCOLaWuYRGYkCvWoVwkrYRHIyeOHsTB8nTjC3wUFwpbOTPV61Cqa+axfwuXgRRVEU5O3b+fzkSf6fOJH5TpzIXg0MIBBra9mvq1cRVI2NxkREsC/GoCB6ekLjYWH0NTAAziUlMb6fH2O5uIATEREIk64u9v3mm415/31jpk0z5u23oQFXV+ipq4vng4OBlyhgJ06A27Nn86yjozGnThlz223sU1ISv3Nz+bu6mn25csWY9naE8OnTCG/Bg95e5h8ezhoCAvjd2Qkf2bYNnDx7lkONnR3rzMtjjLNnjZk5U/nc2bOsPzwcPJk/35j4eGPeesuYSZOM2b+fPV60iP0uL0cJPHyY9QQG0k9+Pniya5cxc+eyRxYL9Ojnp1bqCxeMufdeYz780Jj77kPQe3nxfl0d4+zbx+EiJYVD1/Tpxhw/zng
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGwCAYAAABhDIVPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABINUlEQVR4nO3de3yT9cH+8U+SJm16LpS2EI6CoIBQ5SSCAtqJeJjInOjcBrjhYcWJRTfY9ujj89uGm8q6aSY+PkM257Tz7GQyBAU8MRCtCgywnCnQFuiBhrZpk/z+KBTKsS1pv0l6vV+vviDJ3eaCQnL1vr8HSyAQCCAiIiISJqymA4iIiIg0h8qLiIiIhBWVFxEREQkrKi8iIiISVlReREREJKyovIiIiEhYUXkRERGRsBJlOkCw+f1+9uzZQ0JCAhaLxXQcERERaYJAIMChQ4fo0qULVuuZz61EXHnZs2cP3bp1Mx1DREREWmDXrl107dr1jMdEXHlJSEgA6v/wiYmJhtOIiIhIU1RUVNCtW7eG9/EziZjy4na7cbvd+Hw+ABITE1VeREREwkxThnxYIm1vo4qKCpKSkigvL1d5ERERCRPNef/WbCMREREJKyovIiIiElZUXkRERCSsqLyIiIhIWImY8uJ2u+nfvz/Dhg0zHUVERERakWYbiYiIiHGabSQiIiIRS+VFREREworKi4iIiIQVlRcREREJKyovIiIiElZUXkQkdAQC4NkPZbvA7zt2v9cD1RX1j4tIuxcxu0qLSAjz+9lXuJX91k4NdyVvyqPTlleJvmgiXHoPAIWllbj+0BuAjbd/Sl1sGilxDlxfPQPLHoGLvk3hVU9S6vEC0PnDnxGwxVB76Qw6d+3V5n8sETEjYsqL2+3G7Xbj8/nOfrCItB6/D7yVEJNUf7vmEP7H+5JRe5irqv+EBycAM2yrecC+ir3RXTnQrZwDHi93P7+Wz612LAS4/U9rOEASTruNVZd7SQIqLXFkPbGCqlofVvxsin4Ru8XHlV9czMPfi6djnIPkzS+TtuE5HJmTYdSPzf09iEir0SJ1IhI8ny6AJf8FF90MN/y+4e7a3/QhcPggH457hbTzhwBQVbiOl95ezFe1LjYHugHgtNuYf3smHRPqC05BcSUz8/J5+97RDExzsKGwlGvnf07u5EzO72gnZeOL+A7uYML6q6isrX+u/45ayNSoJey68IeUj36o/s6An+7/uoPqjgPYP/hu/I6E+jM6yc62+7sRkTNqzvt3xJx5EZG2U1h6GP8XeSTsWk7RsAepTagvH4nVsXT3VnJ41xdsLSw/dvzlL/Gjt/bw5vlDGOg6ckbGNYoufS9puAQEnLZQFBRXAvEUlNb/rNUnLZ4BriTofh8A/yqravg6h0q6kf3qILbkp7Lx8w8BON+ym3ej3yNq50dctmoYdUThtNtYPjWd9E4ZkNg52H9FItKKVF5EpGl8tWCzU1hWRda8lfzV4qab9Wse29iBF3xZAMRjo6vlUTbv7Ir/yQ8bfbrT7iAlztHoPley84xnP1LiHDjtNmbm5R/3dWxn/jquQXTveX6jUmSrLqVwO9iqD/JG5tiGMzqWxT+Dko8pvPw3lF5wm87GiIQJlRcRObPKYliUA7vXwswvKfV4qar14R91J8W+HUzt/U1u63DBWb9MS4qBK9nJ0lljmnR25sTPa3xMEvS+C4CMI18j1m5hW9FBUi0WvrPUzo53P8Rpt/GXGxLIOLyZil7XkpyUqDIjEoJUXkTkZEfOsgAU1kSTvv0Toqr2s+2zZRREXQSA85LJpLmSSGvlKGc7O9PSr/nurHGUekaxqeoAbmfHhgHD69/KZVjUEl5d9gq/4F7mf28IHY8706OzMyLmqbyIyDEHt3H4nf8icKiYbTe83PCGPso3hd2BTmx8tQ7IP+Wlm3BzrBQlNdy3dNYYrP/egPerDfQZfie8C1MWrCYRD8OsG3nPfzExdjtLZ41RgRExSOVFRBrs9fhJ3fxP7BYfdz/1GrsDaTjtNr439Z52cfbBleyE8TnwjfsYbLGydFA1pR4vHdf9ic6fPMHejHGM3D6dUo83Iv/8IuFC5UWkvfL7Obgmj5qDuzkw6E4ACvbHsKpuGld/4zrmH5nSHKlF5YysNuC4szO7nBCdiP/88bD9uOPqaiAq2khEkfZM5UWknSrZ8D6d3rmbmoCdG1d0ppgUAJz2b3Bv5sj2V1jO5NJ7IPM7lBXXwLtrKCiuJK7wQ7quyKHysp+SMuoO0wlF2hWVF5F2ZO/eQg744wEo8PYjzncJXQeMYuGIcfgd9fe3yzMtTRGTRHJiVcPU7fn239HLto9//GsJPTrd0C4uq4mEiohZYff47QE2b96sFXZFjlexl+rXsinato5v1PwWL/UziZx2mwafNlPhkQXxLL4aor/4C1NXd2V3bQIAqZTTy7KXdVED9Pcq0kzNWWE3YsrLUdoeQOQUaiqpzc2Ewwf596j/I3nAVYDOEARD4XGr+7pWPEjK5jxy6yaRdU/usdWEReSstD2ASHvn91Oy4X2KUoY23FV28aP84v0ynhpwld5Ug6hhUG8gAInxBCw2VvoG0bO4suEYlUSR4FJ5EYk0vlpqFlxPp8JV3Ov9Bav8/Y884MRpjw/79VlClsUC18+jaNDd/OfZrQ1bGlxp/YxaWxyPzrpHBUYkSFReRCLE8Zcv4uw9SAt8zoOXxhM9ZHTDMToD0Poyuvdl6axulHq82KoPcl5eNtHeUrb/pyeMnGQ6nkhEUHkRiQDFmz5h4l/3UlJbv+ZIIlmk2q/g+Su+pbJiQMOlpGo4eN41bNvwIT7X5aZjiUQMlReRcLf6WTq981MeDIzGMflp+qRpynPIiEliz+W/4eb8ZTx6wEvAVg6BAB02/AUG30qX9NbeGUokMqm8iIS79AEQ8BNj8XJeaowG44aYlDgHFntswxiYb1lX8oRjPts++iOF967ClZpsNJ9IOFJ5EQk3gQBU7IEkV/3tHpdRMGkxP/7bAd626r90qHElO1k6a0zDeKTYfQ48y/7By+WjubbGgstwPpFwpFc6kTBwdDCu1VuBa8UDxBZ/RsGkxficqQAU+LoBB8yGlNM6toM14PoGGzq8y/xnPqPvkenUtsMlpDh8ZPToZzClSPhQeREJcYVlVWQ9sYKqWh9OqnnT8RU9LQeYt+CvvOs/to6L027TNOgwkZTSgWi7nZl5+Vjx8xf7XLpYt7N/0gJSB19jOp5IyFN5EQlxpR4vVbU+cidn0ictHsvBbuz0VXNfp8Hcd9xxGqAbPo6/lGSrLiX9bRv2g3UU21JJNR1OJAyovIiEMq+Hru/NYLz1PPqkja4fjOsabjqVBMGxS0lJrL/pNeY8/SJ31HamprC8/oCAn5T4GBVSkVNQeREJMccvNpf6xXwytrzFb+xx7PHeCWgmUSRKTozn66i+DTOSelr28oz9d9wfuIffzfqBCozICVReRELI8eNbAGxcyO/sI8njGn6b0sFwOmktJ85I6rb0LpK27Wam7wVKK7+r8iJygogpL263G7fbjc/nMx1FpMVKK2u4xJfPt2/5Ln3SE47cO4bfajxLxGs0I+nWZzn4+k+4L38UPy/x1O+bhMY1iRxlCQQCAdMhgqk5W2qLhJoDL99Hx/UL2Tf8Z2Rc+1PTccSgE8/CjbZ+xX9s/Xhr1jUqMBKRmvP+HTFnXkQigTexB/6AhYDFajqKGHb8paS4wo/o8c5jbPS5qDhwCa7kHqbjiRilV0iREHJgwB1c7/0VBwbdaTqKhABXspOBriR6udLxxySzNdAZX7TOKIuovIiYVLSewy/dwfqdJawrLKegxMOGQE/TqSTUuIawZeI/mFV7N+isnIguG4kYU1dD3fM3E1u5h6XrAvyu7mZAK+XKqdXGu6hhGwVHthRI+/QJapLPwzp4ssbASLuj8iJiSlQ0u0c/yt5Fj9Lnhgd5u2v9Fn2aUSKnkhLnwGm3MTMvn1HWr3jB8Xv8AQsTl3qZ9b2b6Hik8Orfj7QHKi8iba22GuwxAFR2G8tttTbe7uqqXz1X5DQarQUTuIwDn+ylwpbC1593Z8qC1Q3
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGwCAYAAACKOz5MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVqElEQVR4nO3de1xUdcLH8c/MwHCViyIgI4qKaeaFEiXL1IoN225WW+i2aW6P7W7lk4tl2W7abR+qtZYtTavd7hctd7d2raikdLNIS6My08Q0FeWiAiMjzMDMPH+QGHkdBM4wfN+v17xWDmcO31li+HLO7/x+Jq/X60VERETEj5mNDiAiIiJyPCosIiIi4vdUWERERMTvqbCIiIiI31NhEREREb+nwiIiIiJ+T4VFRERE/F6Q0QFag8fjYdeuXXTp0gWTyWR0HBERETkBXq+X/fv3k5SUhNl87HMoAVFYdu3aRXJystExREREpAV27NhBz549j7lPQBSWLl26AI0vOCoqyuA0IiIiciLsdjvJyclNv8ePJSAKy8HLQFFRUSosIiIiHcyJDOfQoFsRERHxeyosIiIi4vdUWERERMTvqbCIiIiI31NhEREREb+nwiIiIiJ+T4VFRERE/J4Ki4iIiPg9FRYRERHxeyosIiIi4vdaVFgWLFhASkoKoaGhZGRksGbNmhN63uLFizGZTEyYMKHZdq/Xy5w5c+jRowdhYWFkZmayefPmlkQTERGRAORzYVmyZAk5OTnMnTuXdevWMWzYMLKysigvLz/m87Zt28att97KOeecc9jnHnroIR599FEWLVrE6tWriYiIICsri7q6Ol/jiYiISADyubA88sgjTJs2jalTpzJo0CAWLVpEeHg4Tz/99FGf43a7ueaaa7jnnnvo27dvs895vV7y8vL44x//yGWXXcbQoUN5/vnn2bVrF6+//rrPL0hEREQCj0+FxeVysXbtWjIzMw8dwGwmMzOTwsLCoz7v3nvvJT4+nuuvv/6wz23dupXS0tJmx4yOjiYjI+Oox3Q6ndjt9mYPEenAvN7mH9fXQZ398O0i0mn5VFj27NmD2+0mISGh2faEhARKS0uP+JxVq1bx97//naeeeuqInz/4PF+OmZubS3R0dNMjOTnZl5chIu2kpKqW9SXVzR4lVbWHdvB4YPUTMD8d9v/o5/3bfHggGf4yGFyO9g8uIn4nqC0Pvn//fq699lqeeuop4uLiWu24s2fPJicnp+lju92u0iLiZ0qqasl8eCW19e5m28OCLSyfORZbTBiYTLD+H7C3GD57Gs69s3GnuurG/x1yJVgjKKmqpdLhanac2Ahr4zFEpFPwqbDExcVhsVgoKytrtr2srIzExMTD9t+yZQvbtm3jkksuadrm8Xgav3BQEJs2bWp6XllZGT169Gh2zLS0tCPmCAkJISQkxJfoItLOKh0uauvd5GWnkRofCUBxeQ13LfmYSofrUGE59w+NhSXtl4eePHwKDL4SzJam4hNVX8G0oDf5a8OV7Ce8efERkYDn0yUhq9XK8OHDKSgoaNrm8XgoKChg1KhRh+0/cOBAvvrqK4qKipoel156Keeeey5FRUUkJyfTp08fEhMTmx3TbrezevXqIx5TRDqW1PhIBtuiGWyLZpC1nPdCbqPb+r8f2qHvWBhxPQT/pHiEREJwWFPx+UffZfxP0Nt8kvIUeVcPo7befdhZFxEJXD5fEsrJyWHKlCmkp6czcuRI8vLycDgcTJ06FYDJkydjs9nIzc0lNDSUwYMHN3t+TEwMQLPtM2bM4P7776d///706dOHu+66i6SkpMPmaxGRji3q+3dJNFVSt3ExnHcjBIee8HPr0yYD5URM+Cup9V3aLqSI+CWfC0t2djYVFRXMmTOH0tJS0tLSyM/Pbxo0u337dsxm3+6WnjVrFg6HgxtuuIGqqipGjx5Nfn4+oaEn/mYmIv5vz9DfsPCjEjIy/ode5U7AecT9jjQ+xWEbDSN+3ngZqaS6HdKKiD9p0aDbm2++mZtvvvmIn1uxYsUxn/vss88ets1kMnHvvfdy7733tiSOiPgxU0MdEA1AbGQIr5p/znOv7wB2HPU5Rx2fYjI1/bOfqYS4LxaB7fY2SC0i/qZN7xISkc7NjIc+b2ZD0mmQ9SdsMTEsnzn2mGNPistrmLGk6NDA3COw1FWy1HoPsWtqoEdPOP2atnoJIuInVFhEpM1kmL8hrLwIqr9rvGU5LAZbTNhJ39njDo3lbw0/5ybbZsJPyWqdsCLi11RYRKTNFHpOY+slS+kbdgCiba167AXuy/j5xSM5LaL15ngSEf+lwiIibepA4giwRbfBkU14LT+aj8m5H0J095BIoPJ58UMRkeP69h0sdZXt87XcDbDiAfjLabBva/t8TRFpdyosItK6Kr6FJdfSf+n52Kho+69ntsD3HzdO5//FK23/9UTEELokJCKty+2Crn2pDYmnpLJtx5cUl9cAYE2/m7hTv6XriOw2/XoiYhwVFhFpXYmD4Tcr2fn9Lije0CZfIjbCSliwhRlLipq2hQXHsPyUOq0tJBKgVFhEpPUFheAO7dpmh7fFhDWbz6XZ3C2RFtj+MfQd12ZfX0TanwqLiLSO4gKoq4JBl4OPy3O0xJHmczG77LDoZ7DvO5jxJUQltXkOEWkfGnQrIifP44F3/gBLfw1rnjAuhjUKwrtBaDRUbDIsh4i0Pp1hEZGT56mH0y6HIgcMm2RslssXQUR3sIZTUlV72DIAR1pYUUT8nwqLiJy8oBAYdzuMubXxNmMjxfYGoKSqlsyHV1Jb72726aMurCgifk2FRURaj9Fl5UcqHS76NWzmt1dcQIotETixhRVFxD+psIhIy7kb4J07Yfh1kDDI6DTN9Fj1B5aFvEDpHgeJI+8wOo6InCQVFhFpufVLGwfZrl8KOd80Xhoy0MGJ5ACcoacS7TUTVLvHwEQi0lpUWESk5boPhCFXQeJQQ8vKkSaSC6InKcF/5bnzrzIsl4i0HhUWEWm5pDS48m9GpzhsIrmDdEeQSOBQYRGRgHCkieSa1FZCTTmQ2K6ZRKT1aOI4EfGdYy+s+gs42m58SHF5DetLqpuNS2mRjW/BwwPh3//bOsFExBA6wyIivvv8BVh+N46v/sPWS/91xF1aWjSOvLChhdgIa4uOh+0McNeDqwZzvaNlxxARw6mwiIjP9oYms9vbl+d2Due1x1Yddb+WFI0jjUc5qbEoXRLh5k+ha188u+wtO4aIGE6FRUR8trtHJhc7Q8m7eihTEqKOul9Li8Yxx6O0RLd+rXcsETGECouInJAfr8tz8HJPakIUg23RRsbyjcdNHNVGpxCRFlBhEZHjOrguT2x9GUPMW1nuOYOwYGvLx5UYYdtHnLL0BuZbuwAXGZ1GRHykwiIix1XpcFFb7+aNwV9ySvHfqOp3GY5LnuhYc5zE9ibYsYsBpnBK66qADnRmSER0W7OInLio2DgI70ZM+tUdq6wARPdk24UvcaZzPu7QGKPTiIiPdIZFRE7YnrQbSczKAXPHfOtw2M7Gyapmt1xrNlyRjqFjvuuIiHEMXuDwZPx4jpdgGqgniLBgC8tnjlVpEfFzuiQkIsdldu2nr2mX0TFOmi0mjBWTu7O+73w+7/cUedlp1Na7D1uDSET8jwqLiBxX9HfLeD/kVnp+cIvRUU5aQlwckbs+JrJkFQMia42OIyInSJeEROS4gvfvoMFrpi52gNFRTl5sb7j8Ceg1CveBGGCT0YlE5AToDIuIHFf5iFmMdD7OvlN/aXSU1jFsYmNxEZEOo0WFZcGCBaSkpBAaGkpGRgZr1qw56r7//Oc/SU9PJyYmhoiICNLS0njhhRea7XPddddhMpmaPcaPH9+SaCLSRvYRhSckxugYItJJ+XxJaMmSJeTk5LBo0SIyMjLIy8sjKyuLTZs2ER8ff9j+Xbt25Q9/+AMDBw7EarWybNkypk6dSnx8PFlZWU37jR8/nmeeeabp45CQjnsngkhAcR0wOkHb2P0lPVY9wRSLBRhtdBoROQ6fz7A88sgjTJs2jalTpzJo0CAWLVpEeHg4Tz/99BH3HzduHJd
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Gaussian by ratio of uniforms\n",
    "def f(x):\n",
    "    return np.exp(-x**2/2)/ np.sqrt(2*np.pi)\n",
    "\n",
    "u = 0.7*rng.uniform(size = 10**5)\n",
    "v = -0.6 + 1.2*rng.uniform(size = 10**5)\n",
    "# choose which points to keep\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_distribution(x, f)\n",
    "plot_distribution(x, f, y_log=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68b72c8d",
   "metadata": {},
   "source": [
    "### Rejection method with importance sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e24e2fce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance fraction = 0.16539\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9d3iV15XovSWKGgj13jtNjaaGASEJcKO7ArYTxzZ2JrGNHTuT3mYSl2QyExvHSWxj3G2qG6hQhUSXRFPvvRdQBUnfH7+77sKZZMa5N7nc+/ns59Ej6Zz33WXt1fbaq1iNj4+PG0uzNEuzNEuzNEuztBvUrG/0BCzN0izN0izN0iztq90syoilWZqlWZqlWZql3dBmUUYszdIszdIszdIs7YY2izJiaZZmaZZmaZZmaTe0WZQRS7M0S7M0S7M0S7uhzaKMWJqlWZqlWZqlWdoNbRZlxNIszdIszdIszdJuaJt4oyfwZdrY2JhpamoyU6dONVZWVjd6OpZmaZZmaZZmaZb2Jdr4+Li5fPmy8fHxMdbWf93+8f+EMtLU1GT8/f1v9DQszdIszdIszdIs7X+h1dfXGz8/v7/6/f8TysjUqVONMSzG0dHxBs/G0izN0izN0izN0r5M6+vrM/7+/v9Tjv+19v+EMiJXM46OjhZlxNIszdIszdIs7f+x9t+5WFgcWC3N0izN0izN0izthjaLMmJplmZplmZplmZpN7RZlBFLszRLszRLszRLu6HNooxYmqVZmqVZmqVZ2g1tFmXE0izN0izN0izN0m5osygjlmZplmZplmZplnZDm0UZsTRLszRLszRLs7Qb2izKiKVZmqVZmqVZmqXd0GZRRizN0izN0izN0izthjaLMmJplmZplmZplmZpN7T9zcrIkSNHzG233WZ8fHyMlZWV2b1793/7zqFDh0x8fLyxsbExYWFh5o033vhfmKqlWZqlWZqlWZql/f+x/c3KSH9/v4mJiTEvvfTSl3q+urra3HLLLWbJkiWmsLDQPP744+bBBx80+/fv/5sna2mWZmmWZmmWZmn//2t/c6G8FStWmBUrVnzp51955RUTHBxsXnzxRWOMMdOnTze5ubnmN7/5jVm2bNnfOrylWZqlWZqlWZql/f+s/cN9RvLz801aWtoXPlu2bJnJz8//q+8MDw+bvr6+L/z8o9r4uDFlZfyMj//334+PG1NSYsxrr/Fb3pHnSkv15y/1ef1z138/Pm5Mebn+PzpqzBtv8Psvzam83JixsS++8+d9/Hf/G0MfmZms5c/n+2Vg8+f9l5bS3+jofz/WX5q/jPffre2/g8Wft+vHFtiPjf3X6/syMJR+S0u1v5ISPhsb+89zyMr6z5//+Vijo/qcvDM6yhj79hmzf/9f7uMvvf9freevre3P4fJffV5aynz+0p7+tT377+b01/bkL81d1vrn+PZlaOF6XPvzv6/HEVn79bj9589cT/vXz/2v7fn1n39ZHP7zPd6//4u48Nf2/6/t31+a23+1Z9fTkOC4wP0v7cNf4p3/FQ3/VzQosP/zNf05v/pz/Lh+//4SXv05bK6n3S9DK19mr/58vOvlx1/b878Eu//q/7+1Xc9XvgyN/p9oVuPj/+vTsLKyMrt27TKrVq36q89ERESYBx54wHz3u9/9n5999tln5pZbbjEDAwPGzs7uP73z4x//2PzkJz/5T5/39vYaR0fH/9Xp/sVWWmrMgQPGeHkZY29vTHCwMeHhfFdebkxVlTHV1cYsWGDMlCnG1NbyzpkzxgwPG7NxozHW/0Ol8/Awpr6ejS4qMiY21pjp0+mvosKY0FBjcnKMcXU15sQJY3x8jImK4v3xcf5vajLGysqYo0eN8fdnnKef5r3AQPrKyWE+O3YYExNjTFsbfQQGGuPra0xzszFhYcZkZxuTkGBMSwvvlZUZMzDAOsPDWd+xY8zro4+Mue8+Yzo6GH/pUuDi5sb/Dg4KFyGyykrgM3myMQEBjH/woDGzZ7PedetYz/i4MTU1xjQ00E9XF/Ps6jJm7Vpj8vOBWWGhMcuWMV5REd81NRlTVwcsBgeNaW83Ji2N/mUu/f0w5fR0Y86dMyYpiT6srJjz+Lgxv/oV86ypYW1BQcBtcJC5T5kCzMrL6TcsjDGqq41xdwcu6el85uXFnAMDjcnLY4zOTuDq6ooSuX49fQcHGxMSwp41NBizfDlzTU7WMYQCa2oYKzPTmE2b+O3tDcyOHjVm2jTwy93dGGdnY+6/X/FUmo+PMe+/b0xfH7gVFsaPMex/Xp4xq1YZ88or4Ovatcbs3Ml7VlbALj/fmKEhYD5lCu96ehrz8sv0GRjInGpqmL+LizHd3eypjQ3zKyoyJjHxi3tWXg79pKaCJ6mp4M/YGM/Z2fG97ImDA3MKDQXX5Dl7ez738jLmpZfA2337oLfWVqUFb2/WtmkT73t7g09jY8bs3s3427dDQ3Pm0GdzM30tWsS+tbWxnsJC1unuDv27uNBXQgLrGxxkf9vbjVm5krmPj+v4mZnsueDjc8/x+YULwNfd3ZiCAtbQ3g5+ytoPHABnra2Boaw7IMCY8+eNmTqVfbSyYh8/+gg62rePMcfHmV9bG3tuDDDMywOnt2+Hx1y/FzY2fG8MeGZtbcy2beDS+fM819gIfnzzm8wnMBCcWLGC/mtqmG97uzF33GHMxx8b88gjxpw8ydxraujfygpcOXECPLGyov+cHPYoJ4f9i46G1oOD+f7KFQTq008bM2EC+1JcDB4FBwOLvDzdv5gY6CcqijGMUSXTxgZ+fvasMZGR4L+vr9J+UNAX+fP4OPuwYwdzHx835vhx9tQYYLpgAf20tSkfeOEF9svBge9CQ8Gh48dZ+9gY7yYkwKtPnjRm8WJ44OAgNNnZCfxKSvi/owPYW1n9Z55yfZPvKyvha7Nns6ePPQbey/cCm79X6+vrM9OmTftv5ff/lcrI8PCwGR4e/p//9/X1GX9//3+IMvL558bs2mXMpEnGzJhhjJ8fCoQxxly6BBKHh4Ok6ekg3uOPw3wHBkDWuXNBNnt7kCI315ieHjb8scdgUgkJMHl3d5hEZCQENHkyDPimm3hubAzG39xszHvvwdjKy2EqOTkQzaOPMufoaAh4fBwh6OSk47e1GbN5szF79kAg1tZfVEaMgXD9/Ix5911jbrvNmL17IRhXVxB/zRpjtm5FeIWFKZEIYRQUwBAaG42ZN4/fyckQSVwcY1pZ8ezAAIzByYl3wsPpKzOT769eNSYiAsXJy4uxT5ygz+XLjfnsM8Y1BkKMj+fvwkL6rqsz5vRpY77+dQSxtzewsLMDVgEBrG/ePJjOnDl8L0w6MJD9iYlhziJ8xsZQcKKj2dNNm9iHgQHGt7KCOYgCdOwYMNy+3ZinnoJZfPQRTOe999i/e+5hXXZ2CIfaWvbby0uF5enT7ElxsTHXrhkzaxaCa2wMJe5rX4OpGQM86+pUsT18mDWVlxvzb/8GHI1hra6uxrz+ujEpKaynsRFYXrgAbHp7YaJlZazr/vvBr2eeMebWW405dAgaCAoy5uJFYOHpybizZhljawvez5wJvFNSgG1HBzBzdUURWLWKMeLjUSA8Pdn76mr2zs4OeCYkqNJdWAhe2dqyV+PjqkyNjbGvovj094Ofyckw502bEOoBAfwOCTHmZz9jL1pbVZE/c4a5HzrEnG67zZgPP2Tt587xTEAAytTYGLSUkMAzo6PsQV2dMatXI+ScnVFOb72VMTo7eU8UzPh4eMLHHwOb9nZVeFJT2aP0dGN+9CNjfv974PzGG9DR0aPgzbx50KetLTBLTTXmrbcUl0UJEZgawxxKSuhjyxYU5Lg43Yt9+/QE3tnJHDw8gKmPD3Tr4QEOXrkC7Z4+zbtnzxpz113GjIwY8+mnKCzd3cbcfTf9hoUxzsgIa7Wzgw+6ucGHW1oY9+abmVdMDPttawudDQ3xzKefory1tQHvXbugtxMnjHn4YWDq7g5chU5dXaGfuDg+KyiAFx09Ct9vbKQ/T88vKh9LlgD7qirGj4yEt6akQDutrdDF2bPMfc4c4OHsDKz9/KDLWbMYo6c
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGwCAYAAABhDIVPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMT0lEQVR4nO3deXhU5cH+8e/MJJOEhGyEJCSsGhYxQJSt4ALUKIsbbkXrD3F5sbXBSkFbtK3W921ra114tWnRt3WptUq1ii3uoIAoIlsU2TQICIEkLFlICJlk5vz+CAwEAiRhkmfmzP25rlyQmTOZ+zAkc+c5zzmPw7IsCxEREZEQ4TQdQERERKQlVF5EREQkpKi8iIiISEhReREREZGQovIiIiIiIUXlRUREREKKyouIiIiElAjTAQLN5/Oxc+dOOnbsiMPhMB1HREREmsGyLPbv309GRgZO58nHVmxXXnbu3Em3bt1MxxAREZFW2L59O127dj3pNrYrLx07dgQadj4+Pt5wGhEREWmOyspKunXr5n8fPxnblZfDh4ri4+NVXkREREJMc6Z82GbCbn5+Pv3792fo0KGmo4iIiEgbcthtYcbKykoSEhKoqKjQyIuIiEiIaMn7t21GXkRERCQ8qLyIiIhISFF5ERERkZCi8iIiIiIhxTblRWcbiYiIhAedbSQiIiLG6WwjERERsS2VFxEREQkpKi8iIiISUlReREREJKSovIiIiEhIsc2q0vn5+eTn5+P1ek1HEQl7ReU1lFV7Gt2WFOsmMzHGUCIRsROdKi0iAVVUXkPuo4upqWv8i0RMpIsFM0c1KjDHlhwVHJHw1ZL3b9uMvIiIYd46cEVSVu2hps7LU1ekkRW9n/qYFDYcTGb63AJWbi7FW70Eh+VlS+rF/OAfa6mp85Lt+IaBzi186+zG72f+8EiBWfEXyg/UsrPHVfgiYwFIrd5Eavnn0CkLzhwDNJQga/XfcfjqqDjjMhKSOzd8jb2bYecaSOgG3Ycfybrr84Y/O/UGd4eGv1sWOByNdkkjSCLBSeVFRI5nWVBTBtV7oHOfI7d/8QoHNrzDnsyLqew1HgDXwX30fWkEzvoauH+ff9PBu14mZe3/wcgfEzfsPmIiXdz9SgFfR08D4NKDT0NkAs/fOoyzNn5G6uq/8kJ9LmXVt/rLgfX2LBJ9dYw9GE8JyQD8yP0WP3X+HQbdAGeO8Y/0LHfeT7zjAJMXRrIrohtzJg+m97a3yPj4F1T0HE/VxGePlI6XboDKIgonzudg54EAJBS+TtfFd+PoMxauf9H/dX9s/Z1YDvKcdyzfWBnERLpYeOdgMqLrIC4NnK7j/vk0oiTStlReREJQQN8cd6yE7Z9BxjnQY0TDbZVF8PjZ4IyAX+wGZ8Pc/qotnxG34VXeXFvL7+s7AhCFh03RNQ2Pq6vxf9n6mBRI7AHRCWQmxrBg5ijKqmqpenMEOJy8fPFIEpJSGnJ7BlKxbxwbvs7knKOiVfQaz9KvSrj/qnPp0bUrhaVVzH9lFUU9x+KI7c++ogoKS6uoqfNSkXURuDz85pwR3PqvHUx55jMuce7lZld/PiuM4alHF/sPW3liUiirOMDtczfyjVUJwI2uz/lNZB01dV5iwD+CNCVpJR1qdnHeNXeyliymzy2gfuO78EEedB8Bt75zJPAXr7DX4+CaN6C47sjroUNmIoGl8iISYpqaU9LUm+Nxaqtg1XOwZxNc/sSRQyRf/BM+ewpG3nmkvMSlgcMJUfFwsBw6NIx67O4ymic/2895YyYw/6zzASgsreL8f87m3onD6LG7nsLd1QDsGXQH6RNm+Z8+MzGmId8PGt7s+x+dLfsatifl8o+NSxlWWuW/uTD7d0xfV8D8rl3JzkwgKdbNva5hnPfNYPgG+GCpf/+d1/4f8YkxjAAW9DzzUDE4H/gJPUurqJlbwIot+yhLjaPwOy8xfW4BsyflkJUaB8A3u/ox4tVzeHb4cPodFa383Gl0cJVzZp9sasobfmS66qrA4YL4zMb/xgseoFNlERn1v2LWpO+RlRrH7q9XsOr9l6ktdMCQsaf3GooIoPIiEnIOjwgcfuMtLK1i+twCyqo9R974yrfD5g8gLhX6NhzeKdpfR8b7D+Cw6tnU90fUxWUAEB9/Dgm9xrPffQblRRX+53HcsgnLFQVlQFnD7YWROTzlhcvPOp/szASgYcRgb0QX8l7fCmwFGt6Ik2LdLdqvpFg3MZEups8taHT70V/LP4Jzinko/qJ0kq8dE+liaK/kRtvtohMbq+OoPzSiA7Cv/2QyDu0r5Q3/DiuSL6fs1u/hrD+A7/C/meUjM+074PyKzcUZZKXGkZ2ZQPHnKxkT+Qo7Vu/nyy7fafh3LK1ivPdDrrx4FJ36DKNwT+3xr6GInJDKi0iIOvzmCJDBHvDVHblz3evw/i+p7HEJ38aNZG+1hx++sIqZ1sVUWrH849k17OGbQxvHA5NhA/DO0lM+77HFpKlC0ZpDIK0tJq392sd+3RMVnKP39UQF64hrjntcRJf+/Nt3AR9t684rTzb8+3bgIGujnsL10RwYsh6cDaM/zrpq8HX0H6YTkaapvIiEuJ7zJ/FJ9DK+KUmHbg2HJUpThvGt1Y+3Czvz101HDq30vulJOsW6ueQ0nq+pYtKaQtGUQH2d1nzt5hScExWsYx39uJRzrmBwr4s5o9rDlEP3R1YV4Vl+MTE1JZCQCVUNozdpy38L/3wbLrofBk85wVcXEduUF12kTmzPWw9bl5CyaQWQ7b+5PjYdr+Vg/7YvjhyWONCN6bX3M3tSDvMPzenQhNBTa055au3IT+PHJEDfV47brkPpajiwp+Fw32H7S6Dwfeh3KUW10Tp1WwQblZe8vDzy8vL8F7kRCVVNXVsEGn5b7/vSVaThoDN/9N9+8MKfM27TZRSviIEVRw77NDWnQ4LbN1e+wdm+TZA5+MiNG/4Nb92NZ/kz5Bb9tFkX/xOxO9uUFxE7OHwWSnRdGVe7luLA4i/eS/33/1/kuey2EukYiX9ORXr33vxrZlf9Rm4DlssN3c/3f15UXoOvxk3npL58nXIxNVsPTdROiSFt5R/4Inkct71dpYm+EnZUXkQMOnaU5fA1S/462sHIT/9OfVQiV934q4azfgA4ny7ACwGYxCrBp/Co08QPT7KuqUsB7se5yzoymrZ3GXz+J0a7/04UT5gLLGKIyouIIYdHWTx1HsY5V3CAKD70nUNMpIsewy6DfeOI6H0JZ3fpCJEqJnZ2stPEn791GJ0OjbL5R9NqUqDfZZRFZlC74qhT0jf8B3pdCNE6dC72pvIiYsjh67XMH/wF2eue4GByPwqvnkZSXFTDG9T355qOKO2kuaeJ+3UZCNe/SPGOcljxMYWlVbjLN9P7lclYUfE47yrwX1hQxI5UXkTaW30tHKwEGn5jdg6+CYpfJzr7SrLTYyAi6uSPF1tqzaG/pLgo/4jNYMcmfheZwbcHuxCxw0un2EMXz/PVkdQxVocVxVZUXkTaSFNnDaWXLiHlg59Ct2FwXsNcBV9UAkxbedyKxiKn0njE5nx2Vt3Iz//+EcXPfAZALDW8457FfM7j8jsfI6NzJ7OBRQJE5UWkDTS1dg1ATmQRr7t2UbftM7b2LDlyh4qLtFLjEZsE/jXzUn9pTtrwIplLd3Ox71P2HbTIMBdTJKBUXkTawOH5LM+Oi+LMiN1U9ppw6OwRF5M9s/jsYD88r33dqjWARE6mUZnJuINtMZ351Ztfc/WeWnBWgGWRVlFA5/6jVJolZNmmvOgKu2La0YeJCkurGO7YwOhFv8YRFQ/nXAKZnQ8N8Y/wP0bXYpE25XAQcfZlrHhnMUsOncl0sXMl/+d+jJozxhEz+WUVGAlJtikvusKumNTUYaLYyP7UpQ3EnXKmf9FEXY9F2tuxZzL5lq2ldm0E++OziFFxkRBlm/IiYlJZtYdO9cU8k/0FntG/BIeTpFg37g6jwd3BdDwJc0eX5i9H3MHYVSnk51yKfwWl/SXs/nYDJYnn+B+jUUEJZiovIgHgqD/I6+5f0rmwEs4aCINvNh1J5IS2Wl2wIg+Vasui5vU
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSBElEQVR4nO3dd3xV9eHG8c+9N/dmksHIIIS9RQgkEIMoqFGoE0fFUaFUabXCTwy1FQdUO7BOqqJYraOOgtsWESvRoAzZKCAgQWYwCSObkJvknt8fkUAggdysc8fzfr3uq+bk3Huf0yTkyTnf8/1aDMMwEBERETGJ1ewAIiIi4t9URkRERMRUKiMiIiJiKpURERERMZXKiIiIiJhKZURERERMpTIiIiIipgowO0BDuFwu9u/fT5s2bbBYLGbHERERkQYwDIPi4mI6duyI1Vr/+Q+vKCP79+8nISHB7BgiIiLSCHv37qVTp071ft4rykibNm2A6oMJDw83OY2IiIg0RFFREQkJCTW/x+vjFWXk2KWZ8PBwlREREREvc6YhFhrAKiIiIqZSGRERERFTqYyIiIiIqVRGRERExFQqIyIiImIqlRERERExlcqIiIiImKpRZWTOnDl07dqVoKAgUlJSWLVq1Wn3Lygo4M477yQuLo7AwEB69+7NwoULGxVYREREfIvbk57Nnz+f9PR05s6dS0pKCrNnz2b06NFs27aN6OjoU/Z3Op1cfPHFREdH8+677xIfH8/u3buJjIxsjvwiIiLi5SyGYRjuPCElJYWhQ4fy7LPPAtWL2CUkJDBlyhTuvffeU/afO3cujz32GFu3bsVutzcqZFFRERERERQWFmoGVhERES/R0N/fbl2mcTqdrF27lrS0tOMvYLWSlpbGihUr6nzOf/7zH1JTU7nzzjuJiYlhwIAB/PWvf6Wqqqre9ykvL6eoqKjWQ0RERHyTW2Xk4MGDVFVVERMTU2t7TEwMOTk5dT7nhx9+4N1336WqqoqFCxfy4IMP8sQTT/DnP/+53veZNWsWERERNQ+t2CsiIuK7WvxuGpfLRXR0NP/4xz9ISkpi3Lhx3H///cydO7fe50yfPp3CwsKax969e1s6poiIiJjErQGs7du3x2azkZubW2t7bm4usbGxdT4nLi4Ou92OzWar2davXz9ycnJwOp04HI5TnhMYGEhgYKA70URERMRLuXVmxOFwkJSUREZGRs02l8tFRkYGqampdT7n3HPPJSsrC5fLVbPt+++/Jy4urs4iIiKeKbugjE3ZhbUe2QVlZscSER/g9q296enpTJgwgeTkZIYNG8bs2bMpLS1l4sSJAIwfP574+HhmzZoFwB133MGzzz7LXXfdxZQpU9i+fTt//etf+b//+7/mPRIRaTHZBWWkPbGEsoraA8+D7TYWTxtJfGSwSclExBe4XUbGjRvHgQMHmDFjBjk5OSQmJrJo0aKaQa179uzBaj1+wiUhIYFPP/2Uu+++m4EDBxIfH89dd93FH/7wh+Y7ChFpUfmlTsoqqpg9LpGe0WEAZOWVMHX+BvJLnSojItIkbpcRgMmTJzN58uQ6P5eZmXnKttTUVL7++uvGvJWIeJCe0WEMiI8wO4aI+BitTSMiIiKmUhkRERERU6mMiIiIiKlURkRERMRUKiMiIiJiqkbdTSMiUpfsgjLyS521tkWFOnTrr4iclsqIiDQLTYwmIo2lMiIizUITo4lIY6mMiEizOtPEaA25lKPLPSL+RWVERFpNQy7ltPTlnpOLjkqOiPlURkSk1TTkUk5LXu6pq+hoTIuI+VRGRKRJsvJKav1vQzRkjZuWWAfn5KKjMS0inkFlREQaJSrUQbDdxtT5G2q2BdttRIU6zAvVQFrwT8SzqIyISKPERwazeNpIjb8QkSZTGRGRRouPDFb5EJEmUxkRkeZRWQ4WK9js1R8X5xC5fSGjrbuBETW7XWldTof168F+LUT3AyCafDqsfxayoyH++uOvuWspHC3ETtfWOw4RaXVam0ZEGs4w4PBO2L++9vZXLoU/R8MPmce3HdhGp8y7SQ94t9auP7dlErPmMcjZVLMt1nKYmDWPwrK/137dFc/BvJtos/eL49uK9sNTZ8Pr11TnOTFbC8kuKGNTdmGtR3ZBWYu9n4i/0ZkREWmw8J0fQ8ZvoUM/uPPr458ICKr+35K849vCoimOP581u+0MOuE1Ml2DOKv/ANpGda3Zdpg2HO5zA23bx9R+w5j+UJKLMyz++LbDP0DhHrBawWI5vv0/UyDvOxg1HXpd3ORjPUbT3Iu0PJUREalzxtOOe/9L240vQ8rt0G40AKUdh0NAMNiDweWqLgQAV82BgEAIjjr+AtH92H3pG9z/zFIWnPC6/6y6jKvPH0HbE+5m2WdEs//8R6u3ZRce3/nCB+DCByjJLgSWVm+LS4Rf/Q8qSmsfxK6vIH8Xuw6VURJS/RrtLEXEle+ELiOOZ3WTprkXaXkqIyJ+rr6//B90LOJW61rY/llNGakKagvT9x4fF3JMeFxrxYXAMOiccsrmnKvf4YmXXmPhh05KfyoutzsWca/1X9DnMrjxrSa9rW4HFmk5KiMifi6/1EmPyu280um/FIz6K+VRvcjKK+GJt3O5dORg4obfBMUnPOHkIuIhDtpieMc5vNaEZlvf/Q9VQeHYelxwwp4GVmcRoGIh4ik0gFVEuC/gLTocXEmvve8yID6CntFh7DViOHT2ba171uMMsvJKzjiA9NgZjJ7RYcytupKtv1gLiTfXfH6U9Rv6vJUKmY+06KBXEWk4nRkR8UcVR6vHePw0APQvlTfzZp81RKbeaXKwutU322tDBpAatkBwhNR8fIVtBbaKYjhyuPYAWBExjcqIiL/J/Q7e/RUMmwRDbwVgs9GNfRfcQmSkZ166OHm216YMIP1dxW9IHn0LXYZdDgcqAaov25TbqsejiEirUxkR8Tc7PocDW6rn9Bh8i9lpapxpwb3mmu3VwEpxtzHgCAWq77qJWzYDDq6Fsc9D13Ob/B4i4h6VERF/c85vwVkKyb+CAAdg7uRdZi+4F0EJoTkroTQHqpxnfoKINDuVEREPcfJcH8226NyhHfDVE3DlM2C1Vc+3MeoPTX/dZmL2gnuFhJF17Wf0d30Pte66EZHWojIi4gHqmuujWWb4rDgK/7oKCvdCdH8YPrkZ0jY/sxfccznCIP6EIlJ6EL74C1z8sGmZRPyJyoiIBzh5ls9mm+HTHgSX/BmWPgkDrz/z/lLtvdvghy+qp7c/77kGPaWuWWxb8wyPiDdTGRHxIM02y6erqvqSDMBZY6HfFcc/ljO78AEo3AcXzYAGDCPR+jUiTaMyIuJrNn8Iy5+GWz6EoPDqbSYXkTPdKeNxOiXDnSur/387ca2cemj9GpGmURkR8SXOUlh4D5TmwZp/woi7TY1j9p0yTXJCgetp2Ufn/02Cm/4JgW3qfYrWrxFpHJUREV/iCIXxH8Gal2H4/zX55Zp6RsPsO2WahauKF+xPEb77R/jfg3DFbLMTifgclRERXxPTHy57vEkv0ZxnNFryTpkTS1KLXQKy2phacSf/7r6IsItmtMx7iPg5lRERb3e0CD64HS5+CNr3apaX9PQzGnWVJWi5S0Abje7suvQtBoToEoxIS1AZEfF2n94H2z6Gwz/AHcurJzVrBmbP/XE6dZUlaL3CFJr9FZdY1wAjWvy9RPyByoiIt7toJhTshov/1GxFxBs0Z1k6+RLPaUvNzq/o+sktPGW3s2LHRcBA77lLSMRDqYyIeLuwDjD+P2CxNPgpJ07Q5c+/SE93uafe+UE6p+KMT2XRHjv3LTxIOUtrnuMVdwmJeCCVERFvVJwLBbvJbnP2KZcqzuRQqZPbX197ytTz/viLtK7LPWecH8QWQOCE9zinxOC9IxU1mz1pTI2It1EZEfE2rip4fxLGrqU8U/Ub5jndH7cQbLfx2q+G0e6nAuLPv0gbdbnHHkx8FMRH/fRx/m6I7NLs2UT8hcqIiLepLIfQ9hg2B6uPdqs162dD+VP5aNHbfyud8PHd8O3b8Osl1bdVi4jbVEZEvI0jBK79J1lbvmHHv7I
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Rejection method for exp(-|x^3|)\n",
    "\n",
    "def f(x):\n",
    "    return np.exp(-np.abs(x**3))/(2*0.89298)\n",
    "\n",
    "xmax = 3\n",
    "x = -xmax + 2*xmax*rng.uniform(size = 10**5)    # Generate x values between +-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",
    "print(\"Acceptance fraction = %g\" % (len(x[ind])/len(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\n",
    "xx = np.linspace(-xmax,xmax,1000)\n",
    "plt.plot(xx, f(xx), 'k')\n",
    "plt.xlabel('x')\n",
    "plt.show()\n",
    "\n",
    "plot_distribution(x[ind], f)\n",
    "plot_distribution(x[ind], f, y_log = False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bd8720aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptance fraction = 0.66761\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydd3icV5n2nxn13nvvkqtsx5bkktiWbCekJ0BCOiw1wLKUhWUL+8GysNRlWWCBXSAFCJDEKU5xkUsSF7nLXb333vvMfH/88uSMZMklsdOY57p0SZp53/Oe9p77Pk87FofD4RCXuMQlLnGJS1zikndIrO90BVziEpe4xCUucclft7jIiEtc4hKXuMQlLnlHxUVGXOISl7jEJS5xyTsqLjLiEpe4xCUucYlL3lFxkRGXuMQlLnGJS1zyjoqLjLjEJS5xiUtc4pJ3VFxkxCUucYlLXOISl7yj4v5OV+BSxG63S0tLiwQEBIjFYnmnq+MSl7jEJS5xiUsuQRwOhwwODkpsbKxYrXPrP94TZKSlpUUSEhLe6Wq4xCUucYlLXOKSNyGNjY0SHx8/5/fvCTISEBAgIjQmMDDwHa6NS1ziEpe4xCUuuRQZGBiQhISEN3B8LnlPkBE1zQQGBrrIiEtc4hKXuMQl7zG5mIuFy4HVJS5xiUtc4hKXvKPiIiMucYlLXOISl7jkHRUXGXGJS1ziEpe4xCXvqLjIiEtc4hKXuMQlLnlHxUVGXOISl7jEJS5xyTsqLjLiEpe4xCUucYlL3lFxkRGXuMQlLnGJS1zyjoqLjLjEJS5xiUtc4pJ3VFxkxCUucYlLXOISl7yj4iIjLnGJS1ziEpe45B0VFxlxiUtc4hKXuMQl76i4yIhLXOISl7jEJS55R8VFRlziEpe4xCUucck7Ki4y4pL3nDgcIpWV/H4z3ztfV14usn27iN3Oz44dIjbb3Pc7HCIVFfzY7eb+yUmRRx7hXr2mvJzfWs7Metnt3FtWxs+2bfzW++x2rtffM8ux2ajv1NT59Xau59SUyO9+J3L2rMjLL4v87/+K/OY3XK9is1H/yUmRrVtFfvvb6d9rX23dSj1tNtNG55+zZ0W+8x2Rc+d4rpap9Tt3ju9fftnUV/tx61Y+1/L1sxdf5J4zZ/ju1CmRv/kbkbEx03/l5ZT9u9+Zemudt23jO/29ffv0+p87R3+8/LIpq6yM/3/zG9qk95SX8/lvf3t+XzmPi9ZLx9F5zsw1B5znymzzeea9zjLb3HWebzO/n62sud6bueZtefnF37GLyaU+80rL1S7fJZcvFofj3T8cAwMDEhQUJP39/RIYGPhOV+c9Iw6HSFWVSHq6yEVOb37bnvtm6qQLh4hIRgb3R0eLlJSIFBXx+VzfFxaKVFfP/rzKSoAmLMzcf/31gNYHPyjS2kp5IizAO3eKJCQAGvHxLMj+/tx/5IjIZz4jsnevSH6+yJ/+xOdubiKpqSIpKSJ1dSKRkSJ+ftTn+98XKSgA8FasAFxffVUkN1fktttEOjpEamtFPD1FPvQhkQMHqMvkJHVcsEBkaEhk/36Rn/6Uz0JCaFNUFOWJiOzeLbJpE0A9f75IaanIvHnUb9MmkcREkR/8gLo89RT/BwTwfWysSF6eyHPPiYyOUvbQkIiXF32RlUU99+0TWblSpL6etlVViUREUMcXXhDJzBRpaxPp7hYJDhZxd6fP164VSU4WGR8Xeekl+nPFCpETJ+jngACRiQmu27ZNZNEivluzhvG9806Rvj6R3l6RwEDGLDFRZGBAJCeHNhw6JHL6tIiPj4iHh8iNNzJe6enMrf37RTo7qVd2NmNms4kcO8Y1e/aI3HefKfv//o9yvbxENm6kfkeOiHz848y7l18WGR4WiYtjnnR3c31CgsjJk9S5pYX5cPiwiK8v9erpoe2pqczVtDSRXbsoo6JCpKtLZMkSfm/YMH0+b9vGc44eFfnGN2hzUhLj19JCGzdtgjwtWcKYrVzJdxaLGbOYmOnzXt+TmBhzbW0tdT51ivcrM/PS3uPZ3uviYt6XtrbZnzmzLnrfW13XLlS+S66sXCp+u8jI+1jeqRfuQs99M3WqrGRxFzFAvmMHYOfnx4I0PMwi1dXFArljB+Dp4yOyatXci1p5ucizzwIqIoDFffcBAomJlG23A0Bf/KLI44+zqP/hDyLLlrGj7esT+dSnAJb77xd57DHKfu45keuuE/H25v+ODshIWppIczOfP/aYyOc/L/Ld7wKGFgt1XrKEeywWgDo5WSQ0FCLw/e8DWqdOiQQFUWZ4ONc3NrKwZ2aK1NQAmD4+tOuuu0T+/Ge+O3UKEMjIAEjnzwdEExIA8qQkkaYmkZtvBozXrqU9ExMA/9QUYDsxwTNaWwGotDTqcM01tLe93ZQ1NgYonzvHtQEBIiMjAGt9PW2rqoJINTUxfr6+EKzBQZFrr2W8IiO5NjaWNtfUUP+yMsawvR2SOjTEXGtuBsxHR5kzAwOMYV8ffTU2JtLQQP+GhTEuY2M897XX6K/mZsja3r0iVivXr11L+z08RJYuhZxu3AipiI2FGKxdS52zs3nuvHmMR0uLSH8/de3upv3XXks/hocz/j09PPs3v6GPOjuZl/fdR9sefJD5UV4u8utfU49Fi0SOH4eQtLZSz4QE2l5ayphXVEAirVaIhRJkJfIK8s7aAwX92Fj6Yv9+kcWLIY8z77ucdWK2TYN+J0LZs20idA1REnW5z347N2rv1Kbw3SKXit8uM837WNLTzQv7Tj9XF7a0tOnfOZsT5qLF6eks4G1tAI/dfv73vr4s1vn5vPgtLSzqIjwvLe18U0dVFX9nZbFj9PVlgdcFuqICgPuf/wFMvvAFFs0//lHk7rsBEy8vnvncc+bZutN74AF2615ePG98HO1GdTXgc/q0yPLlEJF16yBStbXszIOC2Nm3tHBvZCRmisOHAbgTJwDcoCDKraszwBQXB4jm5vJdfT1gefSoyL/+K89fuBAisGcPALN7twFJDw/6+UMfEvnZzwC4l18W+dznuCYsDM1GRQXPUtLT3Q0BUW3NqVNoGEZHaUNyMqQoOZlyRkcBwf5+QK2sjLLLyynDYuGaqSnav3WrAVFPT+49dQqQLynhur4+yM+pU2gY7Hau9/OjL1pb0TS88grP8PCgfzIzqevUFPNgeJg65eXRJjc35spDD/F/ejrjde+9lPvss8yNEycYy2PH0IAcOgS5qahAozUyIvKXv9D/np7098iIyD33MLahocyZ0lLq88tfMo+OHWOMk5IgsJmZaOqqqhjD66+n/IYGkX/+Z5GDB+m/ggI+E2GsrFbKsFoZk+pq5uRs715VFXUbHaWsjAz6LzOTfggIMEQkJobfag6a+Y7O9V63tTGPq6tNGVVVzEmLZXbNqsNBvzlrc/RdvlTR9rwd5ODN1vGvTVxk5H0sb+cLd7Hn6gtZXT39u8pKQGh4eO6XVRel2lpA5bHH2OF2dhpth8Uisn69yNNPs4gvWsQutqiIRau4GPBRQlJczK5s3z5U3F5eZtG229mxZmYCJNddh/nkuutEfv97kVtuAezS0jCnbN0KqP3sZ5RbUYHWJDUVEPvgB9nZe3hQ3ssvQ5Ta2iAcixahnZiYgEAoUJw9C+iPjwNKERGAbWAgdfjIR9idhoeziy8rY+fd3Czyn/8JoHl4AP4hIZhMqqvNM8LCAL/aWvpxeJif3l7K/d73aMO2bSKf+ASmoG9/mzr19nL/a6+hQeroMBqQoCCISUAA9T15knYfOoTWo6EBTcayZcacVFfHeFRVUWcRxiQ8nHs7OgD53bupU3c3gBscTF8NDKBJCAvj/xtuAGR7eqjr5CTzLyKC5y9YAIkaHISoVVeLrF7Ncw8dgmy2t9N3Hh706+QkROnzn2cMrrkGopqWBlk9elTkn/6JOfaRj0DILBYAvb9f5AMfgEB87WsAdnU1bU5IgMQ89BDkobmZdjQ2cs+jj1K3nBzmf2oq9S4s5NnJyVy7YQOamf/5H8pMS6MPCgshTvn5EC9fX8a7vp7xF+Fv1VQUF/MOOBym79PSpptzRMy750wOdu6EaO3ceXn
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGwCAYAAABhDIVPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLqklEQVR4nO3deXhU5cH+8e/MZCYZErJASAJhJ8oiS2QVRAENAloVV9xxw6XBygvVqv3VpX1bW2uV1s5b7eJS64JaRUWqZREQRPYoGECDYQuQsGQn68z5/REcQLYkJHlmuT/XlUs4OZO5D2Myd855zvPYLMuyEBEREQkSdtMBRERERBpC5UVERESCisqLiIiIBBWVFxEREQkqKi8iIiISVFReREREJKiovIiIiEhQiTAdoKn5fD527dpF69atsdlspuOIiIhIPViWRWlpKR06dMBuP/m5lZArL7t27aJTp06mY4iIiEgj7Nixg44dO550n5ArL61btwbqDj42NtZwGhEREamPkpISOnXq5H8fP5mQKy/fXyqKjY1VeREREQky9RnyoQG7IiIiElRCprx4PB769OnDkCFDTEcRERGRZmQLtVWlS0pKiIuLo7i4WJeNREREgkRD3r9D5syLiIiIhAeVFxEREQkqKi8iIiISVEKmvGjAroiISHjQgF0RERExTgN2RUREJGSpvIiIiEhQUXkRERGRoKLyIiIiIkFF5UVERESCSsisKu3xePB4PHi9XtNRRMJOXlEFheXV/r8nRLtIjXeHzfOLSMvSrdIi0ngHD7CnIJ+L/7GZAzVOAJIoZIwrmwcvG0zbwVf6dy1c+ndq926hqMflVLXtDUBiTR4pG/4GrRLhgp8f/rpf/AXyN0D6TdBleN224p2w6EmIjIPxvzn8dT/7K/Pmfcx7teew3HcWAO2d5cwbvIqYmNZw4aOHv27OfNiXA52HQYez67Z5a2HnSnDFQEo/sNma4R9KRE6lIe/fIXPmRUSaSNEOKNoOyWeBO568ogqqcj4jad0fqYrtyu6Rh4tDj3d/RMr+r+nn/RlXTJpMWlIMJRs+ZsRyD0Wf9WJD+wsB2F9ejXvePxhqy+bxVRF85NsPwAjnt7zueJGq2K5823Oq/+t2+fq/tN6xkJ2t0ymK6ANA1P4dpK37FzWtktjc72f+fRM2zuNa+wJ6DhqO45yR5BSU8ce3/kPMur/idcWysdf9/n1TV7xBwrfvsGfIQ+yzugMQcbCAXq9NwLLZ+fqOXH95SVr1e9pseY+Ic6fCOffWfQGfF1a8AK2Tofdl4HA2wwsgIqei8iISrgq3UrRuNgdr7Rw4azJw6HLLv66CfZvh5vfIazucjD8sZrB3Ha+6lrJ9xw5+tG6p/0u86fLR1xZFjMPLkG5tSI13k3+wB0s/7893+5N59LnD+97pGkqHXkP5nz7juTexL/vLq3n81UKerbmKwv0x/POIfS+x96WLrR0L59Wy6b912xMp5hrHdRwsjuSVI/Ydb+9D74h4ru87kqTUOBKiXVRFxPJ/tZdRXRvBzCP2vdERyzn2c3hvWS0LP6vb3oF9vOpqjw2LH/15mX/fPzizuMqxk+/27OdgXjFwqOh88jBgg1/sO/xvudwD386D9Bug/7WHt1uWzuSINANdNhIJBxv+Dds+h8F3QHLdmYx9WR+ROPsGNvs6Mq76KQDcTgdLuvyD2JJv2D3iCbIihzBtVhb/N7ET/StXUxPdnvIOIw5/3UNvzj8cY/LDMShw/HEox9uvMerz/A1Vunc7f/z3p2yriWc3bYG6S2KPu/7FmLRY3Le85X+u6A/vIn7LB+we9v/Y3/+uujyuCvjjAEjpR95lb1BY2bBjEAk3umwkEq6qD8K2ZXXjQwbfdnj7utdgywJI6nO4vLRKY6V3KF37DmfOyJHsL6/mnlfXMCTn1rrHzAbIwu10MKDXGaTG9693jNR4d73eiOu7X0M1yddN7UfnrmlHlaCcgjJ+PCuBmf3SScsr9v+bda8dRh97CmuXxLFl8VLcTgevj4ezq0qo3LeVjJnLqaipu5ngFxGv0tO2nee9l7HU18//td1OB/NnjFKBEakHlReRYOatgZqD5FW6KCyvJrLwG85452p8jig2Jv0Iy+ECID51PFGtelAc0YOKQ5dAcspjmFYzjTkjR9I3NQ6A+TNG1euMSbj4YQlKiHbhdjqYNivLv83tdPDgbZNoG+1iMvgLzTUfVtHD9lvaVJeCA165fShto12kvf0oUUU5dJrwP5R2GQlA3tZvyJ77Fypzo+Ds0S16jCLBKGQuGx15q/Q333yjy0YS+r74Cyz8NaX9JzN0xXmHfrO3mO36Bd/6OvLr2hspovVJv4R+22+4+tyWfdJ98rNh+3I46wpo1QaAXfP+TIdlP6e8/TlE3/3J4S/krdGgYAkbDblsFDLl5Xsa8yIhqWwvbPqw7g3PnUBeUQW+L2fR6dP72Zs4jCE772fmpHTSkmIa9GXD+axKIMld8SHZc/5Ep/QM7OfcU7fRV0PPN4ZT0bYveaOexutO1OslIU1jXkSC2PEGm6a9PZ6oohx2llpsSb2Me15dQ0RNNJ1tvyZ7ZxfcTof/bh8JPq6eGfx0jpOKVV5YVXcX1BDbJt6OLKC4fDWX/D0bH3bcTgef3pVGSoeu4NCPbwlfOvMiEkDyiiq49w+vcoG1nJm1VwF1t9lmOmYzzrEKT+3lfOIbitvp4PmbB9E2um5Mi34jD37HlFbLIrLwG1yl2yntMpacgjKmzcpibfvf0bpyFztHP0t5x/P12kvI0JkXkSBVXHiAWfaf47ZVc+mVN1ORPKjuE9YIsNm5D7gPlZVQdNw7pDoOBYYCda95srMC34FcLMq59r0S9lJ3Z9OCqQPp0C4R7FquTsKDyouIQbvyC6jZspjSLhcBkFMEa73ncekZLnp0SIKUOLMBJWCkxrt5d8Yl7CkdTdH+bF5qN8B/Nsa94BHYtwYmPAVnjDUdVaTZqbyIGLI7byut/jqCGCoYWfVH9hyaCM3tnMKYy0cTpzMr8gP+szOdzvdvc1GDe9tCqNrPllIHUUUVOisnIU/lRaQl+Xz+U/v7SaDY15X+cQf51xVdqWw3ANAlIam/hGgXDmcUA4t/z2h7FnPfqsDtXFx3+/uWWVCaX7cuU5TG/0loUXkRaQk1FbD8z7DhPZiyEJxRAPykZiqvTBpP305tDAeUYJQa7z5iYsEMLjp0GWlNzm6SFv4a58EC8mpjKOx901GPU0GWYBcy5eXISepEAtLql6Akj52f/ZOinpPIKShjP3Fgd5hOJkHsyIG+/hmA3/maT+zXcZVjCXfN70jt/Lrbr9tQQiExRDmdmpxQgppulRZpLgUbIam3/6/7V7/Dkx9k8e/qYVjUXTrSDLfS1E64KKVl0W3ONVRXHuSa/FuYeu2PjprUUGdjxDTdKi1iks8Hs++Br2bBLR9A91EA7G4/lneq3UfNhKs3DGlqJ1yU8kAuFG7C7a3Bimh11PpMoCItwUXlRaSp2e0QGYtls5O/8XP2RaYDdSsSA6QlxfgXQhRpMW26wdTV2Het5Z/Jo/1nZ6L2rWeDtwvT3vqKwvJqlRcJCiovIk2hYCO7fAkc8Nb94C/qci/PftGdNZ91gc+W+ndzOx0kHJoVV6TFtU6GnhNIpe4MDfu3wAdXkJI0mNbc7i/YoLOCEthUXkRO19pXsT6azoqa4fxP9d3+zW5nd165/fAU/qA3BAkwBRvBHkGE00WtM+aoS0m6jCSBTOVF5HQlngHeGmKtEv54TR96pNTd9qyiIgGv94/gnqVERUQx32pTdynJV0Pejq3c/UG+LiNJwFJ5EWkoy2J33lb22w7NzeLozb7zZ3HHf2uYk9JG41kkuLTtAXD4UtKS39N76Z8YZ78DGGk0msiJqLyINERVKQf/PRXH5sVMrvpN3Twth7idERrPIsHN54Mti3BUl9DadtB0GpETUnkRaQibA/vebNpQzF/Oq6RV+iX+T+kykQQ9ux1umc32ZW/yztwEbv1+u2WBzWYwmMjRVF5E6uHIib/y05/G8/Fafpl+lS4RSehxOCnpcRmwlJyCMmz
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGwCAYAAAB7MGXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRlElEQVR4nO3deXwU9eHG8c/uZjcJCbkISSCGMxwqSrgNHmiNglfFWovWClJL64EVo7ZgK1Rbi4pHflVarNaj1gO1nmhRiYIiCIqi3BrkDCQEyEFCkk125/dHZCGQQDYk+e5unvfrtS+T2ZndZ8z1MPOd79gsy7IQERERMcRuOoCIiIi0byojIiIiYpTKiIiIiBilMiIiIiJGqYyIiIiIUSojIiIiYpTKiIiIiBgVZjpAU3i9Xnbs2EHHjh2x2Wym44iIiEgTWJbFvn376Nq1K3Z748c/gqKM7Nixg7S0NNMxREREpBm2bdvGCSec0OjzQVFGOnbsCNTtTExMjOE0IiIi0hRlZWWkpaX5/o43JijKyIFTMzExMSojIiIiQeZYQyyaNYB19uzZ9OjRg4iICEaMGMHy5cuPun5JSQk33XQTXbp0ITw8nL59+/Luu+82561FREQkxPh9ZGTu3LlkZ2czZ84cRowYQU5ODqNHj2bDhg0kJSUdsb7b7ea8884jKSmJV199ldTUVLZs2UJcXFxL5BcREZEgZ/P3rr0jRoxg2LBhPPbYY0DdlS5paWncfPPNTJ069Yj158yZw6xZs1i/fj1Op7NJ71FdXU11dbXv8wPnnEpLS3WaRkREJEiUlZURGxt7zL/ffp2mcbvdrFixgqysrIMvYLeTlZXF0qVLG9zmrbfeIjMzk5tuuonk5GQGDBjAX//6VzweT6PvM3PmTGJjY30PXUkjIiISuvwqI7t378bj8ZCcnFxveXJyMgUFBQ1u8/333/Pqq6/i8Xh49913ueuuu3jooYf4y1/+0uj7TJs2jdLSUt9j27Zt/sQUERGRINLqV9N4vV6SkpL45z//icPhYMiQIeTn5zNr1ixmzJjR4Dbh4eGEh4e3djQREREJAH6VkcTERBwOB4WFhfWWFxYWkpKS0uA2Xbp0wel04nA4fMtOPPFECgoKcLvduFyuZsQWERGRUOHXaRqXy8WQIUPIzc31LfN6veTm5pKZmdngNqeffjp5eXl4vV7fsm+//ZYuXbqoiIiIiIj/84xkZ2fzxBNP8Oyzz7Ju3TpuuOEGKioqmDhxIgDjx49n2rRpvvVvuOEG9u7dyy233MK3337LO++8w1//+lduuummltsLERERCVp+jxkZN24cRUVFTJ8+nYKCAjIyMpg/f75vUOvWrVvr3QwnLS2N9957j1tvvZVTTz2V1NRUbrnlFn7/+9+33F6IiIhI0PJ7nhETmnqdsoiIiASOVplnRERERKSlqYyIiIiIUSojIiIiYlSrT3omIsEpv6SS4gq37/P4KBepcZEGE4lIqFIZEZEj5JdUkvXQIiprDt5DKtLpYMFto9qkkBxehEBlSCSUqYyIyBGKK9xU1njIGZdBelI0ebvKmTJ3JcUV7lYvBA0VIWjbMiQibUtlREQalZ4UzYDU2DZ9z8OLENCmZUhE2p7KiIgEJBNFSETM0NU0IiIiYpTKiIiIiBilMiIiIiJGqYyIiIiIUSojIiIiYpSuphGRZtMsrSLSElRGRKRZTM/SKiKhQ2VERJrl8MnJNu0oZPZruezflgBxww+u+MXTsK8AMq6C+B51y3atg5UvQFw3GD7p4LpfPkfizq30tCUfXFa9j4iir0lmb5vsl4i0PY0ZERH/rZtH0hcPEU+Zb3KyQfsWMT98KinL7q2/7rI5sOg+KNl6cNne72HJ3+CbufXX/eIpUj5/gJ62nQeX7fiK9Dcu4QXXYa/7zcvw5XNQmt+y+yYibU5HRkTk6GqqcJV8X39Z7t0k7f6WU+y/9y2q7ZDMHqsjYc4O9dc9+TIo3wUdux5cltALRt4MsWn11+1/EXs79mX7N50PLvO4qYlKYUdZJ+IPXXfJ36BgFVz5AsSm1i2r3gfu/dAxGREJHiojItKoDgWfwzPj6d4hBfjzwSdOvITiwu3sXdXRt6g8bRRnVz9OzqkZpOeXHly3zw11/60GfMu7wsl31H2YX3pw4OtZt7Mjv5Rvv158cPv0LDb8fDnXPLqYnF3lvsXJyWcSaY8hn+7U/PC6qVteJ/79W6D/xXDl8y33P0JEWpXKiIgcVOuGqhIgHICqhH7gqcZeW0Ec+w6ud+508vNLWf3NwdIQH+Ui0ulgytyVfr/tsQa+NvzaZ9Q9vt8CbAHgD64PmWQHumTUf4H9e6FDgt+5RKRtqIyItEOHX5ILELPpf3RZOoPKxIHkZTwCgNcVAzctZ0NlJ0oe+/Sor5kaF8mC20Yd8brH0pQ78jblteteB8687i76pyUdfCJ/BTw1BgZdAxfOArvDr3wi0vpURkTamYYuyQXobytmnmsXFeWfc+eGz4h0RhEf5YK43oecXjm61LjIZl/Wm/fDKZi8Q07FNOe1azskQUTd3X7zSypxrniLJI+bkuIitu8s11woIgFIZUSknTlwSe5jP+nNSfYt7O9ymu+5LYX9qEocwMuO8Db7o93QKZhIp6OuCDXDgTKzp8LN9c+toLJmOKfZ/8imNSkUrllMpNNB7s1D6erZAV1ObYldEJHjpDIi0g71s21l9Ee34PRWw03LIaZL3ROpP2rzLA2dgmlOEWqs1Dz7y+F0ijoDOHhKyLn4AVj1BIyeCadd3yL7ISLNpzIi0g5tsrrgjumGs7YUKooOlpFjOPQUSmOnU5rjeE7vHPoaTSs1Fs6KQrC8EHvCcb2niLQMlRGRENPo/WI8NeBwAuDGydbz/8WJ3bqAM+KYr9nYlTLHczqlNTSt1NjYfs7/Efej30LqkIOLLQtstlbNJyINUxkRCSGN3S/mw+v70+Xd6yDj59D1pwB4IhKaVESg8atZgnYwqM1Wv4hUlcF/Loezfw/pWeZyibRTKiMiIeTw+8UcGCNhX/EM5H8BJVuwX3FBs167JU6lBKxPc2D7cnjjJvjtV+DqcMxNRKTlqIyIhKAD94s5oGjwFJIT4uDES/BWRZkLFqhG/R4qS2DoL1VERAxQGREJVZUlYP3wsc0GZ0yp+7iJc4aEskMH3/pONV38cP2Vat0QFjjjYURCmcqISAhyVO6Gf11Jl6RMbJxnOk7AaOzy3yOmot/9HTz/U7joIY0hEWkDKiMiIShq52ew+zs6VpaTwGnH3qCdOHwgbqNT0S99DIo3w4d/gV4/ArvdTGCRdkJlRCQElfW6GOLC2exIZ8+/803HCShNGoh7wSxwRcPpU1RERNqAyohIKLGsgx+f8lPc+aWAyojfwlww+l7TKUTaDVV+kVCx6lW6vTeRGFpuZlT5QV4uvHtH/bInIi1GR0ZEQoG7At69g5jKvVzlSAbGmE4UOsp2wItXgscNqUNh4DjTiURCjsqISChwRcH4N9i78O888fVFXGI6TyiJ6Qrn3QOFq+Hksb7FjU67LyJ+UxkRCRVdBrLjrAfwfr34iKda6wZ37caIH+7s+8O9axqbdv+IS4RFpElURkSC2XcfQNJJEJva4NPBcoO7gHfYDfSqvl9KZU0tOeMG1Zt2/4hLhEWkSVRGRIJVyVZ45Vqw2eG69yHpxCNWCbkb3JlmWfDGjfT++gV+bJ9MetKZ9abdP/yok/4/izSNyohIsPLUQOd+4HBBYr9GVwvpG9y1kkbHg9hsEN8Dy2bnBFtRvecbOwKlUzcix6YyIhKsOvWGX74PVSWamKsFHXM8yJnZbIw/g7+/WMKFPzzf0BEonboRaTqVEZFg5giDqETTKUJKcYWbyhoPOeMyGh4P4nBSlTgAqD9QWEegRJpP/5wSCSaWBa9MhOVPgNdrOk1IS0+KZkBqLOlJ0Y2u49hfBK/fAGU72zCZSOjRkRGRYLLubVjzGqyfB+nnQkIv04natRMWZcP2RVBbCVc8YzqOSNBSGREJJv0urLutvderIhIACodPo6OtGjInm44iEtRURkSCiSMMhv3KdIp2qaG
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Rejection method for exp(-|x^3|) but now use importance sampling\n",
    "\n",
    "def f(x):\n",
    "    return np.exp(-np.abs(x**3))/(2*0.89298)\n",
    "\n",
    "def p(x):\n",
    "    return np.exp(-x**2/2)/(np.sqrt(2*np.pi))\n",
    "\n",
    "xmax = 2.4\n",
    "x = rng.normal(size = 10**5)    # Generate x values between +-xmax\n",
    "y = 1.5*rng.uniform(size = 10**5)\n",
    "\n",
    "# Reject\n",
    "ind = np.where(y <= f(x)/p(x))\n",
    "ind2 = np.where(y > f(x)/p(x))\n",
    "\n",
    "print(\"Acceptance fraction = %g\" % (len(x[ind])/len(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\n",
    "xx = np.linspace(-xmax,xmax,1000)\n",
    "plt.plot(xx, f(xx), 'k--')\n",
    "plt.plot(xx, p(xx), 'k:')\n",
    "plt.plot(xx, f(xx)/p(xx), 'k')\n",
    "plt.xlabel('x')\n",
    "plt.show()\n",
    "\n",
    "plot_distribution(x[ind], f)\n",
    "plot_distribution(x[ind], f, y_log=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09653426",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}