cmkl_assignments/fall-2024/math/mat-203/00010/01_ExampleDescriptiveStatsStudents.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'numpy'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[1], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m#Numerical Python\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# To plot\u001b[39;00m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyplot\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mplt\u001b[39;00m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'numpy'"
     ]
    }
   ],
   "source": [
    "#Numerical Python\n",
    "import numpy as np\n",
    "\n",
    "# To plot\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "# Gamma function\n",
    "from scipy.special import gamma\n",
    "\n",
    "# To calculate statistics\n",
    "from statistics import mode\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import hmean, trim_mean, iqr, median_abs_deviation, skew, kurtosis\n",
    "from scipy.stats.mstats import gmean, winsorize\n",
    "from statsmodels.distributions import ECDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ALPHA = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = plt.rcParams['axes.prop_cycle'].by_key()['color']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Dat = np.array([0.0, 0.8, 1.0, 1.2, 1.3, 1.3, 1.4, 1.8, 2.4, 4.6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Empirical cumulative distribution function (ECDF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ecdf = ECDF(Dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.05, 4.6499999999999995)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1EAAAINCAYAAADfvvWSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABHhElEQVR4nO3dfVxUdf7//+eAMigKYlx5QYJGmdeCStKapii0ZrldYG6lkrVWympollbatWWuq1um1me9aPuWqK1t24XGQmqprRdoF5aXYXYFyKqDgoLOnN8f/ZwcuZCDwAzwuN9uc7s573mfc17ncGY8z3mfc8ZiGIYhAAAAAECleLm7AAAAAACoSwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAA9cayZctksVicD19fX7Vu3VoJCQn629/+phMnTlRpvps3b9aTTz6p48ePV2/BAIA6iRAFAKh3nn76af3jH//QwoULlZKSIkmaNGmSunbtqi+//NL0/DZv3qynnnqKEAUAkCQ1cncBAABUtxtuuEG9evVyPp82bZoyMzN144036qabbtK3336rJk2auLFCAEBdxkgUAKBBGDhwoJ544gl9//33evPNNyVJX375pcaMGaP27dvL19dXYWFhuueee/S///3POd2TTz6phx9+WJIUGRnpPFXw0KFDkqSlS5dq4MCBCgkJkdVqVadOnbRw4cJaXz8AQO1hJAoA0GDcfffdmj59uj7++GPdd999Sk9P13fffafk5GSFhYVp9+7deu2117R79259/vnnslgsuuWWW7Rv3z69/fbb+utf/6qgoCBJUnBwsCRp4cKF6ty5s2666SY1atRI//73v/Xggw/K4XBo/Pjx7lxdAEANsRiGYbi7CAAAqsOyZcuUnJysbdu2uZzOd74WLVqoffv2ysrK0qlTp0qd1rdixQqNHDlSGzduVL9+/SRJc+bM0cMPP6zs7GxFRES49C9rHomJidq/f78OHjxYfSsHAPAYnM4HAGhQmjVr5rxL3/nh5/Tp08rPz9c111wjScrKyqrU/M6fh81mU35+vvr376/vvvtONputGisHAHgKQhQAoEE5efKkmjdvLkk6evSoJk6cqNDQUDVp0kTBwcGKjIyUpEoHoE2bNik+Pl5+fn5q0aKFgoODNX36dFPzAADULVwTBQBoMH788UfZbDZdccUVkqSkpCRt3rxZDz/8sHr06KFmzZrJ4XAoMTFRDofjovM7ePCgBg0apI4dO2ru3LkKDw+Xj4+PPvzwQ/31r3+t1DwAAHUPIQoA0GD84x//kCQlJCTo2LFjysjI0FNPPaUZM2Y4++zfv7/UdBaLpcz5/fvf/1ZxcbHee+89XX755c72Tz75pJorBwB4Ek7nAwA0CJmZmXrmmWcUGRmpO++8U97e3pKkC++vNG/evFLT+vn5SVKpH9stax42m01Lly6txsoBAJ6GkSgAQL3z0Ucfac+ePTp79qxyc3OVmZmp9PR0tWvXTu+99558fX3l6+ur6667TrNnz9aZM2fUpk0bffzxx8rOzi41v5iYGEnSY489pjvuuEONGzfWsGHDNGTIEPn4+GjYsGEaN26cTp48qddff10hISH65Zdfanu1AQC1hBAFAKh3zp2e5+Pjo5YtW6pr166aN2+ekpOTnTeVkKS33npLKSkpWrBggQzD0JAhQ/TRRx+pdevWLvPr3bu3nnnmGS1atEhr166Vw+FQdna2rrrqKq1evVqPP/64pkyZorCwMD3wwAMKDg7WPffcU6vrDACoPfxOFAAAAACYwDVRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACY0MjdBXg0w5CKin79d9OmksXi3nrgXp6wP9SHGhwOKT//138HBUleVfguxxO2gydgO3gO/hYA4Nmq+XOakaiKFBVJzZr9+ji30dFwecL+UB9qyM+XQkN/fZwLU7VdQ33BdvAc/C0AwLNV8+c0IQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACY0cncB7mYYhk6cOFH2i4WFv/27oECy22unKHgmT9gf6kMN57/fTpyQfH1rv4b6gu3gOfhbAIBnM/k53bx5c1kslnJftxiGYVRXbXVRQUGBAgIC3F0GAAAAAA9hs9nk7+9f7usNPkRVOBKlX0NWeHi4fvjhhwo3JFBb2Cfhadgn4WnYJ+FJ2B/rpouNRDX40/ksFkuldmh/f392fHgU9kl4GvZJeBr2SXgS9sf6hRtLAAAAAIAJhCgAAAAAMIEQdRFWq1UzZ86U1Wp1dymAJPZJeB72SXga9kl4EvbH+qnB31gCAAAAAMxgJAoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhSAem39+vWyWCxav359pfoPGDBAAwYMqLF6IiIiNGbMmBqbv6cwu90ry2Kx6Mknn6zWedaEQ4cOyWKxaNmyZTW+rGXLlslisejQoUPOtoiICN144401vmyp5v7W9dWTTz4pi8Xi7jIAXCJCFIBqc+5grrzH559/7u4S65TTp0/rr3/9q2JjYxUQECBfX19deeWVmjBhgvbt2+fu8mrMhx9+6HFB6fz9uFGjRmrZsqViYmI0ceJEffPNN9W2nFdffbVWgldVeGJtAwYMKPfzZs+ePW6rq6ioSE8++STBEqjHLIZhGO4uAkD9sGzZMiUnJ+vpp59WZGRkqdcTExMVFBRUqzU5HA6VlJTIx8dHXl4X/96opKREkuTj41Mj9URERGjAgAEXPRjNz89XYmKiduzYoRtvvFHx8fFq1qyZ9u7dqxUrVignJ8dZqydav369rr/+en3yySemR/YmTJigBQsWqKz/nk6fPq1GjRqpUaNG1VRp5VgsFg0ePFijRo2SYRiy2Wz64osvtGrVKhUWFurFF19Uamqqs79hGCouLlbjxo3l7e1d6eV06dJFQUFBpg6+7Xa7zpw5I6vV6hzhiIiIUJcuXfT+++9Xej5Vrc3se6w6DRgwQAcPHtSsWbNKvXbTTTfJ39+/Vus5Jz8/X8HBwZo5c2apLwTOnj2rs2fPytfX1y21Aagetfu/EIAG4YYbblCvXr3cXYYkycvLq1IHK0VFRWratGmNhSezxowZo507d2r16tW69dZbXV575pln9Nhjj7mpMvdy54HnlVdeqbvuusul7YUXXtCwYcM0efJkdezYUb///e8l/Rq6arrWwsJC+fn5ydvb21RQq26VfY/VlICAgFJ/F0/mji8BAFQ/TucDUOvOXS8yZ84cLViwQO3bt1fTpk01ZMgQ/fDDDzIMQ88884zatm2rJk2a6Oabb9bRo0dd5nHumo+PP/5YPXr0kK+vrzp16qR//vOfLv3Kul5jwIAB6tKli3bs2KHrrrtOTZs21fTp052vXThycvr0aT355JO68sor5evrq1atWumWW27RwYMHnX3mzJmjuLg4XXbZZWrSpIliYmK0evXqKm2f//73v/rggw80duzYUgFKkqxWq+bMmeOyPmWN9owZM0YRERHO59Wx3cu7Jqky13p9+umnuv3223X55ZfLarUqPDxcDz30kE6dOuVS84IFC5zLOvcoa/mrV6+WxWLRhg0bSi1r8eLFslgs+vrrr51te/bs0W233aaWLVvK19dXvXr10nvvvVdhzRdz2WWXacWKFWrUqJGee+45Z3tZ10Tl5OQoOTlZbdu2ldVqVatWrXTzzTc7r2WKiIjQ7t27tWHDBud6n/u7njtVdsOGDXrwwQcVEhKitm3burx2/jVR51zs/VHe9TkXzrOi2sq7JmrVqlWKiYlRkyZNFBQUpLvuuks//fSTS58xY8aoWbNm+umnnzR8+HA1a9ZMwcHBmjJliux2+0W2/sWVt20q+lz45ptvdP3116tp06Zq06aNZs+eXWq+FX0mHDp0SMHBwZKkp556yrm9zu23ZW3zs2fP6plnnlGHDh1ktVoVERGh6dOnq7i42KXfuc+9zz77TH369JGvr6/at2+vN95445K3FQBz+CoEQLWz2WzKz893abNYLLrssstc2v7f//t/KikpUUpKio4eParZs2crKSlJAwcO1Pr16/XII4/owIEDevnllzVlyhQtWbLEZfr9+/drxIgRuv/++zV69GgtXbpUt99+u9auXavBgwdXWOP//vc/3XDDDbrjjjt01113KTQ0tMx+drtdN954ozIyMnTHHXdo4sSJOnHihNLT0/X111+rQ4cOkqT58+frpptu0p133qmSkhKtWLFCt99+u95//30NHTrU1PY7d2B/9913m5qusi51u1fVqlWrVFRUpAceeECXXXaZtm7dqpdfflk//vijVq1aJUkaN26cfv75Z6Wnp+sf//hHhfMbOnSomjVrppUrV6p///4ur6Wlpalz587q0qWLJGn37t269tpr1aZNGz366KPy8/PTypUrNXz4cL3zzjv6wx/+UOX1uvzyy9W/f3998sknKigoKPcUsltvvVW7d+9WSkqKIiIilJeXp/T0dB0+fFgRERGaN2+eUlJS1KxZM+dI44X75YMPPqjg4GDNmDFDhYWFFdZ1Ke+PC1WmtvOdO7W3d+/emjVrlnJzczV//nxt2rRJO3fuVIsWLZx97Xa7EhISFBsbqzlz5ug///mP/vKXv6hDhw564IEHLlqb3W4v9Xnj6+urZs2amVpHSTp27JgSExN1yy23KCkpSatXr9Yjjzyirl276oYbbnAur6LPhPj4eC1cuFAPPPCA/vCHP+iWW26RJHXr1q3c5d57771avny5brvtNk2ePFn//e9/NWvWLH377bdas2aNS98DBw7otttu09ixYzV69GgtWbJEY8aMUUxMjDp37mx6nQFUkQEA1WTp0qWGpDIfVqvV2S87O9uQZAQHBxvHjx93tk+bNs2QZHTv3t04c+aMs33kyJGGj4+Pcfr0aWdbu3btDEnGO++842yz2WxGq1atjJ49ezrbPvnkE0OS8cknnzjb+vfvb0gyFi1aVGod+vfvb/Tv39/5fMmSJYYkY+7cuaX6OhwO57+LiopcXispKTG6dOliDBw40KW9Xbt2xujRo0vN63x/+MMfDEnGsWPHKuxXXs3njB492mjXrp3zeXVsd0nGzJkzSy3rwvUqa7tfuI0MwzBmzZplWCwW4/vvv3e2jR8/3ijvv6cLlz9y5EgjJCTEOHv2rLPtl19+Mby8vIynn37a2TZo0CCja9euLuvicDiMuLg4IyoqqsxlXbjc8ePHl/v6xIkTDUnGF198YRjGb9t66dKlhmEYxrFjxwxJxksvvVThcjp37lzm3/Lce+t3v/udy7qe/1p2drazrbLvj5kzZ5a5rcuaZ3m1Xfi3LikpMUJCQowuXboYp06dcvZ7//33DUnGjBkznG2jR482JLn8rQzDMHr27GnExMSUWtaFzr2XL3yc2xfLWo+yaj5/Xm+88Yazrbi42AgLCzNuvfVWZ1tlPhOOHDlS7nvlwm2+a9cuQ5Jx7733uvSbMmWKIcnIzMx0tp37u27cuNHZlpeXZ1itVmPy5MnlbygA1Y7T+QBUuwULFig9Pd3l8dFHH5Xqd/vttysgIMD5PDY2VpJ01113uVwzEBsbq5KSklKnArVu3dplBMHf31+jRo3Szp07lZOTU2GNVqtVycnJF12Xd955R0FBQUpJSSn12vmn5DRp0sT572PHjslms6lfv37Kysq66DIuVFBQIElq3ry56Wkr41K3e1Wdv40KCwuVn5+vuLg4GYahnTt3VmmeI0aMUF5enstpWatXr5bD4dCIESMkSUePHlVmZqaSkpJ04sQJ5efnKz8/X//73/+UkJCg/fv3X/I6nhv1OHHiRJmvN2nSRD4+Plq/fr2OHTtW5eXcd999lb7+6VLeH5di+/btysvL04MPPuhyrdTQoUPVsWNHffDBB6Wmuf/++12e9+vXT999912llhcREVHq82bq1KlVqr1Zs2Yu11f5+PioT58+LrVU9jOhsj788ENJcrkxiSRNnjxZkkptr06dOqlfv37O58HBwbrqqqsqvb0AVA9O5wNQ7fr06VOpG0tcfvnlLs/PHdiHh4eX2X7hwecVV1xR6qDlyiuvlPTrNSlhYWHlLrtNmzaVuonEwYMHddVVV130QvD3339fzz77rHbt2uVyHUNVDqrOnQ524sQJl9OeqsulbveqOnz4sGbMmKH33nuv1DxtNluV5pmYmKiAgAClpaVp0KBBkn49la9Hjx7OfeHAgQMyDENPPPGEnnjiiTLnk5eXpzZt2lSpBkk6efKkpPKDr9Vq1YsvvqjJkycrNDRU11xzjW688UaNGjWqwv30QmXd9bI8l/L+uBTff/+9JOmqq64q9VrHjh312WefubT5+vo6ryE6JzAwsNL7nZ+fn+Lj46tYrau2bduW2maBgYH68ssvnc8r+5lQWd9//728vLx0xRVXuLSHhYWpRYsWzu15zoXv33M1Vtf7FEDlEKIAuE1536iX125U4y8ynD8qcqk+/fRT3XTTTbruuuv06quvqlWrVmrcuLGWLl2qt956y/T8OnbsKEn66quvXL5xLo/FYilz25R3YX5NbPeL3QTAbrdr8ODBOnr0qB555BF17NhRfn5++umnnzRmzBg5HI6LLqMsVqtVw4cP15o1a/Tqq68qNzdXmzZt0vPPP+/sc27eU6ZMUUJCQpnzufAA1qyvv/5a3t7eFYacSZMmadiwYXr33Xe1bt06PfHEE5o1a5YyMzPVs2fPSi2nOvdbqfyQXx03daismryzoNn1q43PnvJU9gsXd9YI4DeEKAB11rkRhvMPPs79CO35d6W7FB06dNB///tfnTlzRo0bNy6zzzvvvCNfX1+tW7dOVqvV2b506dIqLXPYsGGaNWuW3nzzzUqFqMDAwDJP5bnwG+zqEBgYqOPHj7u0lZSU6Jdffqlwuq+++kr79u3T8uXLNWrUKGd7enp6qb5mR+9GjBih5cuXKyMjQ99++60Mw3CeyidJ7du3lyQ1bty42kYsznf48GFt2LBBffv2vegpmB06dNDkyZM1efJk7d+/Xz169NBf/vIXvfnmm5KqNnJZnsq8PwIDAyVJx48fdxn1LGvfqWxt7dq1kyTt3btXAwcOdHlt7969ztdrw/nrd75LeW9U5jPBzN+xXbt2cjgc2r9/v66++mpne25uro4fP16r2wtA5XFNFIA66+eff3a5c1VBQYHeeOMN9ejRo9pOVbr11luVn5+vV155pdRr57759fb2lsVicfl2+9ChQ3r33XertMy+ffsqMTFR//d//1fmPEpKSjRlyhTn8w4dOmjPnj06cuSIs+2LL77Qpk2bqrT8inTo0EEbN250aXvttdcuOnJx7tvz878tNwxD8+fPL9XXz89PUukD3/LEx8erZcuWSktLU1pamvr06eMyIhQSEqIBAwZo8eLFZYa987ebWUePHtXIkSNlt9sr/O2uoqIinT592qWtQ4cOat68ucvpn35+fpVe74upzPvj3N0lz/+bFhYWavny5aXmV9naevXqpZCQEC1atMhl3T766CN9++23pu9WeSnKWj+73a7XXnutyvOszGdC06ZNJVVuHz7322Lz5s1zaZ87d64k1er2AlB5jEQBqHYfffSR9uzZU6o9Li7OOSpQHa688kqNHTtW27ZtU2hoqJYsWaLc3NwqjwCVZdSoUXrjjTeUmpqqrVu3ql+/fiosLNR//vMfPfjgg7r55ps1dOhQzZ07V4mJifrjH/+ovLw8LViwQFdccYXLtRRmvPHGGxoyZIhuueUWDRs2TIMGDZKfn5/279+vFStW6JdffnH+VtQ999yjuXPnKiEhQWPHjlVeXp4WLVqkzp07O29SUV3uvfde3X///br11ls1ePBgffHFF1q3bp2CgoIqnK5jx47q0KGDpkyZop9++kn+/v565513yryOIyYmRpL05z//WQkJCfL29tYdd9xR7rwbN26sW265RStWrFBhYaHLb2ids2DBAv3ud79T165ddd9996l9+/bKzc3Vli1b9OOPP+qLL7646Lrv27dPb775pgzDUEFBgb744gutWrVKJ0+edP79K5p20KBBSkpKUqdOndSoUSOtWbNGubm5LusWExOjhQsX6tlnn9UVV1yhkJCQUqM5lVWZ98eQIUN0+eWXa+zYsXr44Yfl7e2tJUuWKDg4WIcPH3aZX2Vra9y4sV588UUlJyerf//+GjlypPMW5xEREXrooYeqtD5V0blzZ11zzTWaNm2ajh49qpYtW2rFihU6e/ZsledZmc+EJk2aqFOnTkpLS9OVV16pli1bqkuXLs5b7p+ve/fuGj16tF577TUdP35c/fv319atW7V8+XINHz5c119//aVsAgA1xQ13BARQT1V0i3Odd7vnc7d/vvB2z+duO7xq1aoy57tt2zZnW7t27YyhQ4ca69atM7p162ZYrVajY8eOpaYt71bGnTt3LnMdyrpdeFFRkfHYY48ZkZGRRuPGjY2wsDDjtttuMw4ePOjs8/e//92Iiopy1rF06dIybx9dmVucn7/cOXPmGL179zaaNWtm+Pj4GFFRUUZKSopx4MABl75vvvmm0b59e8PHx8fo0aOHsW7dunJvcX4p291utxuPPPKIERQUZDRt2tRISEgwDhw4UKlbnH/zzTdGfHy80axZMyMoKMi47777jC+++MJl3zAMwzh79qyRkpJiBAcHGxaLxWUbqpzbRqenpxuSDIvFYvzwww9lbs+DBw8ao0aNMsLCwozGjRsbbdq0MW688UZj9erVZfY/3/n7sZeXl9GiRQujZ8+exsSJE43du3eX6n/hLc7z8/ON8ePHGx07djT8/PyMgIAAIzY21li5cqXLdDk5OcbQoUON5s2bG5Kc+2JZf4tzyrvFeWXeH4ZhGDt27DBiY2MNHx8f4/LLLzfmzp1b5jzLq62sv7VhGEZaWprRs2dPw2q1Gi1btjTuvPNO48cff3TpM3r0aMPPz69UTeXdev1CFb2Xzzl48KARHx9vWK1WIzQ01Jg+fbpzf6nM58KF7yPDqNxnwubNm42YmBjDx8fHZb8ta93OnDljPPXUU875hYeHG9OmTXO5Jb9h/PZ3LWs7lHX7eQA1x2IYXIkIoO6JiIhQly5d9P7777u7FAAA0MBwTRQAAAAAmECIAgAAAAATCFEAAAAAYALXRAEAAACACYxEAQAAAIAJhCgAAAAAMKHBhyjj///RRM5qBAAAAFAZDT5EnThxQgEBATpx4oS7SwEAAABQBzT4EAUAAAAAZhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACY0cncBAAAAAFCTHHaH9u/MkS2/SAFBTRXVM0xe3lUfT/KokaiNGzdq2LBhat26tSwWi959992LTrN+/XpFR0fLarXqiiuu0LJly2q8TgAAAAB1Q1ZmtqYPe1tz739ff388U3Pvf1/Th72trMzsKs/To0JUYWGhunfvrgULFlSqf3Z2toYOHarrr79eu3bt0qRJk3Tvvfdq3bp1NVwpAAAAAE+XlZmtxVPTdSyv0KX9WF6hFk9Nr3KQshiGYVRHgdXNYrFozZo1Gj58eLl9HnnkEX3wwQf6+uuvnW133HGHjh8/rrVr11ZqOQUFBQoICJDNZpO/v/+llg0AAADAAzjsDk0f9napAHW+wFA/Pf/eSNOn9nnUSJRZW7ZsUXx8vEtbQkKCtmzZUu40xcXFKigocHkAAAAAqF/278ypMEBJ0rHcQu3fmWN63nU6ROXk5Cg0NNSlLTQ0VAUFBTp16lSZ08yaNUsBAQHOR3h4eG2UCgAAAKAW2fKLqrXf+ep0iKqKadOmyWazOR8//PCDu0sCAAAAUM0CgppWa7/z1elbnIeFhSk3N9elLTc3V/7+/mrSpEmZ01itVlmt1tooDwAAAICbRPUMU2CI30WviYrqGWZ63nV6JKpv377KyMhwaUtPT1ffvn3dVBEAAAAAT+Dl7aWkKXEV9kmaHFel34vyqBB18uRJ7dq1S7t27ZL06y3Md+3apcOHD0v69VS8UaNGOfvff//9+u677zR16lTt2bNHr776qlauXKmHHnrIHeUDAAAA8CDRAyM1bvZgtQh2PWUvMNRP42YPVvTAyCrN16Nucb5+/Xpdf/31pdpHjx6tZcuWacyYMTp06JDWr1/vMs1DDz2kb775Rm3bttUTTzyhMWPGVHqZ3OIcAAAAqN9OnSzWpAHLJUkp8xPV6Zq2VRqBOsejQpQ7EKIAAACA+q341Bn9ud9SSdLfPk2WtUnjS5qfR53OBwAAAACejhAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJjRydwEAANR1DrtD+3fmyJZfpICgporqGSYvb76nBID6ihAFAMAlyMrM1so5m3Usr9DZFhjip6QpcYoeGOnGygAANYWvyQAAqKKszGwtnpruEqAk6VheoRZPTVdWZrabKgMA1CRGogAAqAKH3aGVczZX2CdtziZd3ac1p/YBgJsVnzpbrfMjRAEAUAX7d+aUGoG60PG8Ik0asLyWKgIA1Ba+GgMAoAps+UXuLgEAYFKH7qHy8b30cSRGogAAqIKAoKaV6pcyP1FR0a1quBoAQGX4+DaSxWK55PkQogAAqIKonmEKDPGr8JS+wFA/dbqmLddEAUA9w6c6AABV4OXtpaQpcRX2SZocR4ACgHqIT3YAAKooemCkxs0erBbBrqf2BYb6adzswfxOFADUUxbDMAx3F+FOBQUFCggIkM1mk7+/v7vLAQDUQadOFjvvwpcyP5FT+ACgnuMTHgCAS3R+YIqKbkWAAoB6jk95AAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmeFyIWrBggSIiIuTr66vY2Fht3bq1wv7z5s3TVVddpSZNmig8PFwPPfSQTp8+XUvVAoD7OOwO7d3+s7auPaC923+Ww+5wd0kAADQIjdxdwPnS0tKUmpqqRYsWKTY2VvPmzVNCQoL27t2rkJCQUv3feustPfroo1qyZIni4uK0b98+jRkzRhaLRXPnznXDGgBA7cjKzNbKOZt1LK/Q2RYY4qekKXGKHhjpxsoAAKj/PGokau7cubrvvvuUnJysTp06adGiRWratKmWLFlSZv/Nmzfr2muv1R//+EdFRERoyJAhGjly5EVHrwCgLsvKzNbiqekuAUqSjuUVavHUdGVlZrupMgAAGgaPGYkqKSnRjh07NG3aNGebl5eX4uPjtWXLljKniYuL05tvvqmtW7eqT58++u677/Thhx/q7rvvrq2yAaBWOewOrZyzucI+aXM26eo+reXl7VHfk9VrxafOursEAEAt8pgQlZ+fL7vdrtDQUJf20NBQ7dmzp8xp/vjHPyo/P1+/+93vZBiGzp49q/vvv1/Tp08vdznFxcUqLi52Pi8oKKieFQCAWrB/Z06pEagLHc8r0qQBy2upIgAAGp46/TXl+vXr9fzzz+vVV19VVlaW/vnPf+qDDz7QM888U+40s2bNUkBAgPMRHh5eixUDwKWx5Re5uwRUoEP3UPn4esz3kwCAGuIxn/RBQUHy9vZWbm6uS3tubq7CwsLKnOaJJ57Q3XffrXvvvVeS1LVrVxUWFupPf/qTHnvsMXl5lc6I06ZNU2pqqvN5QUEBQQpAnREQ1LRS/VLmJyoqulUNV4ML+fg2ksVicXcZAIAa5jEhysfHRzExMcrIyNDw4cMlSQ6HQxkZGZowYUKZ0xQVFZUKSt7e3pIkwzDKnMZqtcpqtVZf4QBQi6J6hikwxK/CU/oCQ/3U6Zq2XBMFAEAN8aj/YVNTU/X6669r+fLl+vbbb/XAAw+osLBQycnJkqRRo0a53Hhi2LBhWrhwoVasWKHs7Gylp6friSee0LBhw5xhCgDqEy9vLyVNiauwT9LkOAIUAAA1yGNGoiRpxIgROnLkiGbMmKGcnBz16NFDa9eudd5s4vDhwy4jT48//rgsFosef/xx/fTTTwoODtawYcP03HPPuWsVAKDGRQ+M1LjZg5X20iYdP/LbNVKBoX5KmszvRAEAUNMsRnnnvTUQBQUFCggIkM1mk7+/v7vLAYBKO3Wy2HkXvpT5iZzCBwBALeF/WwCoo84PTFHRrQhQAADUEv7HBQAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgQiN3FwAAZjnsDu3fmSNbfpECgpoqqmeYvLz5TggAANQOQhSAOiUrM1sr52zWsbxCZ1tgiJ+SpsQpemCkGysDAAANBV/dAqgzsjKztXhqukuAkqRjeYVaPDVdWZnZbqoMAAA0JIxEAagTHHaHVs7ZXGGftDmbdHWf1g3m1L7iU2fdXQIAAA0SIQpAnbB/Z06pEagLHc8r0qQBy2upIgAA0FA1jK9rAdR5tvwid5fgsTp0D5WPL9+JAQBQW/hfF0CdEBDUtFL9UuYnKiq6VQ1X41l8fBvJYrG4uwwAABoMQhSAOiGqZ5gCQ/wqPKUvMNRPna5p22CuiQIAAO7BkQaAOsHL20tJU+Iq7JM0OY4ABQAAahxHGwDqjOiBkRo3e7BaBLue2hcY6qdxswfzO1EAAKBWWAzDMNxdhDsVFBQoICBANptN/v7+7i4HQCWcOlnsvAtfyvxETuEDAAC1iqMOAHXO+YEpKroVAQoAANQqjjwAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJjgcSFqwYIFioiIkK+vr2JjY7V169YK+x8/flzjx49Xq1atZLVadeWVV+rDDz+spWoBmOWwO7R3+8/auvaA9m7/WQ67w90lAQAAmNLI3QWcLy0tTampqVq0aJFiY2M1b948JSQkaO/evQoJCSnVv6SkRIMHD1ZISIhWr16tNm3a6Pvvv1eLFi1qv3gAF5WVma2VczbrWF6hsy0wxE9JU+IUPTDSjZUBAABUnsUwDMPdRZwTGxur3r1765VXXpEkORwOhYeHKyUlRY8++mip/osWLdJLL72kPXv2qHHjxlVaZkFBgQICAmSz2eTv739J9QMoX1ZmthZPTS/39XGzB1c6SBWfOqM/91sqSfrbp8myNqna+x8AAKAqPOZ0vpKSEu3YsUPx8fHONi8vL8XHx2vLli1lTvPee++pb9++Gj9+vEJDQ9WlSxc9//zzstvttVU2gEpw2B1aOWdzhX3S5mzSqZPFKj51phKPs7VUOQAAQGkeczpffn6+7Ha7QkNDXdpDQ0O1Z8+eMqf57rvvlJmZqTvvvFMffvihDhw4oAcffFBnzpzRzJkzy5ymuLhYxcXFzucFBQXVtxIAyrR/Z47LKXxlOZ5XpEkDltdSRQAAAFXnMSNRVeFwOBQSEqLXXntNMTExGjFihB577DEtWrSo3GlmzZqlgIAA5yM8PLwWKwYaJlt+UY3Mt0P3UPn4esx3QQAAoIHwmKOPoKAgeXt7Kzc316U9NzdXYWFhZU7TqlUrNW7cWN7e3s62q6++Wjk5OSopKZGPj0+paaZNm6bU1FTn84KCAoIUUMMCgppWql/K/ERFRbeq9Hx9fBvJYrFUtSwAAIAq8ZiRKB8fH8XExCgjI8PZ5nA4lJGRob59+5Y5zbXXXqsDBw7I4fjtFsn79u1Tq1atygxQkmS1WuXv7+/yAFCzonqGKTDEr8I+gaF+6nRNW1mbNK70gwAFAADcwWNClCSlpqbq9ddf1/Lly/Xtt9/qgQceUGFhoZKTkyVJo0aN0rRp05z9H3jgAR09elQTJ07Uvn379MEHH+j555/X+PHj3bUKAMrg5e2lpClxFfZJmhwnL2+P+kgCAAAok8eczidJI0aM0JEjRzRjxgzl5OSoR48eWrt2rfNmE4cPH5aX128HWeHh4Vq3bp0eeughdevWTW3atNHEiRP1yCOPuGsVAJQjemCkxs0erLSXNun4kd+ukQoM9VPSZH4nCgAA1B0e9TtR7sDvRAG169TJYudd+FLmJ6rTNW0ZgQIAAHUKRy4AatX5gSkquhUBCgAA1DlVPp3vzJkzysnJUVFRkYKDg9WyZcvqrAsAAAAAPJKpr4BPnDihhQsXqn///vL391dERISuvvpqBQcHq127drrvvvu0bdu2mqoVAAAAANyu0iFq7ty5ioiI0NKlSxUfH693331Xu3bt0r59+7RlyxbNnDlTZ8+e1ZAhQ5SYmKj9+/fXZN0AAAAA4BaVPp1v27Zt2rhxozp37lzm63369NE999yjRYsWaenSpfr0008VFRVVbYUCAAAAgCeodIh6++23K9XParXq/vvvr3JBAAAAAODJqnRbrBMnTlR3HQAAAABQJ1QpRPXr1085OTnVXQsAAAAAeLwqhaiePXsqNjZWe/bscWnftWuXfv/731dLYQAAAADgiaoUopYuXaoxY8bod7/7nT777DPt27dPSUlJiomJkbe3d3XXCAAAAAAeo8o/tvvUU0/JarVq8ODBstvtGjRokLZs2aI+ffpUZ30AAAAA4FGqNBKVm5uriRMn6tlnn1WnTp3UuHFjjRkzhgAFAAAAoN6rUoiKjIzUxo0btWrVKu3YsUPvvPOO/vSnP+mll16q7voAAAAAwKNU6XS+JUuW6I477nA+T0xM1CeffKIbb7xRhw4d0oIFC6qtQAAAAADwJFUaiTo/QJ0THR2tzZs3KzMz85KLAgAAAABPVekQdfjw4Yv2iYiI0ObNmyVJP/30U9WrAgAAAAAPVekQ1bt3b40bN07btm0rt4/NZtPq1avVpUsXvfPOO9VSIAAAAAB4kkpfE/XNN9/oueee0+DBg+Xr66uYmBi1bt1avr6+OnbsmL755hvt3r1b0dHRmj17Nj+6CwAAAKBeshiGYZiZ4NSpU/rggw/02Wef6fvvv9epU6cUFBSknj17KiEhQV26dKmpWmtEQUGBAgICZLPZ5O/v7+5ygHqv+NQZ/bnfUknS3z5NlrVJYzdXBAAAYI7pu/M1adJEt912m2677baaqAdAPeewO5z/3p/1izpd01Ze3lW6xw0AAIBbmB6Jqm8YiQJqT1ZmttJe2qTjR4qcbYEhfkqaEqfogZFurAwAAKDyqiVEffPNN/rXv/6lFi1aqHPnzuratasCAwOro74aR4gCakdWZrYWT00v9/VxswcTpAAAQJ1QLefQ3HTTTWratKkKCwv197//XYMGDVKHDh2qY9YA6gGH3aGVczZX2GflXza7nOoHAADgqUxfE1WWsLAwTZw40aXNbrdXx6wB1AP7d+boWF5hhX2O5RZq/84cXdWrdS1VBQAAUDXVMhI1aNAgLV261KXN29u7OmYNoB6w5RddvJOJfgAAAO5kKkQVFZV9gLN9+3Y9+eSTioyMVFJSkp577jn9+9//rpYCAdR9AUFNq7UfAACAO5kKUQEBAcrPzy/V/sEHH+j777/Xl19+qYceekghISHKyMiotiIB1G1RPcMUGOJXYZ/AUD9F9QyrpYoAAACqzlSIstvtcjh+u/C7X79+ys3NdT5v3ry5+vbtq/vuu0/z5s2rtiIB1G1e3l5KmhJXYZ+kyXH8XhQAAKgTLumIZdeuXSosrPhicQCQpOiBkRo3e7BaBLueshcY6sftzQEAQJ1SLXfnA4DKiB4Yqav7tNakAcslSSnzE9XpmraMQAEAgDrF9JHLW2+9paysLJ05c6Ym6gFQz50fmKKiWxGgAABAnWNqJKpfv36aOXOmTpw4ocaNG+vs2bOaOXOm+vXrp549e6pbt26yWq01VSsAAAAAuJ2pELVhwwZJ0v79+7Vjxw5lZWUpKytLjz76qI4fP65GjRrp6quv1hdffFEjxQIAAACAu1XpmqioqChFRUXpjjvucLZlZ2dr+/bt2rlzZ7UVBwAAAACeptpuLBEZGanIyEjdfvvt1TVLAAAAAPA4XNENAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEzwuBC1YMECRUREyNfXV7Gxsdq6dWulpluxYoUsFouGDx9eswUCDZzD7tDe7T9r69oD2rv9ZznsDneXBAAAUKsaubuA86WlpSk1NVWLFi1SbGys5s2bp4SEBO3du1chISHlTnfo0CFNmTJF/fr1q8VqgYYnKzNbK+ds1rG8QmdbYIifkqbEKXpgpBsrAwAAqD0eNRI1d+5c3XfffUpOTlanTp20aNEiNW3aVEuWLCl3GrvdrjvvvFNPPfWU2rdvX4vVAg1LVma2Fk9NdwlQknQsr1CLp6YrKzPbTZUBAADULo8JUSUlJdqxY4fi4+OdbV5eXoqPj9eWLVvKne7pp59WSEiIxo4dWxtlAg2Sw+7QyjmbK+yTNmeTTp0sVvGpMxd5nK2lqgEAAGqGx5zOl5+fL7vdrtDQUJf20NBQ7dmzp8xpPvvsM/3973/Xrl27Kr2c4uJiFRcXO58XFBRUqV6gIdm/M6fUCNSFjucVadKA5bVUEQAAgPt4zEiUWSdOnNDdd9+t119/XUFBQZWebtasWQoICHA+wsPDa7BKoH6w5RdV+zw7dA+Vj6/HfI8DAABQaR5zBBMUFCRvb2/l5ua6tOfm5iosLKxU/4MHD+rQoUMaNmyYs83h+PUuYY0aNdLevXvVoUOHUtNNmzZNqampzucFBQUEKeAiAoKaVqpfyvxERUW3qlRfH99Gslgsl1IWAACAW3hMiPLx8VFMTIwyMjKctyl3OBzKyMjQhAkTSvXv2LGjvvrqK5e2xx9/XCdOnND8+fPLDUZWq1VWq7Xa6wfqs6ieYQoM8avwlL7AUD91uqatvLzr7AA3AABApXhMiJKk1NRUjR49Wr169VKfPn00b948FRYWKjk5WZI0atQotWnTRrNmzZKvr6+6dOniMn2LFi0kqVQ7gEvj5e2lpClxWjw1vdw+SZPjCFAAAKBB8KgQNWLECB05ckQzZsxQTk6OevToobVr1zpvNnH48GF5eXGQBrhD9MBIjZs9WGkvbdLxI79dIxUY6qekyfxOFAAAaDgshmEY7i7CnQoKChQQECCbzSZ/f393lwN4vFMni5134UuZn8gpfAAAoMHhyAeAKecHpqjoVgQoAADQ4HD0AwAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADAhEbuLgBA5TnsDu3fmSNbfpECgpoqqmeYvLz5LgQAAKA2EaKAOiIrM1sr52zWsbxCZ1tgiJ+SpsQpemCkGysDAABoWPgKG6gDsjKztXhqukuAkqRjeYVaPDVdWZnZbqoMAACg4WEkCvBwDrtDK+dsrrBP2pxNurpP61o5ta/41NkaXwYAAIAnI0QBHm7/zpxSI1AXOp5XpEkDltdSRQAAAA0bp/MBHs6WX+TuEsrUoXuofHz5HgYAADQ8HAEBHi4gqGml+qXMT1RUdKsaruY3Pr6NZLFYam15AAAAnoIQBXi4qJ5hCgzxq/CUvsBQP3W6pi23OwcAAKgFHHEBHs7L20tJU+Iq7JM0OY4ABQAAUEs46gLqgOiBkRo3e7BaBLue2hcY6qdxswfzO1EAAAC1yGIYhuHuItypoKBAAQEBstls8vf3d3c5QIVOnSx23oUvZX4ip/ABAAC4AUdfQB1yfmCKim5FgAIAAHADjsAAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJjgcSFqwYIFioiIkK+vr2JjY7V169Zy+77++uvq16+fAgMDFRgYqPj4+Ar7A5fCYXdo7/aftXXtAe3d/rMcdoe7SwIAAIAbNHJ3AedLS0tTamqqFi1apNjYWM2bN08JCQnau3evQkJCSvVfv369Ro4cqbi4OPn6+urFF1/UkCFDtHv3brVp08YNa4D6KiszWyvnbNaxvEJnW2CIn5KmxCl6YKQbKwMAAEBtsxiGYbi7iHNiY2PVu3dvvfLKK5Ikh8Oh8PBwpaSk6NFHH73o9Ha7XYGBgXrllVc0atSoSi2zoKBAAQEBstls8vf3v6T6UT9lZWZr8dT0cl8fN3twrQWp4lNn9Od+SyVJf/s0WdYmjWtluQAAAPiNx4xElZSUaMeOHZo2bZqzzcvLS/Hx8dqyZUul5lFUVKQzZ86oZcuWNVUmGhiH3aGVczZX2CdtziZd3ae1vLxr/uzY4lNna3wZAAAAqJjHhKj8/HzZ7XaFhoa6tIeGhmrPnj2Vmscjjzyi1q1bKz4+vtw+xcXFKi4udj4vKCioWsFoEPbvzHE5ha8sx/OKNGnA8lqqCAAAAO7mcTeWqKoXXnhBK1as0Jo1a+Tr61tuv1mzZikgIMD5CA8Pr8UqUdfY8ovcXUKZOnQPlY+vx3wHAgAA0KB4zFFYUFCQvL29lZub69Kem5ursLCwCqedM2eOXnjhBf3nP/9Rt27dKuw7bdo0paamOp8XFBQQpFCugKCmleqXMj9RUdGtaria3/j4NpLFYqm15QEAAOA3HhOifHx8FBMTo4yMDA0fPlzSrzeWyMjI0IQJE8qdbvbs2Xruuee0bt069erV66LLsVqtslqt1VU26rmonmEKDPGr8JS+wFA/dbqmba1cEwUAAAD386ijvtTUVL3++utavny5vv32Wz3wwAMqLCxUcnKyJGnUqFEuN5548cUX9cQTT2jJkiWKiIhQTk6OcnJydPLkSXetAuoZL28vJU2Jq7BP0uQ4AhQAAEAD4jEjUZI0YsQIHTlyRDNmzFBOTo569OihtWvXOm82cfjwYXl5/XawunDhQpWUlOi2225zmc/MmTP15JNP1mbpqMeiB0Zq3OzBSntpk44f+e0aqcBQPyVN5neiAAAAGhqP+p0od+B3olBZp04WO+/ClzI/kVP4AAAAGiiOAIFKOj8wRUW3IkABAAA0UBwFAgAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBAAAAgAmEKAAAAAAwoZG7C/BkDrtD+3fmyJZfpICgporqGSYvb3InAAAA0JARosqRlZmtlXM261heobMtMMRPSVPiFD0w0o2VAQAAAHAnhlXKkJWZrcVT010ClCQdyyvU4qnpysrMdlNlAAAAANyNkagLOOwOrZyzucI+aXM26eo+rTm1r4EpPnXW3SUAAADAAxCiLrB/Z06pEagLHc8r0qQBy2upIgAAAACehKGUC9jyi9xdAjxch+6h8vHl+wcAAICGiiPBCwQENa1Uv5T5iYqKblXD1cAT+fg2ksVicXcZAAAAcBNC1AWieoYpMMSvwlP6AkP91OmatlwTBQAAADRApIALeHl7KWlKXIV9kibHEaAAAACABookUIbogZEaN3uwWgS7ntoXGOqncbMH8ztRAAAAQANmMQzDcHcR7lRQUKCAgADZbDb5+/u7vHbqZLHzLnwp8xM5hQ8AAAAAI1EVOT8wRUW3IkABAAAAIEQBAAAAgBmEKAAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACYQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADDB40LUggULFBERIV9fX8XGxmrr1q0V9l+1apU6duwoX19fde3aVR9++GEtVQoAAACgIfKoEJWWlqbU1FTNnDlTWVlZ6t69uxISEpSXl1dm/82bN2vkyJEaO3asdu7cqeHDh2v48OH6+uuva7lyAAAAAA2FxTAMw91FnBMbG6vevXvrlVdekSQ5HA6Fh4crJSVFjz76aKn+I0aMUGFhod5//31n2zXXXKMePXpo0aJFlVpmQUGBAgICZLPZ5O/v7/Ja8akz+nO/pZKkv32aLGuTxlVdNQAAAAD1hMeMRJWUlGjHjh2Kj493tnl5eSk+Pl5btmwpc5otW7a49JekhISEcvsDAAAAwKVq5O4CzsnPz5fdbldoaKhLe2hoqPbs2VPmNDk5OWX2z8nJKXc5xcXFKi4udj632WySfh2RKtX31BmV2E85X7eeYSQKAAAAqO+aN28ui8VS7useE6Jqy6xZs/TUU0+Vag8PD69wumVhk2qoIgAAAACepKxLfc7nMSEqKChI3t7eys3NdWnPzc1VWFhYmdOEhYWZ6i9J06ZNU2pqqvO5w+HQ0aNHddlll5WZNgsKChQeHq4ffvihwg0J1Bb2SXga9kl4GvZJeBL2x7qpefPmFb7uMSHKx8dHMTExysjI0PDhwyX9GnAyMjI0YcKEMqfp27evMjIyNGnSJGdbenq6+vbtW+5yrFarrFarS1uLFi0uWp+/vz87PjwK+yQ8DfskPA37JDwJ+2P94jEhSpJSU1M1evRo9erVS3369NG8efNUWFio5ORkSdKoUaPUpk0bzZo1S5I0ceJE9e/fX3/5y180dOhQrVixQtu3b9drr73mztUAAAAAUI95VIgaMWKEjhw5ohkzZignJ0c9evTQ2rVrnTePOHz4sLy8fruhYFxcnN566y09/vjjmj59uqKiovTuu++qS5cu7loFAAAAAPWcR4UoSZowYUK5p++tX7++VNvtt9+u22+/vcbqsVqtmjlzZqlTAAF3YZ+Ep2GfhKdhn4QnYX+snzzqx3YBAAAAwNN5zI/tAgAAAEBdQIgCAAAAABMIUQAAAABgAiEKAAAAAEwgRFVgwYIFioiIkK+vr2JjY7V161Z3l4QGbOPGjRo2bJhat24ti8Wid999190loQGbNWuWevfurebNmyskJETDhw/X3r173V0WGrCFCxeqW7duzh807du3rz766CN3lwU4vfDCC7JYLJo0aZK7S0E1IESVIy0tTampqZo5c6aysrLUvXt3JSQkKC8vz92loYEqLCxU9+7dtWDBAneXAmjDhg0aP368Pv/8c6Wnp+vMmTMaMmSICgsL3V0aGqi2bdvqhRde0I4dO7R9+3YNHDhQN998s3bv3u3u0gBt27ZNixcvVrdu3dxdCqoJtzgvR2xsrHr37q1XXnlFkuRwOBQeHq6UlBQ9+uijbq4ODZ3FYtGaNWs0fPhwd5cCSJKOHDmikJAQbdiwQdddd527ywEkSS1bttRLL72ksWPHursUNGAnT55UdHS0Xn31VT377LPq0aOH5s2b5+6ycIkYiSpDSUmJduzYofj4eGebl5eX4uPjtWXLFjdWBgCeyWazSfr1oBVwN7vdrhUrVqiwsFB9+/Z1dzlo4MaPH6+hQ4e6HFei7mvk7gI8UX5+vux2u0JDQ13aQ0NDtWfPHjdVBQCeyeFwaNKkSbr22mvVpUsXd5eDBuyrr75S3759dfr0aTVr1kxr1qxRp06d3F0WGrAVK1YoKytL27Ztc3cpqGaEKADAJRk/fry+/vprffbZZ+4uBQ3cVVddpV27dslms2n16tUaPXq0NmzYQJCCW/zwww+aOHGi0tPT5evr6+5yUM0IUWUICgqSt7e3cnNzXdpzc3MVFhbmpqoAwPNMmDBB77//vjZu3Ki2bdu6uxw0cD4+PrriiiskSTExMdq2bZvmz5+vxYsXu7kyNEQ7duxQXl6eoqOjnW12u10bN27UK6+8ouLiYnl7e7uxQlwKrokqg4+Pj2JiYpSRkeFsczgcysjI4NxqAJBkGIYmTJigNWvWKDMzU5GRke4uCSjF4XCouLjY3WWggRo0aJC++uor7dq1y/no1auX7rzzTu3atYsAVccxElWO1NRUjR49Wr169VKfPn00b948FRYWKjk52d2loYE6efKkDhw44HyenZ2tXbt2qWXLlrr88svdWBkaovHjx+utt97Sv/71LzVv3lw5OTmSpICAADVp0sTN1aEhmjZtmm644QZdfvnlOnHihN566y2tX79e69atc3dpaKCaN29e6jpRPz8/XXbZZVw/Wg8QosoxYsQIHTlyRDNmzFBOTo569OihtWvXlrrZBFBbtm/fruuvv975PDU1VZI0evRoLVu2zE1VoaFauHChJGnAgAEu7UuXLtWYMWNqvyA0eHl5eRo1apR++eUXBQQEqFu3blq3bp0GDx7s7tIA1EP8ThQAAAAAmMA1UQAAAABgAiEKAAAAAEwgRAEAAACACYQoAAAAADCBEAUAAAAAJhCiAAAAAMAEQhQAAAAAmECIAgAAAAATCFEAAAAAYAIhCgAAAABMIEQBABqst99+W02aNNEvv/zibEtOTla3bt1ks9ncWBkAwJNZDMMw3F0EAADuYBiGevTooeuuu04vv/yyZs6cqSVLlujzzz9XmzZt3F0eAMBDNXJ3AQAAuIvFYtFzzz2n2267TWFhYXr55Zf16aefEqAAABViJAoA0OBFR0dr9+7d+vjjj9W/f393lwMA8HBcEwUAaNDWrl2rPXv2yG63KzQ01N3lAADqAEaiAAANVlZWlgYMGKDFixdr2bJl8vf316pVq9xdFgDAw3FNFACgQTp06JCGDh2q6dOna+TIkWrfvr369u2rrKwsRUdHu7s8AIAHYyQKANDgHD16VHFxcRowYIAWLVrkbB86dKjsdrvWrl3rxuoAAJ6OEAUAAAAAJnBjCQAAAAAwgRAFAAAAACYQogAAAADABEIUAAAAAJhAiAIAAAAAEwhRAAAAAGACIQoAAAAATCBEAQAAAIAJhCgAAAAAMIEQBQAAAAAmEKIAAAAAwARCFAAAAACY8P8B1Pom9LEcV/gAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Divide the figure in two subfigures\n",
    "f, ax = plt.subplots(2, 1, figsize=(10,6), sharex=True, height_ratios=[1, 3]) # The subfigures share the same x axis\n",
    "\n",
    "# In the first subfigure we show our observations\n",
    "ax[0].set_yticks([])                                        # Do not show the ticks in the vertical axis\n",
    "ax[0].set_ylabel(\"\")                                        # Do not show any label in the vertical axis\n",
    "ax[0].spines[['left', 'right', 'top']].set_visible(False)   # Do not plot the left, right and top margins\n",
    "ax[0].text(2.2, 0.75, \"Data\", fontsize=12)\n",
    "\n",
    "for i in range(len(Dat)):\n",
    "    ax[0].axvline(x=Dat[i], ymax=0.2, c='r')                # Plot the observations as a rugplot\n",
    "\n",
    "# In the second subfigure we show the ECDF\n",
    "ax[1].spines[['right', 'top']].set_visible(False)           # Do not plot the right and top margins\n",
    "\n",
    "# Plot the ECDF\n",
    "ax[1].scatter(ecdf.x, ecdf.y, color='rebeccapurple')        \n",
    "ax[1].step(ecdf.x, ecdf.y, where='post', color='rebeccapurple')\n",
    "\n",
    "# Add title and labels to the plot\n",
    "ax[1].set_title(\"Empirical Cumulative Distribution Function\", fontsize=12)\n",
    "ax[1].set_xlabel(r'$x$')\n",
    "ax[1].set_ylabel(r'$F_n(x)$')\n",
    "\n",
    "# Modify the limits of the vertical and horizontal axes\n",
    "ax[1].set_ylim(0, 1.05)\n",
    "ax[1].set_xlim(Dat.min()-0.05, Dat.max()+0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantiles, quartiles, percentiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.792    , 1.3689085, 2.422    ])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Quantiles of specific probabilities\n",
    "p = [0.11, 1-np.exp(-1), 0.89]\n",
    "np.quantile(Dat, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.05, 1.3 , 1.7 ])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Quartiles\n",
    "p = [0.25, 0.5, 0.75]\n",
    "np.quantile(Dat, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Central tendency statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mode = mode(Dat)\n",
    "Median = np.median(Dat)\n",
    "Mean = Dat.mean()\n",
    "GeomMean = gmean(Dat)\n",
    "HarmMean = hmean(Dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "Weights = norm.pdf(Dat, Dat.mean(), Dat.std()) # Weights for the weighted mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGxCAYAAACKvAkXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwPUlEQVR4nO3df1RVdb7/8deBhnP8iTIEB4wCf4wOo0iK8CUzf5HgtFx6b9NS1zQaq2szJN38nltNNAqazmDmeKmrA42z7IdleqcpZ9k4mHMSyxuJQU6Z1soGr5ocQB1BcYTinO8ffj15BjBAOPvAfj7W2ivO53z257x3p+TlZ3/23haPx+MRAACAiQQZXQAAAIC/EYAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIAAAIDpEIBa4fF4VF9fL+4RCQBA70QAasX58+cVGhqq8+fPG10KAADoBgERgDZs2KDY2FjZbDalpKSorKyszb6vv/66kpKSNGjQIPXr10+JiYnavHmzT5/77rtPFovFZ8vIyOjuwwAAAD3EDUYXsG3bNjkcDhUVFSklJUUFBQVKT0/XZ599poiIiBb9w8LC9Itf/EKjRo1SSEiI3nzzTWVmZioiIkLp6enefhkZGXr++ee9r61Wq1+OBwAABD6L0Q9DTUlJ0YQJE7R+/XpJktvtVkxMjB566CE9/vjj7Rpj3Lhxuuuuu7Ry5UpJl2eAzp07p+3bt3eqpvr6eoWGhqqurk4DBw7s1BgAACBwGXoKrKmpSeXl5UpLS/O2BQUFKS0tTaWlpd+6v8fjkdPp1GeffaY77rjD572SkhJFRERo5MiRysrK0pkzZ9ocp7GxUfX19T4bAADovQw9BXb69Gk1NzcrMjLSpz0yMlKffvppm/vV1dVpyJAhamxsVHBwsH7zm9/ozjvv9L6fkZGhf/3Xf1VcXJy++OILPfHEE5o5c6ZKS0sVHBzcYrz8/HytWLGi6w4MAAAENMPXAHXGgAEDdPDgQV24cEFOp1MOh0NDhw7VlClTJEnz5s3z9h0zZowSEhI0bNgwlZSUaPr06S3Gy8nJkcPh8L6ur69XTExMtx8HAAAwhqEBKDw8XMHBwaqurvZpr66ult1ub3O/oKAgDR8+XJKUmJioI0eOKD8/3xuA/tnQoUMVHh6uo0ePthqArFYri6QBADARQ9cAhYSEaPz48XI6nd42t9stp9Op1NTUdo/jdrvV2NjY5vsnT57UmTNnFBUVdV31AgCA3sHwU2AOh0MLFy5UUlKSkpOTVVBQoIaGBmVmZkqSFixYoCFDhig/P1/S5fU6SUlJGjZsmBobG7Vz505t3rxZhYWFkqQLFy5oxYoVuvvuu2W32/XFF1/oscce0/Dhw30ukwcAAOZleACaO3euamtrlZubK5fLpcTERBUXF3sXRh8/flxBQd9MVDU0NOjBBx/UyZMn1adPH40aNUovv/yy5s6dK0kKDg7WRx99pBdffFHnzp1TdHS0ZsyYoZUrV3KaCwAASAqA+wAFIu4DBABA72b4DBCAbzS7PSqrPKua85cUMcCm5LgwBQdZjC4LAHodAhAQIIoPVWnFjsOqqrvkbYsKtSlvVrwyRrOAHwC6UkA8DBUwu+JDVcp6ucIn/EiSq+6Ssl6uUPGhKoMq+0az26PSL87ojwe/VOkXZ9Ts5uw5gJ6LGSDAYM1uj1bsOKzW4oRHkkXSih2HdWe83bDTYcxOAehtmAECDFZWebbFzM/VPJKq6i6prPKs/4q6Sk+YnQKAjiIAAQarOd92+OlMv670bbNT0uXZKU6HAehpOAUG+EHT125tLj2m/z17UbeE9dVPUmMVcsPlv39EDLC1a4z29utKHZmdSh32Xf8VBgDXiQAEdLP8nYe18d1KXT1J8sudR7RoUpxyfhiv5LgwRYXa5Kq71OpMi0WSPfTyJfH+FsizUwBwPTgFBnSj/J2H9dw7vuFHktwe6bl3KpW/87CCgyzKmxUv6XLYudqV13mz4g1ZAB3Is1MAcD0IQEA3afrarY3vVl6zz8Z3K9X0tVsZo6NUeO84RQz0fVyLPdSmwnvHGXal1ZXZqbail0WXrwYzYnYKAK4HAQjoJptLj7WY+flnbs/lfpKUMTpKf3FM9r73QuYE7fv5NEMvMw/k2SkAuB4EIKCb/O/Zix3ud3WQCJTHYATq7BQAXA8WQQPd5Jawvl3az0gZo6M0cXi4xix/S9Ll2alJI24MiIAGAJ3BDBBwHa71eIifpMbq2/JBkOVyv54gEGenAKCzmAECOunbHg8RckOQFk2K03PvtL0QetGkOO/9gAAA/sOfvEAntPfxEDk/jNdP74hrMRMUZJF+esfl+wABAPyPAAR0UEcfD5Hzw3jt+/lUb5//mzZCn66cSfgBAAMRgIAO6szDS68+zfXj/3Mzp70AwGD8KQx0EI+HAICejwAEdBCPhwCAno8ABHQQj4cAgJ6PAAR0EI+HAICejwAEdAKPhwCAno0bIQKdxOMhAKDnYgYIuA48HgIAeiYCEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB2eBQbTaXZ7VFZ5VjXnLyligI1HWACACRGAYCrFh6q0YsdhVdVd8rZFhdqUNyueJ7gDgIlwCgymUXyoSlkvV/iEH0ly1V1S1ssVKj5UZVBlAAB/IwDBFJrdHq3YcVieVt670rZix2E1u1vrAQDobQIiAG3YsEGxsbGy2WxKSUlRWVlZm31ff/11JSUladCgQerXr58SExO1efNmnz4ej0e5ubmKiopSnz59lJaWps8//7y7DwMBrKzybIuZn6t5JFXVXVJZ5Vn/FQUAMIzhAWjbtm1yOBzKy8tTRUWFxo4dq/T0dNXU1LTaPywsTL/4xS9UWlqqjz76SJmZmcrMzNSuXbu8fdasWaNnn31WRUVF2r9/v/r166f09HRdutT2L0D0bjXn2/fdt7cfAKBnMzwArVu3TosWLVJmZqbi4+NVVFSkvn37atOmTa32nzJliv7lX/5F3//+9zVs2DA9/PDDSkhI0L59+yRdnv0pKCjQ0qVLNXv2bCUkJOill17SqVOntH37dj8eGQJJxABbl/YDAPRshgagpqYmlZeXKy0tzdsWFBSktLQ0lZaWfuv+Ho9HTqdTn332me644w5JUmVlpVwul8+YoaGhSklJaXPMxsZG1dfX+2zoXZLjwhQValNbF7tbdPlqsOS4MH+WBQAwiKEB6PTp02publZkZKRPe2RkpFwuV5v71dXVqX///goJCdFdd92l//qv/9Kdd94pSd79OjJmfn6+QkNDvVtMTMz1HBYCUHCQRXmz4iWpRQi68jpvVjz3AwIAkzD8FFhnDBgwQAcPHtSBAwf0y1/+Ug6HQyUlJZ0eLycnR3V1dd7txIkTXVcsAkbG6CgV3jtOEQOtPu32UJsK7x3HfYAAwEQMvRFieHi4goODVV1d7dNeXV0tu93e5n5BQUEaPny4JCkxMVFHjhxRfn6+pkyZ4t2vurpaUVHf/EKrrq5WYmJiq+NZrVZZrdZW30PvkjE6ShOHh2vM8rckSS9kTtCkETcy8wMAJmPoDFBISIjGjx8vp9PpbXO73XI6nUpNTW33OG63W42NjZKkuLg42e12nzHr6+u1f//+Do2J3uvqsMNjMADAnAx/FIbD4dDChQuVlJSk5ORkFRQUqKGhQZmZmZKkBQsWaMiQIcrPz5d0eb1OUlKShg0bpsbGRu3cuVObN29WYWGhJMlisWjJkiVatWqVRowYobi4OC1btkzR0dGaM2eOUYcJAAACiOEBaO7cuaqtrVVubq5cLpcSExNVXFzsXcR8/PhxBQV9M1HV0NCgBx98UCdPnlSfPn00atQovfzyy5o7d663z2OPPaaGhgY98MADOnfunG6//XYVFxfLZuMSZwAAEAABSJKys7OVnZ3d6nv/vLh51apVWrVq1TXHs1gsevLJJ/Xkk092VYkAAKAX6ZFXgQEAAFwPAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADCdG4wuAPg2zW6PyirPqub8JUUMsCk5LkzBQRajywIA9GAEIAS04kNVWrHjsKrqLnnbokJtypsVr4zRUQZWBgDoyTgFhoBVfKhKWS9X+IQfSXLVXVLWyxUqPlRlUGUAgJ6OAISA1Oz2aMWOw/K08t6VthU7DqvZ3VoPAACujQCEgFRWebbFzM/VPJKq6i6prPKs/4oCAPQaBCAEpJrzbYefzvQDAOBqLIJGQIoYYOvSfkBHcOUh0PsRgBCQkuPCFBVqk6vuUqvrgCyS7KGXfzEBXYkrDwFz4BQYAlJwkEV5s+IlXQ47V7vyOm9WPH8rR5fiykPAPAhACFgZo6NUeO84RQy0+rTbQ20qvHccfxtHl+LKQ8BcOAWGgJYxOkoTh4drzPK3JEkvZE7QpBE3MvODLteRKw9Th33Xf4UB6BbMACHgXR12WIyK7sKVh4C5EIAAQFx5CJgNAQgA9M2Vh23NL1p0+WowrjwEegcCEACIKw8BsyEAAcD/x5WHgHlwFRgAXIUrDwFzYAYIAP4JVx4CvR8BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmE5ABKANGzYoNjZWNptNKSkpKisra7Pvxo0bNWnSJA0ePFiDBw9WWlpai/733XefLBaLz5aRkdHdhwEAAHoIwwPQtm3b5HA4lJeXp4qKCo0dO1bp6emqqalptX9JSYnmz5+vPXv2qLS0VDExMZoxY4a+/PJLn34ZGRmqqqrybq+++qo/DgcAAPQAhgegdevWadGiRcrMzFR8fLyKiorUt29fbdq0qdX+r7zyih588EElJiZq1KhR+t3vfie32y2n0+nTz2q1ym63e7fBgwf743AAAEAPYGgAampqUnl5udLS0rxtQUFBSktLU2lpabvGuHjxor766iuFhfk+obmkpEQREREaOXKksrKydObMmTbHaGxsVH19vc8GAAB6L0MD0OnTp9Xc3KzIyEif9sjISLlcrnaN8fOf/1zR0dE+ISojI0MvvfSSnE6nnnrqKe3du1czZ85Uc3Nzq2Pk5+crNDTUu8XExHT+oAAAQMDr0Q9DXb16tbZu3aqSkhLZbDZv+7x587w/jxkzRgkJCRo2bJhKSko0ffr0FuPk5OTI4XB4X9fX1xOCAADoxQydAQoPD1dwcLCqq6t92qurq2W326+579q1a7V69Wq99dZbSkhIuGbfoUOHKjw8XEePHm31favVqoEDB/psAACg9zI0AIWEhGj8+PE+C5ivLGhOTU1tc781a9Zo5cqVKi4uVlJS0rd+zsmTJ3XmzBlFRUV1Sd0AAKBnM/wqMIfDoY0bN+rFF1/UkSNHlJWVpYaGBmVmZkqSFixYoJycHG//p556SsuWLdOmTZsUGxsrl8sll8ulCxcuSJIuXLigRx99VO+//76OHTsmp9Op2bNna/jw4UpPTzfkGAEAQGAxfA3Q3LlzVVtbq9zcXLlcLiUmJqq4uNi7MPr48eMKCvompxUWFqqpqUk/+tGPfMbJy8vT8uXLFRwcrI8++kgvvviizp07p+joaM2YMUMrV66U1Wr167EBAIDAZHgAkqTs7GxlZ2e3+l5JSYnP62PHjl1zrD59+mjXrl1dVBkAAOiNDD8FBgAA4G8EIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoBEYA2bNig2NhY2Ww2paSkqKysrM2+Gzdu1KRJkzR48GANHjxYaWlpLfp7PB7l5uYqKipKffr0UVpamj7//PPuPgwAANBDGB6Atm3bJofDoby8PFVUVGjs2LFKT09XTU1Nq/1LSko0f/587dmzR6WlpYqJidGMGTP05ZdfevusWbNGzz77rIqKirR//37169dP6enpunTpkr8OCwAABDDDA9C6deu0aNEiZWZmKj4+XkVFRerbt682bdrUav9XXnlFDz74oBITEzVq1Cj97ne/k9vtltPplHR59qegoEBLly7V7NmzlZCQoJdeekmnTp3S9u3b/XhkAAAgUBkagJqamlReXq60tDRvW1BQkNLS0lRaWtquMS5evKivvvpKYWFhkqTKykq5XC6fMUNDQ5WSktLmmI2Njaqvr/fZAABA72VoADp9+rSam5sVGRnp0x4ZGSmXy9WuMX7+858rOjraG3iu7NeRMfPz8xUaGurdYmJiOnooAACgBzH8FNj1WL16tbZu3ao33nhDNput0+Pk5OSorq7Ou504caILqwQAAIHmBiM/PDw8XMHBwaqurvZpr66ult1uv+a+a9eu1erVq/WXv/xFCQkJ3vYr+1VXVysqKspnzMTExFbHslqtslqtnTwKAADQ0xg6AxQSEqLx48d7FzBL8i5oTk1NbXO/NWvWaOXKlSouLlZSUpLPe3FxcbLb7T5j1tfXa//+/dccEwAAmIehM0CS5HA4tHDhQiUlJSk5OVkFBQVqaGhQZmamJGnBggUaMmSI8vPzJUlPPfWUcnNztWXLFsXGxnrX9fTv31/9+/eXxWLRkiVLtGrVKo0YMUJxcXFatmyZoqOjNWfOHKMOEwAABBDDA9DcuXNVW1ur3NxcuVwuJSYmqri42LuI+fjx4woK+maiqrCwUE1NTfrRj37kM05eXp6WL18uSXrsscfU0NCgBx54QOfOndPtt9+u4uLi61onBAAAeg/DA5AkZWdnKzs7u9X3SkpKfF4fO3bsW8ezWCx68skn9eSTT3ZBdQAAoLfp0VeBAQAAdAYBCAAAmE6HA9Dx48fl8XhatHs8Hh0/frxLigIAAOhOHQ5AcXFxqq2tbdF+9uxZxcXFdUlRAAAA3anDAcjj8chisbRov3DhAldZAQCAHqHdV4E5HA5Jl6+wWrZsmfr27et9r7m5Wfv372/zTssAAACBpN0B6MMPP5R0eQbo448/VkhIiPe9kJAQjR07Vo888kjXVwgAANDF2h2A9uzZI0nKzMzUM888o4EDB3ZbUQAAAN2pwzdCfP7557ujDgAAAL/pcABqaGjQ6tWr5XQ6VVNTI7fb7fP+3/72ty4rDgAAoDt0OAD927/9m/bu3auf/OQnioqKavWKMAAAgEDW4QD05z//WX/60580ceLE7qgHAACg23X4PkCDBw9WWFhYd9QCAADgFx0OQCtXrlRubq4uXrzYHfUAAAB0u3adArv11lt91vocPXpUkZGRio2N1Xe+8x2fvhUVFV1bIQAAQBdrVwCaM2dON5cBAADgP+0KQHl5ed1dBwAAgN90eA0QAABAT9fhy+AHDx7c6r1/LBaLbDabhg8frvvuu0+ZmZldUiAAAEBX63AAys3N1S9/+UvNnDlTycnJkqSysjIVFxdr8eLFqqysVFZWlr7++mstWrSoywsGAAC4Xh0OQPv27dOqVav0s5/9zKf9ueee01tvvaU//OEPSkhI0LPPPksAAgAAAanDa4B27dqltLS0Fu3Tp0/Xrl27JEk//OEPeSYYAAAIWB0OQGFhYdqxY0eL9h07dnjvEN3Q0KABAwZcf3UAAADdoMOnwJYtW6asrCzt2bPHuwbowIED2rlzp4qKiiRJu3fv1uTJk7u2UgAAgC7S4QC0aNEixcfHa/369Xr99dclSSNHjtTevXt12223SZL+4z/+o2urBAAA6EIdDkCSNHHiRJ4GDwAAeqx2BaD6+noNHDjQ+/O1XOkHAAAQqNoVgAYPHqyqqipFRERo0KBBrd4I0ePxyGKxqLm5ucuLBAAA6ErtCkBvv/229wqvPXv2dGtBAAAA3a1dAejqK7q4ugsAAPR0nXoY6rvvvqt7771Xt912m7788ktJ0ubNm7Vv374uLQ4AAKA7dDgA/eEPf1B6err69OmjiooKNTY2SpLq6ur0q1/9qssLBAAA6GodDkCrVq1SUVGRNm7cqO985zve9okTJ6qioqJLiwMAAOgOHQ5An332me64444W7aGhoTp37lxX1AQAANCtOhyA7Ha7jh492qJ93759Gjp0aJcUBQAA0J06HIAWLVqkhx9+WPv375fFYtGpU6f0yiuv6JFHHlFWVlZ31AgAANCl2v0ojMrKSsXFxenxxx+X2+3W9OnTdfHiRd1xxx2yWq165JFH9NBDD3VnrQAAAF2i3QFo2LBhuuWWWzR16lRNnTpVR44c0fnz53XhwgXFx8erf//+3VknAABAl2l3AHr77bdVUlKikpISvfrqq2pqatLQoUM1bdo0TZs2TVOmTFFkZGR31goAANAl2r0GaMqUKVq+fLlKSkr097//Xbt379b8+fN15MgRLVy4UNHR0frBD37Q4QI2bNig2NhY2Ww2paSkqKysrM2+n3zyie6++27FxsbKYrGooKCgRZ/ly5fLYrH4bKNGjepwXQAAoPfq1J2gbTabpk2bpqVLl2rFihX693//d/Xv31+ffvpph8bZtm2bHA6H8vLyVFFRobFjxyo9PV01NTWt9r948aKGDh2q1atXy263tznuD37wA1VVVXk37lANAACu1qEA1NTUpHfeeUcrVqzQ1KlTNWjQIP3sZz/T3//+d61fv16VlZUd+vB169Zp0aJFyszMVHx8vIqKitS3b19t2rSp1f4TJkzQ008/rXnz5slqtbY57g033CC73e7dwsPDO1QXAADo3dq9BmjatGnav3+/4uLiNHnyZP30pz/Vli1bFBUV1akPbmpqUnl5uXJycrxtQUFBSktLU2lpaafGvOLzzz9XdHS0bDabUlNTlZ+fr5tvvrnN/o2Njd5HekhSfX39dX0+AAAIbO2eAXr33Xf13e9+V9OmTdP06dN15513djr8SNLp06fV3NzcYuF0ZGSkXC5Xp8dNSUnRCy+8oOLiYhUWFqqyslKTJk3S+fPn29wnPz9foaGh3i0mJqbTnw8AAAJfuwPQuXPn9Nvf/lZ9+/bVU089pejoaI0ZM0bZ2dl67bXXVFtb2511ttvMmTN1zz33KCEhQenp6dq5c6fOnTun//7v/25zn5ycHNXV1Xm3EydO+LFiAADgb+0+BdavXz9lZGQoIyNDknT+/Hnt27dPe/bs0Zo1a/TjH/9YI0aM0KFDh9o1Xnh4uIKDg1VdXe3TXl1dfc0Fzh01aNAgfe9732v18R1XWK3Wa64pAgAAvUunrgKTLgeisLAwhYWFafDgwbrhhht05MiRdu8fEhKi8ePHy+l0etvcbrecTqdSU1M7W1YLFy5c0BdffHFdp+sAAEDv0u4ZILfbrQ8++EAlJSXas2eP/ud//kcNDQ0aMmSIpk6dqg0bNmjq1Kkd+nCHw6GFCxcqKSlJycnJKigoUENDgzIzMyVJCxYs0JAhQ5Sfny/p8sLpw4cPe3/+8ssvdfDgQfXv31/Dhw+XJD3yyCOaNWuWbrnlFp06dUp5eXkKDg7W/PnzO1QbAADovdodgAYNGqSGhgbZ7XZNnTpV//mf/6kpU6Zo2LBhnf7wuXPnqra2Vrm5uXK5XEpMTFRxcbF3YfTx48cVFPTNJNWpU6d06623el+vXbtWa9eu1eTJk1VSUiJJOnnypObPn68zZ87oxhtv1O233673339fN954Y6frBAAAvUu7A9DTTz+tqVOn6nvf+16XFpCdna3s7OxW37sSaq6IjY2Vx+O55nhbt27tqtIAAEAv1e4A9NOf/rQ76wAAAPCbdgcgXL9mt0dllWdVc/6SIgbYlBwXpuAgi9FlAQBgOgQgPyk+VKUVOw6rqu6Sty0q1Ka8WfHKGM0VagAA+FOnL4NH+xUfqlLWyxU+4UeSXHWXlPVyhYoPVRlUGQAA5kQA6mbNbo9W7Dis1pZuX2lbseOwmt3XXtwNAAC6DgGom5VVnm0x83M1j6Squksqqzzrv6IAADA5AlA3qznfdvjpTD8AAHD9CEDdLGKArUv7AQCA60cA6mbJcWGKCrWprYvdLbp8NVhyXJg/ywIAwNQIQN0sOMiivFnxktQiBF15nTcrnvsBAQDgRwQgP8gYHaXCe8cpYqDVp90ealPhveO4DxAAAH7GjRD9JGN0lCYOD9eY5W9Jkl7InKBJI25k5gcAAAMwA+RHV4cdHoMBAIBxCEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0DA9AGzZsUGxsrGw2m1JSUlRWVtZm308++UR33323YmNjZbFYVFBQcN1jAgAA8zE0AG3btk0Oh0N5eXmqqKjQ2LFjlZ6erpqamlb7X7x4UUOHDtXq1atlt9u7ZEwAAGA+hgagdevWadGiRcrMzFR8fLyKiorUt29fbdq0qdX+EyZM0NNPP6158+bJarV2yZiS1NjYqPr6ep8NAAD0XoYFoKamJpWXlystLe2bYoKClJaWptLSUr+OmZ+fr9DQUO8WExPTqc8HAAA9g2EB6PTp02publZkZKRPe2RkpFwul1/HzMnJUV1dnXc7ceJEpz4fAAD0DDcYXUAgsFqtbZ5SAwAAvY9hM0Dh4eEKDg5WdXW1T3t1dXWbC5yNGBMAAPQ+hgWgkJAQjR8/Xk6n09vmdrvldDqVmpoaMGMCAIDex9BTYA6HQwsXLlRSUpKSk5NVUFCghoYGZWZmSpIWLFigIUOGKD8/X9LlRc6HDx/2/vzll1/q4MGD6t+/v4YPH96uMQEAAAwNQHPnzlVtba1yc3PlcrmUmJio4uJi7yLm48ePKyjom0mqU6dO6dZbb/W+Xrt2rdauXavJkyerpKSkXWMCAAAYvgg6Oztb2dnZrb53JdRcERsbK4/Hc11jAgAAGP4oDAAAAH8jAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMJiAC0YcMGxcbGymazKSUlRWVlZdfs//vf/16jRo2SzWbTmDFjtHPnTp/377vvPlksFp8tIyOjOw8BAAD0IIYHoG3btsnhcCgvL08VFRUaO3as0tPTVVNT02r/9957T/Pnz9f999+vDz/8UHPmzNGcOXN06NAhn34ZGRmqqqrybq+++qo/DgcAAPQAhgegdevWadGiRcrMzFR8fLyKiorUt29fbdq0qdX+zzzzjDIyMvToo4/q+9//vlauXKlx48Zp/fr1Pv2sVqvsdrt3Gzx4sD8OBwAA9ACGBqCmpiaVl5crLS3N2xYUFKS0tDSVlpa2uk9paalPf0lKT09v0b+kpEQREREaOXKksrKydObMmTbraGxsVH19vc8GAAB6L0MD0OnTp9Xc3KzIyEif9sjISLlcrlb3cblc39o/IyNDL730kpxOp5566int3btXM2fOVHNzc6tj5ufnKzQ01LvFxMRc55EBAIBAdoPRBXSHefPmeX8eM2aMEhISNGzYMJWUlGj69Okt+ufk5MjhcHhf19fXE4IAAOjFDJ0BCg8PV3BwsKqrq33aq6urZbfbW93Hbrd3qL8kDR06VOHh4Tp69Gir71utVg0cONBnAwAAvZehASgkJETjx4+X0+n0trndbjmdTqWmpra6T2pqqk9/Sdq9e3eb/SXp5MmTOnPmjKKiorqmcAAA0KMZfhWYw+HQxo0b9eKLL+rIkSPKyspSQ0ODMjMzJUkLFixQTk6Ot//DDz+s4uJi/frXv9ann36q5cuX64MPPlB2drYk6cKFC3r00Uf1/vvv69ixY3I6nZo9e7aGDx+u9PR0Q44RAAAEFsPXAM2dO1e1tbXKzc2Vy+VSYmKiiouLvQudjx8/rqCgb3Labbfdpi1btmjp0qV64oknNGLECG3fvl2jR4+WJAUHB+ujjz7Siy++qHPnzik6OlozZszQypUrZbVaDTlGAAAQWAwPQJKUnZ3tncH5ZyUlJS3a7rnnHt1zzz2t9u/Tp4927drVleUBAIBexvBTYAAAAP5GAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZzg9EFAAAAc2h2e1RWeVY15y8pYoBNyXFhCg6yGFILAQgAAHS74kNVWrHjsKrqLnnbokJtypsVr4zRUX6vh1NgAACgWxUfqlLWyxU+4UeSXHWXlPVyhYoPVfm9JgIQAADoNs1uj1bsOCxPK+9daVux47Ca3a316D4EIAAA0G3KKs+2mPm5mkdSVd0llVWe9V9RIgABAIBuVHO+7fDTmX5dhQAEAAC6TcQAW5f26yoEIAAA0G2S48IUFWpTWxe7W3T5arDkuDB/lkUAAgAA3Sc4yKK8WfGS1CIEXXmdNyve7/cDIgABAIBulTE6SoX3jlPEQKtPuz3UpsJ7xxlyHyBuhAgAALpdxugoTRwerjHL35IkvZA5QZNG3GjYnaCZAQIAAH5xddgx8jEYEgEIAACYUEAEoA0bNig2NlY2m00pKSkqKyu7Zv/f//73GjVqlGw2m8aMGaOdO3f6vO/xeJSbm6uoqCj16dNHaWlp+vzzz7vzEAAAQA9ieADatm2bHA6H8vLyVFFRobFjxyo9PV01NTWt9n/vvfc0f/583X///frwww81Z84czZkzR4cOHfL2WbNmjZ599lkVFRVp//796tevn9LT03Xpkn9vsgQAAAKT4Yug161bp0WLFikzM1OSVFRUpD/96U/atGmTHn/88Rb9n3nmGWVkZOjRRx+VJK1cuVK7d+/W+vXrVVRUJI/Ho4KCAi1dulSzZ8+WJL300kuKjIzU9u3bNW/ePP8d3D/xeDyyft0oSXJfvCj314b/6+8R3E1fd+m/t64cr71jef7R6O3n+cc/5A5yd3ttXS2Qa+tqZjpWwJ+u/n/L4/Hvs7/+mcVjYAVNTU3q27evXnvtNc2ZM8fbvnDhQp07d05//OMfW+xz8803y+FwaMmSJd62vLw8bd++XX/961/1t7/9TcOGDdOHH36oxMREb5/JkycrMTFRzzzzTIsxGxsb1djY6H1dX1+vmJgY1dXVaeDAgV1yrJJ04Vy9TvyflC4bDwCAnirm/f3qP6jrfsd2lKGnwE6fPq3m5mZFRkb6tEdGRsrlcrW6j8vlumb/K//syJj5+fkKDQ31bjExMZ06nm/TN4S/RQIAIBn/O5HfyJJycnLkcDi8r6/MAHU1S58+GllR3uXjAgDQ01j69DH08w0NQOHh4QoODlZ1dbVPe3V1tex2e6v72O32a/a/8s/q6mpFRUX59Ln6lNjVrFarrFZrq+91JYvFIkvfvt3+OQAA4NoMPQUWEhKi8ePHy+l0etvcbrecTqdSU1Nb3Sc1NdWnvyTt3r3b2z8uLk52u92nT319vfbv39/mmAAAwFwMPwXmcDi0cOFCJSUlKTk5WQUFBWpoaPBeFbZgwQINGTJE+fn5kqSHH35YkydP1q9//Wvddddd2rp1qz744AP99re/lXR5lmXJkiVatWqVRowYobi4OC1btkzR0dE+C60BAIB5GR6A5s6dq9raWuXm5srlcikxMVHFxcXeRczHjx9XUNA3E1W33XabtmzZoqVLl+qJJ57QiBEjtH37do0ePdrb57HHHlNDQ4MeeOABnTt3TrfffruKi4tls9n8fnwAACDwGHoZfKCqr69XaGhol18GDwAAAoPhd4IGAADwNwIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHQIQAAAwHcMfhRGIrtwcu76+3uBKAABARw0YMEAWi+WafQhArTh//rwkKSYmxuBKAABAR7XnUVY8C6wVbrdbp06daleC7Kj6+nrFxMToxIkTPGfMQHwPgYHvITDwPQQGvoeuwwxQJwUFBemmm27q1s8YOHAg/4EHAL6HwMD3EBj4HgID34N/sAgaAACYDgEIAACYDgHIz6xWq/Ly8mS1Wo0uxdT4HgID30Ng4HsIDHwP/sUiaAAAYDrMAAEAANMhAAEAANMhAAEAANMhAAEAANMhAPnRhg0bFBsbK5vNppSUFJWVlRldkum88847mjVrlqKjo2WxWLR9+3ajSzKd/Px8TZgwQQMGDFBERITmzJmjzz77zOiyTKewsFAJCQnem+6lpqbqz3/+s9Flmd7q1atlsVi0ZMkSo0vp9QhAfrJt2zY5HA7l5eWpoqJCY8eOVXp6umpqaowuzVQaGho0duxYbdiwwehSTGvv3r1avHix3n//fe3evVtfffWVZsyYoYaGBqNLM5WbbrpJq1evVnl5uT744ANNmzZNs2fP1ieffGJ0aaZ14MABPffcc0pISDC6FFPgMng/SUlJ0YQJE7R+/XpJl583FhMTo4ceekiPP/64wdWZk8Vi0RtvvKE5c+YYXYqp1dbWKiIiQnv37tUdd9xhdDmmFhYWpqefflr333+/0aWYzoULFzRu3Dj95je/0apVq5SYmKiCggKjy+rVmAHyg6amJpWXlystLc3bFhQUpLS0NJWWlhpYGWC8uro6SZd/+cIYzc3N2rp1qxoaGpSammp0Oaa0ePFi3XXXXT6/J9C9eBiqH5w+fVrNzc2KjIz0aY+MjNSnn35qUFWA8dxut5YsWaKJEydq9OjRRpdjOh9//LFSU1N16dIl9e/fX2+88Ybi4+ONLst0tm7dqoqKCh04cMDoUkyFAATAMIsXL9ahQ4e0b98+o0sxpZEjR+rgwYOqq6vTa6+9poULF2rv3r2EID86ceKEHn74Ye3evVs2m83ockyFAOQH4eHhCg4OVnV1tU97dXW17Ha7QVUBxsrOztabb76pd955RzfddJPR5ZhSSEiIhg8fLkkaP368Dhw4oGeeeUbPPfecwZWZR3l5uWpqajRu3DhvW3Nzs9555x2tX79ejY2NCg4ONrDC3os1QH4QEhKi8ePHy+l0etvcbrecTifn22E6Ho9H2dnZeuONN/T2228rLi7O6JLw/7ndbjU2NhpdhqlMnz5dH3/8sQ4ePOjdkpKS9OMf/1gHDx4k/HQjZoD8xOFwaOHChUpKSlJycrIKCgrU0NCgzMxMo0szlQsXLujo0aPe15WVlTp48KDCwsJ08803G1iZeSxevFhbtmzRH//4Rw0YMEAul0uSFBoaqj59+hhcnXnk5ORo5syZuvnmm3X+/Hlt2bJFJSUl2rVrl9GlmcqAAQNarH/r16+fvvvd77IurpsRgPxk7ty5qq2tVW5urlwulxITE1VcXNxiYTS61wcffKCpU6d6XzscDknSwoUL9cILLxhUlbkUFhZKkqZMmeLT/vzzz+u+++7zf0EmVVNTowULFqiqqkqhoaFKSEjQrl27dOeddxpdGuAX3AcIAACYDmuAAACA6RCAAACA6RCAAACA6RCAAACA6RCAAACA6RCAAACA6RCAAACA6RCAAACA6RCAAACA6RCAAACA6RCAAACA6RCAAPR6tbW1stvt+tWvfuVte++99xQSEiKn02lgZQCMwsNQAZjCzp07NWfOHL333nsaOXKkEhMTNXv2bK1bt87o0gAYgAAEwDQWL16sv/zlL0pKStLHH3+sAwcOyGq1Gl0WAAMQgACYxj/+8Q+NHj1aJ06cUHl5ucaMGWN0SQAMwhogAKbxxRdf6NSpU3K73Tp27JjR5QAwEDNAAEyhqalJycnJSkxM1MiRI1VQUKCPP/5YERERRpcGwAAEIACm8Oijj+q1117TX//6V/Xv31+TJ09WaGio3nzzTaNLA2AAToEB6PVKSkpUUFCgzZs3a+DAgQoKCtLmzZv17rvvqrCw0OjyABiAGSAAAGA6zAABAADTIQABAADTIQABAADTIQABAADTIQABAADTIQABAADTIQABAADTIQABAADTIQABAADTIQABAADTIQABAADT+X9ZQ9o42HBp7wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the weight for each observation\n",
    "ax = plt.subplot(111)\n",
    "\n",
    "ax.stem(Dat, Weights)\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('Weight')\n",
    "\n",
    "ax.spines[['right', 'top']].set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "WeightedMean = np.average(Dat, weights=Weights)\n",
    "TrimmedMean = trim_mean(Dat, 0.1)\n",
    "InterquartileMean = trim_mean(Dat, 0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "masked_array(data=[0.8, 0.8, 1. , 1.2, 1.3, 1.3, 1.4, 1.8, 2.4, 2.4],\n",
       "             mask=False,\n",
       "       fill_value=1e+20)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DatWinsorized = winsorize(Dat, [0.1,0.1]) # Winsorized the data\n",
    "DatWinsorized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "WinsorizedMean = DatWinsorized.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "CentralTendency = [Mode, Median, Mean, GeomMean, HarmMean, WeightedMean, TrimmedMean, InterquartileMean, WinsorizedMean]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(1.75, 0.75, 'Central tendency statistics')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxoAAAFfCAYAAADJSrqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAi4UlEQVR4nO3deXRW9Z3A4W8Ak0ASQGUX2a0KiKmIVq2giAsF6lZRaxWhtpyKCqPFVh13e6jLVCwqamdGrEp1RHGpg8iIyCh2hqqMihtqcBlRBCQQUBFy5w8PGWNAAv7wJeF5zsmR3Pfm3t99740nn9wleVmWZQEAAJBQg1wPAAAAqH+EBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQGwnZo0aVLk5eVVfRQWFka7du3iyCOPjD/+8Y+xcuXKLVrunDlz4rLLLovly5enHTAAdYrQANjOXXHFFXHnnXfGxIkT4+yzz46IiDFjxsRee+0VL7744mYvb86cOXH55ZcLDYDtXKNcDwCA3Bo4cGDsu+++VZ9fcMEFMXPmzBg8eHD8+Mc/jldffTUaN26cwxECUBc5owFADf3794+LL7443nnnnbjrrrsiIuLFF1+M008/Pbp06RKFhYXRpk2bGDFiRCxdurTq6y677LIYO3ZsRER07ty56rKshQsXRkTE7bffHv37949WrVpFQUFBdO/ePSZOnPidbx8AW58zGgBs0KmnnhoXXnhhPP744/GLX/wiZsyYEW+//XYMHz482rRpE/Pnz4/bbrst5s+fH3/7298iLy8vjjvuuHjjjTfiL3/5S1x//fXRokWLiIho2bJlRERMnDgxevToET/+8Y+jUaNG8cgjj8SZZ54ZlZWVMWrUqFxuLgCJ5WVZluV6EAB89yZNmhTDhw+PuXPnVrt06quaN28eXbp0ieeffz4+/fTTGpdQ3XPPPXHyySfH7Nmz4+CDD46IiOuuuy7Gjh0bZWVl0alTp2rzb2gZRx11VCxYsCDeeuutdBsHQM65dAqAjSouLq56+tRXA+Gzzz6LJUuWxA9+8IOIiHj++edrtbyvLqO8vDyWLFkS/fr1i7fffjvKy8sTjhyAXBMaAGxURUVFlJSURETEsmXLYvTo0dG6deto3LhxtGzZMjp37hwRUetIeOaZZ2LAgAFRVFQUzZs3j5YtW8aFF164WcsAoG5wjwYAG/T+++9HeXl5dOvWLSIihg4dGnPmzImxY8dGaWlpFBcXR2VlZRx11FFRWVm5yeW99dZbcdhhh8Uee+wRf/jDH2LXXXeN/Pz8+Pd///e4/vrra7UMAOoOoQHABt15550REXHkkUfGJ598Ek888URcfvnlcckll1TNs2DBghpfl5eXt8HlPfLII/H555/Hww8/HB06dKia/uSTTyYeOQDbApdOAVDDzJkz48orr4zOnTvHKaecEg0bNoyIiK8/P2T8+PE1vraoqCgiosYf7NvQMsrLy+P2229POHIAthXOaABs56ZNmxavvfZarF27Nj766KOYOXNmzJgxIzp27BgPP/xwFBYWRmFhYfTt2zeuueaa+OKLL2KXXXaJxx9/PMrKymosr3fv3hERcdFFF8VJJ50UO+ywQwwZMiSOOOKIyM/PjyFDhsTIkSOjoqIi/vSnP0WrVq1i0aJF3/VmA7CVCQ2A7dz6S6Hy8/Njp512ir322ivGjx8fw4cPr7oRPCJi8uTJcfbZZ8dNN90UWZbFEUccEdOmTYt27dpVW16fPn3iyiuvjFtuuSUee+yxqKysjLKysth9991jypQp8Y//+I/x61//Otq0aRO/+tWvomXLljFixIjvdJsB2Pr8HQ0AACA592gAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJNcr1ADZLlkWsXv3lv5s0icjLy+142Hpyta9zeYx9m3VXVkYsWfLlv1u0iGiwmb9D2B6/t7bHbc4l7zfAlquj/w+tW2c0Vq+OKC7+8mP9m039lKt9nctj7Nuse8mSiNatv/xYHxzf1brrqu1xm3PJ+w2w5ero/0PrVmgAAAB1gtAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByjWozU5ZlsXLlyq09lk1bter//71iRcS6dbkbC1tXrvZ1Lo+xb7Pur35/rlwZUVj43a27rtoetzmXvN8AW24b/X9oSUlJ5OXlbfT1vCzLsk0tZMWKFdGsWbOkAwMAAOqu8vLyaNq06UZfr1VobDNnNOLL6Nl1113jvffe+8YNgy3lGOO74Dhja3OMsbU5xtjUGY1aXTqVl5e3zR1ATZs23ebGRP3iGOO74Dhja3OMsbU5xtgYN4MDAADJCQ0AACC5OhcaBQUFcemll0ZBQUGuh0I95Rjju+A4Y2tzjLG1OcbYlFrdDA4AALA56twZDQAAYNsnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAAAAkhMaAABAckIDAABITmgAAADJCQ0AACA5oQEAACQnNAAAgOSEBgAAkJzQAKjjTj/99OjUqVOuh1ErdWms26qFCxdGXl5eTJo0KcnyZs2aFXl5eTFr1qwkywNYT2gA24233norRo4cGV26dInCwsJo2rRpHHTQQXHDDTfEp59+utXW+8EHH8Rll10W8+bN22rrqI2bb7452Q+nfLNXXnklLrvssli4cOEWL2Py5Mkxfvz4ZGOy/4HvWqNcDwDgu/Doo4/GCSecEAUFBXHaaadFz549Y82aNfH000/H2LFjY/78+XHbbbdtlXV/8MEHcfnll0enTp2itLR0q6yjNm6++eZo0aJFnH766Tkbw/bilVdeicsvvzwOOeSQLT6DM3ny5Hj55ZdjzJgx1aZ37NgxPv3009hhhx02a3kb2/99+/aNTz/9NPLz87donAAbIzSAeq+srCxOOumk6NixY8ycOTPatm1b9dqoUaPizTffjEcffTSHI6xu9erV0aRJk1wPg21UXl5eFBYWJltegwYNki4PYD2XTgH13jXXXBMVFRXxL//yL9UiY71u3brF6NGjq0276667onfv3tG4cePYaaed4qSTTor33nuv2jyHHHJI9OzZM1555ZU49NBDo0mTJrHLLrvENddcUzXPrFmzok+fPhERMXz48MjLy6t2ff36ZTz33HPRt2/faNKkSVx44YUREfHQQw/FoEGDol27dlFQUBBdu3aNK6+8MtatW7fZ70GnTp1i/vz58dRTT1WN4ZBDDql6ffny5TFmzJjYddddo6CgILp16xZXX311VFZWVs2z/t6A6667Lm677bbo2rVrFBQURJ8+fWLu3Lk11vnggw9Gz549o7CwMHr27BlTp07d4NgqKytj/Pjx0aNHjygsLIzWrVvHyJEj45NPPqmxDYMHD46nn3469ttvvygsLIwuXbrEn//85xrLXL58efzDP/xDdOrUKQoKCqJ9+/Zx2mmnxZIlS6KioiKKiopq7POIiPfffz8aNmwY48aN+8b385577onevXtHSUlJNG3aNPbaa6+44YYbIiJi0qRJccIJJ0RExKGHHlr1fq+/B6I2+/WQQw6JRx99NN55552qr19/ZmRD92h8+OGHMXz48Gjfvn0UFBRE27Zt4+ijj666dOub9v/G7tH4r//6r/jRj34UO+64YxQVFUWvXr2qtrE26wRwRgOo9x555JHo0qVLHHjggbWa/3e/+11cfPHFMXTo0DjjjDPi448/jgkTJkTfvn3jhRdeiObNm1fN+8knn8RRRx0Vxx13XAwdOjSmTJkSv/nNb2KvvfaKgQMHxp577hlXXHFFXHLJJfHLX/4yDj744IiIamNZunRpDBw4ME466aT42c9+Fq1bt46IL39gLS4ujnPPPTeKi4tj5syZcckll8SKFSvi2muv3az3YPz48XH22WdHcXFxXHTRRRERVetZvXp19OvXL/73f/83Ro4cGR06dIg5c+bEBRdcEIsWLapxn8DkyZNj5cqVMXLkyMjLy4trrrkmjjvuuHj77berLud5/PHH4/jjj4/u3bvHuHHjYunSpVU/lH7dyJEjY9KkSTF8+PA455xzoqysLG688cZ44YUX4plnnql2idCbb74ZP/nJT+LnP/95DBs2LP71X/81Tj/99Ojdu3f06NEjIiIqKiri4IMPjldffTVGjBgR++yzTyxZsiQefvjheP/996O0tDSOPfbYuPfee+MPf/hDNGzYsGr5f/nLXyLLsjjllFM2+l7OmDEjTj755DjssMPi6quvjoiIV199NZ555pkYPXp09O3bN84555z44x//GBdeeGHsueeeERFV/63Nfr3ooouivLw83n///bj++usjIqK4uHijYzr++ONj/vz5cfbZZ0enTp1i8eLFMWPGjHj33XejU6dO37j/N7aNgwcPjrZt28bo0aOjTZs28eqrr8Zf//rXqkDb1DoBIgOox8rLy7OIyI4++uhazb9w4cKsYcOG2e9+97tq01966aWsUaNG1ab369cvi4jsz3/+c9W0zz//PGvTpk12/PHHV02bO3duFhHZ7bffXmN965dxyy231Hht9erVNaaNHDkya9KkSfbZZ59VTRs2bFjWsWPHTW5bjx49sn79+tWYfuWVV2ZFRUXZG2+8UW36b3/726xhw4bZu+++m2VZlpWVlWURke28887ZsmXLquZ76KGHsojIHnnkkapppaWlWdu2bbPly5dXTXv88ceziKg21v/8z//MIiK7++67q637scceqzG9Y8eOWURks2fPrpq2ePHirKCgIDvvvPOqpl1yySVZRGQPPPBAjW2trKzMsizLpk+fnkVENm3atGqv9+rVa4Pv0VeNHj06a9q0abZ27dqNznPfffdlEZE9+eSTNV6r7X4dNGjQBvfr+v2w/nj65JNPsojIrr322m8c98b2/5NPPlltrGvXrs06d+6cdezYMfvkk0+qzbv+/avtOoHtm0ungHptxYoVERFRUlJSq/kfeOCBqKysjKFDh8aSJUuqPtq0aRO77bZbPPnkk9XmLy4ujp/97GdVn+fn58d+++0Xb7/9dq3HWFBQEMOHD68xvXHjxlX/XrlyZSxZsiQOPvjgWL16dbz22mu1Xv6m3HfffXHwwQfHjjvuWG2bBwwYEOvWrYvZs2dXm//EE0+MHXfcserz9Wdp1m/zokWLYt68eTFs2LBo1qxZ1XyHH354dO/evca6mzVrFocffni1dffu3TuKi4trvN/du3evWl9ERMuWLWP33Xev9n7ff//9sffee8exxx5bY1vz8vIiImLAgAHRrl27uPvuu6tee/nll+PFF1+stj83pHnz5rFq1aqYMWPGN863Man3a+PGjSM/Pz9mzZpV43KzLfHCCy9EWVlZjBkzptrZu4j/f/9SrxOon4QGUK81bdo0Ir78ga42FixYEFmWxW677RYtW7as9vHqq6/G4sWLq83fvn37qh++1ttxxx0364evXXbZZYNP/Jk/f34ce+yx0axZs2jatGm0bNmy6ofg8vLyWi9/UxYsWBCPPfZYje0dMGBARESNbe7QoUO1z9dHx/ptfueddyIiYrfddquxrt13373GusvLy6NVq1Y11l9RUbHJda9f/1ff77feeit69uz5jdvcoEGDOOWUU+LBBx+M1atXR0TE3XffHYWFhVX3V2zMmWeeGd/73vdi4MCB0b59+xgxYkQ89thj3/g1X5V6vxYUFMTVV18d06ZNi9atW0ffvn3jmmuuiQ8//HCzlxXx5fsXEd/4HqZeJ1A/uUcDqNeaNm0a7dq1i5dffrlW81dWVkZeXl5Mmzat2rX76339OvkNzRMRkWVZrcf41d9wr7d8+fLo169fNG3aNK644oro2rVrFBYWxvPPPx+/+c1vqt2k/W1VVlbG4YcfHueff/4GX//e975X7fMU2/zVdbdq1aramYWvatmy5VZb92mnnRbXXnttPPjgg3HyySfH5MmTY/DgwdXOwmxIq1atYt68eTF9+vSYNm1aTJs2LW6//fY47bTT4o477vjGr91a+3XMmDExZMiQePDBB2P69Olx8cUXx7hx42LmzJnx/e9/f4uWuS2uE6hbhAZQ7w0ePDhuu+22ePbZZ+OAAw74xnm7du0aWZZF586da/yAvaW+fsajNmbNmhVLly6NBx54IPr27Vs1vaysLPk4unbtGhUVFVVnML6tjh07RsSXZyu+7vXXX6+x7v/4j/+Igw46aIPBtSW6du1aq7Ds2bNnfP/734+777472rdvH++++25MmDChVuvIz8+PIUOGxJAhQ6KysjLOPPPMuPXWW+Piiy+Obt26bfS93pz9urnHTdeuXeO8886L8847LxYsWBClpaXxT//0T3HXXXdt1vK6du0aEV9eSrapY2JT6wS2by6dAuq9888/P4qKiuKMM86Ijz76qMbrb731VtVjO4877rho2LBhXH755TV+S55lWSxdunSz119UVBQRX/42u7bW/+b+q2NYs2ZN3HzzzZu9/q+OY0NjGDp0aDz77LMxffr0Gq8tX7481q5du1nradu2bZSWlsYdd9xR7VKgGTNmxCuvvFJj3evWrYsrr7yyxnLWrl27We/Zescff3z8z//8zwYfp/v1fXrqqafG448/HuPHj4+dd945Bg4cuMnlf/0YaNCgQfTq1SsiIj7//POI2Pg+35z9WlRUVKtLqVavXh2fffZZtWldu3aNkpKSqvGsX15t3s999tknOnfuHOPHj68x//px13adwPbNGQ2g3uvatWtMnjw5TjzxxNhzzz2r/WXwOXPmxH333Vf115K7du0aV111VVxwwQWxcOHCOOaYY6KkpCTKyspi6tSp8ctf/jJ+/etfb/b6mzdvHrfcckuUlJREUVFR7L///tG5c+eNfs2BBx4YO+64YwwbNizOOeecyMvLizvvvHOLLhFar3fv3jFx4sS46qqrolu3btGqVavo379/jB07Nh5++OEYPHhw1aNiV61aFS+99FJMmTIlFi5cGC1atNisdY0bNy4GDRoUP/zhD2PEiBGxbNmymDBhQvTo0SMqKiqq5uvXr1+MHDkyxo0bF/PmzYsjjjgidthhh1iwYEHcd999ccMNN8RPfvKTzVr32LFjY8qUKXHCCSfEiBEjonfv3rFs2bJ4+OGH45Zbbom99967at6f/vSncf7558fUqVPjV7/6Va3+2vYZZ5wRy5Yti/79+0f79u3jnXfeiQkTJkRpaWnVI2xLS0ujYcOGcfXVV0d5eXkUFBRE//79N2u/9u7dO+69994499xzo0+fPlFcXBxDhgypMd8bb7wRhx12WAwdOjS6d+8ejRo1iqlTp8ZHH30UJ510UrXlbWj/f12DBg1i4sSJMWTIkCgtLY3hw4dH27Zt47XXXov58+fH9OnTa71OYDuXk2ddAeTAG2+8kf3iF7/IOnXqlOXn52clJSXZQQcdlE2YMKHaY0WzLMvuv//+7Ic//GFWVFSUFRUVZXvssUc2atSo7PXXX6+ap1+/flmPHj1qrGdDj5t96KGHsu7du2eNGjWq9mjSjS0jy7LsmWeeyX7wgx9kjRs3ztq1a5edf/75VY9l/epjU2v7eNsPP/wwGzRoUFZSUpJFRLVHna5cuTK74IILsm7dumX5+flZixYtsgMPPDC77rrrsjVr1mRZ9v+PVd3QI00jIrv00kurTbv//vuzPffcMysoKMi6d++ePfDAAxsd62233Zb17t07a9y4cVZSUpLttdde2fnnn5998MEHVfN07NgxGzRoUI2v7devX43Hti5dujQ766yzsl122SXLz8/P2rdvnw0bNixbsmRJja//0Y9+lEVENmfOnG949/7flClTsiOOOCJr1apVlp+fn3Xo0CEbOXJktmjRomrz/elPf8q6dOmSNWzYsNo+q+1+raioyH76059mzZs3r/ZY4K8/3nbJkiXZqFGjsj322CMrKirKmjVrlu2///7Zv/3bv1Ubz8b2/9cfb7ve008/nR1++OFZSUlJVlRUlPXq1SubMGHCZq0T2L7lZdm3+PUYANRxxx57bLz00kvx5ptv5nooAPWKezQA2G4tWrQoHn300Tj11FNzPRSAesc9GgBsd8rKyuKZZ56Jf/7nf44ddtghRo4cmeshAdQ7zmgAsN156qmn4tRTT42ysrK44447ok2bNrkeEkC94x4NAAAgOWc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEiuUa4HsDkqKytj2eo1ERGxU5P8aNBAJ0HEhr83siyL1WvWxqdfrIvGOzSMJvmNIi8vb4Nfn2VZfPrFuoiIaLxDw6r5Njb929gaywQAtj11KjSWrV4T+171RERE/P0fD4sWxYU5HhFsGzb0vfHpF+uix6WPV83zyhVHRpP8DX/Lf/rFuuh+yfQa821s+rexNZYJAGx7nBIAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEhOaAAAAMkJDQAAIDmhAQAAJCc0AACA5IQGAACQnNAAAACSExoAAEByQgMAAEiuUW1myrIsVq5cubXHskkrKz6Lys9Xf/nvFSsiv3JNjkcE24YNfW+sXrO2alpExIoVK2Jt/oa/5b8671fn29j0b2NrLBMA+O6VlJREXl7eRl/Py7Is29RCVqxYEc2aNUs6MAAAoO4qLy+Ppk2bbvT1WoXGtnJGI+LL6Nl1113jvffe+8YNgy3lGOO74Dhja3OMsbU5xtjUGY1aXbOQl5e3zR1ATZs23ebGRP3iGOO74Dhja3OMsbU5xtgYN4MDAADJCQ0AACC5OhcaBQUFcemll0ZBQUGuh0I95Rjju+A4Y2tzjLG1OcbYlFrdDA4AALA56twZDQAAYNsnNAAAgOSEBgAAkJzQAAAAkhMaAABAcnUuNG666abo1KlTFBYWxv777x///d//neshUY/Mnj07hgwZEu3atYu8vLx48MEHcz0k6pFx48ZFnz59oqSkJFq1ahXHHHNMvP7667keFvXIxIkTo1evXlV/qfmAAw6IadOm5XpY1GO///3vIy8vL8aMGZProbANqlOhce+998a5554bl156aTz//POx9957x5FHHhmLFy/O9dCoJ1atWhV777133HTTTbkeCvXQU089FaNGjYq//e1vMWPGjPjiiy/iiCOOiFWrVuV6aNQT7du3j9///vfx3HPPxd///vfo379/HH300TF//vxcD416aO7cuXHrrbdGr169cj0UtlF16u9o7L///tGnT5+48cYbIyKisrIydt111zj77LPjt7/9bY5HR32Tl5cXU6dOjWOOOSbXQ6Ge+vjjj6NVq1bx1FNPRd++fXM9HOqpnXbaKa699tr4+c9/nuuhUI9UVFTEPvvsEzfffHNcddVVUVpaGuPHj8/1sNjG1JkzGmvWrInnnnsuBgwYUDWtQYMGMWDAgHj22WdzODKALVNeXh4RX/4gCKmtW7cu7rnnnli1alUccMABuR4O9cyoUaNi0KBB1X4ug69rlOsB1NaSJUti3bp10bp162rTW7duHa+99lqORgWwZSorK2PMmDFx0EEHRc+ePXM9HOqRl156KQ444ID47LPPori4OKZOnRrdu3fP9bCoR+655554/vnnY+7cubkeCtu4OhMaAPXJqFGj4uWXX46nn34610Ohntl9991j3rx5UV5eHlOmTIlhw4bFU089JTZI4r333ovRo0fHjBkzorCwMNfDYRtXZ0KjRYsW0bBhw/joo4+qTf/oo4+iTZs2ORoVwOY766yz4q9//WvMnj072rdvn+vhUM/k5+dHt27dIiKid+/eMXfu3Ljhhhvi1ltvzfHIqA+ee+65WLx4ceyzzz5V09atWxezZ8+OG2+8MT7//PNo2LBhDkfItqTO3KORn58fvXv3jieeeKJqWmVlZTzxxBOuPQXqhCzL4qyzzoqpU6fGzJkzo3PnzrkeEtuBysrK+Pzzz3M9DOqJww47LF566aWYN29e1ce+++4bp5xySsybN09kUE2dOaMREXHuuefGsGHDYt9994399tsvxo8fH6tWrYrhw4fnemjUExUVFfHmm29WfV5WVhbz5s2LnXbaKTp06JDDkVEfjBo1KiZPnhwPPfRQlJSUxIcffhgREc2aNYvGjRvneHTUBxdccEEMHDgwOnToECtXrozJkyfHrFmzYvr06bkeGvVESUlJjfvKioqKYuedd3a/GTXUqdA48cQT4+OPP45LLrkkPvzwwygtLY3HHnusxg3isKX+/ve/x6GHHlr1+bnnnhsREcOGDYtJkyblaFTUFxMnToyIiEMOOaTa9Ntvvz1OP/30735A1DuLFy+O0047LRYtWhTNmjWLXr16xfTp0+Pwww/P9dCA7VCd+jsaAABA3VBn7tEAAADqDqEBAAAkJzQAAIDkhAYAAJCc0AAAAJITGgAAQHJCAwAASE5oAAAAyQkNAAAgOaEBAAAkJzQAAIDk/g9FqDg+VAR1FgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Divide the figure in two subfigures\n",
    "f, ax = plt.subplots(2, 1, figsize=(10,4), sharex=True)     # The subfigures share the same x axis\n",
    "\n",
    "# In the first subfigure we show our observations\n",
    "ax[0].set_yticks([])                                        # Do not show the ticks in the vertical axis\n",
    "ax[0].set_ylabel(\"\")                                        # Do not show any label in the vertical axis\n",
    "ax[0].spines[['left', 'right', 'top']].set_visible(False)   # Do not plot the left, right and top margins\n",
    "ax[0].text(2.2, 0.75, \"Data\", fontsize=12)\n",
    "\n",
    "for i in range(len(Dat)):\n",
    "    ax[0].axvline(x=Dat[i], ymax=0.2, c='r')                # Plot the observations as a rugplot\n",
    "\n",
    "# In the second subfigure we show the central tendency statistics\n",
    "ax[1].spines[['left', 'right', 'top']].set_visible(False)   # Do not plot the left, right and top margins\n",
    "\n",
    "for i in range(len(CentralTendency)):\n",
    "    ax[1].axvline(x=CentralTendency[i], ymax=0.2)\n",
    "\n",
    "ax[1].set_yticks([])\n",
    "ax[1].set_ylabel(\"\")\n",
    "ax[1].text(1.75, 0.75, \"Central tendency statistics\", fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DispersionStatistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "S2Biased = np.var(Dat)\n",
    "S2Unbiased = np.var(Dat, ddof=1)\n",
    "\n",
    "S1 = np.sqrt(S2Biased)\n",
    "S2 = np.sqrt(S2Unbiased)\n",
    "R = Dat.max() - Dat.min()\n",
    "IQR = iqr(Dat)\n",
    "MAD = median_abs_deviation(Dat)\n",
    "AAD = abs(Dat-Dat.mean()).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def an(n):\n",
    "    return np.sqrt((n-1)/2) * gamma((n-1)/2) / gamma(n/2)\n",
    "\n",
    "def c4(n):\n",
    "    return 1/an(n)\n",
    "\n",
    "# d2 for n in 2,...,25\n",
    "d2 = {2:1.128, 3:1.693, 4:2.059, 5:2.326, 6:2.534, 7:2.704, 8:2.847, 9:2.970, 10:3.078, 11:3.173, 12:3.258, 13:3.336, 14:3.407, 15:3.472, 16:3.532, 17:3.588, 18:3.640, 19:3.689, 20:3.735, 21:3.778, 22:3.819, 23:3.858, 24:3.895, 25:3.931}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate unbiased estimators for the standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(Dat)\n",
    "\n",
    "sigma_1 = S2 / c4(N)\n",
    "sigma_2 = R / d2[N]\n",
    "sigma_3 = IQR / (2 * norm.ppf(0.75))\n",
    "sigma_4 = MAD / norm.ppf(0.75)\n",
    "sigma_5 = AAD * np.sqrt(np.pi/2)\n",
    "\n",
    "sigma_list = [sigma_1, sigma_2, sigma_3, sigma_4, sigma_5]\n",
    "labels = ['Deviation', 'Range', 'IQR', 'MAD', 'AAD']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we assume that the data follows a normal distribution, a prediction interval of prob. $(1-\\alpha)\\times100\\%$ would have limits corresponding with the quantiles of the normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.05, 0.5)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyEAAAH5CAYAAACWHDGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA850lEQVR4nO3deXyNZ/7/8ffJIZtEUmRBqag9ljaKalXVGlUjU2Mrra01tU2NtlQXSztUaYuO0poJulh/xtJ2aq+lQlGqqFgbqoTYt4gl5/794ZszTpMQJNcdyev5eJzHI7nOfV/35z7nFud97uu6b4dlWZYAAAAAwBAvuwsAAAAAkL8QQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQDYYurUqXI4HO6Hr6+vSpQooWbNmumjjz7SuXPnbqvftWvXaujQoTp9+nT2FgwAyDaEEACArd5++2198cUXmjhxovr27StJ6tevn6pVq6atW7fecn9r167VsGHDCCEAkIsVsLsAAED+1rx5cz300EPu3wcNGqTvvvtOTz31lP70pz8pPj5efn5+NlYIAMhunAkBAOQ6DRs21FtvvaUDBw7oyy+/lCRt3bpVXbp0UdmyZeXr66vw8HB169ZNJ06ccK83dOhQvfrqq5KkiIgI91Cv/fv3S5KmTJmihg0bKjQ0VD4+PqpSpYomTpxofP8AIL/jTAgAIFd69tln9frrr2vJkiV64YUXtHTpUv3666/q2rWrwsPD9csvv2jSpEn65Zdf9MMPP8jhcOjpp5/W7t27NWPGDI0ZM0bFihWTJIWEhEiSJk6cqMjISP3pT39SgQIF9PXXX6tXr15yuVzq3bu3nbsLAPmKw7Isy+4iAAD5z9SpU9W1a1dt3LjRYzjW9YKDg1W2bFlt3rxZFy9eTDcsa+bMmerQoYNWr16txx57TJL0/vvv69VXX1VCQoLKlCnjsXxGfURHR2vPnj3at29f9u0cAOCGGI4FAMi1AgIC3FfJuj48pKSk6Pjx43r44YclSZs3b85Sf9f3cebMGR0/flyPP/64fv31V505cyYbKwcA3AghBACQa50/f16BgYGSpJMnT+qll15SWFiY/Pz8FBISooiICEnKcoCIi4tT48aNVahQIQUHByskJESvv/76LfUBALhzzAkBAORKv//+u86cOaNy5cpJktq2bau1a9fq1Vdf1QMPPKCAgAC5XC5FR0fL5XLdtL99+/apUaNGqlSpkj788EOVKlVK3t7e+vbbbzVmzJgs9QEAyB6EEABArvTFF19Ikpo1a6ZTp05p+fLlGjZsmAYPHuxeZs+ePenWczgcGfb39ddf69KlS/rqq69UunRpd/uKFSuyuXIAwM0wHAsAkOt89913eueddxQREaGOHTvK6XRKkv54LZWxY8emW7dQoUKSlO5mhRn1cebMGU2ZMiUbKwcAZAVnQgAAtlq4cKF27typq1ev6ujRo/ruu++0dOlS3Xffffrqq6/k6+srX19f1a9fX6NGjdKVK1dUsmRJLVmyRAkJCen6q1mzpiTpjTfeUPv27VWwYEG1bNlSTZs2lbe3t1q2bKm//vWvOn/+vP71r38pNDRUiYmJpncbAPI1QggAwFZpw6u8vb1VpEgRVatWTWPHjlXXrl3dk9Ilafr06erbt68+/vhjWZalpk2bauHChSpRooRHf7Vq1dI777yjTz75RIsWLZLL5VJCQoIqVqyoOXPm6M0339Qrr7yi8PBw9ezZUyEhIerWrZvRfQaA/I77hAAAAAAwijkhAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCqgN0FII+yLCk5+drP/v6Sw2FvPUjP1HuU27fjcknHj1/7uVgxySsL383kteM7r+2PSbx2AOx0F/8N4kwIckZyshQQcO2R9o8DuYup9yi3b+f4cSks7NojLYzk1LZyq7y2Pybx2gGw0138N4gQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCqwJ12YFmWzp07lx21IC+5cOF/P589K6Wm2lcLMmbqPcrt27n+79e5c5Kvb85tK7fKa/tjEq8dADvl4r9BgYGBcjgcmT7vsCzLupMNnD17VkFBQXfSBQAAAIA85MyZMypcuHCmz99xCLn+TMjZs2dVqlQpHTx48IYbRf7A8YA0HAu4HscD0nAsIA3HQt5zszMhdzwcy+FwpDtYChcuzAEEN44HpOFYwPU4HpCGYwFpOBbyDyamAwAAADCKEAIAAADAqGwNIT4+PhoyZIh8fHyys1vcpTgekIZjAdfjeEAajgWk4VjIf+54YjoAAAAA3AqGYwEAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCqgN0FAAAAIH9ITU3VlStX7C4Dd6BgwYJyOp133A8hBAAAADnKsiwdOXJEp0+ftrsUZIPg4GCFh4fL4XDcdh+EEAAAAOSotAASGhoqf3//O/rwCvtYlqXk5GQlJSVJkooXL37bfRFCAAAAkGNSU1PdAaRo0aJ2l4M75OfnJ0lKSkpSaGjobQ/NYmI6AAAAckzaHBB/f3+bK0F2SXsv72R+DyEEAAAAOY4hWHlHdryXhBAAAAAARhFCAAAAAButXLlSDofjjq8ell39mEAIAQAAADLQpUsXORwOORwOFSxYUGFhYWrSpIkmT54sl8uVbdt55JFHlJiYqKCgoCyv06BBA/Xr1++O+7ELIQQAAADIRHR0tBITE7V//34tXLhQTzzxhF566SU99dRTunr1arZsw9vb+47vu5Gd/ZhACAGAfGro0KEe/1FNnTpVDodD+/fvz5b+s7u/u2XbpuWnfQXs4OPjo/DwcJUsWVJRUVF6/fXXtWDBAi1cuFBTp06VJJ0+fVrPP/+8QkJCVLhwYTVs2FA///yzJGn37t1yOBzauXOnR79jxozR/fffLyn9MKoTJ06oQ4cOKlmypPz9/VWtWjXNmDHDvW6XLl20atUqjRs3zn2mZv/+/RkOx/rPf/6jyMhI+fj4qEyZMvrggw886ihTpoxGjBihbt26KTAwUKVLl9akSZOy+VVMjxACADZI++CY9vD19VWFChXUp08fHT161O7ybsnatWs1dOjQu2IMclbktf252126dEkDBw5UiRIl5Ofnpzp16mjp0qU3XW/jxo3q06ePIiMjVahQIZUuXVpt27bV7t27DVSNrEq+fNXYIzs1bNhQNWrU0Ny5cyVJbdq0UVJSkhYuXKhNmzYpKipKjRo10smTJ1WhQgU99NBDmjZtmkcf06ZN0zPPPJNh/ykpKapZs6b++9//avv27erRo4eeffZZbdiwQZI0btw41a1bVy+88IISExOVmJioUqVKpetn06ZNatu2rdq3b69t27Zp6NCheuutt9zhKc0HH3yghx56SD/99JN69eqlnj17ateuXdnwSmWOmxUCgI3efvttRUREKCUlRWvWrNHEiRP17bffavv27cavqf/ss8+qffv28vHxuaX11q5dq2HDhqlLly4KDg6+4/6yw51sO7P9gT26dOmiOXPmqF+/fipfvrymTp2qJ598UitWrFC9evUyXe+9995TXFyc2rRpo+rVq+vIkSMaP368oqKi9MMPP6hq1aoG9wKZqTJ4sbFt7R/ZIlv7q1SpkrZu3ao1a9Zow4YNSkpKcv/Nef/99zV//nzNmTNHPXr0UMeOHTV+/Hi98847kq6dHdm0aZO+/PLLDPsuWbKkXnnlFffvffv21eLFizV79mzVrl1bQUFB8vb2lr+/v8LDwzOt8cMPP1SjRo301ltvSZIqVKigHTt2aPTo0erSpYt7uSeffFK9evWSJA0cOFBjxozRihUrVLFixTt6jW6EEAIANmrevLkeeughSdLzzz+vokWL6sMPP9SCBQvUoUOHDNe5cOGCChUqlO21OJ3O277zrYn+7pZtZyan3re8bMOGDZo5c6ZGjx7t/kD23HPPqWrVqhowYIDWrl2b6br9+/fX9OnT5e3t7W5r166dqlWrppEjR2b64Q/IKsuy5HA49PPPP+v8+fPp7gZ/8eJF7du3T5LUvn17vfLKK/rhhx/08MMPa9q0aYqKilKlSpUy7Ds1NVUjRozQ7NmzdejQIV2+fFmXLl265S+n4uPj1apVK4+2Rx99VGPHjlVqaqr772T16tXdzzscDoWHhyspKemWtnWrGI4FALlIw4YNJUkJCQmS/jdvY8eOHXrmmWd0zz33uL/9PXTokLp166awsDD5+PgoMjJSkydPzrDfNWvWqFatWvL19dX999+vTz/9NN0ymc0tOHTokLp3764SJUrIx8dHERER6tmzpy5fvqyhQ4fq1VdflSRFRER4jE3OrL+ffvpJzZs3V+HChRUQEKBGjRrphx9+8Fgmbb/37t3rPiMRFBSkrl27Kjk5+aavY0bbzkqfN9qf61+Pm73umb1vc+bMkcPh0KpVq9LV/Omnn8rhcGj79u2SpAMHDqhXr16qWLGi/Pz8VLRoUbVp0+amcz/OnTunfv36qUyZMvLx8VFoaKiaNGmizZs33/R1u5FZs2YpKipKfn5+qly5spYtWybLshQZGanhw4ffUd+ZmTNnjpxOp3r06OFu8/X1Vffu3bVu3TodPHgw03UfeeQRjwAiSeXLl1dkZKTi4+NzpF7cuh1vNzP2yG7x8fGKiIjQ+fPnVbx4cW3ZssXjsWvXLvffk/DwcDVs2FDTp0+XJE2fPl0dO3bMtO/Ro0dr3LhxGjhwoFasWKEtW7aoWbNmunz5crbvhyQVLFjQ43eHw5GtV//KCGdCACAXSfvW7I/fqLVp00bly5fXiBEjZFmWjh49qocfflgOh0N9+vRRSEiIFi5cqO7du+vs2bMel23ctm2bmjZtqpCQEA0dOlRXr17VkCFDFBYWdtN6Dh8+rNq1a+v06dPq0aOHKlWqpEOHDmnOnDlKTk7W008/rd27d2vGjBkaM2aMihUrJkkKCQnJsL9ffvlFjz32mAoXLqwBAwaoYMGC+vTTT9WgQQOtWrVKderU8Vi+bdu2ioiI0LvvvqvNmzfr3//+t0JDQ/Xee+/dysua5T5vtj+38rpL6d+3Fi1aKCAgQLNnz9bjjz/useysWbMUGRnpHia0ceNGrV27Vu3bt9e9996r/fv3a+LEiWrQoIF27NiR6TeiL774oubMmaM+ffqoSpUqOnHihNasWaP4+HhFRUXd1mv25ptvavjw4erSpYu6d++ucePG6bnnntO//vUv/f777+rTp0+6da5cuaIzZ85kqf8iRYrIyyv996I//fSTKlSooMKFC3u0165dW5K0ZcuWDMfBZybt305kZGSW10HO8ve+Oz+Kfvfdd9q2bZv+/ve/695779WRI0dUoEABlSlTJtN1OnbsqAEDBqhDhw769ddf1b59+0yXjYuLU6tWrdSpUydJksvl0u7du1WlShX3Mt7e3kpNTb1hnZUrV1ZcXFy6vitUqGD/2WILAGDclClTLEnWsmXLrGPHjlkHDx60Zs6caRUtWtTy8/Ozfv/9d8uyLGvIkCGWJKtDhw4e63fv3t0qXry4dfz4cY/29u3bW0FBQVZycrK7LSYmxvL19bUOHDjgbtuxY4fldDqt6/8bSKspISHB3fbcc89ZXl5e1saNG9Ptg8vlsizLskaPHp1uvcz6i4mJsby9va19+/a52w4fPmwFBgZa9evXd7el7Xe3bt08+vzzn/9sFS1aNF0tf5TRtrPaZ2b7Y1lZf90ze98sy7I6dOhghYaGWlevXnW3JSYmWl5eXtbbb7/tbrv+PUyzbt06S5L1+eefZ7qvQUFBVu/evTN4VW7P6tWrLUnWwIED3W1z5syxJFlVq1b1aL/eihUrLElZemT0WluWZUVGRloNGzZM1/7LL79YkqxPPvnklvbliy++sCRZsbGxt7Qe7szFixetHTt2WBcvXrS7lFvWuXNnKzo62kpMTLR+//13a9OmTdbw4cOtgIAA66mnnrKuXr1quVwuq169elaNGjWsxYsXWwkJCVZcXJz1+uuve/ztPHv2rOXn52fVqFHDatSokcd20v69nDp1yrIsy/r73/9ulSpVyoqLi7N27NhhPf/881bhwoWtVq1audd54YUXrFq1alkJCQnWsWPHrNTU1HT9bNq0yf23ZdeuXdbUqVMtPz8/a8qUKe5+7rvvPmvMmDEe9dSoUcMaMmRIpq9LdrynDMcCABs1btxYISEhKlWqlNq3b6+AgADNmzdPJUuW9FjuxRdfdP9sWZb+85//qGXLlrIsS8ePH3c/mjVrpjNnzriH3qSmpmrx4sWKiYlR6dKl3X1UrlxZzZrdeHiCy+XS/Pnz1bJlS/e8levd6nXoU1NTtWTJEsXExKhs2bLu9uLFi+uZZ57RmjVrdPbs2Uz3W5Iee+wxnThxIt1yt+J2+7yV1z2zbUnX5iUkJSVp5cqV7rY5c+bI5XKpXbt27jY/Pz/3z1euXNGJEydUrlw5BQcH33BoVXBwsNavX6/Dhw/fcH+yaty4cbrnnnv05ptvutvSvu3du3ev/v73v2e4Xo0aNbR06dIsPTKbWHvx4sUMLy7g6+vrfj6rdu7cqd69e6tu3brq3LlzltcDFi1apOLFi6tMmTKKjo7WihUr9NFHH2nBggVyOp1yOBz69ttvVb9+fXXt2lUVKlRQ+/btdeDAAY8zzoGBgWrZsqV+/vnnGw7Fkq6dfYyKilKzZs3UoEEDhYeHKyYmxmOZV155RU6nU1WqVFFISIh+++23dP1ERUVp9uzZmjlzpqpWrarBgwfr7bff9piUbpe78xwYAOQRH3/8sSpUqKACBQooLCxMFStWzHBYSkREhPvnY8eO6fTp05o0aVKm13JPm1B47NgxXbx4UeXLl0+3TMWKFfXtt99mWtuxY8d09uzZbLuK0LFjx5ScnJzh1VYqV64sl8ulgwcPegyVuT44SdI999wjSTp16lS6ITpZdbt93srrnub69y1NdHS0goKCNGvWLDVq1EjStaFYDzzwgCpUqOBe7uLFi3r33Xc1ZcoUHTp0SJZluZ+70TCnUaNGqXPnzipVqpRq1qypJ598Us8995xH8MuqtODYsmVLBQQEpHu+a9eumQ7ru+eee9S4ceNb3ub1/Pz8dOnSpXTtKSkp7uez4siRI2rRooWCgoLc80yArJg6dWq6y9lmJDAwUB999JE++uijGy43a9YszZo1K117gwYNPP6NFylSRPPnz79hXxUqVNC6des82sqUKePRjyS1bt1arVu3zrSfjOaZbdmy5Ybbzg6EEACwUe3atTM8y/BH13/YSpss2KlTp0y/0b3+Sid3s8w+LP7xP1kTfd7O657Rh2QfHx/FxMRo3rx5mjBhgo4ePaq4uDiNGDHCY7m+fftqypQp6tevn+rWraugoCA5HA61b9/+hhNG27Ztq8cee0zz5s3TkiVLNHr0aL333nuaO3eumjdvfsN9/KNff/1V586dSzeX5NixY5Kk3r17Z7ru5cuXdfLkySxtJyQkJMP3pXjx4jp06FC69sTERElSiRIlbtr3mTNn1Lx5c50+fVrff/99ltYBkPMIIQBwlwkJCVFgYKBSU1Nv+k1zSEiI/Pz8tGfPnnTP3exGVGl3/k27WlNmsjosKyQkRP7+/hlud+fOnfLy8rqlScY5JbP9uZXX/WbatWunzz77TMuXL1d8fLwsy/IYiiVdG6LVuXNnj7sbp6SkZOkmisWLF1evXr3Uq1cvJSUlKSoqSsOHD7/lEJIWNtIm6Kd59913M2y/3tq1a/XEE09kaTsJCQkZTuh94IEHtGLFCp09e9bjLNX69evdz99ISkqKWrZsqd27d2vZsmUek3oB2IsQAgB3GafTqdatW2v69Onavn17uuFSx44dc1/Nyel0qlmzZpo/f75+++0391Ck+Ph4LV5845uEeXl5KSYmRl9++aV+/PHHdGdsrP+7Rn7avS9u9uHY6XSqadOmWrBggfbv3+/+0Hn06FFNnz5d9erVu+0hVtkps/25ldf9Zho3bqwiRYpo1qxZio+PV+3atdMN3XI6nenOzvzzn/+84dVwUlNTdf78eQUFBbnbQkNDVaJEiQyHNd1MWj/XB9Hp06dr9erVkv43LCojaXNCsiKzOSF/+ctf9P7772vSpEnu+4RcunRJU6ZMUZ06dTxCa3Jysn777TcVK1ZMxYoVU2pqqtq1a6d169ZpwYIFqlu3bpZqAWAGIQQA7kIjR47UihUrVKdOHb3wwguqUqWKTp48qc2bN2vZsmUew2CGDRumRYsW6bHHHlOvXr109epV/fOf/1RkZKS2bt16w+2MGDFCS5Ys0eOPP64ePXqocuXKSkxM1P/7f/9Pa9asUXBwsGrWrClJeuONN9S+fXsVLFhQLVu2zLC/f/zjH1q6dKnq1aunXr16qUCBAvr000916dIljRo1KvteoDuQ2f4UKlToll73GylYsKCefvppzZw5UxcuXND777+fbpmnnnpKX3zxhYKCglSlShWtW7dOy5YtS3f55uudO3dO9957r/7yl7+oRo0aCggI0LJly7Rx40aPMyrStTM+jz/+uMcE+T+qXLmyIiIi9NFHH8nf319eXl4aOXKk2rZtq9mzZ2vo0KHq37+/qlWrlm7d7JgTUqdOHbVp00aDBg1SUlKSypUrp88++0z79+9XbGysx7IbNmzQE088oSFDhmjo0KF6+eWX9dVXX6lly5Y6efJkupsTpl36FIA9CCEAcBcKCwvThg0b9Pbbb2vu3LmaMGGCihYtqsjIyHT30KhevboWL16s/v37a/Dgwbr33ns1bNgwJSYm3jSElCxZUuvXr9dbb72ladOm6ezZsypZsqSaN2/uvk9FrVq19M477+iTTz7RokWL5HK53Ddb/KPIyEh9//33GjRokN599125XC7VqVNHX375Zbp7hNgls/0pVKjQLb3uN9OuXTv9+9//lsPhUNu2bdM9P27cODmdTk2bNk0pKSl69NFHtWzZshte1czf31+9evXSkiVLNHfuXLlcLpUrV04TJkxQz5493cudP39e0rVhWzfi5eWluXPnqmfPnnrvvffk7++vl19+We+8844CAgL0+eefq0GDBhmGkOzy+eef66233tIXX3yhU6dOqXr16vrmm29Uv379G66XNrH266+/1tdff53ueUIIYC+HdSez+wAAwF3n22+/1VNPPaWff/45RwMEIF0btpeQkKCIiAj35ZVxd8uO95T7hAAAkM+sWLFC7du3J4AAsA3DsQAAyGdGjx5tdwkA8jnOhAAAAAAwihACAAAAwChCCAAAAJCBLl26yOFwyOFwqGDBgoqIiNCAAQNueI8cZA1zQgAAAIBMREdHa8qUKbpy5Yo2bdqkzp07y+Fw3PJlueGJMyEAAABAJnx8fBQeHq5SpUopJiZGjRs31tKlSyVJJ06cUIcOHVSyZEn5+/urWrVqmjFjhsf6DRo00N/+9jcNGDBARYoUUXh4uIYOHeqxzM6dO1WvXj35+vqqSpUqWrZsmRwOh+bPn+9e5uDBg2rbtq2Cg4NVpEgRtWrVSvv378/hvc85hBAAAADY4/IFc49ssH37dq1du1be3t6Srt0vo2bNmvrvf/+r7du3q0ePHnr22We1YcMGj/U+++wzFSpUSOvXr9eoUaP09ttvu4NMamqqYmJi5O/vr/Xr12vSpEl64403PNa/cuWKmjVrpsDAQH3//feKi4tTQECAoqOjdfny5WzZN9MYjgUAAAB7jChhbltDz9zWat98840CAgJ09epVXbp0SV5eXho/frwkqWTJknrllVfcy/bt21eLFy/W7NmzVbt2bXd79erVNWTIEElS+fLlNX78eC1fvlxNmjTR0qVLtW/fPq1cuVLh4eGSpOHDh6tJkybu9WfNmiWXy6V///vfcjgckqQpU6YoODhYK1euVNOmTW9r3+xECAEAAAAy8cQTT2jixIm6cOGCxowZowIFCqh169aSrp3FGDFihGbPnq1Dhw7p8uXLunTpkvz9/T36qF69usfvxYsXV1JSkiRp165dKlWqlDuASPIIMJL0888/a+/evQoMDPRoT0lJ0b59+7JtX00ihAAAAMAerx+2u4KbKlSokMqVKydJmjx5smrUqKHY2Fh1795do0eP1rhx4zR27FhVq1ZNhQoVUr9+/dINkSpYsKDH7w6HQy6XK8s1nD9/XjVr1tS0adPSPRcSEnIbe2U/QggAAADs4V3I7gpuiZeXl15//XX1799fzzzzjOLi4tSqVSt16tRJkuRyubR7925VqVIly31WrFhRBw8e1NGjRxUWFiZJ2rhxo8cyUVFRmjVrlkJDQ1W4cOHs2yEbMTEdAAAAyKI2bdrI6XTq448/Vvny5bV06VKtXbtW8fHx+utf/6qjR4/eUn9NmjTR/fffr86dO2vr1q2Ki4vTm2++KUnu+R8dO3ZUsWLF1KpVK33//fdKSEjQypUr9be//U2///57tu+jCYQQAAAAIIsKFCigPn36aNSoUXr55ZcVFRWlZs2aqUGDBgoPD1dMTMwt9ed0OjV//nydP39etWrV0vPPP+++Opavr68kyd/fX6tXr1bp0qX19NNPq3LlyurevbtSUlLu2jMjDsuyLLuLAAAAQN6UkpKihIQERUREuD9U48bi4uJUr1497d27V/fff7/d5aSTHe8pc0IAAAAAG82bN08BAQEqX7689u7dq5deekmPPvporgwg2YUQAgAAANjo3LlzGjhwoH777TcVK1ZMjRs31gcffGB3WTmK4VgAAADIMQzHynuy4z1lYjoAAAAAowghAAAAAIwihAAAAAAwihACAAAAwCiujgUgX0t1WdqQcFJJ51IUGuir2hFF5PRy2F0WAAB5GiEEQL61aHuihn29Q4lnUtxtxYN8NaRlFUVXLW5jZQAA5G0MxwKQLy3anqieX272CCCSdORMinp+uVmLtifaVBkAAHkfZ0JskHz5qt0lAPlaqsvSkK9+UUY3SUprG/LVL3q0XDGGZgE28ffmIwqQl/Ev3AZVBi+2uwQAN3H07CVVG7rE7jKAfGv/yBZ2lwCoS5cuOn36tObPny9JOnjwoIYMGaJFixbp+PHjKl68uGJiYjR48GAVLVrUvV6DBg20atUqSZKPj49Kly6trl276rXXXpPDwZdbEsOxAAAAgJv69ddf9dBDD2nPnj2aMWOG9u7dq08++UTLly9X3bp1dfLkSY/lX3jhBSUmJmrXrl0aNGiQBg8erE8++cSm6nMfzoTYYMfbzewuAcjX1v96Ul2nbrzpclO61FKdskUMVAQA+VPylWRj2/Iv6H9H6/fu3Vve3t5asmSJ/Pz8JEmlS5fWgw8+qPvvv19vvPGGJk6c+L/t+fsrPDxcktS1a1eNHz9eS5cuVc+ePe+ojryCEGIDxrkC9qpfIUTFg3x15ExKhvNCHJLCg3xVv0IIc0IAIAfVmV7H2La2dd522+uePHlSixcv1vDhw90BJE14eLg6duyoWbNmacKECemGW1mWpTVr1mjnzp0qX778bdeQ1zAcC0C+4/RyaEjLKpKuBY7rpf0+pGUVAggAQJK0Z88eWZalypUrZ/h85cqVderUKR07dszdNmHCBAUEBMjHx0f169eXy+XS3/72N1Ml53p8JQ8gX4quWlwTO0Wlu09IOPcJAQBj1j+z3u4SbollZXT+/H+8vb3dP3fs2FFvvPGGTp06pSFDhuiRRx7RI488ktMl3jUIIQDyreiqxdWkSjh3TAcAm9zpPA1TypUrJ4fDofj4eP35z39O93x8fLxCQkIUHBzsbgsKClK5cuUkSbNnz1a5cuX08MMPq3HjxqbKztUYjgUgX3N6OVT3/qJq9UBJ1b2/KAEEAJBO0aJF1aRJE02YMEEXL170eO7IkSOaNm2aunTpkun6AQEBeumll/TKK6/c9GxKfkEIAQAAAG5i/PjxunTpkpo1a6bVq1fr4MGDWrRokZo0aaIKFSpo8ODBN1z/r3/9q3bv3q3//Oc/hirO3QghAAAAwE2UL19eGzduVNmyZdW2bVvdd999at68uSpUqKC4uDgFBATccP0iRYroueee09ChQ+VyuQxVnXs5LM4JAQAAIIekpKQoISFBERER8vX1tbucbDVkyBB9+OGHWrp0qR5++GG7yzEmO95TJqYDAAAAt2HYsGEqU6aMfvjhB9WuXVteXgwyyipCCAAAAHCbunbtancJdyXiGgAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAGSgS5cucjgcevHFF9M917t3bzkcDnXp0sWjfd26dXI6nWrRokW6dfbv3y+Hw+F+BAYGKjIyUr1799aePXtyajdyJUIIAAAAkIlSpUpp5syZunjxorstJSVF06dPV+nSpdMtHxsbq759+2r16tU6fPhwhn0uW7ZMiYmJ+vnnnzVixAjFx8erRo0aWr58eY7tR25DCAEAAAAyERUVpVKlSmnu3Lnutrlz56p06dJ68MEHPZY9f/68Zs2apZ49e6pFixaaOnVqhn0WLVpU4eHhKlu2rFq1aqVly5apTp066t69u1JTU3Nyd3INQggAAABs4UpONva4E926ddOUKVPcv0+ePFldu3ZNt9zs2bNVqVIlVaxYUZ06ddLkyZNlWdZN+/fy8tJLL72kAwcOaNOmTXdU692igN0FAAAAIH/aFVXT2LYq74y/7XU7deqkQYMG6cCBA5KkuLg4zZw5UytXrvRYLjY2Vp06dZIkRUdH68yZM1q1apUaNGhw021UqlRJ0rV5I7Vr177tWu8WhBAAAADgBkJCQtzDqyzLUosWLVSsWDGPZXbt2qUNGzZo3rx5kqQCBQqoXbt2io2NzVIISTtj4nA4sr3+3IgQAgAAAFtU3Hz3DD3q1q2b+vTpI0n6+OOP0z0fGxurq1evqkSJEu42y7Lk4+Oj8ePHKygo6Ib9x8dfO1MTERGRjVXnXoQQAAAA2MLL39/uErIsOjpaly9flsPhULNmzTyeu3r1qj7//HN98MEHatq0qcdzMTExmjFjRoaX+U3jcrn00UcfKSIiIt1k97yKEAIAAADchNPpdJ+tcDqdHs998803OnXqlLp3757ujEfr1q0VGxvrEUJOnDihI0eOKDk5Wdu3b9fYsWO1YcMG/fe//03Xd15FCAEAAACyoHDhwhm2x8bGqnHjxhkOuWrdurVGjRqlrVu3utdv3LixJMnf31/33XefnnjiCU2aNEnlypXLueJzGYeVleuGAQAAALchJSVFCQkJioiIkK+vr93lIBtkx3vKfUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAuIF169bJ6XSqRYsWmS4zY8YMOZ1O9e7dO91zK1eulMPhkMPhkJeXl4KCgvTggw9qwIABSkxMzMnScy1CCAAAAHADsbGx6tu3r1avXq3Dhw9nusyAAQM0Y8YMpaSkZLjMrl27dPjwYW3cuFEDBw7UsmXLVLVqVW3bti0ny8+VCCEAAABAJs6fP69Zs2apZ8+eatGihaZOnZpumYSEBK1du1avvfaaKlSooLlz52bYV2hoqMLDw1WhQgW1b99ecXFxCgkJUc+ePXN4L3IfQggAAABsceVSqrHH7Zo9e7YqVaqkihUrqlOnTpo8ebIsy/JYZsqUKWrRooWCgoLUqVMnxcbGZqlvPz8/vfjii4qLi1NSUtJt13g3KmB3AQAAAMifJr20yti2en/S8LbWi42NVadOnSRJ0dHROnPmjFatWqUGDRpIklwul6ZOnap//vOfkqT27dvr5ZdfVkJCgiIiIm7af6VKlSRJ+/fvV2ho6G3VeDfiTAgAAACQgV27dmnDhg3q0KGDJKlAgQJq166dx5mOpUuX6sKFC3ryySclScWKFVOTJk00efLkLG0j7ayKw+HI5upzN86EAAAAwBY9xj1udwk3FBsbq6tXr6pEiRLuNsuy5OPjo/HjxysoKEixsbE6efKk/Pz83Mu4XC5t3bpVw4YNk5fXjb/zj4+PlySVKVMmR/YhtyKEAAAAwBYFfZx2l5Cpq1ev6vPPP9cHH3ygpk2bejwXExOjGTNmqE2bNlqwYIFmzpypyMhI9/OpqamqV6+elixZoujo6Ey3cfHiRU2aNEn169dXSEhIju1LbkQIAQAAAP7gm2++0alTp9S9e3cFBQV5PNe6dWvFxsYqJSVFRYsWVdu2bdMNp3ryyScVGxvrEUKSkpKUkpKic+fOadOmTRo1apSOHz+e6dW08jJCCAAAAPAHsbGxaty4cboAIl0LIaNGjdKmTZvUs2fPDOdztG7dWs8++6yOHz/ubqtYsaIcDocCAgJUtmxZNW3aVP3791d4eHiO7ktu5LD+eI0xAAAAIJukpKS4rxTl6+trdznIBtnxnnJ1LAAAAABGMRwLQP7mSpUOrJXOH5UCwqT7HpG8cu9ESQAA8gJCCID8a8dX0qKB0tnD/2srXEKKfk+q8if76gIAII9jOBaA/GnHV9Ls5zwDiCSdTbzWvuMre+oCACAfyHtnQi5fsLsCALmdK1Va+KqkjK7L8X9tCwdIZRswNAtA5rwL2V3BXYVrIeUd2fFe5r0QMqLEzZcBgJs5lyiNLGV3FQBys6Fn7K7grlCwYEFJUnJyssddxXH3Sk5OlvS/9/Z25L0QAgAAgFzD6XQqODhYSUlJkiR/f/8M76uB3M+yLCUnJyspKUnBwcFyOm9/tEDeu08Iw7EA3Mz+tdL0v9x8uWfmSGUeyfl6ANydGI6VZZZl6ciRIzp9+rTdpSAbBAcHKzw8/I7CZN4LIQBwM65UaWzVa5PQM5wX4rh2lax+25gTAgDZKDU1VVeuXLG7DNyBggUL3tEZkDQMxwKQ/3g5r12Gd/ZzkhzyDCL/961O9EgCCABkM6fTmS0fYHH34xK9APKnKn+S2n4uFS7u2V64xLV27hMCAECOYTgWgPyNO6YDAGAcIQQAAACAUQzHAgAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgVAG7CwAAO6W6UrU5abOOJR9TiH+IokKj5PRy2l0WAAB5GiEEQL617MAyjdwwUkeTj7rbwvzD9Frt19T4vsY2VgYAQN7msCzLsrsIADBt2YFl6r+yvyx5/gl0yCFJ+rDBhwQRAAByCCEEuEskX0m2u4Q8w+VyqdWCVkq6mJTpMqF+oVrQaoG8vJg6d6f8C/rbXQIAIJchhAB3iWqfVbO7BOC2bOu8ze4SAAC5DF/xAQAAADCKMyHAXYLhWNln89HN6rm8502Xm9hooqLCogxUlLcxHAsA8EdcHQu4S/BBLvvULVFXYf5hSkpOSjcxXbo2OT3MP0x1S9Tlcr0AAOQAhmMByHecXk69Vvs1Sf+7GlaatN8H1h5IAAEAIIcQQgDkS43va6wPG3yoUP9Qj/Yw/zAuzwsAQA5jTgiAfI07pgMAYB4hBAAAAIBRDMcCAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGBUAbsLAAA7WampSv5xk64eO6YCISHyf6imHE6n3WUBAJCnEUIA5FtnlyzR0RHv6uqRI+62AuHhCnt9kAo3bWpjZQAA5G0Oy7Isu4sAANPOLlmiQy/1k/74J9DhkCSVHDeWIAIAQA4hhAA5xJWcbHcJyISVmqp9LZ5SalJSpss4Q0N1/3+/YWhWLuTl7293CQCAO0QIAXJIfKXKdpcA5EmVd8bbXQIA4A5xdSwAAAAARnEmBMghDMfKvS78uEm/9+hx0+XunTRJhR6qaaAi3AqGYwHA3Y+rYwE5hA9KuVfAo4+oQHi4rh49mn5iuiQ5HCoQFqaARx9hTggAADmA4VgA8h2H06mw1wf93y+OPzx57few1wcRQAAAyCGEEAD5UuGmTVVy3FgVCAvzaC8QFsbleQEAyGHMCQGQr3HHdAAAzCOEAAAAADCK4VgAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKMIIQAAAACMIoQAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwqoDdBQCAnVwuS4l7TuvC2UsqVNhHxcsHy8vLYXdZAADkaYQQAPnWvp+S9P2sPbpw+pK7rVCwjx5rV173PxhqY2UAAORtDsuyLLuLAADT9v2UpEWfbs/0+ei/ViWIAACQQzgTgiy5cinV7hKAbONyWVo9c/cNl/l+1m7dW6kIQ7OQJxT0cdpdAgB44EwIsuTjF7+zuwQAwG3q/UlDu0sAAA9cHQsAAACAUZwJQZYwHAt5yeE9p/TN+K03Xe6pPtVVovw9BioCchbDsQDkNswJQZbwHxjyklJViqpQsI/HVbH+KOAeH5WqUpQ5IQAA5ACGYwHId7y8HHqsXfkbLlOvbXkCCAAAOYThWADyrYzuExJwj4/qteU+IQAA5CRCCIB8jTumAwBgHiEEAAAAgFHMCQEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYRQgBAAAAYBQhBAAAAIBRhBAAAAAARhFCAAAAABhFCAEAAABgFCEEAAAAgFGEEAAAAABGEUIAAAAAGEUIAQAAAGAUIQQAAACAUYQQAAAAAEYRQgAAAAAYVeBOO7AsS+fOncuOWgAAAADkAYGBgXI4HJk+f8ch5Ny5cwoKCrrTbgAAAADkEWfOnFHhwoUzfd5hWZZ1Jxu4/kzI2bNnVapUKR08ePCGG0X+wPGANBwLuB7HA9JwLCANx0Lek+NnQhwOR7qDpXDhwhxAcON4QBqOBVyP4wFpOBaQhmMh/2BiOgAAAACjCCEAAAAAjMrWEOLj46MhQ4bIx8cnO7vFXYrjAWk4FnA9jgek4VhAGo6F/OeOJ6YDAAAAwK1gOBYAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjciyEDB8+XI888oj8/f0VHBycU5tBLvXxxx+rTJky8vX1VZ06dbRhwwa7S4INVq9erZYtW6pEiRJyOByaP3++3SXBJu+++65q1aqlwMBAhYaGKiYmRrt27bK7LNhk4sSJql69uvvu2HXr1tXChQvtLgu5wMiRI+VwONSvXz+7S0EOy7EQcvnyZbVp00Y9e/bMqU0gl5o1a5b69++vIUOGaPPmzapRo4aaNWumpKQku0uDYRcuXFCNGjX08ccf210KbLZq1Sr17t1bP/zwg5YuXaorV66oadOmunDhgt2lwQb33nuvRo4cqU2bNunHH39Uw4YN1apVK/3yyy92lwYbbdy4UZ9++qmqV69udykwIMfvEzJ16lT169dPp0+fzsnNIBepU6eOatWqpfHjx0uSXC6XSpUqpb59++q1116zuTrYxeFwaN68eYqJibG7FOQCx44dU2hoqFatWqX69evbXQ5ygSJFimj06NHq3r273aXABufPn1dUVJQmTJigf/zjH3rggQc0duxYu8tCDmJOCLLV5cuXtWnTJjVu3Njd5uXlpcaNG2vdunU2VgYgNzlz5oykax88kb+lpqZq5syZunDhgurWrWt3ObBJ79691aJFC4/PD8jbCthdAPKW48ePKzU1VWFhYR7tYWFh2rlzp01VAchNXC6X+vXrp0cffVRVq1a1uxzYZNu2bapbt65SUlIUEBCgefPmqUqVKnaXBRvMnDlTmzdv1saNG+0uBQbd0pmQ1157TQ6H44YPPmgCAG6kd+/e2r59u2bOnGl3KbBRxYoVtWXLFq1fv149e/ZU586dtWPHDrvLgmEHDx7USy+9pGnTpsnX19fucmDQLZ0Jefnll9WlS5cbLlO2bNk7qQd3uWLFisnpdOro0aMe7UePHlV4eLhNVQHILfr06aNvvvlGq1ev1r333mt3ObCRt7e3ypUrJ0mqWbOmNm7cqHHjxunTTz+1uTKYtGnTJiUlJSkqKsrdlpqaqtWrV2v8+PG6dOmSnE6njRUip9xSCAkJCVFISEhO1YI8wNvbWzVr1tTy5cvdE5BdLpeWL1+uPn362FscANtYlqW+fftq3rx5WrlypSIiIuwuCbmMy+XSpUuX7C4DhjVq1Ejbtm3zaOvatasqVaqkgQMHEkDysBybE/Lbb7/p5MmT+u2335SamqotW7ZIksqVK6eAgICc2ixygf79+6tz58566KGHVLt2bY0dO1YXLlxQ165d7S4Nhp0/f1579+51/56QkKAtW7aoSJEiKl26tI2VwbTevXtr+vTpWrBggQIDA3XkyBFJUlBQkPz8/GyuDqYNGjRIzZs3V+nSpXXu3DlNnz5dK1eu1OLFi+0uDYYFBgammxtWqFAhFS1alDljeVyOhZDBgwfrs88+c//+4IMPSpJWrFihBg0a5NRmkQu0a9dOx44d0+DBg3XkyBE98MADWrRoUbrJ6sj7fvzxRz3xxBPu3/v37y9J6ty5s6ZOnWpTVbDDxIkTJSnd3/8pU6bcdJgv8p6kpCQ999xzSkxMVFBQkKpXr67FixerSZMmdpcGwJAcv08IAAAAAFyP+4QAAAAAMIoQAgAAAMAoQggAAAAAowghAAAAAIwihAAAAAAwihACAAAAwChCCAAAAACjCCEAAAAAjCKEAAAAADCKEAIAAADAKEIIAAAAAKP+P6+SuPz319asAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(2, 1, figsize=(10,6), sharex=True, height_ratios=[1, 3])\n",
    "\n",
    "for i in range(len(Dat)):\n",
    "    ax[0].axvline(x=Dat[i], ymax=0.2, c='r')\n",
    "\n",
    "ax[0].set_yticks([])\n",
    "ax[0].set_ylabel(\"\")\n",
    "ax[0].spines[['left', 'right', 'top']].set_visible(False)\n",
    "ax[0].text(2, 0.75, \"Data\", fontsize=12)\n",
    "\n",
    "for i, sigma in enumerate(sigma_list):\n",
    "    ax[1].hlines(\n",
    "        0.1*(4-i), \n",
    "        norm.ppf(ALPHA/2, Median, np.sqrt(1+np.pi/(2*N))*sigma), \n",
    "        norm.ppf(1-ALPHA/2, Median, np.sqrt(1+np.pi/(2*N))*sigma),\n",
    "        color=colors[i],\n",
    "        label=labels[i])\n",
    "    ax[1].scatter(Median, 0.1*(4-i))\n",
    "\n",
    "ax[1].legend()\n",
    "ax[1].set_yticks([])\n",
    "ax[1].set_ylabel(\"\")\n",
    "ax[1].spines[['left', 'right', 'top']].set_visible(False)\n",
    "ax[1].text(1.5, 0.45, r\"Prediction intervals, $\\alpha={}$\".format(ALPHA), fontsize=12)\n",
    "ax[1].set_ylim(-0.05, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Skewness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Empirical central moments and k-statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2 = S2Biased\n",
    "m3 = ((Dat-Mean)**3).mean()\n",
    "\n",
    "k2 = S2Unbiased\n",
    "k3 = N**2 / ((N-1)*(N-2))*m3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coefficient of skewness based on central moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1_byhand = m3 / m2**(3/2)\n",
    "g1 = skew(Dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(1.472), np.float64(1.472))"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g1_byhand.round(4), g1.round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coefficient of skewness based on $k$-statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "G1_byhand = k3 / k2**(3/2)\n",
    "G1 = skew(Dat, bias=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(1.7455), np.float64(1.7455))"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G1_byhand.round(4), G1.round(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: ylabel='Density'>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABW1klEQVR4nO3deXhTZdoG8DtJm3RP97SlhS7sa6EttexIZXMBUQecBeiM6Oc2OtVPxRlhHHUKigyfwoDLKOAG4zqOSxUrlcWytbKvLXShbdKNJm3aJm1yvj/SBCot0vVkuX/XlUs4OTl9ItDcPe/zvq9EEAQBRERERC5EKnYBRERERH2NAYiIiIhcDgMQERERuRwGICIiInI5DEBERETkchiAiIiIyOUwABEREZHLYQAiIiIil+MmdgH2yGw2o6ysDL6+vpBIJGKXQ0RERNdBEATU1dUhIiICUum17/EwALWjrKwMUVFRYpdBREREXVBSUoLIyMhrnsMA1A5fX18Alv+Bfn5+IldDRERE10On0yEqKsr2OX4tDEDtsA57+fn5MQARERE5mOtpX2ETNBEREbkcBiAiIiJyOQxARERE5HIYgIiIiMjlMAARERGRy2EAIiIiIpfDAEREREQuhwGIiIiIXA4DEBEREbkcBiAiIiJyOQxARERE5HIYgIiIiMjlMAARERGRy+Fu8ETUoTPqOuw+V4n9F2qgN7Sgf6AXBql8cee4SCi93MUuj4ioyxiAiOgqNXojXvjyFD7Ou9jm+I8F1QCAdd+dxbLJsbhncgy85Pw2QkSOh9+5iKiN3ecq8ccPfsKlhmZIJMCUQSFIiQtCiI8CxTUNyDyuxhlNHdbuOIsdJzV4a2kSQnwVYpdNRNQpEkEQBLGLsDc6nQ5KpRJarRZ+fn5il0PUZ344W4llWw/B2GLG0DBf/H3BKIzrH9DmHLNZwH+PluFv/z2Jar0RUYGe2JI2HrEhPiJVTURk0ZnPbzZBExEAYNcV4WfmcBU+f2jSVeEHAKRSCebF98PH909A/0AvlNQ04tdv7EdlnUGEqomIuoYBiIhQUtOAB9/Ps4Wf9b8eB7nbtb89RAd74+P7JyA2xBtqXRMefC8PzSZzH1VMRNQ9DEBELq7ZZMYft/2EuqYWjO3vf13hxyrEV4HXf5cIH4UbDhTW4IUvT/VytUREPYMBiMjFrd1xFj8V18LXww2vLBp73eHHamCoD9b+agwAYPOPhfjhbGVvlElE1KMYgIhc2PFSLTb9UAAAWH3HaEQFenXpOjNHhGHphGgAwIr/HEdTs6mnSiQi6hUMQEQuShAE/PXzExAE4NYxEZg7Krxb13ts5mCo/BQoqm7Ahp35PVQlEVHvYAAiclGfHynDoaJL8HSX4em5Q7t9PV8Pd/z11hEAgE0/FCC/or7b1yQi6i0MQEQuqMHYglVfnwYAPDAtDuFKzx657uyRYZg2JATNJgH/2HG2R65JRNQbGICIXNDWnCKUa5sQGeCJZVNie+y6EokET80ZCokE+PJYOU6W6Xrs2kREPYkBiMjFNBpNeHP3eQDAo6mD4eEu69HrDw3zw82t/UTrvuNdICKyT3YRgDZs2IDo6Gh4eHggOTkZBw4c6PDcTz75BImJifD394e3tzfi4+PxzjvvtDln6dKlkEgkbR6zZ8/u7bdB5BA+OFCMqnojIgM8MS8+ole+xqOpgyCVAN+e1ODYRW2vfA0iou4QPQBt374d6enpWLlyJfLy8jBmzBjMmjULFRUV7Z4fGBiIP//5z8jJycHRo0eRlpaGtLQ0fPPNN23Omz17NsrLy22PDz74oC/eDpFda2o24bVdlmnvD0wbCHdZ73wLGBjqi3nx/QAAr35/rle+BhFRd4gegNauXYtly5YhLS0Nw4cPx6ZNm+Dl5YW33nqr3fOnTZuG22+/HcOGDUNcXBweeeQRjB49Gnv27GlznkKhQFhYmO0REHD1nkZEruaj3IvQ6AwIV3rgjoR+vfq1Hpw+EACw45QGxdUNvfq1iIg6S9QAZDQakZubi9TUVNsxqVSK1NRU5OTk/OLrBUFAVlYWzpw5gylTprR5Ljs7G6GhoRgyZAjuv/9+VFdXd3gdg8EAnU7X5kHkbARBwNt7LwAAlk2OhcKtZ3t/fm5gqA+mDg6BIFhWiCYisieiBqCqqiqYTCaoVKo2x1UqFdRqdYev02q18PHxgVwux80334xXX30VN910k+352bNnY+vWrcjKysLq1avxww8/YM6cOTCZ2l+dNiMjA0ql0vaIiorqmTdIZEd+LKhGQaUe3nIZ7kqM7JOv+ftJMQCAfx8qQV1Tc598TSKi6+EmdgFd4evri8OHD6O+vh5ZWVlIT09HbGwspk2bBgBYtGiR7dxRo0Zh9OjRiIuLQ3Z2NmbMmHHV9ZYvX4709HTb73U6HUMQOZ2tOYUAgDsSIuHr4d4nX3PKoGDEhXijoFKPDw9dtAUiIiKxiXoHKDg4GDKZDBqNps1xjUaDsLCwDl8nlUoxcOBAxMfH47HHHsOdd96JjIyMDs+PjY1FcHAw8vPbX55foVDAz8+vzYPImZTWNmLHScu/s9/dMKDPvq5EIkHaREvo2fxjIcxmoc++NhHRtYgagORyORISEpCVlWU7ZjabkZWVhZSUlOu+jtlshsFg6PD5ixcvorq6GuHh3dvriMhRvbevCGYBmBAXhEEq3z792gvG9YOvwg3FNQ3Yd77jXjwior4k+iyw9PR0vPHGG9iyZQtOnTqF+++/H3q9HmlpaQCAxYsXY/ny5bbzMzIysGPHDpw/fx6nTp3Cyy+/jHfeeQe//e1vAQD19fX43//9X+zbtw+FhYXIysrCvHnzMHDgQMyaNUuU90gkpmaTGf8+dBEAsDil7+7+WHnJ3XBr63pD/z5U0udfn4ioPaL3AC1cuBCVlZVYsWIF1Go14uPjkZmZaWuMLi4uhlR6Oafp9Xo88MADuHjxIjw9PTF06FC8++67WLhwIQBAJpPh6NGj2LJlC2praxEREYGZM2fiueeeg0KhEOU9Eolp19lKVNUbEOwjx4xhql9+QS/4VWIU3t9fjK+Pq/FsYzOUnn3Tg0RE1BGJIAgclP8ZnU4HpVIJrVbLfiByePe/m4uvj6vxh0kxeOaW4aLUIAgCZq/bjTOaOjw3f2Sf9iERkevozOe36ENgRNR7LumNyDplWVX9zoS+mfreHolEYpt6/++DHAYjIvExABE5sf8eLYPRZMbwcD8MCxf3buaCcZFwl0lwrFSL02ouNkpE4mIAInJiH+damp/FvPtjFegtx/QhoQCAzw+XiVwNEbk6BiAiJ1VQWY8jF7Vwk0p6bdf3zrqttY7/Hi0D2w+JSEwMQERO6suj5QCAyYOCEeRjHzMgZwxVwUsuQ0lNIw6X1IpdDhG5MAYgIif1xVHLMNPNo+3j7g8AeMpluGm4ZSr+50c4DEZE4mEAInJCZzV1OKuph7tMYgsc9uLW1kD25dFymLg1BhGJhAGIyAl90Tr8NWVQiN0tOjhlsKWmijoD9l/g1hhEJA4GICInIwgCvmwd/rpljP3tfyd3k2L2CMtmx18dKxe5GiJyVQxARE7mtLoOBZV6yN2kSBVp64tfMnuUJQB9e0LDHeKJSBQMQERO5uvjagCW4S9fD/sa/rKaEBcEH4UbKuoMOHyxVuxyiMgFMQAROZkdJzUAgNkjw0SupGMKNxmmD7UsivjNCbXI1RCRK2IAInIiJTUNOFWug1QC3NgaMOzVrBGW4blvjqu5KCIR9TkGICInYr37kxQdiEBvucjVXNu0IaGQu0lRWN2As5p6scshIhfDAETkRKwByN7W/mmPj8INkwYGA+AwGBH1PQYgIidR22DEgcIaAMDM4fbb/3Ml6zCYNbgREfUVBiAiJ/H96QqYzAKGhvmif5CX2OVclxuHWgLQsVItKnRNIldDRK6EAYjISXx3ynGGv6xCfBUYE6kEAOw8UyFyNUTkShiAiJxAs8mM3WerAAAz7HTxw45Yp8N/f5oBiIj6DgMQkRPILbqEOkMLgrzlGN1PKXY5nTKjdRhsz7kqGFpMIldDRK6CAYjICViHj6YODoFUKhG5ms4ZEeGHEF8F9EYTDlyoEbscInIRDEBETiD7dCUAYJqdL37YHqlUghuHcBiMiPoWAxCRgyutbcQZTR2kEmDKoGCxy+mSK/uAuCo0EfUFBiAiB7ez9a7JuP4B8Pey79WfOzJpUDDcZRIUVTegqLpB7HKIyAUwABE5uOzW/p/pDjj8ZeWjcEPCgAAAwK5zlSJXQ0SugAGIyIEZWkzYm18NwNIA7cimtNa/6ywDEBH1PgYgIgeWV1SLxmYTgn0UGB7uJ3Y53TJlkCUA5RRUw9hiFrkaInJ2DEBEDmx363DR5EHBDjf9/eeGh/shyFsOvdGEvOJLYpdDRE6OAYjIge0+Z1n9ebKDzv66klQqsb0PDoMRUW9jACJyUNX1Bhwv0wIAJg10/AAEXNEHxEZoIuplDEBEDmpvQTUEARga5otQPw+xy+kRk1v7gI6X6lBVbxC5GiJyZgxARA5q99nL/T/OIsT3cjP33vwqkashImfGAETkgARBwJ58a/+PY09//7kJcUEAgH3nq0WuhIicGQMQkQMqqKxHubYJcjcpxscEil1Oj5ow0BKAfixgACKi3sMAROSArOEgcUAAPNxlIlfTs5KiAyGTWrbFKK1tFLscInJSDEBEDiinNQClxAaJXEnP8/Vwx+hIJYDL75OIqKcxABE5GLNZsPXHWIeLnI21D+jHAjZCE1HvsIsAtGHDBkRHR8PDwwPJyck4cOBAh+d+8sknSExMhL+/P7y9vREfH4933nmnzTmCIGDFihUIDw+Hp6cnUlNTce7cud5+G0R94oymDpcamuEll2F0pL/Y5fSKlFjLzLacgmoIgiByNUTkjEQPQNu3b0d6ejpWrlyJvLw8jBkzBrNmzUJFRUW75wcGBuLPf/4zcnJycPToUaSlpSEtLQ3ffPON7ZwXX3wRr7zyCjZt2oT9+/fD29sbs2bNQlNTU1+9LaJeY+v/iQ6Eu0z0f8K9ImFAAOQyKcq1TSisbhC7HCJyQqJ/91y7di2WLVuGtLQ0DB8+HJs2bYKXlxfeeuutds+fNm0abr/9dgwbNgxxcXF45JFHMHr0aOzZsweA5e7PunXr8Je//AXz5s3D6NGjsXXrVpSVleGzzz7rw3dG1Ducuf/HylMuw9j+/gDYB0REvUPUAGQ0GpGbm4vU1FTbMalUitTUVOTk5Pzi6wVBQFZWFs6cOYMpU6YAAC5cuAC1Wt3mmkqlEsnJyR1e02AwQKfTtXkQ2SOTWcD+C639P3HOG4AAYEKcZRiMfUBE1BtEDUBVVVUwmUxQqVRtjqtUKqjV6g5fp9Vq4ePjA7lcjptvvhmvvvoqbrrpJgCwva4z18zIyIBSqbQ9oqKiuvO2iHrNyTId6ppa4Ktww4gIP7HL6VXWBm/2ARFRbxB9CKwrfH19cfjwYRw8eBAvvPAC0tPTkZ2d3eXrLV++HFqt1vYoKSnpuWKJelDOecvdkPExgXBz0v4fqzGR/vB0l6Fab8RZTb3Y5RCRk3ET84sHBwdDJpNBo9G0Oa7RaBAWFtbh66RSKQYOHAgAiI+Px6lTp5CRkYFp06bZXqfRaBAeHt7mmvHx8e1eT6FQQKFQdPPdEPU+awN0ipMPfwGA3E2KxOgA7D5XhR8LqjAkzFfskojIiYj6I6RcLkdCQgKysrJsx8xmM7KyspCSknLd1zGbzTAYLDtHx8TEICwsrM01dTod9u/f36lrEtmbZpMZBy/UAABucOIG6CtZ+4DYCE1EPU3UO0AAkJ6ejiVLliAxMRHjx4/HunXroNfrkZaWBgBYvHgx+vXrh4yMDACWfp3ExETExcXBYDDgq6++wjvvvIONGzcCACQSCR599FE8//zzGDRoEGJiYvDMM88gIiIC8+fPF+ttEnXbsVIt9EYTlJ7uth3Tnd2VG6OazAJkUonIFRGRsxA9AC1cuBCVlZVYsWIF1Go14uPjkZmZaWtiLi4uhlR6+UaVXq/HAw88gIsXL8LT0xNDhw7Fu+++i4ULF9rOeeKJJ6DX63HvvfeitrYWkyZNQmZmJjw8PPr8/RH1FOtdkBtiAyF1kSAwIsIPvgo36JpacLJMh1GtW2QQEXWXROD0iqvodDoolUpotVr4+bnGT9pk/3775n7sya/CX28djqUTY8Qup8/cs+UgvjtVgeVzhuK+qXFil0NEdqwzn9/OPY2EyEkYWkw4VGTp/0lp7YtxFSm29YDYB0REPYcBiMgBHCnRoqnZjCBvOQarfMQup09ZV7zOLboEk5k3rImoZzAAETkAW/9PXBAkEtfo/7EaEuYLXw831BtacKqcq7QTUc9gACJyAPvOO//+Xx2RSSVIHBAAANjfugwAEVF3MQAR2Tljixl5xZcAAMkxgSJXI47xMZbgd5ABiIh6CAMQkZ07VqqFocWMQG85Boa6Vv+P1fgYyx2gg4U13BeMiHoEAxCRnTtYaLnrkTggwOX6f6xG9fOHwk2Kar0RBZV6scshIifAAERk56zDPuNddPgLsOwLNra/P4DLgZCIqDsYgIjsmNks2D7wk6JdNwABwPjW93+AfUBE1AMYgIjs2BlNHXRNLfCSyzAiwrVXJU+KYQAiop7DAERkx6x3fxIGBMBN5tr/XMf1D4BMKkFpbSNKaxvFLoeIHJxrf0clsnPWux2uPvwFAN4KN4xsvQvG6fBE1F0MQER2ShDY//Nz1kbwA2yEJqJuYgAislPFNQ3Q6Axwl0lsM6BcXRIboYmohzAAEdkp64f86Eh/eLjLRK7GPlgDUH5FPWr0RpGrISJHxgBEZKc4/HW1AG85Bqssq2FzPSAi6g4GICI7dbDQsv+XdRsIsuAwGBH1BAYgIjtUUdeEC1V6SCRAwgDeAbqStRGad4CIqDsYgIjs0MELlrs/Q8P8oPR0F7ka+2K9A3SiTId6Q4vI1RCRo2IAIrJD1rsb46M5/PVzEf6eiAzwhMksIK/oktjlEJGDYgAiskO2BRBdeAPUa7HeBTrEAEREXcQARGRndE3NOKXWAbi8ASi1lTDAcmeMd4CIqKsYgIjsTG7RJQgCMCDIC6F+HmKXY5cSW4cGfyq+hBaTWeRqiMgRMQAR2Znc1unviZz91aFBob7wVbhBbzThtLpO7HKIyAExABHZmdzWYZ1ENkB3SCaVYKx1GKyYw2BE1HkMQER2pNlkxuGSWgCX+1yofQn9Lf9/DhUyABFR5zEAEdmRU+U6NDab4OfhhoEhPmKXY9esd8hy2QhNRF3AAERkR6wf5gkDAiCVSkSuxr7FR/lDKgFKaxuh1jaJXQ4RORgGICI7cuiKAETX5q1ww7BwPwDAoSJui0FEncMARGQnBEGwzQDj/l/XJ3EAh8GIqGsYgIjsRJm2CWpdE2RSCcZEKcUuxyGMYwAioi5iACKyE4da9/8aEeEHL7mbyNU4hsQrNkZtMHJjVCK6fgxARHYil/0/ndbP3xPhSg+YzAKOlGjFLoeIHAgDEJGdYADqmsvDYGyEJqLrxwBEZAfqDS04VW7ZAJVbYHQOG6GJqCsYgIjswJGSWpgFy5BOmJIboHZGwhUByGwWRK6GiBwFAxCRHThUyOGvrhoW7gdPdxl0TS0oqKwXuxwichB2EYA2bNiA6OhoeHh4IDk5GQcOHOjw3DfeeAOTJ09GQEAAAgICkJqaetX5S5cuhUQiafOYPXt2b78Noi6zLuTHDVA7z10mtS0bcIjDYER0nUQPQNu3b0d6ejpWrlyJvLw8jBkzBrNmzUJFRUW752dnZ+Puu+/Gzp07kZOTg6ioKMycOROlpaVtzps9ezbKy8ttjw8++KAv3g5Rp5nMAg4X1wLgHaCusvZNcWNUIrpeogegtWvXYtmyZUhLS8Pw4cOxadMmeHl54a233mr3/Pfeew8PPPAA4uPjMXToULz55pswm83Iyspqc55CoUBYWJjtERDADxayT2c1dagztMBbLsMQla/Y5TikhNY7Z3nFDEBEdH1EDUBGoxG5ublITU21HZNKpUhNTUVOTs51XaOhoQHNzc0IDGw7cyY7OxuhoaEYMmQI7r//flRXV3d4DYPBAJ1O1+ZB1FeswzZj+wfATSb6zyQOaVyUJQBdqNKjqt4gcjVE5AhE/W5bVVUFk8kElUrV5rhKpYJarb6uazz55JOIiIhoE6Jmz56NrVu3IisrC6tXr8YPP/yAOXPmwGQytXuNjIwMKJVK2yMqKqrrb4qok/K4/k+3Kb3cMVjlA4DT4Yno+jj0evurVq3Ctm3bkJ2dDQ+Py1OHFy1aZPv1qFGjMHr0aMTFxSE7OxszZsy46jrLly9Henq67fc6nY4hiPqMtQGaAah7EgYE4KymHnlFlzBrRJjY5RCRnRP1DlBwcDBkMhk0Gk2b4xqNBmFh1/4GtmbNGqxatQrffvstRo8efc1zY2NjERwcjPz8/HafVygU8PPza/Mg6gsVuiaU1DRCIgHG9vcXuxyHlmBthOYdICK6DqIGILlcjoSEhDYNzNaG5pSUlA5f9+KLL+K5555DZmYmEhMTf/HrXLx4EdXV1QgPD++Ruol6inW4ZojKF74e7iJX49isK0Ifu6iFoaX94W4iIivROy7T09PxxhtvYMuWLTh16hTuv/9+6PV6pKWlAQAWL16M5cuX285fvXo1nnnmGbz11luIjo6GWq2GWq1Gfb1lAbT6+nr87//+L/bt24fCwkJkZWVh3rx5GDhwIGbNmiXKeyTqiPVuBdf/6b4BQV4I8pbDaDLjeCk3RiWiaxM9AC1cuBBr1qzBihUrEB8fj8OHDyMzM9PWGF1cXIzy8nLb+Rs3boTRaMSdd96J8PBw22PNmjUAAJlMhqNHj+K2227D4MGD8Yc//AEJCQnYvXs3FAqFKO+RqCO2AMT9v7pNIpHY+qi4HhAR/RK7aIJ+6KGH8NBDD7X7XHZ2dpvfFxYWXvNanp6e+Oabb3qoMqLe09RswonWOxVsgO4ZCQMC8O1JDWeCEdEvEv0OEJGrOlJSixazgFBfBSIDPMUuxylYhxJziy5BELgxKhF1jAGISCRX9v9IJBKRq3EOI/spIZdJUa03orC6QexyiMiOMQARicS6AOK4/hz+6ikKNxlGRVo2Rs3jMBgRXQMDEJEIzGYBucXWO0BsgO5J1n6qXO4LRkTXwABEJILzVfWobWiGh7sUIyK48GZPst5Ry+VMMCK6BgYgIhFYZymNjvSHOzdA7VHWO0BnK+qgbWwWuRoislf8zkskAus6NYmc/t7jQnwVGBDkBUEADpfUil0OEdkpBiAiEeRyBeheldD/8nR4IqL2MAAR9bEavRHnq/QAOAOst4xrvbPGmWBE1BEGIKI+Zr0rMTDUB/5ecpGrcU7WPqCfii+hxWQWuRoiskcMQER97FBRDYDLwzTU8warfOGrcIPeaMIZTZ3Y5RCRHWIAIupj1mGZBPb/9BqZVIL4/v4AOAxGRO1jACLqQ4YWE45ctGyAyhlgvcu2ICIDEBG1gwGIqA8dL9XB2GJGoLccMcHeYpfj1LgiNBFdCwMQUR+6cv8vboDau+Kj/CGVACU1jajQNYldDhHZGQYgoj5kbYDm+j+9z9fDHUPCLNuMcBiMiH6OAYiojwiCcHkBRPb/9ImEAf4AGICI6GoMQER9pKi6AVX1RshlUozspxS7HJfAPiAi6ggDEFEfOdR6F2JUpBIe7jKRq3ENCf0DAQDHS7VoajaJXA0R2RMGIKI+kmvt/+HwV5+JCvREsI8CzSYBx0q1YpdDRHaEAYioj1h3gE9gAOozEonEFjjZB0REV2IAIuoDtQ1GnKuoB8AA1Ne4ICIRtYcBiKgP5LU24cYEeyPIRyFyNa7lyp3hBUEQuRoishcMQER9gMNf4hnZzw9ymRTVeiOKqhvELoeI7ESXAtD58+d7ug4ip3aI6/+IRuEmw6hIy7IDhzgMRkStuhSABg4ciOnTp+Pdd99FUxOXmCe6FmOLGUdKagFwBWixsA+IiH6uSwEoLy8Po0ePRnp6OsLCwnDffffhwIEDPV0bkVM4Wa6DocUMfy93xAb7iF2OS0q4og+IiAjoYgCKj4/H//3f/6GsrAxvvfUWysvLMWnSJIwcORJr165FZWVlT9dJ5LAOFVrW/0noHwCplBugimFcf0sAOltRB21js8jVEJE96FYTtJubGxYsWIAPP/wQq1evRn5+Ph5//HFERUVh8eLFKC8v76k6iRyWddglgcNfognxVWBAkBcEATjcOhxJRK6tWwHo0KFDeOCBBxAeHo61a9fi8ccfR0FBAXbs2IGysjLMmzevp+okckiCIFzRAB0ocjWuLaH1LlBu6x05InJtbl150dq1a/H222/jzJkzmDt3LrZu3Yq5c+dCKrXkqZiYGGzevBnR0dE9WSuRwympaURlnQHuMglGR3IDVDGNGxCAT34q5caoRASgiwFo48aN+P3vf4+lS5ciPDy83XNCQ0Pxr3/9q1vFETm6Q637f43sxw1QxWadgXe4uBYtJjPcZFwGjciVdSkA7dixA/3797fd8bESBAElJSXo378/5HI5lixZ0iNFEjkqrv9jPwaF+sJX4YY6QwvOaOowIoJ35IhcWZd+BIqLi0NVVdVVx2tqahATE9PtooicRa5tBWj2/4hNJpUgvr8/AE6HJ6IuBqCO9tOpr6+Hh4dHtwoichbaxmacragDwC0w7IX1z4ErQhNRp4bA0tPTAQASiQQrVqyAl5eX7TmTyYT9+/cjPj6+RwskclR5xZcgCEB0kBdCfLkBqj3gitBEZNWpAPTTTz8BsNwBOnbsGORyue05uVyOMWPG4PHHH+/ZCokcFIe/7E98lD+kEuDipUZodE1Q+fGONZGr6tQQ2M6dO7Fz504sWbIEX3/9te33O3fuxDfffIPXXnsNgwYN6nQRGzZsQHR0NDw8PJCcnHzNbTXeeOMNTJ48GQEBAQgICEBqaupV5wuCgBUrViA8PByenp5ITU3FuXPnOl0XUXdYZ4Bx/y/74evhjiFhfgDYB0Tk6rrUA/T222/Dz8+vRwrYvn070tPTsXLlSuTl5WHMmDGYNWsWKioq2j0/Ozsbd999N3bu3ImcnBxERUVh5syZKC0ttZ3z4osv4pVXXsGmTZuwf/9+eHt7Y9asWdy4lfpMs8mMIyVaAJwBZm8SBvgD4DAYkauTCB11NP/MggULsHnzZvj5+WHBggXXPPeTTz657gKSk5ORlJSE9evXAwDMZjOioqLw8MMP46mnnvrF15tMJgQEBGD9+vVYvHgxBEFAREQEHnvsMdtwnFarhUqlwubNm7Fo0aJfvKZOp4NSqYRWq+2xoEeu5ejFWty2fi+Unu746ZmbuAeYHfn0p4v40/YjiI/yx2cPThS7HCLqQZ35/L7uO0BKpRISicT262s9rpfRaERubi5SU1MvFySVIjU1FTk5Odd1jYaGBjQ3NyMw0NJnceHCBajV6jbXVCqVSE5O7vCaBoMBOp2uzYOoOw619v+M6+/P8GNnEvpbvlecKNOiqdkkcjVEJJbrboJ+++232/11d1RVVcFkMkGlUrU5rlKpcPr06eu6xpNPPomIiAhb4FGr1bZr/Pya1ud+LiMjA88++2xnyyfqkHV4JTGaDdD2JirQEyG+ClTWGXD0ohbjY/hnROSKutQD1NjYiIaGBtvvi4qKsG7dOnz77bc9Vtj1WLVqFbZt24ZPP/20W+sPLV++HFqt1vYoKSnpwSrJ1Vg2QLU0QHP9H/sjkUhsfVkHuTEqkcvqUgCaN28etm7dCgCora3F+PHj8fLLL2PevHnYuHHjdV8nODgYMpkMGo2mzXGNRoOwsLBrvnbNmjVYtWoVvv32W4wePdp23Pq6zlxToVDAz8+vzYOoqyxTrA1wk0owJtJf7HKoHUmtd+YOMQARuawuBaC8vDxMnjwZAPDRRx8hLCwMRUVF2Lp1K1555ZXrvo5cLkdCQgKysrJsx8xmM7KyspCSktLh61588UU899xzyMzMRGJiYpvnYmJiEBYW1uaaOp0O+/fvv+Y1iXrKgQuWD9VRkUp4yrkBqj2yBaCiSzCZr2seCBE5mS5thtrQ0ABfX18AwLfffosFCxZAKpXihhtuQFFRUaeulZ6ejiVLliAxMRHjx4/HunXroNfrkZaWBgBYvHgx+vXrh4yMDADA6tWrsWLFCrz//vuIjo629fX4+PjAx8cHEokEjz76KJ5//nkMGjQIMTExeOaZZxAREYH58+d35e0SdYp1WGU8+3/s1rBwX3jLZahrasFZTR2GhfOuL5Gr6VIAGjhwID777DPcfvvt+Oabb/CnP/0JAFBRUdHp4aOFCxeisrISK1asgFqtRnx8PDIzM21NzMXFxW12nd+4cSOMRiPuvPPONtdZuXIl/vrXvwIAnnjiCej1etx7772ora3FpEmTkJmZyX3KqE8caA1ASQxAdstNJsW4AQHYfa4KhwprGICIXNB1rwN0pY8++gi//vWvYTKZMGPGDFvzc0ZGBnbt2oWvv/66xwvtS1wHiLqqqt6AxOe/AwAcXnET/L3kv/AKEsv/fXcO//juLG4dE4FX7x4rdjlE1AM68/ndpTtAd955JyZNmoTy8nKMGTPGdnzGjBm4/fbbu3JJIqdgbaodovJl+LFzSTGtM8Eu1EAQBNs6Z0TkGroUgADLbKufz6oaP358twsicmQHLljW/7F+uJL9GhsVADepBGpdE0prGxEZ4CV2SUTUh7oUgPR6PVatWoWsrCxUVFTAbDa3ef78+fM9UhyRozlQWA0AGB8TJHIl9Es85TKM7KfE4ZJaHCysYQAicjFdCkD33HMPfvjhB/zud79DeHg4bx0TAahrasbJMss2KpwB5hiSogNaA9Al3D42UuxyiKgPdSkAff311/jyyy8xcSI3EiSyyiuuhVmwbLUQpuSMQ0eQGB2IN3Zf4IKIRC6oSwshBgQE2DYfJSKLgxc4/d3RWLfEOKupxyW9UeRqiKgvdSkAPffcc1ixYkWb/cCIXN0BLoDocIJ8FIgL8QZweQNbInINXRoCe/nll1FQUACVSoXo6Gi4u7u3eT4vL69HiiNyFIYWEw6X1AIAkri7uENJig5EQaUeB4tqkDpcJXY5RNRHuhSAuKUEUVvHLmphbDEj2EeO2GBvscuhTkiKDsS2gyW2IUwicg1dCkArV67s6TqIHNr+K/p/OCvSsVh7to6VatHUbIKHOzewJXIFXeoBAoDa2lq8+eabWL58OWpqLN/88/LyUFpa2mPFETmKg9z/y2FFBXpC5adAs0mwDWMSkfPrUgA6evQoBg8ejNWrV2PNmjWora0FAHzyySdYvnx5T9ZHZPdMZgG5hZYG2vHs/3E4EokEia3BldPhiVxHlwJQeno6li5dinPnzrXZYX3u3LnYtWtXjxVH5AhOq3WoM7TAR+HGXcUdVFLrdPgDhZwJRuQquhSADh48iPvuu++q4/369YNare52UUSOxNo8O25AAGRS9v84IuvWJbmFNWgxmX/hbCJyBl0KQAqFAjqd7qrjZ8+eRUhISLeLInIkl9f/4QaojmpomC+Unu7QG004Xnb19zYicj5dCkC33XYb/va3v6G5uRmAZQy9uLgYTz75JO64444eLZDIngmCgP3nLQHohlhugOqopFKJrX9r//lqkashor7QpQD08ssvo76+HiEhIWhsbMTUqVMxcOBA+Pr64oUXXujpGonsVn5FPar1Rni4SzE60l/scqgbklsD0D4GICKX0KV1gJRKJXbs2IG9e/fiyJEjqK+vx7hx45CamtrT9RHZNeuHZeKAQMjduryqBNkB6x28Q4WXYDIL7OcicnKdDkBmsxmbN2/GJ598gsLCQkgkEsTExCAsLAyCIHAROHIp+2zDX5z+7uiGhfvB18MNdU0tOFmmw6hIpdglEVEv6tSPrIIg4LbbbsM999yD0tJSjBo1CiNGjEBRURGWLl2K22+/vbfqJLI7giDY7gCx/8fxyaQS20a2HAYjcn6dCkCbN2/Grl27kJWVhZ9++gkffPABtm3bhiNHjuC7777D999/j61bt/ZWrUR2hf0/zie59U7e/gsMQETOrlMB6IMPPsDTTz+N6dOnX/XcjTfeiKeeegrvvfdejxVHZM/Y/+N8klvXAzpwoQYmsyByNUTUmzr1Xfvo0aOYPXt2h8/PmTMHR44c6XZRRI6A/T/OZ0SEH3wUbtA1teC0musBETmzTgWgmpoaqFSqDp9XqVS4dIlLyZPzY/+Pc3KTSZHYuqClNeASkXPqVAAymUxwc+t44phMJkNLS0u3iyKyd+z/cV7WYTAuiEjk3Do1DV4QBCxduhQKhaLd5w0GQ48URWTv2P/jvKyN0AcKa2A2C5ByPSAip9SpALRkyZJfPGfx4sVdLobIUbD/x3mN6qeEl1yG2oZmnNHUYVi4n9glEVEv6FQAevvtt3urDiKHwf4f5+YukyJhQAB2n6vC/vPVDEBETor37ok6if0/zs8abPdfYCM0kbNiACLqJPb/OD/rxqj7L9RAELgeEJEz4ndvok5i/4/zGx3pDw93KWr0RpyrqBe7HCLqBQxARJ3A/h/XIHeTInGAJeD+mF8lcjVE1BsYgIg64Rz7f1zGhIGWgLu3gOsBETkjBiCiTtjbejcgKZr9P85uYlwwAGBfQTVaTGaRqyGinsbv4ESdsDffcjdg4sBgkSuh3jaynxJ+Hm6oM7TgWKlW7HKIqIcxABFdpxaT2db/M4kByOnJpBKkxFmGwX7kMBiR02EAIrpORy5qUW9ogb+XO4ZzcTyXYL3Tt+ccG6GJnI3oAWjDhg2Ijo6Gh4cHkpOTceDAgQ7PPXHiBO644w5ER0dDIpFg3bp1V53z17/+FRKJpM1j6NChvfgOyFVY+38mxgVzfygXMaG1Dyi3+BKamk0iV0NEPUnUALR9+3akp6dj5cqVyMvLw5gxYzBr1ixUVFS0e35DQwNiY2OxatUqhIWFdXjdESNGoLy83PbYs2dPb70FciF7WgOQdXYQOb+4EG+E+XnA2GLGocJLYpdDRD1I1AC0du1aLFu2DGlpaRg+fDg2bdoELy8vvPXWW+2en5SUhJdeegmLFi3qcEd6AHBzc0NYWJjtERzMfg3qngZjC34qtnwAsv/HdUgkkiumw3MYjMiZiBaAjEYjcnNzkZqaerkYqRSpqanIycnp1rXPnTuHiIgIxMbG4je/+Q2Ki4uveb7BYIBOp2vzILrSgQs1aDYJiAzwRP9AL7HLoT5knQ6/lwsiEjkV0QJQVVUVTCYTVCpVm+MqlQpqtbrL101OTsbmzZuRmZmJjRs34sKFC5g8eTLq6uo6fE1GRgaUSqXtERUV1eWvT87J+uE3aWAwJBL2/7gSayP0sVIttA3NIldDRD1F9CbonjZnzhzcddddGD16NGbNmoWvvvoKtbW1+Pe//93ha5YvXw6tVmt7lJSU9GHF5Aj2tK7/M4HDXy4nTOmBuBBvCAKQc57T4YmchWgBKDg4GDKZDBqNps1xjUZzzQbnzvL398fgwYORn5/f4TkKhQJ+fn5tHkRWVfUGnCq3DItOiGMDtCuy3gXiMBiR8xAtAMnlciQkJCArK8t2zGw2IysrCykpKT32derr61FQUIDw8PAeuya5lpzWRfCGhfsh2Kfj5ntyXtbp8GyEJnIebmJ+8fT0dCxZsgSJiYkYP3481q1bB71ej7S0NADA4sWL0a9fP2RkZACwNE6fPHnS9uvS0lIcPnwYPj4+GDhwIADg8ccfx6233ooBAwagrKwMK1euhEwmw9133y3OmySHd7n/h3d/XFVKbBCkEuB8pR7l2kaEKz3FLomIuknUALRw4UJUVlZixYoVUKvViI+PR2Zmpq0xuri4GFLp5ZtUZWVlGDt2rO33a9aswZo1azB16lRkZ2cDAC5evIi7774b1dXVCAkJwaRJk7Bv3z6EhIT06Xsj5yAIAnafs67/w/4fV6X0cseofkocuajF3vxq3JkQKXZJRNRNEkEQBLGLsDc6nQ5KpRJarZb9QC6uqFqPqS9lw10mweEVM+GtEPVnBhLR6szT2JhdgAVj+2HtwnixyyGidnTm89vpZoER9STr7u9j+wcw/Lg46wKYu/OrYDbz50YiR8cARHQNe/IrAXD1ZwISowPgJZehss6Ak+VcLJXI0TEAEXWgxWS27QI+aRADkKtTuMlsyyD8cLZS5GqIqLsYgIg68FNJLXRNLQjwcseYSH+xyyE7MHWwZTLFD2cYgIgcHQMQUQesH3KTB4VAJuX2FwRMGxIKAMgtvgRtI7fFIHJkDEBEHcg+WwHg8k/9RFGBXogN8YbJLOBHrgpN5NAYgIjaUVHXhOOllkbXKQxAdIVpgy13gbI5DEbk0BiAiNqx+6zlp/tR/ZQI8eX2F3TZ1CGtfUBnK8Fl1IgcFwMQUTuyW2f5cPiLfi45JhAe7lKodU04o6kTuxwi6iIGIKKfMZkF7D5nCUDThjAAUVse7jKkxFqmw3MYjMhxMQAR/cyRi7WobWiGn4cb4qP8xS6H7JD1zmD2mQqRKyGirmIAIvqZ7Cumv7vJ+E+ErmadDn+o8BLqDS0iV0NEXcHv7kQ/80PrT/VTOfxFHYgO9kZ0kBdazAL2cjo8kUNiACK6QnW9AUdLtQCAaWyApmuw3gViHxCRY2IAIrrC7nNVEARgWLgfQv08xC6H7Ji1D2gXp8MTOSQGIKIrWJtaOfuLfskNsUGQu0lRWtuI/Ip6scshok5iACJqZTYL2NW6+zuHv+iXeMovT4fPOs3ZYESOhgGIqNXRUi1q9Eb4KtwwbkCA2OWQA7hpuAoA8O0JtciVEFFnMQARtco6pQEATB4cDHdOf6frYA1AP5XUokLXJHI1RNQZ/C5P1GrHSUsAsn6oEf0SlZ8HxkT5QxCA705xGIzIkTAAEQEoqWnAaXUdZFIJprdObya6HjNbA/OOkxwGI3IkDEBEuHz3Jyk6AP5ecpGrIUcya4QlAO3Nr+aq0EQOhAGICFcOf4WJXAk5mrgQH8QEe8NoMuMHLopI5DAYgMjlaRuacaCwBgBw0zD2/1DnSCQS2zDYtxwGI3IYDEDk8naeqYDJLGCIyhf9g7zELocc0MzWYbDvT1eg2WQWuRoiuh4MQOTydpzi7C/qnvioAAT7KFDX1IL952vELoeIrgMDELk0Q4vJ1reRygBEXSSTSpA6zDJ7kMNgRI6BAYhc2v7zNag3tCDUV4HR/ZRil0MOzDoMtuOkhpujEjkABiByadbZXzOGqSCVSkSuhhzZhLhgeMllKNc24XipTuxyiOgXMACRyxIEAd+19v/M5PAXdZOHuwzThlg20eUwGJH9YwAil3WiTIdybRO85DKkxAWJXQ45AWsjfeZxBiAie8cARC7r29bhrymDQuDhLhO5GnIGNw5VwV0mwbmKepzT1IldDhFdAwMQuaxvT1h+SufsL+opSk93TBlkGQb74mi5yNUQ0bUwAJFLyq+ox2l1HdxlEq7+TD3qljHhAIAvjpZxNhiRHWMAIpf01THLT+eTBgZD6eUucjXkTFKHqSB3k6KgUo8zHAYjslsMQOSSrAFo7qhwkSshZ+Pr4Y6pg1uHwY5wGIzIXokegDZs2IDo6Gh4eHggOTkZBw4c6PDcEydO4I477kB0dDQkEgnWrVvX7WuS68mvqLMNf83k7u/UC24ZbQnWXx4r5zAYkZ0SNQBt374d6enpWLlyJfLy8jBmzBjMmjULFRUV7Z7f0NCA2NhYrFq1CmFh7X9wdfaa5Hq+PGppfp48KITDX9QrZgxTQeEmxYUqPRdFJLJTogagtWvXYtmyZUhLS8Pw4cOxadMmeHl54a233mr3/KSkJLz00ktYtGgRFApFj1yTXM+Xx8oAADdz+It6iY/CzTa78LPDpSJXQ0TtES0AGY1G5ObmIjU19XIxUilSU1ORk5NjN9ck53JOU4ezmnq4yySc/k696vb4fgCAz4+UwWTmMBiRvREtAFVVVcFkMkGlavshpFKpoFZ3bRXVrl7TYDBAp9O1eZBz+rK1+XnKoBAoPTn8Rb1nyuAQBHi5o7LOgB8LqsQuh4h+RvQmaHuQkZEBpVJpe0RFRYldEvWSL1sXp7t5NIe/qHfJ3aS2v2ef/sRhMCJ7I1oACg4Ohkwmg0ajaXNco9F02ODcW9dcvnw5tFqt7VFSUtKlr0/27aymDucq6iGXSTn8RX3i9rGWYbBvjqvRaDSJXA0RXUm0ACSXy5GQkICsrCzbMbPZjKysLKSkpPTpNRUKBfz8/No8yPlYtyaYMjgYfh4c/qLeN65/AKICPaE3mrhDPJGdEXUILD09HW+88Qa2bNmCU6dO4f7774der0daWhoAYPHixVi+fLntfKPRiMOHD+Pw4cMwGo0oLS3F4cOHkZ+ff93XJNckCIJt8UMOf1FfkUgkuH1sJADgo9yLIldDRFdyE/OLL1y4EJWVlVixYgXUajXi4+ORmZlpa2IuLi6GVHo5o5WVlWHs2LG2369ZswZr1qzB1KlTkZ2dfV3XJNd0RlOH/Ip6yN2kSOXeX9SH7kqIxCtZ57AnvwqltY3o5+8pdklEBEAicJnSq+h0OiiVSmi1Wg6HOYmMr0/htR/OY+ZwFV5fnCh2OeRi7n59H3LOVyP9psH444xBYpdD5LQ68/nNWWDk9ExmAZ+1zsJZMC5S5GrIFf0qyfL37sPcEpi5JhCRXWAAIqe3N78KGp0B/l7umD40ROxyyAXNHhEOX4UbSmoasf9CjdjlEBEYgMgFfJJnaT69dXQEFG4ykashV+Qpl+GWMREAgO0Hi0WuhogABiBycvWGFmSesEw/viOBw18knkVJlgVWvzquRo3eKHI1RMQARE7tq2PlaGo2IzbEG2MilWKXQy5sTJQ/RvVTwthixke5XGyVSGwMQOTUrMNfd4yLhEQiEbkacnW/Se4PAHhvfzGboYlExgBETqukpgH7ztdAIgHmt25JQCSm2+Ij4KtwQ1F1A/Zyg1QiUTEAkdOyTn1PiQ3i4nNkF7zkblgwzhLG391XJHI1RK6NAYickiAI+KQ1AN3BtX/IjvzmhgEAgB0nNbh4qUHkaohcFwMQOaW84lpcqNLDSy7D7JFhYpdDZDNY5YuJA4NgFoCtObwLRCQWBiByStbm59kjw+CtEHXLO6Kr/H5iDADggwPF0BtaRK6GyDUxAJHTaTSa8PmRMgAc/iL7NH1IKGKDvVHX1MJd4olEwgBETufLY+Woa2pBVKAnUmKDxC6H6CpSqQRpE6MBAG/vvcAp8UQiYAAip/PBActWA4uS+kMq5do/ZJ/uSIiEn4cbCqsb8O1JtdjlELkcBiByKmfUdcgtugQ3qQR3JXL4i+yXl9wNi1OiAQD/zC6AIPAuEFFfYgAip2K9+5M6TIVQXw+RqyG6trSJ0fBwl+LoRS325HNhRKK+xABETqOp2WSb/XV365YDRPYsyEeBRUmWv6v/3FkgcjVEroUBiJzG54fLoGtqQWSAJyYPDBa7HKLrcu+UWLhJJcg5X43coktil0PkMhiAyCkIgoAtOYUAgN/eMIDNz+QwIvw9bdtjrPvurMjVELkOBiByCnnFtThRpoPCTYqFiVFil0PUKQ/fOAjuMgl2n6vCvvPVYpdD5BIYgMgpbG29+3PrmAgEeMvFLYaok6ICvbAwyRLcX/72DGeEEfUBBiByeJV1Bnx1rBwAsKR1WjGRo3lo+iAo3KQ4WHgJP5ytFLscIqfHAEQO7739RWg2CYiP8seoSKXY5RB1SZjSA79r3Sl+deYZmLg6NFGvYgAih9bUbMK7+yw7av9+UozI1RB1zwPTB8LXww2nynX4KLdE7HKInBoDEDm0zw+XoareiHClB+aMDBO7HKJuCfSW45EZgwAAL31zFvXcKZ6o1zAAkcMSBAH/2nMBALB0QjTcZfzrTI5vcUo0ooO8UFVvwMbsfLHLIXJa/MQgh7UnvwpnNHXwksuwaDxXfibnIHeT4um5wwAAb+y+gPOV9SJXROScGIDIYb2+6zwA4FeJUVB6uotcDVHPuWm4ClMGh8DYYsafPz3OafFEvYABiBzS8VItdp+rgkwqwR/Y/ExORiKR4IX5I+HhLkXO+Wp8nFcqdklETocBiBzSxmzLxpG3jg5HVKCXyNUQ9byoQC88MmMwAOCFL0+iqt4gckVEzoUBiBzOhSo9vjpuWfjwf6bFiVwNUe+5Z3IMhob54lJDM576+BiHwoh6EAMQOZzXdxVAEIAZQ0MxNMxP7HKIeo27TIq1v4qHXCbFd6c02H6QawMR9RQGIHIopbWN+Cj3IgDgft79IRcwPMIPj8+yDIX97YuTKKzSi1wRkXNgACKH8s+d+Wg2CZgQF4TE6ECxyyHqE/dMisUNsYFoMJrwP+/motFoErskIofHAEQOo6y2Ef8+ZBkCsK6WS+QKpFIJ1i0ci2AfOU6r6/DUJ0fZD0TUTQxA5DD+mW25+5MSG4Tk2CCxyyHqU2FKD2z49Ti4SSX4z+Ey2yroRNQ1DEDkEEprG20NoI+k8u4Puabk2CD8+WbLKtEvfHUKXx4tF7kiIsdlFwFow4YNiI6OhoeHB5KTk3HgwIFrnv/hhx9i6NCh8PDwwKhRo/DVV1+1eX7p0qWQSCRtHrNnz+7Nt0C9bN2Os7a7Pzfw7g+5sKUTorE4ZQAEAfjT9sP4saBK7JJ6ldkswNhiRoOxBcYWs9jlkBNxE7uA7du3Iz09HZs2bUJycjLWrVuHWbNm4cyZMwgNDb3q/B9//BF33303MjIycMstt+D999/H/PnzkZeXh5EjR9rOmz17Nt5++23b7xUKRZ+8H+p5+RV1+DjPMvPridlDRK6GSFwSiQQrbx2ByjoDvj6uxrIth/DW0iSHHhbWNTXjZJkOp8p1KKisR1F1A0prG3FJb0RtYzOubHdyl0mg9HRHmNID/fw9MVjli6Fhfhjb3x8R/p7ivQlyOBJB5E665ORkJCUlYf369QAAs9mMqKgoPPzww3jqqaeuOn/hwoXQ6/X44osvbMduuOEGxMfHY9OmTQAsd4Bqa2vx2WefdakmnU4HpVIJrVYLPz+uMyO2/3knF5kn1Jg5XIXXFyeKXQ6RXWhqNuEPWw5ib341PNyleP13iZgyOETssq6LrqkZhwprkFNQjZzz1ThRpkNPfBJFBXpi0sBgzBwehpS4IHi4y7p/UXIonfn8FvUOkNFoRG5uLpYvX247JpVKkZqaipycnHZfk5OTg/T09DbHZs2adVXYyc7ORmhoKAICAnDjjTfi+eefR1BQ+z8hGQwGGAyXl5nX6XRdfEfU0w6X1CLzhBpSCfD4LN79IbLycJfhX0uScP+7udh5phL3bDmEjAWjcEdCpNiltatC14TME2p8ebQcBwtrYP5Z4Onn74nhEX4YovJF/yAvRAZ4IthHAX8vd3i4yyCTSNBiEtDQ3IIavRHltU0oudSAs5o6nCjT4USZDiU1jfjgQAk+OFACH4UbbhkdjjsTIpEwIAASiUScN052S9QAVFVVBZPJBJVK1ea4SqXC6dOn232NWq1u93y1Wm37/ezZs7FgwQLExMSgoKAATz/9NObMmYOcnBzIZFf/RJCRkYFnn322B94R9SRBEPD8FycBALePjcRgla/IFRHZFw93GV77XSIe3f4TvjqmxmMfHsHJch2WzxkKN5n4LZ7l2kZkHlfjq2PlOFR0qc1dnuggL6TEWXr6kmOCEKb0uK5rKuGOcKUnRkQo2xyvN7TgYGENvj9VgW9PqqHRGbDtYAm2HSzBsHA//H5iNG6Lj4DCjXeFyEL0HqDesGjRItuvR40ahdGjRyMuLg7Z2dmYMWPGVecvX768zV0lnU6HqKioPqmVOvbVMTUOFV2Ch7vUthIuEbUld5Ni/d3j8I+Qs3j1+3z8a88F5BVfwkt3jsHAUJ8+r6e0thFfHyvHV8fKkVdc2+a5cf39MXdUOGaNCOvxTYx9FG6YPiQU04eE4tnbRuBAYQ0+yr2IL4+W41S5Dv/70VH8Y8dZPJI6CHeMi7SLgEjiEjUABQcHQyaTQaPRtDmu0WgQFhbW7mvCwsI6dT4AxMbGIjg4GPn5+e0GIIVCwSZpO9PUbMKqzFMAgPumxCFcyeZGoo5IpRI8NnMIhoX74cmPjuKn4lrMfWU3/njjQPxhUiw85b1716O4ugFfHy/HV8fVOFJS2+a5xAEBmDsqHLNHhvVZk7JUKsENrTNG/3LzMHxwoASbf7yAMm0Tnvz4GF7bdR6P3TQEc0aGQSrl0JirEjUAyeVyJCQkICsrC/PnzwdgaYLOysrCQw891O5rUlJSkJWVhUcffdR2bMeOHUhJSenw61y8eBHV1dUIDw/vyfKpF23+sRAlNY1Q+Slw39RYscshcghzR4UjPsofT31yDLvOVmLNt2fxzr4iPHzjICwY1w9e8p75li8IAo6VarHjpAbfntDgjKbO9pxEAiRFB+Lm1js91zu01Vv8veS4f1oc0iZG4919RdiwMx/nK/V48P08jOznh2dvG4mEAQGi1kjiEH0W2Pbt27FkyRK89tprGD9+PNatW4d///vfOH36NFQqFRYvXox+/fohIyMDgGUa/NSpU7Fq1SrcfPPN2LZtG/7+97/bpsHX19fj2WefxR133IGwsDAUFBTgiSeeQF1dHY4dO3Zdd3o4C0xc5dpGzHj5BzQYTVhz1xjcaadNnUT2ShAEfHa4FGu+OYvS2kYAgK+HG+4YF4nZI8OQMCAA7p0YAjKbBZyrqMfBwhrkFl1CTkE11Lom2/MyqQTjowMxd1QYZo0MQ6ivuKHnWuqamvHm7gt4c/d56I0mSCTAr8f3xxOzh0Lp6S52edRNnfn8Fj0AAcD69evx0ksvQa1WIz4+Hq+88gqSk5MBANOmTUN0dDQ2b95sO//DDz/EX/7yFxQWFmLQoEF48cUXMXfuXABAY2Mj5s+fj59++gm1tbWIiIjAzJkz8dxzz13VPN0RBiBxPfR+Hr44Wo5x/f3x0f9M4C1qoi4ytJjw/v5ivL23EMU1Dbbjfh5uGBPlj+ERfhgQ6I1gHzl8PSwf/i1mM2r0RlTVG1Fcrcf5Kj2OlNRC19TS5tpechmmDg7BTcNVmD4kFAHe8j59b91VXW/Aqq9P48NcyxpjwT4KrLx1OG4ZHc4ZYw7M4QKQvWEAEs/e/Cr85s39kEqA/z486aqZHkTUeWazgF3nKvH54TJkn61Ejd7Y6Wt4yWUY298fCQMCkRQdgKToQKdYZyenoBp//uwYzlfqAQDTh4Rg9Z2j7fouFnWMAaibGIDEYWgxYc7/7cb5Sj2WTojGX28bIXZJRE7HZBZwokyLE2U6nCzToVzbiKp6I+oNLZDAMpzl7+WOIG8FogK9EBvsjaHhvhge7ue0M6cMLSZsyj6PDdn5MLaYEegtx6oFozBzRMeTa8g+MQB1EwOQONZ+ewavfJ+PYB8Fsh6byvF4IupTZzV1eGTbYZwqtyyGuygpCs/cMhzeCqdcMcYpdebz2znjPDmc02od/pldAAD427wRDD9E1OcGq3zx2YMTcN/UWEgkwLaDJZj7ym4cvVgrdmnUCxiASHQms4AnPz6GFrOAmcNVmDOSt52JSBwKNxmWzxmG9++5ARFKDxRVN+DOjTl4d18ROGDiXBiASHSv7zqPIyW18FW44W/zRnIGBhGJLiUuCF8/OgU3DVfBaDLjL58dx5+2H4be0PLLLyaHwABEojpZpsPaHWcAAM/cOlz0RdOIiKyUnu54/XcJeHruUMikEnx2uAzzNuzFuSsWfiTHxQBEojG0mJD+78NoNgm4abgKd3HBQyKyMxKJBPdOicMHy25AqK8C+RX1uG39XvzncKnYpVE3MQCRaF7KPIPT6joEecuRsWAUh76IyG6NjwnEl3+cjAlxQWhsNuGRbYfx50+PwdBiErs06iIGIBLFdyc1eHPPBQDAqjtGI9iHm9ESkX0L8VXgnT8k4+EbBwIA3ttfjF+9tg9lrduNkGNhAKI+V1bbiMc/OgIASJsYjZuGX98WJUREYpNJJXhs5hC8nZYEpac7jpTU4tZX9+DHgiqxS6NOYgCiPmVsMeOh9/NQ29CMUf2UeGrOULFLIiLqtOlDQvHfhyZheLgfqvVG/PbN/Xh9VwGnyjsQBiDqU8/+9wTyimvh6+GG9b8eC4Wb4+8lRESuqX+QFz6+fwIWjOsHswD8/avTeOj9n1DPqfIOgQGI+sy2A8V4b38xJBLglUVjMSDIW+ySiIi6xVMuw8t3jcFz80bAXSbBl8fKMX/DXhRU1otdGv0CBiDqE/vPV2PFf04AAB67aTCmDw0VuSIiop4hkUjwu5RobLs3BSo/y1T5eev3IvO4WuzS6BoYgKjXna+sx33v5sJoMmPuqDA8MG2g2CUREfW4hAEB+O/DkzA+JhD1hhb8z7u5WJ15GiYz+4LsEQMQ9arqegN+v/kgahuaER/lj7W/iodUyvV+iMg5hfp64L17knHPpBgAwMbsAix56wBq9EaRK6OfYwCiXlPX1Iylbx9EYXUDIgM88cbiRHi4s+mZiJybu0yKv9wyHK/ePRZechn25Ffh1lf34HBJrdil0RUYgKhXNDWbcM+WQzhWqkWQtxxbfj8eIb5c7JCIXMetYyLw2YMTERPsjdLaRty58Ue8ses8zBwSswsMQNTjmppNuO+dXOy/UAMfhRu2/H484kJ8xC6LiKjPDVb54j8PTcTNo8LRYhbwwlencM/WQxwSswMMQNSjmppNWLb1EH44WwlPdxn+tSQRI/spxS6LiEg0fh7uWP/rsXjh9pGQu0nx/ekKzP2/3dh/vlrs0lwaAxD1mHpDC/6w5SB2n6uCp7sMb6clITk2SOyyiIhEJ5FI8JvkAfjsgYmIDfGGWteEu9/Yh9WZp2FsMYtdnktiAKIeUV1vwN2v78Pe/Gp4yWXYnJaEGxh+iIjaGB7hh/8+NAl3jIuEWbDMEpu/YS/OqOvELs3lMABRt12o0uPOTTk4VqpFoLccHyy7gXd+iIg64K1ww8u/GoNNvx2HAC93nCzX4dZX97BBuo8xAFG3/FhQhfkb9uJClR79/D3x4f+kYEyUv9hlERHZvdkjw/HNn6Zg+pAQGE1mvPDVKdz1Wg7Oang3qC8wAFGXCIKAt/ZcwOJ/HYC2sRlj+/vj0wcncLYXEVEnhPp64K2lSfj77aPgLZcht+gSbn5lN17+9gyamk1il+fUJIIg8H7bz+h0OiiVSmi1Wvj5+Yldjt2pN7TgyY+P4suj5QCAefERWH3HaC5ySETUDWW1jVjxnxP47pQGABAb7I0Xbh+FlDi2FFyvznx+MwC1gwGoY7lFNfjT9iMormmAm1SCP988DEsnREMi4fYWRETdJQgCMo+rsfLzE6ioMwAA5o4Kw/I5wxAV6CVydfaPAaibGICu1tRswqvfn8PG7AKYBaCfvydeuXssEgYEiF0aEZHT0TY246VvTuP9/cUwC4BcJkXapGg8NH0gfD3cxS7PbjEAdRMDUFv7z1dj+afHcL5SDwBYMLYf/jpvBPz4j5CIqFedVuvw/BensCe/CgAQ5C3HwzcOxKLx/dl20A4GoG5iALIorNLjxW9O46tjagBAiK8Cf7ttBOaMChe5MiIi1yEIAr4/XYEXvjpl+0E01FeB+6fF4W4GoTYYgLrJ1QPQJb0Rr3x/Du/uK0KzSYBUAixM6o+n5gyF0pN3fYiIxNBsMmP7wRL8c2c+yrRNACw/mN43JRYLk6I4NAYGoG5z1QBUXW/AO/uK8K89F1DX1AIAmDo4BMvnDsXQMNf5/0BEZM+MLWZ8lHsRG3bmo7S2EQDgLZfhrsQoLE4ZgFgXXo6EAaibXC0AFVTW4197LuDj3IswtO5JMyzcD0/PHYrJg0JEro6IiNpjbDHj47yLeHP3eRS0Do0Blh9c7x7fH9OHhkDh5lrDYwxA3eQKAaip2YSdpyvwYe5FfH+6wnZ8TKQS90yOxdxR4ZBJObWdiMjeCYKAPflV2Ly3EN+fqYD1U93fyx03jwrHgnH9MK5/gEssV8IA1E3OGoDMZgGHii7h058u4suj5dC1DnNJJMCMoSosmxyD8TGBLvGPhIjIGRVV6/He/mJ89lOpbR0hAIgK9MRNw8KQOiwUSTGBcJc550YQDEDd5EwBSNvYjD3nqrDzTAWyz1Siqv7yP4hwpQfmxffDrxIjXXrMmIjI2ZjMAn4sqMKnP5Ui87gaDcbL22r4erhh6uAQTBkUguTYQPQP9HKaH3wZgLrJkQNQdb0BecW1yC26hNyiGuQV18J0xe7CPgo3zB4ZhgVj++GG2CBIOcxFROTUGowt2HW2Ct+d0mDn6QpU641tng/z80BybCDGxwRidD9/DA7zcdjeIYcLQBs2bMBLL70EtVqNMWPG4NVXX8X48eM7PP/DDz/EM888g8LCQgwaNAirV6/G3Llzbc8LgoCVK1fijTfeQG1tLSZOnIiNGzdi0KBB11WPIwQgs1nAxUuNOK3W4aymDqfVdTheqkVhdcNV5w4M9cH0ISGYNiQUidEBDvsXm4iIusdkFnC4pBbfn9Zg3/kaHL1Yi2ZT2xjgLpNgsMoXIyOUGBLmi7hQH8SFeCNC6Wn3PzQ7VADavn07Fi9ejE2bNiE5ORnr1q3Dhx9+iDNnziA0NPSq83/88UdMmTIFGRkZuOWWW/D+++9j9erVyMvLw8iRIwEAq1evRkZGBrZs2YKYmBg888wzOHbsGE6ePAkPD49frMkeAlBTswmVdQZU1htQWWdAeW0jimsaUVzTgJKaBhTXNKCxg52CB6t8MK5/AMYNCEBKbBD3jyEionY1Gk3IK76E/RdqkFd0CcfLtKhtaG73XA93KaKDvNHP3xPh/h4IV3oiwvpfpSdUSoXoP2A7VABKTk5GUlIS1q9fDwAwm82IiorCww8/jKeeeuqq8xcuXAi9Xo8vvvjCduyGG25AfHw8Nm3aBEEQEBERgcceewyPP/44AECr1UKlUmHz5s1YtGjRL9bUWwHofGU9TpXXQdfUjLqmZtQ1taCuqaX19y3QNjajqjXwWNfhuRa5TIqBoT4YEuaLIWG+GBbuh/hIfyi9uBgWERF1niBYRhdOlGlxvFSHcxV1KKjUo6haf9WdovZ4yWUI8JLD38u9zX/9PN3grXCDj8IN3nLLrweG+mBgaM/2n3bm89utR79yJxmNRuTm5mL58uW2Y1KpFKmpqcjJyWn3NTk5OUhPT29zbNasWfjss88AABcuXIBarUZqaqrteaVSieTkZOTk5LQbgAwGAwyGy83BWq0WgOV/ZE/6KKcA63fmX/f57m5ShPjIEeyjQKivApGBXogM8ERkQOt//T3h9vNO/pZG6HSNPVo3ERG5DqUbMKG/Nyb09wZg2fqoxWTGxdpGFFc3oFzXCE2twfJfrQFqXSPKdQY0t5hRbwDq64CS6/g6aZOi8dhNQ3q0duvn9vXc2xE1AFVVVcFkMkGlUrU5rlKpcPr06XZfo1ar2z1frVbbnrce6+icn8vIyMCzzz571fGoqKjreyO96LzYBRAREfWCv60D/tZL166rq4NSqbzmOaIGIHuxfPnyNneVzGYzampqEBQUZHdTA3U6HaKiolBSUmK3DdrOjn8G4uOfgX3gn4P4+GfQliAIqKurQ0RExC+eK2oACg4Ohkwmg0ajaXNco9EgLCys3deEhYVd83zrfzUaDcLDw9ucEx8f3+41FQoFFApFm2P+/v6deSt9zs/Pj3/ZRcY/A/Hxz8A+8M9BfPwzuOyX7vxYiboUpFwuR0JCArKysmzHzGYzsrKykJKS0u5rUlJS2pwPADt27LCdHxMTg7CwsDbn6HQ67N+/v8NrEhERkWsRfQgsPT0dS5YsQWJiIsaPH49169ZBr9cjLS0NALB48WL069cPGRkZAIBHHnkEU6dOxcsvv4ybb74Z27Ztw6FDh/D6668DACQSCR599FE8//zzGDRokG0afEREBObPny/W2yQiIiI7InoAWrhwISorK7FixQqo1WrEx8cjMzPT1sRcXFwMqfTyjaoJEybg/fffx1/+8hc8/fTTGDRoED777DPbGkAA8MQTT0Cv1+Pee+9FbW0tJk2ahMzMzOtaA8jeKRQKrFy58qohO+o7/DMQH/8M7AP/HMTHP4OuE30dICIiIqK+5pzbwRIRERFdAwMQERERuRwGICIiInI5DEBERETkchiAHFRhYSH+8Ic/ICYmBp6enoiLi8PKlSthNBrFLs2pbdiwAdHR0fDw8EBycjIOHDggdkkuJSMjA0lJSfD19UVoaCjmz5+PM2fOiF2WS1u1apVt+RHqW6Wlpfjtb3+LoKAgeHp6YtSoUTh06JDYZTkMBiAHdfr0aZjNZrz22ms4ceIE/vGPf2DTpk14+umnxS7NaW3fvh3p6elYuXIl8vLyMGbMGMyaNQsVFRVil+YyfvjhBzz44IPYt28fduzYgebmZsycORN6vV7s0lzSwYMH8dprr2H06NFil+JyLl26hIkTJ8Ld3R1ff/01Tp48iZdffhkBAQFil+YwOA3eibz00kvYuHEjzp/nFqq9ITk5GUlJSVi/fj0Ay6rlUVFRePjhh/HUU0+JXJ1rqqysRGhoKH744QdMmTJF7HJcSn19PcaNG4d//vOfeP755xEfH49169aJXZbLeOqpp7B3717s3r1b7FIcFu8AORGtVovAwECxy3BKRqMRubm5SE1NtR2TSqVITU1FTk6OiJW5Nq1WCwD8ey+CBx98EDfffHObfxPUdz7//HMkJibirrvuQmhoKMaOHYs33nhD7LIcCgOQk8jPz8err76K++67T+xSnFJVVRVMJpNthXIrlUoFtVotUlWuzWw249FHH8XEiRPbrARPvW/btm3Iy8uzbVFEfe/8+fPYuHEjBg0ahG+++Qb3338//vjHP2LLli1il+YwGIDszFNPPQWJRHLNx+nTp9u8prS0FLNnz8Zdd92FZcuWiVQ5Ud968MEHcfz4cWzbtk3sUlxKSUkJHnnkEbz33ntOsb2QozKbzRg3bhz+/ve/Y+zYsbj33nuxbNkybNq0SezSHIboe4FRW4899hiWLl16zXNiY2Ntvy4rK8P06dMxYcIE24aw1POCg4Mhk8mg0WjaHNdoNAgLCxOpKtf10EMP4YsvvsCuXbsQGRkpdjkuJTc3FxUVFRg3bpztmMlkwq5du7B+/XoYDAbIZDIRK3QN4eHhGD58eJtjw4YNw8cffyxSRY6HAcjOhISEICQk5LrOLS0txfTp05GQkIC33367zaax1LPkcjkSEhKQlZWF+fPnA7D8BJaVlYWHHnpI3OJciCAIePjhh/Hpp58iOzsbMTExYpfkcmbMmIFjx461OZaWloahQ4fiySefZPjpIxMnTrxqCYizZ89iwIABIlXkeBiAHFRpaSmmTZuGAQMGYM2aNaisrLQ9xzsSvSM9PR1LlixBYmIixo8fj3Xr1kGv1yMtLU3s0lzGgw8+iPfffx//+c9/4Ovra+u/UiqV8PT0FLk61+Dr63tVz5W3tzeCgoLYi9WH/vSnP2HChAn4+9//jl/96lc4cOAAXn/9dY4EdAIDkIPasWMH8vPzkZ+ff9UQAFc26B0LFy5EZWUlVqxYAbVajfj4eGRmZl7VGE29Z+PGjQCAadOmtTn+9ttv/+LQMZEzSUpKwqefforly5fjb3/7G2JiYrBu3Tr85je/Ebs0h8F1gIiIiMjlsGmEiIiIXA4DEBEREbkcBiAiIiJyOQxARERE5HIYgIiIiMjlMAARERGRy2EAIiIiIpfDAEREREQuhwGIiIiIXA4DEBEREbkcBiAiIiJyOQxARERE5HL+H8JzOolHZUD3AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.kdeplot(Dat)\n",
    "\n",
    "for i in range(len(Dat)):\n",
    "    plt.axvline(x=Dat[i], ymax=0.05, c='r')                # Plot the observations as a rugplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Deprecated measures of asymmetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sk1 = (Mean - Mode) / S2\n",
    "Sk2 = (Mean - Median) / S2\n",
    "S75 = (np.quantile(Dat, 0.75) + np.quantile(Dat, 0.25) - 2*np.median(Dat)) / IQR\n",
    "S90 = (np.quantile(Dat, 0.9) + np.quantile(Dat, 0.1) - 2*np.median(Dat)) / (np.quantile(Dat, 0.9) - np.quantile(Dat, 0.1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(0.2276),\n",
       " np.float64(0.2276),\n",
       " np.float64(0.2308),\n",
       " np.float64(0.3895))"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Sk1.round(4), Sk2.round(4), S75.round(4), S90.round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kurtosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "m4 = ((Dat-Mean)**4).mean()\n",
    "\n",
    "k4 = N**2*((N+1)*m4 - 3*(N-1)*m2**2) / ((N-1)*(N-2)*(N-3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coefficient of kurtosis based on central moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "g2_byhand = m4/m2**2\n",
    "g2 = kurtosis(Dat, fisher=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(4.8753), np.float64(4.8753))"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g2_byhand.round(4), g2.round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coefficient of kurtosis based on $k$-statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "G2_byhand = k4/k2**2 + 3\n",
    "G2 = kurtosis(Dat, fisher=False, bias=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(7.2795), np.float64(7.2795))"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G2_byhand.round(4), G2.round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Excess of kurtosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "g2_excess_byhand = g2_byhand - 3\n",
    "g2_excess = kurtosis(Dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(1.8753), np.float64(1.8753))"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g2_excess_byhand.round(4), g2_excess.round(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "G2_excess_byhand = G2_byhand - 3\n",
    "G2_excess = kurtosis(Dat, bias=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(4.2795), np.float64(4.2795))"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G2_excess_byhand.round(4), G2_excess.round(4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}