jtrive84/censored-mle.ipynb

## censored-mle.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum Likelihood Estimation with Right-Censored Observations\n",
    "\n",
    "\n",
    "Maximum likelihood estimation (MLE) is a method used in statistics to estimate the parameters of a statistical model. The basic idea behind MLE is to find the values of the model parameters that maximize the likelihood function, which measures how likely the observed data are under the given statistical model. The values of the parameters that maximize the likelihood function are considered the maximum likelihood estimates. These estimates represent the most likely values of the parameters given the observed data.\n",
    "\n",
    "> *Likelihood is a measure of the extent to which a sample provides support for particular values of a parameter in a parametric model.*\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "### Maximum Likelihood Estimation for Non-Censored Data\n",
    "\n",
    "In practice, we do not know the values of the proposed model's parameters, but we do know the data. We use the likelihood function to observe how the function changes for different parameter values while holding the data fixed. Larger values of the likelihood correspond to values of the parameters that are relatively better supported by the data.\n",
    "\n",
    "\n",
    "   \n",
    "The joint density of $n$ independently distributed observations $\\mathbf{x} = (x_{1}, \\cdots, x_{n})^{T}$ is given by:\n",
    "\n",
    "$$\n",
    "f(\\mathbf{x}|\\mathbf{\\theta}) = \\prod_{i=1}^{n} f_{i}(x_{i}|\\mathbf{\\theta}),\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "where $f(\\mathbf{x}|\\mathbf{\\theta})$ is the product of univariate density functions.\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "When this expression is interpreted as a function of unknown $\\theta$ given known data $x$, we obtain the likelihood function:\n",
    "\n",
    "$$\n",
    "L(\\mathbf{\\theta}|\\mathbf{x}) = \\prod_{i=1}^{n} f_{i}(x_{i}|\\mathbf{\\theta})\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "Overall, the goal is to find the model parameters that maximize the likelihood function over the parameter space:\n",
    "\n",
    "$$\n",
    "\\hat {\\theta } = {\\underset {\\theta \\in \\Theta }{\\operatorname {arg\\;max} }}\\, L(\\theta |\\mathbf {x})\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "It is often convenient to work with the natural logarithm of the likelihood function, called the *log-likelihood*:\n",
    "\n",
    "$$\n",
    "\\mathcal{l}(\\mathbf{\\theta}|\\mathbf{x}) = \\mathrm{Ln}(L(\\mathbf{\\theta}|\\mathbf{x})) = \\sum_{i=1}^{n} f_{i}(x_{i}|\\mathbf{\\theta}).\\\\\n",
    "$$\n",
    "             \n",
    "<br>   \n",
    "\n",
    "Because the logorithm is monotonic, the maximum of $\\mathcal{l}(\\mathbf{\\theta}|\\mathbf{x})$ occurs at the same value of $\\theta$ as does the maximum of $L(\\mathbf{\\theta}|\\mathbf{x})$.\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "### MLE Estimation: Numerical Example\n",
    "\n",
    "Assume we have 5 losses with the following values:\n",
    "\n",
    "> 816 580 3288 12494 4587\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "We want to find the best fitting (in the MLE sense) exponential distribution to this data in order to estimate quantiles for some downstream analysis. Recall the density of the exponential distribution is given by:\n",
    "\n",
    "$$\n",
    "f(x; \\theta) = \\frac{1}{\\theta}e^{-x/\\theta}, \\hspace{.50em} x \\geq 0.\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "Expanding the definition of likelihood above, we obtain:\n",
    "\n",
    "$$\n",
    "\\begin{align*} \n",
    "L(\\mathbf{\\theta}|\\mathbf{x}) &= \\prod_{i=1}^{n} f_{i}(x_{i}|\\mathbf{\\theta})\\\\\n",
    "&= \\prod_{i=1}^{5} f_{i}(x_{i}|\\mathbf{\\theta})\\\\\n",
    "&= \\frac{1}{\\theta}e^{-x_1/\\theta} \\cdot \\frac{1}{\\theta}e^{-x_2/\\theta} \\cdot \\frac{1}{\\theta}e^{-x_3/\\theta} \\cdot \\frac{1}{\\theta}e^{-x_4/\\theta} \\cdot \\frac{1}{\\theta}e^{-x_5/\\theta}\\\\\n",
    "&= \\frac{1}{\\theta}e^{-816/\\theta} \\cdot \\frac{1}{\\theta}e^{-580/\\theta} \\cdot \\frac{1}{\\theta}e^{-3288/\\theta} \\cdot \\frac{1}{\\theta}e^{-12494/\\theta} \\cdot \\frac{1}{\\theta}e^{-4587/\\theta}\\\\\n",
    "&= \\Big(\\frac{1}{\\theta}\\Big)^{5} \\cdot e^{-21765/\\theta}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "Taking the natural log of the last expression yields:\n",
    "\n",
    "$$\n",
    "\\mathcal{l}(\\mathbf{\\theta}|\\mathbf{x}) = 5 \\cdot \\mathrm{Ln}\\Big(\\frac{1}{\\theta}\\Big) - 21765/\\theta.\n",
    "$$\n",
    "\n",
    "<br>\n",
    "In order to find the maximum likelihood estimate, we take the derivative of the log-likelihood, set it equal to 0 and solve for theta:\n",
    "\n",
    "$$\n",
    "\\frac{d\\mathcal{l}}{d\\theta} = \\frac{-5}{\\theta} + \\frac{21765}{\\theta^2}\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "Setting equal to 0, rearranging and solving for $\\theta$ results in:\n",
    "\n",
    "$$\n",
    "\\hat \\theta = 4353.\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "Therefore, the best fitting exponential distribution to this data in an MLE sense is:\n",
    "\n",
    "$$\n",
    "f_{MLE}(x) = \\frac{1}{4353}e^{-x / 4353}.\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "Note that some distributions have shortcuts for determining maximum likelihood estimates. For the exponential, the maximum likelihood estimate of $\\theta$ is identical to the average of the underlying data:\n",
    "\n",
    "\n",
    "$$\n",
    "\\hat \\theta = \\frac{816 + 580 + 3288 + 12494 + 4587}{5} = 4353.\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "A general purpose solver can be implemented in Scipy as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "+ MLE scale parameter estimate: 4,353.\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Using Scipy optimization routines to find maximum likelihood estimates.\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "np.set_printoptions(suppress=True, precision=8)\n",
    "\n",
    "\n",
    "x = [816, 580, 3288, 12494, 4587]\n",
    "\n",
    "\n",
    "def loglik(scale):\n",
    "    \"\"\"\n",
    "    Exponential log-likelihood.\n",
    "    \"\"\"\n",
    "    n, sum_x = len(x), sum(x)\n",
    "    return -(n * np.log(1 / scale) - sum_x / scale)\n",
    "\n",
    "\n",
    "# Maximize loglikelihood.\n",
    "results = minimize(loglik, 816, method=\"Nelder-Mead\")\n",
    "\n",
    "scale = results.x.item()\n",
    "\n",
    "print(f\"\\n+ MLE scale parameter estimate: {scale:,.0f}.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the original data with the best fitting exponential distribution:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArIAAAHqCAYAAAD4TK2HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABSiklEQVR4nO3dd3yV5f3/8ffJyTnZe5JAwoawkSV7CYoKauuqu9+qtcNqtcOqrdZq9ddWi639fm3tcKLWgXuAA7eACAiyZSdAEkhyss+6fn8EjhwSIIkk59w5r+fjEeXc83OfK+OdO9d9XTZjjBEAAABgMVGhLgAAAABoD4IsAAAALIkgCwAAAEsiyAIAAMCSCLIAAACwJIIsAAAALIkgCwAAAEsiyAIAAMCSCLIAAACwJIIsAL3wwgvq2bNnqMto5oMPPlD37t2PuU3Pnj31wgsvdE5BFtKa9+5wNptNq1at6riC2uCKK67Q9ddfH+oyAFgAQRZAm2zfvl02m02VlZUdfq7Jkydr9+7dgdfTpk3T/PnzO/y8VtNSmxz53gFAV0SQBQCEDa/XG+oSAFgIQRaIQLt379bs2bOVnJysUaNGad26dUHr77vvPvXr109JSUnq06ePHnjggcC6sWPHSpK6d++uxMREPfHEE6qpqdFZZ52l7OxspaSkaMqUKVq9enWL5/Z6vUpKStKGDRskSS+//LJsNpveeOMNSdKaNWuUmpoqn8+nJUuWKDU1VZJ044036oMPPtAvf/lLJSYmas6cOYFjbtq0SSeffLKSkpI0depU7dq166jXXlpaqosvvljdunVTXl6err/+ejU2NkqSzjvvPF1yySWBbf/4xz9q8ODBqq+vD9z1fOihh9SzZ09lZGTohz/8odxud2D7RYsWaeTIkUpJSdFJJ52kt956K7Duiiuu0FVXXaULL7xQSUlJGjBggJYsWRJY7/F49Jvf/EZ9+vRRRkaG5s2bp5KSksB6m82mBx98UEOGDFFycrLmzZunqqqqo7bJ4e+dJD3++OMaMmSIkpKSVFBQoF//+tcyxhz1fTqcMUb33nuv+vTpo/T0dJ122mnaunWrJOnPf/6zZsyYEbT9008/rYEDBwZeP/XUUxo2bJhSU1M1ZswYffzxx4F106ZN0y9+8QvNnj1bCQkJev3115ud/5JLLlFeXl7g8/Xdd98NvGc5OTlB76MkFRUV6emnn27VtQGwOAMg4kyePNlcdtllpra21qxfv9707NnTFBYWBtY/++yzZufOncbv95t33nnHxMbGmg8//NAYY8y2bduMJFNRURHYvqqqyjz11FOmpqbG1NfXm5/85Cemf//+xu/3t3j+OXPmmP/93/81xhhz/fXXmz59+phf/OIXxhhj5s+fb+bNm2eMMebdd981KSkpgf2mTp1q/vznPwcdq7Cw0AwdOtRs3brV1NfXmzlz5pjLL7+8xfP6/X4zbtw4c8MNN5ja2lpTXl5upk2bZm699VZjjDEVFRWmoKDAPPLII2b58uUmOTnZfPHFF0HXPWfOHFNRUWGKi4vN8OHDze23326MMWbz5s0mNjbWPPfcc8bj8ZhnnnnGxMXFma1btxpjjLn88stNUlKSeffdd43X6zW/+93vgt7zn//852bGjBmmpKTENDY2mhtvvNFMnjw5sF6SmT59utm3b5+pqKgwI0eONLfddttR2+TI9+61114zGzduNH6/36xcudJkZ2ebxx9/POj4K1eubPF9e+SRR0xeXp754osvTH19vbnhhhvMoEGDjMfjMXv37jUOh8Ps3LkzsP0ZZ5xh7rzzTmOMMa+++qrJz883K1asMD6fzzz33HMmPT3dlJeXB9o0KyvLLF261Pj9flNXV2cuv/xyc9111wWO9+9//9tUVlYat9tt/vCHP5j09HTjcrmMMcbceOONQe398ccfm7S0NNPQ0NDitQDoWgiyQITZuXOnkWT27dsXWHbPPfcEhaojnXXWWYFg0lJoOlJFRYWRZHbv3t3i+j/84Q/mvPPOM8YYM2zYMPPoo4+aMWPGGGOMmTdvnpk/f74xpvVB9v/+7/8Crx9//HEzZMiQFs+7bNkyk56ebnw+X2DZokWLTO/evQOvP/jgA5Oammp69eplHnjggcDyQ9e9dOnSwLKnnnrK9OnTxxhjzJ133mlOO+20oPPNmjXL3HXXXcaYpiB7wQUXBNbt3r3bSDLl5eXG7/ebhIQEs2rVqsD6+vp6ExUVFQiIkszrr78eWH/nnXeaM888M6i2YwXZI1133XXmyiuvDLw+VpA95ZRTzD333BN43dDQYJKSksxHH31kjGn6xeTuu+82xhizb98+43Q6zY4dO4wxxpx++umB9jxkwoQJ5tFHHzXGNLXp4aH10Ht15LLDpaamBn6xWrdunUlMTDTV1dXGGGOuvvpq86Mf/eio+wLoWuhaAESYkpISxcbGKjs7O7CssLAwaJsnnnhCJ510ktLT05WamqrXXntN5eXlRz1mfX29fvjDH6pnz55KTk4OjIBwtH2mT5+uJUuWqKysTGVlZbrooou0fft2VVRU6P3332/2p+rjyc3NDfw7ISFB1dXVLW63fft2VVZWBq4rNTVV5557rvbt2xfYZuLEierdu7dcLpeuvPLKZsc4/L0qLCxUcXGxpKbuGkeO/NC7d++gB66OrFOSqqurVV5ertraWk2ZMiVQV25urpxOZ1A3idZeZ0vefPNNTZgwQZmZmUpJSdGDDz54zDY93JHXFhMTo7y8vMC1XXbZZXrsscckSU8++aQmTJiggoICSU3v+c033xy4rtTUVK1atSrwvkkKbNsSv9+vW265Rf369VNycrJSU1NVVVUVqL2oqEhDhgzRs88+q4aGBj399NP6n//5n1a/LwCsjSALRJi8vDw1NDSotLQ0sGznzp1B/7788sv1hz/8QaWlpaqsrNTpp58e6E8ZFdX828a9996rFStW6MMPP5TL5dL27dsl6ah9MEeOHCm3260HHnhAU6dOld1u16RJkzR//nw5HA4NGTKkxf1aOndb9OjRQ9nZ2aqsrAx8VFVVqaamJuhaGhsbVVRUpJtvvrnZMXbs2BH4986dO5Wfny+pqX/qoes+ZPv27a0aAisjI0Px8fFaunRpUG319fWaMGHCcfc/3vvidrv1rW99S9///vdVXFysqqoqXXPNNa3uI3vktbndbpWUlASu7ayzztLu3bu1YsUKPfbYY7r00ksD2/bo0UP33ntv0HXV1tbqpptualX9CxYs0IIFC/Tqq6+qqqpKlZWVSklJCar9e9/7nh5++GEtXLhQhYWFOumkk1p1XQCsjyALRJgePXpo4sSJuummm1RfX6+NGzfq73//e2B9TU2NjDHKzs5WVFSUXnvtNS1atCiwPisrS1FRUfrqq68Cy1wul2JjY5WWlqaampoWA+Dh7Ha7pkyZovnz52v69OmSpBkzZmj+/PmaNm2abDZbi/vl5OQEnbetxowZox49eujWW29VdXW1jDHasWNH4AGjFStW6He/+52efPJJLViwQA8//LDefPPNoGPccccdqqysVElJie6++25dfPHFkqQLLrhAS5Ys0Ysvviiv16vnn39e77//vi688MLj1hUVFaVrrrlGN954Y+AO7P79+1v9wFJLbXK4xsZGNTQ0KCMjQzExMVq6dKkWLFjQqmNLTQ9bPfDAA1q3bp0aGxt16623Kj8/P/CQWVxcnM4991zdcsstWrdunc4777zAvj/60Y/0xz/+UStWrJAxRnV1dXrrrbdaPTSYy+WS0+lUZmam3G637rjjjmZ3oi+44AKtWLFC99xzD3djgQhDkAUi0IIFC7Rr1y5lZ2froosuCvrhP2jQIN1yyy2aMWOGMjIy9PTTT2vevHmB9XFxcbrttts0Z84cpaamasGCBbrhhhtkt9uVk5OjIUOGaPz48cetYfr06XK5XIFuBDNnzgx63ZLrr79eb731llJTU3XmmWe2+brtdrteeeUVFRcXq6ioSCkpKTrjjDO0ZcsW1dTU6Dvf+Y7uuusuDR06VD169NBDDz2kyy+/POju9VlnnaURI0ZoyJAhGjduXCC09+3bV88//7xuu+02paen64477tDChQvVu3fvVtV29913a/z48ZoxY4aSkpI0atSooF8gjqWlNjlcUlKS/va3v+nqq69WcnKy7rrrLl1wwQWtfNeaug5ce+21OvPMM5Wbm6vVq1fr5ZdfVnR0dNA2b775ps4++2wlJSUFls+dO1f33HOPrrrqKqWlpalXr166//775ff7W3Xuyy+/XIMHD1ZhYaF69+6tuLi4Zne5k5KSdN5552nDhg2BXywARAabae3flgAggm3fvl29evVSRUVF0LBWCA933HGHvvjiCz377LOhLgVAJ4o+/iYAAISvsrIyPfTQQ3r44YdDXQqATkbXAgCAZd11113q2bOnzjjjDM2cOTPU5QDoZHQtAAAAgCVxRxYAAACWRJAFAACAJRFkAQAAYEltHrXA7/erpKRESUlJRx20HAAAAGgPY4yqq6uVl5d33JkL2xxkS0pK1KNHj3YXBwAAABzPrl27jjvNd5uD7KEZW7Zt26b09PT2VQZL8ng8WrRokWbPni2HwxHqctCJaPvIRdtHLto+coW67V0ul3r06BE0S+DRtDnIHupOkJSUpOTk5LZXB8vyeDyKj49XcnIy39QiDG0fuWj7yEXbR65wafvWdGHlYS8AAABYEkEWAAAAlkSQBQAAgCURZAEAAGBJBFkAAABYEkEWAAAAlkSQBQAAgCURZAEAAGBJBFkAAABYEkEWAAAAlkSQBQAAgCURZAEAAGBJ0e3d0ePxyOPxnMhaEOYOtTftHnlo+8hF20cu2j5yhbrt23JemzHGtOXgLpdLKSkpWrBggeLj49tcHAAAAHA0dXV1uuiii1RVVaXk5ORjbtvuO7LTp09XRkZGe3eHBXk8Hi1evFizZs2Sw+EIdTnoRLR95KLtIxdtH7lC3fYul6vV27Y7yDocDj6xIxRtH7lo+8hF20cu2j5yhart23JOHvYCAACAJRFkAQAAYEnt7lrQUXw+v/xte/4MnSDKZgt1CQCACHL++efrzTff1MaNG5WbmxvqchCmwi7I+o3RntIqeTy+UJeCgxwOu7plp4S6DABAhFi8eLHWrFmja665RjfccIMWLFgQ6pIQpsIuyEqSx+OTmyALAEDEaWxs1PXXX6/HHntMI0aM0Pjx4/XOO+9oxowZoS4NYSgsgywAAIhMMTEx+vLLLwOvly9fHsJqEO542AsAAACWRJAFAABhZfXq1Zo9e7YSEhKUn5+vv/zlL6Euqd1qamrUvXt32Ww2ffbZZ8fd/rXXXtPUqVOVlZWlmJgY9e7dWzfccIOqqqoC2zz88MOy2WzNPm666aY2H8vq6FoAAADCxieffKJZs2bpqquu0s0336xXXnlF1113nUaMGKEpU6aEurw2+93vfiev19vq7Q8cOKBx48bpJz/5iTIyMrR27VrdfvvtWrt2rRYtWhS07RtvvKGUlK8fxs7Pz2/3sayKIAsAAMJCY2OjLrroIn3/+9/XvffeK0maNm2aFi5cqKeeespyQXbDhg3629/+pnvvvVfXXHNNq/a55JJLgl5PmzZNMTExuvrqq1VSUqK8vLzAulGjRikzM/OEHMuq6FoAAADCwn/+8x+Vl5frtttuC1reo0cP7dixI0RVtd+1116ra665RgMGDPhGx8nIyJAkud3ub1zTiTxWOCDIAgCAsPDYY49p3rx5io+Pl9frDXzU1tbK4XB0Wh3GmKDzH+3jWJ599lmtWbNGv/nNb9pVg8/nU0NDgz7//HPdcccdmjdvnnr27Bm0zeDBg2W329W7d2/dfffd8vlaHrq0NceyKoIsAAAIuerqai1dulQLFiyQw+EI+vjss89UUFDQabU88sgjzWpo6WP79u0t7l9XV6cbbrhBv//975WcnNyuGgoLCxUXF6dRo0apW7duQZNCdOvWTb/97W/16KOP6vXXX9fpp5+uW2+9Vdddd12bj2V19JEFAAAh98UXX8jn8+mZZ54Julu4bds2nX/++Ro5cmSn1TJ37txWjV97tD6md955p3JycvTd73633TW89tprqq2t1Zdffqk777xTc+fO1eLFi2W323Xqqafq1FNPDWw7e/ZsxcXF6c9//rNuueUWdevWrdXHsjqCLAAACLldu3ZJanog6fAHmD777DPZbDbNnDlTkvTKK69o7ty5+sMf/qCf//znkqT9+/erb9++6tu3r5YvX64XXnhB8+fP15IlS4LO8corr+jSSy9VYWFhYNm3v/1t/frXvw7aLj09PWg0gKOJjm4eo3bs2KF7771XCxcuDAxzVVNTE/h/TU2NEhMTj3vsYcOGSZLGjx+vMWPGaMSIEVq4cKHOPffcFrc///zz9ac//UmrVq1qFmTbeiwrIcgCAICQO9Tn9Mi7hE888YSmTZsW6FqwcuVKjRw5UuvWrQts89vf/lY9e/bUiBEjJEmrVq0K/PtwK1eu1Pe//33dc889x6zlkUceadXd1G3btjXra7pt2za53W6dccYZzbafPn26xo0bp08//fS4xz7csGHD5HA4tGXLljbt19HHCgcEWQAAEHKH7pKuW7dOEydOlNR0B/Wjjz7Se++9F9hu1apVuuyyy/Tkk09KkjZv3qyPPvpII0aMCHQ/WLVqlc4+++xm51i1apUuvPDC49byTboWjBgxQu+++26z8/70pz/Vgw8+qDFjxhz3uEdaunSpPB6PevfufdRtnnrqKdnt9uN2wWjNsayEIAsAAEJu/PjxKigo0HXXXae7775bmzdv1k033aTbb79dkydPDmy3cuVK3X333YHZvn7xi1/o97//vW6++WZdddVVkpqC4+23397sHCtXrtT69et11113SZL69Omj5557rtl2GRkZgWGq2io1NVXTpk1rcd2oUaN00kknBV4vWbJE06dP13/+8x9dccUVkqRvfetbGj16tIYNG6a4uDitXr1af/zjHzVs2LBAOD/11FM1Y8YMDR06VJL00ksv6R//+Ieuu+465ebmBo7fmmNZHUEWAACEXHR0tJ5//nldddVVmjt3rnr37q37778/6E/8VVVVKisrU9++fVVYWKgnnnhCdXV1mjlzptavX6/hw4ersrJSe/bs0aBBg4KOX1VVpf3794fV9Ky1tbWSFBQ+x44dq6efflr33HOP/H6/evbsqauuuko/+9nP5HQ6JUkDBw7Uv/71L+3evVt+v1/9+/fX/Pnzde211wYdvzXHsjqCLAAACAujRo3S559/ftT1q1at0pAhQxQVFaVhw4bpxz/+sd577z2tX79ePXr0UEJCgpYsWaKBAwc2C2qrVq3SwIEDO/oSWjRt2jQZY5ot//TTTzV48OCgEQhuuukm3XTTTcc83v3336/777//uOdtzbGsjiALAAAs4fCHuC699FKNGTNGw4YN02OPPRbUP7alB71WrVqloqKiTqz2+D766CPdfPPNstlsoS7FsgiyAADAElauXKnx48dLkkaPHq3Ro0dLklavXh00YsGLL74YNJrAq6++qpUrV+r1118PCrl33nmnzjzzzM4qv5l33nknZOfuKgiyAADAEh5++OEWl//pT3867jZHWw5rY4paAAAAWBJBFgAAAJZEkAUAAIAlEWQBAABgSQRZAAAAWBJBFgAAAJZEkAUAAIAlEWQBAABgSQRZAAAAWBJBFgAAAJZEkAUAAIAlEWQBAABgSQRZAAAAWBJBFgAAAJZEkAUAAIAlEWQBAABgSQRZAAAAWBJBFgAAAJZEkAUAAIAlEWQBAABgSQRZAAAAWBJBFgAAAJZEkAUAAIAlEWQBAABgSQRZAAAAWBJBFgAAAJZEkAUAAIAlEWQBAABgSQRZAAAAWBJBFgAAAJZEkAUAAIAlEWQBAABgSQRZAAAAWFJ0e3f0eDzyeDwnspaDbDLGyBjTAcdGexhjJGPk9XolqYPaHeHsUJvT9pGHto9ctH3kCnXbt+W8NtPGxOhyuZSSkqIFCxYoPj6+zcUdS2JiosaMGauvdpapocF9Qo+N9ouNdapPQZaWL1+mmpqaUJcDAAC6sLq6Ol100UWqqqpScnLyMbdt9x3Z6dOnKyMjo727H4NNSYmJionxdcCx0R5Oh11xcXGaMGGCFi1apFmzZsnhcIS6LHQij8ejxYsX0/YRiLaPXLR95Ap127tcrlZv2+4g63A4OuTiPF6fbDabbDbbCT822sdms0k2m6Kjmz5dOqrtEf5o+8hF20cu2j5yhart23JOHvYCAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJRFkAQAAYEkEWQAAAFgSQRYAAACWRJAFAACAJUW3d0ePxyOPx3MiaznIJmOMjDEdcGy0hzFGMkZer1eSOqjdEc4OtTltH3lo+8hF20euULd9W85rM21MjC6XSykpKVqwYIHi4+PbXNyxJCYmasyYsfpqZ5kaGtwn9Nhov9hYp/oUZGn58mWqqakJdTkAAKALq6ur00UXXaSqqiolJycfc9t235GdPn26MjIy2rv7MdiUlJiomBhfBxwb7eF02BUXF6cJEyZo0aJFmjVrlhwOR6jLQifyeDxavHgxbR+BaPvIRdtHrlC3vcvlavW27Q6yDoejQy7O4/XJZrPJZrOd8GOjfWw2m2SzKTq66dOlo9oe4Y+2j1y0feSi7SNXqNq+LefkYS8AAABYEkEWAAAAlkSQBQAAgCURZAEAAGBJBFkAAABYEkEWAAAAlkSQBQAAgCURZAEAAGBJBFkAAABYEkEWAAAAlkSQBQAAgCURZAEAAGBJBFkAAABYEkEWAAAAlkSQBQAAgCURZAEAAGBJBFkAAABYEkEWAAAAlkSQBQAAgCURZAEAAGBJ7Q6ypdUNJ7IOAAAAoE3aHWSvemyl9lYRZgEAABAa7Q6yOw/U64J/fKLiyvoTWQ8AAADQKt+oj+yO/XW64O+faNeBuhNVDwAAANAq7Q6yBelxkqTdFfW68B+faud+wiwAAAA6T7uD7EOXjlTvrARJUnFlUzeD7eW1J6wwAAAA4FjaHWSzk2L11NUnq192oiRpT1WDLvjHJ9paVnPCigMAAACO5hv1kc1OitWTV5+sATlJkqR9rkbdu2jTCSkMAAAAOJZvPCFCZmKMnrz6ZBV1S9bowjT9v3OHnYi6AAAAgGOKPhEHSU9wasGV4xRttykx5oQcEgAAADimEzZFbVqCU0mxjqBlpa4GfbG78kSdAgAAAAg4YUH2SAdq3br4n0v1nX98qqVb93fUaQAAABChOizI3v/WJm0urVGt26fL/7NM720q66hTAQAAIAJ1WJD91elFmj4gS5LU4PHrykeW6421ezrqdAAAAIgwHRZkYx12/f3S0Tp9aK4kyeMz+tGClVq4cndHnRIAAAARpMOCrCQ5o6P0lwtH6tsndZck+fxGN/x3tZ5YuqMjTwsAAIAI0KFBVpKi7VH647nDdOnJhZIkY6RbFq7VQ+9v7ehTAwAAoAvr8CArSVFRNt1x1mBdM7VPYNldr63XuxtLO+P0AAAA6II6JchKks1m0y9PG6Cfze4vSTp3VHdN7ZfVWacHAABAF9Op03DZbDb9eEY/DcpL1tT+2YqKsnXm6QEAANCFdNod2cPNGJgj+xEhdsNelxo8vlCUAwAAAAsKSZA90sa91Tr/wU902b+WyVXvCXU5AAAAsICQB1mf3+gHT6yQq8GrZdsP6OJ/LdOBOneoywIAAECYC3mQtUfZNP+CEUpPcEqSNuyt1nUL12uPqzHElQEAACCchTzIStKw7ql65prxyk+NkySVuBp148sbtG1/XYgrAwAAQLgKiyArSX2yEvXsD8arb3aiJKmi3qufv7pJa/dWh7gyAAAAhKOwCbKS1C0lTguuHKuinARJUq3bp1te36ylOytDWxgAAADCTlgFWUlKi3fqj3MHalT3ZEmS22d051tbVVbDA2AAAAD4WtgFWUmKc9h126w+mto7TZJ05bjuykp0hrgqAAAAhJNOndmrLRz2KP1iei9N65OukwtTQ10OAAAAwkxY3pE9JMpmazHELttVxSxgAAAAES6sg2xLPtlRqd8u2qKbXtusSmYBAwAAiFiWCrINXr/++uEO+Y20saxWN768USVMnAAAABCRLBVkY6OjdNdp/ZQR75DUNHHCDS9t0May2hBXBgAAgM5mqSArSb0y4nXfvIEqSI2VJFU1ePXLVzdp2a6qEFcGAACAzmS5ICtJ2YlO/WnuAA3JbZoFrNHr128XbdGbG8tDXBkAAAA6S7uH3/J4PPJ4OuJhK5uMMTLGHHOrRKddd57aV396f7s+3FYpv5Hmf7BDJa4GXTYqT1E2WwfUFpmMMZIx8nq9ktRB7Y5wdqjNafvIQ9tHLto+coW67dtyXps5XmI8gsvlUkpKihYsWKD4+Pg2F3csiYmJGjNmrL7aWaaGhtbN5OU3Ro+v3q/XtzR1LYh3ROmeU7orK8FxQmuLZLGxTvUpyNLy5ctUU1MT6nIAAEAXVldXp4suukhVVVVKTk4+5rbtDrJ79uxRRkbGNyr0KCVpR/EBuds4TuxLX5bqX8uLdfusPhqZf+yLRts4HXYV5qfL6/Vo0aJFmjVrlhwOflGIJB6PR4sXL6btIxBtH7lo+8gV6rZ3uVzKzMxsVZBtd9cCh8PRIRfn8fpks9lka2PXgLOG5GhSrzRlJDCV7Ylms9kkm03R0U2fLh3V9gh/tH3kou0jF20fuULV9m05pyUf9jqaI0OsMUYPfrJLq0pcIaoIAAAAHaVLBdkjPb16r178slS3vr6ZEQ0AAAC6mC4bZP3GaENp00QJvoMjGjy8vFj+tnUJBgAAQJjqskE2ymbTr0/po7MGZweWPb16r/7fO9vU6PWHsDIAAACcCF02yEqSPcqma8b30DXjeyjq4LNj72+r0C9f3agDdYyLBwAAYGVdOsgectbgbP1mVh/FRjdd7sayOv3khfXaXF4b4soAAADQXhERZCVpXEGq7p07IDBRwv46j37+8katLqkOcWUAAABoj4gJspLUOyNe959VpEHZCZKkzESnemfEhbgqAAAAtEe7J0SwqrR4h+4+o7/+uXS35g3OVlJMxL0FAAAAXUJEpjinPUo/nFDQbPmBOo98fqOsRGYHAwAACHcR1bXgWNxev3731le67sX1Wl9aE+pyAAAAcBwE2YMeW1GiDaW1qqj36pevbtLbm/eHuiQAAAAcA0H2oHOH52pobqIkyeMz+tN72/WPT3fJ52cmMAAAgHBEkD0oJTZad83pp9MHZgaWLVxbqlve2KyqBm8IKwMAAEBLCLKHcdij9OOJBfrxxAJFH5wKbHVJtX7ywnptKa8LcXUAAAA4HEH2CDabTWcUZeme0/srLa5pUIfSGrdufHmD3tlCv1kAAIBwQZA9isG5ifrL2UUakNU0eYLbZ7grCwAAEEYIsseQmeDUH87sr9MGZGp4XpK+N7Z7qEsCAADAQRE5IUJbOO1R+smkArl9RvaD/WYPqXX7lOC0h6gyAACAyMYd2Vaw2WyKiQ5+q7aU1+mKp9cw3iwAAECIEGTboarBq9+99ZVqGn3603vb9bePdsrt84e6LAAAgIhCkG2HOEeUTspPDrx+ZX2ZfvHKRpXWuENYFQAAQGQhyLaD0x6l6yYX6rpJhXLYm/rNbiyr07UL12nF7qoQVwcAABAZCLLfwGkDM3Xf3IHKTXJKklyNPv36jS164vMS+Q1T2wIAAHQkguw31DczXn85u0jjClIkSUbS45/v0W/e3MLUtgAAAB2IIHsCJMVE6zez+uiK0Xk6NELX6pJq7a1uDG1hAAAAXRhB9gSJstl0wYhu+v2c/kqNjdb3T+4emBUMAAAAJx4TIpxgw/OS9PdzByspJniiBJ/fyOPzK9bBBAoAAAAnAndkO0BybLRstuBZwB5bUaJrX1ivbfvrQlQVAABA10KQ7QTLdlbp6dV7tbuqUde9tEGvrS+TYVQDAACAb4Qg2wnyU2LUJyNOkuTxGf31o526551tqnX7QlwZAACAdRFkO0F+SqzumzdQ8wZlBZa9v61CP164TpvKakNYGQAAgHURZDuJ0x6lH0wo0K0zeyvB2fTA195qt258eaNeWLuPrgYAAABtRJDtZBN7pemBc4o0ICtekuT1G/390926Y/FXdDUAAABoA4JsCOQmxeiPZw7Qt4fmBJbtr/PIYbcdYy8AAAAcjiAbIg57lK4c112/nd1XuUlO3TSjt5x2mgMAAKC1mBAhxMYWpGhU9yGyRwXfjd1RUa94h11Zic4QVQYAABDeuAUYBo4MsQ1ev+56a6t+uHCdPtxWEaKqAAAAwhtBNgw9uXKPdlU1qKbRp7ve3qr7P9ihBg8PggEAAByOIBuGzh2Wo8m90gKv39hYrh+/sF6byxlzFgAA4BCCbBhKionWr2b00k+nFCo2uqmJiqsadcNLG/XM6r3yM+YsAAAAQTZc2Ww2ze6fqQfOKVK/zK/HnP338mLd8vpmlde6Q1whAABAaBFkw1x+SqzunTtA5w/P1aFHwlaVVOsnL6xXPf1mAQBABCPIWoDDHqXvjsnX3af3V0a8Q5J09pAcxTnsIa4MAAAgdBhH1kKG5yXp/741SC+tKw2aFUySjDGy2ZgZDAAARA6CrMUkxUbr4pPymi1/bEWJ6r1+XTE6XzHR3GgHAABdH0G2C1i/r0ZPr94rv5FW7Hbp59N6ql9mQqjLAgAA6FDcuusCil2NgdnBdlU26KcvbtCClXvk8zNMFwAA6LoIsl3AKf0y9Nezi9QnI06S5DNNXQ1ufHmjdlc1hLg6AACAjkGQ7SIK0+L053kD9Z0RuTp4c1Yby2r14+fX6eV1pUyiAAAAuhyCbBfisEfpstH5+tPcAcpPjpEkNfqM/vfjXfrVa5tUWe8JcYUAAAAnDkG2CyrKTtQD5xTpzKKswLJat0+JMTzbBwAAug6CbBcV67DrRxML9Ps5/ZSfHKMbp/RUdBTjzAIAgK6DINvFjcxP1t/PHaxeGfFBy7eU1+nFtfSdBQAA1sXfmiOA/Yg7sR6fX/e9v13bDtTrw+0V+unkQuWlxIaoOgAAgPbhjmwEWllcrW0H6iVJa/fW6IfPr9PCtfsYdxYAAFhKu+/IejweeTwd8RS8TcYYGf7k3WHG9EjWPaf30/wPdmhvtVuNPqN/fLpbH26r0PWTC9X9iLuzxhjJGHm9XknqoHZHODvU5rR95KHtIxdtH7lC3fZtOa/NtDExulwupaSkaMGCBYqPjz/+Dm2QmJioMWPG6qudZWpocJ/QY6O5Bq9fT609oDe3VAWWOaJs+lZRms4ckBp4OCw21qk+BVlavnyZampqQlUuAACIAHV1dbroootUVVWl5OTkY27b7iC7Z88eZWRkfKNCj1KSdhQfkNvj64BjoyVr9lRr/gc7tae6MbCsZ1qsbpzaU30y4uV02FWYny6v16NFixZp1qxZcjgcIawYnc3j8Wjx4sW0fQSi7SMXbR+5Qt32LpdLmZmZrQqy7e5a4HA4OuTiPF6fbDabbDaGiuosw/KS9b/fKtLjn+/RwrX75DfSzsqGQDvYbDbJZlN0dNOnS0e1PcIfbR+5aPvIRdtHrlC1fVvOyagFkNQ07uyV47prap80zf9gh0Z3T1GfjBPbdQQAAOBEIsgiSL/MBN1/VlGz8WU9Pr/uf2uzCujzDwAAwgRBFs00PeQV3LXjnx9s098/2KaEaLscBSU6b0wB3T8AAEBIMY4sjsvt9euJpTslSbVem37x/Fpd9u9l2rm/LsSVAQCASEaQxXE5o6P0wo8m6IwhuYFlH2wu16w/v6e/vr1ZjV5GmAAAAJ2PIItWyUyM0fwLhumqgT7lJsdIkhq9ft27eJPm3P+BPt5SHuIKAQBApCHIok2GpBm9/pOJunJSL9kPTpiwtaxWF/1zqX72zGpmZAMAAJ2GIIs2S4yJ1q1nDtLLP56kkwpSA8tzk2N5AAwAAHQagizabVBesp69ZoLu/tZQDeueoh/P6BvqkgAAQARh+C18I1FRNn1nbIEuHNOj2d3Yv769Wftr3bphdn8lxzIrDAAAOLEIsmiTxMTEFpcfGWK3ldfqr+9ukdvr16tr9uiW04t01og8uh4AAGABMTExoS6hVQiyOC57lE02m2SMTWPGjJVkk+c4Q26t2V2pg8+Cqay6Udc/vUqPf7pDvz6zSIO6JXd80REiymaT3U4PIQDW5fP5m80miVCzaezYcaEuolUIsjiuqKgoGb+0p9Sl/RWVSkpMPO6d1cFpTv37gqF64MMd+nh7pSTpsx0VOud/P9aZg7L13bHdlRzLp9834XDY1S07RfZQFwIA34DfGO0prZLHw5jk4cIRbVdqkjW6BJIk0Gpuj1cNDW7FxPha1UUgLTZavz6lj5bvqtLfP9mlYlej/EZ66ctSvbtlvy4fna/TBmQGhvECAEQmj8cnN0E2bDQNpUmQBSRJY3qkaHhekl5YW6onV+5Rg9ev6kafHvhop2w26fSBWaEuEQAAWBCd69ApnPYonT88Vw+dN1jT+qRJkvKSY3RKv4wQVwYAAKyKO7LoVJkJTv1yem+dPrBaRk0B93Bf7KlWUXaCHDzABAAAjoMgi5AY2i2p2bKdFfW6+bVNyk2K0ZXjumtcQQrDdQEAgKPithfCxj+W7pbPSMWuRv128Vf61eubtXV/XajLAgAAYYogi7Bxxeh8Dcn9esKF1SXV+vHC9br/gx2qqPOEsDIAABCOCLIIG30z4/WHM/rr1pm9lZvklCQZSW9sLNf3nlmr/67eK7fXH9oiAQBA2CDIIqzYbDZN7JWmv587WN8bm694R9OnaL3Hr/8sL9bVz36pEldjiKsEAADhgCCLsOS0R+ncYbn61/lDdPrAzMB0t3GOKOUkOkNbHAAACAuMWoCwlhrn0LWTCjV3ULYeWrpb3x6a02wmsKoGr1KY7hYAgIjDHVlYQs/0ON01p59O6p4ctHxjWa0uffIL/ePTXXI1eENUHQAACAWCLCzLGKN/L9stj89o4dpS/c9/mx4Ia+SBMAAAIgJBFpblM1JRdqKc9qauBrVun/6zvFhXPrNWizaVy+c3Ia4QAAB0JIIsLCs6yqYrxuTrn+cP0ez+GYEHwsprPfrz+zv044XrtGxXlYwh0AIA0BURZGF5WQlO/XRKT/3tnEEa2yMlsHx7RYNue3OLbnptk6ob6T8LAEBXQ5BFl9EzPU6/PbWv/t8Z/dU/Kz6w3O0zSnTaQ1gZAADoCARZdDnDuiVp/ryB+tWMXuqWHKPvjc2XzRY8ZBdT3gIAYH0MvokuyWazaUrvdE3smdZs3NmVxS795s0tOqMoSxcMz1VavCNEVQIAgG+CO7Lo0o4MscYYPfJZibx+oxe/LNV3/7tWDy8vpg8tAAAWRJBFRPH6jUbkJSkmuulTv9Hr19Or9+p/nl6rp1ftUYPHF+IKAQBAaxFkEVEc9ihdMSZf/z5/iOYNylL0wTu2NW6fHv6sRN99eq1eXFsqt49JFQAACHcEWUSk9HiHfjChQP88b3DQGLSVDV49+OkuXfnftapmylsAAMIaQRYRLScpRj+d0lMPfnuwpvRKCywvSItTUizPQgIAEM4IsoCkHqmx+tXM3vrr2UUa0yNZF4/sFrTeGKP3tx6Ql2lvAQAIG9xyAg7TNzNed5zar9nyT3ZU6e53tiknsVgXjuimU/pnBPrXAgCA0OCOLHAcxhg9sbJEkrSvxq37P9yhK/+7Vm9uLOcOLQAAIUSQBY7DZrPpxxMKNKp7cmDZvhq35n9AoAUAIJQIskArFOUk6s7T+um+uQOOGmhfXV8mt5dhuwAA6CwEWaANDgXae+cO0En5wYH2gY92qsTVGMLqAACILARZoB0G5STqrjlNgfbQHdrxhSnqmR4XtJ0xdDkAAKCjMGoB8A0MOniHdlNZrWKjg38v9PmNfvrSBp2Un6yzh2QrNc4RoioBAOiaCLLACdA/K6HZsg+2VWhzeZ02l9fphbX7dNrALH17WI6yEpwhqBAAgK6HIAt0kL3VjYqOssnrN2r0Gb34ZaleXV+mU/pl6LzhucpLjgl1iQAAWBpBFuggF47oppl9M/Tcmn16Y0OZGn1GXr/RGxvLtWhTuSb1TNO5w3PUL7P53VwAAHB8POwFdKCsRKeuGd9DD184VBcMz1W8o+lLzm+k97dV6CcvbNDza/aFuEoAAKyJIAt0gtQ4h64Yk69HLhyqy0blKTW26Y8hNknjClJCWxwAABZF1wKgEyXGROs7I7vp20Nz9NaW/SqubFB+SmzQNm9t3q/qBq9OHZCpeKc9RJUCABD+CLJACDijo3T6wKxmy71+o8dWlKi0xq0FK/fojKIszRucrfR4hu4CAOBIBFkgjKzdW63SGrckqcbt09Or9+q5Nft0Sr8MnT0kW4Vpccc5AgAAkYM+skAYGZGXrAe/PUiz+mcoOsomSYGRDq55bp1ueX2zPttVxYxhAACIO7JA2ClMi9MNU3rqslF5evHLUr22vkx1Hr8k6fNilz4vdql/Vrz+9u3BIa4UAIDQ4o4sEKYyE5z63tjuevQ7w3T1yd2Vm/T1jGADshIUZbOFsDoAAEKPIAuEuQSnXecMydE/zxuiW0/prSG5iTprcHbQNlX1Hv3smdX6YndlaIoEACAE2t21wOPxyOPxnMhaDrLJGEMfwDBiZAL/lUTbhEiUTZpQmKoJhamSDraDMfJ4PFrw6XY9u2K3nl2xW6MLU3XF+ELNHJilaPuJ+V310Nd6x3zNI5zR9pGr89qen/vh5lBTeL3ekJy/LZ9z7Q6y7777ruLj49u7e4sSExM1ZsxYVdfUqKHBfUKPjW/AJMiYZNXV1UqSXC5XiAuCJMXGOlVfH6fly5fpiY/r1TS9gvTZjkp9tqNSqU6jiTl+jc8xSjpBo3ctXrz4xBwIlkPbR66ObHt+7oen2FincjLi9PHHH6umpqbTz19XV9fqbW2mjb8CuVwupaSkaM+ePcrIyGhzca0oSTuKD8jt8XXAsdEeCfFO5WWnaHvxAZWXH1BycrJs9M8MOafDrsL8dElG9W6fXvpijx7+eIe2lNUGbeew23TGkFxdcnKBhndv3yxiHo9Hixcv1qxZs+RwMKZtJKHtI1fntT0/98ONI9qunIw4ORzRio7u/HEBXC6XMjMzVVVVpeTk5GNu2+7qHA5Hh3xie7w+2Ww2glIYsckW+K8k2idM2Gw2yWaTIzpaDodDl4zvpYtP7qkPNpfr0U+26+0NpTJG8viMXli9Ry+s3qP7zh+ub53Uvd3n7Kive4Q/2j5ydXTb83M//BxqiuiDP186W1vOyfBbQBdis9k0pX+WpvTP0q4DdXr80x16+rNdqqzzKM5h18yinKDtjTH88AAAWBZBFuiieqTH61enF+mns/rrpdUlOlDrVkpc8G+5t7/0pYorG3T5hEJN6ptJqAUAWApBFujiYh12nT+6R7PlrgaPnlmxW3Vun95av0+9MxP0nbEF+vao7kpPcLZwJAAAwgvjyAIRaltZrZJjv75Du7W8Vne9tl4n//5tXffUSn26dT/D4QAAwhp3ZIEINbxHqj745XQtXrdPj36yXZ9uPSBJcvv8enFViV5cVaI+WU13aS8cnR/iagEAaI47skAEc9ijdPrQbnrq6vF6+8apumpyL6XFf32X9quyWj36yQ45oug7CwAIPwRZAJKkPlmJuuWMQfrkVzN1/4UjNLZXuiTpwrE9FHVEkH1tzR5V1THTEwAgtOhaACBIrMOus0bk66wR+dpSWqOMIx782lpWox8+8blioqN0xtBuOm90D43rld4s7AIA0NEIsgCOqm92oqTgea+fXLZTktTo9ev5lcV6fmWxCtLjdd6o7jp3dHd1S4kLSa0AgMhDkAXQJheM6SGPz+j5z3fL1eCVJO08UKd7F2/SfW9t0pR+WTp/dA+dMihbMdH2EFcLAOjK6CMLoE36Zifp9nmDteyWU/SX74zU5H6ZgekMjZHe21SmHy34XDc8vTq0hQIAujzuyAJol1iHXfOG52ne8DztrqjTcyuK9cyKXdpdUS9JOmNYt6DtvT6/aht9Sonv/Hm7AQBdE0EWwDfWPS1e153ST9fO6KtPt+7XS6tLNLMoO2ib9zaV6QdPfK5TB+fqWyfla3LfTEXb+aMQAKD9CLIATpioKJsm9M3UhL6Zzdb997Ndcnv9enl1iV5eXaLMxBidNSJP54zM1+C8ZNlsjHoAAGgbgiyADmeMUc/MBKXFO1RxcPzZ8ppG/evDbfrXh9s0ICdJ55yUr7NH5Cs3JTbE1QIArIK/6wHocDabTb+aU6SlN5+if1w6SnOG5Mp5WLeCjfuqdc/rGzT+nrf1xNIdIawUAGAl3JEF0Gmc0VGaPThXswfnqqrOo1fWlOj5z4u1YkeFpKZRD0YVpgXt0+DxyWGPkp0JFwAARyDIAgiJlHiHLh5XqIvHFWp7ea0WrizW+j0uDcxNDtru4Y+36z8fbdOZw/J01og8Dc1PoT8tAEASQRZAGOiZmaCfzurfbLkxTRMv7HN93Z+2Z0Z807BfI/LUNzspBNUCAMIFQRZA2Kpz+1SYkaBt5bXy+Iwkafv+Ov3lnS36yztbVNQtWfOG52nu8G7qnhYf4moBAJ2NIAsgbCXEROuhy0arss6tN9bu1UurS/TJ1v0yTZlW6/e4tH6PS//vjQ361+WjNbMoJ7QFAwA6FUEWQNhLjXfqwrEFunBsgfa5GvTKF3v00uoSrd5VKUmKiY7S2F7pQfvsrWpQnNOulDhmEgOAroogC8BScpJj9b1JvfS9Sb20Y3+tXl5douoGr5JigwPrH9/cqJdWF2tyvyydPrSbZhXlMD0uAHQxBFkAllWYkaAfz+jXbHmDx6dFX+6Vx2f0zoZSvbOhVA67TRP7Zur0od00e1COUuOdIagYAHAiEWQBdDn1bp/OHd1dr6/Zq72uBkmSx2e0ZGOZlmws080Hp9I9Y2iuThvSje4HAGBRBFkAXU5aglO3zR2sX58xSCt3Vei1NXv1+po9KqlqCrVev9H7m8r0/qYyFXVL1rDuqaEtGADQLgRZAF1WVJRNowrTNaowXbecXqRVuyv12hd79PravSqurFf3tDgNzU8J2ueNtXtUXuPW7EE5yk6ODVHlAIDWIMgCiAhRUTadVJCmkwrSdMsZRVq9u0oHahubzRL24HtbtWpXpX794lqN7JGq2YNzdergXPXKTAhR5QCAoyHIAog4NptNI3qkNlteUlmvVQeH9DJG+nxnpT7fWal7Xt+g/jmJOnVwrmYPytWQ/GSmyQWAMECQBYCDuqXE6pVrJ+nNL/fqzS/3atO+msC6TftqtGnfFv31nS3KS4nV/10ySsNbCMMAgM5DkAWAg2w2m4bkp2hIfopunD1A28prtehgqF25qzIwo9i+6kYVZgRPiVtW3ajEmGjFOe0hqBwAIhNBFgCOoldmgr4/tY++P7WPSl0NWrx+nxZ9uU82m5qNQ/unNzfqhVXFmtg3UzOLsjVjYLa6pcSFqHIAiAwEWQBohezkWF08rlAXjyuU32+C1vn8Rm+t36dGrz8wAYMkDc5L1syB2ZpZlKOh+SmKiqJfLQCcSARZAGijIwNpTaNXswfn6O31pSqtbgws/7LEpS9LXPrLO1uUlRSjGQOy9b3JvdQ/J6mzSwaALokgCwDfUEqcQ3d/a5j8fqMvS1x6e8M+vb2+VGuKqwLblFU36unPdunCsT2C9vX7DXdqAaCdCLIAcIJERdk0tHuKhnZP0fWn9Nc+V4Pe2VCqt9fv04dbypXgjNbwI2YR+/dH2/Tsit2aPjBbU/tnaVRhmhz2qNBcAABYDEEWADpITnKsvjO2QN8ZW6AGj0/bymub3X19e32pNuyt1oa91fq/JV8pMSZaE/tmaNqApmCbl8oDYwBwNARZAOgEsQ67irolBy3z+Y3cPn/QsppGr978cp/e/HKfJKl/TqKmDcjWOSPzm+0PAJGOIAsAIWKPsum5H0xQeU2jPthcpiUby/T+pjJV1HkC2zRNxFCjftmJQUHWGMPsYgAiHkEWAEIsMzFG54zsrnNGdpfPb7SmuEpLNpbqvU1lWnVwIoap/bOC9lmysUy/e2Wdpg7I0pR+WRrbK10JMXxLBxBZ+K4HAGHEHmXTiB6pGtEjVdef0l8VtW6t3FWh7OTYoO3e21SmreW12lpeq/98tF0Ou00jC9I0uW+mJvXL1ND8FEXz0BiALo4gCwBhLC3BqRkDc5otd9V7ZI+yyXdwcgaPz2jZtgNatu2A7l28SUmx0ZrQJ0Nnj8jXnKHdOrtsAOgU/LoOABZ03wUjtPI3s/TgJSfpkpML1DMjPmh9dUPTQ2Ordlc227eyzt1JVQJAx+KOLABYVHKsQ6cN6abThjTdcd11oE4fbSnXB1vK9fGWclXUeTSpb2bQPjv312nqn97V4LxkTeybqYl9MjW6Z5rinfw4AGA9fOcCgC6iR3q8LhxboAvHFgRmGeuXkxi0zYdbymWMtLbYpbXFLv39va1y2G0a3j1V4/tkaHzvDJ1UmKZYhz1EVwEArUeQBYAu6NAsY0eKd9o1qFuy1u1xBZZ5fEaf7ajQZzsq9Nd3tshpj9KU/ln65+WjO7NkAGgzgiwARJCzR+br7JH52l/TqI++2q9PvtqvT7fu17by2sA2bp9ffmOa7fviqmLlpcZpePdUOaN5xAJA6BFkASACZSTGaN7wPM0bnidJ2lvVoE+3NgXbT7bu1/jeGUHbe3x+/er5Napz+xTnsGt0zzSd3DtDJ/dO19B8gi2A0CDIAgCUmxIbuFsrKTCs1yFri12qc/skSfUenz7YXK4PNpdLkmIdURrRI1Vje6ZrbK8Mje5JH1sAnYNfoQEAzdijgqe/LcyI1//79lCdMzJfuUdMztDg8evTrQf0l3e26JJ/LdX+2uDhvfz+5t0UAOBE4I4sAOC40hOcumBMgS4YUyBjjLbvr9MnX+3X8u1NkzAUV9ZLkvJT45SfGhe0792vr9eSjWUa0ytd43qla0zPdOUdsQ0AtAdBFgDQJjabTb0yE9QrM0EXjSuQJJVU1mv59gNq9Pqbbb902wFtLq3R5tIaLVi6U1JT4B3bK12je6ZpVGGa+mUnNbsLDADHQ5AFAHxjealxOmtEfrPlPr9RdJQtaDpdSSqurNfClcVauLJYkpQUE61fnzlI54/p0Wk1A7A+giwAoMPYo2x6/ocTVdvo1cqdlVq2bb+WbT+glTsrg+7eVjd6lZnkDNp35/46/d97X2lUYZpOKkhVr8wE2WzctQXwNYIsAKDDJcREa1K/TE3q1zRlrtvr15riKq3cWaEVOyr0+c4KjeyRFrTP0m379eSynXpyWVN3hPQEp04qSNVJhWkaVZCmYd1TFedkdAQgkhFkAQCdzhkdpVGFTf1jr5zc8jaf76wIen2g1q231pfqrfWlkqToKJuKuiVrav8s/ezUAR1dMoAw1O4g6/F45PF4TmQtB9lkjJFpYVYZhIaRCfxXEm0TJowxkjEd9HUY7NA5OuNcCC+hbPubTu2nM4bk6POdlfp8Z6VW7apSZf3XdXj9RmuKq5Se4GhW39Of7VZGglPDu6coKymms0vvEjqv7fm5H24ONYXX6w3J+dvyOdfuIPvuu+8qPj6+vbu3KDExUWPGjFV1TY0aGtzH3wGdwyTImGTV1TVNYelyuY6zAzpDbKxT9fVxWr58mWpqajrlnIsXL+6U8yD8hLLtCyUVZkhnpUtlDdK2apu2V9u0tdqmffU2xdft02uvvRbY3m+kO5bZ5fY39adNdRoVJn790SNRiqFHQqt1ZNvzcz88xcY6lZMRp48//rjTfr4crq6urtXbtjvITp8+XRkZGcffsM1sSkpMVEyMrwOOjfZIiHfKZrMpPj5BdXWNSk5O5oGLMOB02BUXF6cpU6Z0+Lk8Ho8WL16sWbNmyeFwdPj5ED7Cve2rGzzy+aXU+K9r27i3Wu5PPwm8rnTbVHnAptUHml5H2aR+2Yka1j1F10zppYL0E3tTpqvovLbn5364cUQ3/aY3YcIERUd3fi/Uttwwa3d1DoejQz6xPV6fbDYbQSmM2GQL/FcS7RMmbDabZLPJ0YnfZDrq6x7hL1zbPr2FmgqzkvTgJSdp1a4qrdpVoTW7q1Tr/jok+Y20cV+NNu6r0XWn9A+6rlW7KrVxr0tD8lPUPydJDjsTYHZ02/NzP/wcaoro6OiQfN235Zw87AUA6FKSYh06bUg3nTakm6SmsWy3lNZo9a5KrdpdqVU7K7VxX7XS4p3NZiF7aVWJ/v3RNklND6QV5SZpSH6KhuanBMKtM5pwC4QLgiwAoEuzR9k0IDdJA3KTAhMu1Lt92l1R1+wu4NriqsC/3V6/Vu+u0urdXy9z2qM0sFuSzhmZr+9O7NU5FwDgqAiyAICIE+e0q19OUrPlvzhtgFburNSa4iqtLa7S1vLaoPVun19f7K7ShD6ZQcuNMfrty+vUNztRRd2SNTA3SQkx/IgFOhpfZQAAHDS6Z7pG90wPvK5u8OjLEpfWFldpzcGPrWW1GpqfErTfjv11evjj7YHXNpvUKyNBRXnJGtTt4EdesrKTYugLCpxABFkAAI4iKdahk3tn6OTeX4/SU93gUXRUcD/ZNYd1SZCaxuHcWl6rreW1evWLPYHlGQlOvXH9FMa2BU4QgiwAAG2QFNv8ieqZRdl69prxWrfHpXUlLq3b49KGvdVye/1B27m9fmUmOoOW3bdoo97dWKaibkkH79ymaEBuklLiwm+UCCDcEGQBAPiG4p3RzboleH1+bSuvDQq3Cc7oZl0LVu6qDHRbOFxeSqz6H3xIbWBukkb2SFPPzIROuR7AKgiyAAB0gGh7lPrlJKlfTpLOGpF/1O0aPX7ZbF9PC3pISVWDSqoatGRjmSTp+1N761dzigLr/X6jxev3aUBOkgrS4xUVRd9bRB6CLAAAIfTfa8arzu3Vxr3Vgbu3m/ZVa8PealU3fD3X/cDc4FEWdlfU6/uPrZAkxTns6peTqAE5SYGhxgbkJikrkYfL0LURZAEACLF4Z7RGFqRpZEFaYJkxRnuqGrRxX7U27q3W6ML0oH027qsO/Lve49MXu6v0xe7g7gnpCU71y07UQ5ePVnILfXsBqyPIAgAQhmw2m/JS45SXGqfpA7Kbre+ZEa+fzOirDXurtWlftXYcqGvWPeFArVtri6uUdMSYtvPf2qT3N5Wpb3ai+mUnqW9OovpmJSo/NY4uCrAUgiwAABbULydJN8weEHhd5/Zq876awB3cQ90TuqXENutesHpXpT7f2fRxuDiHXX2yE9QvO0m9M+LkqxYQ1giyAAB0AfHOaA3vkarhPVKDljd4fM22rWn0NlsmNXVRWFvs0tpilyRpRrfg8XL9fqP7396sXpkJ6pWZoN5ZCS0ORwZ0FoIsAABdWKzD3mzZM9dMUG2jV1+V1WjzvhptLq3RltIabSmt1s4DdfIf7KKQEx/cV6Gkql73v705aFlWUox6ZSaoT1aCemcmBgJuYUaC7HRTQAcjyAIAEIESYqI1rHuqhnVPDVre4PFpW3mtNpRUqvKrlUHrtpbVNjtOWXWjyqobtWzbgaDlH/xiunqkxwdef1VWo/LqRvXOSlRmopPRFHBCEGQBAEBArMOuom7J6psZp9d2BwfZofkpeuiy0dpWXqOtZbVNH+W1Kq9pDNouJjpK+alxQcv++9ku/f29rZKkpNho9c5MUO+sRBWkxynRZpSV4FBeUoySYokmaD0+WwAAQKukJTg1a1COpJyg5VX1Hm0rrw0E3Hq3r9noB4ffza1u8Gr17iqtPmK4MEma0itNv5rZO2jZV/vrlBHvUEps85nRENkIsgAA4BtJiXNoRI9UjTjiQbPDnTmsm/JSYrW1vOlObklVfbPhwiQpPSH44TFjjH728kY1eP2Kc0QpLzlGecmxykuOUbfkmMD/0+MdiiLkRhyCLAAA6HBnjcgPmqr3UF/crWU1Wr21XLsq61XialS/jPig/SrqvWrw+iVJ9R6/vtpfr6/21zc7fozdpt+d1k9Du309A1p1g1eVDV7lJDrljI5qtg+sjyALAAA6XaAvblaCBqY45G5hmDCp6Y7smUVZKnE1qsTVoNIad2BUhcM1+ozS44Pv5n66s1L3vb9DkpQe71BuklO5STHKSYoJ/Ds3KUYZ8Q5GWLAogiwAAAhbGQlO/WhiQeC1129UWt14MNg2fexxNWpvdaOyE51B++6tdgf+faDOowN1Hq3b13zkhe4pMXrovCFBy9bvq1GUzabc5Bglx9jpmxumCLIAAMAyoqNsykuJVV5K7HG37Zkep2l90rS32q191Y2qqG95IoisBGezZQ9+ukubyuokSXGOKOUkxig70amcJKeyE53KSnAqO8mp7imxzaYARufhnQcAAF3S5F5pmtwrLfC6wePTvhq39lY3BsLt3upG9ctMaLbv4Xdz6z1+ba+o1/aK5n1zrx7XXecM/XoUh5pGr15YW6qsxKbAm5PoVGaiU047fXQ7AkEWAABEhFiHXYVpcSpMizvmdsYYfXtozsGg2xR8S2vc8rbQOffI7gzFrkY9sXJPs+3S4qKVnRijnERnIOTOHpCpWB5C+0YIsgAAAIex2Ww6f3hu0DKf36iy3qPSGrf21bhVVuNWaY1bPdODQ3FpjVstqaj3qqLeq42Hjad72oDMoG1eXleqT3ZUKjPBefDDocyEpm4MmQkOJTrpq3skgiwAAMBx2KNsykhwKiPBqaKco283NDdRt83qo9LDwu6+GrdKa4L76KbFRTcbEmxzeZ1WFlcf9dgx0VHKTHBofGGqvje2e9C63VUNSo6JVlKEPZhGkAUAADhBUuMcOrkwtcV1bq9f5bVNwbbx4Ni4h6s6ysNohzR6/SquapSrofl2N768Ua4Gr2LstmZ3dDMSHMqIdyoj3qEeqbGKd9rbdW3hiCALAADQCZzRUcccceH22X1U4/apvNaj8lq3ymrdgX+X17pVXuNRWa272SgLbq8/EG4bfUbFrkYVuxpbPMdts/oEBe2dFfV68ctSpcc3hd70eIdyk2OUmOhU3LG7EocFgiwAAEAYsNlsSoqJVlJMtHqlt5wijTHNJoRw+/w6tX+GyusOBuAat+o8ze/4Smo2acSOyga9tqG82XZ/OW+w5gxrPppDuCHIAgAAWITNZpP9iC6wiTHRun5Kz6BltW6f9te6VVbr0YG6pju7B+o8ykmKCdpuf23LD6cdORpDuCLIAgAAdDEJTrsSnHEqOM5QY6f0y9CgnEQdqPMcDLtuVTR41a0VE06EA4IsAABAhEqMiVb/rOA46IiOUtoRXRDCFaPwAgAAwJIIsgAAALAkgiwAAAAsiSALAAAASyLIAgAAwJIIsgAAALAkgiwAAAAsiSALAAAASyLIAgAAwJIIsgAAALAkgiwAAAAsiSALAAAASyLIAgAAwJIIsgAAALAkgiwAAAAsiSALAAAASyLIAgAAwJIIsgAAALAkgiwAAAAsiSALAAAASyLIAgAAwJIIsgAAALAkgiwAAAAsiSALAAAAS4pu6w7GGEnSgQMHTngxTWxqqK+V2+ProOOjraLklsslNTTUyettVGN9nWQLdVXwe+1yVdklmQ4/l9frlc1m04EDBxQd3eZvG7Aw2j5ydV7b83M/3Pii7XK5PKqvrwvJ1311dbWkrzPnsbS5ukMH79+/f1t3BQAAAFqlurpaKSkpx9zGZloTdw/j9/tVUlKipKQk2WzclgMAAMCJY4xRdXW18vLyFBV17F6wbQ6yAAAAQDjgYS8AAABYEkEWAAAAltSmIPvLX/5SkydP1qWXXiqPx9NRNaGTLVu2TOPHj9eUKVP0ne98J6htfT6f/ud//keTJ0/W9ddfH1h+//33a+LEiZo3b55cLpck6cMPP9SECRM0adIkrVmzprMvA+305JNPKisrK2gZ7R4ZlixZopkzZ2r69OlauHBh0LqWvt8/88wzmjBhgmbOnKndu3dLkjZs2KApU6ZowoQJevvttzv9GtB2fr9fV1xxhSZPnqxJkyZpw4YNQetp+66nqqpKY8eOVWJiotauXSup5TY9XGu/3+/du1ezZ8/WxIkT9fjjj3feRR1iWmnVqlXm4osvNsYYc+edd5oFCxa0dleEuZKSElNXV2eMMeamm24yzzzzTGDdCy+8YG655RZjjDFXXnml+fjjj01ZWZmZMWOG8fv95vHHHze///3vjTHGTJkyxRw4cMDs2LHDzJkzp/MvBG3m9XrNOeecY0aOHBm0nHbv+urq6syZZ55pGhsbm61r6fu9x+MxJ598smlsbDQffvihufrqq40xxpxzzjlm06ZNpqqqykyYMKFTrwHts2LFCnPhhRcaY4x5//33zVVXXRVYR9t3TW6325SWlprLL7/crFmz5qhtekhbvt9ff/31ZvHixYFj1tfXd+q1tfqO7Mcff6zZs2dLkk477TR99NFHHRau0bm6deumuLg4SZLT6Qx6QrCldl++fLmmTp0qm80WWFZfXy+73a60tDQVFBR04DjDOJGefPJJnXfeec2eCqXdu75PPvlEcXFxmjt3rs455xzt3bs3sK6l9t+8ebOKiorkdDo1ceJEffHFF5KkkpIS9evXT8nJyUpPT1d5eXlIrget1717dxljZIxRRUWFMjMzA+to+67J4XAE/eXtaG16SFu+3y9btkwzZsxQdHS0Ro8eHbjj21laHWQrKiqUnJwsSUpJSeEHVhe0Y8cOLVq0SHPnzg0sa6ndj7dMkqKjo+V2uzv3AtAmPp9P//3vf3XBBRc0W0e7d3379u3Tli1b9PLLL+uqq67S7bffHljXmrb2+ZoGr/f7/YFl/GywhszMTDkcDg0cOFDXXnutfvjDHwbW0faR4Wht2tL6432/93g8gZshofg8aHWQTU1NDfSRqKqqUnp6eocVhc7ncrl06aWX6uGHH5bD4Qgsb6ndj7dMapoRxul0du5FoE0ef/xxnX/++S2O0Ue7d32pqamaOHGinE6nZs6cqS+//DJo3fHa2m63S1LQ5w8/G6xh0aJFio6O1saNG/Xcc8/pxhtvDKyj7SPD0dq0pfXH+37vcDgCv9SE4vOg1UF2woQJeuuttyRJb775piZOnNhhRaFzeb1eXXjhhbrttts0YMAASVJxcbGMMS22+5gxY/T+++8HLYuPj5fX61VlZaV27drFNzQLWLdunR599FGddtpp2rx5s37yk5/Q7hFkzJgxWr9+vYwxWrVqlXr37q3S0lK53e4W279fv35av3693G63Pv74Yw0bNkxSU9ekr776StXV1Tpw4EDQn6kRnowxysjIkNR0d7aqqoq2jzBHa9MDBw6orq6uTd/vx4wZoyVLlsjr9WrFihUaPHhw515MWzrU/uxnPzOTJk0yF110UYsPCMCaHn30UZOenm6mTp1qpk6dap566ikzbdo009DQYDwej7n88svNpEmTzLXXXhvY57777jMTJkwwZ5xxhqmsrDTGGPPee++Z8ePHmwkTJphVq1aF6nLQDqNGjTLGGNo9wjzwwANm8uTJZsqUKWbLli3m4osvNlu2bDHGtPz9/qmnnjLjx48306dPNzt37jTGGPPll1+aSZMmmfHjx5tFixaF7FrQeh6Px5x//vlmypQpZty4ceajjz6i7SPAnDlzTLdu3czJJ59s/vOf/7TYprfccot56623jDGt/35fUlJiTjnlFDN+/HjzyCOPdPp1MbMXmvH5fPrRj36kBx98MNSloBPR7rjqqqv00EMPhboMhABtD0n6wQ9+oL/+9a+Kjo4OdSmtRpAFAACAJTGzFwAAACyJIAsAAABLIsgCAADAkgiyAAAAsCSCLAAAACyJIAsAAABLIsgCAADAkgiyAAAAsCSCLAAAACyJIAsAAABL+v8qhRfixxIGTwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import expon\n",
    "\n",
    "# Bind reference to best fitting exponential distribution.\n",
    "scale = 4353\n",
    "dist = expon(scale=scale)\n",
    "\n",
    "xx = np.linspace(0, max(x), 1000)\n",
    "yy = dist.pdf(xx)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(7, 5), tight_layout=True) \n",
    "\n",
    "ax.set_title(\"data with exponential overlay\", fontsize=9, weight=\"normal\")\n",
    "ax.hist(\n",
    "    x, 6, density=True, alpha=.80, color=\"#ccd3e1\", \n",
    "    edgecolor=\"#FFFFFF\", linewidth=1.0\n",
    "    )\n",
    "ax.plot(xx, yy, linewidth=2.0, linestyle=\"--\")\n",
    "ax.tick_params(axis=\"x\", which=\"major\", direction='in', labelsize=6)\n",
    "ax.tick_params(axis=\"y\", which=\"major\", direction='in', labelsize=6)\n",
    "ax.get_xaxis().set_major_formatter(mpl.ticker.FuncFormatter(lambda x, p: format(int(x), ',')))\n",
    "ax.set_yticklabels([])\n",
    "ax.set_xlim(0)\n",
    "ax.xaxis.set_ticks_position(\"none\")\n",
    "ax.yaxis.set_ticks_position(\"none\")\n",
    "ax.annotate(\n",
    "    r\"$\\hat \\theta_{MLE} = $\" + f\"{scale:,.0f}\", xy=(.80, .90), \n",
    "    xycoords=\"axes fraction\", ha=\"left\", va=\"bottom\", fontsize=11, \n",
    "    rotation=0, weight=\"normal\", color=\"#000000\"\n",
    "    ) \n",
    "ax.grid(True)   \n",
    "ax.set_axisbelow(True) \n",
    "plt.show()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Maximum Likelihood Estimation for Censored Data\n",
    "\n",
    "Right-censoring occurs when the event of interest has not occurred for some of the observations by the end of the study period. Their data is censored at the last known point. One example is for loss data, when a particular loss is at the limit: If we do not incorporate censoring within our likelihood estimator, the distributional assumption will not accurately reflect the underlying process, specifically by underestimating the tail of the distribution.\n",
    " \n",
    "\n",
    "For censored observations, instead of including the density we use the survivial function. Recall that the CDF for the exponential distribution is given by:\n",
    "\n",
    "\n",
    "$$\n",
    "F(x) = 1 - e^{x/\\theta}, \\hspace{.50em} x \\geq 0.\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "The survivial function is defined as $S(x) = 1 - F(x)$. To incorporate right censoring, the likelihood function is modified in the following way: \n",
    "\n",
    "$$\n",
    "L(\\mathbf{\\theta}|\\mathbf{x}) = \\prod_{i=1}^{n} f_{i}(x_{i}|\\mathbf{\\theta})^{\\delta_{i}} \\cdot [1 - F(x_{i}|{\\theta})]^{1 - \\delta_{i}},\\\\\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\delta_{i} = \n",
    "\\begin{dcases}\n",
    "1 & \\mathrm{if} \\hspace{.25em} x_{i} \\hspace{.25em} \\text{is not censored} \\\\\n",
    "0 & \\mathrm{if} \\hspace{.25em} x_{i} \\hspace{.25em} \\text{is censored} \\\\\n",
    "\\end{dcases}\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's looks at a second example:\n",
    "\n",
    "\n",
    "> uncensored =  [361., 433., 716., 753., 772., 2007., 3563.] \n",
    "  \n",
    "> censored = [10000., 10000.]\n",
    "\n",
    "\n",
    "We will fit two distributions: The first not considering right-censored data, the second considering losses at the limit as censored. Adpating our code from the first example, we obtain the maximum likelihood estimate as follows:\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "+ Uncensored MLE scale parameter estimate: 3,178.\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Using Scipy optimization routines to find maximum likelihood estimates.\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "np.set_printoptions(suppress=True, precision=8)\n",
    "\n",
    "\n",
    "# uncensored = [361., 433., 716., 753., 772., 2007., 3563.] \n",
    "# censored = [10000., 10000., 15000.,]\n",
    "\n",
    "uncensored = [361., 433., 716., 753., 772., 2007., 3563.] \n",
    "censored = [10000., 10000.]\n",
    "vals = uncensored + censored\n",
    "\n",
    "\n",
    "def loglik(scale):\n",
    "    \"\"\"\n",
    "    Exponential log-likelihood.\n",
    "    \"\"\"\n",
    "    n, sum_x = len(vals), sum(vals)\n",
    "    return -(n * np.log(1 / scale) - sum_x / scale)\n",
    "\n",
    "\n",
    "# Maximize uncensored loglikelihood.\n",
    "results = minimize(loglik, 361, method=\"Nelder-Mead\")\n",
    "\n",
    "scale_u = results.x.item()\n",
    "\n",
    "print(f\"\\n+ Uncensored MLE scale parameter estimate: {scale_u:,.0f}.\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we implement a function which represents the censored loglikelihood for the exponential distribution. The censored likelihood is:\n",
    "\n",
    "$$\n",
    "L(\\mathbf{\\theta}|\\mathbf{x}) = \\Big(\\frac{1}{\\theta}\\Big)^{7} \\cdot e^{-8605/\\theta} \\cdot e^{-35000/\\theta}.\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "The loglikelihood is then:\n",
    "\n",
    "$$\n",
    "\\mathcal{l}(\\mathbf{\\theta}|\\mathbf{x}) = 7 \\cdot \\mathrm{Ln}\\Big(\\frac{1}{\\theta}\\Big) - 43605/\\theta\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "+ Censored MLE scale parameter estimate: 4,086.\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "def loglik2(scale):\n",
    "    \"\"\"\n",
    "    Exponential log-likelihood considering right-censored observations.\n",
    "    \"\"\"\n",
    "    n_censored = len(censored)\n",
    "    n, sum_x = len(x), sum(x)\n",
    "    return -((n - n_censored)  * np.log(1 / scale) - sum_x / scale)\n",
    "\n",
    "\n",
    "# Maximize uncensored loglikelihood.\n",
    "results2 = minimize(loglik2, 361, method=\"Nelder-Mead\")\n",
    "\n",
    "scale_c = results2.x.item()\n",
    "\n",
    "print(f\"\\n+ Censored MLE scale parameter estimate: {scale_c:,.0f}.\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Let's visually compare the uncensored and censored maximum likelihood estimates:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAAHqCAYAAADBKxzcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAACNsElEQVR4nOzdd3xUVf7G8c+dkt4rSQi9g/QuVcSCgmBBQUXXgg1XXVdd0RXXws/eZe2iLtiwIyoooDRBERAFpEgNgZDeJ1Pu74/IwJAASZgU4Hm/XpGZO/ee+52TIT7cnHuOYZqmiYiIiIhIA2Wp7wJERERERI5EgVVEREREGjQFVhERERFp0BRYRURERKRBU2AVERERkQZNgVVEREREGjQFVhERERFp0BRYRURERKRBU2AVERERkQZNgVVE6s3q1asxDKO+yzhhREVFsXDhwvouo9acffbZTJs2rb7LEJF6oMAqInISmD59OoZhcNlll/ls37NnD3a7naioKO+2K6+8kltvvbXSdq688koCAgIICwvz+crIyPBrvZXV8NVXX3HjjTf69Tz7TZ8+na5du9ZK2yJy7BRYRcQvnE5nfZdwwqitvmzatClz5swhPz/fu+3tt9+mVatW1WrnxhtvpLCw0OcrISHB3+WKiHgpsIrUEcMwWL16tff5M888w5AhQ3xef+mll+jUqRMRERGMGjWKvLw87+ubNm1i1KhRxMfHExMTw/nnn+99bcuWLYwcOZL4+HiaNm3KQw89hMfjAQ5cOXrwwQdJSEggMTGRZ555xnvsL7/8Qt++fYmIiCAuLo6RI0d6X9u8eTNnnnkmMTExtGzZ0ue4/e1OmTKFRo0acckllwDw3nvv0blzZ6KioujVqxdLly71HpObm8vYsWOJioqiXbt2/PDDD4ftr08++YSWLVv6bFu+fDlRUVGUlpaydetWTj/9dCIjI4mJieHUU0+luLj4yN+Ev+Tn5zNp0iSaNm1KREQEvXr1YufOnQAUFhYyadIkmjRpQkJCAhMmTPB+H7Zt24ZhGLzzzju0atWKqKgorrzySm/AzM7OZsyYMURHRxMVFUWPHj3Yvn07AAUFBUycOJGkpCSSkpK4/vrrKSoq8mn3zTffpFWrVjRu3Nj7vRk6dCgxMTG0atWKV1991fsePB4P//73v0lMTCQ5OZkXX3zxqO87KiqKM888k/fee8+77c033+Rvf/tblfqtJo7Unw6Hg6uuuoq4uDgiIyPp1KkTP/30E8899xwzZsxg2rRphIWF0bFjRwCGDBni/QwuXLiQqKgopk2bRkpKCtHR0TzzzDNs2LCBPn36EBERwejRo719DHDZZZeRnJxMREQEPXr0YMGCBQCsWrWK66+/nrVr13qvGO/YsQM48ud5xowZtG7dmvDwcFJSUnjwwQdrrR9FTnqmiNQJwFy1apX3+dNPP20OHjzY5/WhQ4eae/fuNXNycsxu3bqZU6ZMMU3TNAsLC80mTZqYd999t1lYWGg6HA5z/vz5pmmaZlFRkdm0aVPz6aefNh0Oh7l9+3azY8eO5muvvWaapmm++eabps1mM5944gmzrKzMXLBggWmz2czNmzebpmma/fr1Mx966CHT7XabpaWl5vfff2+apmk6nU6zbdu25h133GGWlJSYa9asMZOSkswZM2Z427VareYDDzxgOhwOs6ioyPzyyy/NlJQUc+XKlabb7TY/+ugjMyYmxszMzDRN0zQvv/xyc/jw4WZOTo6ZlpZm9ujRwzzcjyGHw2HGxMSYixcv9m676aabzGuuucY0TdMcN26ced1115llZWVmWVmZuWTJEtPhcFTpezFmzBjzzDPPNNPS0ky3223+8ssv5r59+0zTNM2LLrrIHDdunJmTk2MWFhaal1xyiXnZZZeZpmmaW7duNQFz3LhxZn5+vpmWlmY2btzYfPPNN03TNM27777bPPfcc82ioiLT5XKZq1atMrOyskzTNM2//e1v5tChQ83MzExz37595uDBg81rr73Wp93Ro0ebOTk5ZlFRkZmenm7GxMSY77//vulyucy1a9eaSUlJ5rfffmuapmm+/vrrZuPGjc3169ebRUVF5pVXXmlaLBZzwYIFlb7nN9980+zSpYv5zTffmH369DFN0zSXLl1qdujQwVywYIEZGRnp3feKK64wb7nllkrbOdJrlTlSf7788stm9+7dzZycHNPj8Zh//PGHuWPHjsOeZ/DgwebTTz9tmqZpLliwwLRYLOYdd9xhOhwOc968eabVajXPOeccc8eOHWZubq7ZsWNH88knn/Qe/8Ybb5i5ublmWVmZ+dhjj5kxMTFmfn6+T/8c7Eif58LCQtNms3n/vuTk5JgrVqyocr+ISPUosIrUkaoE1q+++sr7/KGHHjLPPfdc0zRN87333jNbtmxpejyeCu1+8MEHZteuXX22vfLKK+Zpp51mmmb5/4gbNWrk83qrVq3MWbNmmaZpmoMGDTKvvfZac+fOnT77LF682IyIiPAJgQ8//LA5fPhwb7sxMTGm2+32vj5ixAjzmWee8Wmnf//+5ttvv226XC4zICDAXL58ufe1995777CB1TRN84YbbjCvu+460zRNs6yszIyLizN/+OEH0zRNc8KECeaoUaPMjRs3Hvb4yuzZs8cEzO3bt1d4LSMjw7RYLGZ2drZ328aNG0273W66XC5vsFy/fr339WuuucacNGmSaZqmed9995n9+vUzV69e7dOu2+02AwICzB9//NG7bcmSJWZgYKDpdru97R78+XjsscfM0aNH+7QzefJk86qrrjJN0zRPO+0089FHH63wvo4WWN1ut5mammquW7fOvPbaa83HH3+82oE1ICDAjIyM9H61adOm0n2P1p9vvPGG2bp1a3Pp0qU+n6PD1VBZYC0uLva+Hh8fb7700kve53fccYd56aWXVlqbaZpmVFSU9x9ElQXWI32eCwsLzeDgYPOll14y8/LyDnsOEfEPDQkQaUAaNWrkfRwaGkpBQQEA27dvp2XLlpXeUb9t2zZ+++03oqKivF+33347e/bs8e6TmJjoc8zBbb/xxhuUlpbSo0cP2rVrxwsvvADArl27SE5OJiAgwHtcixYt2LVrl/d5SkoKFsuBHyPbtm1j8uTJPrWsXr2atLQ0MjMzKSsro2nTpt79D35cmQkTJvDBBx/gcDiYM2cO4eHhDBgwAIDHH3+clJQUTj/9dJo1a8b999/vHQZxJNu3bycwMJAmTZpU2pcej4fmzZt76+/VqxcWi8WnPw/3fbrjjjsYOHAgY8eOpVGjRtxyyy2UlJSwb98+ysrKaNasmU9fOhwOMjMzvdsOrmnbtm3MmTPHpy+fe+450tPTAdi9e7dP/yUmJhIYGHjU92+xWJgwYQIvvvgiH330EZdffvlRjznUDTfcQG5urvfrjz/+qHS/o/Xn5ZdfzpVXXsn1119PXFwcV155pU9/HE14eDjBwcHe5yEhIT6f9ZCQEAoLC4HyIRT33HMPrVu3JiIigqioKPLy8o54viN9nkNDQ/niiy/47LPPSE1NZcCAAd4hBiLif7b6LkDkZBEaGuozxnJ/8KiKpk2bsmXLFkzTrBBaU1NT6dGjBz/++GON6mrZsiVvv/02pmmyZMkSTj/9dPr160fjxo3ZvXs3TqcTu90OlP8PfP/4SsAnrO6v5eabb+b666+vcB63243dbmf79u3eULF/nODh9O3bl7i4OGbPns27777LZZdd5n3/CQkJ3imO1q5dy/DhwznllFO44IILjthm06ZNcTgc7Ny5k9TU1Ar1WywWdu/eTUhISIVjt23bdsS2w8LCePTRR3n00UfZunUrI0eOZNq0adx2220EBASwbds273vftm0bgYGBxMXFefvh4P5MTU1lzJgxPuNND5acnOwdHwuQkZGBw+E4Yn37XXnllbRt25ZzzjmHxMRE1q9fX6Xjquto/QkwefJkJk+ezN69exk3bhz/+c9/eP755yt8to7VzJkzmTlzJt988w2tW7fGMAyio6MxTROo+FneX//hPs8Aw4YNY9iwYTidTqZNm8bo0aPJycnxe+0iopuuROpM9+7deeedd3C5XKxevZp33nmnyseec845OBwO7rvvPoqKiigrK/NezTn33HPZu3cv06ZNo7S0FLfbzR9//FHl+Tjffvtt9u7di2EYREVFYbFYsFqt9O7dm8TERO677z4cDge//fYbzz//PFdcccVh27rpppt4/PHHWblyJaZpUlxczLfffsuuXbuwWq2MHTuW++67j9zcXHbv3s3jjz9+1Pouv/xynn/+eb788ksmTJjg3f7BBx+wY8cOTNMkKioKq9WKzXb0f4MnJiZy3nnncf3115Oeno7H42HVqlVkZWXRqFEjRo8ezaRJk7xX3vbs2cMnn3xShZ6E2bNns3HjRjweDxEREdjtdmw2GxaLhfHjx3PPPfeQnZ1NVlYWkydP5vLLLz9suLn88suZP38+H330EU6nE6fTyerVq/npp58AGDduHC+++CJ//PEHJSUl3H333VUOSq1ateL777/3Xk2vjNvtprS01Oeruo7Wn/Pnz2f16tW4XC5CQ0MJCgryfg8TExP5888/vYHyWOXn5xMQEEBcXBxlZWU88MAD3ivj+8+Xnp5OSUmJd9uRPs979+7lk08+oaCgAJvNRkRERJU+fyJSMwqsInXk+eefZ9myZURFRXHXXXcdMfgdKiwsjG+//ZaVK1fSpEkTkpKSvHeF73/tu+++o1mzZsTGxjJ+/HifX2EfybfffkuXLl0ICwvjvPPO4/HHH6dr167Y7XZmz57NypUradSoEaNGjeIf//gH48ePP2xbI0eO5JFHHuHaa68lOjqa5s2b8+yzz3p/Vf/8888TFhZG06ZNOe2006r06+jLL7+cH374gW7duvlMv7Ry5Ur69+9PWFgY/fr14+qrr2bUqFEAXH/99Ye9Kgbw1ltvkZqaSs+ePYmKiuL666/3BpXp06d7f3UdERHBwIEDWblyZZX6cvPmzZx11lmEh4fToUMH+vXrxw033ADAs88+S7NmzejQoQMdO3akVatWPPXUU4dtKyUlhW+++YaXX36ZpKQkEhMTuemmm7xTUl111VVcdtllDBw4kBYtWtCtWzfCw8OrVCfAgAEDKh0Wsd8LL7xAcHCwz9d+++/eP/hr1apVlbZzpP7cf1U1KiqK5s2bExkZyZQpUwC45pprSEtLIyYmhs6dO1f5fR3OFVdcQceOHWnatCktWrQgODjY57cFp512Gn379iUlJYWoqCh27NhxxM+zx+Ph2WefJTU1lcjISF588UVmzZqlq6sitcQw/fXPVxERERGRWqB/CoqIiIhIg6bAKiIiIiINmgKriIiIiDRoCqwiIiIi0qApsIqIiIhIg6bAKiIiIiINWrVnOfZ4POzevZvw8PBKl4kUERERETka0zQpKCggOTn5qHMYVzuw7t69u8JyhiIiIiIiNbFz506fhTwqU+3Aun8lla1btxITE1OzysTL6XQyd+5czjjjDO967VJz6k//Un/6l/rTv9Sf/qX+9C/159Hl5+eTmppapVX6qh1Y9w8DCA8PJyIiovrViQ+n00lISIh33XE5NupP/1J/+pf607/Un/6l/vQv9WfVVWWIqW66EhEREZEGTYFVRERERBq0ag8JEBEREalvpmnidDpxuVz1XUqlnE4ndrud4uLik35IgM1mw263H9PsUgqsIiIiclxxOBxs27aNwsLC+i7liBITE9m8eXN9l9EghIWF0axZMwIDA2t0vAKriIiIHDc8Hg/r1q3DZrPRvHlzAgMDNS98A2aaJg6Hg7S0NNatW0eXLl2OOudqZRRYRURE5LhRWlqKx+OhefPmhIWF1Xc5UgWhoaEEBATwxx9/4HA4CA4OrnYbuulKREREjjs1uUon9Wf/98s0zZod789iRERERET8TYFVRERERI7Z9OnTOf3002ulbQVWEREREWnQFFhFREREpIKGNMetAquIiIiIHxiGwa5du7zPhwwZwv/+9z8ArrzySv7+978zbNgwwsPDOeOMM8jOzvbuO3/+fHr27ElERAStW7dm0aJFAGRnZzN+/HgSEhJo0aIFb731lk/7U6ZM8R538cUX43A4ANi4cSMDBgwgIiKCxMRE7rjjDu9x06ZNo0WLFsTHx3PZZZeRl5cHwMKFC2nVqhVTpkwhLi6OKVOmUFJSwqRJk0hOTqZx48Y88sgj3naKioq49NJLiYqKonv37mzatKkWerVcjae1cjqdOJ1Of9ZyUtrfh+pL/1B/+pf607/Un/6l/vSv46U/G3p9R/LBBx8wd+5c2rRpwznnnMOzzz7Lf/7zH/78809Gjx7NzJkzOfvss0lLS6OsrAyAyy+/nPbt27Nz5062bt3KaaedRrdu3ejcubO3za+//prIyEhOPfVUZs6cyd/+9jfuu+8+zjnnHBYtWkRxcTG///47APPmzePBBx/ku+++o2nTpkyYMIFbbrmF6dOnA7Bt2zasVivp6em4XC7++c9/kpOTw8aNG8nPz2f48OF06tSJc889l//85z/s3buXHTt2sHv3boYPH07btm2P2AcH58fqfC9rHFgXLFhASEhITQ+XQ8ybN6++SzihqD/9S/3pX+pP/1J/+ldD70+73U5iYmKF7Z9++imfffbZUY9v2bIl9957r8+2hx56iC1bthz2mPPOO4/Ro0dXu9ZDXXTRRd6gecEFFzB37lwA3n33XUaOHMm5554LQJMmTQDYs2cPCxcu5NNPP8Vut9OuXTvGjx/Pxx9/7G3nmmuuoWnTpgCcc845rFmzBijvp+3bt7Nnzx6SkpLo3bs3AO+99x4TJ06kQ4cOAEydOpWuXbvy5ptvAhAYGMjkyZOx2WzYbDbefPNNNm/eTFhYGGFhYdxwww3MmjWLc889lw8//JA333yTiIgIIiIiuOKKK/jxxx+P2AdLlizxBtXi4uIq912NA+vQoUOJjY2t6eHyF6fTybx58xg+fPhJv9awP6g//Uv96V/qT/9Sf/rX8dKfxcXFlS53WlxcTFZW1lGPj4uLq7AtLy/viMdWJ1gdycFBOyQkxLu07K5du2jevHmF/Xfs2EFpaSnx8fHebW63m0svvfSwbe4fZvDYY49xzz330LVrV5KTk3nggQcYOXIku3fvpn///t5jmjZtSmlpqfe4Ro0aYbOVx8N9+/ZRUlLiDbdQvtLYqaeeCkB6ejqpqane11JTU48aWE899VTvBc/8/Pwj7nuwGgdWu93eoD/Qxxv1p3+pP/1L/elf6k//Un/6V0Pvz8PVFhISUqULaZGRkZVuO9KxVf2NckhICCUlJd7ne/furdJxqamprFu3rsL2lJQUwsLCyMnJqfbys0lJSbzxxhuYpsnnn3/O2LFjycnJITk5mR07dnj327FjB0FBQcTExAD4nCcuLo7AwED+/PNP7+uHnmPnzp20bNkSgJ07dx61roM/X9X5nGlpVhERETnujR49usa/tj90iEBNdenShffee4/JkyczY8aMSq8EV2bcuHF07dqVOXPmcNZZZ3nHsLZs2ZJ+/fpx7733cs899xAQEMCvv/5KUFCQz1XPysyaNYv+/fuTnJxMVFQUhmFgGAYXX3wxV155JZdccglNmjThnnvuYezYsZUGYovFwhVXXMHtt9/O008/TUREBH/88QcFBQX07t2bCy+8kKlTp9K9e3fS09N5++23adOmTY367mg0S4CIiIiIHzz99NPMmDGDmJgYVq5c6fOr9yNp3rw5H330Effccw+RkZEMGzaM9PR0AGbMmMGuXbto0aIFCQkJ3HrrrT5XcQ9nxYoV9OjRwzvu9N133yUwMJAzzjiDu+++mxEjRtC0aVPsdjvPPPPMEd9TZGQkp5xyCjExMUyYMIGcnBwApkyZQmxsLKmpqYwbN47LL7+8Su+3Jgyzmou65ufnExkZSWZmZo3GsLrdHjw1XEf2RGMxDDweN3PmzGHEiBEN+lcwxwun06n+9CP1p3+pP/1L/elfx0t/FhcXs379etq3b6+bv48jlX3f9mfKvLw8IiIijnh8nQ8J8Jgm6Rl5OJ3uuj51g2K3W0lKqDiORkRERER81csYVqfTTdlJHlhFREREpGo0hlVEREREGjQFVhERERFp0BRYRURERKRBU2AVERERkQZNgVVEREREGjQFVhERERFp0BRYRURERASAbdu2YbPVy6ynR6TAKiIiIiINmgKriIiIyAnK5XLVdwl+ocAqIiIi4idbt27lnHPOITY2lqSkJJ577jncbjdTpkyhadOmJCYmcvvtt3uD5P3338+ll17KRRddRHh4OH369GHr1q0AlJSUMG7cOGJiYoiJiWHgwIHe8yxatIhu3boRFRXF4MGDWb9+vfc1wzB44YUXaN68OUOHDgXgxRdfpHXr1sTFxXHFFVdQVFTk3X/q1KkkJibSrFkzPv/887ropmpreIMURERERKop/62F5L+98Kj7BXRoTMLz1/hsy7j5NcrW7TrsMREThhBxxZCjtu1yuTjnnHMYO3YsH330EWVlZWzatImnnnqKRYsW8fPPP2O32xkzZgwvvfQSkyZNAuCTTz7hyy+/5N133+Wqq67i/vvv56233uKtt96iqKiItLQ07HY7y5YtAyArK4tRo0bx+uuvM3LkSJ555hlGjRrF+vXrveNP582bx5o1a7Db7Xz44Ye89NJLfPvttyQkJHD11VczZcoUnnjiCebMmcOLL77IokWLSEhI4MILLzzq+6wPusIqIiIixz1PUSnujLyjf2UXVjjWnV14xGM8RaVVqmH58uUUFBRw3333ERQUREREBD169OD111/noYceIj4+nqioKG6//XZmzZrlPW7YsGEMHToUm83GJZdcwpo1awCw2+1kZWXx559/YrPZvFdY58yZQ+fOnTn//POx2+3cfvvtFBcX89NPP3nb/Ne//kVERATBwcG8/vrr3H333TRt2pTg4GAmT57sPf+HH37ItddeS5s2bYiKiuJf//pXjb8HtUlXWEVEROS4ZwkNwpoQedT9rDFhlW470rGW0KAq1bBr1y6aNm2KxeJ7PXDHjh2cffbZGIYBgGmapKSkeF9PTEz0Pg4JCaGwsDxUX3755Wzfvp0xY8bgcDi4/vrrufvuu9m9ezdNmjQ5UJ/FQmpqKrt37/Zua9y4sc/5r7vuOm688UbvNqfTCUB6ejoDBgzwbk9NTa3Se61rCqwiIiJy3Iu4omq/tq/MoUMEaio1NZXt27djmqY3nAKkpKTw/vvv071792q1FxAQwAMPPMADDzzAhg0bGDJkCP369SM5OZk5c+Z49zNNk507d5KcnOzdduj5H3roIc4///wK50hKSmLnzp3e5wc/bkg0JEBERETED3r37k14eDgPPvggpaWl5Ofns3LlSq666iruvfde0tPTMU2Tbdu28f333x+1vQULFvD777/j8XiIiIjAZrNhtVo5++yzWbNmDZ999hkul4unn36a4OBgevbsWWk7V111FVOnTmXLli1A+VXVr7/+GoALL7yQ1157jU2bNpGXl8djjz3mvw7xIwVWERERET+w2WzMnj2bpUuXkpSURNu2bVm2bBl33HEH/fr149RTTyUyMpKRI0dW6Upmeno6o0ePJiIigl69enHdddcxcOBA4uLi+PTTT5kyZQqxsbF88sknfPrpp9jt9krbGTduHFdffTXnnHMOERERDB48mHXr1gFwzjnncN1113HqqafSuXNnzj33XL/2ib8Ypmma1TkgPz+fyMhIMjMziY2NrfYJnS43O9KyKXO6q33siSTAbqVJSgyYHubMmcOIESMO+0GTqnM6nepPP1J/+pf607/Un/51vPRncXEx69evp3379oSEhNR3OVJFlX3f9mfKvLw8IiIijni8rrCKiIiISIOmwCoiIiIiDZoCq4iIiIg0aAqsIiIiItKgKbCKiIjIccfj8dR3CVIN+79fB88PWx1aOEBERESOG0FBQVgsFrZu3UpKSgqBgYE1DkFS+0zTxOFwkJaWhsViITAwsEbtKLCKiIjIccNisdChQwe2bdvG1q1b67scqaKwsDDatGlTYdnaqlJgFRERkeNKYGAgbdq0wel04nK56rucSjmdTpYsWcKpp57aoOe1rQs2mw273X5MV8IVWEVEROS4YxgGAQEBBAQE1HcplXI6nTidTkJCQk76wOoPuulKRERERBo0BVYRERERadAUWEVERESkQVNgFREREZEGTYFVRERERBo0BVYRERERadAUWEVERESkQVNgFREREZEGrcaBtaSkxJ91iIiIiIhUqsaBdeXKlf6sQ0RERESkUjUOrEuXLvVnHSIiIiIilapxYP3tt9/IycnxZy0iIiIi9ebJJ59kzpw59V2GVKLGgdXj8fDDDz/4sxYRERGRevHhhx9y//33c9lll5GVlVXf5cghjmmWgO+//95fdYiIiIjUC4fDwV133cXrr7/O8OHDmTJlSn2XJIc4psC6efNmdu7c6a9aREREROrcM888Q6dOnRg7dizPPfccH3zwAevWravvsuQgtmNt4Pvvv+eyyy7zRy0iIiIide6uu+7yPk5MTCQjI6Meq5HK1PgKq8VSfujChQvxeDx+K0hERERE5GA1DqydOnUCICMjg/Xr1/utIBERERGRg9U4sA4aNIj+/fszefJk2rRp48+aREREROqMy+Xi4Ycfpnnz5oSEhDB48GA2btxY32XJQWocWPv27cu//vUv+vbti91u92dNIiIiInXC7XZz/vnn8/TTT3PNNdfw8MMPs27dOkaOHInL5aqXmn7//XcuuugiWrRoQUhICHFxcQwaNIgvvviiSscXFhYyZcoUzjrrLGJiYjAMg+nTp1e675VXXolhGIf9SktL8+67adMmLrnkEho3bkxISAjt2rXjgQceoLi42B9v+4iO+aYrERERkePVE088wXfffceKFSvo2LEjUH7j1aWXXsrChQs5/fTT67ym7du3U1BQwBVXXEFycjLFxcV89NFHjBo1ipdffpmJEyce8fjMzEweeOABmjRpQpcuXVi4cOFh973uuusqvEfTNLn++utp1qwZKSkpAOzcuZPevXsTGRnJpEmTiImJYdmyZUyZMoWVK1fy2WefHfP7PhIFVhERETkp5eXlMXXqVG699VZvWAXo378/AGvWrKmXwDpixAhGjBjhs23SpEn06NGDp5566qiBNSkpifT0dBo1asTPP/9Mr169Drtvv3796Nevn8+2xYsXU1xczKWXXurd9s4775Cbm8vixYu9fTVx4kQ8Hg9vv/02OTk5REdHV/etVtkxzcMKUFZWxtKlS3nkkUcoLCz0R00iIiIitW7GjBkUFBRUCID7hzoWFBTUR1mVslqtpKamkpube9R9AwMDadSoUY3PNXPmTAzDYPz48d5t+fn5QPnV54MlJSVhsVgICAio8fmq4pivsM6YMYNPPvkEgO7du3PGGWccc1EiIiIite3jjz+mQ4cOhIaGkpmZ6d2+f1Gk0NDQarfpdDrJy8vD6XSSn59PZmbmYe/1iYmJ8U4TWpmioiJKSkrIy8vj888/56uvvuLiiy+udk3V4XQ6+eCDD+jfvz/NmjXzbh8yZAiPPvooV199Nf/5z3+IjY1l6dKl/Pe//+Xvf/97jfqqOo45sA4YMMAbWOfPn6/AKiIiIg2e2+3mxx9/pKioiPj4+Er3ad68ebXbXbJkCUOHDq3Svlu3bvUJhYe6/fbbefnll4Hy+e/PP/98XnjhhWrXVB3ffPMNWVlZPsMBAM466ywefPBBpk6dyueff+7dfs899/DQQw/Vak1wDIHV6XTidDpp2rQpjRs3ZteuXaxbt46dO3ce5TK0gWmamKZZ01OfEEzTBNP03oHodDrruaITw/5+VH/6h/rTv9Sf/qX+9K+TrT83btxIUVERt99+e4VxqtOnT+f999+nffv21e6PDh068NVXX+Fyufjll1/o3r07NlvlcSs2NvaI7d90002MHj2a3bt389FHH+F0OikqKiIsLKzK9ezPGW63u0rv5X//+x92u50xY8ZU2L9x48YMHDiQMWPGEBMTw1dffcXUqVOJj4/nxhtvrHJN+1Wnb2scWBcsWEBISAhQPn5h165dALz66qv06NGj0mPCwsLo1as3BYWFlJaW1fTUJ4SgoABKSoL56acVAMybN6+eKzqxqD/9S/3pX+pP/1J/+tfJ0p+rVq0CICQkBIfD4fPa+vXriYqKYvPmzWzevLnG5+jSpQtutxu3213p6/Pnz69SO7GxsUycOJEpU6Zw2mmn8dhjj2EYRpWO3V//mjVrmDNnzhH3LSkp4dNPP6VLly4sX77c57VFixbx/PPPM23aNOLi4gAYPXo0O3bs4K677iImJoaIiIgq1bRfdabDqnFgHTp0KLGxsUD5nKwrV67E4/Gwa9cu7rnnniOMyTAIDwsjMLDyb97JIsBuJTg4mP79+zN37lyGDx+u+Wz9wOl0Mm/ePPWnn6g//Uv96V/qT/862fqzrKz8wtmQIUMYOHCgd3teXh4bNmzgyiuvrHCnflXbzc7Oxul08sMPPzBo0KDD9md8fDxWq7XKbe/evZsbb7yRVq1a0bZt2yods3LlSqA8PB/t/cyYMQOHw8Ett9xSYd8nnniCHj16MGHCBJ/tZWVlzJ8/n/j4eIYNG1bl9wIHbuSqihoHVrvd7v0GJCYm0q1bN1auXElmZiYbNmygS5culR7ndLm9k9GezAzDAMPw/prg4P6UY6f+9C/1p3+pP/1L/elfJ0t/7p+CqbS01Of9zpw5k7KyMm666Sbv9ptvvpmEhAT+/e9/A/Dll1/y2GOP8f3331do159jWA+1P2QXFxdX+Xu0P2dYrdajHvP+++8TFhbG+eefX2HfjIwMoqOjK2w/eIhndT831dnfb/OwDhs2zJvi58+ff9jAKiIiIlLfOnfujMViYcGCBd6ribt27eLBBx9kwoQJdO7c2bvvr7/+yj/+8Q/v87Vr1/q8frAuXbowb948XC4XK1asoHfv3ocdw3q4e34yMjJISEjw2eZ0Onn77bcJDg6mQ4cO1XqvVbFv3z6+/fZbxo0b5x3yebA2bdowd+5cNm7cSJs2bbzb3333XSwWy2H7w1/8Flh79+5NWFgYhYWFLFmyhOuuu67SNywiIiJS3xISEhg9ejTPPvssISEhREZG8swzz5CSksLzzz/vs++hAXXt2rUMHjy40najo6M5/fTTcTqdOBwOhg0bVu0rj9dddx35+fkMGjSIlJQU9uzZw4wZM9iwYQNPPvlkhZuuDMNg8ODBPitavfDCC+Tm5rJ7924AvvjiC+/9RjfffDORkZE+bbz//vu4XK4KswPsd8cdd/DVV18xcOBAJk2aRGxsLLNnz+arr77immuuITk5uVrvsbr8FlgDAgIYNGgQc+bMoaysjDVr1lRYOUFERESkoXjttde45pprvCFw7NixPPzww4SHh3v32bVrFy6Xy+dX92vXrmXSpEm1VtfFF1/M66+/zn//+1+ysrIIDw+nR48ePProo4waNcpn3/2LNiUlJflsf+KJJ9i+fbv3+ccff8zHH38MwGWXXVYhsM6YMYOEhITDruw1aNAgli5dyv3338+0adPIysqiefPmPPzww9x5553H/J6Pxq9Lsw4fPpzg4GCGDRtG48aN/dm0iIiIiF9FR0fz0UcfHXGfdevW0bFjR++9N3v27GHdunV06tSp1uq65JJLuOSSS6q07w8//IBhGEyePNln+7Zt26p1zmXLlh11n969ex91poHa4tfA2rJlS1q2bOnPJkVERETqxJNPPkn79u197pA3DIOCggI8Ho/3ZqzU1FSfq7D1acGCBVxyySWccsop9V1KrfJrYBURERE5Hn344Yfcf//92O12Nm3a5J26c/DgwaSkpNC+fXuaNWtGp06dGtTiCo8//nh9l1AnajWw7p/q4GSfwkpEREQaLofDwV133cXrr7/ORx99xJQpU7xLoAYEBPDNN9/Uc4VyuNn9j0lWVhYffPABN9xwA7/++mttnEJERETEL5555hk6derE2LFjee655/jggw9Yt25dfZclB6mVK6zr16/nf//7H6A5WUVERKRhu+uuu7yPExMTycjIqMdqpDK1coW1d+/ehIaGAuUrPlRnrVgRERERkYPVSmANCAjwTqhbVlbG4sWLa+M0IiIiInISqJXACuVLte43d+7c2jqNiIiIyDFxuVw8/PDDNG/enJCQEAYPHszGjRvruyw5SK0F1latWtG8eXMANm7cWO0JbEVERERqm9vt5vzzz+fpp5/mmmuu4eGHH2bdunWMHDkSl8tV3+Xx8MMPYxhGtRYq2D/rQXJyMsHBwfTp04d58+ZVuu+mTZu45JJLaNy4MSEhIbRr144HHnig0uGcv/zyC6NGjSImJoaQkBA6derEc889V+P3Vh21Nq2VYRgMHz6cV155BYB58+Zx7bXX1tbpRERERKrtiSee4LvvvmPFihV07NgRKL/x6tJLL2XhwoWHXaq0LuzatYupU6d67wuqqiuvvJJZs2Zx66230rp1a6ZPn86IESNYsGABAwYM8O63c+dOevfuTWRkJJMmTSImJoZly5YxZcoUVq5cyWeffebdd+7cuYwcOZJu3brx73//m7CwMLZs2cKuXbv89n6PpFbnYR0yZAjTp0+nrKyMBQsWcMUVV2BYrLV5ShEREZEqycvLY+rUqdx6663esArQv39/ANasWVOvgfWf//wnffv2xe12k5mZWaVjVqxYwXvvvcfjjz/OP//5TwAmTJhAp06duPPOO1m6dKl333feeYfc3FwWL17sff8TJ07E4/Hw9ttvk5OTQ3R0NPn5+UyYMIFzzjmHWbNmYbHU2i/oD6tWzxgWFub9phcWFlZpnVoRERGRujBjxgwKCgqYOHGiz3a73Q5AQUFBfZQFwA8//MCsWbN45plnqnXcrFmzsFqtPu8pKCiIq6++mmXLlrFz507v9vz8fKD8ivLBkpKSsFgsBAQEADBz5kz27t3Lww8/jMVioaioCI/HU8N3VjO1HpHPOOMMoLyzcnNza/t0IiIiIlXy8ccf06FDB0JDQ8nMzPR+7Q911f1VvNPp9GknPz/f5/nBX0cKfG63m5tvvplrrrmGU045pVo1rFq1ijZt2hAREeGzvXfv3gCsXr3au23IkCEAXH311axevZqdO3fy/vvv89///pe///3v3vf/7bffEhERQVpaGm3btiUsLIyIiAhuuOEGSktLq1VfTdXqkACAjh078o9//IPevXsTEhKC0+Wu7VOKiIiIHJHb7ebHH3+kqKiI+Pj4SvfZf/N4VS1ZsoShQ4dWad+tW7fSrFmzSl976aWX2L59O99++221zg+Qnp5OUlJShe37t+3evdu77ayzzuLBBx9k6tSpfP75597t99xzDw899JD3+aZNm3C5XJx33nlcffXV/N///R8LFy7k+eefJzc3l3fffbfadVZXrQdWwzC8CV5ERESkIdiyZQtFRUXceeedDB8+3Oe1N954g3fffZfOnTtXq80uXbp478Z3uVysWLGC3r17Y7NVjFuNGjWqtI2srCzuu+8+/v3vfx82SB9JSUkJgYGBFbYHBQV5Xz9Ys2bNGDRoEBdccAGxsbF8+eWXTJ06lUaNGjFp0iSgfFhncXEx119/vXdWgPPPP5+ysjJefvllHnjgAVq3bl3tWquj1gOriIiISEOzf7rNIUOGVLix6pFHHiExMZE2bdpUq83o6GhvW06nE4fDwbBhw7xjYqvi3nvvJSYmhptvvrla594vODgYh8NRYfv+X90HBwd7t7333ntMnDiRjRs30rhxY6A8iHo8Hu666y7GjRtHbGys95hx48b5tDl+/Hhefvllli1bVuuBtc5v8youLiY7q2p3uomIiIjUhqKiIqDiONW8vDwWLVrEmDFjqt1mWVkZe/bs8X7l5OT4PD/4y+2uOERy06ZNvPLKK/z9739n9+7dbNu2jW3btlFaWorT6WTbtm1kZ2cfsYakpCTS09MrbN+/LTk52btt2rRpdOvWzRtW9xs1ahTFxcWsWrXK55hDb85KSEgAICcn52hdc8zqLLAWFhby7LPPcs3VV/HxrJl1dVoRERGRCsLDw4HyfHKwt956i7KyMm644QbvNtM0eeWVV2jZsiWhoaF07tzZ5277/ZYuXUpSUhJJSUk0adKEv/3tbzRp0sS77eCvyo5PS0vD4/Hw97//nebNm3u/li9fzsaNG2nevDkPPPDAEd9X165d2bhxo3cGgP2WL1/ufX2/vXv3VhqcnU4ngHfhhB49enjrO9j+8bA1GbpQXXU2JCA4OJg1a9bgcDhY9/uv5OXmEBkVXVenFxEREfHq3LkzFouFBQsWMGLECKB8ov4HH3yQCRMm+IxfffDBB/nyyy+ZM2cOrVq1YvHixcTGxlZo81jHsHbq1IlPPvmkwvZ7772XgoICnn32WVq2bHnE93XhhRfyxBNP8Morr3jnYXU4HLz55pv06dOH1NRU775t2rRh7ty5bNy40Wf4w7vvvovFYvH2wdixY3nkkUd4/fXXOe2007z7vfbaa9hstjq5V6nOAqvVauX000/nvffewzRNVv60lNOGn1NXpxcRERHxSkhIYPTo0Tz77LOEhIQQGRnJM888Q0pKCs8//7x3vz179vDkk0+ycuVKWrVqBcDgwYMrbfNYx7DGxcUxevToCtv3z8Va2WuGYTB48GAWLlwIQJ8+fbjooou4++67ycjIoFWrVrz11lts27aN119/3efYO+64g6+++oqBAwcyadIkYmNjmT17Nl999RXXXHONdyhAt27duOqqq3jjjTdwuVze83344YfcfffdPsMMakud3nQ1bNgw3n//fUzT5KflSxgy7Ox6WS1BRERE5LXXXuOaa67hySefJCwsjLFjx/Lwww97hwtA+dLyffv29YbVhmT/cIZDp7F6++23+fe//80777xDTk4OnTt3Zvbs2QwaNMhnv0GDBrF06VLuv/9+pk2bRlZWFs2bN+fhhx/mzjvv9Nn3pZdeokmTJrz55pt88sknNG3alKeffppbb721Vt/jfnUaWBMTE+ncpQtrVq8mOzuTTRvX07Zdx6MfKCIiIuJn0dHRfPTRR0fcJzs7m6ioqLop6DD2Xz091A8//IBhGEyePNlne1BQEI8//jiPP/74Udvu3bs3c+bMOep+drudKVOmMGXKlCrV7G91fnnzzDPO9D5evvT7uj69iIiISJV17tyZBQsWsGnTJjweD7/88kuld+HXhwULFnDJJZdUezWs41Gdz8Pas1cvIiOjyMvLZd3va8jNzSFKN1+JiIhIAzR06FBuuukmBg0aRGFhIe3bt6/SFcm6UJUrqCeKOr/CarVa6du/fAyFaZr89OOiui5BREREpMqmTJlCeno6BQUFrFixgri4uPou6aRTL3c89e030Huz1YofF+N2u+qjDBERERE5DtRLYI2KjqF9h87ExiUwcMjpeNye+ihDRERERI4DdT6Gdb8LL7mCoKBgTWslIiIiIkdUb4E1JCT06DuJiIiIyElPlzdFREREpEFrEIF1757dfPn5LJxOZ32XIiIiIiINTL0NCdjvu3lfMverzwBISm5M955967kiEREREWlI6v0Ka4uWbbyPf9TKVyIiIiJyiHoPrM2at6JRo2QAtm/bQvruXfVckYiIiIg0JPUeWA3DoE//wd7nusoqIiIiIger98AK0L1nXwICAgH45ecfKSkprueKRERERKShaBCBNSgo2HuzVVmZg59XLK3nikRERESkoWgQgRWg/4Ch3sfLlizE49FyrSIiIiLSgAJrYqNkWrVuB0BWZgYbN/xezxWJiIiISEPQYAIrQP+Bp3kfr1n9Uz1WIiIiIiINRb0vHHCw9h060617H07p0p32HbvUdzkiIiIi0gA0qMBqsVi45LKr67sMEREREWlAGtSQABERERGRQzX4wFpW5qjvEkRERESkHjWoIQEH27B+LUsWzSc3J5vb7piCxdLgs7WIiIiI1IIGmwK/n/8NGzf8TsbedDZvWl/f5YiIiIhIPWmwgfXUg6a4WrJofj1WIiIiIiL1qcEG1vYduxAZFQ3AH+t/Iyszo54rEhEREZH60GADq9Vqpd+pQwAwTZMlixbUb0EiIiIiUi8abGAF6N13IHa7HYCfViympKS4nisSERERkbrWoANraGgYPXr1B6DM4WDFj4vruSIRERERqWsNOrACDBg0zPt4yaLvcLtd9ViNiIiIiNS1Bh9Y4xMa0b5DZwDycnNY++sv9VyRiIiIiNSlGi8c4HQ6cTqdNTjSwDRNTNOs8hEDBp/O+nW/0rxFa8LDI6t1bENlmiaYJi5X+RXjmvWlHGp/P6o//UP96V/qT/9Sf/qX+tO/1J9HV52+Mcxqpr/8/HwiIyOZOXMmISEh1SosLCyMXr16s2XHPkpLy6p8nGma7MvYQ0JiUrXO15AFBQXQskk8P/20gsLCwvouR0RERKROFRcXM378ePLy8oiIiDjivjUOrOnp6cTGxtagPIPtadmUOd01OPbEEWC30jQlBpfLydy5cxk+fLh3RgSpOafTybx589SffqL+9C/1p3+pP/1L/elf6s+jy8/PJy4urkqBtcZDAux2e42+AU6XG8MwMAyjpqc+IRiGAYaBzVb+Lahpf0rl1J/+pf70L/Wnf6k//Uv96V/qz8OrTr80+JuuDmWaJhv/WMfyZYvquxQRERERqQM1vsJaHzweD9Oee4SdO7YREBhI5649CA6u3jhaERERETm+1PgKq+mo+7veLBYLKY2bAlpIQERERORkUePAWvjeEn/WUWUDBp3ufayFBEREREROfDUOrEX/+wHX7mx/1lIl8QmJtO94YCGBNat+qvMaRERERKTu1HxIQKmTnCc+82ctVTZ46Jnex9/P/+aEWEhARERERCp3TLMEFM/7lZKlf/irlipr1rwVTZq1AGDPnt38sf63Oq9BREREROrGMU9rlf3Ix5jOuh1HahgGQw66yrpw/td1en4RERERqTs1Dqz2Tk0AcG3NIH9G3c+J2r5jF+ITGgGw9c9NbN/2Z53XICIiIiK1r8aBNeLWEfDXalV5//0GV0ae34qqCovFwuChZ3ifb9ygYQEiIiIiJ6IaB9aAtimEXdQPALPYQe5TX/itqKrq1qMPvfsO5ObbJjP8rFF1fn4RERERqX3HNIY16uYRWCJDwGbBGheO6fH4q64qsdnsXDD2chqnNqvT84qIiIhI3TmmpVmtUaHETh2PLSWWgJaN/FWTiIiIiIjXMQVWgJBBHf1Rh194PB6yMjO8N2OJiIiIyPHvmANrQ7H6lxXM/3YOubnZ3P3vRwgODqnvkkRERETED455HtaDmWUu8l6ZR9aU9/zZbJVs2fwHe/fsxlFayvKlP9T5+UVERESkdvgtsJqmyd5rp5H7/BwKP15e5ytgDR56BsZf02wt+n4ezrKyOj2/iIiIiNQOvwVWwzAIPa+393n2w7MwHU5/NX9UcfGJnNKlBwCFhQWsWL64zs4tIiIiIrXHr0MCwkb3JrB7cwBcOzLJe/07fzZ/VKedPsL7+PsF3+By1V1gFhEREZHa4dfAalgsxNx7EdjKm8177Vuc2zL8eYojSkpuTIeOXcrPnZvDyp9+rLNzi4iIiEjt8GtgBQhonUTEhCHlT5xush+ahWma/j7NYR18lXXhd1/hdrvr7NwiIiIi4n9+D6wAkdedgTU5GoDS5Zso+vKX2jhNpVKbNqd12w4AZGdnsmbVijo7t4iIiIj4X60EVktIIDH3XOB9nvP4p7jzimvjVJUaNvwc7+MF331dp1d4RURERMS/aiWwQvkKWCGndwbAk11I/vQFtXWqCpq3aE2Llm3o2q03l06Y6J3uSkRERESOP7W60lX0XWMo/Wkz4RefSsS1p9fmqSq45vpbsVpPmIW8RERERE5atZrobI2iSPnm31hCg2rzNJVSWBURERE5MdTakADvCeohrFbG4/FoLKuIiIjIcajWA+uhnFv2QFp2nZ3P5XLx0/LFPPHIfWzauK7OzisiIiIi/lFngdV0usj97zfsG/c01mdng6durnZuWLeWWe+/TVZmBvO++UJXWUVERESOM3UXWF0eir74GZxujHW7sM1dXSfn7dCpC4mNkgHYse1PNv6hq6wiIiIix5M6C6yW4ABi7x/rfW5/eyFGZn7tn9di4fQzR3qfz/v6M11lFRERETmO1OkY1qDerQkZ3RsAo6SMgJfnQh2Ex06ndCMpqTEAO3ds44/1v9X6OUVERETEP+r8pquI287FjA4FwPrzFqyL19f6OStcZf3mc11lFRERETlO1HlgtUSE4LnhTO/zgNe+g/ySWj9vx1O6kpySCsCundtZv+7XWj+niIiIiBy7Og+sAGb/drj6tgHAyC8m4M35tX5OwzAOGcuqq6wiIiIix4N6CawAZdcOxwwNBMC28Dcsq/6s9XN26NiFlMZNANidtpN1v62u9XOKiIiIyLGpt8BKTBhlVw4FwJMSA2HBtX5KwzAYfuYoLBYLvfoMIOmvIQIiIiIi0nDZ6vPk7mGdKXN5cA3tBIH2Ojlnuw6ncOfkh4mOia2T84mIiIjIsanXwIph4DqrWx2f0lBYFRERETmO1N+QgMNxe8DhrNNT6uYrERERkYarQQVWIy2LwMkzCHh1Xp2cz+VysmTRfF564THcblednFNEREREqqd+hwQczOEk6J6ZGHnFsHE3rn5t8fRoWaunnPX+26xauRyAn5YvoW//wbV6PhERERGpvoZzhTXQTtn4gd6nAdO+hqLSWj1lv1OHeh9/N/dLnGVltXo+EREREam+hhNYAffwLri7NgPAkl1IwBvf1er5mjZrQYdOXQDIz89l6eIFtXo+EREREam+BhVYMQzKbjobMyQAANv837D8vKVWT3nm2aMxDAOAhfO/pqSkuFbPJyIiIiLV07ACK2DGRVD2t2He5wHTvobC2hsa0CgphW49+gBQXFzEDwvr5oYvEREREamaBhdYAdzDTsHdvQUAlpzCWp81YPiZo7BarQAs/v5bCgrya/V8IiIiIlJ1DTKwYhiU3XgWZkggALYf1mFdsr7WThcTG0fvvuU3fJWVOVjw7ZxaO5eIiIiIVE/DDKyAGRtO2cTh5Y8NMNKya/V8w4afgz2gfOzsj0t/ICc7q1bPJyIiIiJV03DmYa2Ee1AHnFv24O7dGk+nJrV6rvCISAYMHMYPC+fRt/9gAv4KryIiIiJSvxp0YMUwcF417Oj7+cng086kT79BRMfE1tk5RUREROTIGnZgPRyPCRbD780GB4cQHBzi93ZFREREpOYa7BjWSnlMbF/8TOA9M8DprpNTmqZZJ+cRERERkcodV4HV/uo8At74DuuGNOzvL67Vc5WVOZg/70teeuFxPB5PrZ5LRERERA7vuAqsrtM7Y1rLS7Z9shzLul21dq6Z77zKN199xratm1m18sdaO4+IiIiIHNlxFVjNlo1wjhsAgOExCXhuNhQ7auVcAwed7n38zVef4Swrq5XziIiIiMiR1fimK6fTidPprMGRBqZp1nhsqPO83lh/3oJ1QxqWvXnYX/uWsptH1KitI2nRqi3t2p/ChvVrycvNYdEP3zJ02Nl+a980TTBNXC4XQA37Ug61vx/Vn/6h/vQv9ad/qT/9S/3pX+rPo6tO3xhmNZNjfn4+kZGRzJw5k5CQ6t1RHxYWRq9evdmyYx+lpTW/Ymndl0/8fbOwlJa/0ZyJp1HSr3WN2zucfRl7eP3lpzFNk4CAQG64+S5CQsP80nZQUAAtm8Tz008rKCws9EubIiIiIseL4uJixo8fT15eHhEREUfct8aBNT09ndjYmsxXarA9LZuyY7zL3/r97wQ9+yUAZnAAJU9egdko+pjarMxHH7zDT8vLb/DqP2Aoo8Zc4pd2A+xWmqbE4HI5mTt3LsOHD8dut/ul7ZOZ0+lk3rx56k8/UX/6l/rTv9Sf/qX+9C/159Hl5+cTFxdXpcBa4yEBdru9Rt8Ap8uNYRgYxrHNo+oZ0gnXmm3YFv6OUVJG4FOzcUy9FOzWY2r3UGecPYrVq1bgLCvjx6Xfc+rA04iLTzzmdg3DAMPAZiv/FtS0P6Vy6k//Un/6l/rTv9Sf/qX+9C/15+FVp1+Oq5uuDlU2cTiepPKrqkZRKUaO/3+1HhERxaAhZwDg8Xj48otZfj+HiIiIiBzecR1YCQ6k7B8jcZ12CqVPXoGZEFkrpxk89AzCI8rbXvfbGjZv2lAr5xERERGRio7vwAp4WiWVzxIQHFhr5wgMDOLsc87Hbrdz+pkjadK0ea2dS0RERER81XgMa4Pm9oClfIyov3Tr0YdWrdsRGeX/G7tERERE5PCO+yushzIy8gi8dya22Sv92q7FYlFYFREREakHJ9YV1vwSgm6fjlFYimXzHtwdUzFbHPsd/YdTUlJMcHD15qIVERERkeo5sa6wRgTjGtYZAMPlJvCJz6DI/0u3OkpL+Wr2x0x94C4y9qb7vX0REREROeDECqyA89JBuFslAWBJzyHgxa+ghsvAHs7SJQtYOP9ryhwOvvz8Q7+2LSIiIiK+TrjAit1K2T9HYYaWzxpgW/YHtjm/+PUUpw48zTuedcP63/hj/W9+bV9EREREDjjxAitgJkZRdvM53uf26fOxbPLfr+4DAgI5+5zzvc+/+OwD3G6X39oXERERkQNOyMAK4O7TGud5vQAwXB4CHv8UCkv91n7X7r1p0rQFAPsy9rBk0Xy/tS0iIiIiB5ywgRXAedlg3G2TAbDsyyfguS/9Np7VMAzOO/8SjL/mep33zRfk5+X6pW0REREROeCEDqzYrJT98zzM8GAAzORo8PjvBqzGqc3o3XcgQPkNWF/M8lvbIiIiIlLuxA6sgBkXgeO2c3H863ycV54GVv++5TNHjCYkJBSA1b+sYMvmP/zavoiIiMjJ7oQPrACebi1w92ldK22HhoZx1jljvM81llVERETEv06sla6qwdixDzMl1i9XXHv1GcCaVT/Rpl1HBgw63Q/ViYiIiMh+J2VgtX67hoCX5+Ea1Qvn5YOPuT2LxcK1N/zDewOWiIiIiPjPSTEk4GBGWjYB//0Gw+XG/vGPWJf5Z8ypwqqIiIhI7TjpAquZElN+89VfAp6bg7Ery+/nSd+9i107t/u9XREREZGTzUkXWAFc5/bANagDAEZpGYGPfAzFDr+07Swr44tP3+e5px7i/Zlv4HJpBSwRERGRY3FSBlYMg7IbzsTTNB4AS1o2Ac/N8cuiAlabjW1bN+PxeMjYm86ihfOOuU0RERGRk9nJGVgBggJw3DUGMzQQANvyjdg+Xn7MzVosFs6/6DLvmNZv584mKzPjmNsVEREROVmdvIEVMJOicdw2EvOv+6XsM3/AsnrrMbeb0rgppw4aBoDL5eSTj2Zi+mlJWBEREZGTzUkdWAE8PVrivHgAAIbHJOD5OVB27ONOzzhrFJFR0QBs+mMda1b9dMxtioiIiJyMTvrACuC6qD/uni3xxIZTdtcYCDj26WkDA4MYff447/MvPn2f4uKiY25XRERE5GSjwApgMXDcei6lT1yBp02y35rt0KkrHU/pBkBhYQFfzf7Yb22LiIiInCwUWPcLDYKoUL83O2rMJQQElt/YteLHRWzbutnv5xARERE5kSmwHo7bg336/GNeCSsqKpqzzh4NQFJyY+x2ux+KExERETl5HPtgzRORw0ngo59gXbUVM3A1pY2iMZsn1Li5fgOGEhAURPcefbFarX4sVEREROTEpyuslQmwYUaEAGA4nAT+30eQW/MbpiwWC716n6qwKiIiIlIDCqyVMQzKbjwLd+skACz78gl8/FNwuv12Crfbhdvtv/ZERERETlQKrIcTYKPsX2PwRIcBYF23i4BX5/ll+da0XTt4+omH+OqrOcfcloiIiMiJToH1CMyYcMr+NQbTXv6rfNu8NdhmrzymNnNzsnnx2f8jbddOZs6YwZ49e/xRqoiIiMgJS4H1KDxtkim76Wzvc/v0+Vh+qvnUVFHRMfTuOxAAh8PBSy+9pGVbRURERI5AgbUK3IM74rywH1C+fGvgU59jbN1b4/bOPud8oqNjAFi7di0bN270S50iIiIiJyIF1ipyjhuIa0A7ADyN4zCPYZGBwKAgxo67wvt8+fLlZGVlHXONIiIiIiciBdaqshiUTRqB86L+OB4aB3/djFVT7dp3YujQ0wBwOp288sorGhogIiIiUgkF1uoItOMcPxAC/bNa1ZV/+xtRUVEA/PzzzyxYsMAv7YqIiIicSBRYj1VRKda5q2t0aFhYGBMnTvQ+f/XVV9m3b5+fChMRERE5MSiwHgNjby5Bd/+PwP9+g+2rX2rURp8+fWjVqhUAJSUl/P777/4sUUREROS4Z6vvAo5nlt93YtlZfrOU/bVv8cRH4unZstrt9O3bl9DQUCZMmEDr1q39XaaIiIjIcU1XWI+B+7RTcI7uDfw13dUTn2HZlF7tdgIDA7nvvvsUVkVEREQqocB6jJyXD8F1avl0V4bDSeBDszDSc+q5KhEREZEThwLrsbIYlP39HNwdUwEw8osJfOADyC2qcZNut5v333+fzz77zF9VioiIiBy3FFj9IcCG4+7z8aTGAWDZk0vgwx9BaVm1m3K5XPzrX/9ixowZvPXWW2zfvt3f1YqIiIgcVxRY/SU0CMd9F+GJLV9QwLo5ncAnPgO3p1rN2Gw2OnToAJSH16effhqXy+X3ckVERESOFwqsfmTGReD491jMkMDy5+EhUIPVqy699FJSU8uHGPz555+89957fq1TRERE5HiiwOpnZtN4HHefj/Oi/pT9fQTYrNVuIyAggNtuuw2rtfzYDz/8kN9++83fpYqIiIgcFxRYa4GnU5PyJVwNo8ZttGrVivHjxwNgmiZPPfUUhYWF/ipRRERE5LihwFpHjO37sC77o1rHnH/++ZxyyikAZGZmMm3aNMwaDDEQEREROZ4psNYBy4Y0gu6ZQcCTn2P55c8qH2e1WrntttsICyu/kWvx4sXMnz+/tsoUERERaZAUWOuAdfF6jCIHhttD4KOfYNmwq8rHxsXFMWnSJO/zuXPn6iqriIiInFRs9V3AycD5t9MwsguxLfsDo8xF4EOzcD9yOaTEVOn4/v37c8YZZxAUFMSECRMwjmFsrIiIiMjxRoG1LlgtlN12LkaxA+uabRhFDqz3vYvrrZuxNa5aaL3xxhuxWHRBXERERE4+SkB1xW7DcdcY3G2SADByisi68VXc+/KrdLjCqoiIiJyslILqUnAAjnsu8i7h6t6VRdZNr2Etdla7qR07dnDfffeRl5fn7ypFREREGpQaDwlwOp04ndUPWmBgmubJe+NQeBClUy4iePJMjIw8XJv30PKtIsrOObPKTaxatYrHHnuMsrIynnzySe655x5dgf3L/s9kzT6bcij1p3+pP/1L/elf6k//Un8eXXX6xjCrmRzz8/OJjIxk5syZhISEVKuwsLAwevXqzZYd+ygtLavWsSeakNxioh/6FDOrgJ0jW5PVJ6XKxxYXF/PJJ59QWloKQM+ePenSpUttlSoiIiLid8XFxYwfP568vDwiIiKOuG+Nr7AOHTqU2NjYGhxpEB4WRmCgu6anPiEExMUQ9/J1lP62nSzbPoYPH47dbq/y8a1bt+bBBx/ENE1++eUXzjvvPDp06FCLFR8fnE4n8+bNq3Z/SuXUn/6l/vQv9ad/qT/9S/15dPn5VbuPB44hsNrt9hp9A5wuN4ZhnPRTMxmGgb1lI+wtEmDOHJ/+NE3zqP3Ts2dPLr74Yt577z08Hg/PPPMMzzzzDJGRkXVRfoNX08+nVE796V/qT/9Sf/qX+tO/1J+HV51+0cDHBqZo7mr23fYmptN11H0vvvhi79KtWVlZPP3003g8ntouUURERKROKbA2IEVzfiHzjrcp+W4tmXe9g+k68rAJq9XK7bff7r2q+ssvv/Dxxx/XRakiIiIidUaBtQGxxkdg2MtHaRTP+5XMe2Ziuo98xTQmJobbb7/dO4Tgf//7H+vWrav1WkVERETqigJrAxLUqxXxz10FdisAxXN+IWvK+0cNrV27duXiiy/2Pt+2bVttlikiIiJSp7Q0awMT3L8d8U9dyb7b3gSXh6LPVoBpEvvAJRjWw//74uKLLyY9PZ0zzzyTTp061WHFIiIiIrVLV1gboJAhnYh7bALYyr89RZ//RNa/3z3ildb941kVVkVEROREo8DaQIUO70L8E1ccCK1f/Ezm5BlHvRHrUNWZ40xERESkIVJgbcBChnUm/skrwVY+prVk/lqcf+6t0rEej4cPPviAiRMnsmPHjlqsUkRERKR2KbA2cCGnnUL801dihAcR/8K1BLRJrtJxc+bM4X//+x/FxcVMnTqVoqKiWq5UREREpHYosB4HQoZ0ovHX9xHcp3WVjxk+fDjNmzcHYPfu3TzzzDNaVEBERESOSwqsDUBYWNhR97FEBPs8N02Toq9+OeyKWIGBgdx9993etpcvX86sWbOOvVgRkQYkMDCwvksQkTqgaa3qidViYBhgmga9evUGDJzVuKEq/6VvKHx5HoED2xPz2ASMoIrr8cbGxXPrrbfx8MMPYZomM2bMoGmzZnTv3sOP78R/LIaB9QhTd4kIuN0ePKZZ32U0EAa9e/ep7yJEpA4osNYTi8WC6YH0jHyycnIJDwvzrlZ1VHtysL65AANwLFpP2nUv4bn3QgipeKUhNrEZZ404j6++/BTTNHnyiSe59fZ7SGyU5N83dIzsditJCZFY67sQkQbOY5qkZ+ThdFZvxpATkd1mJSq84j/WReTEo8Baz8qcLkpLywgMdFc9sMZGYPn3RQQ+/BFGaRmWX7dj3jsTx7/HQlhQhd0HDT2LnTu289vaVZSWlvDaK89x0y13ExIS6ud3IyJ1wel0U6bAimmagAKryMlAv389Tnk6NcHxn4sx/wqo1o3pBP17JuRWnA3AYrEwdvzfSEpqDEDmvgw+/vB/dVqviIiISE0psB7HPG2SKX1wHGZU+ZVSy7Z9BN0zEyOz4mIBgYFBXHH1TYSGhROf0IizRoyu42pFREREakaB9ThnNkug9KHxeOLCAbDsziZw8gyM9JwK+0bHxHLNdbdy0y3/Ii4+sa5LFREREakRBdYTgJkSg+PhS/EkRQNg2ZdPwHNfQiV3EienpBIcHFLXJYqIiIjUmALrCcJMiCy/0tokDk90GGV/PweqcBOXs6yMTz+aSXZ2Zh1UKSIiIlJ9miXgRBITRumD4zEKijH/utp6JAX5eUx//UV27dzGn1s2csPNd+rqq4iIiDQ4usJ6ookIxkyJ9d1W5sLy85YKu1ptNkpLiwHYu2c3/3vrZdzuylfOEhEREakvCqwnOreHgKe+IOjhWdg+We4zrjUkJJS/XXMzIaHlswxs3riejz+c8dfchiIiIiINgwLrCc764x/Ylm8EIODthdjfnA+eA4E0Lj6RK666CZutfHTIzyuWsOC7r+qlVhEREZHKKLCe4Nz921F26SDvc/sXPxPw7Gw4aJWcZs1bMXbc37zPv5nzKat/WVGndYqIiIgcjgLric4wcF3YD8eNZ2FaymcNsP2wjsCps6CkzLtbl269OGvEGO/zD96dztY/N9V5uSIiIiKHUmA9SbiHd6HsrjGYAeW/+reu3kbQvTMxsgu8+wwZdha9+gwo39/t4u03ppG5L6Ne6hURERHZT4H1JOLu3RrHlIsxQwMBsPy5l8C73sHYvg8AwzAYc+F4WrdpD0Bqk+aEh0fUW70iIiIioMB60vF0aEzp/12GJ748iFoyC7B/tMz7utVq49IrruPMs8/jiqtvIjAoqL5KFREREQEUWE9KZmocpY9ejrtVIzwtEim74Syf14ODQzht+DlYrdZ6qlBERETkAK10dbKKDsPx4DgodUJwwFF3z8nJ4scl33PmiNFYLPp3joiIiNQdBdaTWVBA+ddBjMx8bJ+uwHnFELCXfzz2pKfx+ivPkp+Xi9PpZOTosRiGUQ8Fi4iIyMlIl8rkgCIHgQ9+iP3LlQTe/wHkly/bmpOTRWFBPgBLFn3H9/O/qc8qRURE5CSjwCpelm0ZGHtyAbCu20nQne9g7NhH+w6dOX/s5d79vvryY35asaSeqhQREZGTjQKreHk6puJ4cBxmdCgAlr25BP3rf1h+3kKv3qdy1jkHFhb4+IN3+O3XVfVVqoiIiJxEFFjFh6dNMqWPTcDTIhEAo6SMwKmzsH22giFDz+TUgcPK9/N4mPnOq2z84/f6LFdEREROAgqsUoEZF0Hpw+Nx9WsDgGFCwPQFBL74NeeOGEP3nn2BA6thaQlXERERqU0KrFK5oADK/jka59j+3k22+WsJ/s+HXHj+ZXQ6pRsATqeTN199npzsrPqqVERERE5wCqxyeBYD57iBOG4fhRlQPsWVp2Mq1qAAxl1+DW3adQSg/8ChREXH1GelIiIicgLTPKxyVO4B7XE0isL63Vqc4wYCYLPZufzK6/nt11XeIQIiIiIitUGBVarE0yoJT6skn20BAYH0DG2Mx+0B64GL9R6PR6thiYiIiN8oVUiNWdbvIvDf7xL4nw8gr3yRgS2bNvDskw+Sm5Ndz9WJiIjIiUKBVWrG6Sbg6S8w3B6sa7cTdMdbpH+/gjdfe5496Wm8Mu1JcnNz6rtKEREROQEosErN2K2U3TrywCID+/Jp8cIPDHLGA5CVtY9Xpj1JnkKriIiIHCMFVqkxT4fGlD5xBe62yQAYLjejtwZzeWEqNhOyMjN4WaFVREREjpECqxwTMyYcx4PjcZ7Z1but974Abt/XjBinjazMDF7575Pk5+XWW40iIiJyfFNglWNnt+K8/kwcN52NabcC0LjIyl3pqbQtCSZzX/mVVoVWERERqQkFVvEb9+mdKf2/y/AkRgEQbLFjxIQDkLlvLy9Pe0KzB4iIiEi1KbCKX5ktG1H65BW4+rSh7LozGP3P24iJiQMgOyuT9N276rlCEREROd4osIr/hQZRdtdo3MM6Ex0dy8Sbbic+IZFLLr6SjlYt4SoiIiLVo8AqtcMwvA+jo2O59Z9T6PlrIYGTZ2B7fzG4PfVYnIiIiBxPtDSr1ImADbuxz15Z/vi9JVh/28nygXHEtm1FatPm9VydiIiINGS6wip1wtMhlbLLBmFayq+8Wn/bQZeXf2LZ4/9l65+b6rk6ERERacgUWKVuWAxcF/TD8eA4PHHlMweEeqxcmRZL7gNvs2HVqnouUERERBoqBVapU54OqZQ+fRXOPq28207NCyP2gc/48d0v6rEyERERaagUWKXuhQXhvOt8Sq47HVf5OgMkOQNIenwB8154p35rExERkQZHgVXqh2FgntUDx1NXkRNtB2BHgIMXv/mQt956C9M067lAERERaSgUWKVeGU3iCfjvzWzpGstb8Rl4DPjoo494/vnncbvd9V2eiIiINAA1ntbK6XTidDprcKSBaZon/RU0E9P7X+Dk7o8AG00enMiFa1fw2muvYpomK1as4Pz2fQlctIWIf5yLJTSoSk3t/0zW7LMph1J/+pd/+lM/Q/fb3wUul6t+CzlB6O+7f6k/j646fWOY1fypl5+fT2RkJDNnziQkJKRahYWFhdGrV2+27NhHaWlZtY490URGhNKqWSPWb95FcbGjvsupd0FBAbRsEs/MmTP47rvvGDH0dE59fxuB2aU4ooPYcUE7ippF1XeZIvVKP0N97f+58dNPKygsLKzvckSkmoqLixk/fjx5eXlEREQccd8aB9b09HRiY2NrUJ7B9rRsypwn9697Q0MCSE6IZFtaNpmZ2URERGActDrUySbAbqVpSgxgUlxcjPWPvWRPeh2z5K//KRsGoZcPIuLGMzECDv+LAafTybx58xg+fDh2u71uij+BqT/9yz/9qZ+h+9ltVhJjg7HbbdhsWgfnWOnvu3+pP48uPz+fuLi4KgXWGv8Nt9vtNfoGOF1uDMM4qcMZgIHh/S9w0veJYRhgGNhtNiIjI6F3JIEf3UHWvTNx/LIVTJOit7+nbOlG4h65lIC2KUdsr6afT6mc+tO/jqU/9TP0gP1dYLPZ9Pn0I/199y/15+FVp19005U0WPbUOBLfmMSmAYk4/xrt69ycTvolT5P36jxMXWESERE5KSiwSoPm8rhZGFvIEym7SAv4a6yvy03uc3PYc+kzlG3ZU78FioiISK1TYJUGzW63M2XKFNqdNYAnkncxNzIHz19XW8u27NGvRUVERE4CCqzS4Nntdm655RbGXjqeL2KyeTI5jd32Mn5uY6MsMay+yxMREZFapsAqxwXDMLjkkku47bbb2B3q5vGUnbxT+Dt33nknu3fvBsBTWkbB9IUYGtsqIiJyQlFglePK0KFDuf/++wkKD8NjwM6dO/nnP//Jli1byJv2DQXPzaHtiz/jWL2tvksVERERP1FgleNO586defLJJ0lNTQUgMTGRpJAoCt5dBEBQZglZV00j68EP8eSX1GepIiIi4gcKrHJcSkpK4vHHH2fYsGFMnjyZkKRYGr17G/ZOqd59Cj9Yyu7z/o+ir1dpGUsREZHjmAKrHLdCQkK45ZZbiI+PByCgVRJxb9zIltMaQ1D5ZMTuzAIy73ibfTe9iistuz7LFRERkRpSYJUTSqmzjLdyf+XpNll4ejXzbi9ZtJ7dox8hf8YP9VeciIiI1IgCq5xQpk+fTk5ODn8W7OP2nO/ZdUU3rAmRAJilTsxSZz1XKCIiItWlwConlIsvvpiEhAQAXG4Xj/4wi8/PjiLk4v7Y2yYTMWFI/RYoIiIi1abAKieUmJgYRowYwYgRI7zbvv5hPg9nLsF4ahyG3eqzf+4LX1Hw/hJMt6euSxUREZEqUmCVE47VauXqq6/m9ttvJzAwEIA///yT2/91JytXrvTuV/ZHGnmvziP7oVnsGfeU5m4VERFpoBRY5YQ1ePBgHn/8cZKSkgAoLCzkgQce4P333wegZPEG8JRPd1W2Po09lz9L5r3v4s4qqLeaRUREpCIFVjmhNWvWjKeeeoo+ffoAYJomdnv5lFeRVw8j8a2bsbdJ9u5f9NkK0kZOJX/GD5guLfEqIiLSECiwygkvNDSUu+++mwkTJtC7d29Gjx7tfS2oewuS3v8H0XefjxEeBIBZUErOI5+QftGTlCz9o56qFhERkf0UWOWkYLFYuPDCC5k8eTIWi+/Hfs1vawkZ24+ULyYTOrq3d7tzczoZ171EyaL1dV2uiIiIHESBVU4qh4bVX3/9lSlTpjB58mSy3CXEPTiORu/cQkDH8iVe7e1SCOrftj5KFRERkb8osMpJy+l08txzz2GaJhs2bOCWW25h2bJlBHZtRqOZtxI7dTwxd5+PYfX9a1KyaD2m01VPVYuIiJx8FFjlpGW327nzzju9Cw0UFRXxf//3fzz//POUOhyEjexFUPcWPsc41m4n48ZX2H3+YxQv/A3TNOujdBERkZOKAquc1Nq0acMzzzxD//79vdvmzZvHLbfcwoYNG3z2NU2TnMc/A8C1bR/7bn6dvddMw/HbjjqtWURE5GSjwConvbCwMO666y4mTZpEUFD5TAF79uzhX//6F++++y5ud/n0VoZhEH3XaAK7N/ce61ixmT3jnmbfP6bj3JZRL/WLiIic6BRYRSgPo2eccQbPPPMMbdq0AcDj8fDuu+/y8ssve/cL7NiExOk3E/fkFdhS47zbi+etYffoR8l68ENc+/LqvH4REZETmQKryEGSk5N59NFHGTduHBaLhcDAQM477zyffQzDIPSMriR/dhcxky/AEhNW/oLbQ+EHS9l9zlScW3W1VURExF9s9V2ASENjtVoZN24c3bt3JyMjg5SUFJ/XPR4PFosFw24jfNwAQs/rRf5bC8mfvgCz2EFAuxRszeLrqXoREZETj66wihxG27ZtGThwoM82h8PB7bffztdff+2dIcASEkjUDWeSMucewi8dSNSt52IYhvcY0zQp/vZXTYUlIiJSQwqsItXwzjvvsGXLFqZNm8Z9991HRsaBX/1bY8OJ+df5FabCKl22kX23vUnayP+j8JPlmC53XZctIiJyXFNgFaki0zQpLS31Pl+zZg0333yzz9XWyo7JffErANxp2WTd9x67z3uEwi9+xnR76qRuERGR450Cq0gVGYbBpEmTmDJlCnFx5TMElJSUeK+27tmzp9JjYu4a47O8q2tHJlmTZ7B7zKMUfb0K06PgKiIiciQKrCLV1KNHD55//nmGDx/u3bZmzRomTZrERx99hMvlO1Y1sHNTEl++nsS3biawdyvvdtfWDDLveJv0C56gaN4aBVcREZHDUGAVqYHQ0FBuvvlmn6utZWVlvPXWW0yePBlPJeEzqHsLGr1+E4mv30hgtwOLDzg3p5P5j+mUbUirs/pFRESOJwqsIsegR48evPDCC4wcOdI7M0C/fv2wWA7/Vyuod2sS37qZhJevI6Bz0/Jt/doQ2CG1TmoWERE53mgeVpFjFBISwrXXXsuQIUOYPXs2o0aN8nnd7XZjtVp9thmGQXD/dgT1a0vp4g0HFh/4i+nxsG/SawQNaE/Y+X2wBAXU+vsQERFpqBRYRfykdevW3HbbbRW2T58+nT179nD11VfTqFEjn9cMwyB4YPsKx5Qs/J2SRespWbSevFfmEXHlUMLH9scSElhr9YuIiDRUCqwitWjz5s188cUXeDweVq1axQUXXMD5559PYOCRg6dj9TbvY09WAblPfk7+a98SPm4A4eMGYj3kiqyIiMiJTGNYRWpRfn4+UVFRQPlNWe+++y433XQTP/7442HnbgWI/sdIkj78JyFndIG/xsZ68orJe2kuaWc8QNaDH+Lcvq8u3oKIiEi9U2AVqUXdu3dn2rRpjB492juONSMjg6lTp/LAAw+we/fuwx4b0C6F+CevJPnTuwgd1RNs5X9dTYeTwg+Wsnvk/5E/44c6eR8iIiL1SYFVpJaFhIRw1VVX8eyzz9K5c2fv9pUrVzJp0iTefvttSkpKDnu8vUUicQ9fSspX9xI+YQjG/nGspklQr1aHPU5EROREocAqUkeaNGnCgw8+yJ133umdu9XlcjFr1iyWLl161ONtjaKJueM8Gs+bQtSt5xI6ujcBbZJ99imau5qCD5fiKSmrlfcgIiJSH3TTlUgdMgyDAQMG0LNnTz744AM+/fRTUlNTGTJkSJXbsEQEE3n1sArbTbeH3Ofm4Nq+j9xnvyTsgr6EjxuArVG0H9+BiIhI3VNgFakHQUFBTJgwgdNPP52SkpIK87R+9tln9OjRg8aNG1e5Tccvf+L660YsT14x+W/MJ/+thYScdgrhlw0isFtz7+IGIiIixxMNCRCpR8nJybRs2dJn28aNG3n99deZNGkSL7/8Mvn5+VVqK6hXKxq99w9CR/YE218B2O2heN4a9l7xPHsufpLCz1ZgOpz+fhsiIiK1SoFVpIH5+OOPAfB4PHz55ZdMnDiR999//4g3Zu0X2DGVuKmX0njefUTeeCaW2HDva2Xr08i6913SRv0fptNda/WLiIj4mwKrSANz2223MX78eIKCggAoLi5mxowZTJw4kdmzZ+N0Hv0KqTUugqgbzqLxvPuInXopAR1Tva8F922DYbce4WgREZGGRYFVpIEJDAzkkksu4aWXXuKMM87AYin/a5qXl8crr7zCDTfcwIIFC3C7j36V1LDbCBvZk0bv3kajd24h5KxuhI8f6LOPp7CU9LFPkDd9Ae7colp5TyIiIsdCN12JNFAxMTFMmjSJMWPGMGPGDBYvXgyULzzw9NNPExkZSffu3avUlmEYBHZtRnzXZhVeK/pyJWXr0yhbn0bu83MIPaMLYWP7E9i1uT/fjoiISI0psIo0cCkpKdx5551ccMEFvP3226xatYp27drRrVs3v7Rftin9oCcuimavpGj2Suytkgi5oA8Wu8sv5xEREakpBVaR40TLli35z3/+w9q1awkKCvKZoso0TV599VV69uxJt27dqjV9Vey9FxJx2SAKZy2j8NMVePKKAXBuTifv0U/paLeQs6aEyPGDCOzUxO/vS0RE5GgUWEWOM6ecckqFbWvXrmX27NnMnj2btm3bcskll9C9e/cqB1d7swSi/3kekZPOpnjeGgo/WIpj9TYArE4PJV+sJLBNigKriIjUCwVWkRPAwoULvY//+OMP/vOf/9CmTRvGjRtXreBqCQogbGQvwkb2ouyP3eS9v5iCL1ZgdZmEndvTZ1/n9n2UbdxNyJCOGHb9KBERkdqj/8uInAAmTZpE9+7def/999m+fTtQvgDB/uB64YUX0rt3b++MA1UR0DaZqLvHsKyDjdNSOmCNC/d5veC9xRT87wcs0aGEntuTsDF9CGid5Nf3JSIiAgqsIicEi8XCgAED6N+/P8uWLeO9997zCa5Tp04lNTWViRMn0qVLl2q1bdqtBPb0XY3LdJbfnAXgySmi4J3vKXjne+ztUgg7tychI7phi4/0z5sTEZGTngKryAnEYrFw6qmn0q9fP3788Ufee+89tm3bBsDOnTsJCAjwz4msFuIevZzCT5dT/O2v8NfKWc4NaeRsSCPnqc8J6tOa0JE9CRnWGUtIoH/OKyIiJyUFVpETkMVioX///vTr14+ff/6Zjz76CMMwaN++vc9+27ZtIyYmhoiIiGq1b1gsBPdvS3D/trjziij68heKZv9M2dod5Tt4TEqXbaR02Uas/w0jeED7IzcoIiJyBAqsIicwwzDo1asXvXr1ori42Oc10zR5+umnSU9PZ9iwYYwcOZLk5ORqn8MaGUrE+IFEjB+Ic2sGRV+upGj2z7jSsrHEhhPUt43P/o7fy0NtQIfUak2/JSIiJ68aB1an01mlNc0rMjBNE9M0a3rqE4KJ6f0voP4wTTDNGn6mDth//LG2cyKy2+0+/bJq1Sq2bt0KwJdffsmcOXPo0aMH5557Lp06dcIwjOr3Z+NoQq87nZCJwyhbsx3PvnxcpgecHu8uOc9+iWPZRqyNYwge3oXgM7tga510UoRX/3w+9TN0v/1d4HJpcQt/0M9P/1J/Hl11+sYwq/lTLz8/n8jISGbOnElISEi1CgsLC6NXr95s2bGP0tKyah17oomMCKVVs0as37yL4mJHfZdT74KCAmjZJJ6fflpBYWFhfZdzUigsLGTt2rVs3Lixwv/wY2Ji6NixIy1btsRqtfrtnLYCBx0fW4ZxyE+d0rhgcjslkHtKAqWJoX4734lGP0N96eeGyPGtuLiY8ePHk5eXd9ShaTUOrOnp6cTGxtagPIPtadmU/XWTxskqNCSA5IRItqVlk5mZTURExElxhelwAuxWmqbEAMd21cjpdDJv3jyGDx+O3W73T3EnuIKCAr777jvmzJlDVlaWz2uRkZG0atWKO+64wy/96Sl2UPL1akrmrqHs5y3gqfj9trVMJPiMLoRc0BdrTNgxn7Mh8c/nUz9D97PbrCTGBmO327DZNMLtWOnnp3+pP48uPz+fuLi4KgXWGv8Nt9vtNfoGOF1uDMM4qcMZgIHh/S9w0veJYRhgGNj99D+dmn4+T0YxMTFcdNFFjBkzhmXLlvHZZ5+xceNGAPLy8sjJyfFff0baCbx4AFEXD8CdWUDxt2so+mY1jpV/en+/69qyl4L/ziViTF9sJ+j38Fj6Uz9DD9jfBTabTX/f/Ug/P/1L/Xl41ekX/ZNURIDy/+kPHDiQgQMHsmHDBj7//HOWLl1Kx44dffZzOBwsXLiQQYMGERwcXOPzWePCCb9kAOGXDMCVkUfxvDUUf70Kx+ptBHZphi0p2mf//LcW4M4tJmTYKQR01A1bIiInEwVWEamgXbt2tGvXjoyMDJYtW+bz2qJFi3jxxRd54403GDp0KCNGjKBJkybHdD5bQiQRlw4i4tJBuNJzcOf4jkc0TZP8mYtw784h/7VvsTaKIuS0Uwg5vTOB3Zpj2Pw3zlZERBoeBVYROazo6OgK2+bMmQNASUkJc+bMYc6cOXTo0IGzzz6bfv36HfPiBLak6ApXV13b9+Hek+t97t6TS8HMRRTMXIQlKpTgIR0JOb0zwX3bYATqV28iIicaBVYRqZabbrqJr776iu+//x6Ho3yGi3Xr1rFu3TrCwsIYMmQIp59+Oi1atPDbOe3NEmg8/z8UL/ydku9+peTHjd7VtTy5RRR9uoKiT1dgBAfQ6H+3ENCm+vPJiohIw2Wp7wJE5PjSsmVLJk2axPTp05k4cSKpqane1woLC5k9eza33norX3/9tV/Pa40NJ/yCviRMm0jqDw8R99jlhJzRBSP4oCu6Vgv25ok+xzl+34Hj952at1RE5DimK6wiUiOhoaGce+65nHPOOfz+++/MmzePJUuWUFZWhsVioVevXj77OxwO7HY7Fsux/zvZEhZE6NndCT27O57SMkp/3Ejxd2sxAu0Ydt/xrHnTvqHkh3VYEyIJHtie4MEdCerbBkvwsQ1dEBGRuqPAKiLHxDAMOnXqRKdOnZg4cSKLFi2qdJ7md999l8WLFzN48GCGDBnic2X2WFiCAggZ0omQIZ0qvOYpKaN0+SYA3Bl5FH70I4Uf/YgRaCewVyuCB7Qj+NR22JrGa9YBEZEGTIFVRPwmNDSUs846q8J2l8vF/Pnzyc3N5cMPP+TDDz+kZcuWDBkyhIEDBxITE1M7BZkm0XeOpuT73yldvgnTUb4MoOlwUrp4PaWL15MDWFNiiHt4PEE9WtZOHSIickwUWEWk1uXn59OyZUtWrVqFx+MBYMuWLWzZsoU333yTLl26MGTIEPr27XtMc7seyhISSPjY/oSP7V9+tXXFJkq+/52S79fhzsjz7udOy64wM4FzZyae/GIC2jfG8MMwBhERqTkFVhGpdTExMUyZMoWcnBwWLVrEwoUL2bx5MwAej4dVq1axatUqAgMDeeSRR2jZ0v9XOi3BAYQM7kjI4I6YpolzYzolSzdQumQD7pwibMm+V3kLP1hK/vQFWKJDCerXluBT2xHUpzW2xCi/1yYiIkemwCoidSY6OppRo0YxatQodu3axcKFC1m4cCEZGRkABAYG0rRpU59jsrKyiIiI8OvShoZhENA2mYC2yUT+7TRMt6fCPiWLNwDgySmieM4vFM/5BQBbswSC+rQmqHcrgnq3xhoV6re6RESkcgqsIlIvGjduzGWXXcall17K+vXrWbhwIREREdhsvj+WnnvuOTZu3EifPn0YMGAAXbp08fu63IbV91f+pmkSPm4AJUs2ULp8I2aRw/uaa1sGhdsyKHx/CRgGUf8YSeSVQ/1aj4iI+FJgFZF6ZRgGHTp0oEOHDhVey8/PZ82aNXg8HubPn8/8+fMJDQ2lb9++9O/fny5duhzzylqHq2n/2FfT6caxZhulP26kdMUmHGu3g+uvK7KmSUAL33lfXXtyKZy1jKDerQjo3BRLkKbPEhE5VgqsItJglZaWMmjQIJYvX05JSQkARUVFfPfdd3z33XcEBQXRvXt3+vTp4/cbtvYz7FaCerYkqGdL4Gw8RaU4ftlK6YpNlP60mcAevit6lS77g7yX55L38lywWQk8pQmBPVoS1KMFlo6N/V6fiMjJQIFVRBqshIQE/vGPf1BWVsaqVatYvHgxK1as8IbX0tJSli5dytKlS3nzzTdrJbAeyhIaVL4AwcD2lb6+f95XAFxuHKu24li1lfzXAItBm0ah5K1zEXpqe4IHVbyqLCIiFSmwikiDFxAQQJ8+fejTpw9lZWX88ssvLF++nBUrVlBQUEDr1q0rLFTw+eefU1xcTJ8+fWjWrFmdLQwQddu5BJ3aFsfPWyhd+Seu7fsOvOgxCdldSNHMxXh2ZlUIrO7sQqwxYXVSp4jI8USBVUSOKwEBAfTt25e+ffvidrtZv349LpfLZx/TNPn888/JyMhg5syZxMbG0qNHD3r27Ennzp0JCQmptfpsiVGEjexF2MjypWld+/JwrPwTx8o/Kfl5M67NewAIPGSRAtPpJu3MB7BEhBDYpRmBXZsR2KVZ+TywAfpRLSInN/0UFJHjltVqpVOnikuypqWleafKgvKpsebOncvcuXOx2Wx07NjRG2BTUlJq9eqrLT4S21ndCD2rG06nk29mfcbg+NYEtU3x2a9s427MUifu0jyK562heN6a8hcCbAR2aFweXrs0w9rRP0vaiogcTxRYReSE07hxY15//XWWL1/OypUrWbt2LWVlZUD5MrFr1qxhzZo1vPHGGzzyyCOVzlBQW9whdoIGd6g4NZfbQ1Cf1jjW7sAsPjCNFmUuHKu34Vi97cC26TdDZO1dJRYRaWgUWEXkhBQfH8+5557Lueeei8PhYO3atfz888/8/PPP3quvwcHBtGnTxue4JUuWsGPHDrp27Urr1q0rzAtbWwI7NyXxtRsx3R6cm9NxrNlW/rV6G64dmd79LImREBcOTrd3m336fCy/78LTqhGe1kl4WjXCTIkFq5aUFZETgwKriJzwAgMD6dmzJz179sQ0TdLS0vj5558pLS2tEEi/++47fv75Z959912Cg4M55ZRT6Nq1K127dq314QNQvohBQNsUAtqmED72VKD8ZizHr+Xh1bRbKTvkGMu6XVg3p2PdnA5frwLADArA0yIRT+tGeFol4WmVhJkYCXV085mIiD8psIrIScUwDBo3bkzjxhXnRHW5XPz222/e5yUlJaxYsYIVK1YAEBcXR5cuXejcuTNdu3YlOjq6Tmq2xoQRMqQTIUM64XS5yUnLPvCiaWK43BWOMUrLsK7biXXdTu+2svEDcV3U3+fY8p0VYkWkYVNgFRH5i81m48UXX2T16tXeca55eXne1zMzM72LFkyaNIkzzjjD+5rH48EwjDqbPsvLMCh96m9Q5MCyZQ+Wzel/fe3Bsi/fZ1ezWYLvoTsyCbp3Jp7miXhaJJT/2TxBwwlEpMFRYBUROUh8fDzDhw9n+PDheDwetm/fzurVq1m9ejW///679+atQ2cnWLVqFS+88AKdOnWiY8eOdOrUqU6GEHiFBuLp3BRP56YHtuUWYdm8B+tfIdbdqpHPIZY/92IUlmJdux3r2u3e7WaADU/TeDzNEzFblIdYT+skXYkVkXqjwCoichgWi4XmzZvTvHlzxowZg9PpZMOGDWzcuJGkpCSffdeuXUtWVhbff/8933//PQDR0dF06NCB9u3b065du0qHIdSqqFA8PVvi6dmy0peNMhdmVChGblGF7dZN6Vg3pQNgRoVS8uYk33127IOQQMzYcAVZEal1CqwiIlVkt9s55ZRTOOWUUyq85nK5CAwMxOE4MCVVTk4OS5YsYcmSJUD5ogctWrRgxIgRdVbzkbjO7IrrzK6QU4jlz71YtmaU/7ktA0t6jnc/T/OECscG/PcbrBvSMEMD8TSJx9M0HrNp/F+P4yA0qA7fiYic6BRYRUT84JprruHKK69k8+bN/P777/z222+sW7eOkpIS7z5lZWVYrdYKx77yyis0adKEdu3akZqaWuk+tSo6DE+PMDwHr75V7CgPrn/uxTx0uVjTxLKjfMlZo8iBdf0urOt3+eziiQvHbBKP87xeeDo3q+U3ICInOgVWERE/sdlstGvXjnbt2nHBBRfgdrvZvn07GzZsYP369WzYsIHExESfY/bt28fs2bO9z4OCgmjRogWtW7emdevWtGrViqSkpLq/mSskEE+HVDwdKllZy+nGdWY3LNv3YWzfhyWroMIulswCyCzANbyLz3ZjZyb2d77HbByLp3Fs+Z+psRAcWFvvREROAAqsIiK1xGq10qJFC+8wAKfTyZdffumzz4YNG3yel5aWsm7dOtatW+fdFhYWRqtWrbj99tsJCT3kamd9CLDhnDDkwPPCUiw7MrHs+CvA/vVlFDvwNI33OdSyNQPbT5vhp80+2z2x4b4htnEsno6pGh8rIoACq4hInTr0Smnfvn157LHH2LBhAxs2bGDz5s3elbj2KywsZMOGDYSFheExD2xfs+onMjL2kJySSnJKKlFRMXV/JRYgLAhPh8Z4Ohx0U5lpYmQVYMaE++xqScuqtAlLVgFkFWBds6388OhQSt7wvdHLsnILhsuNJzkGs1EU2AL8+S5EpAFTYBURqUd2u907jGC/vLw8Nm3axKZNm9i8eTObNm0iOTkZq9WK56BFAn5Z+SMb1q31Pg8ODvkrvDYhOaUxScmpJCQ2wmqthx/1hoEZF1Fhs/PiAbiGdcbYmYllVxaWXVkY+/8sLPXu52kcV+FY+6xlWDekAWAaQHwkBU3jcDRPJLB5AvamCdiaxmFLjsGw1fE4YBGpVQqsIiINTGRkpHcpWQDTNH1u3tq/bdeObT7bSkqK2bL5D7Zs/sO7zWazccZZ5zH4tDNrve4qsRiYCZGYCZG+N3mZJuQVe0OsGR5c8dDdB2YuMEwgIw9XRh6un7ZQfNB+UX8fQeS1ww80Xeai9Jc/safGYU2MVJgVOQ4psIqINHCGYRASElJh+w0338nutJ3s3r2T9LSdpKXtpCA/z2cfl8tFaJjvr+VzcrKY9tyjNGqUQmJSMomNkssfJyYREFhPNz8ZRvm8sVGheDo1qfi6aVJ2zTCM3dlYdudg7M7Bkp6NUeSosKvtkHGzzq0ZZFz7379etGBLisHWOAZb41hsqXHYGsdibxyLrXEslkqCsojUPwVWEZHjkGEYxMUnEhefSOeuPb3bCwry2Z1WHmB3797J7rQdJKf43um/N303+Xm55OflsvGP333ajI6JLQ+vjcqD7CldemCzNYD/VRgG7oEdfDbZrQZxFhNbRgGk5eDcvg/X9gwCWif77OfcftCYYJcH185MXDszKz1N6o//h+WgOWTL/kjDnVNUHm4bRenqrEg9aQA/hURExF/CwyNo264jbdt1POw+JSXFBAeHUFJS7LPdNE2yszLJzspk3e9rsNvtdOnWy2ef9b//SkFBPvEJiSQkNKpw9bZOGQaW6BACk2Ox92p92N3sqXGETxiCKy0L164sXDuzMIsrXpm1xIT5hFWAgncXU/jRj+VPbBasiVHYGkVjS47GlhSNNan8sb15IrakaL++PRE5QIFVROQk061HH7p2701BQR570nezd0+a98+9e9IpKysPcwmJyVgsFp9jf1z2vc+NXiGhocTHNyI+oREJCY2IT0gkPiGJmNi4ul8A4TAC2jcmpv2BGQxM08STW1QeXg/6wl7xf4muXQfNauDy4E7Lxp2WjWOl736ho3oS9/ClPttyX/wKS1QotqS/wm1yNJaIkPqZyUHkOKfAKiJyEjIMg4iIKCIiomjT9sCv2j0eD7m52exN313pcfsy9vo8Ly4qYnvRFrZv2+KzfdCQMzhn1IU+7W7843fi4hKIjomtn5kL/mIYBtboMKzRYQSe0vSI+4ad35eADo0PBNvdOXjyiivsZ0uK8XluOpzkvTS34rlDArElRWFNLP+yJUYSNroPtpSYCvuKyAEKrCIi4mWxWIiJiSMmpuK0UgDnX3gpe/emsy9jD/sy9pCRsYf8vNwK+8Un+K7olZebw5uvPu89R1RUDLFx8cTGJRz0ZwKxMXHYAxrO/KqhI7oTOqK7zzZPUSmu9Bzcu3Nw7cnBtTuHoD6+QxJce3Mrbc8sduDcshfnlgPBP3hwR5/AWrzgN3Ke+AxrQiS2xCisiZHecGvd/zwmHMNqqewUIickBVYREamyVm3a06pNe59tjtJS9u3b6xNiUxr73umflXngxiePx0N2dibZ2Zls2ri+wjnuvu9RoqIOjAfdu2c3+fm5REfHERUdU+83gVlCgwholQStkg67jzU+ksQ3bsK1OxtXei7u9GxcfwVcd3oupsN5YN/ESJ9jXWlZuHZk4tqRScWRtvsPshDQJpmkD2732Vz602Y8JQ6ssRFY4yOwxoTpRjE5ISiwiojIMQkMCqJxalMapx7+1+uRUdEMPf1ssjL3kZWZQWZmBo7S0gr7Wa02IiJ8A9xPK5awaOE8oPzX+eERkcTExBITG0dqSiOSk5NJSkoiOTmZhIQE/765GrIEBxDUq1Wlr5mmiSevGPfeXNx787AeshqY6XRjhAdhFlTsHy+3ByoZCpv38lxKl286sMEwsESHlofXuAjvn8H92x62PpGGSIFVRERqXXxCI84aMcb73DRNiosKyczcR3bWPjIzM8jO2ofb7a5wo1duTpbPcfun5Nq2dQu//Hxgvz59+nDPPff4HPvee+8REBBAbGws8fHxxMbGEhMTg91ur503WgWGYWCNCsUaFQptUyq8Hvm304j822l4ikpx783DtTcXd0beIY9zsVdyhde9L993g2niyS7Ek12I848D45KNAJtPYPUUO0gfdj/tAg0yZ23DGhuONSYMa0w4lpiw8sex4ViiQ7E3iccIUHyQuqVPnIiI1DnDMAgNCyc0LJymzVoccd8uXXsRExNPTk4WOdmZ5GRnUVhYUGG/Q6+umqbJrFmzKCsrq3DuqKgoYmNjiYuLIy4ujjPOOINmzZr5HFvfd/NbQoOwtAjC3iLx6Dv/JeKaYbh2ZeHel487swB3Zr73MQct62uN872q687MxywsJagQyrK2HfEcSR/dQUCbA3Pdliz7g8JPVmCNLQ+2lpi/wm50GJbYMKxRoRhhQfXen3J8U2AVEZEG7ZQuPTilSw+fbWVlDgrzc/A4C8nLyyUrK4uOHX3nni0oKKgQVqE8jObk5JCTk8PmzZsB6NWrl09gXb16NY899hjR0dFER0cTExNT4fH+P0NCGs5UVWEje1W63TsMYV95gLW3OCTclzixNomjLCMHa6m70jb2O3QIg3NjOsVf/XLkwmwWAjqkkjTjVp/NhV/8hDurEGtUKJaoECx/XXm2RIViiQjGsOjGMimnwCoiIsedgIBAEhslkxATTHBwUKW/4g8ODuaRRx4hKyuLzMxMMjMzfR7n5ORgmiYAcXG+syJkZ2dTVFREUVERu3btOmwdgYGBfPDBBz7blixZwu7du4mMjCQqKorIyEjv46CgoMO0VLt8hiG0rjiUIKBtMomf3smcOXM4+/QzsBQ4cGcX4MkuxP3Xlye7AHdWIZYo32WC3dkVr3ZX4PJUurlw1jIcv2yt/BiLgSWiPMSGjxtAxPiB3pdMl5vCT5aXvx4R7PtneJCC7glIgVVERE5IdrudDh06HPZ1l8tFTk4OWVlZNGrUyOc1i8VCo0aNyMnJweE47L36REdHV7i6+sMPP7Bs2bJK9w8MDPSG11NPPZUxY8b4vL527VrCwsKIiIggMjKyXmZEMAJs2BoFY2sUVaX9I687g7AL++HJKsSdXfBXuP3rcVYhntwiPHlFlQ5tcOcUHb5hT/kCD57cogork3kKSsh+4MPDvAEDS1gQRkQw1ohgYqdeWj6rw1+cWzMo/WnzgYAbeVDgDQvWdGENlAKriIiclGw2G/Hx8cTHx1d4bejQoQwdOhTTNCkpKSE7O5ucnBzvn/sfR0REVDg2JyfnsOd0OBxkZGSQkZFB+/a+04N5PB7+/e9/4/EcuBoZFhZGeHg44eHhREREeP8866yzSEk5cMNWaWkpBQUFREREEBgYWJPuqDFLSCCWkEBIrXzu3iOJvf9i3Bl55cMVcsqDrSe36K/HxXhyCnHnFWOJCvU5zpNfcvhGTRNPQQkUlOBOo8LV1tJf/iT7wcOEXcAIC8ISEYw9NY7E1270ea3o61W4dmVhCQ/GCA3CEh6EJSwIS1hweUgOK3+u0Ot/CqwiIiKHYRgGISEhhISE0Lhx46MfAFx33XXs3buX/Px8cnNzycvL837tf15QUEBUVJTPcYWFhT5hdf+2wsJC0tPTfbb37dvXJ7CuXr2aqVOnAhAQEOANuYeG3aioKEaOHFnhHEC9zJwQ1P3IN9ztt3/oxn6WqFBi7r8YT34xnvySw/9ZUIIlItjnWE9+xZXKfM5VWIq7sLTSmRCKPv+JkkUV5w4+VPj4gYT/07ef9931DpYgO8Zf4Xb/l/FX4DdCA7GEBmJLjSv/B4D4UGAVERHxo5YtW9KyZcsj7uN2uyuEU4vFwtixY31Cbm5uLgUFBRQV+f7qPDzc98angoID40jLysrIysoiKyuLQ4WHh1cIrK+++ioLFizAMAzsdjtffPEFoaGhhIWFERoa6n3crl07Tj31VJ9j09LSCAwMJCQkhKCgoApTkvnLocMurJEhhF/Q96jHHRp0AYIHtscaFYonvxh33iEht+DAc2t0WIVjPUVHmBv34HoDfcO/6XRTPOcoN6b9JeGl6wg+tZ33ecnyTWTdMxMjJABLaFD5nyFBWEIDK2wLnzD4/9u7/5i2yjUO4N/S9vDDbXSFuwWuboZcnEbvEq7i1gJlwLKAiHf8Ac4hmTEyHQvLoiQuWW7cjNH4hxqjJvPuD3UhbtmyqDExG5vJnAzNzCJjsk3HsqCObGiA/uC0PT09z/2jt0dOe0pbBrTA80ma0vd9T3vOQw88ffue99XESr49DhKlUDKckwlDtjBve385YWWMMcbmmNFohNGoXYFqyZIlePrpp3XbB4NBuN1u9bZypXY8aH5+Pmw2m6aNy+WCLMuadnfdpf1qHYCaDBMRJEnCyMhIVBsAEEUxKmHt7OxUtzcYDMjOzkZ2djZycnLU+5ycHDQ2NmLNmjXqduPj47h48aKm3eT2M9Xbqzd7g/CPAs2Y1mQs7/x3aAiD2wdlwgfF7QN5vFA8vr9ubi9Mq7TDIxJNdAHAcJe2d1VxTiAYY6lfzXaZZizbtkFT5vzvKXiO9ka1M+QIMGQLyMjJhDF/GVYe3JHw/qUKJ6yMMcZYmjMajbBYLFHDCMJKSkpQUlKiKQuPv52cxOopKipCIBCAx+PB7du3AYSS2GBQO71VZLKrKApE8a+v14kIoihCFMWo3t2amhrN46GhIbz11lsxj9dkMqnJ6wcffKBJYHt6etDf34+srKwpbxaLBYWFhTFfYzoy/xl7NbdIgcBfy+9mLM3G30/+J5ToerxQJnwgtw/KhB+K6AeJfigToXvTSovmeQwZGTCuyA1tI0qATq8xABhyhKiyyIvVAID8gdDSwGMTCCLOeOA0wgkrY4wxtgBNHn8b2SM72datWwGEEqyvvvoKjz32GEwmE/x+vzq1l8fjiUqWZVlGTU0NPB4PvF4vvF6vmrCGH4dlZ2vHkU6u0yPLMlwuFzweT9RMCZcvX8aJEyfiHn9JSQn279+vKXvxxRcxMjKim+AKgoDMzEwIgoCKigo89NBDmv39/vvvNW0m34dvgiDAbDZH9ewajBkwFVrj7rOenI1rkbNxLYDQhwLySiBRgiKGElhlwgea8OsOf8j8VxFgMISS4QkfFK/0/+396s+GeTJelhNWxhhjjGkYDAY1kcvLy9NtIwgCdu3aFfM5FEWBz+eDKIpRsymsXr0abW1tmgR38n345/C+TObzJfb1ut6ct06nEy6XCy6XS2cL7f5NTljHxsbwzjvvJPS6Bw4c0PTsfvvtt/jiiy80Ca7ZbIYgCBAEASaTCYIgwGq1YvPmzZrn6u/vh9PphNlsVrcJ/2w2myFYBJj/Zgl9MIn4UAAAS5vsWNpkn3J/9RLddMQJK2OMMcZmXEZGhtrDG6mgoCDq4q9EtbS0oL6+Hn6/H16vN+b96tXRX9+vWLECRqNRbRcr+Y2cGmyquXjjbTs2NoYbN2IsjjDJPffcE5WwHj9+HD/++GPcbRsaGtDW1qY+JiK0trbCZDLpJrpms1mta25uRlFRYrM1pBInrIwxxhibN/Ly8mL2+sbzxhtvaB4riqImrpIkQZIk+P3+qLl5rVYrnn/+efj9frVN+D58Cz+O7NlVFAUmkynqArhIeotETB4HO5XIi9TCQyoSUVdXl1C7VOOElTHGGGOLUkZGhjqzwVRyc3NRX18/rdfYvHkzmpqaEAwG1aQ2EAggEAhAkiT1Z72ZERoaGrBu3bqotpGPI3uTg8EgCgoKotrrJc2pWE1tOubHXjLGGGOMzWNGozGh5Hgym802rdfKysrChx9+GFWuKApkWYYsy2oiq7daWzpKOmEND84dHR2d5ksa4PNOQAoE4zddwDIgweUCfD4RsuyH3ysC0dPFLRqKbITLaQRwZ4O/ZVmGwWDA6OjovPnUmM44njNrZuLJf0PDgiYjXK4AvF6R358zIBAIqFNSpWLVq4VmPsTTYDDEnO5sLoRfO5ELv5I+w8NPft999yW7KWOMMcYYYxputxu5ublTtjFQkvMZKIqC4eFhLF26VHcFCcYYY4wxxuIhIrjdbhQWFsZd1jfphJUxxhhjjLG5NHU6yxhjjDHGWIpxwsoYY4wxxtJaUgnryy+/jIqKCrS2tiY8me1idP78edhsNjgcDjz11FOaWAWDQTz77LOoqKjA7t271fJ3330XZWVleOKJJ9TJfnt6emC321FeXo5Lly7N9WGkncOHD0dN5szxnJ4zZ86gpqYGVVVV+OyzzzR1euf5sWPHYLfbUVNTg99//x0AcPXqVTgcDtjtdnz99ddzfgzpRFEUPPPMM6ioqEB5eTmuXr2qqeeYxud0OvHoo49iyZIl+OmnnwDox2iyRM/zW7duYdOmTSgrK0NXV9fcHVQKRcbT7XajuroaDocD1dXVGBoaitqG4xmb3vsTAIaGhpCZmakpC+N4zjBKUF9fH7W0tBAR0WuvvUaffvppopsuOsPDwySKIhER7dmzh44dO6bWff7557R3714iInruueeot7eX/vjjD6quriZFUairq4tef/11IiJyOBw0OjpKQ0NDVFdXN/cHkkZkWabGxkYqKSnRlHM8kyeKIj3++OPk9/uj6vTO80AgQOvXrye/3089PT20fft2IiJqbGykX375hZxOJ9nt9jk9hnRz4cIF2rJlCxERnT17ltra2tQ6jmliJEmikZER2rZtG126dClmjMKSOc93795Np06dUp/T6/XO+fHNtch4er1eunnzJhERnThxgnbu3Klpz/GcWmQ8w9rb26mqqkpTRsTxnA0J97D29vZi06ZNAIDa2lqcO3du1pLo+a6goECdGFgQBM2Vb3px/OGHH1BZWQmDwaCWeb1eGI1GLF++HKtWrbqDeW8XhsOHD6OpqSnqKkKOZ/K+++47ZGdno6GhAY2Njbh165ZapxfPa9eu4YEHHoAgCCgrK0N/fz8AYHh4GMXFxVi2bBmsViv+/PPPlBxPOrj77rtBRCAijI2NIT8/X63jmCbGbDZrvkGJFaOwZM7z8+fPo7q6GiaTCY888ohub9hCExnPrKwsFBYWAoj+vwRwPOOJjCcA3LhxAwaDAatWrYpqz/GceQknrGNjY+pqCLm5uYv6H36ihoaG0N3djYaGBrVML47xyoDQ0mmSJM3tAaSJYDCIo0eP4sknn4yq43gm7/bt2xgcHMSXX36JtrY27Nu3T61LJHbBYGjCekVR1LLF/jchPz8fZrMZ999/Pzo6OtDe3q7WcUynJ1aM9OrjneeBQEBN0BZ7XCVJwr59+9DR0aEp53gm780330RnZ6duHcdz5iWcsFosFnUMhtPphNVqnbWdWghcLhdaW1vx8ccfa1a40ItjvDIgtEKOIAhzexBpoqurC83NzbpztHE8k2exWFBWVgZBEFBTU4OBgQFNXbzYGY1GAND8Phb734Tu7m6YTCb8/PPPOH78OF566SW1jmM6PbFipFcf7zw3m83qh4HFHtft27ejvb0dxcXFmnKOZ3KuX78OALj33nt16zmeMy/hhNVut+P06dMAgJMnT6KsrGzWdmq+k2UZW7ZswSuvvII1a9YAAG7evAki0o1jaWkpzp49qynLycmBLMsYHx/Hb7/9tqjfwJcvX8ahQ4dQW1uLa9euYdeuXRzPO1BaWoorV66AiNDX14eioiKMjIxAkiTdeBYXF+PKlSuQJAm9vb1Yu3YtgNDQl+vXr8PtdmN0dFTzNfhiQ0TIy8sDEOptdTqdHNM7FCtGo6OjEEUxqfO8tLQUZ86cgSzLuHDhAh588MGUHVcq7d+/H0VFRZpvqzie03Px4kUMDAygtrYWp06dwgsvvACfz8fxnE3JDHjt7Oyk8vJy2rp1q+4FGyzk0KFDZLVaqbKykiorK+nIkSO0YcMG8vl8FAgEaNu2bVReXk4dHR3qNm+//TbZ7Xaqr6+n8fFxIiL65ptvyGazkd1up76+vlQdTlp5+OGHiYg4nnfo/fffp4qKCnI4HDQ4OEgtLS00ODhIRPrn+ZEjR8hms1FVVRX9+uuvREQ0MDBA5eXlZLPZqLu7O2XHkg4CgQA1NzeTw+GgdevW0blz5zim01BXV0cFBQW0fv16+uijj3RjtHfvXjp9+jQRJX6eDw8P08aNG8lms9Enn3ySmoNLgcnxfPXVV8loNKr/l/bs2UNEHM9kRL4/wyZfiMXxnD280tUcCAaD2LlzJw4cOJDqXVkQOJ4zr62tDQcPHkz1biwoHNPZsWPHDrz33nswmUyp3pUFgeM5szies4cTVsYYY4wxltZ4pSvGGGOMMZbWOGFljDHGGGNpjRNWxhhjjDGW1jhhZYwxxhhjaY0TVsYYY4wxltY4YWWMMcYYY2mNE1bGGGOMMZbWOGFljDHGGGNpjRNWxhhjjDGW1jhhZYwxxhhjae1/IhO/mvWV+mAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import expon\n",
    "\n",
    "\n",
    "# Bind reference to best fitting exponential distribution.\n",
    "dist_u = expon(scale=scale_u)\n",
    "dist_c = expon(scale=scale_c)\n",
    "\n",
    "\n",
    "xx = np.linspace(0, 1.5 * max(vals), 1000)\n",
    "yy_u = dist_u.pdf(xx)\n",
    "yy_c = dist_c.pdf(xx)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(7, 5), tight_layout=True) \n",
    "\n",
    "ax.set_title(\"uncensored vs. censored MLE estimates\", fontsize=9, weight=\"normal\")\n",
    "ax.hist(\n",
    "    vals, 6, density=True, alpha=.80, color=\"#ccd3e1\",\n",
    "    edgecolor=\"#FFFFFF\", linewidth=1.0\n",
    "    )\n",
    "\n",
    "ax.plot(xx, yy_u, linewidth=2.0, linestyle=\"--\", color=\"#4d4d4d\", label=\"uncensored\")\n",
    "ax.plot(xx, yy_c, linewidth=2.0, linestyle=\"--\", color=\"#E02C70\", label=\"censored\")\n",
    "ax.tick_params(axis=\"x\", which=\"major\", direction='in', labelsize=6)\n",
    "ax.tick_params(axis=\"y\", which=\"major\", direction='in', labelsize=6)\n",
    "ax.get_xaxis().set_major_formatter(mpl.ticker.FuncFormatter(lambda x, p: format(int(x), ',')))\n",
    "ax.set_yticklabels([])\n",
    "ax.set_xlim(0)\n",
    "ax.xaxis.set_ticks_position(\"none\")\n",
    "ax.yaxis.set_ticks_position(\"none\")\n",
    "ax.legend(loc=\"best\", fancybox=True, framealpha=1, fontsize=\"small\")\n",
    "ax.annotate(\n",
    "    r\"$\\hat \\theta_{u} = $\" + f\"{scale_u:,.0f}\", xy=(.85, .80), \n",
    "    xycoords=\"axes fraction\", ha=\"left\", va=\"bottom\", fontsize=12, \n",
    "    rotation=0, weight=\"normal\", color=\"#000000\"\n",
    "    ) \n",
    "ax.annotate(\n",
    "    r\"$\\hat \\theta_{c} = $\" + f\"{scale_c:,.0f}\", xy=(.85, .75), \n",
    "    xycoords=\"axes fraction\", ha=\"left\", va=\"bottom\", fontsize=12, \n",
    "    rotation=0, weight=\"normal\", color=\"#000000\"\n",
    "    )\n",
    "ax.grid(True)   \n",
    "ax.set_axisbelow(True) \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>q</th>\n",
       "      <th>uncensored</th>\n",
       "      <th>censored</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.00</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.05</td>\n",
       "      <td>163.0</td>\n",
       "      <td>210.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.10</td>\n",
       "      <td>335.0</td>\n",
       "      <td>431.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.15</td>\n",
       "      <td>517.0</td>\n",
       "      <td>664.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.20</td>\n",
       "      <td>709.0</td>\n",
       "      <td>912.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.25</td>\n",
       "      <td>914.0</td>\n",
       "      <td>1176.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.30</td>\n",
       "      <td>1134.0</td>\n",
       "      <td>1458.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.35</td>\n",
       "      <td>1369.0</td>\n",
       "      <td>1760.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.40</td>\n",
       "      <td>1624.0</td>\n",
       "      <td>2087.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.45</td>\n",
       "      <td>1900.0</td>\n",
       "      <td>2443.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.50</td>\n",
       "      <td>2203.0</td>\n",
       "      <td>2832.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0.55</td>\n",
       "      <td>2538.0</td>\n",
       "      <td>3263.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>0.60</td>\n",
       "      <td>2912.0</td>\n",
       "      <td>3744.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.65</td>\n",
       "      <td>3337.0</td>\n",
       "      <td>4290.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>0.70</td>\n",
       "      <td>3827.0</td>\n",
       "      <td>4920.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>0.75</td>\n",
       "      <td>4406.0</td>\n",
       "      <td>5665.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>0.80</td>\n",
       "      <td>5115.0</td>\n",
       "      <td>6577.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.85</td>\n",
       "      <td>6030.0</td>\n",
       "      <td>7752.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>0.90</td>\n",
       "      <td>7318.0</td>\n",
       "      <td>9409.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>0.91</td>\n",
       "      <td>7653.0</td>\n",
       "      <td>9840.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>0.92</td>\n",
       "      <td>8028.0</td>\n",
       "      <td>10321.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>0.93</td>\n",
       "      <td>8452.0</td>\n",
       "      <td>10867.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>0.94</td>\n",
       "      <td>8942.0</td>\n",
       "      <td>11497.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>0.95</td>\n",
       "      <td>9521.0</td>\n",
       "      <td>12242.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>0.96</td>\n",
       "      <td>10231.0</td>\n",
       "      <td>13154.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>0.97</td>\n",
       "      <td>11145.0</td>\n",
       "      <td>14329.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>0.98</td>\n",
       "      <td>12434.0</td>\n",
       "      <td>15986.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>0.99</td>\n",
       "      <td>14637.0</td>\n",
       "      <td>18819.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       q  uncensored  censored\n",
       "0   0.00         0.0       0.0\n",
       "1   0.05       163.0     210.0\n",
       "2   0.10       335.0     431.0\n",
       "3   0.15       517.0     664.0\n",
       "4   0.20       709.0     912.0\n",
       "5   0.25       914.0    1176.0\n",
       "6   0.30      1134.0    1458.0\n",
       "7   0.35      1369.0    1760.0\n",
       "8   0.40      1624.0    2087.0\n",
       "9   0.45      1900.0    2443.0\n",
       "10  0.50      2203.0    2832.0\n",
       "11  0.55      2538.0    3263.0\n",
       "12  0.60      2912.0    3744.0\n",
       "13  0.65      3337.0    4290.0\n",
       "14  0.70      3827.0    4920.0\n",
       "15  0.75      4406.0    5665.0\n",
       "16  0.80      5115.0    6577.0\n",
       "17  0.85      6030.0    7752.0\n",
       "18  0.90      7318.0    9409.0\n",
       "19  0.91      7653.0    9840.0\n",
       "20  0.92      8028.0   10321.0\n",
       "21  0.93      8452.0   10867.0\n",
       "22  0.94      8942.0   11497.0\n",
       "23  0.95      9521.0   12242.0\n",
       "24  0.96     10231.0   13154.0\n",
       "25  0.97     11145.0   14329.0\n",
       "26  0.98     12434.0   15986.0\n",
       "27  0.99     14637.0   18819.0"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# Quantiles.\n",
    "import pandas as pd\n",
    "\n",
    "q0 = np.arange(0, .91, .05)\n",
    "q1 = np.arange(.91, 1.0, .01)\n",
    "q = np.append(q0, q1)\n",
    "\n",
    "dfq = pd.DataFrame({\n",
    "    \"q\": q, \n",
    "    \"uncensored\": np.round(dist_u.ppf(q), 0), \n",
    "    \"censored\": np.round(dist_c.ppf(q), 0), \n",
    "    })\n",
    "\n",
    "    \n",
    "dfq.head(30)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Takeaways:\n",
    "\n",
    "\n",
    "- We should account for censoring for AGs/loss causes in which limits/TIVs are readily available. When limits/TIVs are not readily available, just treat losses as uncensored. \n",
    "\n",
    "- Not accounting for censored observations results in underestimation of the tail of the distribution. \n",
    "\n",
    "- If the number of losses is large, and the number of losses at the limit are relatively small, censoring will not change the model parameterization significantly.\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py312",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}