{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e83a02ca-6265-4d13-9a29-7e08e7ac7568",
   "metadata": {},
   "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": {},
   "source": [
    "### Problem 1 - MLE [10 points]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f280efc4-5207-46b4-9a16-d2a656933b2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the following libraries\n",
    "# numpy\n",
    "# matplotlib.pyplot\n",
    "# scipy.stats\n",
    "# scipy.optimize\n",
    "# pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2cdf87eb-7f11-4d8d-b5ae-3b22a41cf18b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate normally distributed data\n",
    "np.____.____(123) # Use a random seed of 123\n",
    "data = np.____.____(loc=3, scale=2, ____=___) # Use random.normla to generate normally distributed data of size 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "347eb869-2121-4234-b709-593a64f844ae",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.82873879  4.99469089  3.565957   -0.01258943  1.8427995   6.30287307\n",
      " -1.85335849  2.14217474  5.53187252  1.2665192 ]\n"
     ]
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5106ab7f-4174-4bd6-8ed1-cc425f4234c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood(x, mu, sigma):\n",
    "  #your code here\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": {},
   "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": 6,
   "id": "2ce3d1b4-0904-4100-8d08-2ae05b81197f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_likelihood(x, mu, sigma):\n",
    "  #your code here\n",
    "  p = _____\n",
    "  log_lhood = np.sum(p)\n",
    "  return log_lhood"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ceba8dcd-f48a-4bfc-a998-413405902e47",
   "metadata": {},
   "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": 7,
   "id": "d086a093-9ff9-4794-af4c-12afe122587c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#maximum likelihood estimation can be defined as:\n",
    "def max_likelihood_estimation(x):\n",
    "    # defining the negative log-likelihood function\n",
    "    neg_log_likelihood = lambda params: -np.____(np.____(likelihood(x, *params)))\n",
    "    \n",
    "    # Set the initial guess for the parameters and optimize\n",
    "    mu_init = np._____(x) # Set the mean\n",
    "    sigma_init = np.____(x) # Set the intial 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",
   "execution_count": 8,
   "id": "f00392df-500a-464b-85d8-b768a33ba9c1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "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"
     ]
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1fb798a9-8086-465c-b4c0-66b7a7edd739",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAGwCAYAAAC99fF4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABtmElEQVR4nO3deVxU9f7H8dcwbIIsAgKigGi5oyLmgllZpmmbZaWWS6WV3crU22bWzayrv8rKNrXFtUWttLIyl8o1dxKXJLVcUBYRVBCQbZjfH5NTXFAZBA/L+/l4nMedOXznnPfMJefD93zP92uyWq1WRERERKQYJ6MDiIiIiFRFKpJERERESqEiSURERKQUKpJERERESqEiSURERKQUKpJERERESqEiSURERKQUzkYHqK6KiopISkrCy8sLk8lkdBwREREpA6vVyunTpwkJCcHJ6fx9RSqSyikpKYnQ0FCjY4iIiEg5HDlyhEaNGp23jYqkcvLy8gJsH7K3t7fBaURERKQsMjMzCQ0NtX+Pn4+KpHI6e4nN29tbRZKIiEg1U5ahMhq4LSIiIlIKFUkiIiIipVCRJCIiIlIKjUkSERFDWSwWCgoKjI4hNYSLiwtms7lCjqUiSUREDGG1WklJSeHUqVNGR5EaxtfXl+Dg4Iuex1BFkoiIGOJsgRQYGIiHh4cm5pWLZrVaycnJITU1FYAGDRpc1PFUJImIyCVnsVjsBZK/v7/RcaQGqVOnDgCpqakEBgZe1KU3DdwWEZFL7uwYJA8PD4OTSE109vfqYse6qUgSERHD6BKbVIaK+r1SkSQiIiJSChVJIiIiIqXQwG0REalS3ly575Keb8z1zS7p+aT6UE+SiIiIA+69915MJhMmkwkXFxeCgoK4/vrrmTVrFkVFRWU+zpw5c/D19a28oHLR1JMkIlLRrFbIOw25GZB7Cs6csj2v4wtewVA3GFx1V1d1dsMNNzB79mwsFgvHjh1j2bJlPP7443z55ZcsWbIEZ2d9vdYE+n9RRORiFeZB8g5I2ARHNtv+Nyft/K9x94F6jaFxd2jSA8JjVDhVI25ubgQHBwPQsGFDOnToQJcuXbjuuuuYM2cOI0aM4I033mD27NkcOHAAPz8/br75Zl599VXq1q3L6tWrue+++4C/78R64YUXmDBhAp988glTp05l7969eHp6cu211zJ16lQCAwMNe7+1leGX26ZNm0ZERATu7u5ER0ezbt26c7ZNTk7m7rvvpnnz5jg5OTF69OgSba655hp7N+g/txtvvNHeZsKECSV+fvaXXUSkTKxWOPQLLH4Q/i8cZl4PK5+H37/7u0Ayu4JnIAQ0g4bR4NcEnG0T3ZGbYSusNr4Ln/aHV8Jhzk2w9SPIyzLufUm5XXvttbRr147FixcD4OTkxNtvv83u3buZO3cuP//8M0899RQAMTExTJ06FW9vb5KTk0lOTuaJJ54AID8/n5deeokdO3bw9ddfc/DgQe69916j3latZmhP0sKFCxk9ejTTpk2jW7duvP/++/Tp04c9e/YQFhZWon1eXh7169dn/PjxvPnmm6Uec/HixeTn59ufp6en065dO+68885i7Vq3bs2PP/5of15Ri+GJSA2XnQ47PoPYuZC+/+/9Hv4Q2gXCOtv+N7gNuHjA/87XYrVCXiacToGUXXBgtW3LOAKH1tm2HydChyHQ6UGoF34p351cpBYtWrBz506AYn/IR0RE8NJLL/Hwww8zbdo0XF1d8fHxKfWP9Pvvv9/+uEmTJrz99tt06tSJrKws6tate0neh9gYWiS98cYbDB8+nBEjRgAwdepUli9fzvTp05k8eXKJ9o0bN+att94CYNasWaUe08/Pr9jzBQsW4OHhUaJIcnZ2dqj3KC8vj7y8PPvzzMzMMr9WRGqA/GzY8C788hYUZNv2uXhCZH/ocC807FCyICqNyWS71ObuA/WbQ+QdtsLpxAHY+wNsm2l7vPFd2DQNmveF6/5jaytVntVqtV8+W7VqFZMmTWLPnj1kZmZSWFhIbm4u2dnZeHp6nvMY27dvZ8KECcTFxXHixAn7YPCEhARatWp1Sd6H2Bh2uS0/P5/Y2Fh69epVbH+vXr3YsGFDhZ1n5syZDBw4sMQv5P79+wkJCSEiIoKBAwdy4MCB8x5n8uTJ+Pj42LfQ0NAKyygiVVhREcR9Bu90hNWTbAVScCTcNBWe2Au3vAONostWIJ2LyQT+TSHmUXg0Fu7+3DZOyVpku3w3PQZWPG8b/C1VWnx8PBERERw+fJi+ffvSpk0bFi1aRGxsLO+99x5w/qUysrOz6dWrF3Xr1uWTTz5h69atfPXVVwDFrpLIpWFYkZSWlobFYiEoKKjY/qCgIFJSUirkHFu2bGH37t32nqqzOnfuzLx581i+fDkffvghKSkpxMTEkJ6efs5jjRs3joyMDPt25MiRCskoIlXY0Vj44Gr4+mE4nQQ+YdB/Jjy0DjreB25eFX9OJydo1huGfg3/2gzN+kBRIWx421ao7fzC1vMkVc7PP//Mrl276N+/P9u2baOwsJDXX3+dLl260KxZM5KSkoq1d3V1xWKxFNv3+++/k5aWxv/93//RvXt3WrRoYV/RXi49w+9u+9/1Vf7ZVXmxZs6cSZs2bejUqVOx/X369LE/joyMpGvXrjRt2pS5c+cyduzYUo/l5uaGm5tbheQSkSquyALr34TVk20Fips3dP83dB4JLu6XLkdgC7h7AexbDj88DScPwuIREPcJ3PYBeAVd+BhSKfLy8khJSSk2BcDkyZO56aabGDp0KLt27aKwsJB33nmHm2++mV9++YUZM2YUO0bjxo3Jysrip59+ol27dnh4eBAWFoarqyvvvPMOI0eOZPfu3bz00ksGvUsxrEgKCAjAbDaX6DVKTU0t0btUHjk5OSxYsICJEydesK2npyeRkZHs37//gm1FpIbLOAqLH4LD623PW98OfV8DzwDjMjXrDRFX28YprZ1iG+g940ro/xE0udq4XJWkOsyAvWzZMho0aICzszP16tWjXbt2vP322wwbNgwnJyfat2/PG2+8wSuvvMK4ceO46qqrmDx5MkOHDrUfIyYmhpEjRzJgwADS09PtUwDMmTOHZ599lrfffpsOHTowZcoUbrnlFgPfbe1lslqN67ft3Lkz0dHRTJs2zb6vVatW3HrrraUO3P6na665hvbt2zN16tRSfz5nzhxGjhxJYmIi/v7+5z1WXl4eTZs25cEHH+Q///lPmbJnZmbi4+NDRkYG3t7eZXqNiFRxe5bAksdsE0C6eMKNU6DdoIsbb1TRju+DL4ZB6h4wOcHVz8BVT4BT9bpDNzc3l4MHD9qngBGpSOf7/XLk+9vQeZLGjh3LRx99xKxZs4iPj2fMmDEkJCQwcuRIwDYO6J9VN0BcXBxxcXFkZWVx/Phx4uLi2LNnT4ljz5w5k379+pVaID3xxBOsWbOGgwcPsnnzZu644w4yMzMZNmxY5bxREanarFZY/Qp8PsRWIIVEwch10P7uqlUgAdRvBiN+gqghtoHdqyfBx7dB9gUmrxQRhxk6JulsF+PEiRNJTk6mTZs2LF26lPBw27wgycnJJCQkFHtNVFSU/XFsbCyfffYZ4eHhHDp0yL5/3759rF+/nhUrVpR63qNHjzJo0CDS0tKoX78+Xbp0YdOmTfbzikgtYimAbx+HuE9tz7s+Cte9AM6uxuY6H1cPuPVdCO8G34+Fg2tgVm8YvFjzKolUIEMvt1VnutwmUgPkZsDnQ21jfExOcOPr0PH+C76sSkmNh0/vtE1GWTcYBi+yTWRZxelym1SmGnG5TUTEMBlHYVYfW4Hk4gmDFla/AgkgsCUMXwGBrSArBWb3hUPrjU4lUiOoSBKR2ifjqK2YSP3N1vty31Jo1uvCr6uqvENs7yEsBvIy4OPbIf5bo1OJVHsqkkSkdslMsi0ke+ow1IuAET9CSHujU128OvVgyGJocRNY8uDzYRD/ndGpRKo1FUkiUnucToG5N9smZfQNh3u/A98atMSQSx24c65t2gKrBb68D/74yehUItWWiiQRqR2yUmHuLZD+B/iEwrBvwaeR0akqntkZbnkXWt0KlnxYcA8crrj1MEVqExVJIlLz5ZyAebdC2l7wbmgrkGryrfJmZ7j9I7i8FxSegU/vgsRfjU4l5zFnzhx8fX2NjuGQ6pjZUSqSRKRmK8yDhUNsM1TXDbYVSH4RRqeqfM6ucNc8aNwd8k/DJ7fDsZIT74rj7r33XkwmU4nthhtuKNPrGzduXGK1iAEDBrBv375KSFvcpS5s/vn5eHl50bFjRxYvXmz/+YQJE+w/d3Z2JiAggKuuuoqpU6eSl5dX7FjXXHNNqZ97YWFhpeVXkSQiNZfVCktG2dZhc/WyDWz2b2p0qkvHpQ4Mmg8NO8KZk/DZAMg6bnSqGuGGG24gOTm52DZ//vxyH69OnToEBgZWYMKqY/bs2SQnJ7N161batWvHnXfeycaNG+0/b926tX3y6FWrVnHnnXcyefJkYmJiOH36dLFjPfDAAyU+d2fnypsXW0WSiNRca16FnQvAZIa75kBQa6MTXXpuXnDPF+DXFDISYMHdUJBrdKpqz83NjeDg4GJbvXr17D+fMGECYWFhuLm5ERISwqhRowBbb8jhw4cZM2aMvScESvbwTJgwgfbt2zNr1izCwsKoW7cuDz/8MBaLhVdffZXg4GACAwP573//WyzXG2+8QWRkJJ6enoSGhvKvf/2LrKwsAFavXs19991HRkaG/dwTJkwAID8/n6eeeoqGDRvi6elJ586dWb16dbFjz5kzh7CwMDw8PLjttttIT08v02fl6+tLcHAwLVq0YMaMGbi7u7NkyRL7z52dnQkODiYkJITIyEgee+wx1qxZw+7du3nllVeKHcvDw6PE516ZVCSJSM20Y6FtXTOwzaR9WU9j8xjJww/uXgjuPnB0C3w7ytbLVtVYrZCfbcxWgZ/Hl19+yZtvvsn777/P/v37+frrr4mMjARg8eLFNGrUyL4cV3Jy8jmP8+eff/LDDz+wbNky5s+fz6xZs7jxxhs5evQoa9as4ZVXXuG5555j06ZN9tc4OTnx9ttvs3v3bubOncvPP//MU089BUBMTAxTp07F29vbfu4nnngCgPvuu49ffvmFBQsWsHPnTu68805uuOEG9u/fD8DmzZu5//77+de//kVcXBw9evTg5ZdfdvizcXFxwdnZmYKCgvO2a9GiBX369Cl2ac4Ihq7dJiJSKQ79AksetT3u9jh0vM/YPFVBwOW26QE+6Q87F0JAM7jqCaNTFVeQA5NCjDn3s0ng6lnm5t999x1169Yttu/pp5/m+eefJyEhgeDgYHr27ImLiwthYWF06tQJAD8/P8xmM15eXhfsBSkqKmLWrFl4eXnRqlUrevTowd69e1m6dClOTk40b96cV155hdWrV9OlSxcARo8ebX99REQEL730Eg8//DDTpk3D1dUVHx8fTCZTsXP/+eefzJ8/n6NHjxISYvv8n3jiCZYtW8bs2bOZNGkSb731Fr179+aZZ54BoFmzZmzYsIFly5aV+TPLy8vjtddeIzMzk+uuu+6C7Vu0aFFiDdZp06bx0Ucf2Z8/9NBDvP7662XO4CgVSSJSs2Qchc+H2G5/b3kLXDfB6ERVR9Me0Pc126K4P79kK5xa3Wp0qmqpR48eTJ8+vdg+Pz8/AO68806mTp1KkyZNuOGGG+jbty8333yzw2NnGjdujJeXl/15UFAQZrMZJyenYvtSU1Ptz1etWsWkSZPYs2cPmZmZFBYWkpubS3Z2Np6epReBv/76K1arlWbNmhXbn5eXh7+/PwDx8fHcdtttxX7etWvXMhVJgwYNwmw2c+bMGXx8fJgyZQp9+vS54OusVqv9cuRZ99xzD+PHj7c/r+xB6CqSRKTmKMyHL+6FnHQIjoTb3gcnjSoo5orhcHwvbHkfFj8E/pdDUCujU9m4eNh6dIw6twM8PT257LLLSv1ZaGgoe/fuZeXKlfz444/861//4rXXXmPNmjW4uLiUPdL/tDWZTKXuKyoqAuDw4cP07duXkSNH8tJLL+Hn58f69esZPnz4eS9vFRUVYTabiY2NxWw2F/vZ2d4y60VcjnzzzTfp2bMn3t7eDg1Oj4+PJyKi+J2oPj4+5/zcK4OKJBGpOVY+D0e32sbe3PUxuDr2xVdr9J4EafvgwCr4Yhg8sArc6l74dZXNZHLokldVVqdOHW655RZuueUWHnnkEVq0aMGuXbvo0KEDrq6uWCyWCj/ntm3bKCws5PXXX7f3Nn3++efF2pR27qioKCwWC6mpqXTv3r3UY7dq1arY2CegxPNzCQ4Odriw+f3331m2bBnjxo1z6HUVTUWSiNQMuxfB5hm2x7e9XzvmQiovszP0/whmXGkrlr4fa/vM/ufShpxbXl4eKSkpxfadnednzpw5WCwWOnfujIeHBx9//DF16tQhPNw2gWnjxo1Zu3YtAwcOxM3NjYCAgArJ1LRpUwoLC3nnnXe4+eab+eWXX5gxY0axNo0bNyYrK4uffvqJdu3a4eHhQbNmzbjnnnsYOnQor7/+OlFRUaSlpfHzzz8TGRlJ3759GTVqFDExMbz66qv069ePFStWODQe6XwKCwtJSUmhqKiI9PR0Vq9ezcsvv0z79u158sknK+Qc5aUiSaSaeXNl5U84V9HGXN/swo0uxvG9tvmQAK4cA80vPN6h1vMMgDtm2Rb73bkQwrtB9DCjU1Uby5Yto0GDBsX2NW/enN9//x1fX1/+7//+j7Fjx2KxWIiMjOTbb7+1j++ZOHEiDz30EE2bNiUvL++iLmX9U/v27XnjjTd45ZVXGDduHFdddRWTJ09m6NCh9jYxMTGMHDmSAQMGkJ6ezgsvvMCECROYPXs2L7/8Mv/+979JTEzE39+frl270rdvXwC6dOnCRx99ZG/fs2dPnnvuOV566aWLzv3bb7/RoEEDzGYzPj4+tGrVinHjxvHwww/j5uZ20ce/GCZrRf2/U8tkZmbi4+NDRkYG3t7eRseRWkRF0v/Iy4IPr7UtOdK4Owz52tZTImWz7g346UVwdocRP0Fwm0ty2tzcXA4ePEhERATu7u6X5JxSe5zv98uR72+NaBSR6m35OFuBVDfY1jOiAskx3Ub/tcZbrm18Ut7pC75EpLZQkSQi1Vf8d/DrPMBkG2NTt2Yu61CpnJyg3wzbwr/pf8D3VWzuJBEDqUgSkerp9DHbzNEAMY9BROl35UgZePrDHbPB5GRbxmXPN0YnEqkSVCSJSPVjtcI3j9jmQwqKhGufMzpR9RfWGa4ca3v87WhbESpSy6lIEpHqZ+tH8MdKMLtB/w/B2dg7YGqMq5+G4LZw5sQlW99N9w5JZaio3ysVSSJSvRzfByuetz2+fiIEtjQ2T03i7Aq3f2ArPvct+2u8V+U4O3N0Tk5OpZ1Daq+zv1eOzHBeGt0GIiLVh6UQvnoQCs9Akx7Q6UGjE9U8gS3huudhxXOw/FmIuKpSJuY0m834+vra1x3z8PAosU6XiKOsVis5OTmkpqbi6+tbYpkVR6lIEpHqY+O7kLQd3H2h3zSty1ZZujwCe5fB4fXw9cNw7/fgdHFfNqU5uxL9PxdoFakIvr6+9t+vi6EiSUSqh7Q/YPVk2+Pek8A7xNg8NZmTk60InR4DCRthy4fQZWSFn8ZkMtGgQQMCAwPPuwCriCNcXFwuugfpLBVJIlL1FRXZBhIX5tous7W/2+hENV+9cOj1Enw3Bn6aCC1uBN/QSjmV2WyusC81kYqkvmoRqfp+nQOHfwEXT7j5LS3Eeql0uBfCYqAg27YIru5Ek1pGRZKIVG0ZibDiP7bH1/3H1sMhl4aTk60oNbvC/hWwe5HRiUQuKRVJIlJ1Wa22yz35p6FRJ+j0gNGJap/6zaD7X0uVLHsGck4Ym0fkElKRJCJV1+5FsH+5rSfjlncq5Q4rKYMrx0D9FpB9/O85qkRqARVJIlI15WbC8vG2x92fgMAWxuapzZz/KlIxQdwncGC10YlELgkVSSJSNa2eDFkp4NcUrhxtdBoJ7QRXjLA9/m4MFOYZm0fkElCRJCJVT8pu2Py+7XHf17Q2W1Vx3X+gbjCcOAAb3jE6jUil0zxJIlLp3ly5r+yNrUXctetfNLRa2Od/Hd8fDIWDDry+goy5vtklP2eV5+5tmztp8QOwdgq0HXDRcyc59LtRReh3o/ZQT5KIVCmtUr+n4ekd5DvVYU3EGKPjyP+KvNM2d1LhGdvabiI1mIokEaky3Aoz6X7YdhlnU+gDZLkFGZxISjCZbJdATWaIXwJ//mx0IpFKoyJJRKqMboen4VFwkvQ6EWwPGWR0HDmX4DZ/z1m19CkozDc2j0glMbxImjZtGhEREbi7uxMdHc26devO2TY5OZm7776b5s2b4+TkxOjRo0u0mTNnDiaTqcSWm5tb7vOKSOULyN5H25TFAPzc9GmKnDRkskq7Zhx41of0/bBpmtFpRCqFoUXSwoULGT16NOPHj2f79u10796dPn36kJCQUGr7vLw86tevz/jx42nXrt05j+vt7U1ycnKxzd3dvdznFZFKZrVyzYE3MGFlb8D1HPWJNjqRXEgdX7h+ou3xmldty8eI1DCGFklvvPEGw4cPZ8SIEbRs2ZKpU6cSGhrK9OnTS23fuHFj3nrrLYYOHYqPj885j2symQgODi62Xcx5RaRyNT2xmtDMWAqd3FgX/pjRcaSs2g60LRdTkA0/v2R0GpEKZ1iRlJ+fT2xsLL169Sq2v1evXmzYsOGijp2VlUV4eDiNGjXipptuYvv27Rd93ry8PDIzM4ttInLxzEX5XHXoLQC2hdzDafcGBieSMnNygj7/Z3u8Yz4k/mpsHpEKZliRlJaWhsViISio+N0rQUFBpKSklPu4LVq0YM6cOSxZsoT58+fj7u5Ot27d2L9//0Wdd/Lkyfj4+Ni30NCLmxtERGzaJy3ENzeRLJcAtjUaZnQccVTDaNt8SWBbRsZqNTaPSAUyfOC2yWQq9txqtZbY54guXbowePBg2rVrR/fu3fn8889p1qwZ77xTfHZYR887btw4MjIy7NuRI0fKnVFEbDzy0+l8dCYAvzR+hAKzh8GJpFyu+w8414GEDRD/rdFpRCqMYUVSQEAAZrO5RO9NampqiV6ei+Hk5MQVV1xh70kq73nd3Nzw9vYutonIxema8D5ulmyOebZkT/2+RseR8vJpBDF/jSVb+bzWdZMaw7AiydXVlejoaFauXFls/8qVK4mJiamw81itVuLi4mjQoMElPa+InF9A9j7aHPsGgNVNxoLJ8I5tuRjdHret63byEGz5wOg0IhXC0IlIxo4dy5AhQ+jYsSNdu3blgw8+ICEhgZEjRwK2S1yJiYnMmzfP/pq4uDjANjj7+PHjxMXF4erqSqtWrQB48cUX6dKlC5dffjmZmZm8/fbbxMXF8d5775X5vCJS+a46+BZOFLHXvydJ3u2NjiMXy60uXPc8fPMIrHkN2t0Nnv5GpxK5KIYWSQMGDCA9PZ2JEyeSnJxMmzZtWLp0KeHh4YBt8sj/nbsoKirK/jg2NpbPPvuM8PBwDh06BMCpU6d48MEHSUlJwcfHh6ioKNauXUunTp3KfF4RqVxhJzcRnrEFi8mZ9Y0fNTqOVJR2g2DzDEjZBasnw41TjE4kclFMVqtuRSiPzMxMfHx8yMjI0PgkuaSq46rpxViLuGfHEAKz9/Frg0GsaTLW6ESl0krv5XRwHcy9yba226Nbwb/peZtXx99n/W5Ub458f2sQgIhcUi2OLyMwex95Zk82h95vdBypaBHd4fLeYLVogkmp9lQkicglYy7KIyZhBgBbGw0j18XX2EBSOXq+AJjgt68gMdboNCLlpiJJRC6Zdslf4pOXzGnXQLY3GGR0HKksQa2h3UDb45UvaIJJqbZUJInIJeFWeJpOR2cDsDHsQQrN7hd4hVRrPZ4FsyscWgd//mR0GpFyUZEkIpfEFUfnUKcwg/Q6EewJvNHoOFLZfMOg04O2xz9OgKIiQ+OIlIeKJBGpdJ55qUQlLwRgfeNHsZoMnX1ELpXu/wY3b9uUALsXGZ1GxGEqkkSk0nU+OgvnojwSvdpxoF53o+PIpeLhZ5uJG+DniVquRKodFUkiUql8co/S5tjXAPwS/i+4iAWspRrq8rBtuZJTCRA7x+g0Ig5RkSQilapLwoeYrRYO+XYh0aeD0XHkUnP1hKufsj1eOwXyc4zNI+IAFUkiUmn8c/6k5fEfANgQ9rDBacQwUUNsA7mzU2Hrh0anESkzFUkiUmm6JryPCSv7/a7hmFcro+OIUZxd4Zpxtsfr34TcTGPziJSRiiQRqRSBWfFcnr4KKyY2ho80Oo4YLfIu8L8czpyETdONTiNSJiqSRKRSxBy2fRH+Xv8G0j3Ov8ip1AJmZ+jxV2/Sxnch54SxeUTKQEWSiFS4hhnbiTi1EYvJzMbQB42OI1VFq9sgsDXkZcKGd4xOI3JBKpJEpGJZrcQk2HqRfgu8lYw6jQwOJFWGkxNcO972ePMMyEo1No/IBahIEpEK1SgjlkaZ2yk0ubA59H6j40hV07wvhHSAghzbIG6RKkxFkohUHKuVrkc+AGB3UD+y3IIMDiRVjskE1z5ne7x1Jp55x43NI3IeKpJEpMKEZmz9qxfJlS2N7jU6jlRVTa+F0C5gyeOKxLlGpxE5JxVJIlIxrFa6Jth6kXYF30a2W6DBgaTKMpngmmcAiEz5Sr1JUmWpSBKRChGWsYWGp3dQaHJla8NhRseRqq7JNRDaBWdrvnqTpMpSkSQiF69EL1J9gwNJlafeJKkGVCSJyEULO7WZkNM7KXRyY6vGIklZNbmGRK92OFvz6Zg4z+g0IiWoSBKRi/OPO9p2Bt1OtmuAwYGk2jCZ2BT2AABtj32FZ36awYFEinM2OoCIVG/hpzYRcnrXX71IQ42OU2HeXLnP6AgOG3N9M6MjOCzBpxOJXu1oeHoHHY/OZU2TfxsdScROPUkiUn5WK52PzARsvUg56kUSR6k3SaowFUkiUm6NMmLtd7RtazjE6DhSTZ3tTXIuyqPjUY1NkqpDRZKIlFvno7ZepN1Bt+iONik/k4lNoSMAiDy2mDr5JwwOJGKjIklEyiUkcwdhGduwmJzZ1kjzIsnFSfDtTHLd1rgU5RGd9KnRcUQAFUkiUk6dj3wEwJ7AmzjtFmxwGqn2TCa2/LUgcrvkL3EryDA4kIiKJBEph6DTv9H41CaKMLNVvUhSQQ7U606q5+W4FuUQlbzQ6DgiKpJExHFn72iLD7yBDPdGBqeRGsNkYkuj4QBEJS3AtTDL4EBS26lIEhGH1M/aS9OT6yjCiS2N7jM6jtQw+/17kF4nAnfLadqlfGF0HKnlVCSJiEM6H50FwL6A6zlVJ9zgNFLjmP4uvjskfoaz5YzBgaQ2U5EkImXml3OQy9JXAbAlVL1IUjn21r+eU+6N8Cg8RduUxUbHkVpMRZKIlNkVR+diwsp+v2tI92hqdBypoawmZ3tvUsfEjzEX5RmcSGorFUkiUibeuUm0OL4MgK2N7jU2jNR48fX7kOkWjGdBOq2PfWd0HKmlVCSJSJlEJ36MExYO+3bmmFdro+NIDVfk5GJf6qZj4jxM1kKDE0ltpCJJRC7IMz+NNseWALBZd7TJJbI78BayXfzwyUui+fEVRseRWkhFkohcUIekz3C25pPk1ZZE7w5Gx5FawmJ2Z3vIIAA6HZ0D1iJjA0mtY3iRNG3aNCIiInB3dyc6Opp169ads21ycjJ33303zZs3x8nJidGjR5do8+GHH9K9e3fq1atHvXr16NmzJ1u2bCnWZsKECZhMpmJbcLCWVRApjVtBBm2TFwHYBtOaTAYnktpkR/Ad5Jk98T9zkKYn1hodR2oZQ4ukhQsXMnr0aMaPH8/27dvp3r07ffr0ISEhodT2eXl51K9fn/Hjx9OuXbtS26xevZpBgwaxatUqNm7cSFhYGL169SIxMbFYu9atW5OcnGzfdu3aVeHvT6QmiEr+HNeiHFI9L+dgvW5Gx5FaJt+5LnEN7gLgiqNzwGo1NpDUKoYWSW+88QbDhw9nxIgRtGzZkqlTpxIaGsr06dNLbd+4cWPeeusthg4dio+PT6ltPv30U/71r3/Rvn17WrRowYcffkhRURE//fRTsXbOzs4EBwfbt/r161f4+xOp7lwKs4lKXgDAVvUiiUG2NxhIgZMbDbJ+IzRjq9FxpBYxrEjKz88nNjaWXr16Fdvfq1cvNmzYUGHnycnJoaCgAD8/v2L79+/fT0hICBEREQwcOJADBw6c9zh5eXlkZmYW20RqushjX+NemMlJ9zD2+19rdByppc64+rE7qB8AnY7ONjaM1CqGFUlpaWlYLBaCgoKK7Q8KCiIlJaXCzvPMM8/QsGFDevbsad/XuXNn5s2bx/Lly/nwww9JSUkhJiaG9PT0cx5n8uTJ+Pj42LfQ0NAKyyhSFZmL8olO+hSAbQ2HYDWZDU4ktdm2hoOxmMyEZWwj+LSGR8il4XCRdO+997J2bcUNnjP9T/e91Wotsa+8Xn31VebPn8/ixYtxd3e37+/Tpw/9+/cnMjKSnj178v333wMwd+7ccx5r3LhxZGRk2LcjR45USEaRqqrF8WXUzT9Olmt94gP7Gh1Harkst2Di69t+DzsdnWNsGKk1HC6STp8+Ta9evbj88suZNGlSiQHRZRUQEIDZbC7Ra5Samlqid6k8pkyZwqRJk1ixYgVt27Y9b1tPT08iIyPZv3//Odu4ubnh7e1dbBOpqUxWCx2P2v5oiA25B4uTq8GJRGBbw6FYMdH0xFr8cs4/REKkIjhcJC1atIjExEQeffRRvvjiCxo3bkyfPn348ssvKSgoKPNxXF1diY6OZuXKlcX2r1y5kpiYGEdjFfPaa6/x0ksvsWzZMjp27HjB9nl5ecTHx9OgQYOLOq9ITdE0fQ1+uQnkOnuz66+xICJGO+nRmD/8rwFss3CLVLZyjUny9/fn8ccfZ/v27WzZsoXLLruMIUOGEBISwpgxY87bI/NPY8eO5aOPPmLWrFnEx8czZswYEhISGDlyJGC7xDV06NBir4mLiyMuLo6srCyOHz9OXFwce/bssf/81Vdf5bnnnmPWrFk0btyYlJQUUlJSyMrKsrd54oknWLNmDQcPHmTz5s3ccccdZGZmMmzYsPJ8HCI1i9XKFYlzAIgLvpMCZ09j84j8w9aGtn+nWxxfhldexY1fFSnNRQ3cTk5OZsWKFaxYsQKz2Uzfvn357bffaNWqFW+++eYFXz9gwACmTp3KxIkTad++PWvXrmXp0qWEh4fbj/+/cyZFRUURFRVFbGwsn332GVFRUfTt+/d4iWnTppGfn88dd9xBgwYN7NuUKVPsbY4ePcqgQYNo3rw5t99+O66urmzatMl+XpHaLDRjK8FZ8RQ4uREXMsDoOCLFHPNqTYJPR8xWCx0SPzU6jtRwJqvVsZm5CgoKWLJkCbNnz7aP9xkxYgT33HMPXl5eACxYsICHH36YkydPVkroqiAzMxMfHx8yMjI0PkkuqTdX7qvU49+++xHCM7awvcEAVjd5olLPJRVrzPXNjI7gsPL8Poed2kz/3x6lwMmdjzp+S66Lb8UHO4/q+DnL3xz5/nZ29OANGjSgqKiIQYMGsWXLFtq3b1+iTe/evfH19XX00CJisKDTewjP2ILFZCY25B6j44iUKsGnE8c8mxOUvZf2yV+wKewBoyNJDeXw5bY333yTpKQk3nvvvVILJIB69epx8ODBi80mIpfY2cGwewNu4LS7bmSQKspkYlsj29ik9skLcbacMTiQ1FQOF0mrVq0q9S627Oxs7r///goJJSKXnu+Zw1ye/jMA2xoNMTiNyPnt97+WU+6NqFOYQeSxr42OIzWUw0XS3LlzOXOmZNV+5swZ5s3TLZki1VXHxE8wYeXPet1J92hqdByR87KazGxraCvmoxM/wamo7FPQiJRVmYukzMxMMjIysFqtnD59utgaZidPnmTp0qUEBgZWZlYRqSQe+Wm0TLXNPL+t0dALtBapGvYE3ki2iz9e+ak0T1thdBypgco8cNvX1xeTyYTJZKJZs5Ij+00mEy+++GKFhhORSyMqaQHO1gKSvNqS5N3e6DgiZWJxcmN7yECuPPweHRPnEV+/D5gMW5JUaqAyF0mrVq3CarVy7bXXsmjRIvz8/Ow/c3V1JTw8nJCQkEoJKSKVx7Uwi7YpiwDY2lC9SFK97AzuzxVH5xCQc4CIkxs46Hel0ZGkBilzkXT11VcDcPDgQcLCwipsEVoRMVZkyle4W7JIrxPBAb/uRscRcUiesxe7gm6jY9IndEycpyJJKlSZiqSdO3fSpk0bnJycyMjIYNeuXedse6HFZEWk6jAX5dMh6TMAtjUcrEsVUi39GjKIqOQFNMrcToPMnSR763tIKkaZiqT27duTkpJCYGAg7du3x2QyUdpE3SaTCYvFUuEhRaRytDi+jLoFaZx2DeT3+n2MjiNSLtlugcTX70Ob1G/pmPgx33q/ZnQkqSHKVCQdPHiQ+vXr2x+LSA1gLbJPHvlryCCKnFwMDiRSfrENh9Am9VuanlhDvZxDnPRobHQkqQHKVCT9c+FXLQIrUjM0PbEWvzOHyTXbxnSIVGcnPCL40+8qmp5YS3TSJ/x42XNGR5IaoFyTSX7//ff250899RS+vr7ExMRw+PDhCg0nIpXEarX3Iu1s0J8CZ0+DA4lcvLN3Z7ZMXYpnfprBaaQmcLhImjRpEnXq1AFg48aNvPvuu7z66qsEBAQwZsyYCg8oIhUv5PQOQk7votDkyvYGA4yOI1Ihkr3bkejVDmdrAVFJ842OIzWAw0XSkSNHuOyyywD4+uuvueOOO3jwwQeZPHky69atq/CAIlLxOiZ+DEB8YF9yXAMMTiNScc7OGN82ZRGuhVkGp5HqzuEiqW7duqSnpwOwYsUKevbsCYC7u3upa7qJSNXil3OQpifWYsVEbMg9RscRqVAH6l1Jep0I3CzZRB77yug4Us05XCRdf/31jBgxghEjRrBv3z5uvPFGAH777TcaN25c0flEpIJFJ34CwB9+V+sOIKl5TE7ENhwM2Jbb0cK3cjEcLpLee+89unbtyvHjx1m0aBH+/v4AxMbGMmjQoAoPKCIVxzPvOC2PLwW0kK3UXL/Xv4EslwC88lNpkbbc6DhSjZV5WZKzfH19effdd0vs1+K2IlVfVPICzNZCjnpHkeIVaXQckUphcXJle8gguh9+h+jEj9lTv69mk5dycbhIAjh16hRbtmwhNTWVoqIi+36TycSQIUMqLJyIVJx/LmR79nKESE21M/h2Oh2dpYVv5aI4XCR9++233HPPPWRnZ+Pl5VVsoVsVSSJVV+Sxr3CzZNsWsq2nLwyp2fKd69oXvo1O/FhFkpSLw/2P//73v7n//vs5ffo0p06d4uTJk/btxIkTlZFRRC6SU1EBHf6aN0YL2UptsT1kIBaTM6GZvxJ8erfRcaQacvhfysTEREaNGoWHh0dl5BGRStDi+DLq5h8nyyWAvfVvMDqOyCWR5RbE73/9vp+dG0zEEQ4XSb1792bbtm2VkUVEKoO1yP4FsT1kEBYnV4MDiVw6Z8ffXZa+Ct8zCQankerG4TFJN954I08++SR79uwhMjISF5fiK4ffcsstFRZORC5exMkN+J85SJ7Zk53BtxsdR+SSSvdoyoF6V9Lk5Ho6JH3Gz02fMTqSVCMOF0kPPPAAABMnTizxM5PJhMViufhUIlJhzi5kuyvoNvKd6xqcRuTS29ZwME1Orqf1sW/ZGPogZ1z9jI4k1YTDl9uKiorOualAEqlagk/vplHmdiwmZ7aHDDQ6joghEr07kFy3Nc7WfNonf250HKlGLuoWl9zc3IrKISKV4OxYpN/r9ybLLcjgNCIGMZmIbWibnqZ9yhe4WHIMDiTVhcNFksVi4aWXXqJhw4bUrVuXAwcOAPD8888zc+bMCg8oIuXjc+YIl6WvAiA2RJNHSu32h/81nHJvhHthJq2PLTE6jlQTDhdJ//3vf5kzZw6vvvoqrq5/3yUTGRnJRx99VKHhRKT8opM+xYSVA/W6ke55mdFxRAxlNZmJDbkHgA5Jn2GyFhqcSKoDh4ukefPm8cEHH3DPPfdgNpvt+9u2bcvvv/9eoeFEpHzq5J+g9bFvAdjWULPgiwD8FngTOS718MlLplnaT0bHkWrA4bvbEhMTueyykn+VFhUVUVBQUCGhROTitE/+HGdrPil1W5Ho3cHoOHKJvLlyn9ERqjSL2Z244DuJOfIB0YkfszegF/xjaS2R/+VwT1Lr1q1Zt25dif1ffPEFUVFRFRJKRMrPxZJD+5QvgL96kfQlIGK3o8GdFDi5E5S9l7CMLUbHkSrO4Z6kF154gSFDhpCYmEhRURGLFy9m7969zJs3j++++64yMoqIA1ofW4J7YSan3Bvxh38Po+OIVCm5Lr7sDrqVqOSFdEz8mATfzkZHkirM4Z6km2++mYULF7J06VJMJhP/+c9/iI+P59tvv+X666+vjIwiUkYmayEdkj4DIDbkHqwm8wVeIVL7/BoyiCKcCD+1mfpZe42OI1WYwz1JYFu/rXfv3hWdRUQuUrO0n/DJSybHpR6/Bd5kdByRKinTvSH7AnrSIm0F0Ykfs6z5y0ZHkirqoiaTFJEqxGq1L0ES1+AuLGZ3gwOJVF1nJ5dsnvYj3rlJBqeRqqpMPUn16tXDVMbBnydOnLioQCJSPmEZWwjM3keBkzs7gu8wOo5IlZZatwWHfToRnrGFDkmfsbrJE0ZHkiqoTEXS1KlT7Y/T09N5+eWX6d27N127dgVg48aNLF++nOeff75SQorIhZ1dgmR30K3kuvgaG0akGtjWaCjhGVtoc+wbNoWO0H83UkKZLrcNGzbMvv3yyy9MnDiR+fPnM2rUKEaNGsX8+fOZOHEia9ascTjAtGnTiIiIwN3dnejo6FKnFzgrOTmZu+++m+bNm+Pk5MTo0aNLbbdo0SJatWqFm5sbrVq14quvvrqo84pUdfWz9hJ+ajNFmPk15G6j44hUCwk+nUj1bIZLUS5tUxYZHUeqIIfHJC1fvpwbbrihxP7evXvz448/OnSshQsXMnr0aMaPH8/27dvp3r07ffr0ISEhodT2eXl51K9fn/Hjx9OuXbtS22zcuJEBAwYwZMgQduzYwZAhQ7jrrrvYvHlzuc8rUtWd7UXaF9CTTPcQg9OIVBMmE9saDgUgKnkhZosWbZfiHC6S/P39S+2Z+frrr/H393foWG+88QbDhw9nxIgRtGzZkqlTpxIaGsr06dNLbd+4cWPeeusthg4dio+PT6ltpk6dyvXXX8+4ceNo0aIF48aN47rrrit2ydDR84pUZd65STRLs/2BoiVIRByzL+A6Mtwa4FFwktap3xsdR6oYh6cAePHFFxk+fDirV6+2j0natGkTy5Ytc2iB2/z8fGJjY3nmmWeK7e/VqxcbNmxwNJbdxo0bGTNmTLF9vXv3thdJ5T1vXl4eeXl59ueZmZnlzihSkTokfYYTFhJ8ruB43eZGxxGpVqwmZ34NuZseB18nOukTdgX30/xiYudwT9K9997Lhg0b8PX1ZfHixSxatAgfHx9++eUX7r333jIfJy0tDYvFQlBQULH9QUFBpKSkOBrLLiUl5bzHLO95J0+ejI+Pj30LDQ0td0aRiuJecIo2x74GYGujYcaGEammdgfdyhlnH3xzj3JZ+iqj40gVUq7JJDt37synn35aIQH+d2oBq9Va5ukGLuaYjp533LhxjB071v48MzNThZIYrl3yF7gU5XHMszkJPp2MjiNSLRWa67Aj+A66HJ1Jx8SP2e9/ndY8FMDAySQDAgIwm80lem9SU1NL9PI4Ijg4+LzHLO953dzc8Pb2LraJGMnZkkv75M8BbINP9Y+6SLnFhQyg0MmN4Kw9NMqINTqOVBGGFUmurq5ER0ezcuXKYvtXrlxJTExMuY/btWvXEsdcsWKF/ZiVdV6RS6116hI8Ck+R4RbC/oBrjY4jUq2dcanH7sCbAewz14uU63JbRRk7dixDhgyhY8eOdO3alQ8++ICEhARGjhwJ2C5xJSYmMm/e37+wcXFxAGRlZXH8+HHi4uJwdXWlVatWADz++ONcddVVvPLKK9x666188803/Pjjj6xfv77M5xWp6kzWQqITbZe8Yxveg9Vk6H/KIjXCrw3voW3KYiJObSQgex9pns2MjiQGM/Rf1gEDBpCens7EiRNJTk6mTZs2LF26lPDwcMA2eeT/zl0UFRVlfxwbG8tnn31GeHg4hw4dAiAmJoYFCxbw3HPP8fzzz9O0aVMWLlxI586dy3xekaru8rSf8clLIsfZl98CbzE6jkiNkOHeiP0B19E8bSUdEz9mWbOXjI4kBjNZrVarIy+YM2cOd911Fx4eHpWVqVrIzMzEx8eHjIwMjU+SS+rNFXu5e8cQgrL3siH0QTaHPWB0JJEao37WXgbvGEwRZmZHLy51ctYx16uHqTpz5Pvb4TFJ48aNIzg4mOHDh1/UfEYiUj5hGVsIyt5rW8i2wZ1GxxGpUY7Xbc5h3844YaFD0mdGxxGDOVwkHT16lE8++YSTJ0/So0cPWrRowSuvvHJRcxuJSNl1PGobo6eFbEUqx9mZ69sc+xr3glPGhhFDOVwkmc1mbrnlFhYvXsyRI0d48MEH+fTTTwkLC+OWW27hm2++oaioqDKyikjSdsIztlCEmdiQe4xOI1IjJfh04phnC1yK8mif/IXRccRAFzUFQGBgIN26daNr1644OTmxa9cu7r33Xpo2bcrq1asrKKKI2K2fCsDv9Xtz2r2BsVlEaiqTyd6b1D55Ic6WMwYHEqOUq0g6duwYU6ZMoXXr1lxzzTVkZmby3XffcfDgQZKSkrj99tsZNkxLJIhUqPQ/Yc83gBayFals+wOu5ZR7I+oUZtDm2BKj44hBHC6Sbr75ZkJDQ5kzZw4PPPAAiYmJzJ8/n549ewJQp04d/v3vf3PkyJEKDytSq/3yFmDlQL0rSfe8zOg0IjWa1eRMbMhgAKKTPsGpqNDgRGIEh+dJCgwMZM2aNXTt2vWcbRo0aMDBgwcvKpiI/MPpFNgxH9BCtiKXym9BN9HlyAd456XQPG0F8YF9jY4kl5jDPUlXX301HTp0KLE/Pz/fPjO2yWTSxIwiFWnTNLDkQ2gXkrzbG51GpFawOLmxPWQQAB0T54JVNyXVNg4XSffddx8ZGRkl9p8+fZr77ruvQkKJyD/kZsC22bbHV442NIpIbbMj+A7yzJ4E5Bygycn1F36B1CgOF0lWqxVTKauNHz16FB8fnwoJJSL/sHUm5GVC/ZZweW+j04jUKvnOddkZ3B+AK47OAccWqZBqrsxjkqKiojCZTJhMJq677jqcnf9+qcVi4eDBg9xwww2VElKk1irIhU3TbY+7PQ5OFzVrh4iUw/aQQUQlLSDk9C4aZsYBzY2OJJdImYukfv36ARAXF0fv3r2pW7eu/Weurq40btyY/v37V3hAkVot7lPITgXvRhB5h9FpRGqlbNcA9gTeRNtji+mYOA8YYHQkuUTKXCS98MILADRu3JgBAwbg7u5eaaFEBLAU/nXbPxDzGJhdjM0jUottaziYNse+to1LStkNwW2MjiSXgMN998OGDVOBJHIp/LYYTh0GD3/oMNToNCK1WkadUPYHXGd7sv5NY8PIJVOmIsnPz4+0tDQA6tWrh5+f3zk3EakARUV//0Pc5WFw9TA2j4iwteFfc5T9thhOHDA2jFwSZbrc9uabb+Ll5WV/XNrdbSJSgfYvh9Q94OoFVzxgdBoRAY7Xbc7BejFEnNwAv7wNN081OpJUsjIVSf9ch+3ee++trCwiArZbjNe9YXt8xf1Qx9fQOCLyty2N7rMVSXGfwtVPg7cWmq7JylQkZWZmlvmA3t7e5Q4jIsDhX+DoFjC7QZdHjE4jIv+Q5N0ewrpCwkbY9B70etnoSFKJylQk+fr6XvAS29lJJi0WS4UEE6m1zvYiRQ0GryBjs4hISd3/DZ/eAVtnwZVjwUPjcWuqMhVJq1atquwcIoZ4c+U+oyMUE5j1O/f8+RNFmJltvZnMKpZPRIDLekJwJKTsgi0fwDXPGJ1IKkmZiqSrr766snOICH8tewDsrX89me4NjQ0jIqUzmWw9SF/eB5tnQNdHwa3uhV8n1U6ZiqSdO3fSpk0bnJyc2Llz53nbtm3btkKCidQ2fjkHuTz9Z+AftxqLSNXU6lbwawon/oTYORDzqNGJpBKUqUhq3749KSkpBAYG0r59e0wmE9ZSFvnTmCSR8rvi6FxMWPnD7xrSPS8zOo6InI+TGa4cDUseg43vQqcHwNnN6FRSwcpUJB08eJD69evbH4tIxfLOTaTF8WWA7RZjEakG2g6EVZPhdJJtSoCO9xudSCpYmYqk8PDwUh+LSMXomPgxTlg45NuFY16tjI4jImXh7ArdHodlT9tmyI8aojUWaxiH124D2Lt3L48++ijXXXcdPXv25NFHH2Xv3r0VnU2kVvDMS6X1sSUAbGmkv0RFqpUOQ8GzPpxKgF1fGp1GKpjDRdKXX35JmzZtiI2NpV27drRt25Zff/2VNm3a8MUXX1RGRpEaLTrpU5ytBRz1jiLRJ8roOCLiCFcP291tAOtehyKNy61JynS57Z+eeuopxo0bx8SJE4vtf+GFF3j66ae58847KyycSE1Xp+AkbVMWAxqLJFJtXTHcdrktfT/s+Qba3G50IqkgDvckpaSkMHTo0BL7Bw8eTEpKSoWEEqktopLm41KUS0rdlhz27WJ0HBEpDzcv6PKw7fHaKVBUZGweqTAOF0nXXHMN69atK7F//fr1dO/evUJCidQGboWnaZ/8OfDXWKQLLP0jIlVYpwfB1QtSf4N9y4xOIxWkTJfblixZYn98yy238PTTTxMbG0uXLra/fDdt2sQXX3zBiy++WDkpRWqgdsmf42bJJs2jCX/6XWV0HBG5GB5+0GmE7bLb2tegeR/94VMDmKylzQr5P5ycytbhVJsmk8zMzMTHx4eMjAy8vb2NjiPlZNTabS6F2YyIvQX3wky+b/Zf9tXvZUgOEXHcmOublf6DrOMwNRIKz8DgxXDZdZc2mJSJI9/fZap+ioqKyrTVlgJJ5GK1T/kC98JMTtQJZ3+A/iEVqRHq1oeOf92AsfY1uHAfhFRx5ZonSUTKz8WSQ4fETwHY3Oh+rCazwYlEpMLEjAKzGyRshEMlx+9K9eLwFAAA2dnZrFmzhoSEBPLz84v9bNSoURUSTKSmapuyCI/CU5x0D2WvLrOJ1CzeDWwTTG79EFa/AhEab1idOVwkbd++nb59+5KTk0N2djZ+fn6kpaXh4eFBYGCgiiSR83C25BKd+AlgmxfJairX3ykiUpVdOQZ+nQuH18Oh9dD4SqMTSTk5fLltzJgx3HzzzZw4cYI6deqwadMmDh8+THR0NFOmTKmMjCI1RuSxxXgWnCDDLYTf6/cxOo6IVAafhhA12PZ4zSvGZpGL4nCRFBcXx7///W/MZjNms5m8vDxCQ0N59dVXefbZZysjo0iNYLbkcsXReYCtF6nISb1IIjXWlWPByQUOroXDG41OI+Xk8L/SLi4umP6a+yEoKIiEhARatmyJj48PCQkJFR5QpKZok7oEz4J0Mt2C2RN4o9FxRKScyjp1SM/6NxJ57GsOf/UCi1u/W8mpLuycUxfIOTnckxQVFcW2bdsA6NGjB//5z3/49NNPGT16NJGRkQ4HmDZtGhEREbi7uxMdHV3qbN7/tGbNGqKjo3F3d6dJkybMmDGj2M+vueYaTCZTie3GG//+UpowYUKJnwcHBzucXaSszEV5XHF0LgBbGw6jyMnF4EQiUtm2NLoPi8lM+KnNNMjcaXQcKQeHi6RJkybRoEEDAF566SX8/f15+OGHSU1N5YMPPnDoWAsXLmT06NGMHz+e7du30717d/r06XPOHqmDBw/St29funfvzvbt23n22WcZNWoUixYtsrdZvHgxycnJ9m337t2YzeYSC++2bt26WLtdu3Y5+EmIlF2bY9/glZ/KaddAfgu6xeg4InIJZLqHEF/f9gd65yMfGZxGysPhy20dO3a0P65fvz5Lly4t98nfeOMNhg8fzogRIwCYOnUqy5cvZ/r06UyePLlE+xkzZhAWFsbUqVMBaNmyJdu2bWPKlCn0798fAD8/v2KvWbBgAR4eHiWKJGdnZ4d6j/Ly8sjLy7M/z8zMLPNrpXYzF+XR6egc4K+/LJ1cjQ0kIpfMlkb30Sr1eyJObST49G5SvNoYHUkcUO7JJFNTU1m3bh3r16/n+PHjDr8+Pz+f2NhYevUqPk9Mr1692LBhQ6mv2bhxY4n2vXv3Ztu2bRQUFJT6mpkzZzJw4EA8PT2L7d+/fz8hISFEREQwcOBADhw4cN68kydPxsfHx76FhoZe6C2KABCZ8hV184+T6RqkXiSRWiajTiPiA213snZNeN/gNOIoh4ukzMxMhgwZQsOGDbn66qu56qqrCAkJYfDgwWRkZJT5OGlpaVgsFoKCgortDwoKIiUlpdTXpKSklNq+sLCQtLS0Eu23bNnC7t277T1VZ3Xu3Jl58+axfPlyPvzwQ1JSUoiJiSE9Pf2ceceNG0dGRoZ9O3LkSFnfqtRiZkvu371IoferF0mkFtocOpwizDQ+tYkGmTuMjiMOcLhIGjFiBJs3b+a7777j1KlTZGRk8N1337Ft2zYeeOABhwOY/meVZKvVWmLfhdqXth9svUht2rShU6dOxfb36dOH/v37ExkZSc+ePfn+++8BmDt37jnP6+bmhre3d7FN5ELaHluMZ0E6GW4N+C3wZqPjiIgBMtwb8VvQTQB0TXBs7K4Yy+Ei6fvvv2fWrFn07t0bb29vvLy86N27Nx9++KG92CiLgIAAzGZziV6j1NTUEr1FZwUHB5fa3tnZGX9//2L7c3JyWLBgQYlepNJ4enoSGRnJ/v37y5xf5EKcLbn2O9ps8yLpjjaR2mpLo/uxmJwJz9hCw4xfjY4jZeRwkeTv74+Pj0+J/T4+PtSrV6/Mx3F1dSU6OpqVK1cW279y5UpiYmJKfU3Xrl1LtF+xYgUdO3bExaX4F9Dnn39OXl4egwcPvmCWvLw84uPj7XftiVSEtimL7LNr71Evkkitlukewu6gWwH1JlUnDhdJzz33HGPHjiU5Odm+LyUlhSeffJLnn3/eoWONHTuWjz76iFmzZhEfH8+YMWNISEhg5MiRgG0c0NChQ+3tR44cyeHDhxk7dizx8fHMmjWLmTNn8sQTT5Q49syZM+nXr1+JHiaAJ554gjVr1nDw4EE2b97MHXfcQWZmJsOGDXMov8i5OFvO0DHRNrv25tD7Nbu2iLCl0b0UmlwIzYyl0altRseRMijTv9xRUVHFxvzs37+f8PBwwsLCAEhISMDNzY3jx4/z0EMPlfnkAwYMID09nYkTJ5KcnEybNm1YunQp4eHhACQnJxebMykiIoKlS5cyZswY3nvvPUJCQnj77bftt/+ftW/fPtavX8+KFStKPe/Ro0cZNGgQaWlp1K9fny5durBp0yb7eUUuVruUL/EsOMEp94b2eVJEpHbLcgtmd1A/2qd8QdcjH/CFTzScZwyuGM9kPTvy+TxefPHFMh/whRdeuKhA1UVmZiY+Pj5kZGRoEHc1VtblBRzhWpjF/bH9qFOYwfLLnmePbvsXkb945qVyf+xtOFvzWdT6XRJ8O1+yc2tZEhtHvr/L1JNUWwofkYoQlTSfOoUZnKgTTnxgX6PjiEgVku0WyM7g2+mQvICYwzNI8Omk3qQqrNwDJWJjY4mPj8dkMtGqVSuioqIqMpdIteRWkEF00qcAbAx9EKtJY5FEpLitjYYReewrGmTtpsnJdRzwu8roSHIODv8LnpqaysCBA1m9ejW+vr5YrVYyMjLo0aMHCxYsoH79+pWRU6Ra6Jj4MW6WbI57XM6+gJ5GxxGRKijHNYC4BgO4InEeMYdncKDelWAq9wIYUokc/n/lscceIzMzk99++40TJ05w8uRJdu/eTWZmJqNGjaqMjCLVgkd+GlHJCwHYEPaQ/tETkXPa2nAoeWZP6ufsp1naj0bHkXNw+F/xZcuWMX36dFq2bGnf16pVK9577z1++OGHCg0nUp10OjoHl6Jckuu2Vve5iJxXnosPsQ1t8/jFJLyPyVpocCIpjcNFUlFRUYmJGwFcXFwoKiqqkFAi1U3dvBQiUxYDsCH8YQ3EFJEL+rXBIHKcfamXm0Cr1LKvWCGXjsNF0rXXXsvjjz9OUlKSfV9iYiJjxozhuuuuq9BwItVFlyMzcbYWcMS7g+1uFRGRCyhw9mRro3sB6JLwIeaifGMDSQkOF0nvvvsup0+fpnHjxjRt2pTLLruMiIgITp8+zTvvvFMZGUWqNN8zCbQ+9i0AG8L/pV4kESmzHcH9Oe0aiHf+MXtvtFQdDt/dFhoayq+//srKlSv5/fffsVqttGrVip49dSeP1E4xh6fjhIUD9bqR5N3O6DgiUo1YzO5sDh1Ozz8n0+nobHYH3UqhuY7RseQvDhVJhYWFuLu7ExcXx/XXX8/1119fWblEqoXArHiap/+IFRO/hD9idBwRqYZ+C7yFjokf45t7lA5J89kSer/RkeQvDl1uc3Z2Jjw8HIvFUll5RKqVKw+9C8Dv9W8gzfNyg9OISHVU5OTMhjDbwu4dE+fhXnDK2EBi5/CYpOeee45x48Zx4sSJysgjUm2EndpMeMYWLCZn27xIIiLltDfgeo55NsfNkk2no7ONjiN/cXhM0ttvv80ff/xBSEgI4eHheHp6Fvv5r7/+WmHhRKosq5UrD70HwM7g/mS6NzQ4kIhUayYn1oc/Sv89j9Eu+QviGgwg0z3E6FS1nsNF0q233opJd+9ILXd5+k8EZceT7+TB5kYaPyAiFy/BtzMJPlcQlrGVrgnvs7zZi0ZHqvUcLpImTJhQCTFEqg+nokK6HZ4OQGzDezjj6mdwIhGpEUwm1oc/yt07h9Hy+A/ENryHNM9mRqeq1co8JiknJ4dHHnmEhg0bEhgYyN13301aWlplZhOpklqnLqFebgI5LvWIDbnH6DgiUoMc82rF3oDrMWHlysPvGR2n1itzkfTCCy8wZ84cbrzxRgYOHMjKlSt5+OGHKzObSJXjUphN14T3Adjc6H4KnD0v8AoREcf8EvYwFpOZiJMbaHRqm9FxarUyX25bvHgxM2fOZODAgQAMHjyYbt26YbFYMJvNlRZQpCrpmPQJngUnOOXeiJ3B/Y2OIyI1UEadUHYF3U77lC/ofvgd5vvMBpPDN6NLBSjzp37kyBG6d+9uf96pUyecnZ2LreEmUpN55h0nOvETANaHP0qRU8mFnkVEKsLm0OHkO3kQnLWH5mkrjY5Ta5W5SLJYLLi6uhbb5+zsTGFhYYWHEqmKuia8j0tRLklekez3v9boOCJSg+W4+rO10VAAuh1+D3NRnsGJaqcyX26zWq3ce++9uLm52ffl5uYycuTIYnMlLV6sBfqk5vHP/oM2qUsAWNt4tBaxFZFK92vIPbRNWYxPXjJRSQvZ9lfRJJdOmYukYcOGldg3ePDgCg0jUlV1P/Q2Jqzs87+OZO+2RscRkVqg0OzOL+H/4ob9E+h0dBa/Bd3MGZd6RseqVcpcJM2erWnSpXYKO7mJiFMbsZicWa9FbEXkEoqv34eopPkEZe+l85GPWN3kSaMj1SoaLi9yHiarhasOvQ3AjuA7yKgTanAiEalVTE62S/xA25RF1Ms5ZGic2kZFksh5tD72LfVz9pNrrsvm0OFGxxGRWuiob0f+rNcds9VC98PvGB2nVlGRJHIOroVZxCTYlh/ZFPoAuS6+xgYSkVprXeNRFGGm6Ym1mmDyElKRJHIOnY/OwrPgBCfcw9jR4E6j44hILXbSozE7g28H4JqDb2CyWgxOVDuoSBIphc+ZI0QlzQdgbcQYTRwpIobbGPYguc7e1M/ZT5tj3xgdp1ZQkSRSiqsOvYXZWsgh3y4crNfN6DgiIuS6+LIx9EEAuh2ehlthpsGJaj4VSSL/I+zUZi47sYYizKyJGKOJI0WkytjRoD9pHk2oU5hBl4QPjY5T46lIEvkHk7WQqw++AcCOBndwwqOJwYlERP5mNTmzJmIsAO2Tv8Av54DBiWo2FUki/xCZ8hUBOQc44+zDxtAHjI4jIlJCgm9n/vC7GicsXH3wTbBajY5UY6lIEvmLe8EpuiXMAGwDJPNcfAxOJCJSurWNR1NocqHxqU00ObnO6Dg1lookkb9cefg93AszOe5xuf1WWxGRqiijTiN+DbkbgKsPvom5KM/gRDWTiiQRIOj0b/Zbalc1eRKrqczLGoqIGGJLo/vIcgnAN/coHRM/NjpOjaQiScRaxLUHXsWElT31+5LoE2V0IhGRCypw9mRtxGgAOh2dg3duorGBaiAVSVLrtTn2DcFZe8gze7Ku8WNGxxERKbO9Ab1I8OmIc1Ee1xx43eg4NY6KJKnV3AtOceXh9wDbYO0c1wCDE4mIOMBkYlWTp7CYzDQ9uY4mJ9YanahGUZEktVq3w9OpU5hBmkdT4hrcZXQcERGHnfCI4NeQewC45sDrmC25BieqOQwvkqZNm0ZERATu7u5ER0ezbt35b2Vcs2YN0dHRuLu706RJE2bMmFHs53PmzMFkMpXYcnOL/9I4el6pgRJjiTz2FQA/N3lKg7VFpNraHDqc066B+OQl0enoHKPj1BiGFkkLFy5k9OjRjB8/nu3bt9O9e3f69OlDQkJCqe0PHjxI37596d69O9u3b+fZZ59l1KhRLFq0qFg7b29vkpOTi23u7u7lPq/UQJZC+HY0JqzE17+BRJ8ORicSESm3ArMHqyP+DUDHxHn4ntH3WUUwWa3GTdXZuXNnOnTowPTp0+37WrZsSb9+/Zg8eXKJ9k8//TRLliwhPj7evm/kyJHs2LGDjRs3AraepNGjR3Pq1KkKOy9AXl4eeXl/z0ORmZlJaGgoGRkZeHt7l/k9SxWx8T1Y/iy5zt7MifqCM65+RicSEbk4Viu37RlF41ObOOzTicWt3y229uSY65sZGK7qyMzMxMfHp0zf34b1JOXn5xMbG0uvXr2K7e/VqxcbNmwo9TUbN24s0b53795s27aNgoIC+76srCzCw8Np1KgRN910E9u3b7+o8wJMnjwZHx8f+xYaGlrm9ypVzKkj8PN/AVgX/pgKJBGpGUwmfm7yFIVOboRnbKHF8R+MTlTtGVYkpaWlYbFYCAoKKrY/KCiIlJSUUl+TkpJSavvCwkLS0tIAaNGiBXPmzGHJkiXMnz8fd3d3unXrxv79+8t9XoBx48aRkZFh344cOeLwe5YqwGqFpU9CQTaEdWV30C1GJxIRqTAZdULZ1Gg4YJuJ273glLGBqjnDB26b/tEVCGC1Wkvsu1D7f+7v0qULgwcPpl27dnTv3p3PP/+cZs2a8c4771zUed3c3PD29i62STX0+3ew7wdwcoGbpoLJ8P8EREQqVGzDwaR5NMWj8BRXHXrL6DjVmmHfEAEBAZjN5hK9N6mpqSV6ec4KDg4utb2zszP+/v6lvsbJyYkrrrjC3pNUnvNKDZGbCUufsj3u9jgEtjA2j4hIJShycuHHps9ixUTr1O8IPbXV6EjVlmH3PLu6uhIdHc3KlSu57bbb7PtXrlzJrbfeWuprunbtyrffflts34oVK+jYsSMuLi6lvsZqtRIXF0dkZGS5zys1xKr/wukkqBcBVz1hdBoRkUqT7N2WHcH9aZ/yJdf9OZmP23/Gmyv3GR3LYUYPNjf0WsPYsWP56KOPmDVrFvHx8YwZM4aEhARGjhwJ2MYBDR061N5+5MiRHD58mLFjxxIfH8+sWbOYOXMmTzzx9xfeiy++yPLlyzlw4ABxcXEMHz6cuLg4+zHLcl6pgRI2w+b3bY9vegNc6hibR0Skkv0S/ghZLgHUyz1C56OzjY5TLRk6e96AAQNIT09n4sSJJCcn06ZNG5YuXUp4eDgAycnJxeYuioiIYOnSpYwZM4b33nuPkJAQ3n77bfr3729vc+rUKR588EFSUlLw8fEhKiqKtWvX0qlTpzKfV2qYglz45hHACu0HQ9NrjU4kIlLp8p3rsqrpk9z8+9N0TJzLvoCepHlebnSsasXQeZKqM0fmWRCD/TgB1r8JdYPgkc1Qp579R9Wx+1lEpMysVm7+/SkuO7GaY54tWNB2NkVO1Wd1gcq43FYt5kkSuSSStsMvb9se3/RmsQJJRKTGM5n4qenT5Dp7E5T9Ox0T5xmdqFpRkSQ1V2E+fPMoWC3Q+nZocaPRiURELrkc1wBW/bVkSecjH+Gf86fBiaoPFUlSc61/E47tBg9/6Pua0WlERAzze/0+HKh3Jc7WAnrtn4jJWmh0pGpBRZLUTMd+g7V/FUZ9XgXPAGPziIgYyWTix6bjyDXXJThrD9GJnxqdqFpQkSQ1T2E+LH4IigqgeV9o0//CrxERqeGy3QJZEzEWgK4JH1Av55CxgaoBFUlS86z5Pzi2C+r4/bX0yLmXmxERqU32BN7EwXoxOFvzuWH/CzgV6bLb+ahIkprlyBbbWCSAm6eCl5aaERGxM5n4semz5Jq9CM7aQydNMnleKpKk5sjPhq8eAmsRtB0IrbTMjIjI/8pyC+Lnpk8D0PnITIJO/2ZwoqpLRZLUHCv/AycOgHdD6POK0WlERKqsvfV783tAb5yw0Gfff3C2nDE6UpWkIklqhj9+hK0f2R7f+h7U8TU0johIVfdz06c47RpIvdwErjr0ltFxqiQVSVL9ZafbJo0E6PQQNO1hbB4RkWogz9mbFZf/B4B2KYtofPIXgxNVPSqSpHqzWm2L155OhoBm0HOC0YlERKqNBN/O/NpgIAC99r9EnYKTBieqWlQkSfW25UPY9wOYXeGOWeDqYXQiEZFqZX34I6TXicCzIJ1e+1+0/fEpgIokqc5SdsOK52yPe70MwZHG5hERqYYsZne+bz6JQpMrTU7+QlTSfKMjVRkqkqR6ys+BL+8HSx40uwE6PWh0IhGRaivd8zLWRIwBoPvhdwg6vcfgRFWDiiSpnpaPg7S9UDcYbp2mWbVFRC7SzuD+7PfvgdlaSN9943EtzDI6kuFUJEn189vXEDsHMMHt74Onv8GBRERqAJOJlZc9R6ZbML65R7n2z/+r9eOTVCRJ9ZL2x9+3+185BppcY2gcEZGaJM/Zm6XNXqYIMy3TltMq9VujIxlKRZJUH/k58PlQyD8N4d2gx3ijE4mI1DjJ3u3YEP4QANcdeJWA7H0GJzKOiiSpHqxW+H4spP4GnoG22/3NzkanEhGpkbY2HMbBejE4F+Vx8+9P41Z42uhIhlCRJNXDr/Ngx3wwOdkKJK9goxOJiNRcJieWXf4iGW4N8M09Su/9E2yLh9cyKpKk6kuKg6VP2h5f+zxEdDc0johIbZDr4st3LV6h0ORC0xNruSJxrtGRLjkVSVK15ZywjUM6Ox9St9FGJxIRqTVS67ZkVRPbH6kxh2cQemqLwYkuLRVJUnVZCuGLe+HUYfANh37TwUm/siIil9LuoH7sDrwZJ4rou+856ualGB3pktE3jlRdK5+Hg2vAxRMGzQcPP6MTiYjUPiYTPzd5imOezfEoOMkt8U/gbDljdKpLQkWSVE3bP4VN02yPb5sOQa2NzSMiUotZzO582+I1clzqEZS9l977J9aKiSZVJEnVc3QbfDfa9vjqp6HVrYbGEREROO3egO+av4LF5Eyz9B/pfHSm0ZEqnYokqVoyk2HBPWDJh+Y3wtXPGJ1IRET+kugTxc9NngYgJuF9Lkv/2eBElUtFklQd+dkwfyBkpUD9lrZ12TRQW0SkStkd3I/tDQYAcMO+F2r0jNz6BpKqocgCXw6H5Djw8IdBn4Gbl9GpRESkFGsiRnPYtzMuRbncumcsnnnHjY5UKVQkSdWw/FnY9wOY3WDgfPBrYnQiERE5B6vJme+bT+JEnXC884/RL340LoXZRseqcCqSxHibZsDmGbbHt78PYZ2NzSMiIheU5+zN162mku3iR2D2Pm7cOw6nokKjY1UoFUlirN+XwrK/Bmf3nACtbzM0joiIlF2GeyO+afkGBU5uRJzayLUH/q9GTQ2gIkmMc2QrLBoOWKHDMC05IiJSDR3zas3S5pMowonIY99wxdE5RkeqMCqSxBip8fDpHVCQA5f1hBtfB5PJ6FQiIlIOB/yuYlWTJwC4MmEarY59a3CiiqEiSS69k4fh49sg9xQ0ugLumgdmF6NTiYjIRdjZ4E62NhwKwPV/vEzT9FUGJ7p4KpLk0so6Dh/3g9PJtrmQ7v4cXD2NTiUiIhVgffij7Aq81bYY7t7xhJ3abHSki6IiSS6d3Ez45HY4cQB8wmDIYi1aKyJSk5hM/HTZOPb5X4eztYCb458k+PQuo1OVm+FF0rRp04iIiMDd3Z3o6GjWrVt33vZr1qwhOjoad3d3mjRpwowZM4r9/MMPP6R79+7Uq1ePevXq0bNnT7Zs2VKszYQJEzCZTMW24ODgCn9v8g95WfDZXZCyEzwCYOjX4B1idCoREalgVpOZZc0mcsi3C65FZ7jtt8fxz/7D6FjlYmiRtHDhQkaPHs348ePZvn073bt3p0+fPiQkJJTa/uDBg/Tt25fu3buzfft2nn32WUaNGsWiRYvsbVavXs2gQYNYtWoVGzduJCwsjF69epGYmFjsWK1btyY5Odm+7dpVfSvdKi8/Gz4bAAkbwc0HBi8C/6ZGpxIRkUpicXLl2xavkuTVFnfLafr/9gh+OQeNjuUwk9Vq3IQGnTt3pkOHDkyfPt2+r2XLlvTr14/JkyeXaP/000+zZMkS4uPj7ftGjhzJjh072LhxY6nnsFgs1KtXj3fffZehQ20DyiZMmMDXX39NXFxcubNnZmbi4+NDRkYG3t7e5T5OjZefA/MHwMG14OYNQ76GRtFGp7J7c2XNXXNIRMRoboWnuWP3SAKz95Ht4scXbd7npEfjMr9+zPXNKjyTI9/fhvUk5efnExsbS69evYrt79WrFxs2bCj1NRs3bizRvnfv3mzbto2CgoJSX5OTk0NBQQF+fsXHvuzfv5+QkBAiIiIYOHAgBw4cOG/evLw8MjMzi21yAQW5sOBuW4HkWtfWg1SFCiQREalcec5eLGr9Hqmel+NZcII7dz9EvZxDRscqM8OKpLS0NCwWC0FBQcX2BwUFkZKSUuprUlJSSm1fWFhIWlpaqa955plnaNiwIT179rTv69y5M/PmzWP58uV8+OGHpKSkEBMTQ3p6+jnzTp48GR8fH/sWGhpa1rdaOxXkwsLBcGAVuHjCPV9CaCejU4mIyCWW6+LLotbTqmWhZPjAbdP/TCBotVpL7LtQ+9L2A7z66qvMnz+fxYsX4+7ubt/fp08f+vfvT2RkJD179uT7778HYO7cuec877hx48jIyLBvR44cufCbq63ODtL+YyU414F7PofwrkanEhERg5RWKFWHMUqGFUkBAQGYzeYSvUapqakleovOCg4OLrW9s7Mz/v7+xfZPmTKFSZMmsWLFCtq2bXveLJ6enkRGRrJ///5ztnFzc8Pb27vYJqU4c8o2UeTBNbZLbPd8AY2vNDqViIgY7H8Lpbt2PUDQ6T1Gxzovw4okV1dXoqOjWblyZbH9K1euJCYmptTXdO3atUT7FStW0LFjR1xc/p6x+bXXXuOll15i2bJldOzY8YJZ8vLyiI+Pp0GDBuV4J2KXnQZzb4KjW8DdF4Z+AxHdjU4lIiJVxNlCKaVuK+oUZtD/t3/RMCPW6FjnZOjltrFjx/LRRx8xa9Ys4uPjGTNmDAkJCYwcORKwXeI6e0ca2O5kO3z4MGPHjiU+Pp5Zs2Yxc+ZMnnjiCXubV199leeee45Zs2bRuHFjUlJSSElJISsry97miSeeYM2aNRw8eJDNmzdzxx13kJmZybBhwy7dm69pMpNgdh9I2QWe9eHe76HRhQtUERGpXXJdfPmy9TQSfDriZsnm9j2PE3Hi/HMkGsXQImnAgAFMnTqViRMn0r59e9auXcvSpUsJDw8HIDk5udicSRERESxdupTVq1fTvn17XnrpJd5++2369+9vbzNt2jTy8/O54447aNCggX2bMmWKvc3Ro0cZNGgQzZs35/bbb8fV1ZVNmzbZzysOOr4XZvaCtH3g3RDu+wGC2xidSkREqqgCZ0++bjWVP/2uwrkoj5t/f5Lmx5cZHasEQ+dJqs40T9JfDm+A+QMhNwP8msKQr6Be9Sk2NU+SiIhxnIoK6fXHRFoe/wGAdeGPsa3hEPjrZqxaO0+S1AC/fQXzbrUVSI06wfCV1apAEhERYxU5ObPs8glsbzAAgO6H3+HaA69gshYanMxGRZKUz8b34It7wZIPLW6CYUvA0/+CLxMRESnG5MTqJk+wOmIMVky0S1nErfH/xqUw2+hkKpLEQYX58O1oWP6s7Xmnh+CueeBSx9BYIiJSvW0PuZtvW7xCgZMbESc3cNfuByEz2dBMKpKk7LLTbJfXYmcDJuj1MvR5BZzMRicTEZEa4E//HnzZZgY5LvUIzN4Hn94BRUWG5XE27MxyXlVtQHFA9j5uiX8Cn7xk8syeLG32Xw5ld4Mfzz0Bp4iIiKNSvNqwoO0sbo5/ivp9XgEn4/pz1JMkF3R52k8M3Dkcn7xkTrqHsqDtbA75dTM6loiI1FAZ7o34tP3Hhq/YoJ4kOSenogK6H3qbDskLADjs25nvm08iz7kWT3kgIiKXhNVk/FAOFUlSKq+8FG78fRwNsnYDsK3hENaH/wurSb8yIiJSO+gbT0qIOLGe3vsnUKcwg1yzF8svf4ED/lcbHUtEROSSUpEkduaiPLodnkZ00mcApNRtxffNJ5PpHmJwMhERkUtPRZIAtrvX+uz7DwE5fwKwvcFdrGv8OBYnV4OTiYiIGENFUi1nslrokPgZMQnTcbYWkO3ix8rLnuegn7F3FIiIiBhNRVIt5nsmgev/eJlGmdsB+MPvan68bDxnXOoZnExERMR4KpJqIZO1kOjET+l65EOci/LId6rD6ib/5rfAW+wrL4uIiNR2KpJqmfpZe7n+j5cIyt4L2OY++rHpODLdGxqcTEREpGpRkVRLuBZm0eXIh0QlLcQJC7nO3qyJGMOe+jeq90hERKQUKpJqOmsRrY4v5cpD7+BZcAKAff49WdXkCXJc/Q0OJyIiUnWpSKrBArPi6XHgNUJO7wLghHsYq5v8m8P1YgxOJiIiUvWpSKqBvHOTiEmYTsvjywDId6rD5tAR/BoyiCInF4PTiYiIVA8qkmoQ94JTdDo6m3bJX+BsLQAgvv4NrAt/jGy3QIPTiYiIVC8qkmoAt8LTtE9eSHTiJ7hZsgFI8LmCdY0fI7VuS4PTiYiIVE8qkqoxt8JMopIWEJU0H3dLFgCpnpezLnwUCb6dddeaiIjIRVCRVA3VyT9B++SFRCUvtPccpdeJYFPoCPYF9ASTk8EJRUREqj8VSdVIvZxDdEj6jFap3+NszQcgzaMpm0JHsN//WhVHIiIiFUhFUlVntRKasY2opPk0PbnOvju5bmu2NRzKH/7XqDgSERGpBCqSqii3wkxapX5H25TF+J05DIAVE3/6XUVsw8EkebXTmCMREZFKpCKpqknZDZum88DOL3ApygMg38mD+MA+bG8wkJMejY3NJyIiUkuoSKpqjm6BuE9wAY57XM7O4P7E17+BAmdPo5OJiIjUKiqSqprIuyAxlgWF15Ds1VaX1ERERAyiEb9VjVtduPU9kr015khERMRIKpJERERESqEiSURERKQUKpJERERESqEiSURERKQUKpJERERESqEiSURERKQUKpJERERESqEiSURERKQUKpJERERESmF4kTRt2jQiIiJwd3cnOjqadevWnbf9mjVriI6Oxt3dnSZNmjBjxowSbRYtWkSrVq1wc3OjVatWfPXVVxd9XhEREaldDC2SFi5cyOjRoxk/fjzbt2+ne/fu9OnTh4SEhFLbHzx4kL59+9K9e3e2b9/Os88+y6hRo1i0aJG9zcaNGxkwYABDhgxhx44dDBkyhLvuuovNmzeX+7wiIiJS+5isVqvVqJN37tyZDh06MH36dPu+li1b0q9fPyZPnlyi/dNPP82SJUuIj4+37xs5ciQ7duxg48aNAAwYMIDMzEx++OEHe5sbbriBevXqMX/+/HKdtzSZmZn4+PiQkZGBt7e3Y2+8DN5cua/CjykiIlKdjLm+WYUf05Hvb+cKP3sZ5efnExsbyzPPPFNsf69evdiwYUOpr9m4cSO9evUqtq93797MnDmTgoICXFxc2LhxI2PGjCnRZurUqeU+L0BeXh55eXn25xkZGYDtw64MudlZlXJcERGR6qIyvmPPHrMsfUSGFUlpaWlYLBaCgoKK7Q8KCiIlJaXU16SkpJTavrCwkLS0NBo0aHDONmePWZ7zAkyePJkXX3yxxP7Q0NBzv0kREREpt2cr8dinT5/Gx8fnvG0MK5LOMplMxZ5brdYS+y7U/n/3l+WYjp533LhxjB071v68qKiIEydO4O/vf97XVWWZmZmEhoZy5MiRSrlkWFvpc608+mwrjz7byqPPtvKU57O1Wq2cPn2akJCQC7Y1rEgKCAjAbDaX6L1JTU0t0ctzVnBwcKntnZ2d8ff3P2+bs8csz3kB3NzccHNzK7bP19f33G+wGvH29tZ/uJVAn2vl0WdbefTZVh59tpXH0c/2Qj1IZxl2d5urqyvR0dGsXLmy2P6VK1cSExNT6mu6du1aov2KFSvo2LEjLi4u521z9pjlOa+IiIjUPoZebhs7dixDhgyhY8eOdO3alQ8++ICEhARGjhwJ2C5xJSYmMm/ePMB2J9u7777L2LFjeeCBB9i4cSMzZ86037UG8Pjjj3PVVVfxyiuvcOutt/LNN9/w448/sn79+jKfV0RERASrwd577z1reHi41dXV1dqhQwfrmjVr7D8bNmyY9eqrry7WfvXq1daoqCirq6urtXHjxtbp06eXOOYXX3xhbd68udXFxcXaokUL66JFixw6b22Rm5trfeGFF6y5ublGR6lR9LlWHn22lUefbeXRZ1t5KvuzNXSeJBEREZGqyvBlSURERESqIhVJIiIiIqVQkSQiIiJSChVJIiIiIqVQkSQcOnSI4cOHExERQZ06dWjatCkvvPAC+fn5RkerlqZNm0ZERATu7u5ER0ezbt06oyNVe5MnT+aKK67Ay8uLwMBA+vXrx969e42OVeNMnjwZk8nE6NGjjY5SYyQmJjJ48GD8/f3x8PCgffv2xMbGGh2r2issLOS5556zf281adKEiRMnUlRUVKHnMXxZEjHe77//TlFREe+//z6XXXYZu3fv5oEHHiA7O5spU6YYHa9aWbhwIaNHj2batGl069aN999/nz59+rBnzx7CwsKMjldtrVmzhkceeYQrrriCwsJCxo8fT69evdizZw+enp5Gx6sRtm7dygcffEDbtm2NjlJjnDx5km7dutGjRw9++OEHAgMD+fPPP2vMag1GeuWVV5gxYwZz586ldevWbNu2jfvuuw8fHx8ef/zxCjuPpgCQUr322mtMnz6dAwcOGB2lWuncuTMdOnRg+vTp9n0tW7akX79+TJ482cBkNcvx48cJDAxkzZo1XHXVVUbHqfaysrLo0KED06ZN4+WXX6Z9+/ZMnTrV6FjV3jPPPMMvv/yi3uRKcNNNNxEUFMTMmTPt+/r374+Hhwcff/xxhZ1Hl9ukVBkZGfj5+Rkdo1rJz88nNjaWXr16Fdvfq1cvNmzYYFCqmikjIwNAv6MV5JFHHuHGG2+kZ8+eRkepUZYsWULHjh258847CQwMJCoqig8//NDoWDXClVdeyU8//cS+ffsA2LFjB+vXr6dv374Veh5dbpMS/vzzT9555x1ef/11o6NUK2lpaVgslhILJQcFBZVYUFnKz2q1MnbsWK688kratGljdJxqb8GCBfz6669s3brV6Cg1zoEDB5g+fTpjx47l2WefZcuWLYwaNQo3NzeGDh1qdLxq7emnnyYjI4MWLVpgNpuxWCz897//ZdCgQRV6HvUk1WATJkzAZDKdd9u2bVux1yQlJXHDDTdw5513MmLECIOSV28mk6nYc6vVWmKflN+jjz7Kzp07i63ZKOVz5MgRHn/8cT755BPc3d2NjlPjFBUV0aFDByZNmkRUVBQPPfQQDzzwQLHL8VI+Cxcu5JNPPuGzzz7j119/Ze7cuUyZMoW5c+dW6HnUk1SDPfroowwcOPC8bRo3bmx/nJSURI8ePeyL/opjAgICMJvNJXqNUlNTS/QuSfk89thjLFmyhLVr19KoUSOj41R7sbGxpKamEh0dbd9nsVhYu3Yt7777Lnl5eZjNZgMTVm8NGjSgVatWxfa1bNmSRYsWGZSo5njyySd55pln7N9xkZGRHD58mMmTJzNs2LAKO4+KpBosICCAgICAMrVNTEykR48eREdHM3v2bJyc1MnoKFdXV6Kjo1m5ciW33Xabff/KlSu59dZbDUxW/VmtVh577DG++uorVq9eTUREhNGRaoTrrruOXbt2Fdt333330aJFC55++mkVSBepW7duJaaq2LdvH+Hh4QYlqjlycnJKfE+ZzWZNASAVLykpiWuuuYawsDCmTJnC8ePH7T8LDg42MFn1M3bsWIYMGULHjh3tPXIJCQmMHDnS6GjV2iOPPMJnn33GN998g5eXl723zsfHhzp16hicrvry8vIqMa7L09MTf39/jfeqAGPGjCEmJoZJkyZx1113sWXLFj744AP11FeAm2++mf/+97+EhYXRunVrtm/fzhtvvMH9999fsSeySq03e/ZsK1DqJo577733rOHh4VZXV1drhw4drGvWrDE6UrV3rt/P2bNnGx2txrn66qutjz/+uNExaoxvv/3W2qZNG6ubm5u1RYsW1g8++MDoSDVCZmam9fHHH7eGhYVZ3d3drU2aNLGOHz/empeXV6Hn0TxJIiIiIqXQwBMRERGRUqhIEhERESmFiiQRERGRUqhIEhERESmFiiQRERGRUqhIEhERESmFiiQRERGRUqhIEhERESmFiiQRERGRUqhIEhERESmFiiQRERGRUqhIEhEBjh8/TnBwMJMmTbLv27x5M66urqxYscLAZCJiFC1wKyLyl6VLl9KvXz82bNhAixYtiIqK4sYbb2Tq1KlGRxMRA6hIEhH5h0ceeYQff/yRK664gh07drB161bc3d2NjiUiBlCRJCLyD2fOnKFNmzYcOXKEbdu20bZtW6MjiYhBNCZJROQfDhw4QFJSEkVFRRw+fNjoOCJiIPUkiYj8JT8/n06dOtG+fXtatGjBG2+8wa5duwgKCjI6mogYQEWSiMhfnnzySb788kt27NhB3bp16dGjB15eXnz33XdGRxMRA+hym4gIsHr1aqZOncrHH3+Mt7c3Tk5OfPzxx6xfv57p06cbHU9EDKCeJBEREZFSqCdJREREpBQqkkRERERKoSJJREREpBQqkkRERERKoSJJREREpBQqkkRERERKoSJJREREpBQqkkRERERKoSJJREREpBQqkkRERERKoSJJREREpBT/D1zQPf81fSi+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ce44cb9-35ec-4dec-bb5d-aa9bc0d75f03",
   "metadata": {},
   "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": "code",
   "execution_count": 10,
   "id": "1c3036f5-f214-48ee-a2ff-93009d010435",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.425"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#your code here\n",
    "\n",
    "#define function for Bayes' theorem\n",
    "def bayesTheorem(pA, pB, pBA):\n",
    "    return _______\n",
    "\n",
    "#define probabilities\n",
    "pRain = ____\n",
    "pCloudy = ____\n",
    "pCloudyRain = _____\n",
    "\n",
    "#use function to calculate conditional probability\n",
    "bayesTheorem(pRain, pCloudy, pCloudyRain)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dd0cfc0-be5c-4987-babc-93735aba1f32",
   "metadata": {},
   "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",
   "execution_count": 11,
   "id": "9890d5d6-2632-411a-b38a-dfae9e9d6c79",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "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",
      "\n",
      "\n",
      "Posterior probability of having the drug given a positive test: 0.16102\n"
     ]
    }
   ],
   "source": [
    "# Defining the prior probabilities and test accuracies\n",
    "prior_drug = _____  # What's the chance someone from the population has the drug?\n",
    "prior_no_drug = 1 - _____  # 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 = _____  # If a person has the drug, what's the chance the test detects it?\n",
    "p_pos_given_no_drug = _____  # 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 * _____  # Multiply the prior and the likelihood for people who have the drug\n",
    "unnormalized_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_no_drug  # Total probability of getting a positive test result\n",
    "\n",
    "# Calculate the posterior probabilities\n",
    "posterior_drug = unnormalized_drug / _____  # Given a positive test, what's the updated chance of truly having the drug?\n",
    "posterior_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": {
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}