{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "e83a02ca-6265-4d13-9a29-7e08e7ac7568",
      "metadata": {
        "id": "e83a02ca-6265-4d13-9a29-7e08e7ac7568"
      },
      "source": [
        "## Computational\n",
        "\n",
        "The solutions to problems in this section should be in the form of code. You might be required to state your observations, in which case, an empty markdown cell will be provided."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "204ccbde-07b7-467b-a0eb-6357c24dab7c",
      "metadata": {
        "id": "204ccbde-07b7-467b-a0eb-6357c24dab7c"
      },
      "source": [
        "### Problem 1 - MLE [10 points]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "f280efc4-5207-46b4-9a16-d2a656933b2c",
      "metadata": {
        "id": "f280efc4-5207-46b4-9a16-d2a656933b2c"
      },
      "outputs": [],
      "source": [
        "# Import the following libraries\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import scipy.stats as stats\n",
        "import scipy.optimize as optimize\n",
        "import pandas as pd\n",
        "\n",
        "from scipy.stats import norm\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "2cdf87eb-7f11-4d8d-b5ae-3b22a41cf18b",
      "metadata": {
        "id": "2cdf87eb-7f11-4d8d-b5ae-3b22a41cf18b"
      },
      "outputs": [],
      "source": [
        "# generate normally distributed data\n",
        "np.random.seed(123) # Use a random seed of 123\n",
        "data = np.random.normal(loc=3, scale=2, size=100) # Use random.normal to generate normally distributed data of size 100"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# printing the first 10 numbers to check the data generated\n",
        "print(data[:10])\n",
        "\n",
        "# NOTE: It can be different from the output below\n",
        "# It should give you a general idea about how your output should look like"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "iJcFb2Sh_BTz",
        "outputId": "f4dc7ffa-85d5-4e39-c67a-7bcc95611704"
      },
      "id": "iJcFb2Sh_BTz",
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[ 0.82873879  4.99469089  3.565957   -0.01258943  1.8427995   6.30287307\n",
            " -1.85335849  2.14217474  5.53187252  1.2665192 ]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "5106ab7f-4174-4bd6-8ed1-cc425f4234c9",
      "metadata": {
        "id": "5106ab7f-4174-4bd6-8ed1-cc425f4234c9"
      },
      "outputs": [],
      "source": [
        "def likelihood(x, mu, sigma):\n",
        "  p = 1 / np.sqrt(2 * np.pi * sigma**2) * np.exp(-0.5 * (x - mu)**2 / sigma**2)\n",
        "  l_hood = np.prod(p)\n",
        "  return l_hood"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5fc327e7-c384-4783-aae1-9fb06dfccca8",
      "metadata": {
        "id": "5fc327e7-c384-4783-aae1-9fb06dfccca8"
      },
      "source": [
        "Formula for the probability density function (PDF) of the normal (or Gaussian) distribution for a given value x.\n",
        "\n",
        "$$ p = -0.5 \\times \\log(2 \\times \\pi \\times \\sigma^2) - \\frac{0.5 \\times (x - \\mu)^2}{\\sigma^2} $$\n",
        "\n",
        "Use this formula to fill the blank below, function parameters are provided for you\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "2ce3d1b4-0904-4100-8d08-2ae05b81197f",
      "metadata": {
        "id": "2ce3d1b4-0904-4100-8d08-2ae05b81197f"
      },
      "outputs": [],
      "source": [
        "def log_likelihood(x, mu, sigma):\n",
        "  p = -0.5 * np.log(2 * np.pi * sigma**2) - 0.5 * (x - mu)**2 / sigma**2\n",
        "  log_lhood = np.sum(p)\n",
        "  return log_lhood"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ceba8dcd-f48a-4bfc-a998-413405902e47",
      "metadata": {
        "id": "ceba8dcd-f48a-4bfc-a998-413405902e47"
      },
      "source": [
        "### Formula for Maximum Likelihood Estimation:\n",
        "\n",
        "$$ \\text{NLL}(\\theta) = -\\sum_{i} \\log \\left( L(x_i; \\theta) \\right) $$\n",
        "\n",
        "Where:\n",
        "\n",
        "- $ \\text{NLL}(\\theta) $: Represents the negative log-likelihood for parameters $ \\theta $.\n",
        "- $ \\sum_{i} $: Indicates the summation over all data points.\n",
        "- $ \\log $: Denotes the natural logarithm.\n",
        "- $ L(x_i; \\theta) $: Stands for the likelihood of observing data point $ x_i $ given the parameters $ \\theta $.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "d086a093-9ff9-4794-af4c-12afe122587c",
      "metadata": {
        "id": "d086a093-9ff9-4794-af4c-12afe122587c"
      },
      "outputs": [],
      "source": [
        "# Maximum likelihood estimation can be defined as:\n",
        "\n",
        "from scipy.optimize import minimize\n",
        "\n",
        "def max_likelihood_estimation(x):\n",
        "    # Defining the negative log-likelihood function\n",
        "    neg_log_likelihood = lambda params: -np.sum(np.log(likelihood(x, *params)))\n",
        "\n",
        "    # Set the initial guess for the parameters and optimize\n",
        "    mu_init = np.mean(x) # Set the mean\n",
        "    sigma_init = np.std(x) # Set the initial standard deviation\n",
        "    params_init = [mu_init, sigma_init]\n",
        "    res = minimize(neg_log_likelihood, params_init, method='Nelder-Mead')\n",
        "\n",
        "    # Return the estimated parameters and negative log-likelihood\n",
        "    return res.x, -res.fun"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Running the funtions\n",
        "\n",
        "mean_mle, std_dev_mle = max_likelihood_estimation(data)\n",
        "mu_mean = np.mean(data)\n",
        "sigma_mean = np.var(data, ddof=1)\n",
        "print(\"Maximum likelihood estimate of mu:\", mean_mle)\n",
        "print(\"Sample mean of x:\", mu_mean)\n",
        "print(\"Maximum likelihood estimate of sigma^2:\", std_dev_mle)\n",
        "print(\"Sample variance of x:\", sigma_mean)\n",
        "\n",
        "# Note - mu and sigma are used here to avoid confusion.\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jlzgFZpJCYiB",
        "outputId": "67f3f4f7-6455-4ec1-8bf9-2c71fed3b0dd"
      },
      "id": "jlzgFZpJCYiB",
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Maximum likelihood estimate of mu: [3.05421815 2.25648094]\n",
            "Sample mean of x: 3.05421814698072\n",
            "Maximum likelihood estimate of sigma^2: -223.2745027168046\n",
            "Sample variance of x: 5.143137613027599\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Estimate the maximum likelihood parameters\n",
        "params, neg_log_likelihood = max_likelihood_estimation(data)\n",
        "mu_ml, sigma_ml = params\n",
        "\n",
        "# Generate the PDF of the estimated distribution\n",
        "x = np.linspace(np.min(data), np.max(data), 100)\n",
        "pdf = norm.pdf(x, loc=mu_ml, scale=sigma_ml)\n",
        "\n",
        "# Plot the data and the estimated distribution\n",
        "plt.hist(data, bins=10, density=True, alpha=0.5, label='Data')\n",
        "plt.plot(x, pdf, label='Estimated PDF')\n",
        "plt.xlabel('x')\n",
        "plt.ylabel('Probability density')\n",
        "plt.legend()\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 449
        },
        "id": "OTwEX3p-C3wz",
        "outputId": "1da937d8-9eee-4a3d-9cda-e83a82bf6065"
      },
      "id": "OTwEX3p-C3wz",
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAGwCAYAAABSN5pGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABvTklEQVR4nO3dd3wUdf7H8ddueoAkQEiDQCjSewtVRKMgoqIcAuKBiFiOHvUO/HnC2UBFjCiCcoCoICgiKofxIAgondCR3hISUmgJSUjb3d8fK+uFmoSESXk/H495OJmdnXnvumQ/+c53vl+TzWazISIiIiIOZqMDiIiIiJQ0KpBERERErqACSUREROQKKpBERERErqACSUREROQKKpBERERErqACSUREROQKzkYHKK2sVivx8fFUqlQJk8lkdBwRERHJB5vNxsWLFwkKCsJsvn47kQqkQoqPjyc4ONjoGCIiIlIIsbGx1KhR47qPq0AqpEqVKgH2N9jLy8vgNCIiIpIfqampBAcHO77Hr0cFUiFdvqzm5eWlAklERKSUuVn3GHXSFhEREbmCCiQRERGRK6hAEhEREbmC+iCJiIihLBYLOTk5RseQMsLFxQUnJ6dbPo4KJBERMYTNZiMhIYELFy4YHUXKGB8fHwICAm5pnEIVSCIiYojLxZGfnx+enp4adFdumc1mIyMjg6SkJAACAwMLfSwVSCIicttZLBZHcVS1alWj40gZ4uHhAUBSUhJ+fn6FvtymTtoiInLbXe5z5OnpaXASKYsuf65upW+bCiQRETGMLqtJcSiKz5UKJBEREZErqEASERERuYI6aYuISIny/spDt+1c4+6tf9vOJaWLWpBEREQK4Mknn8RkMmEymXBxccHf3597772XuXPnYrVa832czz77DB8fn+ILKrdELUgiIkXNkguZKZB5wb5cugCWHKjkDxUDoEI1cNKv39KsZ8+ezJs3D4vFQmJiIpGRkYwZM4YlS5bwww8/4Oys/7+lnf4PiojcqvQzELsZYjbZl9M7wZJ9/f1NZqjgBwFNoc5dUKc7+DcB3dFVari5uREQEABA9erVad26NR06dOCee+7hs88+4+mnn2batGnMmzePY8eOUaVKFR588EHeeecdKlasyJo1axg6dCjw5x1XEydOZNKkSXzxxRd88MEHHDx4kAoVKnD33XcTERGBn5+fYa+3PDL8EtuMGTMICQnB3d2d0NBQtmzZct199+3bR9++fQkJCcFkMhEREXHVPpcfu3IZMWKEY5+77rrrqsefe+654nh5IlJWZaVB9HyYfTe8WxcWPQ4bpsOpLX8WR66VwKsG+DeFwJZQKQhMTmCzQloCHFkF/30FZnWGqXfA0mfg+K9gsxn60qRw7r77blq0aMHSpUsBMJvNTJ8+nX379jF//nxWr17N3//+dwA6depEREQEXl5enD59mtOnT/Piiy8C9rF7Xn/9dXbt2sWyZcs4ceIETz75pFEvq9wytAVp8eLFhIeHM2vWLEJDQ4mIiKBHjx4cPHjwmpVyRkYGderUoV+/fowbN+6ax9y6dSsWi8Xx8969e7n33nvp169fnv2GDx/Oa6+95vhZg5WJSL7E74Doz2DPEshO+3N7tYYQHAo1O9j/61Pr2pfRrBZIT4aUOIjZCMfWwMn19m27F9sX/6YQ+iw06wcuHrfrlUkRaNiwIbt37wZg7Nixju0hISG88cYbPPfcc3z88ce4urri7e2NyWRytERd9tRTTznW69Spw/Tp02nXrh1paWlUrFjxtrwOMbhAmjZtGsOHD3c0M86aNYv//Oc/zJ07l/Hjx1+1f7t27WjXrh3ANR8HqFatWp6fp0yZQt26denWrVue7Z6enld9KG8kKyuLrKwsx8+pqan5fq6IlAHJh2DlP+FQ5J/bqtSFNkOgxUComM/LH2YnqBRgX2q0gU4jITcLTm2Fvd/CrkWQuBd+GAUrJ0Loc9BpFLjqj7jSwGazOS6ZrVq1ismTJ3PgwAFSU1PJzc0lMzOTjIyMG/5RHh0dzaRJk9i1axfnz593dPyOiYmhcePGt+V1iIGX2LKzs4mOjiYsLOzPMGYzYWFhbNy4scjO8eWXX/LUU09dNarmggUL8PX1pWnTpkyYMIGMjIwbHmvy5Ml4e3s7luDg4CLJKCIlXPoZ+M8L8HEHe3FkcoKmfWHIchgVDZ3H5L84uh5nNwjpAr3fh/Df4d7XwDsYLp2DNW/BjFDY/6MuvZUC+/fvp3bt2pw4cYLevXvTvHlzvv32W6Kjo5kxYwZg/266nvT0dHr06IGXlxcLFixg69atfPfddzd9nhQ9w1qQzpw5g8Viwd/fP892f39/Dhw4UCTnWLZsGRcuXLjq2u3jjz9OrVq1CAoKYvfu3fzjH//g4MGDjuvG1zJhwgTCw8MdP6empqpIEinLrFbYPBPWTIGsP1qMG/SyFy++dxTfeT0q24uuDiPg92X2VqSUGFj8BNS9G+5/p3jPL4W2evVq9uzZw7hx44iOjsZqtfLee+9hNtvbIr7++us8+7u6uubpEgJw4MABzp49y5QpUxzfMdu2bbs9L0DyKNN3sc2ZM4f777+foKCgPNufeeYZx3qzZs0IDAzknnvu4ejRo9StW/eax3Jzc8PNza1Y84pICZEaD989C8fX2X8ObAH3vQG177x9GZycodlfoMH98Nv7sP4DOLoaPu4I3V+GzmPBbPh9NuVWVlYWCQkJeW7znzx5Mr1792bw4MHs3buXnJwcPvzwQx588EHWr1/PrFmz8hwjJCSEtLQ0oqKiaNGiBZ6entSsWRNXV1c+/PBDnnvuOfbu3cvrr79u0Kss52wGycrKsjk5Odm+++67PNsHDx5se+ihh276/Fq1atnef//96z5+4sQJm9lsti1btuymx0pLS7MBtsjIyJvue1lKSooNsKWkpOT7OSJSCvz+o802pZbNNtHLZnsjwGbbOtdms1iMTmWznTlis33Zz55ropfN9uVfbLa0M0anKrRLly7Zfv/9d9ulS5eMjlJgQ4YMsQE2wObs7GyrVq2aLSwszDZ37lyb5X8+K9OmTbMFBgbaPDw8bD169LB9/vnnNsB2/vx5xz7PPfecrWrVqjbANnHiRJvNZrMtXLjQFhISYnNzc7N17NjR9sMPP9gA244dO27vCy3FbvT5yu/3t8lmM+6idmhoKO3bt+fDDz8EwGq1UrNmTUaOHHndTtiXhYSEMHbs2Dx3CfyvSZMm8cknnxAbG3vTAbvWr19Ply5d2LVrF82bN89X9tTUVLy9vUlJScHLyytfzxGREiznEvz8Mmyba/85sAX0nVOyLmfZbLDjC1jxEuRmgld1+Mtc+51zpUxmZibHjx+ndu3auLu7Gx1Hypgbfb7y+/1t6CW28PBwhgwZQtu2bWnfvj0RERGkp6c77mobPHgw1atXZ/LkyYC9g9rvv//uWI+Li2Pnzp1UrFiRevXqOY5rtVqZN28eQ4YMuao4Onr0KAsXLqRXr15UrVqV3bt3M27cOO688858F0ciUsakJcNXAyDuj74enUbD3f8EZ1djc13JZILWgyGoNXwzBM4egXm9IGyS/U43DTQpUmQMLZD69+9PcnIyr776KgkJCbRs2ZLIyEhHx+2YmBhH5zaA+Ph4WrVq5fh56tSpTJ06lW7durFmzRrH9lWrVhETE5NnLInLXF1dWbVqlaMYCw4Opm/fvrzyyivF90JFpOQ6cxgW/AXOn7B3kP7LXHtn6JIsoCk8swZ+HAt7l9iHHzh/Anq9ax9GQERumaGX2EozXWITKQNOboRFA+HSeagcAoOWlKxLajdjs8GWT+GnfwA2aPQQPDobXEr+JStdYpPiVBSX2HQLhIiUT3uXwucP24uj6m1h2KrSVRyB/ZJa6LPQbx44ucL+H+DLvvaJckXklqhAEpHyZ+dXsOQpsGRBw94w5EeoWO3mzyupmjwCT3xrn/vt5G8w7wG4mGB0KpFSTQWSiJQvu7+B7/8G2KDtU/DY52VjGo/ad8LQ/0AFP0jcA/MftI8CLiKFogJJRMqPfd/Bd8+AzQqth0Cv98pWp+bAFjDsv+BVA84cgi/6wKULRqcSKZVUIIlI+bD/R1gyzF4ctXwCekeUzZGoq9SGwd/bW5IS9tjv0MtKMzqVSKlTBn87iIhc4dB/4ZuhYLNA8wHw0PSyWRxd5lsPBi8Ddx84tdU+xlPOJaNTyQ189tln+Pj4GB2jQEpj5oIow78hRESA+J32QRWtOdC0L/T5uGxdVrse/ybw16X2jtsnfoWvh4Alx+hUZcKTTz6JyWS6aunZs2e+nh8SEkJERESebf379+fQoUPFkDav213U/O/74+3tTefOnVm9erXj8f99L11cXPD39+fee+9l7ty5WK3WPMcKCQm56j2vUaNGsWVXgSQiZVfKKVjYH3IyoE53eOST8lEcXVa9DTy+GJw94PDP9vGSNPRdkejZsyenT5/Os3z11VeFPp6Hhwd+fn5FmLDkmDdvHqdPn2b9+vX4+vrSu3dvjh075nj88nt54sQJfvrpJ7p3786YMWPo3bs3ubm5eY712muv5XnPd+zYUWy5VSCJSNmUmQoLHoO0BKjWCB6bD04uRqe6/UI620cHxwTb5tgHlpRb5ubmRkBAQJ6lcuXKANhsNiZNmkTNmjVxc3MjKCiI0aNHA3DXXXdx8uRJxo0b52gFgatbdiZNmkTLli2ZO3cuNWvWpGLFivztb3/DYrHwzjvvEBAQgJ+fH2+++WaeXNOmTaNZs2ZUqFCB4OBg/va3v5GWZu+DtmbNGoYOHUpKSorj3JMmTQIgKyuLF198kerVq1OhQgVCQ0PzzFBxOWPNmjXx9PTkkUce4ezZs/l6r3x8fAgICKBp06bMnDmTS5cusXLlyqvey+rVq9O6dWtefvllvv/+e3766Sc+++yzPMeqVKlSnve8WrXiG55DBZKIlD2WXFgyFJL2QUV/GPQ1uHsbnco4DXvBvf+yr0eOh8OrjM1zPTYbZKff/qWIW9W+/fZb3n//fT755BMOHz7MsmXLaNasGQBLly6lRo0aeVpCrufo0aP89NNPREZG8tVXXzFnzhweeOABTp06xdq1a3n77bd55ZVX2Lx5s+M5ZrOZ6dOns2/fPubPn8/q1av5+9//DkCnTp2IiIjAy8vLce4XX3wRgJEjR7Jx40YWLVrE7t276devHz179uTw4cMAbN68mWHDhjFy5Eh27txJ9+7deeONNwr83nh4eAD2+VRv5O6776ZFixYsXbq0wOcoKobOxSYiUuRsNvjpJTiyyn5paeAi8KlpdCrjdRoNyYdg55f24nHYSvBraHSqvHIy4K2g23/el+PBtUKBnrJ8+XIqVqyY9zAvv8zLL79MTEwMAQEBhIWF4eLiQs2aNWnfvj0AVapUwcnJydESciNWq5W5c+dSqVIlGjduTPfu3Tl48CArVqzAbDbToEED3n77bX755RdCQ0MBGDt2rOP5ISEhvPHGGzz33HN8/PHHuLq64u3tjclkynPumJgY5s2bR0xMDEFB9vf/xRdfJDIyknnz5vHWW2/xwQcf0LNnT0exVb9+fTZs2EBkZGS+37OMjAxeeeUVnJyc6Nat2033b9iwIbt3786z7R//+EeeuVPfeustR+tcUVOBJCJlS/Q82PbHJaW+/4bqrY1OVDKYTNB7Gpw7BjEb4Kv+8PRqqFDV6GSlUvfu3Zk5c2aebVWqVAGgX79+REREUKdOHXr27EmvXr148MEHcXYu2FduSEgIlSpVcvzs7++Pk5NTnknc/f39SUpKcvy8atUqJk+ezIEDB0hNTSU3N5fMzEwyMjLw9Lz2gKh79uzBYrFQv379PNuzsrKoWtX++di/fz+PPPJInsc7duyYrwJp4MCBODk5cenSJapVq8acOXNo3rz5TZ9ns9kclyAve+mll3jyyScdP/v6+t70OIWlAklEyo647X9M3Arc8yo06m1snpLG2Q36fwmzu8P5E/a7+wZ/X3I6rrt42ltzjDhvAVWoUIF69epd87Hg4GAOHjzIqlWrWLlyJX/729949913Wbt2LS4u+e8Hd+W+l+/0unLb5bu9Tpw4Qe/evXn++ed58803qVKlCr/99hvDhg0jOzv7ugVSWloaTk5OREdH4+SU97NwZStZYbz//vuEhYXh7e1doD5D+/fvp3bt2nm2+fr6Xvd9L2oqkESkbMg498et7NnQ4AHoMs7oRCVTharw+Ncw+2777f9rpsDd/2d0KjuTqcCXukoqDw8PHnzwQR588EFGjBhBw4YN2bNnD61bt8bV1RWLxVLk54yOjsZqtfLee+85Wpm+/vrrPPtc69ytWrXCYrGQlJRE165dr3nsRo0a5enrBLBp06Z85QoICChwUbN69Wr27NnDuHHG/TtWgSQipZ/VCt89CykxULm2fayjK5rm5X/4NYQHP4ClT8O6d6FmB6h3j9GpSpWsrCwSEvJOCOzs7Iyvry+fffYZFouF0NBQPD09+fLLL/Hw8KBWrVqA/dLZunXrGDBgAG5ubkV2mahevXrk5OTw4Ycf8uCDD7J+/XpmzZqVZ5+QkBDS0tKIioqiRYsWeHp6Ur9+fQYNGsTgwYN57733aNWqFcnJyURFRdG8eXMeeOABRo8eTefOnZk6dSoPP/wwP//8c4H6H93I5ffSYrGQmJhIZGQkkydPpnfv3gwePLhIzlEYKpBESpn3Vxb/YHJFbdy99W++06347T04/F9wdrdPPuvhU7znKwua94OTv0H0Z7D0GXjuN/AKNDpVqREZGUlgYN73q0GDBhw4cAAfHx+mTJlCeHg4FouFZs2a8eOPPzr687z22ms8++yz1K1bl6ysLGxFdBddixYtmDZtGm+//TYTJkzgzjvvZPLkyXmKjE6dOvHcc8/Rv39/zp49y8SJE5k0aRLz5s3jjTfe4IUXXiAuLg5fX186dOhA7972y9QdOnRg9uzZTJw4kVdffZWwsDBeeeUVXn/99VvOffm9dHZ2pnLlyrRo0YLp06czZMiQPP2tbjeTraj+z5QzqampeHt7k5KSgpeXl9FxpBxRgXSFo7/AF48ANnjoI2j91+I7V1mTcwn+fS8k7oFanWHwD+B0e/5uzszM5Pjx49SuXRt3d/fbck4pP270+crv97fGQRKR0iv9jP3SGjZo9VcVRwXl4mEfQNO1EpxcD2veMjqRSImhAklESiebDX4cA2mJUK0h9HrX6ESlU9W69sl7AX59z94iJyIqkESklNrxJRxYDmYXeHS2vTVECqfpo9D2Kfv69yPg0gVD44iUBCqQRKT0OXfsz/GO7n4FAm8+6JzcxH1vQJW6kBoHP/3d6DQihlOBJCKliyXXftdVTjrU6gKdRhmdqGxwrQCPfAImM+xeDPuW3ZbT6j4hKQ5F8blSgSQipctv0+DUVnDzhkdmlZxRoMuC4HbQ9QX7+vJxcDHhxvvfgssjQmdkZBTbOaT8uvy5KsjI5VfSOEgiUnrEbbeP/AzwwFTwCTY2T1l059/h0M+QsBt+GA2PLy6WQTednJzw8fFxzCPm6el51bxbIgVls9nIyMggKSkJHx+fq6ZOKQgVSCJSOuRmw/cjwWaBJo9As35GJyqbnF3h0U/hk25w+GfYPh/aPFksp7o8o/z/TrYqUhR8fHwcn6/CUoEkIqXD+g8gaR94VIFeUzWVSHHya2Sf7Pe//wc/vwL17gXv6kV+GpPJRGBgIH5+fuTk5BT58aV8cnFxuaWWo8tUIIlIyZd8ENa9Y1+//x2oUDRzV8kNdPgb7P8BYjfDihdhwMJiK0qdnJyK5AtNpCipk7aIlGxWi/3SmiUb7ugBzf5idKLywWyGB6fbx5k6uAJ+/97oRCK3lQokESnZtv4bTm2xT4fRe5ourd1Ofg2ha7h9fcVLcOm8sXlEbiMVSCJScp0/Cav+ZV+/91/gXcPYPOVR1xfAtz6kJ8HKV41OI3LbqEASkZLJZrOPxZOTbp9pvs1QoxOVT85u9kttANs/h+O/GptH5DZRgSQiJdPv38PRKHD64wvarF9XhqnV8c+52n4cAzmXjM0jchvoN46IlDxZaRA5wb7eZSz41jM0jgBhk6BSIJw7Cr+9b3QakWKnAklESp61b8PFePCpBV3GGZ1GANy9oecfo5j/FmGfMFikDNM4SCJS7N5feSjf+1bNOMqgnTNwApYFjuX4mtjiC3YD4+6tb8h5S7TGD0Odu+DYGnsL3+OLb/mQBflslBT6bJQPakESkZLDZqP70Xdwslk4UqUbx6t0MTqR/C+TyT5Qp9kZDkXCwUijE4kUGxVIIlJiNEyOJDh1OzlmN9bUfsHoOHIt1RrYR9kGiPwH5GQam0ekmKhAEpESwTU3jTtPRACwucYwLroHGhtIrq/b3+0dts+fsM+RJ1IGGV4gzZgxg5CQENzd3QkNDWXLli3X3Xffvn307duXkJAQTCYTERERV+0zadIkTCZTnqVhw4Z59snMzGTEiBFUrVqVihUr0rdvXxITE4v6pYlIAYTG/psKOec4516T7dUHGR1HbsStEtz3hn39t2n2AT1FyhhDC6TFixcTHh7OxIkT2b59Oy1atKBHjx4kJSVdc/+MjAzq1KnDlClTCAgIuO5xmzRpwunTpx3Lb7/9lufxcePG8eOPP/LNN9+wdu1a4uPjefTRR4v0tYlI/vlciqHVaXuH3zV1XsBidjU4kdxU074Q0hVyM+Hnl41OI1LkDC2Qpk2bxvDhwxk6dCiNGzdm1qxZeHp6Mnfu3Gvu365dO959910GDBiAm5vbdY/r7OxMQECAY/H1/XPm75SUFObMmcO0adO4++67adOmDfPmzWPDhg1s2rSpyF+jiNzcnSc+wMmWy7HKnTlZuZPRcSQ/TCbo9S6YnODAcji21uhEIkXKsAIpOzub6OhowsLC/gxjNhMWFsbGjRtv6diHDx8mKCiIOnXqMGjQIGJiYhyPRUdHk5OTk+e8DRs2pGbNmjc8b1ZWFqmpqXkWEbl1NS9spu65dVhxYl3IWKPjSEH4NfpzhO2f/w+sFmPziBQhwwqkM2fOYLFY8Pf3z7Pd39+fhISEQh83NDSUzz77jMjISGbOnMnx48fp2rUrFy9eBCAhIQFXV1d8fHwKdN7Jkyfj7e3tWIKDgwudUUTsTLZcuh23j8q8K/AvnPcMMTaQFNxdE8DNGxL3wM6FRqcRKTKGd9Iuavfffz/9+vWjefPm9OjRgxUrVnDhwgW+/vrrWzruhAkTSElJcSyxscYMXidSljRN/AHfjKNkOnuxMXi40XGkMCpUhW4v2ddXv26fJkakDDCsQPL19cXJyemqu8cSExNv2AG7oHx8fKhfvz5HjhwBICAggOzsbC5cuFCg87q5ueHl5ZVnEZHCc81No9PJmQBsDH6GLBdvgxNJobV/BirXhrREWB9hdBqRImFYgeTq6kqbNm2IiopybLNarURFRdGxY8ciO09aWhpHjx4lMNA+pkqbNm1wcXHJc96DBw8SExNTpOcVkRsLjZ2DZ+4FznqEsDugr9Fx5FY4u8G9r9nXN3wIKaeMzSNSBAydiy08PJwhQ4bQtm1b2rdvT0REBOnp6QwdOhSAwYMHU716dSZPngzYO3b//vvvjvW4uDh27txJxYoVqVfPPtv3iy++yIMPPkitWrWIj49n4sSJODk5MXDgQAC8vb0ZNmwY4eHhVKlSBS8vL0aNGkXHjh3p0KGDAe+CSPnjlRlHq9OLAFhXeyxWs6aFLPUaPQi1OsPJ9RD1Gjz6qdGJRG6Job+V+vfvT3JyMq+++ioJCQm0bNmSyMhIR8ftmJgYzOY/G7ni4+Np1aqV4+epU6cydepUunXrxpo1awA4deoUAwcO5OzZs1SrVo0uXbqwadMmqlWr5nje+++/j9lspm/fvmRlZdGjRw8+/vjj2/OiRYTOJ2fiZMvlpE8oJyp3NjqOFAWTCXq8CZ92h92Lof2zUKON0alECs1ks9lsRocojVJTU/H29iYlJUX9keS2Ko2zn/8vv7T9DNo1GIAvW3xJcsUGBie6Ns3YXkjfPQ+7FkKtLvDkcnvhdAOl8fOsz0bplt/v7zJ3F5uIlGA2G11PfAjA/mo9S2xxJLfg7v8DJzc4+RscWWV0GpFCU4EkIrdNrQubqJmylVyTCxtqPm90HCkO3jUg9Bn7+sqJGjxSSi0VSCJye9isdDlpbz3aFdiPVPcggwNJsekSbh88Mmkf7PnG6DQihaICSURui0bJkfilHybTqSJbagw1Oo4UJ88q0HWcfX31m5CbZWwekUJQgSQixc7JmkWnGPugkNtqDCHTxcfYQFL8Qp+DSkGQEgNb5xidRqTAVCCJSLFrfvpbvLISuOjqx/bAAUbHkdvBxQPuGm9fX/cuZKYYm0ekgFQgiUixcslNp/2peQBsCh6Oxcnd4ERy27QcBL714dI5WD/d6DQiBaICSUSKVevTX+GZe4Hz7jXZ59/b6DhyOzk5wz0T7esbZ8DFxBvvL1KCqEASkWLjlpNCm7gvAdhQ81lsJk0pUu40fACqt4XcS/DbNKPTiOSbCiQRKTZt477AzZJOsucdHPINMzqOGMFkgrtfsa9vm6uJbKXUUIEkIsXCM/uMY0LaDTWfBZN+3ZRbde6CkK5gyYa17xidRiRf9BtLRIpF+1Of4WLN4nTFphyrcqfRccRIJhN0/z/7+o4v4exRY/OI5IMKJBEpcpUyT9MsYSkA62s9f9MJS6UcqNUR6oWBzaJWJCkVVCCJSJHrEPtvnG05xHi3JdanvdFxpKS43Iq0ezEkHTA2i8hNqEASkSLlcymGxkn/AdCEtJJX9dbQsDdggzVvGZ1G5IZUIIlIkQqNnYsZC8cqd+a0V3Oj40hJ0/1lwAS/fw+ndxmdRuS6VCCJSJHxuXSShsk/AfZRs0Wu4t8Emva1r/8y2dgsIjegAklEikyH2DmYsXKschcSKzUxOo6UVHeNtw/7cOgn/NLUF0lKJhVIIlIkKmecoEHyzwBsrPmMwWmkRPO9A5r+BYAOsbMNDiNybSqQRKRIhJ6ytx4drdyVpIqNjI4jJd2dL4HJTN1z69SKJCWSCiQRuWX21qP/ArBJrUeSH9XqqxVJSjQVSCJyyzrE/hszVo5U6UZSxYZGx5HS4s6XsGJvRaqWdtDoNCJ5qEASkVtSJeMYDc780XqkO9ekIKrV52C1+wC1IknJ42x0ABEp3UJj52DCxuEqd5FcsYHRcYrM+ysPGR2hwMbdW9/oCAW2ucYwGiT/l3rn1lIt7WCZ+gxJ6aYWJBEptMoZJ2hwZiUAm4OfNjiNlEbnPUPUiiQlkgokESm09qc+w4SNo5W76i9/KbTNNYZhxexoRRIpCVQgiUiheGeeomFyJACbg4cZnEZKs/OeIRzyvReA9qfmGZxGxE4FkogUSrtTn2HGwgmfDho1W27ZluChANxxdjVVMo4ZnEZEBZKIFEKlrAQaJ/0HgE3qeyRF4KxnXQ5X7Y4Jm1qRpERQgSQiBdb21HycbLnEeLfltFcLo+NIGbG5hv1SbYPk/+J9KdbgNFLeqUASkQKpkJVM08QfAN25JkUruWIDjlXughkr7U99ZnQcKedUIIlIgbSN+wJnWzZxXi055dXa6DhSxmwOfgqARsn/oVLmaYPTSHmmAklE8s0j+xzNEpcCsKnGMDCZDE4kZU1CpWac9G6Pk81Cu7j5RseRckwFkojkW+v4hbhYszhdsQkxPqFGx5Ey6nIrUpPEH6iQlWxwGimvVCCJSL645V6kRcISALbUGKrWIyk2cd5tOOXVCmdbDm3jvjA6jpRTKpBEJF9anP4GN0s6ZzzrcqxKV6PjSBm3uYa9FalZ4ne451wwNoyUSyqQROSmnC2XaBX/FQBbajwJJv3qkOIV4xNKYoVGuFgzaRW/yOg4Ug7pt5yI3FSzxO/wzL3ABfcaHPINMzqOlAcmE1uCnwSg5emvcc1NMzaPlDsqkETkhpys2bSJWwDA1uqDsZmcDU4k5cWRKndx1iMEd8tFmicsNTqOlDOGF0gzZswgJCQEd3d3QkND2bJly3X33bdvH3379iUkJASTyURERMRV+0yePJl27dpRqVIl/Pz86NOnDwcP5p0d+q677sJkMuVZnnvuuaJ+aSJlQqOkFVTKTiLNtRr7/R4wOo6UJyYz22oMAaB1/AKcLJkGB5LyxNACafHixYSHhzNx4kS2b99OixYt6NGjB0lJSdfcPyMjgzp16jBlyhQCAgKuuc/atWsZMWIEmzZtYuXKleTk5HDfffeRnp6eZ7/hw4dz+vRpx/LOO+8U+esTKe1MtlzHWDTbgp7AYnY1OJGUNwd8e5LqFkCFnHM0SVpudBwpRwwtkKZNm8bw4cMZOnQojRs3ZtasWXh6ejJ37txr7t+uXTveffddBgwYgJub2zX3iYyM5Mknn6RJkya0aNGCzz77jJiYGKKjo/Ps5+npSUBAgGPx8vIq8tcnUtrVPxOFT+YpLjl7syfgEaPjSDlkNTuzrfpfAfso7mZrrsGJpLwwrEDKzs4mOjqasLA/O3yazWbCwsLYuHFjkZ0nJSUFgCpVquTZvmDBAnx9fWnatCkTJkwgIyPjhsfJysoiNTU1zyJSptlstPtjPqztQQPJdfIwNo+UW3v9HiLdpQreWfE0OPNfo+NIOWFYgXTmzBksFgv+/v55tvv7+5OQkFAk57BarYwdO5bOnTvTtGlTx/bHH3+cL7/8kl9++YUJEybwxRdf8MQTT9zwWJMnT8bb29uxBAcHF0lGkZIq5PwGqmUcIdvsya7AfkbHkXLM4uTO9qDHAexFu81qbCApFwpcIE2cOJGTJ08WR5YiN2LECPbu3cuiRXnH0HjmmWfo0aMHzZo1Y9CgQXz++ed89913HD169LrHmjBhAikpKY4lNja2uOOLGOpy36PdAY+S5axL0GKs3QF9yXSqSNVLx6l7bq3RcaQcKHCB9P3331O3bl3uueceFi5cSFZWVqFO7Ovri5OTE4mJiXm2JyYmXrcDdkGMHDmS5cuX88svv1CjRo0b7hsaap9T6siRI9fdx83NDS8vrzyLSFkVlLqLGqk7yDW5OP5yFzFStnNFR0tmu1PzwWYzOJGUdQUukHbu3MnWrVtp0qQJY8aMISAggOeff56tW7cW6Diurq60adOGqKgoxzar1UpUVBQdO3YsaCwHm83GyJEj+e6771i9ejW1a9e+6XN27twJQGBgYKHPK1KWXO57tN/vAdLdqhkbRuQPOwP7k2t2IzBtHzVStxsdR8q4QvVBatWqFdOnTyc+Pp45c+Zw6tQpOnfuTPPmzfnggw8cHaNvJjw8nNmzZzN//nz279/P888/T3p6OkOHDgVg8ODBTJgwwbF/dnY2O3fuZOfOnWRnZxMXF8fOnTvztPyMGDGCL7/8koULF1KpUiUSEhJISEjg0qVLABw9epTXX3+d6OhoTpw4wQ8//MDgwYO58847ad68eWHeDpEypWr6Eeqc/w0bJsfdQyIlQYZrVfb6PQj8WcSLFJdb6qRts9nIyckhOzsbm81G5cqV+eijjwgODmbx4sU3fX7//v2ZOnUqr776Ki1btmTnzp1ERkY6Om7HxMRw+vRpx/7x8fG0atWKVq1acfr0aaZOnUqrVq14+umnHfvMnDmTlJQU7rrrLgIDAx3L5Tyurq6sWrWK++67j4YNG/LCCy/Qt29ffvzxx1t5K0TKjHZxnwNwuOo9XPCoaXAakbyiqz+BFSdCLmyiWtrBmz9BpJAKNWdAdHQ08+bN46uvvsLNzY3BgwczY8YM6tWrB8CHH37I6NGj6d+//02PNXLkSEaOHHnNx9asWZPn55CQEGw3ue58s8eDg4NZu1Yd/ESuxSszngbJ9tuot9YYbHAakaululfnoG8Yjc78TNu4z/mpwZtGR5IyqsAtSM2aNaNDhw4cP36cOXPmEBsby5QpUxzFEcDAgQNJTk4u0qAiUvxaxy/AjIWTPqEkVWxkdByRa9r2R/Fe/8wqvC+dMjiNlFUFLpAee+wxTpw4wX/+8x/69OmDk5PTVfv4+vpitWqcCpHSxCP7HM0SvwdgS/UnjQ0jcgNnKtTneOVOmLHSJv5Lo+NIGVXgAulyX6MrXbp0iddee61IQonI7dfy9GKcrVkkVGzMKe82RscRuaGt1e2T2DZJ/BHP7DMGp5GyqMAF0r/+9S/S0tKu2p6RkcG//vWvIgklIreXiyWDlglLANhaYwiYTAYnErmxOK9WxFdqhrMtm9bxXxkdR8qgQrUgma7xy3PXrl1XzXcmIqVD08RluOemct69JkerdDM6jsjNmUyOVqTmCd/imnv1H+4ityLfd7FVrlwZk8mEyWSifv36eYoki8VCWloazz33XLGEFJHiY7bm0jpuIQDbqj+BzXR1v0KRkuhYla6c9ahN1UvHaZbwHdE1NG6XFJ18F0gRERHYbDaeeuop/vWvf+Ht7e14zNXVlZCQkFsaAVtEjNHgzM94ZSeS7lKV/X69jI4jkn8mM9uqP0GPI6/T+vRX7Azqj8XsanQqKSPyXSANGWJvyqxduzadOnXCxcWl2EKJyG1is9L2j4EhtwcNxGJ2MziQSMEcrNaTzjGzqJidTMPkSPb5P2R0JCkj8tUHKTU11bHeqlUrLl26RGpq6jUXESk9ap9fj2/GMbKcKrAn4FGj44gUmMXsyvaggQC0ifsCbBpiRopGvgqkypUrk5SUBICPjw+VK1e+arm8XURKj7ZxXwCwJ+BRspwrGZxGpHD2+D9CplNFql46QZ1zvxodR8qIfF1iW716teMOtV9++aVYA4nI7RGYupsaqTvINbmwPXCg0XFECi3buSK7A/rSPm4+7eI+51hV3Ykpty5fBVK3bt2uuS4ipdflvkf7q/Ui3a2awWlEbs3OoAG0jl9I0MXdBKXuIt6rhdGRpJQr8DhIkZGR/Pbbb46fZ8yYQcuWLXn88cc5f/58kYYTkeJROeME9c6txYaJ6OpPGB1H5Jalu/qy3+8BANqemm9wGikLClwgvfTSS47O2Hv27CE8PJxevXpx/PhxwsPDizygiBS9NnH2+auOVunGec8QY8OIFJFt1Z/Ahom653+lasZRo+NIKVfgAun48eM0btwYgG+//ZYHH3yQt956ixkzZvDTTz8VeUARKVoVss/QKHkFYP9CESkrLnjU4kjVuwBoHbfA2DBS6hW4QHJ1dSUjIwOAVatWcd999wFQpUoV3eYvUgq0jF+Esy2HuEotOK1+GlLGbKs+GIBGyT9RISvJ4DRSmhW4QOrSpQvh4eG8/vrrbNmyhQcesF/zPXToEDVq1CjygCJSdFxz02ie8C0A26prWgYpexIqNeWUVyucbLm0Pr3I6DhSihW4QProo49wdnZmyZIlzJw5k+rVqwPw008/0bNnzyIPKCJFp2niMtwtaZz1COFYla5GxxEpFpeL/2YJSzWJrRRavqcauaxmzZosX778qu3vv/9+kQQSkeJhtubQOv4rAKKr/xVMBf77SKRUOF65syaxlVtW4AIJwGq1cuTIEZKSkrBa8w7rfueddxZJMBEpWg3P/Eyl7CTSXHw5UE2tvVKGXTGJ7Y6gAVjNmj9UCqbABdKmTZt4/PHHOXnyJDabLc9jJpMJi8VSZOFEpIjYbPZ5qoAdQQM047mUeVdOYvu7/4NGR5JSpsBt7M899xxt27Zl7969nDt3jvPnzzuWc+fOFUdGEblFIec3aFJaKVcsZlfHFDptNYmtFEKBW5AOHz7MkiVLqFevXnHkEZFicHlakT3+fTQprZQbewIeIfTUHKpeOk7t8+s5rhsTpAAK3IIUGhrKkSNHiiOLiBQD/4v7CE7djsXkxI4gTUor5Ue2c0VHi2nbPy4xi+RXgVuQRo0axQsvvEBCQgLNmjXDxSVvx7fmzZsXWTgRuXWXvxgOVOtJmpu/wWlEbq/tgQNoFf8VNVJ3EHBxDwmVmhkdSUqJAhdIffv2BeCpp55ybDOZTNhsNnXSFilhvC/FcsfZ1QBEB2laESl/0t38OFCtJ02SltM27guWN3zH6EhSShS4QDp+/Hhx5BCRYtAmfgEmbByr3JmzFdRvUMqn6OpP0CRpOfXOrsHnUgwXPGoaHUlKgQIXSLVq1SqOHCJSxDyyz9EkyT6oa7SmFZFy7KxnXY5V7kKd87/RJm4BUfUmGB1JSoFCDaX7xRdf0LlzZ4KCgjh58iQAERERfP/990UaTkQKr2XCNzhbs0io2JhTXq2NjiNiqMvTjzROWo5n9lmD00hpUOACaebMmYSHh9OrVy8uXLjg6HPk4+NDREREUecTkUJwtlyixelvgD++GEwmgxOJGCvOqxWnKzbB2ZZNi9NfGx1HSoECF0gffvghs2fP5v/+7/9wcnJybG/bti179uwp0nAiUjhNE3/AIzeFC+7VOVK1u9FxRIxnMrGt+mAAWiYswcWSYXAgKekK1Um7VatWV213c3MjPT29SEKJSOGZbLm0jl8IQHTQIGwmp5s8Q8qK91ceMjpCiXa0ajfOuwdTOTOWJok/sDNogNGRpAQrcAtS7dq12blz51XbIyMjadSoUVFkEpFbcMeZ1XhnxZPh7MPvfpp/SuQym8mJ7UGDAGgdvxCTLdfgRFKSFbgFKTw8nBEjRpCZmYnNZmPLli189dVXTJ48mX//+9/FkVFE8stmc0wrsjPwMXKd3A0OJFKy7PN7gI6xn+CddZr6Z1ZxsFpPoyNJCVXgAunpp5/Gw8ODV155hYyMDB5//HGCgoL44IMPGDBAzZUiRgpO2Yp/+kFyzG7sCuxndByREsfi5M7OwMfoFPMJbeO+4KBvD93EINdUqNv8Bw0axOHDh0lLSyMhIYFTp04xbNiwos4mIgV0eVqRvf4Pk+niY2wYkRJqV8BfyDG745d+iJopW4yOIyVUoQqkyzw9PfHz8yuqLCJyC6qlHSTkwiasmNke9LjRcURKrEwXH/b6PwxA21OfG5xGSqp8XWJr1aoVpnw2QW7fvv2WAolI4bT5o/XosO89pLpXNziNSMm2PehxWpxeQq2ULfilHSCpYkOjI0kJk68WpD59+vDwww/z8MMP06NHD44ePYqbmxt33XUXd911F+7u7hw9epQePXoUd14RuYZKmadpcGYV8OeIwSJyfanuQRz0DQOgTdyXBqeRkihfBdLEiRMdS3JyMqNHj2bjxo1MmzaNadOmsWHDBsaOHUtiYmKBA8yYMYOQkBDc3d0JDQ1ly5brXw/et28fffv2JSQkBJPJdN2Ru292zMzMTEaMGEHVqlWpWLEiffv2LVR2kZKidfxCzFiI8W5HUkUNtyGSH5fnKKx/ZhVemXEGp5GSpsB9kL755hsGDx581fYnnniCb7/9tkDHWrx4MeHh4UycOJHt27fTokULevToQVJS0jX3z8jIoE6dOkyZMoWAgIBCH3PcuHH8+OOPfPPNN6xdu5b4+HgeffTRAmUXKSnccy7QLHEZAFurX/1vU0SuLbliA074dMCMhdbxXxkdR0qYAhdIHh4erF+//qrt69evx929YGOuTJs2jeHDhzN06FAaN27MrFmz8PT0ZO7cudfcv127drz77rsMGDAANze3Qh0zJSWFOXPmMG3aNO6++27atGnDvHnz2LBhA5s2bSpQfpGSoHnCt7hYM0mqcAcxPqFGxxEpVS5fkm6auAz3nAvGhpESpcDjII0dO5bnn3+e7du30759ewA2b97M3Llz+ec//5nv42RnZxMdHc2ECRMc28xmM2FhYWzcuLGgsfJ9zOjoaHJycggLC3Ps07BhQ2rWrMnGjRvp0KHDNY+dlZVFVlaW4+fU1NRCZRQpSk6WTFqdXgxAdNATGs9FpIBivduRWKEB/ukHaXH6GzbXHG50JCkhCtyCNH78eObPn090dDSjR49m9OjRbN++nXnz5jF+/Ph8H+fMmTNYLBb8/f3zbPf39ychIaGgsfJ9zISEBFxdXfHx8SnQeSdPnoy3t7djCQ4OLlRGkaLUJOk/eOacJ9UtgEO+9xkdR6T0+d9JbE9/jbMl0+BAUlIUuAUJ4LHHHuOxxx4r6iwl2oQJEwgPD3f8nJqaqiJJDGWyWWgTb7/7JjpoEFZzof45i5R7h33vJuVkEN5Z8TRO+pHdGoVeuMWBIm+Fr68vTk5OV909lpiYeN0O2EVxzICAALKzs7lw4UKBzuvm5oaXl1eeRcRI9c7+gk/mKS45ezsGvRORgrOZnImubp/Etm3cl5rEVgADCyRXV1fatGlDVFSUY5vVaiUqKoqOHTsW2zHbtGmDi4tLnn0OHjxITExMoc8rctvZbLT7YwTgXYH9yHXyMDiQSOm2z+8hMpx98M6K544zq42OIyWAoW3y4eHhDBkyhLZt29K+fXsiIiJIT09n6NChAAwePJjq1aszefJkwN4J+/fff3esx8XFsXPnTipWrEi9evXydUxvb2+GDRtGeHg4VapUwcvLi1GjRtGxY8frdtAWKWlqpETjn76fHLMbOwPL1+VukeKQe3kS29hPaRv3OYd879VND+WcoQVS//79SU5O5tVXXyUhIYGWLVsSGRnp6GQdExOD2fxnI1d8fDytWrVy/Dx16lSmTp1Kt27dWLNmTb6OCfD+++9jNpvp27cvWVlZ9OjRg48//vj2vGiRItAubj5g/6v3kktlg9OIlA27AvvRLu5z/NMPUjNli4bNKOdMNpvNVpAn/PLLL3Tv3r248pQaqampeHt7k5KSov5Iclt9sWw5f905CCtmPmvzLSnuNYyOJFJm3HVsKq1OL+akd3uWNp1xzX3G3Vv/NqeSopTf7+8C90Hq2bMndevW5Y033iA2NvaWQopIwbX9Y1LaQ75hKo5Eilh00CCsOP0xie1+o+OIgQpcIMXFxTFy5EiWLFlCnTp16NGjB19//TXZ2dnFkU9E/tf5kzRIXgngGLtFRIrORfdADla7F/jzjxEpnwpcIPn6+jJu3Dh27tzJ5s2bqV+/Pn/7298ICgpi9OjR7Nq1qzhyigjAxo8wY+GkTyjJFRsYnUakTLr8x8cdZ6LwvnTK4DRilFu6zb9169ZMmDCBkSNHkpaWxty5c2nTpg1du3Zl3759RZVRRADSz8B2+1+0W6sPMTiMSNl1psIdHPfpiBmrYzBWKX8KVSDl5OSwZMkSevXqRa1atfj555/56KOPSExM5MiRI9SqVYt+/TQSqUiR2vwJ5F4ioWIjYr3bGp1GpEzbWsP+R0iTpOV4Zp81OI0YocAF0qhRowgMDOTZZ5+lfv367Nixg40bN/L0009ToUIFQkJCmDp1KgcOHCiOvCLlU1YabPkUgK3Vn9T4LCLFLM6rNacrNsXZmkXLPyaElvKlwAXS77//zocffkh8fDwRERE0bdr0qn18fX355ZdfiiSgiADb50PmBahSl6NVuxmdRqTsM5kcrUgtTn+Da26awYHkditwgTRx4kT69euHm5tbnu25ubmsW7cOAGdnZ7p10y9xkSKRmw0bPrKvdx6DzeRkbB6RcuJolTs56xGCuyWNZgnfGR1HbrMCF0jdu3fn3LlzV21PSUnRAJIixWHPN3AxHioGQIsBRqcRKT9MZscdbW3iF+BkzTI4kNxOBS6QbDYbpmv0fzh79iwVKlQoklAi8gerFdZH2Nc7/g2c3W64u4gUrQPVenLR1Y8KOWdplLTC6DhyG+V7LrZHH30UAJPJxJNPPpnnEpvFYmH37t106tSp6BOKlGcHV8CZQ+DmDW2GGp1GpNyxml2Irj6Iu46/T9u4L9jn/5DRkeQ2yXeB5O3tDdhbkCpVqoSHh4fjMVdXVzp06MDw4cOLPqFIeWWzwW/T7OvtngJ3zfknYoS9/n0IjZ1L5cxY6p1dAzQyOpLcBvkukObNmwdASEgIL774oi6niRS34+sgLhqc3aHD34xOI1Ju5Th5sjPwMTrGzqbdqflge05DbZQDhbqLTcWRyG1wufWo1V+hop+xWUTKuZ2Bj5Fjdsc/fT8cXW10HLkN8tWC1Lp1a6KioqhcuTKtWrW6Zifty7Zv315k4UTKrbjtcGwNmJyg0yij04iUe5kuPuzx70Pr04vgt/eh3j1GR5Jilq8C6eGHH3Z0yu7Tp09x5hER+LP1qFk/qFzL2CwiAkB09UG0SFiC04lfIXYLBLc3OpIUI5PNZrMZHaI0Sk1Nxdvbm5SUFLy81HlWilDyIZjRHrDB3zaBX94Ooe+vPGRMLhEh7PAbNEv6HurfD48vMjqOFEJ+v78LNVmtiBSj9RGADRo8cFVxJCLG2lbjr4AJDv0EifuMjiPFKF+X2CpXrnzDfkf/61qjbItIPl2Ihd1/TIzZNdzYLCJylQsetaBJH9j3nb0vUt9/Gx1Jikm+CqSIiIhijiFijJJ2uequY1NpZc0lxrst3+73gv0lK5+IAF3G2Qukvd9C95ehSh2jE0kxyFeBNGTIkOLOIVLueeScp2niMgC21njS0CwicgOBLaDevXBkJayfDg9GGJ1IikG++iClpqbmWb/RIiKF0zp+IS7WLBIqNiLGW3fHiJRoly+B71wAqaeNzSLFIl8FUuXKlUlKSgLAx8eHypUrX7Vc3i4iBeeWe5EWp78BYHONYRqlV6Skq9UJanYESzZs/MjoNFIM8nWJbfXq1VSpUgWAX375pVgDiZRHLU5/jZslnTOedTlWpavRcUQkP7qEw8J+sG2efb1CVaMTSRHKV4HUrVu3a66LyK1zsWTQOv4rALbUGAomjb4hUirccS8ENIeE3bB5Jtz9itGJpAjle7La/3X+/HnmzJnD/v37AWjcuDFDhw51tDKJSP41S1iKR24K591rcsg3zOg4IpJfJhPc+RJ8/VfY/Kl9WiB3b6NTSREp8J+q69atIyQkhOnTp3P+/HnOnz/P9OnTqV27NuvWrSuOjCJllpM1i7ZxXwKwtcYQbCYngxOJSIE07A3VGkJWCmyZbXQaKUIFLpBGjBhB//79OX78OEuXLmXp0qUcO3aMAQMGMGLEiOLIKFJmNUn8kQo5Z0l19Wd/tfuNjiMiBWU2Q9cX7OsbZ0B2urF5pMgUuEA6cuQIL7zwAk5Of/6l6+TkRHh4OEeOHCnScCJlmdmaS9u4zwHYVmMwVrOLwYlEpFCaPAqVa8Olc/YO21ImFLhAat26taPv0f/av38/LVq0KJJQIuVBw+Sf8M46TbpLFfb6PWR0HBEpLCfnP8dF2jAdcjKNzSNFIl+dtHfv3u1YHz16NGPGjOHIkSN06NABgE2bNjFjxgymTJlSPClFyhiTzUK7U58BEB30BBYnd2MDicitaT4A1rwNqadgxxfQfrjRieQWmWw2m+1mO5nNZkwmEzfb1WQyYbFYiixcSZaamoq3tzcpKSl4eXkZHUcKyai52BokR9Lr0D+55OzNnLY/kOPkaUgOESm4cffWv/YDW2bDihfBOxhG7wAnXTYvifL7/Z2vFqTjx48XWTCR8s5ksxAaOxeA6KBBKo5EyopWT8C6dyElFnZ9Ba0HG51IbkG+CqRatWoVdw6RcuOOs6upeuk4mU6V2BXYz+g4IlJUXDyg02j47//BuqnQYqBakUqxQg0UCfD7778TExNDdnZ2nu0PPaTOpiLXZbMSGjsHgO1BA8l2rmhwIBEpUm2Hwm/vw4WTsPtraDXI6ERSSAUukI4dO8YjjzzCnj178vRLMv0xuWZ56YMkUhj1zq3BN+MoWU4V2Bk0wOg4IlLUXCtA59Gw8lX4dSo072+/y01KnQLf5j9mzBhq165NUlISnp6e7Nu3j3Xr1tG2bVvWrFlTDBFFygibjdDYfwOwI3AAWc6VDA4kIsWi7TDwrArnjsHeJUankUIqcIG0ceNGXnvtNXx9fTGbzZjNZrp06cLkyZMZPXp0cWQUKRPqnFuHX/phss2ebA8aaHQcESkubhWh40j7+rp3waorK6VRgdv9LBYLlSrZ//L19fUlPj6eBg0aUKtWLQ4ePFjkAUXKBJuNDn/0PdoZ9BhZLprQUqS0ys/wIC653RnmHIHH2SOsWDSDg9V63oZk13fdoQnkugrcgtS0aVN27doFQGhoKO+88w7r16/ntddeo06dOoUKMWPGDEJCQnB3dyc0NJQtW7bccP9vvvmGhg0b4u7uTrNmzVixYkWex00m0zWXd99917FPSEjIVY9roEspLrXPr8c/fT85Zne2Bz1udBwRKWY5zhUc/9ZDY+disqkVqbQpcIH0yiuvYLVaAXjttdc4fvw4Xbt2ZcWKFUyfPr3AARYvXkx4eDgTJ05k+/bttGjRgh49epCUlHTN/Tds2MDAgQMZNmwYO3bsoE+fPvTp04e9e/c69jl9+nSeZe7cuZhMJvr27ZvnWK+99lqe/UaNGlXg/CI3ZbPRIfZTAHYF9uOSS2WDA4nI7bAz8DEynSpR9dJx7jgTZXQcKaB8jaR9M+fOnaNy5cqOO9kKIjQ0lHbt2vHRRx8BYLVaCQ4OZtSoUYwfP/6q/fv37096ejrLly93bOvQoQMtW7Zk1qxZ1zxHnz59uHjxIlFRf35AQ0JCGDt2LGPHjs1XzqysLLKyshw/p6amEhwcrJG0S7nbMZJ2nXPreHj/C2SbPZjb9nsVSCLlSGjMbDrFfspZj9p80eorbCanmz+pGOgS25/yO5J2gVuQ/ldsbCyxsbFUqVKlUMVRdnY20dHRhIWF/RnIbCYsLIyNGzde8zkbN27Msz9Ajx49rrt/YmIi//nPfxg2bNhVj02ZMoWqVavSqlUr3n33XXJzc6+bdfLkyXh7ezuW4ODg/LxEKe9sNjrEzAbUeiRSHu0MGqBWpFKqwAVSbm4u//znP/H29iYkJISQkBC8vb155ZVXyMnJKdCxzpw5g8Viwd/fP892f39/EhISrvmchISEAu0/f/58KlWqxKOPPppn++jRo1m0aBG//PILzz77LG+99RZ///vfr5t1woQJpKSkOJbY2Nj8vEQp5+qcW4d/+gGyzZ5sq/5Xo+OIyG2W5VyJ7dXtfZE6xM5WX6RSpMB3sY0aNYqlS5fyzjvv0LFjR8DeqjNp0iTOnj3LzJkzizzkrZg7dy6DBg3C3T3vbOnh4eGO9ebNm+Pq6sqzzz7L5MmTcXNzu+o4bm5u19wucl02Gx3/6Hu0M+gxMl18jM0jIobYETiA1vFfUfXSCRok/5cDfvcbHUnyocAF0sKFC1m0aBH33//n/+DmzZsTHBzMwIEDC1Qg+fr64uTkRGJiYp7tiYmJBAQEXPM5AQEB+d7/119/5eDBgyxevPimWUJDQ8nNzeXEiRM0aNAg369B5HrqnluDX/ohspwqEB2k6QZEyqts54psq/4EXU5+TIfYf3Ow2r3YTBpdu6Qr8CU2Nzc3QkJCrtpeu3ZtXF1dC3QsV1dX2rRpk6fztNVqJSoqytE6daWOHTvm2R9g5cqV19x/zpw5tGnThhYtWtw0y86dOzGbzfj5+RXoNYhck81Kxz/6Hu0M7K/WI5FybmfAY2Q4+1A5M4ZGyZFGx5F8KHCBNHLkSF5//fU8d3RlZWXx5ptvMnLkyAIHCA8PZ/bs2cyfP5/9+/fz/PPPk56eztChQwEYPHgwEyZMcOw/ZswYIiMjee+99zhw4ACTJk1i27ZtV507NTWVb775hqeffvqqc27cuJGIiAh27drFsWPHWLBgAePGjeOJJ56gcmV1opVbV+/sGqplHP6j9UjjHomUdznOFdhWfTAAobH/xmy9/k1BUjLkq43vyg7Oq1atokaNGo6WmV27dpGdnc0999xT4AD9+/cnOTmZV199lYSEBFq2bElkZKSjI3ZMTAxm8591XKdOnVi4cCGvvPIKL7/8MnfccQfLli2jadOmeY67aNEibDYbAwdePaWDm5sbixYtYtKkSWRlZVG7dm3GjRuXp1+SSKHZrI6+RzsCB2jUbBEBYFfgX2gT/yU+mXE0Sl7BPv+HjI4kN5CvcZAut+bkx7x5824pUGmR33EUpGQrjnGQGiT/TK9Dr5DpVJG5bb8ny1mfDxGxax23gG4nIkhxC+Kz1kuwml1uy3k1DtKf8vv9na8WpPJS9IjcKpMtl44x9taj6OpPqDgSkTx2B/SlbdwXeGfF0zTxB3YH9r35k8QQhR4oMjk5md9++43ffvuN5OTkoswkUmo1TlpB5cwYMpx92BE4wOg4IlLC5Dq5s6WG/apM6Kk5OFkyDU4k11PgAik9PZ2nnnqKwMBA7rzzTu68806CgoIYNmwYGRkZxZFRpFRwsmbTIdZ+59rWGk+S41zB4EQiUhLtCXiEVLcAKmYn0yLhW6PjyHUUuEAKDw9n7dq1/Pjjj1y4cIELFy7w/fffs3btWl544YXiyChSKjRNXIZXVgJprtXYFaBmcxG5NovZlU3BwwFod+ozXHLTDU4k11LgAunbb79lzpw53H///Xh5eeHl5UWvXr2YPXs2S5YsKY6MIiWesyWT0Ni5AGyu8RQWJ/ebPENEyrPf/Xpxzr0mnrkXaH36K6PjyDUUuEDKyMi4ai40AD8/P11ik3KrxelvqJBzlhS3IPb6P2x0HBEp4WwmZzbWfBaANnFf4paTYnAiuVKBC6SOHTsyceJEMjP/7Fh26dIl/vWvf1139GuRssw1N412cfMB2BT89G27bVdESrdDvmEkVbgDN0s67eI+NzqOXKHAk8FERETQs2fPqwaKdHd35+effy7ygCIlXev4hXjkpnDOoxb7NQmliOSXycyGms/TZ384LU8vZkfQQNJdfY1OJX8ocIHUrFkzDh8+zIIFCzhw4AAAAwcOZNCgQXh4eBR5QJGSzCP7HG3iFgCwoeZzmoBSRArkeOUuxFdqRtDFPbSPncsvdf9udCT5Q4F+m+fk5NCwYUOWL1/O8OHDiyuTSKkRemourtYMEio24nDVgk+1IyLlnMnE+pp/o9++52mW+B3bqz9OinsNo1MJBeyD5OLikqfvkUh55p15iuZ/jGHyW62RYDIZnEhESqNTPm054dMBJ1suHWM+MTqO/KHAnbRHjBjB22+/TW6uZiKW8q1jzCc42XI56RNKrE97o+OISCn2W62RADRKjqRa2kGD0wgUog/S1q1biYqK4r///S/NmjWjQoW8owUvXbq0yMKJlFS+6YdomGy/KeG3WiMMTiMipV1yxQYc8O1BwzM/0+XkR3zX5EOjI5V7BS6QfHx86NtXowRL+db55MeYsHGwahhJFRsZHUdEyoANtZ7jjrNRhFzYRPCFLWqZNliBC6R58+YVRw6RUqN6ynbqnF+PxeTEhlrPGx1HRMqIFPca7A7oS6vTi+l64iMWtvgMTIWeU15uUb7feavVyttvv03nzp1p164d48eP59KlS8WZTaTksdnocvIjAPb69+GCR02DA4lIWbK5xlNkmz3xT99P/bNRRscp1/JdIL355pu8/PLLVKxYkerVq/PBBx8wYoT6Xkj5csfZ1QRd3EOO2Z3NwU8bHUdEyphLrlXYVv0JwH4p32zVDVFGyXeB9Pnnn/Pxxx/z888/s2zZMn788UcWLFiA1WotznwiJYbZmuNoPYqu/oRGvBWRYrG9+iDSXargk3mKZom68cko+S6QYmJi6NWrl+PnsLAwTCYT8fHxxRJMpKRpkbAEn8xTpLtUYVv1vxodR0TKqBwnTzb90ULdIWY2rrlpBicqn/JdIOXm5uLu7p5nm4uLCzk5OUUeSqSkccu9SGjsHAA21nyWHCdPgxOJSFm21/8RznnUwjP3Au1P6eYoI+T7LjabzcaTTz6Jm5ubY1tmZibPPfdcnrGQNA6SlEXtT83DIzeFsx612ev/kNFxRKSMs5qd+TVkNA/vf4FW8YvYFfAXLroHGh2rXMl3gTRkyJCrtj3xxBNFGkakJPLKjKdl/GIAfg0ZrQlpReS2OFa5KzHebamZso3OJ2cQ2eANoyOVK/n+Ta/xj6S86nTyY5xt2cR4t+V45c5GxxGR8sJk4teQMTy+azCNzvzMjqCBJFZqYnSqckMjUIncgP/F32l0xj6lyK8hYzQhrYjcVkkVG7K/mv0GqTtPfAA2m8GJyg8VSCLXY7PR7fg0APZXu5+kig0NDiQi5dH6Ws+TY3ajRuoO6p5bY3ScckMFksh11D+zkuoXd5FjdteEtCJimDQ3f6KDBgHQ9cSHmK26e/x2UIEkcg1Olky6nrTPpr21+mDS3PwNTiQi5dm26oNJd6lC5cxYWp5ebHScckEFksg1tI3/Eq+sBFJd/YnWoJAiYrAc5wqOluwOsf/GM/uswYnKPhVIIleokJVEu1PzAftt/blO7jd5hohI8fvdrzcJFRvhZkmn08mZRscp81QgiVyh68mPcLFmElepBYd87zU6joiIncnMmtovAtA06Qf80vYbHKhsU4Ek8j8CLu6hUfJPAKytHa7b+kWkRDnt1Zz91XpiwsZdx97Tbf/FSAWSyGU2K3cds9/Wv8+vN4mVGhscSETkar/VGkmO2Z3qF3fR4Mx/jY5TZqlAEvlD46TlBKbtJdvswXrd1i8iJVSamz9bajwJQNcT03G2XDI2UBmlAkkEcMtNpevJjwDYFDycdFdfgxOJiFxfdNAgUtyCqJSdRLtTnxkdp0xSgSQCdIz5BM+c85z1CGFH0ACj44iI3JDFyZ11IWMAaBv3BT6XYgxOVPaoQJJyr1raQVqcXgLAL3Vewmp2MTiRiMjNHananRM+HXC25dD92LvqsF3EVCBJ+Waz0v3Yu5ixcrBqGLE+7Y1OJCKSPyYTv9R5iVyTCyEXNlHv3C9GJypTVCBJudY4eQXVL+4i2+zButpjjY4jIlIgFzxqsu2P0f67HZumDttFSAWSlF+XLtDlhH2+tc3BwzTfmoiUSltrDCXFLRCv7ERCY+cYHafMKBEF0owZMwgJCcHd3Z3Q0FC2bNlyw/2/+eYbGjZsiLu7O82aNWPFihV5Hn/yyScxmUx5lp49e+bZ59y5cwwaNAgvLy98fHwYNmwYaWlpRf7apAT75S0q5JzjnEcttgc9bnQaEZFCyXVyZ03tFwBoE7+AyhknjA1URhheIC1evJjw8HAmTpzI9u3badGiBT169CApKema+2/YsIGBAwcybNgwduzYQZ8+fejTpw979+7Ns1/Pnj05ffq0Y/nqq6/yPD5o0CD27dvHypUrWb58OevWreOZZ54pttcpJUxcNGz5FFDHbBEp/Y5V7caxyl1wsuVy97G31WG7CJhsNmPfxdDQUNq1a8dHH9nHoLFarQQHBzNq1CjGjx9/1f79+/cnPT2d5cuXO7Z16NCBli1bMmvWLMDegnThwgWWLVt2zXPu37+fxo0bs3XrVtq2bQtAZGQkvXr14tSpUwQFBV31nKysLLKyshw/p6amEhwcTEpKCl5eXoV+/WIASy7M7g4Ju9lfrSeR9V83OpGIyC3zzjzF4B0DcLZm8dMdr3HA737HY+PurW9gspIlNTUVb2/vm35/G9qClJ2dTXR0NGFhYY5tZrOZsLAwNm7ceM3nbNy4Mc/+AD169Lhq/zVr1uDn50eDBg14/vnnOXv2bJ5j+Pj4OIojgLCwMMxmM5s3b77meSdPnoy3t7djCQ4OLvDrlRJiyyeQsBvcfVgbMs7oNCIiRSLFvQabazwFQLcT7+Oec8HYQKWcoQXSmTNnsFgs+Pvn7Rzr7+9PQkLCNZ+TkJBw0/179uzJ559/TlRUFG+//TZr167l/vvvx2KxOI7h5+eX5xjOzs5UqVLluuedMGECKSkpjiU2NrbAr1dKgAuxsPpN+/q9r3HJtYqxeUREitC26n/ljGcdPHPO0/XEdKPjlGrORgcoDgMG/DkScrNmzWjevDl169ZlzZo13HPPPYU6ppubG25ubkUVUYxgs8GKlyAnHWp2hFZ/hagjRqcSESkyVrMLq+q+zIA9T9M06Uf2V+vFKZ+2N3+iXMXQFiRfX1+cnJxITEzMsz0xMZGAgIBrPicgIKBA+wPUqVMHX19fjhw54jjGlZ3Ac3NzOXfu3A2PI6XcgeVw6Ccwu0DvCDAbfo+CiEiRO+3Vgl0BfQEIOzoZJ2vWTZ4h12JoC5Krqytt2rQhKiqKPn36APZO2lFRUYwcOfKaz+nYsSNRUVGMHTvWsW3lypV07Njxuuc5deoUZ8+eJTAw0HGMCxcuEB0dTZs2bQBYvXo1VquV0NDQonlxUrJkpsKKv9vXO48Bv4bG5hERKUbra42g7tm1VM6MoX3sPN5f+ZzRkQrM6I7lhv8JHR4ezuzZs5k/fz779+/n+eefJz09naFDhwIwePBgJkyY4Nh/zJgxREZG8t5773HgwAEmTZrEtm3bHAVVWloaL730Eps2beLEiRNERUXx8MMPU69ePXr06AFAo0aN6NmzJ8OHD2fLli2sX7+ekSNHMmDAgGvewSZlwKqJcDEeKteGO180Oo2ISLHKcq7Emjr233Xt4uZTJeOYwYlKH8MLpP79+zN16lReffVVWrZsyc6dO4mMjHR0xI6JieH06dOO/Tt16sTChQv59NNPadGiBUuWLGHZsmU0bdoUACcnJ3bv3s1DDz1E/fr1GTZsGG3atOHXX3/N04dowYIFNGzYkHvuuYdevXrRpUsXPv3009v74uX2OL4Ots21rz/0Ibh4GJtHROQ2OFz1bo5W7oqTLZd7j7yByWYxOlKpYvg4SKVVfsdREINlp8PMTnD+BLR9Cnq/n+fh91ceMiaXiMhtUDErgcE7BuBmSWdNyDh2VC89swYU1yW2UjEOkkixW/2mvTjyqgFh/zI6jYjIbZXmFsCvIaMB6BLzMT6XYgxOVHqoQJKyK3YLbPrYvv7gB+Culj4RKX/2+D/CSe/2OFuzuPfIG2CzGh2pVFCBJGVTTiZ8PwKwQYvH4Y6wmz5FRKRMMplYVe//yDZ7UCN1By1Pf2N0olJBBZKUTWvfhjOHoIIf9HjT6DQiIoZKdQ/i15BRAHQ5+RHel04ZnKjkU4EkZU/sFlgfYV/vPQ08NZ2IiMjugL7EerXBxZqpS235oAJJypbsdPjuWfs//Ob9odGDRicSESkZTGZW3vEKOWZ3glOjaXn6a6MTlWgqkKRsWfkqnDsGXtXh/neMTiMiUqKkuNdwXGrrevIjqmQcNzhRyaUCScqOI6tg67/t6w/PAA8fQ+OIiJREuwL6ccKnA87WLHoeehWzNcfoSCWSCiQpGy6dh+//mL+v/bNQt7uxeURESiqTif/We5VLzt74px+gQ+y/jU5UIqlAkrJhxUtw8TRUrQdhk4xOIyJSoqW7VSOqrn2e03anPiMwdZfBiUoeFUhS+u1ZAnu+AZMTPPIpuHoanUhEpMQ77HsPv1frhRkrPQ9PxMWSYXSkEkUFkpRu547D8nH29TtfghptjM0jIlKK/FLnJVLdAvDJjKPbsWlGxylRVCBJ6WXJgW+fhqxUCO5gL5BERCTfsp0r8vMdk7BholnS99Q/s9LoSCWGCiQpvdZMhrht4O4NfWeDk7PRiURESp1T3m3YUmMoAGFH3sQrM87gRCWDCiQpnY6thV//aA5+cDr41DQ2j4hIKbax5nDiKzXHzZJOr4P/h9maa3Qkw6lAktIn/SwsfQawQesh0KSP0YlEREo1m8mZFfXfINOpEoFp++gUM9PoSIZTgSSli9UK3/8N0hLAtwH0nGJ0IhGRMuGieyAr73gFgHZxn1Pz/CaDExlLBZKULhs+gEOR4OQGf5mjW/pFRIrQkap3szPgLwD0PDyRCtlnDE5kHBVIUnoc/xWiXrOv93oXApoZm0dEpAxaV3ssyZ53UCHnHL0OvozJVj77I6lAktLhYgIseQpsVmjxOLQebHQiEZEyyWJ24z8N3iLLqQI1UnfQ5cQMoyMZQgWSlHyWXHtxlJ4Efk3ggffAZDI6lYhImXXeM4T/1nsVgLbxX1LvzGqDE91+KpCk5Fv9GpxcD66V4LHP1e9IROQ2OOJ7N9uCngDgviOvUTnjhLGBbjMVSFKy7f8R1n9gX3/4Q/CtZ2weEZFy5LeQEZzyaoWbJZ0HD/y9XM3XpgJJSq7EfbD0Wft66PPQ5BFj84iIlDM2kzP/afAWaS6+VL10nLAjb4LNZnSs20IFkpRMGefgq4GQkw6174T7Xjc6kYhIuZTh6st/Gk7GYnKi4Zn/0i5uvtGRbgsVSFLyWHLhmyfhwknwqQX95oOTi9GpRETKrXivlqyp/SIAnU9+TJ1z6wxOVPxUIEnJs/KfcHwtuFSAgV+BZxWjE4mIlHu7A//CroC+mLBx/8F/UjXjqNGRipUKJClZdi6ETR/b1x+ZCf5NjM0jIiIOa2q/SKxXG1ytGTy0/wXccy4YHanYqECSkuPEevhxjH292z+g8cPG5hERkTysZmeWN5xCilsQPplx9D4wHrO1bI60rQJJSoYzh2HR42DJhkYPQrfxRicSEZFryHTx4fvG08g2exKcGs3dx6aUyTvbVCCJ8dLPwIK/QOYFqN4WHp0NZn00RURKqrOedVnR4A2smGmW+D3tT80zOlKR07eQGCvnkv12/vMn7HesDVwELh5GpxIRkZs4XqUra+q8AEDnmJk0TFphcKKipQJJjGO1wnfPwakt4O4Ng76BitWMTiUiIvm0K/AxtlX/KwD3HXmdGhe2GZyo6KhAEmPYbPbb+X9fBmYX6L8AqjUwOpWIiBTQr7VGctD3XpxsuTx44KUyc/u/CiQxxm/vw8aP7OsPfwS1uxqbR0RECsdk5uc7JnLKqxXuljT67BtDxawEo1PdMhVIcvtFfwZR/7Kv3/cmtBhgaBwREbk1FrMbPzR8l7MeIXhlJ9J330g8ss8ZHeuWqECS2+v372H5OPt6l3DoNNLYPCIiUiSyXLxZ2uRDUt0CqHLpJI/8PhrX3DSjYxWaCiS5fY6tgW+fBpsVWg+Be141OpGIiBShNLcAvm0yg3SXKvinH+Th/eE4WzKNjlUoJaJAmjFjBiEhIbi7uxMaGsqWLVtuuP8333xDw4YNcXd3p1mzZqxY8eethTk5OfzjH/+gWbNmVKhQgaCgIAYPHkx8fHyeY4SEhGAymfIsU6ZMKZbXJ0DMJvjq8kCQD0Hv98FkMjqViIgUsQseNVna5EMynSpSI3UHvQ/+A7M1x+hYBWZ4gbR48WLCw8OZOHEi27dvp0WLFvTo0YOkpKRr7r9hwwYGDhzIsGHD2LFjB3369KFPnz7s3bsXgIyMDLZv384///lPtm/fztKlSzl48CAPPfTQVcd67bXXOH36tGMZNWpUsb7Wcit2K3zZF3LSoc5d0PffYHYyOpWIiBSTMxXq833jCHLMbtQ+v4Feh14pdVOSmGw2Y8cHDw0NpV27dnz0kf2OJqvVSnBwMKNGjWL8+Kunm+jfvz/p6eksX77csa1Dhw60bNmSWbNmXfMcW7dupX379pw8eZKaNWsC9haksWPHMnbs2ELlTk1Nxdvbm5SUFLy8vAp1jHIhLho+7wNZqRDSFR7/Glw9jU7l8P7KQ0ZHEBEps2qd38hD+1/A2ZbDoar38FP9N7CanfP13HH31i+WTPn9/ja0BSk7O5vo6GjCwsIc28xmM2FhYWzcuPGaz9m4cWOe/QF69Ohx3f0BUlJSMJlM+Pj45Nk+ZcoUqlatSqtWrXj33XfJzb1+dZuVlUVqamqeRW4ifgd88Yi9OKrVGR5fXKKKIxERKV4nK3fkx4bvkmtyof7ZKO4vRS1JhhZIZ86cwWKx4O/vn2e7v78/CQnXHkMhISGhQPtnZmbyj3/8g4EDB+apFEePHs2iRYv45ZdfePbZZ3nrrbf4+9//ft2skydPxtvb27EEBwfn92WWT6d321uOMlMguMMfxVEFo1OJiMhtdqJK51JZJBneB6k45eTk8Nhjj2Gz2Zg5c2aex8LDw7nrrrto3rw5zz33HO+99x4ffvghWVlZ1zzWhAkTSElJcSyxsbG34yWUTrFbYH5v++SzNdrZpxBxq2R0KhERMciJKp1Z3vCdK4qkkt1x29ACydfXFycnJxITE/NsT0xMJCAg4JrPCQgIyNf+l4ujkydPsnLlypv2EwoNDSU3N5cTJ05c83E3Nze8vLzyLHINx9bmbTl64ltw13slIlLeHa/SJU+R9NCBF0v0EACGFkiurq60adOGqKgoxzar1UpUVBQdO3a85nM6duyYZ3+AlStX5tn/cnF0+PBhVq1aRdWqVW+aZefOnZjNZvz8/Ar5aoSDkbCg3x93q3WHvy61T0IrIiKCvUj6odF7jrvbSvJgkoZfYgsPD2f27NnMnz+f/fv38/zzz5Oens7QoUMBGDx4MBMmTHDsP2bMGCIjI3nvvfc4cOAAkyZNYtu2bYwcaR+ROScnh7/85S9s27aNBQsWYLFYSEhIICEhgezsbMDe0TsiIoJdu3Zx7NgxFixYwLhx43jiiSeoXLny7X8TyoK938LiQWDJggYPqM+RiIhc08nKHVna5CPHOEl/2fscHjnnjY51lfzda1eM+vfvT3JyMq+++ioJCQm0bNmSyMhIR0fsmJgYzOY/67hOnTqxcOFCXnnlFV5++WXuuOMOli1bRtOmTQGIi4vjhx9+AKBly5Z5zvXLL79w11134ebmxqJFi5g0aRJZWVnUrl2bcePGER4efntedFmz+VP46e+ADZr1gz4zwcnF6FQiIlJCxXu1ZEnTWTz6+yj80w/y2J7hLG3yERfdrt29xgiGj4NUWmkcJMBqhVUTYcN0+89tn4JeU0vVIJAaB0lExDg+l07Sd+8IvLITSXOtxrJG75NcsQFQzsdBklIsNwuWPv1ncXTPq/DAtFJVHImIiLEueNRicfN/c9ajNhWzk3lszzOEnF9vdCxABZIUxqUL8MWj9n5HZmd45BPo+oLmVhMRkQJLcwtgcfM5xHi3xdWawcO/v0CzhKVGx1KBJAV05gj8OwxO/gaulWDQEmgxwOhUIiJSimU5V+K7xtP5vdoDmLEQdnQyrJxo78phEBVIkn9HVsHsu+HsYfCqDk9FQt3uRqcSEZEywGp24ec7JrIx+Bn7hvURsGmGYXkMv4tNrlbiOg7bbLSOX0jXE9MxYyW+UnN+bPAOGXtcYU8JyyoiIqWXycSmmsNJdQ+kx6UV0GaoYVFUIMkNOVkyuefYFJok/QeAvX4PsbruP7CYXQ1OJiIiZdXvfr3pcc8YQ2/8UYEk1+VzKYYHDo7HL/0wVpxYW3ssOwP7qzO2iIgUP4PvilaBJNd0x5lV3HvkDdws6WS4VGZF/TeI9WlvdCwREZHbQgWS5OFkzabriem0Or0YgFNerVhR/03S3aoZnExEROT2UYEkDpUzTnD/oVfxT98PwJbqQ9hQ6zlsJn1MRESkfNE3n4DNRouEb+h6Yjou1iwuOXvz8x0TOV6lq9HJREREDKECqZyrkH2Gew+/Ru0LGwE44dOB/9Z7VZfURESkXFOBVF7ZbDQ48zPdj03FIzeFXLMbv9Yaxc7AfmDS+KEiIlK+qUAqhypmJXDP0SnU+WNCwMQKDYis/xrnPOsYnExERKRkUIFUntistEhYQpcTM3C1ZpBrcmFL8FNsrT4Eq9nF6HQiIiIlhgqkcsIv7QDdj71D0MU9AMRXas7Keq9wzrO2wclERERKHhVIZZx7zgU6n5xJs8TvMGEj2+zJbyEj2BXwF/U1EhERuQ4VSGWUyZZLs4Tv6RQzE4/cFAAO+PZgXcho0t38DE4nIiJSsqlAKmtsNuqc/5UuJz6i6qXjACR71uOXOi8S593G4HAiIiKlgwqkMiTg4h66nviQGqk7ALjk7M2m4OHsCuyr0bBFREQKQN+aZYBv+iE6xMzmjnNrAMg1u7EjsD9bazxJlnMlY8OJiIiUQiqQSjG/tAOExv6beufWAmDFzO9+vdlY8xnS3PwNTiciIlJ6qUAqhQIu7qF97Dzqnv8VABsmDvrey+bgYRrsUUREpAioQColTDYLdc6to23clwRd3A3YW4wOVruPLTWe0nhGIiIiRUgFUgnnmptGo+QVtIpfROXMWAByTS4cqNaTbdUHc94zxNiAIiIiZZAKpBLKN/0QLU5/S8PkSFytGQBkOnuxK6AvuwIfI93V1+CEIiIiZZcKpJLEZoPdX9N/90eOy2gAZz1C2BXYj31+D5Lr5GFgQBERkfJBBVJJYjLB1tkEXdyNxeTE0Srd2RXQl1PebeyPiYiIyG2hAqmk6TyW9Rt/Za//w2ToMpqIiIghVCCVNI16s+VUfaNTiIiIlGuazl1ERETkCiqQRERERK6gAklERETkCiqQRERERK6gAklERETkCiqQRERERK6gAklERETkCiqQRERERK6gAklERETkCiWiQJoxYwYhISG4u7sTGhrKli1bbrj/N998Q8OGDXF3d6dZs2asWLEiz+M2m41XX32VwMBAPDw8CAsL4/Dhw3n2OXfuHIMGDcLLywsfHx+GDRtGWlpakb82ERERKX0ML5AWL15MeHg4EydOZPv27bRo0YIePXqQlJR0zf03bNjAwIEDGTZsGDt27KBPnz706dOHvXv3OvZ55513mD59OrNmzWLz5s1UqFCBHj16kJmZ6dhn0KBB7Nu3j5UrV7J8+XLWrVvHM888U+yvV0REREo+k81msxkZIDQ0lHbt2vHRRx8BYLVaCQ4OZtSoUYwfP/6q/fv37096ejrLly93bOvQoQMtW7Zk1qxZ2Gw2goKCeOGFF3jxxRcBSElJwd/fn88++4wBAwawf/9+GjduzNatW2nbti0AkZGR9OrVi1OnThEUFHTT3KmpqXh7e5OSkoKXl1dRvBUO7688VKTHExERKW3G3Vs885Lm9/vb0Mlqs7OziY6OZsKECY5tZrOZsLAwNm7ceM3nbNy4kfDw8DzbevTowbJlywA4fvw4CQkJhIWFOR739vYmNDSUjRs3MmDAADZu3IiPj4+jOAIICwvDbDazefNmHnnkkavOm5WVRVZWluPnlJQUwP5GF7XMdF3qExGR8q04vl//97g3ax8ytEA6c+YMFosFf3//PNv9/f05cODANZ+TkJBwzf0TEhIcj1/edqN9/Pz88jzu7OxMlSpVHPtcafLkyfzrX/+6antwcPD1Xp6IiIgU0svFfPyLFy/i7e193ccNLZBKkwkTJuRpubJarZw7d46qVatiMpkMTFY4qampBAcHExsbW+SXCMs7vbfFR+9t8dF7W3z03hafwry3NpuNixcv3rQ7jaEFkq+vL05OTiQmJubZnpiYSEBAwDWfExAQcMP9L/83MTGRwMDAPPu0bNnSsc+VncBzc3M5d+7cdc/r5uaGm5tbnm0+Pj43foGlgJeXl/7BFhO9t8VH723x0XtbfPTeFp+Cvrc3ajm6zNC72FxdXWnTpg1RUVGObVarlaioKDp27HjN53Ts2DHP/gArV6507F+7dm0CAgLy7JOamsrmzZsd+3Ts2JELFy4QHR3t2Gf16tVYrVZCQ0OL7PWJiIhI6WT4Jbbw8HCGDBlC27Ztad++PREREaSnpzN06FAABg8eTPXq1Zk8eTIAY8aMoVu3brz33ns88MADLFq0iG3btvHpp58CYDKZGDt2LG+88QZ33HEHtWvX5p///CdBQUH06dMHgEaNGtGzZ0+GDx/OrFmzyMnJYeTIkQwYMCBfd7CJiIhI2WZ4gdS/f3+Sk5N59dVXSUhIoGXLlkRGRjo6WcfExGA2/9nQ1alTJxYuXMgrr7zCyy+/zB133MGyZcto2rSpY5+///3vpKen88wzz3DhwgW6dOlCZGQk7u7ujn0WLFjAyJEjueeeezCbzfTt25fp06ffvhduMDc3NyZOnHjVZUO5dXpvi4/e2+Kj97b46L0tPsX53ho+DpKIiIhISWP4SNoiIiIiJY0KJBEREZErqEASERERuYIKJBEREZErqEASTpw4wbBhw6hduzYeHh7UrVuXiRMnkp2dbXS0UmnGjBmEhITg7u5OaGgoW7ZsMTpSqTd58mTatWtHpUqV8PPzo0+fPhw8eNDoWGXOlClTHEOlyK2Li4vjiSeeoGrVqnh4eNCsWTO2bdtmdKxSz2Kx8M9//jPPd9brr79+07nVCsrw2/zFeAcOHMBqtfLJJ59Qr1499u7dy/Dhw0lPT2fq1KlGxytVFi9eTHh4OLNmzSI0NJSIiAh69OjBwYMHr5r/T/Jv7dq1jBgxgnbt2pGbm8vLL7/Mfffdx++//06FChWMjlcmbN26lU8++YTmzZsbHaVMOH/+PJ07d6Z79+789NNPVKtWjcOHD1O5cmWjo5V6b7/9NjNnzmT+/Pk0adKEbdu2MXToULy9vRk9enSRnUe3+cs1vfvuu8ycOZNjx44ZHaVUCQ0NpV27dnz00UeAfWT44OBgRo0axfjx4w1OV3YkJyfj5+fH2rVrufPOO42OU+qlpaXRunVrPv74Y9544w1atmxJRESE0bFKtfHjx7N+/Xp+/fVXo6OUOb1798bf3585c+Y4tvXt2xcPDw++/PLLIjuPLrHJNaWkpFClShWjY5Qq2dnZREdHExYW5thmNpsJCwtj48aNBiYre1JSUgD0GS0iI0aM4IEHHsjz2ZVb88MPP9C2bVv69euHn58frVq1Yvbs2UbHKhM6depEVFQUhw4dAmDXrl389ttv3H///UV6Hl1ik6scOXKEDz/8UJfXCujMmTNYLBbHKPCX+fv7c+DAAYNSlT1Wq5WxY8fSuXPnPCPoS+EsWrSI7du3s3XrVqOjlCnHjh1j5syZhIeH8/LLL7N161ZGjx6Nq6srQ4YMMTpeqTZ+/HhSU1Np2LAhTk5OWCwW3nzzTQYNGlSk51ELUhk2fvx4TCbTDZcrv7jj4uLo2bMn/fr1Y/jw4QYlF7m+ESNGsHfvXhYtWmR0lFIvNjaWMWPGsGDBgjxTMcmts1qttG7dmrfeeotWrVrxzDPPOOb/lFvz9ddfs2DBAhYuXMj27duZP38+U6dOZf78+UV6HrUglWEvvPACTz755A33qVOnjmM9Pj6e7t2706lTJ8fkv5J/vr6+ODk5kZiYmGd7YmIiAQEBBqUqW0aOHMny5ctZt24dNWrUMDpOqRcdHU1SUhKtW7d2bLNYLKxbt46PPvqIrKwsnJycDExYegUGBtK4ceM82xo1asS3335rUKKy46WXXmL8+PEMGDAAgGbNmnHy5EkmT55cpK1zKpDKsGrVqlGtWrV87RsXF0f37t1p06YN8+bNyzNBsOSPq6srbdq0ISoqij59+gD2vyKjoqIYOXKkseFKOZvNxqhRo/juu+9Ys2YNtWvXNjpSmXDPPfewZ8+ePNuGDh1Kw4YN+cc//qHi6BZ07tz5qqEoDh06RK1atQxKVHZkZGRc9R3l5OSE1Wot0vOoQBLi4uK46667qFWrFlOnTiU5OdnxmFo+CiY8PJwhQ4bQtm1b2rdvT0REBOnp6QwdOtToaKXaiBEjWLhwId9//z2VKlUiISEBAG9vbzw8PAxOV3pVqlTpqn5cFSpUoGrVqurfdYvGjRtHp06deOutt3jsscfYsmULn376qVrni8CDDz7Im2++Sc2aNWnSpAk7duxg2rRpPPXUU0V7IpuUe/PmzbMB11yk4D788ENbzZo1ba6urrb27dvbNm3aZHSkUu96n8958+YZHa3M6datm23MmDFGxygTfvzxR1vTpk1tbm5utoYNG9o+/fRToyOVCampqbYxY8bYatasaXN3d7fVqVPH9n//93+2rKysIj2PxkESERERuYI6moiIiIhcQQWSiIiIyBVUIImIiIhcQQWSiIiIyBVUIImIiIhcQQWSiIiIyBVUIImIiIhcQQWSiIiIyBVUIImIiIhcQQWSiIiIyBVUIImIiIhcQQWSiAiQnJxMQEAAb731lmPbhg0bcHV1JSoqysBkImIETVYrIvKHFStW0KdPHzZs2ECDBg1o2bIlDz/8MNOmTTM6mojcZiqQRET+x4gRI1i1ahVt27Zlz549bN26FTc3N6NjichtpgJJROR/XLp0iaZNmxIbG0t0dDTNmjUzOpKIGEB9kERE/sfRo0eJj4/HarVy4sQJo+OIiEHUgiQi8ofs7Gzat29Py5YtadCgAREREezZswc/Pz+jo4nIbaYCSUTkDy+99BJLlixh165dVKxYkW7duuHt7c3y5cuNjiYit5kusYmIAGvWrCEiIoIvvvgCLy8vzGYzX3zxBb/++iszZ840Op6I3GZqQRIRERG5glqQRERERK6gAklERETkCiqQRERERK6gAklERETkCiqQRERERK6gAklERETkCiqQRERERK6gAklERETkCiqQRERERK6gAklERETkCiqQRERERK7w/wvv7vNNx+R4AAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0ce44cb9-35ec-4dec-bb5d-aa9bc0d75f03",
      "metadata": {
        "id": "0ce44cb9-35ec-4dec-bb5d-aa9bc0d75f03"
      },
      "source": [
        "### Problem 2 - Bayes Theorem [5 points]\n",
        "\n",
        "**Solve the following question using python.**\n",
        "\n",
        "**Suppose the probability of the weather being cloudy is 40%. Also suppose the probability of rain on a given day is 20%. Also suppose the probability of clouds on a rainy day is 85%. If it’s cloudy outside on a given day, what is the probability that it will rain that day?**"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Given:\n",
        "\n",
        "To solve this problem using Bayes' theorem, we have the following information:\n",
        "\n",
        "- $P(Cloudy)$: Probability of it being cloudy, which is 40%.\n",
        "- $P(Rain)$: Probability of it raining on any given day, which is 20%.\n",
        "- $P(Cloudy∣Rain)$: Probability of it being cloudy given that it is raining, which is 85%.\n",
        "\n",
        "We're trying to find:\n",
        "\n",
        "$P(Rain∣Cloudy)$\n",
        "\n",
        "Using Bayes' theorem, we can express this as:\n",
        "\n",
        "$$ P(\\text{Rain}|\\text{Cloudy}) = \\frac{P(\\text{Cloudy}|\\text{Rain}) \\times P(\\text{Rain})}{P(\\text{Cloudy})} $$\n"
      ],
      "metadata": {
        "id": "ha3FiIRoEYCm"
      },
      "id": "ha3FiIRoEYCm"
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "1c3036f5-f214-48ee-a2ff-93009d010435",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1c3036f5-f214-48ee-a2ff-93009d010435",
        "outputId": "638d31d4-3b85-42ef-c090-b903ea6b1e42"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.425"
            ]
          },
          "metadata": {},
          "execution_count": 9
        }
      ],
      "source": [
        "# your code here\n",
        "\n",
        "# define function for Bayes' theorem\n",
        "def bayesTheorem(pA, pB, pBA):\n",
        "    return (pBA * pA) / pB\n",
        "\n",
        "# define probabilities\n",
        "pRain = 0.20\n",
        "pCloudy = 0.40\n",
        "pCloudyRain = 0.85\n",
        "\n",
        "# use function to calculate conditional probability\n",
        "bayesTheorem(pRain, pCloudy, pCloudyRain)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2dd0cfc0-be5c-4987-babc-93735aba1f32",
      "metadata": {
        "id": "2dd0cfc0-be5c-4987-babc-93735aba1f32"
      },
      "source": [
        "### Problem 3 - 10 points (Bonus)\n",
        "\n",
        "**Solve the following problem using Python by filling in the code blocks using comments.**\n",
        "\n",
        "Suppose a drug test is 95% accurate in detecting the presence of a certain drug in a person's system. That is, if a person has the drug in their system, there is a 95% chance that the test will come back positive. On the other hand, if a person does not have the drug in their system, there is a 5% chance that the test will falsely report a positive result.\n",
        "\n",
        "Suppose a random person is selected from a population where the prevalence of the drug is 1%. If that person takes the drug test and it comes back positive, what is the probability that they actually have the drug in their system?"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Defining the prior probabilities and test accuracies\n",
        "prior_drug = 0.01  # What's the chance someone from the population has the drug?\n",
        "prior_no_drug = 1 - prior_drug  # If we know the chance of having the drug, what's the chance of NOT having it?\n",
        "\n",
        "# Likelihoods\n",
        "p_pos_given_drug = 0.95  # If a person has the drug, what's the chance the test detects it?\n",
        "p_pos_given_no_drug = 0.05  # Even if a person doesn't have the drug, the test isn't perfect. What's the chance it falsely reports a positive?\n",
        "\n",
        "# Calculate the unnormalized probabilities (prior x likelihood)\n",
        "unnormalized_drug = prior_drug * p_pos_given_drug  # Multiply the prior and the likelihood for people who have the drug\n",
        "unnormalized_no_drug = prior_no_drug * p_pos_given_no_drug  # Now, do the same for those who don't\n",
        "\n",
        "# Calculate the evidence (probability of a positive test)\n",
        "p_pos = unnormalized_drug + unnormalized_no_drug  # Total probability of getting a positive test result\n",
        "\n",
        "# Calculate the posterior probabilities\n",
        "posterior_drug = unnormalized_drug / p_pos  # Given a positive test, what's the updated chance of truly having the drug?\n",
        "posterior_no_drug = unnormalized_no_drug / p_pos  # And if the test is positive, what's the updated chance they don't have the drug?\n",
        "\n",
        "# Printing the results\n",
        "print(\"Bayes Table:\")\n",
        "print(f\"{'Situation':<30} {'Prior':<10} {'Likelihood':<15} {'Unnormalized':<20} {'Posterior':<10}\")\n",
        "print('-' * 90)\n",
        "print(f\"{'Has Drug':<30} {prior_drug:<10.2f} {p_pos_given_drug:<15.2f} {unnormalized_drug:<20.5f} {posterior_drug:<10.5f}\")\n",
        "print(f\"{'Doesnt Have Drug':<30}\", end=' ')\n",
        "print(f\"{prior_no_drug:<10.2f}\", end=' ')\n",
        "print(f\"{p_pos_given_no_drug:<15.2f} {unnormalized_no_drug:<20.5f} {posterior_no_drug:<10.5f}\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "M-fU5t_zHkFl",
        "outputId": "4fd22907-889e-4110-8df2-78572c4de9e0"
      },
      "id": "M-fU5t_zHkFl",
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Bayes Table:\n",
            "Situation                      Prior      Likelihood      Unnormalized         Posterior \n",
            "------------------------------------------------------------------------------------------\n",
            "Has Drug                       0.01       0.95            0.00950              0.16102   \n",
            "Doesnt Have Drug               0.99       0.05            0.04950              0.83898   \n"
          ]
        }
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.13"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}