alinaselega/lognormal_cdf.ipynb

## lognormal_cdf.ipynb
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "lognormal_cdf.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1YNdXJuElQpk",
        "colab_type": "text"
      },
      "source": [
        "This notebook was intended to simulate some lognormal data and evaluate likelihoods. It became a numerical evaluation of whether lognormal PDF on exponents on MC samples is equal to the normal PDF on those samples."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CQYyWEXGlQ3-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import numpy.random as rn\n",
        "import scipy.special as sc\n",
        "import matplotlib.pyplot as plt\n",
        "plt.style.use('ggplot')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GDDKyRHyllRI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Set parameters\n",
        "sigma = 0.25\n",
        "mean = 0"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zRI84YoxxF4g",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Lognormal PDF\n",
        "def lognormal_pdf(x, mean, sigma):\n",
        "    c = 1/(x*sigma*np.sqrt(2*np.pi))\n",
        "    return c * np.exp(-(np.log(x) - mean)**2 / (2*sigma**2)) "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1ePbGvfKym84",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Normal PDF\n",
        "def normal_pdf(x, mean, sigma):\n",
        "    c = 1/(np.sqrt(2*np.pi)*sigma)\n",
        "    return c * np.exp(-(x-mean)**2 / (2*sigma**2))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zaXOzINGlp31",
        "colab_type": "text"
      },
      "source": [
        "Check that log-taking was not causing problems."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q_Fahd3v5seP",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def lognormal_logpdf(x, mean, sigma):\n",
        "    return -np.log(x*sigma*np.sqrt(2*np.pi)) - (np.log(x) - mean)**2 / (2*sigma**2)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8Kki8i5Kqof9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def lognormal_cdf(x, mean, sigma):\n",
        "    return 0.5 + 0.5 * sc.erf((np.log(x) - mean) / np.sqrt(2 * sigma**2))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gEyErqcSrhVM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def normal_cdf(x, mean, sigma):\n",
        "    return 0.5 * (1 + sc.erf((x - mean) / np.sqrt(2 * sigma**2)))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n70gTbY3lt_u",
        "colab_type": "text"
      },
      "source": [
        "Plot lognormal and normal PDFs to validate my functions."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "14_KoUjOmflH",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        },
        "outputId": "0e4a802e-9658-40f9-f3f7-102cb719d72d"
      },
      "source": [
        "x = np.linspace(0.1, 3, 100)\n",
        "plt.axes(xlabel=\"x\", ylabel=\"Lognorm PDF\")\n",
        "plt.plot(x, lognormal_pdf(x, mean, sigma),'r-', lw=5, \n",
        "         alpha=0.6)\n",
        "plt.show()"
      ],
      "execution_count": 28,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEJCAYAAABlmAtYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3hTdZ4/8Pc5J7e2KaVJoKXQAi3l\nUhChlFsZWGgLUx1XWS/w7DgrgrrDOK6LPqvCPjyjo4tWEdZBQAetHR3HHXYcZ3Sd0Z92EBEqDFAK\ncm3LtUChtClt0zbN5ZzfH4HYk3PSpG2Sk5x+Xs/DQ8/3nCSfLyn55HyvjCAIAgghhJAesEoHQAgh\nJPpRsiCEEBIQJQtCCCEBUbIghBASECULQgghAVGyIIQQEpBG6QDC6fLly96fLRYLGhsbFYwm9NRW\nJ7XVB1BfndRWH0B9depPfdLS0vyeozsLQgghAVGyIIQQEhAlC0IIIQFRsiCEEBIQJQtCCCEBqXo0\nFIk+TGcn2KYmCAYDeItF6XAIIUGiZEHCTnPiBPQVFeDq68G0tnrLnRMnonPJEghxcQpGRwgJBjVD\nkbDSVlUhoawMmlOnRIkCALTHjsH4+utgr1xRKDpCSLAoWZCw4erqEP+HP/R4DdvYCOOWLdAcOxah\nqAghfUHJgoQF09KChHffBZzOwNd2dSHh/feBs2cjEBkhpC8oWZDQczqR8N57kmYnAADLAhqZrjK3\nG+zvfw/wfPjjI4T0GiULEnL6nTvB1dVJyt0jRqDll79E2xNPgB86VHKeqauD7ttvIxEiIaSXKFmQ\n0HI6oZf5wBcSE9H+wAOAXg9+yBC0PfYY3MOGSa4z/L//B6atLRKREkJ6gZIFCSltVRUYm01cyLJo\nf+ABCIMHf19mMKDz7rslj2fsdhj+8pcwR0kI6S1KFiR0BAH63bslxY4pU+AeOVJS7h45Eo7p0yXl\nuspKcNTZTUhUoWRBQoY7cwZcfb2k3DFnjt/H2G+/HUJ8vKRc//XXIY2NENI/lCxIyOj37JGUuUeO\nhDs93e9jhIQE2G+7TVKuPXECrNUa0vgIIX1HyYKEBGO1Qiszsa7rBz8I+FjH9Ongu/dnAIAg0Mgo\nQqIIJQsSEvpvvwUEQVTGJyXBOWlS4AezLByzZ0uKdX//O+BwhCpEQkg/ULIg/ScI0B45Iil2zJ4N\ncFxQT+GYMUMyWY/p7ISuqiokIRJC+oeSBek39upVsM3NPoWsJwEESUhIgCA3MmrPHskdCyEk8ihZ\nkH7TnjwpKXONGgXBaOzV8/Dz5knKuPp6GkZLSBSIyH4WW7duRWVlJZKSkrBhwwbJ+WPHjuGVV17B\n0BtLQMycORP33nsvAKCqqgplZWXgeR6FhYVYvHhxJEImvaCRSxbjx/f+iUaOhHvkSHDnz4uKdfv3\nozMzs6/hEUJCICLJYv78+SguLsaWLVv8XjNhwgSsXr1aVMbzPEpLS7F27VqYzWasWbMGeXl5GDFi\nRLhDJkFiOjqgOXdOUu6cMKFPz9eVn494n2ShPXYMnU4noNX26TkJIf0XkWaonJwcGHvZJAEAtbW1\nSE1NRUpKCjQaDfLz87F///4wREj6SlNdLVkpljeZZBcKDIZz0iQIBoOojLHboTl1qs8xEkL6L2r6\nLKqrq/HUU0/hxRdfRN2NFUutVivMZrP3GrPZDCtN1Ioq2hMnJGXOCRMAhunjE2rhnDhRUqw7fLhv\nz0cICYmo2IN79OjR2Lp1KwwGAyorK7F+/Xps2rSp189TXl6O8vJyAEBJSQksFov3nEajER2rgeJ1\n4nlwFy4APncCutmzYexDXDfrw8yfD9Zngp/h7FnEJyYCen2/Qo40xd+jEFNbfQD11Slc9YmKZBHf\nbW2g3NxclJaWorW1FSaTCU1NTd5zTU1NMJlMfp+nqKgIRUVF3uPGxkbvzxaLRXSsBkrXiTt7FsZu\n7w8ACDodWpOTgT7E5a2P2YxBLAumo+P7k3Y7OnbvhvPWW/sbdkQp/R6FmtrqA6ivTv2pT1pamt9z\nUdEMdf36dQg3xtLX1taC53kkJiYiKysL9fX1aGhogMvlQkVFBfLy8hSOltwk1wTlGjOm/x3RGo3s\nzG8tTdAjRDERubN47bXXcPz4cbS1tWHlypVYsmQJXC4XAGDRokXYu3cvvvjiC3AcB51Oh1WrVoFh\nGHAchxUrVmDdunXgeR4LFixAeg+L0pHI0lRXS8pcOTkheW7nlCme5T660Z46BXR2AnFxIXkNQkjw\nIpIsVq1a1eP54uJiFBcXy57Lzc1Fbm5uOMIi/WG3yy5H7hw3LiRP78rMhJCYKN41z+WC9vhxOKdN\nC8lrEEKCFxXNUCT2aOrqpAsHmkwQkpJC8wIsC+ctt0iKtd99F5rnJ4T0CiUL0ieczEQ81+jRIX0N\nh0xntramhlaiJUQBlCxIn8jN2pbbOrU/3CNHSteXcjqhqa0N6esQQgKjZEF6j+fB3Zg42Z0rIyO0\nr8Oysn0gcqOwCCHhRcmC9Bp75QoYu11UJhgM4FNTQ/5acqOrNCdP0rLlhEQYJQvSaxqfhf4AwJ2R\nAbCh/3VyZmdLNlBiW1rAXb4c8tcihPhHyYL0mu8S4gDgCnF/hZfBAJfM8uSa48fD83qEEFmULEiv\nyXZujxoVtteTW+6c+i0IiSxKFqRXmNZWsL4r/zIMXGGcWS/Xb8FdvAimpSVsr0kIEaNkQXqFu3BB\nUuYeNkyy8mwo8SYT+JQUSbmW9rggJGIoWZBe0cjshx3OJqib5JqiqN+CkMihZEF6RSNzZxG2zu1u\nZJPF6dPAjQUpCSHhRcmCBM/tBnfpkqQ4EsnCPXIkBJ/VZpmuLtlmMUJI6FGyIEFjGxok3+SFhAQI\nyckReHHWs1eGD63MMumEkNCjZEGCJjcRzp2W1vf9tnvJNXaspExuTw1CSOhRsiBBk2uCcg8fHrHX\nd2ZnS8q4S5fAtLdHLAZCBipKFiRofu8sIkQwmcAPGeJTKEBTUxOxGAgZqChZkOAIgnyyiOCdBQA4\nqSmKEEVQsiBBYZubpSvN6vXgzeaIxiHbb1FTQ6vQEhJmlCxIUOT6K/jU1LCsNNsTV2YmoBFvHc+2\ntIC9ejWicRAy0FCyIEGJhiYoAIBeD5fMjHFqiiIkvChZkKDIjoSKYOd2d3JNUTTfgpDw0gS+pP+2\nbt2KyspKJCUlYcOGDZLz33zzDT7++GMIgoC4uDg8/PDDGHXj2+PPf/5zGAwGsCwLjuNQUlISiZCJ\nD7a+XlKmVLJwZmfDd9lCzdmzgNMJaLWKxESI2kUkWcyfPx/FxcXYsmWL7PmhQ4fiueeeg9FoxKFD\nh7Bt2za8+OKL3vPPPvssBg0aFIlQiQzGZgPruxw4y8Idhm1Ug8GnpUEwGsHYbN8XOp3gzp+HW2aW\nNyGk/yLSDJWTkwOj0ej3/Lhx47zns7Oz0dTUFImwSJBkm6BSUiQdzRHDMPJLf5w+rUAwhAwMCv1v\n92/Hjh2YOnWqqGzdunUAgIULF6KoqMjvY8vLy1FeXg4AKCkpgcVi8Z7TaDSiYzWIVJ2YgwfB+uxX\nIYwbB0OIX7s39WGmTQN78qSoTH/5MoxR9h6r7fdObfUB1FencNUnqpLF0aNH8dVXX+H555/3lr3w\nwgswmUxoaWnBf/3XfyEtLQ05MjunAUBRUZEomTQ2Nnp/tlgsomM1iFSd4k+ehNZnjkVnYiIcIX7t\n3tSHGTIEg3xiwsmTaK2rk6xOqyS1/d6prT6A+urUn/qk9dAPGTWjoc6fP49f//rXeOqpp5CYmOgt\nN5lMAICkpCRMnz4dtbW1SoU4YCm9JpQcwWQCf+N34/tCAdyZM8oERIjKRUWyaGxsxKuvvorHHntM\nlNnsdjs6Ozu9Px85cgQZGRlKhTkwdXWBlelDcg8bpkAwYi6ZhQU19GWCkLCISDPUa6+9huPHj6Ot\nrQ0rV67EkiVL4LqxL8KiRYvw4Ycfwmaz4e233wYA7xDZlpYWvPrqqwAAt9uNH/zgB5gyZUokQiY3\ncNeuSZbS4JOTgSho6nGNGQPdvn2iMm1NDex+rieE9F1EksWqVat6PL9y5UqsXLlSUp6SkoL169eH\nKywSBNn5FVFwVwFAdkQU29AApqUFQlKSAhERol5R0QxFohfX0CAp44cOVSASKSEhQXZiIDVFERJ6\nlCxIj7grVyRlSk3Gk0P9FoREBiUL0iO51VzdKSkKRCJPrilKU1tLS5YTEmKULIh/nZ1gr18XlzFM\n1DRDAfCsQMtxojK2pQXstWvKBESISlGyIH7J9ldYLNG1WJ9eD9fIkZJiaooiJLQoWRC/ZPsroqgJ\n6iZXVpakjJIFIaFFyYL4xUZ55/ZNsv0Wp09TvwUhIdRjsvjkk09Ex9d926+JqnEyndt8FN5ZuNPT\nIeh0ojKms1N2dz9CSN/0mCz++Mc/io6feOKJsAZDootcsojGZihoNHCPHi0tpqYoQkKmx2Qh0G38\ngMW0t4NpaxMXsqyngzsK+R1CSwgJiR6TBcMwkYqDRBnZ/oqhQ5Xb8CgAuWTBnT0L3FiDjBDSPz3+\nz3c4HNi8ebP3uKurS3QMAI899lh4IiOKipX+ipvcw4ZBiIsDc2OVYgBgHA5wdXWyTVSEkN7pMVnc\nfffdouN/+qd/CmswJHrETH/FTSwLV1YWtEePioo1tbWULAgJgR6TxX333RepOEiUkV3mIwqHzXbn\nGjNGmixOn0bXwoUKRUSIegTVAF1ZWYmTJ0/CZrPBaDRiwoQJkn2yiYoIguyEvGhuhgL8dHKfPw90\ndQF6vQIREaIePSYLl8uFl156CdXV1cjMzERycjIuXbqEzz77DNnZ2fjP//xPaKK0w5P0HWOzgeno\nEBdqNODNZmUCChI/ZAj4pCSwLS3fF7rd0Jw7B9e4ccoFRogK9PhJ/+mnn6KtrQ3//d//DUu3IZON\njY1Yv349Pv30UyxevDjsQZLIkm2CGjIEYKN8wj/DeHbPO3hQVKypraVkQUg/9fi/f9++fXjwwQdF\niQIALBYLli1bhm+//TaswRFlyC4gGOVNUDfJrhNVU6NAJISoS4/Jor6+HmNk2oEBYMyYMbgi065N\nYp/c8t7uIUMUiKT3ZOdb1NdLm9UIIb0ScAa3zmfNnZv8lZPYF2tzLLoTBg8G75vYBMGzsCAhpM8C\ndnB/9dVXfpf9cLvdYQmKKEv2ziKKNjwKxJWdDZ1PHTQ1NXDecotCERES+3pMFtnZ2di1a1eP54O1\ndetWVFZWIikpCRs2bJCcFwQBZWVlOHToEPR6PR599FFkZmYCAHbu3ImPPvoIgGei4Pz584N+XdJL\nnZ3i0USAZ3e8KF0TSo5zzBjoKipEZbROFCH902OyeO6550L2QvPnz0dxcTG2bNkie/7QoUO4cuUK\nNm3ahJqaGrz99tt48cUXYbPZ8OGHH6KkpAQAsHr1auTl5cFoNIYsNvI9Tuaugjebo3ZNKDnuzEyA\nYUT7WbCNjWCuX4cweLCCkRESu3rss+jq6sIHH3yAl19+Gf/7v/8Lp9PZ5xfKycnp8QP+wIEDmDdv\nHhiGwdixY9He3o7m5mZUVVVh8uTJMBqNMBqNmDx5MqqqqvocB+kZKzMSKpaaoABAiI+He/hwSTmN\niiKk73r8ulhaWorTp09j6tSp2LdvH2w2G1asWBGWQKxWq2iIrtlshtVqhdVqhbnbZDCTyQSr1Sr7\nHOXl5SgvLwcAlJSUiJ5Po9FIhgDHunDUibHbwRoMojLdmDGIj8C/XSjrw+Tmgm1sFJXpr16NeHOa\n2n7v1FYfQH11Cld9ekwWVVVVePnll5GcnIzi4mI8++yzYUsWoVBUVISioiLvcWO3DwuLxSI6VoNw\n1Cm+thZau11U1qHXwxmBf7tQ1keTkoIEn3oIhw6h9Uc/8jRRRYjafu/UVh9AfXXqT33S0tL8ngvY\nDJWcnOwNoCOMY9VNJpOogk1NTTCZTDCZTGhqavKWW61WmEymsMUx0MXyhLzuXCNHSvpZmLY22dnp\nhJDAekwWbrcbR48e9f7heV50fNRnhc/+yMvLw65duyAIAqqrqxEfH4/k5GRMmTIFhw8fhs1mg81m\nw+HDhzFlypSQvS7pxuUC2y0x3xQrE/JEdDpPwvBB/RaE9E2PzVBJSUl44403vMdGo1F0zDCMZDMk\nf1577TUcP34cbW1tWLlyJZYsWQLXjV3MFi1ahKlTp6KyshKPP/44dDodHn30Ue9r3nPPPVizZg0A\n4N5776WRUGHCNjaKRhABAJ+UBMTFKRRR/7iysyWT8bQ1NXDMnatQRITErh6Thb9hrn2xatWqHs8z\nDIOHH35Y9lxBQQEKCgpCFguRJ9sEFYt3FTe4xo4FPv9cVMadOePZajWGhgITEg2ifBlREkmyq83G\nYH/FTe60NAjx8aIyxuEAd/68QhERErsoWRAv2Ql5MXxnAZaVXVhQS/0WhPQaJQviJXtnEWMT8ny5\nZJakoU5uQnqPkgXxEARwMmOzY3HYbHdOmWTBXbxIS5YT0ku96uXr6OiA3WeiE815UAe2uRnwWc5F\nMBggxPjIM8FkAm+xiGdz31iynFahJSR4QSWLI0eOYNu2bbgm06a9ffv2kAdFIk9uTSg+JSWis53D\nxTV2LHQ+d02a6mpKFoT0QlDJ4s0338Q999yDOXPm0KZHKiW7gGAsd253I7tkeXW1Z06JCpIhIZEQ\nVJ+F0+nEggULYDAYwLKs6A9RB9mRUDHeuX2TKysL8PldZZubZWerE0LkBfVp/6Mf/Qgff/yx3x3z\nSOxTw9LkfsXFwZ2eLinWnDqlQDCExKagmqFmzpyJdevW4c9//jMSExNF54Jd7oNEN9nZ22pJFgCc\nY8dKJuNpqqvhmDNHoYgIiS1BJYuNGzdi/PjxmD17NvVZqBDT3g6mvV1cqNGAV9FIN9f48cCXX4rK\nNKdPe0aAabUKRUVI7AgqWTQ0NODll1+mPgqVkm2Cslgk7fyxzD18OIT4eNH8CsbhgObcOdmJe4QQ\nsaA+DfLy8kK6HDmJLmpbQFAWy3oWFvShOXlSgWAIiT1B3Vk4nU688sormDBhApKSkkTnHnvssbAE\nRiJHdo6FivorbnKOGwetz/7t2upq2P1cTwj5XlDJIj09Hekyo0mIOsgNm1XLHIvu5O4s2KtXwTQ3\nQ7ixIyQhRF7AZMHzPMxmM+bOnQstdQSqkt/Z2yojJCbCPXw4uEuXROXa6mo4Zs5UKCpCYkPAPguW\nZfHee+9RolArp9OzLpQPt8WiQDDh5xo3TlJG8y0ICSyoDu5p06bhwIED4Y6FKIC9dk26lergwYBe\nr1BE4eWU6+SurQXcbgWiISR2BN3BvXHjRowdOxZmsxlMt/V0qIM7tql5mQ857lGjIBgMYLqtnszY\n7eDOnYM7K0vByAiJbtTBPcDJbnikws5tL5aFKzsb2u++ExVrT5ygZEFID4JKFvfdd1+44yAKkb2z\nUGHndnfO8eNlk4X9jjsUioiQ6Bf05kfHjh3D119/jebmZiQnJ2PevHmYNGlS0C9UVVWFsrIy8DyP\nwsJCLF68WHT+N7/5DY4dOwYAcDgcaGlpwW9+8xsAwNKlS5GRkQEAsFgseOaZZ4J+XdIzNS9N7o9r\nwgTP0uTd+mrYa9fANjaCV2nHPiH9FVSy+Nvf/ob/+Z//QUFBAbKzs9HY2Ihf/epXWLp0KYqKigI+\nnud5lJaWYu3atTCbzVizZg3y8vIwYsQI7zUPPvig9+fPPvsMZ8+e9R7rdDqsX7++F9UiQeF5+a1U\nVdxnAQCC0Qh3RoZ0YcHjx+GYN0+hqAiJbkEli08++QRr167FqFGjvGX5+fnYsGFDUMmitrYWqamp\nSLnRvJGfn4/9+/eLkkV3e/bswZIlS4IJjfQDe/26dCvVuLiY30o1GM7x4yXJQnviBCULQvwIKlm0\ntbVJPtjT0tJgs9mCehGr1Qqz2ew9NpvNqKmpkb322rVraGhoEDVxOZ1OrF69GhzH4a677sKMGTNk\nH1teXo7y8nIAQElJCSzdmhQ0Go3oWA36WyfmyhWwBoOoTBg9GjqFmqEi+h7NmQPu66/FZfX1iIuP\nB+LjQ/Yyavu9U1t9APXVKVz1CSpZjB8/Hu+99x7uv/9+6PV62O12fPDBBxgrM2a9v/bs2YNZs2aJ\nVrjdunUrTCYTrl69iueffx4ZGRlITU2VPLaoqEh0p9PYrYnFYrGIjtWgv3XSnzoFg128MpIjIQGd\nCv07RfQ90uuRGBcnmZDY8e23cN56a8heRm2/d2qrD6C+OvWnPmlpaX7PBTUp75FHHsH58+fx4IMP\n4pFHHsHy5ctx/vx5/Ou//mtQAZhMJjR128KyqakJJj97JVRUVGCOz4Y0N69NSUlBTk4Ozp07F9Tr\nkp7JDZtV+0goL4bxdHT70Jw4oUAwhES/oO4skpOT8ctf/hJNTU3e0VDdm5UCycrKQn19PRoaGmAy\nmVBRUYHHH39cct2lS5fQ3t4uumOx2WzQ6/XQarVobW3FqVOncNdddwX92sQ/uaXJVbOVahCcEyZA\nV1EhKtOeOIFOnlfVXh6EhELQQ2cBT1tYYmIiXC4Xrt74VpoSxDdRjuOwYsUKrFu3DjzPY8GCBUhP\nT8f27duRlZWFvLw8AJ4mqPz8fNEM8UuXLmHbtm1gWRY8z2Px4sV+O8ZJLwiCuvfdDoIrMxOCXg+m\nq8tbxnR2emZzZ2YqGBkh0SeoZFFVVYU33ngD169fl5zbvn17UC+Um5uL3NxcUdnSpUtFx3IjoMaN\nG4cNGzYE9RokeExLi2jJCwAQ9PqBtVS3VgvXmDHQ3pjf4y0+epSSBSE+gkoWpaWluOeeezB//nza\ng1sl/O6O1+2ubiBwTpwoTRbHjsH+j/844P4tCOlJUA2zNpsNCxcupEShIrJrQg2Uzu1uXDk5kv4J\ntrlZsucFIQNdUMmioKAAX331VbhjIREke2cxAJOFEB8Pl8wCglrac54QkaCaoWpqavDZZ5/h448/\nxuDBg0XnfvnLX4YlMBJenNydxQDq3O7OOWkSND6TRLXffQf7D39ITVGE3BBUsigoKEBBQUG4YyGR\nIggDe46FD+fEiYj785+lCws2NAzYfxNCfAWVLObPnx/mMEgkMTYbmM5OcaFGA34gjYTqRhg0CO6R\nI8H5TPbUfvcduihZEAIgyGSxY8cO2XKtVguz2Yzs7GzaozuGyDZBDRkyoCeiOSZNQpxvsjh6FF1B\nLJRJyEAQVLLYtWsXqqurkZSUBLPZjKamJrS0tCArKwsNNzpKn376aWTRTmMxgZqgpFyTJgGffioq\n4y5fBmu1gvezNA0hA0lQyWLEiBGYMWMGbr/9dm/Z559/jkuXLuH555/HRx99hHfeeQfr1q0LW6Ak\ndGSX+RjgyYI3meBOSwN3+bKoXHvkCLqoGZaQ4IbO7tmzB8XFxaKyRYsWYffu3WAYBnfeeScuXrwY\nlgBJ6NEcC3nOW26RlGkPHVIgEkKiT1DJIikpCQcPHhSVVVZWYtCgQQA8+01oNL1aZoooSHaOxQAd\nNtud3NLkXH29bHIlZKAJ6hN++fLl2LhxIzIyMrx9FhcuXMCTTz4JwDMPw/fOg0Qnpr0djO+mVRwH\nvherCKsVb7HAnZ4Orq5OVK49fBhdixYpFBUh0SGoZHHrrbfi9ddfR1VVFaxWK6ZOnYrc3FwkJiZ6\nz98awg1jSPjINkENGQJwnALRRB/HrbcizidZ6Kqq0LVwIU3QIwNa0GMlBw0ahJycHOTk5GDixIne\nREFiC3fliqSMmqC+57z1VklSYBsbaa0oMuAFdWfR3NyM1157DTU1NTAajWhra8PYsWPx7//+7353\nvCPRiauvl5S5hw1TIJLoJCQlwTV6NDRnzojKtYcOwU37qJABLKg7i7feegsjR47EO++8g23btqGs\nrAyjRo3CW2+9Fe74SIjJ3Vm4ZfYzH8icU6dKyrRHjgA8r0A0hESHoJLFqVOn8MADD8BgMAAADAYD\nfvKTn6C6ujqswZEQEwSwcs1QdGch4rzlFumy5S0t4M6eVSgiQpQXVLJISEiQzKO4fPky4uPjwxIU\nCQ+2uVm0hSgACAbDgF0Tyh8hPh7OceMk5Tqf4eOEDCRB9VnceeedeOGFF1BQUIAhQ4bg2rVr2Llz\np2RbVBLdWJn+Cj41lUb5yHBOnQrtiROiMu1336HzzjuBG3fYhAwkQSWLoqIipKamYvfu3bhw4QKS\nk5Px+OOP4xaZGa8kesl2blN/hSznxIkQDAbRPuVMVxd0330Hx/TpCkZGiDKCnnY9adIkTJo0yXvM\n8zy2b99OdxcxhJJFL2i1cE6ZAt3evaJi3f79lCzIgNTnNTrcbjc++uijoJNFVVUVysrKwPM8CgsL\nsXjxYtH5nTt34re//a13KG5xcTEKCwu95z766CMAwN133037a/SR7Ego6tz2yzF9uiRZcOfOgb12\nDfyQIQpFRYgyIrKgE8/zKC0txdq1a2E2m7FmzRrk5eVhhM+49fz8fDz00EOiMpvNhg8//BAlJSUA\ngNWrVyMvLw9GozESoauHwwG2sVFSTHcW/rlHjIA7NVWSZHUHDsB+220KRUWIMiKy201tbS1SU1OR\nkpICjUaD/Px87N+/P6jHVlVVYfLkyTAajTAajZg8eTKqqqrCHLH6cA0Nom1DAXhGQcXFKRRRDGAY\nOPLyJMXaykqac0EGnB7vLI4ePer3nMvlCvpFrFYrzN0WqjObzaipqZFct2/fPpw4cQLDhg3DsmXL\nYLFYJI81mUywWq2yr1NeXo7y8nIAQElJCSwWi/ecRqMRHatBb+rEVFeD9RnFI4wZA30U/ZtE5Xu0\ncCG4HTvEyaGrC4aGBgjd+vD8ico69YPa6gOor07hqk+PyeKNN97o8cGhDGjatGmYM2cOtFotvvzy\nS2zZsgXPPvtsr56jqKgIRd22wWzs1uxisVhEx2rQmzoZTp6EvtvIHgDoSkyEPYr+TaL1PYrPzITW\n54uT8/PP0RFEE1601qmv1FYfQH116k990tLS/J7rMVls2bKlTy/oy2QyoampyXvc1NQkWVOq+8KE\nhYWFeP/9972PPX78uPec1WpFTk5OSOIaSKhzu+8ceXmSZKE9cYK2XCUDSkT6LLKyslBfX4+Ghga4\nXC5UVFQgz6ctuLm52fvzgfBWIKAAABbJSURBVAMHvJ3fU6ZMweHDh2Gz2WCz2XD48GFMmTIlEmGr\nhyBItgsFqHM7WK7x46Wz3AUBum+/VSYgQhQQkdFQHMdhxYoVWLduHXiex4IFC5Ceno7t27cjKysL\neXl5+Oyzz3DgwAFwHAej0YhHH30UAGA0GnHPPfdgzZo1AIB7772XRkL1EtPWBqajQ1yo0dDwz2Cx\nLByzZ8Pw17+KinV//zvsCxcCOp1CgRESORHbCzU3Nxe5ubmisu5zNH784x/jxz/+sexjCwoKUFBQ\nENb41Ex2Mt7QoZLF8oh/junTYfjiC6DbwA6msxO6w4dpkh4ZEOjTYgCQ27iH+it6R0hIgEOm+VO3\nZ49kSDIhakTJYgDgfFYMBkAb+fSBIz9fUsZdvgzu/HkFoiEksihZDACULELDPWIE3BkZknL9N98o\nEA0hkUXJQuWYtjaw16+LC1kW7h7GUxP/umTuLrRHj4K9dk2BaAiJHEoWKid7V5GSAmi1CkQT+5yT\nJ4NPShIXCgL0X3+tTECERAglC5Xj6uokZe70dAUiUQmNBl1z50qKdQcPgmlpUSAgQiKDkoXKyY6E\nov6KfnHMmAHBdwFGtxv63buVCYiQCKBkoWaCAA3dWYSewSA7Mkq3dy+Yzk4FAiIk/ChZqBjT0gLG\nZhMXajSePgvSL11z5kj6fZiuLs+8C0JUiJKFisl2bg8bBmgiNnFftQSjEV0zZkjK9bt2SZdWIUQF\nKFmomEYuWVATVMh0zZ0rWTKFsdtpZBRRJUoWKiY3Eso1fLgCkaiTYDLBIXd3sXs3mNZWBSIiJHwo\nWaiVIMg3Q9GdRUjZCwulc1acThh27FAmIELChJKFSrFNTZKROYJOB37oUIUiUichKUl2Vrdu3z6w\nfrb/JSQWUbJQKdm7iuHDaVnyMOiaPx+Cz/7mcLth+PxzZQIiJAzok0OlNGfPSspoMl54CAkJsrO6\ntVVVQE2NAhEREnqULFRKc+aMpMw1erQCkQwMXXPnQpDZwZH9wx8AnlcgIkJCi5KFCjFtbWCvXpWU\nuzMzFYhmgDAYYC8ulhQzFy/SXt1EFShZqJDcXYV72DAI8fEKRDNwOKZPl23qM3zxhXQmPSExhpKF\nCsk2QdFdRfgxDDoXL5YWd3bC8Je/KBAQIaFDyUKFZJNFVpYCkQw87owMOPLyJOW6gwehOXFCgYgI\nCY2ILRJUVVWFsrIy8DyPwsJCLPb5Bvbpp5/ib3/7GziOw6BBg/Czn/0MQ4YMAQAsXboUGTe2s7RY\nLHjmmWciFXbMYWw26q9QmP2226A9ehSM3S4qj/voI9iefFK6vDkhMSAiyYLneZSWlmLt2rUwm81Y\ns2YN8vLyMKJb++6oUaNQUlICvV6PL774Au+//z6eeOIJAIBOp8P69esjEWrMo/4K5QmJibD/6EeI\n++MfReVsSwsMn3yCzqVLFYqMkL6LSDNUbW0tUlNTkZKSAo1Gg/z8fOzfv190zaRJk6DX6wEA2dnZ\nsNLs1z6h/oro4JgxA66xYyXluoMHoTl+XIGICOmfiNxZWK1WmM1m77HZbEZND5OVduzYgSlTpniP\nnU4nVq9eDY7jcNddd2GGzOJtAFBeXo7y8nIAQElJCSwWi/ecRqMRHauBXJ24+nrAZzaxbto0GGOg\n7qp7j376U7Dr1sFnbjcM//d/cN9yC5CcrEhY/aG69wjqq1O46hN1Gxvs2rULZ86cwXPPPect27p1\nK0wmE65evYrnn38eGRkZSE1NlTy2qKgIRUVF3uPGxkbvzxaLRXSsBr51Ymw2DJKZud2anAwhBuqu\nxvdoyN13w/H22+JCux3uTZtg++lPY25vETW+R2qrU3/qk5aW5vdcRJqhTCYTmpqavMdNTU0wmUyS\n644cOYI//elPePrpp6HttpLnzWtTUlKQk5ODc+fOhT3mWKQ5fVpS5h42DEJCggLREAAQZs+Ga/x4\nSTl3/jwMf/2rAhER0jcRSRZZWVmor69HQ0MDXC4XKioqkOczvPDs2bN466238PTTTyMpKclbbrPZ\n4HQ6AQCtra04deqUqGOcfE9z8qSkjPorFMYw6FiyBHy33+mb9Lt3Q3vkiAJBEdJ7EbkH5jgOK1as\nwLp168DzPBYsWID09HRs374dWVlZyMvLw/vvvw+73Y6NGzcC+H6I7KVLl7Bt2zawLAue57F48WJK\nFnJcLmhlOk5d2dkKBEO6E4xGdNx/P4xvvilZJyp++3bYBg+G+8bQcEKiFSMIgqB0EOFy+fJl789q\na5cExHXSnDqFhNJS0XlBr0frL34h3ZwnSqn9PdLt2oW4Tz+VXCMYjbD9/Ofguw0CiVZqf4/UIKb7\nLEj4aY8elZS5xo+PmUQxEDjmzoVz0iRJOWOzIeGdd8B0dCgQFSHBoWShBjwP7bFjkmLnLbcoEAzx\ni2HQsXSpZxMqH+y1a0goLZXsbkhItKBkoQLc2bPSVU21WjjHjVMmIOKfXo/25cvBy8yx4OrqkPD2\n25QwSFSiZKECck1Qzuxs4MaMeBJdhEGD0L58uXQrVtxIGG+9RQmDRB1KFrFOEOSTBTVBRTU+NRXt\ny5bJ9ilxFy8i4c03wVy/rkBkhMijZBHjuAsXwLa0+BRycOXkKBMQCZo7KwvtDz4onzDq65G4eTO4\nS5ciHxghMihZxDhdZaWkzDVmDC2DHSNc2dloX75cNmEwra1IePNNaGTuHAmJNEoWsay9HboDByTF\ncsMzSfRyjRnj6cOQ6WNiurqQ8N57MHzyCeByKRAdIR6ULGIYs2sXcGMplJsEgwGOW29VKCLSV64x\nY2BbuRLCoEGy5/W7d8P4xhtgu62xRkgkUbKIVU4n2K+/lhQ7ZsyQLFFOYgM/fDjaHnsM7mHDZM9z\ndXVI3LgR+l27JMuGEBJulCxilK6qCmhtFReyLLrmzFEmIBISwuDBsP3sZ3D6G6DgdMLw6acwbt4M\n7uLFyAZHBjRKFrFIEDzfLn04J0+GEIMb6hAfBgM6li1D5x13AKz8f1Hu4kUYN21C3O9/D4Z2lSQR\nEFs7rxAAgObECbBXr0qam7rmzVMoIhJyDAPHvHlwjxqF+N/9Dmxzs+xluspK6I4cQdfMmeiaN4++\nLJCwoTuLWGO3I+7jjyXFrsxMuGnpdtVxZ2Sg7YknPM2LDCN/kcsF/Z49GPTyy4jbvh0szc0gYUB3\nFjEm7q9/lf2W2TV3rgLRkIgwGGC/6y44p0xB/Icfeu4q5fA8dAcPQnfwINzp6eiaNQvOyZNp2RcS\nEpQsYoimuhq6vXsl5a7Ro2nG9gDgHjkSbatWQbdvHwxffgmmvd3vtVxdHeLr6oA//xnOCRPgmDIF\nrrFjAZ0ughETNaFkESOYzk7EffihpFzQ6dC5ZIn/JgqiLhwHR34+HLm5MOzcCd2ePWC6uvxf73RC\ne+SIZ/tWjQaurCw4x4+HKzsb/JAh9HtDgkbJIgYwLS1IKC0FK7OwnP3222NihzUSYgYD7MXF6Jo3\nD7qKCuh37w68eZLLBc2pU9CcOgUAEBIT4crMhCsjA+70dM8+G7RZFvGDkkWUY69e9ZsohHHj4Jg9\nW4GoSLQQ4uPRVVSErrlzoTt0CLq9e8F12064J0xbG7SHD0N7+LCngGXhHjoUfGoq3Kmpnp+HDPF8\nGdHQR8VAR78BUYrp7ITum2883xjtdsl5wWAA/5OfAOrdQp30hl4Px6xZcMycCe7iRej274f2yJHe\nbdXK8+CuXAF35QpE9xcMA37wYPDJyWAzMmDgOPBJSeAHDYIwaBB4oxGC0UgJReXo3Y0iTEsLNOfO\nQXPmDLSHDskmCQAQ9Hp0PPAAdGYzoKKN5kkIMAzc6enoTE9H5113QVNb67l7OHlSuptisAQBbHMz\n2OZmMJcvQ+/v99JggBAf//2fuDhPmcEAQa/3LJRoMEDQaiHodIBO5/lZqwW0WggajedvjvMkHo2G\n+lSiSMSSRVVVFcrKysDzPAoLC7F48WLReafTic2bN+PMmTNITEzEqlWrMHToUADAn/70J+zYsQMs\ny2L58uWYMmVKyOPTVVSAa2jo/xN1/6bv+62f5z1/BAGMwwG4XGC6usC2toJpa+u5o/LmUxqNaF+x\nguZUkMA4Dq5x4+AaNw6dggDu0iVoTpyA5vRpaC5cCPkqtozd7vmCE8oZ5SzrSSIsC3AcwHGeZMIw\nnmOWhcCynuOb1zCMp+zG473nbiQeofuseJYFazQivr3dc757crr5M8NA6PaztOIyj+nLNf6u7yU2\nMRFxbW1wp6fDMX16n5/HV0SSBc/zKC0txdq1a2E2m7FmzRrk5eVhRLcPvB07diAhIQGvv/469uzZ\ng9/97nd44okncPHiRVRUVGDjxo1obm7GCy+8gF/96ldg/SyD0FfakyehOXkypM8ZarzZjPaHH6YO\nbdJ7DAP3iBFwjxiBroULAacT3IUL0NTVgbtwAdzFi7L9Yorjec8Xq25Cfa/BGAzQ+rlbikWMwQCd\n3Q5nZ2fsJYva2lqkpqYiJSUFAJCfn4/9+/eLksWBAwdw3333AQBmzZqFd955B4IgYP/+/cjPz4dW\nq8XQoUORmpqK2tpajB07NhKhRweOgyMvD/bbboMQH690NEQNtFq4s7LgzsryFjHt7WDr6z39Flev\ngm1s9Pzx3YmRDEgRSRZWqxXmbt+GzWYzampq/F7DcRzi4+PR1tYGq9WK7Oxs73UmkwlWP7e55eXl\nKC8vBwCUlJTAYrF4z2k0GtGxLzYhAUy0Le2t0UDIzwe/aBG0JhMSJKd7rlOsUVt9gBirk8UCjBwp\nLXc4gKYmMFYrmOvXoWtq8uwPfv060NYGpq0NaGuL2cEWLMvCEG3/9/vhZn30RiPiQ/i7p6oO7qKi\nIhQVFXmPG7t1/losFtGxr4T2dmiUvhXVaOAeOtSzztPo0XBlZXnuJHhetiM7UJ1ijdrqA6ioTlot\nkJICy8SJ8vXheU+fRXu754/dDqaz0/Onq8vTH3fjb8bh8DQtORxgnE7A6fT87XZ7/na5wLjdEdsZ\n0GAwwK70//0Qulkfp82Gjl7+7qWlpfk9F5FkYTKZ0NRth6+mpiaYTCbZa8xmM9xuNzo6OpCYmCh5\nrNVqlTw2FLpmzYJz/PiQP6+oo4plPR1ygGckiEYD6HSeoYeDBnkSA43+ILGIZb2joDBkSGieUxAA\nt9uTRHjekzx4/vvjm3/cbjCC8P31guBJXjfLbhwD8A4wufn8jCBAN3gwOpqbv7+2++t3u1a2HH76\nUHoa6BKovKd/jyDokpLQ2dIS8r7NiCSLrKws1NfXo6GhASaTCRUVFXj88cdF10ybNg07d+7E2LFj\nsXfvXkycOBEMwyAvLw+bNm3CHXfcgebmZtTX12PMmDEhj5HWViIkyjCMdwit78dkKBu8BIsFTjXc\n/d0gWCxwhKE+EUkWHMdhxYoVWLduHXiex4IFC5Ceno7t27cjKysLeXl5KCgowObNm/Fv//ZvMBqN\nWLVqFQAgPT0ds2fPxpNPPgmWZfHQQw+FfCQUIYSQnjGCEKO9UkG43G3ZA9W0HXejtjqprT6A+uqk\ntvoA6qtTf+rTU58FfUUnhBASECULQgghAVGyIIQQEhAlC0IIIQGpuoObEEJIaAyYO4vVq1crHULI\nqa1OaqsPoL46qa0+gPrqFK76DJhkQQghpO8oWRBCCAmIe+65555TOohIyczMVDqEkFNbndRWH0B9\ndVJbfQD11Skc9aEObkIIIQFRMxQhhJCAKFkQQggJSFWbHwFAVVUVysrKwPM8CgsLsXjxYtF5p9OJ\nzZs348yZM0hMTMSqVaswdOhQhaINLFB9du7cid/+9rfePT6Ki4tRWFioRKhB2bp1KyorK5GUlIQN\nGzZIzguCgLKyMhw6dAh6vR6PPvpo1LcnB6rTsWPH8Morr3h/z2bOnIl777030mEGrbGxEVu2bMH1\n69fBMAyKiopw++23i66JpfcpmPrE2nvkcDjw7LPPwuVywe12Y9asWViyZInompB/1gkq4na7hcce\ne0y4cuWK4HQ6hf/4j/8Q6urqRNd8/vnnwq9//WtBEARh9+7dwsaNG5UINSjB1Oerr74S3n77bYUi\n7L1jx44Jp0+fFp588knZ8wcPHhTWrVsn8DwvnDp1SlizZk2EI+y9QHU6evSo8NJLL0U4qr6zWq3C\n6dOnBUEQhI6ODuHxxx+X/N7F0vsUTH1i7T3ieV7o7OwUBEEQnE6nsGbNGuHUqVOia0L9WaeqZqja\n2lqkpqYiJSUFGo0G+fn52L9/v+iaAwcOYP78+QCAWbNm4ejRoxCitI8/mPrEmpycHBiNRr/nDxw4\ngHnz5oFhGIwdOxbt7e1obm6OYIS9F6hOsSY5Odl7lxAXF4fhw4dL9r2PpfcpmPrEGoZhvPuGu91u\nuN1uMD67bIb6s05VzVBWqxXmblsJms1m1NTU+L2G4zjEx8ejra0NgwYNimiswQimPgCwb98+nDhx\nAsOGDcOyZctgCeEm7ZFmtVpF8ZvNZlitViQnJysYVf9VV1fjqaeeQnJyMv7lX/4F6enpSocUlIaG\nBpw9e1ayO2Wsvk/+6gPE3nvE8zyeeeYZXLlyBT/84Q+RnZ0tOh/qzzpVJYuBaNq0aZgzZw60Wi2+\n/PJLbNmyBc8++6zSYZFuRo8eja1bt8JgMKCyshLr16/Hpk2blA4rILvdjg0bNuDBBx9EfHy80uH0\nW0/1icX3iGVZrF+/Hu3t7Xj11Vdx4cIFZGRkhO/1wvbMCjCZTGhqavIeNzU1eTt+5a5xu93o6OhA\nYmJiROMMVjD1SUxMhFarBQAUFhbizJkzEY0x1Ewmk2iXL7k6x5r4+Hhvk0Fubi7cbjdaW1sVjqpn\nLpcLGzZswNy5czFz5kzJ+Vh7nwLVJxbfo5sSEhIwceJEVFVVicpD/VmnqmSRlZWF+vp6NDQ0wOVy\noaKiAnl5eaJrpk2bhp07dwIA9u7di4kTJ0ra+qJFMPXp3k584MABjBgxItJhhlReXh527doFQRBQ\nXV2N+Pj4qG/aCOT69evetuLa2lrwPB+1X1AAz0inN998E8OHD8cdd9whe00svU/B1CfW3qPW1la0\nt7cD8IyMOnLkCIYPHy66JtSfdaqbwV1ZWYl3330XPM9jwYIFuPvuu7F9+3ZkZWUhLy8PDocDmzdv\nxtmzZ2E0GrFq1SqkpKQoHbZfgerzwQcf4MCBA+A4DkajEQ8//LDklyaavPbaazh+/Dja2tqQlJSE\nJUuWwOVyAQAWLVoEQRBQWlqKw4cPQ6fT4dFHH0VWVpbCUfcsUJ0+//xzfPHFF+A4DjqdDg888ADG\njRuncNT+nTx5Er/4xS+QkZHh/XD553/+Z++dRKy9T8HUJ9beo/Pnz2PLli3geR6CIGD27Nm49957\nw/pZp7pkQQghJPRU1QxFCCEkPChZEEIICYiSBSGEkIAoWRBCCAmIkgUhhJCAKFkQQggJiJIFIYSQ\ngChZEEIICYiSBSERcOXKFSxfvty7dpfVasVDDz2EY8eOKRwZIcGhZEFIBKSmpuL+++/H66+/jq6u\nLrzxxhv4h3/4B0ycOFHp0AgJCi33QUgEvfzyy2hoaADDMHjppZe8KwYTEu3ozoKQCCosLERdXR2K\ni4spUZCYQsmCkAix2+149913UVBQgD/84Q+w2WxKh0RI0ChZEBIhZWVlyMzMxMqVK5Gbm4tt27Yp\nHRIhQaNkQUgE7N+/H1VVVXjkkUcAAMuWLcPZs2fxzTffKBwZIcGhDm5CCCEB0Z0FIYSQgChZEEII\nCYiSBSGEkIAoWRBCCAmIkgUhhJCAKFkQQggJiJIFIYSQgChZEEIICej/A8xTGMwoP4rZAAAAAElF\nTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jJuPh4y2rUAX",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        },
        "outputId": "6a77e6a6-3e69-4dc1-b65b-1e7519b39d19"
      },
      "source": [
        "x = np.linspace(0.1, 3, 100)\n",
        "plt.axes(xlabel=\"x\", ylabel=\"Lognorm PDF\")\n",
        "plt.plot(x, lognormal_cdf(x, mean, sigma),'r-', lw=5, \n",
        "         alpha=0.6)\n",
        "plt.show()"
      ],
      "execution_count": 29,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXAU55k/8G/PrQuhGYEESFwCYUAY\nI2TOGKPDmLK9MYmNib2xcdikKks5GLuyzpLCC1kva3yAvTHYiR1FZTvZXZWDlxBc5ueIwxhkMEKI\n+5CAcEoIaRDomtHMdP/+EJ5o1KPRCM1MT/d8P1UqTT9vj+Z5a6R51G93v68gSZIEIiIiADqlEyAi\noujBokBERF4sCkRE5MWiQEREXiwKRETkxaJAREReBqUT6K+rV6/6bKempqKhoUGhbEJPa/0BtNcn\nrfUH0F6ftNYfoH99Gjp0aI9tPFIgIiIvFgUiIvJiUSAiIi8WBSIi8mJRICIir4hcffTuu++isrIS\nycnJWLdunaxdkiSUlJTg0KFDMJvNWLp0KUaPHh2J1IgoFCSp80sUAY8HkCQIkuR97NMuSRAA72N8\n+/jbn9P1MdC5b5ftHr93f9x1+8YN6G/c6L0PdyqE84oKve/SqaEBhqYmuDMzAbM5ZK8fkaIwd+5c\nzJ8/Hxs3bvTbfujQIdTV1eHXv/41qqur8bvf/Q7/+Z//GYnUiLRLFCG0t0Noa4PgcHi/4HBA6OiA\n4HRCcLmAb7+7XJ3f3W7oTCYk3rwJuN2AxwPh9nd4PBBuf4codj4Wxb9/kEcpvcWCRIdD6TRCSm+x\nIMHhQPPPfw5x8OCQ/dyIFIUJEyagvr6+x/aKigrMmTMHgiAgOzsbra2tuHHjBlJSUiKRHpG6tLdD\nd/MmdE1N0N28CeHmTeiamyE0N0PX2gqhpQVCa2tnAbhDgsUCvcY+RCk4UXHzmt1uR2pqqnfbZrPB\nbrf7LQplZWUoKysDAKxdu9bneQBgMBhkMTXTWn8A7fUpLP1xu4G6Ogh1dZ3f6+uB69chXL8OtLUF\n9zMsljt+eZ1OB0s/nh9ttNYf4O99MlqtQAh//6KiKPRFUVERioqKvNvd7+jT2p2LWusPoL0+9bs/\nHg/0V69Cf/Ei9Jcvdz6+dk3RIRmLxQKHho4UtNYf4O99arbbIRr69lEe6I7mqCgKVqvV54+qsbER\nVqtVwYyIwsjjgf7SJRiqq2E4exaGS5cAl0vprPpPEAC9HhAESDodoNN1xr79/m1cuH0qtevjb/cD\nIN3e1+fndn3c9Tn+4t1zAiDFx8Pd3u7/Z0ajIPKTEhLgbm2FZDSG9KWjoijk5eVh27ZtmD17Nqqr\nqxEfH8/zCaQpQns7DKdOwXjiBAynT/drvL8vJIsFUlxc55fF4v2C2QzJbO78QLn9XTIagdvfjYMG\noaW5uXNbrwcMBsBg6PxQ1+sBvf7vj7sWgChlSU1Fq4aOTgEgLkx9ikhRePvtt3HixAk0Nzfjpz/9\nKZ544gm43W4AwLx58zBlyhRUVlZi2bJlMJlMWLp0aSTSIgqvjg4YT56E8dAhGE+dCt1wkF4PMTkZ\n4sCBEFNSIA0YADEpqfN7YiKkhARIiYmQ4uK8/333WWoqPBr7EKXgRKQoLF++PGC7IAj48Y9/HIlU\niMJOV1cHc3k5jFVV/ToikAYMgCctDZ7BgyGmpUFMTYXHaoU0cOCdf9gT9SIqho+IVE+SYDhxAuav\nvoLh3Lm+Pz0hAe4RI+AZPhyeYcPgGTYMUmJiGBIlCoxFgag/RBHCwYNI/PRT6Gtrg36aZLHAPWaM\n90scNCiqx+QpdrAoEN0hw+nTsHz2GXRNTUHd6CUmJ8OVkwP3hAlwjxrVefKWKMrwt5Koj3TXriFu\n61YYTp/uDAS4KUoym+G6+264pkyBe/RongugqMeiQBQstxvmHTtg2bGj1yuJxLQ0OGfOREdubr/u\nLCaKNBYFoiDoL1xA/J/+BN21awH384wcCUdhIdzZ2TxHQKrEokAUiCjCvHMnLF98EXB6ZM+IEWh/\n8EF4srJYDEjVWBSIeiA0NyP+f/4HhpqaHvcRU1IgPvUUWoYPZzEgTWBRIPJDf+ECEj76CEJzs/8d\ndDo4CgrgzM+HecgQgHf/kkawKBB1Yzx8GPGlpZ3TV/vhychA2+OPQwww0ySRWrEoEH1LkjrPH2zb\n5r9dEODMz4dj3jxeWkqaxaJABACSBMuWLTDv3eu/OTERbT/4QedVRUQaxqJAJEmI+/RTmPbv99vs\nychA6+LFkJKTI5wYUeSxKFBsE0XE/elPMFVU+G125eSgbdEiwGyOcGJEymBRoNglSYjbtKnHguCc\nMweOhx/mpaYUU1gUKGZZ/t//g+nAAb9tjnnz4OyyFjhRrGBRoJhk2rsX5h07/LY5HnoIzrlzI5sQ\nUZRgUaCYYzx6FHFbtvhta3/kEXTMmRPhjIiiBy+2ppiiu3oV8f/7v37nMXIWFrIgUMxjUaCYIbS1\nIeHDDwGXS9bWce+9nTelEcU4FgWKDaKI+D/+EbobN2RNrgkT0P7YY7zKiAgsChQjLNu2wVBdLYt7\nhgxB25NPctoKotv4l0CaZ6iuhnnXLllciotD2zPP8MY0oi5YFEjThNbWzhlPZQ0C2p56CqLNFvmk\niKIYiwJplyQh7pNPINy6JWtyPPgg3OPGKZAUUXRjUSDNMn3zDYwnTsji7jFj4MzPVyAjoujHokCa\nJNjtsGzdKotLcXGdE9zxSiMiv1gUSHskCfGbN0NwOmVN7Y8/zimwiQJgUSDNMR4+DMOpU7J4x733\nwjVpkgIZEakHiwJpitDW5ndeIzE5Ge2PPKJARkTqwqJAmmL57DMILS2yePujjwJxcQpkRKQuEZsl\ntaqqCiUlJRBFEYWFhViwYIFPe0NDAzZu3IjW1laIooinnnoKubm5kUqPNEB/4YLf9RFcOTlw5+Qo\nkBGR+kSkKIiiiOLiYqxcuRI2mw0rVqxAXl4eMjIyvPts2rQJM2fOxLx583D58mW8+uqrLAoUPEny\nO2wkWSydRwlEFJSIDB/V1NQgPT0daWlpMBgMmDVrFg50+49OEAS0tbUBANra2pCSkhKJ1EgjjJWV\n0F+6JIs75s/n1UZEfRCRIwW73Q5bl+kEbDYbqrtNTrZw4UL8x3/8B7Zt2wan04mXX37Z788qKytD\nWVkZAGDt2rVITU31aTcYDLKYmmmtP0AY+uR0Qv/ll4DF4hOWMjJgfOSRsE92x/co+mmtP0D4+hQ1\nK6/t3bsXc+fOxT/8wz/gzJkzeOedd7Bu3Trouv1BFxUVoajL2rkNDQ0+7ampqbKYmmmtP0Do+2T+\n4gtY6upk8ZaCAnjs9pC9Tk/4HkU/rfUH6F+fhg4d2mNbRIaPrFYrGhsbvduNjY2wWq0+++zYsQMz\nZ84EAGRnZ8PlcqG5uTkS6ZGKCU1NsHz5pSzumjgRnjFjFMiISN0iUhSysrJQW1uL+vp6uN1ulJeX\nIy8vz2ef1NRUHDt2DABw+fJluFwuDBgwIBLpkYpZtm+Xr6Sm18PBexKI7khEho/0ej2WLFmCNWvW\nQBRF5OfnIzMzE6WlpcjKykJeXh6eeeYZ/Pa3v8Vnn30GAFi6dCkEzk9DAegaG/1egur8znc4JTbR\nHYrYOYXc3FzZJaaLFi3yPs7IyMArr7wSqXRIA8x//Ssgij4xyWKBs6BAoYyI1I93NJMq6a5dg+nQ\nIVncef/9kHjnMtEdY1EgVbJ88QUgST4xKSEBzu98R6GMiLSBRYFUR3flCoxHj8rijvx8rrdM1E8s\nCqQ6lp07ZTExORkdty9pJqI7x6JAqqK7ft3vUYKzoAAwGhXIiEhbWBRIVcy7dsnOJYjJyei4915l\nEiLSGBYFUg2hqQmmykpZvOO++wBD1MzYQqRqLAqkGuavvgI8Hp+YFBcH57RpCmVEpD0sCqQKQlsb\nTN98I4s7Z8+WzY5KRHeORYFUwfT11xCcTp+YZDKhg/clEIUUiwJFP7cb5vJyWbhj2jRI8fEKJESk\nXSwKFPWMR45A6D6Nuk4H5333KZMQkYaxKFB0kySY9+yRhV05OZC4ZCtRyLEoUFTTX7gA/eXLsjjn\nOCIKDxYFimrmvXtlMU9GBjwjRiiQDZH2BSwKW7Zs8dluamoKazJEXQlNTf6ntPjOdwAuwEQUFgGL\nwqZNm3y2X3jhhbAmQ9SVed8++SI6SUlw3X23QhkRaV/AoiB1m2OGKGLcbv9Lbc6cySktiMIoYFHg\nGsmkFOOJE34vQ+2YPl2ZhIhiRMB/uTo6OrBhwwbvttPp9NkGgOeeey48mVFMM+3bJ4u5Jk2ClJSk\nQDZEsSNgUfj+97/vs/29730vrMkQAYCuoQGGmhpZnEcJROEXsCgsXLgwUnkQeZn275fFxNRUuLOy\nFMiGKLYEdcausrISp06dQktLCxITEzF+/HhMmTIl3LlRLOrhBHPHjBm8DJUoAgIWBbfbjVdffRVn\nzpzB6NGjkZKSgitXruDzzz/H2LFj8ctf/hIGXglCIWQ8ehRCW5tv0GBAx9SpyiREFGMCfqJv3boV\nzc3NeOutt5CamuqNNzQ04I033sDWrVuxYMGCsCdJscPfUYJr0iRICQkKZEMUewJekrp//348++yz\nPgUBAFJTU7F48WJ8/fXXYU2OYotgt/s9wezkCWaiiAlYFGprazFmzBi/bWPGjEFdXV1YkqLYZDp4\nUBYTU1PhGTVKgWyIYlOvdzSbTCa/bT3Fie6IJPktCh15eTzBTBRBvZ5o3rlzZ4/TXXi6LaJOdKf0\n585BZ7f7BgUBHbm5yiREFKMCFoWxY8di9+7dAduJQsFUUSGLuceOhTRwoALZEMWugEVh9erVIXuh\nqqoqlJSUQBRFFBYW+r1qqby8HJ988gkEQcCIESPw/PPPh+z1KYo5HH6nyO7Iy1MgGaLYFrAoOJ1O\nbNq0CZcuXcKoUaPwve99D0ajsc8vIooiiouLsXLlSthsNqxYsQJ5eXnIyMjw7lNbW4vNmzfjlVde\nQWJiIm7evNn33pAqmY4ehdDR4ROTLBa4Jk5UKCOi2BXwRHNxcTEOHjyIYcOGYf/+/fj444/v6EVq\namqQnp6OtLQ0GAwGzJo1Cwe6XY++fft2PPjgg0hMTAQAJCcn39FrkfoY/Zxgdt1zD3AH/4AQUf8E\nPFKoqqrCa6+9hpSUFMyfPx+rVq3CkiVL+vwidrsdNpvNu22z2VBdXe2zz9WrVwEAL7/8MkRRxMKF\nC3HPPffIflZZWRnKysoAAGvXrpXdQ2EwGGQxNdNaf4BufbLbob96FbBYfPYxzpuHBJX0W/PvkQZo\nrT9A+PrU6/BRSkoKgM4b1tq6Tz8QQqIoora2FqtWrYLdbseqVavw5ptvIqHbnaxFRUUoKirybjc0\nNPi0p6amymJqprX+AL59Mu/aBYvD4dMuDhqE5oQEQCX91vp7pAVa6w/Qvz4NHTq0x7aARcHj8eDY\nsWPebVEUfbYBICcnp9cErFYrGhsbvduNjY2wWq2yfcaOHQuDwYDBgwdjyJAhAW+eI20wVlbKYh1T\npvDeBCKFBCwKycnJeO+997zbiYmJPtuCIMgW3fEnKysLtbW1qK+vh9VqRXl5OZYtW+azz7Rp07Bn\nzx7k5+fj1q1bqK2tRVpaWl/7Qyqiu3oVej93xbs4Ay+RYgIWhY0bN4bkRfR6PZYsWYI1a9ZAFEXk\n5+cjMzMTpaWlyMrKQl5eHiZPnozDhw/jhRdegE6nww9/+EMkcZUtTTMdOiSLeUaMgNjl/BMRRVbE\n5r3Ozc1Fbre7UxctWuR9LAgCFi9ejMWLF0cqJVKSJMFYVSULd/AogUhRAS9JJQoX/dmz0HW/F0Wn\ng+vuu5VJiIgAsCiQQvwNHbnGjYN0+z4VIlIGiwJFntsNY7er2ACeYCaKBn06p9DW1gZHt2vKu19a\nStQb4eRJCO3tPjHJbIZrwgSFMiKibwVVFI4cOYL3338f169fl7WVlpaGPCnSNsHPtBbu8eMBrtFB\npLigisJvfvMbPPbYY5g9ezYX16H+cbkgHD4sC3dMnqxAMkTUXVBFweVyIT8/HzodT0FQ/xhPnQKc\nTp+YZLHAPW6cQhkRUVdBfco//PDD+POf/9zjCmxEwTL6OUpw5eQAhojdMkNEAQT1lzh9+nSsWbMG\nmzdvlt1lHMw0F0QAAKcTxpMnAb3eJ+zi0BFR1AiqKKxfvx533XUXZs6cyXMKdMeMJ08CLpdPUZDi\n4+HmpIdEUSOoolBfX4/XXnuN5xSoX/xNa+GaNEl25EBEygnqUz4vL082ZTZRnzgcMJ45IwvzqiOi\n6BL01Uevv/46xo8fL1sm87nnngtLYqQtxhMnALfbJyYlJcEzerRCGRGRP0EVhczMTGRmZoY7F9Iw\n49GjspgrJwfgkCRRVOm1KIiiCJvNhvvuuw9GLqROd8LhgPH0aVm4Y9IkBZIhokB6/TdNp9Pho48+\nYkGgO2Y8dUo+dJSYyKEjoigU1LH71KlTUVFREe5cSKOMR47IYhw6IopOQZ9oXr9+PbKzs2Gz2SB0\nWVSdJ5opIKfT79ARF9Mhik480UxhZTx1qvOGta4SE+Hm0BFRVAqqKCxcuDDceZBG+Rs6ku65h0NH\nRFEq6FnIjh8/ji+//BI3btxASkoK5syZg5ycnHDmRmrndHYeKXQjcoU1oqgV1L9r27dvx1tvvYWB\nAwdi2rRpSElJwX/913+hrKws3PmRihnPnJENHUnx8QCnySaKWkEdKWzZsgUrV67EyJEjvbFZs2Zh\n3bp1KCoqCldupHJ+b1ibOBEmDh0RRa2g/jqbm5uRkZHhExs6dChaWlrCkhRpgNsNg5+hIxdvWCOK\nakEVhbvuugsfffQRnLdXzHI4HPj444+RnZ0d1uRIvQxnzkBwOHxiksXCabKJolxQw0c/+clP8Pbb\nb+PZZ59FYmIiWlpakJ2djeeffz7c+ZFKGf3MquuaMIErrBFFuaD+QlNSUvCrX/0KjY2N3quPbDZb\nuHMjtfJ4YDx+XBZ28Wo1oqjXpzN+BoMBSUlJcLvduHbtGq5duxauvEjFDOfOQWhv94lJZjPcvOqI\nKOoFdaRQVVWF9957D01NTbK20tLSkCdF6ubvqiP3uHEAJ1UkinpBFYXi4mI89thjmDt3LtdopsBE\n0f/5BF51RKQKQRWFlpYWPPDAAz4T4RH5o//b3yB0v1TZYICLQ0dEqhDUOYWCggLs3LmzXy9UVVWF\n559/Hj/72c+wefPmHvfbt28fnnjiCZw9e7Zfr0fK8HvDWnY2YLEokA0R9VVQRwrV1dX4/PPP8ec/\n/xkDBw70afvVr37V6/NFUURxcTFWrlwJm82GFStWIC8vT3ZDXHt7Oz7//HOMHTu2D12gqCFJvOqI\nSOWCKgoFBQUoKCi44xepqalBeno60tLSAHROkXHgwAFZUSgtLcWjjz6KLVu23PFrkXL0ly9D1/1i\nBJ0O7okTlUmIiPosqKIwd+7cfr2I3W73ua/BZrOhurraZ59z586hoaEBubm5AYtCWVmZdyK+tWvX\nIjU11afdYDDIYmqmpv4Ie/ZA122YSBo/HrZua3GoqU/B0Fp/AO31SWv9AcLXp6CKwo4dO/zGjUYj\nbDYbxo4d2681nEVRxEcffYSlS5f2um9RUZHPJHwNDQ0+7ampqbKYmqmmP5KEpPJy6LpNbdE+ahQ6\n+B6pjtb6pLX+AP3r09ChQ3tsC6oo7N69G2fOnEFycjJsNhsaGxtx8+ZNZGVlob6+HgDw0ksvISsr\ny+/zrVYrGhsbvduNjY2wWq3ebYfDgUuXLnnPTzQ1NeH1118P+DMpuuhqa6Hr/gsqCHBx6IhIVYIq\nChkZGZg2bRoeeughb2zbtm24cuUK/v3f/x2ffvopfv/732PNmjV+n5+VlYXa2lrU19fDarWivLwc\ny5Yt87bHx8ejuLjYu7169Wo8/fTTLAgq4u/eBPeoUZCSkhTIhojuVFCXpO7duxfz58/3ic2bNw97\n9uyBIAj47ne/i8uXL/f4fL1ejyVLlmDNmjV44YUXMHPmTGRmZqK0tBQVFRX96wFFBb83rPEogUh1\ngjpSSE5OxsGDB3Hvvfd6Y5WVlRgwYAAAwOVywdDL7Je5ubnIzc31iS1atMjvvqtXrw4mLYoSuuvX\noa+rk8V5FzOR+gRVFH70ox9h/fr1GD58uPecwsWLF/Hiiy8C6LyPofuRBMUOfzeseTIzIXW7p4WI\nol9QRWHy5Ml45513UFVVBbvdjilTpiA3NxdJt8eLJ0+ejMmTJ4c1UYpefu9i5lECkSoFveLJgAED\nMGHCBNjtdlitVm9BoNims9uhv3JFFmdRIFKnoIrCjRs38Pbbb6O6uhqJiYlobm72rrzW9dJSij1+\nh46GDoXIRZiIVCmoq48++OADjBgxAr///e/x/vvvo6SkBCNHjsQHH3wQ7vwoynGabCJtCaoonD59\nGs888wwst6cwsFgs+OEPf4gzZ86ENTmKbkJTE/QXLsjiLApE6hVUUUhISJDdh3D16lXEx8eHJSlS\nB39HCWJaGsTBgxXIhohCIahzCt/97nfxyiuvoKCgAIMGDcL169exa9euHu8zoNjgryh08CiBSNWC\nKgpFRUVIT0/Hnj17cPHiRaSkpGDZsmWYxA+AmCU0N8Nw/rwszrUTiNQt6EtSc3JykNPlD14URZSW\nlvJoIUYZjx4FJMknJqamQhwyRKGMiCgUgjqn4I/H48Gnn34aylxIRYxHjshirrvvBriON5Gq3XFR\noNgl3Lrld+iI5xOI1I9Fgfqsx6GjAAt3EJE6BDyncMzP1SXfcrvdIU+G1MHvXEeTJ3PoiEgDAhaF\n9957L+CTtbbmKfVOuHmTQ0dEGhawKGzcuDFSeZBK+B06GjSIVx0RaQTPKVCfmPxcddTBoSMizWBR\noKAJTU3Q/+1vsrjr7rsjnwwRhQWLAgXN31GCmJYGMS1NgWyIKBxYFChoxkOHZLGOSZM4dESkISwK\nFBRdfb3/FdamTFEgGyIKFxYFCoqxqkoW82RkQBw0SIFsiChcWBSod5IEk5+i0HHPPQokQ0ThxKJA\nvdJfvgxdQ4NvUBA672ImIk1hUaBe+TvB7M7KgpScrEA2RBROLAoUmCjCdPiwLMwTzETaxKJAARnO\nnoXQ3Owb1Ou5whqRRrEoUECmigpZzDVuHKS4OAWyIaJwY1GgnrW3w+hn+nTX1KkKJENEkcCiQD0y\nHTkCuFw+MSk+Hq7x4xXKiIjCLeDU2aFUVVWFkpISiKKIwsJCLFiwwKd969at2L59O/R6PQYMGIB/\n/ud/xiDeGKUof0NHHVOmAIaI/doQUYRF5EhBFEUUFxfjl7/8Jd566y3s3bsXly9f9tln5MiRWLt2\nLd58803MmDEDf/jDHyKRGvVAV18P/YULsnhHXp4C2RBRpESkKNTU1CA9PR1paWkwGAyYNWsWDhw4\n4LNPTk4OzGYzAGDs2LGw2+2RSI16YDp4UBbzDBnCdZiJNC4i4wB2ux02m827bbPZUF1d3eP+O3bs\nwD09TKFQVlaGsrIyAMDatWtlS4IaDAZNLROqSH9EEfqTJwGLxTdcWAhLCIb0+B5FP631SWv9AcLX\np6gbHN69ezfOnTuH1atX+20vKipCUVGRd7uh2/QLqampspiaKdEfw6lTSKir8w3qdLg1ejSkEOTC\n9yj6aa1PWusP0L8+DQ1wxB+R4SOr1YrGxkbvdmNjI6xWq2y/I0eO4P/+7//w0ksvwWg0RiI18sO0\nb58s5ho/HlJiogLZEFEkRaQoZGVloba2FvX19XC73SgvL0detxOW58+fxwcffICXXnoJyZxTRzGC\n3Q7jyZOyeMe99yqQDRFFWkSGj/R6PZYsWYI1a9ZAFEXk5+cjMzMTpaWlyMrKQl5eHv7whz/A4XBg\n/fr1ADoPjX7xi19EIj3qwvzNN4Ak+cTElBS477pLoYyIKJIidk4hNzcXubm5PrFFixZ5H7/88suR\nSoV64nbD9M03snDH9OmAjvc5EsUC/qWTl/HYMQgtLb5BvR4d06YpkxARRRyLAnn5PcE8aRJPMBPF\nEBYFAgDoamthOHdOFnfOnKlANkSkFBYFAgCYv/pKFvOkp8MzcmTkkyEixbAoEISmJpj8LLnZMWMG\nIAgKZERESmFRoM6jBI/HJybFx6OD6yYQxRwWhRgntLf7vQzVOWsWcHuCQiKKHSwKMc5UXg7B6fQN\nGo3omD1bmYSISFEsCrHM5YJ5715Z2DltGqSEBAUSIiKlsSjEMNM338hvVtPp4LzvPmUSIiLFsSjE\nKocDlu3bZWHX5MmQ/MxgS0SxgUUhRpm/+kp+lADAcf/9CmRDRNGCRSEGCS0tMO/eLYu7pkzhcptE\nMY5FIQaZt2+XX3Gk18Px4IPKJEREUYNFIcbo7HaY/Ux855wxAyLPJRDFPBaFWCJJiNu8WX73stkM\nZ2GhQkkRUTRhUYghxqoqGE6dksWdc+ZwemwiAsCiEDOE1lbEbdkii4vJybwvgYi8WBRihOUvf4HQ\n2iqLty9YAFgsCmRERNGIRSEGGI4fh6myUhZ33X033BMnKpAREUUrFgWN0zU0IL60VBaX4uLQ/uij\nCmRERNGMRUHLHA4kfPghBIdD1tT+yCOQkpIUSIqIohmLglZJEuI/+QS6a9dkTe677oIrL0+BpIgo\n2rEoaJEkwfKXv8B49KisSbRa0fbkk1xmk4j8MiidAIXY7RvUTF9/LW8ymdC6eDGkuDgFEiMiNWBR\n0BJRRNymTTAdOOC3uX3hQohDhkQ4KSJSExYFjRDsdsSXlsJw/rzfdkdREVyTJ0c4KyJSGxYFtZMk\nGI8cQdymTX6vMgI6C4LzgQcinBgRqRGLglpJEgw1NbB88QX0Fy70uJtj/nw4CwoimBgRqRmLgsoI\nN2/CePQoTFVV0F+8GGBHAe0PP4yOOXMilxwRqR6LQrSSpM65ilpaYDx5EoYrV6C/eBH6K1cASQr4\nVHHgQLT94AfwjB4doWSJSCsiVhSqqqpQUlICURRRWFiIBQsW+LS7XC5s2LAB586dQ1JSEpYvX47B\ngweHNAfBbofFzzKUIdP1w7r7B3fXbVHs/JIkCB4P4HZDcLsBlwtCezsEhwO6lhbA5YLeYkF8D+cK\n/HFNmoT2xx6DFB/fz84QUUZiloQAAAi8SURBVCyKSFEQRRHFxcVYuXIlbDYbVqxYgby8PGRkZHj3\n2bFjBxISEvDOO+9g7969+OMf/4gXXnghpHnoWlthKi8P6c+MFp4hQ+B44IHOCe54YxoR3aGI3NFc\nU1OD9PR0pKWlwWAwYNasWTjQ7Vr6iooKzJ07FwAwY8YMHDt2DFIvwyTUWQzann4aLcuXw52Tw4JA\nRP0SkSMFu90Om83m3bbZbKiuru5xH71ej/j4eDQ3N2PAgAE++5WVlaGsrAwAsHbtWqSmpvq0GwwG\nWcyrpQV6la0doNPpYOmec1oaxClTIE2dCuPQobCorBAEfI9USGv9AbTXJ631Bwhfn1R3ormoqAhF\nRUXe7YaGBp/21NRUWexb+hs3kNiH8XmlSWYzLOnpaDWZ4BkyBO7MTHgyMiAlJ//9iKCxUdkk70Cg\n90iNtNYfQHt90lp/gP71aejQoT22RaQoWK1WNHb58GpsbITVavW7j81mg8fjQVtbG5JCPLWzOHBg\n50pjkdL1P/hujyWdDrj9JRmNgF4PGI0Q4+KAuDhIcXGQLBaYBg1Cq8Z+mYkoekWkKGRlZaG2thb1\n9fWwWq0oLy/HsmXLfPaZOnUqdu3ahezsbOzbtw8TJ06EEOJhESkpCR2zZoX0ZxIRaUlEioJer8eS\nJUuwZs0aiKKI/Px8ZGZmorS0FFlZWcjLy0NBQQE2bNiAn/3sZ0hMTMTy5csjkRoREXURsXMKubm5\nyM3N9YktWrTI+9hkMuHFF1+MVDpEROQHF9khIiIvFgUiIvJiUSAiIi8WBSIi8hIkziVBRES3ae5I\n4V//9V+VTiGktNYfQHt90lp/AO31SWv9AcLXJ80VBSIiunMsCkRE5KVfvXr1aqWTCLXRGltxTGv9\nAbTXJ631B9Ben7TWHyA8feKJZiIi8uLwERERebEoEBGRl+oW2flWVVUVSkpKIIoiCgsLsaDbOgku\nlwsbNmzAuXPnkJSUhOXLl2Pw4MEKZdu73vqza9cufPzxx951KObPn4/CwkIlUg3Ku+++i8rKSiQn\nJ2PdunWydkmSUFJSgkOHDsFsNmPp0qVRPebbW3+OHz+O119/3fs7Nn36dDz++OORTrNPGhoasHHj\nRjQ1NUEQBBQVFeGhhx7y2UdN71Mw/VHb+9TR0YFVq1bB7XbD4/FgxowZeOKJJ3z2CflnnaRCHo9H\neu6556S6ujrJ5XJJP//5z6VLly757LNt2zbpt7/9rSRJkrRnzx5p/fr1SqQalGD6s3PnTul3v/ud\nQhn23fHjx6WzZ89KL774ot/2gwcPSmvWrJFEUZROnz4trVixIsIZ9k1v/Tl27Jj06quvRjir/rHb\n7dLZs2clSZKktrY2admyZbLfOzW9T8H0R23vkyiKUnt7uyRJkuRyuaQVK1ZIp0+f9tkn1J91qhw+\nqqmpQXp6OtLS0mAwGDBr1iwcOHDAZ5+KigrMnTsXADBjxgwcO3YMUpSeUw+mP2ozYcIEJCYm9the\nUVGBOXPmQBAEZGdno7W1FTdu3Ihghn3TW3/UKCUlxftff1xcHIYNGwa73e6zj5rep2D6ozaCIHjX\naPd4PPB4PLLFx0L9WafK4SO73Q6bzebdttlsqK6u7nEfvV6P+Ph4NDc3Y8CAARHNNRjB9AcA9u/f\nj5MnT2LIkCFYvHixqhcit9vtPvnbbDbY7XakpKQomFX/nDlzBv/yL/+ClJQUPP3008jMzFQ6paDV\n19fj/PnzGDNmjE9cre9TT/0B1Pc+iaKIX/ziF6irq8ODDz6IsWPH+rSH+rNOlUUhFk2dOhWzZ8+G\n0WjEX//6V2zcuBGrVq1SOi26bdSoUXj33XdhsVhQWVmJN954A7/+9a+VTisoDocD69atw7PPPov4\n+Hil0+m3QP1R4/uk0+nwxhtvoLW1FW+++SYuXryI4cOHh+/1wvaTw8hqtaKxsdG73djY6D0B628f\nj8eDtrY2JCUlRTTPYAXTn6SkJBiNRgBAYWEhzp07F9EcQ81qtaKhocG77a/PahIfH+89zM/NzYXH\n48GtW7cUzqp3brcb69atw3333Yfp06fL2tX2PvXWH7W+TwCQkJCAiRMnoqqqyice6s86VRaFrKws\n1NbWor6+Hm63G+Xl5cjLy/PZZ+rUqdi1axcAYN++fZg4caJsLC5aBNOfruO4FRUVyMjIiHSaIZWX\nl4fdu3dDkiScOXMG8fHxUT8kEUhTU5N3HLempgaiKEbtPyHfkiQJv/nNbzBs2DA88sgjfvdR0/sU\nTH/U9j7dunULra2tADqvRDpy5AiGDRvms0+oP+tUe0dzZWUlPvzwQ4iiiPz8fHz/+99HaWkpsrKy\nkJeXh46ODmzYsAHnz59HYmIili9fjrS0NKXT7lFv/fnv//5vVFRUQK/XIzExET/+8Y9lvxzR5O23\n38aJEyfQ3NyM5ORkPPHEE3C73QCAefPmQZIkFBcX4/DhwzCZTFi6dCmysrIUzrpnvfVn27Zt+OKL\nL6DX62EymfDMM89g3LhxCmcd2KlTp/Bv//ZvGD58uPdD5Mknn/QeGajtfQqmP2p7ny5cuICNGzdC\nFEVIkoSZM2fi8ccfD+tnnWqLAhERhZ4qh4+IiCg8WBSIiMiLRYGIiLxYFIiIyItFgYiIvFgUiIjI\ni0WBiIi8WBSIiMiLRYEoROrq6vCjH/3IOy+V3W7HP/3TP+H48eMKZ0YUPBYFohBJT0/HP/7jP+Kd\nd96B0+nEe++9h/vvvx8TJ05UOjWioHGaC6IQe+2111BfXw9BEPDqq696Z7clUgMeKRCFWGFhIS5d\nuoT58+ezIJDqsCgQhZDD4cCHH36IgoICfPLJJ2hpaVE6JaI+YVEgCqGSkhKMHj0aP/3pT5Gbm4v3\n339f6ZSI+oRFgShEDhw4gKqqKvzkJz8BACxevBjnz5/HV199pXBmRMHjiWYiIvLikQIREXmxKBAR\nkReLAhERebEoEBGRF4sCERF5sSgQEZEXiwIREXmxKBARkdf/B+pLjfR55C8mAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4slV2qd7zqxh",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        },
        "outputId": "7c4e6ee8-07db-4422-d25b-3cf08508983f"
      },
      "source": [
        "x = np.linspace(-5, 5, 100)\n",
        "plt.axes(xlabel=\"x\", ylabel=\"Normal PDF\")\n",
        "plt.plot(x, normal_pdf(x, mean, sigma),'r-', lw=5, \n",
        "         alpha=0.6)\n",
        "plt.show()"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEJCAYAAABlmAtYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1wU570/8M9c9gIsorvooogSVjRe\nqhbXG5pGhJPanNZwElNT2/yamJOk0Vw0rx6riU2b5JiQGs1pNclpciixsT2H3GxzVYPWmIhGNMH7\nBQQ1CEJkvYDsws7O/P5At+zOAAvszuzl+369eOl+d4b9Pqz43XmeZ56HkSRJAiGEENIFVusECCGE\nhD8qFoQQQrpFxYIQQki3qFgQQgjpFhULQggh3aJiQQghpFu81gmESm1trdYp9FhycjIuXLigdRqq\nojbHBmpzZBgyZEinz9GVBSGEkG5RsSCEENItKhaEEEK6RcWCEEJIt1QrFuXl5XjsscfwyCOP4G9/\n+1unx+3Zswc//vGPcerUKW9s06ZNeOSRR/DYY4+hvLxcjXQJIYR0oEqxEEURhYWFeOKJJ/DSSy9h\n165dqKmpkR3ndDrxySefIDMz0xurqalBaWkp1q5diyeffBKFhYUQRVGNtAkJHpcL3KlT4CsrgePH\nwdXUAPTvmEQQVabOVlZWIiUlBVarFQCQnZ2NsrIyDB061Oe44uJi3HbbbXj//fe9sbKyMmRnZ0On\n02HQoEFISUlBZWUlRo4cqUbqhPRNayuMW7bAsGcPIAgAAM5ohMnlgmQywfmv/wp3VhbAMBonSkjX\nVCkWDocDFovF+9hisaCiosLnmKqqKly4cAFZWVk+xcLhcPhcaZjNZjgcDtlrlJSUoKSkBABQUFCA\n5OTkYDcj5Hiej8i8+yKq23z8OLi//hW4cAHg+fYvACzLwmg0AoKAuL//HVJlJcQFC4ABAzROOHSi\n+n3uRLS1OSxuyhNFEX/+85+xaNGiXn+PvLw85OXleR9H2s0wQGTexNNX0dpmQ0kJjFu3wq3wnNFo\nhMvl+mdg/35IR47g6j33wJORoVqOaorW97krkdhmzW/KM5vNaGxs9D5ubGyE2Wz2Pna5XPjmm2/w\n9NNPY/HixaioqMDvfvc7nDp1Snauw+HwOZeQcKMrL4dx69YencO4XEj485/BKFw1ExIOVCkWNpsN\ndXV1aGhogCAIKC0thd1u9z4fHx+PwsJCvPzyy3j55ZeRmZmJZcuWwWazwW63o7S0FG63Gw0NDair\nq8OIESPUSJuQHmPPn0fcO+8oP6nTQcjIgJServg009KChI0bAbfS9Qgh2lKlG4rjOCxcuBCrVq2C\nKIrIyclBWloaiouLvQWhM2lpaZg+fToef/xxsCyL++67DyxLt4eQMOR0IuHNN8G0tcmeEjIy4Jw3\nD2JyMuKSk3H1s88Qv2kTmCtXfI7jamoQ9/77cN5xh1pZExIQJlr34KaFBCND1LRZkhD/5pvQHT4s\ne8o9YQJaFizwzni63mamuRmmP/wB7KVLsnNafvxjuLv4EBVpouZ97oFIbLPmYxaERDv+2DHFQiFa\nrWiZN09xaqxkMqHl7ru9s6Q6ivvgAzBOZ0hyJaQ3qFgQ0leiCOPmzbKwZDTi6v/7f4DB0OmpnrQ0\nOOfOlcUZpxOGzz4LapqE9AUVC0L6SFdeDu78eVncefvtEAcO7Pb8tqlT4Z44URY3fP65bEyDEK1Q\nsSCkLwQBxi1bZGFPejrcEyYE9j0YBs5bb5V3R7ndMG7bFoQkCek7KhaE9IH+yy/BXrwoizvnzOnR\nEh5S//5onTFD+ftH2CApiU5ULAjprdZWxU/+wo039upO7NacHEhGo29QFHt8gx8hoUDFgpBe0u/f\nD6a5WRZ3zpnTq+8nxcej9eabZXHdgQNg6c5uojEqFoT0hiTBsHu3LOz+7nchdjFXvTutM2dCSkyU\nvZb+yy97/T0JCQYqFoT0AlddDba+XhZ35eT07RsbDGidOVMW1u/d613inBAtULEgpBeUriqEjAyI\nKSl9/t5tkyfLZkYxV69Cd/Bgn783Ib1FxYKQHmKamqA7dEgWb5s+PSjfXzKZ0DZ+vCxu2LMnKN+f\nkN6gYkFID+n37pVtiSolJsI9dmzQXqNt2jRZjDt9GmwErnlGogMVC0J6QhQVB5vbpkxRXOOptzzD\nh8MzeLAsTlcXRCtULAjpAf7YMfkqsQyD1ilTgvtCDKPYraX7+mug4y57hKiEigUhPaD/+mtZzD16\nNKQQ7J/d9t3vym7SY1pboTt6NOivRUh3qFgQEiiXC7pjx2ThtqlTQ/N6BoPiAoO6AwdC83qEdEGV\nnfIAoLy8HEVFRRBFEbm5ucjPz/d5fuvWrdiyZQtYloXRaMSDDz6IoUOHoqGhAUuXLvVuypGZmYkH\nHnhArbQJ8dIdPSrb8lSKj4cwcmTIXrNt4kTo/cYpdCdOgGlpgRQfH7LXJcSfKsVCFEUUFhZi5cqV\nsFgsWLFiBex2O4YOHeo9ZubMmbjlllsAAPv27cOGDRvw5JNPAgBSUlKwevVqNVIlpFNKn+jd48cD\nHBey1/TccAPEpCSwly//MyiK0B06FLorGkIUqNINVVlZiZSUFFitVvA8j+zsbJSVlfkcE9/hU5LL\n5QLTgxU7CQk1pqUFuhMnZPG2QJch7/ULM+0FyY+uvDy0r0uIH1WuLBwOBywWi/exxWJBRUWF7LjN\nmzfjo48+giAIeOqpp7zxhoYGLFu2DHFxcbjrrrswevRo2bklJSUoKSkBABQUFCA5OTkELQktnucj\nMu++iJQ2M198AVav9w327w/d5MkA27PPXD1u8+zZ4Pw+XKG2FnE8D/Tv36PX1kqkvM/BFG1tVm3M\nIhBz5szBnDlz8MUXX+Ddd9/Fww8/jAEDBuCVV15BYmIiqqqqsHr1aqxZs8bnSgQA8vLykJeX530c\naRulA5G5wXtfRUqbEz77DLzflNXWzEy4erEabI/bHB+PRJNJtq+F8x//QNtNN/X49bUQKe9zMEVi\nm4d0sQimKt1QZrMZjY2N3seNjY0wm82dHt+xm0qn0yHx2iqcGRkZsFqtqKurC23ChHTAXLkCvqpK\nFleaqRSaBBjFXff0NCuKqEiVYmGz2VBXV4eGhgYIgoDS0lLY7XafYzoWgK+++gqDr929euXKFYjX\nllaor69HXV0drFarGmkTAgDtC/hJkk9MtFjg6TBBI9TaFAoTd/Ys7XNBVKNKNxTHcVi4cCFWrVoF\nURSRk5ODtLQ0FBcXw2azwW63Y/PmzTh06BA4joPJZMLixYsBAEePHsVbb70FjuPAsizuv/9+mEwm\nNdImBAAUb4JzT5zYo21T+0q0WuEZPBic31U1f+RIxHRFkcim2phFVlYWsrKyfGLz58/3/v3ee+9V\nPG/atGmYprCoGiFqYJxOxS6otnHjVM/F/Z3vyIqF7tgxKhZEFXQHNyFd4E+ckK0wK/bv36fd8HrL\nPWaMLMZXVYFxOlXPhcQeKhaEdEF35IgsJowZo2oX1HXi4MEQ/degEkXwx4+rnguJPVQsCOmMILRf\nWfgJ5r4VPcIwilcXtLAgUQMVC0I6wVdXg/G7t0IyGiHccINGGV27qvHDnzhB+3OTkKNiQUgneIVP\n7MKoUUHd5KinhBtukC9b7nKBr67WKCMSK6hYEKJEkpSnzCosNaMqnm8vWP5h6ooiIUbFghAFbF0d\n2IsX/YIshBtv1CahDjodt/C7cZCQYKJiQYgCpU2OhBtuCIs9JIQbb5QtXshevAiWlsEhIUTFghAF\nSsuRK32i14IUF6c4yK47eVKDbEisoGJBiB/G6QR35owsHg5dUNcJCmMndL8FCSUqFoT44Ssr5QsH\nms0Qw2hvArfCVq78mTOA31RfQoKFigUhfpQ+oQujRmly13ZnRKsVYlKSb9DjAX/qlDYJkahHxYKQ\njiQJvMIujkqf5DXFMBAUctIp5E5IMFCxIKQDtr4e7KVLfkEWgs2mTUJdULzfgsYtSIhQsSCkA15h\nRpGQng743TUdDoTMTFnXGOtwyLZfJSQYqFgQ0oHSlNlwmgXVkRQXB8/w4bK40uKHhPSVaovclJeX\no6ioCKIoIjc3F/n5+T7Pb926FVu2bAHLsjAajXjwwQcx9Nq2lZs2bcL27dvBsizuvfdeTFRr72MS\nW9raFNdYcmdmapBMYNyZmeBOn/aJ8SdPom3GDG0SIlFLlSsLURRRWFiIJ554Ai+99BJ27dqFmpoa\nn2NmzpyJNWvWYPXq1bjtttuwYcMGAEBNTQ1KS0uxdu1aPPnkkygsLPTuyU1IMPFVVbLVW6XERE02\nOgqU0lUPf+oUrUJLgk6VYlFZWYmUlBRYrVbwPI/s7GyUlZX5HBPfYRkFl8sF5lpfbFlZGbKzs6HT\n6TBo0CCkpKSgsrJSjbRJjFHcu2LkyLCaMuvPk5oqW4KEaWsD73e1QUhfqdIN5XA4YLFYvI8tFgsq\nFKb4bd68GR999BEEQcBTTz3lPTezQzeA2WyGw+GQnVtSUoKSkhIAQEFBAZLD6AaqQPE8H5F590U4\ntZmrrZUNZOunTkVCkPMLdpvZrCww+/b5xPT19ZDCaO/6cHqf1RJtbdZuYX4Fc+bMwZw5c/DFF1/g\n3XffxcMPPxzwuXl5ecjLy/M+vhCBM0KSk5MjMu++CJc2M5cvo5/CeMWV5GRIQc4v2G3WDR6MeL87\ntz3796M5jMYtwuV9VlMktnlIF12uqnRDmc1mNDY2eh83NjbCbDZ3enzHbir/cx0OR5fnEtIbSnc+\newYPhmQyaZBNzwgjRshi3LlzYJxODbIh0UqVYmGz2VBXV4eGhgYIgoDS0lLY7XafY+o6LK/81Vdf\nYfDgwQAAu92O0tJSuN1uNDQ0oK6uDiMUfjkI6Qulu7aFMJ4F1ZE0YIB83SpJAkdLf5AgUqUbiuM4\nLFy4EKtWrYIoisjJyUFaWhqKi4ths9lgt9uxefNmHDp0CBzHwWQyYfHixQCAtLQ0TJ8+HY8//jhY\nlsV9990HlqXbQ0gQSVL74oF+lD6xhythxAjo/bo8dBUVEMaN0ygjEm1UG7PIyspCVlaWT2z+/Pne\nv997772dnnv77bfj9ttvD1luJLaxFy6AvXzZN8hxintGhCshMxP6PXt8YkoFkJDeoo/oJOYpdkEN\nHw4YDBpk0zuCzSZf+uPbb8H4r3NFSC9RsSAxT7ELKgwXDuyKFB8PT2qqLK5UCAnpDSoWJLaJonKx\niJDB7Y6UxlhofwsSLFQsSEzjzp0D43ePgmQwwJOWplFGvadU4HQVFbJd/wjpDSoWJKYpjlfYbADH\naZBN3wjp6QDvO2eFaWoCW1+vTUIkqlCxIDEtGsYrvHS69oF5PzRuQYKBigWJXW43+DNnZOFIur/C\nn1JXFI1bkGCgYkFiFnf2LOB2+8QkkwliSopGGfWd0lURX1UF0LL+pI+oWJCYpfSJW+l+hUjiGToU\nkt/9IYzL1b6iLiF9QMWCxKxOi0Uk4zh4FO48p64o0lddFotLdPcniVZtbeDPnpWFI75YoJOuKCoW\npI+6LBaPPfaYz+MXX3wxpMkQohb+zBnA4/GJiUlJ8tVbI5DikuXV1bL2EtITXRYLye9mniNHjoQ0\nGULU0umU2Qger7jOM3gwpLg4nxjT2grOb997Qnqiy2LBRMEvDiFKonK84jqWhZCRIQvTKrSkL7pc\notzj8eDw4cPex6Io+jwGgHG0Xj6JNC6X4qfsqCkWaG+Lzq8ngD91Cq25uRplRCJdl8UiKSkJr776\nqvexyWTyecwwDNavXx+67AgJAb66WnbfgWg2Q4qi7XoVFxU8c6b9vhKdToOMSKTrsli8/PLLQXuh\n8vJyFBUVQRRF5ObmIj8/3+f5Dz/8ENu2bQPHcejXrx8eeughDBw4EED7JknDhg0D0L4J+q9+9aug\n5UViT1Qt8dEJ0WqFZDKBaW7+Z9DtBnf2LDxR1laijoB3yqutrUVzczNMJhOGDBnSoxcRRRGFhYVY\nuXIlLBYLVqxYAbvdjqFDh3qPSU9PR0FBAQwGA7Zu3YqNGzdi6dKlAAC9Xo/Vq1f36DUJ6UxUj1dc\nxzDtXVEHDviE+VOnqFiQXum2WHz22WfYuHEjrly54o0lJSVhwYIFmDVrVkAvUllZiZSUFFitVgBA\ndnY2ysrKfIpFx7GPzMxMfP7554G2gZCAMS0t4OrqZPGoKxaAcrGoqkKrRvmQyNZlsTh48CAKCwtx\n5513YurUqRgwYAAcDge+/PJLFBUVwWw2Y/z48d2+iMPhgMVi8T62WCyo6GIlzO3bt2PixInex263\nG8uXLwfHcbjtttswZcoU2TklJSUoKSkBABQUFCA5AufL8zwfkXn3hdptZg4cAOu/XarVCp2KxUK1\nNk+eDO7jj31j9fWI69cP0OtD//od0L/tyNdlsfjkk09w11134dZbb/XGrFYr5s6dC71ej48//jig\nYtETO3fuRFVVFX772996Y6+88grMZjPq6+vxzDPPYNiwYUjxW+wtLy8PeXl53scXLlwIal5qSE5O\njsi8+0LtNhv37YPBb7OjtsGD4VQxB9XazLLop9eD6dArAABX9+9XfSdA+rcdGboaYujyPotTp05h\nxowZis9Nnz4dpwJcQsBsNqOxsdH7uLGxEWaFmScHDx7Epk2bsGzZMug6zNi4fqzVasWYMWNw+vTp\ngF6XEH98VZUsFo1dUAC84xb+aOkP0htdFovW1lYkJSUpPpeUlITW1sB6P202G+rq6tDQ0ABBEFBa\nWgq73e5zTHV1NV5//XUsW7bM5zWbm5vhvraM9JUrV3DixAmfsQ5CAsU0NyuPVyjcwBYtqFiQYOl2\ngFuSJNmyH9cFeoc3x3FYuHAhVq1aBVEUkZOTg7S0NBQXF8Nms8Fut2Pjxo1wuVxYu3YtgH9OkT13\n7hxee+01sCwLURSRn59PxYL0Cl9dLYuJViukxEQNslGHUrHgvvkGaG0F/MduCOlCl8XC5XLhrrvu\nCsoLZWVlISsryyc2f/58799//etfK543atQorFmzJig5kNgWE1Nm/YhmM8T+/cF2XEFaFMGfPg1h\n1CjtEiMRp8tiQXdnk2gSi8Xi+riFfv9+nzB/6hQVC9IjXRaLgQMH4ty5c6ipqcHw4cNlM5AIiRRM\nUxPY+npZPJrHK64TMjIUiwUhPdFlsdixYwf++Mc/IiEhAS0tLXjkkUcwffp0tXIjJGiU/nP0DB4M\nKSFBg2zUpThuUVMDOJ2A31LmhHSmy2Lx97//HY8//jgmT56MvXv34t1336ViQSJSTHZBXSOZzRDN\nZrAOR4egBL66GsKYMdolRiJKl1NnHQ4HJk+eDACYPHlyxN1gQsh1sVwsAJpCS/quy2LREcMwEP2W\ndSYkEjCXLoH1/6DDMPDEwHjFdVQsSF91O3X2oYce8j5uaWnxeQzAZ38LQsKR4njFkCGyrUejmeK4\nRV0dmJYWSPHxGmREIk2XxeI3v/mNWnkQEjKKXVAKmwNFMykpCeLAgWC//bZDUAJXVQWBdrskAeiy\nWIyhwS8S6SQp5scrrhNsNug7Fgtcu9+CigUJQMBjFoREIvbiRbAXL/oFWQjp6ZrkoyWlqymdwq6B\nhCihYkGimtIWqp60NMBo1CAbbSndgMjW14NpatIgGxJpqFiQqKa433aMjVdcJ5lM8AweLIvTrCgS\nCCoWJHrReIUMTaElvdXpAPe6desCWoL84YcfDmpChAQL29Ag72LheQjDh2uTUBgQbDYYvvjCJ6Z0\n9UWIv06LBS0aSCKd4lXFsGFAh10YY40nIwNgGKDDHjVsYyOYixchDRigYWYk3HVaLO6880418yAk\n6Gi8Qk6Ki4MnNbV9IcEO+FOn4PbbvZKQjrrdKe86QRBQW1uLK36bv48LcI52eXk5ioqKIIoicnNz\nkZ+f7/P8hx9+iG3btoHjOPTr1w8PPfQQBg4cCKB99dv33nsPAHD77bdj1qxZgaZNYpUoKu+3HePF\nAmj/GVCxID0VULE4fvw41q5dC7fbDafTibi4OLhcLlgsloA2SBJFEYWFhVi5ciUsFgtWrFgBu93u\nsz1qeno6CgoKYDAYsHXrVmzcuBFLly5Fc3Mz3nnnHRQUFAAAli9fDrvdDpPJ1Msmk1jA1daCaWnx\niUl6PTy0JW/7uMWOHT4xvrKyvWsqwK2SSewJaDbUhg0bMHfuXBQVFSEuLg5FRUW44447cMsttwT0\nIpWVlUhJSYHVagXP88jOzkZZWZnPMePGjYPh2p7AmZmZcFxbTrm8vBzjx4+HyWSCyWTC+PHjUV5e\n3pM2khikeH9FRgbAB3wxHbWE9HSA43xi7OXLvkuBEOInoN+c2tpa3HrrrT6x/Px8LF68GHPnzu32\nfIfDAYvF4n1ssVhQUVHR6fHbt2/HxIkTFc81m83eQtJRSUkJSkpKAAAFBQVITk7uNq9ww/N8RObd\nF6FqM1tbC8bvxjtx0iTEhcHPNxzeZ3bMGDB+v4P6hgZIIVriJxzarLZoa3NAxSI+Ph5OpxMJCQno\n378/ampqYDKZ4HK5gp7Qzp07UVVVhd/+9rc9Oi8vLw95eXnex5G490ZycnJE5t0XIWmz242kI0cA\nt9sn3DRwIMQw+PmGw/tsGDIExkOHfGLuffvQEqJ1osKhzWqLxDYPGTKk0+cC6oaaOnUqvv76awBA\nTk4Onn76aSxfvhzTpk0LKAGz2YzGxkbv48bGRpjNZtlxBw8exKZNm7Bs2TLork1v9D/X4XAonkvI\ndfyZM7JCIZlMEBXuXo5VSgP9fFUVQHvWkE4EdGVxzz33eP8+d+5cZGZmwuVyYcKECQG9iM1mQ11d\nHRoaGmA2m1FaWopHH33U55jq6mq8/vrreOKJJ5CUlOSNT5w4Ef/7v/+L5uZmAMCBAwewYMGCgF6X\nxKZOp8zS4K2XJy0NktEIpkPvAONygaupgWfYMA0zI+GqV6N9o0eP7tHxHMdh4cKFWLVqFURRRE5O\nDtLS0lBcXAybzQa73Y6NGzfC5XJh7dq1ANov4X71q1/BZDLhjjvuwIoVKwAA8+bNo5lQpEu8wniY\nkJmpQSZhjGUhZGRAd/SoT5ivqKBiQRQFVCwuXLiAt99+G6dPn5aNU/z+978P6IWysrKQlZXlE5s/\nf77377/+9a87PXf27NmYPXt2QK9DYhvjdMruIQAAdwyvB9UZITNTsVi05uZqlBEJZwEVi7Vr12LI\nkCH48Y9/DL1eH+qcCOk17tQpn6UsAEC0WCDROJeM4rjF2bNAWxtAv+fET0DF4ty5c/jP//xPsCwt\nUkvCm9JmPtQFpUwcNAhSv35gOq7KIAjgT5+GMHKkdomRsBTQ//6TJk3CUb/LVULCEY1X9ADDwK3w\ns1H6GRIS0JXFwoULsXLlSlitVp+ZSgCwaNGikCRGSE8xFy/K70JmmJjev6I7wogR0O/f7xPjT54E\n/vVfNcqIhKuAisUrr7wClmWRmppKYxYkbOlOnpTFPKmpkOLjNcgmMiiNW3B1dWCamiAlJmqQEQlX\nARWLw4cP449//CPi4uJCnQ8hvcYrFAvqe++alJQET0oKuPPnfeL8yZNwT5qkUVYkHAU0ZjF8+HA0\n0abuJJyJomJfu5uKRbeUCqrSVRqJbQFdWYwdOxarVq3CrFmzZGMWdP8DCQfc2bM+dyMDgGQ0whPD\nW6gGShg1CoadO31i/MmTtGQ58RFQsThx4gTMZjMOHjwoe46KBQkHil1QGRmypbiJnJCe3r7VbIf1\ntJirV8GdO0f7fxCvbouFJEn4xS9+geTkZHD0i0fClFK3iTBqlAaZRCCdDkJGBvgTJ3zC/MmTVCyI\nV7djFgzD4Je//CUYuhwlYYppaQH3zTeyOBWLwCmN7ShdrZHYFdAAd3p6Ourq6kKdCyG9wldUyJf4\nSE6GSEt8BExpkJs/fRoIwZ41JDIFPMD93HPP4eabb5bt/ERjFkRrNGW278RBgyD27w/20qUOQRH8\nqVMQxo7VLjESNgIe4B40aBCOHTsme46KBdGUJCkWC5oy20MMA2HkSOj37vUJ606epGJBAARYLH7z\nm9+EOg9CeoU9fx7s5cu+QY6jJT56QRg1SlYs+OPHaQotAdCDzY+am5uxf/9+77amkyZN6tEmROXl\n5SgqKoIoisjNzUV+fr7P80ePHsWGDRtw5swZLFmyxGfL1vnz52PYtQ1Zrm+KRAgA6I4fl8WE4cMB\ng0GDbCKbMGIEwLI+W6uyFy+CbWiAaLVqmBkJBwEVi5MnT+L5559HamoqkpOT8dVXX+GNN97AihUr\nMDKAy31RFFFYWIiVK1fCYrFgxYoVsNvtGNphWl5ycjIWLVqEDz74QHa+Xq/H6tWre9AsEit0Cl2j\nwpgxGmQS+aS4OAjp6e17cXegO3YMrVQsYl5AxeKNN97Av//7v2PGjBneWGlpKYqKivD88893e35l\nZSVSUlJgvfYPLjs7G2VlZT7FYtCgQQBAU3RJwJiWFnBnzsji7htv1CCb6CDceKOsWPDHj6N11ixt\nEiJhI6BiUVdXh+nTp/vEpk2bhtdffz2gF3E4HLBYLN7HFosFFT1YM9/tdmP58uXgOA633XYbpkyZ\nIjumpKQEJSUlAICCggLZrK1IwPN8RObdF31pM7N3L1j/7qaBA6EbPTqs+9jD+n2eMQPc9u2+sfPn\n2xcRTUjo9bcN6zaHSLS1OaBikZKSgtLSUsycOdMb2717t/dKIdReeeUVmM1m1NfX45lnnsGwYcOQ\nkpLic0xeXh7y8vK8jy9cuKBKbsGUnJwckXn3RV/aHL9nD3R+9wG0Dh8OV2NjMFILmbB+n3keifHx\nYB0On3DL7t1wT5zY628b1m0OkUhs85AhQzp9LqBicc8996CgoACffPIJkpOT8e2336Kurg7Lly8P\nKAGz2YzGDr/AjY2NMPfghqnrx1qtVowZMwanT5+WFQsSY0RRtjwFAAijR2uQTBRhGLhHj4Zh1y6f\nMH/8eJ+KBYl8Ad3BPWrUKKxbtw7f//73kZGRgTlz5mDdunUYFeByCjabDXV1dWhoaIAgCCgtLYXd\nbg/o3ObmZrivLXB25coVnDhxwmesg8Qm7vRpME6nT0wyGCDccINGGUUPQWHMR3f8uM8sKRJ7Ap46\nazKZ8L3vfa9XL8JxHBYuXKBz8QsAABjgSURBVIhVq1ZBFEXk5OQgLS0NxcXFsNlssNvtqKysxIsv\nvoirV69i//79eOutt7B27VqcO3cOr732GliWhSiKyM/Pp2JBlGdBZWYCfMD/pEknhIwMSAYDmNZW\nb+z6ZAIPFeOYxUiS36I6HTz99NNdn8wweOqpp4KeVDDU1tZqnUKPRWIfZ1/1ts2Ja9aAra/3ibXc\neSfckycHK7WQiYT3OX7DBuiOHPGJtebkwPWDH/Tq+0VCm4MtEtvc6zGLm266STHucDjwySefoLXD\nJw9C1MI6HLJCASh3n5DeEcaMkRUL/tgxoJfFgkS+LouF/7pPTU1N2LRpE7Zt24bs7GzMmzcvpMkR\nokR3+LAs5klLg5SYqEE20ck9ahTi/GLc+fNgv/0W4sCBmuREtBVQB29LSwvef/99bNmyBVlZWXjh\nhRdoNhLRjFKxcNNid0El9esHz7Bh4M6e9YnrjhyhG/RiVJfFoq2tDR999BE+/PBDjBkzBs888wzS\n0tLUyo0QGebKFXCnT8viVCyCzz1unLxYHDpExSJGdVksFi9eDFEUMXfuXNhsNly+fBmX/Vb4HDdu\nXEgTJKQj/3504NpeDLR2UdC5x46F8eOPfWLcN9+AuXQJUv/+GmVFtNJlsdDr9QCArVu3Kj7PMAzW\nr18f/KwI6YRiFxR9YAkJceBAeFJSwJ0/7xPXHTmCtg7rxJHY0GWxePnll9XKg5BuMS0t4E+dksWp\nWISOe9w4ebE4dIiKRQwK6A5uQsIBf+yY7C5iccAAeFJTNcoo+ikVYr66GkxzswbZEC1RsSARQ3fo\nkCzmHjs2rFeYjXTi4MEQ/ddxkyTFO+hJdKNiQSJDayt0CsvaUxdUiDEM3N/5jiysVLhJdKNiQSKC\n7tgx4NqCktdJJhM86enaJBRDFLuiKirAtLRokA3RChULEhF05eWymHvs2PY9o0lIeYYNg5iU5Bf0\n0NVFjKHfNBL2mJaW9iWy/bRNmKBBNjGIYeAeP14W1n39tQbJEK1QsSBhT3fokHwWVFISPBkZGmUU\ne9zf/a4sxldXg7l0SYNsiBaoWJCwp/QJ1j1+PHVBqciTmgrRfz9pSYL+4EFtEiKqo982EtaYS5fA\nV1fL4kqfdEkIMQzaFLZVVRpLItFJtW3FysvLUVRUBFEUkZubi/z8fJ/njx49ig0bNuDMmTNYsmQJ\npk2b5n1ux44deO+99wAAt99+O2bRQmYxQ3/gAOC3P5c4cCDdiKcB93e/C2NJiU+Mq6mhZctjhCpX\nFqIoorCwEE888QReeukl7Nq1CzU1NT7HJCcnY9GiRZg5c6ZPvLm5Ge+88w6ee+45PPfcc3jnnXfQ\nTHePxgylLqi2iRPpRjwNiAMHwqOwpTENdMcGVYpFZWUlUlJSYLVawfM8srOzUVZW5nPMoEGDMHz4\ncDB+/wmUl5dj/PjxMJlMMJlMGD9+PMrp0jcmsPX14BS2x3UrdIcQdSh1RenLy2VXfyT6qNIN5XA4\nYLFYvI8tFgsqFO7GDeRcs9kMh8MhO66kpAQl1y6RCwoKkOw/GBcBeJ6PyLz7oqs2M59/DtZo9IlJ\nw4dDN3q0GqmFTES/z7Nng9u2zbc4NDfDcPkyMGJEp6dFdJt7KdrarNqYRajl5eUhLy/P+zjSNkoH\nInOD977qtM2CgH7/+AcYl8sn7MzMRFuE/4wi/X1OSE0FX1npE2vbsgXOLva4iPQ290YktnnIkCGd\nPqdKN5TZbEZjY6P3cWNjI8z+i5MFeK7D4Qj4XBK5dMeOyVc25Xm4J03SJiHi1TZ5siymP3gQjNOp\nQTZELaoUC5vNhrq6OjQ0NEAQBJSWlsJutwd07sSJE3HgwAE0NzejubkZBw4cwETqs456+r17ZTH3\n2LGQEhI0yIZ05B43DlJcnF/QTdNoo5wq3VAcx2HhwoVYtWoVRFFETk4O0tLSUFxcDJvNBrvdjsrK\nSrz44ou4evUq9u/fj7feegtr166FyWTCHXfcgRUrVgAA5s2bB5PJpEbaRCPMxYvgT56UxdumTNEg\nGyKj06EtKwuGXbt8wvq9e9E2fbpGSZFQYyQpOqcx1CrMogl3kdjH2VdKbTZ8+imMn37qExMHDEDT\n8uVRMWU2Gt5ntrYWif/1X7J482OPKd4DEw1t7qlIbLPmYxaEBEwUod+3TxZumzIlKgpFtBCHDIEn\nLU0W1/tNiSfRg4oFCSv8yZNgL170DTIM2gIc4yLqUeoW1H31FdDaqkE2JNSoWJCwYvjiC1nMPXo0\nJP/9FIjm2iZMgGQw+MQYlwv6/fs1yoiEEhULEjbY8+dpYDuSGI1wK+wpYvjiC7qjOwpRsSBhQ+mq\nQkxOhhDhd2xHs1a/tdwAgL1wAfzRoxpkQ0KJigUJC0xzM/RffSWLt86YQQPbYUxMSYGQmSmLGz7/\nXINsSChRsSBhQb9nDyAIPjEpLo4GtiNA6003yWJ8VRXYc+c0yIaEChULoj1BgGH3blm4bcoUwG8A\nlYQfYdQoiIMGyeJK3YokclGxIJrTf/01mKYm3yDLtndBkfDHMIpjF/ryctqjO4pQsSDa8nhg2LZN\nFnZ/5zuQuljFlISXtqwsSPHxvkGPB8Z//EObhEjQUbEgmmL27AGrsD9J6/e+p0E2pNf0erRmZ8vD\ne/eC8b/JkkQkKhZEO4IA9pNP5OEbb1RcSoKEt7aZMyH5bVZFVxfRg4oF0Yx+3z5A4arCdcstGmRD\n+kqKj1ceu9i7F+iwJw2JTFQsiDYEAYbt22Vh95gx8AwdqkFCJBjabrpJfnUhimA3b9YmIRI0VCyI\nJvS7d4NVmCnT+i//okE2JFikuDjF+y6Y3bvBNjRokBEJFioWRHVMc7Nsvwrg2lWFwl4IJLK0zpwp\n30lPFBH3wQe0ZlQEU2WnPAAoLy9HUVERRFFEbm4u8vPzfZ53u91Yv349qqqqkJiYiCVLlmDQoEFo\naGjA0qVLvZtyZGZm4oEHHlArbRICxi1bwLhcvkGGgYuuKqJDXBxab74ZRr+uJ/7ECfDHjkEYM0aj\nxEhfqFIsRFFEYWEhVq5cCYvFghUrVsBut2Noh77p7du3IyEhAevWrcOuXbvwl7/8BUuXLgUApKSk\nYPXq1WqkSkKMq6lR3F+7zW6HSFcVUaP1ppug//JL2d4kce+/j6bMTECn0ygz0luqdENVVlYiJSUF\nVqsVPM8jOzsbZX47au3btw+zZs0CAEybNg2HDx9GlO74GrskCXF//7usK0IyGuH6wQ80SoqEhE4H\n1w9/KAuzDgctMhihVLmycDgcsFgs3scWiwUVFRWdHsNxHOLj49F0bQmIhoYGLFu2DHFxcbjrrrsw\nWmHJ6pKSEpSUlAAACgoKkJycHKrmhAzP8xGZd6CY3bvB1tcDHWbLsCwLw+23Q5+erl1iKov299lr\n1iywBw+COXECLMvCeO19N5aWwjN7NhDlP4Noe59VG7PorQEDBuCVV15BYmIiqqqqsHr1aqxZswbx\nfksL5OXlIS8vz/s40jZKByJzg/dAMQ4HEt98UzZWYbzhBnw7diwQpe1WEs3vsz82NxeJhw7BqNfD\ndf29d7kg/Pd/4+qDD0b18vOR+D5fHxtWoko3lNlsRmOHm3IaGxthNps7Pcbj8aClpQWJiYnQ6XRI\nTEwEAGRkZMBqtaKurk6NtEmwSBLi335bPqgNwDNvHsBxGiRF1CCmpKB1+nRZnK+qgp66oyKKKsXC\nZrOhrq4ODQ0NEAQBpaWlsPvtUzBp0iTs2LEDALBnzx6MHTsWDMPgypUrEEURAFBfX4+6ujpYrVY1\n0iZBot+1C/ypU7J4W1YWQLvgRT3XLbcAAwbI4nGbN7d3S5KIoEo3FMdxWLhwIVatWgVRFJGTk4O0\ntDQUFxfDZrPBbrdj9uzZWL9+PR555BGYTCYsWbIEAHD06FG89dZb4DgOLMvi/vvvh8lkUiNtEgRs\nfT3iFNZ/EpOS4LrtNiRokBNRWVwcPHffDfzud75xQUB8cTGaFy0C+LDvEY95jBSlU45qa2u1TqHH\nIrGPsyuM0wnT+vVgv/1W9tzV+++HkJkZdW0ORKy2uflPf1LcEKlt8mQ4582LuvGLSHyfNR+zIDFI\nFBH/178qForWGTMU920m0c31gx8o7qinLyuDvrRUg4xIT1CxICFh/Phj8CdOyOLiwIF0T0Ws0unQ\nMn++4oSGuA8+AF9ZqUFSJFBULEjQ6ffsgWHnTllcMhpx9ec/B/R6DbIi4cCTlgan31I/ANqvRN98\nEyzNdAxbVCxIUOnLyhD33nvyJxgGLQsWKHZDkNjSNnUq2hR21WOcTphef51mSIUpKhYkaHT79iHu\nnXcUn3P94AcQbrxR5YxIuHL+6EcQMjJkcaa5GabXXqPlzMMQFQsSFPrSUsS//bbiEtRtkyah9eab\nNciKhC2OQ8vdd0PssAzQdUxTE0x//CO4mhoNEiOdoWJB+kYUEfe3vyHub39TLBTu73wHzjvvjLpp\nkaTvpIQEND/wAESFG/aYpiaYXn0VuoMHNciMKKFiQXqNaW5GQmFhp9Me3WPHomXBAoClf2ZEmTRg\nAJoffFCxYMDtRvzGjTBs3QpcW8WBaId+i0mv8IcPI3HtWvB+qwdf5x4zBi0//Smt+0S6JZnNuPrg\ngxD791d83lhSAtOrryres0PUQ8WC9AjT1IS4//s/JPz5z2CamxWPaZs6FS13301LOJCAiWYzmhct\ngqeTO4i5M2eQ+F//Bf3OnYAgqJwdASJgiXISJlwuGHbuhOHzz8G0tiofwzBw/uhHaJsxg8YoSI9J\n/fuj+aGHEF9cDN3hw/ID3G7EffghDLt3wzVnDtzjx9O/MxVRsSBdYpqaoN+9G4Y9ezq9kgAAKT4e\nLT/5CYRRo1TMjkQdgwEtd98Nw6efwrhtm+KkCbaxEfF/+Qs8//gHWmfOhHvCBNqmVQVULIicIICv\nqID+66+hO3QI8Hi6PNw9Zgycd9wB6dq+I4T0CcOg9ZZbIGRmIv6tt8B22AunI662FvFvvQXp44/R\nZrejbcIEiEOG0NVGiFCxIAAA5soV8JWV4CsqoDt6FIzT2e05ktEI59y5cE+aRL+gJOg8N9yApiVL\nEPfJJ10uNMg0N8OwYwcMO3ZATE6Ge+xYCCNGQEhPBwwG9RKOclQsYo0kgWluBtvQAO7cOXC1teC/\n+aZnM004Dq3Tp6M1NxdSAu1IQULIYIAzPx+tkycj7uOPO519dx174QIMn30Gw2efARwHT2oqPEOH\nQkhNhTh4MDzJyT57wJPAUbGIBpIECAIYl6v9y+kEc/Vqe1G4ehXM5ctgL14Ee/ky2MZGxe1NA8Ky\naJs4Ea233ALRb1tcQkJJTE3F1fvvB19RAeOWLeDOnu3+JI8H3Nmz4M6eRcelK6V+/SCazRD794eY\nlASxXz9ICQmQEhMhxcdDiouDZDRCMhho6ncHqhWL8vJyFBUVQRRF5ObmIt9v5Um3243169ejqqoK\niYmJWLJkCQZdW3Ru06ZN2L59O1iWxb333ouJEycGPT9DSQnYpqbADg50v6iOx0nSPx8r/MmIItiE\nBMRfudIeF0Uwoth+M9K1L0YQ2ouCxwO43WDc7vY/W1sDz6kXJJMJbVOnonXaNEhJSSF7HUK6I2Rm\nojkzE9zZszDs2gXdgQM9vmGPuXIF3JUrCKgMcBwkvR6SXg/wPCSdrn1KOMdB4rj2YsKy7X8H2m9A\nvfbFJiYi7urVf3bRKv3Zsfs2kK7cHnT3ujMzIYwbF/Dx3VGlWIiiiMLCQqxcuRIWiwUrVqyA3W7H\n0KFDvcds374dCQkJWLduHXbt2oW//OUvWLp0KWpqalBaWoq1a9fi4sWLePbZZ/H73/8ebJDvCtYf\nOKD5apeM0Qhdbz/1B5lkMEAYPRptEyZAGDmSZpuQsOIZNgwtw4aB+dGPoDt0CLoDB8BXVwf/Q5PH\n036lHsAYnj/GaIRew99nyWCIvGJRWVmJlJQUWK1WAEB2djbKysp8isW+fftw5513AgCmTZuGP/3p\nT5AkCWVlZcjOzoZOp8OgQYOQkpKCyspKjBw5Uo3UY4pn8OD2gcHMTAg2GxUIEvYkkwlt06ejbfr0\n9kkaFRXtEzUqK8Fevqx1elFFlWLhcDhg6bC6pMViQYXfQFXHYziOQ3x8PJqamuBwOJDZYQtOs9kM\nh8Mhe42SkhKUlJQAAAoKCpCcnNyjHLn4eM0HvliWhVGNHHQ6SFYrYLVCSkuDNGwYMGwYdPHxoX9t\nPzzP9/i9inTU5hBJTgYyMoDvf7/9CsPhAHP2LHD2LJiaGjANDcCFCyHtsu1Itd/nTuj79YMpiD/z\nqBngzsvLQ15envdxTzdKT2xpAatxF5DRaISrtzlwHCSDoX1w7tqXmJjoHbjzDuYNGNA+7uDf99nS\n0v6lskjc1L6vqM0qSktr/7pOEMA6HO2TPS5dAnPpUvskkKam9kkhTuc/J4n0cSywT7/PQdB65Qpc\nPfyZD+lkuRVApWJhNpvR2OHGmsbGRpj9ZtNcP8ZiscDj8aClpQWJiYmycx0Oh+zcYHDl5oLpyX+W\nXQ00+T/nP4jFMJA6PmZZgGGgHzAALZcvtw+YXY93HEC79iXxfPvVAccBBsM/B90IIV3jeYiDBgW2\nY+P1WYatre1/CgLQ1tY+weTaF+Px+E5CuTY5BaIIfVISWi5e9E5uYTp+X//JL/6v21k+PeBJTe3R\n8d1R5X8Ym82Guro6NDQ0wGw2o7S0FI8++qjPMZMmTcKOHTswcuRI7NmzB2PHjgXDMLDb7fjDH/6A\nH/7wh7h48SLq6uowYsSIoOfoDsEMq56SkpPhjrFPnISELYZp/1B2beyup9cY0fb7rEqx4DgOCxcu\nxKpVqyCKInJycpCWlobi4mLYbDbY7XbMnj0b69evxyOPPAKTyYQlS5YAANLS0jB9+nQ8/vjjYFkW\n9913X9BnQhFCCOkaI0kqjfaorLa2VusUeoz6smMDtTk2RGKbuxqzoI/ohBBCukXFghBCSLeoWBBC\nCOkWFQtCCCHditoBbkIIIcFDVxZhZPny5VqnoDpqc2ygNkc+KhaEEEK6RcWCEEJIt6hYhJGOCyHG\nCmpzbKA2Rz4a4CaEENIturIghBDSLSoWhBBCukWbIISpDz74AG+++Sb+53/+B/369dM6nZB68803\nsX//fvA8D6vVikWLFiEhIUHrtIKuvLwcRUVFEEURubm5yM/P1zqlkLpw4QJefvllXLp0CQzDIC8v\nD7feeqvWaalCFEUsX74cZrM5aqbQUrEIQxcuXMDBgwdjZuvN8ePHY8GCBeA4Dhs3bsSmTZvws5/9\nTOu0gkoURRQWFmLlypWwWCxYsWIF7Ha7zz700YbjONx9993IyMiA0+nE8uXLMX78+Khu83Uff/wx\nUlNT4XQ6tU4laKgbKgxt2LABP/3pT8F0tRtfFJkwYQI4jgMAjBw5UnGP9UhXWVmJlJQUWK1W8DyP\n7OxslJWVaZ1WSA0YMAAZGRkAgLi4OKSmpkble+uvsbERX331FXJzc7VOJaioWISZsrIymM1mpKen\na52KJrZv346JYbBrYbA5HA5YLBbvY4vFEhP/cV7X0NCA6urqkOxyGW7eeOMN/OxnP4u6D3vUDaWB\nZ599FpcuXZLF77rrLmzatAkrV67UIKvQ6qrNkydPBgC899574DgON910k9rpkRByuVxYs2YN7rnn\nHsTHx2udTkjt378fSUlJyMjIwJEjR7ROJ6joPoswcvbsWTzzzDMwGAwA2i9nBwwYgOeffx79+/fX\nOLvQ2rFjBz799FM89dRT3vZHk5MnT+Ltt9/Gk08+CQDYtGkTAODf/u3ftEwr5ARBwAsvvIAJEybg\nhz/8odbphNxf//pX7Ny5ExzHoa2tDU6nE1OmTMGjjz6qdWp9J5GwtWjRIuny5ctapxFyX3/9tbRk\nyZKobqsgCNLixYul+vp6ye12S7/85S+ls2fPap1WSImiKK1bt04qKirSOhVNHD58WHr++ee1TiNo\nqBuKaK6wsBCCIODZZ58FAGRmZuKBBx7QOKvg4jgOCxcuxKpVqyCKInJycpCWlqZ1WiF14sQJ7Ny5\nE8OGDcN//Md/AAB+8pOfICsrS+PMSG9QNxQhhJBu0WwoQggh3aJiQQghpFtULAghhHSLigUhhJBu\nUbEghBDSLSoWhBBCukXFghBCSLeoWBBCCOkWFQtCVHD+/Hnce++9qKqqAtC+Cu19990XdYvNkehF\nxYIQFaSkpOCnP/0p1q1bh9bWVrz66qu4+eabMXbsWK1TIyQgtNwHISp64YUX0NDQAIZh8Pzzz0On\n02mdEiEBoSsLQlSUm5uLb775BnPmzKFCQSIKFQtCVOJyubBhwwbMnj0bb7/9Npqbm7VOiZCAUbEg\nRCVFRUXIyMjAL37xC2RlZeG1117TOiVCAkbFghAVlJWVoby8HPfffz8A4Oc//zmqq6vx+eefa5wZ\nIYGhAW5CCCHdoisLQggh3aJiQQghpFtULAghhHSLigUhhJBuUbEghBDSLSoWhBBCukXFghBCSLeo\nWBBCCOnW/wcjw1l8rrkGGAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XRJ4IK3kmMon",
        "colab_type": "text"
      },
      "source": [
        "Sample normal observations and take exponents (lognormally distributed)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "g-eK0jPw0Cg3",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Sample normally distributed energies\n",
        "energies = sigma * rn.randn(20) + mean\n",
        "p_hat = np.exp(energies)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zy87yk_QmTpX",
        "colab_type": "text"
      },
      "source": [
        "Compute ML estimates of the mean and variance (equal for both normal and lognormal distribution)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6pw-pntx2C4j",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# ML estimates of parameters\n",
        "mean_ml = np.sum(np.log(p_hat)) / p_hat.size\n",
        "var_ml = np.sum((np.log(p_hat) - mean_ml)**2) / p_hat.size"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZhhYa0VenJKn",
        "colab_type": "text"
      },
      "source": [
        "Compute likelihood of the data under lognormal and normal. PDFs are not equal:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ozC6We5T21jW",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "39d40382-4d7f-4198-f67a-5dcaed87c344"
      },
      "source": [
        "# Compute lognormal likelihood of n probabilities\n",
        "lognormal_lik = np.sum(np.log(lognormal_pdf(p_hat, mean_ml, np.sqrt(var_ml))))\n",
        "# Compute normal likelihood of n energies\n",
        "normal_lik = np.sum(np.log(normal_pdf(energies, mean_ml, np.sqrt(var_ml))))\n",
        "print(f\"PDFs: \", lognormal_lik, normal_lik)"
      ],
      "execution_count": 33,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "PDFs:  -2.731334262422567 -0.8733253276697781\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d0Npe_mVsJ82",
        "colab_type": "text"
      },
      "source": [
        "But CDFs are!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CgMAa-cvsBvL",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "4dfa60b0-6e93-4ab1-9cae-1a9119af9394"
      },
      "source": [
        "# Compute lognormal likelihood of n probabilities\n",
        "lognormal_lik = np.sum(np.log(lognormal_cdf(p_hat, mean_ml, np.sqrt(var_ml))))\n",
        "# Compute normal likelihood of n energies\n",
        "normal_lik = np.sum(np.log(normal_cdf(energies, mean_ml, np.sqrt(var_ml))))\n",
        "print(f\"CDFs: \", lognormal_lik, normal_lik)"
      ],
      "execution_count": 34,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "CDFs:  -20.075537356794374 -20.075537356794374\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}