dfm/emcee_pz_dist.ipynb

## emcee_pz_dist.ipynb
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "emcee-pz-dist.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FgBLtvgYhh7x",
        "outputId": "a6a2584e-ee9e-4224-e1c7-de4b934c8f13"
      },
      "source": [
        "import sympy as sm\n",
        "\n",
        "a, z, u, zp = sm.symbols(\"a, z, u, zp\", real=True, positive=True)\n",
        "\n",
        "# The density\n",
        "g = 1 / sm.sqrt(z)\n",
        "\n",
        "# Normalize as a probability within the valid range\n",
        "norm = sm.integrate(g, (z, 1 / a, a))\n",
        "print(\"norm =\", norm)\n",
        "g /= norm\n",
        "\n",
        "# Compute the cumulative density and solve for z in: u = cdf(z)\n",
        "soln = sm.solve(sm.Eq(sm.integrate(g, (z, 1 / a, zp)), u), zp)\n",
        "print(\"z =\", soln[0])\n",
        "\n",
        "# Let's check that this is the same as the version used in the code,\n",
        "# this should give 0\n",
        "assert sm.expand(soln[0] - ((a - 1) * u + 1) ** 2 / a) == 0"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "norm = 2*sqrt(a) - 2/sqrt(a)\n",
            "z = (a*u - u + 1)**2/a\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        },
        "id": "kHiAFbt0hjIO",
        "outputId": "f80a7056-b211-470f-81dc-f1595383da80"
      },
      "source": [
        "# Let's check the result numerically\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "a0 = 2.0\n",
        "u0 = np.random.rand(100000)\n",
        "zz = ((a0 - 1.0) * u0 + 1) ** 2.0 / a0\n",
        "plt.hist(zz, 100, density=True, histtype=\"step\", color=\"k\")\n",
        "\n",
        "x0 = np.linspace(1/a0, a0, 500)\n",
        "plt.plot(x0, (1 / np.sqrt(x0)) / (2 * np.sqrt(a0) - 2 * np.sqrt(1/a0)))\n",
        "plt.xlabel(\"z\")\n",
        "plt.ylabel(\"p(z)\");"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxV9Z3/8dcnCSRkIZCFLZCEXUFUFgUR69oKzFi07rjVcaTWpZ0ZbWvb+TmOnY61rdXWwYVaq22tjFptsa4dBVdQgigIgmwJhC0h7IQly+f3x73QS8hO7r0J5/18PPLwnnu/995PMPe+z/f7Ped7zN0REZHgSoh3ASIiEl8KAhGRgFMQiIgEnIJARCTgFAQiIgGXFO8CWionJ8cLCwvjXYaISIeyYMGCLe6eW99jHS4ICgsLKSoqincZIiIdipmVNPSYhoZERAJOQSAiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgEXtSAwsyfMrMzMPmvgcTOzX5nZSjNbZGajolWLiIg0LJo9gieBiY08PgkYHP6ZBjwSxVpwd0q3VUbzLUREOqSoBYG7vwNsbaTJFOB3HjIP6GZmvaNVz8NzVjHpl+/yybrt0XoLEZEOKZ5zBHnAuojt0vB9RzCzaWZWZGZF5eXlrXqzC0fm0T21M9c8/iELSra16jVERI5FHWKy2N1nuPsYdx+Tm1vvUhlNyuvWhZnTxpGd3plrf/Mh84sb66yIiARHPINgPdAvYrtv+L6o6dOtCzOnnUbPrilc98RHfLi6IppvJyLSIcQzCGYB14aPHhoH7HD3jdF+016ZKcycNo7emSl8/bfz+WDVlmi/pYhIuxbNw0efAeYCQ82s1MxuMLObzOymcJNXgNXASuDXwM3RqqWuHl1TmDntNPp278L1v53Pm59vjtVbi4i0O1Fbhtrdr2zicQduidb7N+WUEUNZV7aNHpf+Jzf8tor7Lx/FxaP7xqscEZG46RCTxdFQUlJCTeUOVsy4FcpXcvtzn9L1lCmYWYM/uiCOiByLAhsEB6UnJ7H8199m4vBeZJ07jZ++9jm1tbW4+xE/JSUNXtdBRKTDCnwQACQnJTL9qlFceWo/ps9exfdfWEx1TW28yxIRiYkOd6nKaElMMP77ohFkpyXzP7NXsmHHPqZPHUlGSqd4lyYiElWB6hEUFhYeGu8vKCg44nEz447zh3LfxSP4YOUWLn10Lhu2741DpSIisROoICgpKTk03l9cXNxgu8tPyefJ609l/ba9XDj9fRaX7ohdkSIiMRaoIGiJCYNz+NPN4+mUmMBlj83ljSWb4l2SiEhUKAgaMaRnBi/eMp4hPdOZ9vsFZI6/gtpaj3dZIiJtSkHQhB4ZKfzvN07ja6Py6HbG1Xzz6QXs3l8d77JERNqMgqAZUjolcv+lJ8HHz/Ha4g0MuWUGhSecEu+yRETahIKgmcyM4jee5I/TxtMzfxA1593BnOVl8S5LROSoKQha6PRBOcy6dQI1O8q4/sn5PPC3L6jRvIGIdGAKglbol5XKpqe/w0Uj8/jlmyu49okPKd+1P95liYi0ioKglbxqP/dfehI/vfhE3l+2kZHff5aU/BO1MJ2IdDgKglYqKCggISGBy0/Nh7/9jCH9+9Hnqnth2ETMErRaqYh0GFprqJXqnpm8Z381P3xxMX/2C5j6ta/z80tPok+31PgUJyLSAuoRtJG05CQeuPxk7rt4BB+XbGfig++SOvT0eJclItIkBUEbMjMuPyWfV759BoU5aeRe+H1uf/ZTdu2rOqxd5OJ3Gj4SkXiz0BUjO44xY8Z4UVFRq55rZsTq962qqWXwlNvwYedTvbOcTkVPU1z01hF1xLImEQkuM1vg7mPqe0w9gijplJhA8V+n88ItExgwoD+c8690P/ufSOiUXO8S2CIi8aIgiLLRBVm88q0zuHJcIZljL+ase1/nT29/HO+yREQOURDEQEZKJ+792gj+cMNY9lfXcsmjc7nnpaXsPVAT79JERHT4aCxNGJzD6//6Je57dRlPvL+GN5dtJrnv8HiXJSIBpx5BjKUnJ/GjC0/gmRvH4Q69rrqP7Im3kdglQ0cQiUhcKAji5LSB2bz2L2dw4xn96TZqEif/4AUoPFWHlYpIzCkI4ii1cxI//IdhvHTrBPKzU2HcdVzx2FxWlu2ipKQk3uWJSEAoCNqBYX268qebxvPfF41gyYYdTHrwXTInXMW+Kk0mi0j0KQjaiYQEY+rYfN68/Sz+4cTedDv9Ss69/21eXrSxwRPOdIayiLQFBUE7k5uRzAOXnwxvPsCa5Uu45Y8f0/+GX7Jkw44j2paUlODuuPthQ0kKCBFpCR0+2k4Vz/8/amqdmfPXcufTO7jgofe4/JR87vjKELLTkxt97sGAgNASFiIijVGPoB1LTDCuGltA4qs/YvtHf+aPc1cz+q5ZPP7uavZXa/5ARNqGegQdQPGKzwFYWbaL0299gP96OZWn5haTevyZ1NY6CQna6xeR1otqj8DMJprZcjNbaWZ31vN4vpnNNrOFZrbIzCZHs56OblCPDMqe+w+e+qdTSU/uRO5Xv8NXp7/Heyu2xLs0EenAohYEZpYITAcmAcOAK81sWJ1m/w486+4jgSuAh6NVz7HkzCG5vHzbBLa89HO27ani6t98SI/L7uGz9UdOKIuINCWaPYJTgZXuvtrdDwAzgSl12jjQNXw7E9gQxXqOKQkJxp6lc3jrjjP59384ns69BvGPD73Ht55ZSFJWXrzLE5EOJJpBkAesi9guDd8X6W7gajMrBV4BbqvvhcxsmpkVmVlReXl5NGrtsJKTEvnnMwaw/tF/5uazBvK3pZvpc8PD/Nv/fkLxlj3xLk9EOoB4HzV0JfCku/cFJgO/N7MjanL3Ge4+xt3H5ObmxrzIjiC/dy7fm3Q8X/ziCmz5bF75bCPn/uJtsid9m7UVlfEuT0TasWgeNbQe6Bex3Td8X6QbgIkA7j7XzFKAHKAsinUdk4qLiw/bLtu1j0fnrObx/V/inPvncMnovtxy9iD6ZaUCoZPODp6EVlBQcMTzRSQ4otkjmA8MNrP+ZtaZ0GTwrDpt1gLnApjZ8UAKoLGfNtAjI4W7LhjGhhk3cvW4Al74eD1n/3wOdzz3KSvLdjd4VrKIBE9UL14fPhz0QSAReMLdf2xm9wBF7j4rfBTRr4F0QhPH33X3Nxp7zY5y8fpoaemF7w+22bhjL4+9vZqZ89eyv7qW3cve582HvsuJfbsdE/8uItK4xi5eH9UgiAYFQeuC4KCK3fv57fvF/OrVT0hISeeMwTk8d8+N7C3+tE2Wo9CQk0j71FgQxHuyWFqooKDg0IJyBQUFLX5+dnoyd5w/lNJHrufOScfx+cZd9LryXi56+ANeX7KJwv79j2rBOg05iXQ8CoIOpri4+NAXbXP2tiODI/KL3Q/s5aYzB/Le986m4vXpVOzZzzd+v4Cqr/yA3763mt37qgC0iqlIACgIjnGRwQEc0ZtI6ZTI7k9eZfbtZzF96ihqKndw90tLGXfvm0x7+FVKt1U2usy1QkKk49McgRwx77CgZCu/eW8Nr322CYCJJ/Ti8e9MZf/6ZUe0r+/5LZnDEJHYaGyOQKuPyqHho4O3R+V3Z9TU7qzfvpenPijmmY/W0vua++l97f3s+vhlCvoPjHPFItKW1COQJu3ZX81zRev43dwSVm/ZQ/fUTlw2ph9XjS0gPztVPQKRDkA9AjkqaclJfP30/lw3vpAPVlXwh3klPP7eGma8u5ozh+TSZeCp1NQ6ibougkiHpCCQZjMzTh+Uw+mDcti0Yx/PfLSWmfPX0uOSu/jST2czdWw+CWnd4l2miLSQhobkqFTV1NJt2BlMueN+PlhVgdfW8OXhfbhsTF/OP7EvXlMd7xJFBJ1ZfIiCIDoO/ruuKt/NmEtvZdCXr6Z8135qdm/j5smjuXRMPwb1SI93mSKBpiAIUxBER90J4gPVNby9vJyp//4QGceNp6bWGV3QncvH9GPyib1JT279iKSWsBBpHQVBmIIgOhr6ci4sLGTdlh2kDz+H7qMmQ2YvUjsnMumE3lw0Mo/TBma3eIJZRyWJtI6OGpKoamivPPJ+M6OoeCvPFa3j5UUb+dPHpfTISGbKyX24cGQew3p3PapF79RTEGk99QgkJiL/7fdV1fDWsjJeXLieOcvLqKpxhvRM58KReUw5OY+8bl2a9To6f0Gk+TQ0FKYviPhp6N++cMhwKtIKSBt+Fil9hwOwb+1i9iydQ+UXc6ndu/OwPfy2CgL1ICRoNDQk7VbJiqW4LwFgbUUlf/lkPS9+ksbq/BH0mPwtxg/M5t2nf0lil67U7tvVqqW3633f8HLZQJtch0GkI1OPQGKioX/7+u53d5Zu3MnLizby8uKNlFRUkpgQOpntH0f05ivDe9IttfNR9Qg0lCRBox6BdChmxvA+mQzvk8l3zh/Kkg07eXnxRl5etJHv/mkRP3jRmDA4h7QR57G98gDdUjvHu2SRDk09AomJlvQIGuLufLZ+J39dvIFXFm9k3da9JCYYY/tn8ddHf8znbz5LXrcu6hGI1EM9Aom7uktdt2Zy1swY0TeTEX0zuXPicfQfdSZbMwqZUzaOrC/fxOk/eYsT8rqSOf4Klm3aydCeGYeN/9edIBaREPUIJOaisTe+unw3f1u6mTeWbqaouAKzBPKzUvnKsJ58ZXgvRhd0Jykx4ah7JSIdlQ4fDdMHvn2I9rBMYnp3fv/mp7yxdBMfrKzgQE0tWWmdWfvhazzx49v50uCcw+YV9HchQaAgCNMHvn2IdhBEvuaufVXMWV7OW8vKeP79pSSmZpJgMLqgO2cf14Ozh/ZgWJ9M/V3IMU9BEKYgaB9iGQSH3Z+QSNGaLcxeVsbs5WUs2bATgOqdZVz75TGcM7QH4wdlk9pZU2dy7NFksbQrdSeOY8ZrGV3QndEF3bnj/KFs3rmP2cvK+Na9j/GXhb3544dr6ZyUwNj+WXxpcC4TBudwXK/DJ5xFjkXqEcgxp6WHqpoZ+6qqmb9mG28tK+OdFeWsLNsNQG5GMmcMymHC4NBPj4yUZi9PoWUspD1Rj0CkCclJiYe+7AE27tjLuyu28O6KLcxeXsYLC9cDcFyvDHYMOJe33/gJp/bPoksjw0haxkI6CgWBHHMih57q3t9cvTO7cNmYflw2ph+1taElL95ZUc67X2yh66gLuPaJj+iclECPy3/Ew3NWctqAbC46czQlxWta9X4i8aShIQm8lg7hJHRO4c3Fa3lvxRYeeXEOnXNDX/i1+ys558QCThuQzWkDsxneJ/PQhXf0tyfxpqEhkUbUvYBOU7xqP2cPDR16etcFwynftY95qyu47js/Zm3fS5mzvByAjJQkxvbPYtyAbDr1GEBtrZPQwiuytYbmJqSl1CMQiRD5NxL5hRqpqesjlO3cx9zVFcxbXcHcVRUUV1QCkNmlE2P7ZzF2QDanFmZxfO8MkhITovo76G9eDlKPQKSZ6h7a2pov0R5dU5hycuhqaxCaeB48fjKX/eeDzF1dwRtLNwOQ2jmRUfndOaUwi1MKuzP1/PGUrF5x6L3r25M/mr199RSkIVHtEZjZROCXQCLwuLv/pJ42lwF3Aw586u5TG3tN9QikPWnu3vdhX8LHncRDM19hfvFW5hdvY9mmnbiD11RzcmEOpxR05/F772TdwrebfYW25tSknkKwxaVHYGaJwHTgy0ApMN/MZrn70og2g4HvA6e7+zYz6xGtekTiqe48xAUn9eGCk/oAsGNvFR+v3cbXbryD5IHf5HfzSuCMb9DvjG8wMDeNT/7vBZ4tWseo/G7A3+cY2mJFVxGIYo/AzE4D7nb388Pb3wdw93sj2vwU+MLdH2/u66pHIO1Ja/ayGxqiOfj8/dU1fLZ+Bx+t2UZR8VbeWLiKxC4ZANTu282ZIwoZld+dkfndOLlft0av1tbUezXWRo4t8ZojyAPWRWyXAmPrtBkCYGbvExo+utvdX4tiTSJxF/lFW1hYeMRyG8lJiYwuyGJ0QRYwELt+LCs272Th2m3c9MP72NK/Dw+9tYLacOYMyE0je/K/8Id5JYzM7waWUO97NUQnvkk0ewSXABPd/Z/D29cAY9391og2fwWqgMuAvsA7wAh3317ntaYB0wDy8/NH13ckRzNrUo9A2lQsxt3re489+6tZVLqDj9duY+Ha7bxetIzEtO4A1B7Yy/jj8hiZ352T+nbjxL6Z9M5MOexLXnMHwROvHsF6oF/Edt/wfZFKgQ/dvQpYY2ZfAIOB+ZGN3H0GMANCQ0NRq1ikheK1gF5achKnDQyduAZgXz+FtRV7+HjtNq6/4x72DryWx99dTVVN6OOSk96ZEXmZjOjbjRPzMg+FRmtoKOnYE80eQRLwBXAuoQCYD0x19yURbSYCV7r7dWaWAywETnb3ioZeV3MEEjTN2WOv78t5X1UNyzbtYnHpdhaV7mDx+h18sXnXoSGlnl2TGZHXjecf+zl//u2vGJGXSU5Giq73fIyK2/UIzGwy8CCh8f8n3P3HZnYPUOTusyy0K3U/MBGoAX7s7jMbe00FgQRNW37xVh6oZumGnZxzyfXc8J3/ZFHpdlaW7cLC8wrVOzZzwekncUJeJsP7dGV4n0xyM5KjVo/Eji5ME6Y/WumIovHFG/k6CcmpvL90HYvXb+eH989g2ISJlITPhobQUtzDendleJ+uDOvTlYvOOoX9W0pJSDB9pjoQBUGY/milI4p2ENR3e0dlFSedNZmtnkrnHgNI73scSdn9qA6PK6V2TuT43l1Z8H9/pmLlpxwoW82B8hKoqQI0d9AeaYkJkQ4slhPSdd9rZ8R5B/uqqlmxeTfjJl3KHf/1AEs37CT9hHNh8JkAJCYYg3LTGdanKy8+8RBdBoymqryYvtkZCoV2Tj0CkQBqaS+jofa1tc66bZUs3bCTJRt2snTjTpZu2MmmnfsOPbdm707GDyvkuF4ZDO3VlaG9MhjSM52MlE5R+u2kPkfdI7DQTNJJQB9gL/CZu5e1XYkiEkst7WU01D4hwSjITqMgO41JI3ofun/rngMs37SL5Zt2cvuPfsGBofk8v6CUPQdqDrXJ69YlHA5//xmQk07npLZfkVUa12iPwMwGAt8DzgNWAOVACqEzgiuBx4Cn3L02+qWGqEcg0rEc/NzV1jrrt+9l2aZdfLF5F8s27eKlt4vwjB5YYmiftFOiMSAn/VCvYVCPDAb1SKcgO5VOUViyO0haPVlsZs8AjwDvep2GZtYTuBLY5u5PtWG9jVIQiHQsjX3uDs49rC7fwynnTeH/3f9IuCexi/Xb9x5q1ynRKMxOY3DPdAblpjOoZwaDe6TTPyeNlE6JbVZrc06W66gn1OmooTAFgUjs1b3AT3OX1d6zv5pV5btZsXk3K8P/XVW+m5KKPYdOikswyM9KZVCPUO9hcI90br9xKmuXFOEH9rb4i/polvlu7446CMxsFfAzd3804r6/uvs/tl2ZzaMgEOnYjmado8LCQkpKN9Cpex865eTTKacfnbL7kdZrIJ2y8w4tqQHQJzOFVZ98wDev+hoDctPon5PGgNx0endNYcCA/k2uyhqkIGju4aNVwNlmNhb4hrsfILS6qIhIzESulBrJzKiqrqFkayUjTv8yP3vs96ws282a5V15rmjdYZPUyUkJHDjnDr557mn0z0lj+r13kZx3HNVb18d0vaj2pLlBUOnul5vZd4F3zexSQlcUExFpkWidF5GUmMDA3HT2rpjHLWcPAuDBK0ayobaW8l37WVW+hzVb9rC6fDe/enIen2/cxetLNsO46+g9LvQaWWmdufiRDxiQk0bXsZfw2mebGJCbRn5War1zEcfKxYGaOzS00N1Hhm+fB/wPkOXuMb+imIaGRI5NdecSDmpo6Kah50a2b+rCPFU1tazdWsma8j2s3rKbNVv2HAqM8l37//4GXkvfrDQKs9N49bnfce8P/y182GwqBVlpdOmceNTfL9GehG6LOYIL3P2liO0C4Dp3v6ftymweBYFIsLTVmHxLX2fnviqKt4RC4brbvss1N98RWup7RemhK8Yd1LNrMiVLFnDNhRND4ZAdCo387FQyuzTvxLlozz0czeGjhe5e3MjjBuS5e+lRV9lMCgKRYGmrPeWj+aKt+9xte/ZTUlFJydZKSrbsoWRrJc+/NofqlO4kZWQf9tzuqZ3Iz06jMDuVgqxwSOSkkp+VRk5650NDS9G+fOjRBMFzQALwF2ABfz+hbBBwNqFrDfyHu//tqKtsJgWBiLRGWwZBY8+tPFBNt7xBzJo9j7Vb91BcUcnaikqKK/awYfveQ4e+AnTplEi/rC70657Ki394nJ/dfSf5WamH7ktLTmrLhQZbPzRkZsOAq4DTgV6Elpj4HHgFeN7d9zXy9DanIBCR1ohVEDTW5kB1LaXb/t6TWLdtL+u2VrJ2ayVLSzaTkJx6WPvstM5sWLGYSyadTb/uXZh4Qi9O7Nut2XXXqan1h4+6+1Iz+y/gZmACoaOF5hOHEBARaa14XVa0vmGeAbnpMPTwdmZGxe79h4Jh3bZK1m2t5IklH7KodDuvLt5IYXZaq4OgMc2dLH4W2Ak8Hb5rKpDp7pe1eUVNUI9ARGKtpWP2rTlprskjoiyBgsL+FK9e2arfoS1OKDvB3YdFbM82s6WtqkZEpINp6WRtW/Y+It/74Gu2teYGwcdmNs7d54WLGQu0brdcROQY11BwRAZEfY/FS3ODYDTwgZmtDW/nA8vNbDHg7n5iVKoTETmGtNczj5sbBBOjWoWIiMRNs4LA3Y8871tERI4JuuSPiEjAKQhERAJOQSAiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIREQCTkEgIhJwUQ0CM5toZsvNbKWZ3dlIu4vNzM2s3osmiIhI9EQtCMwsEZgOTAKGAVeGr39ct10G8G3gw2jVIiIiDYtmj+BUYKW7r3b3A8BMYEo97X4E3Afo+sciInEQzSDIA9ZFbJeG7zvEzEYB/dz95cZeyMymmVmRmRWVl5e3faUiIgEWt8liM0sAfgHc3lRbd5/h7mPcfUxubm70ixMRCZBoBsF6oF/Edt/wfQdlACcAc8ysGBgHzNKEsYhIbEUzCOYDg82sv5l1Bq4AZh180N13uHuOuxe6eyEwD/iquxdFsSYREakjakHg7tXArcDrwOfAs+6+xMzuMbOvRut9RUSkZZp78fpWcfdXgFfq3HdXA23PimYtIiJSP51ZLCIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJOAWBiEjAKQhERAJOQSAiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIREQCTkEgIhJwCgIRkYBTEIiIBJyCQEQk4BQEIiIBpyAQEQk4BYGISMApCEREAk5BICIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJOAWBiEjAKQhERAIuqkFgZhPNbLmZrTSzO+t5/N/MbKmZLTKzN82sIJr1iIjIkaIWBGaWCEwHJgHDgCvNbFidZguBMe5+IvA88NNo1SMiIvWLZo/gVGClu6929wPATGBKZAN3n+3uleHNeUDfKNYjIiL1iGYQ5AHrIrZLw/c15Abg1foeMLNpZlZkZkXl5eVtWKKIiLSLyWIzuxoYA/ysvsfdfYa7j3H3Mbm5ubEtTkTkGJcUxddeD/SL2O4bvu8wZnYe8EPgTHffH8V6RESkHtHsEcwHBptZfzPrDFwBzIpsYGYjgceAr7p7WRRrERGRBkQtCNy9GrgVeB34HHjW3ZeY2T1m9tVws58B6cBzZvaJmc1q4OVERCRKojk0hLu/ArxS5767Im6fF833FxGRprWLyWIREYkfBYGISMApCEREAk5BICIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJOAWBiEjAKQhERAJOQSAiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIREQCTkEgIhJwCgIRkYBTEIiIBJyCQEQk4BQEIiIBpyAQEQk4BYGISMApCEREAk5BICIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJuKgGgZlNNLPlZrbSzO6s5/FkM/vf8OMfmllhNOsREZEjRS0IzCwRmA5MAoYBV5rZsDrNbgC2ufsg4AHgvmjVIyIi9Ytmj+BUYKW7r3b3A8BMYEqdNlOAp8K3nwfONTOLYk0iIlJHUhRfOw9YF7FdCoxtqI27V5vZDiAb2BLZyMymAdPCm7vNbHlri4pCzuRQp952pr3XB+2/xvZeH7T/Gtt7fdBBajSz1tZY0NAD0QyCNuPuM4AZ8a6jPmZW5O5j4l1HQ9p7fdD+a2zv9UH7r7G91wfBrjGaQ0PrgX4R233D99XbxsySgEygIoo1iYhIHdEMgvnAYDPrb2adgSuAWXXazAKuC9++BHjL3T2KNYmISB1RGxoKj/nfCrwOJAJPuPsSM7sHKHL3WcBvgN+b2UpgK6Gw6Gja5ZBVhPZeH7T/Gtt7fdD+a2zv9UGAazTtgIuIBJvOLBYRCTgFgYhIwCkImqGppTLCbS4zs6VmtsTM/tjeajSzfDObbWYLzWyRmU2OcX1PmFmZmX3WwONmZr8K17/IzEa1s/quCte12Mw+MLOTYllfc2qMaHeKmVWb2SWxqi3ivZus0czOMrNPwp+Vt9tTfWaWaWYvmdmn4fquj3F9/cKf04PfJd+up03bf1bcXT+N/BCa6F4FDAA6A58Cw+q0GQwsBLqHt3u0wxpnAN8M3x4GFMe4xi8Bo4DPGnh8MvAqYMA44MN2Vt/4iP+/k2JdX3NqjPhbeAt4BbikvdUIdAOWAvnh7Vh/Vpqq7wfAfeHbuYQOYukcw/p6A6PCtzOAL+r5LLf5Z0U9gqY1Z6mMG4Hp7r4NwN3L2mGNDnQN384ENsSwPtz9HUIfqoZMAX7nIfOAbmbWOzbVNV2fu39w8P8vMI/QeTEx1Yx/Q4DbgD8Bsf4bBJpV41TgBXdfG24f0zqbUZ8DGeGlbtLDbatjURuAu29094/Dt3cBnxNagSFSm39WFARNq2+pjLr/Y4YAQ8zsffTDmQ0AAAMXSURBVDObZ2YTY1ZdSHNqvBu42sxKCe0t3hab0pqtOb9De3EDoT2ydsXM8oCLgEfiXUsjhgDdzWyOmS0ws2vjXVAd/wMcT2hHaTHwbXevjUch4dWYRwIf1nmozT8rHWKJiQ4gidDw0FmE9hTfMbMR7r49rlUd7krgSXe/38xOI3T+xgnx+iPvqMzsbEJBMCHetdTjQeB77l7bjtduTAJGA+cCXYC5ZjbP3b+Ib1mHnA98ApwDDAT+ZmbvuvvOWBZhZumEenb/Eov3VhA0rTlLZZQSGqerAtaY2ReEgmF+bEpsVo03ABMB3H2umaUQWmQrLkMI9WjO7xBXZnYi8Dgwyd3b41IoY4CZ4RDIASabWbW7/zm+ZR2mFKhw9z3AHjN7BziJ0Fh4e3A98BMPDcavNLM1wHHAR7EqwMw6EQqBp939hXqatPlnRUNDTWvOUhl/JtQbwMxyCHV/V7ezGtcS2gvDzI4HUoDyGNbYlFnAteEjIsYBO9x9Y7yLOsjM8oEXgGva0d7rYdy9v7sXunshoWXdb25nIQDwF2CCmSWZWSqhFYk/j3NNkSI/Jz2BocTwsxyem/gN8Lm7/6KBZm3+WVGPoAnevKUyXge+YmZLgRrgO7HcY2xmjbcDvzazfyU0Ifb18F5PTJjZM4TCMic8T/EfQKdw/Y8SmreYDKwEKgntmcVMM+q7i9AS6Q+H97irPcYrVTajxrhrqkZ3/9zMXgMWAbXA4+7e6OGwsawP+BHwpJktJnRUzvfcPZZLU58OXAMsNrNPwvf9AMiPqLHNPytaYkJEJOA0NCQiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIRI6Smd0UXl//EzNbY2az412TSEvohDKRNhJeI+Yt4Kfu/lK86xFpLvUIRNrOL4G3FALS0WitIZE2YGZfBwqAW+NcikiLaWhI5CiZ2WjgKeCMiKuYiXQYGhoSOXq3AlnA7PCE8ePxLkikJdQjEBEJOPUIREQCTkEgIhJwCgIRkYBTEIiIBJyCQEQk4BQEIiIBpyAQEQm4/w83mD8NjkhfswAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yrLOQaGXk4F6"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}
	{
	"nbformat": 4,
	"nbformat_minor": 0,
	"metadata": {
	"colab": {
	"name": "emcee-pz-dist.ipynb",
	"provenance": []
	},
	"kernelspec": {
	"name": "python3",
	"display_name": "Python 3"
	},
	"language_info": {
	"name": "python"
	}
	},
	"cells": [
	{
	"cell_type": "code",
	"metadata": {
	"colab": {
	"base_uri": "https://localhost:8080/"
	},
	"id": "FgBLtvgYhh7x",
	"outputId": "a6a2584e-ee9e-4224-e1c7-de4b934c8f13"
	},
	"source": [
	"import sympy as sm\n",
	"\n",
	"a, z, u, zp = sm.symbols(\"a, z, u, zp\", real=True, positive=True)\n",
	"\n",
	"# The density\n",
	"g = 1 / sm.sqrt(z)\n",
	"\n",
	"# Normalize as a probability within the valid range\n",
	"norm = sm.integrate(g, (z, 1 / a, a))\n",
	"print(\"norm =\", norm)\n",
	"g /= norm\n",
	"\n",
	"# Compute the cumulative density and solve for z in: u = cdf(z)\n",
	"soln = sm.solve(sm.Eq(sm.integrate(g, (z, 1 / a, zp)), u), zp)\n",
	"print(\"z =\", soln[0])\n",
	"\n",
	"# Let's check that this is the same as the version used in the code,\n",
	"# this should give 0\n",
	"assert sm.expand(soln[0] - ((a - 1) * u + 1) ** 2 / a) == 0"
	],
	"execution_count": 1,
	"outputs": [
	{
	"output_type": "stream",
	"text": [
	"norm = 2*sqrt(a) - 2/sqrt(a)\n",
	"z = (au - u + 1)*2/a\n"
	],
	"name": "stdout"
	}
	]
	},
	{
	"cell_type": "code",
	"metadata": {
	"colab": {
	"base_uri": "https://localhost:8080/",
	"height": 279
	},
	"id": "kHiAFbt0hjIO",
	"outputId": "f80a7056-b211-470f-81dc-f1595383da80"
	},
	"source": [
	"# Let's check the result numerically\n",
	"import numpy as np\n",
	"import matplotlib.pyplot as plt\n",
	"\n",
	"a0 = 2.0\n",
	"u0 = np.random.rand(100000)\n",
	"zz = ((a0 - 1.0) * u0 + 1) ** 2.0 / a0\n",
	"plt.hist(zz, 100, density=True, histtype=\"step\", color=\"k\")\n",
	"\n",
	"x0 = np.linspace(1/a0, a0, 500)\n",
	"plt.plot(x0, (1 / np.sqrt(x0)) / (2 * np.sqrt(a0) - 2 * np.sqrt(1/a0)))\n",
	"plt.xlabel(\"z\")\n",
	"plt.ylabel(\"p(z)\");"
	],
	"execution_count": 6,
	"outputs": [
	{
	"output_type": "display_data",
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxV9Z3/8dcnCSRkIZCFLZCEXUFUFgUR69oKzFi07rjVcaTWpZ0ZbWvb+TmOnY61rdXWwYVaq22tjFptsa4dBVdQgigIgmwJhC0h7IQly+f3x73QS8hO7r0J5/18PPLwnnu/995PMPe+z/f7Ped7zN0REZHgSoh3ASIiEl8KAhGRgFMQiIgEnIJARCTgFAQiIgGXFO8CWionJ8cLCwvjXYaISIeyYMGCLe6eW99jHS4ICgsLKSoqincZIiIdipmVNPSYhoZERAJOQSAiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgEXtSAwsyfMrMzMPmvgcTOzX5nZSjNbZGajolWLiIg0LJo9gieBiY08PgkYHP6ZBjwSxVpwd0q3VUbzLUREOqSoBYG7vwNsbaTJFOB3HjIP6GZmvaNVz8NzVjHpl+/yybrt0XoLEZEOKZ5zBHnAuojt0vB9RzCzaWZWZGZF5eXlrXqzC0fm0T21M9c8/iELSra16jVERI5FHWKy2N1nuPsYdx+Tm1vvUhlNyuvWhZnTxpGd3plrf/Mh84sb66yIiARHPINgPdAvYrtv+L6o6dOtCzOnnUbPrilc98RHfLi6IppvJyLSIcQzCGYB14aPHhoH7HD3jdF+016ZKcycNo7emSl8/bfz+WDVlmi/pYhIuxbNw0efAeYCQ82s1MxuMLObzOymcJNXgNXASuDXwM3RqqWuHl1TmDntNPp278L1v53Pm59vjtVbi4i0O1Fbhtrdr2zicQduidb7N+WUEUNZV7aNHpf+Jzf8tor7Lx/FxaP7xqscEZG46RCTxdFQUlJCTeUOVsy4FcpXcvtzn9L1lCmYWYM/uiCOiByLAhsEB6UnJ7H8199m4vBeZJ07jZ++9jm1tbW4+xE/JSUNXtdBRKTDCnwQACQnJTL9qlFceWo/ps9exfdfWEx1TW28yxIRiYkOd6nKaElMMP77ohFkpyXzP7NXsmHHPqZPHUlGSqd4lyYiElWB6hEUFhYeGu8vKCg44nEz447zh3LfxSP4YOUWLn10Lhu2741DpSIisROoICgpKTk03l9cXNxgu8tPyefJ609l/ba9XDj9fRaX7ohdkSIiMRaoIGiJCYNz+NPN4+mUmMBlj83ljSWb4l2SiEhUKAgaMaRnBi/eMp4hPdOZ9vsFZI6/gtpaj3dZIiJtSkHQhB4ZKfzvN07ja6Py6HbG1Xzz6QXs3l8d77JERNqMgqAZUjolcv+lJ8HHz/Ha4g0MuWUGhSecEu+yRETahIKgmcyM4jee5I/TxtMzfxA1593BnOVl8S5LROSoKQha6PRBOcy6dQI1O8q4/sn5PPC3L6jRvIGIdGAKglbol5XKpqe/w0Uj8/jlmyu49okPKd+1P95liYi0ioKglbxqP/dfehI/vfhE3l+2kZHff5aU/BO1MJ2IdDgKglYqKCggISGBy0/Nh7/9jCH9+9Hnqnth2ETMErRaqYh0GFprqJXqnpm8Z381P3xxMX/2C5j6ta/z80tPok+31PgUJyLSAuoRtJG05CQeuPxk7rt4BB+XbGfig++SOvT0eJclItIkBUEbMjMuPyWfV759BoU5aeRe+H1uf/ZTdu2rOqxd5OJ3Gj4SkXiz0BUjO44xY8Z4UVFRq55rZsTq962qqWXwlNvwYedTvbOcTkVPU1z01hF1xLImEQkuM1vg7mPqe0w9gijplJhA8V+n88ItExgwoD+c8690P/ufSOiUXO8S2CIi8aIgiLLRBVm88q0zuHJcIZljL+ase1/nT29/HO+yREQOURDEQEZKJ+792gj+cMNY9lfXcsmjc7nnpaXsPVAT79JERHT4aCxNGJzD6//6Je57dRlPvL+GN5dtJrnv8HiXJSIBpx5BjKUnJ/GjC0/gmRvH4Q69rrqP7Im3kdglQ0cQiUhcKAji5LSB2bz2L2dw4xn96TZqEif/4AUoPFWHlYpIzCkI4ii1cxI//IdhvHTrBPKzU2HcdVzx2FxWlu2ipKQk3uWJSEAoCNqBYX268qebxvPfF41gyYYdTHrwXTInXMW+Kk0mi0j0KQjaiYQEY+rYfN68/Sz+4cTedDv9Ss69/21eXrSxwRPOdIayiLQFBUE7k5uRzAOXnwxvPsCa5Uu45Y8f0/+GX7Jkw44j2paUlODuuPthQ0kKCBFpCR0+2k4Vz/8/amqdmfPXcufTO7jgofe4/JR87vjKELLTkxt97sGAgNASFiIijVGPoB1LTDCuGltA4qs/YvtHf+aPc1cz+q5ZPP7uavZXa/5ARNqGegQdQPGKzwFYWbaL0299gP96OZWn5haTevyZ1NY6CQna6xeR1otqj8DMJprZcjNbaWZ31vN4vpnNNrOFZrbIzCZHs56OblCPDMqe+w+e+qdTSU/uRO5Xv8NXp7/Heyu2xLs0EenAohYEZpYITAcmAcOAK81sWJ1m/w486+4jgSuAh6NVz7HkzCG5vHzbBLa89HO27ani6t98SI/L7uGz9UdOKIuINCWaPYJTgZXuvtrdDwAzgSl12jjQNXw7E9gQxXqOKQkJxp6lc3jrjjP59384ns69BvGPD73Ht55ZSFJWXrzLE5EOJJpBkAesi9guDd8X6W7gajMrBV4BbqvvhcxsmpkVmVlReXl5NGrtsJKTEvnnMwaw/tF/5uazBvK3pZvpc8PD/Nv/fkLxlj3xLk9EOoB4HzV0JfCku/cFJgO/N7MjanL3Ge4+xt3H5ObmxrzIjiC/dy7fm3Q8X/ziCmz5bF75bCPn/uJtsid9m7UVlfEuT0TasWgeNbQe6Bex3Td8X6QbgIkA7j7XzFKAHKAsinUdk4qLiw/bLtu1j0fnrObx/V/inPvncMnovtxy9iD6ZaUCoZPODp6EVlBQcMTzRSQ4otkjmA8MNrP+ZtaZ0GTwrDpt1gLnApjZ8UAKoLGfNtAjI4W7LhjGhhk3cvW4Al74eD1n/3wOdzz3KSvLdjd4VrKIBE9UL14fPhz0QSAReMLdf2xm9wBF7j4rfBTRr4F0QhPH33X3Nxp7zY5y8fpoaemF7w+22bhjL4+9vZqZ89eyv7qW3cve582HvsuJfbsdE/8uItK4xi5eH9UgiAYFQeuC4KCK3fv57fvF/OrVT0hISeeMwTk8d8+N7C3+tE2Wo9CQk0j71FgQxHuyWFqooKDg0IJyBQUFLX5+dnoyd5w/lNJHrufOScfx+cZd9LryXi56+ANeX7KJwv79j2rBOg05iXQ8CoIOpri4+NAXbXP2tiODI/KL3Q/s5aYzB/Le986m4vXpVOzZzzd+v4Cqr/yA3763mt37qgC0iqlIACgIjnGRwQEc0ZtI6ZTI7k9eZfbtZzF96ihqKndw90tLGXfvm0x7+FVKt1U2usy1QkKk49McgRwx77CgZCu/eW8Nr322CYCJJ/Ti8e9MZf/6ZUe0r+/5LZnDEJHYaGyOQKuPyqHho4O3R+V3Z9TU7qzfvpenPijmmY/W0vua++l97f3s+vhlCvoPjHPFItKW1COQJu3ZX81zRev43dwSVm/ZQ/fUTlw2ph9XjS0gPztVPQKRDkA9AjkqaclJfP30/lw3vpAPVlXwh3klPP7eGma8u5ozh+TSZeCp1NQ6ibougkiHpCCQZjMzTh+Uw+mDcti0Yx/PfLSWmfPX0uOSu/jST2czdWw+CWnd4l2miLSQhobkqFTV1NJt2BlMueN+PlhVgdfW8OXhfbhsTF/OP7EvXlMd7xJFBJ1ZfIiCIDoO/ruuKt/NmEtvZdCXr6Z8135qdm/j5smjuXRMPwb1SI93mSKBpiAIUxBER90J4gPVNby9vJyp//4QGceNp6bWGV3QncvH9GPyib1JT279iKSWsBBpHQVBmIIgOhr6ci4sLGTdlh2kDz+H7qMmQ2YvUjsnMumE3lw0Mo/TBma3eIJZRyWJtI6OGpKoamivPPJ+M6OoeCvPFa3j5UUb+dPHpfTISGbKyX24cGQew3p3PapF79RTEGk99QgkJiL/7fdV1fDWsjJeXLieOcvLqKpxhvRM58KReUw5OY+8bl2a9To6f0Gk+TQ0FKYviPhp6N++cMhwKtIKSBt+Fil9hwOwb+1i9iydQ+UXc6ndu/OwPfy2CgL1ICRoNDQk7VbJiqW4LwFgbUUlf/lkPS9+ksbq/BH0mPwtxg/M5t2nf0lil67U7tvVqqW3633f8HLZQJtch0GkI1OPQGKioX/7+u53d5Zu3MnLizby8uKNlFRUkpgQOpntH0f05ivDe9IttfNR9Qg0lCRBox6BdChmxvA+mQzvk8l3zh/Kkg07eXnxRl5etJHv/mkRP3jRmDA4h7QR57G98gDdUjvHu2SRDk09AomJlvQIGuLufLZ+J39dvIFXFm9k3da9JCYYY/tn8ddHf8znbz5LXrcu6hGI1EM9Aom7uktdt2Zy1swY0TeTEX0zuXPicfQfdSZbMwqZUzaOrC/fxOk/eYsT8rqSOf4Klm3aydCeGYeN/9edIBaREPUIJOaisTe+unw3f1u6mTeWbqaouAKzBPKzUvnKsJ58ZXgvRhd0Jykx4ah7JSIdlQ4fDdMHvn2I9rBMYnp3fv/mp7yxdBMfrKzgQE0tWWmdWfvhazzx49v50uCcw+YV9HchQaAgCNMHvn2IdhBEvuaufVXMWV7OW8vKeP79pSSmZpJgMLqgO2cf14Ozh/ZgWJ9M/V3IMU9BEKYgaB9iGQSH3Z+QSNGaLcxeVsbs5WUs2bATgOqdZVz75TGcM7QH4wdlk9pZU2dy7NFksbQrdSeOY8ZrGV3QndEF3bnj/KFs3rmP2cvK+Na9j/GXhb3544dr6ZyUwNj+WXxpcC4TBudwXK/DJ5xFjkXqEcgxp6WHqpoZ+6qqmb9mG28tK+OdFeWsLNsNQG5GMmcMymHC4NBPj4yUZi9PoWUspD1Rj0CkCclJiYe+7AE27tjLuyu28O6KLcxeXsYLC9cDcFyvDHYMOJe33/gJp/bPoksjw0haxkI6CgWBHHMih57q3t9cvTO7cNmYflw2ph+1taElL95ZUc67X2yh66gLuPaJj+iclECPy3/Ew3NWctqAbC46czQlxWta9X4i8aShIQm8lg7hJHRO4c3Fa3lvxRYeeXEOnXNDX/i1+ys558QCThuQzWkDsxneJ/PQhXf0tyfxpqEhkUbUvYBOU7xqP2cPDR16etcFwynftY95qyu47js/Zm3fS5mzvByAjJQkxvbPYtyAbDr1GEBtrZPQwiuytYbmJqSl1CMQiRD5NxL5hRqpqesjlO3cx9zVFcxbXcHcVRUUV1QCkNmlE2P7ZzF2QDanFmZxfO8MkhITovo76G9eDlKPQKSZ6h7a2pov0R5dU5hycuhqaxCaeB48fjKX/eeDzF1dwRtLNwOQ2jmRUfndOaUwi1MKuzP1/PGUrF5x6L3r25M/mr199RSkIVHtEZjZROCXQCLwuLv/pJ42lwF3Aw586u5TG3tN9QikPWnu3vdhX8LHncRDM19hfvFW5hdvY9mmnbiD11RzcmEOpxR05/F772TdwrebfYW25tSknkKwxaVHYGaJwHTgy0ApMN/MZrn70og2g4HvA6e7+zYz6xGtekTiqe48xAUn9eGCk/oAsGNvFR+v3cbXbryD5IHf5HfzSuCMb9DvjG8wMDeNT/7vBZ4tWseo/G7A3+cY2mJFVxGIYo/AzE4D7nb388Pb3wdw93sj2vwU+MLdH2/u66pHIO1Ja/ayGxqiOfj8/dU1fLZ+Bx+t2UZR8VbeWLiKxC4ZANTu282ZIwoZld+dkfndOLlft0av1tbUezXWRo4t8ZojyAPWRWyXAmPrtBkCYGbvExo+utvdX4tiTSJxF/lFW1hYeMRyG8lJiYwuyGJ0QRYwELt+LCs272Th2m3c9MP72NK/Dw+9tYLacOYMyE0je/K/8Id5JYzM7waWUO97NUQnvkk0ewSXABPd/Z/D29cAY9391og2fwWqgMuAvsA7wAh3317ntaYB0wDy8/NH13ckRzNrUo9A2lQsxt3re489+6tZVLqDj9duY+Ha7bxetIzEtO4A1B7Yy/jj8hiZ352T+nbjxL6Z9M5MOexLXnMHwROvHsF6oF/Edt/wfZFKgQ/dvQpYY2ZfAIOB+ZGN3H0GMANCQ0NRq1ikheK1gF5achKnDQyduAZgXz+FtRV7+HjtNq6/4x72DryWx99dTVVN6OOSk96ZEXmZjOjbjRPzMg+FRmtoKOnYE80eQRLwBXAuoQCYD0x19yURbSYCV7r7dWaWAywETnb3ioZeV3MEEjTN2WOv78t5X1UNyzbtYnHpdhaV7mDx+h18sXnXoSGlnl2TGZHXjecf+zl//u2vGJGXSU5Giq73fIyK2/UIzGwy8CCh8f8n3P3HZnYPUOTusyy0K3U/MBGoAX7s7jMbe00FgQRNW37xVh6oZumGnZxzyfXc8J3/ZFHpdlaW7cLC8wrVOzZzwekncUJeJsP7dGV4n0xyM5KjVo/Eji5ME6Y/WumIovHFG/k6CcmpvL90HYvXb+eH989g2ISJlITPhobQUtzDendleJ+uDOvTlYvOOoX9W0pJSDB9pjoQBUGY/milI4p2ENR3e0dlFSedNZmtnkrnHgNI73scSdn9qA6PK6V2TuT43l1Z8H9/pmLlpxwoW82B8hKoqQI0d9AeaYkJkQ4slhPSdd9rZ8R5B/uqqlmxeTfjJl3KHf/1AEs37CT9hHNh8JkAJCYYg3LTGdanKy8+8RBdBoymqryYvtkZCoV2Tj0CkQBqaS+jofa1tc66bZUs3bCTJRt2snTjTpZu2MmmnfsOPbdm707GDyvkuF4ZDO3VlaG9MhjSM52MlE5R+u2kPkfdI7DQTNJJQB9gL/CZu5e1XYkiEkst7WU01D4hwSjITqMgO41JI3ofun/rngMs37SL5Zt2cvuPfsGBofk8v6CUPQdqDrXJ69YlHA5//xmQk07npLZfkVUa12iPwMwGAt8DzgNWAOVACqEzgiuBx4Cn3L02+qWGqEcg0rEc/NzV1jrrt+9l2aZdfLF5F8s27eKlt4vwjB5YYmiftFOiMSAn/VCvYVCPDAb1SKcgO5VOUViyO0haPVlsZs8AjwDvep2GZtYTuBLY5u5PtWG9jVIQiHQsjX3uDs49rC7fwynnTeH/3f9IuCexi/Xb9x5q1ynRKMxOY3DPdAblpjOoZwaDe6TTPyeNlE6JbVZrc06W66gn1OmooTAFgUjs1b3AT3OX1d6zv5pV5btZsXk3K8P/XVW+m5KKPYdOikswyM9KZVCPUO9hcI90br9xKmuXFOEH9rb4i/polvlu7446CMxsFfAzd3804r6/uvs/tl2ZzaMgEOnYjmado8LCQkpKN9Cpex865eTTKacfnbL7kdZrIJ2y8w4tqQHQJzOFVZ98wDev+hoDctPon5PGgNx0endNYcCA/k2uyhqkIGju4aNVwNlmNhb4hrsfILS6qIhIzESulBrJzKiqrqFkayUjTv8yP3vs96ws282a5V15rmjdYZPUyUkJHDjnDr557mn0z0lj+r13kZx3HNVb18d0vaj2pLlBUOnul5vZd4F3zexSQlcUExFpkWidF5GUmMDA3HT2rpjHLWcPAuDBK0ayobaW8l37WVW+hzVb9rC6fDe/enIen2/cxetLNsO46+g9LvQaWWmdufiRDxiQk0bXsZfw2mebGJCbRn5War1zEcfKxYGaOzS00N1Hhm+fB/wPkOXuMb+imIaGRI5NdecSDmpo6Kah50a2b+rCPFU1tazdWsma8j2s3rKbNVv2HAqM8l37//4GXkvfrDQKs9N49bnfce8P/y182GwqBVlpdOmceNTfL9GehG6LOYIL3P2liO0C4Dp3v6ftymweBYFIsLTVmHxLX2fnviqKt4RC4brbvss1N98RWup7RemhK8Yd1LNrMiVLFnDNhRND4ZAdCo387FQyuzTvxLlozz0czeGjhe5e3MjjBuS5e+lRV9lMCgKRYGmrPeWj+aKt+9xte/ZTUlFJydZKSrbsoWRrJc+/NofqlO4kZWQf9tzuqZ3Iz06jMDuVgqxwSOSkkp+VRk5650NDS9G+fOjRBMFzQALwF2ABfz+hbBBwNqFrDfyHu//tqKtsJgWBiLRGWwZBY8+tPFBNt7xBzJo9j7Vb91BcUcnaikqKK/awYfveQ4e+AnTplEi/rC70657Ki394nJ/dfSf5WamH7ktLTmrLhQZbPzRkZsOAq4DTgV6Elpj4HHgFeN7d9zXy9DanIBCR1ohVEDTW5kB1LaXb/t6TWLdtL+u2VrJ2ayVLSzaTkJx6WPvstM5sWLGYSyadTb/uXZh4Qi9O7Nut2XXXqan1h4+6+1Iz+y/gZmACoaOF5hOHEBARaa14XVa0vmGeAbnpMPTwdmZGxe79h4Jh3bZK1m2t5IklH7KodDuvLt5IYXZaq4OgMc2dLH4W2Ak8Hb5rKpDp7pe1eUVNUI9ARGKtpWP2rTlprskjoiyBgsL+FK9e2arfoS1OKDvB3YdFbM82s6WtqkZEpINp6WRtW/Y+It/74Gu2teYGwcdmNs7d54WLGQu0brdcROQY11BwRAZEfY/FS3ODYDTwgZmtDW/nA8vNbDHg7n5iVKoTETmGtNczj5sbBBOjWoWIiMRNs4LA3Y8871tERI4JuuSPiEjAKQhERAJOQSAiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIREQCTkEgIhJwUQ0CM5toZsvNbKWZ3dlIu4vNzM2s3osmiIhI9EQtCMwsEZgOTAKGAVeGr39ct10G8G3gw2jVIiIiDYtmj+BUYKW7r3b3A8BMYEo97X4E3Afo+sciInEQzSDIA9ZFbJeG7zvEzEYB/dz95cZeyMymmVmRmRWVl5e3faUiIgEWt8liM0sAfgHc3lRbd5/h7mPcfUxubm70ixMRCZBoBsF6oF/Edt/wfQdlACcAc8ysGBgHzNKEsYhIbEUzCOYDg82sv5l1Bq4AZh180N13uHuOuxe6eyEwD/iquxdFsSYREakjakHg7tXArcDrwOfAs+6+xMzuMbOvRut9RUSkZZp78fpWcfdXgFfq3HdXA23PimYtIiJSP51ZLCIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJOAWBiEjAKQhERAJOQSAiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIREQCTkEgIhJwCgIRkYBTEIiIBJyCQEQk4BQEIiIBpyAQEQk4BYGISMApCEREAk5BICIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJOAWBiEjAKQhERAIuqkFgZhPNbLmZrTSzO+t5/N/MbKmZLTKzN82sIJr1iIjIkaIWBGaWCEwHJgHDgCvNbFidZguBMe5+IvA88NNo1SMiIvWLZo/gVGClu6929wPATGBKZAN3n+3uleHNeUDfKNYjIiL1iGYQ5AHrIrZLw/c15Abg1foeMLNpZlZkZkXl5eVtWKKIiLSLyWIzuxoYA/ysvsfdfYa7j3H3Mbm5ubEtTkTkGJcUxddeD/SL2O4bvu8wZnYe8EPgTHffH8V6RESkHtHsEcwHBptZfzPrDFwBzIpsYGYjgceAr7p7WRRrERGRBkQtCNy9GrgVeB34HHjW3ZeY2T1m9tVws58B6cBzZvaJmc1q4OVERCRKojk0hLu/ArxS5767Im6fF833FxGRprWLyWIREYkfBYGISMApCEREAk5BICIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJOAWBiEjAKQhERAJOQSAiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIREQCTkEgIhJwCgIRkYBTEIiIBJyCQEQk4BQEIiIBpyAQEQk4BYGISMApCEREAk5BICIScAoCEZGAUxCIiAScgkBEJOAUBCIiAacgEBEJuKgGgZlNNLPlZrbSzO6s5/FkM/vf8OMfmllhNOsREZEjRS0IzCwRmA5MAoYBV5rZsDrNbgC2ufsg4AHgvmjVIyIi9Ytmj+BUYKW7r3b3A8BMYEqdNlOAp8K3nwfONTOLYk0iIlJHUhRfOw9YF7FdCoxtqI27V5vZDiAb2BLZyMymAdPCm7vNbHlri4pCzuRQp952pr3XB+2/xvZeH7T/Gtt7fdBBajSz1tZY0NAD0QyCNuPuM4AZ8a6jPmZW5O5j4l1HQ9p7fdD+a2zv9UH7r7G91wfBrjGaQ0PrgX4R233D99XbxsySgEygIoo1iYhIHdEMgvnAYDPrb2adgSuAWXXazAKuC9++BHjL3T2KNYmISB1RGxoKj/nfCrwOJAJPuPsSM7sHKHL3WcBvgN+b2UpgK6Gw6Gja5ZBVhPZeH7T/Gtt7fdD+a2zv9UGAazTtgIuIBJvOLBYRCTgFgYhIwCkImqGppTLCbS4zs6VmtsTM/tjeajSzfDObbWYLzWyRmU2OcX1PmFmZmX3WwONmZr8K17/IzEa1s/quCte12Mw+MLOTYllfc2qMaHeKmVWb2SWxqi3ivZus0czOMrNPwp+Vt9tTfWaWaWYvmdmn4fquj3F9/cKf04PfJd+up03bf1bcXT+N/BCa6F4FDAA6A58Cw+q0GQwsBLqHt3u0wxpnAN8M3x4GFMe4xi8Bo4DPGnh8MvAqYMA44MN2Vt/4iP+/k2JdX3NqjPhbeAt4BbikvdUIdAOWAvnh7Vh/Vpqq7wfAfeHbuYQOYukcw/p6A6PCtzOAL+r5LLf5Z0U9gqY1Z6mMG4Hp7r4NwN3L2mGNDnQN384ENsSwPtz9HUIfqoZMAX7nIfOAbmbWOzbVNV2fu39w8P8vMI/QeTEx1Yx/Q4DbgD8Bsf4bBJpV41TgBXdfG24f0zqbUZ8DGeGlbtLDbatjURuAu29094/Dt3cBnxNagSFSm39WFARNq2+pjLr/Y4YAQ8zsffTDmQ0AAAMXSURBVDObZ2YTY1ZdSHNqvBu42sxKCe0t3hab0pqtOb9De3EDoT2ydsXM8oCLgEfiXUsjhgDdzWyOmS0ws2vjXVAd/wMcT2hHaTHwbXevjUch4dWYRwIf1nmozT8rHWKJiQ4gidDw0FmE9hTfMbMR7r49rlUd7krgSXe/38xOI3T+xgnx+iPvqMzsbEJBMCHetdTjQeB77l7bjtduTAJGA+cCXYC5ZjbP3b+Ib1mHnA98ApwDDAT+ZmbvuvvOWBZhZumEenb/Eov3VhA0rTlLZZQSGqerAtaY2ReEgmF+bEpsVo03ABMB3H2umaUQWmQrLkMI9WjO7xBXZnYi8Dgwyd3b41IoY4CZ4RDIASabWbW7/zm+ZR2mFKhw9z3AHjN7BziJ0Fh4e3A98BMPDcavNLM1wHHAR7EqwMw6EQqBp939hXqatPlnRUNDTWvOUhl/JtQbwMxyCHV/V7ezGtcS2gvDzI4HUoDyGNbYlFnAteEjIsYBO9x9Y7yLOsjM8oEXgGva0d7rYdy9v7sXunshoWXdb25nIQDwF2CCmSWZWSqhFYk/j3NNkSI/Jz2BocTwsxyem/gN8Lm7/6KBZm3+WVGPoAnevKUyXge+YmZLgRrgO7HcY2xmjbcDvzazfyU0Ifb18F5PTJjZM4TCMic8T/EfQKdw/Y8SmreYDKwEKgntmcVMM+q7i9AS6Q+H97irPcYrVTajxrhrqkZ3/9zMXgMWAbXA4+7e6OGwsawP+BHwpJktJnRUzvfcPZZLU58OXAMsNrNPwvf9AMiPqLHNPytaYkJEJOA0NCQiEnAKAhGRgFMQiIgEnIJARCTgFAQiIgGnIBARCTgFgYhIwCkIRI6Smd0UXl//EzNbY2az412TSEvohDKRNhJeI+Yt4Kfu/lK86xFpLvUIRNrOL4G3FALS0WitIZE2YGZfBwqAW+NcikiLaWhI5CiZ2WjgKeCMiKuYiXQYGhoSOXq3AlnA7PCE8ePxLkikJdQjEBEJOPUIREQCTkEgIhJwCgIRkYBTEIiIBJyCQEQk4BQEIiIBpyAQEQm4/w83mD8NjkhfswAAAABJRU5ErkJggg==\n",
	"text/plain": [
	"<Figure size 432x288 with 1 Axes>"
	]
	},
	"metadata": {
	"tags": [],
	"needs_background": "light"
	}
	}
	]
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "yrLOQaGXk4F6"
	},
	"source": [
	""
	],
	"execution_count": null,
	"outputs": []
	}
	]
	}