exercises/Week-1-intro-23_sol.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 1 - Normalization and data types\n",
    "\n",
    "This first week is going to be a warm-up week where we conider different statistical tools data types.\n",
    "\n",
    "Lets first import the libraries that we are going to need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: Sample statistics and normalization\n",
    "\n",
    "Let $X$ and $Y$ be two normal distributions of discrete variables with ($\\mu_X=25,\\sigma_X=5$) and ($\\mu_Y=60,\\sigma_Y=10$) that describe the *Age* and the *Weight* of a population. From $X$ and $Y$ sample the *Age* and the *Weight* of $n=500000$ individuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.random.seed(1)\n",
    "# TODO: use np.random.normal to initialize X and Y and round() function for discretization of the data\n",
    "X = np.random.normal(25,5,500000)\n",
    "Y = np.random.normal(60,10,500000)\n",
    "X = X.round()\n",
    "Y=Y.round()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)** Find the mean, the median and the mode for $X,Y$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean: 25.005264\n",
      "Median: 25.0\n",
      "Mode: 25.0\n",
      "\n",
      "Mean: 60.002704\n",
      "Median: 60.0\n",
      "Mode: 60.0\n"
     ]
    }
   ],
   "source": [
    "def count_dictionary(X):\n",
    "    \"\"\"\n",
    "    Input: certain list of data   \n",
    "    Output: counts of list elements as a dictionary \n",
    "    Example: Input: [0,1,0,2,3,2,2] => Output: {0:2, 1:1, 2:3, 3:1}\n",
    "    \"\"\"\n",
    "    # TODO: the calculation process of counts\n",
    "    # We need counts also for later when we are computing the empirical PMF\n",
    "    count={}\n",
    "    for i in X:\n",
    "        if i not in count:\n",
    "            count[i]=1\n",
    "        else:\n",
    "            count[i]+=1\n",
    "    return count\n",
    "\n",
    "def compute_mean(X):\n",
    "    \"\"\"\n",
    "    Input: certain list of data   \n",
    "    Output: mean of this list   \n",
    "    \"\"\"\n",
    "    # TODO: the calculation process of mean\n",
    "    return sum(X)/len(X)\n",
    "\n",
    "def compute_median(X):\n",
    "    \"\"\"\n",
    "    Input: certain list of data   \n",
    "    Output: median of this list\n",
    "    Tip:Note: If the number of data values is odd, returns the exact middle value. \n",
    "    If the number of data values is even, returns the average of the two middle values.   \n",
    "    \"\"\"\n",
    "    # TODO: the calculation process of median\n",
    "    if len(X) % 2 == 1:\n",
    "        median=sorted(X)[int(len(X)//2)]\n",
    "    else:\n",
    "        median=(sorted(X)[len(X)//2]+sorted(X)[len(X)//2-1])/2\n",
    "    return median\n",
    "def compute_mode(X):\n",
    "    \"\"\"\n",
    "    Input: certain list of data   \n",
    "    Output: element with highest frequency \n",
    "    Hint: You can use function count_dictionary   \n",
    "    \"\"\"\n",
    "    # TODO: the calculation process of mode  \n",
    "    count=count_dictionary(X)\n",
    "    mode=max(count, key=count.get)\n",
    "    return mode\n",
    "\n",
    "\n",
    "x_mean = compute_mean(X)\n",
    "x_median = compute_median(X)\n",
    "x_mode = compute_mode(X)\n",
    "\n",
    "print(f\"Mean: {x_mean}\\nMedian: {x_median}\\nMode: {x_mode}\\n\" )\n",
    "\n",
    "y_mean = compute_mean(Y)\n",
    "y_median = compute_median(Y)\n",
    "y_mode = compute_mode(Y)\n",
    "\n",
    "print(f\"Mean: {y_mean}\\nMedian: {y_median}\\nMode: {y_mode}\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b)** Using the sampled data of a) find their:\\\n",
    "i) Probability Mass Function, ii) Expected Value and iii) Variance.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "For X:\n",
      "Expected Value of X is 25.005264000000004\n",
      "Variance of X is 25.050060290304\n",
      "\n",
      "For Y:\n",
      "Expected Value of Y is 60.002703999999994\n",
      "Variance of Y is 99.99015668838402\n"
     ]
    }
   ],
   "source": [
    "# In order to make our life a lot easier, we will compute the probability of each element\n",
    "# Because that will be useful for the next couple of tasks\n",
    "# NOTE: We can reuse the dictionary from earlier (Task A)\n",
    "\n",
    "\n",
    "def compute_pmf(X):\n",
    "    \"\"\"\n",
    "    Input: certain list of data   \n",
    "    Output: pmf as a dictionary \n",
    "    Example: Input: [0,1,0,2,3,2,2] => Output: {0:2/7, 1:1/7, 2:3/7, 3:1/7}\n",
    "    Hint: you can use count_dictionary function from a)\n",
    "    \"\"\"\n",
    "    # TODO: return pmf of X\n",
    "    count=count_dictionary(X)\n",
    "    for val in count:\n",
    "        count[val] = count[val]/len(X)\n",
    "    return count\n",
    "\n",
    "def compute_ev(X):\n",
    "    \"\"\"\n",
    "    Input: certain list of data   \n",
    "    Output: expected value of X \n",
    "    Hint: you can use compute_pmf function \n",
    "    μ = Σx * P(x) \n",
    "    \"\"\"\n",
    "    # TODO: return the expected value of X\n",
    "    pmf = compute_pmf(X)\n",
    "    expected_value = 0\n",
    "    for val, prob in pmf.items():\n",
    "        expected_value += val * prob\n",
    "    return expected_value\n",
    "\n",
    "def compute_variance(X):  \n",
    "    \"\"\"\n",
    "    Input: certain list of data   \n",
    "    Output: variance of X \n",
    "    Hint: you can use compute_pmf and compute_ev functions\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    # TODO: return the variance of X\n",
    "    pmf = compute_pmf(X)\n",
    "    ev = compute_ev(X)\n",
    "    variance = 0\n",
    "    for val, prob in pmf.items():\n",
    "        variance += prob * ((val - ev)**2)\n",
    "    return variance\n",
    "\n",
    "\n",
    "    \n",
    "print(\"\\nFor X:\")\n",
    "x_pmf = compute_pmf(X)\n",
    "x_ev = compute_ev(X)\n",
    "print(f'Expected Value of X is {x_ev}')\n",
    "x_variance = compute_variance(X)\n",
    "print(f'Variance of X is {x_variance}')\n",
    "\n",
    "print(\"\\nFor Y:\")\n",
    "y_pmf = compute_pmf(Y)\n",
    "y_ev = compute_ev(Y)\n",
    "print(f'Expected Value of Y is {y_ev}')\n",
    "y_variance = compute_variance(Y)\n",
    "print(f'Variance of Y is {y_variance}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c)** Suppose the only information we had is the sampled data $X$ and $Y$ for the *Age* and *Weight* respectively. Plot the distributions of $X$ and $Y$ and consider if the data seem to fit a normal distribution.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABnW0lEQVR4nO3deXxU5dnG8d+ZyR4gYZGEQCCgKCAIgoBBK1JT4S1tpVpEqqKItlJRFEXFBbFW4wYFBUUQFKsUxCoqIopRECSC7KDs+5KFLStkmznvH7NAJEAmJDkzyfXtZz5JJs/M3HOoyZX7ec5zDNM0TURERET8mM3qAkRERETORYFFRERE/J4Ci4iIiPg9BRYRERHxewosIiIi4vcUWERERMTvKbCIiIiI31NgEREREb8XZHUBlcHpdHLw4EHq1q2LYRhWlyMiIiLlYJomubm5xMXFYbOdvYdSIwLLwYMHiY+Pt7oMERERqYB9+/bRrFmzs46pEYGlbt26gOsN16tXz+JqREREpDxycnKIj4/3/h4/mxoRWDzTQPXq1VNgERERCTDlWc6hRbciIiLi9xRYRERExO8psIiIiIjfqxFrWERERM7ENE1KSkpwOBxWl1Ir2e12goKCznvbEQUWERGpsYqKikhLS+P48eNWl1KrRURE0KRJE0JCQir8HAosIiJSIzmdTnbt2oXdbicuLo6QkBBtLlrNTNOkqKiIQ4cOsWvXLlq3bn3ODeLORIFFRERqpKKiIpxOJ/Hx8URERFhdTq0VHh5OcHAwe/bsoaioiLCwsAo9jxbdiohIjVbRv+il8lTGv4H+FUVERMTvKbCIiIiI31NgERER8VOpqanY7Xb69u1rdSmWU2ARERHxU9OmTeP+++/n+++/5+DBg1aXYykFFhERET+Ul5fH7NmzGTp0KH379uXdd98t9f3PPvuM1q1bExYWRq9evZgxYwaGYZCVleUds3TpUn7zm98QHh5OfHw8DzzwAPn5+dX7RiqJAosEjkOH4NZbYdYsqysRkUBlmpCfX/030/S51A8//JA2bdpwySWXcNtttzF9+nRM9/Ps2rWLv/zlL/Tr149169bx97//nSeffLLU43fs2EGfPn246aabWL9+PbNnz2bp0qUMGzasUg5lddM+LBIYnE6+79aba3avIe+jT7hmmZPVr/3V6qpEJNAcPw516lT/6+blQWSkTw+ZNm0at912GwB9+vQhOzubxYsXc+211/LWW29xySWX8MorrwBwySWXsHHjRp5//nnv45OTk7n11lt58MEHAWjdujWvvfYaPXv25M0336zwfihWqVCHZdKkSSQkJBAWFkb37t1ZsWLFWcfPmTOHNm3aEBYWRocOHZg/f36p7+fl5TFs2DCaNWtGeHg47dq1Y/LkyRUpTWqql17imt1rAKhTdIJ7f/zI4oJERKrOli1bWLFiBQMHDgQgKCiIAQMGMG3aNO/3u3btWuox3bp1K/X1unXrePfdd6lTp4731rt3b+8OwIHG5w7L7NmzGTFiBJMnT6Z79+6MHz+e3r17s2XLFho3bnza+GXLljFw4ECSk5P5wx/+wMyZM+nXrx+rV6+mffv2AIwYMYJvv/2W999/n4SEBL7++mv+8Y9/EBcXx5/+9Kfzf5cS2JYuhaefBuB/7X/LTRu/5Y7V82D/fmjWzOLiRCSgRES4uh1WvK4Ppk2bRklJCXFxcd77TNMkNDSUiRMnlus58vLy+Pvf/84DDzxw2veaN2/uUz1+wfRRt27dzPvuu8/7tcPhMOPi4szk5OQyx998881m3759S93XvXt38+9//7v360svvdT85z//WWpM586dzSeffLJcNWVnZ5uAmZ2dXd63IYHi0CHTbNrUNMH836W9zBaPfm7+2OxS0wTTPOX/QyIiv3bixAnzl19+MU+cOGF1KT4pLi42Y2JizLFjx5obNmwodbvwwgvNN99803zsscfMDh06lHrcU089ZQLmsWPHTNM0zb/+9a/mddddZ8E7ON2Z/i18+f3t05RQUVERq1atIikpyXufzWYjKSmJ1NTUMh+TmppaajxA7969S43v0aMHn332GQcOHMA0Tb777ju2bt3K9ddfX+ZzFhYWkpOTU+omNZDTCXfeCQcOwMUX89T1/wDD4NVrbnd9f9o02LHD0hJFRCrbvHnzOHbsGEOGDKF9+/albjfddBPTpk3j73//O5s3b+axxx5j69atfPjhh96ziDwXeHzsscdYtmwZw4YNY+3atWzbto1PP/00YBfd+hRYDh8+jMPhICYmptT9MTExpKenl/mY9PT0c45//fXXadeuHc2aNSMkJIQ+ffowadIkrrnmmjKfMzk5maioKO8tPj7el7chgeLdd+GLLyA0FD78kOMh4QD8FN+eRS27QEkJjBljaYkiIpVt2rRpJCUlERUVddr3brrpJlauXElubi4fffQRH3/8MZdddhlvvvmm9yyh0NBQAC677DIWL17M1q1b+c1vfsPll1/O6NGjS00zBRK/OEvo9ddf58cff+Szzz6jRYsWfP/999x3333ExcWd1p0BGDVqFCNGjPB+nZOTo9BSwyQ8/gWT5k6nLzChy438+7/7S33/1Wtu59pdq+CDD+Cxx8C9HkpEJNB9/vnnZ/xet27dvKc2X3bZZaXWeT7//PM0a9as1Nk/Xbt25euvv666YquRT4GlUaNG2O12MjIySt2fkZFBbGxsmY+JjY096/gTJ07wxBNP8Mknn3i3Hr7ssstYu3Ytr776apmBJTQ01JsgpeZql7kTgBXxl572vY2xF8FNN8H//udakPvJJ9VdnoiIpd544w26du1Kw4YN+eGHH3jllVcCdrqnPHyaEgoJCaFLly6kpKR473M6naSkpJCYmFjmYxITE0uNB1i4cKF3fHFxMcXFxaddetput+N0On0pT2qQyMLjtDyWBsAvjVuWOSYp+jrXJ3Pn0uX+D6qrNBERv7Bt2zZuuOEG2rVrx3PPPcfDDz/MmBo8Te7zlNCIESO44447uOKKK+jWrRvjx48nPz+fwYMHAzBo0CCaNm1KcnIyAMOHD6dnz56MHTuWvn37MmvWLFauXMmUKVMAqFevHj179mTkyJGEh4fTokULFi9ezHvvvce4ceMq8a1KIGl7yLVHwMG6jTgWcfo8LsD2Rs3Z3KgFbQ7vofu+jdVZnoiI5f7973/z73//2+oyqo3PgWXAgAEcOnSI0aNHk56eTqdOnViwYIF3Ye3evXtLdUt69OjBzJkzeeqpp3jiiSdo3bo1c+fO9e7BAjBr1ixGjRrFrbfeytGjR2nRogXPP/889957byW8RQlE7TJc00Fn6q54pLa4jDaH95C4d311lCUiIhap0KLbYcOGnXGebNGiRafd179/f/r373/G54uNjeWdd96pSClSQ13qDiw/x1x41nE/Nu/A4FWfk7hHgUVEpCbTxQ/FL3kW3P7SuNVZx/0Y3wEnBhcd3Q9padVRmoiIWECBRfxPcTEXH94DwM8xZw8s2eF12eSZNiqjuyciIjWDAov4n02bCHWUkBMSwf6omHMOT23ewfXJd99VcWEiImIVBRbxP2vXArApphW4t5g+m9QWl7k+UWAREamxFFjE/6xZA8DP51i/4rEivj0Owwbbt7uu4CwiIqUsWrQIwzDIysoq92PGjBlDp06dqqwmXymwiP9xd1h+Ocf6FY/c0Eg2es4mUpdFRALc5MmTqVu3LiUlJd778vLyCA4O5tprry011hNEdpzjQrA9evQgLS2tzOsTnY9rr72WBx98sFKf80wUWMS/mKY3sJxrwe2ptI5FRGqKXr16kZeXx8qVK733LVmyhNjYWJYvX05BQYH3/u+++47mzZtz4YVn3wIiJCSE2NhY75WcA5ECi/iXPXsgK4siWxDbG5b/gpY/Ntc6FhGpGS655BKaNGlSal+zRYsWccMNN9CyZUt+/PHHUvf36tULp9NJcnIyLVu2JDw8nI4dO/LRRx+VGvfrKaGpU6cSHx9PREQEf/7znxk3bhzR0dGn1fOf//yHhIQEoqKiuOWWW8jNzQXgzjvvZPHixUyYMAHDMDAMg927d1f24fBSYBH/4u6ubGvUnGJ7cLkf9lOzdmC3w+7drpuISBlM0yS/KL/ab54rLJdXr169+O6UP8C+++47rr32Wnr27Om9/8SJEyxfvpxevXqRnJzMe++9x+TJk/n555956KGHuO2221i8eHGZz//DDz9w7733Mnz4cNauXcvvfvc7nn/++dPG7dixg7lz5zJv3jzmzZvH4sWLefHFFwGYMGECiYmJ3HPPPaSlpZGWlkZ8fPn/0PRVhXa6FakyngW3PkwHAeSHRkDXrvDjj64ui/vaViIipzpefJw6yXWq/XXzRuURGRJZ7vG9evXiwQcfpKSkhBMnTrBmzRp69uxJcXExkydPBiA1NZXCwkKuvfZa2rVrxzfffOO9sHCrVq1YunQpb731Fj179jzt+V9//XX+7//+j0ceeQSAiy++mGXLljFv3rxS45xOJ++++y5169YF4PbbbyclJYXnn3+eqKgoQkJCiIiIIDY2tkLHxRfqsIh/8Sy4LecZQqX06uX6qGkhEQlw1157Lfn5+fz0008sWbKEiy++mAsuuICePXt617EsWrSIVq1akZeXx/Hjx/nd735HnTp1vLf33nvvjItxt2zZQrdu3Urd9+uvARISErxhBaBJkyZkZmZW7pstJ3VYxL+4OyzlPUOolF69IDnZFVhMs1x7uIhI7RIRHEHeqDxLXtcXF110Ec2aNeO7777j2LFj3i5JXFwc8fHxLFu2jO+++47f/va35OW53s8XX3xB06ZNSz1PaGjoedUdHFx6at4wDJxO53k9Z0UpsIj/OHIE9u0DOLndvi+uugqCg117sezYARddVMkFikigMwzDp6kZK/Xq1YtFixZx7NgxRo4c6b3/mmuu4csvv2TFihUMHTqUdu3aERoayt69e8uc/inLJZdcwk8//VTqvl9/XR4hISE4HA6fH1cRmhIS/+GeDqJVK3JDff+BkvDP71jbyBV0hj4yrRILExGpfr169WLp0qWsXbu2VBDp2bMnb731FkVFRfTq1Yu6devyyCOP8NBDDzFjxgx27NjB6tWref3115kxY0aZz33//fczf/58xo0bx7Zt23jrrbf48ssvfT7tOSEhgeXLl7N7924OHz5cpd0XBRbxH57AcvnlFX6KnQ1c7dCELF25WUQCW69evThx4gQXXXQRMTEnr6vWs2dPcnNzvac/Azz33HM8/fTTJCcn07ZtW/r06cMXX3xBy5Zld6uvuuoqJk+ezLhx4+jYsSMLFizgoYceIiwszKcaH3nkEex2O+3ateOCCy5g7969FX/D52CYvp5r5YdycnKIiooiOzubevXqWV2OVNTtt8P778Nzz5GQV7HQMnzpTB76YSazLrueW9Z9VckFikggKSgoYNeuXbRs2dLnX8S10T333MPmzZtZsmRJpT/3mf4tfPn9rQ6L+IWEx79g01dLARiytqjCz7OnvuuvjYRjByulLhGRmurVV19l3bp1bN++3Tt9dMcdd1hd1hlp0a34BZvTQUt3yNjWsHmFn2dPtCuwtDimKSERkbNZsWIFL7/8Mrm5ubRq1YrXXnuNu+++2+qyzkiBRfxCXM4hwkqKKLQHsT+qcYWfZ7e7w9Ik7wgcPw4Rvp1KKCJSW3z44YdWl+ATTQmJX7joyH4AdtVvitNmr/DzHAuvR47nDKOdOyujNBER8QMKLOIXLjzqCiw7GjY7vycyDG+Xhe3bz7MqERHxFwos4hcudHdYdjQ4z8DCyXUsnGFLahGpXWrAybABrzL+DRRYxC9UWocF2F0/zvWJOiwitZpnW/njx49bXIl4/g1+vdW/L7ToVvxCK0+HpeH5X5p8j6aERASw2+1ER0d7L9YXERHh806ucn5M0+T48eNkZmYSHR2N3V7xNYoKLGK9o0e54HgWcHKn2vOhNSwi4hEbGwtg2RWGxSU6Otr7b1FRCixivS1bADhYtxHHQ8LP++n2RLunhPbuhcJCOM+rlYpI4DIMgyZNmtC4cWOKi4utLqdWCg4OPq/OiocCi1hv82agchbcAhyKjCY/OIzI4gLYvRsuuaRSnldEApfdbq+UX5piHS26Fet5AkslLLgFwDDYG+1uPWpaSESkRlBgEetVdmBBZwqJiNQ0CixivUqeEoJTzhTSXiwiIjWCAotYq7jYu4V+pXZYonWmkIhITaLAItbasQNKSsgLCSejTsNKe9o9mhISEalRFFjEWu7poJ0NmkIlbujk3Ytl1y4oKam05xUREWtUKLBMmjSJhIQEwsLC6N69OytWrDjr+Dlz5tCmTRvCwsLo0KED8+fPL/V9wzDKvL3yyisVKU8CSRWsXwFIr9vQtf9KSYlrPxYREQloPgeW2bNnM2LECJ555hlWr15Nx44d6d279xl3EVy2bBkDBw5kyJAhrFmzhn79+tGvXz82btzoHZOWllbqNn36dAzD4Kabbqr4O5PAUAVnCAGYhg0uvND1haaFREQCns+BZdy4cdxzzz0MHjyYdu3aMXnyZCIiIpg+fXqZ4ydMmECfPn0YOXIkbdu25bnnnqNz585MnDjROyY2NrbU7dNPP6VXr160atWq4u9MAkMVdVgAuOgi10cFFhGRgOdTYCkqKmLVqlUkJSWdfAKbjaSkJFJTU8t8TGpqaqnxAL179z7j+IyMDL744guGDBlyxjoKCwvJyckpdZMAZJpV1mEB1GEREalBfAoshw8fxuFwEBMTU+r+mJgY0tPTy3xMenq6T+NnzJhB3bp1ufHGG89YR3JyMlFRUd5bfPz5X+FXLJCRAdnZYLOdPKunMqnDIiJSY/jdWULTp0/n1ltvJSws7IxjRo0aRXZ2tve2b9++aqxQKo27u0LLlhQGhVT+83sCizaPExEJeD5d/LBRo0bY7XYyMjJK3Z+RkXHGy0bHxsaWe/ySJUvYsmULs2fPPmsdoaGhhOoKvIHPE1jatKma5z81sDidYPO7fC4iIuXk00/wkJAQunTpQkpKivc+p9NJSkoKiYmJZT4mMTGx1HiAhQsXljl+2rRpdOnShY4dO/pSlgSqqg4szZtDUBAUFsKBA1XzGiIiUi18/pNzxIgRTJ06lRkzZrBp0yaGDh1Kfn4+gwcPBmDQoEGMGjXKO3748OEsWLCAsWPHsnnzZsaMGcPKlSsZNmxYqefNyclhzpw53H333ef5liRgVHVgCQqCli1dn2sdi4hIQPNpSghgwIABHDp0iNGjR5Oenk6nTp1YsGCBd2Ht3r17sZ3Seu/RowczZ87kqaee4oknnqB169bMnTuX9u3bl3reWbNmYZomAwcOPM+3JAHj1MCyPbtqXuOii2DbNldg6dWral5DRESqnGGapml1EecrJyeHqKgosrOzqVevntXlSHkUFEB4uOvzzEwSxp59t+SKenbhm9yx+gveuPIvvNzzTna/2LdKXkdERHzny+9vrUIUS1wz/H0ATgSFkvDq8ip7nfS6jQCIyTtaZa8hIiJVT4FFLBGbdwSA9LoNKvWih7+WUacBAI1zFVhERAKZAotYIjb3MHCyA1JVMuo0BNRhEREJdAosYonYXHeHxR0oqoqnwxLj7uiIiEhgUmARS3gCS0bdqg0sme7AElWYT1hxQZW+loiIVB0FFrGEZ0oorYqnhHJCIzkR5NoVuXHesSp9LRERqToKLGIJ76LbKp4SwjBcC3vRtJCISCBTYBFLxLjP2qnqKSGATC28FREJeAosUv0cDm+3I60aAsvJhbcKLCIigUqBRapfZiZBphOHYeNwZP0qfznvXiwKLCIiAUuBRaqf+8rJhyKjcdjsVf5yOrVZRCTwKbBI9du/H6j6TeM8tIZFRCTwKbBI9XN3WNKrYf0KaEpIRKQmUGCR6ucJLFV9SrObFt2KiAQ+BRapftU+JeQKLHWKTkBubrW8poiIVC4FFql+1TwldDwknJyQiFKvLSIigUWBRaqfOzRkVNOUEJzssnDwYLW9poiIVB4FFqlepumdEqqOTeM8MuoqsIiIBDIFFqleOTmQnw9U35QQnLLAV4FFRCQgKbBI9XJPB2WHRlIQHFZtL5upwCIiEtAUWKR6ec8Qqr7uCpw8tVmBRUQkMCmwSPXyniFUPac0eyiwiIgENgUWqV7VvGmcR4amhEREApoCi1Qvi6aEMk89S8g0q/W1RUTk/CmwSPWyaEooM9IdWAoL4dixan1tERE5fwosUr2qeZdbj6KgYI6G13N9oWkhEZGAo8Ai1cs9JVSdu9x6aOGtiEjgUmCR6lNYCIcOAdW7y62HtucXEQlcCixSfdLSXB9DQznmmZ6pRuqwiIgELgUWqT7u6SDi4sAwqv3ltT2/iEjgUmCR6uNecEuzZpa8fGZdBRYRkUClwCLVxxNYmja15OU1JSQiErgqFFgmTZpEQkICYWFhdO/enRUrVpx1/Jw5c2jTpg1hYWF06NCB+fPnnzZm06ZN/OlPfyIqKorIyEi6du3K3r17K1Ke+CvPlJDVgcUTnEREJGD4HFhmz57NiBEjeOaZZ1i9ejUdO3akd+/eZGZmljl+2bJlDBw4kCFDhrBmzRr69etHv3792Lhxo3fMjh07uPrqq2nTpg2LFi1i/fr1PP3004SFVd/VfKUaWDwl5A0saWngdFpSg4iIVIxhmr7tU969e3e6du3KxIkTAXA6ncTHx3P//ffz+OOPnzZ+wIAB5OfnM2/ePO99V155JZ06dWLy5MkA3HLLLQQHB/Of//ynQm8iJyeHqKgosrOzqVev+s8+kXK6+mr44Qf48EMSVkVU+8vbnQ52vNrPtTV/ejrExFR7DSIicpIvv7996rAUFRWxatUqkpKSTj6BzUZSUhKpqallPiY1NbXUeIDevXt7xzudTr744gsuvvhievfuTePGjenevTtz5849Yx2FhYXk5OSUukkAsHhKyGGznwwpWsciIhJQfAoshw8fxuFwEPOrv0xjYmJIT08v8zHp6elnHZ+ZmUleXh4vvvgiffr04euvv+bPf/4zN954I4sXLy7zOZOTk4mKivLe4uPjfXkbYgWn82RIsGhKCHCdUg0KLCIiAcbys4Sc7rUEN9xwAw899BCdOnXi8ccf5w9/+IN3yujXRo0aRXZ2tve2b9++6ixZKuLwYSgudu2/0qSJdXUosIiIBKQgXwY3atQIu91ORkZGqfszMjKIjY0t8zGxsbFnHd+oUSOCgoJo165dqTFt27Zl6dKlZT5naGgooaGhvpQuVvNMBzVuDMHB1tWhwCIiEpB86rCEhITQpUsXUlJSvPc5nU5SUlJITEws8zGJiYmlxgMsXLjQOz4kJISuXbuyZcuWUmO2bt1KixYtfClP/NhdL3wKwAbqkPD4F9YVosAiIhKQfOqwAIwYMYI77riDK664gm7dujF+/Hjy8/MZPHgwAIMGDaJp06YkJycDMHz4cHr27MnYsWPp27cvs2bNYuXKlUyZMsX7nCNHjmTAgAFcc8019OrViwULFvD555+zaNGiynmXYrmYvKMApFtw0cNSPIFFe7GIiAQUn9ewDBgwgFdffZXRo0fTqVMn1q5dy4IFC7wLa/fu3Uua5yJ3QI8ePZg5cyZTpkyhY8eOfPTRR8ydO5f27dt7x/z5z39m8uTJvPzyy3To0IG3336b//3vf1x99dWV8BbFH8TkHQHgUGQDS+u482tXUNm4crO1nR4REfGJzx0WgGHDhjFs2LAyv1dWV6R///7079//rM951113cdddd1WkHAkAjd0dFu/mbRbJdF8A0ROgREQkMFh+lpDUDjF+Elg8U1IN87MJcpRYWouIiJSfAotUC2+HxeI1LMfC61JkC8KGSeP8o5bWIiIi5afAItXC02HJtLjDYho2bw0xuQosIiKBQoFFql5JCY3yswDItHjRLZycltI6FhGRwKHAIlUvIwMbJiWGjSMR1l+c0rOOxdP1ERER/6fAIlXPvUnbocj6OG12i4uBDPeZQrG56rCIiAQKBRapeu59eTLqWj8dBCfr0JSQiEjgUGCRqufusHj2QLFauvZiEREJOAosUvU8gSWyvsWFuGR6p4S0hkVEJFAosEjVcwcWqzeN8/Asum2sDouISMBQYJGq5w0s/jEl5AlOdYtOQG6uxdWIiEh5KLBI1XMvurV60ziP4yHh5IREuL7QVZtFRAKCAotUPU+HxU/OEoJTLhHgrk1ERPybAotUreJiyMwE/GOXWw/vehp1WEREAoICi1St9HQAim12jvrBLrce3g6LAouISEBQYJGq5T2luQGm4T//d/N2WDQlJCISEPznN4jUTH624NbDs3mcOiwiIoFBgUWqlh8uuAUtuhURCTQKLFK1TpkS8icZ6rCIiAQUBRapWn62y62Hd0ooLQ2cTmuLERGRc1JgkarlvfChfwWWw5HRODGgpAQOHbK6HBEROQcFFqla7kW3/tZhKbEHcTgy2vWFpoVERPyeAotULe+iW/+4jtCp0rXwVkQkYCiwSNUpLITDhwHIjKxvcTGn0263IiKBQ4FFqo57l1uCgzkW7j+73HroTCERkcChwCJVxzPVEhcHhmFtLWXQbrciIoFDgUWqzqmBxQ+l63pCIiIBQ4FFqo77DCF/DSyZmhISEQkYCixSdQKlw6IpIRERv6fAIlXHEwSaNLG2jjPw7nZ75AgUFFhbjIiInJUCi1QdP++wZIfVgdBQ1xee6SsREfFLCixSdfw8sGAY0LSp63OtYxER8WsVCiyTJk0iISGBsLAwunfvzooVK846fs6cObRp04awsDA6dOjA/PnzS33/zjvvxDCMUrc+ffpUpDTxJ36+6BY4WZsCi4iIX/M5sMyePZsRI0bwzDPPsHr1ajp27Ejv3r3JzMwsc/yyZcsYOHAgQ4YMYc2aNfTr149+/fqxcePGUuP69OlDWlqa9/bf//63Yu9I/ENBARw96vrcnwOLp8OihbciIn7N58Aybtw47rnnHgYPHky7du2YPHkyERERTJ8+vczxEyZMoE+fPowcOZK2bdvy3HPP0blzZyZOnFhqXGhoKLGxsd5b/fr+t5W7+MDTXQkNhehoS0s5K00JiYgEBJ8CS1FREatWrSIpKenkE9hsJCUlkZqaWuZjUlNTS40H6N2792njFy1aROPGjbnkkksYOnQoR44cOWMdhYWF5OTklLqJn/HzXW69PN0fdVhERPyaT4Hl8OHDOBwOYmJiSt0fExNDuue6Mb+Snp5+zvF9+vThvffeIyUlhZdeeonFixfzf//3fzgcjjKfMzk5maioKO8tPj7el7ch1cHfF9x6qMMiIhIQgqwuAOCWW27xft6hQwcuu+wyLrzwQhYtWsR111132vhRo0YxYsQI79c5OTkKLf4mEBbcghbdiogECJ86LI0aNcJut5ORkVHq/oyMDGJjY8t8TGxsrE/jAVq1akWjRo3Yvn17md8PDQ2lXr16pW7iPxIe/4I3Zy4B4J2dBSQ8/oXFFZ3FqR0W07S2FhEROSOfAktISAhdunQhJSXFe5/T6SQlJYXExMQyH5OYmFhqPMDChQvPOB5g//79HDlyhCZ+ukOqnFvjfNcZQpmeKyL7qUveWOf6pKCAyx760NpiRETkjHw+S2jEiBFMnTqVGTNmsGnTJoYOHUp+fj6DBw8GYNCgQYwaNco7fvjw4SxYsICxY8eyefNmxowZw8qVKxk2bBgAeXl5jBw5kh9//JHdu3eTkpLCDTfcwEUXXUTv3r0r6W1Kdbsg7xgAmZH+HVgKg0M5FlYXgCa5hy2uRkREzsTnNSwDBgzg0KFDjB49mvT0dDp16sSCBQu8C2v37t2LzXYyB/Xo0YOZM2fy1FNP8cQTT9C6dWvmzp1L+/btAbDb7axfv54ZM2aQlZVFXFwc119/Pc899xyhnm3TJeA0zvN0WPz/9PT0ug2pX5BLbO6Zz0wTERFrVWjR7bBhw7wdkl9btGjRaff179+f/v37lzk+PDycr776qiJliB9rnO/usPj5lBC4AkvbQ7uJVYdFRMRv6VpCUumCHcU0OOHaGycz0v87LGl1GwGowyIi4scUWKTSXeDurhTZgjgW7v9ncGXUaQigDouIiB9TYJFK19i94PZQZH3/3uXWzdthyVOHRUTEXymwSKXzLLg9FAALbgEy6no6LAosIiL+SoFFKl0gLbgF16JbUGAREfFnCixS6S7wnNIcAAtu4eSUUP2CXDhxwuJqRESkLAosUulO7sESGB2WnNBITgS59/zRNYVERPySAotUOu+UUIB0WDAM0tzTQgosIiL+SYFFKl2gdVjg5MJbBRYREf+kwCKVLtAW3cLJdSzs329tISIiUiYFFqlcDgeN8rOAAJoS4uTmceqwiIj4JwUWqVyHDmE3nTgxOBIZbXU15eY5tVkdFhER/6TAIpUrLQ2AI5FROGx2i4spv3StYRER8WsKLFK53IElMzJw1q8ApGtKSETErymwSOXyBJYA2Zbfw7voNi0NHA5rixERkdMosEjlCtAOy+HIaEoMmyusZGRYXY6IiPyKAotULm+HJbACi9NmP1mzFt6KiPgdBRapXAE6JQQ6tVlExJ8psEjlCtApIUDb84uI+DEFFqlc7sByKBA7LNqLRUTEbymwSOUxTUhPBwJvDQtoLxYREX+mwCKVJysLCgsBOBRA2/J76HpCIiL+S4FFKo97OigrrA6FQSEWF+M7LboVEfFfCixSeQJ4wS2c0mE5cMA1vSUiIn5DgUUqTwCf0gyQ4Vl3c/y4a3pLRET8hgKLVJ4A3TTOozA4FBpqWkhExB8psEjl8U4JBWaHBYCmTV0ftfBWRMSvKLBI5fHuwRKYHRbgZGBRh0VExK8osEjlqQkdlmbNXB/VYRER8SsKLFJ5AnwNC6AOi4iIn1JgkcpTEwKLp8OiwCIi4lcUWKRy5OdDbi4Q4FNCWnQrIuKXFFikcri7K0REkB8Sbm0t50MdFhERv1ShwDJp0iQSEhIICwuje/furFix4qzj58yZQ5s2bQgLC6NDhw7Mnz//jGPvvfdeDMNg/PjxFSlNrOIJLE2agGFYW8v58HRYjhyBEyesrUVERLx8DiyzZ89mxIgRPPPMM6xevZqOHTvSu3dvMjMzyxy/bNkyBg4cyJAhQ1izZg39+vWjX79+bNy48bSxn3zyCT/++CNxcXG+vxOx1qmBJZBFR0O4u0N08KClpYg1CksKmbF2Bvuy91ldioicwufAMm7cOO655x4GDx5Mu3btmDx5MhEREUyfPr3M8RMmTKBPnz6MHDmStm3b8txzz9G5c2cmTpxYatyBAwe4//77+eCDDwgODq7YuxHr1JTAYhg6tbkWM02TIZ8N4c5P7yRxWiIHcjQ1KOIvfAosRUVFrFq1iqSkpJNPYLORlJREampqmY9JTU0tNR6gd+/epcY7nU5uv/12Ro4cyaWXXnrOOgoLC8nJySl1E4vVkMCS8PgXLC8IBeD+sV9YXI1Ut4nPRvDBhg8AOJB7gL7jmpNbmGtxVSICPgaWw4cP43A4iImJKXV/TEwM6enpZT4mPT39nONfeuklgoKCeOCBB8pVR3JyMlFRUd5bfHy8L29DqkINCSxw8iKIjfOPWVyJVIW16WtZvn/5afcv3buUERQAMNIMobFpsM5wcvNHN1PsKK7uMkXkVyw/S2jVqlVMmDCBd999F6OcizVHjRpFdna297Zvn+aaLVeDAovntOwL8o5aXIlUtq1HttJtajeunHYlv//g9/yc+TMAablp9J/TnxIDbjGDeIlQ5hFBuAkLti/gvvn3YZqmxdWL1G5Bvgxu1KgRdrudjIyMUvdnZGQQGxtb5mNiY2PPOn7JkiVkZmbSvHlz7/cdDgcPP/ww48ePZ/fu3ac9Z2hoKKGhob6ULlUo4fEv+Hr1Fi4Gbv1yPyQ0srqk85JRx3XF5hgFlhpn1OuXUmyUAPDl9i/5asdX3H353Ww8tJH0vHTamzbeJhwDg67Y+S/h/JkCpq6eis2wcWWzK4mJjKFxZGNaRLegUURg/39dJJD41GEJCQmhS5cupKSkeO9zOp2kpKSQmJhY5mMSExNLjQdYuHChd/ztt9/O+vXrWbt2rfcWFxfHyJEj+eqrr3x9P2IF0yQu9xAAaXUD/wd4Zh1Xh6VxvgJLTbJs3zI+NkqwmTDXDOdGMwin6WTK6iks27eMqNAoPiacSE52em8gmAl9JgDw1qq3GPzpYH4/8/dcMfUKGr/SmDvn3qmziUSqiU8dFoARI0Zwxx13cMUVV9CtWzfGjx9Pfn4+gwcPBmDQoEE0bdqU5ORkAIYPH07Pnj0ZO3Ysffv2ZdasWaxcuZIpU6YA0LBhQxo2bFjqNYKDg4mNjeWSSy453/cn1aBeYT51ilx7lhysVxMCi3sNS57WsNQUpmnyyNePAHAXwdzgvi0xS3iEAtbj5IOCIlpz+hmK93/5FPUJYyEOMnCSgUkmJgcNkxnrZjBr4yyGdx/OqN+MIjosuprfmUjt4XNgGTBgAIcOHWL06NGkp6fTqVMnFixY4F1Yu3fvXmy2k42bHj16MHPmTJ566imeeOIJWrduzdy5c2nfvn3lvQuxVJPcwwAcDa9HQXCYxdWcP09g0ZRQzfHJ5k9I3Z9KuAnPcnI6+TcEsZw6FGESwpnX0N1GCLf96r4VpoORFPC9o5CXl73M22veZka/Gfzh4j9U0bsQqd0MswasJMvJySEqKors7Gzq1atndTm1zp39x/DuR8/yc+NW9B38mtXlnLe6hflsGD/A9cXx4yc3kpOAVOwo5tI3LmXb0W08ZYbwHJUXqk1MvqCERylkk+GkoWmwjTrUH6OtFkTKw5ff35afJSSBL87dYakJ00EAuSERnAhy/xXuOftJAtbU1VPZdnQbF0RcwKNU7mJ9A4M/EMw6IrnUtHHEMHmGwkp9DRFxUWCR89YkxxVY0upeYHEllcQwvHuxaHv+wJZVkMWYRWMAGHPtGOqeZdrnfARjMMHduXmDIjZmnn7pERE5Pwosct7iclzXkUqrIR0WOHmmkDosgeuHvT/Q+a3OHDp+iNYNWnNP53uq9PWuI4h+ZhAOAx5c8KD2bRGpZAosct68U0I14JRmj8xIdVgCVZGjiCdSnuCa6VezK2sXLUyD/x45SPBzVf//z7GEEWpCyq4U5m6eW+WvJ1KbKLDIefNOCdWrIVNCnDxTSB2WwLLtyDaufPtKkpcm4zTgDjOY9dShC/Zqef1W2HiYEAAe/vphCkoKquV1RWoDBRY5P6bpPa35oAKLWKjIUcQNs25gTfoaGoY35CMznHcJp14VrVs5k1GEElc3jl1ZuxiXOq5aX1ukJlNgkfNz6BChjmKcGKTXaXju8QHCu4ZFU0IBY8Jz9dh0eBMXmAbrjxdyUxmbwFWHOhi8nPQyAC8seYHM/ExL6hCpaRRY5Py4Lzx5qE59Suw+70PotzzXE1KHJTDsz9nPs+7TiV8mlDiLf7T99eOhXGHayC/O5+VXdDV5kcqgwCLnxx1YasI1hE7luWKzAktgeOTrR8g3oIdpZ5BFnZVTGRj8073nyxsUkZ6XbnFFIoFPgUXOz969QM06Qwggo667w3L0KBRo4aQ/+3bXt8z+eTY2EyYRhq2a16ycSR+CuNK0c8KAl5a+ZHU5IgFPgUXOj6fDUoMW3ALkhEZSaHf/pZ6uv479VZGjiGHzhwEwlGA6VdPZQOVxapflzZVvcjBX66FEzocCi5wfd2CpaR0WDOPkmUJaeOu3Xlv+mmuhbcQFlXqNoMqShJ2rTDuFjkKSlyRbXY5IQFNgkfPjCSw1rMMCnNyeX+tY/NKR40d4dvGzALyU9BL1/WQq6FSndlmmrJ7Cvux9FlckErgUWOT81NApIdDCW3/32stx5BXl0cm0ccenw60u54x6Yadni54UOYpIXqoui0hFKbBIxZWUeKdLatyUEGhKyI/lFubyGkUAPEmo3yy0LYuBwbPXujpBb69+m91Zu60tSCRAKbBIxaWlgcNBsc3O4choq6updNrt1n9NXjmZLAPamDZuxP/3/+n57p+4zrRT7Czm7vEX4RxTz+qSRAKOAotUnHs6KKNOQ5w2/zk7o7IosPingpICxqaOBeBxQvy6u3KqSYQRbkKK4eANiq0uRyTgKLBIxXkX3Na86SA4ZQ2LpoT8yvQ108nIz6C5afBXP9gkrrwuwc7L7jOZHqWArUe2WlyRSGBRYJGK8+5yW/MW3MIpm8epw+I3ih3FvPyD6zo9jxJKcIB0Vzz+QTBJ7s3kBn0yiBJnidUliQQMBRapuBp8SjOc0mE5fBiKiqwtRgD478b/sid7D40jG3NXAHVXPGwYTCecKBOWH1iuHXBFfKDAIhVXw6eEjoXXg2DtdusvnKbTe1rwiCtHEB5g3RWPeGxMdE8NjVk8hjVpayyuSCQwKLBIxbmvI1RTp4QwDGjSxPW5poUs98mmT9h8eDNRoVEM7TrU6nLOy60Ec2PbGylxlnDvF/dimqbVJYn4PQUWqTjvpnE1s8MCnAwsWnhrKeeYejz74c0APFBwgnrJ8RZXdH4MDCb98jURJqw4sILPno20uiQRv6fAIhVTWAiZmUDN3DTOKy7O9VEdFkt9TAkbDCf1THjIvdV9oIvFxgOEAPA0hThNp8UVifg3BRapmP37XR/DwlxrPWoqTQlZzmk6eZZCAB4kxC+vGVRRIwklyoQNhpPZG2dbXY6IX1NgkYpxTwcRH+9a61FTaUrIcv/75X9sNJxE1aDuikcDDB5xv6fRi0ZT7NCGciJnosAiFeMJLM2bW1tHVdOUkKWcptN7ReaHCCW6BnVXPIYTQiPTYPvR7cxYN8PqckT8lgKLVIz7DCHiA3vx4zlpSshSc36ew8+HfibKdP1ir4nqYvCE+739c/E/KSwptLgiEf+kwCIVc+qUUE2mKSHLOJwO/vn9PwEYUUO7Kx5DCaFp3absy9nHW6vesrocEb+kwCIVU1sCi2dK6NAhKNE26tVpzi9z+OXQL0SHRdfY7opHGAZPX/M0AM8veZ4TxScsrkjE/yiwSMXUlsDSqBEEBYFpQkaG1dXUGqZp8sKSFwB4OPFhompwd8XjrnkjaWEaZOZn8p/nG1hdjojfUWCRiqktgcVmg5gY1+eaFqo23z1bhw2ZG4g04b5vX7a6nGoRjMGD7k7SOIq0L4vIr1QosEyaNImEhATCwsLo3r07K1asOOv4OXPm0KZNG8LCwujQoQPz588v9f0xY8bQpk0bIiMjqV+/PklJSSxfvrwipUk1uPShOZCVBUD7t3+xtpjqoDOFqt14XBebvJPgGrXvyrkMIYQoE7YYTuZvm3/uB4jUIj4HltmzZzNixAieeeYZVq9eTceOHenduzeZ7l1Pf23ZsmUMHDiQIUOGsGbNGvr160e/fv3YuHGjd8zFF1/MxIkT2bBhA0uXLiUhIYHrr7+eQ4cOVfydSZWJy3H9u2SHRpIXGmFxNdVAC2+r1bYj25iHa73QAzV87cqv1cXgb+73PDZ1rMXViPgXnwPLuHHjuOeeexg8eDDt2rVj8uTJREREMH369DLHT5gwgT59+jBy5Ejatm3Lc889R+fOnZk4caJ3zF//+leSkpJo1aoVl156KePGjSMnJ4f169dX/J1JlWma4wqnB6IaW1xJ1Ut4/Ave3+vazOv19xaR8PgXFldU8722/DVMA/qaQVyM3epyqt39hBBkwqLdi1idttrqckT8hk+BpaioiFWrVpGUlHTyCWw2kpKSSE1NLfMxqamppcYD9O7d+4zji4qKmDJlClFRUXTs2LHMMYWFheTk5JS6SfVp6u6wHKhXQ6/S/Cv7o1xrWDxBTapOVkEW76x9B4CHall3xSMeGzcTBKjLInIqnwLL4cOHcTgcxHgWIbrFxMSQnp5e5mPS09PLNX7evHnUqVOHsLAw/v3vf7Nw4UIaNSr7onrJyclERUV5b/E1feGnn/F2WOrV/A4LwD53YInP0llCVW3a6mnkF+fT3rTx21rYXfF42L1d/+yNs9mXvc/iakT8g9+cJdSrVy/Wrl3LsmXL6NOnDzfffPMZ18WMGjWK7Oxs723fPv0HXZ3ivB2W2hFY9runvpplK7BUpRJnCa+veB1wXeTQqEWLbX+tM3auTbgWh+ngteWvWV2OiF/wKbA0atQIu91Oxq/2o8jIyCA2NrbMx8TGxpZrfGRkJBdddBFXXnkl06ZNIygoiGnTppX5nKGhodSrV6/UTapP02xXYDlYS6aE9kW7/r8ak3eUkBJdnK6qzN08lz3Ze2gU0Yi/Emx1OZZ7OPFhAKasnkJOoaa9RXwKLCEhIXTp0oWUlBTvfU6nk5SUFBITE8t8TGJiYqnxAAsXLjzj+FOft7BQ19TwR7Vp0S3A0fB6HA8OxYZJnNaxVJnxP44H4N4u9xJei7srHr9v/XvaNGpDTmEO/1n3H6vLEbGcz1NCI0aMYOrUqcyYMYNNmzYxdOhQ8vPzGTx4MACDBg1i1KhR3vHDhw9nwYIFjB07ls2bNzNmzBhWrlzJsGHDAMjPz+eJJ57gxx9/ZM+ePaxatYq77rqLAwcO0L9//0p6m1JpiouJyTsK1J5FtxjGyXUsmhaqEuvG1OGHfT8QZMLQ78dbXY5fsD1bn3sP7QLgnfkPwJgoiysSsVaQrw8YMGAAhw4dYvTo0aSnp9OpUycWLFjgXVi7d+9ebLaTOahHjx7MnDmTp556iieeeILWrVszd+5c2rdvD4Ddbmfz5s3MmDGDw4cP07BhQ7p27cqSJUu49NJLK+ltSqU5cAC76aTQHsThyGirq6k2+6NiuOTwXpplq8NSFabhmmq7gSDi/GdpneVuJZiRZiGrDCcbTAcdrC5IxEI+BxaAYcOGeTskv7Zo0aLT7uvfv/8ZuyVhYWF8/PHHFSlDrLB3LwBpdS/ANGrPLxbPwtv47LLPhpOKKygp4H33zrZ319JTmc+kETb+QBCfUMK7FKOTnKU2qz2/caRy7NkDwIGoWjId5OaZElKHpfLN3TyXYwbEmwa/q8WnMp/JYPcC5PcpptihRd9SeymwiG/cHZaDdWvHgluP/d7AojUslW3aGtfZgHcSjF2LbU/ThyAamwaZhsmX27+0uhwRyyiwiG9qeYdFi24r165ju/hm5zcYJtyl6aAyBWNwu7vL4tkFWKQ2UmAR37g7LLXmDCE3T4flgvwsOHHC2mJqEM8v4Ouwk6AfR2d0pzuwzNs6j0P5uiis1E76CSG+8XRYaskutx7ZYXXICXFfmXr3bktrqSkcToc3sAxRd+Ws2mPnCtNGibOEDzZ8YHU5IpZQYJHyM82Ta1hqWYcFwzi5UZ4CS6X4esfX7M/ZT4PwBvSr2AmLtcpgd6h7Z+07mKZpcTUi1U+BRcrvyBE4fhyAtNoWWDg5LcSuXdYWUkN4Ftve1uE2wrTY9pwGEkyIPYT1GetZm77W6nJEqp0Ci5Sfu7tyKDKawqDa18L3LLxVh+X8Hco/xGdbPgNgSOchFlcTGOpj0K9NP0CLb6V2UmCR8qul61c81GGpPNPWTKPYWUzXuK5cFnOZ1eUEjMEb5wHw3+UTKRqji75K7aKJYym/WnqGkMe+aHVYKkPRmHq8Th4YcN+BDbpGjg+SsBNrGqQbJgvMEv5kdUEi1UgdFim/Wt9hcb9vdVjOyyyKOWiYNDENBrpP15XyCcLgr+5j9h7a9VZqFwUWKT9PhyWqtgYWd4flyBHIzbW2mABlmiZj3dcNeoAQQrTY1meD3IHlc0o4duKYxdWIVB8FFik/d4el1p3S7JYbGklWWB3XF+5jIb75Zuc3rDecRJrwd+29UiEdsXOZaaPIgA9//tDqckSqjQKLlJ93DUvt7LCAFt6er7GprusNDyGE+uquVJiny/Le+vcsrkSk+iiwSPmcOAGZrisV19ZFt6BTm8/HhowNfLXjK2wmPKjuynn5K8HYTFi2bxnbj263uhyRaqHAIuXj7q5Qpw7ZnmmRWkgLbytu3I/jALiJIFrqR895aYKN32EH4P3171tcjUj10E8NKR9PYGneHIza28pXh6ViDuYe5IP1rmvgPKzuSqUY5D6O7617T1v1S62gwCLl41lk2qKFtXVYTGtYKmbiiokUO4u5Kv4qumv7p0rRjyDqhNRhV9Yuftj3g9XliFQ5BRYpn1M7LLWYOiy+KxhTjylLXgTg4b2rLa6m5ojA4C/t/gLAf9b9x+JqRKqeAouUjzoswCl70GRluW5yTh9RzBHDJN40+KO6K5Vq0GWDAJj982zyi/ItrkakaimwSPmowwLA8ZBwaNTI9YW6LOXypntH1r8RQpBOZa5UPWfcwIWmQXZhNtNfaGR1OSJVSoFFykcdlpNatnR91DqWc1qXvo5lhoMgE+7WNvyVzobBI4QCMJZCSpwlFlckUnUUWOTcHA7Yv9/1eS3vsACQkOD6qA7LOb258k0AbiSIWP24qRJ3EMwFpsEew2TOz3OsLkekyugniJxbejoUF4PdDnFxVldjPU+HRYHlrHIKc7x7hAzVqcxVJhyDB9zH9+VlL+sUZ6mxFFjk3DzTQU2bQpAWTXo7LDt3WlqGv/vPuv+QX5xPW9NGT/cmZ1I1/kEIESasTV/LNzu/sbockSqhwCLn9GDyxwCkEkXC419YXI31blt0BIAdy9bqeJyBaZre6aChhGBosW2VaoDBPad0WURqIgUWOaeWRw8CsKu+poMAdjRsBkDzrDSCHFrkWJale5fy86GfiQiO8F6oT6rWQ4RgN+x8s/MbVqdpvxupeRRY5JxaZLkCy576TSyuxD+k121IfnAYwU4HzbPSrS7HL72x8g0A/tr+r0Spu1ItWmDjlva3APDKslcsrkak8imwyDklHHMFlt3R6rAAmIaNXQ2aAnDh0f0WV+N/MvIy+N8v/wNgaNehFldTu4zsMRKAD3/+kJ3HtMZKahYFFjmnhGNpAOxuoMDisaOBa1rowiMKLL/21qstKHYWc6Vpp/Nbvawup1bpOPkaept2nKaTlye0gTFRVpckUmkUWOTsjh4luiAPgL1RsRYX4z8861haqcNSSpGjiMkUAXhPtZXq9aR7I7npFLMPp8XViFQeBRY5u23bAEiv04ATIWEWF+M/dnqmhNRhKeXjTR+TZpjEmgY36bpBlvgNQfQ07RQb8DKFVpcjUmkqFFgmTZpEQkICYWFhdO/enRUrVpx1/Jw5c2jTpg1hYWF06NCB+fPne79XXFzMY489RocOHYiMjCQuLo5BgwZx8ODBipQmlW37dgB26wyhUjwdlguP7gdt1OX12vLXALiXEEK02NYyo91dlqkUk5abZnE1IpXD58Aye/ZsRowYwTPPPMPq1avp2LEjvXv3JjMzs8zxy5YtY+DAgQwZMoQ1a9bQr18/+vXrx8aNGwE4fvw4q1ev5umnn2b16tV8/PHHbNmyhT/96U/n986kciiwlGlX/TicGK7pssOHrS7HL6w6uIrU/akEm/B3ncpsqV7Y6WHaKTR0xpDUHIbp4z7O3bt3p2vXrkycOBEAp9NJfHw8999/P48//vhp4wcMGEB+fj7z5s3z3nfllVfSqVMnJk+eXOZr/PTTT3Tr1o09e/bQvBzXrsnJySEqKors7Gzq1avny9uRc7ntNvjgA17qeQdvXtnf6mr8ytI376JZTiZ8/z385jdWl2O5O+feyYx1M7jVDOZ9wq0up9ZbQAn/ZxwnPCic3Q/upnFkY6tLEjmNL7+/feqwFBUVsWrVKpKSkk4+gc1GUlISqampZT4mNTW11HiA3r17n3E8QHZ2NoZhEB0d7Ut5UhXcHRZtGnc6z7QQmzdbW4gfyMzP5L8b/wvA/equ+IXe2Olq2jhRcoJxqeOsLkfkvPkUWA4fPozD4SAmJqbU/TExMaSnl72BVnp6uk/jCwoKeOyxxxg4cOAZ01ZhYSE5OTmlblJF3IFFm8adznNqM1u2WFuIH5i6aipFjiK6Ne1Gdy229QsGBk+717JMXDGRI8ePWFyRyPnxq7OEiouLufnmm13XIXnzzTOOS05OJioqynuLj4+vxiprkWPH4Ijrh9yeaAWWX9upDgsAxY5i73WD7u92v8XVyKn+QBCdYjuRX5zP+B/HW12OyHnxKbA0atQIu91ORkZGqfszMjKIjS17j47Y2NhyjfeElT179rBw4cKzzmWNGjWK7Oxs723fvn2+vA0pL3d3JaNOA46HaE3Cr6nD4jJ381wO5B6gcWRj+rfTOid/YmDwZJrr/59vff88RWO0xk8Cl0+BJSQkhC5dupCSkuK9z+l0kpKSQmJiYpmPSUxMLDUeYOHChaXGe8LKtm3b+Oabb2jYsOFZ6wgNDaVevXqlblIFPGcIqbtSpu2eDsvOnVBYS/e7GBPFa3MGAvD3vCxC/6WFnf7mBoJoYhocMkw+QxfrlMDl85TQiBEjmDp1KjNmzGDTpk0MHTqU/Px8Bg8eDMCgQYMYNWqUd/zw4cNZsGABY8eOZfPmzYwZM4aVK1cybNgwwBVW/vKXv7By5Uo++OADHA4H6enppKenU1RUVElvUyrEvWmc1q+U7VBkfXJDwsHp9Ia72mYNDpYaDoJMGKqdbf1SMAZ3uRdCT0E/UyVw+RxYBgwYwKuvvsro0aPp1KkTa9euZcGCBd6FtXv37iUt7eRGRT169GDmzJlMmTKFjh078tFHHzF37lzat28PwIEDB/jss8/Yv38/nTp1okmTJt7bsmXLKultSoVoD5azM4yTZwrV0mmh192/APsTRBP/WhInpxhCCIYJCw2HLoooAatCy/mHDRvm7ZD82qJFi067r3///vTvX/bcdkJCAj5uBSPVRYHlnHY0aEantG21cuHtofxDzKQY0HWD/F1LbFyPna9w8Pbqt3nhuhesLknEZ/qTSM5MU0LntLNB7T1T6O3Vb1NowBWmje7YrS5HzuFv7lA5fc10ih3FFlcj4jsFFilbVpZ3y3ktuj2z2jolVOwo5o2VbwCu7oqh6wb5vT8SRIxpkJGfwedbP7e6HBGfKbBI2XbscH2MiSE/NMLaWvzYjlM7LLVoanPu5rnsz9lPY9PgZu1sGxBKLb5dNcXiakR8p8AiZXNPB3HRRdbW4ef21I8Dmw1ycuAMuzfXRK+veB1wXeQwVN2VgHG3e1ro6x1fsztrt7XFiPhIgUXK5jlNt3Vra+vwc0VBwdCypeuLWjIttCZtDUv2LiHIFsS9WmwbUFphI6lVEiYm01ZPs7ocEZ8osEjZPIFFHZZza9PG9bGWLLx9bcVrAPyl3V+I04+QgPO3zn8DYNqaaRQ5tC+LBA79tJGyeaaE1GE5t0sucX2sBR2WXWPq8v6adwEYvuEza4uRCrlhzhCamAZpeWnMeC7K6nJEyk2BRcqmDkv51aIOywsUUmLA9aadK3VV5oAUgsGj7qm8FyjUKc4SMBRY5HQ5OZCZ6fpcgeXcaklg2XVsF++6N4p7hlCLq5Hz8TdCiDENdhsm/1n/H6vLESkXBRY5nae70rgx6MKS53T5HPfVwnfvpt2Ij6wtpgq9sOQFSgz4nWmnh7orAS0Cg5HuLsvzS55Xl0UCggKLnE7TQT45FhFFRp0GALQ5tNvaYqrI7qzdvLvuXUDdlZriXkK4wDTYeWwnH2z4wOpyRM5JgUVOt3Wr66MCS7ltviABgLaZu6wtpIq8sOQFSpwl/M60c5W6KzVC5K+6LCXOEosrEjk7BRYpJeHxL1jwn/kA/GtfEAmPf2FxRYFhU2PXXiw1MbDsztrNO2vfAdRdqWmGEkKjiEZsP7qd/274r9XliJyVAouc5rI01ynN65tcbHElgeMXd2CpiVNC3u5Kq9+pu1LD1MHg4cSHAfjXkn/hcDosrkjkzBRYpJQL8o4Rl3sYJwY/N25ldTkBY9MFpwQWp9PaYirR/pz9J7srPZ+xuBqpCvd1vY8G4Q3YemQrH/78odXliJyRAouU0iHd1V3Z3jBeFz30wc6GzSi0B1On6ATsqjnTQm+Mu5ASZwnXmHaumv57q8uRKlA3uRnDj+cDMO5/t2GO0ZmB4p8UWKSUy9yBZUMTLbj1hcNmZ2uj5q4v1q2ztphKcrz4OG+59115SNcMqtGGEkyYCSsNJ0vRtJD4JwUWKcWzfmVdrLbk95Vn4W1NCSzvr3+fo4ZJS9Pgj1q7UqNdgI1BBAMwDl1fSPyTAoucZJp0SHftwbJBgcVnNSmwmKbJ+B/HA/AAIdgxrC1IqtyD7i7ap5Sw/eh2i6sROZ0Ci5y0fz8XHM+i2Gb3nvUi5VeTAsvCnQvZdHgTdU24S9NBtUJb7PzeDMI0YMKPE6wuR+Q0Cixy0k8/AbC1UQsKg7Xfhq9+8ZxVtXs3ZGdbWsv58nRXBhNCPXVXao0R7nA6fe10jp44anE1IqUpsMhJK1cCsD5WC24rIiesDgfqXuD6Yv16a4s5D5sPb+bL7V9iYHC/uiu1ym+xc5lp43jxcaasmmJ1OSKlKLDISe4OizaMq7jNjRNcnwTwtNBry18D4I+X/JGL9COiVjEwvF2W11e8TpFDC3DFf+inkbiYpjoslSDQ17EcPXGUGetmAPBg9wetLUYsMZBgmtRpwsHcg9pITvyKAou47NwJWVkU2oPZekELq6sJWJ4dbwN1Sujfqf/mePFxLou5jGsTrrW6HLFACAbDug0D4NnFz1JQUmBxRSIuCizi4p4O2tS4JcX2YIuLCVzeDsuGDeAIrA24do2pyyvf/wuA0enbMJ6NtrYgscx9KS/SxDTYfnQ7L/0r2upyRAAFFvFwTweta6L9V87H7vpNIDwcTpyA7YG1l8UjFFBoQC/Tzo3aKK5Wi8JgPGEAvEAh245ss7giEQUW8XB3WLRh3Plx2uzQoYPriwBax5KyM4WPjRLsJkwgDEOnMtd6/Qmit2mnyIB/zP8HpmlaXZLUcgos4pq6WL0a0Jb8laJjR9fHAAksJc4Shi8YDriuKdMBu8UViT8wMJhEOKEmfLPzG2ZtnGV1SVLLKbAIbN0KeXkQEcGOhs2sribwBVhgefOnN/n50M80NA2edU8DiABciI2ncG0i+dBXD5FVkGVtQVKrKbCIdzqIzp1dUxpyfgIosBw+fpjRi0YD8C9CaaCpIPmVkYRwScNLyMjP4MmUJ60uR2oxBRbxLrila1dr66gpLrvM9XH/fjjq39ubj/5uNFkFWXSM6cg96OwwOV0oBm/2fROAN1e+yU8HfrK4IqmtFFhquYTHv2D9nC8BeGC7zgypDAkvLGFPdCwAd9w70eJqzmzHmLpM/cn1i2hC+jZdkVnOqNeMftxmBmNi8o+pV+JwBtYp+1IzVCiwTJo0iYSEBMLCwujevTsrVqw46/g5c+bQpk0bwsLC6NChA/Pnzy/1/Y8//pjrr7+ehg0bYhgGa9eurUhZUgHRJ3Jon74DgNTmHSyupuZYHt8egB57/HdaaAyFlBjwf2YQPXUas5zDK4RSz4SVhlPXGRJL+BxYZs+ezYgRI3jmmWdYvXo1HTt2pHfv3mRmZpY5ftmyZQwcOJAhQ4awZs0a+vXrR79+/di4caN3TH5+PldffTUvvfRSxd+JVEiPPeuxYbK5UQsO1WlgdTk1xg8tXOtYrvLTwPJz5s98QDHgWrsici6x2HjevSj7iW+fIDO/7J/5IlXF58Aybtw47rnnHgYPHky7du2YPHkyERERTJ8+vczxEyZMoE+fPowcOZK2bdvy3HPP0blzZyZOPNkqv/322xk9ejRJSUkVfydSIVfvXgPADwmdrC2khlnmDiztMnbCkSMWV3O6p797GtOAm8wgOus0ZimnoQRzuWkjqyCLRxc+anU5Usv4FFiKiopYtWpVqWBhs9lISkoiNTW1zMekpqaeFkR69+59xvHlUVhYSE5OTqmbVMzVu9cCsESBpVIdqtOALY2aY8OE776zupxSVh5cySebP8Ew4Z/qrogP7Bi84e6yzFg3g6V7l1pckdQmPgWWw4cP43A4iImJKXV/TEwM6enpZT4mPT3dp/HlkZycTFRUlPcWHx9f4eeq1XbsoHl2BkW2IFa411xI5fF0WUhJsbaQX3nq26cAuJ1g2qm7Ij66kiDu6XwPAEO/GEqxo9jiiqS2CMizhEaNGkV2drb3tm/fPqtLCkwLFwKwpmkbjoeEW1xMzfNDi06uT/wosHy/53u+2vEVQbYgnlF3RSoo+bpkGoY3ZGPmRi3AlWrjU2Bp1KgRdrudjIyMUvdnZGQQGxtb5mNiY2N9Gl8eoaGh1KtXr9RNKuCbbwBNB1WV5c3b4zBssG0b+EGodppOnkh5AoC7L7+bVoH594r4gYYvt+Kfx/MBeG7+/eSP0c9gqXo+/cQKCQmhS5cupJzyF6PT6SQlJYXExMQyH5OYmFhqPMDChQvPOF6qicMB334LnNIJkEqVGxrJes+1mfygyzL62XB+2PcDYSY8tfJ9q8uRAHc3wbQ0DTIMk4kUWV2O1AI+/4k1YsQIpk6dyowZM9i0aRNDhw4lPz+fwYMHAzBo0CBGjRrlHT98+HAWLFjA2LFj2bx5M2PGjGHlypUMGzbMO+bo0aOsXbuWX375BYAtW7awdu3a81rnIuewahUcO0ZOaCTrm+iCh1XlhwT/WMfy7tp3ed5w/VJ5kzCaqrsi5ykEg2fd04ovUajrDEmV8/mn1oABA3j11VcZPXo0nTp1Yu3atSxYsMC7sHbv3r2kpaV5x/fo0YOZM2cyZcoUOnbsyEcffcTcuXNp3/7kIs/PPvuMyy+/nL59+wJwyy23cPnllzN58uTzfX9yJu7poNTmHXDo+kFVxrMfC998A6ZpSQ2Ldi/ib5//DYAnzBDuJMSSOqTm+SvBtDNtHDNg7LKxVpcjNZxhmhb9FK1EOTk5REVFkZ2drfUs5dWrFyxaxFO/G8r7nftaXU2NFVpSxJZJf4WCAvj5Z2jXrlpff8vhLSROS+RYwTEGmEHMJBybtuCXSvQJxdxonCAyOJKdw3fSOLKx1SVJAPHl97f6wrVRfj788AMAS7XgtkoVBoXAVVe5vqjmaaHDxw/Td2ZfjhUcI7FZIu8orEgV6EcQV5g28ovzSV6SbHU5UoMpsNRGS5ZAcTE0b87u+nFWV1PzXXed62N1BpYxUQx7uQk7ju2gpWnw6b4NhCusSBUwMHjBvZncGyvfYG/2XosrkppKgaU2cu+/wu9+B4Z+iVU5z07PixZBSUm1vORnFDPbKMFuwhwiuED/qUsVSsLOtQnXUuQo4uGvH8ZpOq0uSWog/RSrjTyBRdduqh6dO0N0NGRnw+rVVf5y2QXZDKUAgIcJoYt2s5UqZmDw0q6fsJvw0S8fMepZbUQplU+BpZb5zb3TYMMGHIaNzsv0V1C1sNvh2mtdn8+fX+Uv9+jCRzlomFxk2hij3WylmnTDztvuqaGXjSJeW/6axRVJTaPAUsv8YfMSwHU689GIKIurqR0SHv+CB0suBGDXa1NIeGxelb3Wot2LmLLatVX624Rp3YpUqzsJ4XnTFZIfXPAgH/78ocUVSU2iwFLL/HHT9wB83vYaiyupXb5ufSXHg0NpeSyNTmlbq+Q1jhcf557PXRel+7sZTE+CquR1RM5mFCHcZwZjYnL7J7ezaPciq0uSGkKBpTbZtIl2mbsottlZcHEPq6upVY6HhPNVa9flKPr9/F2VvMbj3zzO9qPbaVq3KS+5W/Mi1c3AYAJh3NT2JoocRdww6wZWHVxldVlSAyiw1CazZwOwJOFyssPrWlxM7TP30l6Au8tVXFypz/3KmDBeX/E6AG/mHCVKU0FiITsG7//yNdeYdnIKc7h+Slc2Zm60uiwJcAostYVpegOLpoOssTShE4cio2l4Ige++qrSnnfa6mk8ahQC8JIZyh8JrrTnFqmoMAzmEUE308ZRwyTpvSS2Hqma6VCpHRRYaov162HzZgrtwSxsfaXV1dRKDpudz9r2dH3xfuVcLfnjTR/zt3mu6wQ9aobwqM4KEj9SF4MFRNLRtJGRn8F1713H7qzdVpclAUqBpbaYNQuAby/sSl5ohMXF1F6fuKeF+PRTyMk5r+f6dte3DPzfQJymkyFmMC8qrIgfqo/B10TQplEb9ufs57czfsv+nP1WlyUBSIGlNjBNb2DRdJC1NsZcyPYGzVwXQ/z44wo/z08HfuKGWTdQ5CjixrY3MpkwDK1bET/VGBsph/ZzoWmwK2sX14xrzq4xWkcnvlFgqQ1++gl274bISL698Aqrq6ndDONkl6WC00I/Z/5Mn6ndySvK47emnQ9+WUiQwor4uThsfEskF5k2dhkm15CvNS3iEwWW2sDdXeFPf6IgWKe7Wu3Tdu51LN9+CwcO+PTYXcd2cf3713PUMOlm2phLBGEKKxIgmmNjMRG0NW3sN0yueecanT0k5abAUtM5nd6zg7jlFmtrEQD2R8fC1Ve7pupmziz349Jy00j6TxIHcw9yqWnjSyKpq7AiASbOHVo8C3GvffdaVqdV/TW2JPApsNR0n34KBw9CVBT07m11NeL2eGQnAA4+9zKXPHzutSxHjh+h9/u92XlsJy2jW/I1ETRQWJEAdQE2viOSbk27ceTEEZLeS2JN2hqryxI/p8BSkxUWsueufwAwsV1vEp75xuKCxOOT9r/lQN0LiMs9zKDVX5x17JbDW7hy2pVsyNxAbJ1YFt6+kDj9pysBrj4GC/dvItG0c6zgGElvdWFt+lqryxI/pp96Ndn48bTISiejTgPeuLK/1dXIKQqDQvj3b24F4L4fP4SsrDLHLdq9iMSJbdh+dDstTIOU3FwufK1zNVYqUnXqYbCACK407Rw1TK577zrWpa+zuizxUwosNVVaGvzrXwC82PNOjoeEW1yQ/NrHl/ZiS6PmRBfkwUsvnfb9d9a8w+/+8zuOGXClaWc5kbTDbkGlIlXHE1q6m3aOnjjKde9dx/qM9VaXJX5IgaWmeuIJyMtjbZOLmXvptVZXI2Vw2uy83PMO1xfjx3vPGNqXvY+7P7ubuz67ixJnCQPMIL4lghj95yo1VBQGXxHhXdPS892efLD+A0zTtLo08SP6CVgT/fQTvPsuAM9e9zdMQ//M/irlwm6saNYOCgo4+Nyj3D//fi76d3OmrZkGwFNmCDMJJ1wLbKWGi8LgK/ealqyCLG775DZu/PBGMvIyrC5N/IRh1oAIm5OTQ1RUFNnZ2dSrV8/qcqxlmnDVVZCaCrffTkLcAKsrknO47MA6Lsl4kondoMB93cKepp3nCOU3BFlbnEg1K8bkJYr4J4UUG9AwvCGTfj+Jmy+9GcNQcK9pfPn9rT+9a5pp01xhJSICkpOtrkbOodDYwpYLn+bVq1xh5apM+NYZziIiFVakVgrG4ClC+YlIOsZ05MiJI9zyv1u4YdYN7MnaY3V5YiEFlhqk750TKLzXdRrzi11vJuH1tdYWJGdkUkJW0EzSQ0ey1XAS5zCYNxOWvAG9fnRaXZ6I5TpiZ0X6TkabIQSb8PnWz2k3PoGXf3iZYkex1eWJBRRYaoojR5g8N5lQRzHftrqCt7rfaHVFcgZFxg7SQx8lO3gmGE5uMYPYYKtD35ahrpUqCwthZ4nVZYpYLgSDZwljLZH8xrRz3IDHvnmMzlM68/2e760uT6qZAktN4HDArbcSn53BnuhYHvzjI1po64dKjEyuDvszaWHDKbJtJdqEmWY4//XsWts9BDoGgwl8dAKOqdMiAtAOO4uJ4B0zjIbhDdmYuZGe7/bk5jk3sztrt9XlSTXRb7WaYMwY+OorTgSF8vc/P0lOWB2rK5JTOMjhWNB0DoT+jQ8MVyv7VjOYDdRhIMEnBxoG/CEMmtrghAmzjkNRwK+JF6kUBgZ3EsKW44X83QzGZsKcX+bQZmIbnkx5ktzCXKtLlCqms4QC3ezZ3osaDv/Dw3x6aS+LCxKPQmM7uUHzOG7/HtMoAuC3pp1XCKPz2TaAy3HC1HzIM6FNEPwlHOw6O0LkVOtx8CAFfGc4AGgQ3oBBlw3i7s53c2njSy2uTsrLl9/fCiyBqqCAd66+mcGrPgfgnS5/5Nmkv1tclDjI4YR9BXn2BRTaN3vvv9y08Tyh9CEIozx7quwrgXePgxOItcGfw6GxdrkVOZWJyWeU8AiFbDdOTqEmNkvkrsvvomeLnlzU4CKdDu3HFFhqul9+gYEDYb1r++ppV9xA8rWDKbHrNNjqZmJSbOzlhP0nTthWUGjbDJ4fnGYQfwXuI4RE7OULKqfaWgxzC1zTQ3bgulC4MsQ1dSQiXiWYfEUJb1PM5zYTh+nwfi86LJor4q6ga1xXrk24lmtaXENYUJiF1cqpFFhqqqwsePttGD0aTpzgcEQUj/z+IRZdeIXVldUKJiZOsiiy7aDQtoVC2xaKbFtxGnmlxnU0bfyFYO4mmNjzXSaW54TPCmCb+6yhBDv8JtT10abgIvJr6Th5l2LmUsJaHBT+6j+TiOAIklol8fuLfs/vLvwdLaNbqgNjoSoPLJMmTeKVV14hPT2djh078vrrr9OtW7czjp8zZw5PP/00u3fvpnXr1rz00kv8/ve/937fNE2eeeYZpk6dSlZWFldddRVvvvkmrVu3Llc9NTqwmKZrI7gpUzjxwSzCSwoB+D7hch7uO4JDdepbXGDN4wkmxbb9FBv7KbbtpcjYTbFtD04j57TxYSb0Iog/EsQfCCK+steymyasLoavCsCz/URdAzoEu84q0lSRSJmKMPkZJz/h4EccfEUJB43Sv/IahDegS5MudGnShc5NOtP2grZcWP9CwoN1wdjqUKWBZfbs2QwaNIjJkyfTvXt3xo8fz5w5c9iyZQuNGzc+bfyyZcu45pprSE5O5g9/+AMzZ87kpZdeYvXq1bRv3x6Al156ieTkZGbMmEHLli15+umn2bBhA7/88gthYedu3dWowGKasHMn/PADLF0KixfD1q3eb29u1IIZXf7IrI7X69TlCjJx4OQ4TiMfh3GMEiPNeyu2pVFi7Mdp5Jf5WMOE1tjojp0rsdMdO5dhI7g6rvVz1AnLCuHnYig45f6GNmgd5Lo1t0OQ/loUKYuJyTqczKeE+ZSwAgfFZfznYmDQPKo5Fze8mBZRLWhWrxlN6zV1fazblKb1mlI/rL46M5WgSgNL9+7d6dq1KxMnTgTA6XQSHx/P/fffz+OPP37a+AEDBpCfn8+8efO891155ZV06tSJyZMnY5omcXFxPPzwwzzyyCMAZGdnExMTw7vvvsst7jNgKusNW8bhgBMn4PhxyMuD3FzXx8OHYds2VyjZts21PiUzs9RDTwSFMq/Nb5jZqQ9r4i6pdWsYTEzAgUkJJsWYFINRgkmhO3gUYFKAaRS4x5S4xxfhMI5SYhzBYRzBYRzGYeRgGsfP+ZqGCS0xaIOdttjogI322GmHzfoLEZaYrimidcWuj6du1xIMxNihjuG+2U7/PMKAIGrd/49Efq0Qk404WYWDlThYi4OtOMkux38a4UHhNK3XlKZ1mxJXN44mdZq4PtZtQr3QeoTaQwkLCiMsKIxgezB2w47NsGG32bEbdsKCwggNco0JtYcSbA/GVgv/CPXl97dPqzSLiopYtWoVo0aN8t5ns9lISkoiNTW1zMekpqYyYsSIUvf17t2buXPnArBr1y7S09NJSkryfj8qKoru3buTmppaZmApLCyksLDQ+3V2djbgeuNV5uab4ejR0veZJjidrpvD4boVF7tuRUWuW2GhK6iUlH/n0iLDzs8xrVjdtA1r49rwU7N25IVGur95ohLf1JmVGJkcDhkLp/1yNsq4z/TeTByAExMneG8mri1/7BjYAAMTxynBo8Q9xgBsGBjukFLs/l6x+/tnYP7qYzlFAA2AlthoiY1W7o+tsXEhNsLKCCbF7qos1yrIdSswYXcJ7Chx7Y6bD+x3nPPhgOufJAjXgl7PR7vh+uj6Z8L7b21wyvF1f3JVKFx0yj4yIgGotfu/+VvceyKZmBzBZDtOtuNkPyYHcXIQkwM4ScPkGHCCE2zP2852tldqPXabnSBbEEG2IOyG3RVw3PfZDXup+z0/K8/Wd/A8V5AtCLvNfs5QdOrzmaaJiUnrhq2Z+seplfo+PTy/t8vTO/EpsBw+fBiHw0FMTEyp+2NiYti8eXOZj0lPTy9zfHp6uvf7nvvONObXkpOTefbZZ0+7Pz4+vnxvxN+ZDkjf5rq5T1uuTaorDhx33/bjZAm1cFdZJ1D06zt/ffTP8q/xUQGl56ZE5Hw53P8rpPDcg6vJGtbwIR9W6Wvk5uYSFRV11jEBeR7sqFGjSnVtnE4nR48epWHDhj7NKebk5BAfH8++ffv8dyopQOnYVh0d26qh41p1dGyrTqAfW9M0yc3NJS4u7pxjfQosjRo1wm63k5GRUer+jIwMYmNjy3xMbGzsWcd7PmZkZNCkSZNSYzp16lTmc4aGhhIaGlrqvujoaF/eSin16tULyH/oQKBjW3V0bKuGjmvV0bGtOoF8bM/VWfHwaYVPSEgIXbp0ISUlxXuf0+kkJSWFxMTEMh+TmJhYajzAwoULveNbtmxJbGxsqTE5OTksX778jM8pIiIitYvPU0IjRozgjjvu4IorrqBbt26MHz+e/Px8Bg8eDMCgQYNo2rQpycnJAAwfPpyePXsyduxY+vbty6xZs1i5ciVTpkwBwDAMHnzwQf71r3/RunVr72nNcXFx9OvXr/LeqYiIiAQsnwPLgAEDOHToEKNHjyY9PZ1OnTqxYMEC76LZvXv3YrOdbNz06NGDmTNn8tRTT/HEE0/QunVr5s6d692DBeDRRx8lPz+fv/3tb2RlZXH11VezYMGCcu3Bcj5CQ0N55plnTptekvOnY1t1dGyrho5r1dGxrTq16djWiK35RUREpGarfbvUiIiISMBRYBERERG/p8AiIiIifk+BRURERPxerQ4skyZNIiEhgbCwMLp3786KFSusLimgJCcn07VrV+rWrUvjxo3p168fW7ZsKTWmoKCA++67j4YNG1KnTh1uuumm0zYSlHN78cUXvVsAeOjYVtyBAwe47bbbaNiwIeHh4XTo0IGVK1d6v2+aJqNHj6ZJkyaEh4eTlJTEtm3bLKzY/zkcDp5++mlatmxJeHg4F154Ic8991ypa8TouJbP999/zx//+Efi4uIwDMN77T2P8hzHo0ePcuutt1KvXj2io6MZMmQIeXl51fguqoBZS82aNcsMCQkxp0+fbv7888/mPffcY0ZHR5sZGRlWlxYwevfubb7zzjvmxo0bzbVr15q///3vzebNm5t5eXneMffee68ZHx9vpqSkmCtXrjSvvPJKs0ePHhZWHXhWrFhhJiQkmJdddpk5fPhw7/06thVz9OhRs0WLFuadd95pLl++3Ny5c6f51Vdfmdu3b/eOefHFF82oqChz7ty55rp168w//elPZsuWLc0TJ05YWLl/e/75582GDRua8+bNM3ft2mXOmTPHrFOnjjlhwgTvGB3X8pk/f7755JNPmh9//LEJmJ988kmp75fnOPbp08fs2LGj+eOPP5pLliwxL7roInPgwIHV/E4qV60NLN26dTPvu+8+79cOh8OMi4szk5OTLawqsGVmZpqAuXjxYtM0TTMrK8sMDg4258yZ4x2zadMmEzBTU1OtKjOg5Obmmq1btzYXLlxo9uzZ0xtYdGwr7rHHHjOvvvrqM37f6XSasbGx5iuvvOK9LysrywwNDTX/+9//VkeJAalv377mXXfdVeq+G2+80bz11ltN09RxrahfB5byHMdffvnFBMyffvrJO+bLL780DcMwDxw4UG21V7ZaOSVUVFTEqlWrSEpK8t5ns9lISkoiNTXVwsoCW3Z2NgANGjQAYNWqVRQXF5c6zm3atKF58+Y6zuV033330bdv31LHEHRsz8dnn33GFVdcQf/+/WncuDGXX345U6dO9X5/165dpKenlzq2UVFRdO/eXcf2LHr06EFKSgpbt24FYN26dSxdupT/+7//A3RcK0t5jmNqairR0dFcccUV3jFJSUnYbDaWL19e7TVXloC8WvP5Onz4MA6Hw7s7r0dMTAybN2+2qKrA5nQ6efDBB7nqqqu8uxinp6cTEhJy2oUpY2JiSE9Pt6DKwDJr1ixWr17NTz/9dNr3dGwrbufOnbz55puMGDGCJ554gp9++okHHniAkJAQ7rjjDu/xK+vng47tmT3++OPk5OTQpk0b7HY7DoeD559/nltvvRVAx7WSlOc4pqen07hx41LfDwoKokGDBgF9rGtlYJHKd99997Fx40aWLl1qdSk1wr59+xg+fDgLFy6s8ktU1DZOp5MrrriCF154AYDLL7+cjRs3MnnyZO644w6LqwtcH374IR988AEzZ87k0ksvZe3atTz44IPExcXpuEqlqJVTQo0aNcJut592RkVGRgaxsbEWVRW4hg0bxrx58/juu+9o1qyZ9/7Y2FiKiorIysoqNV7H+dxWrVpFZmYmnTt3JigoiKCgIBYvXsxrr71GUFAQMTExOrYV1KRJE9q1a1fqvrZt27J3714A7/HTzwffjBw5kscff5xbbrmFDh06cPvtt/PQQw95L4Sr41o5ynMcY2NjyczMLPX9kpISjh49GtDHulYGlpCQELp06UJKSor3PqfTSUpKComJiRZWFlhM02TYsGF88sknfPvtt7Rs2bLU97t06UJwcHCp47xlyxb27t2r43wO1113HRs2bGDt2rXe2xVXXMGtt97q/VzHtmKuuuqq006/37p1Ky1atACgZcuWxMbGljq2OTk5LF++XMf2LI4fP17qwrcAdrsdp9MJ6LhWlvIcx8TERLKysli1apV3zLfffovT6aR79+7VXnOlsXrVr1VmzZplhoaGmu+++675yy+/mH/729/M6OhoMz093erSAsbQoUPNqKgoc9GiRWZaWpr3dvz4ce+Ye++912zevLn57bffmitXrjQTExPNxMREC6sOXKeeJWSaOrYVtWLFCjMoKMh8/vnnzW3btpkffPCBGRERYb7//vveMS+++KIZHR1tfvrpp+b69evNG264QaffnsMdd9xhNm3a1Hta88cff2w2atTIfPTRR71jdFzLJzc311yzZo25Zs0aEzDHjRtnrlmzxtyzZ49pmuU7jn369DEvv/xyc/ny5ebSpUvN1q1b67TmQPb666+bzZs3N0NCQsxu3bqZP/74o9UlBRSgzNs777zjHXPixAnzH//4h1m/fn0zIiLC/POf/2ympaVZV3QA+3Vg0bGtuM8//9xs3769GRoaarZp08acMmVKqe87nU7z6aefNmNiYszQ0FDzuuuuM7ds2WJRtYEhJyfHHD58uNm8eXMzLCzMbNWqlfnkk0+ahYWF3jE6ruXz3Xfflfmz9Y477jBNs3zH8ciRI+bAgQPNOnXqmPXq1TMHDx5s5ubmWvBuKo9hmqdsQygiIiLih2rlGhYREREJLAosIiIi4vcUWERERMTvKbCIiIiI31NgEREREb+nwCIiIiJ+T4FFRERE/J4Ci4iIiPg9BRYRERHxewosIiIi4vcUWERERMTvKbCIiIiI3/t/38IJLQNRG4EAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "# TODO: use plt.hist() and bins = len(np.unique()) in order to visualize the distribution\n",
    "n,x,_ = plt.hist(X, bins = len(np.unique(X)), density=True)\n",
    "bin_centers = 0.5*(x[1:]+x[:-1])\n",
    "plt.plot(bin_centers, n, color = 'red', label = 'Age')\n",
    "\n",
    "n,x,_ = plt.hist(Y, bins = len(np.unique(Y)), density=True)\n",
    "bin_centers = 0.5*(x[1:]+x[:-1])\n",
    "plt.plot(bin_centers,n, color = 'green', label = 'Weight')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**d)** What is the probability of observing an age of 80 or higher?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HINT: To calculate the probability we have to use the probability density function of normal distribution for $X$, in range $[80,\\infty)$. Calculate the sample PDF and the expected PDF given that the distribution is normal with ($\\mu_Y=60,\\sigma_Y=10$). How close are these two values?\n",
    "\n",
    "**Answer**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sampled PDF is 0.02016\n",
      "Expected PDF is 0.022759418529548325\n",
      "0.022759418529548325\n",
      "0.0025994185295483245\n"
     ]
    }
   ],
   "source": [
    " # TODO: find the probability that age is >80 (PDF), then find the expected PDF of the normal distribution (μ=60,σ=10) and compare the results\n",
    "Y_greater_80 = Y[Y>80]\n",
    "pdf = (len(Y_greater_80)/len(Y))\n",
    "print(f'Sampled PDF is {pdf}')\n",
    "#to get any range ex {80-X}\n",
    "y1 = scipy.stats.norm(y_mean, y_variance**(1/2)).cdf(float('inf'))\n",
    "y2 = scipy.stats.norm(y_mean, y_variance**(1/2)).cdf(80)\n",
    "epdf = y1-y2\n",
    "print(f'Expected PDF is {epdf}')\n",
    "#to get 80-inf\n",
    "print(1-scipy.stats.norm(y_mean, y_variance**(1/2)).cdf(80))\n",
    "print(np.abs(pdf-epdf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**e)** Find the 2-dimensional mean $\\hat \\mu$ and the covariance matrix $\\widehat \\Sigma$ for these two variables.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean:\n",
      " [[25.005264 60.002704]]\n",
      "Cov: \n",
      " [[ 2.50500603e+01 -5.56962339e-02]\n",
      " [-5.56962339e-02  9.99901567e+01]]\n"
     ]
    }
   ],
   "source": [
    "# TODO: Compute 2-Dimensional mean is the vector (E[X], E[Y])\n",
    "Z = np.c_[X, Y]\n",
    "n, d = Z.shape\n",
    "mu = np.mean(Z, axis=0, keepdims=True)\n",
    "cov = (Z - mu).T @ (Z - mu)/n\n",
    "\n",
    "# TODO: The diagonal of the Covariance matrix must equal the initial covariance values of X and Y. Check your results.\n",
    "print(\"Mean:\\n\", mu)\n",
    "print(\"Cov: \\n\", cov)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**f)** Normalize the data with _range normalization_ to the range $[0, 1]$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**range normalization**:\n",
    "\n",
    "$x_i^{'} = \\frac{x_i-\\min_i{x_i}}{\\max_i{x_i} - \\min_i{x_i}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Range normalized: \n",
      " [[0.67391304 0.56666667]\n",
      " [0.43478261 0.27777778]\n",
      " [0.43478261 0.43333333]\n",
      " ...\n",
      " [0.56521739 0.42222222]\n",
      " [0.32608696 0.62222222]\n",
      " [0.56521739 0.52222222]]\n"
     ]
    }
   ],
   "source": [
    "# TODO: Compute the range normalization: (x_i - min)/(max - min)\n",
    "min_z = Z.min(0, keepdims=True)\n",
    "max_z = Z.max(0, keepdims=True)\n",
    "Z_normalize = (Z - min_z) / (max_z - min_z)\n",
    "print(\"Range normalized: \\n\", Z_normalize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**g)** Normalize the data with _standard score normalization_ such that is has mean 0 and standard deviation 1.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**z-score**:\n",
    "\n",
    "$x_i^{'} = \\frac{x_i-\\hat{\\mu}}{\\hat{\\sigma}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.283844142108137e-16 1.0\n",
      "-7.653255806872039e-17 1.0\n"
     ]
    }
   ],
   "source": [
    "# TODO: #  Compute the standard score normalization:, x_i = (x_i - mu)/(sigma)\n",
    "X_ssn = [(x - x_mean)/np.sqrt(x_variance) for x in X]\n",
    "# Cheating and computing mean of Y with Numpy\n",
    "Y_ssn = [(y - np.mean(Y))/np.sqrt(y_variance) for y in Y]\n",
    "print(np.mean(Y_ssn), np.std(Y_ssn))\n",
    "print(np.mean(X_ssn), np.std(X_ssn))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2: Robustness\n",
    "Determine whether the following statements are true or false and explain why.\n",
    " - Mean is robust against outliers\n",
    " - Median is robust against outliers\n",
    " - Standard deviation is robust against outliers\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "\n",
    "\n",
    "\n",
    " - Mean is robust against outliers \n",
    " \n",
    "      **False**: We say that a statistic is robust if it is not affected by extreme values (such as outliers) in the data. The sample mean is unfortunately not robust because a single large value (an outlier) can skew the average.\n",
    " \n",
    " <em>Data Mining and Analysis Page 47</em>\n",
    " - Median is robust against outliers\n",
    " \n",
    "      *True：* It is not affected by extreme values\n",
    "\n",
    "- Standard deviation is robust against outliers\n",
    "\n",
    "    **False**: Standard deviation $\\sigma = \\sqrt\\frac{\\sum(x_i-\\mu)^2}{N}$, it could influenced by each $x_i$ \n",
    "      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3: More on robustness\n",
    "Provide an informal definition of when a measure is robust. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "\n",
    "It is not affected by extreme values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4: Independence analysis\n",
    "Please explain what the idea of independence analysis is using contingency tables; for which kind of data is this particularly relevant? \n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "\n",
    "Determine if two categorical variables are independent or if they are in fact related to one another. If two categorical variables are independent, then the value of one variable does not change the probability distribution of the other. If two categorical variables are related, then the distribution of one depends on the level the other. This test measures the differences in the observed conditional distribution of one variable across levels of the other, and compares it to the marginal (overall) distribution of that variable.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5: Contingency table and $\\chi^2$ statistics\n",
    "In the table below,  assuming that $X_1$ is discretized into three bins, as follows: $[-2, -0.5]$, $[-0.5, 0.5]$, and $[0.5, 2]$.\n",
    "\n",
    "![Table 3.11](graphics/3.11.png)\n",
    "\n",
    "Answer the following questions\n",
    "\n",
    "**a)** Construct the contingency table between the discretized $X_1$ and $X_2$ attributes, including the row and column counts.<br>\n",
    "**b)** Construct the expected table between the discretized $X_1$ and $X_2$ attributes.<br>\n",
    "**c)** Compute the $\\chi^2$ statistic between them.<br>\n",
    "**d)** Determine whether they are dependent or not at the 5% significance level. Use the $\\chi^2$ critical values from Table 3.10.\n",
    "\n",
    "![Table 3.10](graphics/3.10.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "**a)**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X1: [1, 1, 1, 0, 1, 2, 1, 0, 2, 0]\n",
      "X2: [0, 1, 0, 0, 0, 1, 0, 1, 1, 0]\n",
      "Our table before row counts: [[2, 1], [4, 1], [0, 2]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "X_1 = [0.3, -0.3, 0.44, -0.60, 0.40, 1.20, -0.12, -1.60, 1.60, -1.32]\n",
    "X_2 = [\"a\", \"b\", \"a\", \"a\", \"a\", \"b\", \"a\", \"b\", \"b\", \"a\"]\n",
    "# For sanity's sake, use code to bin X_1\n",
    "# Converting everything to numbers allow us to index in them for easier computation\n",
    "X_tmp = []\n",
    "for x in X_1:\n",
    "    if x >= -2 and x <= -0.5:\n",
    "        X_tmp.append(0)\n",
    "    elif x >= -0.5 and x <= 0.5:\n",
    "            X_tmp.append(1)\n",
    "    elif x > 0.5:\n",
    "            X_tmp.append(2)\n",
    "X_1 = X_tmp\n",
    "X_2 = [0 if x == \"a\" else 1 for x in X_2]\n",
    "print(f\"X1: {X_1}\")\n",
    "print(f\"X2: {X_2}\")    \n",
    "\n",
    "# Our table looks like\n",
    "## Bins: a b\n",
    "## bin1\n",
    "## bin2\n",
    "## bin3\n",
    "b1_row = [0, 0]\n",
    "b2_row = [0, 0]\n",
    "b3_row = [0, 0]\n",
    "count_row = [0, 0, 0] # THIS ONE is for column counts\n",
    "c_table = [b1_row, b2_row, b3_row]\n",
    "# TODO: Compute contingency table c_table \n",
    "\n",
    "for i in range(len(X_1)):\n",
    "    x1_index = X_1[i]\n",
    "    x2_index = X_2[i]\n",
    "    c_table[x1_index][x2_index] += 1\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "print(f\"Our table before row counts: {c_table}\")   "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Table: Contingency table\n",
    "\n",
    "|  | [-2,-0.5] | [-0.5,0.5] | [0.5, 2] |row counts|\n",
    "|---|-----------|------------|----------|----|\n",
    "| a | 2         |4          |0      | 6  |\n",
    "| b | 1         |1          |2         | 4  |\n",
    "| column counts  | 3         | 5          | 2        | 10 |"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "**Answer:**\n",
    "**b)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The expected table is: [[1.8, 1.2], [3.0, 2.0], [1.2, 0.8]]\n"
     ]
    }
   ],
   "source": [
    "# Now add row counts\n",
    "# Remember these are reference variables, so modifying these\n",
    "# will also modify our table\n",
    "\n",
    "# TODO: Compute expected table\n",
    "E = [[0, 0] for _ in range(3)]\n",
    "b1_row.append(sum(b1_row))\n",
    "b2_row.append(sum(b2_row))\n",
    "b3_row.append(sum(b3_row))\n",
    "c_table.append([b1_row[i] + b2_row[i] + b3_row[i] for i in range(3)])\n",
    "\n",
    "\"\"\"\n",
    "print(f\"Our contingency table:\")\n",
    "print(\".......a, b, counts\")\n",
    "print(f\"bin1: {b1_row}\")\n",
    "print(f\"bin2: {b2_row}\")\n",
    "print(f\"bin3: {b3_row}\")\n",
    "print(f\"count: {c_table[3]}\")\n",
    "\"\"\"\n",
    "\n",
    "# Construct matrix of e_{i,j} needed for Chi^2\n",
    "E = [[0, 0] for _ in range(3)]\n",
    "for i in range(3):\n",
    "    for j in range(2):\n",
    "        E[i][j] = (c_table[i][2] * c_table[3][j])/(len(X_1))\n",
    "        \n",
    "\n",
    "\n",
    "print(f\"The expected table is: {E}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "**c)**\n",
    "|  | [-2,-0.5] | [-0.5,0.5] | [0.5, 2] |row counts|\n",
    "|---|-----------|------------|----------|----|\n",
    "| a | 1.8         |3          |1.2       | 6  |\n",
    "| b | 1.2         |2          |0.8         | 4  |\n",
    "| column counts  | 3         | 5          | 2        | 10 |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi_squared is: 0\n",
      "Chi_squared is: 3.8888888888888884\n"
     ]
    }
   ],
   "source": [
    "# TODO: Compute Chi_squared\n",
    "Chi_squared = 0\n",
    "print(f\"Chi_squared is: {Chi_squared}\")\n",
    "for i in range(3):\n",
    "    for j in range(2):\n",
    "        Chi_squared += ((c_table[i][j]-E[i][j])**2)/(E[i][j])\n",
    "print(f\"Chi_squared is: {Chi_squared}\")\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "**d)**\n",
    "\n",
    "null hypothesis $H_0$, $X_1$ and $X_2$ are independent\n",
    "\n",
    "**Chi-square**\n",
    "\n",
    "$X^2 = \\sum_{i_1}^{m_1}\\sum_{j_1}^{m_2} \\frac{(n_{ij}-e_{ij})^2}{e_{ij}}$\n",
    "\n",
    "$X^2 = 0.022+0.333+0+0.333+0.5+1.8 = 3.886$\n",
    "\n",
    "**degree of freedom**\n",
    "\n",
    "$q = (2-1)*(3-1) = 2$\n",
    "\n",
    "**cirtical value** by table $= 5.991$\n",
    "\n",
    "$X^2 \\le 5.991 $ it accept the null hypothesis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 6: Working with Metrics\n",
    "\n",
    "Consider the following situation: \n",
    "\n",
    "We know some distances between data points: $d(p_1,p_2)=1$, $d(p_1,p_3)=2$, $d(p_3, p_4)=1$. \n",
    "We also know that $d$ is a metric. \n",
    "\n",
    "1. What do we know about the remaining distances?\n",
    "2. Do we need to compute further distances if we want to find the two points that are most similar to $p_1$? \n",
    "3. Can $p_4$ be closer to $p_2$ than $p_3$ is to $p_2$?\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answers:**\n",
    "\n",
    "    1. dist(x,z)  ≤ dist(x, y) + dist(y, z)\n",
    "A1. $d(p_1,p_4)\\leq d(p_1,p_3)+d(p_3,p_4)$ ->$d(p_1,p_4)\\leq 3$\\\n",
    "$d(p_1,p_4)\\leq d(p_1,p_2)+d(p_2,p_4)$ ->$d(p_1,p_4)\\leq 1+d(p_2,p_4)$\\\n",
    "$d(p_2,p_3)\\leq d(p_2,p_1)+d(p_1,p_3) $->$d(p_2,p_3)\\leq 3$\\\n",
    "$d(p_2,p_4)\\leq d(p_2,p_1)+d(p_1,p_4)\\leq d(p_1,p_2)+d(p_1,p_3)+d(p_3,p_4)$->$d(p_1,p_2)+d(p_1,p_4)\\leq 4$\\\n",
    "\n",
    "\n",
    "$2.$\n",
    "Yes, since $p_4$ can be closer to $p_1$ than $p_3$ but also the opposite scenario\n",
    " \n",
    "\n",
    "$3.$\n",
    "Yes $p_4$ can be equal to $p_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optional Exercises\n",
    "\n",
    "## Exercise 7: Mean absolute deviation \n",
    "Define a measure of deviation called _mean absolute deviation_ for a random variable $X$ as follows:\n",
    "\n",
    "$$\n",
    "\\frac{1}{n} \\sum_{i=1}^n |x_i - \\mu|\n",
    "$$\n",
    "\n",
    "Is this measure robust? Why or why not?"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer:**\n",
    "\n",
    "It is not robust, as a single outlier can skew the mean absolute deviation. (Similar to mean and variance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 8: Correlation\n",
    "\n",
    "Consider the table below. Assume that both the attributes $X$ and $Y$ are numeric, and the\n",
    "table represents the entire population. And we know that the correlation between $X$\n",
    "and $Y$ is zero.\n",
    "1. What can you infer about the values of $Y$? (Find a relationship between $a$, $b$ and $c$)\n",
    "2. If we know there is a missing row, what can we infer about it?\n",
    "\n",
    "|$X$ | $Y$ |\n",
    "|:---:|:---:|\n",
    "|$1$|$a$|\n",
    "|$0$|$b$|\n",
    "|$1$|$c$|\n",
    "|$0$|$a$|\n",
    "|$0$|$c$|"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answer**\n",
    "\n",
    "    1. since it is independent, we have\n",
    "\n",
    "$E[XY] = E[X]E[Y]$\n",
    "\n",
    "and\n",
    "\n",
    "$E[XY] = \\frac{a+c}{5}$\n",
    "\n",
    "$E[X] = \\frac{2}{5}$\n",
    "\n",
    "$E[Y] = \\frac{2a+b+2c}{5}$\n",
    "\n",
    "We have $\\frac{a+c}{5}=\\frac{2}{5}\\cdot\\frac{2a+b+2c}{5}$\n",
    "\n",
    "solve this we got, $2b=a+c$\n",
    "\n",
    "\n",
    "\n",
    "    2.\n",
    "|  | 1 | 0 |\n",
    "|-----|---|---|\n",
    "| a   | 1 | 1 |\n",
    "| b   | - | 1 |\n",
    "| c   | 1 | 1 |\n",
    "\n",
    "from the contigency table we could figure $[1,b]$ are now empty and could add one more line.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 9: 3-way contingency table\n",
    "\n",
    "Consider the 3-way contingency table for attributes $X$,$Y$,$Z$ shown in the table below.\n",
    "Compute the $\\chi^2$ metric for the correlation between $Y$ and $Z$. Are they dependent\n",
    "or independent at the 95% confidence level? See Table 3.10 above for $\\chi^2$ values.\n",
    "\n",
    "![Table 3.9](./graphics/3.9.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.957264957264956\n",
      "Chi2ContingencyResult(statistic=0.6515384615384604, pvalue=0.4195631759021463, dof=1, expected_freq=array([[37.14285714, 27.85714286],\n",
      "       [42.85714286, 32.14285714]]))\n"
     ]
    }
   ],
   "source": [
    "import scipy.stats.contingency\n",
    "from scipy.stats import chi2_contingency\n",
    "# TODO: Compute contingency using chi2_contingency\n",
    "xy = np.array([[40,25],[40,35]])\n",
    "chi=chi2_contingency(xy)\n",
    "chi2 = np.sum((xy - chi[3])**2/chi[3])\n",
    "print(chi2)\n",
    "#print(x)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We obtain $\\chi^2 = 0.957$. With one degree of freedom, we have $\\chi^2$ being at left of the $1-0.95=0.05$ critical $\\chi^2$. Thus we cannot reject the null hypothesis, and we concludethat the two variables are independent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 10: Mixed data\n",
    "\n",
    "Consider the \"mixed\" data given in the table below. Here $X_1$ is a numeric attribute and\n",
    "$X_2$ is a categorical one. Assume that the domain of $X_2$ is given as $dom(X_2) = \\{a, b\\}$.\n",
    "Answer the following questions.\n",
    "\n",
    "**a)** What is the mean vector for this dataset?  \n",
    "**b)** What is the covariance matrix?  \n",
    "\n",
    "![Table 3.11](graphics/3.11.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean vector:\n",
      " [0.144 0.6   0.4  ]\n"
     ]
    }
   ],
   "source": [
    "fine = np.array([[0.3,1,0],[-0.3,0,1],[0.44,1,0],[-0.6,1,0],[0.4,1,0],[1.2,0,1],[-0.12,1,0],[-0.16,0,1],[1.6,0,1],[-1.32,1,0]])\n",
    "# TODO: compute the mean vector\n",
    "mu = np.mean(fine, axis=0)\n",
    "print(\"mean vector:\\n\", mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cov: [[ 0.71873778 -0.196       0.196     ]\n",
      " [-0.196       0.26666667 -0.26666667]\n",
      " [ 0.196      -0.26666667  0.26666667]]\n",
      "Cov1: [[ 0.71873778 -0.196       0.196     ]\n",
      " [-0.196       0.26666667 -0.26666667]\n",
      " [ 0.196      -0.26666667  0.26666667]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "n, d = fine.shape\n",
    "# TODO: compute the covariance matrix\n",
    "cov = (fine - mu).T @ (fine - mu)/(n-1)\n",
    "print(f\"Cov: {cov}\")\n",
    "#print(f\"Cov1: {np.cov(fine.T)}\")"
   ]
  }
 ],
 "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.10.9 (tags/v3.10.9:1dd9be6, Dec  6 2022, 20:01:21) [MSC v.1934 64 bit (AMD64)]"
  },
  "vscode": {
   "interpreter": {
    "hash": "80c6980e3f02e16e7aee110d1bf85cdb98492d897883fc820eec5f24d661294b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}