{ "cells": [ { "cell_type": "markdown", "id": "8119837f", "metadata": {}, "source": [ "# Tutorial exercises: Sampling\n", "\n", "In these exercises we again work with the Brexdex data\n", "\n", "We are going to investigate how the sampling distribution of the mean depends on $n$, the relationship between SEM and $\\sqrt{n}$, and how we assess whether a distribution, such as the sampling distribution of the mean, is Normal." ] }, { "cell_type": "markdown", "id": "79456fc9", "metadata": {}, "source": [ "### Set up Python libraries\n", "\n", "As usual, run the code cell below to import the relevant Python libraries" ] }, { "cell_type": "code", "execution_count": 1, "id": "5cd335fe", "metadata": {}, "outputs": [], "source": [ "# Set-up Python libraries - you need to run this but you don't need to change it\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as stats\n", "import pandas as pd\n", "import seaborn as sns\n", "sns.set_theme(style='white')\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "import warnings \n", "warnings.simplefilter('ignore', category=FutureWarning)" ] }, { "cell_type": "markdown", "id": "fcf4f280", "metadata": {}, "source": [ "## Import and plot the data\n", "\n", "Let's remind ourselves of the dataset we are working with" ] }, { "cell_type": "code", "execution_count": 13, "id": "466ab767", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ID_codescore
018664053
158814090
297739030
394847042
456436084
.........
999585178081
999669834045
999769358051
999887273078
999938564288
\n", "

10000 rows × 2 columns

\n", "
" ], "text/plain": [ " ID_code score\n", "0 186640 53\n", "1 588140 90\n", "2 977390 30\n", "3 948470 42\n", "4 564360 84\n", "... ... ...\n", "9995 851780 81\n", "9996 698340 45\n", "9997 693580 51\n", "9998 872730 78\n", "9999 385642 88\n", "\n", "[10000 rows x 2 columns]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "UKBrexdex=pd.read_csv('https://raw.githubusercontent.com/jillxoreilly/StatsCourseBook_2025/main/data/UKBrexdex.csv')\n", "UKBrexdex" ] }, { "cell_type": "markdown", "id": "0e76cd91-aded-467b-8e5b-a171d90e8cc2", "metadata": {}, "source": [ "* What variables are in the dataset?\n", "* How many people are in the sample?\n", "\n", "and let's plot the data" ] }, { "cell_type": "code", "execution_count": 14, "id": "9f1bcab5-ffd8-4aa2-bd6a-569c348d59b8", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAG1CAYAAAAWb5UUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA39klEQVR4nO3deVyVZd7H8S/LAKKi4iD4OJO5hIqiomCUC0qZTWql1pSKk0qmaZqlaC6TWmqWC4l77s+U6aTmlFmNltlkZmqLe+ZGliKGCm6AwPX80cMZDmB69MAB7s/79eL1kuu+zzm/c0ecL9dyX27GGCMAAACLcHd1AQAAAMWJ8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACzF09UFlDTh4eHKzMxUQECAq0sBAAA36MyZM/Ly8tLOnTuvey7hJ5+MjAxlZ2e7ugwAAOCArKws3eimFYSffKpVqyZJ+uSTT1xcCQAAuFH33HPPDZ/LnB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGApnq4uAACA4rDus8NKOntZkhTk76uH29Z1cUVwFcIPAMASks5e1k9JF1xdBkoAhr0AAIClEH4AAIClMOwFAMA15J0nJDFXqKwg/AAAcA3MEyqbGPYCAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWwsamAAD8v7y7uIfc7u/ialBUCD8AAPy/vLu4B/r7urgaFBWGvQAAgKUQfgAAgKUQfgAAgKUQfgAAluPu7ubqEuBCTHgGAFhOQOVydiu7JFZ3WQnhBwBgSXlXdkms7rIShr0AAIClEH4AAICluDz8nD9/Xi+++KLatGmjZs2aqXv37tq5c6ft+KhRo1SvXj27rzZt2tiO5+TkKCEhQa1bt1aTJk3Ut29fJSYmuuKtAACAUsDlc36ef/55paSkaMaMGfL399eKFSsUGxurtWvXqk6dOvrhhx80YMAAxcTE2B7j4eFh+/fcuXO1cuVKvfLKKwoMDNTUqVPVr18/rV+/Xl5eXq54SwAAoARzac9PYmKitm7dqnHjxik8PFy1a9fWmDFjFBgYqPXr1ys7O1uHDx9WaGioAgICbF/+/r/NyM/MzNSSJUs0ePBgRUVFqX79+oqPj9fp06e1ceNGV741AABQQrk0/FSpUkVvvPGGGjVqZGtzc3OTMUapqak6fvy4MjIyVKdOnUIff/DgQV26dEmRkZG2Nj8/P4WEhGjHjh1FXj8AACh9XDrs5efnp6ioKLu2Dz/8UD/99JNatWqlQ4cOyc3NTcuXL9fnn38ud3d3RUVFaejQoapYsaKSkpIkSdWrV7d7jmrVqunUqVPF9j4AANbAzRHLBpfP+clr165dGj16tO655x5FR0crISFB7u7uqlGjhubPn6/ExES9+uqrOnTokJYvX64rV65IUoG5Pd7e3kpNTXXFWwAAlGH5b44Y5O+rh9vWdXFVcFSJCT+bNm3S8OHD1aRJE82YMUOSNHjwYPXu3Vt+fn6SpODgYAUEBOixxx7Tnj175OPjI+m3uT+5/5akjIwMlStXrvjfBACgzMt/c0SUPi5f6i5Jb775pgYPHqw2bdpo4cKFtiDj5uZmCz65goODJUlJSUm24a7k5GS7c5KTkxUUFFQMlQMAgNLG5eFnxYoVevnll9WzZ0+9/vrrdkNYw4YNU2xsrN35e/bskSTVrVtX9evXV4UKFbR9+3bb8bS0NO3fv1/h4eHF8wYAAECp4tJhr2PHjmny5Mlq3769+vfvr5SUFNsxHx8fderUSU8//bTmzZunjh076tixY3rppZfUqVMn2wqwmJgYTZs2Tf7+/qpRo4amTp2qoKAgtW/f3lVvCwAAlGAuDT8ff/yxrl69qo0bNxa4L0+XLl00ZcoUzZw5U/Pnz9f8+fNVsWJFde7cWUOHDrWdN2TIEGVlZWns2LFKT09XRESEFi9ezA0OAQBAoVwafgYMGKABAwb87jkdOnRQhw4drnncw8NDcXFxiouLc3Z5AACgDHL5nB8AAIDiRPgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAACWQvgBAOAmubu7uboE3ARPVxcAAEBpFVC5nNZ9dlhJZy/b2oL8ffVw27ourArXQ/gBAOAWJJ29rJ+SLri6DDiAYS8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAGAphB8AAJzI3d3N1SXgOjxdXQAAAGVJQOVyWvfZYSWdvSxJCvL31cNt67q4KuRF+AEAwMmSzl7WT0kXXF0GroFhLwAAYCn0/AAAypy8w06SFHK7vwurQUlD+AEAlDn5h50C/X1dWA1KGoa9AACApRB+AACApRB+AACApbh8zs/58+c1Y8YMffbZZ7p48aLq1aunYcOGKTw8XJJ04MABTZo0SXv37lXlypXVq1cvxcbG2h6fk5Oj2bNn65133lFaWpqaN2+ucePGqWbNmq56S0CRyT+Jk/uHAIDjXN7z8/zzz+v777/XjBkztHr1ajVs2FCxsbE6cuSIzp07pz59+uj222/XmjVrNHjwYM2cOVNr1qyxPX7u3LlauXKlJk6cqFWrVsnNzU39+vVTZmamC98VUDRyJ3HmfuUNQgCAG+PSnp/ExERt3bpVb7/9tpo1ayZJGjNmjD7//HOtX79ePj4+8vLy0vjx4+Xp6ak6deooMTFRCxcuVLdu3ZSZmaklS5YoLi5OUVFRkqT4+Hi1bt1aGzduVMeOHV359gAAQAnk0p6fKlWq6I033lCjRo1sbW5ubjLGKDU1VTt37lRERIQ8Pf+b0SIjI3Xs2DGlpKTo4MGDunTpkiIjI23H/fz8FBISoh07dhTrewEAAKWDS3t+/Pz8bD02uT788EP99NNPatWqleLj4xUcHGx3vFq1apKkkydPKikpSZJUvXr1AuecOnWqCCsHbl3e+Tsht/vrbFo683kAoBi4fM5PXrt27dLo0aN1zz33KDo6Wunp6fLy8rI7x9vbW5KUkZGhK1euSFKh52RkZBRP0cBNyjt/J+X/gw/zeQCg6JWY8LNp0ybFxsaqcePGmjFjhiTJx8enwMTl3FDj6+srHx8fSSr0nHLlyhVD1QAAoLQpEeHnzTff1ODBg9WmTRstXLjQFmqCgoKUnJxsd27u94GBgbbhrsLOCQoKKobKAQBAaePy+/ysWLFCL7/8snr16qXRo0fL3f2/eSwiIkIrV65Udna2PDw8JEnbtm1TrVq1VLVqVVWsWFEVKlTQ9u3bddttt0mS0tLStH//fsXExLjk/QAAihebmMJRLg0/x44d0+TJk9W+fXv1799fKSkptmM+Pj7q1q2bFi1apDFjxujJJ5/U7t27tXz5ck2YMEHSb3N9YmJiNG3aNPn7+6tGjRqaOnWqgoKC1L59e1e9LQBAMWITUzjKpeHn448/1tWrV7Vx40Zt3LjR7liXLl00ZcoULVq0SJMmTVKXLl0UEBCgESNGqEuXLrbzhgwZoqysLI0dO1bp6emKiIjQ4sWLC0yCBgAAkFwcfgYMGKABAwb87jmNGzfWqlWrrnncw8NDcXFxiouLc3Z5AACgDHL5nB+gLMo7B4H79QBAyUL4AYpA/jkIAICSo0QsdQcAACguhB8AAGAphB8AAGAphB8AAGAphB8AAGAphB+giLm7u7m6BABAHix1Bxzk6D18AiqXY+8hAChBCD+Ag27mHj7sPQQAJQfDXsAtcPWQlqtfHwBKI3p+gFvg6iGt/K/PVhoAcH2EH+AWuXpIi600AMAxDHsBAABLIfwAAABLIfwAAABLIfwAAABLIfwAAABLIfwAAABLIfwAAABLIfwAAABLIfwAAABLIfwAAABLYXsL4He4ct8uAEDRIPwAv8PV+3Y5il3eAeD6CD9ACXUzQaawXebZ6R1lTd6fcXpjcTMIP0AJlT/I3OgveWft8k6IQkmV92e8pPfGomQi/AAlmCt/yTsrRAFASePwaq/09PSiqAMAAKBYOBx+7r77bo0dO1bffPNNUdQDAABQpBwOPwMGDNC3336rHj16qEOHDlqwYIGSkpKKojYAAACnc3jOz1NPPaWnnnpKu3fv1tq1a7VkyRIlJCQoMjJS3bp107333isvL6+iqBW4YXkn6zJRFwCQ101PeG7cuLEaN26sMWPGaOvWrVq8eLGGDRumihUr6qGHHtLf/vY3/fnPf3ZmrcANY7IuAOBabml7i5MnT2rx4sWKj4/Xjh07VKtWLXXr1k1fffWVOnbsqPXr1zurTgAAAKdwuOfn4sWL+vjjj7Vu3Trt2rVLPj4+uv/++zVu3Dg1a9ZMkjRy5Ej1799fU6ZMUadOnZxeNAAAwM1yOPy0bNlSGRkZatq0qV566SU98MAD8vUteP+R0NBQ7d+/3ylFAgAAOIvD4adnz5565JFHVLt27d89r0+fPnr66advujAAAICi4PCcnxEjRiglJUWzZs2yte3du1fPPPOMdu/ebWsrX768PDw8nFMlUEzWfXZY89fu1vy1u/X5Nz+7uhwAQBFwOPxs3rxZvXv31ldffWVr8/T01MmTJ9WzZ0/t2LHDqQUCxSl3ldhPSReUklY27mbOTu8AYM/h8DN79mw9+OCDeuutt2xt9evX19q1a9WpUyfNmDHDqQUCuDW5G6Tm9mit++ywq0sCAJdyeM7P0aNHFRcXV+ixBx98UAMHDrzlogA4143c9+hmdpAHgNLI4Z4fPz8/HT16tNBjiYmJKl++/C0XBaD4lcUhPwAojMPh5/7779fMmTP12Wef2bVv2bJFCQkJuu+++5xVGwAAgNM5POz17LPPavfu3RowYID+8Ic/qHLlyjp//ryysrLUpEkTPf/880VRJwAAgFM4HH58fX21YsUKbdmyRTt37lRqaqoqVqyo8PBwtW3bVu7ut7RjBgAAQJG6qY1N3dzc1LZtW7Vt29bJ5QDFJ+8EX4lJvgBgFTcVfrZu3arNmzfrypUrysnJsTvm5uamyZMnO6U4oCjlXwEV6F9wmxYAQNnjcPhZtGiRpk2bJm9vb/n7+8vNzf4Gavm/BwAAKEkcDj9vvfWWOnfurEmTJsnLy6soagIAACgyDs9OTklJ0SOPPELwAQAApZLD4SckJEQ//vhjUdQCFAn2trLH9QBgdQ4Pe40ePVpDhw6Vr6+vmjRponLlyhU453/+53+cUhzgDLl7W7Gy6zdcD6B48QdHyeNw+OnevbtycnI0evToa05uPnDgwE0VM3fuXG3btk3/+Mc/bG2jRo3S2rVr7c4LDAzU559/LknKycnR7Nmz9c477ygtLU3NmzfXuHHjVLNmzZuqAWUTK7vscT2A4lPYHxz/U9VXD0bVdWFV1uZw+Jk4cWJR1KFly5YpISFBERERdu0//PCDBgwYoJiYGFubh4eH7d9z587VypUr9corrygwMFBTp05Vv379tH79euYlAQBKhML+4MgbiIL8ffVwW8JQcXE4/HTp0sWpBZw+fVpjxozRrl27VKtWLbtj2dnZOnz4sAYOHKiAgIACj83MzNSSJUsUFxenqKgoSVJ8fLxat26tjRs3qmPHjk6tFQAAZ8kfiFB8bmoviszMTK1YsULPPPOMHnvsMR05ckRvv/22du/e7fBz7du3T5UqVdJ7772nJk2a2B07fvy4MjIyVKdOnUIfe/DgQV26dEmRkZG2Nj8/P4WEhGjHjh0O1wIAAMo+h3t+zp49qyeeeEJHjx5V7dq1dfjwYaWnp2vLli2aMmWKli1bprCwsBt+vujoaEVHRxd67NChQ3Jzc9Py5cv1+eefy93dXVFRURo6dKgqVqyopKQkSVL16tXtHletWjWdOnXK0bcGwEH55zHQdQ+gNHC45+e1117TpUuXtGHDBr377rsyxkiSZs6cqdDQUCUkJDituB9//FHu7u6qUaOG5s+fr5EjR2rLli0aOHCgcnJydOXKFUkqMLfH29tbGRkZTqsDQOFyu+1zv/IGIQAoqRzu+dm8ebNGjx6tmjVrKjs729bu7e2tvn376oUXXnBacYMHD1bv3r3l5+cnSQoODlZAQIAee+wx7dmzRz4+PpJ+G4bL/bckZWRkFLoEHwAAwOGen4yMDFWuXLnQYx4eHrp69eqt1mTj5uZmCz65goODJUlJSUm24a7k5GS7c5KTkxUUFOS0OgAAQNnhcPgJDQ3VihUrCj32/vvvq1GjRrdcVK5hw4YpNjbWrm3Pnj2SpLp166p+/fqqUKGCtm/fbjuelpam/fv3Kzw83Gl1AACAssPhYa9nn31WvXv31kMPPaSoqCi5ublp/fr1mjVrlr744gstWrTIacV16tRJTz/9tObNm6eOHTvq2LFjeumll9SpUyfbCrCYmBhNmzZN/v7+qlGjhqZOnaqgoCC1b9/eaXUAAICyw+HwEx4erqVLl2r69OlatGiRjDFatmyZQkJCtGDBArtl57eqXbt2mjlzpubPn6/58+erYsWK6ty5s4YOHWo7Z8iQIcrKytLYsWOVnp6uiIgILV68mBscAgCAQjkcfiQpIiJCK1euVHp6ulJTU1WhQgWVL1/+louZMmVKgbYOHTqoQ4cO13yMh4eH4uLiFBcXd8uvDwAAyr6bCj+5fHx87FZZASi72JwRQFnhcPipX7/+NTc0zXWzG5sCKLnyb87ITvAASiuHw8+gQYMKhJ9Lly7pm2++0U8//aThw4c7rTgAJUvevYjYCR5AaeVw+Bk8ePA1j40cOVJ79+5Vt27dbqkoAACAonJTG5tey8MPP6wNGzY48ykBAACcyqnh5/jx48rKynLmUwIAADiVw8Nes2fPLtCWk5OjU6dOacOGDdfcoR0AAKAkcEr4kaQKFSqoffv2GjVq1C0XBQAAUFQcDj8HDx4sijoAAACKhVPn/AAAAJR0Dvf8ODKs5ebmpsmTJzv6EgAAAEXG4fCTlJSk/fv3KzU1VTVq1FBgYKDOnz+vxMREGWMUFBRkO/d6d4IGAAAobg6HnwceeEA//vijVqxYoWbNmtnajx49qqefflo9evTQE0884dQiAQAAnMXhOT/z5s3T8OHD7YKPJNWuXVtDhw7V4sWLnVYcAACAszkcfs6dOyc/P79Cj7m5uenChQu3XBQAAEBRcTj8NGnSRHPmzNG5c+fs2k+fPq2EhAS1atXKacUBAAA4m8Nzfl544QXFxMQoOjpaYWFhqlKlin799Vd99913qlq1qkaPHl0UdQLXtO6zw0o6e9n2fcjt/i6sBgBQ0jnc81O/fn198MEHevzxx3Xx4kXt3btXGRkZ6tu3r9auXavq1asXRZ3ANSWdvayfki7YvlLS0l1dEgCgBHO450eSAgMDNXLkSGfXAgAAUORuKvxkZmZq9erV+vLLL3XmzBlNnjxZX3/9tRo2bKjGjRs7u0YApYS7O/f2AlDyORx+zp49qyeeeEJHjx5V7dq1dfjwYaWnp2vLli2aMmWKli1bprCwsKKoFUAJF1C5nN0crCB/Xz3ctq6LqwIAew7P+Xnttdd06dIlbdiwQe+++66MMZKkmTNnKjQ0VAkJCU4vEkDpkXcOVt6J6ABQUjgcfjZv3qxnn31WNWvWtNu+wtvbW3379tW+ffucWiAAAIAzORx+MjIyVLly5UKPeXh46OrVq7daEwAAQJFxOPyEhoZqxYoVhR57//331ahRo1suCgAAoKg4POH52WefVe/evfXQQw8pKipKbm5uWr9+vWbNmqUvvvhCixYtKoo6AZQhTIoG4EoO9/yEh4dr6dKlKleunBYtWiRjjJYtW6YzZ85owYIFioyMLIo6AZQhTIoG4EoO9/x8+eWXatq0qVauXKn09HSlpqaqQoUKKl++fFHUBwAA4FQO9/yMGDFCn3zyiSTJx8dHgYGBBB8AAFBqOBx+vLy85O3tXRS1AAAAFDmHh7369++vF198UQcPHtQdd9yhP/7xjwXOiYiIcEpxAAAAznZD4ScjI8PW2zNu3DhJ0ty5cyXJ7kaHxhi5ubnpwIEDzq4TAADAKW4o/ERHR2v27NkKCwtTRESEHn30UQUFBRV1bYDdkmiJZdEAgFt3Q+HnwoULSk5OliTt3LlTcXFx7N6OYpG7JBoAAGe5ofDTuHFjDRs2TK+++qqMMRo0aJC8vLwKPdfNzU2bNm1yapFALnd3t+ufBADA77ih8DN9+nQtW7ZM58+f17vvvquQkBD5+/sXdW1AAQGVy9kNhYXczs8hAMAxNxR+AgMDNXLkSEnS9u3b9dxzz6l+/fpFWhhwLXmHwgL9fV1cDQCgtHF4qfunn35aFHUAAAAUC4fDDwAAxSX/ik+GuuEMhB8AQImVf8UnQ91wBoe3twAAACjN6PkBAJQYrOZEcSD8AHAp7t2EvFjNieJA+AHgUvnv3SSxjQmAokX4AVBkbrRXh21MABQnwg+AIlNYrw7zOAC4GuEHQJFiqTKAkoal7gAAwFIIPwAAwFIIPwAAwFIIPwAAwFIIPwAAwFJKVPiZO3euevXqZdd24MABxcTEqGnTpmrbtq0WL15sdzwnJ0cJCQlq3bq1mjRpor59+yoxMbE4ywYAAKVIiQk/y5YtU0JCgl3buXPn1KdPH91+++1as2aNBg8erJkzZ2rNmjW2c+bOnauVK1dq4sSJWrVqldzc3NSvXz9lZmYW91sAAAClgMvv83P69GmNGTNGu3btUq1ateyO/fOf/5SXl5fGjx8vT09P1alTR4mJiVq4cKG6deumzMxMLVmyRHFxcYqKipIkxcfHq3Xr1tq4caM6duzoircE4Bax3xeAouTy8LNv3z5VqlRJ7733nubMmaNffvnFdmznzp2KiIiQp+d/y4yMjNSCBQuUkpKiX375RZcuXVJkZKTtuJ+fn0JCQrRjxw7CD1BK5b8z9LX2+rqRc4DSjJ/xouHy8BMdHa3o6OhCjyUlJSk4ONiurVq1apKkkydPKikpSZJUvXr1AuecOnWqCKoFUFxuZL8v9gRDWcfPeNEoMXN+CpOeni4vLy+7Nm9vb0lSRkaGrly5IkmFnpORkVE8RQIAgFKlRIcfHx+fAhOXc0ONr6+vfHx8JKnQc8qVK1c8RQIAgFKlRIefoKAgJScn27Xlfh8YGGgb7irsnKCgoOIpEgAAlColOvxERERo165dys7OtrVt27ZNtWrVUtWqVVW/fn1VqFBB27dvtx1PS0vT/v37FR4e7oqSAQBACVeiw0+3bt108eJFjRkzRocPH9batWu1fPly9e/fX9Jvc31iYmI0bdo0ffLJJzp48KCee+45BQUFqX379i6uHgAAlEQuX+31e6pWrapFixZp0qRJ6tKliwICAjRixAh16dLFds6QIUOUlZWlsWPHKj09XREREVq8eHGBSdAoHfIu6wy53d/F1QAAyqISFX6mTJlSoK1x48ZatWrVNR/j4eGhuLg4xcXFFWVpKCZ5l3UG+vu6uBoAQFlUooe9AAAAnI3wA6DEY7sLAM5Uooa9AKAw+be7kJgTBuDmEX4AlAr5b/PPnDAAN4thLwAAXIyh3eJFzw9chmEMFIf8t084m5Zu93PHTtkoCRjaLV6EH7gMwxgoDvlvn3CaXbJRQvE7sfgw7AXA0hhuAKyHnh8AZcLNhpj8ww0MgwFlH+EHQJlwK3Mm8g83ACjbCD8AygzmTAC4Ecz5AQAAlkL4AQAAlkL4AQAAlsKcHxSJ/BNPWUEDID9u6gdXIfygSLB6BsD1MEEdrsKwFwAAsBTCDwAAsBTCDwAAsBTCDwAAsBTCDwAAsBTCDwAAsBTCDwAAsBTCDwAAsBTCDwDk4e7u5uoSABQx7vAMAHkEVC53Q9uzsIULUHoRfgAgn/zbLhTWG8QWLkDpRfgBgOvI3xvEBpxA6Ub4gcPo7ocV5e3pYQNOoHQj/MBhdPcDAEozVnuhWLCCBgBQUtDzg2LBnAkAQElB+EGxYc4EYG38AYSSgvADACgW/AGEkoI5PwAAwFIIPwAAwFIIPwAAwFIIPwAAwFIIPwAAwFJY7QUAcLr82+CwtB0lCeEHAOB0+bfBYWk7ShKGvQAAgKUQfgAAgKUQfgAAgKUQfgDACdzd3VxdAoAbxIRn3DJ+6QNSQOVydiucgvx99XDbui6uCkBhCD+4Zfl/6bOkFVaVf4UTgJKJ8AOnYLdmAEBpwZwfAABgKYQfAABgKYQfAABgKaVizs8vv/yi6OjoAu0TJ07Uo48+qgMHDmjSpEnau3evKleurF69eik2NtYFlQLA72NFGOB6pSL8/PDDD/L29tamTZvk5vbfZdUVK1bUuXPn1KdPH917772aMGGCvvvuO02YMEGVK1dWt27dXFg1ABTEijDA9UpF+Dl06JBq1aqlatWqFTi2fPlyeXl5afz48fL09FSdOnWUmJiohQsXEn4AAEABpWLOzw8//KC6dQvvGt65c6ciIiLk6fnfHBcZGaljx44pJSWluEoEAAClRKkIP4cOHVJKSop69Oihu+++W927d9d//vMfSVJSUpKCgoLszs/tITp58mSx1woAAEq2Ej/slZmZqePHj6tcuXIaMWKEfH199d5776lfv35aunSp0tPT5eXlZfcYb29vSVJGRoYrSgYAACVYiQ8/Xl5e2rFjhzw9PW0hp1GjRjpy5IgWL14sHx8fZWZm2j0mN/T4+nKnYQAAYK9UDHv5+voW6N0JDg7W6dOnFRQUpOTkZLtjud8HBgYWW40AABQlNpF2nhLf83Pw4EF1795dCxcuVHh4uK197969qlu3rho0aKCVK1cqOztbHh4ekqRt27apVq1aqlq1qqvKBgDAqfJvIi1xr6ibVeJ7foKDg3XHHXdowoQJ2rlzp44cOaJXXnlF3333nQYMGKBu3brp4sWLGjNmjA4fPqy1a9dq+fLl6t+/v6tLBwDAqXLvE5X7lTcI4caV+J4fd3d3zZ8/X9OmTdPQoUOVlpamkJAQLV26VPXq1ZMkLVq0SJMmTVKXLl0UEBCgESNGqEuXLi6uHAAAlEQlPvxIkr+/vyZPnnzN440bN9aqVauKsSIAAFBalfhhLwAojZicCpRcpaLnBwBKm8Imp4bc7u/CigDkIvwAQBHJv4lpoH/ZvfdY3qBHyENJR/gBANyyvEGvLIe8kobh1ZtD+AGAUihvTwv3erGu/MOr/CzcGMIProvubKBo3Ohf7YXNHco/pAbr4mfBcYQfXBfd2UDRuNak6LNp6XZ/cFhp7hBQHAg/AOBChQWb0/zBARQp7vMDAKUck14Bx9DzAwClHBteAo4h/ABAGcCkV+DGEX4AwCLoHQJ+Q/gBAIugdwj4DeEHAOAQ9ixDaUf4AQA4hPsOobQj/AAAbJgXBCsg/AAAbJgXBCvgJocAgGviBoooi+j5AQCLupFgk/8GikxuRllA+AEAi7rRYMPmxihrCD8AYGEEG1gRc34AAICl0PMDO9y8DABQ1hF+YIeblwEAyjqGvQAAgKUQfgAAgKUQfgAAgKUQfgAAKCO4I/eNYcIzAJRBfAhaU/4bV0psTlsYwg8AlEFsS2FdbE57fYQfACijuHszUDjm/AAAAEuh58dCGAcGAIDwYymMAwMAwLAXAACwGMIPAACwFMIPAACwFMKPhXETNAAo+/hdXxATni2Mm6ABQNmX/3c9K30JP5bHTdAAoOxjta89hr0AAIClEH4AAIClEH4AALAQJkAz5wcAAEvJPwFast4kaMIPAAAWY/UJ0Ax7AQAASyH8AAAASyH8AAAASyH8AAAASyH8AABgcVZb/s5qLwAALM5q+38RfgAAgKWWv5eJ8JOTk6PZs2frnXfeUVpampo3b65x48apZs2ari7NpdixHQCAgspE+Jk7d65WrlypV155RYGBgZo6dar69eun9evXy8vLy9XlFYv8d+sMud2fHdsBAE51M0NjJfFu0qU+/GRmZmrJkiWKi4tTVFSUJCk+Pl6tW7fWxo0b1bFjRxdXWDzyd1cSdgAAN+taE6BvZmisJA6nlfrwc/DgQV26dEmRkZG2Nj8/P4WEhGjHjh1lNvwwpAUAKCqF7f+V/7PmWgGpNHw+uRljjKuLuBX//ve/NXjwYH3//ffy8fGxtT/77LNKT0/XggULHHq+0NBQZWdnq3r16s4u1akup19VVvZv/+m8/uChnJwc2/eFtXHO75/j6tfnHM7hHP7/Lo3nuEnKzhMjPD3clZWV/bvP4+nhJl+fP8jZTp06JQ8PD+3Zs+e655b6np8rV65IUoG5Pd7e3kpNTXX4+by9vZWZmemU2opSwR8cj0LOyt/GOSX79TmHcziH/79L2zkFeXnmv4Xg9R/jDJ6enjc8z7fUh5/c3p7MzEy7np+MjAyVK1fO4efbuXOn02oDAAAlT6m/w3Pu8FRycrJde3JysoKCglxREgAAKMFKffipX7++KlSooO3bt9va0tLStH//foWHh7uwMgAAUBKV+mEvLy8vxcTEaNq0afL391eNGjU0depUBQUFqX379q4uDwAAlDClPvxI0pAhQ5SVlaWxY8cqPT1dERERWrx4sWVucAgAAG5cqV/qDgAA4IhSP+cHAADAEYQfAABgKYQfAABgKYQfAABgKYQfAABgKYQfAABgKYQfAABgKYSfYpCTk6OEhAS1bt1aTZo0Ud++fZWYmOjqskq98+fP68UXX1SbNm3UrFkzde/e3W5j2gMHDigmJkZNmzZV27ZttXjxYhdWWzYcO3ZMYWFhWrt2ra2N6+w869at0wMPPKDQ0FB17NhRH374oe0Y19l5rl69qvj4eLVt21ZhYWHq0aOHvvnmG9txrvWtmzt3rnr16mXXdr3rWqyflQZFbtasWeauu+4yn332mTlw4IDp27evad++vcnIyHB1aaVanz59zIMPPmh27Nhhjhw5Yl5++WXTuHFjc/jwYXP27Flz5513mjFjxpjDhw+b1atXm9DQULN69WpXl11qZWZmmq5du5rg4GCzZs0aY4zhOjvRunXrTIMGDcyyZcvM8ePHzezZs039+vXNN998w3V2spkzZ5qWLVua//znP+b48eNmzJgxplmzZiYpKYlr7QRLly419erVMzExMba2G7muxflZSfgpYhkZGSYsLMysWLHC1paammoaN25s1q9f78LKSrfjx4+b4OBgs2vXLltbTk6Oad++vXn99dfN/PnzTevWrc3Vq1dtx6dPn246dOjginLLhOnTp5tevXrZhR+us3Pk5OSYdu3amSlTpti19+3b18yfP5/r7GQPPvigeeWVV2zfX7hwwQQHB5uPPvqIa30LkpKSTGxsrGnatKm5//777cLP9a5rcX9WMuxVxA4ePKhLly4pMjLS1ubn56eQkBDt2LHDhZWVblWqVNEbb7yhRo0a2drc3NxkjFFqaqp27typiIgIeXr+d/u6yMhIHTt2TCkpKa4ouVTbsWOHVq1apVdffdWunevsHEePHtUvv/yizp0727UvXrxY/fv35zo7WeXKlbV582b9/PPPys7O1qpVq+Tl5aUGDRpwrW/Bvn37VKlSJb333ntq0qSJ3bHrXdfi/qwk/BSxpKQkSVL16tXt2qtVq6ZTp065oqQywc/PT1FRUXab13744Yf66aef1KpVKyUlJSkoKMjuMdWqVZMknTx5slhrLe3S0tI0YsQIjR07tsDPMdfZOY4fPy5Junz5smJjY3XXXXfp0Ucf1aeffiqJ6+xsY8aMkaenp+655x6FhoYqPj5er7/+um677Tau9S2Ijo7W9OnT9ec//7nAsetd1+L+rCT8FLErV65IUoEd5r29vZWRkeGKksqkXbt2afTo0brnnnsUHR2t9PT0Qq+5JK67g8aPH6+mTZsW6JWQxHV2kosXL0qSRo4cqU6dOmnJkiVq2bKlBg4cqG3btnGdnezIkSPy8/PTnDlztGrVKnXt2lUjR47UwYMHudZF5HrXtbg/Kz2vfwpuhY+PjyQpMzPT9m/pt//Y5cqVc1VZZcqmTZs0fPhwNWnSRDNmzJD023XPzMy0Oy/3fyBfX99ir7G0WrdunXbu3Kn333+/0ONcZ+f4wx/+IEmKjY1Vly5dJEkNGjTQ/v37tXTpUq6zE/3yyy+Ki4vTsmXLFB4eLkkKDQ3V4cOHNWvWLK51EbnedS3uz0p6fopYbhdecnKyXXtycnKBLkA47s0339TgwYPVpk0bLVy40PY/TVBQUKHXXJICAwOLvc7Sas2aNUpJSbEtCQ4LC5MkjRs3Th07duQ6O0nu74Lg4GC79rp16+rnn3/mOjvR7t27dfXqVYWGhtq1N2nSRMePH+daF5HrXdfi/qwk/BSx+vXrq0KFCtq+fbutLS0tTfv377f91YGbs2LFCr388svq2bOnXn/9dbvu0oiICO3atUvZ2dm2tm3btqlWrVqqWrWqK8otlaZNm6YNGzZo3bp1ti9JGjJkiN544w2us5OEhISofPny+v777+3aDx06pNtuu43r7ES5H7I//PCDXfuhQ4dUs2ZNrnURud51LfbPSqevH0MBM2bMMC1atDCbNm2y3bvgvvvu4z4/t+Do0aOmYcOGZtCgQSY5OdnuKy0tzfz6668mIiLCjBw50vz4449mzZo1JjQ01Kxdu9bVpZd6eZe6c52dZ86cOSYsLMy8//77JjEx0cydO9fUr1/ffPXVV1xnJ8rOzjY9evQw999/v9m2bZs5duyYiY+PNw0aNDDffvst19pJRo4cabfU/Uaua3F+VhJ+ikFWVpZ57bXXTGRkpGnatKnp16+fOXHihKvLKtXmzZtngoODC/0aOXKkMcaY77//3vz1r381jRo1Mu3atTP/+Mc/XFx12ZA3/BjDdXamJUuWmOjoaNOwYUPz4IMPmo0bN9qOcZ2d5/z582b8+PGmbdu2JiwszDz22GNm+/bttuNc61uXP/wYc/3rWpyflW7GGOP8/iQAAICSiTk/AADAUgg/AADAUgg/AADAUgg/AADAUgg/AADAUgg/AADAUgg/AADAUtjYFIBlrV27VqNGjbJr8/b2Vo0aNfTQQw+pX79+8vDwKJZa6tWrZ/e9h4eHKlasqJCQEPXu3VtRUVHFUgdgBYQfAJY3e/ZsBQQEyBijK1eu6JtvvlFCQoLS09M1dOjQYqvjkUce0aOPPipJunr1qs6cOaPVq1frqaee0t///nfFxMQUWy1AWUb4AWB5DRo00J/+9Cfb93fffbdOnDihlStXFmv4CQoKUtOmTe3a/vKXv2jQoEGaMmWK2rZta1cngJvDnB/AQvbt26cnnnhCzZs3V1hYmHr37l1gJ/GtW7eqZ8+eCgsLU6tWrfTiiy8qNTXVdvz48eMaMmSIWrZsqaZNm6pXr17atWuX7fjPP/+sevXqaenSpfrLX/6iFi1aaO3atZJ+2zm7f//+atasmZo1a6ZBgwbpxIkT161769at6tGjh5o3b64777xTw4YN06lTp2zH165dq5CQEH3//fd67LHHFBoaqrZt22rhwoU3fa0qVaokNze3Aq/xzjvvqFWrVmrTpo1+/PFHSdKmTZvUtWtXhYaGqmXLlpo4caIuX74sSbp48aKio6N1//33KzMzU5JkjFHfvn1111136ddff/3dOtzc3DRs2DBdvXpVq1evtrVnZGTotddeU1RUlBo1aqTOnTtrw4YNtuOffPKJ6tWrp1mzZtnajh07pqZNm2rkyJE3fV2AsoDwA1jExYsX9eSTT6pKlSpKSEhQfHy8rly5otjYWF24cEGStGXLFj355JOqXLmy4uPjFRcXp08//VRDhgyRJB0+fFhdu3bViRMnNHbsWE2bNk1ubm564okn9PXXX9u9Xnx8vGJjYzVx4kRFRkbq2LFjevzxx5WSkqIpU6Zo0qRJOnHihLp3766UlJRr1v2vf/1Lffv2VWBgoGbMmKFRo0bp22+/1WOPPWb3uJycHA0dOlQPPPCA3njjDTVv3lzTpk3Tf/7zn+tem5ycHGVlZSkrK0sXL17U559/rn/961/q2bOn3XnZ2dmaP3++Jk6cqKFDh6pu3bp6//33NWjQINWuXVtz5szRM888o/fee08DBw6UMUYVKlTQpEmTdPz4cc2fP1+StGLFCm3dulWTJk3SH//4x+vWV6dOHVWvXt0WMo0xGjRokFauXKk+ffpo3rx5CgsL03PPPad169ZJku655x49/PDDWrBggY4cOaLs7Gy98MIL8vf319///vfrviZQphXJdqkASpxvv/3WBAcHm507d9raEhMTzauvvmpOnjxpjDGma9eu5uGHH7Z73EcffWTuu+8+k5SUZJ599lnTokULk5aWZjt+9epV06FDB/PII48YY4w5ceKECQ4ONsOGDbN7nueff97cdddd5sKFC7a2c+fOmebNm5spU6YUWnN2drZp2bKl6d27t117YmKiadiwoXnttdeMMcasWbPGBAcHm3/+85+2czIyMkxoaKh56aWXrnlNch9X2Fe3bt3s3mdhr5GTk2PatGljYmNj7Z73yy+/NMHBwWbz5s22tgkTJpiGDRuaLVu2mKZNm5q///3vdo8JDg42CQkJ16y1W7du5v777zfGGPPFF1+Y4OBg88EHH9idM3z4cNOyZUtz9epVY4wxqamppnXr1iYmJsYsWLDANGjQwOzYseOarwFYBT0/gEXccccd8vf319NPP61x48bp008/VUBAgEaMGKHq1asrPT1d+/bt07333mv3uA4dOujjjz9WYGCgvv76a7Vr104VK1a0Hff09FTHjh21Z88eXbp0ydYeHBxs9zxfffWV7rzzTvn4+Nh6WSpUqKDw8HB9+eWXhdZ87NgxnTlzRp07d7Zrv+222xQWFqbt27fbtYeFhdn+7eXlJX9/f9vw0++ZN2+eVq9erdWrV+utt97S+PHjlZycrMcff1wXL160Ozfv+zp69KiSkpIUHR1te09ZWVmKiIhQhQoVtHXrVtu5w4cPV/Xq1dW/f39Vq1atwCqzG5E7DLdt2za5ubkpKirK7nWjo6N15swZ23Ccn5+fJk6cqK+//lrx8fHq16+fwsPDHX5doKxhwjNgEeXLl9dbb72lefPmacOGDVq5cqXKlSunBx98UGPGjFFqaqqMMapateo1nyM1NbXQYZo//vGPMsbYBYX8550/f14bNmywm5eSy9/fv9DXO3/+fKHPldu2f/9+uzYfHx+7793d3WWMKfzN5BEcHGw3kTg8PFzBwcHq0aOH3nnnHfXp08d2LO/1ya1vwoQJmjBhQoHnTU5Otv3b19dXHTp00MKFCxUZGaly5cpdt668Tp8+rTvuuMP2usYYNWvWrNBzk5OT1aBBA0nSXXfdperVq+vUqVOKjo526DWBsorwA1hI7dq1NXXqVGVnZ2v37t3617/+pbffflt/+tOf1LNnT7m5uens2bN2j8nMzNS2bdvUuHFjVapUqdAJumfOnJEkValSxe4DP6+KFSvq7rvvtgsSuTw9C/9VVLlyZUm65mtWqVLld9/vrcgND8ePH7/mOX5+fpKkESNGqEWLFgWOV6pUyfbvw4cPa/ny5WrQoIH++c9/qnPnzjfcC3PkyBElJyerR48ekn67lr6+vvrf//3fQs+vWbOm7d9z5szRr7/+qjp16mjs2LFas2aNvLy8buh1gbKKYS/AIj766CNFRkbqzJkz8vDwUFhYmMaPHy8/Pz8lJSWpfPnyatCggT755BO7x33xxRd66qmnlJSUpIiICG3evNk2QVr6bRLwBx98oNDQ0N/9UG3RooUOHz6sBg0aKDQ0VKGhoWrUqJGWLVumjRs3FvqYWrVqKSAgQO+//75d+4kTJ/Tdd99ds+fDGb777jtJ0u23337Nc2rXrq2qVavq559/tr2n0NBQBQUFafr06baeqaysLI0cOVI1atTQ22+/rUaNGmnUqFE3NCQnSQkJCfLx8VGXLl0k/XYtL1++LGOM3ev++OOPmjNnjrKysiRJe/bs0cKFCzVgwABNnz5dR48etVv9BVgVPT+ARTRr1kw5OTkaNGiQnnrqKZUvX14ffvihLly4oPvuu0+SNGTIED399NMaOnSounbtqrNnz2r69Olq166dGjRooGeeeUaff/65/va3v+mpp56Sl5eX3nzzTZ04cUKLFi363dcfOHCgHn/8cfXv31/du3eXt7e3Vq1apU2bNikhIaHQx7i7u+v555/XqFGj9Nxzz+nhhx/WuXPnNHv2bFWqVKnQXqSbceDAAVvvkjFGR44cUUJCggICAmyBozAeHh567rnn9OKLL8rDw0Pt2rVTWlqa5s6dq9OnT6thw4aSpAULFmjfvn168803Va5cOb388svq1q2bpk2bphdffNH2fElJSbbQlZWVpdOnT+vdd9/VF198oZdeeklBQUGSpKioKEVERGjgwIEaOHCg6tSpo927d2vWrFlq1aqV/P39lZmZqRdeeEG1atWy/bf629/+psWLF+vee+9VkyZNnHLtgFLJlbOtARSv77//3vTt29e0aNHChIaGmq5du5p///vfdud89tlnplu3bqZRo0amdevWZtKkSebixYu24/v37zdPPvmkadq0qQkLCzNPPPGE3Qqi3NVea9asKfD6e/fuNbGxsSYsLMw0bdrU/PWvfzWbNm26bt0fffSR6dKli2nYsKG58847zfDhw20r1Iz570qsEydO2D2uXbt2ZuTIkdd83sJWe4WEhJiWLVua559/3iQmJl73NYwx5oMPPjBdunQxjRo1Mi1atDADBgwwBw8eNMYYc+DAAdOwYUMzbtw4u8e89tprpl69eubLL780xpgCdTRs2NC0a9fODB48uNAVWpcuXTKTJ082bdq0MQ0bNjTR0dFm+vTpJj093RhjzKuvvmrq169vvv32W9tjLl++bKKjo02HDh1s5wFW5GbMDcwGBAAAKCOY8wMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACyF8AMAACzl/wCDwfeW3+SgPgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.histplot(UKBrexdex.score, bins=range(101))\n", "plt.xlabel('score on BrexDex')\n", "plt.ylabel('frequency')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "aadea440", "metadata": {}, "source": [ "How would you describe this data?\n", "* write some text, including descritive statatistics, to describe the distribution of Brexdex scores in the UK population" ] }, { "cell_type": "markdown", "id": "8f3bca1d-f652-46f8-be41-51057d6be0ba", "metadata": {}, "source": [ "## Drawing a sample - recap\n", "\n", "We can draw a sample from our dataframe using `df.sample()`. The output of `df.sample()` is another dataframe (with $n$ rows, where $n$ is the sample size)\n", "\n", "If we like, we can save this dataframe with a new name, like `mySample`\n", "\n", "* Create a new dataframe called `mySample` containing 10 people randomly sampled from the Brexdex dataset\n", "* Get the mean score for this sample\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "6240f953-9b53-484a-a868-1c0ea279242e", "metadata": {}, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "markdown", "id": "756e9c2e-8cf0-46c2-946d-2329c2f62906", "metadata": {}, "source": [ "## Sampling distribution of the mean\n", "\n", "In general we are not interested in any specific sample, but just in working out how much random variation to expect between samples. To quantify this random variation we will need to generate a lot of samples (in a loop) but we don't need to save each sample as a new dataframe, we can just save the sample mean from each sample into an array.\n", "\n", "Let's do this in steps:\n", "\n", "First, get a sample of 10 people from the Brexdex distribution, find their mean score, and save it as a variable `m` - this should be a single line of code and if we print m, we should get a single number (somewhere near 50)" ] }, { "cell_type": "code", "execution_count": 6, "id": "264b4f92-bb49-425d-a2f1-24a211d9b969", "metadata": {}, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "markdown", "id": "bf5d1d79-f06b-4441-ade4-219a11a0e8d2", "metadata": {}, "source": [ "Now let's make a loop to do this 50 times and save the mean of each sample into an array called `m`. When we print out `m`, we should now get an array of 50 numbers." ] }, { "cell_type": "code", "execution_count": 8, "id": "856ab1fa-be87-4550-ac02-abc0e9aa8bed", "metadata": {}, "outputs": [], "source": [ "# I've included part of the code to get you started\n", "# You will need to uncomment it\n", "\n", "#m = np.zeros(50)\n", "#for i in range(len(m)):\n", " # Your code here\n", "\n", "#m" ] }, { "cell_type": "markdown", "id": "b693c02c-9bdc-43af-b65a-0fa4825c3ada", "metadata": {}, "source": [ "Change the code above to get\n", "* 100 samples of size 10\n", "* 10 samples of size 100" ] }, { "cell_type": "markdown", "id": "dfe366ba-f71a-4e76-9ae6-27a91177faea", "metadata": {}, "source": [ "Finally, let's get the sample means for 10000 samples of size 50 and plot them on a histogram" ] }, { "cell_type": "code", "execution_count": 9, "id": "d74f7f65-3d83-49b1-8115-9575cf6f5e8c", "metadata": {}, "outputs": [], "source": [ "# Your code here\n", "\n", "#sns.histplot(m)\n", "#plt.xlabel('sample mean')\n", "#plt.show()" ] }, { "cell_type": "markdown", "id": "1e2e2dd7-c779-468d-ac56-7808bb328476", "metadata": {}, "source": [ "## Standard error of the mean\n", "\n", "The standard error of the mean is the standard deviation of the sampling distribution of the mean.\n", "\n", "Let's work out the SEM from the distribution of sample means we just generated:" ] }, { "cell_type": "code", "execution_count": 31, "id": "5faeec75-1ea8-4b13-a6c7-93734518134a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.521772474810944" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# m.std()" ] }, { "cell_type": "markdown", "id": "a45efb4d-17d8-42c3-a6d4-0093d2efa101", "metadata": {}, "source": [ "The SEM can also be calculated as $$ SEM = \\frac{\\sigma}{\\sqrt{n}} $$" ] }, { "cell_type": "code", "execution_count": 10, "id": "6728bb36-c299-4591-b69d-27ae728da0e1", "metadata": {}, "outputs": [], "source": [ "# calculate the standard deviation of the Brexdex data\n", "# sigma = \n", "\n", "# n is the SAMPLE size\n", "n=50\n", "\n", "# Fill in the equation for the SEM\n", "# Your code here" ] }, { "cell_type": "markdown", "id": "558da4d2-6db4-4bd1-ba52-312aa3e68d47", "metadata": {}, "source": [ "Hopefully the two values match well." ] }, { "cell_type": "markdown", "id": "650efbe0", "metadata": {}, "source": [ "## Shape of the sampling distribution\n", "\n", "Here we explore the shape of the sampling distribution of the mean\n", "\n", "### Sampling distribution of the mean when $n$ is small\n", "\n", "Let's start by simulating what the sampling distribution of the mean looks like for small values of $n$ - starting with $n=1$\n", "\n", "#### $n=1$\n", "\n", "Write some code to draw 10,000 samples with n=1, obtain the mean of each sample, and plot the means.\n", "\n", "**Think**\n", "* Before you run it, think: what will this 'sampling distribution of the mean' look like?\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "8f8ef0c3", "metadata": {}, "outputs": [], "source": [ "# Your code here!\n", "# Note this can be copied from examples earlier in the chapter\n", "nSamples = 10000 # we will draw 10,000 samples\n", "samplesize=1 # each sample contains n people\n", "\n", "# make an empty array 'm' to store the means\n", "\n", "# Make a loop to get the samples and store their means in m\n", "\n", "# sns.histplot(m)\n", "# plt.xlabel('sample mean')\n", "# plt.show()" ] }, { "cell_type": "markdown", "id": "7c7549dd", "metadata": {}, "source": [ "In fact when $n=1$, the sample mean is simply the value of the (one) person in the sample's score, so the sampling distribution of the mean is exactly the sample data distribution\n", "\n", "#### $n=2$\n", "\n", "Write some code to draw 10,000 samples with n=2, obtain the mean of each sample, and plot the means." ] }, { "cell_type": "code", "execution_count": 12, "id": "aa81f9b2", "metadata": {}, "outputs": [], "source": [ "# Your code here!\n", "# Should be the same as the previous one but with a different n" ] }, { "cell_type": "markdown", "id": "fd102df1", "metadata": {}, "source": [ "Hopefully you have noticed a middle peak emerging.\n", "\n", "**Think**-\n", "\n", "A simple summary of the Brexdex distribution is that people are either pre-Brexit (belonging to the lower mode of the distribution, the hump of scores below 50%), or they are against Brexit (belonging to the upper mode).\n", "\n", "If we draw a sample of $n=2$, there are four possible outomes:\n", "\n", "* pro-pro\n", "* pro-against\n", "* against-pro\n", "* against-against\n", "\n", "Case 1 yields low scores, case 4 yields high scores, and cases 2 and 3 yield intermediate scores.\n", "\n", "How does this relate to the simulated sampling distribution of the mean you plotted?" ] }, { "cell_type": "markdown", "id": "f9302f17", "metadata": {}, "source": [ "#### $n = 3,4,5$\n", "\n", "As $n$ increases, we rapidly see a unimodal, bell-curve-like shape emerging\n", "\n", "Write some code to simulate the sampling distribution of the mean for 10,000 samples each of $n=3,4,5$ and plot a histogram of the sample means for each value of $n$. Organise these plots as subplots next to or above each other (you decide which is more informative)\n", "\n", "To do this, you will need to use a double `for`-loop. I have completed some parts of the code to get you started.\n", "\n", "To compare the sampling distributions effectively, students should plot them above eachother and match the scales on the x-axis" ] }, { "cell_type": "code", "execution_count": 15, "id": "9e3efd2b", "metadata": {}, "outputs": [], "source": [ "# Note this can be copied from examples in chapter 0 (\"prepwork\")\n", "nSamples = 10000 # we will draw 10,000 samples\n", "samplesize=[3,4,5] # each sample contains n people\n", "\n", "for j in range(len(samplesize)):\n", " m=np.empty(nSamples) # make an array to store the means\n", "\n", " for i in range(nSamples):\n", " pass # This is a placeholder, delete it!\n", " # Your code here to fill in m with sample means\n", " \n", " #plt.subplot(3,1,j+1)\n", " #sns.histplot(m, bins=range(0,100,5))\n", " #plt.xlabel('sample mean')\n", "\n", "#plt.tight_layout()\n", "#plt.show()" ] }, { "cell_type": "markdown", "id": "cf481a80-12da-4579-8134-f95261a61e02", "metadata": {}, "source": [ "## SEM is proportional to $\\frac{1}{\\sqrt{n}}$\n", "\n", "The standard error of the mean, which is the width of the sampling distribution of the mean, is inversely proportional to the square root of the sample size $n$\n", "\n", "Using the double loop, as above, we can plot the SEM calculated in two ways:\n", "* as the standard deviation of the simulated sampling distribution of the mean\n", "* From the equation\n", "\n", "... for a range of values of n\n", "\n", "Complete the code block below:" ] }, { "cell_type": "code", "execution_count": 16, "id": "51d81422-53e1-45b4-8d0e-4db056c93e72", "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (1398808157.py, line 13)", "output_type": "error", "traceback": [ "\u001b[0;36m Cell \u001b[0;32mIn[16], line 13\u001b[0;36m\u001b[0m\n\u001b[0;31m simulatedSEM[j] = # your code here\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "nSamples = 100 # we will draw 100 samples at each samplesize\n", "samplesizes=range(1,100) # try all values of n from 1 to 100\n", "simulatedSEM=np.zeros(len(samplesizes))\n", "calculatedSEM=np.zeros(len(samplesizes))\n", "\n", "for j in range(len(samplesizes)):\n", " m=np.zeros(nSamples) # make an array to store the means\n", "\n", " for i in range(nSamples):\n", " # Your code here to fill in m with sample means\n", " pass # This is a placeholder, delete it!\n", " \n", " simulatedSEM[j] = # your code here\n", " calculatedSEM[j] = # your code here\n", "\n", "#plt.plot(samplesizes, calculatedSEM, 'r')\n", "#plt.plot(samplesizes, simulatedSEM, 'k.-')\n", "\n", "#plt.xlabel('sample mean')\n", "#plt.ylabel('SEM')\n", "\n", "#plt.tight_layout()\n", "#plt.show()" ] }, { "cell_type": "markdown", "id": "9ba12ced", "metadata": {}, "source": [ "## When does the CLT apply?\n", "\n", "The Central Limit Theorem states that the sampling distribution of the mean is estimated by $\\mathcal{N}(\\bar{x},\\frac{s}{\\sqrt{n}})$ when $n$ is large\n", "\n", "But how large is large enough?\n", "\n", "A good rule of thumb is that the Central Limit Theorem applies for $n>50$, and a larger $n$ is required for a roughly normal sampling distribution when the data distribution is grossly non-normal (such as the bimodal Brexdex distribution). \n", "\n", "In reality, the normal distribution becomes a closer and closer fit tot the sampling distribution of the mean as $n$ gets larger\n", "\n", "### $n=100$\n", "\n", "Let's start with a value of $n$ for which the central limit theorem should definitely apply, $n=100$\n", "\n", "Now, we work out the mean and SEM that would be predicted for the sampling distribution of the mean, if the central limit theorem applied.\n", "\n", "Finally we compare the predicted normal distribution to the simulated sampling distribution of the mean in a plot\n", "\n", "**Note -**\n", "The code to get the normal curve and histogram to match in scale is a bit fiddly, this was covered in the tutorial exercises last week" ] }, { "cell_type": "code", "execution_count": 7, "id": "96ed1083", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAG1CAYAAAAWb5UUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABgx0lEQVR4nO3deXxU9f398ddMksm+74QtLAlbWGR1ARQXWsWFarVudUGrpdUvtVVrbSv+WhWrrXVpRC2KtbYqSBHXVkBRERAQkC1AAglLNrJP1klm7u+PmEhIgCQkuTOZ83w8InDvnZtzGe/wzud+FothGAYiIiIiXsJqdgARERGRnqTiR0RERLyKih8RERHxKip+RERExKuo+BERERGvouJHREREvIqKHxEREfEqKn5ERETEq/iaHcDdTJgwAYfDQWxsrNlRREREpJ2OHj2KzWZj06ZNpzxWxc9x6urqcDqdZscQERGRDmhoaKC9i1ao+DlOXFwcAKtWrTI5iYiIiLTX+eef3+5j1edHREREvIqKHxEREfEqKn5ERETEq6j4EREREa+i4kdERES8ioofERER8SoqfkRERMSrqPgRERERr6LiR0RERLyKih8RERHxKip+RERExKuo+BERERGvouJHREREvIqKHxEREfEqKn5ERETEq6j4ERFpp/Sl20hfus3sGCJymnzNDiAi4ins1Q6zI4hIF1DLj4iIiHgVFT8iIiLiVVT8iIiIiFdR8SMiIiJeRcWPiIiIeBUVPyIiIuJVVPyIiIiIV1HxIyIiIl5FxY+IiIh4FRU/IiIi4lVU/IiIiIhXUfEjIiIiXkXFj4iIiHgVFT8iIiLiVVT8iIiIiFdR8SMiIiJeRcWPiIiIeBVfswMcKz09nXXr1vHaa68BcOONN/LVV1+1eezjjz/OFVdcwZEjR5gxY0ar/X/84x/54Q9/2K15RURExPO4TfGzePFinnnmGSZOnNi87dlnn6W+vr7Fcb/97W85ePAgF1xwAQB79uzB39+flStXYrFYmo8LDQ3tmeAiIiLiUUwvfgoKCnjwwQfZvHkzycnJLfZFRES0+PN7773HF198wbJlywgJCQFg7969JCcnExcX11ORRURExIOZ3udn586dhIeHs2LFCsaMGXPC46qrq/nTn/7ETTfdRGpqavP2PXv2MGTIkJ6IKiIiIr2A6S0/M2bMaLPPzvHeeOMNqqqq+OlPf9pi+969e4mNjeW6664jOzubAQMGMHfuXKZOndpdkUVERMSDmd7y0x5Op5PXXnuN6667rkVfHofDQXZ2NpWVlcybN48XX3yRtLQ0br/9dtatW2diYhEREXFXprf8tMdXX31Fbm4uV199dYvtNpuNjRs34uvri81mA2DUqFFkZWWxaNEizjzzTDPiioiIiBvziJaflStXMnr0aPr169dqX1BQUHPh0yQlJYWCgoKeiiciIiIexCOKn82bNzNlypRW2zMyMhg3bhybNm1qsX3Hjh3qBC0iIiJtcvvHXk6nk8zMTObMmdNqX0pKCkOHDuXhhx/moYceIjIykrfeeoutW7eydOlSE9KKSHdKX7oNe7WD0CAbc6868ejQtl4HdOg1ItJ7uX3xU1ZWRn19fas5fwCsVisLFy7kySefZN68eVRUVDBixAheeeWVFsPhRaR3sFc7KK90dOp1IiJN3Kr4WbBgQatt0dHR7Nmz54SviYqK4tFHH+3OWCLiZkKDbGrNEZFOc6viR0SkvdSaIyKd5REdnkVETiV96bbm1iAzXi8inkMtPyLSK5xuS5BakkS8h1p+RERExKuo+BERERGvouJHREREvIqKHxEREfEqKn5ERETEq6j4EREREa+i4kdERES8iub5EZFeo2nZi84sfioi3kPFj4j0Kp1d/FREvIcee4mIiIhXUfEjIiIiXkXFj4iYQguJiohZ1OdHREyhhURFxCxq+RERERGvouJHREREvIoee4lIr9Q0509IoB8AFmcDrvp6rH5+JicTEbOp+BGRXsfqaiA8eycBh/YRWXSQ5Go7VpeTdSusBCQkEDygP9Fnn0X05IlYbTaz44pID1PxIyK9hqWhnsTMTYzPWE9gXWXrA1wuanNzqc3NpXjdenyCg+lz6SUk/eCKHs8qIuZR8SMivUJISS4pm98noKIYAEdQKEeThlHXbyjF/uE0+Plz1xUjqTl8mPIdOylc/SmOoiIOvfEWhas/IXLIdEoTh5h8FSLSE1T8iIhHMwyDw0uXMebTf2MxDOoCQtgzaAqMnURFnYvQIBt13w6r94+Owj86iogxo+n/o6sp+nI92a+8Sl3hUUYWLqVg1DmkL/EDi0Xrgon0Yip+RMRzuVxkPvc8hStXYQFKB4xkz+gLKHb40NfHFzjxXEIWHx9ip55N1MTx5Lz2OnnvfUD8ji/wK8ojY/LlPXYJItLzNNRdRDySxdnA8PX/oXDlKrBayRw3k4PTrqTBFtih8/gEBDDo9jnsmTgLl48vUflZjPhyKc66um5KLiJmU/EjIp7HMEjZ9B7Refuw2mwMu/9X5A8ad1qnPNp/FFkX3kiDr42IowfZ/cfHVACJ9FIqfkTEsxgGfTb+l9jDGbgsVoY/+Guip0zuklNXx/Zj59lX4/Txo/yb7ex7+lkMw+iSc4uI+1DxIyIeZWDOFmL3fAXA3omziBjbtR2T7TF92Xn2D7H4+lK8dh2H3lzSpecXEfOp+BERjxFafJjhez8D4EDaeRT1G9Et36citj+D77wdgEP/fpOiL9d1y/cREXNotJeIeIT68nKGbXgHq+GidOBIjgyd1OlzpS/dhr3aQWiQ7YRD2uMvvIDqg4fIXfEemc+mEzJ4cKe/n4i4F7X8iIjbMwyDfU8/i3+NncrgSA5PmQUWS6fPZ692UF7pwF594qHwAANv/jGhw1JxVlez76/PgOHq9PcUEfeh4kdE3F7hqtWUbt6Cy+rD5jGX4vLz75Hva/HxIeWe/8MnMJCKXbvpu2d9j3xfEeleKn5ExK3VFRVzYNFiAHJGTKUyNKZHv39AfDyDvu3/03/XFwSVF/bo9xeRrqfiR0TclmEYZKU/j7O6mpCUoRxJ6Xw/n9MRO30aUVMmYzVcDPn6vxguPf4S8WQqfkTEbRWvW0/p5i1Y/PwYevfPwWLOR5bFYmHQ7XNo8LURVnKE/P9+bEoOEekablX8pKenc+ONN7bY9sADD5Camtria9q0ac37XS4XzzzzDFOnTmXMmDHceuut5OTk9HR0EemE0CAb6Uu3kb50W6t9zro6sl9eDEDfH1xBUL++PZyuJf+YaHJGNn725Lz2TxylpabmEZHOc5viZ/HixTzzzDOttu/Zs4c777yTL774ovlr+fLlzfvT09N54403+OMf/8ibb76JxWLh9ttvx+E4+SgOEXEP9uq2R10dXrqMuqNF+MfGkHTlbBOStZY3+AzskQk4q6o5+PobZscRkU4yvfgpKCjgtttu4+mnnyY5ObnFPqfTSWZmJmlpacTGxjZ/RUVFAeBwOHj55Ze56667mD59OsOGDeOpp56ioKCAjz9Ws7SIp6rNz+fIf94BIHnOLfj4d8/orqaWp3+8v6t9L7BY2T/mAgAKVq6i6kB2t+QSke5levGzc+dOwsPDWbFiBWPGtJxsLDs7m7q6OgafYHKxjIwMqqqqmDJlSvO2sLAwRowYwcaNG7s1t4h0n5zX/41RX0/46DSiumjdrhOxVzuorKlv//HRfYk++ywwDA68vFhrf4l4INNneJ4xYwYzZsxoc9/evXuxWCy8+uqrfPbZZ1itVqZPn868efMIDQ0lPz8fgMTExBavi4uLIy8vr9uzi0jXe+XFj0j57AsABt7yYyynMZlhdxl40w2UbPiK8m+2U7r5a6ImjDc7koh0gOktPyezb98+rFYrSUlJLFy4kPvvv581a9Ywd+5cXC4XNTU1ANhsthav8/f3p66uzozIInKaYjetBKCw3whCBg0yOU3bAuLj6XPpJQB8/dwi0pdsNTeQiHSI6S0/J3PXXXdx8803ExYWBkBKSgqxsbFcc801bN++nYCAAKCx70/T7wHq6uoIDAw0JbOIdF7ZN9uJLDiAy2Ll4IipZsc5qaQfzCb/o/8RWFqAbd92YKzZkUSkndy65cdisTQXPk1SUlIAyM/Pb37cVVjYcsbVwsJCEhISeiakiLRb+tJtPP6PjW12MDYMg4Ov/xuAg33TqA2J7Ol4HeIXFkqfy2YBMGDX5xhOp8mJRKS93Lr4+eUvf8mcOXNabNu+fTsAQ4YMYdiwYYSEhLBhw4bm/RUVFezatYsJEyb0aFYRObWmBUXb6mBcvn0H9ow9uKw+ZA7q3k7OXaXPZZfS4OdPkL2Yo5+vNTuOiLSTWxc/s2bNYu3atTz//PMcPHiQNWvW8Jvf/IZZs2YxePBgbDYbN9xwA08++SSrVq0iIyODX/ziFyQkJHDhhReaHV9EOuDQm0sAyE8eQ11ASIt9HRqO3oN8Q4I5nNJYqB16c4laf0Q8hFv3+TnvvPN4+umnWbhwIQsXLiQ0NJRLL72UefPmNR9z991309DQwG9/+1tqa2uZOHEiixYtatUJWkTcV9jRg1Ts2InF15fDKVPguKWz2poE0V3kDR5P370bqM3NpXj9V8ScfabZkUTkFNyq+FmwYEGrbTNnzmTmzJknfI2Pjw/33nsv9957b3dGE5Fu1C/jSwDizp+BIygMKt2r2Elfuo2QQL829zn9/MkdPJ7+GV9y+O3/EH3WFLccni8i33Hrx14i0vsFl+UTWZgNVit93WQZi+OdaiLE3CETsNpsVGVlUb7tmx5MJiKdoeJHREyVtPcrAGLOOYuA+DiT03ROg38Q8Rc1LntxeOky4LuRbW0t2ioi5lLxIyKm8asqJ/bwbgCSrrjc5DSnJ+mKy7H4+DSOWtuX2TyyzZ37K4l4K7fq8yMivdOJ+szE7l6PxTAoix1AyOCOz+bctDDpifrj9JTQIBuL1hyhf/8RRB7YTt6770PMWaZmEpETU8uPiHS7tvrM+NTXErVvC0DzcPGuOrcZ7NUOspMb1/gq+mItthq7yYlE5ERU/IiIKeKzt+PT4KAqLIay+GSz43SJqsgEwkYMx3A6Scz62uw4InICKn5EpOcZBolZm4HGeXLoRUPDm5a8SDiwBavT/BYpEWlNfX5EpMeFHtlHYFUZDbYACvuPNDtOK019iezVDhKjgzv02qhJE/GPj4OCQpJyd1MRrqV2RNyNWn5EpMfF7NkIQMngsbh83XM29pOtQ3YyFh8fEi/+PgADDm0Dw+iOeCJyGlT8iEiPCq4qJSw3CwMoTp1odpxuETfjPJxWX8LtRwktyTU7jogcR8WPiPSo/ocaZ0AuTRiMIzTS5DTdwy8slKJ+wwFI3K+OzyLuRsWPiPQYi7OBpNzG1dnzk8eaG6ab5Q06A4CYwxnUV1SYnEZEjqUOzyLSY6Lz9uFfX0N9YCglCYMJNTtQJ7S3M3RlVCJlYfFEVBRQ8PGqNtctazpPaJCNuVeN6c7YInIMtfyISI+JP9C4zlXJ4DFg9dyPn/Z2hj7YbzQABStXYbTR8VlLYIiYw3M/fUTEo9jspY2rtwMlQ8aZG6aH5Cak4vTxozY3D/vuDLPjiMi3VPyISI+IytoKwNHo/r22o/PxnL42ivoOA6Bg5WqT04hIE/X5EZFu07zoqMtFVOZWAA71TcPnuOPcZYHS7lAwcDTxOdspWvslybfdim9QoNmRRLyeWn5EpNs0LToaWbAfvxo79bZACuIGn/TY3qYiui8Bffrgqq2l+MsvzY4jIqj4EZFukL50G/94f1fznxMObAWgcMAoXFYva3C2WIi/YAagR18i7kLFj4h0uWNbcXyr7UTlZwGQP9A7h3PHnXcuWK3Yd2dQffiI2XFEvJ6KHxHpVlFZW7EYBlWxfakJizE7jilsUZFEjm8c4Va4Sq0/ImZT8SMi3ccwmkd5FQ89w9QoTZ2qj30c15Pizz8fgMJPPsVwOk3JICKNVPyISLcJLcnF316K08eP8v4jzI5jaqfqyAln4BceRn1pGaVfbzElg4g0UvEjIt0m9uBOAIqTUnD52UxOYy6rnx+x504H1PFZxGwqfkSkW1icTmIP7wagsP8ok9O4h7jzG0d9lW7cRH15uclpRLyXih8R6RahuZn4OWqoDwyhLHaA2XHcQvCA/oQMGYzhdFK0dp3ZcUS8loofEekWkQe2A1A2cKRHL2La1WKmTQXg6JrPTE4i4r30iSQiXc6nvpaww3sBKE0e3alzmD06q7vETj0HLBbsGXvwryozO46IV1LxIyJdLvrIXqzOBqpDo6mJSuj0eXrjkhe2qEjC0xr7QMUe2m1yGhHv5GXzzItIT4j7dpRXYf+RYLGYnKbntHeB1tjpUyn/Zjuxh3ZCkrnzH4l4I7X8iEiXqisqJvxoDgBH+5k/t09Pa09rVfSZU7D4+RFcUURoZVEPJRORJip+RKRLHf3scyxAZVx/6oIjzI7jlnyDg4ma0Nji0ycv46THpi/dxuP/2Ej60m09EU3EK6j4EZEu1TSKqTQ5zeQk7q1p1FdSXgYYxgmPs1c7KK90YK929FQ0kV5PxY+IdJmq7Gyqs3NwWX0oH+B9j7w6ImrCeBp8/QmstRNWfNjsOCJeRcWPiHSZo599AUBpwiCc/oEmp3FvVpuN4qQUAGIP9q7h/CLuzq2Kn/T0dG688cYW21avXs2VV17JuHHjmDFjBo8//ji1tbXN+48cOUJqamqrryVLlvR0fBGvZhgGxWu/BLyzo/OJNI0Aa6vPTmG/kQDEHNmNq753DekXcWduM9R98eLFPPPMM0ycOLF526ZNm/j5z3/OvHnzmDlzJjk5Ofz+97+nrKyMxx57DIA9e/bg7+/PypUrsRwzpDY0NLTHr0HEm1Vl7ac2vwCrzUZJwmCCzQ7kRk7UX6c8rj+1tmACHFWUbd1G1MQJPZxMxDuZ3vJTUFDAbbfdxtNPP01ycnKLfW+88QZTpkzhJz/5CQMGDGDatGn84he/YMWKFTgcjR8me/fuJTk5mbi4OGJjY5u/AgICzLgcEa9V9MVaACInjsfl690ruLebxUpuYirQOEpORHqG6cXPzp07CQ8PZ8WKFYwZM6bFvltvvZX77ruv1WsaGhqorKwEGlt+hgwZ0iNZRaRthmFw4H+fArAzUIuYdkReQmPxU/LVJpx1dSanEfEOpj/2mjFjBjNmzGhz34gRLfsNOBwOXnnlFUaOHElUVBTQ2PITGxvLddddR3Z2NgMGDGDu3LlMnTq127OLSKPKfZnYqspp8PGjMDYZHCceui0tlYUnUBsURkB1BWVbthI9ZbLZkUR6PdNbftqroaGB++67j8zMTB566CGgsRjKzs6msrKSefPm8eKLL5KWlsbtt9/OunXrTE4s4j2Kvu3oXBA7CMP35Es7yHEsFoqTGlt/mv4eRaR7md7y0x5Nxc2GDRt45plnmh+P2Ww2Nm7ciK+vLzZbYx+DUaNGkZWVxaJFizjzzDPNjC3iFQyXi6IvGv/RzktIQaVPxxUlDSNp38bmR18+/v5mRxLp1dy+5aewsJDrr7+eLVu28NJLL7V6RBYUFNRc+DRJSUmhoKCgJ2OKeC373n04iopo8LVxNCb51C+QVuxRfbDFxOCqraVsy1az44j0em5d/JSXl3PTTTdRUlLCv/71L6ZMmdJif0ZGBuPGjWPTpk0ttu/YsUOdoEV6SFOrT0niUFw+HtGY7H4sFmLOavx8K1qrR/Yi3c2tP6kee+wxDh06xN///neioqI4evRo876oqChSUlIYOnQoDz/8MA899BCRkZG89dZbbN26laVLl5qYXKR3SV+6DXu1g9AgG3Ov+m5UpuFyUfxlY/FT1HeYWfE8QtNkh0CLv8Mm0WefRe6K9yjduAmXw4HVpukCRLqL2xY/LpeLDz74gPr6em666aZW+1etWkXfvn1ZuHAhTz75JPPmzaOiooIRI0bwyiuvkJqaakJqkd6paXHNVtsz9uAoLsEnKIjS+GSocZmQznOcbHHS0JSh2GJicBQVUbplK9GTJ/VgMhHv4lbFz4IFC5p/b7Va+eabb075mqioKB599NHujCUiJ9D0yCtq8iQMH19AK493lsVqJeasKeSueI/itetU/Ih0I7fu8yMi7stwuSj6srF/Ssw5Z5mcpneIPrvx77Hkq424HCokRbqLW7X8iIjnsO/ZS31pKT5BQUSMGQ27tpodySM09f0JCWw9KUBoylBs0dE4iov519+W40oZZUJCkd5PLT8i0inF69YDEDVxAlY/ze7TEfZqB5U1rVdxt1itRJ/VOD9ZcNaONo8RkdOn4kdEOswwjObiJ/rMKac4Wjoi5uzG4icqbx8WZ4PJaUR6JxU/ItJhVfsPUFd4FKvNRsQZY82O06uEpqZgi47Gt8FBSN5+s+OI9EoqfkSkw5pafSLHj9NSDF3MYrUSPaVxpFf4wQyT04j0Tip+RKTDvnvkpfXzukPUtyu7hx/eCy7NnSTS1VT8iEiHVB86TM3hI1h8fYmccIbZcXql8JEjqLcF4FtXTVjxYbPjiPQ6Kn5EpEOaWn0ixqThGxxscpreyeLjQ0li4/qE0bl7TU4j0vuo+BGRDtEor55R3CcF+Lb4MQyT04j0LprkUETazb+qjKr9B8BqJWrSRIATTtgnp6csPhmnjx8B1RWEVRRCaL/mfSdbIFVETk0tPyLSbtFHGh/BhI0Yjl94OHDiCfvk9Lh8/LAnDQYgoTCzxT57teOki6SKyMmp+BGRdmvqfxJzlh559YTyfsMASCjIPMWRItIRKn5EpF386yqbRx5FTZ5schrvYE8aistiJbSqmAB7idlxRHoNFT8i0i7xBVlYgJCUofjHRJsdxys4/QMpj+0PaNSXSFdS8SMi7dLU7yR6ilp9elKLUV8i0iVU/IhIK+lLtzWPKALwqa8luuQQAFGTJ5kVyyuV9BkKQFhJLnXFxSanEekdVPyISCvHjyaKzN+P1XBRHRpNUN8kE5N5H0dgKKURiQCUbNhochqR3kHFj4icUnTuPgCKE4eanMQ75cc1zvZcvH6DyUlEegcVPyJyUq76eiLz9wNQO2hEq0dicvpCg2wn/XttKn4qduykobKyJ6OJ9Eqa4VlETqp8x058G+qotQVTHZOkyfW6ycn+XquDI6kKiyW44iglGzcBWlNN5HSo5UdETqpkw1cAFMQNAovF5DTeq/jbjs/F678yOYmI5+tU8bNx40aqqqra3FdRUcH7779/WqFExD0YLhclXzV2si2IG2xyGu9WnNQ45L3s6y1YG7SciMjp6FTx8+Mf/5isrKw29+3atYsHHnjgtEKJiHuozNqPo7iEBl8bxVH9m7c39VH5x/u7TEznXarC4/GPi8XlcBBRcMDsOCIerd19fu6//37y8vIAMAyD+fPnExIS0uq47OxsYmJiui6hiJim6ZFXWXwyLp+WHxfq+9PDLBaiJk8m7933iM7dS8m3LUEi0nHtbvmZOXMmhmFgGEbztqY/N31ZrVbGjh3LY4891i1hRaRnNT3yaupvIuaKPrNxgsmovEwsLqfJaUQ8V7tbfmbMmMGMGTMAuPHGG5k/fz6DB6sPgEhvFVBZSnXOQbBaKUkYAmroMV3YsGH4hYdBeQVhRYeAKWZHEvFInerz89prr6nwEenlor6d2DB81EictgCT00hokI3n/7OD4vjGOX+aJp4UkY7r1Dw/NTU1LFy4kE8++YSamhpcLleL/RaLhZUrV3ZJQBExR3Re40KaUZMnQYnJYQRo7GdV3GcoYXu/JipvH4ZhYNH0AyId1qni55FHHuHtt99m0qRJDB8+HKtV0wWJ9Ca+ddWEFR0BIGrSBPgox+RE0sSeOAinjx8B1RVUHThAyKBBZkcS8TidKn7+97//8Ytf/IKf/OQnXZ1HRNxAVF4mFgyCByUTEBcHqPhxF4avH2XxA4nO3UfJho0qfkQ6oVNNNg0NDYwePbqrs4iIm2jqTxI1eZLJSbxLe+dPKu7TOMy9aSoCEemYThU/55xzDp999llXZxERN+CsqyOisHESvWgVPz3OXu2gsubkMziXJAzGwELVgWxqCwp6KJlI79Gpx14XX3wxDz30ECUlJYwZM4bAwMBWx1xxxRWnm01ETFC2dRs+zgYcweG8uqmMkF2axdndNPgHURXfn5CCHFa8+A5X/05dEEQ6olPFz7x58wBYvnw5y5cvb7XfYrGo+BHxUCUbGic2LO+Xir2mXouZuqnyfqmEFOQQfiiD9KXbsFc7CA2yMfeqMWZHE3F7nSp+Vq1a1dU5AEhPT2fdunW89tprzdt2797NI488wo4dO4iIiODGG29kzpw5zftdLhfPPfccS5YsoaKigvHjx/PQQw8xYMCAbsko0psZTiclGzcBjf+4ivuq6JtK0qb/EVx4kJrScsrrfcyOJOIxOtXnJykp6ZRfHbV48WKeeeaZFttKS0u55ZZbGDhwIG+//TZ33XUXTz/9NG+//XbzMenp6bzxxhv88Y9/5M0338RisXD77bfjcGg6WpGOqsjYQ0NFBfV+AVTF6QcId+YIjaQqPBaLYRCVn2l2HBGP0qmWn+eee+6Ux/z85z9v17kKCgp48MEH2bx5M8nJyS32vfXWW9hsNubPn4+vry+DBw8mJyeHl156iSuvvBKHw8HLL7/Mvffey/Tp0wF46qmnmDp1Kh9//DGXXHJJxy9OxIs1reVVmjgYNH+X2yvuk0Jw+dHG2bij1VIn0l5dXvyEhIQQFxfX7uJn586dhIeHs2LFCv72t79x5MiR5n2bNm1i4sSJ+Pp+F3PKlCm88MILFBcXc+TIEaqqqpgy5bv1bcLCwhgxYgQbN25U8SPSAYZhULK+cei0FjL1DMWJQ+m/ey2RBQewjqgHbGZHEvEInSp+MjIyWm2rrq5m8+bNzJ8/n9/97nftPtexC6YeLz8/n5SUlBbb4uLiAMjNzSU/Px+AxMTEVsfk5eW1O4OIQM2hQ9Tm52Px86M0fhDBZgeSU6qKiMcRFIatuoKY4oPUhQ83O5KIR+iydu2goCCmTp3Kz372M/70pz91yTlra2ux2Vr+JOPv7w9AXV0dNTU1AG0eU1dX1yUZRLxF8betPhFjRuPyVQuCR7BYqPi2Y3p8YZbJYUQ8R5c/1E9MTCQrq2tuwoCAgFYdl5uKmqCgIAICGleabuuYtuYeEpETa+rvo1mdPUvTqLz4o1lguE5xtIhAFxY/hmGQm5vLSy+91KnRXm1JSEigsLCwxbamP8fHxzc/7mrrmISEhC7JIOIN6oqLqdyXCRZL40Km4jEq4wfQ4OePv6OGsOIjp36BiHSuz8+wYcOwnGDiM8Mwuuyx18SJE3njjTdwOp34+DTOYbFu3TqSk5OJjo4mNDSUkJAQNmzYQP/+/QGoqKhg165d3HDDDV2SQcQbNLX6hKakYIuIMDeMdIzVh5KEIcQd2tk46ktETqlTxc/PfvazNoufkJAQzj33XAYOHHi6uQC48sor+fvf/86DDz7IbbfdxjfffMOrr77Kww8/DDT29bnhhht48skniYqKIikpiSeeeIKEhAQuvPDCLskg4g2aZnX+xi+RradYVFPcT3GfocQd2kl07l4Mw8BisZC+dBuAZnwWaUOnip+77rqrq3O0KTo6mr///e888sgjzJ49m9jYWO677z5mz57dfMzdd99NQ0MDv/3tb6mtrWXixIksWrSoVSdoEWlbQ1UV5dt3AJATMZDIUyyqKe6nLD4Zp9WHwKoyag4dIqh/f+zVmuhV5EQ6VfxAYyfjZcuWsWHDBioqKoiMjGTChAnMnj27eURWRy1YsKDVttGjR/Pmm2+e8DU+Pj7ce++93HvvvZ36niLervTrrRgNDVSHRFEVHEWk2YGkw5x+/hRH9Seu6ADFGzYS9G03ABFpW6c6PFdUVHD11Vczf/58tm3bRmVlJV9//TXz58/nqquuwm63d3VOEekmJV81TWyYcoojxZ3lxw0GoGTDVyYnEXF/nSp+/vznP5Ofn88///lPVq9ezZtvvsnq1av55z//SXFxMU8//XRX5xSRbuCqr6d009cAlGhWZ49WGDcYA6jcl0ldcbHZcUTcWqeKn1WrVjFv3jwmTGg5JHbChAncfffd/O9//+uScCLSvSp27sJZXY1fRAT2qD5mx5HTUOcfjD2qcZqRptF7ItK2ThU/VVVV9OvXr819/fr1o6ys7HQyiUgPKf72EUnUpAlwgukrxHM0rcnWtEabiLStU8XPoEGD+OSTT9rct2rVKgYMGHBaoUSk+xmG0TzEXbM69w5Njy7Ld+zEp77W5DQi7qtTo73mzJnDPffcg8Ph4NJLLyUmJoaioiLeffddlixZwvz587s4poh0taqs/TiKi7EGBBAxOg12bTM7kpymmtBoAvsmUXP4CJH5+ynqN8LsSCJuqVPFz8UXX0x2djYLFy5kyZIlzdv9/Pz42c9+xjXXXNNlAUWkezQ98oocNxar5sXqNaImTeTI4SNE5+5T8SNyAp0qfqqrq5k7dy433HADW7dupby8nLy8PK655hrCw8O7OqOIdIOmIdFRkyeanES6UvSUyRxZtpzI/Cwszgaz44i4pQ71+dm9ezdXXHEFixcvBiAsLIxp06Yxbdo0/vrXv3Ldddd12YruItJ9agsKqM45CFYrkRPGmx1HulDI0CH4RUbg2+Ag/OhBs+OIuKV2Fz+HDh3i5ptvpry8nCFDhrTYZ7PZ+M1vfkNVVRXXXXcd+fn5XR5URLpOU0fn8JEj8AsNNTmNdCWL1UrUpMbWvOg8LXQq0pZ2Fz8vvvgikZGR/Oc//+Giiy5qsS8wMJAbbriBt99+m6CgIBYuXNjlQUWk6xTrkVevFv3t6L2o3H0YLpfJaUTcT7uLn3Xr1nHbbbcRERFxwmOio6O55ZZbWLduXVdkE5FuUF9hp2LXbgCiJmmIe28UPjqNBl8b/rWVVGaqK4LI8drd4fno0aPtmr8nJSVFj71E3Fjpps3gchGcPJCA+DjSl24jJNDP7FjShax+fpTGDyL2SAYlG74iNKXtpUvSlzZObzD3qjE9GU/EdO1u+YmKiqKwsPCUx5WUlJy0dUhEzPXdrM6Nj7zs1Q4qa+rNjCTdoGnCw+KTLHRqr3Zgr3b0VCQRt9Hu4mfixIksW7bslMctX76c4cOHn1YoEeka6Uu3Nf90D+Csq6Nsy1YAoqbokVdvVpIwGJfFSs2hw9Tk5podR8SttLv4ufHGG9mwYQMLFiygrq6u1X6Hw8Hjjz/O559/zvXXX9+lIUWkc47/yb78m+246urwj40hODnZxGTS3Zy2AMpj+wPfje47vhgW8Vbt7vOTlpbGAw88wKOPPso777zDmWeeSd++fXE6neTm5rJhwwZKS0v5v//7P6ZOndqdmUWkk4rXf/fIy6KFTHu9kj5DiSzMpnjDVyTNvlyPuES+1aEZnq+//nqGDRvGokWLWLVqVXMLUHBwMOeccw633norY8ao45yIOzKcTko3bgK0kKm3KE4cyuCtH2PP2IOjrNzsOCJuo8PLW4wfP57x4xtnhC0tLcVqtWpJCxEPYN+7j/rycnyCgwkbqTWfvIEjKIzgwYOpysqidONGQJ/VItDJtb2aREZGdlUOEelmzWt5TRiP1ddXQ9y9QGiQjcyQfiSSxbYVq2HibLMjibiFDq3tJSKeyTCMVrM6a4i7d8iLHQRA4OFMrA3q8yMCKn5EvELN4SPU5uZh8fUlYtw4s+NID6oOi6UqMByry0lkwQGz44i4BRU/Il6g6ZFXxJg0fIMCTU4jPcpioSBuMNC41peInGafHxHxDCVfNc7zEjVpkvr6eKGCuCEMyvmaqLxM0EKnImr5Eent/Goqse/ZC0DkxAnq6+OFSiP60OAfhF99LeFFh8yOI2I6FT8ivVxUXiYAISlD8Y+OMjmNmMGwWqno27jWV1SeHn2JqPgR6eWicxtbfaI1saFXK+83DGj8/8EwDJPTiJhLxY9IL+bjqCWiMBuAqCmTzQ0jprInDsLp40tAdQXV2TlmxxExlYofkV4sKj8Lq+EisF9fgvommR1HTGT4+lEW17iYbdOcTyLeSsWPSC8WfWRP469nTjE5ibiD4j6N/X5K1n9X/IQG2Uhfuo3H/7FRK76L19BQd5FeyllbS2TBfkDFjzQqSRyCgYWqAweoLSxs3m6vdlBeqdmfxXuo5UeklyrbshUfZwO1QeEEJw80O464gQb/ICpi+gJQsmGjyWlEzKPiR6SXKvpyfeOvSalYLBaT04i7aH709ZWKH/FeKn5EeiFXfT2lmzYDUJyUanIacSfFiSkAlO/Yia+jxuQ0IuZQnx+RXqhs2zc4q6upCwjBHtXH7DjSQ5o6L59s+ZK6kAiCBvSnOucgkXlZ1Eac0YMJRdyD2xc/GzZs4Mc//nGb+/r27cuqVat44IEHWLZsWYt98fHxfPbZZz0RUcTtFK/b0PhrnxTQIy+vYq8+dcflqMmTqM45SHTePo4MV/Ej3sfti59x48bxxRdftNi2d+9efvKTn3DnnXcCsGfPHu68805uuOGG5mN8fHx6NKeIuzCczuZV3IuTUkxOI+4oevIkDr+1lMj8/eQ2aJ038T5uX/zYbDZiY2Ob/1xfX89jjz3GRRddxA9/+EOcTieZmZnMnTu3xXEi3qp85y4a7HZ8Q0Mpj+lvdhxxQ8GDB2GLicFRVERo3n4qopPNjiTSozyuw/Prr79OXl4eDzzwAADZ2dnU1dUxePBgk5OJuIeS9Y2PvKImTwSrx93i0gMsFkvz3E/hB3ebnEak53nUJ2NdXR0LFy7kpptuIi4uDmh8BGaxWHj11VeZMWMGF1xwAX/4wx+w2+0mpxXpeYbL1dzfRxMbysnEnH0mAOGH9mBxNpicRqRneVTx884771BXV8eNN97YvG3fvn1YrVaSkpJYuHAh999/P2vWrGHu3Lm4XC4T04r0vMp9mThKSvAJCiJizGiz44gbahoRtjyrgbqAEHzq64go1EKn4l3cvs/PsZYvX85FF11EZGRk87a77rqLm2++mbCwMABSUlKIjY3lmmuuYfv27YwZM8asuCI9rujLdQBEThiP1e/Ew53FuzWOCLNR3CeFPvu/JuZIBtnD+pkdS6THeEzLT0lJCVu2bOHiiy9usd1isTQXPk1SUhpHuOTn5/dYPhGzGYbR3N8n+szJJqcRT1Dct3ECzKi8fVhcTpPTiPQcjyl+vv76aywWC5MmTWqx/Ze//CVz5sxpsW379u0ADBkypMfyiZit6kA2tfkFWG02Is8YZ3Yc8QDlMf2o9w/Cz1FLdMlhs+OI9BiPKX4yMjLo168fgYGBLbbPmjWLtWvX8vzzz3Pw4EHWrFnDb37zG2bNmqURYOJVitc1ruUVccY4fAICTE4jHsFipbz/cAASCvaaHEak53hM8VNUVERERESr7eeddx5PP/00//vf/7j00kt58MEHueiii3j00Ud7PqSIiZoeeW3xSTQ5iXiS8gHfFj+FmaBBIuIlPKbD8/z580+4b+bMmcycObPnwoi4merDR6g+eAiXxUpu1ECz44gHqYwfQL0tAH9HDeFFhwD1F5Pez2OKHxFpW/rSbYR//SkDgbK4gThteuQlHWD1obhPCgnZ3xB9ZI/ZaUR6hMc89hKRttmrHUTl7AK+G70j0hHFSY3/38Tk7sHQoy/xAip+RDxcYEUxYZVFGBZr4yruIh1UFjeQel9/bLVV2DPU+iO9n4ofEQ8Xc6RxbSZ7n0E02AJPcbRIa4bVh4K4xtGxRWu/NDmNSPdT8SPi4WIOZwBQNmCkyUnEk+XFDwWgaO06DKcmPJTeTcWPiAerPniQ4IoinBYfyvupv490XlHMABr8/KkvLaVid4bZcUS6lUZ7iXiwoi8aH1EUxQzAZQuABkerY5oWsrRXO0iMDu7piOIhXFZfivukEJ+znaIv1hI+Si2J0nup5UfEQxmGQdEXawHITTh5q4+92kF5pYPKmvqeiCYe6mi/EQAUf6lHX9K7qfgR8VDVOTnUHMnFZfWhMG6Q2XGkFyiLHYBvWBj15RWUb99hdhyRbqPiR8RDFX3e2OpTkjCYBl9/k9NIr2C1EnPWFACOfv6FyWFEuo+KHxEPZBhG85Dkor7DTE4jvUVokI0vnAkA5H72Jc+/+bXJiUS6h4ofEQ9Utf8AtXn5WG02ShKHmB1HepH80ATqA0PwddTik60JD6V3UvEj4oGaOjpHThiPy9dmchrpVSzW5jmjYg/vNjmMSPdQ8SPiYY595LUzqL/JaaQ3KhvYOOorKncfzro6k9OIdD0VPyIexp6xh7qCQpw+fhTGapSXdL3qmL7UBoXh2+CgdLP6/Ujvo+JHxMMcXfM5AMVJKRi+fiankV7JYqGo73Dgu0esIr2Jih8RD+JqaGh+5FXYTzPwSvc5+m3xU7pxM86aGpPTiHQtFT8iHqRs6zYaKirwCw+nLG6g2XGkF6uKiKcmJBKXw0Hxhq/MjiPSpVT8iHiQo2s+AyDmnLPBqttXupHF0ty6ePTTz0wOI9K19Okp4iGcNTWUbNgIwEa/vianEW9wtH9j8VO27RscpaUmpxHpOip+RDxE8YavcNXVURMcQXFIvNlxxAvUhkQSmpoKLhdHP9NyF9J7qPgR8RBNo7yO9h8JFovJacQbhAbZ2B3eOJ1C0yNXkd5AxY+IB3CUlVO2dRugUV7Ssw7HDcVlsVKVtZ/qg4fMjiPSJVT8iHiAoi/WgstFyJDB1IZGmR1HvEiDfxD2pMb141Y89ybpS7eZnEjk9Kn4EfEATY8cYqdPMzmJeKPSQaMBiM7Zib1Ky12I51PxI+LmavLyqNy7D6xWYqaebXYc8UIVSUNp8PUnqLaCsCI9+hLPp+JHxM01dXSOGJ2GLTLS5DTijQxfP4r6pgIQd3CnyWlETp+v2QFE5MTSl2wl9Z2PCABiz5verteEBtlIX7qNkECt+yVd52i/kSRkf0PMkQxcDgdWm83sSCKdppYfETdmOZhFQHU5Db42os+c0u7X2asdVNbUd2My8Tblsf2pCQjBt75OK72Lx1PxI+LG4nO+AaCo73B8/P1NTiNezWIhN2EYAIWffGpuFpHTpOJHxE01VNcQc3gPAAUD00xOIwKHkxrnmCrd9DWOsjJzw4icBhU/Im6q+Mt1+DjrqQyKxB6VZHYcESpDoqmI6oPhdHL0kzVmxxHpNBU/Im6qcNVq4NuftrWchbiJggGNc/4UrFyNYRgmpxHpHBU/Im6oJi+Pil27MbBwuM9ws+OINCvqNxyrzUbN4cON80+JeCCPKH6OHDlCampqq68lS5YAsHv3bm644QbGjh3Lueeey6JFi0xOLHJ6Cld/CkBZ/EDqAkLNDSNyDKefP9FnnwVAwcpVJqcR6RyPmOdnz549+Pv7s3LlSizHNP+HhoZSWlrKLbfcwgUXXMDDDz/M1q1befjhh4mIiODKK680MbVI5xhOZ3Px0/SIQcSdxF8wg6OffErR52tJnnMLPgEBZkcS6RCPKH727t1LcnIycXFxrfa9+uqr2Gw25s+fj6+vL4MHDyYnJ4eXXnpJxY94pPLtO3AUFeETHExxn6FQ4zI7kkgLYSNHEJCYQG1ePsVfriNuxnlmRxLpEI947LVnzx6GDBnS5r5NmzYxceJEfH2/q+OmTJnCgQMHKC4u7qmIIl2mYNUnAMROOwfDxyN+PhEvY7FYiDt/BtDY8VnE03hE8bN3716Ki4u57rrrOOuss7j22mv5/PPG9Y7y8/NJSEhocXxTC1Fubm6PZxU5HQ2VlZSs3wDQ/I+LiDuKm3EuWK1U7NxFjT5rxcO4ffHjcDjIzs6msrKSefPm8eKLL5KWlsbtt9/OunXrqK2txXbcGjP+386EW1dXZ0ZkkU4r/GQNLoeDoAH9CRky2Ow4IifkHx1N5LixABR+21op4incvk3dZrOxceNGfH19m4ucUaNGkZWVxaJFiwgICMDhcLR4TVPRExQU1ON5RTrLMAzyP/ofAAnfm9mic7+IO4q7YAalm7+mcPWn9L/uRzz/nx3Yqx2EBtmYe9UYs+OJnJDbt/xAYxFzfOtOSkoKBQUFJCQkUFhY2GJf05/j4+N7LKPI6arYtZuaw4exBgQQe+40s+OInFLUxAn4hoXhKCmhdPPX2KsdlFc6sFc7Tv1iERO5ffGTkZHBuHHj2LRpU4vtO3bsYMiQIUycOJHNmzfjdDqb961bt47k5GSio6N7Oq5IpzW1+sROPQff41otQ4NspC/dxj/e32VGNJE2Wf38Gvv+APkf/dfcMCId4PbFT0pKCkOHDuXhhx9m06ZNZGVl8dhjj7F161buvPNOrrzySiorK3nwwQfJzMxk2bJlvPrqq9xxxx1mRxdpt/qKCoq/XAdAwvcuavMYe7WDypr6nowl0kpTIZ6+dBvQ+IgWoPTrrfhXlpmYTKT93L74sVqtLFy4kLS0NObNm8fs2bPZtm0br7zyCqmpqURHR/P3v/+dAwcOMHv2bJ577jnuu+8+Zs+ebXZ0kXYrXPUJRkMDIUMGq6OzuD179XePtgITE4gYNxYMg8QDW8wNJtJObt/hGSAqKopHH330hPtHjx7Nm2++2YOJRLqO4XKR/9/GR17xM9tu9RFxZwnf/x5lW7YSn/0N1v6TAdspXyNiJrdv+RHp7cq376A2Lx+foCBip55tdhyRDouacAb+sTH4OWpILNhrdhyRU1LxI9LDju0vAcd0dJ4+DZ/AQLNiiXTIsX1/LD4+xF90IQADDm07xStFzKfiR6SHHdtfwlFaSsmGrwBI+N6FZsYS6bBj/1+Ov/B8XBYrkWV5BJcVmJxM5ORU/IiYqODjVRhOJ6HDUgkeOLBVq5CIp7BFRlKclApAwn51fBb3puJHxCSu+nryP2ycGyXh+43DhY/9SVrE0+QNGgdA3MGdNFRVmZxG5MRU/IiYpPjL9ThKSqgPDGFxtr8mMBSPVxHTD3tIND7Oego/WWN2HJETUvEjYgbDIHfFuwDkJo+jrNqpCQzF81ks5PRrXNMr/4MPMVwukwOJtE3Fj4gJQkuOUJmZhcXPr/lRgYgnOn7plcN9RtDg50/NkVxKN39tcjqRtqn4ETFBn32Na9XFTp9Gg3/QKY4WcW/HLr3i9LWRnzwWgNx33jUxlciJqfgR6WH+VeXE5O4BYHPEMJPTiHS93MHjsfj4UL59B5X795sdR6QVFT8iPSxx/9dYDIOy2AGUBEa32q8V3MXTOYLCiD77LECtP+KeVPyI9CBnTQ0JB7YCcGToxBMepxXcxdP1uWwWAEWfr6WuuNjkNCItqfgR6UGFqz/Ft76OutAoShO0erv0XqFDhxA2cgSG00neex+YHUekBRU/Ij3EcLnIfe99AI4OnwwWi8mJRLpXn8svAyD/vx/jrKkxOY3Id1T8iPSQkg1fUZubR4OfP6WDxpgdR6TbRU0cT0CfRJxVVRSs+sTsOCLNVPyI9ADDMDi0ZBnQOBLG5WczOZFI97NYrc19f/LefQ/D6TQ5kUgjFT8iPaBsy1aqsrKw+vuTO2SC2XFEekzcjPPwDQ2hNr+A4vUbzI4jAqj4EekRh5c2tvokzLxQkxqKV/Hx9yfx4u8DcOitpVryQtyCih+RblaxazcVO3dh8fWlzxWXmR1HpNs1zVWVvnQbAImXXoJPYCDV2TmUbNxkcjoRFT8i3e7w0rcBiJtxLv7RrSc1FOmN7NUO7NUOAPxCQ0m85NvWnzeXYhiGmdFEVPyIdKfKrP2Ubt4CVitJP5htdhwR0/S5bBZWf3+qsrK04KmYTsWPSDdq6usTc87ZBCYmmJxGxDx+4eEkfH8mAIffUuuPmEvFj0g3qT58mKIv1wGwNmykyWlEzJd0xWVYbTbse/by2nPLO/TaY/sQiZwuFT8i3eTwkrexAPlxgznqH2F2HBHT2SIjib/oAgBitn7aodce24dI5HSp+BHpBlXZORxd8zkAmYMmm5xGxH0k/eAKXFYfwosOU75jp9lxxEv5mh1AxFM0NbnPverUS1Mc/Ne/wTAoSkqlPDyB8O4OJ+Lm0pduw17tIDE6mKCBo0ncv4V1z7zM/gtvZO4Px5odT7yMWn5E2qm9ze72PXsp2bARrFZyRkztgWQi7s9e7aC80kFlTT2HU6fgsvoQWpCNT/Zes6OJF1LxI9IFmjpjGoZBzmuvA1AxZCw1YTEmJxNxP3VB4RSnNi7zMnD7J5r1WXqcih+RLtDUKlS+7RvKt+/AZfXh8Ci1+oicSEHaVBp8/QkpL+ToZ1+YHUe8jIofka5yTKtP3qBx1Aerp4/IiTj9gzicOgWAg6//G1d9vcmJxJuo+BHpItG5e6nMzMIaEMDh1DPNjiPi9nKHTKAuIIS6wkLyP/yv2XHEi6j4EemA4xdsbOZyMWDnZ0DjNP71AcEmpBPxLC5fPw5+Oyjg0FtLaaiqMjmReAsVPyId1Naor4QDWwmyF+MbGkKSVm4XabeCAWkE9k2iwW7nyLLlZscRL6HiR+Q01VfYGbCrsdWn/7U/wje4datPU4vRP97f1dPxRNyb1cqAG68HIHfFe9QVl7TYrftGuoOKH5HTdPBfb+DnqKUqLJaE7110wuPs1Y1znIhIS1GTJxE6LBWXw0HOP15rsU/3jXQHty9+ysrK+P3vf8+0adM444wzuPbaa9m0aVPz/gceeIDU1NQWX9OmTTMxsXiTqgPZ5P/3fwDsH3sBFh8fkxOJeB6LxULybbeCxcLRTz+jfGfby140taA+/o+NWuRUTovbFz/33HMP27Zt4y9/+QtLly5l5MiRzJkzh6ysLAD27NnDnXfeyRdffNH8tXz5cnNDi1cwDIP9Ly0Cl4ujfYdRHjvA7EgiHit06BDiL7oQgP0v/B1XQ0ObxzXNFK1FTuV0uHXxk5OTw9q1a3nooYeYMGECgwYN4sEHHyQ+Pp733nsPp9NJZmYmaWlpxMbGNn9FRUWZHV28QNEXX1KxcxdWm43stPPMjiPi8VaFj6beFkh1zkGWPrpIfX2k27h18RMZGcmLL77IqFGjmrdZLBYMw6C8vJzs7Gzq6uoYPHiwiSnFG1kbHGS/8ioASVfOpi5IExqKnK4ypy+7h54DQPy2T6k9rvOzSFdx61Xdw8LCmD59eottH374IQcPHuScc85h7969WCwWXn31VT777DOsVivTp09n3rx5hIaGmpRaeoOmFahDg2xtruLed896HMXF+MfFkjT7cnjzm+b+CCGBfiYkFnE/nbknDiWNYnD+ToKLc+mz+WNKx8/qxoTirdy65ed4mzdv5je/+Q3nn38+M2bMYN++fVitVpKSkli4cCH3338/a9asYe7cubi0UJ6chpP1KwiqOErfPesBSL71Znz8/Vu8TiNTRL7T4XvCYuHI5IsxgMjsHYQfzem2bOK93Lrl51grV67kV7/6FWPGjOEvf/kLAHfddRc333wzYWFhAKSkpBAbG8s111zD9u3bGTOm9U/sIqfFMBiy+SOshovIiROImjLZ7EQivU5NdB/yB40jcf8WBm/5H5l9k82OJL2MR7T8/POf/+Suu+5i2rRpvPTSSwQEBACN/X+aCp8mKSkpAOTn5/d4Tun9ovduIqzkCA2+NgbfcTsWi8XsSCIe7UQTgOaMnEZ9QDBB9mIStn1qTjjptdy++PnXv/7FH/7wB66//nr++te/YrPZmvf98pe/ZM6cOS2O3759OwBDhgzp0ZzS+/lXlZO4ZRUAOSOn4x8bY3Iikd6hrUdjDbZADk9p7O8Tu2sdocVHzIgmvZRbFz8HDhzg0Ucf5cILL+SOO+6guLiYo0ePcvToUex2O7NmzWLt2rU8//zzHDx4kDVr1vCb3/yGWbNmaQSYdCnD5WLo5g/wqXdQHt2XvMHjzI4k0utV9EulsP9ILIbB0E3vY2n4rkA64SLDIu3g1n1+/vvf/1JfX8/HH3/Mxx9/3GLf7NmzWbBgAU8//TQLFy5k4cKFhIaGcumllzJv3jxzAkuvlf/hf4k4moPTx499Ey4Bi1v/3CDSa2SNuZDIowcJqiwhcctqKkZ+N6eWJjqUznLr4ufOO+/kzjvvPOkxM2fOZObMmT2USLxRTW4u2a82rjeUd8b51IZEmpxIxHs4bQEcOnMWg1b/m5iMDeTFDaY8MMHsWOLh9OOryElYXE72/vmvuOrqKIsdQHHqRLMjiXgde9JQ8geOxgKNj58b1OIjp0fFj8hJ9N/5OZWZWfiGhLB3wiWg0V0ipjgw+nwcQWEEVpUxMuMTs+OIh1PxI3IC0cUH6bu3cTLDIT//KY6gsFO8QkS6i9PPn4NnX4EB9Duyk4j9282OJB5MxY9IG2w1dsZ98wEWIH7mhUSfOcXsSCJeryphIIeGnw1A3w3vE2DX2l/SOSp+RI7jamggdcM7+DuqqQyPI3nOLWZHEpFvHRx+NsWRffFpcDB8/X+wqv+PdIKKH5Hj5Lz2OuHFh6n3tZEx+YoWa3eJiMksVraMvpj6gGCCK44y5OuPMAzD7FTiYVT8iByj8NM15C5fAcA3Iy+iNjTK5EQicry6gBBypl2Fy2Il7tAu8t593+xI4mHcep4fke6UvnQb9moHoUE25l41Bvu+TDKfex6AQ6lnkp+QQt9vZ5ENCfQzOa2IHKsqfgAHRs9g8LaV7H95MSuzHfz47h+YHUs8hFp+xGvZqx2UVzqwVzuoKyom47HHMerriZw4npyR01ocd/y6QyJivrzB4ykZNAaLYTBgzVKqDx02O5J4CBU/4vV86uvY9YdHcBSXENi3Lym/+D/N5yPiCSwWDk+5hPKYvvg21LH7j4/iKCs3O5V4ABU/4tUsLifD1v+H6uwc/CIiGPH7B/ENDjY7loi0k+Hjy+4pP6AmOILa/AJ2/+ERnDU1zfu1+Km0RcWPeC/DxZjtHxFZmI3V35/hv32AgPg4s1OJSAc1+Aex8+yr8Q0LozIzi4wFT+Cqb3xUba92aAFUaUUdnsXrpC/dRkiAL4O3fkxi/h5cFisj7v8VoUOHtPscod92hLZXO0iMVkuRiNlqQ6MY8bvfsOO3D1G2dRt7n3qa1F/+wuxY4qbU8iNex15VR+jaD0ncvwUD2DvxUiLHn9Hx83zbYVqdoUXcQ2jKUIb9+l4svr4Ur13Hvmf+BpoDSNqglh/pFZqe6c+9asxJjzMMg4HbPyFu31cAbB9xARX9hrc4j4a1i3ie0OZpKfyxnv0DBn62hKOfrmHogBL2jf++2fHEzajlR3qF9jzXNwyD7JcX0/fbwidz7EUc6je61XnUkiPimZru3yMxg9iSdjGGxUJ8znZSNr6H4XSaHU/ciIof6ZXSl27j8X9sbG4RMpxO3v3lH8ld8R4AhyddTP7gjj/qEhHPkJeYSs7UK5tngc74059x1tWZHUvchIof6ZWOncDQWVdHxp/+THTWVgyLhb3jL6Y4dYLZEUWkm5UPGEHGlNm4rD6UrN/Azt89TH255gESFT/Sy0Va6lk1915K1m/AZfUhe9oPKRw4+tQvFJFeoaTPUHaccw2+ISHY9+zhm/seoPrwyWeCPr7lWHofFT/Sa4XYixjy0SKCi47QYAtgxznXUNF/mNmxRKSHVcT2J+3xRwhIiKc2v4Dt9z9I+Y6dJzz+2JZj6Z002ks82olGZ0UfzmDopvfxddZTExxB9gXXU+EbSmgb5wjV4qUivVpokI3F64vxmXoDZ37zLvY9e9n50P9j0O1ziJ95IRYtZ+N11PIjHu340Vkuh4P9f3+F4RuW4+usx56QzLbzfowjLLpD5xGR3sVe7aDM5cfIP8wn+uwzMRoayHr+BfY88WcaqqrMjic9TMWP9BpB5UfZ9qv7yXu3cURX1sDx7D//ehr8g0xOJiLuwsffn9Rf3cPAm3+MYbFSvHYdW3/xK+x797Xr9VorrHdQ8SOezzBIzNzE2NWLqc45iF94GDvPuoqM1Olg1f/iItKSxWolafblbDv3BmqDwqkrKGT7rx/k8LLlp5wPSGuF9Q76l0E8mn9lGYNWvc7gbSuxupxEjj+Dsc88RWli+9fpEhHvVBnVhy3n39L4GMzpJOfV1/jmvgcILs03O5p0M3V4Fo/krKvjyH/eYfzHS7G6nDitvmSPPo8bf3eHOi+KSLsFRYSxOvwCopzRDNy+msrMLMZm7id8wBlUnXlhu5fOgfYvsyPmU/EjHqdk4yb2v7SIuoJCrIA9IZm9oy+gJjRahY+IdJi9ph570ih+cOcVHFj0CkWfr2VQzmYcR/exP+08ipLaN0WGHod5DhU/4vaahqG7Du4nYeunhOYfAMAWHcW2IdNwDE2jRiO1ROQ02SIjSf3VPaxx9SH56/8SVF3BsA3vUBmxntI0PyLGjW3zdZoqw/Oo+BG3ZxzOJnbPWsJyswCw+PrS57JZ9Lv6KlYv2UGoWntEpAuVJgzm4Nk3Ma7wG2J2rSOkrIBdD/+RsFEjCY0Zhz26b4vj1eLjeVT8iKlO9IzccLko/2Y7uSveZezmLY3bLBZKBo+lIG0qn1uCSFx9oMVrNFmhiJxM02eEvdpBYnTwSY91+fhRMGY62f3H0G/POvpmb6Vix07GsJOK6CSKBjcQPWUSFh+fDmVQvyD3oOJHTHX8T0wN1dUUrv6U/A8+pOZILtBY9JQOGkPWkCnY4uKap54PDbKd8nwiIsc62edHWxr8gzgw+nwuv/82Dr25hLyVqwkrPsKePz2Jf1wsibMuxscRCe08nz6j3IOKH+kSp/PTTGiAL/949j9EZO8k9vBuXLW1APgEBhJ3/nm8W98Pv7h46qodtO/jRUSkY07VcuwfG8OQn/+UFbZhJGZ9zaAj26krPEr2y68y2epDRb9UjiSNoCw+uYeTS2eo+JEu0dGfZgyXi8p9mSRvW0nckQz8aioBcAGBffuSeMn3iT13Or5BgdT+YyN6kCUi3a09n2P1ASEcHDmNq//fzzm65jPy3vuA6pyDROTsIiJnFw7/IPY7dhI1ZRJhI4Zj9dU/s+5I74r0mBdeX49vzl76lB4k8OBe/GqrSPp2X71fABUDR1CaPIr8kERCq/2ZGxRoal4REWi7VcjH35+Eiy4k/sIL+Nuz79M3dxfh+7djq6sm7/0PyHv/A3yCg4g8YxxRkyay7KAvgZHhXZ6tqQ9TaJBN/Yg6QMWPdJv68nIqdmdQsWs3FbsyGJmZicUwmvc7/WyUxA+mcshociP7ExIS2Pw8Ho3gEhE3cqJWIYvFQlVEPLl9+rFn2HT6lB0keP9OIvOzsFVVU/T5Woo+X8sooCYijsDovpTH9MdRloIt4vSLoebPTOmQXlH8uFwunnvuOZYsWUJFRQXjx4/noYceYsCAAWZH8xq+ddUElxVyeNkRqg4coCprf3OH5SYWwB4cTc2AFAqiB2L0G0RFnZPQIBuGOgGKiIc44agxqxV73xRyowaC4eKnkyMo+WojpZs2U51zkMCyQgLLCumT9TUbNyzHPy6WkMGDCBkyhODBgwhOHohfeLgma+0BvaL4SU9P54033uCxxx4jPj6eJ554gttvv5333nsPm01dZLtKQ3U1jqIiaguPUnMk99uvI9QcOcKU0jIAcr5o+Zqg/v0IGzGc0OHD+VeGg0JXIH3jQhqbaX18gJMvIigi4o5OOWrMYiVsWCphw1IZ+OMb+PNLa0ioyMN2ZD/hRw8SXHGUusLGr+J1G5pf5hsaQmDfvgT160tg3yQC4uLxj4/FPzYW35AQFUZdxOOLH4fDwcsvv8y9997L9OnTAXjqqaeYOnUqH3/8MZdcconJCd2bYRg4q6qoL6+gvqKC+vJyPv5kF761VYxL9Keu6Ci5mYfxqyrH11F70nPVhUZSE5VIWUgstv4DKAmNJzAijMqaekKq/KgLqoJ2Ns9qzh4R8WTHtw41+AdRPmA49tjBAPzyyhFU7T9AZdZ+KjMzqczaT01uHg32Suy7M7Dvzmh1Tp/AQPzjYhu/YmLwi4ggIauUQMOGX0Q4NUeS8AsPxyc4SEXSKXh88ZORkUFVVRVTpkxp3hYWFsaIESPYuHGjRxU/hsuF4XKBy4XhdDb+2en69lcnuFy4GupxOepxORwYDQ0sX7kbi7OBiyf3a9xe7+DTDQewOJ2cOSwaZ3U1DdU1OGtqcNZU46yuwVldQ0lRGb71Dqx1NVgMV4sc/b79NX9r46/Hdjuu9wvAFRpOTUgURnQc9qBIrLHxHHIFERQe2vzTUN+YxtYdZ019p+e10HwYIuLJTtQ6FBpk48UPM78tjIZQOXAAIcO/R0FhGYH2Eq4ZG071ocPs2LALn4pSAmvt+NZU4qypoTrnINU5B5vPNeSY83695vXG31it+AQG4hsc9O2vwfgEBeETFIhvUFDj7wMCsPj54eNvw+Jnw2rzw+pnw+pvw+rnh9X23a8WPz8sPlYsVp/GX318sFit0PRn6zHb8IyJHD2++MnPzwcgMTGxxfa4uDjy8vI6fL7CwkKcTifnn39+l+RrZhg4SksbixgAo/k/YICBcaJXtsvCv7fe9nSHzmDBsFjAasXAAhYrNn9fLD4+1NS7wOqDCwsuwNdqxWUYWC2WNn91GcYJjznZPnc7RjmUwxNyeFJW5Tj1PoD/vNNYLFVUOY45xoXF5SLY3wfD6aS21tE4gMQwsLhcgIHVMDCO+2G2p1mwNP9r9q+HLY2dPY/Z2/SLb0gIPgEBXfq98/Ly8GnnjNseX/zU1NQAtOrb4+/vT3l5eYfP5+/vj8PRDS0OFgu2qKiuP28PCNGIcxGRHhcWfOI+q8GB+mA+nq+vb7v7+Xp88RPwbeXocDiafw9QV1dHYCf+59i0aVOXZRMRERH3YzU7wOlqetxVWFjYYnthYSEJCQlmRBIRERE35vHFz7BhwwgJCWHDhu+GClZUVLBr1y4mTJhgYjIRERFxRx7/2Mtms3HDDTfw5JNPEhUVRVJSEk888QQJCQlceOGFZscTERERN+PxxQ/A3XffTUNDA7/97W+pra1l4sSJLFq0SBMcioiISCsWwzBOb4y1iIiIiAfx+D4/IiIiIh2h4kdERES8ioofERER8SoqfkRERMSrqPgRERERr6LiR0RERLyKih8RERHxKip+usCBAwcYN24cy5Yta962e/dubrjhBsaOHcu5557LokWLTnmeDz/8kIsvvpi0tDQuvfRSPvvss+6M3SFtXePq1au58sorGTduHDNmzODxxx+ntrb2pOeZMWMGqampLb5+9atfdXf8U2rr+h544IFWWadNm3bS83jSe3jjjTe2ur6mr+XLl5/wPO72Hh45cqTNa1iyZAnQO+7FU12jp9+Lp7q+3nAvnuwae8u9CLB8+fLmv/dLLrmEDz/8sHmfW92LhpwWh8Nh/OAHPzBSUlKMt99+2zAMwygpKTEmT55sPPjgg0ZmZqaxdOlSIy0tzVi6dOkJz7Nu3Tpj5MiRxmuvvWZkZmYaCxYsMEaNGmVkZmb21KWcUFvXuHHjRmP48OHGCy+8YGRnZxtr1qwxpk+fbvz6178+4XnsdruRmppqfPLJJ0ZhYWHzV0VFRU9dSpvauj7DMIzZs2cbf/nLX1pkLS4uPuF5PO09LC0tbXFthYWFxk9+8hPje9/7nmG329s8jzu+h6tWrTLS0tKMgoKCFplqamp6zb14smvsDffiya7PMHrHvXiya+wt9+Ly5cuN4cOHG4sXLzays7ON5557zhg2bJjx9ddfu929qOLnNP35z382brzxxhb/qCxcuNCYOnWqUV9f3+K4mTNnnvA8t956qzFv3rwW26655hrjd7/7XfcE74C2rvGXv/ylccstt7Q4bvny5caIESOMurq6Ns+zefNmIyUlxSgvL+/2zB3R1vU1NDQYaWlpxscff9zu83jae3i8d9991xgxYoSRkZFxwvO443v4/PPPG5dddlmb+3rLvXiya+wN9+LJrq+33Isnu8bjeeK96HK5jPPOO89YsGBBi+233nqrsXDhQre7F/XY6zRs3LiRN998k8cff7zF9k2bNjFx4kR8fb9bOm3KlCkcOHCA4uLiVudxuVx8/fXXTJkypcX2yZMns2nTpu4J304nusZbb72V++67r9XxDQ0NVFZWtnmuPXv2EBsbS1hYWLdk7YwTXV92djZ1dXUMHjy4XefxxPfwWNXV1fzpT3/ipptuIjU19YTHueN7uGfPHoYMGdLmvt5yL57sGnvDvXiy6+st9+LJrvFYnnov7t+/nyNHjnDppZe22L5o0SLuuOMOt7sXVfx0UkVFBffddx+//e1vSUxMbLEvPz+fhISEFtvi4uIAyM3NbfNc1dXVbb4mLy+vi5O338muccSIEQwbNqz5zw6Hg1deeYWRI0cSFRXV5vn27t1LUFAQd911F+eccw6XXXYZixcvxuVydet1nMjJrm/v3r1YLBZeffVVZsyYwQUXXMAf/vAH7Hb7Cc/lae/hsd544w2qqqr46U9/etLzudt72JSpuLiY6667jrPOOotrr72Wzz//HOg99+LJrrE33Isnu77eci+e7BqP5an3YnZ2NtBYvM2ZM4czzzyTH/7wh6xevRpwv3tRxU8nzZ8/n7Fjx7aqcgFqa2tbrSjv7+8PQF1dXZvHA22+pq3je8rJrvFYDQ0N3HfffWRmZvLQQw+d8Lh9+/Zht9u5+OKLWbRoEddccw1PP/00zz77bFdHb5eTXd++ffuwWq0kJSWxcOFC7r//ftasWcPcuXPb/HDx5PfQ6XTy2muvcd111xEaGnrS87nbe+hwOMjOzqayspJ58+bx4osvkpaWxu233866det6xb14qms8lifei6e6vt5wL7b3PfTke7GplfH+++9n1qxZvPzyy5x99tnMnTvXLe9F31MfIsdbvnw5mzZt4t13321zf0BAAA6Ho8W2pjcrKCio1fFN/wO09ZrAwMCuiNxhp7rGJk0384YNG3jmmWcYM2bMCY995ZVXqKurIyQkBIDU1FSqqqp4/vnnueuuu7Bae64WP9X13XXXXdx8883NTcopKSnExsZyzTXXsH379lbX6cnv4VdffUVubi5XX331Kc/pTu8hNH4wbty4EV9f3+YPyVGjRpGVlcWiRYt6xb14qms888wzAc+9F091fS+99JLH34vtfQ89+V708/MDYM6cOcyePRuA4cOHs2vXLl555RW3uxfV8tMJb7/9NsXFxZx77rmMGzeOcePGAfDQQw9xySWXkJCQQGFhYYvXNP05Pj6+1fkiIiIICgpq8zXHN/n1lFNdY1O+66+/ni1btvDSSy8xY8aMk57Tz8+v+UZtkpKSQnV1NeXl5d1zISdwquuzWCytnqWnpKQAjc23x/PU9xBg5cqVjB49mn79+p3ynO70HjYJCgpq9dNhSkoKBQUFveJehJNfI3j2vQgnv77ecC/Cqd9D8Ox7senvtum9aTJkyBAOHz7sdveiip9OePLJJ/nggw9Yvnx58xfA3XffzYsvvsjEiRPZvHkzTqez+TXr1q0jOTmZ6OjoVuezWCycccYZfPXVVy22b9iwgfHjx3frtZzIqa6xvLycm266iZKSEv71r3+16pR2PJfLxYwZM3j++edbbN++fTsxMTFERkZ216W06VTX98tf/pI5c+a0ygq02WnRE9/DJps3bz7l+wfu9x4CZGRkMG7cuFYdIHfs2MGQIUN6xb14qmv09HvxVNfXG+7FU11jE0++F0eMGEFwcDDbtm1rsX3v3r3079/f/e7F0x4vJoZhGC2GEBcVFRkTJ0407r//fmPfvn3G22+/baSlpRnLli1rPr6ioqLFPBWff/65MXz4cOPll182MjMzjccff9wYPXq0W8wR0+TYa7z//vuNkSNHGuvWrWs1P0VDQ4NhGK2vccGCBcYZZ5xhfPDBB0ZOTo7xxhtvGKNHjzbefPNNU67neMde3+rVq43U1FQjPT3dyMnJMT799FNjxowZxj333NN8vKe/h4bROIx45MiRxooVK9o83t3fQ6fTafzwhz80Zs2aZWzcuNHIzMw0Hn30UWPUqFFGRkZGr7gXT3WNnn4vnur6esO9eKprNAzPvxcNwzD+9re/GePGjTPeffddIycnx0hPTzeGDRtmrF+/3u3uRRU/XeT4f1S2bdtmXH311caoUaOM8847z3jttddaHH///fcb5513Xott//nPf4wLL7zQSEtLM2bPnm18+eWXPZK9vZqu0el0GmlpaUZKSkqbX4cOHTIMo/U11tfXG+np6cb5559vjBw50pg5c6bbFD6G0fo9/Oijj4wrrrjCGD16tHH22WcbCxYsMGpra5v3e/J72KSoqMhISUkxPvvsszaP94T3sLi42HjggQeMs88+20hLSzOuueYaY+PGjc37e8O9eKJr7C334qnew95wL57qGnvDvWgYhvHyyy8bM2bMMEaOHGlcdtllLeZncqd70WIYhnH67UciIiIinkF9fkRERMSrqPgRERERr6LiR0RERLyKih8RERHxKip+RERExKuo+BERERGvouJHREREvIqKHxHpNQ4fPkxqairLli0zO4qIuDEVPyIiIuJVVPyIiIiIV1HxIyLttnPnTm666SbGjx/PuHHjuPnmm1ut4rxkyRJ+8IMfMHbsWEaPHs3ll1/OBx980Lx/2bJlpKWlsXnzZq688krS0tKYOXMmq1evZv/+/dx0002MGTOGCy+8kPfff7/F61JTU9m2bRuzZ89m9OjRXHrppS3O3Zbc3FzuueceJk2axJgxY7jpppvYtWvXSV/z7LPP8r3vfY+VK1cya9Ys0tLSuPzyy9myZQtbt27lhz/8IaNHj2bWrFmsW7euxWv37t3LHXfcwRlnnMEZZ5zBz372Mw4dOtTimIyMDH7+858zZcoURo4cydSpU/njH/9IbW1t8zGpqam8/vrrPPjgg0yaNIlx48Zx9913U1RUdNLsInJqKn5EpF0qKyu57bbbiIyM5JlnnuGpp56ipqaGOXPmYLfbAXj99df5/e9/z/nnn88LL7zAE088gZ+fH/feey+5ubnN52poaOCee+7hRz/6Eenp6fj7+/OrX/2KO++8k3PPPZenn36a2NhY7r//fvLz81vkuOOOOzj//PN57rnnSE5O5p577mHVqlVtZi4pKeFHP/oRO3fu5He/+x1//vOfcblcXH/99WRlZZ30evPz83nssce48847+etf/0p5eTl3330399xzD1dffTV/+ctfcLlc/OIXv2guWg4cOMCPfvQjiouLWbBgAY888giHDh3i2muvpbi4GIDCwkKuv/56ampqWLBgAS+99BLf//73ee2111i8eHGLDE899RQul4u//OUv3HfffXz66ac8+uijHXrfRKQNXbI8qoj0elu2bDFSUlKMTZs2NW/LyckxHn/8cSM3N9cwDMN47LHHjD/96U8tXrdjxw4jJSXFePfddw3DMIy3337bSElJMf71r381H/Pee+8ZKSkpxl//+tfmbdu3bzdSUlKaV4Vuet2zzz7bfIzL5TIuv/xy4wc/+IFhGIZx6NChFivX/+UvfzHS0tKMw4cPN7+mrq7OOP/884277rrrhNf6zDPPGCkpKcaaNWuat73wwgtGSkqKsWTJkuZtH330kZGSkmLs2rXLMAzDuOeee4wzzzzTsNvtzceUlpYa48ePNxYsWGAYhmF8/vnnxvXXX9/iGMMwjFmzZhm33npr859TUlKMa6+9tsUxv/71r42xY8eeMLeItI+v2cWXiHiGoUOHEhUVxU9/+lO+//3vM336dM4880zuu+++5mN+/etfA2C328nOziY7O7v5sVB9fX2L840bN6759zExMQCMHTu2eVtERAQAFRUVLV53+eWXN//eYrFw4YUX8uyzz1JTU9Mq87p16xg+fDjx8fE0NDQAYLVamTZtGitWrDjlNZ9xxhkdyrh+/XomT55MQEBA8/cLCQlhwoQJfPnllwCcc845nHPOOdTX13PgwAGys7PZs2cPJSUlzedrcuz3AkhISGjzOkWkY1T8iEi7BAcH8/rrr/P888/zwQcf8MYbbxAYGMhll13Ggw8+iL+/PwcPHuT3v/8969evx9fXl0GDBpGamgqAYRgtzhcSEtLqewQEBJwyR3x8fIs/R0dHYxhG86O3Y5WVlZGTk8PIkSPbPFdNTQ2BgYEn/F4dzVhWVsYHH3zQZj+kqKgogObHWK+//jrV1dUkJiYyevRo/P39W73m+GxWq7XV36OIdJyKHxFpt0GDBvHEE0/gdDr55ptveOedd/j3v/9N3759ue222/jJT36Cn58fb731FiNGjMDX15fMzMx2tbK0V2lpaYsCqKioCB8fHyIiIigsLGxxbGhoKJMmTWrROnUsm83WZbmavt9ZZ53FLbfc0mqfr2/jx+2LL77I4sWLmT9/PjNnziQ0NBSAq666qkuziMiJqcOziLTLRx99xJQpUzh69Cg+Pj6MGzeO+fPnExYWRn5+PqWlpRw4cICrrrqK0aNHN/9j/9lnnwGNLR5dYfXq1c2/NwyD//3vf4wfP77NQmbSpEkcOHCA5ORk0tLSmr9WrFjBkiVL8PHx6ZJMx36/zMxMhg8f3vy9Ro0axeLFi/n4448B2Lx5M0OGDOGqq65qLnwKCgrYu3dvl/0dicjJqeVHRNrljDPOwOVy8bOf/Yyf/OQnBAcH8+GHH2K327nooouIjo4mKSmJ119/nYSEBMLCwvjiiy949dVXAbqsr8oTTzyBw+EgOTmZJUuWkJWV1fw9jnfzzTfzzjvvcPPNN3PrrbcSGRnJBx98wFtvvcUDDzzQJXmONXfuXH70ox9xxx13cO211+Lv78+bb77JypUreeaZZwAYPXo06enpvPjii4wdO5acnBxeeOEFHA6H+vOI9BAVPyLSLnFxcfz973/n6aef5sEHH6SmpoahQ4fy7LPPMmXKFADS09N55JFH+PWvf43NZmPIkCE8//zzPProo2zatIkbb7zxtHPMnz+fF154gUOHDjFixAhefvllJkyY0Oax8fHxvPHGG/z5z39m/vz51NXVMXDgQB555JFuecw0bNgwXn/9dZ566inuu+8+DMMgJSWFv/3tb5x//vlA41D90tJS/vGPf/C3v/2NxMRELr/8ciwWCy+88ALl5eWEh4d3eTYR+Y7FUO85EfEAy5Yt44EHHmDVqlX07dvX7Dgi4sHU50dERES8ioofERER8Sp67CUiIiJeRS0/IiIi4lVU/IiIiIhXUfEjIiIiXkXFj4iIiHgVFT8iIiLiVVT8iIiIiFdR8SMiIiJeRcWPiIiIeBUVPyIiIuJV/j/1UG+NZT/gyAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# first simulate the sampling distribution of the mean for 10,000 samples\n", "nSamples = 10000\n", "samplesize = 100\n", "m=np.empty(nSamples) # make an array to store the means\n", "binwidth=0.1\n", "\n", "for i in range(nSamples):\n", " sample = UKBrexdex.sample(n=samplesize, replace=False)\n", " m[i]=sample.score.mean()\n", " \n", "sns.histplot(m, bins=np.arange(40,60,binwidth))\n", " \n", "# now work out the expected normal distribution\n", "mu = UKBrexdex.score.mean()\n", "sigma = UKBrexdex.score.std()\n", "SEM = sigma/(samplesize**0.5)\n", "\n", "x = np.arange(40,60,binwidth) # x axis values - you may wat to change these once you have tried plotting it\n", "p = stats.norm.pdf(x,mu,SEM)\n", "freq = p*nSamples*binwidth # exected frequency in each ibn is probability of the bin * total nSamples * bin width\n", "\n", "plt.plot(x,freq,'r')\n", "plt.xlabel('sample mean')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "8c05c451", "metadata": {}, "source": [ "This is clearly quite a good fit.\n", "\n", "Now try changing the `samplesize`, $n$, in the code block above, to $n=4$ (hint: change the histogram bins and x-axis values to `np.arange(0,100,binwidth)`, and change `binwidth` to 1)\n", "\n", "Hopefully, you can see that although the histogram on its own looked quite normal, it is actually not a great fit to the normal distribution we would expect if the Central Limit Theorem applied - the peak is too flat and there are fewer sample means out in the tails than we would expect - the distribution looks like a piece of Toblerone\n", "\n", "\n", "\n", "### Q-Q plot\n", "\n", "The differences in the peak and tails of the distribution can be hard to spot on a histogram/Normal plot as above.\n", "\n", "A type of plot designed to make these clearer exists - it is called a Q-Q (quantile-quantile) plot. In the Q-Q plot, we plot the quantiles of the data distribution (in this case our 10,000 simulated sample means) against the quantiles of the normal distribution.\n", "\n", "If our data distribution was normal, the points would all fall on a straight line, but here we see the deviation at the tails of the distribution, reflecting the difference between the triangular tails of the simulated sampling distribution as opposed to the finely tapered tails of the normal distribution." ] }, { "cell_type": "code", "execution_count": 86, "id": "499c1583", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAHJCAYAAACG+j24AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAACEe0lEQVR4nO3dd3RUxd/H8feW9N4DhF4ChN5baIIoTcVekCYoHZWmYEHFSlG6AiKoiApIRxRRaYEkFAUpofcU0nuy5fmDh/25bnZTyKbsfl/neDBzZyczu0A+3Jk7o9Dr9XqEEEIIIeyUsrw7IIQQQghRniQMCSGEEMKuSRgSQgghhF2TMCSEEEIIuyZhSAghhBB2TcKQEEIIIeyahCEhhBBC2DUJQ0IIIYSwaxKGhBCiGCrLPrWVpZ9CVAQShoQQRTJ48GBCQ0ON/mvSpAndu3dn1qxZpKamlsr36dmzJ9OnT7/ndjZu3EhoaCjXr183W+f69euEhoaycePGAl8zffp0evbsaaj/22+/MW3atHvq193v+e//GjZsSMuWLRk0aBDr168327+iWrp0KStXrrynfgphT9Tl3QEhROXRuHFj3nrrLcPX+fn5/PPPP8ybN4/Tp0/z3XffoVAoyrGHxRMYGMj3339PjRo1Crw+ZswYnn/+ecPXX331Val979GjR9O9e3fgzl2czMxMfvzxR2bMmIFGo+Gpp54qcduffvop48aNK6WeCmH7JAwJIYrM3d2dFi1aGJW1bduWzMxMFixYwF9//WVyvSJzdHS02F9zIak01KhRw+R7d+rUiTNnzvDVV1/dUxgSQhSPTJMJIe5ZkyZNALh58yZwZ0pt8uTJTJgwgVatWjFq1CgA0tPT+eCDD+jVqxdNmzalf//+RtNCd+Xn5/Pee+/Rtm1b2rZty7Rp00hKSjKq8+OPPzJo0CBatGhBs2bNeOihh9ixY4dJW0ePHuXhhx+madOmDBgwwKhOYdNQ/54mGzx4MJGRkURGRhIaGsrBgwfp0qULr776qsnrHnzwQV577bWivHVGlEoljRo1MryPBbl8+TITJkygc+fOtGjRgsGDB3PkyBHD9dDQUAAWLVpk+H8hhGUShoQQ9+zSpUsAVK9e3VC2c+dOHBwcWLx4Mc8//zw5OTk888wzbNmyheHDh7NkyRJat27NjBkzWLZsmVF7O3fu5OTJk3z44YdMnTqVP/74gzFjxhiuf/vtt7z55pvcd999fP7553zyySc4ODgwZcoUkyDxxhtv8MADD7B48WLq1avHyy+/zP79+4s9xrfeeovGjRvTuHFjvv/+e5o1a8bDDz/M7t27ycjIMNT766+/uHjxIoMGDSr294A776W5O1Lnz59n0KBBXLt2jZkzZzJnzhwUCgVDhgwhMjISgO+//x6Axx57zPD/QgjLZJpMCFFker0ejUZj+Do1NZXIyEiWLl1KixYtDHeI4M5djnfffRdXV1cA1q5dS0xMDGvXrqV169YAhIeHo9FoWLJkCU899RTe3t4AeHp6smLFCtzd3QHw8fFh7Nix7N+/ny5dunDt2jWGDx/O2LFjDd8vJCSEQYMGcfToUapWrWooHzt2rOHOVNeuXbl8+TKLFi2iS5cuxRp7vXr1DP25O7316KOPsnz5cnbt2sWjjz4KwE8//USNGjVo06aNxfZ0Op3hvdTpdMTFxfH1119z5swZo3VZ/7Zo0SIcHBxYs2YNHh4eAHTv3p3+/fvzySef8OOPPxr6FhwcXKmmLIUoTxKGhBBFFhUVRVhYmFGZUqmkY8eOvPvuu0aLp0NCQgxBCCAyMpJq1aoZgtBdAwcOZP369fz1119069YNgG7duhmCB9x5wszBwcEwNXX3abP09HQuX77M5cuXiYiIAO5Msf3bgw8+aPR1r169WLhwIZmZmSV9Gwxq165N69at2bx5M48++ih5eXns2LGDIUOGFLqQfMaMGcyYMcOozN3dnZdeesnseqHIyEh69OhhCEIAarWafv36sXjxYjIzM3Fzc7vncQlhbyQMCSGKLCwsjFmzZgGgUChwcnKiSpUqRsHlLn9/f6OvU1NTTcr+XS8tLc3sa5VKJd7e3oY6V69e5c033+TQoUOo1Wrq1KljWB/z3/11AgICjL728/NDr9cbTW3di8cee4zXX3+dmzdv8tdff5GWlsYjjzxS6OvGjRtneJpMqVTi4eFBtWrVUKlUZl9j6T28OyYJQ0IUn4QhIUSRubm50bRp0xK91svLiytXrpiUJyQkAHemwu76dzAC0Gq1JCcn4+fnh06nY9SoUTg4OPDDDz/QuHFj1Go158+fZ8uWLSbtp6am4uzsbPj69u3bqFQqvLy8uH37donG8m8PPPAA7733Hrt27eLYsWN07NjRaJrOnGrVqhX7vTTX54LeQyFE0ckCaiFEmWjbti03btwwevIJYMuWLTg4ONCsWTND2cGDB43WJu3atQuNRkP79u1JTk7m0qVLPPbYYzRr1gy1+s6/6fbu3QvcWX/zb/v27TP8v06n4+eff6Z58+ZGAamolErTvzJdXV3p27cv27ZtY9++fUW6K1RSbdu25ffffyc9Pd1QptVq2b59O02bNsXR0dFsP4UQ5smdISFEmRg0aBBr165l3LhxTJgwgerVq7Nnzx42bNjAuHHj8PT0NNS9ffs248ePZ/DgwVy+fJl58+bRuXNnOnbsiEKhoFq1anz77bcEBwfj6enJ/v37Wb16NQDZ2dlG3/fTTz9Fq9VSpUoVvvvuOy5dusSqVatKNAZPT0+OHTtGREQEjRs3xsvLC7gzVfbkk0/i7u7O/fffX8J3qHDjxo1j7969PP/884waNQpHR0e++eYbrl27xooVK0z6GRUVRZs2bSrVRphClAf554MQoky4uLjw9ddf07NnTxYsWMDo0aM5cuQIs2fPZvz48UZ1n3jiCfz9/Rk7diyfffYZAwYMYNGiRYYf6kuWLCEoKIjp06czadIkjh8/ztKlS6lTpw7R0dFGbc2ePZs1a9YwZswY4uLiWL58Oe3atSvRGJ599lkcHBwYOXKk4U4U3Hm6zMfHh379+pXojlNR1a9fn7Vr1+Lv78/rr7/OlClT0Ov1rFmzhk6dOhnqvfTSS5w4cYKRI0dy69Ytq/VHCFuh0MtpfkIIcU/+/vtvHn/8cTZs2GC0vYAQonKQaTIhhCihw4cPc/jwYTZt2kSHDh0kCAlRSck0mRBClFBycjKrVq3Cz8+PDz74oLy7I4QoIZkmE0IIIYRdkztDQgghhLBrEoaEEEIIYdckDAkhhBDCrkkYEkIIIYRdk0fri0iv16PTFW2tuVKpKHLdyk7GapvsaaxgX+OVsdomGav5ukXZgV3CUBHpdHqSkjILradWK/HxcSMtLQuNRldo/cpMxmqb7GmsYF/jlbHaJhmreb6+bqhUhYchmSYTQgghhF2TMCSEEEIIuyZhSAghhBB2TcKQEEIIIeyahCEhhBBC2DUJQ0IIIYSwaxKGhBBCCGHXJAwJIYQQwq5JGBJCCCGEXasQYWjTpk307duXpk2b0q9fP3bu3Gm4dvr0aZ577jlatGhB9+7dWblyZaHt7dy509DegAED2Lt3rzW7L4QQQogS0On0nLmSzKFTsZy5klxuR4qU+3Ecmzdv5vXXX2fatGl0796dbdu28corrxAcHEytWrUYNmwYvXr1YtasWRw/fpxZs2bh7e3No48+WmB7hw4dYsqUKUyfPp2OHTuyfv16xo4dy6ZNm6hbt24Zj04IIYQQBTlyNp61u8+RnJ5rKPPxcOKZXvVpHRpYpn0p1ztDer2ezz77jCFDhjBkyBBq1qzJ2LFj6dSpE5GRkfzwww84Ojry9ttvU7duXR599FGGDh3K8uXLzba5fPlyevfuzXPPPUfdunWZNm0aYWFhrF69ugxHJoQQQghzjpyNZ/FPJ42CEEByei6LfzrJkbPxZdqfcg1DFy9e5MaNGwwYMMCofOXKlbz44otER0fTtm1b1Or/3cDq0KEDly5dIjEx0aQ9nU7H0aNH6dChg1F5+/btiY6Ots4ghBBCCFFkOp2etbvPWazz3e5zZTplVq7TZJcvXwYgKyuLESNGcOrUKUJCQhg9ejQ9e/YkNjaWBg0aGL0mMPDOrbObN2/i5+dndC0tLY2srCyCg4NNXnPr1q177q9aXXh2VKmURr/aMhmrbbKnsYJ9jVfGapsq21hPX04yuSP0X0npuVy4mUqjWr5G5dYaa7mGoYyMDACmTZvGuHHjmDx5Mrt27WLMmDGsWrWKnJwcHB0djV7j5OQEQG6u6RuZk5MDUOBrCqpfHEqlAh8ftyLX9/R0uafvV5nIWG2TPY0V7Gu8MlbbVFnGmn8puWj19OZ/7pb2WMs1DDk4OAAwYsQIHnnkEQAaNWrEqVOnWLVqFc7OzuTl5Rm95m6ocXV1NWnvblAq6DUuLvf2xul0etLSsgqtp1Ip8fR0IS0tG61Wd0/fs6KTsdomexor2Nd4Zay2qbKN1UFRtOkvB4We5ORMo7LijtXT06VId5HKNQzdnc7671RYvXr1+OOPP6hWrRrx8caLqO5+HRQUZNKet7c3rq6uBb7mv1NnJaHRFP03mVarK1b9ykzGapvsaaxgX+OVsdqmyjLWulW98PFwsjhV5uvhRN2qXmbHU9pjLdcJxsaNG+Pm5sZff/1lVB4TE0ONGjVo27YtR44cQavVGq5FRERQu3Ztk/VCAAqFglatWhEZGWlUfvjwYVq3bm2dQQghhBCiyJRKBc/0qm+xztO96qNUKsqoR+UchpydnXnhhRdYvHgx27Zt4+rVqyxdupQDBw4wbNgwHn30UTIyMpgxYwbnz59n48aNrF69mhdffNHQRnp6OklJSYavhw0bxvbt21m1ahUXLlzg448/5vTp0wwZMqQ8hiiEEEKI/2gdGsjYR5rg4+FkVO7r4cTYR5qU+T5D5b7p4pgxY3BxcWH+/PnExcVRt25dFi5cSPv27QFYsWIFs2fP5pFHHiEgIICpU6ca1hcBzJ49m8jISPbs2QNAly5deP/991myZAnz58+nXr16LFu2TDZcFEIIISqQ1qGBtKwfQMy1FFIyc/F2c6JBde8yvSN0l0Kv15fP3teVjFarIykps9B6arUSHx83kpMzK8Xc7b2Qsdomexor2Nd4Zay2ScZqnq+vW5EWUFeOTQmEEEIIIaxEwpAQQggh7JqEISGEEEKUC21mJvHrviXp5x3l2o9yX0AthBBCCPuTfS6GW8s/R5OUiNrfH98H+pZbXyQMCSGEEKLM6LVakrZvJXHrZtDrcQgIpMqLY8q1TxKGhBBCCFEm8hMTiV3xOdnnYgDw7NiZwGefQ+lcvueqSRgSQgghhNWlR0cRt2YVuqwslM7OBA4egmf7juXdLUDCkBBCCCGsSJebS/y6b0nbtxcA59p1CB71Eo4BZbvLtCUShoQQQghhFTlXrxD7xTLyYm+BQoHvg/3wG/gwCnXFih8VqzdCCCGEqPT0ej0pu3/h9oYf0Ws0qLy9qTJiFK6NGpd31wokYUgIIYQQpUaTlkbslyvIOvk3AG4tWhI8ZDgqD49y7pl5EoaEEEIIUSoy/zlJ7Mov0KaloVCrCXjyaby690ShKPvDV4tDwpAQQggh7oleo+H2xvUk//IzAI5Vq1Fl1Es4hVQv554VjYQhIYQQQpRYXmwst75YSu7VKwB4de9JwBNPoXR0LOeeFZ2EISGEEEIUm16vJ+3gfuLXfoM+NxelmxvBQ4fj3rJ1eXet2CQMCSGEEKJYtFlZxH+zmvTIwwC4hDYkaPhILmYoSTkVi7ebEw2qe6NUVuy1QndJGBJCCCFEkWVfOM+t5cvQ3L4NSiV+Dz3C5Tpt+XTdGZLTcw31fDyceKZXfVqHVpzNFc2RMCSEEEKIQul1OpJ2bCNxyybQ6XDwDyB45Iuc0niy+KeTJvWT03NZ/NNJxj7SpMIHIglDQgghhLAoPymJWys+JyfmLAD6sFaEjByBytWVtUsPWnztd7vP0bJ+QIWeMpMwJIQQQgiz0o8e4eaXK1DkZJOnUPNLQHtO5tTBZ/VxujWvYjQ1VpCk9FxirqXQsKZPGfW4+CQMCSGEEMKELjeXhB/Wkfrn7yiAW05+bAkKJ9nRE7gzDbZp/+UitZWSaTkwlTcJQ0IIIYQwknv9Gre+WErezZsAHPIOY69fC3QKVYna83ZzKs3ulToJQ0IIIYQA/v+A1d9/4/YP69BrNOjdPfjesz2XXauWuE1fjzuP2VdkEoaEEEIIgTY9ndivVpL513EA3Jo243rngVzeffWe2n26V/0KvXgaJAwJIYQQdi/r9ClurfgCbWoKCrUa/8eewPu+3qRfTQEKD0MPd6nNn3/dNFpM7evhxNOyz5AQQgghKjK9RsPtzT+R/PMO0OtxDK5C8KiXcAypwdmrKSRl5ODh4kB6dr7ZNnw9nOjfqRb9O9Ui5loKKZm5sgO1EEIIISomnU5vCCyeOWmoNn2D4uadOz+6lh2I7dKXv69q+XPrwUIfmb/r39NgFfnxeUskDAkhhBA26m74ScrI4fSlZI6fv01mjoaw9IvcH38YtT6fHKUjOwM7cja9Juw8X+S2K9M0WGEkDAkhhBA2RqfTs+3gZX6NvkZmjsZQ7qjLo39CJE3SLwJwzTmQrUFdSHNwL1K7Hi4OPHVffXw8Ktc0WGEkDAkhhBA25MjZeL7aecYoBAFUybnNwNi9+Ggy0KHggG8zDvo0Ra9QFrnt9Ox8fDycKu10mDkShoQQQggbEXUmnqWb/nNoql5Ph5SThCceR4WeVLUbW4O6cN0lqETfo6LvJl0S5R6Gbty4Qc+ePU3K33vvPbZs2UJkZGSBr/voo494+OGHC7zWs2dPbty4YVQ2YMAA5syZc8/9FUIIISoSnU7PifO32RN1hd3R142uuWuy6B+3n1rZsQCcdq/JzwEdyFWVfEfoir6bdEmUexg6e/YsTk5O7N69G4Xif3OPHh4e9O7dm/x848f5Zs6cydWrV+nVq1eB7WVkZHDz5k0+//xzwsLCDOXOzs7WGYAQQghRTqLOxPPNL2dJzzJ99L1e5jX6xh3EVZdLnkLN7oC2/O1RDxQlX+dTGXaTLolyD0MxMTHUrl2bwEDT1ej/DTDbtm1j//79bNy4EXf3ghd7xcTEoNfradWqFZ6enlbpsxBCCFFedDo9Z64ks3HfRS7eTDO5rtJp6ZkYTevUswDEOfqwObgrSY5e9/y9K8Nu0iVR7mHo7Nmz1KtXr9B6WVlZfPzxxwwZMoTQ0FCL7QUEBFglCKnVhS8yU6mURr/aMhmrbbKnsYJ9jVfGWnnpdHpOX07ityPXOX7uNhqdvsB6fnkpPBS7l8C8FAAivRrxp38rtCU8YPUuX08nnr0/lLYNy/cxemt9rgq9Xl/wO1pG+vXrR0BAAHl5eVy+fJmaNWsyZswYwsPDjep9+eWXLF68mD/++AMPDw+z7c2aNYsDBw4QGhrKsWPH8PX1ZdCgQTz//PMolSV/8/R6vdE0nhBCCGFNWp2ek+dvs/3gJaJOxaHR6sxX1utpkRbDfbejcdBryVQ5sz2wMxfdqpXoe/t5OdOnQ02q+rvj6+lM4zp+qGzwjtBd5Xpn6G4AcnFxYerUqbi6urJlyxZGjhzJqlWr6NixIwBarZavv/6aZ555xmIQAjh37hzp6en07duXcePGER0dzZw5c0hNTWXixIkl7qtOpyctLavQeiqVEk9PF9LSstFa+o1rA2Sstsmexgr2NV4Za+Vx+J9Ylm89RZ6m8L47a3N5MD6C0Mw7O0lfdKnC9qAuZKpdzL7GzUVN7zbVCa3hQ1pmHp6uDqBQkJaZh7e7I6E1fIymw9JSC//5VxaK+7l6eroU6S5SuYYhR0dHoqKiUKvVODo6AtCkSRMuXLjAypUrDWEoMjKSmzdv8sQTTxTa5qpVq8jNzTWsKQoNDSUzM5OlS5cyfvz4e7o7pCnCb8q7tFpdsepXZjJW22RPYwX7Gq+MtWK6uxZoza4zxKfkFOk11bNjGRC3H09NFlqU/OnXkkjvxgUuku7UJIiwWn5F2jBRp9OjMzMVVxGU9uda7muGXF1dTcoaNGjA/v37DV/v3r2bZs2aUb169ULbc3BwwMHBwaS9rKwsUlNT8fGxrY2ihBBCVH7mNko0R6HX0SXpLzoln0ABJDl4sDmoK3HOfgXWf/GhMNo3Ktm+QvagXFeWnTlzhpYtWxIdHW1UfvLkSaNF1UeOHKFDhw6FtqfT6ejZsydLly41Kj9x4gT+/v4ShIQQQlQ4kafiWPzTySIHIa/8DJ69sYvO/x+E/vaoy6rq/c0GoQfaVZcgVIhyvTPUoEED6tevz6xZs3jrrbfw8fHhhx9+4Pjx46xfvx64s17o/PnzjBgxosA20tPTyc/Px9fXF6VSSZ8+fVixYgW1atUiLCyMiIgIVqxYwYwZM8pyaEIIIUSh1v0Wwy9R1wuv+P8apV+iT8IhnHX55Cgd2BXQgdMetQus6+HiwHN9GtC2oQShwpRrGFIqlSxbtow5c+YwadIk0tLSaNy4MatWrTI8Pp+SkkJ+fj7e3t4FtjF79mwiIyPZs2cPAK+++iqenp7MnTuX2NhYQkJCmDFjRpHWGwkhhBBl5bMf/+KvC4lFquugy6d3QiTN0i8AcN05gK1BXUh1MH6oyEGtpFuLqrSqH2BTB6laW7k/Wl9ZaLU6kpIyC62nVivx8XEjOTmz0izaKykZq22yp7GCfY1XxlpxzP/+GCcuJRepblBOIg/F7cU3Px09cNCnKft9mxsdsOqgVtC/U236dahp0wGouJ+rr69bxX+aTAghhLA3b395mKvxhf/jGr2edimn6JZ4DBU60tSubA3qwjWXYEMVB5WCfp1rMXRAU9JSsypk8KsMJAwJIYQQZeStlYe4llD4nj1ummz6xe2nTvYtAM661WBnYEdyVE6olNC8nj89W4XQsIYPjo4qm94QsSxIGBJCCCGsSKfTc+pSEks2nyAnr/A7N3Uyb9Av/gBu2hzyFSp+82/Lcc/6BPq4MO6BhjT8z4aI4t5JGBJCCCGsQKfTs3nfRbYfukJR9i9U6bV0v32UtqmnAYh39GZLUFeSXXx4cUAj2jcKLqQFUVIShoQQQohSduRsPJ9v+QeNtmjPKPnmpfJQ7F6C8u4sqj7iFcoevza0ahTEhw81lTtBViZhSAghhChFUWfiWbrpZNEq6/U0Sz9Pr4QoHPUaspRO7AjqxCWPGowa2Fj2CCojEoaEEEKIUhJ1Jo6lm/4pUl0nbR4PJETQKOMKAJddgtkW1IWAkEA+f76t3A0qQxKGhBBCiFJw5Gx8kYNQtex4Bsbtw0uTiRYFe/1actg7jBpB7rwxtJ2Veyr+S8KQEEIIcY90Oj3f/HK20HoKvY5OySfonPQ3SvQkq93ZEtyVW87+VA9w5e3h7cugt+K/JAwJIYQQ9yjmWgqpmfkW63jmZzAgbj/Vc+IBOOlRh18C2pGndKR6gCuzRhR+ILmwDglDQgghxD3QaHTM/+G4xTqhGVd4MD4CZ10euQo1vwR24B+POgDUCHTl7eEShMqThCEhhBCihL7fc45dkdfMXnfQ5XPf7WhapJ0D4KaTP1uCw0n5/wNW729Tjad6hZZJX4V5EoaEEEKIEliw/i+Onzd/6nxgbhIDY/fhn5+KHjjk3YR9fi3QKZR0bBLEsAcaoVYXfoiosD4JQ0IIIUQxRZyKNR+E9HrapJ6h++0jqNGRrnJhW1AXrrhWwc1JyWcTu8lj8xWMhCEhhBCiGA6cuMXK7acLvOaqyaZv/EHqZd0A4JxrCDuCOpGtcgZg/viuEoQqIAlDQgghRBFoNDomLdpPVo6mwOu1sm7SP+4A7tpsNAole/zacNQrFBR3wk/vNiEyLVZBSRgSQgghzNBodPwceZltB6+Qpyn4nDGlXkvXxON0SLmz4WKCoxdbgrqS4ORjqOPl6sjTvRqUSZ9F8UkYEkIIIQrww55z/GzhSTEAn7w0Bsbto0runfVDRz0bsMe/DRql8Y/XueM6W62f4t5JGBJCCCH+o7BH5tHraZJ+kfsTDuOo15CtdGRnYCdi3GuYVB37SBNZJ1TBSRgSQggh/iXiVKzFIOSkzeP+hEOEZVwG4KpzEFuDu5CudjOqpwDGPNKE1qGBVuytKA0ShoQQQgggJ0fDGysjSEw3f6xG1ZwEBsbuw1uTgQ4F+32bE+HTBL3CeGF0rSA3Zg5pJ3eEKgkJQ0IIIexaXp6WyUv2k5GjNVtHodfRIfkfwpOOo0RPitqdLUHh3HQJMKnr4aLmzWFy4GplImFICCGEXdLp9Ly3JorLsRkW67lrshgQt4+a2XEAnHKvxa6ADuSqHE3rOqv4bGJXq/RXWI+EISGEEHYn6kw8SzedLLRe/Yyr9I0/iIsujzyFml8D2nHCo65h76B/69WqCs/c38ga3RVWJmFICCGEXfludwy/Rl+3WEet09DzdjSt0mIAiHXyZXNQV5IdPU3qOihh8SvdZUPFSkzCkBBCCLug0eiYvHQ/aZkF7yB9l39uMg/F7SMgLwWAw96N+dOvJTqFyqSu3A2yDRKGhBBC2Lx1u8/xS7TlDRTR62mVepaeidGo9ToyVM5sD+zMJbdqJlXbNw5gRN8wuRtkIyQMCSGEsFl5eVomLNhLloUnxQBctDn0jT9I/cw702cXXKuxPbATWWoXk7pjZe8gmyNhSAghhE2ateIQ0afjCq1XIyuWAXH78NBmo0HJH/6tiPZqZLJIWgEsn9pD9g6yQRKGhBBC2Jxx8/4kLcv85okASr2OLknH6Zh8EgWQ6ODJluBw4pz8TOqG+DvzzgudrNRbUd4kDAkhhLAZGo2OUXP+KLSeV346A2P3US33NgDHPevxm39b8pUORvWqB7ry2jNtcHaWH5e2rNw/3Rs3btCzZ0+T8vfee4/HH3+c1157jY0bNxpdCwoKYu/evWbb3LlzJwsXLuTatWvUqlWLKVOm0LWrbIIlhBC27Oufz/L78RuF1mucfpE+8Ydx0ueTo3RgZ2BHzrrXMqpTM9CFN4Z2kCkxO1HuYejs2bM4OTmxe/duFP+an/Xw8DBcf+mll3juuecM11Qq08cb7zp06BBTpkxh+vTpdOzYkfXr1zN27Fg2bdpE3bp1rTcQIYQQ5Wbkx3vQ6izXcdTl0zvhME3TLwJwzTmArUHhpDm4G9Vb9ko3HB3N/5wRtqfcw1BMTAy1a9cmMNB0Zb5Wq+X8+fOMGTOGgADT818Ksnz5cnr37m0IT9OmTePYsWOsXr2ad955p1T7LoQQovwN/3BPoXWCc24zMG4fvvnp6FBw0KcpB3ybmRyw+uV005kKYfvKPQydPXuWevXqFXjt8uXL5ObmFvmOjk6n4+jRo0yfPt2ovH379vz666/33Nei7CehUimNfrVlMlbbZE9jBfsary2O9fn3dluuoNfTPuUfuiYeQ4WeVLUrW4PCue4SZFStmp8zH4zuYsWeWo8tfq7mWGus5R6GYmJiCAgI4JlnnuHy5cvUrFmTMWPGEB4eTkxMDAqFgtWrV7N3716USiXdunVj0qRJhmm0f0tLSyMrK4vg4GCj8sDAQG7dunVP/VQqFfj4uBW5vqen6d4UtkrGapvsaaxgX+O1lbEOeHWzxetumiz6xx2gdvadv//PuNVkZ2AHclVORvXWze6Hmw0skLaVz7UoSnus5frp5+XlcfnyZVxcXJg6dSqurq5s2bKFkSNHsmrVKs6dO4dSqaRatWosW7aMK1eu8NFHHxETE8Pq1atRKo2TYU5ODgCOjsYnCTs5OZGbm3tPfdXp9KSlZRVaT6VS4unpQlpaNtrCJrArORmrbbKnsYJ9jddWxpqXp+WFj3+3WKdu5nX6xR3AVZdLvkLFr/7t+NuzntHeQS4OCj6fdh952bnkZd/bz4jyZCufa1EUd6yeni5FuotUrmHI0dGRqKgo1Gq1IcA0adKECxcusHLlSpYvX87QoUPx9LxzMF6DBg0ICAjgySef5MSJEzRv3tyoPSenO2k/Ly/PqDw3NxcXl3tPkRpN0X+TabW6YtWvzGSstsmexgr2Nd7KPNb53x/nxKUks9dVOi09Eo/QJvUMAHGOPmwJDifR0duo3oJxXXB3d6y070NBKvPnWlylPdZyvy/o6upqUtagQQP279+PQqEwBKF/XwOIjY01CUPe3t64uroSHx9vVB4fH28ydSaEEKJyeWnOH+RZ+AHol5fCQ7H7CMxLBiDKqxF/+LVCqzR+MkwWSYv/KtfVVmfOnKFly5ZER0cblZ88eZJ69erx6quvMmLECKNrJ06cAChw0bVCoaBVq1ZERkYalR8+fJjWrVuXcu+FEEKUlfHzLQQhvZ7mqTEMvbadwLxkMlXO/FClJ78FtDUKQo1reEgQEgUq1ztDDRo0oH79+syaNYu33noLHx8ffvjhB44fP8769eu5efMmo0ePZunSpfTr149Lly7xzjvv0L9/f8MTZunp6eTn5+Pr6wvAsGHDGDVqFI0bN6Zr165s2LCB06dPM3v27PIcqhBCiBJ6+bM/ycwtOAg5a3N5MD6C0MyrAFxyqcK2oM5kqo1nHWTvIGGJQq/X68uzA0lJScyZM4e9e/eSlpZG48aNmTx5Mm3atAFg165dLFu2jIsXL+Lh4cGAAQOYNGmSYX3Q9OnTiYyMZM+e/+0zsWnTJpYsWUJsbCz16tVjypQpdOzY8Z76qdXqSErKLLSeWq3Ex8eN5ORMm5+7lbHaJnsaK9jXeCvjWF/4aA86Mz+lqmfHMSBuH56aLLQo+dOvJZHejY0WSTeo6sb059uXUW/LR2X8XEuquGP19XUr0gLqcg9DlYWEIVMyVttkT2MF+xpvZRuruc0UFXodXZL+pmPyCZToSXLwYEtQOLHO/kb1nB2VfDG1Z6UY672obJ/rvbBWGCr3BdRCCCHEv2Vk5DFh0f4Cr3nlZzAgbh8hOQkAnPCoy68B7cj7zwGrCuDHDwaQnFz4P2KFkDAkhBCiwpjw6V4ycjQFXmuYfpkHEiJw1t05YHVXQAdOe9Q2qVfV15EPx8jh3KLoJAwJIYSoEMxNizno8umVEEXz9PMA3HDyZ0twOKkOpicRLJnUFWcb2E1alC35HSOEEKLcmQtCQTmJDIzbh19+Gnr4/wNWm6NTmK4DkcfmRUlJGBJCCFFu8vK0vDTvT9MLej1tU07TPfEoKnSkqVzZFtSFq64Fb6ArQUjcCwlDQgghypxGo2PGFwdISMs3ueaqyaZf/AHqZt0EIMatOjsCO5Kjci6wLQlC4l5JGBJCCFGmvt9zjl2R1wq8VjvzBv3jD+CmzSFfoWKPfxuOeTYw2jvo3yQIidIgYUgIIUSZ+eTbo5y+lmJSrtJr6ZZ4jHYppwCId/RmS1A4t518zLYlQUiUFglDQgghyoS5RdK+eakMjNtHcO6d0+iPeIXyu19rNErzP6IkCInSJGFICCGE1RUYhPR6mqWfp1dCFI56DVlKJ3YEdeK8W3Wz7TSq7s6UZ9tZsafCHkkYEkIIYTU3b6Yzc02USbmTNo8HEiJolHEFgCsuwWwN6kLGfw5Y/Tc5bFVYi4QhIYQQVmFuWqxadjwD4/bhpclEh4K9vi047BOGvoC9g+6SaTFhTRKGhBBClCqNRseoOX+YlCv0Ojomn6RL0l8o0ZOidmdzcDi3nAPMtuXlqmT+hO7W66wQSBgSQghRilZvP82fJ26ZlHvkZzIgbj81cuIAOOlem18C25OndDTb1qIJ4bi6Opi9LkRpkTAkhBCiVJibFmuQcYUH4yNw0eWRq1DzS0B7/vGsa7EtmRYTZUnCkBBCiHty9Woqb689YlKu1mm473YULdPOAXDTyY8tQeGkOHqabeud59oQEmL+uhDWIGFICCFEiZm7GxSYm8TA2H3456eiBw55h7HPrwU6hfmnweRukCgvEoaEEEKUiLm9g1qnnqFH4hHUeh3pKhe2BXXhimsVi21JEBLlScKQEEKIYok5l8iHG/4yKXfR5tAv7iD1sq4DcN41hO1Bncg2c8AqgAJYKUFIlLNSCUMajYaMjAy8vb1LozkhhBAVlLlpsZpZN+kfdwAPbTYahZI9fm046hVq9oBVgCH9GtCtaYi1uipEkRU7DGk0GpYtW0aNGjUYOHAgERERTJw4kfT0dNq1a8eCBQvw8vKyRl+FEEKUE3N7Byn1WromHqd9yj8ogNsOXmwO7kqChQNWAVZM7YFSaT4oCVGWzG/3acbChQtZunQp6enpALz//vv4+Pjw2muvcfXqVebOnVvqnRRCCFF+lqz/q8Ag5J2XxuDrP9Ph/4PQMc8GfFW9n8Ug1DDEjS+n95QgJCqUYt8Z2rZtG6+88grPPvssFy9e5Ny5c3z44Yc8/PDDeHt78/HHH/POO+9Yo69CCCHKmLlpsbC0C9yfcBgnvYZspSM7AzsS417TYltytpioqIodhuLj42nevDkAe/fuRalU0rVrVwCCg4MNd4yEEEJUXuamxRx1edwff5gmGZcAuOocxNagLqQ7uFlsT54WExVZscNQYGAg169fp02bNvz66680atQIX19fAI4dO0ZwcHCpd1IIIUTZ+XrnGX7/66ZJeZWcBAbG7sNHk4EOBft9mxPh08TiAau+rjBnggQhUbEVOwwNHDiQDz74gK1bt3LkyBHefPNNAGbPns13333HSy+9VOqdFEIIUTZGfLgH/X/KFHod7ZP/ITzpOCr0pKrd2BIUzg2XQIttLZnUFWdn2cFFVHzF/l06YcIEnJ2diYqK4tVXX+WZZ54B4MSJEwwfPpwxY8aUeieFEEJYX0Hrg9w1WQyI20/N7FgATrnXYldAB3JV5g9YBZkWE5VLscOQQqHgxRdf5MUXXzQqX7duXal1SgghRNnJy9Py0rw/TcrrZV6jb9xBXHW55CnU/BrQjhMedWXvIGFzSnT/Mi8vj/Xr13Pw4EESEhJ4//33iYyMJCwsjGbNmpV2H4UQQljJgvV/cfx8olGZWqehR+IRWqeeBSDWyZfNQV1JtnDAKsjeQaLyKnYYSkpKYsiQIVy8eJE6depw/vx5cnJy+OOPP/jwww/56quvaNmypTX6KoQQohS9vzqa87fSjMr8c1MYGLeXwLwUAA57N2avX0u0Fg5YffqBuvRuYfmxeiEqsmJvuvjxxx+TmZnJjh07+Omnn9Dr7yy1W7BgAU2bNmXBggWl3kkhhBCla/iHe4yDkF5Py9SzDLm+ncC8FDJUznxf5T5+929jMQitmNpDgpCo9Ip9Z+j333/n9ddfp2bNmmi1WkO5k5MTw4cPZ/r06cVq78aNG/TsabrQ7r333uPxxx9nz549LF68mIsXL+Lj40OfPn2YOHEizs7mD/7r2bMnN27cMCobMGAAc+bMKVbfhBDC1pw9n8TkpfuMypy1OfSNj6BB5jUALrhWZXtgZ7LULhbbkkXSwlYUOwzl5uaaPZBVpVKRn59frPbOnj2Lk5MTu3fvRvGvRXkeHh5ER0czbtw4Jk2aRJ8+fbhy5QpvvvkmKSkpfPDBBwW2l5GRwc2bN/n8888JCwszlFsKT0IIYQ+ef2+3SVmNrFj6x+3HU5uFFiV/+LciyquRxUXSQV4qPhjdzZpdFaJMFTsMNW3alLVr19Ktm+kfhK1bt9KkSZNitRcTE0Pt2rUJDDTdr2LdunV06NCBUaNGAVCzZk1efvllXn/9dWbNmoWjo+mjnTExMej1elq1aoWnp+XFfkIIYQ/On0/i/fXHjcqUeh1dkv6iY/IJFECigydbgsOJc/Kz2JbsHSRsUbF/R0+cOJGhQ4fy0EMP0a1bNxQKBdu2bWPhwoXs37+fFStWFKu9s2fPUq9evQKvDR8+HKXSdFmTRqMhIyPDsPP1f9sLCAiwShBSqwtfYqVSKY1+tWUyVttkT2MF2x9vQXeDvPLTGRi7j2q5twH4y7Meu/3bkq90sNjWmpm9rNJHa7D1z/XfZKz3TqG/uwK6GKKiopg7dy5///03Op0OhUJB48aNeeWVV+jcuXOx2urXrx8BAQHk5eVx+fJlatasyZgxYwgPDzepm5eXxxNPPIFSqWTjxo0Ftjdr1iwOHDhAaGgox44dw9fXl0GDBvH8888XGKyKSq/XG03jCSFERXbtZjpj5ppuotgo/RJ9Eg7hrMsnR+nAzwEdOeNRq9D2ts59yAq9FKJiKNG9zrZt27Ju3TpycnJITU3F3d0dNzfLh/QV5G4AcnFxYerUqbi6urJlyxZGjhzJqlWr6Nixo6GuRqNh6tSpnD9/nm+//dZsm+fOnSM9PZ2+ffsybtw4oqOjmTNnDqmpqUycOLEkwwVAp9OTlpZVaD2VSomnpwtpadlotboSf7/KQMZqm+xprGCb4y3obpCjLp/eCZE0Tb8AwHXnALYEhZPm4G6xrcf61GZg27okJ2dapa/WYoufqzkyVvM8PV2KdBfpniZ+nZ2d72lhsqOjI1FRUajVasP6nyZNmnDhwgVWrlxpCEMZGRlMmjSJw4cPs2DBApo3b262zVWrVpGbm4u7+50/4KGhoWRmZrJ06VLGjx9/T3eHNJqi/ybTanXFql+ZyVhtkz2NFWxjvPHxmUz/8rBJeXDObQbG7cM3Px0dCg76NOWAbzOLB6zC/zZRrMzviy18rkUlYy25Yoehhg0bFjpddPr06SK35+rqalLWoEED9u/fD0B8fDwjR47k+vXrLF++nA4dOlhsz8HBAQcH43nvBg0akJWVRWpqKj4+PkXumxBCVBYFnSuGXk+7lH/olngcFTpS1a5sDQrnuktQoe3JY/PCnhQ7DI0dO9YkDGVmZnL06FGuXr3K5MmTi9zWmTNnePrpp1m+fDlt2rQxlJ88eZJ69eqRmprKkCFDyMjIYO3atYSGhlpsT6fT0atXLx5//HFGjx5tKD9x4gT+/v4ShIQQNicrK59xC/aZlLtpsugfd4Da2bcAOONWg52BHclVORXapgQhYW+KHYbGjx9v9tq0adM4efIkjz76aJHaatCgAfXr12fWrFm89dZb+Pj48MMPP3D8+HHWr1/PBx98wLVr11ixYgW+vr4kJCQYXuvr64tKpSI9PZ38/Hx8fX1RKpX06dOHFStWUKtWLcLCwoiIiGDFihXMmDGjuEMVQogK7ZUF+0nJyjMpr5N5nX7xB3HT5pCvULHbvy1/eda3uHfQXRKEhD0q0dNk5kRERBjW9hRVUlISc+bMYe/evaSlpdG4cWMmT55Mq1ataNGiBbm5uQW+7rfffiMkJITp06cTGRnJnj13bhFrNBqWL1/Ohg0biI2NJSQkhOHDh/PEE0/c09i0Wh1JSYUvIFSrlfj4uJGcnGnzc7cyVttkT2OFyjvegqbFVHot3W8fpW3qnaUKcY4+bAkOJ9HRu9D23nu+LVWrepR2N8tNZf1cS0LGap6vr5v1F1D/1+XLl9FoNMV6ja+vL++//36B1/7+++9CX//hhx8afa1Wqxk9erTRNJkQQtgKjUbHqDl/mJT75aUwMHYfQXnJAER7NeR3v9ZolebPFQNoVcuJcU8Vb0sUIWxNscPQokWLTMp0Oh23bt1ix44dBZ4zJoQQ4t4t/+kkEWfjjQv1epqnnaPX7Sgc9FqylE5sD+rMBbeQQtv7YnL3Im0mK4StK5UwBODu7k7v3r157bXX7rlTQgghjBU0LeakzeXB+EM0zLwCwCWXKmwL6kym2vQp3f+StUFC/E+xw9CZM2es0Q8hhBAFSErKZvIXESblIdlxDIjbh5cmCy0K/vRrSaR3WKGLpFe81htHld7m15YIURxy2p4QQlRQL3y0B91/HnFR6HV0TvqbTsknUKInycGDLUHhxDr7F9rempm98PFxrXS7SQthbUUKQ88//3yRG1QoFKxevbrEHRJCCFHwtJhnfgYD4vZRPefONiMnPOrya0A78go5YPXTMZ3x9Cx8fyEh7FWRwlBxnr4vxSf1hRDC7pw6m8Ccn06YlIdmXObB+AicdfnkKhzYFdieUx51Cm1P1gYJUbgihaGvv/7a2v0QQgi7V9DdIAddPr1uR9E87TwAN5z82RIcTqqD5T2B5r3UCW/vkp8dKYQ9KdU1Q1lZWURHR9O1a9fSbFYIIWzawb+usWLnOZPyoNxEBsbuwy8/DT0Q4dOU/b7N0RVywKrcDRKieIodhm7cuMGbb75JVFQU+fn5BdYpzkGtQghhz8wdsNom9TTdbx9FjY50lQtbg8K56hpcaHsShIQovmKHoQ8++IBjx47xxBNPcPToUVxcXGjRogUHDhwgJiaGhQsXWqOfQghhU3Q6PS98/LtJuasmm37xB6mbdQOAGLfq7AjsSI7K8pTXuAGNaBVWxSp9FcLWFXvr0aioKCZNmsTMmTN59NFHcXR0ZMqUKWzYsIG2bdvy22+/WaOfQghhM/YevVFgEKqVdZPh17ZSN+sG+QoVuwLasTG4e6FB6MvpPSUICXEPin1nKDMzk0aNGgFQt25dw50glUrFs88+a3JWmBBCiP8Z8dEe/vvQrVKvpVviMdqnnAIgwdGbzUHh3HbyKbQ9mRYT4t4VOwwFBgaSkHBnj4uaNWuSmppKfHw8gYGBeHl5kZiYWOqdFEKIyi4rK59xC/aZlPvkpfFQ3F6Cc5MAOOIVyu9+rdEoLf/1/OHw9gQGulmlr0LYm2KHoW7duvHZZ58RHBxMq1atCA4O5ssvv2TcuHFs2LCBoKAga/RTCCEqrUmf7SUtW2NcqNfTNP0CvRMicdRryFY6siOwE+fcaxTantwNEqJ0FXvN0IQJE/D09GTBggUAvPzyy6xZs4a2bduydetWhg0bVuqdFEKIymr4h3tMgpCTNo+BcfvoF38QR72GKy7BrKwxUIKQEOWkSHeGnnzySR577DH69euHj48PP/74I/Hx8QAMHDiQqlWrcvz4cZo1a0a7du2s2mEhhKgMzO0kXTU7gYFx+/DWZKBDwV7fFhz2CUNfyN5Bbz7Vilq1vK3UWyHsW5HCUE5ODm+88QYffPABffv25bHHHqNFixaG623atKFNmzbW6qMQQlQqBe0dpNDr6Jh8ki5Jf6FET4ranS3B4dx0DrDYVngzF4b17WitrgohKGIY2rx5M2fOnGHTpk1s27aNDRs2ULduXR577DEGDhyIr6+vtfsphBAV3s2b6cxcE2VS7qHJZEDsfmrkxAHwj3ttfgloT67K0WJ7K6b2QKlUWKWvQoj/UeiLebKqTqdj3759/PTTT/z+++/odDp69uzJE088QefOna3Vz3Kn1epISsostJ5arcTHx43k5Ew0Gl0Z9Kz8yFhtkz2NFUpvvAXuJA00yLjKg/EHcdHlkatQ82tAe0561AGF5ZBjjbVB9vTZylhtU3HH6uvrhkpV+PLoYj9NplQq6datG926dSMjI4Pt27ezefNmXnjhBapUqcKgQYMYN25ccZsVQohKKSdHw5hP95qUq3Uaet6OplVaDAC3nPzYHBROiqOnxfYmPRRGs0byVK4QZanYT5P9m7u7O08++SRr165lzZo1ODo6snjx4tLqmxBCVGgzvzhUYBAKyE1myPXthiB0yDuMr0MeKDQIfTm9pwQhIcrBPZ1aHxcXx/bt29m6dStnzpyhWrVqjB8/vrT6JoQQFZa5A1ZbpZ6lZ2I0ar2ODJUL24I6c9m1aqHtySPzQpSfYoehjIwMdu3axdatW4mKikKtVtOrVy+mTp1Kx47yxIMQwvYVFIRctDn0jTtI/azrAJx3rcb2oM5kF3Ku2PRHm9Ogvp9V+imEKJoihSGNRsOff/7Jli1b+OOPP8jNzaVx48a8/vrrDBw4EA8PD2v3Uwghyl3kiRss237WpLxm1i36x+3HQ5uNBiW/+7fmiFfDclkkLYQoviKFoc6dO5OWloanpyePP/44jz32GA0bNrR234QQosIo6G6QUq8jPPE4HVJOogBuO3ixJTiceKfCtxuRICRExVGkMBQWFsZjjz1Gr169cHS0vC+GEELYkvj4TKZ/edik3Ds/nYGx+6iaexuAY571+c2/rRywKkQlVKQw9OWXX1q7H0IIUeGY2zsoLP0i98cfxkmfT47SkZ2BHTnrXrPQ9uRukBAV0z09TSaEELaqoCDkqMvj/oRImqRfBOCacyBbg7qQ5uBusa3XH2tBvXqyU78QFZWEISGE+Jej/9xi0dbTJuVVcm4zMHYvPv9/wOoB32Yc9Gla6AGrcjdIiIpPwpAQQvw/c3sHdUg5SXjicVToSVW7sSUonBsugRbben9oO4KDLd8xEkJUDBKGhBB271Z8JqM+2m1S7q7Jon/cfmplxwJw2r0mPwd0LPSAVbkbJETlUqQwdPPmzWI1WrVq4but3nXjxg169jT9i+O9997j8ccf5/Tp08yePZuTJ0/i7e3N4MGDGTFihMU2d+7cycKFC7l27Rq1atViypQpdO3atVhjEELYh+ffMw1BAPUyr9E37iCuulzyFGp2B7Tlb496sneQEDaoSGGoZ8+eKAr5C+DfTp82nW835+zZszg5ObF7926j7+Hh4UFycjLDhg2jV69ezJo1i+PHjzNr1iy8vb159NFHC2zv0KFDTJkyhenTp9OxY0fWr1/P2LFj2bRpE3Xr1i1yv4QQtq+gaTGVTkvPxGhap97ZXDHO0YfNwV1JcvQqtD0JQkJUTkUKQ++//74hqKSmpjJnzhw6duzIgw8+SEBAACkpKezZs4c//viD6dOnF6sDMTEx1K5dm8BA0/n31atX4+joyNtvv41araZu3bpcuXKF5cuXmw1Dy5cvp3fv3jz33HMATJs2jWPHjrF69WreeeedYvVNCGGbTp6JZ96mkyblfnkpPBS7l8C8FAAivRrxp38rtAqVxfbmvdQJb2/Lx24IISquIoWhQYMGGf5/7NixPPLII7z77rtGdQYMGMDs2bPZuXMnTz75ZJE7cPbsWerVq1fgtejoaNq2bYta/b9udujQgc8//5zExET8/IzP89HpdBw9etQkkLVv355ff/21yH0yR622/NQIgEqlNPrVlslYbZOtj7XAaTG9nhZpMdx3OxoHvZZMlTPbAjtzya1aoe2tmdnLCr20Dlv/bP9NxmqbrDXWYi+gPnDgAIsXLy7wWvfu3fnhhx+K1V5MTAwBAQE888wzXL58mZo1azJmzBjCw8OJjY2lQYMGRvXv3kG6efOmSRhKS0sjKyuL4OBgk9fcunWrWP36L6VSgY9P0XeN9fR0uafvV5nIWG2TrY313MVkXlm816TcWZvLg/ERhGZeBeCia1W2BXYmS215/O8M6UDLZkFW6au12dpna4mM1TaV9liLHYZ8fHw4fvw4nTt3Nrl26NAhgoKK/pdDXl4ely9fxsXFhalTp+Lq6sqWLVsYOXIkq1atIicnx+T4DycnJwByc3NN2svJyQEo8DUF1S8OnU5PWlpWofVUKiWeni6kpWWj1eru6XtWdDJW22SLYzW3SLp6diwD4vbjqclCi5I//FoR5d2o0EXSd+8GJSdnlnpfrckWP1tzZKy2qbhj9fR0KdJdpGKHoccff5wlS5aQnZ1Nz5498fX15fbt2/z888989913vP7660Vuy9HRkaioKNRqtSHANGnShAsXLrBy5UqcnZ3Jy8szes3dUOPq6mrS3t2gVNBrXFzuPUVqNEX/TabV6opVvzKTsdomWxlrQYukFXodXZL+olPyCRRAooMnW4LCiXP2M23gP76c3rPSvy+28tkWhYzVNpX2WIsdhkaPHk16ejpfffUVK1euBECv1+Ps7MzEiRN59tlni9VeQaGmQYMG7N+/n+DgYOLj442u3f26oDtQ3t7euLq6Fvia/06dCSFsW1JSNpO/iDAp98rPYEDcPkJyEgD426Muvwa0I1/pYLE92URRCNtV7DCkUCiYNm0aY8aM4fjx46SmpuLj40PLli0LDDaWnDlzhqeffprly5fTpk0bQ/nJkyepV68ejRo1Yt26dWi1WlSqO09zREREULt2bZP1Qnf71qpVKyIjI3n88ccN5YcPH6Z169bFHaoQopIa8dEe9HrT8kbpl+iTcAhnXT45Sgd+DujAGY/ahbYnj8wLYdtKvAO1m5sbAQEB6PV6mjdvTl5eXrHDUIMGDahfvz6zZs3irbfewsfHhx9++IHjx4+zfv16/P39WbFiBTNmzOCFF17g77//ZvXq1cyaNcvQRnp6Ovn5+fj63jkEcdiwYYwaNYrGjRvTtWtXNmzYYNi4UQhh2y5fTuGddUdNyh10+fROiKRZ+gUArjsHsDUonNRCDlj9cHh7AgOL/uCEEKJyUuj1Bf37ybLNmzczd+5cEhISUCgU/PjjjyxcuBAHBwfmzp1rsoDZkqSkJObMmcPevXtJS0ujcePGTJ482XCn6O+//2b27NmcOnWKgIAAhg8fbthDCGD69OlERkayZ8//1gVs2rSJJUuWEBsbS7169ZgyZQodO3Ys7jCNaLU6kpIKXyypVivx8XEjOTnT5uduZay2qbKOtcBzxYCgnEQeituLb346OhRE+DRlv28zuzxgtbJ+tiUhY7VNxR2rr69bkRZQFzsM7dixg1deeYWBAwfSo0cPXn75ZTZs2MDZs2eZNWsWw4YNY9KkScVpslKQMGRKxmqbKttYf9xznJ2RSaYX9HrapZyiW+IxVOhIU7uyNSicay6FP/Fqi0EIKt9ney9krLbJWmGo2NNky5Yt46mnnuLtt99Gq9UaygcNGkRiYiI//PCDTYYhIUTFY+5ukJsmm35x+6mTfWd/sbNuNdgZ2JEclZPF9t57vi1Vq3qUej+FEBVbsbdwvHTpEr179y7wWvPmzYmLi7vnTgkhhCVpablmg1CdzBsMv7aVOtm3yFeo+DmgAz8Fdys0CK2Z2UuCkBB2qth3hvz8/Lhw4UKBmy5euHChwKe8hBCitIz6ZA8arWm5Sq+l++2jtE29c1B0vKM3m4O7kujobbG9qYOaEd65dqXbQFEIUXqKHYb69u3LggULCAwMpFu3bsCdR9pPnjzJkiVL6N+/f6l3UgghdDo9L3z8e4HXfPNSeSh2L0F5yQBEezXkd7/WaJWWD1j9cnrPIp05KISwbcUOQ5MmTSImJoZJkyahVN75S2Tw4MFkZWXRpk0bJk6cWOqdFELYt29+PsWe47GmF/R6mqWfp1dCFI56DVlKJ7YHdeaCW0ihbdrqImkhRPEVOww5OjqyYsUKDhw4wKFDh0hJScHDw4N27drRrVs3FIWc6SOEEMVhbm2QkzaPBxIiaJRxBYDLLsFsC+pChtryfmfvPNeGkBDPUu+nEKLyKnYYeumll3j++efp3LlzgeuGhBCitJgLQtWy4xkYtw8vTSZaFOz1a8lh77BCD1iVu0FCiIIUOwxFRUUxbNgwa/RFCCEAiDxxg2Xbz5qUK/Q6OiWfoHPS3yjRk+zgwZagcG45+1tsb/gD9ejSooa1uiuEqOSKHYY6d+7Mjz/+SIsWLQynxAshRGkxdzfIMz+DAXH7qZ5z5yDmkx51+CWgPXmFHLAqd4OEEIUpdhhycnJi586d/Prrr4SEhJg8Sq9QKFi9enWpdVAIYR+ysvIZt2BfgddCM67wYHwEzro8chUO7ApszymPOhbbWzCuC+7uRT8aSAhhv4odhmJjY2nZsqXh6/+e5lGCo86EEHZu4md7Sc/WmJQ76PK573Y0LdLOAXDTyZ8tweGkOFjeHFHuBgkhiqPYYejrr7+2Rj+EEHYoL0/LS/P+LPBaYG4SA2P34Z+fih6I8GnCft8W6CwcsPrmU62oVcvbOp0VQtisYoehu1JTU4mOjiY+Pp4+ffqQkpJC7dq15dF6IUSRvLsqkktxGaYX9HrapJ6h++0jqNGRrnJhW1AXrrhWsdie3A0SQpRUicLQ0qVL+fzzz8nJyUGhUNCsWTPmz59PSkoKX375JZ6esoeHEMI8c4ukXTXZ9I0/SL2sGwCccw1hR1AnslXOZtsaN6ARrcIsByUhhLCk2PvQf/PNNyxcuJBhw4bxww8/GNYIDRkyhGvXrvHZZ5+VeieFELbDXBCqlXWT4de2US/rBhqFkl/827GhSg+LQejL6T0lCAkh7lmxw9DXX3/NqFGjmDhxImFhYYby8PBwJk2axJ49Bf9FJ4Swb0lJ2QUGIaVeS/fbR3jq5m7ctdkkOHqzOqQfR70bWtxEUabFhBClpdjTZDdv3qRdu3YFXqtTpw63b9++504JIWzLiA/3UNBzpj55aQyM20eV3EQAjno2YI9/GzRK83811Q5Q88aIrlbqqRDCHhU7DFWpUoVjx47RqVMnk2snT56kShW5ZS2E+J8Cp8X0epqkX+T+hMM46jVkKx3ZEdiJc+6Wd4le9ko3HB0tn0QvhBDFVeww9Nhjj7Fw4UKcnZ3p3r07AFlZWezatYvPP/9cjuoQQgBwJuY2H2/826TcSZvH/QmHCMu4DMBV5yC2BnchXe1msT2ZFhNCWEuxw9DIkSO5fv06c+bMYc6cOQA8//zzAAwYMIAXX3yxdHsohKh0zC2SrpqTwMDYfXhrMtChYJ9vcw75NEFvYe+g955vS9WqljdZFEKIe1HsMKRQKHjnnXcYNmwYhw4dIjU1FQ8PD9q1a0f9+vWt0UchRCVSUBBS6HV0SP6H8KTjKNGTonZnS1A4N10CLLYld4OEEGWhxJsu1q5dm9q1a5dmX4QQldjlyym8s+6oSbm7JosBcfuomR0HwD/utfgloAO5KvPnhskp80KIslSkMPTaa68Vq9EPPvigRJ0RQlRO5qbF6mdcpW/8QVx0eeQp1PwS0J6THnXkkXkhRIVSpDB0+PBho6/j4+PRaDRUrVqVgIAAUlJSuHbtGo6OjjRs2NAqHRVCVEwFBSG1TkPP29G0SosB4JaTH1uCwkl2tLw7vQQhIUR5KFIY+vdGilu3bmXOnDksXLiQZs2aGcrPnz/P2LFjefDBB0u/l0KICic+PpPpXx42KffPTeahuH0E5KUAcNi7MX/6tUSnsPxIvAQhIUR5Kfaaofnz5/Pqq68aBSGAevXqMXHiRD744AOGDBlSah0UQlQ85vYOapV6lp6J0aj1OjJUzmwL6sJl16oW2/r4hQ74+7taqadCCFG4Yoeh5ORkPDwKfsxVrVaTlZV1z50SQlRMKSk5vLLsoEm5izaHvvEHqZ95HYDzrtXYEdiJLLWLxfbkbpAQoiIodhhq0aIFixYtokWLFvj4+BjK4+PjWbhwIe3bty/VDgohKoYXPt6DTmdaXiMrlgFx+/DQZqNByR/+rYn2snyuGEgQEkJUHMUOQ9OmTeO5556jZ8+etGzZEh8fHxITEzl27BheXl4sXbrUGv0UQpSTrKx8xi3YZ1Ku1OvoknScjsknUQCJDp5sDu5KvJOvxfY+HdMZT08nK/VWCCGKr9hhqGHDhmzfvp2vvvqKo0ePcuPGDXx8fBg+fDhDhgzB29vbCt0UQpSHCZ/tJSNbY1LulZ/OwNh9VMu9czDzcc96/Obflnylg8X25G6QEKIiKnYYWrZsGffddx/Tpk2zRn+EEBXE8+/tLrC8cfpF+sQfxkmfT47SkZ2BHTnrXtNiW4GeSj4c090KvRRCiHtX7DC0YsUKGjdubJWjNy5dusSgQYN44403GDRoEIMHDyYyMrLAuh999BEPP/xwgdd69uzJjRs3jMoGDBhgOEtNCGHe+fNJvLPONAg56vLpnXCYpukXAbjmHMjWoC6kObhbbG/JpK44O5d4s3shhLC6Yv8NVatWLc6dO0fXrl1LtSP5+flMnjzZ6Gm0hQsXkp+fb1Rv5syZXL16lV69ehXYTkZGBjdv3uTzzz8nLCzMUO7s7Fyq/RXCFpnbSTo45zYD4/bhm5+ODgUHfJtx0KepxQNWQabFhBCVQ7HDUPfu3Zk/fz6///479evXx8/Pz+i6QqFg7Nixxe7IwoULcXNzMyr77/qjbdu2sX//fjZu3Ii7e8H/Go2JiUGv19OqVSs8PS3vdiuE+B9zewe1T/mHronHUKEnVe3G1qAuXHcJKrQ9CUJCiMqi2GFo0aJFAERHRxMdHW1yvSRhKCoqiu+//55NmzbRvXv3AutkZWXx8ccfM2TIEEJDQ822dfbsWQICAqwShNRqy/8KBlCplEa/2jIZq224fTurwL2D3DRZ9I87QO3sWwCccavJzsAO5KosPwn28QsdCA62PHVWkdjyZ/tfMlbbJGO9d8UOQ2fOnCnVDqSlpTF16lRmzpxJlSpVzNZbt24dmZmZjB492mJ7MTExuLq6Mn78eI4dO4avry+DBg3i+eefR6ks+ZunVCrw8XErvOL/8/S0vNmcLZGxVl4DJ29Grzctr5t5nX5xB3DV5ZKnULPbvy1/e9YrdO+grXMfslJPrc/WPltLZKy2ScZaciVe1Zieno5er7/nOzBvv/02LVq0YMCAAWbraLVavv76a5555hmzu1/fde7cOdLT0+nbty/jxo0jOjqaOXPmkJqaysSJE0vcT51OT1pa4btrq1RKPD1dSEvLRqstYIc6GyJjrdwKelpMpdPSI/EIbVLv/KMnztGHzcFdSXL0stjW9Eeb07hRAMnJmVbpqzXZ4mdrjozVNslYzfP0dCnSXaRihaELFy6wfPlyfvvtNzIyMgBwdXWlV69eDB8+3OL0VUE2bdpEdHQ0W7dutVgvMjKSmzdv8sQTTxTa5qpVq8jNzTWsKQoNDSUzM5OlS5cyfvz4e7o7pNEU/TeZVqsrVv3KTMZa+RS0PsgvL4WHYvcRmJcMQJRXI/7wb4W2iAesVvb3xVY+26KQsdomGWvJFTkM7dixg9deew2lUkmnTp2oUaMGarWaa9eusWfPHnbu3Mn7779P//79i/zNN2zYQGJiosk6obfeeouVK1eyfft2AHbv3k2zZs2oXr16oW06ODjg4GC88VuDBg3IysoiNTXV6AgRIezNmZjbfLzxb+NCvZ7maefodTsKB72WTJUz2wM7cdEtxGJb817qhLe3PKUphKj8ihSGLly4wGuvvUa3bt145513TJ7yysjI4K233mLmzJk0atSIunXrFumbz5kzh5ycHKOy+++/nwkTJtC3b19D2ZEjRwgPDy+0PZ1OR69evXj88ceN1hadOHECf39/CULCrhV0N8hZm8uD8RGEZl4F4JJLFbYFdSFTDlgVQtiRIoWhr776inr16jF//nxUKtNb5u7u7nzyySc888wzrF69mnfeeadI3zwoqODHc/38/KhWrRpwZ73Q+fPnGTFiRIF109PTyc/Px9fXF6VSSZ8+fVixYgW1atUiLCyMiIgIVqxYwYwZM4rUJyFsUUFBqHp2HAPi9uGpyUKLkj/9WhLp3djiIulB99Wgf9t61uyqEEKUuSKFoYiICEaPHl1gELpLqVTy1FNPGR69Ly0pKSnk5+ebPfNs9uzZREZGsmfPnb/sX331VTw9PZk7dy6xsbGEhIQwY8aMIq03EsLWXL2ayttrjxiVKfQ6uiT9TcfkEyjRk+TgweagrsQ5+5lp5Y4VU3ugVFp+mkwIISqjIoWh+Ph4ata0fPYQQEhICAkJCffUobNnzxp97efnZ1L2bx9++KHR12q1mtGjRxf6CL4Qtq6gu0Fe+RkMiNtHSM6dP6d/e9Tl14B2csCqEMKuFSkMeXp6Eh8fX2i9hIQEfH1977lTQoiSS0nJKXATxYbpl3kgIQJnXT45Sgd2BXTgtEdti229P7RdpdpAUQghSqJIYahVq1Zs3LjRaFFzQX766SdatWpVKh0TQhTfqI9/R6Mz3kXRQZdPr4QomqefB+CGkz9bgsNJdbC8Z9eamb3s5jFdIYR9K9KmO0OGDOHAgQMW1wPNnz+fAwcOMGTIkFLrnBCiaHQ6PcM/3GMShIJyEhl6bTvN08+jBw74NOXbkAcKDUKVeSdpIYQoriLdGWrdujUvv/wy8+bNY8eOHfTo0YOQkBDUajU3btzg119/5eLFi0ybNo1mzZpZu89CiH/Zffgqa38/b1yo19M25TTdE4+iQkeaypVtQV246hpssa3HetdiYHt5WkwIYV+KvOniqFGjqF+/PosWLWLlypVG11q0aMHy5cvp0qVLqXdQCGFeQYukXTXZ9Is/QN2smwDEuFVnR2BHclSWN0iUp8WEEPaqWMdx9OjRgx49epCcnMyNGzfQ6/VUq1ZNFk0LUcaO/nOLRVtPm5TXzrxB//gDuGlzyFeo+M2/Dcc9GxR6wKo8LSaEsGclOqjVx8dHdnMWohxcvpzCO+uOmpSr9Fq6JR6jXcopAOIdvdkS1JXbTt6FtilBSAhh70p8ar0Qouzcvp3F1BWHCrzmm5fKwLh9BOcmAXDEK5Tf/VqjUVr+4921iR9D+zcv9b4KIURlI2FIiAruhY/28J+HxO7Q62mWfp5eCVE46jVkKZ3YEdSJ826FH2j8xeTuqNVFephUCCFsnoQhISqopKRsJn8RUeA1J20eDyRE0CjjCgCXXYLZFtSFDLVroe3KtJgQQhiTMCREBTTy4z1ozex3WC07noFx+/DSZKJFwT6/Fhz2DkOvsHynZ8G4Lri7O1qht0IIUblJGBKiginocXm4c8Bqx+STdEn6CyV6ktXubAnuyi1n/0LblLtBQghhnoQhISqIvDwtL837s8BrHvmZDIjbT42cOABOetThl4B25Ckt3+mZ91InvL0t7y8khBD2TsKQEBXAnG+PcupaSoHXGmRc4cH4CFx0eeQq1PwS0J5/POsW2qbcDRJCiKKRMCREOTM3LabWabjvdhQt084BcNPJjy1B4aQ4elpsz9sZ5k2SICSEEEUlYUiIcnLxYjLv/XCswGuBuUkMjN2Hf34qeuCQdxj7/FqgU6jMtufvreadoZ1wdpY/1kIIURzyt6YQZczSBoro9bROPUOPxCOo9TrSVS5sC+rCFdcqFtuUc8WEEKLkJAwJUUYsLZAGcNHm0C/uIPWyrgNw3jWE7UGdyLZwwOqcUR3x9XUp9b4KIYQ9kTAkRBmYt+44Jy8nmb1eM+sm/eMO4KHNRqNQssevDUe9Qi0esCoLpIUQonRIGBLCyswepwEo9Vq6Jh6nfco/KIDbDl5sDu5KgpP5g5DlbpAQQpQuCUNCWElKSg6vLDto9rp3XhoPxe2jSm4iAMc8G/CbfxuLB6zK3SAhhCh9EoaEKGU6nZ4X5/xu9jgNgLC0C9yfcBgnvYZspSM7AzsS417TYrsShIQQwjokDAlRio6cjWfxTyfNXnfU5XF//GGaZFwC4KpzEFuDupDu4Gb2NbKLtBBCWJeEISFKyeF/Yvl86ymz16vkJDAwdh8+mgx0KNjv25wInyYWD1iVu0FCCGF9EoaEuEcZGXm8vGg/WjPXFXod7ZP/ITzpOCr0pKrd2BIUzg2XQLNtPv1AXXq3sDxtJoQQonRIGBLiHkz4bC8Z2Rqz1901WQyI20/N7FgATrnXYldAB3JVBR+w2rGZLyMeaC4bKAohRBmSMCRECeh0el74+HeLdeplXqNv3EFcdbnkKdT8GtCOEx51ze4dJLtICyFE+ZAwJEQxHfz7Fit2nDZ7Xa3T0CPxCK1TzwIQ6+TL5qCuJJs5YPXTMZ3x9HSySl+FEEIUTsKQEEWUp9Hxwod7yNOYf2bePzeFgXF7CcxLAeCwd2P2+rVEa+aAVVkgLYQQ5U/CkBBFsHr7aX47dsN8Bb2elmkx9LwdjYNeS4bKme2BnbnkVq3A6kFeKj4Y3c1KvRVCCFEcFSoMXbp0iUGDBvHGG28waNAgAF577TU2btxoVC8oKIi9e/eabWfnzp0sXLiQa9euUatWLaZMmULXrl2t2ndhuywdpwHgrM2hb3wEDTKvAXDBtSrbAzuTpS74yIwlk7ri7Fyh/ugJIYRdqzB/I+fn5zN58mSysrKMys+ePctLL73Ec889ZyhTqQqecgA4dOgQU6ZMYfr06XTs2JH169czduxYNm3aRN26da3Wf2Gbhn+4x+L1Glmx9I/bj6c2Cy1K/vBvRZRXowIXSS8Y1wV394KfIhNCCFF+zO/2VsYWLlyIm5vxLrxarZbz58/TtGlTAgICDP/5+vqabWf58uX07t2b5557jrp16zJt2jTCwsJYvXq1tYcgbIylIKTU6+iaeIynb/6CpzaLRAdP1lR/kCjvxgUGoS+n95QgJIQQFVSFCENRUVF8//33fPTRR0blly9fJjc3t8h3dHQ6HUePHqVDhw5G5e3btyc6OrrU+itsn6Ug5JWfzrPXf6ZT8gkUwF+e9fiqej/inPwKrC+LpIUQomIr92mytLQ0pk6dysyZM6lSpYrRtZiYGBQKBatXr2bv3r0olUq6devGpEmT8PDwKLCtrKwsgoODjcoDAwO5devWPfdVrS48O6pUSqNfbZktjlWj0VkMQo3SL9En4RDOunxylA78HNCRMx61Cqw776VO+Pu7Wqmn1mOLn6sl9jReGattkrHeu3IPQ2+//TYtWrRgwIABJtfOnTuHUqmkWrVqLFu2jCtXrvDRRx8RExPD6tWrUSqN34ycnBwAHB2NpyOcnJzIzc29p34qlQp8fMwfpvlfnp4FL561RbYw1vSsfEZ/8AupWQXvJu2oy6d3QiRN0y8AcN05gC1B4aQ5uJvU7d+pOi8+2sqq/S0LtvC5Foc9jVfGaptkrCVXrmFo06ZNREdHs3Xr1gKvjx8/nqFDh+LpeWezugYNGhAQEMCTTz7JiRMnaN68uVF9J6c7G9fl5eUZlefm5uLicm9vnE6nJy0tq9B6KpUST08X0tKy0WrN70djC2xlrK8u2k9CSo7Z68E5txkYtw/f/HR0KDjo05QDvs0KPGD1y+k9UauVJCdnWrPLVmUrn2tR2dN4Zay2ScZqnqenS5HuIpVrGNqwYQOJiYl0797dqPytt95i5cqVbN++3RCE7mrQoAEAsbGxJmHI29sbV1dX4uPjjcrj4+NNps5KQmNhs73/0mp1xapfmVXmsb405w/zmyjq9bRL+Yduicf+/4BVV7YGhXPdJcikqhr44v/XBlXW9+K/KvPnWhL2NF4Zq22SsZZcuYahOXPmGKa27rr//vuZMGECffv25dVXXyUlJYWVK1carp84cQKAevXqmbSnUCho1aoVkZGRPP7444byw4cP07p1ayuNQlRWM784YDYIuWmy6B93gNrZd9aanXGrwc7AjuSqTI/NkEfmhRCicivXMBQUZPovbAA/Pz+qVatG//79GT16NEuXLqVfv35cunSJd955h/79+xueMEtPTyc/P9/wuP2wYcMYNWoUjRs3pmvXrmzYsIHTp08ze/bsMhuXqPjmrjvKzaSC15HVzbxOv7gDuOpyyVeo2O3flr8865s8Mt+oujtTnm1XFt0VQghhReW+gNqSHj168Nlnn7Fs2TKWLVuGh4cHAwYMYNKkSYY6s2fPJjIykj177jwB1KVLF95//32WLFnC/PnzqVevHsuWLZMNF4XB1CUHuJ1mGoRUOi09Eo/QJvUMAHGOPmwJDifR0duk7n2tqvLs/Q2t3VUhhBBlQKHX6y0cNCDu0mp1JCUVvihWrVbi4+NGcnKmzc/dVsaxvjTnd/I0pr/l/fJSGBi7j6C8ZACivRryu19rtErT3c5rBrnx1rD2Vu9reamMn+u9sKfxylhtk4zVPF9ft4q/gFqIsjTiwz2YxCC9nuZp5+h1OwoHvZYspRPbgzpzwS2kwDZqBXvw5tC2Vu+rEEKIsiNhSNiFgoKQkzaXB+MjaJh5FYBLLlXYFtSZTHXBGyW++nQrmtf1tfl/eQkhhL2RMCRsmk6n54WPfzcpD8mOY0DcPrw0WWhR8KdfSyK9wwo8V0ytgBWv3Yefn3ul3j9ICCFEwSQMCZsVdSaepZtOGpUp9Do6J/1Np+QTKNGT5ODBlqBwYp39C2zj7mPzSqVpSBJCCGEbJAwJm/TtL2f57egNozLP/AwGxO2jek4CACc86vJrQDvylA4FtiEHrAohhH2QMCRsik6nZ9z8P8nJN17X0zD9Mg8kROCsyydX4cCuwPac8qhjth0JQkIIYT8kDAmbEXkqjmVb/jEqc9Dl0+t2FM3TzgNww8mfLcHhpDp4FNiGClguQUgIIeyKhCFR6el0emZ/Hc2lW+lG5UG5iQyM3Ydffhp6IMKnKft9m6Mr4IBVAB93NXPHdS2DHgshhKhIJAyJSkuj0fHljlMcOmV8MC96PW1TT9Pt9lHU6EhXubA1KJyrruYP621Sy5tXnmpl5R4LIYSoiCQMiUrpu90x/Bp93aTcVZNNv/gD1M26CUCMW3V2BHYkR+Vstq1mdXyY9ERLq/VVCCFExSZhSFQqGo2OSYv2k5WjMblWO/MG/eIP4K7NIV+hYo9/a455hha4d9BdfdqG8OR9DazZZSGEEBWchCFRaZi7G6TSa+mWeIx2KacASHD0ZnNQOLedfCy298Xk7qjVhZ9ZI4QQwrZJGBKVwrSlB0lIzTEp98lL46G4vQTnJgFwxCuU3/1ao1Ga/62tVsAX0+SJMSGEEHdIGBIV3oRP/yQjR2tcqNfTNP0CvRMicdRryFY6siOwE+fca1hsq0ktL155qrUVeyuEEKKykTAkKrQ3V0SYBCEnbR59Eg7ROOMyAFdcgtka1IUMMwes3rXslW44Oqqs1VUhhBCVlIQhUeHodHpOXUrii20nyMg23km6WnY8A+P24aXJRIeCvb4tOOwTht7M3kEg+wcJIYSwTMKQqFCOnI1nxbZT5P7nOA2FXkfH5JN0SfoLJXpS1O5sCQ7npnOAxfZ6t67K070bWrPLQgghKjkJQ6LCOHI2nsU/nTQp98jPZEDcfmrkxAHwj3ttfgloT67K0WxbQT7OvDuigzwtJoQQolAShkSFoNPp+WLrPyblDTKu8mD8QVx0eeQq1Pwa0J6THnXM7h0U4u/KzOfbytogIYQQRSZhSJQ7nU7PmysOk6/RG8rUOg333Y6mZVoMALec/NgcFE6Ko2eBbXi5qflkdBe5EySEEKLYJAyJcnXkbDxLN51E978cREBuMgPj9hKQlwrAIe8w9vq1QKco+G7PyIGN6Ni4Sll0VwghhA2SMCTKjckaIb2e1qln6JF4BLVeR4bKhW1BnbnsWtVsG/K4vBBCiHslYUiUi/+uEXLR5tAv7iD1su4ct3HetRrbgzqTbeGA1QfaVZcgJIQQ4p5JGBJlTqPR8doXEYY1QjWzbtE/bj8e2mw0KPndvzVHvBpaPGD1gXbVeaJn/bLqshBCCBsmYUiUqe/3nGNX5DUAlHod4YnH6JDyDwrgtoMXW4LDiXfyNfv6Kr7OzBouj8wLIYQoPRKGhNXpdHrOXElmza4zxKfcOWzVOz+dgbF7qZqbCMAxz/r85t/W4gGrslBaCCGENUgYElYVdSaeVTtOk5P3v/PFwtIucH/CYZz0GnKUjuwM7MhZ95oW2xn9cBPaNgy0dneFEELYIQlDwmr+PSUG4KjL4/74wzTJuATANedAtgZ1Ic3B3WI7ox8OkyAkhBDCaiQMCatY91sMv0RdN3xdJSeBgbH78NFkoEPBAd9mHPRpavGAVZA7QkIIIaxPwpAoNXfXBq3/8zyXYzPuFOr1dEg5SXjicVToSVW7sSUonBsuhQccuSMkhBCiLEgYEqXiyNl4vtp5hswcjaHMXZNF/7j91MqOBeC0e01+Duho8YDVu+SOkBBCiLJSocLQpUuXGDRoEG+88QaDBg0CYM+ePSxevJiLFy/i4+NDnz59mDhxIs7O5jfj69mzJzdu3DAqGzBgAHPmzLFq/+3F3TtAMddTcHJ24GZ8OgdOxBrVqZd5jb5xB3HV5ZKnULM7oC1/e9SzuHfQXXJHSAghRFmqMGEoPz+fyZMnk5WVZSiLjo5m3LhxTJo0iT59+nDlyhXefPNNUlJS+OCDDwpsJyMjg5s3b/L5558TFhZmKLcUnkTR6HR6th28zM7DV8jN1xVYR63T0CPxCK1TzwIQ5+jD5uCuJDl6Fdq+m7OaoQ82pHWoBCEhhBBlp8KEoYULF+Lm5mZUtm7dOjp06MCoUaMAqFmzJi+//DKvv/46s2bNwtHRdLolJiYGvV5Pq1at8PQs+IRzUXwFPSL/X/65KQyM20tgXgoAkV6N+NO/FVozB6zeFejtxPMPNKJhDR+UysLvHAkhhBClqUKEoaioKL7//ns2bdpE9+7dDeXDhw9HqTR92kij0ZCRkYGvr+lOxWfPniUgIECCUCn6Yc85fv7XI/Im9HpapMVw3+1oHPRaMlXObAvszCW3aoW23adtCE/e16AUeyuEEEIUT7mHobS0NKZOncrMmTOpUsV4d+HGjRsbfZ2Xl8eqVasICwsrMAjBnTtDrq6ujB8/nmPHjuHr68ugQYN4/vnnCwxWxVGUIyBUKqXRr5WVTqfn7NVkos7Eszv6utl6ztoc+sZH0CDzTli66FqVbYGdyVK7WGy/c7NgRvRtXGmO1bCVz7Uo7GmsYF/jlbHaJhnrvSv3MPT222/TokULBgwYYLGeRqNh6tSpnD9/nm+//dZsvXPnzpGenk7fvn0ZN24c0dHRzJkzh9TUVCZOnFjifiqVCnx83Aqv+P88PS2HgYrs4N83+WLTCRJTcyzWq5EVS/+4/Xhqs9Ci5A+/VkR5N7K4SNrZUcnLT7emU7Oqpd3tMlGZP9fisqexgn2NV8Zqm2SsJVeuYWjTpk1ER0ezdetWi/UyMjKYNGkShw8fZsGCBTRv3txs3VWrVpGbm4u7+51djUNDQ8nMzGTp0qWMHz++xHeHdDo9aWlZhdZTqZR4erqQlpaNVlvwIuOKLPJUHIs2nrBYR6nX0TnpLzoln0ABJDp4siUonDhnP4uva9sokLGPNEWpVJCcnFmKvba+yv65Foc9jRXsa7wyVtskYzXP09OlSHeRyjUMbdiwgcTERKN1QgBvvfUWK1euZPv27cTHxzNy5EiuX7/O8uXL6dChg8U2HRwccHBwMCpr0KABWVlZpKam4uPjU+L+ajRF/02m1eqKVb+86XR6th64xJYDly3W88pPZ2DsPqrl3gbgb4+6/BrQjnylg9nXODuqGNa3IW0bBqHT6dHp9KXZ9TJV2T7Xe2FPYwX7Gq+M1TbJWEuuXMPQnDlzyMkxnoq5//77mTBhAn379iU1NZUhQ4aQkZHB2rVrCQ0NtdieTqejV69ePP7444wePdpQfuLECfz9/e8pCNkinU5PzLUUjp1LYN/ftyw+KQbQKP0SfRIO4azLJ0fpwM8BHTjjUdviawZ2rsXAzrXlKTEhhBAVVrmGoaCgoALL/fz8qFatGtOnT+fatWusWLECX19fEhISDHV8fX1RqVSkp6eTn5+Pr68vSqWSPn36sGLFCmrVqkVYWBgRERGsWLGCGTNmlNWwKoUjZ+NZu/scyem5hdZ10OXTOyGSZukXALjuHMDWoHBSLRyw6u7iwJAHQmXPICGEEBVeuS+gNken07Fjxw7y8/MZMmSIyfXffvuNkJAQZs+eTWRkJHv27AHg1VdfxdPTk7lz5xIbG0tISAgzZszgiSeeKOshVFhHzsaz+KeTRaoblJPIQ3F78c1PR4eCCJ+m7PdtZjhgtUF1L1o2CCA9I4+k9Fz8vJxpVNNH9gwSQghRaSj0en3lXcBRhrRaHUlJhS/6VauV+Pi4kZycWSHnbjUaHa8uPkB6dr7lino97VJO0S3xGCp0pKld2RoUzjWX/93NG/NIEx7sUrfCjrU0VfTPtTTZ01jBvsYrY7VNMlbzfH3dKv4CalG2jpyNZ/XPZ8koJAi5abLpF7efOtm3ADjrVoOdgR3JUTkZ6ox+OIwOYcFW7a8QQghRFiQM2YmiTo3VybxOv/iDuGlzyFeo+M2/Lcc96xv2DvL1cOLpXvVlLZAQQgibIWHIDuh0etbuPmexjkqvpfvto7RNPQ1AvKM3m4O7kujoDUCvNiG0qh9Ag+reshZICCGETZEwZAdirqVYfGrMNy+Vh2L3EpSXDEC0V0N+92uNVqmSp8KEEELYPAlDdiAl00wQ0utplnaeXrejcNRryFI6sT2oMxfcQnBzVtO7TQj9O8keQUIIIWybhCE74O3mZFLmpM3lwfhDNMy8AsBll2B2V+9Gy1Z1eVSmw4QQQtgRCUN2oEF1b3w8nAxTZSHZcQyI24+XJhMtCvb6teRUlRbMHdel0pwiL4QQQpQW+clnB5RKBc/0qo/i/w9YfebGL3hpMkl28OCbkAc57NOE5x9sKEFICCGEXZI7Q3aimb+KV3P2o066DMBJjzr8EtAedy93xsqj8kIIIeyYhCEbp9Ppidn1B2z7HnVuDgpnZ7QPPEbVOk2Y5OYka4OEEELYPQlDNuzIietc/+ZrGiWeBeCmkz9/1rmP/vWb0UHuBAkhhBCArBmyWcf2Hid3ySc0SjyLHjjo04RvQh7gSp4Ti386yZGz8eXdRSGEEKJCkDtDNkav15P86y84//g9bnod6SoXtgV14YprFaN63+0+R8v6ATJFJoQQwu5JGLIheSmpXFq2DMX506iAc64h7AjqRLbK2aRuUnouMddSaFjTp+w7KoQQQlQgEoZsxPGf98Omb3HVZKNRKNnj14ajXqGGA1YLYnZnaiGEEMKOSBiq5PQaDf+sXINr1F4AEhy92RIUToJT4Xd8CtqZWgghhLA3EoYqsby4WG59sQzHK5cBOOrZgD3+bdAoC/9YfT3uPFYvhBBC2DsJQ5WQXq8n7eB+4td+gz43l2ylIzsCO3HOvUaR23i6V31ZPC2EEEIgYahS0en0xJy7Rf6W73E6+xcA+SG1+VLdinS1W5Ha8PVw4mnZcVoIIYQwkDBUSRw5G8/uzfvpcXEP3poMdCiIrtIaz54PkH7waqGv79+xJo1r+cqO00IIIcR/SBgqJzqdnphrKaRk5uJdyLEYR07HcnTVOh5K+gslelLU7mwJCuemSwAcvIqbs5rMHI3Z7+Xr4cTD4XUkBAkhhBAFkDBUDo6cjWft7nMkp//v0XYfDyeeKWD6Ki8xkbTPP6Nbxi0A/nGvxS8BHchVORrqFBZxZH2QEEIIYZ4cx1HGjpyNZ/FPJ42CEEByeq7JMRkZx45w6e03qJZxizyFmm2BndkaFG4UhAAycjQ83KU2Ph7Gj8r7ejgx9pEmsj5ICCGEsEDuDJUhnU7P2t3nLNb5bvc5mtf0IvHHdaT++TsK4JaTH1uCwkl29DT7ukBfFz4Z3anIU29CCCGEuEPCUBmKuZZickfov1S3Y7kw6y0UCbEA6Dv24Ov4qugUKouv83ZzQqlUyPEaQgghRDFJGCpDFo+/0OtplXqWnonRKPQ6VJ6eBI8YhUujMLyWHrQYomQDRSGEEKLkZM1QGTJ3/IWLNodHY3/n/tuRqPU69PUaUfPt93ALa4JSqeCZXvUttisLpIUQQoiSkzBUhhpU9zZZ5Fwz6xbDr26lfuZ1NCg5WK0j9aZMQe35v/VBrUMDGftIE1kgLYQQQliBTJOVobt3eRb/dBKAdskn6ZF4FAWQ6ODJ5uCuPP5UV1Qq04zaOjSQlvUDZIG0EEIIUcokDJWxu3d51u4+R9tLp1EAxz3rcbRWFx7v09jiXR5ZIC2EEEKUPglD5cBwlyfam7S0DJrVb8BjcpdHCCGEKBcVas3QpUuXaNmyJRs3bjSUnT59mueee44WLVrQvXt3Vq5cWWg7O3fupG/fvjRt2pQBAwawd+9ea3a7RJRKBQ3bhdGmV3sa1vSRICSEEEKUkwoThvLz85k8eTJZWVmGsuTkZIYNG0atWrXYsGED48eP57PPPmPDhg1m2zl06BBTpkzhmWeeYdOmTXTp0oWxY8dy4cKFshiGEEIIISqZChOGFi5ciJubm1HZDz/8gKOjI2+//TZ169bl0UcfZejQoSxfvtxsO8uXL6d3794899xz1K1bl2nTphEWFsbq1autPQQhhBBCVEIVIgxFRUXx/fff89FHHxmVR0dH07ZtW9Tq/y1t6tChA5cuXSIxMdGkHZ1Ox9GjR+nQoYNRefv27YmOjrZO54UQQghRqZX7Auq0tDSmTp3KzJkzqVKlitG12NhYGjRoYFQWGHjnaaubN2/i5+dn0lZWVhbBwcEmr7l169Y991WtLjw73n0svqDH422NjNU22dNYwb7GK2O1TTLWe1fuYejtt9+mRYsWDBgwwORaTk4Ojo7GJ7Q7Od3ZeDA31/R4ipycHIACX1NQ/eJQKhX4+LgVXvH/eXq63NP3q0xkrLbJnsYK9jVeGattkrGWXLmGoU2bNhEdHc3WrVsLvO7s7ExeXp5R2d1Q4+rqalL/blAq6DUuLvf2xul0etLSsgqtp1Ip8fR0IS0tG61Wd0/fs6KTsdomexor2Nd4Zay2ScZqnqenS5HuIpVrGNqwYQOJiYl0797dqPytt95i5cqVVK1alfj4eKNrd78OCgoyac/b2xtXV9cCX/PfqbOS0GiK/ptMq9UVq35lJmO1TfY0VrCv8cpYbZOMteTKNQzNmTPHMLV11/3338+ECRPo27cv27dvZ926dWi1WlQqFQARERHUrl3bZL0QgEKhoFWrVkRGRvL4448byg8fPkzr1q2tOxghhBBCVErlutoqKCiImjVrGv0H4OfnR7Vq1Xj00UfJyMhgxowZnD9/no0bN7J69WpefPFFQxvp6ekkJSUZvh42bBjbt29n1apVXLhwgY8//pjTp08zZMiQMh+fEEIIISq+Cr303M/PjxUrVnDp0iUeeeQRFi1axNSpU3nkkUcMdWbPns1jjz1m+LpLly68//77fPfddzzyyCMcOnSIZcuWUbdu3fIYghBCCCEqOIVer9eXdycqA61WR1JSZqH11GolPj5uJCdn2vzcrYzVNtnTWMG+xitjtU0yVvN8fd2KtIC6Qt8ZEkIIIYSwNglDQgghhLBrMk1WRHq9Hp2uaG+VSqW0+b0e7pKx2iZ7GivY13hlrLZJxlowpVKBQqEotJ6EISGEEELYNZkmE0IIIYRdkzAkhBBCCLsmYUgIIYQQdk3CkBBCCCHsmoQhIYQQQtg1CUNCCCGEsGsShoQQQghh1yQMCSGEEMKuSRgSQgghhF2TMCSEEEIIuyZhSAghhBB2TcKQEEIIIeyahCEhhBBC2DUJQ2UgOjqaRo0acfjw4fLuilVcvXqV0aNH06ZNG9q0acPLL79MbGxseXfLKm7dusUrr7xC586dadu2LSNGjODcuXPl3a0yMWPGDKZPn17e3Sg1Op2OBQsWEB4eTvPmzRk+fDhXrlwp725Z3ZIlSxg8eHB5d8NqUlJSePPNN+natSutWrXi6aefJjo6ury7ZRWJiYlMmTKFDh060LJlS0aNGsX58+fLu1tWd+nSJVq2bMnGjRtLrU0JQ1aWnp7O1KlT0el05d0Vq8jNzWXo0KEAfPfdd3z99dckJCTw4osvotfry7dzpSwvL49Ro0aRmJjI559/ztq1a/Hw8GDIkCEkJSWVd/esRqvV8tFHH7F+/fry7kqpWrJkCevWreO9997j+++/R6FQMHLkSPLy8sq7a1bz1VdfsWDBgvLuhlW98sor/PXXX8ybN4/169cTFhbGiBEjuHDhQnl3rdSNHj2aa9eusXz5ctavX4+zszNDhw4lOzu7vLtmNfn5+UyePJmsrKxSbVfCkJW9/fbbVK9evby7YTU3b96kadOmzJ49m/r169OoUSOGDh3KmTNnSE5OLu/ularo6GhiYmL4+OOPadKkCfXr1+fjjz8mKyuLPXv2lHf3rOLChQs8/fTTbNq0iapVq5Z3d0pNXl4eX375JePHj6dbt240bNiQ+fPnExcXx6+//lre3St1cXFxvPDCC3z22WfUrl27vLtjNVeuXOHAgQO89dZbtGnThjp16jBjxgyCgoLYtm1beXevVCUnJxMSEsK7775L06ZNqVu3LmPGjCEhIcGm71YvXLgQNze3Um9XwpAVbd68mWPHjvH666+Xd1espnbt2nz22Wf4+voCcP36ddauXUtYWBg+Pj7l3LvSVb9+fb744guCgoKMyvV6PampqeXUK+uKjIykUaNGbNu2jZCQkPLuTqk5c+YMmZmZdOjQwVDm6elJ48aNiYqKKseeWcc///yDl5cXW7ZsoXnz5uXdHavx8fHhiy++oEmTJoYyhUJhk39GfXx8mDdvHvXr1wfg9u3brFy5kuDgYOrVq1fOvbOOqKgovv/+ez766KNSb1td6i0K4E4omD17NkuWLLFKiq2Ihg8fzoEDB/Dy8mL16tUoFIry7lKpCggIoFu3bkZla9asITc3l86dO5dTr6zr6aefLu8uWMXdNW1VqlQxKg8MDOTWrVvl0SWr6tmzJz179izvblidp6enyZ/RnTt3cvXqVbp06VJOvbK+N954gx9++AFHR0eWLl2Kq6treXep1KWlpTF16lRmzpxp8ue2NEgYKoHr169z3333mb2+d+9epk6dypNPPkmbNm24fv16GfaudBU21v379xMQEADAlClTmDhxIkuXLmXo0KFs2rTJKr9praU4YwX45ZdfmD9/PoMHD6Zhw4Zl0cVSVdzx2pK7ayocHR2Nyp2cnGzuDoI9O3LkCK+//jr33XefTYfBIUOG8OSTT/Ldd98xduxYw915W/L222/TokULBgwYYJX2JQyVQFBQEDt27DB7/ccffyQrK4vx48eXYa+so7Cx3p0eA2jUqBEA8+fPp3v37mzYsIFx48ZZvY+lpThj/e6773j33Xfp27cvr732Wll0r9QVZ7y2xtnZGbizduju/8OdBwJcXFzKq1uiFO3evZvJkyfTvHlz5s2bV97dsaq702Lvvvsux48f55tvvuGDDz4o516Vnk2bNhEdHc3WrVut9j0kDJWAg4MDdevWNXt948aNxMfH0759ewDDU1UjR46kXbt2rFixokz6WRoKG+uNGzc4efIkffr0MZS5uLgQEhJCfHx8WXSx1BQ21rvmzJnD8uXLGTx4MDNmzKi004FFHa8tunvHMj4+nho1ahjK4+PjK+VdPmHsm2++Yfbs2fTu3Zs5c+aY3AG0BYmJiURERPDggw+iUqkAUCqV1K1bt9L93VuYDRs2kJiYSPfu3Y3K33rrLVauXMn27dvv+XtIGLKCr7/+Go1GY/g6Li6OwYMH89577xkCkq04ffo0EyZM4NdffzX8UElLS+PSpUsMHDiwnHtX+j755BNWrFjB1KlTGTFiRHl3R5RQw4YNcXd35/Dhw0a/b0+dOsVzzz1Xzr0T92Lt2rW8++67DB48mNdffx2l0jafE4qPj+fVV1/Fz8+Pjh07AnceOz916pTNTQnOmTOHnJwco7L777+fCRMm0Ldv31L5HhKGrKBatWpGX99N7UFBQSZPIlV2Xbt2JTQ0lKlTp/LGG2+g1+v55JNP8PHx4dFHHy3v7pWqw4cPs2LFCgYPHszAgQNJSEgwXHN1dbWbhfK2wNHRkeeee445c+bg6+tLtWrV+OSTTwgODqZ3797l3T1RQpcuXeL999+nd+/evPjiiyQmJhquOTs74+HhUY69K10NGzakS5cuzJo1i/feew9PT0+WLVtGWlqaYe83W2Hu56afn5/Jz9uSkjAk7omjoyMrVqzgo48+YsSIEeTl5dGlSxc+/PBD3N3dy7t7peruPiVff/01X3/9tdG1cePG2cQaMXsyYcIENBoNM2fOJCcnh7Zt27Jy5UqbnFKxF7t27SI/P59ff/3VZL+oRx55hA8//LCcelb6FAoFn376KXPnzmXSpEmkp6fTpk0bvv32W5vaE6ysKPS2tk2wEEIIIUQx2OZkqhBCCCFEEUkYEkIIIYRdkzAkhBBCCLsmYUgIIYQQdk3CkBBCCCHsmoQhIYQQQtg1CUNCCCGEsGsShoQQlYJsiWaZvD9ClJyEISFsyPTp0wkNDbX4391zi6ZPn15pzjD68ccf+eijjwxfb9y4kdDQUK5fv15q3+P69euEhoaycePGUmuzLMTGxvLiiy9y48YNQ1nPnj2ZPn06UHnHJURZkuM4hLAhY8aM4amnnjJ8vWTJEk6dOsWiRYsMZZXxuImlS5fSrl07w9fdu3fn+++/JzAwsBx7VTEcPHiQP/74gzfeeMNQtmjRIps7DkcIa5IwJIQNqVGjhuEUdgBfX18cHR1p0aJF+XXKCnx9ffH19S3vblRYjRs3Lu8uCFGpyDSZEHZu48aN9OnTh6ZNmzJw4ED27t1rdP3mzZu88sortGvXjubNmzNkyBBOnTplVCc9PZ0PPviAXr160bRpU/r378/69euN6vTs2ZP333+fIUOG0KpVK958800AUlJSePPNN+nUqRNNmzbliSeeICIiwuh1N27c4KeffjJMjRU0TXbgwAGeffZZWrZsSZcuXXjzzTdJTU01XI+KimLEiBG0bduWJk2a0LNnTxYuXIhOpyvW+7V27Vr69OlDs2bNePbZZzl48CChoaEcPnzY8H4WNIX376krgKSkJGbNmkWPHj1o0qQJ7dq1Y+zYsUavGzx4MDNmzOCLL76ge/fuNG3alKeeeoq//vrL8L1ee+01AO677z5D+//9Xv9VlM90x44dDBw4kGbNmtGhQwcmT55MfHx8sd4rISoLCUNC2LFbt27xxRdfMHHiRBYsWIBer2f8+PEkJiYCd35gP/XUU/zzzz+88cYbzJ07F51Ox7PPPsuFCxcAyMnJ4ZlnnmHLli0MHz6cJUuW0Lp1a2bMmMGyZcuMvt+3335LaGgoCxcu5KGHHiI3N5chQ4bw22+/8fLLL7No0SKCg4N54YUXDIFo0aJFBAQE0K1bN7NTY3/++ScvvPAC3t7ezJ8/nylTprBnzx4mTJgAwJkzZxg6dKjh+tKlS2nVqhWLFi1i+/btRX6/vv76a2bNmkV4eDiLFy+madOmvPzyy8V+3/V6PS+++CIHDhzg1VdfZeXKlYwZM4aDBw8aQuJdu3bt4rfffmPmzJnMmzeP27dvM2HCBLRaLd27d2f06NGG92nMmDGFfu+ifKZHjhxh8uTJ3H///SxfvpzXXnuNQ4cO8eqrrxZ7rEJUBjJNJoQd0+l0LF68mLp16wLg5OTEsGHDOH78OPfddx+rV68mJSWF7777jmrVqgHQtWtX+vbty2effcaCBQvYuHEjMTExrF27ltatWwMQHh6ORqNhyZIlPPXUU3h7ewMQGBjI9OnTUSrv/Dvshx9+4MyZM/zwww80b97c0P7gwYOZM2cOGzZsoHHjxjg6OuLr62t2um/BggU0bNiQxYsXG8qcnZ2ZN28ecXFxnDlzhk6dOvHJJ58Yvnfnzp35448/iIqKYsCAAUV6r5YuXUqfPn2YOXOmYZwZGRn8+OOPxXrf4+PjcXFxYdq0abRp0waA9u3bc/36ddatW2dUV6PRsHLlSsMaoMzMTKZNm8bp06dp0qSJYVq0UaNGhISEFPq9i/KZHjlyBCcnJ0aOHImTkxMA3t7enDhxAr1ej0KhKNZ4hajo5M6QEHbMx8fHEIQAqlevDtyZ9gKIiIigUaNGBAUFodFo0Gg0KJVKunbtysGDBwGIjIykWrVqhiB018CBA8nNzTVM6QDUrVvXEEbuth8QEEBYWJihfa1WS48ePTh58qTRNJc5OTk5/PPPP/Tq1cuovE+fPuzatYugoCAefvhhli9fTn5+PufOnWP37t0sXLgQrVZLfn5+kd6rS5cukZiYyH333WcyzuIKCgpizZo1tGnThps3bxIREcE333zD0aNHTfpTr149o8XQQUFBAGRnZxf7+0LRPtO2bduSk5PDgAEDmD9/PkeOHKFLly6MGzdOgpCwSXJnSAg75urqavT13R90d9fRpKSkcOXKFcLCwgp8fXZ2Nqmpqfj7+5tcu1uWlpZmUnZXSkoKCQkJZttPSEjAy8vL4hhSU1PR6/X4+fmZrZOTk8O7777L5s2b0Wg0hISE0LJlS9RqdZH350lJSQEwWbh9N5wU15YtW5g3bx63bt3C29ubhg0b4uzsbFLPxcXF6Ou7YbK4a53uKspn2rJlS7744gu++uorVq5cybJlywgICGDkyJEMGTKkRN9XiIpMwpAQwiwPDw/atWvH1KlTC7zu6OiIl5cXV65cMbmWkJAA3Ln7ZKn9WrVqMWfOnAKvF2Xax93dHYVCQVJSklF5Xl4eERERNGvWjHnz5rFr1y4+/fRTOnXqZAiBHTt2LLT9u+6O4/bt20bld0PSXf8NlHdlZmYa/j86Oppp06bx3HPPMWLECIKDgwH4+OOPOXLkSJH7VBJF+UzhzhRgeHg42dnZHDp0iDVr1vD+++/TokULw5SmELZCpsmEEGa1a9eOS5cuUbt2bZo2bWr4b8uWLfz444+oVCratm3LjRs3TH6Ib9myBQcHB5o1a2ax/Vu3buHn52fUfkREBCtWrEClUgEYTa39l5ubG40aNeK3334zKt+/fz+jRo0iNjaWI0eO0L59e3r16mUIQidPniQpKanId1hq165NlSpV2LFjh1H5nj17jL6+O6V169YtQ9nFixeNQtOxY8fQ6XRMmDDBEIS0Wq1hmqo4d30svTcFKcpn+tFHH/HYY4+h1+txcXGhR48eTJs2zWRcQtgKCUNCCLOGDh2KTqdj6NCh7Nixg4iICN544w3WrFlDnTp1ABg0aBD16tVj3LhxfPfdd+zfv5933nmHDRs28OKLL+Lp6Wm2/UGDBlG1alWGDRvGTz/9xKFDh5g3bx7z588nMDAQBwcHADw9PTl16hSRkZHk5OSYtDNhwgT++ecfJk2axN69e9m0aRNvvfUWPXr0oFGjRjRr1oz9+/fz3XffERkZyZo1axg5ciQKhaLIa28UCgVTp05l3759vPbaa+zbt48lS5bw5ZdfGtXr0KEDLi4ufPjhh/z555/s2LGDcePGGRaRA4aA+M4773Do0CF++eUXhg0bxpkzZwDIysoqUp/uvjcAv/76q+FpMEuK8pl27NiRkydPMn36dA4cOMAff/zBe++9h7e3Nx06dChy34SoLGSaTAhhVlBQEOvWrWPu3Lm8/fbb5ObmUqtWLWbPns1jjz0G3FnT8vXXXzN37lwWLFhARkYGderUMapjjqurK99++y1z587lk08+IT09nWrVqvHqq68yfPhwQ73hw4fz/vvvM2LECFatWmXSTo8ePfj8889ZuHAhY8eOxcfHhwcffJCJEycCd44eyc/P59NPPyUvL4+QkBBGjx7N+fPn2bNnD1qttkjvR9++fVGpVCxcuJCtW7fSqFEjXn31VT744ANDHQ8PDxYsWMDcuXMZO3Ys1apVY9y4cWzatMlQp3379rz55pusWrWKn3/+GX9/f9q3b8+iRYsYO3YsR44coVu3bkXqU/v27enUqRNz584lIiKCL774wmL9onymXbt2Zc6cOXz55ZeGRdOtW7dmzZo1RqFOCFuh0MvpfkIIUWKHDx/m+eefZ82aNbRv3768uyOEKAGZJhNCCCGEXZMwJIQQQgi7JtNkQgghhLBrcmdICCGEEHZNwpAQQggh7JqEISGEEELYNQlDQgghhLBrEoaEEEIIYdckDAkhhBDCrkkYEkIIIYRdkzAkhBBCCLv2fzJvWHcL16w/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# first simulate the sampling distribution of the mean for 10,000 samples\n", "nSamples = 10000\n", "samplesize = 100\n", "m=np.empty(nSamples) # make an array to store the means\n", "binwidth=0.1\n", "\n", "for i in range(nSamples):\n", " sample = UKBrexdex.sample(n=samplesize, replace=False)\n", " m[i]=sample.score.mean()\n", "\n", "# Now make the Q-Q plot\n", "stats.probplot(m, plot=plt)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d4c35390", "metadata": {}, "source": [ "* Try changing `samplesize` ($n$) back to smaller values ($n=2,3,4$) or larger value such as ($n=10,100$) in the code block above and remaking the Q-Q plot.\n", "\n", "You should see the tails of the sampling distribution (in both the histogram and the Q-Q plot) start to match the normal distribution\n", "\n", "* Try setting $n=2$ in the code block above and remaking the Q-Q plot. \n", "\n", "You should see the funny three-peak histogram - how is the shape of the histogram reflected in the Q-Q plot?" ] }, { "cell_type": "code", "execution_count": null, "id": "6c98deeb", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.7" } }, "nbformat": 4, "nbformat_minor": 5 }