{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4304eb40-80fd-46b5-a6ef-6e7a84d33dbe",
   "metadata": {},
   "source": [
    "## Monte Carlo Tests\n",
    "To motivate resampling methods, suppose we wish to infer from measurements whether the weights of adult human males in a medical study are normally distributed [[1](https://www.jstor.org/stable/2333709https://www.jstor.org/stable/2333709)]. This may be prerequisite to performing other statistical tests, many of which are developed under the assumption that samples are drawn from a normally distributed population (although some tests are quite robust to deviations from this assumption, as we will see later). The weights are recorded in the array `x` below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "02824778-802d-4d9d-8cd5-0a69cc028d13",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236])  # weights (lbs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe14279d-496f-45fe-9ce0-937ff5844a3b",
   "metadata": {},
   "source": [
    "One way of testing for departures from normality, chosen based on its simplicity rather than its sensitivity, is the [Jarque-Bera test [2]](https://www.sciencedirect.com/science/article/abs/pii/0165176580900245) implemented in SciPy as `scipy.stats.jarque_bera`. The test, like many other hypothesis tests, computes a *statistic* based on the sample and compares its value to the distribution of the statistic derived under the *null hypothesis* that the sample is normally distributed. If the value of the statistic is extreme compared to this *null distribution* - that is, if there is a low probability of sampling such data from a normally distributed population - then we have evidence to reject the null hypothesis.\n",
    "\n",
    "The statistic is calculated based on the skewness and kurtosis of the sample as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8475312e-7810-4f61-8fed-7fe8297afe28",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.982848237344646"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "\n",
    "def statistic(x):\n",
    "    # Calculates the Jarque-Bera Statistic\n",
    "    # Compare against `scipy.stats.jarque_bera`:\n",
    "    # https://github.com/scipy/scipy/blob/4cf21e753cf937d1c6c2d2a0e372fbc1dbbeea81/scipy/stats/_stats_py.py#L1583-L1637\n",
    "    n = len(x)\n",
    "    mu = np.mean(x, axis=0)\n",
    "    x = x - mu  # remove the mean from the data\n",
    "    s = stats.skew(x)  # calculate the sample skewness\n",
    "    k = stats.kurtosis(x)  # calculate the sample kurtosis\n",
    "    statistic = n/6 * (s**2 + k**2/4)\n",
    "    return statistic\n",
    "\n",
    "stat1 = statistic(x)\n",
    "stat2, _ = stats.jarque_bera(x)\n",
    "np.testing.assert_allclose(stat1, stat2, rtol=1e-14)\n",
    "stat1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a38c293-1414-4d5c-8c6f-84d9d2e77e93",
   "metadata": {},
   "source": [
    "Note that the value of the statistic is unaffected by the scale and location of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "138b5da5-d953-4630-97c4-ebd8ad6420c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "old_location = np.mean(x)\n",
    "old_scale = np.std(x)\n",
    "x_new = (x - old_location) / old_scale  # make location 0 and scale 1\n",
    "stat3 = statistic(x_new)\n",
    "np.testing.assert_allclose(stat1, stat3, rtol=1e-14)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69fad794-1d58-49c9-9c75-09f642cdfbf0",
   "metadata": {},
   "source": [
    "Consequently, it can be shown that large samples drawn from a normal distribution with any mean and variance will produce statistic values that are distributed according to the chi-squared distribution with two degrees of freedom. We can check this numerically by drawing 1000 samples of size 500 from a standard normal distribution and computing the value of the statistic for each sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9e2c2e75-1f3b-477c-9b6a-20d6bbfd6efd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/RElEQVR4nO3deXhTZfbA8e9J0tKWfZdNWWRvaQuF4iCbC5sIo6yOIoiI6CDuA+MGLqOOMC64DoLAb0ZBRUHGDVRAQdlaKJvsilJQtkJZSpc07++PpDEtaZuWpmna83mePEnufe/NyaXk3Pvee88rxhiUUkpVXJZAB6CUUiqwNBEopVQFp4lAKaUqOE0ESilVwWkiUEqpCs4W6ACKqk6dOqZp06aBDkMppYJKYmLicWNMXW/zgi4RNG3alISEhECHoZRSQUVEfslvnnYNKaVUBaeJQCmlKjhNBEopVcEF3TmC4sjKyiI5OZn09PRAh6KUUn4VFhZG48aNCQkJ8XmZCpEIkpOTqVq1Kk2bNkVEAh2OUkr5hTGGEydOkJycTLNmzXxerkJ0DaWnp1O7dm1NAkqpck1EqF27dpF7P/yaCESkn4jsFpF9IjLFy/xeIpIqIkmuxxN+jMVfq1ZKqTKjOL91fusaEhEr8DpwLZAMbBSRpcaYH/M0XW2MGejreu0OLZutlFIlyZ9HBF2AfcaYn4wxmcBCYPDFrnTPkTME6xgKixcvRkTYtWtXQD7/wIEDvPfee4W2S0pK4vPPP3e/X7p0Kc8//7zPn9O0aVOGDBnifr9o0SLGjBlT6HJVqlRxxxkZGXnB/AMHDhAeHk5sbCxt27alS5cuzJ8/3+c4836vvBISEpg0aRIA06ZNY8aMGYXG7Onll18mLS3N/X7AgAGcOnWqSOtQKhD8mQgaAQc93ie7puV1hYhsEZEvRKS9txWJyHgRSRCRhGyHYc+Rs/6I1+8WLFjAlVdeycKFCwPy+cVNBIMGDWLKlAt69gqUkJDAjh07ihxjYVq0aMHmzZvZuXMnCxcu5KWXXmLu3Lk+xVlQIrDb7cTFxTFz5sxix5Y3EXz++efUqFGj2OtTqrT4MxF466jKuyu/CbjMGBMNvAos8bYiY8wsY0ycMSYOYN1PJ0oyzlJx9uxZvv/+e+bMmZMrEfz222/06NGDmJgYIiMjWb16NXPmzOH+++93t3n77bd54IEHOHDgAG3atGHcuHFERkZy88038/XXX9OtWzdatmzJhg0bAOfe7KhRo7jqqqto2bIlb7/9NgBTpkxh9erVxMTE8NJLL5Gens5tt91GVFQUsbGxrFy5kszMTJ544gnef/99YmJieP/995k3bx4TJ04E4MiRI9xwww1ER0cTHR3NDz/84PX7PvTQQzz77LMXTM+7px0ZGcmBAweKtU2bN2/Oiy++6P7x9ozzww8/JDIykujoaHr06OH1e02bNo3x48fTp08fbr31VlatWsXAgX/0Um7ZsuWCbZi3zcSJE5k3bx4zZ87k8OHD9O7dm969ewPOI6Pjx48D8OKLLxIZGUlkZCQvv/wy4EzMbdu25Y477qB9+/b06dOH8+fPF2tbKHUx/Hn5aDLQxON9Y+CwZwNjzGmP15+LyBsiUscYczy/lYZYLaz76QSj/9S0WEE9+b8d/Hj4dOENi6Bdw2pMvd7rwYzbkiVL6NevH61ataJWrVps2rSJjh078t5779G3b18effRRsrOzSUtLo2PHjnTo0IEXXniBkJAQ5s6dy7///W8A9u3bx4cffsisWbPo3Lkz7733HmvWrGHp0qU8++yzLFmyBICtW7eybt06zp07R2xsLNdddx3PP/88M2bM4NNPPwXgX//6FwDbtm1j165d9OnThz179vDUU0+RkJDAa6+9Bjh/YHNMmjSJnj17snjxYrKzszl71vvR2fDhw3njjTfYt2/fxWzaQnXs2NFrV9tTTz3FsmXLaNSoEadOnSI0NPSC7zVt2jQSExNZs2YN4eHhrFq1Ktc6vG3D/EyaNIkXX3yRlStXUqdOnVzzEhMTmTt3LuvXr8cYQ3x8PD179qRmzZrs3buXBQsW8PbbbzN8+HA++ugjbrnllovfMEoVgT+PCDYCLUWkmYiEAiOBpZ4NROQScZ3iFpEurngK3N2vUsnG+p9TcATZSeMFCxYwcuRIAEaOHMmCBQsA6Ny5M3PnzmXatGls27aNqlWrUrlyZa666io+/fRTdu3aRVZWFlFRUQA0a9aMqKgoLBYL7du35+qrr0ZEiIqKyrVnPXjwYMLDw6lTpw69e/d2Hy14WrNmDaNGjQKgTZs2XHbZZezZs6fA77FixQruuusuAKxWK9WrV/fazmq18vDDD/Pcc88VbUMVUX7ni7p168aYMWN4++23yc7Oznf5QYMGER4e7nWeL9vQF2vWrOGGG26gcuXKVKlShRtvvJHVq1cDzn/PmJgYADp16lTsoyOlLobfjgiMMXYRmQgsA6zAO8aYHSIywTX/LWAocJeI2IHzwEhTyJngypWspJzLZO/Rs7S+pGqR4ypsz90fTpw4wYoVK9i+fTsiQnZ2NiLCCy+8QI8ePfjuu+/47LPPGDVqFA8//DC33nor48aN49lnn6VNmzbcdttt7nVVqlTJ/dpisbjfWywW7Ha7e17eS8i8XVLm75Puo0aN4rnnnqN9+z+2uc1mw+FwuN9f7N3emzdvpm3bthdMf+utt1i/fj2fffYZMTExJCUleV2+cuXK+a7b2zYsTvwFbWfPf0+r1apdQyog/HofgTHmc2NMK2NMC2PMP1zT3nIlAYwxrxlj2htjoo0xXY0x3jucPVQOdeautfvz7T0qcxYtWsStt97KL7/8woEDBzh48CDNmjVjzZo1/PLLL9SrV4877riD22+/nU2bNgEQHx/PwYMHee+997jpppuK/JmffPIJ6enpnDhxglWrVtG5c2eqVq3KmTNn3G169OjBu+++C8CePXv49ddfad269QXtPF199dW8+eabAGRnZ3P6dP7dbCEhIdx///3uPnFw9pvnfMdNmzbx888/F/m75Thw4AAPPfQQ99xzzwXz9u/fT3x8PE899RR16tTh4MGDBX4vb7xtw8suu4wff/yRjIwMUlNT+eabb9zt81t/jx49WLJkCWlpaZw7d47FixfTvXv34n1ppfwg6O4sDrVZaFQjnHU/pQQ6FJ8tWLCAG264Ide0IUOG8N5777Fq1SpiYmKIjY3lo48+4t5773W3GT58ON26daNmzZpF/swuXbpw3XXX0bVrVx5//HEaNmxIhw4dsNlsREdH89JLL3H33XeTnZ1NVFQUI0aMYN68eVSqVInevXvz448/uk+qenrllVdYuXIlUVFRdOrUqdArg26//fZcRypDhgwhJSWFmJgY3nzzTVq1alWk77V//3735aPDhw/nnnvuyXXElOPhhx8mKiqKyMhIevToQXR0dIHfyxtv27BJkyYMHz6cDh06cPPNNxMbG+tuP378ePr37+8+WZyjY8eOjBkzhi5duhAfH8+4ceNyLadUoEmwXZMfFxdnek2ewzc7j5D42LVYLIXfRbdz506v3Qdl3cCBA7n//vu5+uqri7TctGnTqFKlCg899JCfIlNKlWXefvNEJDHnysu8grLoXNfmtVmUmMyeo2doc0m1P2ZM837ikhHrSiewEnLq1Cm6dOlCdHR0kZOAUkoVVVAmgvhmtQBYu/9E7kRQTtSoUaPQq3cKMm3atJILRilV7gXdOQKAJrUiaFwzPChvLFNKqbImKBMBwBXNawfl/QRKKVXWBG0i6Nq8NqfSsth9xPfLAZVSSl0oaBNBfPM/zhMopZQqvqBNBI1rRtCkVjHPE0yrXrIPH4iIu5wDOKtd1q1bN1cBs6I4deoUb7zxRpGXO3v2LHfeeSctWrSgffv29OjRg/Xr1xdpHb169SIhIaFI7ePi/rhqLSEhgV69ehW6nGfRtpwS1XlZrVZiYmJo37490dHRvPjii+47fz3LSntTWDXWw4cPM3ToUCB3QTtfzZs3j8OH/yivNW7cOH78Me9wHMWzZMkSnnrqKffn1K1bl5iYGGJiYpg9e7a73fz582nZsiUtW7bMVbL7559/Jj4+npYtWzJixAgyMzN9/uzibAt/WLJkSa7t+cQTT/D111+X+Of069ePQ4cOlfh6C+JLGfTXXnvNXXn3YgVtIoDgOk9QuXJltm/f7i4h8NVXX9Gokbeq3L4pbiIYN24ctWrVYu/evezYsYN58+a5f2x9UVDdnoIcPXqUL774oljLFiQ8PJykpCR27NjBV199xeeff86TTz4JUGhZ6YISgd1up2HDhixatKjYseVNBLNnz6Zdu3bFXp+nF154gbvvvtv9fsSIESQlJZGUlMS4ceMASElJ4cknn2T9+vVs2LCBJ598kpMnTwIwefJk7r//fvbu3UvNmjWZM2dOicRVXJ43HfoqbyJ46qmnuOaaa0oyLM6fP09KSspF/V/1l7Fjx15U2XRPQZ0IujavTer5LHb9HhznCfr3789nn30GOO829iwdkZKSwp///Gc6dOhA165d2bp1K+DcMxg7diy9evWiefPm7n/4KVOmsH//fmJiYnj44YcBmD59Op07d6ZDhw5MnTr1gs/fv38/69ev55lnnsFicf7TN2/e3F1V889//jOdOnWiffv2zJo1y71clSpVeOKJJ4iPj2ft2rW51rlgwQL3HbyTJ0/O97s//PDDPPPMMxdMz7t3OXDgwAuqgPqqXr16zJo1i9deew1jTK6S0d9++617jzk2NpYzZ85cUJZ73rx5DBs2jOuvv54+ffpcMEDOwYMH6devH61bt3Ynm7xtZsyYwbRp01i0aBEJCQncfPPNxMTEcP78+VxHUvlttypVqvDoo48SHR1N165dOXLkyAXfc8+ePVSqVOmCKqd5LVu2jGuvvZZatWpRs2ZNrr32Wr788kuMMaxYscJ9tDN69Gh31VpP3sqUF7Qtzp07x3XXXUd0dDSRkZHuu7cTExPp2bMnnTp1om/fvvz222+A80jxkUceoWfPnvzjH/+gadOm7qO5tLQ0mjRpQlZWFm+//TadO3cmOjqaIUOGkJaWxg8//MDSpUt5+OGHiYmJYf/+/YwZM8aduL/55htiY2OJiopi7NixZGRkAM6jzKlTp9KxY0eioqLclWu9/X2As+x4ztHrlClTaNeuHR06dHDfrPm///2P+Ph4YmNjueaaa9z/XtOmTWP06NH06dOHpk2b8vHHH/O3v/2NqKgo+vXrR1ZWljueyZMn06VLF7p06eK1Wu/+/fvp168fnTp1onv37u6YIyIiaNq0abGLIXoK6kQQ37w2EDzjE4wcOZKFCxeSnp7O1q1biY+Pd8+bOnUqsbGxbN26lWeffZZbb73VPW/Xrl0sW7bMvVeXlZXF888/T4sWLUhKSmL69OksX76cvXv3smHDBpKSkkhMTOS7777L9fk7duwgJiYGq9XqNb533nmHxMREEhISmDlzJidOOLfruXPniIyMZP369Vx55ZXu9ocPH2by5MmsWLGCpKQkNm7c6PUHBeCKK66gUqVKuX5M/KF58+Y4HA6OHj2aa/qMGTN4/fXXSUpKYvXq1YSHh/P888/TvXt3kpKS3OM/rF27lvnz57NixYoL1r1hwwbeffddkpKS+PDDDwvsHhs6dChxcXHu9p4VTgvabufOnaNr165s2bKFHj16uMdB8PT999/TsWPHXNM++ugjOnTowNChQzl40Dke1KFDh2jS5I9K8I0bN+bQoUOcOHGCGjVqYLPZck3P6/XXXwecZcoXLFjA6NGj3UX2vG2LL7/8koYNG7Jlyxa2b9/u/sG75557WLRoEYmJiYwdO5ZHH33U/RmnTp3i22+/ZerUqURHR/Ptt98Czh/Yvn37EhISwo033sjGjRvZsmULbdu2Zc6cOfzpT39i0KBBTJ8+naSkJFq0aOFeZ3p6OmPGjOH9999n27Zt2O12d30sgDp16rBp0ybuuusud/eLt78PgC+++IJ+/fqRkpLC4sWL2bFjB1u3buWxxx4D4Morr2TdunVs3ryZkSNH8sILL7g/Z//+/Xz22Wd88skn3HLLLfTu3Ztt27YRHh7u3iEEqFatGhs2bGDixIncd999F/w7jB8/nldffZXExERmzJiR60gwLi7OXcn2YgR1ImhUI5xLa0WwNkgSQYcOHThw4AALFixgwIABueZ5loS+6qqrOHHiBKmpqQBcd9117j3AevXqed1LXL58OcuXLyc2NtZdo3/v3r1Fim/mzJnuPdGDBw+6l7darbmGnsyxceNGevXqRd26dbHZbNx8880XJB9Pjz32mNejgpLmrWxKt27deOCBB5g5cyanTp1y/wjmlbMHnd+82rVrEx4ezo033siaNWuKFV9B2y00NNR9FJNfWerffvuNunXrut9ff/31HDhwgK1bt3LNNdcwevRowPt2EJF8p+dVUJlyb9siKiqKr7/+msmTJ7N69WqqV6/O7t272b59O9deey0xMTE888wzJCcnuz9jxIgRuV7nHEUsXLjQPW/79u10796dqKgo3n333ULrW+3evZtmzZq561iNHj0619/ljTfeCOTevvn9fXz//fdceeWVVKtWjbCwMMaNG8fHH39MREQEAMnJyfTt25eoqCimT5+eK7b+/fsTEhJCVFQU2dnZ9OvXD+CCkvE5PQM33XTTBUfcZ8+e5YcffmDYsGHExMRw5513uo+owHkU7Nn9WFxBnQjAeZ5gQ5CcJwBn/fuHHnrogoqiBf3nzFuq2Ft/qjGGv//97+5+4n379nH77bfnatO+fXu2bNmSq4xyjlWrVvH111+zdu1atmzZQmxsrHvvLywszOtRRFHrVF111VWkp6ezbt0fJT9Kuiz1Tz/9hNVqpV69ermmT5kyhdmzZ3P+/Hm6du2a77jRgS5LHRIS4v6c/P6tw8PDc31O7dq13X8jd9xxB4mJiYBzTz/n6ACcP1oNGzakTp06nDp1yr3unOlFidPbtmjVqhWJiYlERUXx97//naeeegpjDO3bt3f/XW7bto3ly5e7l/Pc3oMGDeKLL74gJSWFxMRErrrqKgDGjBnDa6+9xrZt25g6dWqh27iwv8ucbeW5fb39ffz00080adKE0NBQbDYbGzZsYMiQIe5BpgDuueceJk6cyLZt2/j3v/+dKzbPEvGe/64FlYzPu10dDgc1atRwb7+kpCR27tzpnp+enp7veBpFEfSJoGuLWqSez2Ln7yU76pi/jB07lieeeMI90EwOz5LQq1atok6dOlSrln/5jLwlj/v27cs777zjHjHs0KFDF3SPtGjRgri4OKZOner+z7J3714++eQTUlNTqVmzJhEREezatSvXj3V+4uPj+fbbbzl+/DjZ2dksWLCAnj17FrjMo48+muvwuWnTpiQlJeFwODh48OBF9XceO3aMCRMmMHHixAv+Q+3fv5+oqCgmT55MXFwcu3btKnJZ6q+++oqUlBTOnz/PkiVL6NatG/Xr1+fo0aOcOHGCjIwM9+hvkH9Z6uJsN09t27bN1ZfsuYe4dOlSd7Gxvn37snz5ck6ePMnJkydZvnw5ffv2RUTo3bu3uz99/vz5DB48+ILPya9MeX7b4vDhw0RERHDLLbfw0EMPsWnTJlq3bs2xY8fce7pZWVn57tFXqVKFLl26cO+99zJw4ED3zseZM2do0KABWVlZ7ngK2r5t2rThwIED7m30n//8p9Dt6+3vI6dbCJx75qmpqQwYMICXX37ZPb5Famqq+0Sy51VZRZFzFPT+++9zxRVX5JpXrVo1mjVrxocffgg4k9yWLVvc8/fs2ZPrHFVxBWWtIU/xzXLOE6Tg85Az01L9Fk9hGjdunKvUdI5p06Zx22230aFDByIiIgr9o6pduzbdunUjMjKS/v37M336dHbu3On+Q6pSpQr//e9/L9gznj17Ng8++CCXX345ERER1K5dm+nTp9OhQwfeeustOnToQOvWrenatWuh36VBgwY899xz9O7dG2MMAwYM8PqD4mnAgAG5ujW6devmHnUtMjLygr7vwpw/f56YmBiysrKw2WyMGjWKBx544IJ2L7/8MitXrsRqtdKuXTv69++PxWJxl+UeM2ZMoeW+r7zySkaNGsW+ffv4y1/+4r4kNudEerNmzWjTpo27/ZgxY5gwYQLh4eG5DvmLs9089ejRgwcffBBjDCLCzJkzWbp0KTabjVq1armHFq1VqxaPP/44nTt3dseZ0+31z3/+k5EjR/LYY48RGxt7wdEjwN13382ECROIiorCZrO5y5Tnty2WLVvGww8/7N4DfvPNNwkNDWXRokVMmjSJ1NRU7HY79913X67BijyNGDGCYcOG5bpg4OmnnyY+Pp7LLruMqKgo94//yJEjueOOO5g5c2auq7vCwsKYO3cuw4YNw26307lzZyZMmFDgNvX29zF06FBeffVVwJmMBg8eTHp6OsYYXnrpJcD5/3bYsGE0atSIrl27Fmt8jYyMDOLj43E4HO6RCz29++673HXXXTzzzDNkZWUxcuRIoqOjAWfXlbcLQ4oqKMtQ5z1J13P6SlrWq8rsn71X6tw5Yl1QlqFWKj/33nsv119/fYlfLqmcMjIy6NatW5HulymOpk2bkpCQUOgVYN5s3ryZF198kf/85z8XzCtqGeqg7xqCnPMEJ8g2hY9NoFR58Mgjj5CWlhboMMqtSpUq+T0JXKzjx4/z9NNPl8i6gq9r6PDmC+7m7ZrdjYVZf2Vn6KVEyi8BCkyp0lO/fn0GDRoU6DDURfJ2VZivrr322hKLo1wcEcRbnGfR1znyv2sz2LrAlFKqOIrzW1cuEkEDSaGp/J5vIggLC+PEiROaDJRS5ZoxhhMnThAWFlak5YKvaygf3SzbWZLdjQxjo5Lkvva6cePGJCcnc+zYsQBFp5RSpSMsLIzGjRsXaZlykwh6W5J4N/saEhyt6WbNfZ1ySEgIzZo1C1BkSilVtpWLriGAP1l2EEoWKx0xgQ5FKaWCSrlJBBGSQbxlpyYCpZQqonKTCMDZPbTfNOJXR73CGyullALKYSIAWOWIDmwgSikVRMpVImhm+Z2m8rt2DymlVBGUq0QA0MuSxA+O9qSbkECHopRSQaHcJYLeliQyCGVtAXcZK6WU+kO5SwTxlp2EkcEq7R5SSimflLtEECZZdLPsYKUjBq0ooZRShfNrIhCRfiKyW0T2iciUAtp1FpFsERlaEp/by5LEr6Y+P5kGJbE6pZQq1/yWCETECrwO9AfaATeJyAUd9652/wSWldRn93JdRqpXDymlVOH8eUTQBdhnjPnJGJMJLAS8jcd3D/ARcNTLvGJpYjlOS0nmW72fQCmlCuXPRNAIOOjxPtk1zU1EGgE3AG8VtCIRGS8iCSKScCzNt47/XpYk1jvacs5UKlrUSilVwfgzEXgbNzLvr/jLwGRjTHZBKzLGzDLGxBlj4upG+DYcZW9LEpmE8IPD5yHtlVKqQvJnGepkoInH+8bA4Txt4oCFIgJQBxggInZjzJKL/fA4y24qc56VjhhKbkA3pZQqf/yZCDYCLUWkGXAIGAn8xbOBMcY9SICIzAM+LYkkABAq2Vxp2caq7BiMMbiSjVJKqTz81jVkjLEDE3FeDbQT+MAYs0NEJojIBH99rqfeliQOU4c9R86WxscppVRQ8usIZcaYz4HP80zzemLYGDOmpD+/l3UL2GHl7qO0vqRqSa9eKaXKhXJ3Z7GnS+QkbeUAK3eV2JWpSilV7pTrRADQ27KFhF9Ocjo9K9ChKKVUmVT+E4E1iWyHYc3e44EORSmlyqRynwhiZS/Vw0P4eueRQIeilFJlkl9PFpcFNnFwdcY3fL2pE5k7riJU8ty7Ni01MIEppVQZUe6PCACus67nNJX53hEV6FCUUqrMKTQRiMgMEQnqOg1XWrZRlXN85ogPdChKKVXm+HJEsAuYJSLrXTeDVfd3UCWtkti51pLI8uw4Mo010OEopVSZUmgiMMbMNsZ0A24FmgJbReQ9Eent7+BKknYPKaWUdz6dI3ANHtPG9TgObAEeEJGFfoytRDm7h9K0e0gppfLw5RzBizi7hwYAzxpjOhlj/mmMuR6I9XeAJUW7h5RSyjtfjgi2A9HGmDuNMRvyzOvih5j85jrrOlf3UGSgQ1FKqTLDl0RwszEmzXOCiHwDYIwJqovw/+ge6hroUJRSqszINxGISJiI1ALqiEhNEanlejQFGpZahCVIu4eUUupCBR0R3Akk4jxBvMn1OhH4BHjd/6H5xwD31UPaPaSUUlBAIjDGvOIaQewhY0wzj0e0Mea1UoyxRHW3bKUqaXyuVw8ppRRQcNfQVa6Xh0TkxryPUoqvxFUSO9do95BSSrkVVHSuJ7ACuN7LPAN87JeISsF11vUsdnTne0ckQXVXnFJK+UG+icAYM9X1fFvphVM6PLuHNBEopSo6X24ou1dEqonTbBHZJCJ9SiM4f/HsHsrKdgQ6HKWUCihf7iMYa4w5DfQB6gG3Ac/7NapSMMC6nlSq8P0+HblMKVWx+ZIIxPU8AJhrjNniMS1odbdsowppfL7tt0CHopRSAeVLIkgUkeU4E8EyEakKBH1/SphkcY1lE8t2HNHuIaVUheZLIrgdmAJ0dpWaCMXZPRT0rrOuI/V8lnYPKaUqNF/GI3AAR4B2ItIDaA/U8HNcpaK7ZRvVwmws2Xwo0KEopVTAFDp4vYj8ExgB/AjkjPxugO/8GFepCJMsBsc04oOEgzx5Povq4SGBDkkppUqdL11DfwZaG2MGGGOudz0G+TmuUjM8rgkZdgf/23I40KEopVRA+JIIfgLK7a5yZKNqtLmkKh8mHAx0KEopFRCFdg0BaUCSawyCjJyJxphJfouqFIkIw+Oa8NSnP7Lr99O0uaRaoENSSqlS5csRwVLgaeAH/ihFnejPoErbn2MbEWIVPkxIDnQoSilV6go9IjDGzBeRcOBSY8zuUoip1NWqHMq17eqzePMhJvdrQ6jNl/yolFLlgy+1hq4HkoAvXe9jRGSpn+MqdcPimpByLpMVu44EOhSllCpVvuz6TsM5SP0pAGNMEtDMl5WLSD8R2S0i+0Rkipf5g0Vkq4gkiUiCiFzpc+QlZVp1mFadHu+14hJO8MF7s93TlFKqIvAlEdi9DFJvCltIRKw4h7TsD7QDbhKRdnmafQNEG2NigLHAbB/i8QurGIZYV7PKEcMRUyNQYSilVKnzJRFsF5G/AFYRaSkir+I8cVyYLsA+Y8xPxphMYCEw2LOBMeasMSYnqVTGhwTjT8Os3+LAwkfZ3QMZhlJKlSpfEsE9OMtKZAALgNPAfT4s1wjwvDg/2TUtFxG5QUR2AZ/hPCoImKaWI3SRnXyY3QsT0JSklFKlx5daQ2nGmEeNMZ2NMXGu1+k+rNtbqeoLfl6NMYuNMW1w3sH8tNcViYx3nUNIOJbm31/o4bZV/GwakGBa+/VzlFKqrCgwEYjIaNeIZOdcjwQRudXHdScDTTzeNwbyreNgjPkOaCEidbzMm+VKQnF1I/w7FMIAywYqc54Psnv69XOUUqqsyDcRuH7w7wMeBBri7Nb5G3Cvj8lgI9BSRJqJSCgwEufNaZ6fcbmIiOt1R5wlrk8U43uUmAjJ4HrrWj7L7srZDHsgQ1FKqVJR0BHB3cANxpiVxphUY8wpY8wKYIhrXoGMMXZgIrAM2Al8YIzZISITRGSCq9kQnCejk3BeYTTC4+RxwAyzfksaYXy+VUcvU0qVf5Lf766I/GiMyXu5Z6Hz/C2uodUkjK/i188wBq7JnE6NSyP56K4/+fWzlFKqNIhIojEmztu8go4IzhdzXtATgeHWb0n85ST7jp4JdDhKKeVXBSWCtq67fvM+tgFtSivAQLnR+h2hVgvzf/gl0KEopZRfFVR0rm2pRVEG1ZXTDIppyKLEZB7s04oaEaGBDkkppfwi30RgjKnwu8Jjt93CoqznWfDsWO6y/e+PGdPyVtxQSqngpfWWC9DO8it/smxnvr0PWcYa6HCUUsovNBEU4nbrF/xObT53dAl0KEop5Re+jEcwUEQqbMLobUmiuRzmHfsArT+klCqXfPmBHwnsFZEXRKTCnUC2iOE265dsMS1INK0CHY5SSpU4X4rO3QLEAvuBuSKy1lUErqrfoysjhlhXU52zzLH3D3QoSilV4godsxjAGHNaRD4CwnHWH7oBeFhEZhpjXvVjfGVChGRwk3UFs7IHctBRlyYFjV6mVxQppYKML+cIBonIYmAFEAJ0Mcb0B6KBh/wcX5lxq205gmFedt9Ah6KUUiXKl3MEQ4GXjDEdjDHTjTFHwTlOAQEeSKY0NZQUBljW8352L86Y8ECHo5RSJcaXRPCba6wANxH5J4Ax5hu/RFVG3W77grNE8KGOVaCUKkd8SQTXeplWIc+axlj200l2Mze7H9nGvwPkKKVUaSloYJq7cgrM5Sk69zOwtfRCLFtut33BQVOPrxydAh2KUkqViIKuGnoP+AJ4DpjiMf2MMSbFr1GVYX0sCTTiGO/Y+9PPmhDocJRS6qIV1DVkjDEHgL8CZzweiEgt/4dWNtnEwW22ZWwwbUl0tAx0OEopddEKSgTvuZ4TgQTXc6LH+wrrJus31CaVl+xDAx2KUkpdtHwTgTFmoOu5mTGmues559G89EIseypLBhNs/2ONI4r1jnI/Ro9Sqpwr6GRxx4IepRlkWXSL9WvqcpKX7EMCHYpSSl2Ugk4W/6uAeQa4qoRjCSrhksndtqU8aR/ND9nt+JP1R+eM/MpPaOkJpVQZVdAIZb1LM5BgdJN1BW/Zr+cl+1CusDyF6K0FSqkglG8iEJGrjDErRORGb/ONMR/7L6zgECZZTLQt4XH7WNY4Iulu3R7okJRSqsgKumoop47C9V4eA/0cV9AYbl1FQ47zon2YDlyjlApKBXUNTXU931Z64QSfSmJnom0Jj9jHscoRTW/rlkCHpJRSReJLGeraIjJTRDaJSKKIvCIitUsjuGAx1PotjeUoL9uH6lGBUiro+FJ0biFwDBiCsyT1MeB9fwYVbEIlm0nWxWwxLVjhiA10OEopVSS+JIJaxpinjTE/ux7PADX8HFfQucG6hsvkd17UowKlVJDxJRGsFJGRImJxPYYDn/k7sGATItlMsi1mh2nGckdcoMNRSimfFXRn8RkROQ3cibPuUKbrsRC4v3TCCy6DLd/TXA7zkn2ojleglAoaBdUaqmqMqeZ6thhjbK6HxRhTrTSDDBY2cfCg7UN2mUtZmK334ymlgkNBJSbcRKQm0BIIy5mWd/hK5TTAsp54+ZEZ9hEMtK6nupwLdEhKKVUgXy4fHQd8BywDnnQ9T/NvWMFLBKaFzCeVylqQTikVFHw5WXwv0Bn4xVV/KBbnJaSFEpF+IrJbRPaJyBQv82/2GALzBxGJLlL0ZVRby0Futn7Nf7KvZZejSaDDUUqpAvmSCNKNMekAIlLJGLMLaF3YQiJiBV7HOdB9O+AmEWmXp9nPQE9jTAfgaWBWUYIvyx6wLaIqaTxpv1UvJ1VKlWm+JIJkEakBLAG+EpFPgMM+LNcF2GeM+ckYk3O10WDPBsaYH4wxJ11v1wGNfQ28rKspZ3nQ9iFrHe350tE50OEopVS+Ck0ExpgbjDGnjDHTgMeBOcCffVh3I+Cgx/tk17T83A584W2GiIwXkQQRSTiWFjy71zdZV9BGfuGZrFs4n5kd6HCUUsorX44IckYrmwR0AJJde/iFLuZlmtdfcRHpjTMRTPY23xgzyxgTZ4yJqxsRPNfn28TBtJD5HKIu//5uf6DDUUopr3y5augJYD5QG6gDzBWRx3xYdzLgeaa0MV66lESkAzAbGGyMOeFL0MGkq2UX11nW8uaq/SSfTAt0OEopdQExhZzJFJGdQKzHCeNwYJMxpm0hy9mAPcDVwCFgI/AXY8wOjzaXAiuAW40xP/gScFxDq0kYX8WXpmXGIVObqzNmcLVlM6+HzrywgQ5jqZTyMxFJNMZ4rX/jS9fQATxuJAMqAYX2cxhj7MBEnPcd7AQ+MMbsEJEJIjLB1ewJnEcab4hIkogk+BBP0GkkJ7jbtpTPHF35ITvvhVNKKRVY+R4RiMirOPv0L8V5H8FXrvfXAmuMMSNLK0hPwXhEAJBuQrgmczphZPJp6KOESdYfM/WIQCnlZwUdERRUYiJn7zwRWOwxfVUJxVWhhEkW/7C9w+isKbxsH8KUkIWBDkkppYCCh6qcn/NaREKBVq63u40xWd6XUgXpad3KTY4VzMoeSB9rAh0t+wIdklJK+XTVUC9gL867hN8A9ohID/+GVX49YnuXBpzgoawJpJuQQIejlFI+nSz+F9DHGNPTGNMD6Au85N+wyq+qcp4XQmbxk2nIDPvwQIejlFI+JYIQY8zunDfGmD2A7spehG7WHYyyLmdOdn82OAot26SUUn7lSyJIFJE5ItLL9Xgb5wlkdRGm2BbQWI7zcNadpGXaAx2OUqoC8yURTAB2AJNwlqT+0TVNXYTKksH0kLf4xVzCC1/uLnwBpZTykwJHKBMRC5BojIkEXiydkCqOrpZd3Gb9grk/9KdP+/r8qUWdQIeklKqACjwiMMY4gC2uUhDKD/5me5+mtSP426KtnM3QLiKlVOnzpWuoAbBDRL4RkaU5D38HVlGESyYzhkVz6NR5nly6o/AFlFKqhPkyeP2Tfo+igotrWouJvS/n1RX7iGtakxGd9QBMKVV68k0EIhKG86Tw5cA2YI6rkJzyg/uuaUXSwVM8/skO2jesTmSj6oEOSSlVQRTUNTQfiMOZBPrjvLFM+YnVIrwyMpY6lUOZ8N9ETqX5MvaPUkpdvIKqj24zxkS5XtuADcaYjqUZnDfBWn3UV5sdLRieOZXulq3MfuYxLJbgGZFNKVV2Fbf6qLuwnDHGLqI/SKUh1rKfJ2z/x+P2sbz+xBjusS3x3lBLVyulSkhBiSBaRE67XgsQ7novgDHGVPN7dBXULdav2eRoyYv2ocTIPrpbtwc6JKVUOZbvOQJjjNUYU831qGqMsXm81iTgRyLwj5B3aCXJTMq6h0OmdqBDUkqVY77cR6ACIEIyeCvkZexYuTvzXjKML1f6KqVU0WkiKMOaWX5neshbbDGX8/esO8jnvL5SSl0UTQRlXD9rAg/aPuBjR3em20cEOhylVDmk/Q1BYKJ1Cb+ZWryRPZj6cpLRtuWBDkkpVY5oIggCIvC0bS7HTA2m2W+lvpykX6CDUkqVG9o1FCSsYpgZ8hqxso9JWX9lw88pgQ5JKVVOaCIIIuGSyZzQGTSW44ybv5G9R84EOiSlVDmgiSDI1JSzzA95nkohVka/s4HfUs8HOiSlVJDTRBCEmliOM++2zpxOtzPmnY2kpmUVvpBSSuVDE0GQat+wOv8e1Ymfj5/jprfXkXJOq5UqpYpHE0EQ63Z5HWbd2on9x87yl7fXcfxsRqBDUkoFIU0EQa5X63q8M6YzB06cY+SsdRw9nR7okJRSQUYTQTnQ7fI6zL+tC7+dOs+IWev0BLJSqkg0EZQT8c1r83+3x3P8TAbD/72WgylpgQ5JKRUkNBGUI50uq8l/x8WTmpbFyFnr+OXEuUCHpJQKAn5NBCLST0R2i8g+EZniZX4bEVkrIhki8pA/Y6koopvUYMH4rqRl2hn+77XsOKwjmSmlCua3RCAiVuB1nAPftwNuEpF2eZqlAJOAGf6KoyJqP+tSFtrvx3L6MMNmfs1Xj/eGadUDHZZSqozy5xFBF2CfMeYnY0wmsBAY7NnAGHPUGLMRj/GRVclobUnmk0qP01IOMT7rAd62D8DogAZKKS/8mQgaAQc93ie7phWZiIwXkQQRSTiWpj9mvqonp1gY+jT9LRv5h/0WHlm8jaxsR6DDUkqVMf4sQy1ephXrV9wYMwuYBRDX0KqZAHzu6gmXTF4LmclL9qG8uuEGfkn4gjdCXqGG5DmRPE3PJShVUfkzESQDTTzeNwYO+/HzVD4sYngw5EOaWX5jStYd3Jj5JHNCZtDM8vsfjfJLLJoglCr3/Nk1tBFoKSLNRCQUGAks9ePnqULcaF3Df0Of5aSpyuDMp/kyu3OgQ1JKlQF+SwTGGDswEVgG7AQ+MMbsEJEJIjIBQEQuEZFk4AHgMRFJFpFq/opJQRfLbpaGPkYz+Z0JWffzWNZtpJuQQIellAogCbYrSeIaWk3C+CqBDiPoZRorM+wjmJU9kDbyK6+FzORyi5eeO+0aUqpcEJFEY0yct3l6Z3EFFSrZPBLyHnND/slRU4PrM5/hA3tPgmy/QClVAjQRVHC9rVv4otIUYi37+Jv9Tu7N+itnTHigw1JKlSJNBIr6cor/hDzHQ7b3+czRlX4Zz7MyOzrQYSmlSokmAgWAVQwTbZ/wQehThEkmt2VN5t7Mv+pgN0pVAJoIVC6dLHv5PPTv3Gv9iM8d8Vzz4rcsSkzW8hRKlWOaCNQFKomd+0M+4rPQv9OibhUe+nALo+Zs0LLWSpVTmghUvlpZDvHhnVfw9OD2JB08Rd+Xv+PVb/ZyPjM70KEppUqQJgJVIItFGHVFU756oAc9W9XlX1/tofeMVSxKTCbbod1FSpUHekOZKpL1jjY8m3UzW0wL2jWoxiMD2nJlyzqBDkspVQi9oUyVmHjLLhaHPsErIa+Sej6LW+asZ8zcDez+/UygQ1NKFZMmAlVkFjEMtq7lmwd78siANiT+cpL+r3zHfQs3s/eIJgSlgo0/y1Crci4sxMr4Hi0Y1qkJb6zax3/X/cqSpMP0a38JE6+6nMhGOjymUsFAzxGo4stTkC7lXCbvrPmZ+T8c4EyGnV6t6zKx9+XENa0VoACVUjkKOkegiUAVXz6VSU+nZ/Gftb8wZ83PpJzLJL5ZLcZe2Yyr29TDZtXeSKUCQROB8o9CSlSnZdpZsOEgs1f/xG+p6TSsHsbNXS9jROcm1KlSqZSCVEqBJgLlLz6OVWDPdvD1zqP8Z90Bvt93glCrhQFRlzDqiqZ0vLQGIt6Gt1ZKlaSCEoGeLFZ+Z7Na6Bd5Cf0iL2Hf0bP8d90vfJSYzJKkw7RrUI0hnRozKLohdavqUYJSgaBHBKr4LmL0snMZdhZvPsTCjb+y/dBprBahO5u5wbqaPpZEwiWzRD5HKeWkRwSqzKlcycYtXS/jlq6XsffIGRZvPsSSVce413EPVUijn3Ujf7Z8T7xlJzqislL+pUcEqvhKeE/dMbUG6x1tWOy4ki+y4zlDBNU5y9WxrenT/hJ6tKpDRKjuuyhVHHpEoIKCRQxXWHdyhXUnT9nm8a0jmmXZcXyzqyYfbz5EWIiF7i3r0qddfa5uW59alUMDHbJS5YIeEaiSV9wjhWne70TOevwkG39OYdmO31n+4xF+S01HBKIaVad7yzp0b1mXjpfWJNSm9ygolR+9fFSVrhJOBJ7rM8aw7VAqq3YfY/XeY2z69RTZDkNEqJUrmteme8s6XNGiDi3rVcFi0ctSlcqhXUOq3BAROjSuQYfGNZh0dUtOp2exbv8JVu89zuq9x/hm11EAakSEEHdZTTo3rUXnZrWIbFhdjxiUyocmAlXy8tuzhxI/wVzt+Tr0Afq43h8Mrcs6R1sS2v2DjQdS+HqnMzGEhViIaVKD2Etr0qFRdTo0qUHD6mF6M5tSaCJQ5UwTyzGaWI4xbGgHAI6dySDhQAobDqSQcOAkb3/3E3bXyGp1qoS6ji6q06Fxddo1qE79apU0OagKRxOBKtfqVq1E/6gG9I9qAEB6VjY7fzvN1uRUtiSfYmtyKit3HyXnVFmNiBDaXFKVNpdUo20D53Or+lUJD7UG8Fso5V+aCFT5lE/3VBgQOy2V2EtruqedSc/ix8On2fX7GXb9fpqdv53h/Y0HOZ+V7W7TqEY4LepVoUXdyrSoW8X5qFeZulX0CEIFP00Equwr6JxDCagaFkJ889rEN6/tnuZwGH5NSWPX72fY/fsZfloxl/2pDdm4pwHnCXO3q0Ialza4hMtqR3BprQgudT1fVqsyDWuEadltFRQ0ESjlhcUiNK1TmaZ1KtMv8hJY8zoADiP8Tk32Oxqy3zTkZ9OAX6vdwu4jZ/hm51Eysx1/rEOgfrUwGtYIdz3CaOx6Xb9aGPWrhVG7cqhe5qoCThOBUkVgEUNDUmhoTaE7250Tb5sJOI8ifj+dzi8n0vg15RyHTp7n0Kl0Dp86z9bkUyzbnp4rUQDYLEK9qpWoVy2M+tUqUX/Xf6gjqdQhlTqSSm05TV3X64gnj5b211UVhN5QplSOgi5tLU73VJ71ORyG4+cyOHTyPEfeHsYRU9P1qMFRarrfp+L97zs8xEqtyqHUrBxCzYhQakaEOt9HOKdVDw+hWrjzuXp4CDVc70O0e0qhN5QpVSZYLEK9qmHUqxoG1o35tss0Vk5QneMm51GN41TnePwjnEzL5OS5TE6mZfFrShop5zI5k24v8HMjQq1UDbNRNSwk13O1MBtVKtmoXMn5HBFqo3Ilq3ta5VAb4aFWItwPm96UV075NRGISD/gFcAKzDbGPJ9nvrjmDwDSgDHGmE3+jEmpUlPMk9yhkk0DUmggKblnDPyv1/ZZ2Q5OpWWRet75OH3+j9c5j7Ppds5kZHEm3U7q+SyST6ZxJt3O2XR7rqujCmOzCOGhVsJDrO7nsBArYSEW97RKNuf7SjYrlWwWKoW4nnNeWy1UCrEQ6n62EmqzEGqzEGIVQq0WQqw5753tbFYhxOqcr1dplTy/JQIRsQKvA9cCycBGEVlqjPnRo1l/oKXrEQ+86XpWqvT5+eqki5ZPfCFAXdejaOtzdl1lOwznMu2kZWRzNsPOude6c86EcY4w0qjEeVOJNNfrtCsfIS0zm/SsbM5n5Tw7SN/7HScI5TyVyCCEdBPqfCaUDEq2SqwNOyHYsZFNCNmuZzu22pdhswg2izNx2KwWbBbBapE8z67pVsEqf8zzfFgkzzQRLDnP4jy6c7YDi/zR3vL5g1hwYMFgxYFgsIrzvQx5293WIs5yKTnLSwHPQu724noteEzD2RY8lvOYXtj1CP48IugC7DPG/IQz4IXAYMAzEQwG/s84T1SsE5EaItLAGPObH+NSSnmwWoRqYSFUC3MNAWT5Kf/G/f7P+/Rp1+a7iJl6igy7gwy7g0y7gwx7Npl2B5mvdiWDEDIJIdPYyMRGJiFkYXU+u6fZsGMjCytZxkYWNuxYXdOt2LGSZazYm0Rjdxjs2Q7s2Qa7w5DtMNgdDrKyHZzPcr7PyjY4XNOzHYZsY8j2aJ9tnM8Oj9fZDoPDp9OpY/OftTDJlxUEhD8TQSPgoMf7ZC7c2/fWphGQKxGIyHhgvOtthjx5envJhlrq6gDHAx3ERQj2+CH4v8PFx/9kMbpYirVMvucVSvjfYHrJrco3wfY3dFl+M/yZCLz9xeTNqb60wRgzC5gFICIJ+Z35DhbB/h2CPX4I/u8Q7PFD8H+HYI/fkz8vAUgGmni8bwwcLkYbpZRSfuTPRLARaCkizUQkFBgJLM3TZilwqzh1BVL1/IBSSpUuv3UNGWPsIjIRWIbz8tF3jDE7RGSCa/5bwOc4Lx3dh/Py0dt8WPUsP4VcmoL9OwR7/BD83yHY44fg/w7BHr9b0N1ZrJRSqmTpbYJKKVXBaSJQSqkKLqgSgYj0E5HdIrJPRKYEOp6iEJEmIrJSRHaKyA4RuTfQMRWHiFhFZLOIfBroWIrDddPiIhHZ5fq3uCLQMRWViNzv+hvaLiILRCSs8KUCR0TeEZGjIrLdY1otEflKRPa6nmsWtI5Ay+c7THf9HW0VkcUiUiOAIV6UoEkEHiUr+gPtgJtEpF1goyoSO/CgMaYt0BX4a5DFn+NeYGegg7gIrwBfGmPaANEE2XcRkUbAJCDOGBOJ80KMkYGNqlDzgH55pk0BvjHGtAS+cb0vy+Zx4Xf4Cog0xnQA9gB/L+2gSkrQJAI8SlYYYzKBnJIVQcEY81tOQT1jzBmcP0CNAhtV0YhIY+A6YHagYykOEakG9ADmABhjMo0xpwIaVPHYgHARsQERlPF7b4wx3wF5KugxGJjvej0f+HNpxlRU3r6DMWa5MSan9Os6nPdBBaVgSgT5laMIOiLSFIgF1gc4lKJ6Gfgb4CikXVnVHDgGzHV1b80WkcqBDqoojDGHgBnArzhLsaQaY5YHNqpiqZ9zz5DruV6A47lYY4EvAh1EcQVTIvCpHEVZJyJVgI+A+4wxpwMdj69EZCBw1BiTGOhYLoIN6Ai8aYyJBc5R9rskcnH1pQ8GmgENgcoicktgo6rYRORRnF2/7wY6luIKpkQQ9OUoRCQEZxJ41xjzcaDjKaJuwCAROYCzW+4qEfFeIL/sSgaSjTE5R2KLcCaGYHIN8LMx5pgxJgv4GPhTgGMqjiMi0gDA9RyU43CKyGhgIHCzCeKbsoIpEfhSsqLMcg3CMwfYaYx5MdDxFJUx5u/GmMbGmKY4t/0KY0xQ7YkaY34HDopIa9ekq8ldFj0Y/Ap0FZEI19/U1QTZCW+XpcBo1+vRwCcBjKVYXANvTQYGGWPSAh3PxQiaROA6KZNTsmIn8IExZkdgoyqSbsAonHvSSa7HgEAHVQHdA7wrIluBGODZwIZTNK6jmUXAJmAbzv/DZbrUgYgsANYCrUUkWURuB54HrhWRvTgHr3q+oHUEWj7f4TWgKvCV6//zWwEN8iJoiQmllKrgguaIQCmllH9oIlBKqQpOE4FSSlVwmgiUUqqC00SglFIVnCYCFRREZJWI9M0z7T4ReaOQZfw6uLir+udWEbnfh7YxvlwynLediAwqqNpuUdsrlZcmAhUsFnBhlc2RrukBISKXAH8yxnQwxrzkwyIxOIdmLVI7Y8xSY0xB19kXtb1SuWgiUMFiETBQRCqBu3BfQ2CNiLwpIgmuGv1PeltYRM56vB4qIvNcr+uKyEcistH16OZl2TARmSsi21zF6nq7Zi0H6rluJuqeZ5lhrvECtojId6674Z8CRrjajxCRLiLyg2udP4hI63zajRGR14qwXs/29V218re4HsFYjkL5md8Gr1eqJBljTojIBpw14T/BeTTwvjHGiMijxpgU15gV34hIB2PMVh9X/QrwkjFmjYhcivPO9bZ52vzVFUOUiLQBlotIK2AQ8KkxJsbLep8A+hpjDolIDWNMpog8gXMcgYnwR1lsY4xdRK4BnjXGDPHSbkwR1+vZfibwrTHmBtf2qeLjdlEViCYCFUxyuodyEsFY1/ThIjIe599zA5wDF/maCK4B2jnL9gBQTUSqusaMyHEl8CqAMWaXiPwCtAIKqh77PTBPRD7AWRjOm+rAfBFpibOSbogP8fqyXk9XAbe6Ys8GUn1YRlUwmghUMFkCvCgiHYFwY8wmEWkGPAR0NsacdHX5eBu60bOWiud8C3CFMeZ8AZ/rrQR6gYwxE0QkHudAPkkiEuOl2dPAStfeelNgVQmtV6ki0XMEKmgYY87i/LF8hz9OElfDOa5AqojUxzmUqTdHRKStiFiAGzymL8dZzBBwXoHjZdnvgJtd81sBlwK7C4pVRFoYY9YbY54AjuMsoX4GZ5GyHNWBQ67XYzym521X1PV6+ga4y7Ws1dUdpVQumghUsFmAc6zhhQDGmC3AZmAHzgTxfT7LTQE+BVbgHNkrxyQgznUJ6I/ABC/LvgFYRWQb8D4wxhiTUUic010nl7fjTCRbgJU4u6GSRGQE8ALwnIh8j3Ps4Rx52xV1vZ7uBXq7Yk8E2hcSt6qAtPqoUkpVcHpEoJRSFZwmAqWUquA0ESilVAWniUAppSo4TQRKKVXBaSJQSqkKThOBUkpVcP8PUob1Gf3FeJgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "n_observations = 500  # number of observations\n",
    "n_samples = 1000  # number of samples\n",
    "# standard normal distribution can be used, since the\n",
    "# statistic is unaffected by location and scale\n",
    "norm_dist = stats.norm()\n",
    "# Draw 1000 samples, each with 500 observations\n",
    "y = norm_dist.rvs(size=(n_observations, n_samples))\n",
    "\n",
    "# calculate the value of the statistic for each sample\n",
    "# we'll call this the \"Monte Carlo null distribution\"\n",
    "null_dist_mc = statistic(y)\n",
    "\n",
    "# the asymptotic null distribution is chi-squared with df=2\n",
    "null_dist = stats.chi2(df=2)\n",
    "y_grid = np.linspace(0, null_dist.isf(0.001))\n",
    "pdf = null_dist.pdf(y_grid)\n",
    "\n",
    "# compare the two\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(y_grid, pdf)\n",
    "plt.hist(null_dist_mc, density=True, bins=100)\n",
    "plt.xlim(0, np.max(y_grid))\n",
    "plt.xlabel(\"Value of statistic\")\n",
    "plt.ylabel(\"Probability Density\")\n",
    "plt.legend(['Asymptotic Null Distribution', 'Monte Carlo Null Distribution (500 observations/sample)'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87d87faf-f1e1-462f-a94c-cf09d9445767",
   "metadata": {},
   "source": [
    "As we can see, the *Monte Carlo null distribution* <a name=\"cite_ref-2\"></a>[<sup>**†**</sup>](#cite_note-1) of the test statistic when samples are drawn according to the null hypothesis (from a normal distribution) appears to follow the *asymptotic null distribution* (chi-squared with two degrees of freedom). \n",
    "\n",
    "<a name=\"cite_ref-1\"></a>[<sup>**†**</sup>](#cite_note-1) Named after the Monte Carlo Casino in Monaco, apparently [[3]](https://en.wikipedia.org/wiki/Monte_Carlo_method#Historyhttps://en.wikipedia.org/wiki/Monte_Carlo_method#History).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4195c3fd-12fe-46b4-b681-83e0f57550ae",
   "metadata": {},
   "source": [
    "Note that the originally observed value of the statistic, 6.98, is located in the right tail of the null distribution. Random samples from a normal distribution usually produce statistic values less than 6.98, and only occasionally produce higher values. Therefore, there is rather low probability of observing such an extreme value of the statistic under the null hypothesis that the sample is drawn from a normal population. This probability is quantified by the *survival function* of the null distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0f2fad67-7dc3-4bf4-af5f-12bad97321e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Under the null hypothesis, the chance of drawing a sample that produces a statistic value greater than \n",
      "6.982848237344646\n",
      "is\n",
      "0.03045746622458189\n"
     ]
    }
   ],
   "source": [
    "pvalue = null_dist.sf(stat1)\n",
    "message = (\"Under the null hypothesis, the chance of drawing a sample \"\n",
    "           f\"that produces a statistic value greater than \\n{stat1}\\n\"\n",
    "           f\"is\\n{pvalue}\")\n",
    "print(message)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54c6b2aa-98ce-4463-be49-78732a50083d",
   "metadata": {},
   "source": [
    "This is the `pvalue` returned by `stats.jarque_bera`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bead807c-e722-4269-bfbe-3be8fccb2da8",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(pvalue, stats.jarque_bera(x).pvalue, rtol=1e-14)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78720622-ff97-4248-85a1-6e13b786ea06",
   "metadata": {},
   "source": [
    "When the $p$-value is small, we take this as evidence against the null hypothesis, since samples drawn under the null hypothesis have a low probability of producing such an extreme value of the statistic. For better or for worse, a common \"confidence level\" used for statistical tests is 0.99, meaning that the threshold for rejection of the null hypothesis is $p \\leq 0.01$. If we adopt this criterion, then the Jarque-Bera test was inconclusive; it gives insufficient evidence to conclude the null hypothesis is false. Although this should *not* be taken as evidence that the null hypothesis is *true*, the lack of evidence against the hypothesis of normality is often considered sufficient to proceed with tests that assume the data is drawn from a normal population."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8feb92d9-99b7-4863-885c-8876a64905f6",
   "metadata": {},
   "source": [
    "There are a few shortcomings with the test procedure outlined above. The Monte Carlo null distribution agreed with the asymptotic null distribution when the number of observations per sample was 500, but our original sample `x` had only 11 observations. If we generate a Monte Carlo null distribution of the statistic for sample size of only 11 observations, there is marked disagreement between the Monte Carlo null distribution and the asymptotic null distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1fe77f87-b4cc-4b0f-baee-cf79ca2088ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABAX0lEQVR4nO3dd3wUdf748dc7m0AIJbSAAlJVENJAmqJURZooFiyIoiKH/Syc3Hkqlp/nCWc79TxEgdMDCwrytStFwEOkSAeVpjQFQhdC2vv3x0zWTdhsNmWz2eT9fDzmsTuzU96z2ex7Zz4z74+oKsYYYyqvqHAHYIwxJrwsERhjTCVnicAYYyo5SwTGGFPJWSIwxphKLjrcARRV/fr1tXnz5uEOwxhjIsry5cv3qWqCv9ciLhE0b96cZcuWhTsMY4yJKCLyU0Gv2akhY4yp5CwRGGNMJWeJwBhjKrmIayMoTZmZmezYsYP09PRwh2KMMaUiNjaWJk2aEBMTE/QylToR7Nixg5o1a9K8eXNEJNzhGGNMiagqaWlp7NixgxYtWgS9XKU+NZSenk69evUsCRhjKgQRoV69ekU+y1FoIhCRZSJyu4jUKXZ05ZglAWNMRVKc77RgjgiuBhoBS0XkLRG5SOzb0xhjKoxCE4GqblLVB4EzgWnA68DPIvKoiNQNdYCVwcyZMxERNm7cGJbtb9u2jWnTphU638qVK/n444+947Nnz+app54KejvNmzfn8ssv947PmDGDESNGFLpcjRo1vHEmJiae9Pq2bduoVq0a7du356yzzqJz585MnTo16Djz71d+y5Yt46677gJg3LhxTJgwodCYfT333HMcO3bMOz5gwAAOHjxYpHUYE0pBtRGISDLwD2A88B5wBXAYmBu60EJgXPzvQzkyffp0zjvvPN56662wbL+4iWDw4MGMHTu2SNtatmwZ69atK3KMhWnVqhXfffcdGzZs4K233uLZZ59l8uTJQcUZKBFkZWXRsWNHXnjhhWLHlj8RfPzxx9SuXbvY6zOmtAXTRrAceBZYCiSr6l2qukRV/wFsCXWAFd3Ro0f5+uuvee211/Ikgt27d9O9e3dSU1NJTExk4cKFvPbaa9xzzz3eeV599VXuvfdetm3bRps2bRg5ciSJiYkMGzaML7/8km7dunHGGWfw7bffAs6v2eHDh9O7d2/OOOMMXn31VQDGjh3LwoULSU1N5dlnnyU9PZ0bb7yRpKQk2rdvz7x588jIyODhhx/m7bffJjU1lbfffpspU6Zwxx13APDrr78yZMgQUlJSSElJ4X//+5/f/b3//vt58sknT5qe/5d2YmIi27ZtK9Z72rJlS5555hnvl7dvnO+++y6JiYmkpKTQvXt3v/s1btw4Ro0aRd++fbn++uuZP38+gwYN8q5/1apVJ72H+ee54447mDJlCi+88AK7du2iV69e9OrVC3COjPbt2wfAM888Q2JiIomJiTz33HOAk5jPOussbrnlFtq1a0ffvn05fvx4sd4LY4IRzOWjV6qq3y98Vb2slOMJm0f/bx3rdx0u1XW2bVSLRy5uF3CeWbNm0a9fP84880zq1q3LihUr6NChA9OmTeOiiy7iwQcfJDs7m2PHjtGhQweSk5N5+umniYmJYfLkyfz73/8GYNOmTbz77rtMnDiRTp06MW3aNBYtWsTs2bN58sknmTVrFgCrV6/mm2++4bfffqN9+/YMHDiQp556igkTJvDhhx8C8I9//AOANWvWsHHjRvr27csPP/zAY489xrJly3jxxRcB5ws211133UWPHj2YOXMm2dnZHD161O/+Dh06lJdffplNmzaV5K0tVIcOHfyeanvsscf47LPPaNy4MQcPHqRKlSon7de4ceNYvnw5ixYtolq1asyfPz/POvy9hwW56667eOaZZ5g3bx7169fP89ry5cuZPHkyS5YsQVXp0qULPXr0oE6dOvz4449Mnz6dV199laFDh/Lee+9x3XXXlfyNMcaPYE4NjRSR2rkjIlJHRJ4IXUiVy/Tp07n66qsBuPrqq5k+fToAnTp1YvLkyYwbN441a9ZQs2ZNqlevTu/evfnwww/ZuHEjmZmZJCUlAdCiRQuSkpKIioqiXbt29OnTBxEhKSkpzy/rSy65hGrVqlG/fn169erlPVrwtWjRIoYPHw5AmzZtaNasGT/88EPA/Zg7dy633norAB6Ph/h4/6ffPB4PY8aM4W9/+1vR3qgiKqgv7m7dujFixAheffVVsrOzC1x+8ODBVKtWze9rwbyHwVi0aBFDhgyhevXq1KhRg8suu4yFCxcCzt8zNTUVgLPPPrvYR0fGBCOYI4L+qvqX3BFVPSAiA4C/hi6sslfYL/dQSEtLY+7cuaxduxYRITs7GxHh6aefpnv37ixYsICPPvqI4cOHM2bMGK6//npGjhzJk08+SZs2bbjxxhu966patar3eVRUlHc8KiqKrKws72v5L/jydwFYQV+ipWX48OH87W9/o12739/z6OhocnJyvOMlvdv7u+++46yzzjpp+iuvvMKSJUv46KOPSE1NZeXKlX6Xr169eoHr9vceFif+QO+z79/T4/HYqSETUsEcEXhExPupFJFqQNUA85sgzZgxg+uvv56ffvqJbdu2sX37dlq0aMGiRYv46aefaNCgAbfccgs333wzK1asAKBLly5s376dadOmcc011xR5mx988AHp6emkpaUxf/58OnXqRM2aNTly5Ih3nu7du/Pf//4XgB9++IGff/6Z1q1bnzSfrz59+vCvf/0LgOzsbA4fLvg0W0xMDPfcc4/3nDg4581z93HFihVs3bq1yPuWa9u2bdx///3ceeedJ722efNmunTpwmOPPUb9+vXZvn17wP3yx9972KxZM9avX8+JEyc4dOgQc+bM8c5f0Pq7d+/OrFmzOHbsGL/99hszZ87k/PPPL95OG1MCwSSCN4E5InKziNwEfAFMLWQZE4Tp06czZMiQPNMuv/xypk2bxvz580lNTaV9+/a899573H333d55hg4dSrdu3ahTp+j3+HXu3JmBAwfStWtXHnroIRo1akRycjLR0dGkpKTw7LPPctttt5GdnU1SUhJXXXUVU6ZMoWrVqvTq1Yv169d7G1V9Pf/888ybN4+kpCTOPvvsQq8Muvnmm/McqVx++eXs37+f1NRU/vWvf3HmmWcWab82b97svXx06NCh3HnnnXmOmHKNGTOGpKQkEhMT6d69OykpKQH3yx9/7+Fpp53G0KFDSU5OZtiwYbRv3947/6hRo+jfv7+3sThXhw4dGDFiBJ07d6ZLly6MHDkyz3LGlBUJ5jSAiPQH+gACfK6qn4U6sIJ07NhRi90xje9lo+MOsWHDBr+nD8q7QYMGcc8999CnT58iLTdu3Dhq1KjB/fffH6LIjDHlgb/vNhFZrqod/c0fVNE5Vf0E+KTk4ZmSOHjwIJ07dyYlJaXIScAYYwpSaCIQkcuAvwMNcI4IBFBVrRXi2Ew+tWvXLvTqnUDGjRtXesEYYyqMYI4IngYuVtUNoQ7GGGNM2QumsfhXSwLGGFNxBXNEsExE3gZmASdyJ6rq+6EKyhhjTNkJJhHUAo4BfX2mKWCJwBhjKoBgylDf6Ge4qSyCK3O+1UlLYwiCiHjLOYBT7TIhISFPAbOiOHjwIC+//HKRlzt69Ch/+MMfaNWqFe3ataN79+4sWbKkSOvo2bMnRbm0t2fPnnTs+PvVbMuWLaNnz56FLudbtC23RHV+Ho+H1NRU2rVrR0pKCs8884z3zl/fstL+FFaNddeuXVxxxRVA3oJ2wZoyZQq7du3yjo8cOZL169cXaR0FmTVrFo899hgACxYsoEOHDkRHRzNjxow88/Xr14/atWsX63NW1L9zqOQvXnjuueeW+jZ2795N3759C5+xlPl+xgtywQUXcODAgVLZXjDVR88UkTkistYdTxaRClVeIpyqV6/O2rVrvSUEvvjiCxo3blzs9RU3EYwcOZK6devy448/sm7dOqZMmVLoB9FXoLo9gezZs4dPPin9K5OrVavGypUrWbduHV988QUff/wxjz76KEChZaUDJYKsrCwaNWp00hdrUeRPBJMmTaJt27bFXp+vp59+mttuuw2Apk2bMmXKFK699tqT5hszZgxvvPFGqWyzpFQ1T3mOYOVPBAVVvC2JTz/9lIsuuqjU11sahg8fXqz/dX+CaSx+FfgzkAmgqqtxei2LHOWwDwJf/fv356OPPgKcu419S0fs37+fSy+9lOTkZLp27crq1asB51LQm266iZ49e9KyZUvvF9vYsWPZvHkzqampjBkzBoDx48fTqVMnkpOTeeSRR07a/ubNm1myZAlPPPEEUVHOR6Jly5beqpqXXnopZ599Nu3atWPixIne5WrUqMHDDz9Mly5dWLx4cZ51Tp8+3XsH7wMPPFDgvo8ZM4Ynnji5hmH+X9qDBg06qQposBo0aMDEiRN58cUXUdU8JaO/+uorUlNTvXdxHzly5KSy3FOmTOHKK6/k4osvpm/fvid1kLN9+3b69etH69atvckm/zwTJkxg3LhxzJgxg2XLljFs2DBSU1M5fvx4nl/YBb1vNWrU4MEHHyQlJYWuXbvy66+/nrSfP/zwA1WrVvVWOW3evDnJycnev6mvPn36ULNmzYDv28qVK+natSvJyckMGTIkz6/PN998k3PPPZfExERv0T1/7yX4//zlltq+7bbb6NChA48//jh/+tOfvOufMmWKt0SIv8/f2LFjOX78OKmpqQwbNsz7HoGTWMaMGUNiYiJJSUneu8Xnz59Pz549ueKKK2jTpg3Dhg3z1nsaO3Ysbdu2JTk5Oc8Nl59++in9+/f3WxYe4NZbb6Vjx460a9cuz/9W8+bN+ctf/sI555xDx44dWbFiBRdddBGtWrXilVde8cbTvXt3hgwZQtu2bRk9erTfhPjmm2/SuXNnUlNT+cMf/uD90TV48GBvkcoSU9WAA7DUffzOZ9rKwpYL1XD22WdrkT1S6+RBVdevX1/4fCUZglC9enVdtWqVXn755Xr8+HFNSUnRefPm6cCBA1VV9Y477tBx48apquqcOXM0JSXFCfWRR/Scc87R9PR03bt3r9atW1czMjJ069at2q5dO+/6P/vsM73llls0JydHs7OzdeDAgfrVV1/lieGDDz7QSy+9tMAY09LSVFX12LFj2q5dO923b5+qqgL69ttve+fr0aOHLl26VHfu3KmnnXaa7tmzRzMzM7VXr146c+bMk9abO3+vXr107ty5unTpUu3Ro4eqqk6ePFlvv/1277wDBw7UefPmqapqs2bNdO/evd73r6D3Nb/atWvrL7/8kuf9HTRokC5atEhVVY8cOaKZmZl5Xs+NpXHjxt73wfc9njx5sp5yyim6b98+7/uzdOnSk/4O48eP10ceeSTPfhflfQN09uzZqqo6ZswYffzxx0/av9dff13vvffek6bfcMMN+u677540Pf9+5peUlKTz589XVdWHHnpI7777bm+8I0eOVFXVr776yruf/t7Lgj5/W7duVRHRxYsXq6rqnj17tFWrVt5t9+vXTxcuXKiqBX/+8v+Nc8dnzJihF1xwgWZlZekvv/yip512mu7atUvnzZuntWrV0u3bt2t2drZ27dpVFy5cqGlpaXrmmWdqTk6OqqoeOHBAVVWzsrK8/28TJkzQJ554wjv98OHDeWLLysrSHj166KpVq1TV+Yy+/PLLqqr6xz/+UZOSkvTw4cO6Z88eTUhI8L7/VatW1c2bN2tWVpZecMEF3r9T7md8/fr1OmjQIM3IyFBV1VtvvVWnTp3q3efTTz/d+374Oum7TVWBZVrA92owRwT7RKQVTgMxInIFsLt00lCIlfMjgVzJycls27aN6dOnM2DAgDyv+ZaE7t27N2lpaRw6dAiAgQMHen8BNmjQwO+vxM8//5zPP/+c9u3be2v0//jjj0WK74UXXvD+Et2+fbt3eY/Hk6fryVxLly6lZ8+eJCQkEB0dzbBhw1iwYEGB6//rX//q96igtKmfcirdunXj3nvv5YUXXuDgwYNER/u/fuLCCy+kbl3/PbNeeOGF1KtXj2rVqnHZZZexaNGiYsUX6H2rUqWK9yimoLLUu3fvJiEhoVjbzu/QoUMcPHiQHj16AHDDDTfk+RvmHrV2796dw4cPc/DgQb/vZaDPX7NmzejatSsACQkJtGzZkm+++Ya0tDS+//57unXrBhT8+SvIokWLuOaaa/B4PDRs2JAePXqwdOlSwKkT1aRJE6KiokhNTWXbtm3UqlWL2NhYRo4cyfvvv09cXBwAS5YsoUuXLoD/svAA77zzDh06dKB9+/asW7cuT1vP4MGDAUhKSqJLly7UrFmThIQEYmNjvV2Vdu7cmZYtW+LxeLjmmmtO+uzMmTOH5cuX06lTJ1JTU5kzZw5btvzePUyDBg3ynGYsrmCuGrodmAi0EZGdwFbAesgoZYMHD+b+++9n/vz5pKWleaf7+/LKLYOcv1SxbxE33+X//Oc/84c//KHAbbdr145Vq1aRk5Nz0mmE+fPn8+WXX7J48WLi4uLo2bOnt8RybGwsHo/H7zaLonfv3jz00EN888033mmlXZZ6y5YteDweGjRowIYNv98WM3bsWAYOHMjHH39M165d+fLLL/0uH+6y1DExMd7tFPS3rlatmvdHQqj522d/72VBn79t27ad9J5eddVVvPPOO7Rp04YhQ4YgIgE/fwUJ9D76+5+Jjo7m22+/Zc6cObz11lu8+OKLzJ07l08++YR+/foB+C0Lf/755zNhwgSWLl1KnTp1GDFiRJ7YfEvB5y8Tn/v3K6wsvKpyww03FNh/R3p6eoH9ZhRFMFcNbVHVC4AEoI2qnqeq20q8ZZPHTTfdxMMPP+ztaCaXb0no+fPnU79+fWrVKri6R/6SxxdddBGvv/66t8ewnTt3smfPnjzLtGrVio4dO/LII494/4l+/PFHPvjgAw4dOkSdOnWIi4tj48aNeb6sC9KlSxe++uor9u3bR3Z2NtOnT/f+sizIgw8+yNNPP+0db968OStXriQnJ4ft27cXu/MXgL179zJ69GjuuOOOk/7RNm/eTFJSEg888AAdO3Zk48aNRS5L/cUXX7B//36OHz/OrFmz6NatGw0bNmTPnj2kpaVx4sQJb+9vUHBZ6uK8b77OOuusUuv5LT4+njp16njPhb/xxht5Ysk9775o0SLi4+OJj4/3+14G8/nLddlllzFr1iymT5/OVVddBRDw8xcTE0NmZuZJ6+nevTtvv/022dnZ7N27lwULFtC5c+cC9/Xo0aMcOnSIAQMG8Nxzz3n7qJgzZ463ppe/svCHDx+mevXqxMfH8+uvvxbroodvv/2WrVu3kpOTw9tvv815552X5/U+ffowY8YM73u2f/9+fvrpJ8BJEr/88gvNmzcv8nbzC6bW0MP5xnODeKzEWy9vxpXNryl/mjRpkqfUdK5x48Zx4403kpycTFxcHFOnBq4AXq9ePbp160ZiYiL9+/dn/PjxbNiwgXPOOQdwGtTefPNNGjRokGe5SZMmcd9993H66acTFxdHvXr1GD9+PMnJybzyyiskJyfTunVr76F8IKeeeip/+9vf6NWrF6rKgAEDuOSSSwIuM2DAgDynNbp16+btdS0xMZEOHToUul1fuQ2JmZmZREdHM3z4cO69996T5nvuueeYN28eHo+Htm3b0r9/f6KiorxluUeMGFFoue/zzjuP4cOHs2nTJq699lrvJbG5DektWrSgTZs23vlHjBjB6NGjqVatWp5G9uK8b766d+/Offfdh6oiIixdutTbyPt///d/PPLII97y4Oeffz4bN27k6NGjNGnShNdee+2kq2OmTp3K6NGjOXbsGC1btmTy5Mne1+rUqcO5557L4cOHef311wt8L6tWrer38+fvSLJOnTq0bduW9evXe7+4+/XrV+Dnb9SoUSQnJ9OhQwfvjyWAIUOGsHjxYlJSUrwdPZ1yyil+uy4FOHLkCJdccgnp6emoKs8++yx79+4lNjbW+6Nr/vz5jB8/npiYGGrUqMF//vMfWrRoQfv27WnXrh0tW7b0nsoqinPOOYexY8eyZs0ab8Oxr7Zt2/LEE0/Qt29fcnJyiImJ4aWXXqJZs2YsX76crl27Fng6sygKLUMtIvf5jMYCg4ANGsS9BCLSD3ge8ACTVPWpfK/H4/R30BQnKU1Q1cknrchHkcpQB2ofiOAy1MYU5O677+biiy/mggsuCHcoEe3NN99kx44djB07NmTbmD9/fp6+wovq7rvvZvDgwX4rEZd6GWpV/Ue+lU0AZhe2nIh4gJeAC4EdwFIRma2qvnfO3A6sV9WLRSQB+F5E/quqGYWt3xhzsr/85S9FvhHQnOy668p/M2hiYmKplaMP5qqh/OKAlkHM1xnY5LYxZABvAfmPcxWoKc75phrAfuDkVjBjTFAaNmzovVrFlG89e/Ys9tEAwC233FJqsQTTRrAG99JRnFM8CUAw7QONge0+4zuALvnmeRHn6GIXUBO4SlVPuqNCREYBo8C5W7I05Z5PNcaYiqCoV+1BcJeP+hYjycIpSx3Mr3Z/3675I7wIWAn0BloBX4jIQlXN0/O5qk7EuYSVjh07Fn0vCxAbG0taWhr16tWzZGCMiXiqSlpaGrGxsUVaLphEkP86t1q+X5qqur+A5XYAp/mMN8H55e/rRuAp9663TSKyFWgDFP9aQQj6JrImTZqwY8cO9u7dW6LNGWNMeREbG0uTJk2KtEwwiWAFzhf6AZxf+bWBn93XlILbC5YCZ4hIC2AnTn2i/NWvfgb6AAtFpCHQGthCGYmJiaFFixZltTljjCmXgmks/hSnq8r6qloP51TR+6raQlULbDR2Tx/dAXwGbADeUdV1IjJaREa7sz0OnOu2Q8wBHlDV4EtelkSElJ8wxphQC+aIoJOq5n5xo6qfiMjjwaxcVT8GPs437RWf57vI2+GNMcaYMhZMItjn9j/wJs6poOuAtMCLGGOMiRTBnBq6BueS0ZnukOBOM8YYUwEEc2fxfuBuEamhqkfLICZjjDFlKJiuKs8VkfXAenc8RURKp380Y4wxYRfMqaFncW78SgNQ1VVA91AGVabs6iFjTCUXVK0hVd2eb1Lxeio3xhhT7gRz1dB2ETkXUBGpAtyFc1+AMcaYCiCYI4LROOWiG+OUjUh1x40xxlQAAY8I3D4FnlPVYWUUjzHGmDIW8IhAVbOBBPeUkDHGmAoomDaCbcDXIjIb+C13oqo+E6qgjDHGlJ1gEsEud4jC6TzGGGNMBVJgIhCRKao6QlUfFZEbVHVqWQZmjDGmbARqI0jxeX53qAMxxhgTHoESQal1CWmMMab8CtRG0EREXsDplSz3uZeq3hXSyIwxxpSJQIlgjM/zZaEOxBhjTHgUmAiscdgYYyqHoIrOGWOMqbgsERhjTCVnicAYYyq5QDeU/ZMAl5DaVUPGGFMxBDoiWAYsB2KBDsCP7pCKdUxjjDEVRqFXDYnICKCXqma6468An5dJdMYYY0IumDaCRuQtNlfDnWaMMaYCCKb66FPAdyIyzx3vAYwLWUQlYZ3QG2NMkRWaCFR1soh8AnRxJ41V1V9CG5YxxpiyUuipIRER4AIgRVU/AKqISOeQR1bWxsXbEYUxplIKpo3gZeAc4Bp3/AjwUsgiMsYYU6aCaSPooqodROQ7AFU9YH0YG2NMxRHMEUGmiHhwby4TkQQgJ6RRGWOMKTPBJIIXgJlAAxH5f8Ai4MmQRmWMMabMBHPV0H9FZDnQB6eTmktVdUPIIzPGGFMmCk0EIvI88Laqlt8GYrvaxxhjii2YU0MrgL+KyCYRGS8iHUMdlDHGmLJTaCJQ1amqOgDoDPwA/F1Efgx5ZMYYY8pEUfojOB1oAzQHNgazgIj0E5Hv3aOJsQXM01NEVorIOhH5qgjxGGOMKQXBtBH8HbgM2Ay8DTyuqgeDWM6Dc+PZhcAOYKmIzFbV9T7z1Ma5Ya2fqv4sIg2KsxPGGGOKL2AicMtLHAXOUdV9RVx3Z2CTqm5x1/UWcAmw3meea4H3VfVnAFXdU8RtGGOMKaGAp4ZUVXEuFy1qEgBoDGz3Gd/hTvN1JlBHROaLyHIRub4Y2zHGGFMCwZSY+EZEOqnq0iKuW/xMy9/1ZTRwNs49CtWAxSLyjar+kGdFIqOAUQBNmzYtYhjGGGMCCaaxuBdOMtgsIqtFZI2IrA5iuR3AaT7jTYBdfub5VFV/c486FgAp+VekqhNVtaOqdkxISAhi08YYY4IVzBFB/2Kueylwhoi0AHYCV+O0Cfj6AHhRRKKBKjh9HjxbzO0ZY4wphmDuI/gJ55d9b/f5sSCXywLuAD4DNgDvqOo6ERktIqPdeTYAnwKrgW+BSaq6trg7Y4wxpuiCuXz0EaAj0BqYDMQAbwLdCltWVT8GPs437ZV84+OB8cGHbIwxpjQF00YwBBgM/AagqrvI25m9McaYCBZMIshwLyPN7Y+gemhDMsYYU5aCSQTviMi/gdoicgvwJfBqaMMyxhhTVoLpj2CCiFwIHMZpJ3hYVb8IeWThklvSetyh8MZhjDFlJJjG4urAXFX9QkRaA61FJEZVM0MfnjHGmFAL5tTQAqCqiDTGOS10IzAllEEZY4wpO8EkAlHVYzgVSP+pqkOAtqENyxhjTFkJKhGIyDnAMOAjd1owdyQbY4yJAMEkgj8CfwZmuncGtwTmhTQqY4wxZSaYq4a+Ar4SkVoiUtPtX+Cu0IdmjDGmLBR6RCAiHUVkDU49oLUiskpEzg59aMYYY8pCMOf6XwduU9WFACJyHk7NoeRQBmaMMaZsBNNGcCQ3CQCo6iLgSOhCMsYYU5YKPCIQkQ7u02/dEhPTceoNXQXMD31ohci9A9gYY0yJBDo19I9844/4PM/f5aQxxpgIVWAiUNVeZRmIMcaY8AjYRiAiiSIyVUSWichS93lSWQXnz54jJ8K5eWOMqXAKTAQicgkwE/gKuAkY6T5/330tLH49nM6WvUfDtXljjKlwArURPAZcqKrbfKatEpG5OJ3OfxDKwAoiwKRFW3kyHBs3xpgKKNCpoZh8SQAAd1pMqAIqTJ3qVZixfAf7tFa4QjDGmAolUCLIFJGm+SeKSDMgK3QhBVa/RlUys3P4T1bfcIVgjDEVSqBE8AjwpYiMEJEkt+H4RuBz4OGyCe9kVaOjuOCshvwn+0KOadVwhWGMMRVGgYlAVWcBVwK9cTqi+Q/QCxjqvhY2o3u05CA1eSe7RzjDMMaYCiFgrSFVXQVcX0axBO3sZnU5W75nUvYArvN8SbTkhDskY4yJWMHUGiqXRkV/xA5twCc5ncMdijHGRLSITQQXRi2npexiYtYg1ApeGGNMsQW6oewaEalXlsEURZQoIz0fs0ZbsjjHulA2xpjiCnRE0Ax4V0QWisg4EekiIlJWgQXjMs9C6nOIf2cPCncoxhgTsQJdNfSUqvYGBgCrcMpMrBCRaSJyvYg0LKsgCxIrmYyI/pSvclLZmHNauMMxxpiIVGgbgaoeUdWZqvoHVW0PPAEk4FxOGnbXeb4kjnQmZg0s3RWPi/99MMaYCqzIjcWqul5V/6GqF4UioKKqLb8x1DOf2TnnslvrhjscY4yJOBF71ZCvmz0fowivZ/ULdyjGGBNxKkQiOC1qHxdHLebN7AvYo3YqxxhjiqLQRCAiE0SkXVkEUxJ3R79PBjG8lHVpuEMxxpiIEswRwUZgoogsEZHRIlIuf3K3iPqFoZ75TMvuw/achHCHY4wxESOYq4YmqWo3nJpDzYHV7iWk5a5P47uj30fI4dmsy8MdijHGRIyg2ghExAO0cYd9OPcV3Csib4UwtiI7RQ5wg+dzZuacxw85jcMdjjHGRIRg2giewTk9NAB4UlXPVtW/q+rFQPtClu0nIt+LyCYRGRtgvk4iki0iVxQa8a7vAl7bf2v0/1GddCZkDS10VcYYY4I7IlgLpLg3lH2b77UCS3+6RxEvAf2BtsA1InJSUSB3vr8DnwUddQB15Qi3RH/E5zmdWJnTqjRWaYwxFVowiWCYqh7znSAicwBU9VCA5ToDm1R1i6pmAG8Bl/iZ707gPWBPcCEX7mbPJ9TjEOOzriqtVRpjTIUVqPporIjUBeqLSB0RqesOzYFGQay7MbDdZ3yHO813G42BIcArgVYkIqNEZJmILNt7rPCa0zUknduiP+DrnES+zi73V74aY0xYBToi+AOwHKeBeIX7fDnwAc4pn8L4q1Sa/1v8OeABVc0OtCJVnaiqHVW1Y0JccAVQh3nm0Ih9PJ11lfVXYIwxAQSqPvq8qrYA7lfVFj5Diqq+GMS6dwC+JUGbALvyzdMReEtEtgFXAC+LyKVF2oMCxEomf4x+j1V6Op/ldCyNVRpjTIUU6NRQb/fpThG5LP8QxLqXAmeISAsRqQJcDcz2ncFNLM1VtTkwA7hNVWcVa0/8uMyzkFayk39kDSVby1VXCsYYU24EOjXUw3282M9QaE8wqpoF3IFzNdAG4B1VXefenTy6RFEHKVpyuC/6XX7UJryffX5ZbNIYYyKOaISdQO/YyKPLRtUIev4cFYZkPMYurcucqvdTS44XfaPjAl0cZYwx5Z+ILFdVv+fJg7mh7G4RqSWOSSKyQkT6ln6YoRElymMxk9lHPM9mFX6/mjHGVDbB3Edwk6oeBvoCDYAbgadCGlUpS4nawrWeuUzNvoj1OU3DHY4xxpQrwSSC3FbWAcBkVV2F/0tDy7Ux0W9Tm6M8lHkjOdZwbIwxXsEkguUi8jlOIvhMRGoCOaENq/TVlt8YGz2d5dqaGdndwx2OMcaUG8EkgpuBsUAnt9REFZzTQxHnCs8COsr3PJV1DQe1evALWif2xpgKLJj+CHKAX4G2ItIdaAfUDnFcIRElyuMxr3OI6jxtdYiMMQaA6MJmEJG/A1cB64HcUhAKLAhhXCFzVtR2bvB8xuTsfgz1fEVq1OZwh2SMMWEVzKmhS4HWqjpAVS92h8Ehjiuk7ol+jwQO8dfMm+yOY2NMpRdMItgCxIQ6kLJUU47z15g3WastmJbdJ9zhGGNMWBV6agg4Bqx0+yA4kTtRVe8KWVRl4OKoxbwV1Yuns66in+dbEuRwuEMyxpiwCOaIYDbwOPA/fi9FvTyUQZUFEXgsegrpVOWhzJusVLUxptIq9IhAVaeKSDWgqap+XwYxlZnTo3ZxX/Q7PJV1LTOyu3NldES2fxtjTIkEU2voYmAl8Kk7nioiswMuFEFu8XxEZ9nAo1nXsz0nIdzhGGNMmQvm1NA4nP6HDwKo6kqgRcgiKmMeUZ6p8i8A7s281a4iMsZUOsEkgiw/ndRXqDPqTWQfj8ZMZam24d/ZhXa1YIwxFUowiWCtiFwLeETkDBH5J07DcYVyWdRCBkQt4dmsK1mb0yzc4RhjTJkJJhHciVNW4gQwHTgM/DGEMYWFCPy/mNeowxH+mHk76Vqhbp0wxpgCBVNr6JiqPqiqnVS1o/s8vSyCK2t15CgTYl5hkzbhqaxrwh2OMcaUiYCJQERucHsk+80dlonI9WUVXDh096xhhOdTpmT3Y2F2YrjDMcaYkCswEbhf+H8E7gMaAY2BPwF3V/Rk8ED0W7SSndyfOZoDGnz/yMYYE4kCHRHcBgxR1XmqekhVD6rqXOBy97UKq5pk8HzMSxygJndm3kmWBtOUYowxkSnQN1wtVd2Wf6I7rVaoAiovEqO28UT06yzKSeLprKvDHY4xxoRMoBITx4v5WoUxNPor1moLJmYPol3UNi4Jd0DGGBMCgRLBWSKy2s90AVqGKJ5y56HoN9iYcxp/yhxFq52HSGxsXVYaYyoW0QLKbopIwLuqVPWnkERUiI6NPLpsVNk24O7TWlx84v8RVbsJs+/oRr0aVct0+8YYU1IislxVO/p7rcA2AlX9KdAQunDLn/pymIlVnmHf0RPcMe07MrNzwh2SMcaUGrscJkhJUVv522VJLN6SxpMfbwh3OMYYU2qC6aHMuC7r0IR1uw7z2qKttGsUzxVnNwl3SMYYU2LB9EcwSETsyAFgXDx/XtqNc1vV4y8z1/DNlrRwR2SMMSUWzBf81cCPIvK0iJwV6oDKu2jJ4aVrO9C0bhy3TF3G2p35K3QbY0xkCabo3HVAe2AzMFlEFovIKBGpGfLoyqk61avwxs2dqVUthhte/5Yte4+GOyRjjCm2oE75qOph4D3gLeBUYAiwQkTuDGFs5dqp8dV44+bOAAx/7Vt2H6oU99gZYyqgYNoIBovITGAuEAN0VtX+QApwf4jjK9daJtRg6k2dOXw8k+smLWH/bxnhDskYY4osmCOCK4BnVTVZVcer6h5w+ikAbgppdBEgsXE8k27oyI4Dxxkx+VuOnsgKd0jGGFMkwSSC3aq6wHeCiPwdQFXnhCSqCNOlZT1eurYD63YdZtR/lpGemR3ukIwxJmjBJIIL/UzrX9qBRLoL2jZkwpXJ/G9zGndN/46MLLv72BgTGQJ1THOriKwB2ojIap9hK+CvGJ2/dfQTke9FZJOIjPXz+jCf9f5PRFKKvyvhN6R9Ex4d3I7P1//KLf9ZxvEMOzIwxpR/ge4sngZ8AvwN8P0SP6Kq+wtbsYh4gJdwjih2AEtFZLaqrveZbSvQQ1UPiEh/YCLQpYj7UK7ccG5zqkZH8eeZa7j+9SW8NqITtWJjwh2WMcYUKNCpIXU7obkdOOIzICJ1g1h3Z2CTqm5R1QycS0/zlPRX1f+p6gF39BugQtRsuLpzU168pgMrtx/kmonfsO/oiXCHZIwxBQqUCKa5j8uBZe7jcp/xwjQGtvuM73CnFeRmnCOQk7g3sC0TkWV7j/kvm13eDEw+lVev78jmvUcZ+u/F7Dpo9xkYY8qnQGWoB7mPLVS1pfuYOwTTMY34W63fGUV64SSCBwqIZaKqdlTVjglx/lZbPvVs3YA3bu7C3sMnuPKVxWzd91u4QzLGmJMEaizuEGgIYt07gNN8xpsAu/xsJxmYBFyiqhWuilun5nWZPqorxzOzufKVxazfdTjcIRljTB6BeiibF2A5VdXeAVcsEg38APQBdgJLgWtVdZ3PPE1x7li+XlX/F0zA4eihrEDjgi84t2nPUYa/toTDxzN57ur2XNi2YQgDM8aYvAL1UFbgVUOq2qskG1XVLBG5A/gM8ACvq+o6ERntvv4K8DBQD3hZRACyCgo00p3eoAbv33Yuo/6znFFvLOO+C8/k9l6n4+63McaETaAjgt6qOldELvP3uqq+H9LIChCpRwS50jOzGfveamat3MXApFMZf2UycVWsfyBjTGgV64gA6IFz2uZiP68pEJZEEOliYzw8e1UqbRvV4qlPNrJl32+8ev3ZNKkTF+7QjDGVVIFHBOVVpB8R+Jr//R7unP4dMZ4oXh7Wga4t65VSYMYYk1egI4JgylDXE5EXRGSFiCwXkedFxL6xSkHP1g2YdXs3asfFcN2kJby2aCs5OZGVmI0xkS+YonNvAXuBy3FKUu8F3g5lUJVJq4QazLq9Gz1bJ/D4h+sZMWUpew6nhzssY0wlEkwiqKuqj6vqVnd4Aqgd4rgqlVqxMbx6fUcevzSRb7em0e/5hXy+7pdwh2WMqSSCSQTzRORqEYlyh6HAR6EOrLIREYZ3bcaHd57HqfGxjHpjOX9+fzXHMqyjG2NMaAW6fPQIztVBAlQHcgvsRwFHVbVWmUSYT0VqLC5IRlYOz3zxA/9esJkW9arz3NWpJDepHZJtGWMqh0CNxXbVUGkIUUJYvDmNe99Zyd4jJ7i1Zytu73U6sTGekGzLGFOxleiqIXcFdUSks4h0zx1KN0Tjzzmt6vHp3d0ZnNKIf87dRN9nFzDv+z3hDssYU8EEc/noSGABTqmIR93HcaENy+SKj4vhmatSmXZLF6I9wo2Tl3Lbf5fzyyG7ssgYUzqCOSK4G+gE/OTWH2qPcwmpKUPntqrPJ3efz/19z2TOhj30+cd8Ji3cQla29Y1sjCmZYBJBuqqmA4hIVVXdCLQObVgRZly8M4RY1WgPd/Q+gy/u6UGnFnV54qMNXPzi13y9aV/It22MqbiCSQQ7RKQ2MAv4QkQ+wE+/AqbsNK0Xx+QRnfjXsA4cPp7JsElLGP7aEtbsCE2jtTGmYivSVUMi0gOIBz51+yEuc+XyqqFcIbp6KJD0zGze/OYnXpq3iQPHMhmYfCr3921Ni/rVyzwWY0z5Vdzqo74r6ACch3NfwdfhSgLmZLExHkae35KhnU5j0oItTFq0lU/X/sJVnU7j7j5n0LBWbLhDNMaUc8FcNfQwMBWnA5n6wGQR+WuoAzNFUys2hnv7tuarMb0Y1qUp7yzdTven5/HQrLX8nHYs3OEZY8qxQk8NicgGoL1Pg3E1YIWqnlUG8Z3ETg0F56e033hp3iZmfreT7BxlQNKpjO7RisTGoW/UNsaUPyW9oWwb4Ht+oSqwuRTiMiHUrF51nr4ihUUP9OaW81sy//u9DPrnIoa/toRFP+4j0u4oN8aETqBaQ//EaRNoinMfwRfu+IXAIlW9uqyC9GVHBMVzOD2T/37zM69/vZW9R07Q9tRaXNe1GYNTG1GjqnWVaUxFV6xaQyJyQ6CVqurUUoityMp1IvCnnCWH9MxsZn63k6n/28bGX45QvYqHS9o35trOTe20kTEVWImLzolIFeBMd/R7Vc0sxfiKxBJB6VBVvtt+kGlLfubD1btIz8whpUk813ZpysUpjYirYkcJxlQkJUoEItIT56qhbTglqU8DblDVBaUaZZAiLhHkKqcJAeDQsUxmfreDad/+zA+/HiWuiocL2zZkcEojzj8jgSrRQdUmNMaUYyVNBMuBa1X1e3f8TGC6qp5d6pEGwRJB6Kgqy386wPvf7eTjNbs5eCyT2nEx9E88lcEpjejSoi5RURLuMI0xxVDSRLBaVZMLm1ZWIj4R5NYkKueJISMrh0Wb9jJ75S4+X/8rxzKyaVirKv0TT+WCsxrSuUVdO1IwJoKU9M7i5SLyGvCGOz4MWF5awZnyqUp0FL3bNKR3m4Ycy8hizoY9zF61i+nf/syU/22jZtVoerRO4IKzGtKzdQK146qEO2RjTDEFkwhGA7cDd+G0ESwAXg5lUKZ8iasSzcUpjbg4pRHHM7JZtGkfX67/lTkb9/Dh6t14ooROzevQ/cwEzju9Pu0axeOxU0jGRIyAp4ZEJApYraqJZRdSYBF7aii/cn5qKBg5OcqqHQf5csOvzNmwh42/HAEgvloM57SsR7cz6tOtVT1a1K+OiCUGY8Kp2KeGVDVHRFaJSFNV/Tk04ZlIaTfILypKaN+0Du2b1mHMRW3YcySdxZvT+HrTPr7elMan634B4NT4WDo1r0vH5nXo0LQObU6pSbTH2heMKS+CaSyei3Nn8bfAb7nTVXVwaEPzr8IcEfgTYYkgEFXlp7RjLNq0j8Wb01j2035+PXwCgOpVPKQ2rc3ZzepydrM6JDWOp251a2MwJpRK2lj8aCnHYyoBEaF5/eo0r1+d67o2Q1XZefA4y386wLJtB1j+0wFenPsjOe7vkMa1q5HYuBZJjeNJbBxPUuN46tWoGt6dMKaSCFRiIhanofh0YA3wmqpmlWFsflXoI4L8KtARgj9H0jNZs+MQa3Y6w9qdh9jmUzL71PhYzmxYkzan1OTMhjVpfUpNTm9Qg9gYTxijNiYyFfeIYCqQCSwE+gNtcTqyN6ZU1IyN4dzT63Pu6fW90w6nZ7Ju52HW7jzE+t2H+f6XIyzekkZGVg4AUeJUVj29QQ1aJlSnZf3qtKjvPK9XvYo1ShtTDIESQVtVTQJw7yP4tmxCMiETAY3StWJjOKdVPc5pVc87LSs7h5/2H+OHX46w8ZcjfP/LETbvPcpX3+8lIzvHZ9loWiTUoFndOJrWjeO0utU4rW4cp9WJ49T4WGugNqYAgRKBt7CcqmbZL60wGhdEVdBy/OVeUtGeKFol1KBVQg36J53qnZ6do+w8cJwt+46ydd9vbNn7G1v2HeW77Qf4aM1usnN+P+0ZHSU0ql2NRrVjaRRfjVNrx3JqfDVOjXceG9WOJb5ajB1RmEopUCJIEZHD7nMBqrnjAqiq1gp5dJVdMAkg/7wRVsqiJDxRQtN6cTStF0fP1nlfy8rOYfehdLbvP8bP+4+x/cAxft5/nN0Hj7Nk635+OZyeJ1GAczd1g5pVaVCzKg1rxTrPa8WSULMq9WtUoV71qtRzH6tVsXYKU3EUmAhU1T7pFVVBScI38UR4Aon2RDmnherGca6f17NzlL1HTrDr0HF2H0xn96Hj7Dlygj2H09lz5AQ/7jnK15v2cTjd//URcVU81KtRhbrVq1InLoY6cVWo7T7WiYuhtjseXy2GWrHOY83YaDs9ZcolKzpf0ZTkKKK480QgT5RwSnwsp8THOn3wFSA9M5u9R06Q9lsGaUdPkHY0g32/OY9pR3OnZ7Bpz1EOHsvk6InAF9ZVr+Jxk0IMNWKjqRkbTY2qzmPN2BhqVI2metVoqlfxOI9VPVSv4kyLq+Ihrko01ap4iKviIcaSiiklQXVMU+yVi/QDngc8wCRVfSrf6+K+PgA4BoxQ1RWB1lmpLh8tT4JJBPmTUFGSRwVJOBlZORw8nsHBY5kcOp7JoWOZHE53nh8+nuV9fiTdSRpH07M4kp7FEff58czsoLcV4xGqxfyeHGJjPMTGRBEb7XHHo9xpHqpGR1E12n2M+f15lego9zXneRWPx3mMjqKKJ4oq0UIVj4eYaCHGE0WMx5ke7RGio8TaVCJISW8oK+5GPcBLOH0c7wCWishsVV3vM1t/4Ax36AL8y3005U1RjjTyL+Pvyz3Y9fk7XVWUhFPYabCSJp5863HaGWJpUDO2WKvLzM7hWEY2v53I4lhGFkdPZHPsRBZHT2RxLCPbHbI4npHN8Uxn/HhGNscys0n3GfYeySI905knPTOHE1nZnMjK8V6GW1piPEJ0lJMYfk8QvyeKGE8Unigh2hNFTJS4zwVPVBTRueNRQpT76BFn2kmDOz3KfR7lneaUOokSZ1wE7zIiQpTgzO++FiVCVJT76B3wzus7n79HIe/8Th78fZrw+7z4PM9d1lmXs4wz7fflIO983vECXvem4PzLFLDOQEJ5aqgzsElVtzjByFvAJYBvIrgE+I86hyXfiEhtETlVVXeHMC5T1krrdFVB6wnqqqpizFOcxBNoe4Wtb1w8MUC8OwRef4D3p6CkFw05HiHjr2mcyMrhRKaTHE5k5ZDx4jlkEEMG0WSo+0g0mVe8QWZ2DpnZOWRkK5lZOWR8+ley8JDR469kZeeQufB5MnOiyUodRWaWkrniv2TjIQsPmW0uITtHyfzhC7LwkN2sOycyc8jKySY7R8nKUbJzcsjKUXLcce/jkT1kEUU2HrKJIjs6jhxVsnOUnNCdyKiUQnZqSESuAPqp6kh3fDjQRVXv8JnnQ+ApVV3kjs8BHlDVZfnWNQoY5Y4mAmtDEnTZqQ/sC3cQJRDp8UPk70Okxw+Rvw+RFn8zVU3w90Iojwj8HYvkzzrBzIOqTgQmAojIsoLOc0WKSN+HSI8fIn8fIj1+iPx9iPT4fYXysoMdOB3d52oC7CrGPMYYY0IolIlgKXCGiLQQkSrA1cDsfPPMBq4XR1fgkLUPGGNM2QrZqSG3LMUdwGc4l4++rqrrRGS0+/orwMc4l45uwrl89MYgVj0xRCGXpUjfh0iPHyJ/HyI9foj8fYj0+L1Ceh+BMcaY8s9uTTTGmErOEoExxlRyEZUIRKSfiHwvIptEZGy44ykKETlNROaJyAYRWSciEdnJj4h4ROQ79x6QiOPetDhDRDa6f4tzwh1TUYnIPe5naK2ITHd7Eyy3ROR1EdkjImt9ptUVkS9E5Ef3sU44YyxMAfsw3v0crRaRmSJSO4whlkjEJAKfkhW5vaVdIyJtwxtVkWQB96nqWUBX4PYIiz/X3cCGcAdRAs8Dn6pqGyCFCNsXEWkM3AV0VNVEnAsxrg5vVIWaAvTLN20sMEdVzwDmuOPl2RRO3ocvgERVTQZ+AP5c1kGVlohJBPiUrFDVDCC3ZEVEUNXduQX1VPUIzhdQ4/BGVTQi0gQYCEwKdyzFISK1gO7AawCqmqGqB8MaVPFE4/QPEg3EUc7vvVHVBcD+fJMvwekOF/fx0rKMqaj87YOqfu7Tj/s3OPdBRaRISgSNge0+4zuIsC/SXCLSHGgPLAlzKEX1HPAnoHQrl5WdlsBeYLJ7emuSiFQPd1BFoao7gQnAz8BunHtvPg9vVMXSMPeeIfexQZjjKambgE/CHURxRVIiCKocRXknIjWA94A/qurhwuYvL0RkELBHVZeHO5YSiAY6AP9S1fbAb5T/UxJ5uOfSLwFaAI2A6iJyXXijqtxE5EGcU7//DXcsxRVJiSDiy1GISAxOEvivqr4f7niKqBswWES24ZyW6y0ib4Y3pCLbAexQ1dwjsRk4iSGSXABsVdW9qpoJvA9+O2Er734VkVMB3Mc9YY6nWETkBmAQMEwj+KasSEoEwZSsKLfcTnheAzao6jPhjqeoVPXPqtpEVZvjvPdzVTWifomq6i/AdhHJ7eG4D3nLokeCn4GuIhLnfqb6EGEN3q7ZwA3u8xuAD8IYS7G4HW89AAxW1WPhjqckIiYRuI0yuSUrNgDvqOq68EZVJN2A4Ti/pFe6w4BwB1UJ3Qn8V0RWA6nAk+ENp2jco5kZwApgDc7/cLkudSAi04HFQGsR2SEiNwNPAReKyI84nVc9FWgd4VbAPrwI1AS+cP+fXwlrkCVgJSaMMaaSi5gjAmOMMaFhicAYYyo5SwTGGFPJWSIwxphKzhKBMcZUcpYITNiIyHwRuSjftD+KyMuFLBPSDsPdip6rReSefNMvLW6hQBFJLcrlwm6V1NuKOp+INBKRGaU1v6kcLBGYcJrOyZUzr3anh4WInAKcq6rJqvpsvpcvxal8WxypON2yBqs2UGgiyD+fqu5S1StKcX5TCVgiMOE0AxgkIlXBW4yvEbBIRP4lIsvcuvuP+ltYRI76PL9CRKa4zxNE5D0RWeoO3fwsGysik0VkjVuArpf70udAA/cGofN95j8XGAyMd19r5Q6fishyEVkoIm3cea90+wpYJSIL3DvhHwOucpe9Kl8s7UTkW/e11SJyBs4NVq3caeNFpIaIzBGRFW7MuZV388/XPLdmfpDr9Z3fIyIT3PWvFpE7g/5LmsimqjbYELYB+Ai4xH0+FhjvPq/rPnqA+UCyOz4fpxY/wFGf9VwBTHGfTwPOc583xSnrkX+79wGT3edtcEo3xALNgbUFxDoFuMJnfA5whvu8C07ZDXDu+G3sPq/tPo4AXixgvf/EqVUDUAWolj8OnIJ5tdzn9YFNOIUY88/nHQ9yvb7z34pTCyva929gQ8UfogvJE8aEWu7poQ/cx5vc6UNFZBTOF+CpOKdkVge5zguAtk4pHgBqiUhNdfqByHUezhclqrpRRH4CzgSCqggrThXZc4F3fbZT1X38GpgiIu/gFIUrzGLgQXH6e3hfVX/0Wad3k8CTItIdpwx4Y6BhKazX1wXAK+rW2FfV/H0ImArKEoEJt1nAMyLSAaimqitEpAVwP9BJVQ+4p3z8dcfoWx/F9/Uo4BxVPR5guwG/EYMQBRxU1dSTglIdLSJdcDrxWSkiJ82Tb/5pIrLEnf8zERkJbMk32zAgAThbVTPFqQIbsIvKINfrS4jA0u6m5KyNwISVqh7FOd3zOr83EtfC6SvgkIg0xOme1J9fReQsEYkChvhM/xynQCHgXLHjZ9kFOF+uiMiZOKeQvi8k3CM4RcZQpy+JrSJypbsOEZEU93krVV2iqg8D+3DKp3uXzU9EWgJbVPUFnKqcyX7mj8fpDyLTbc9olj+mYq7X1+fAaHF6PkNE6hbyfpgKwhKBKQ+m4/Qf/BaAqq4CvgPW4SSIrwtYbizwITAXp7euXHcBHd0Gz/XAaD/Lvgx4RGQN8DYwQlVPFBLnW8AYt3G5FU4iuVlEVrmx5jbgjncbXNfiJJxVwDyc01UnNRYDVwFrRWQlTnvFf1Q1DfjabXQej9PpSUcRWeZud6P7XuWfr6jr9TUJp61ktbtP1xbyfpgKwqqPGmNMJWdHBMYYU8lZIjDGmErOEoExxlRylgiMMaaSs0RgjDGVnCUCY4yp5CwRGGNMJff/AUSnK4ZrDKtYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Draw 10000 samples, each with 11 observations\n",
    "n_observations = 11\n",
    "n_samples = 10000\n",
    "y = norm_dist.rvs(size=(n_observations, n_samples))\n",
    "\n",
    "# calculate the value of the statistic for each sample\n",
    "null_dist_mc = statistic(y)\n",
    "\n",
    "# compare the MC and asymptotic distributions\n",
    "plt.plot(y_grid, pdf)\n",
    "plt.hist(null_dist_mc, density=True, bins=200)\n",
    "plt.xlim(0, np.max(y_grid))\n",
    "plt.xlabel(\"Value of test statistic\")\n",
    "plt.ylabel(\"Probability Density / Observed Frequency\")\n",
    "plt.legend(['Asymptotic Null Distribution', 'Monte Carlo Null Distribution (11 observations/sample)'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c15ba7a-53c4-42b4-833c-39d86e2ba969",
   "metadata": {},
   "source": [
    "This is because the asymptotic null distribution was derived under the assumption that the number of observations approaches infinity (hence the name \"asymptotic\" null distribution). Apparently, it is quite different from the null distribution of the test statistic when the number of observations is 11.\n",
    "\n",
    "The true theoretical null distribution when the number of observations is 11 may not be possible to calculate analytically (in a way that can be expressed in terms of a finite number of common functions). So rather than comparing a test statistic against a theoretical null distribution to determine the $p$-value, the approach of the *Monte Carlo test* is to compare the test statistic against the Monte Carlo null distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "17e06351-479d-49b6-989f-7a30bfde43ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Asymptotic p-value 0.03045746622458189\n",
      "Monte Carlo p-value: 0.0085\n"
     ]
    }
   ],
   "source": [
    "res = stats.jarque_bera(x)\n",
    "stat, p_asymptotic = res\n",
    "count = np.sum(null_dist_mc >= stat)\n",
    "p_mc = count / n_samples\n",
    "print(f\"Asymptotic p-value {p_asymptotic}\")\n",
    "print(f\"Monte Carlo p-value: {p_mc}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e27fd4f9-64eb-42ca-9131-b6f1cefa9698",
   "metadata": {},
   "source": [
    "These $p$-values are substantially different, so we might draw different conclusions about the validity of the null hypothesis depending on which test we perform. Under the 1% threshold used above, the Monte Carlo test would suggest that there is enough evidence for rejection of the null hypothesis whereas the asymptotic test performed by `stats.jarque_bera` would not. In other cases, the opposite may be true. In any case, it seems that the Monte Carlo test should be preferred when the number of observations is small.\n",
    "\n",
    "`stats.monte_carlo_test` simplifies the process of performing a Monte Carlo test. All we need to provide is the obverved data, a function that generates data sampled under under the null hypothesis, and a function that computes the test statistic. `monte_carlo_test` returns an object with the observed statistic value, an empirical null distribution of the statistic, and the corresponding $p$-value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4869f8dd-4d59-432e-80be-b6ec321aef26",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed values of the statistic: 6.982848237344646\n",
      "Monte Carlo p-value: 0.0076\n",
      "Empirical null distribution shape: (10000,)\n"
     ]
    }
   ],
   "source": [
    "res = stats.monte_carlo_test(sample=x, rvs=norm_dist.rvs, statistic=statistic, \n",
    "                             n_resamples=10000, alternative='greater')\n",
    "print(f\"Observed values of the statistic: {res.statistic}\")\n",
    "print(f\"Monte Carlo p-value: {res.pvalue}\")\n",
    "print(f\"Empirical null distribution shape: {res.null_distribution.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46763f1a-a696-41b9-ace3-9124760b7e65",
   "metadata": {},
   "source": [
    "Note that the $p$-value here is slightly different than `p_mc` above because the algorithm is stochastic. Nevertheless, `res.pvalue` is much closer to `p_mc` than to `p_asyptotic`: for small samples, the Monte Carlo null distribution generated from a sufficiently large number of random samples is often more accurate than the asymptotic null distribution, despite the error inherent in random sampling.\n",
    "\n",
    "Here, we also passed optional arguments to control the behavior of `monte_carlo_test`. As one might expect, the parameter `n_resamples` controls the number of samples to draw from the provided *random variate sample* function. Perhaps less obvious is the meaning of `alternative`, which controls which *alternative hypothesis* we are testing, that is, which tail of the null distribution the observed statistic value should be compared against. In this case, we are assessing the null hypothesis that the sample was drawn from a normal population against the alternative that the sample is drawn from a population which tends to produce a *greater* value of the test statistic. Perhaps more intuitively, the argument `'greater'` corresponds with the sign of the comparison against the null distribution (i.e., `null_dist_emperical >= statistic(x)` from above). Another options for `alternative` is `'less'` (i.e., `null_dist_emperical <= statistic(x)`). In some cases, a `'two-sided'` alternative is desired, which is twice the minimum of the \"one-sided\" $p$-values. (More on the choice of this convention below.)\n",
    "\n",
    "We can improve performance of `monte_carlo_test` by ensuring that our test statistic is \"vectorized\". That is, instead of requiring a one-dimensional sample array as input, the statistic should accept an $n$-dimensional array of samples in which each *axis-slice* (e.g. row, column) is a distinct sample. Our `statistic` function is already vectorized in some sense. Above, we wrote `null_dist_emperical = statistic(y)`, and `statistic` computed the statistic for each column of the two-dimensional `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6d0f8f55-5497-428c-a54e-8664c1cfd6d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11, 10000)\n",
      "(10000,)\n"
     ]
    }
   ],
   "source": [
    "print(y.shape)  # 10000 samples, each with 11 observations\n",
    "print(statistic(y).shape)  # statistic for each column of y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c4da9ef-f7b2-4540-a195-b322dd089692",
   "metadata": {},
   "source": [
    "However, `monte_carlo_test` requires that the statistic function accept an `axis` argument to compute the statistic along *any* axis. Only minor modifications are required:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b7d6797e-7dab-46ac-a097-ec41716cd8b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def statistic_vectorized(x, axis=0):\n",
    "    # Calculates the Jarque-Bera Statistic\n",
    "    # Compare against https://github.com/scipy/scipy/blob/4cf21e753cf937d1c6c2d2a0e372fbc1dbbeea81/scipy/stats/_stats_py.py#L1583-L1637\n",
    "    n = x.shape[axis]\n",
    "    mu = np.mean(x, axis=axis, keepdims=True)\n",
    "    x = x - mu  # remove the mean from the data\n",
    "    s = stats.skew(x, axis=axis)  # calculate the sample skewness\n",
    "    k = stats.kurtosis(x, axis=axis)  # calculate the sample kurtosis\n",
    "    statistic = n/6 * (s**2 + k**2/4)\n",
    "    return statistic\n",
    "\n",
    "np.testing.assert_allclose(statistic_vectorized(y, axis=0), statistic_vectorized(y.T, axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f5001f8-c828-42f2-848c-8d79157bd159",
   "metadata": {},
   "source": [
    "But `monte_carlo_test` becomes much faster:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "198b8734-c516-4394-af01-614712138c10",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.15 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# Before\n",
    "%timeit -r1 -n1 stats.monte_carlo_test(sample=x, rvs=norm_dist.rvs, statistic=statistic, n_resamples=10000, alternative='greater')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "20212a36-30ab-4486-81d8-fbe8683ecaa6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.8 ms ± 583 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "# After\n",
    "%timeit stats.monte_carlo_test(sample=x, rvs=norm_dist.rvs, statistic=statistic_vectorized, n_resamples=10000, alternative='greater', vectorized=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f376518b-23b7-45bf-9921-4f717378f0aa",
   "metadata": {},
   "source": [
    "When a statistical test is already implemented in SciPy (like `stats.jarque_bera`), it becomes even easier to perform a Monte Carlo version of the test. We simply need to \"wrap\" it in another function which only returns the test statistic rather than the full result object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4fc52618-5476-4042-9bea-8b463c2767f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.008"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def statistic_scipy(x):\n",
    "    # `jarque_bera` returns a result obeject\n",
    "    res = stats.jarque_bera(x)  \n",
    "    # Our wrapper returns only the statistic, as required by `monte_carlo_test`\n",
    "    return res.statistic\n",
    "\n",
    "res = stats.monte_carlo_test(sample=x, rvs=norm_dist.rvs, statistic=statistic_scipy, n_resamples=10000, alternative='greater')\n",
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ac75775-40e2-44c7-882b-009b2adef689",
   "metadata": {},
   "source": [
    "Of course, besides enabling more accurate tests for small sample sizes, `monte_carlo_test` makes it easy to perform hypothesis tests *not* implemented in SciPy. For instance, suppose we want to assess whether our data is distributed according to a [`rayleigh` distribution](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rayleigh). Just as the \"normal distribution\" is really a *family* of distributions parameterized by mean and standard deviations, so `rayleigh` is a family of distributions rather than one specific distribution. In contrast with the normal distribution, however, there are no tests in SciPy specifically designed to determine whether a sample is drawn from a Rayleigh distribution. \n",
    "\n",
    "The closest options are tests like [`ks_1samp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ks_1samp.html) and [`cramervonmises`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.cramervonmises.html), which can assess whether a sample is drawn from a *specific* distribution. If we want to use these tests to assess whether the sample is distributed according to *any* Rayleigh distribution, one approach is to fit a Rayleigh distribution to the data, and then apply `ks_1samp` and `cramervonmises` to test whether the data were drawn from the fitted Rayleigh distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7ecf4ac3-3430-4508-b42e-625e3d7f08cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEGCAYAAACtqQjWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABHpElEQVR4nO3dd3hUVfrA8e+bkEIPLUhP6BAglNCUJjZAFMvae2OxrLru+hN0Xcuq6+ruolgXFUUXxYIFV7CACFIlQXpNAkggQEgglPTk/f0xNziGlAnJZDLJ+3meeWbm3nPuvMfgvHPvOfccUVWMMcaYyhDg6wCMMcbUHJZUjDHGVBpLKsYYYyqNJRVjjDGVxpKKMcaYSlPH1wH4UvPmzTUiIsLXYRhjjF+Ji4s7pKotittXq5NKREQEsbGxvg7DGGP8iojsLmmfXf4yxhhTaSypGGOMqTSWVIwxxlSaWt2nYoypGrm5uSQlJZGVleXrUEw5hIaG0rZtW4KCgjyuY0nFGON1SUlJNGzYkIiICETE1+EYD6gqqampJCUlERkZ6XE9u/xljPG6rKwsmjVrZgnFj4gIzZo1K/fZpVeTioiMEZFtIhIvIpOL2S8iMs3Zv15E+jvb24nIIhHZIiKbROQ+tzpNReQ7EdnhPDdx2zfFOdY2EbnAm20zxpSPJRT/czp/M69d/hKRQOAV4DwgCVgtInNVdbNbsbFAF+cxGHjNec4D/qSqa0SkIRAnIt85dScDC1X1WSdRTQYeEpGewNVAFNAaWCAiXVU131ttNKcn9Xg2KxPT2JV6gvCGIbQOq0urxqG0alyXusGBvg7PGFMB3uxTGQTEq2oigIjMBiYA7kllAvCuuhZ1WSkiYSLSSlWTgWQAVT0mIluANk7dCcAop/5M4AfgIWf7bFXNBnaKSLwTwwovttF44GhWLj8lprE8IZXlCYfYuv9YiWXD6gXRvmk9LuvXhisHtqNesHX7mcqxf/9+7r//flavXk1ISAgRERG88MILdO3atULH3bVrF+PHj2fjxo3Exsby7rvvMm3atFPKFd5s3bx58xKP9cwzz/Dwww+X6/NHjRpFcnIyoaGhNGjQgBkzZtCtW7eT20NCQsjJyeHcc8/lqaeeIiwsDIDAwEB69+598jiff/45lTHDiDf/j20D7HF7n4TrLKSsMm1wEgqAiEQA/YBVzqaWTtJBVZNFJNztWCuLOdZviMhEYCJA+/bty9UgUz4Hj2bxxJebmb8xmQKFkDoBxEQ04cELujG0UzO6tWzIoePZ7DuSRXJ6JsnpWew7ksnGfUd5/MvNTF2wgxuGdODGMzsQ3jDU180xfkxVufTSS7npppuYPXs2AGvXruXAgQMVTiruYmJiiImJOe36p5NUAGbNmkVMTAzTp0/nwQcfZO7cub/ZnpOTw5QpU5gwYQKLFy8GoG7duqxdu/a0Yy2JN5NKcRfjii4zWWoZEWkAzAHuV9WjlfB5qOp0YDpATEyMLXvpBarKJ3FJ/O1/m8nKK+CO4R0Z1S2cfu3DCA367eWt+iF16NCs/inHiNudxvQlibzyQzzTlyRyWf823D48ks7hDauqGaYGWbRoEUFBQUyaNOnktr59+wJw/PhxJkyYwOHDh8nNzeWpp55iwoQJ7Nq1i7FjxzJs2DCWL19OmzZt+OKLL6hbty5xcXHceuut1KtXj2HDhp085g8//MA///lP/ve//5Gamso111xDSkoKgwYNwn2V3UsuuYQ9e/aQlZXFfffdx8SJE5k8eTKZmZn07duXqKgoZs2axX//+1+mTZtGTk4OgwcP5tVXXyUwsORLxCNGjOCFF144ZXtwcDDPPfccnTt3Zt26dURHR1f8P2oJvJlUkoB2bu/bAvs8LSMiQbgSyixV/dStzIHCS2Qi0go4WI7PM162Jy2Dhz/bwI87DjEwogn/uLwPHVs0KPdxBnRoyn9uaEpiynHeWrqTT+KSmL16DzefGcGUcd0JqWN9L/7qiS83sXlfWb8Ry6dn60Y8dlFUifs3btzIgAEDit0XGhrKZ599RqNGjTh06BBDhgzh4osvBmDHjh188MEHvPHGG1x55ZXMmTOH66+/nltuuYWXXnqJkSNH8uCDDxZ73CeeeIJhw4bx17/+la+++orp06ef3DdjxgyaNm1KZmYmAwcO5PLLL+fZZ5/l5ZdfPnn2sGXLFj788EOWLVtGUFAQd911F7NmzeLGG28ssZ1ffvnlby5puQsMDCQ6OpqtW7cSHR19MoEBREZG8tlnn5V43PLwZlJZDXQRkUhgL65O9GuLlJkL3OP0twwG0p1kIcBbwBZV/XcxdW4CnnWev3Db/r6I/BtXR30X4KfKb5YpTkGB8t7K3fzj660I8OSEKK4f3IGAgIqN+OnYogFPX9qbB87rykvfx/PO8l38/MthXr62P+2a1quc4E2tpqo8/PDDLFmyhICAAPbu3cuBAwcA15dt4RfvgAED2LVrF+np6Rw5coSRI0cCcMMNNzB//vxTjrtkyRI+/dT1e/jCCy+kSZOTA1WZNm3ayS/xPXv2sGPHDpo1a/ab+gsXLiQuLo6BAwcCkJmZSXh4OMW57rrrqFu3LhEREbz00kultrWQ313+UtU8EbkH+AYIBGao6iYRmeTsfx2YB4wD4oEM4Ban+lnADcAGEVnrbHtYVefhSiYfichtwC/AFc7xNonIR7g68/OAu23kV9XIzMnn9ndXsyw+lRFdW/DMpb1o26Ryv/CbNQjh8YujGNKxGQ9+vI7xLy3l31dGc06PlpX6Ocb7Sjuj8JaoqCg++eSTYvfNmjWLlJQU4uLiCAoKIiIi4uS9GSEhISfLBQYGkpmZiap6PNS2uHI//PADCxYsYMWKFdSrV49Ro0YVey+IqnLTTTfx97//vczPKew7KU1+fj4bNmygR48eHsV+urx6n4qqzlPVrqraSVWfdra97iQU1OVuZ39vVY11ti9VVVHVPqra13nMc/alquo5qtrFeU5z+7ynnWN1U9VTfzqYSpedl8/E92JZkZDKM5f2ZuYtAys9obgb0+sM/nfvMNo2qcttM2N5dv5W8vILvPZ5pmYYPXo02dnZvPHGGye3rV69msWLF5Oenk54eDhBQUEsWrSI3btLnNUdgLCwMBo3bszSpUsB1xd6cUaMGHFy3/z58zl8+DAA6enpNGnShHr16rF161ZWrvx1fFFQUBC5ubkAnHPOOXzyySccPOi6wp+WllZmbCXJzc1lypQptGvXjj59+pzWMTxld9Sb05aXX8B9H6zlxx2HePayPlw7uH2V3ODWoVl95tx5JtcObs/rixO49o1VHDhqc0qZkokIn332Gd999x2dOnUiKiqKxx9/nNatW3PdddcRGxtLTEwMs2bNonv37mUe7+233+buu+9m6NCh1K1bt9gyjz32GEuWLKF///58++23J0ebjhkzhry8PPr06cOjjz7KkCFDTtaZOHEiffr04brrrqNnz5489dRTnH/++fTp04fzzjuP5OTkYj+rJNdddx19+vShV69enDhxgi+++KLsShUk7tfYapuYmBi1RbpOT0GB8uAn65mzJolHx/fktmGezw1UmT77OYmHP91I84bBfPT7obRqXPz/4Ma3tmzZ4vXLLsY7ivvbiUicqhZ7vc3OVEy5qSpPfLmJOWuS+OO5XX2WUAAu7deW2ROHcORELte+sYqDdsZijE9ZUjHl9q9vtzNzxW7uGB7Jved09nU4RLcL451bB3LgaBbXvrmKQ8ezfR2SMbWWJRVTLv9ZnMDLi+K5ZlA7Hh7Xo9pMEjigQ1Nm3DyQpMMZXP/mKg6fyPF1SMbUSpZUjMeWbE/h7/O3Mr5PK566pHe1SSiFhnRsxps3DiTx0AlumLGK9MxcX4dkTK1jScV4JD0zl4fmrKdzeAP+eUU0gRW8qdFbhnVpzn+uH8C2/ce4acZPHMuyxGJMVbKkYjzy5JebOXgsm39dEX3K/F3Vzdndw3nl2v5s3JvObTNjycmz+1iMqSqWVEyZvtt8gDlrkrhrVCei24X5OhyPnB91Bv+6Mpqfdqbxt/9tLruCqfGSkpKYMGECXbp0oVOnTtx3333k5Lj63t555x3uueceH0d4qgYNip83LzAw8OTEk9HR0fz73/+moKD0H0+7du3i/fff90aYv2FJxZTq8Ikcpny6gR6tGvGH0V18HU65TOjbht+P6Mh7K3fz0eo9ZVcw1cPSF2Dnkt9u27nEtf00qSqXXXYZl1xyCTt27GD79u0cP36cRx55pEKhliYvL89rxy6ct2vTpk189913zJs3jyeeeKLUOpZUTLXw6BcbSc/M4d9XRhNcx//+uTx4QTeGdW7OXz7fyNo9R3wdjvFEm/7w8c2/JpadS1zv2/Q/7UN+//33hIaGcsstrukFAwMDmTp1KjNmzCAjIwNwTew4ZswYunXrdvIL+sSJE1x44YVER0fTq1cvPvzwQwDi4uIYOXIkAwYM4IILLjh5p/uoUaN4+OGHGTlyJE8//TQREREnzyAyMjJo164dubm5JCQkMGbMGAYMGMDw4cPZunWrq6k7dzJ06FAGDhzIo48+6lHbwsPDmT59Oi+//DKqyq5duxg+fDj9+/enf//+LF++HIDJkyfz448/0rdvX6ZOnVpiuQpT1Vr7GDBggJqSfblur3Z46H/68vc7fB1KhaQdz9aznl2og59eoAePZvk6nFpp8+bN5auQuFj1H5GqC59yPScurtDnv/jii3r//fefsr1v3766bt06ffvtt/WMM87QQ4cOaUZGhkZFRenq1av1k08+0dtvv/1k+SNHjmhOTo4OHTpUDx48qKqqs2fP1ltuuUVVVUeOHKl33nnnyfIXX3yxfv/99yfL3XbbbaqqOnr0aN2+fbuqqq5cuVLPPvtsVVW96KKLdObMmaqq+vLLL2v9+vWLbU9x28PCwnT//v164sQJzczMVFXV7du3a+H33KJFi/TCCy88Wb6kckUV97cDYrWE71X/++lpqsTBY1k8+vlGotuF8fsRHX0dToU0qR/Mf24YwJHMHO5+fw25NgFl9Rc5AmJugyXPuZ4jR1TocFrCzMLu28877zyaNWtG3bp1ueyyy1i6dCm9e/dmwYIFPPTQQ/z44480btyYbdu2sXHjRs477zz69u3LU089RVJS0sljXnXVVb95XXh2M3v2bK666iqOHz/O8uXLueKKK+jbty+///3vT57pLFu2jGuuuQZwTalf3jaCa/LIO+64g969e3PFFVeweXPxfYqelisvWwDcnEJVefjTjWTk5POvK6KpE+j/vz2iWjfmH5f34b7Za3n6qy08fnHVT79uymHnEoh9C0b8n+s5cniFEktUVBRz5sz5zbajR4+yZ88eOnXqRFxc3ClJR0To2rUrcXFxzJs3jylTpnD++edz6aWXEhUVxYoVK4r9rPr1f13J9OKLL2bKlCmkpaURFxfH6NGjOXHiBGFhYSWuZXI6938lJiYSGBhIeHg4TzzxBC1btmTdunUUFBQQGlr8UtxTp071qFx5+f+3hal08zfuZ8GWAzx4QTc6h5d/1cbqakLfNtw2LJJ3lu9iTlxS2RWMbxT2oVzxDox+xPXs3sdyGs455xwyMjJ49913AdfaIn/605+4+eabqVfPtVTDd999R1paGpmZmXz++eecddZZ7Nu3j3r16nH99dfz5z//mTVr1tCtWzdSUlJOJpXc3Fw2bdpU7Oc2aNCAQYMGcd999zF+/HgCAwNp1KgRkZGRfPzxx4DrR9y6desAOOuss5g9ezZQ8pT6RaWkpDBp0iTuueceRIT09HRatWpFQEAA7733Hvn5rmWlGjZsyLFjx07WK6lcRVlSMb+Rl1/AP7/ZRteWDbjlLN9NFOktU8Z2Z0jHpjz82Qa2HzhWdgVT9faucSWSwjOTyBGu93vXnPYhC6e+//jjj+nSpQtdu3YlNDSUZ5555mSZYcOGccMNN9C3b18uv/xyYmJi2LBhA4MGDaJv3748/fTT/OUvfyE4OJhPPvmEhx56iOjoaPr27VtqJ/dVV13Ff//7399cFps1axZvvfUW0dHRREVFnZyS/sUXX+SVV15h4MCBpKenl3hM97Xszz33XM4//3wee+wxAO666y5mzpzJkCFD2L59+8kzpz59+lCnTh2io6OZOnVqieUqyqtT34vIGOBFXCs/vqmqzxbZL87+cbhWfrxZVdc4+2YA44GDqtrLrc6HQDfnbRhwRFX7ikgEsAXY5uxbqaqTSovPpr4/1eyffmHypxuYfsMAzo86w9fheEXKsWzGvLCElo1C+fzus/xyVJu/sanv/Ve1mfpeRAKBV4CxQE/gGhHpWaTYWFxryXcBJgKvue17BxhT9LiqepU6q0ECc4BP3XYn6K8rRZaaUMypsnLzeWHBDvq1D+O8njV3md4WDUP4+2W92Zx8lBcWbPd1OMbUKN78iTYIiFfVRFXNAWYDE4qUmQC864xSWwmEiUgrAFVdAqRRAucs50rgA69EXwu9t2I3+49m8X8XdK92k0VWtvOjzuCqmHa8vjiB2F0l/jMzxpSTN5NKG8D9NuYkZ1t5y5RkOHBAVXe4bYsUkZ9FZLGIDC+ukohMFJFYEYlNSUnx8KNqvqNZubzyQzwjurZgaKdmvg6nSjx6UU/aNKnLAx+t43i29+5+Ni7evNRuvON0/mbeTCrF/dQtGqEnZUpyDb89S0kG2qtqP+AB4H0RaXTKwVWnq2qMqsa0aNHCw4+q+d5cksiRjFz+74JuZReuIRqE1OHfV/Zlz+EMnrL5wbwqNDSU1NRUSyx+RFVJTU0t91Bjb96nkgS0c3vfFth3GmVOISJ1gMuAAYXbVDUbyHZex4lIAtAVsJ74MqQcy+bNpTu5sE8rerVp7OtwqtTAiKZMGtmJ135I4JweLWt0X5IvtW3blqSkJOzqgH8JDQ2lbdu25arjzaSyGugiIpHAXuBq4NoiZeYC94jIbGAwkK6qyR4c+1xgq6qevNlARFoAaaqaLyIdcXX+J1ZCO2q8VxbFk51XwJ/O6+rrUHzij+d25YdtKUyes55+7UfQvEGIr0OqcYKCgoiMrHlD1M2pvHb5S1XzgHuAb3AN9f1IVTeJyCQRKRyZNQ/XF3888AZwV2F9EfkAWAF0E5EkEbnN7fBXc2oH/QhgvYisAz4BJqmq9cCWYU9aBrNW7ebKmLZ0bFFzbnQsj+A6AbxwVV+OZecxec4Gu0RjTAV49T6V6s7uU4EHPlrLV+uT+eHBUbRqXNfX4fjUmz8m8tRXW/jXFdFcPqB8p/zG1CY+uU/FVH/b9h/js5/3ctOZEbU+oQDcelYkAzo04amvNpN2IsfX4Rjjlyyp1GKv/hBP/eA63Dmyk69DqRYCAoS/X9ab49l5PPWVjQYz5nRYUqml9qdn8dX6ZK6MaUeT+sG+Dqfa6NqyIZNGduLTNXtZuuOQr8Mxxu9YUqml/rtyN/mq3HxmhK9DqXbuPrszkc3r8/BnG8jMqZyZW42pLSyp1EJZufnMWrWb83q0pH2zer4Op9oJDQrk6Ut78UtaBtO+31F2BWPMSZZUaqHPf97L4Yxcbh1m9w2U5MxOzbliQFumL0lkS/JRX4djjN+wpFLLqCozlu2kZ6tGDI5s6utwqrWHx/UgrG4Qkz/dQH5B7R16b0x5WFKpZZbFp7L9wHFuOSuixs9EXFFN6gfz14t6sm7PEd5bscvX4RjjFyyp1DJvL9tJ8wbBXBTd2teh+IWLo1szvEtznv9mG/uOZPo6HGOqPUsqtcjOQydYuPUg1w3uQGhQoK/D8QsiwtOX9CZflSe/tHtXjCmLJZVa5J1lOwkODOC6Ie19HYpfad+sHvec3ZmvN+3nxx02y64xpbGkUkukZ+bycVwSF0W3Jrxh+dZHMHD78I50aFaPx+duIievwNfhGFNtWVKpJT5avYeMnHxuOSvC16H4pdCgQB67qCcJKSd4Z/lOX4djTLVlSaUWyMsv4J3luxgU2bTWLcJVmUZ3b8k53cN5ccEODhzN8nU4xlRLllRqgQVbDrD3SCa3nmU3O1bUXy/qSW6+8vd5W3wdijHVkleTioiMEZFtIhIvIpOL2S8iMs3Zv15E+rvtmyEiB0VkY5E6j4vIXhFZ6zzGue2b4hxrm4hc4M22+ZOZy3fTtkldWyq3EnRoVp/fj+zI52v3sSox1dfhGFPteC2piEgg8AowFugJXCMiPYsUG4tr2d8uwETgNbd97wBjSjj8VFXt6zzmOZ/XE9eKkFFOvVedGGq1PWkZrEhM5eqB7QgMsJsdK8NdozrTJqwuj83dRF6+ddob486bZyqDgHhVTVTVHGA2MKFImQnAu+qyEggTkVYAqroEKM9ywBOA2aqarao7cS1RPKjCrfBzc9YkIQKX9reVDCtL3eBA/nJhD7buP8asVb/4OhxjqhVvJpU2wB6390nOtvKWKc49zuWyGSLSpDzHEpGJIhIrIrEpKTX7noOCAmXOmiTO6tScNmG2smNlGtPrDIZ1bs6/vt3GoePZvg7HmGrDm0mluGstRWfl86RMUa8BnYC+QDLwr/IcS1Wnq2qMqsa0aNGijI/ybz/tSmNPWia/s/XWK52I8PjFPcnIyef5r7f5Ohxjqg1vJpUkoJ3b+7bAvtMo8xuqekBV81W1AHiDXy9xlftYNd0ncUk0CKnDBVFn+DqUGqlzeENuOSuCj+L2sHFvuq/DMaZa8GZSWQ10EZFIEQnG1Yk+t0iZucCNziiwIUC6qiaXdtDCPhfHpUDh6LC5wNUiEiIikbg6/3+qjIb4oxPZeczbkMz4Pq2oG1zrxyt4zR/O6ULTesE8+eVmVG16fGPKTCoi8k8RiSrvgVU1D7gH+AbYAnykqptEZJKITHKKzQMScXWqvwHc5fa5HwArgG4ikiQitzm7nhORDSKyHjgb+KPzeZuAj4DNwNfA3apaa9eCnb9xPxk5+Xbpy8sahQbxwPld+WlXGvM37vd1OMb4nJT160pEbgduAeoAbwMfqGqNONePiYnR2NhYX4fhFVdPX8H+9CwW/XmUrZviZfkFyoXTfuR4dh4LHhhpM0CbGk9E4lQ1prh9ZZ6pqOqbqnoWcCMQAawXkfdF5OzKDdNUlj1pGaxMTON3A9paQqkCgQHCX8f3JOlwJm8ttXnBTO3mUZ+KcxNhd+dxCFgHPCAis70YmzlNdm9K1Tuzc3PO79mSVxbFc9DmBTO1mCd9Kv8GtgLjgGdUdYCq/kNVLwL6eTtAUz52b4rvPHJhD3LzC3juGxtibGovT85UNgLRqvp7VS06mqrW37Fe3di9Kb7ToVl9bj0rkk/iklifdMTX4RjjE54kletUNcN9g4gsBKgpHfY1id2b4lv3jO5M8wY2xNjUXiUmFREJFZGmQHMRaSIiTZ1HBNC6yiI0HrN7U3yvYWgQfzq/G7G7D/O/9aXecmVMjVTamcrvgThcnfNrnNdxwBe4Zh821Yzdm1I9XBnTjh6tGvHs/K1k5dbaW6VMLVViUlHVF1U1Evizqka6PaJV9eUqjNF46JO4PUQ0q8eADk3KLmy8pnCI8d4jNsTY1D6lXf4a7bzcKyKXFX1UUXzGQ/vTs1iZmMal/ezelOpgaKdmnNezJa8uiufgMRtibGqP0i5/jXSeLyrmMd7LcZlymr/Rdf3+wj6tyihpqsrD43qQnVfA1O+2+zoUY6pMnZJ2qOpjzvMtVReOOV3zNiTTrWVDOoc38HUoxhHZvD43Do3gneU7uXFoBD1aNfJ1SMZ4nSc3P94nIo2cmYTfFJE1InJ+VQRnPHPgaBaxuw8zrredpVQ3957TmYahQTz1lQ0xNrWDJ/ep3KqqR4HzgXBck0s+69WoTLnM35CMKlzYx+5NqW7C6gVz/7ldWBafyqJtB30djjFe50lSKez1HQe8rarrKH6VReMj8zbsp2vLBnQOb+jrUEwxrh/SgY7N6/PUV1vIzS/wdTjGeJUnSSVORL7FlVS+EZGGgP2fUU0cPJrF6t1pdumrGgsKDGDKuB4kppzg/VW/+DocY7zKk6RyGzAZGOhM1xKM6xKYqQbmb9zvuvRlSaVaO7dHOGd2asYLC7aTnpHr63CM8RpP1lMpAA4APUVkBBAFhHlycBEZIyLbRCReRCYXs19EZJqzf72I9HfbN0NEDorIxiJ1nheRrU75z0QkzNkeISKZIrLWebzuSYz+7qsNyXQJb0CXlnbpqzoTER65sAdHMnN56fsdvg7HGK/xZPTXP4BlwF+AB53Hnz2oF4hrOpexQE/gGhHpWaTYWFxryXcBJgKvue17BxhTzKG/A3qpah9gOzDFbV+CqvZ1HpOKqVujHDyaxepddunLX0S1bswVA9oyc8Uudh064etwjPEKTy5/XQJ0U9VxqnqR87jYg3qDgHhVTVTVHGA2MKFImQnAu+qyEggTkVYAqroESCt6UFX9VlXznLcrgVo70dXXm5xLX3bDo9/48/ndCAoM4O/zt/g6FGO8wpOkkggEncax2wB73N4nOdvKW6Y0twLz3d5HisjPIrJYRIYXV0FEJopIrIjEpqSklOOjqp+v1ifTObwBXe3Sl98IbxTKnSM78c2mA6xISPV1OMZUOk+SSgawVkT+4/R/TBORaR7UK27YcdG7vzwpU/zBRR4B8oBZzqZkoL2q9gMeAN4XkVNuYVbV6aoao6oxLVq08OSjqqWDx7L4yS59+aU7RnSkdeNQnvpqMwUFdkOkqVk8SSpzgb8By/l1+vs4D+olAe3c3rcF9p1GmVOIyE245h+7Tp3blFU1W1VTnddxQALQ1YM4/dI3mw7YqC8/FRoUyENju7Np31HmrEnydTjGVCpPRn/NBD4CVqrqzMKHB8deDXQRkUgRCQauxpWg3M0FbnRGgQ0B0lW11JWNRGQM8BBwsfuKlCLSwhkcgIh0xNX5n+hBnH5p3vpkOrWoT9eWNteXP7o4ujV924Xx/DfbOJGdV3YFY/yEJ6O/LgLWAl877/uKSNHkcAqnM/0e4BtgC/CRqm4SkUkiUjgyax6uL/544A3gLrfP/QBYAXQTkSQRuc3Z9TLQEPiuyNDhEcB6EVkHfAJMUtVTOvprgkPHs1m1M5ULe7eyae79lIjw6PieHDyWzX+W1NjfPqYWKnGWYjeP4xrJ9QOAqq4VkUhPDq6q83AlDvdtr7u9VuDuEupeU8L2ziVsnwPM8SQuf/f1xv0UKIyzUV9+bUCHJozv04rpSxK4ZlA7WjWu6+uQjKkwT/pU8lQ1vcg26130oXkbkunYoj7dbNSX35s8tjsFCs99vc3XoRhTKTxJKhtF5FogUES6iMhLuDrtjQ+kHs9mZaJd+qop2japx+3DIvns572s23PE1+EYU2GeJJU/4JqaJRv4ADgK3O/FmEwpFm49SIHCBVE2zX1NcdfZnWneIJi//c/WXDH+z5PRXxmq+oiqDnTu73hEVW3RbR9ZuOUArRqHEtXaVhGsKRqE1OFP53cjdvdh5m3Y7+twjKmQUpOKiNzkrPR4wnnEisiNVRWc+a2s3Hx+3HGI0d3D7dJXDXNlTDt6tGrEM/O2kJWb7+twjDltJSYVJ3ncD/wJaI1r+pT/A+6zxOIbKxNTycjJ59weLX0diqlkgQHCX8f3ZO+RTN6wIcbGj5V2pnIXcKmqLlLVdFU9oqrfA5fjdj+JqToLtxykblAgQzs183UoxguGdmrG2F5n8OoPCexPtyvMxj+VllQaqequohudbXZBv4qpKgu3HGB4l+aEBgX6OhzjJQ+P60G+Kv/4equvQzHmtJSWVDJPc5/xgs3JR9mXnmWXvmq4dk3rccdw1xDjNb8c9nU4xpRbaUmlh7O6YtHHBqB7VQVoXBZuOYgInN093NehGC+7a1RnwhuG8MSXNoux8T+lTdPSo8qiMGVauOUAfduF0aJhiK9DMV5WP6QOD43pzp8+Xsfna/dyWf9auw6d8UMlnqmo6u7SHlUZZG134GgW65LS7dJXLXJpvzZEtwvj2flbbRZj41c8uaPe+Nj3Ww8CcE4Pu/RVWwQECI9d5JrF+LUfEnwdjjEes6TiBxZuOUCbsLo2gWQt0799Ey7t14bpPyayJy2j7ArGVAOerKcyXkQs+fhIZk4+S+MPcW4Pu4u+NnpoTHcCRXjqq82+DsUYj3iSLK4GdojIcyJSrs57ERkjIttEJF5EJhezX5w17+OdkWX93fbNEJGDIrKxSJ2mIvKdiOxwnpu47ZviHGubiFxQnlirq2Xxh8jKLeDcntafUhud0TiUe0Z35ptNB1iyPcXX4RhTJk8mlLwe6Idrzfe3RWSFiEwUkVKvxThL+74CjAV6AteISM8ixcbiWva3CzAReM1t3zvAmGIOPRlYqKpdgIXOe5xjX41rRuUxwKuFywv7s4VbD9AgpA6DI+0u+trq9uGRRDSrx+NfbiInr8DX4RhTKo8ua6nqUVyrKs4GWgGXAmtE5A+lVBsExKtqoqrmOHUnFCkzAXhXXVYCYSLSyvnMJUBxywFPAGY6r2cCl7htn62q2aq6E9cSxYM8aV91VVCgLNxykBFdmxNcx65A1lYhdQJ57OIoElNOMGPZTl+HY0ypPOlTuVhEPgO+B4KAQao6FogG/lxK1TbAHrf3Sc628pYpqqWqJgM4z4VDok7nWNXaxn3pHDyWzTnd7dJXbXd2t3DO7dGSaQt3kJxuE1qY6suTn7+/A6aqah9VfV5VD4JrnRXg1lLqFderXPT2YE/KeMqjYzmX7mJFJDYlpXpfo16w5SABdhe9cTx2UU/yCpRn5tm8YKb68iSpJDuXok4SkX8AqOrCUuolAe3c3rcF9p1GmaIOFF4ic54PludYqjrdWWwspkWLFmV8lG8t2HyAAR2a0LR+sK9DMdVAu6b1uHNkJ75ct48VCam+DseYYnmSVM4rZttYD+qtBrqISKSIBOPqRJ9bpMxc4EZnFNgQIL3w0lYp5gI3Oa9vAr5w2361iISISCSuzv+fPIizWtp3JJPNyUc5x+6iN27uHNWJtk3q8tjcjeTmW6e9qX5KW6TrzsLJI4tMKLkTWF/WgVU1D7gH+AbYAnykqptEZJKITHKKzQMScXWqv4HbOi0i8gGwAugmIkkicpuz61ngPBHZgSvhPet83ibgI2Az8DVwt6r67RJ6i7Y5d9HbpS/jJjQokL+O78n2A8d5d4XNlmSqH1EtvgtDRBoDTYC/4wzbdRxT1eJGZfmdmJgYjY2N9XUYxZr4biyb9h1l6UNn202P5jdUlZvfXs2a3YdZ+OeRhDcM9XVIppYRkThVjSluX2mXv9RZkOtu4JjbAxFpWtlBml/l5BWwPCGVkd1aWEIxpxARHr84iuy8Ap61TntTzZSWVN53nuOAWOc5zu298ZI1vxzmeHYeI7pU74EExncim9fnjhGRfPrzXpYnHPJ1OMacVNrU9+Od50hV7eg8Fz46Vl2Itc+S7SnUCRDO7Gx30ZuS/WF0F9o3rcdfPttIdp7fdh+aGqa0jvr+pT2qMsjaZvH2FPp3aEKj0CBfh2KqsdCgQP52SS8SD52w6fFNtVHayo//KmWfAqMrORYDHDyWxaZ9R3nwgm6+DsX4gZFdW3BRdGteXZTARdGt6dSiga9DMrVciUlFVc+uykCMy4/bXdfHR3a1/hTjmUfH9+CHbQf5y2cbef+OwTa4w/hUaZe/RjvPlxX3qLoQa5clO1Jo3iCYnq0a+ToU4yfCG4by0JjurEhM5dM1e30djqnlSrv8NRLXJJIXFbNPgU+9ElEtll+gLNmewtndwgkIsF+bxnPXDmrPnDVJPD1vC6O7h9PEpvYxPlLa5a/HnOdbqi6c2m3j3nQOZ+Qywi59mXIKCBCeubQ3419ayt/nb+G530X7OiRTS3ky9X0zZ3XGNSISJyIvioiNdfWCJdtTEIHhXZr7OhTjh3q0asTtwyP5KDaJVYk24aTxDU8mlJwNpACX45oGPwX40JtB1VaLt6fQu01jmjUI8XUoxk/dd04X2oTV5eHPNti9K8YnPEkqTVX1b6q603k8BYR5Oa5aJz0jlzW/HLa76E2F1Auuw9OX9iIh5QTTFu7wdTimFvIkqSwSkatFJMB5XAl85e3AaptlCYcoUBjZzZKKqZhR3cK5vH9bXl+cyIakdF+HY2qZ0oYUHxORo8Dvcc0DluM8ZgN/rJrwao8l21NoGFqHfu3CfB2KqQH+Or4nzeoH8+An68jJs3VXTNUpbe6vhqrayHkOUNU6ziNAVe0mikqkqizensKwzs2pE+jJyaMxpWtcL4hnLu3N1v3HeHlRvK/DMbWIR99gItJERAaJyIjCh7cDq012HDxOcnqWDSU2lercni25pG9rXl0Uz6Z9dhnMVA1PhhTfDizBtYLjE87z454cXETGiMg2EYkXkcnF7BdnuHK8s6pk/7LqisiHIrLWeewSkbXO9ggRyXTb97onMVYHi7elAFhSMZXusYuiCKsXzIMfr7flh02V8ORM5T5gILDbmQ+sH65hxaUSkUDgFVzr2fcErhGRnkWKjcW1lnwXYCLwWll1VfUqVe2rqn2BOfz2zv6Ewn2qOgk/sWRHCl3CG9AmrK6vQzE1TJP6wTx1SRSbk4/yus1kbKqAJ0klS1WzAEQkRFW3Ap5MoTsIiFfVRFUt7OCfUKTMBOBddVkJhIlIK0/qimvWvCuBDzyIpdrKyMljVWKanaUYrxnTqxUX9mnFtO93sG3/MV+HY2o4T5JKkoiEAZ8D34nIF8A+D+q1Afa4H8fZ5kkZT+oOBw6oqvtg/EgR+VlEFovI8OKCEpGJIhIrIrEpKWWecHndqsQ0cvILbFZi41VPXhxFw9AgHvxkHXl2Gcx4UZlJRVUvVdUjqvo48CjwFnCJB8cubkZE9bCMJ3Wv4bdnKclAe1XtBzwAvC8ip4xSU9XpqhqjqjEtWvj+i3zx9hRC6gQwKLKpr0MxNVizBiE8OSGK9UnpvLLILoMZ7/F09Fd/EbkX6AMkOZekypIEtHN735ZTz3BKKlNqXRGpA1yG23QxqpqtqqnO6zggAejqQZw+tSz+EIMimxIaFOjrUEwNN75Pay7p25pp3+9gzS+HfR2OqaE8Gf31V2Am0AxoDrwtIn/x4NirgS4iEikiwcDVwNwiZeYCNzqjwIYA6aqa7EHdc4GtqprkFmcLp4MfEemIq/M/0YM4febA0Sx2HDzOsM42gaSpGk9e0oszGoXyxw/Xcjw7z9fhmBrIkzOVa4CBqvqYMx3+EOC6siqpah5wD64hyFuAj1R1k4hMEpHCkVnzcH3xxwNvAHeVVtft8Fdzagf9CGC9iKwDPgEmqWqaB+3zmWXxrlUez7KkYqpIo9Agpl7Vlz1pGTz55aayKxhTTqUt0lVoFxAKZDnvQ3BdWiqTqs7DlTjct73u9lqBuz2t67bv5mK2zcE1xNhvLI0/RJN6QbbKo6lSgyKbcueoTryyKIGzu4UztncrX4dkapASk4qIvISrczwb2CQi3znvzwOWVk14NZeqsiz+EGd2bm6rPJoqd/+5XflxxyEmf7qBfu2bcEbjUF+HZGqI0i5/xQJxwGfAw8Ai4AfgEWC+1yOr4RJSjnPgaLb1pxifCAoM4IWr+pKTV8CfPl5LQUHRwZXGnJ7SlhOeWfja6SwvHEm1TVVzvR1YTbd0h6s/xZKK8ZWOLRrw14t6MuXTDcxYtpPbh3f0dUimBvBk9NcoYAeuaVNeBbbbhJIVtywhlfZN69GuaT1fh2JqsasHtuO8ni157uttNumkqRSejP76F3C+qo5U1RHABcBU74ZVs+XlF7AyIdVGfRmfExH+cXkfmtQP4u5ZaziaZRchTMV4klSCVHVb4RtV3Q4EeS+kmm/93nSOZefZpS9TLTStH8zL1/Znz+FM/u/j9bgGZRpzejxJKnEi8paIjHIeb+DqwDenadmOQ4jA0E7NfB2KMQAMjGjK5DHd+XrTft5autPX4Rg/5klSmQRsAu7FNQ3+ZmebOU1L4w8R1boRTesH+zoUY066fXgk5/dsybPztxK7q1rfN2yqsVKTiogEAHGq+m9VvcyZXHKqqmZXUXw1TkZOHmt+OWz9KabaERGevyKaNk3qcs/7P3PouP1vbsqv1KSiqgXAOhFpX0Xx1Hg/7UwjN1+tP8VUS43rBvHqdf05nJHD/bPXkm/3r5hy8uTyVytcd9QvFJG5hQ9vB1ZTLYs/RHCdAAZG2FT3pnqKat2Yv03oxdL4Q7y4cEfZFYxx48ncX094PYpaZGl8KjEdmthU96Zau3JgO1bvSuOl73fQr30YZ3cL93VIxk+UeKYiIqEicj9wBdAdWKaqiwsfVRVgTXLoeDZbko9af4rxC09O6EX3Mxpx7wc/E3/wuK/DMX6itMtfM4EYYAMwFtdNkKYCliekAjY1i/EPdYMDeePGAYTUCeC2mas5fMKTtflMbVdaUumpqter6n+A3+FaE95UwLIdh2gUWodebRr7OhRjPNK2ST3+c0MMyelZTPpvHDl5tr69KV1pSeXkfA3OolmmAlSVpfGHOLNTcwJtqnvjRwZ0aMJzl/dh1c40Hv18o91xb0pVWlKJFpGjzuMY0KfwtYgc9eTgIjJGRLaJSLyITC5mv4jINGf/ehHpX1ZdEXlcRPaKyFrnMc5t3xSn/DYRucCz/wRVY3dqBnuPZHJWF7v0ZfzPJf3a8IfRnfkwdg9v/mh33JuSlTb1fYWGJznrxb+Ca1GvJGC1iMxV1c1uxcbiWku+CzAYeA0Y7EHdqar6zyKf1xPXMsNRQGtggYh0VdX8irSjsiyNt6nujX/747ldSUg5zjPzt9CxRX3O6dHS1yGZasiT+1RO1yAgXlUTVTUHmA1MKFJmAvCuuqwEwkSklYd1i5oAzFbVbFXdiWvd+0GV2aCKWJ5wiDZhdYloZlPdG/8UECD864q+9GrdmHs/+JktyR5dsDC1jDeTShtgj9v7JGebJ2XKqnuPc7lshog0KcfnISITRSRWRGJTUlLK057TVlCgLE9IZWinZohYf4rxX64RYTE0CK3Dre+sZu+RTF+HZKoZbyaV4r49i/bwlVSmtLqvAZ2AvkAyvw519uTzUNXpqhqjqjEtWrQopkrl27L/KEcycjnTZiU2NcAZjUOZcfNAjmflccNbq0i1OcKMG28mlSSgndv7tsA+D8uUWFdVD6hqvjMv2Rv8eonLk8/ziRXO/Sk21b2pKaJaN+atmwey93AmN7+9mmO2uJdxeDOprAa6iEiks8b91UDROcPmAjc6o8CGAOmqmlxaXafPpdClwEa3Y10tIiEiEomr8/8nbzWuPFYkpNKxeX1aNa7r61CMqTSDIpvy2vX92ZJ8lDvejSUrt1qMiTE+5rWk4tzbcg/wDbAF+EhVN4nIJBEpXI9lHpCIq1P9DeCu0uo6dZ4TkQ0ish44G/ijU2cT8BGu9V6+Bu6uDiO/8vILWLUzjSF2lmJqoNHdW/LPK6JZmZjGPe//TF6+3RxZ20ltvpEpJiZGY2NjvfoZa/cc4ZJXlvHytf0Y36e1Vz/LGF+ZuXwXj83dxOX92/L87/oQYDf41mgiEqeqMcXt82SWYlMByxNc96cM6WhnKqbmuunMCI5k5DJ1wXYa1w3i0fE9bKRjLWVJxctWJKTSrWVDmjcI8XUoxnjVved05khmDjOW7SRA4JELLbHURpZUvCgnr4DVu9K4eqAtnGlqPhHhr+N7ogpvLt1Jdl4BT1wcZZfCahlLKl60ds8RsnILbCixqTVEhMcu6klIUAD/WZxITl4Bz1zW2yZRrUUsqXjR8oRDiMCQSEsqpvYQESaP6U5InUCmLdxBTn4Bz/+uD3UCvXkHg6kuLKl40YqEVHq1bkzjekG+DsWYKiUiPHBeV0LqBPD8N9vIySvghav7EmSJpcazpOIlWbn5/PzLEW4+K8LXoRjjM3ef3ZmQOgE89dUWsvMKePnafoQGVWgCdFPN2c8GL4nbfZicfOtPMeb24R3524QoFmw5wPVvriLNliWu0SypeMnyhEPUCRAGRjT1dSjG+NwNQyN45dr+rN+bzmWvLiMx5bivQzJeYknFS5YnpNKnbWMahNgVRmMALuzTig/uGMKxrDwue205P+1M83VIxgssqXjB8ew81ielc2YnW+XRGHcDOjThs7vOomn9YK5/cxWf/7zX1yGZSmZJxQtW70wjv0CtP8WYYrRvVo9P7zyTfu3DuP/Dtby4YAe1eQ7CmsaSihcsTzhEcGAAAzo0KbuwMbVQWL1g3rttMJf1b8PUBdu55/2fbU2WGsKSihesSEylf4cwGzppTCmC6wTwryuieWhMd77etJ+LX17G5n227r2/s6RSyY5k5LBp31GGdrT+FGPKIiLcOaoT798+mBPZeVz66jI+XP2LXQ7zY15NKiIyRkS2iUi8iEwuZr+IyDRn/3oR6V9WXRF5XkS2OuU/E5EwZ3uEiGSKyFrn8bo321aSlYlpqMKZna0/xRhPDe7YjHn3DScmogkPzdnAnz9eT0ZOnq/DMqfBa0lFRAKBV4CxQE/gGhHpWaTYWFzL/nYBJgKveVD3O6CXqvYBtgNT3I6XoKp9ncckfGBlYip1gwKJbhvmi483xm81bxDCu7cO5r5zuvDpz0lc8soy4g8e83VYppy8eaYyCIhX1URVzQFmAxOKlJkAvKsuK4EwZw36Euuq6rfOcsMAK4G2XmxDuS1POERMRBOC69iVRWPKKzBA+ON5XXn31kGkHs9h3LSl/GdxAvkFdjnMX3jzm68NsMftfZKzzZMyntQFuBWY7/Y+UkR+FpHFIjK8uKBEZKKIxIpIbEpKimct8dCh49lsP3DchhIbU0HDu7Rg/v3DObtbC/4+fyuXvbacHQfsrMUfeDOpFLeAQtGfGyWVKbOuiDwC5AGznE3JQHtV7Qc8ALwvIo1OOYjqdFWNUdWYFi1alNGE8lmZmArAUFs62JgKC28YyuvXD+Cla/rxS+oJLpy2lFd/iCcvv8DXoZlSeDOpJAHt3N63BfZ5WKbUuiJyEzAeuE6dYSKqmq2qqc7rOCAB6FopLfHQysRU6gcH0rtN46r8WGNqLBHhoujWfPfASM7tGc5zX2/jsteWs22/nbVUV95MKquBLiISKSLBwNXA3CJl5gI3OqPAhgDpqppcWl0RGQM8BFysqhmFBxKRFk4HPyLSEVfnf6IX23eKFQmpDIxsaosRGVPJmjcI4dXrBvDKtf3ZeziTcdN+5PG5mziSYTMeVzdem+1QVfNE5B7gGyAQmKGqm0RkkrP/dWAeMA6IBzKAW0qr6xz6ZSAE+E5EAFY6I71GAE+KSB6QD0xS1Sqbse7gsSwSUk5wZUy7sgsbY07LhX1aMbRTM/757TbeXbGLz9fu5YHzunLtoPb2Y66akNp8k1FMTIzGxsZWyrHmrtvHvR/8zBd3n0V0u7BKOaYxpmRbko/y5JebWZGYSteWDXh0fE+Gd6ncflJTPBGJU9WY4vZZaq8kKxNTaRhSh6jWp4wNMMZ4QY9WjXj/jsG8fv0AsnILuOGtn7h95mqb6sXHLKlUkpXWn2JMlRMRxvQ6g2//OIKHxnRnVWIa46b9yB3vxrJxb7qvw6uV7BuwEhw4mkXioRM2lNgYHwkNCuTOUZ1YOnk095/bhVWJqYx/aSm3z1zN+qQjvg6vVrGkUgkK708ZYknFGJ9qXDeI+8/tytLJo/nTeV1ZveswF7+8jFve/onl8YdsosoqYGvdVoKViak0DK1DT+tPMaZaaBQaxB/O6cLNZ0Xw3srdvPnjTq59cxVdwhtw49AOXNq/rS317SV2plIJViSkMjiyKYEBxU0EYIzxlYahQdw1qjPLJ4/mn1dEExoUyKNfbGLIMwt5fO4mElKO+zrEGsdSdQUlp2eyKzWD64d08HUoxpgShAYF8rsBbbm8fxvW7jnCuyt28/6qX3hn+S4GRTRlQr/WXNi7FWH1gn0dqt+zpFJBqxJd91daf4ox1Z+I0K99E/q1b8LD43rwUewePl2TxCOfbeTxuZsY1S2cS/u1YXT3cFu59TRZUqmgFQmpNK4bRM9W1p9ijD9p0TCEu8/uzF2jOrFp31E+/3kvc9ft47vNB2gYUofzerbknB4tGdG1OQ1Dg3wdrt+wpFJBK3emMiiyKQHWn2KMXxIRerVpTK82jZkyrgcrElL5fO1eFmw5wKc/7yUoUBgc2YxzeoRzbo+WtGtaz9chV2uWVCpg35FMdqdmcNPQCF+HYoypBIEBwrAuzRnWpTl5+QWs+eUIC7ccYMGWAzzx5Wae+HIzncMbMLRjM4Z0bMbgjk1p3iDE12FXK5ZUKsDuTzGm5qoTGMCgyKYMimzKlHE92HXoBAu3HmTJ9hQ+XZPEeyt3A9C1ZQOGdGzGoMim9G0XRpuwujiT3dZKllQqYEVCKmH1guh+RkNfh2KM8bKI5vW5bVgktw2LJDe/gI1701mZmMbKxFQ+iUvi3RWuJNOsfjB92jamT9swotu5nmvT2YwllQpYudN1f4r1pxhTuwQFBpwcRXbnqE7k5hewed9R1icdYV1SOuuTjvDD9hQKb+Bv0TCEri0b0LVlw5OPLi0b0KgGDgCwpHKakg5nsCctk9vOivR1KMYYHwsKDCC6XRjR7cK4wdl2IjuPjXvT2bA3nW37j7H9wDFm/7SHzNz8k/VaNAwholk9OjSrT0SzerR3nts1qUdYvSC/vIxmSeU0rSy8P6WT9acYY05VP6QOgzs2Y7Bbn2tBgbL3SCbb9h9j24Fj7Dp0gt2pGfy4I4VP4rJ/Uz80KIDWjetyRuNQWjWuS+uwUM5oHErzBiE0bxBCiwYhtGgYQt3g6nU/jVeTirP074u4Vm98U1WfLbJfnP3jcK38eLOqrimtrog0BT4EIoBdwJWqetjZNwW4DdfKj/eq6jfeatuKhFSa1g+ma3gV9qcsfQHa9IfIEb9u27kE9q6BYfdXbv3iyn55Hyhw8Yu/7gdYNg3Outf1uvBYxR23op9fnrZWhuoQQ3WMxZSP298uIEBo17Qe7dJjOTdoDVxx/8limTn5/JKWwc5DJ9h7JJPkI5kkp2exLz2TZfGHOHgsi4Ji5sOsHxxI84YhhNULpkm9IMLqBhFWL5iwekE0qRdMw9A6NAwNcp7r0Cg0iAYhrtfeWKrDa0nFWS/+FeA8IAlYLSJzVXWzW7GxuNaS7wIMBl4DBpdRdzKwUFWfFZHJzvuHRKQnrrXso4DWwAIR6aqq+XjBykQf9Ke06Q8f3wxXvOP6ctm55Nf3lV2/uLIbP3Xt6325a//s61zvRz706+urZ5V83Ip+fnnaWhmqQwzVMRZTPh7+7eoGB9LtjIZ0K2HgT15+ASnHszl0LIdDx7Ndr49nk3Ism9TjORzOyCH1eA4JKcc5ciKXY9l5pYY1ttcZvHb9gEppojtvnqkMAuJVNRFARGYDEwD3pDIBeFdd81GvFJEwEWmF6yykpLoTgFFO/ZnAD8BDzvbZqpoN7BSReCeGFZXdsD1pGew9ksnvR3as7EOXLnKE6x/ixzdDzG0Q+9av/1Aru35xZa+e5dpXuK1QlttiSDt/LPm4Ff388rS1MlSHGKpjLKZ8KulvVycwgFaN69KqcV2PyufmF5CemcuxrDyOZRU+//q6vZdu4vRmUmkD7HF7n4TrbKSsMm3KqNtSVZMBVDVZRMLdjrWymGP9hohMBCYCtG/fvhzN+VV2Xj4XRLXkzE7NT6t+hUSOcP3DXPIcjPi/8n+plKd+SWXdt0Hxr0s6bmV8flWqDjFUx1hM+fjgbxcUGHCy/6UqeXPq++KuCxW9IlhSGU/qns7noarTVTVGVWNatGhRxiGL1zm8If+5IYbO4Q1Oq36F7Fzi+qUz4v9czzuXeK9+cWXdt6163fUo+rq041b086tadYihOsZiyqc2/e1U1SsPYCjwjdv7KcCUImX+A1zj9n4b0Kq0uoVlnNetgG3FHR/4BhhaWowDBgxQv5K4WPUfka7n4t5XZv3iyj7T1vVIXPzb98teOnVfccet6OeXp62VoTrEUB1jMeVTA/92QKyW9N1f0o6KPnBdWksEIoFgYB0QVaTMhcB8XGcZQ4CfyqoLPA9Mdl5PBp5zXkc55UKceolAYGkx+l1S+XFq8V/UP06t/PrFlZ17r+oX9/52f+Ji1fcu//V14bGKO25FP788ba0M1SGG6hiLKZ8a+LcrLamIenHNZhEZB7yAa1jwDFV9WkQmAajq686Q4peBMbiGFN+iqrEl1XW2NwM+AtoDvwBXqGqas+8R4FYgD7hfVeeXFl9MTIzGxsZWapuNMaamE5E4VY0pdp83k0p1Z0nFGGPKr7SkYmvUG2OMqTSWVIwxxlQaSyrGGGMqjSUVY4wxlaZWd9SLSAqw28sf0xw45OXP8AVrl/+pqW2zdlW9Dqpa7N3jtTqpVAURiS1plIQ/s3b5n5raNmtX9WKXv4wxxlQaSyrGGGMqjSUV75vu6wC8xNrlf2pq26xd1Yj1qRhjjKk0dqZijDGm0lhSMcYYU2ksqVSAiMwQkYMisrGYfX8WERWR5m7bpohIvIhsE5ELqjba8impbSLyByf+TSLynNt2v2hbce0Skb4islJE1opIrIgMctvnL+1qJyKLRGSL87e5z9neVES+E5EdznMTtzrVvm2ltOt5EdkqIutF5DMRCXOr47ftctvvv98fJc2Jbw+P1owZAfQHNhbZ3g7XImG7gebOtp78dr2XBMpY76W6tQ04G1gAhDjvw/2tbSW061tgrPN6HPCDH7arFdDfed0Q2O7E/xy/XX/oH/7UtlLadT5Qx9n+j5rSLue9X39/2JlKBajqEiCtmF1Tgf/jt8sZTwBmq2q2qu4E4oFBxdStFkpo253As6qa7ZQ56Gz3m7aV0C4FGjmvGwP7nNf+1K5kVV3jvD4GbAHa4GrDTKfYTOAS57VftK2kdqnqt6qa5xRbCbR1Xvt1u5zdfv39YUmlkonIxcBeVV1XZFcbYI/b+yR+/UfkL7oCw0VklYgsFpGBznZ/b9v9wPMisgf4J66lqcFP2yUiEUA/YBXQUlWTwfVFBoQ7xfyubUXa5e5WXCvIgp+3qyZ8f9TxdQA1iYjUAx7BdWp+yu5itvnbeO46QBNcSz8PBD4SkY74f9vuBP6oqnNE5ErgLeBc/LBdItIAmINr5dOjrsVViy9azLZq27ai7XLb/giulV5nFW4qprpftAtXO/z++8POVCpXJ1zXO9eJyC5cp+RrROQMXL8s2rmVbcuvl1n8RRLwqbr8BBTgmvTO39t2E/Cp8/pjfr2s4FftEpEgXF9Qs1S1sD0HRKSVs78VUHjJ0m/aVkK7EJGbgPHAdep0PODf7aoZ3x++7tTx9wcQQZGOerd9u/i1oy2K33a0JVJNO9pKahswCXjSed0V1+m4+FvbimnXFmCU8/ocIM7f/mbO3+Fd4IUi25/ntx31z/lT20pp1xhgM9CiyHa/bleRMn75/eHzAPz5AXwAJAO5uH5J3FbSPwrn/SO4Rm1swxltVF0fxbUNCAb+C2wE1gCj/a1tJbRrGBDn/E+7Chjgh+0ahutyyHpgrfMYBzQDFgI7nOem/tS2UtoVj+tHTeG212tCu4qU8cvvD5umxRhjTKWxPhVjjDGVxpKKMcaYSmNJxRhjTKWxpGKMMabSWFIxxhhTaSypGOMhEZkqIve7vf9GRN50e/8vEXmglPpPisi5ZXzG4yLy52K2h4nIXaXUq+tMnRMoIhGFszCLyM0i8nIx5YNFZImI2KwaplJZUjHGc8uBMwFEJADXbAJRbvvPBJaVVFlV/6qqC07zs8OAEpMKrvmvPlXVfE8Opqo5uO5bueo04zGmWJZUjPHcMpykgiuZbASOiUgTEQkBegA/i8gA56whzjmbKZwm5R0R+Z3zepyzHshSEZkmIv9z+5yeIvKDiCSKyL3OtmeBTs6aL88XE9t1wBclxN1ORL521uF4zG375049YyqNnfoa4yFV3ScieSLSHldyWYFrptihQDquu6MVeAmYoKopInIV8DSuMwkARCQU+A8wQlV3isgHRT6qO661axoC20TkNVxTrPRS1b5F4xKRYKCjqu4qIfRBQC8gA1gtIl+paiyupDiwhDrGnBZLKsaUT+HZypnAv3EllTNxJZXlQDdcX+DfOTMEB+KaFsZddyBRXetigGvqmIlu+79S15o12SJyEGhZRkzNgSOl7P9OVVMBRORTXFOExKpqvojkiEhDda3pYUyFWVIxpnwK+1V64/qlvwf4E3AUmIFrosBNqjq0lGOUOB+9I9vtdT5l/3+aCYSWsr/oXEzu70OArDKOb4zHrE/FmPJZhmu69TRVzVfVNFyd6ENxXQ7bBrQQkaHgmt5cRKKKHGMr0NFZnAk86yw/huty2ClU9TAQ6FxWK8554lqrvi6ulR+XObE1A1JUNdeDzzfGI5ZUjCmfDbguN60ssi1dVQ85o6p+B/xDRNbhmn32TPcDqGomrpFcX4vIUuAArstnJXIuXy0TkY0ldNR/i+uyVnGWAu85scxx+lPA1W8zr7TPNaa8bJZiY3xARBqo6nFxdby8AuxQ1akVOF4/4AFVvaEcdT4FpqjqttP9XGOKsjMVY3zjDhFZC2wCGuMaDXbaVPVnYJGIBHpS3hkx9rklFFPZ7EzFGGNMpbEzFWOMMZXGkooxxphKY0nFGGNMpbGkYowxptJYUjHGGFNp/h9AUSF6/lgbyAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dist_family = stats.rayleigh\n",
    "params = dist_family.fit(x)\n",
    "dist_specific = dist_family(*params)\n",
    "\n",
    "z = np.linspace(dist_specific.ppf(0), dist_specific.isf(0.001))\n",
    "plt.plot(z, dist_specific.pdf(z))\n",
    "plt.plot(x, np.zeros_like(x), 'x')\n",
    "plt.legend(('Candidate PDF', 'Observed Data'))\n",
    "plt.xlabel('Weight (lb)')\n",
    "plt.ylabel('Probability Density')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf116179-6679-4ca4-b938-549886b92b0d",
   "metadata": {},
   "source": [
    "To the eyes of the author, this does not look like a terrific fit. The mode of the Rayleigh distribution is too far to the right compared to the cluster of observations around 160 lb. Also, according to this Rayleigh distribution, there is zero probability that any weights could be less than ~135 lb, which does not seem realistic. However, the `ks_1samp` and `cramervonmises` tests are both inconclusive, with relatively large $p$-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "09860ff5-3a29-4d41-bf9f-d69770ab9aa3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KstestResult(statistic=0.26884627441317677, pvalue=0.3412228239139401)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.ks_1samp(x, dist_specific.cdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "09867af4-cd95-4286-8eb7-3e79e5b97b7a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CramerVonMisesResult(statistic=0.17536330558707267, pvalue=0.32395064743536117)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.cramervonmises(x, dist_specific.cdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7666ae1d-8c0b-48df-b40d-06050da64edc",
   "metadata": {},
   "source": [
    "A much more powerful test of the null hypothesis that the data is distributed according to *any* Rayleigh distribution is the [Anderson-Darling Test](https://en.wikipedia.org/wiki/Anderson%E2%80%93Darling_test). [`scipy.stats.anderson`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.anderson.html) implements the test for some families of distributions, but not for the Rayleigh distribution. A simple implementation is included below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "e9eb2df7-4560-4b17-98f5-e9611a9223a5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0425042504250425"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def statistic(x):\n",
    "    \"\"\"Compute the Anderson-Darling statistic A^2\"\"\"\n",
    "    # fit a distribution to the data\n",
    "    params = dist_family.fit(x)\n",
    "    dist = dist_family(*params)\n",
    "    \n",
    "    # compute A^2\n",
    "    x = np.sort(x)\n",
    "    n = len(x)\n",
    "    i = np.arange(1, n+1)\n",
    "    Si = (2*i - 1)/n * (dist.logcdf(x) + dist.logsf(x[::-1]))\n",
    "    S = np.sum(Si)\n",
    "    return -n - S\n",
    "\n",
    "params = dist_family.fit(x)\n",
    "dist = dist_family(*params)\n",
    "res = stats.monte_carlo_test(x, rvs=dist.rvs, statistic=statistic, alternative='greater')\n",
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db20efd8-6b9f-41d2-a44b-3b26e847ab22",
   "metadata": {},
   "source": [
    "Although this does not meet the threshold for significance used above (1%), it does begin to cast doubt on the null hypothesis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb114e39-ea07-451d-b6d1-1c965fca311e",
   "metadata": {},
   "source": [
    "As we can see, `monte_carlo_test` is a versatile tool for comparing a sample against a distribution by means of an arbitrary statistic. Provided a statistic and null distribution, it can replicate the $p$-value of any such tests in SciPy, and it may be more accurate than these existing implementations, especially for small samples:\n",
    "\n",
    "- [`skewtest`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skewtest.html)\n",
    "- [`kurtosistest`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kurtosistest.html)\n",
    "- [`normaltest`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.normaltest.html)\n",
    "- [`shapiro`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html)\n",
    "- [`anderson`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.anderson.html)\n",
    "- [`chisquare`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html)\n",
    "- [`power_divergence`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.power_divergence.html)\n",
    "- [`cramervonmises`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.cramervonmises.html)\n",
    "- [`ks_1samp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ks_1samp.html)\n",
    "- [`binomtest`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binomtest.html)\n",
    "\n",
    "In addition, `monte_carlo_test` can be used to perform tests not yet implemented in SciPy, such as [the Lilliefors Test](https://www.tandfonline.com/doi/abs/10.1080/01621459.1967.10482916) for normality.\n",
    "\n",
    "However, there are other types of statistical tests that do not test whether a sample is drawn from a particular distribution or family of distributions, but instead test whether  multiple samples are drawn from the same distribution. For these situations, we turn our attention to [Permutation Tests](https://nbviewer.org/github/scipy/scipy-cookbook/blob/main/ipython/ResamplingAndMonteCarloMethods/resampling_tutorial_2.ipynb)."
   ]
  }
 ],
 "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}