{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "57c13277-a239-452b-a202-e499e0a068b0",
   "metadata": {},
   "source": [
    "### Independent-Sample Tests\n",
    "#### Two-sample Test\n",
    "In [Individual Comparisons by Ranking Methods](https://www.jstor.org/stable/3001968#metadata_info_tab_contents), Wilcoxon considers two sprays designed to kill flying insects. A subset of the data, the percentage of flies killed in repeated trials of two treatments, is recorded below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fe264441-ebf5-47f0-83fb-81b824751b2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "x = np.array([61, 62, 67, 63, 56, 58])\n",
    "y = np.array([60, 68, 59, 72, 64])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8525c260-6699-4f90-b20c-85e0d948bb3a",
   "metadata": {},
   "source": [
    "In the paper, Wilcoxon describes a test to assess whether the two samples are drawn from the same population that is now commonly described as a *nonparametric* version of the independent sample t-test - that is, a version of the t-test that does not make a normality (or any particular distributional) assumption. \n",
    "\n",
    "Given samples `x` and `y`, Wilcoxon introduces a statistic which is proportional to an empirical estimate of the probability that a random observation from the distribution underlying `x` will be less than a random observation from the distribution underlying `y`. Suppose we want to test the null hypothesis that the samples are drawn from the same distribution against the alternative that they are drawn from different distributions which tend to produce samples that give lower values of the statistic. Under certain assumptions, this can be argued as evidence that the location of the distribution underlying `x` is less than the location of the distribution underlying `y`. To perform this test, we pass the data into [`scipy.stats.mannwhitneyu`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mannwhitneyu.html) with `alternative='less'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9fc1c550-67a2-4483-b760-e119b16e6b2d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.16450216450216448"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "_, pvalue = stats.mannwhitneyu(x, y, alternative='less')\n",
    "pvalue  # p-value is greater than our threshold; test is inconclusive"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b709400-3601-4ce5-9a36-538498abab45",
   "metadata": {},
   "source": [
    "Like the mean comparison test in Efron's example from the previous tutorial on [Permutation Tests](https://nbviewer.org/github/scipy/scipy-cookbook/blob/main/ipython/ResamplingAndMonteCarloMethods/resampling_tutorial_2.ipynb), this is an example of an \"independent sample\" test of the null hypothesis that group labels (`x`, `y`) are entirely random. In fact, because `mannwhitneyu` claims to produce an exact $p$-value, we would expect `permutation_test` to return precisely the same $p$-value (using `mannwhitneyu` only to compute the statistic)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f526ac5a-b0a8-4c60-a47d-324c22215431",
   "metadata": {},
   "outputs": [],
   "source": [
    "def statistic(x, y):\n",
    "    # return just the Mann-Whitney U statistic\n",
    "    return stats.mannwhitneyu(x, y, alternative='less').statistic\n",
    "\n",
    "# \"independent\" is the default `permutation type`, so we are not required to pass it here\n",
    "# We pass `alternative='less'` because lesser values of the statistic are more extreme\n",
    "res = stats.permutation_test((x, y), statistic, permutation_type='independent', alternative='less')\n",
    "np.testing.assert_allclose(res.pvalue, pvalue, atol=1e-15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c346fa8-1b34-4443-8220-c018b6d0f5a8",
   "metadata": {},
   "source": [
    "Just as with `monte_carlo_test`, vectorizing the `statistic` function can greatly improve the speed of the test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4f2dae73-2bf1-4c66-89b3-ac1ff0487e4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "272 ms ± 4.61 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# Before\n",
    "%timeit stats.permutation_test((x, y), statistic, alternative='less')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b5c8a5c2-e958-4164-bd19-0b9fffc342b5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56.5 ms ± 859 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "# After \n",
    "def statistic_vectorized(x, y, axis=0):\n",
    "    # return just the Mann-Whitney U statistic\n",
    "    return stats.mannwhitneyu(x, y, axis=axis, alternative='less').statistic\n",
    "\n",
    "%timeit stats.permutation_test((x, y), statistic_vectorized, alternative='less', vectorized=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a751fc04-e796-4d8e-bf5a-3c47da9e4323",
   "metadata": {},
   "source": [
    "Although `mannwhitneyu` provides an exact $p$-value for the data above, `permutation_test` comes in handy when there are ties in the samples. As the [`mannwhitneyu` documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mannwhitneyu.html) states,\n",
    "> `'exact'`: computes the exact p-value by comparing the observed statistic against the exact distribution of the statistic under the null hypothesis. **No correction is made for ties.**\n",
    "\n",
    "The complete data set in Wilcoxon's original paper had ties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ccadcdf2-1341-48c6-bcc3-92be7b80b725",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.014763014763014764 0.01351981351981352\n"
     ]
    }
   ],
   "source": [
    "x = [60, 67, 61, 62, 67, 63, 56, 58]\n",
    "y = [68, 68, 59, 72, 64, 67, 70, 74]\n",
    "res1 = stats.mannwhitneyu(x, y, method='exact', alternative='two-sided')\n",
    "# By default, only 9,999 random permutations are used. \n",
    "# We pass n_resamples=np.inf to ensure that all 12,870 possible permutations are used\n",
    "res2 = stats.permutation_test((x, y), statistic_vectorized, alternative='two-sided', vectorized=True, n_resamples=np.inf)\n",
    "print(res1.pvalue, res2.pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c85b266-314c-490f-9a22-a5c1e1400496",
   "metadata": {},
   "source": [
    "The two $p$-values are similar despite the ties, but only `permutation_test` is truly \"exact\" in this case. Either way, our 1% threshold for statistical significance is not met, and the test is inconclusive."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cfb7ac5-6ec6-43a4-ae3b-f5ef32050f17",
   "metadata": {},
   "source": [
    "#### Multi-sample Test\n",
    "`scipy.stats.kruskal` is a many-sample extension of the Mann-Whitney U test, but SciPy provides only an approximate (asymptotic) $p$-value. It is possible to perform an exact version of the test using `permutation_test` for very small samples, and a randomized test using a subset of the possible permutations may yield more accurate results than the approximation implemented by `kruskal`, especially if there are ties or the sample size is small. Using the (artificial) data for milk cap production from [Kruskal and Wallis' original paper](https://www.tandfonline.com/doi/abs/10.1080/01621459.1952.10483441), we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d8879fa7-9576-4d84-acfb-85962ebb0b02",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KruskalResult(statistic=5.656410256410254, pvalue=0.059118869289796136)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = [340, 345, 330, 342, 338]\n",
    "y = [339, 333, 344]\n",
    "z = [347, 343, 349, 355]\n",
    "stats.kruskal(x, y, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "246840a9-f7fe-4af6-942e-ad97540e5730",
   "metadata": {},
   "source": [
    "At the expense of some time, the exact p-value for this data is given by `permutation_test`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5b274945-72d2-4597-a2dd-e9d7efbe293e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.048629148629148626"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def statistic(x, y, z, axis=0):\n",
    "    return stats.kruskal(x, y, z, axis=axis).statistic\n",
    "\n",
    "res = stats.permutation_test((x, y, z), statistic, vectorized=True, alternative='greater', n_resamples=np.inf)\n",
    "res.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc877591-9f75-4990-a035-a6a48b4e4597",
   "metadata": {},
   "source": [
    "Note that we passed `alternative='greater'` into `permutation_test` but not into `kruskal`. This is because the `kruskal` test is inherently one-sided: data generated under the null hypothesis tends to generate small positive values of the statistic with greater values always being more exceptional. This raises the point that setting up a permutation test requires some study of both the underlying statistic and SciPy's implementation. Another example of this is shown in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eb247ee-5dc9-4951-9124-11bd3a88c510",
   "metadata": {},
   "source": [
    "#### Gotchas\n",
    "Suppose that we wish to perform the two-sample Kolmogorov-Smirnov test to test the null hypothesis that two samples were drawn from the same distribution against the alternative that the distribution $X$ underlying sample `x` is [stochastically greater](https://en.wikipedia.org/wiki/Stochastic_ordering) than the distribution $Y$ underlying sample `y`. Roughly speaking, this is the alternative that $X$ \"tends to be\" greater than $Y$.\n",
    "\n",
    "Here, we'll use randomly-generated data that best illustrates some confusing (but important) points. We choose shapes of the samples to generate the `RuntimeWarning` reported in [gh-14019](https://github.com/scipy/scipy/issues/14019)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "08329d40-cd24-482b-9a40-24e2ed90b23e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABDv0lEQVR4nO3dd3hcV5n48e87TSNp1HuxLfca98ROcSokdnohpAdIQja7S4clQGgLywK7/AIBkg2BQEggPSEJ6dWkufduy5Ysq/feZ87vjzs2ku7Ilm1JI43ez/PMY2veO3deXc28c+acc88VYwxKKaVGP0e4E1BKKTU4tKArpVSE0IKulFIRQgu6UkpFCC3oSikVIbSgK6VUhNCCPghE5EER+d4g7Wu8iDSLiDP480oRuWMw9h3c32si8pnB2l8kE5FCEflEuPM4GhE5V0SKh3D/g/baVkNPC/oxBN/UbSLSJCL1IvKxiNwlIkeOnTHmLmPMjwe4r6MWCGNMkTHGZ4zxD0LuPxSRv/TZ/wpjzJ9Pdt8hnusREekMfhjVishbIjKjRx5dwWPYJCJ7ReS3IpLV4/Hnikgg+PjDt78Pdp6DRURyReQ5EakWkQYR2SYinw13Xscj+Df7rz735YmIEREXDO5rWw09LegDc5kxJg6YAPwMuBt4eLCf5PCbaBT7H2OMD8gFKoFHesSeCh7DZOAqIBPY0LOoA6XBD7PDt8uGK/ET8BhwCOs1kQLcClSENaMIFQHvi2GjBf04GGMajDEvAdcBnxGROdC7pSMiqSLycrA1XysiH4iIQ0QeA8YDfw+2Pr/ZozV0u4gUAe/2bSEFTRaRtcGW4Isikhx8LtvX7cMtJRFZDnwHuC74fFuC8SNdOMG8visiB0WkUkQeFZGEYOxwHp8RkaJgS/SeAR6nVuBxYE6IWJcxZkfwGFYBXz/W/kTkNBFZLyKNIlIhIvf2s11S8NhXiUhd8P+5PeIrReTHIvJR8JvCmyKS2iN+S/BY1Azgdz0VeMQY02KM6TbGbDLGvNZjX8+ISHnwb/a+iMzuEXtERB4Qq/urOZhPpoj8Kpj3bhFZ0GP7QhH5tojsDMb/JCLefo5BdvCbQ5WIFIjIl451fI/mRF/bwe0vF5Edwe1XisjMHvtdKCKbgn+HZ0TkqR7Pc66IFIvI3SJSDvxpgH/b/xLrG3SziPxdRFJE5K/B1806EckLbisi8svga75BRLZK8L082mlBPwHGmLVAMbAsRPjrwVgakIFVVI0x5hagCKu17zPG/E+Px5wDzAQu6ucpbwVuA7KBbuDXA8jxdeC/sVrGPmPMvBCbfTZ4Ow+YBPiA3/bZ5ixgOnAB8P2eb8r+iIgPuAnYdJT8/MCLhD6Gfd0H3GeMiQcmA0/3s50D+BNWq3k80Ib997kR+ByQDniAbwRzngX8H3AL1nFOwfqm0Z/VwP0icr2IjA8Rfw2YGnyejcBf+8Q/DXwXSAU6gFXB7VKBZ4G+H1o3Yb0+JgPTgo/tRaxuwL8DW4AcrL/ZV0Skv9fV8Rrwa1tEpgFPAF8Jbv8qVsH3iIgH+BvWN7jk4HZX9XmuzGBsAnAnA/vbXo/198vBOk6rgo9JBnYBPwhudyFwNtZxTMRqXNSc8FEZQbSgn7hSrBdKX11AFjAh2Br9wBx7wZwfBlt6bf3EHzPGbDfGtADfAz4twUHTk3QTcK8x5oAxphn4NnC99P528J/GmDZjzBasQhHqg+Gwb4hIPZCP9eHw2WM8f99jmB1szR2+fTp4fxcwRURSjTHNxpjVoXZmjKkxxjxnjGk1xjQBP8H6sOzpT8aYvcFj/TQwP3j/p4CXjTHvG2M6sI5z4Ci5Xwt8ENyuQEQ2i8ipPXL5ozGmKbivHwLzJPjtJ+hvxpgNxph2rOLWbox5NPhB9xSwgN5+a4w5ZIypDf5eN4TI6VQgzRjzI2NMpzHmAPB7rELXn2/0PObA1qNsezyv7euAV4wxbxljuoBfANHAGcBSwAX8Orif54G1fR4fAH5gjOkIvv4G+rfdb4xpwPpA3W+MedsY0w08wz+PaRcQB8wAxBizyxhTdpTfe9TQgn7icoDaEPf/L1ZBe1NEDojItwawr0PHET8IuLFacicrO7i/nvt2YbW+Divv8f9WrELdn18YYxKNMZnGmMuNMfuP8fx9j2Fp8PGHb4db4rdjtaZ2B786XxpqZyISIyK/C3abNALvA4l9Pvz6+32y6XGcgx+e/bbajDF1xphvGWNmYx2vzcALwa/zThH5mYjsD+ZRGHxYz79Zz/72thA/9z3OfV8D2SHSmkCfD0WsVnRGiG0P+0XPYw7MPcq2x/Pa7vXaMsYEgr9DTjBW0ufDoO97oCr4YQcM+G87oGNqjHkXq3V/P1AhIg+JSPxRfpdRQwv6CQi2xHKAD/vGgq2yrxtjJgGXAV8TkQsOh/vZ5bFa8ON6/H88VgujGmgBYnrk5cT6ejvQ/ZZiFYGe++5mGAb3gt0Dl2G1co/KGLPPGHMDVvfFz4FnRSQ2xKZfx+oeWhLsnjn78NMNIKUyehxnEYnB6nY5JmNMNVYLNBvrG8eNwBXAJ4AEIO848uhP39dAaYhtDgEFfT4U44wxF5/E8x5xnK/tXq8tEZHg71CCdaxzgvcdNq73w237O5m/bajf5dfGmEXAbKzGwn+cyH5GGi3ox0FE4oOtwyeBvxhjtoXY5lIRmRJ8sTYC/uANrEI56QSe+mYRmRUsMj8Cng1+Nd8LeEXkEhFxY/WrRvV4XAWQJz2mWPbxBPBVEZkY7Pc+3OfefQI5DoiIuIP98E9g9ZOGHODs85ibRSQt2MqrD94dalpnHFZLrF6sgeMfhNimP88Cl4rIWcE+3h9xlPeHiPxcROaIiEtE4oB/BfKNMTXBPDqwWvgxWMf1ZP27WFMlk7Fa3U+F2GYt0BgcTIwOflOY07Mr6GQc52v7aeASEbkg+Nr8OtYx+Rirb9sPfCF4/K4ATjvG05/M37bv73GqiCwJ5tUCtBP69TTqaEEfmL+LSBNWC+gerCL0uX62nQq8DTRjvXAfMMasDMZ+Cnw3+HX4G8fx/I9hDSCVA17gS2DNugH+DfgDVsunBWvQ6rBngv/WiMjGEPv9Y3Df7wMFWC/sLx5HXsfjOhFpxirIL2EVu0XGmFAtzb6WAzuCj78PuL7n1/EefoXVT1uNNWj5+kCTM9bMm3/Hmp1TBtTR+1j2FYPV910PHMBqjV4ejD2K1d1QAuwM5nKyHgfeDD7XAeC/+m4Q/JC/DGtcoADrOPwB61vCYBjwa9sYswe4GfhNMI/LsAZNO40xncDVWF1p9cHtXsYq+P35FSf4tw0hHmtsoQ7r71SD9Q1r1JNjj9cppcJJRAqBO4wxb4c7l6EiImuAB40xfwp3LqOZttCVUsNORM4Ra+69S6ylKOZycq1uhTWjQSmlhtt0rH52H7Af+FSkTB0MJ+1yUUqpCKFdLkopFSHC1uWSmppq8vLywvX0Sik1Km3YsKHaGJMWKha2gp6Xl8f69evD9fRKKTUqicjB/mLa5aKUUhFCC7pSSkUILehKKRUhdB66UmrM6erqori4mPb2UCtIjAxer5fc3FzcbveAH6MFXSk15hQXFxMXF0deXh69F30cGYwx1NTUUFxczMSJEwf8OO1yUUqNOe3t7aSkpIzIYg4gIqSkpBz3Nwgt6EqpMWmkFvPDTiQ/LehKKRUhtKArpVSE0EFRNTq899P+Y+d9e/jyUGoE0xa6UkoNs3Xr1jF37lza29tpaWlh9uzZbN++/aT3qy10pdSY9p9/38HO0sZB3ees7Hh+cNnsfuOnnnoql19+Od/97ndpa2vj5ptvZs6cOSf9vFrQlVIqDL7//e9z6qmn4vV6+fWvfz0o+9SCrpQa047Wkh5KtbW1NDc309XVRXt7O7GxsSe9Ty3oamw72mAr6ICrGjJ33nknP/7xjykoKODuu+/mt7/97UnvUwu6UkoNs0cffRSXy8WNN96I3+/njDPO4N133+X8888/qf1qQVdKqWF26623cuuttwLgdDpZs2bNoOxXpy0qpVSE0IKulFIRQgu6UkpFCC3oSikVIbSgK6VUhNCCrpRSEUILulJKRQgt6EopFSG0oCul1DD73ve+x3333Xfk53vuuWdQFujSM0WVOhq9sEbke+1bUL5tcPeZeQqs+Fm/4dtvv52rr76aL3/5ywQCAZ588knWrl170k+rBV0ppYZZXl4eKSkpbNq0iYqKChYsWEBKSspJ71cLulJqbDtKS3oo3XHHHTzyyCOUl5dz2223Dco+tQ9dKaXC4KqrruL1119n3bp1XHTRRYOyT22hK6VUGHg8Hs477zwSExNxOp2Dsk8t6Eodw6oDNSHvX929l69+ctowZ6MiRSAQYPXq1TzzzDODtk/tclFKqWG2c+dOpkyZwgUXXMDUqVMHbb/aQlfqJPzyrb1HjWsLXoUya9YsDhw4MOj71Ra6UmpMMsaEO4WjOpH8tKArpcYcr9dLTU3NiC3qxhhqamrwer3H9TjtclFKjTm5ubkUFxdTVVUV7lT65fV6yc3NPa7HDKigi8hy4D7ACfzBGBNyJr6InAqsBq4zxjx7XJkopdQwcbvdTJw4MdxpDLpjdrmIiBO4H1gBzAJuEJFZ/Wz3c+CNwU5SKaXUsQ2khX4akG+MOQAgIk8CVwA7+2z3ReA54NRBzVCp4dJaC6095pxH+cCXCUDAQEVXNIc6Y2nyuxnnaaHLH8Dt1GEoNXIMpKDnAId6/FwMLOm5gYjkAFcB53OUgi4idwJ3AowfP/54c1Vq8HW3Q+Uua7W9plJ72JtMcdcSftO6nCKT1ismxftJjYvitLxkJqfFIiLDlbVSIQ2koId6lfYdGv4VcLcxxn+0F7Ux5iHgIYDFixePzOFlNXYcWAnrfg+dLRCTCpPOg4TxIEJzt/D+wTZS6zZzreM1roh6i7e9F1Hrm0acq5uiDh8feM8hv6KZV7aVkZMYzdlTU0mPP75ZCUoNpoEU9GJgXI+fc4G+TZnFwJPBYp4KXCwi3caYFwYjSaUGlb8bVv4UPvh/EJMMs66G+GwINkZK2tzcumUi+1u8XJW1jEbnGhY3vcPF7a9QLbs4kLyMdF873eNTWJKXzPbSBlYfqOXJdYc4b3o6p+QmhPkXVGPVQAr6OmCqiEwESoDrgRt7bmCMOTJcLCKPAC9rMVcjUlcbPH4dFPwD5t8MvnRweo6E9zRFceuGSbT6HTxx6n5OT25h1YEodnlXkN20hfEN6/FVVLMj7VIAHA5hbm4i0zPjeH17Oe/uqaSt28+pE5K0C0YNu2OO6BhjuoEvYM1e2QU8bYzZISJ3ichdQ52gUoPG3w3P3gYF78Plv4Er7+9VzNfXxXDt2skAPHOaVcyPEKE0fj470i/F5W9jZtWrRHU1HAlHuZxcOjebGZlxrNpfw/v7qkfsSSsqcg1oHrox5lXg1T73PdjPtp89+bSUGmSBALz0RdjzKlz8C1h4a69wUauH2zfmkerx8+jiA+RGd4XcTVNUJntSL2Rm1etctfMrPDfnfrqcMQA4HcKFszLwupxsPlSP16UzYNTw0lecGhve/gFseRzO/Q6c9vleoZZuB3dumgDAnxYV9FvMD2v0ZrM39Xwymndy6e67EdN9JCYinD0tlRmZcawuqOXd3RWD/7so1Q8t6Cry7X4FPv41LL4dzvlmr5Ax8B/bc9nb7OU384qYENM5oF3WRefx9pTvkFe/miWH/tgrJiKcPyOdNF8UX35yM4XVLf3sRanBpQVdRbamCqurJfMUWP7TIzNZDnuwII1XKxL51rQyzk5tPq5d78i4gp1pF7Pk0MNkNW7tFXM7HVw6NwunQ/iXxzbQ0tHdz16UGjxa0FXkMgZe/Ddrnvk1D4Mrqld4b0UT9+ZncElGPZ/Pqz6hp3hv0n/QFJXJ8n3fx93duyUeH+3mtzcsZG9lE//7xp4T/jWUGigt6CpyrX0I8t+GC/8L0qb3CvkDhruf24rPFeBHs0r6NtwHrNPl4/VpPyK+vYzzCn5hi581NZVbl07gz6sK2Xyo/sSeRKkB0oKuIlN9Ebz1fZh6IZx6hy382KpCNhXV8/0ZpaR4/Cf1VKXx81ib+1lmV77MhLqPbfFvXDSdjDgv33puK13+wEk9l1JHowVdRaY3vwsIXPpLW795SX0b//PGHs6elsaVWfWD8nRrxt1BnXc85xbciyPQe5ZMnNfNj6+cw+7yJn7/weBfdkypw/QCFyryFLwPO1+E8+6BhFx476e9wj/YmIfx+/hJ5vsn3NXSV8Dh5h8Tv8qVu77K/LKn2Jhzc6/4J2dlsGJOJve9vY+L52SRt+3XR9/hed8enMTUmKItdBVZ/N3w2t2QOB7O+KIt/HFNLG9XxfOlSRWMO8Z88+NVkHwWBUlnsPTQH4jprLHFf3j5bFwO4eev7x7U51XqMC3oKrJs+BNU7rQGQt3RvULGwM/3ZpHl7eRzE05sVsux/CPvq7gC7Zx58AFbLCPey+fPnsRr28vZVB8d4tFKnRwt6CpytNXDez+BvGUw83Jb+NWKBLY0xvDVKRV4nUOzzkpdTB6bsq5nduXfSWu2T1W8Y9kkUn0efr43C13qRQ02Legqcqz+P2irg4t+YhsI7QrA/+7LZLqvjWuy64Y0jTXj7qDDFcfpRb+zxXxRLr54/lRW1/n4R7VvSPNQY48OiqrI0FYHqx+AGZdC1jxb+MniFApbo/jjwgKcfQZCVx2w93efjE6Xjw3ZN3Jm0YM8/vwLVMT1vgSvP2BId7fxg+1pvHduMw5dZVcNEi3oKjKsegA6GuFc++yQdr9w3/50Tktq5rzUpkF7yqVFD/UbcwY6aXMlsPTQ73lx1i97xxzCp1MO8Nvy2fy9LJErsusHLSc1tmmXixr9Wmut7paZl0PmHFv4mZJkqjvdfG1KxaBNUzwWv8PDhuybmFT3IRlNO2zx032VjPM089sD6QS0L10NEi3oavRb/QB0NsE5d9tCXf4ADxaksTCxhSVJw7vq4Zasa4+00vtyCFyeVMS+Fi9vV8UPa14qcmlBV6NbVzusfhBmXRGydf7S5lJK2j38+8TKYWudH9bp8rEh5yYm1X0UupUeV8m46A4eOJCmM17UoNCCrka3sk1W6/zs/7CFAgHD//1jP+M9zUQ3FbLqQE3I21DanPlp2p1xLC55zBZziuHOvCo2N8SyqjZ2SPNQY4MWdDV6BfxQsgEmnWetd97HmzsryK9s5vLkg8PeOj+syxXLtsyrmFLzHvHtJbb4tTl1pHq6eKAgPQzZqUijBV2NXpU7obMZzviCLWSM4YGV+UxIiWGpryoMyf3TpqzrMCIsLH3CFvM6DbdPqObDmji2NujZo+rkaEFXo5MxULwWYtNg8gW28PqDdWwtbuDzyybhlPB2ULdEpbMn9SJmV7xEVHejLX7z+BriXH5+V5gWhuxUJNGCrkanukJoqYLcU21nhQI88lEhCdFurl6YM/y5hbAh5yY8gTbmlj9vi8W5AlyXU8vrFQmUtbvDkJ2KFFrQ1ehUvBY8sZA+yxYqrW/j9R3lXH/qOGI8I+PcuerYaRxMOI35ZU/b1ksHuHV8DQEDfz2UHIbsVKTQgq5Gn5YqqCuA7EXgsBfsx1YfxBjDLadPCENy/duQcxO+ziqmV7915L7DM21KystYGFvNnw8m8Y/8On751t4wZqpGKy3oavQp3QTihOz5tlB7l58n1hZx4axMcpNihj+3oziYeDq10ROYW/5syPjyxBKa/B5WNeuMF3VitKCr0aW7Ayq2Q/pMcNsL9oubS6hv7eKzZ+YNf27HIsLWzGvIbtpGbKd9PfbZ0XXkeFp4vT4Xo2caqRMwMjoYlRqoyh3g74TshUfuOtw9YYzh8bVFpPo8fJxfzar91klDS8OSaGg70y/lzIP3k9G8kwPJZ/eKicDyxGIerpxOWUN7mDJUo5m20NXoYYzV3eLLgLgsW7i0oZ3q5k7mjUtEwnUm0TF0uOLYnbac1NZ8nIEOW/ysuHJiHV1sOVQ//MmpUU8Luho9GoutAdHshSGnKm4vacDjdDA9Iy4MyQ3c1sxP4TR+0lr22WJeR4Bl8eXkVzVT02wv+EodjRZ0NXqUbgJnVMipiu1dfvZVNjMjMw63c2S/rCt9M2jypJPZvJNQq3KdH19KwMBzG4vDkJ0azUb2K1+pwzpboGq3tWaL037yza6yRvwBw5ychDAkd/zKfbOI7m4gvqPUFhsX1UpWgpcn1h7SwVF1XLSgq9GhfBuYAGQvsIWMMWwvaSQz3ktaXFQYkjt+NTET6XJEkdG8O2T8lJwECqpbWH2gdpgzU6OZFnQ18hkD5VshIRdiUmzh0vp2als7mZMzei4UYcRFdcxUktsKcfntM1qmpvuI97p4Ym1RGLJTo5UWdDXyFa2CtlrItF/8GWB7qTUYOm2ED4b2VRk7HQcBUlvtg6Mup4OrF+by+vZyals6w5CdGo20oKuRb+Oj1mBo2nRbqL7TOWoGQ/tq9STT5Ekjo2VPyMHRG04bT6c/wHMbdHBUDcyATiwSkeXAfYAT+IMx5md94lcAPwYCQDfwFWPMh4Ocqxrt3vtp/7Hzvh36/vYG2PGCdWao02MLv1iWOKoGQ/uqjJ3O5LoP8XVW0RzV+5T/6ZlxLJqQxJPrirhj2cQRO7dejRzHbNKIiBO4H1gBzAJuEJG+88beAeYZY+YDtwF/GOQ81Vi17Vnobuu3u+WZ0iTSfFGjZjC0r5qYyfjFRXrLnpDxaxflsr+qhU16opEagIF8Rz0NyDfGHDDGdAJPAlf03MAY02z+Ob8qFtC5VmpwbHoMMuZAXKYttKvJy/bGGGZlj57B0L78Dg810ZNIbd0fclndS+Zm4XU7eGa9druoYxtIQc8BDvX4uTh4Xy8icpWI7AZewWql24jInSKyXkTWV1WF97JgahQo32adTLTglpBnhj5TkoRHAkzPHF2DoX1V+qbjNF2ktB6wxeK8bi6ek8XLW0pp6/SHITs1mgykoIfquLO1wI0xfzPGzACuxOpPtz/ImIeMMYuNMYvT0vRyW+oYNv3V6jef+2lbqDMgvFCaxCfSG4l2O8OQ3OBp8mTQ5kogvTX0GuifWpxLU0c3b+woH+bM1GgzkIJeDIzr8XMuYD+9LcgY8z4wWURSTzI3NZb5u2DbMzBtOcTYr+LzblUctV0urs2pC0Nyg0yEqtipxHeUh7zm6NKJKeQmRfPMhkMhHqzUPw1klss6YKqITARKgOuBG3tuICJTgP3GGCMiCwEPUDPYyaoI1ncGTPU+aK22WughZsc8U5JMelQXy1Ka2BwYphyP09Kihwa8bVXMVMY1rCetZR/FCYt6xRwO4VOLcrnvnX0U17WOuAt3qJHjmC10Y0w38AXgDWAX8LQxZoeI3CUidwU3uwbYLiKbsWbEXGd0EQp1Miq2WRewSJpoC1V2uFhZHcc12XW4RtfU8351unw0ROWQ1rI35Jz0axbmYgw8t6EkDNmp0WJA89CNMa8Cr/a578Ee//858PPBTU2NWV1tUJMfvGaovX/8xbJE/Ea4JjsCult6qIqdytTalcR32PvKxyXHcMbkFJ7bWMyXLpiic9JVSBHSvlERpXKntRBX5pyQ4edLk5iX0MoUX2StF14bnYdf3KT1Mzh6zcJcimpb2XAwsj7I1ODRS9Cpkad8G8SmW1cm6mNno5ddTdH8aGbkdT0EHG5qYiaS0lqAy9925NJ6h3V2B3A5hB++tIMLZmbw1U9OC1OmaqTSFroaWVqqoLncWvc8hL+VJeGWAJdm1g9vXsOkMmYaTtPF1Jr3bDGPy8GUdB97K5vp9o/QkWAVVtpCVyNLxXZAQl6VqDsAL5Qmcl5aE8meyDzJpikqk3ZnHDMrX2FX+sW2+IzMOHaXN1FQ3WJrwfelLfixR1voauQwAajYCcmTwBNrC39Y46Oq083VETYY2osI1bFTGN+wjtiOSlt4XHIMsVFOdpU3hSE5NdJpQVcjR30RdDZZa7eE8HxpEonubs5Li+xiVhUzFcEwo/oNW8whwozMeA7WtNDa2R2G7NRIpgVdjRwVO6wTiVKm2EJN3Q7eqEzgssx6ohyRfYpDuzuBMt8cZla+FjI+IzOOgIG9Fc3DnJka6bSgq5HB3wXVe6yLWIS4CPRrFQl0BBxclV0//LmFwa70FaS17iO1xX41o9TgcsG7yuzLBKixTQdF1chQsw/8nZAeurvlhdJE8mI6WJDQaosdzyn2o8We1As5p+BeZla9xgexU23xGZlxfLCvmtqWTpJj7Rf+UGOTttDVyFCxA6LiIXG8LVTW7mZVrY8rs+pCraIbkdrdiRQmnsGMqtcRY5/RMz0jDgH26OCo6kELugq/zhaoPWBNVQxRsV8qS8QgXJlVP/y5hdGu9IvxdVYxrmG9LRYb5WJccgy7yxvRZZPUYVrQVfhV7gJMv7Nb/laayIKEFvJiO4c3rzA7kLyMDmfsUQdHG9u7KWtoH+bM1EilBV2FX+UO6zT/WPsS+ruavOxujh4zg6E9+R1R7E35BFNq3sPltxftyWk+XA5ht3a7qCAt6Cq8WmuhqQzSZ4cMv1CaiEtMxJ7qfyy701fgCbQyqfZ9W8zjcjA5zcfeiib8Ae12UVrQVbhV7rD+TZ9pCwWMtVTuuamRe6r/sRTHL6DJk87Mqv67XTq6AxTWtAxzZmok0oKuwscYa3ZLYh5E2S/0vLo2lvIOD1dG8qn+xyIOdqctZ0L9KqK77MdhfHIM0W6ndrsoQAu6CqfGUmivh4x+ulvKkvA5/XwibWyfQLMrbQVO42da9Vu2mMMhTM+Io6C6hY6usfktRv2TFnQVPpU7wOGCVPuqgO1+4bWKBJZnNOB1ju3+4ZrYKVTFTGFG1esh49Oz4vAHDPlVuhTAWKcFXYVHwA9VuyBlKriibOF3q+Jp6nZy5Ric3RLK7rQVZDdtI6Gt2BbLiIsiMdqt3S5KC7oKk7oC69qh/XS3/K0skfSoLk5P1lYnwO60CzEIM6rtrXQRYUZmHMV1bTS1d4UhOzVS6FouKjwqdoArGpIm2kL1nU5WVsXxmfE1OIMnjq46UDPMCY4szVGZFMcvZGbla6zJvd12Ru30zDhWF9Syt6KZRROSwpSlCjdtoavh191hLcaVPhMcTlv4lYoEuoxjbM9uCWFX+gqS2ovIaN5piyXGeMiM97K7fGwPII91WtDV8KveC4Hufk8merEskamx7cyO01Pae8pPOZ9u8Rx1Tnp1cyfVzR3DnJkaKbSgq+FXuQO8iRCfbQsdanOzts7HldljZ2XFgepwxXEgeRnTqt9CjP1qRVMzfIigg6NjmBZ0Nbw6mqDu4FFXVgS4fIytrDhQu9OWE9tVy/j6tbZYjMfFhOQY9pQ36QqMY5QOiqrhdWRlxd7dLasO1GAM/PXgZGZ46ykuK8c+QU8VJp1BtyOKZYW/Iatpuy0ecKfzm47ZlNS3hSE7FW7aQlfDq3IHxGVBTIotVNjho7QrlrPiy8OQ2Ojgd3iojp5IclshjoB9iuKi2GrcTl2BcazSgq6GT0s1NFf0Oxj6QVMmLgmwxFc1zImNLtWxU3GabpLbCm2xKEeAKek+9lU0065LAYw52uWihsUv39rLtXvWk4OwoS2Drj7zyv1G+LgpnYUx1fic9gE/9U9NngzanT7SWvOpDnm90Xh2lTXxzq5KLpmbFYYMVbhoC10NDxMgrTWfem8OXc4YW3h7axIN/ijOiq8IQ3KjjAjVMVNIaC/B7bdfNDs3KZrYKCd/21QShuRUOGlBV8Mit3ETUf5mqmPsLUqAD5syiHV0MT9mbJ8ROlBVsVMRDKmt+baYQ6wVGFfuqaS2ZWxdtm+s04KuhsXMylfxi5va6DxbrD3gZF1zGkt9lbgdOt1uINrdiTR70khrsRd0sLpdugOGV7aWDnNmKpy0oKuh19XG1Jp3qInOI+CwD9usa06lwzi1u+U4VcVMIbarhujOWlss1edhekacdruMMVrQ1dDb8ypR/haqQgzgAXzYlEmaq41p3oZhTmx0q46ZTAAhrXWfLSYiXLkgh41F9RRW6+Xpxgot6GrobXmKJk86jVH2GRe13R62tSaxLL4Ch57qf1y6ndHUe8eR2rofTMAWv3JBNiLwvLbSx4wBFXQRWS4ie0QkX0S+FSJ+k4hsDd4+FpF5g5+qGpWaqyD/bXanLQexv9w+asrAIJwVpycTnYjq2ClE+VuI7yizxbISojlrSirPbywmENCxibHgmAVdRJzA/cAKYBZwg4jM6rNZAXCOMWYu8GPgocFOVI1S258D42dX2gpbyBh4vzGTad4Gsjx6qvqJqPNOoFvcpLXYu10Arl6YQ3FdG+sK7f3sKvIMpIV+GpBvjDlgjOkEngSu6LmBMeZjY8zhxatXA7mDm6YatbY8AZmnUBM7xRYq7PBR3OljmZ7qf8ICDhc1MZNI6WcpgItmZxLjcfL8Ru12GQsGUtBzgEM9fi4O3tef24GQCzaLyJ0isl5E1ldV6endEa9iJ5Rthvk3hQwfPtV/qa9yePOKMFWxU3GarpBLAcR4XKyYk8Ur28p0KYAxYCAFPdRQVcgOORE5D6ug3x0qbox5yBiz2BizOC0tbeBZqtFpy+PgcMEp19pC3Ub4qCmDRbF6qv/JavJk0u6MI71lb8j4NYtyaO7o5s2dOi000g2koBcD43r8nAvYzlYQkbnAH4ArjDF6ut9Y5++GrU/D1IsgNtUW3tqSTKPfwzIdDD15IlTFTiO+oxRPt32VxaUTU8hO8PL8Rl2QONINpKCvA6aKyEQR8QDXAy/13EBExgPPA7cYY0I3E9TYsv9da2XF+TeEDL/flEmcs5N5sTpYNxispQAIOTjqcAhXLczh/b1VVDbqZf0i2TELujGmG/gC8AawC3jaGLNDRO4SkbuCm30fSAEeEJHNIrJ+yDJWo8OWxyE62Wqh99Hkd7GhJZWz4ipwiU6nGwwdrjgaorKsk4xCXK3omoW5BIzOSY90A5qHbox51RgzzRgz2Rjzk+B9DxpjHgz+/w5jTJIxZn7wtngok1YjXFsd7H7V6jt3eWzhj5oy6DYOzom3z51WJ64qdhrR3Y1kNW21xSal+Vg8IYln1h/Sy9NFMD1TVA2+7c+DvwPm3xgyvLIxi4lRTUyI0lPSB1NN9ET84mJ25csh49cuzmV/VQubDtUPb2Jq2GhBV4Nv8+PWRaCz7CcMF3b4ONgRp63zIRBwuKmJnsi06rdw+e195ZfMzSba7eSZ9YdCPFpFAi3oanBV7oKS9dbcc7HPeP1HozX3/Mw4nUI3FKpipxHlb2Fqzbu2mC/KxYpTMvn7ljLaOnVOeiTSS9CpwbXxMXC4Yd71tlC3ET5s1LnnA7G06MRWz2iMyqLem8vsiheBr9ji1y4ax/MbS/jKk5uYkRXf736++slpJ/T8Kry0ha4GT3cnbH0Spq8IOfd8Q3MqzQEP52p3y9ARYXv65Yxr3Ag1+23hJROTGZ8cw46yxjAkp4aaFnQ1ePa8Cq01sPDWkOGVjVkkOTuYG6Nzz4fSzvRLCOCATX+xxRwO4VOLcimua6Ohzb72ixrdtMtFDZqCtx4k1ZPOwwW5mMLe55c1tXexpTWZK5IO6rrnQ6wlKp3CpDOYtPlxOO8ecPZ+m39qUS73vrWXnaWNnD45JUxZqqGgLXQ1OBqKyatfzY70SzHitIV3lDZiEM5L0O6W4bA943JoLof8t22x7MRo8lJi2FHWoOukRxgt6GpwbH4cwbAz4zJbKGAMO0obOSWmlnS3nno+HAqSlkFsGmx6LGR8Tk4CLR1+Cmv0XIBIol0u6uQFArDpMRqisplZ+aotvLklmeaOeZyfqVegHy4Bhwvm3QCrH4CmCojL6BXPS4kl1uNke2kjk9J8YcpSDTZtoauTt/9dqC+iwjcjZPjdhmzinZ0s9lUPc2Jj3MLPQKA7ZCvd6RBmZcdTWN1CU7sOjkYKbaGrk7f+jxCbRm10ni1U3+1hY0sKK5KKdSGuYfbLTQGuSTiVhI/+wJ/aL7GNbczOTmBdYR07SxtZMkkHRyOBttDVyWkogb2vwYKbQw6G/qMxEz8Ozo/X7pZw2Jp5DQkdZUyoW22LJUS7GZ8cw/bSRgK6YFdE0IKuTs7GR63lWhd+xhYKGKu7ZVZ0nV4EOkz2J59DizuZueXPhYzPyY6nuaObgzWtw5yZGgpa0NWJ83fDxj/DlAsgeaItvLk1hcruaD6ZoGtwh0vA4WJ7xhVMrPuIuA771aEmpfmI8TjZWlw//MmpQacFXZ24va9DUxks+lzI8Fv1OSQ6O3QwNMy2Z1yJYJhT/oIt5nQIc7ITKKxp1TNHI4AWdHXiNvwJ4rJg2nJbqKLLy5bWZC5IKNXB0DBr9GZTkHQGcypexBGwL4p2Sk4CIrCtuCEM2anBpAVdnZia/dZZiAs/Yzu1HODthhwEw/kJOhg6EmzLvBpfVzWTa1faYj6vi8mpPnaUNtDtDwx/cmrQaEFXJ2bt761lchfbu1s6Aw5WNmRxqq+aZFdnGJJTfRUknUm9N4f5ZU+HjM/NTaC9O8DeiuZhzkwNJi3o6vh1NFkr+c2+CuIybeFVzek0B9w6GDqCGHGyJfNachs3kda8xxbPTYomOdbD1pL64U9ODRot6Or4bX4COptgyV22kDHwZn0OOZ4WZkXXD39uql87Mi6ny+FlftlTtpiIMDcngYrGDsobdL2d0UoLujo+gQCs/R3kLIbcRbbwnvYEDnTEc1FCcagr0Kkw6nDFsTP9EmZUvUF0V50tPjMrHo/TwWa9iPSopQVdHZ/970JNfsjWOcBrdbn4HF2cHW+f86zCb3PWp3GZzpBTGD0uB7Nz4tlX2URZg54INhppQVfHZ82D4MuAWVfYQhVdXta1pHFBQilRDp0tMRLVxkziYOIS5pU/G3IK4/zcRIyBP398MAzZqZOlBV0NXNUeyH8LFt8GLo8t/Hp9Lg4MFyYWhyE5NVCbsq4jrrOSqTXv2GLx0W6mpPt4fM1BWjr0Qt6jjRZ0NXAf/wZcXjj1Dluosb2LlQ1ZnB5XqVMVR7iCpDOpic5jUclj1ih2HwvGJ9LY3s2zG/SDebTR5XNVb+/9NPT9HU2w+XFr3nlsqi381NpDtBsXFyceGuIE1bEsLXqo39jq8XeCONiQczMX5v8X4xvWUpS4pNc2WQnRLByfyJ8+KuDmpRNw6kVgRw1toauBKdkAJgCn/7st1OUP8MjHhcyIrmeiV09MGQ12p62gxZ3C4pLQl6i7Y9kkCmtaeXtXxTBnpk6GFnR1bN0dULoJ0qZD8iRb+MXNpZTUt3F5kg6kjRZ+h4dN2dczoX5NyBONLpyVwbjkaB5YuR+ja6WPGlrQ1bGVbQZ/B+QusYX8AcMDK/OZmRXP/Jja4c9NnbCtmdfQ6Yix+tL7cDkd3HXOZLYcqmfV/powZKdOhBZ0dXQBP5Ssh4TxEJ9lC7+5o5wDVS3827mT9USiUabDFcfWzKuYXv02ce1ltvg1C3NJj4vi/pX5YchOnQgdFFVHV7HNGhCdtsIWMsbwwMr95KXEcPEpWaxdH4b81HHpO2AqJgAYLt19NwXJZ1mDpkFet5PPL5vET17dxaaiOhaMTxrmbNXx0ha66l/AD0WrrDXPk+xXJPpgXzXbShr413Mn60yIUarT5aMydjrpLXvwdNsHtG9cMp6EaDf3v7c/DNmp46UtdNW/iu3Q3gBTPkmo/pT738snK8HLVQtyw5CcGiwl8fNIb9lDTtOWI/f98q29R/4/IzOOt3dVcM/ftpHqiwLgq5+cNux5qmPTFroK7XDr3JcJyZNt4Y/zq1lTUMvnl03C49KX0WjW6YqjKnYa6c278XXYpynOH5eI2ymsLdBB75FuQO9EEVkuIntEJF9EvhUiPkNEVolIh4h8Y/DTVMOucge018OEM22tc2MM/++tvWTGe7lxyfjw5KcGVUn8fMCwuORRW8zrdrJgXBL7KpupauoY9tzUwB2zoIuIE7gfWAHMAm4QkVl9NqsFvgT8YtAzVMPPBIKt8wxImWILr9xbxYaDdXzxgil43c4wJKgGW4crjurYqZxS/gKxHVW2+MLxiUS5HKw6oFMYR7KBtNBPA/KNMQeMMZ3Ak0CvpfaMMZXGmHWAXjY8EpRvg7a6flvn9765l3HJ0Vy7aFyYElRDoTh+AYKfJcUP22JRbicLJyRRUN2iS+uOYAMZFM0Bei7QUQzYzzAZABG5E7gTYPx4/ao+Ivm7oPBDiMuGlKm28Js7K9hW0sAvrp2nfecRpsMVz7aMq5hb/jc2Zt9AffSEXvH5uYlsLqpn1YGaXoOmfemAafgM5B0Zaj7aCZ0LbIx5yBiz2BizOC0t7UR2oYZayQbr8nKTzrW1zv0G7n1zL5NSY7lyfnZ48lNDas242+l2eDjz4IO2mMflYHFeEodq2yiuaw1DdupYBlLQi4Ge361zgdKhSUeFVWut1XeePBkS7d+gnitJYk9FE1/95DRcTm2dR6JWTyobc25iWs3bZDTtsMXn5iTgi3LxYX61rvEyAg2ky2UdMFVEJgIlwPXAjUOalQqPD++11myZdK4t1NLt4L93p5MZ72VfRVPIr9xLhyFFNfQ2ZN/E3PLnOOvgb3lu9gO9vqm5nA5On5zCWzsr2FPRxIzM+DBmqvo6ZjPLGNMNfAF4A9gFPG2M2SEid4nIXQAikikixcDXgO+KSLGI6F96NKkvgjUPQcYpEGvvDvtdQRr1/ijOnpaK6KItEa3T5WNN7u2Mb1hPXv3HtvjMzDjS46L4KL+GLr9eanAkGdCZosaYV4FX+9z3YI//l2N1xajR6s3vgjggb5ktVNrm5qHCNE73VXBVw3vQEIb81LDamnk188ue4pyCX1KUcBoBh/tITEQ4e2oaz24sZmNRHUsmpoQxU9WTdoQqOLASdr4Iy74OXvsXq1/kZxIAbkg9MOypqfAIONysnPh1ktsOsqD0SVs8JymayWmxrC+so1mvPTpiaEEf6/xd8Oo3ISkPzviiLbyhLobnS5O4bUI1ae724c9PhU1h8pnsT1rG0kN/CHmy0VlTUjHGWgZCjQxa0Me6Nb+D6j2w/Ofg9vYKdQXgOztzyPJ28sVJlWFKUIXTPyZ+DYfpZtnBX9tiiTEeFoxPZFd5k05jHCG0oI9lTeWw8mcw9UKYvtwWfvhgGnuao/nhjFJiXTr4NRY1ROeyPudmZla9Tk7DRlv8tInJJES7eWd3Jd0BfY2Emy6fG4ne++nR4+d9G4yBV74O/k5Y/jPbJofa3PwqP4NPpDVwUUbjECWqRoN1uZ9jZtVrfGL/f/OX+X/F74g6EnM7HZw3PY0XNpdSvuVtPpVSCO+FGCQ979vDl/AYpi30sWrnC7D7ZeuNltJ7eVxj4Ie7cnCI4T9n6jlkY12308vbk79DcttBlhb93hafkBLLtAwfL9ZNoLQzOgwZqsO0oI9FLTXwyjcgaz6cbh8IfakskXeq4vnqlApyonW9NQVFSUvZnn45i0v+QnrzLlv87KlpRImf31fMwK8nkIaNFvSx6PVvWVciuvIBcPbudStrd/O9XdksTGzhc+N19oL6p/cnfoVWdxIX7vsxjkDvD/rYKBe3pOWzuz2RhwtTw5Sh0j70saZ6L+x4Hs79NmTM7hW69809vLs2gw4/3JKwjXWFukyq+qcOVxzvTP4WV+z+BksOPcyqCXf1ip8dV86G5lT+Z18m8e3FTIhqORJb3f3PpSJ0Ncahoy30saS9Efa8Cplz4ayv2cJbihvY3pbMzWn7yfRoMVd2B1LOYWfaJZxW/CdyGjb0ionAHel7iHV0cX/5LLoCukTEcNOCPlaYAOx+yfr32kfA5ekV3lvRxIf51cyPqeGCeB0IVf17d9J/0ODNZcXe7+Ptqu8Vi3d1cWfGHg51+ni6ZlJ4EhzDtKCPFQc/goZia855n1ktzR3d3PWXDXicDu7M2N13GXSleulyxfLq9J8Q3VXHhft+ZE2L6mFhbA2fSCjh5frxrG/W/vThpAV9LKg7CAc/how51q0HYwx3P7uVgzWtXHxKJkmuzjAlqUaTSt8MPsj7EpPrPmBB6RO2+C2p+UyMauSBipmU61TGYaMFPdK11VtzzmOSrdZ5H3/8qJBXtpXxzYumk5sUM+zpqdFrc9Z15Cefy9mFv2Zc/dpeMY8jwFeztuPEcG/ZHF1md5hoQY9k3R2w4znAwOxrwNm733zV/hp++uouLpyVwZ1na3+nOk4ivD71h9TG5HHpnm/j7eq9rnKau4MvZO6kuDOWd3ZV6hWOhoFOW4xUxsDuv0NLNcy9zmqh97Cvool/eWw9eamx/OLT8/SiFQqApUUP9RtbPf5O231drlhenPn/uHHLZ5he/SbbM67A7/hnw2FebC2fTingqYpJJES7OX2yrp0+lLSgRyJj4MC7UJMPUz5pLY0btOpADe+07OSp9YfwBwzLpqTy8AcF4ctVjXqN3hz+PuPnfGr7vzKt+m12p12EEeeR+BVJB9nmmc/awlrivFpyhpJ2uUSioo+heB3kLILshb1CbQEnL24ppb3LzxXzsomPdvezE6UGriRhEQeSlpHYUcLUmves6bFBInD+jHQmpMTw7p5K3tutSzEPFS3okWbN76DwA+ss0Mmf6HWB35ZuB/9TOpfq5g4unpNFerz3KDtS6vhU+aZTmLiUlLYCJtd90Gs6o9MhXDwni1RfFP/61w18vF+XlRgKWtAjycbH4LVvQspUmH6JrZh/bmMee9viWT47k7zU2DAmqiJVWdwpHIpfSHrLXvLqV/Uq6h6XgyvmZTMuKYbbHlmnVzoaAlrQI8Xq/4OXvgCTzoNZV1gXfA46XMw31MfyhcxdTMuIC2OiKtIVxy+k1HcKWc07mFz3fq/ul9goF0/cuZTxyTHc9ud1fKRFfVDpCMVINZCLVIDVAlr5U/jHz2HmZXDNw/DBvUc2q+xwccfGPLY3RnPf3CJS2ypZPYRpK4UIBxOX4He4Gde4EWegk3Xjbjsy+yXVF8UTn1/Kjb9fw22PrONX181nRfUj/e/vZC6OMdD3UYTQFvpo5u+CV75mFfP5N8GnHgHXP68mk98cxdVrprCvxcvvFhzksqyG/vel1GASoThhEQWJS0lpK+TKnV8mqsc89RRfFE/cuZTZ2fH82+Mbeaggte8KAuoEaEEfrZqr4NErYP0f4cwvw+W/7bW2+Uc1Pq5eM5l2v/DUqfv5ZLpeRk4Nv/K4U8hPPoecxs3cuPUzpLTkH4klx3p4/PNLuXhOFv+9N5vv7sqhU1doPCna5TLKrDpQQ2xnNXmrzyKmq463pv4nu7kY3rHeKAFjKN0cx/O1eWR7Wrk7eystNe2sqglz4mrMqoqdxvqcW7hs991cv/U2mPYQzLocAK/byW9uWMC45i08WJDOjkYvv5lXxDi9UtYJ0Rb6aGICZDduZk7Fi2AMT53ye3anX3wk3NLRzQubS3iudiJnxlXwX+M2kOZuD2PCSlnK4ufx13mPUhMzCZ6+Bf7+ZehoBsDhEL41rZwH5h1kf4uXSz6eyhsV8WHOeHTSgj5atNXD5seZ0LCOuugJ/HX+X6j0zQSsFRN3lTXyl9UHKa1v58703fxbxi68Dn94c1aqh5aodJ455SGri3DDn+HBs6BozZH4xZkNvHz6PibEdPIvm/P42rZc6jqdR9mj6ku7XMLpWCPwAIFuOLTWOvtTnOxLPpfqmCm0uxMBaGjr4t3dlRTVtpKV4OWCGemcV/vuUXd5tPU6lBoKvV5zk1Jg3g2w5xX444XWxconng3uGCbEdPLskv38Zn86Dxak84/qeL6fUcLl87J1vaEB0II+UhkDNftg/7vQXg+p02Hy+VSXdgPQ3uVnXWEtWw414HDAudPSmJubYL3oa8ObulLHlDgeFt0GhR9CyXqo2g15yyBrPlEOJ9+YWsGlmQ3cvSOXLz+5mUdXHeQ7F89g0YTkY+97DNOCPtIYA7X74eCH0FQOMSkw9/ojC2y1+ut5qyGHvxUU0tEdYGZWHKdPSiHOq2uyqFHGFQVTLoCsuZD/NuS/BcVrYfwZkDGHGXHtPL8kn6dib+KXb+/lmv9bxUWzM/jSBVOZnZ0Q7uxHJC3oI0Wg22qlFK+H5nLwJsC0FdYVhhxOKtpd/PFgKo8WzaIt4GJCipczJ6eSFhd17H0rNYKsOtB3ypUDfJ8k0VXMuMYN+Pa+Rvv+DyiPm01l7HQq8tq5dlEum4rq+Si/hjd2VHDWlFT+5ZxJnDUlVbtietCCHm4t1VCxHcq3QlcrRCcfKeRdOHmvOp6ni5N5rzoOY2CJr5LLkoqomHZjuDNXavCIUB89jnpvLonth8hp3Exe/RrGNawnuruBnemX4c47hQdvWcTja4r400cF3PLwWiamxnLt4lw+tTBXF5sDJFxXEVm8eLFZv359WJ47rIyB6n2w51VYfT80VwICKVMgZxGd8Xk8vNOwrjmN9S2pNPo9JDo7ODu+nPPjS8nw6DRENTbEdNaQ2byD1NZ8nMZPu9OHN3sOpE6jIzaLl8uTeKokmbV1PpwOYemkZNxOB5PTfPiirLZq3wkAp0/qc4GNUXjqv4hsMMYsDhXTFvpwaKqwZqkUfgj73oL6g9b9cdkEJn+CAzHz+KApk48K4lhdG0uz34lXulkQW8MZcZXMj63BJXpetBpbWj0pHEg+m8LEpSS3HSS1dT/eQ2vg0Gqi3NFckzyJa8ZPoOiae3gq38lr2ys4UN3Cyj1VpMVFMT4phujuZKZHN4yZKbzaQh9MxlhdKFW7oGyLdSvZaA1yAsYdQ2vOmRxIPJOPHQv5cE8ZmxtiaOq25tpOiO7gzJRmcvwlzImpw+PQC+sq1dPp42KgrgBq9kPdAehqswK+TEzuIt6pz2Jj13hWN2eypTEGv3EgGHI9LUzxNnJBboAZvnam+dqJdwciroU+oIIuIsuB+wAn8AdjzM/6xCUYvxhoBT5rjNl4tH2OyoJujDWFsKkcmsqgsRTqDlot7toCTPVepL3+yOY1zjTynZPZxAw+7p7GmvZxdASs4i0CEzxNTPE2MsXbyMzoetL1rE6ljqpXl4kx0FptXcylaBWUbsbU5CNYNa3T4aVB4ikhg/3+DLZ351AQSKfSJFFhkvBERTE+O5u8lFjGp8SQleAlKyGazAQvqT4PvijXiBxwPakuFxFxAvcDnwSKgXUi8pIxZmePzVYAU4O3JcD/Bf8dfMZY6yv3vQX8mOAt4O/GGD/+7m7wdxEI3kxXB4HuTgLdHQS62jHdHQQ62zCdrQQ6WzGdrdDRBJ3NSEcTjo4GHJ2NuDrq8XTU4e2qw2m6e6Xjx0G1pFBk0tnTvYj9Jpt8k8OOQB61xON0CInRbhJj3cxJ9ZDis27JMR7OKvnDkBwipcYEEYhNg1Nvt27A/a9tIq11HymtB0huLWRS3QdM7z7IXHZwjav3N94AQlN5PDVlcVQFfDSaGEqJYbeJpZUo2iUGR5QPR1QMzqgY3FGxuDzRuKO8eKK81r8eLx5PFG6PB4/bg9vjxu1y43a5cLrduJ0unC6X9bPTeeQWE+UmJmrwpxoPpA/9NCDfGHMAQESeBK4Aehb0K4BHjdXcXy0iiSKSZYwpG+yEN77+CAvXfCVkTIK3w+sZnMjh8huhBS/NRFNrYmkklkYTR43Jpo446iSBRlcaLVGptHoz6I7NJi42msQYN6m+KKbERVG1v4aLolz4vC6i3c4R+SmvVCTqcsVSGj+f0vj5AHQU+ayAMXj8LSzKdFlryHQ24ehsJSFtOgmtNUxoqaG7tR7TVomjsxFnd6vVeOvGurUMbp6rsm7m9H+5f3B3ysAKeg5wqMfPxdhb36G2yQF6FXQRuRO4M/hjs4jsOa5s/ykVGMJLnZzwuuFDnNdJGam5aV7HR/Oy+c7RgsfI6/VBzmWgHkjlrgdO9HhN6C8wkIIeqnnZt+N9INtgjHkIOOmFRERkfX99SOE0UvOCkZub5nV8NK/jM9byGshqi8XAuB4/5wKlJ7CNUkqpITSQgr4OmCoiE0XEA1wPvNRnm5eAW8WyFGgYiv5zpZRS/Ttml4sxpltEvgC8gTVt8Y/GmB0iclcw/iDwKtaUxXysaYufG7qUgUHothkiIzUvGLm5aV7HR/M6PmMqr7CdWKSUUmpw6RWLlFIqQmhBV0qpCDHqC7qIfENEjIikhjsXABH5sYhsFZHNIvKmiGSHOycAEflfEdkdzO1vIpIY7pwARORaEdkhIgERCfv0MhFZLiJ7RCRfRL4V7nwOE5E/ikiliGwPdy6Hicg4EXlPRHYF/4ZfDndOACLiFZG1IrIlmNd/hjunnkTEKSKbROTlwd73qC7oIjIOa0mConDn0sP/GmPmGmPmAy8D3w9zPoe9BcwxxswF9gIjZVWi7cDVwPvhTqTHMhcrgFnADSIyK7xZHfEIsDzcSfTRDXzdGDMTWAr8+wg5Xh3A+caYecB8YHlw9t1I8WVg11DseFQXdOCXwDcJcRJTuBhjGnv8GMsIyc0Y86YxRxaiWY11rkDYGWN2GWNO9IzhwXZkmQtjTCdweJmLsDPGvM8Iu1qsMabs8CJ8xpgmrCKVE96swFiagz+6g7cR8T4UkVzgEmBIFnIatQVdRC4HSowxW8KdS18i8hMROQTcxMhpofd0G/BauJMYgfpbwkIdg4jkAQuANWFOBTjSrbEZqATeMsaMiLyAX2E1QodkbewRfYELEXkbyAwRugdrAYcLhzcjy9HyMsa8aIy5B7hHRL4NfAH4wUjIK7jNPVhflf86HDkNNK8RYkBLWKjeRMQHPAd8pc831LAxxviB+cGxor+JyBxjTFjHH0TkUqDSGLNBRM4diucY0QXdGPOJUPeLyCnARGBLcCXDXGCjiJxmjCkPV14hPA68wjAV9GPlJSKfAS4FLjDDeALCcRyvcNMlLI6TiLixivlfjTHPhzufvowx9SKyEmv8IdwDymcCl4vIxYAXiBeRvxhjbh6sJxiVXS7GmG3GmHRjTJ4xJg/rjbhwOIr5sYjI1B4/Xg7sDlcuPQUvUnI3cLkxpjXc+YxQA1nmQgUFL2zzMLDLGHNvuPM5TETSDs/iEpFo4BOMgPehMebbxpjcYM26Hnh3MIs5jNKCPsL9TES2i8hWrC6hETGVC/gtEAe8FZxS+WC4EwIQkatEpBg4HXhFRN4IVy7BQePDy1zsAp42xuwIVz49icgTwCpguogUi8jt4c4Jq8V5C3B+8DW1Odj6DLcs4L3ge3AdVh/6oE8RHIn01H+llIoQ2kJXSqkIoQVdKaUihBZ0pZSKEFrQlVIqQmhBV0qpCKEFXSmlIoQWdKWUihD/H52ICZ0jDBEgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Indeed, the distribution $X$ is stochastically greater the distribution $y$\n",
    "X = stats.norm(loc=+0.2)\n",
    "Y = stats.norm(loc=0)\n",
    "x = X.rvs(size=801)\n",
    "y = Y.rvs(size=399)\n",
    "\n",
    "grid = np.linspace(-4, 4, 100)\n",
    "plt.plot(grid, X.pdf(grid), 'C0')\n",
    "plt.plot(grid, Y.pdf(grid), 'C1')\n",
    "plt.hist(x, density=True, color='C0', bins=30, alpha=0.5)\n",
    "plt.hist(y, density=True, color='C1', bins=30, alpha=0.5)\n",
    "plt.title('Distribution PDFs and Sample Histograms')\n",
    "plt.legend(['x', 'y'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d10f5fe3-43f8-46e9-8347-01b926170eaa",
   "metadata": {},
   "source": [
    "Our first difficulty is determining the correct value of `alternative` to pass into `ks_2samp`. From its [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ks_2samp.html), we see that the alternatives are expressed not in terms of the values of the samples or the location of the underlying distributions, but in terms of _cumulative density functions_ of the underlying distributions. \n",
    "\n",
    "> - `two-sided`: The null hypothesis is that the two distributions are identical, $F(x)=G(x)$ for all $x$; the alternative is that they are not identical.\n",
    "> - `less`: The null hypothesis is that $F(x) >= G(x)$ for all $x$; the alternative is that $F(x) < G(x)$ for at least one $x$.\n",
    "> - `greater`: The null hypothesis is that $F(x) <= G(x)$ for all $x$; the alternative is that $F(x) > G(x)$ for at least one $x$.\n",
    "\n",
    "Note that if a distribution $X$ tends to be greater than $Y$, we find that the cumulative distribution function of $X$ lies _below_ the cumulative distribution function of $Y$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5ebe1d2f-f067-46f1-a415-986b19b933be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwa0lEQVR4nO3dd3xUZdbA8d/JpJGEmoQaekeaCojYwIp9dVcXGxZc1lVsq6vuWnfdte+ua2VZe0VRfAVFsNBUpPdODCSETkggIX3mvH/cAZNMIAEmuZPJ+fqZD3Pn3LlzGGcOzzz3uc8jqooxxpi6L8LtBIwxxgSHFXRjjAkTVtCNMSZMWEE3xpgwYQXdGGPChBV0Y4wJE1bQTa0RkbEi8nCQjtVORPJExOPfnikiNwfj2P7jfSUi1wfreMbUBivoJihEZJOIFIhIrojkiMgcEblFRA5+xlT1FlV9vJrHOvtw+6hqhqomqKo3CLk/JiLvVTj++ar69rEe+xCvN0hEpvjfpz0iMl9EbvTHhoqIz/+PVZ6IZIrIxyIysMIxVET2l9kvpyZyNXWLFXQTTBerakOgPfAUcD/werBfREQig33M2iIiJwPTgVlAFyAR+ANwfpndtqpqAtAQGAysBb4XkbMqHK6f/x+1BFVtUuPJm5BnBd0EnaruVdVJwG+B60WkN4CIvCUif/ffTxKRL8q0Ur8XkQgReRdoB0z2tzzvE5EO/hbpKBHJAKaXeaxsce/sb+3uFZHPRaSZ/7WGikhm2RwP/AoQkeHAX4Df+l9vmT9+sAvHn9dDIpIuIjtF5B0RaeyPHcjjehHJEJHdIvLgYd6eZ4G3VfVpVd2tjkWqemUl76OqaqaqPgK8Bjxd1Xt/qPe1queZ8GD/o02NUdX5QCZwWiXhe/yxZKAFTlFVVb0OyMBp7Seo6jNlnnMG0BM47xAvORK4CWgNlAIvVCPHqcATwEf+1+tXyW43+G/DgE5AAvBShX1OBboDZwGPiEjPigcRkTjgZOCTqvKqxETgBBGJr2K/St/Xo3g9UwdZQTc1bSvQrJLHS4BWQHtVLVHV77XqiYUeU9X9qlpwiPi7qrpSVfcDDwNXHjhpeoyuAf6lqmmqmgf8GRhR4dfBX1W1QFWXAcuAyv5haIrzndt2FDlsBQRoUuaxxf6WeI6IHPjH62jeVxMmrKCbmtYG2FPJ488CqcDXIpImIg9U41ibjyCeDkQBSdXK8vBa+49X9tiROC3gA7aXuZ+P04qvKBvw4RTcI9UGp6WdU+axE1S1if92h/+xo3lfTZiwgm5qjH9kRhvgh4oxVc1V1XtUtRNwMfDHMif9DtWirKql2bbM/XY4rdXdwH4grkxeHpwuieoedyvOid6yxy4FdlTxvHJUNR/4Cfj1kTzP7zJgsf/Xx+Fe43DvqwlzVtBN0IlIIxG5CBgPvKeqKyrZ5yIR6SIiAuwDvP4bOIWy01G89LUi0svfV/034BP/sMb1QKyIXCgiUcBDQEyZ5+0AOhzm5OGHwN0i0lFEEvilz730KHK8D7hBRP4kIokAItJPRMZX3FEcbUTkUeBmnP7ww6rifTVhzgq6CabJIpKL0/XxIPAv4MZD7NsV+BbIw2m1vqKqM/2xJ4GH/H3D9x7B678LvIXT/REL3AHOqBvgVpyRIltwWuxlR71M8P+ZJSKLKznuG/5jzwY2AoXA7UeQ10GqOgc4039LE5E9wDhgSpndWotIHs57swDoAwxV1a+r8RKHe19NmBM7X2KMMeHBWujGGBMmrKAbY0yYsIJujDFhwgq6McaECdcmOUpKStIOHTq49fLGGFMnLVq0aLeqJlcWc62gd+jQgYULF7r18sYYUyeJSPqhYtblYowxYcIKujHGhAkr6MYYEyZCauWXkpISMjMzKSwsdDuVQ4qNjSUlJYWoqCi3UzHGmHJCqqBnZmbSsGFDOnTogDO3UGhRVbKyssjMzKRjx45up2OMMeVU2eUiIm/4l91aeYi4iMgLIpIqIstF5ISjTaawsJDExMSQLOYAIkJiYmJI/4IwxtRf1elDfwsYfpj4+TgzvHUFRgOvHktCoVrMDwj1/Iwx9VeVXS6qOltEOhxml0uBd/zLXM0VkSYi0kpVj2aZLWOMqR5V8JZAST6UFEBpIZQWgbcISotRbxHFRcUUFxdSUlJMSXExpaXFeEtL8XpL8JZ68XlL8fmcP9Xnw6c+1OdFfT5UFVUf6vOBHthW53VVcf7z54E6y6Soz0mtXJ6g+CgbSOh2Gr1Pvyzob0kw+tDbUH7pr0z/YwEFXURG47TiadeuXRBe2hhT56lCYQ7s2wZ5O2D/LueWvwcK9uDLz6Zkfzal+XuhcC9Ssh9PyX6ivAVEHGbtDsFZxSTmkHvULp/+8ut+nrcEQrSgV9YHUekk66o6DmcyfwYMGGATsRtTX3hLIXsT7F4HWanO/exNkJMB+7Y6reyKTyGCfSSQ7YtnH3Hs0zj2kcx+bUc+MRRHxKJR8RDVAE90HBHRDYiIiiUy2n+LiiYqOoaoqBgio6Od7ahoPJFRREZG/nLzRBHh8RAVGUWEJwJPRAQRER48kR4iRIiIiDh480gE4hEE53ER575ECCIgRCARTkkUfumide47fy8R4eQaepuDUdAzKb+WYwrOGox1zoIFCxg1ahTz58/H6/UyaNAgPvroI3r37u12asbUHaXFsH05bFkE25bD9mWwax14iw/uUhLdhKzoVmz2tSbV14vUkkbs0Kbs1CZkSWOiG7WgadMkWjWNp02TWFo0jqVFw1jaNowhMSGaxPgYGkR7XPxLhqZgFPRJwBj/mognAXuD0X/+18mrWL113zEnV1av1o149OLjDhkfOHAgl1xyCQ899BAFBQVce+21VsyNqUppEWyeD2kzIH0ObF3i9GcDxCVRlNyb9A7XsLSwJTOymvJjTlP2FcYTIdCleQI92jaie8uGnNQ8gU7JCbRt1oCYSCvWR6PKgi4iHwJDgSQRyQQeBaIAVHUszlqIFwCpQD6HXkOyTnjkkUcYOHAgsbGxvPDCC26nY0xo2r8b1k2BtV/CxtlOl4l4oPXx6IBRpMf3ZlpOWyalKavW5QLQKDaSQR0TGTOkKf3bNuW41o2IjwmpS2HqvOqMcrmqirgCtwUtI7/DtaRr0p49e8jLy6OkpITCwkLi4+NdycOYkFOUC6s/h2XjIf1HZ0RH43bQ/xroPIyf4/ozcU0uXyzfRnpWPhGyjxPbN+W+4d05vWsyPVs1whNhw35rkv3zWMHo0aN5/PHH2bhxI/fffz8vvfSS2ykZ467MhTD/f04xLy2AZp3htHuh50UUNDuOL1ZsY/yMzSxKX4onQhjSOZHbhnbh3ONa0CQu2u3s6xUr6GW88847REZGcvXVV+P1ehkyZAjTp0/nzDPPdDs1Y2qXtxRWTYS5r8LWxRDdEPqNgP5XQ8pAduYV8c6cdN6dO529BSV0So7noQt78qvj25CUECoDBesfK+hljBw5kpEjRwLg8XiYN2+eyxkZU8u8JU6XyvfPOcMKk7rBBc85xTymIVtzCnjp/1byyaJMSrw+zuvVkhtP6cCgjs3sKuoQYAXdGONc3LPyU/jub5CTDq36w4gPodtwiIggK6+IV75Zzbtz01FVrhjQlt+d1omOSXaOKZRYQTemvsuYB9P+AlsWQss+cPUE6HoOiFDi9fHODxt5/pv17C8u5dcnpHDn2V1JaRrndtamElbQjamvCrLh64dgyXuQ0BIufcXpWolwxoDPTcvi0c9XsW5HLmd0S+bhi3rSpXlDl5M2h2MF3Zj6RhVWfQZf3Q/5WXDKnXDG/RDtdJ/sLyrlya/W8N7cDNo0acC4607knF4trI+8DrCCbkx9UpADX/7R6S9v1R+u/QRa9TsYnr9xD/dOWMbm7HxGndqRe8/tbpfY1yFW0I2pLzLmwac3w74tMOwhOPVu8DglwOtTXpy+gf98t4G2TeP4aPTJDOrYzOWEzZGygm5MuFOFOS/Ct49B4xS4aRq0HXgwvDuviLs/Wsr3G3Zz+fFtePxXve2S/DrK/q8ZE86K82HyHbBiAvS6FC55CWIbHQwvz8xh9DuLyM4v5ulf9+HKAW2tr7wOs4JexsMPP0xSUhJ33nknAA8++CAtWrTgjjvucDkzY47C3kz48CrYvgLOfBhOu+eXSbmBqSu3cddHS0lKiGHirUM4rnVjF5M1wRC6Bf2rB5wPYjC17APnP3XI8KhRo7j88su588478fl8jB8/nvnz5wc3B2Nqw47V8N6voTgPrv4Iup13MKSqjJ2VxtNT13J8uyb8b+QAu1w/TIRuQXdBhw4dSExMZMmSJezYsYPjjz+exMREt9My5shs+hHGXwWRDeDGr6DlL3P6+3zK375YzVtzNnFR31Y8d0U/YqNsFEu4CN2CfpiWdE26+eabeeutt9i+fTs33XSTKzkYc9TWfQUfXw9N2sF1E50//Uq9Ph6YuIJPFmUy6tSOPHhBTyJsOtuwEuF2AqHmsssuY+rUqSxYsIDzzjuv6icYEyrWfgkfXQctjnNGspQp5sWlPu4Yv4RPFmVy19ldeehCK+bhKHRb6C6Jjo5m2LBhNGnSBI/HfoqaOmLNZJhwg3Ox0HUTIfaXE5ylXh+3f7iYaat28NCFPbn5tE6upWlqlhX0Cnw+H3PnzmXChAlup2JM9az90inmrY+Haz8tV8y9PuWeCcuYtmoHj17cixtP6ehenqbGWZdLGatXr6ZLly6cddZZdO3a1e10jKnaxu9hwo3+y/jLt8x9PuUvE1fw+dKt3D+8hxXzesBa6GX06tWLtLQ0t9Mwpnq2LXPGmTfrCNdMKHfBEMBTU9fy0cLN3HFWV/4wtLNLSZraFHItdGfN6dAV6vmZeiLrZ2eceWxjp2UeV37elbfnbGLc7DRGntyeu8+2X5v1RUgV9NjYWLKyskK2aKoqWVlZxMbGup2Kqc/y98D7V4DPC9d9Bo3blAt/vWo7j01exTm9WvDoxcfZpfz1SEh1uaSkpJCZmcmuXbvcTuWQYmNjSUlJcTsNU195S2DC9bB3M1w/GZK7lQsv3ZzDHeOX0DelCS+MOB6PDU2sV0KqoEdFRdGxo524MaZSqjDlXtg4G341FtoNLhfeua+Q37+7kKSEGF6/foDNY14PhVRBN8Ycxrz/wqK3nHnM+19VLlRU6uWW9xaxr6CUibcOsblZ6ikr6MbUBelznIWcu18IZz5SLqSqPPr5KhZn5PDy1SfQs1WjQxzEhLuQOilqjKlE7g5nrHnT9nDZqxBR/mv7wfwMxi/YzK1DO3Nh31YuJWlCgbXQjQll3lL4dBQU7g24ChRg1da9/HXyak7rmsQ953Z3KUkTKqygGxPKZvwdNn3vnAQtMw0uQF5RKWM+WELTuCie/21/G9FirKAbE7J+ngE//BtOGBlwElRVefCzFaRn7eeD3w0m0U6CGqwP3ZjQtD8LPrsFkrrB8KcDwhMWZvL50q3cdXY3BneyRViMo1oFXUSGi8g6EUkVkQcqiTcWkckiskxEVonIjcFP1Zh6QhUm3Q4Fe+DXr0N0XLlwetZ+Hpu8ipM7JXLbsC4uJWlCUZUFXUQ8wMvA+UAv4CoR6VVht9uA1araDxgK/FNEooOcqzH1w6I3Yd2XcPZj0KpvuZDXp/zx42V4IoR/XtnP+s1NOdVpoQ8CUlU1TVWLgfHApRX2UaChOJNGJAB7gNKgZmpMfbAnDaY9CJ3PhJP+EBAeO+tnFqVn8/ilvWndpIELCZpQVp2C3gbYXGY70/9YWS8BPYGtwArgTlX1VTyQiIwWkYUisjCU52sxxhU+H3w+BiKi4NKXA8abr9q6l+e/Xc+FfVtxaf/WLiVpQll1Cnplv+kqTod4HrAUaA30B14SkYDL1VR1nKoOUNUBycnJR5iqMWFu/jhI/xGGPwmNyhfs4lIf93y8jKZx0fzjV71tBkVTqeoU9EygbZntFJyWeFk3AhPVkQpsBHoEJ0Vj6oGsn+Hbx6DredD/6oDw2Fk/s3Z7Lk9c1ocmcXZ6ylSuOgV9AdBVRDr6T3SOACZV2CcDOAtARFoA3QFb+seY6vD54PPbwBMNFz8PFVrfG3bk8uL0DVzSrzVn92rhTo6mTqjywiJVLRWRMcA0wAO8oaqrROQWf3ws8DjwloiswOmiuV9Vd9dg3saEj8VvQcZPTr95ha4Wr0+579PlJMRE8ujFFQeXGVNeta4UVdUpwJQKj40tc38rcG5wUzOmHsjdDt88Bh1Og/7XBITfnrOJJRk5/GdEf7sa1FTJrhQ1xk1TH4DSQrjo+YCulq05BTz39TqGdU/mkn42qsVUzQq6MW5Z/zWs+gxOvxeSAq/4fPyL1fhU+dulNqrFVI8VdGPcUJwPU+6BpO5wyl0B4RnrdvLVyu3cfmZX2jaLC3y+MZWw2RaNccMP/4acDLhhCkSWH4ZYWOLl0c9X0Tk5nt+d1smlBE1dZAXdmNq2Jw1+/A/0uRI6nBIQfmVGKhl78vngdycRHWk/ok312afFmNo29c/giYJz/hYQysjKZ+zsNC7t35ohnZNcSM7UZVbQjalN66bC+qkw9AFoFLj+59+/XE1khPCXC3q6kJyp66ygG1NbSgph6v3OidCTbgkI/7BhN1+v3sFtw7rQolGsCwmaus760I2pLfPGQvYmuO4zp8uljFKvj79OXkW7ZnGMOrWjO/mZOs9a6MbUhrxdMPs56Dbcmeu8gvfmprNhZx4PXtiT2CiPCwmacGAF3ZjaMOMfUFoA5/49IJSTX8y/v93AqV2SONcm3zLHwAq6MTVtxypY/DYMvBmSugaEX5yeSm5hCQ9d1NOuCDXHxAq6MTVJ1VlSLqYRnHF/QDgjK593ftrEFSe2pUfLgDVhjDkiVtCNqUmp30LaDKeYxzULCD89bS2RERH88dxuLiRnwo0VdGNqis8L3zwKTTs63S0VLErP5svl2xh9eicbpmiCwoYtGlNTln8MO1fBb94ImK9FVXlyyhqSG8Yw+nSbr8UEh7XQjakJJYXOyJZW/aHXZQHhb1bvYGF6Nnef3Y34GGtXmeCwT5IxNWH+ONi72VlWLqJ8u8nrU56dto5OyfFcOSDFpQRNOLIWujHBVpAN3/8TupwNnc4ICE9cnMmGnXn86dzuRHrsK2iCxz5NxgTbjy9A4V44+7GAUGGJl39/s55+KY0Z3rtl7edmwpoVdGOCKXeHM2dLn99Ayz4B4ffmprN1byH3D+9hFxGZoLOCbkwwff9PKC2CoX8OCOUWlvDyjFRO65rEkC4217kJPivoxgRLTgYsfANOuA4SOweE3/hhE9n5JfzpvO4uJGfqAyvoxgTLzKdBIuD0+wJCOfnFvPZ9Gucd14K+KU1qPzdTL1hBNyYYdq2HZR/AoN9B4zYB4XGz08grLuXuc+wSf1NzrKAbEwyznoLIBnDq3QGh3XlFvPnjJi7u29om4DI1ygq6Mcdqx2pYORFO+j3EB57sfHXmzxSVernr7MCpc40JJivoxhyrWU9BdAIMuT0gtH1vIe/OTefXJ6TQKTnBheRMfWIF3ZhjsX0FrP4cTr610ulxX52Zis+n3HGWtc5NzbOCbsyxmPkUxDSGwbcGhLbtLeDD+Zu5YkAKbZvFuZCcqW+qVdBFZLiIrBORVBF54BD7DBWRpSKySkRmBTdNY0LQ1iWw9gsYMgYaNAkIvzLjZ3yq3Dq0S+3nZuqlKmdbFBEP8DJwDpAJLBCRSaq6usw+TYBXgOGqmiEizWsoX2NCx8ynIbYJnHRLQGhrTgEfLdjMFQPaWuvc1JrqtNAHAamqmqaqxcB44NIK+1wNTFTVDABV3RncNI0JMVuXwvqv4OQxEBs4FPGVmakoym3DAq8YNaamVKegtwE2l9nO9D9WVjegqYjMFJFFIjKysgOJyGgRWSgiC3ft2nV0GRsTCmY9A7GN4aTRAaGyrfOUptY6N7WnOgW9sinhtMJ2JHAicCFwHvCwiARcEqeq41R1gKoOSE5OPuJkjQkJ25bBui9h8G1OUa/g1Zk/A3DbMOs7N7WrOisWZQJty2ynAFsr2We3qu4H9ovIbKAfsD4oWRoTSg62zn8fENq+t5CPFmzmNye2pU2TBi4kZ+qz6rTQFwBdRaSjiEQDI4BJFfb5HDhNRCJFJA44CVgT3FSNCQHbVzgjWwbfWunIlrGzDoxssb5zU/uqbKGraqmIjAGmAR7gDVVdJSK3+ONjVXWNiEwFlgM+4DVVXVmTiRvjilnPQEyjSke27NxXyIfzM7j8hDY2ssW4olqLRKvqFGBKhcfGVth+Fng2eKkZE2J2roE1k+D0P1XaOh83O41Sn1rfuXGNXSlqTHXNfs6Zs6WSq0J35xXx3rx0Lu3fmvaJ8S4kZ4wVdGOqZ/cGWDURBo6qdM6W177fSFGpz1rnxlVW0I2pju//CZ4YODlwRsXs/cW8+9MmLurbms42o6JxkRV0Y6qyZyMs/xgG3AQJgddPvDlnE/uLvYyx1rlxmRV0Y6ryw78hIrLS+c73FZbw5o8bGX5cS7q3bOhCcsb8wgq6MYezNxOWfgAnXAeNWgWE35mzidzCUsacaa1z4z4r6MYczo8vAAqn3BkQ2l9Uyus/bOTMHs3p3SZwCgBjapsVdGMOJXcHLH4b+o2AJu0Cwh/MyyA7v8RGtpiQYQXdmEP56SXwFsOpfwwIFZZ4Gfd9Gqd0SeTE9k1dSM6YQFbQjalM/h5Y8Dr0/jUkBs7L8vHCzezKLWLMMFsr1IQOK+jGVGbuq1CyH067JyBUXOrjv7PSOLF9UwZ3CrzIyBi3WEE3pqLCvTDvv9DjImjeMyD8f0u2sCWngDFndkGksuUCjHGHFXRjKpr/Pyja60zCVUGp18crM1Pp3aYRQ7vZIi0mtFhBN6as4v0w9xXocg607h8Q/nLFNjZl5TNmmLXOTeixgm5MWYvegvysSlvnPp/y8oxUujZP4NxeLWs/N2OqYAXdmANKCmHOi9DhNGh3UkD4mzU7WL8jj9uGdSEiwlrnJvRYQTfmgKXvQ+42OP3egJCq8tL0VNonxnFR38ApAIwJBVbQjQHwlsCPz0ObAdDxjIDwrPW7WLFlL384ozORHvvamNBkn0xjAFZMgJwMp++8wsnOA63zVo1jufyEFJcSNKZqVtCN8XmdBSxa9IFu5wWE523cw8L0bG45ozPRkfaVMaHLPp3GrP4/yEp1+s4rGYr40vRUkhJi+O3AtrWfmzFHwAq6qd98Pmfx56Tu0POSgPDijGx+SN3N6NM7EhvlcSFBY6rPCrqp39Z/BTtXO3O2RAR+HV6ankqTuCiuOam9C8kZc2SsoJv6SxVmPwtNOzizKlawcstepq/dyahTOhIfE1n7+RlzhKygm/or9VvYusSZ79wTWLBfnpFKw9hIrj+lQ+3nZsxRsIJu6idVmPUMNG4L/a4KCK/fkctXK7dzw5AONIqNciFBY46cFXRTP22cBZnznbVCI6MDwi/PSCUu2sONp3R0ITljjo4VdFM/zXoWGraC468LCP28K4/Jy7Zy3cntaRYfWOyNCVVW0E39kz4H0n9wWudRsQHhl2ekEh0Zwe9O6+RCcsYcPSvopv6Z9QzEJ8MJ1weENu3ez+dLt3LtSe1JSohxITljjp4VdFO/bJ4PaTNgyB0QHRcQfmVmKpERwujTrXVu6p5qFXQRGS4i60QkVUQeOMx+A0XEKyK/CV6KxgTRzKcgLgkGjgoIbd6Tz8TFW7hqUDuaNwrsijEm1FVZ0EXEA7wMnA/0Aq4SkV6H2O9pYFqwkzQmKDYvgJ+/gyG3Q3R8QPiVmalEiHDLGZ1dSM6YY1edFvogIFVV01S1GBgPXFrJfrcDnwI7g5ifMcEz6ymIS4SBNweENu/JZ8LCTH47sC0tG1vr3NRN1SnobYDNZbYz/Y8dJCJtgMuAsYc7kIiMFpGFIrJw165dR5qrMUcvc6FzZejJYyAmISD88gyndX7rMGudm7qrOgW9ssUTtcL288D9quo93IFUdZyqDlDVAcnJydVM0ZggmPkUNGgKg34XENq8J59PFmVy1aC2tGrcwIXkjAmO6sw4lAmUnQg6BdhaYZ8BwHhx5pJOAi4QkVJV/b9gJGnMMdm8AFK/gbMegZiGAeGXpqcSESH8YWgXF5IzJniqU9AXAF1FpCOwBRgBXF12B1U9eH20iLwFfGHF3ISMmU84feeDRgeEMrLy+XRxJtcObm9956bOq7LLRVVLgTE4o1fWAB+r6ioRuUVEbqnpBI05Juk/wc/TnatCK2udz9jgb51b37mp+6o1ybOqTgGmVHis0hOgqnrDsadlTJDMfALim8PAwL7zjbv38+niLYw8uT0tbNy5CQN2pagJXxu/h42z4dS7K70q9D/frifaE2GtcxM2rKCb8KQKM56AhJYw4MaA8PoduXy+bCvXD+lA84bWOjfhwQq6CU8/T4eMOc5aoVGBQxGf/3Y98dGR/N7mbDFhxAq6CT+q8N3foHE7ODFwRsWVW/YyZcV2Rp3akaY237kJI1bQTfhZMxm2LYWhD0Bk4BS4//pmPY0bRDHqNFuNyIQXK+gmvPi8MP3vkNQd+o0ICC/YtIfpa3cy+vROtlaoCTvVGrZoTJ2x/CPYvQ6ufAciPOVCqsrTX62lecMYbrK1Qk0Ysha6CR+lRTDzSWjVH3peEhCevnYnC9OzueOsrjSI9gQ+35g6zlroJnwsfANyMuCif4OUn1PO61OembqODolx/HZg20McwJi6zVroJjwU7nXWCu14BnQ+KyA8adkW1u3I5Z5zuxPlsY+9CU/2yTbh4cf/QMEeOOevAa3zolIv//x6Pce1bsSFfVq5lKAxNc8Kuqn79m2Fn16B3r+B1scHhN+Zk05mdgEPnN+DiIjKpvc3JjxYQTd138wnwVcKZz0cEMrJL+bF6Rs4o1syp3W1RVVMeLOCbuq2nWtgyXvOOqFNOwSEX5yeSl5RKX+5oGft52ZMLbOCbuq2aQ9CdEM4/U8BofSs/bzz0yauOLEt3VsGzoVuTLixgm7qrg3fws/fwRn3QXxiQPiZaeuIjIjgj+d2cyE5Y2qfFXRTN3lL4esHoVmnSpeWW7BpD18u38bo0zvZ4hWm3rALi0zdtOhN2LUWfvs+RJafMdHrUx6btIpWjWP5/Rk2Pa6pP6yFbuqeghxnZEuH06DHhQHhTxZtZtXWfTxwfg/ioq3NYuoPK+im7pn5JBRkw3lPBFxEtK+whGenrWNA+6Zc0q+1Swka4w4r6KZu2b4S5o+DATdBq74B4Re/20DW/mIevfg4ROwiIlO/WEE3dYcqTPkTxDaBYQ8GhFN35vLWnE1ccWIKfVIa135+xrjMOhhN3bHiE2ed0Iv/A3HNyoVUlYf/bxUNojzcN7yHSwka4y5roZu6oSgXvn7Imavl+OsCwp8v3cpPaVncN7wHSQmBy84ZUx9YC93UDdP/AXk7YMT7ASsR7S0o4e9frqFf2yZcNaidSwka4z4r6Cb0bVkM8/8LA0dByoCA8L++Xsee/UW8ecNAPDaboqnHrMvFhDZvKUy+A+Kbw1mPBISXbc7h3bnpXDe4vZ0INfWetdBNaJv3KmxfAVe8DbHlC3ZxqY/7P11O84ax3HNed5cSNCZ0WEE3oSsnA2Y8Ad2GQ69LA8JjZ/3M2u25vDZyAI1io1xI0JjQYl0uJjSpwqTbQSLggucCrgjdsCOXl6ancnG/1pzdq4VLSRoTWqpV0EVkuIisE5FUEXmgkvg1IrLcf5sjIv2Cn6qpVxa9CWkz4dzHoUnbciGvT7n/0+XEx3h49OJe7uRnTAiqsqCLiAd4GTgf6AVcJSIVv0UbgTNUtS/wODAu2ImaeiQ7Hb5+GDoNhRNvDAi//kMaizNyeOTiXjbm3JgyqtNCHwSkqmqaqhYD44FyHZqqOkdVs/2bc4GU4KZp6g2fDyaNAQQueTGgq2Xd9lyem7ae845rwa/6t3EnR2NCVHUKehtgc5ntTP9jhzIK+KqygIiMFpGFIrJw165d1c/S1B8LX4eNs/1dLeUvEiou9XHXR0tp1CCSJy7rY5NvGVNBdQp6Zd8arXRHkWE4Bf3+yuKqOk5VB6jqgORkW4HdVLBjtXN5f5dz4MQbAsL/+W49a7bt48nL+5JoXS3GBKjOsMVMoOxZqRRga8WdRKQv8BpwvqpmBSc9U2+UFMKnoyCmIfzq1YCulgWb9vDqzJ+5ckAK59ioFmMqVZ0W+gKgq4h0FJFoYAQwqewOItIOmAhcp6rrg5+mCXvfPgo7VzvFPKH8r7fs/cXc8eES2jaL4+GLbFSLMYdSZQtdVUtFZAwwDfAAb6jqKhG5xR8fCzwCJAKv+Ps1S1U1cNINYyqzfhrMGwsn/QG6nlMupKrcO2EZWXnFTLx1CA3tAiJjDqlaV4qq6hRgSoXHxpa5fzNwc3BTM/VCdjp89nto0QfOfiwg/PoPG/lu7U4eu7gXvdvYXC3GHI5dKWrcU1IIH490hipe+TZExZYLL8nI5umpazm3VwuuH9LBnRyNqUNsLhfjnqkPwLalMOIDSOxcLrQzt5Bb3ltEi0axPPubfjZE0ZhqsIJu3LHkfefy/lPugh4XlgsVl/q49b3F7C0oYeIfTqFxnPWbG1MdVtBN7cuYC1/cBR1PhzMfDgj/dfIqFqZn8+JVx9OrdaPaz8+YOsr60E3tyk6H8ddA47bOHOee8m2K9+am8/68DH5/Ricu7tfapSSNqZusoJvaU7gPPhwBvhK4+iOIa1YuPH3tDh75fCXDuidz33k9XErSmLrLulxM7Sgthgk3wK51cO2nkNS1XHhF5l7GfLCEXq0b8dLVJ9jaoMYcBSvopub5fPD5rfDzd84Mip2HlQtnZudz09sLaBoXzRvXDyQ+xj6WxhwN++aYmqUK0/4MKybAWY/CCSPLhXfuK+Ta1+ZRWOLl/ZtPonmj2EMcyBhTFetDNzVr9rPOZf2Db4NT7y4X2rO/mGtfn8fO3CLeunEQ3Vo0dClJY8KDFXRTc2Y9CzP+Af2ugnP/Xm4Gxb0FJYx8Yx7pWfm8dv0ATmzf1MVEjQkP1uViasasZ5xi3ncEXPoyRPzSdsjeX8wNb85n3fZcxo0cwJDOSS4makz4sIJugksVZj4Fs55yWuaXvgwRnoPhnbmFXPfafDZm7efVa05kWPfmLiZrTHixgm6Cx+eFKffCwjeg/zXOiJYyxXxLTgHX/G8uO3OLePOGgZzSxVrmxgSTFXQTHAdWHFr7hTM/y9mPleszX7llLze9tYCCEi/vjjrJ+syNqQFW0M2xy9sJH10Hm+fC8Kdg8B/Khb9etZ07xy+lWXw07446ie4tbTSLMTXBCro5NluXOHOz5O+B37wJvS8/GFJV/vd9Gk9+tZa+KU3438gTad7QxpkbU1OsoJujt+wjmHwHxCfDqGnQqt/B0N6CEu77ZBnTVu3gwj6teO6KfjSI9hzmYMaYY2UF3Ry5olz46n5Y+j60P9VZbSj+lxOcK7fs5bYPFrMlu4CHLuzJqFM72gIVxtQCK+jmyGxZBJ/eDNmb4Iz74fT7Dk6BW+r18d/ZaTz/7XoS42MYP3owAzo0O/zxjDFBYwXdVE9JAcx6Gn58ARq2ghu+hPZDDobTduVxz4RlLMnI4cI+rXj8V71pFh/tYsLG1D9W0E3VNs6GyXfCnjQ4/lrnMv4GzrDDgmIvr85MZeysNOJiPLx41fG2MIUxLrGCbg5tTxp88yismQRNO8LISdDpDMAZwfL16h08/sVqMrML+FX/1vzlgp42W6IxLrKCbgLl7YIf/g3zx4EnGoY9CENuh6gGAMxNy+KZqWtZnJFDtxYJjB89mMGdEl1O2hhjBd38IncHzHnBuXS/pMDpXjnzIWjYElVl7s9ZvDIzle837KZlo1ieuKwPVwxIIcpjk3YaEwqsoBvYvgLm/ddZhMJbDH2uhNPugeRulHh9TFu+lf/NTmNZ5l6SEqL5ywU9GHlyB2KjbFy5MaHECnp9VZTnzLuy+B1I/xEiG0C/ETDkDkjszOY9+Xw4dS0fL8xkd14RHZPi+cdlvfn1CSlWyI0JUVbQ65OSQtg4C1Z9BqsnQcl+aNoBznkcjr+Wnd44pizfxuTlc1iUnk2EwJk9WnDVoLYM7d7cFm42JsRZQQ93+7ZB2gzY8DVs+AaK8yCmMfT5Db6+I1gd2Yvp63Yx/c01LMvMQRV6tGzIn87rzuUntKFV4wZu/w2MMdVkBT2cqEJOBmye78x8uOkH2LXWicUn4z3ucjYln8n3pb2YsymX+e/sISf/R0Sgb0oT7jqrGxf2bUmX5jYbojF1kRX0uqo4H7I2wK71sGMlbF8O25ZD/m4AfFHx7Es6gQ1dh/MTfZm+J5lV8/Mo8SqQSrtmcZzbqwWDOyVyerdkkhJi3P37GGOOWbUKuogMB/4DeIDXVPWpCnHxxy8A8oEbVHVxkHOtP1ShINuZZzx3G+zbAvu2Qk46ZKej2Rth7xYEBcArkexq0Im0yIEsi+/Ad3kdWZzbGl+uM5ywSVwUPVtGM+rUTvRLaUz/dk2sK8WYMFRlQRcRD/AycA6QCSwQkUmqurrMbucDXf23k4BX/X+GL58PfKXgKwFviXPfWwzeYrS0GF9pEaVFBfhKCvAW5eMtzsdXlI+vKA9fUR5alIcW7oOifUjhXjxFOXiKcoguziGmeA8e9Qa85G5pSoY2Z5O3A5t8g0nV1qRqGzZpSyiJoU2TBqQ0a8Bx3eK5MCmeTskJ9GjZkOYNY2y2Q2Pqgeq00AcBqaqaBiAi44FLgbIF/VLgHVVVYK6INBGRVqq6LdgJL5/5CY1nP3pw+0ArFUC0zP2D9/TgPlJmu/x9Jx6BD1EQfAjqbKN4ymx78OHBS0SZ161IcH7KVDW4L1cbkEsDcjWOHBLI0cZkaxuyaMQempAX1Yz86GTyG7TAG9+C+PgEEuOjSYyPoXlCNH0axdKiUQwtG8WSlBBDhI1CMaZeq05BbwNsLrOdSWDru7J92gDlCrqIjAZGA7Rr1+5IcwUgOr4Ju+O6lHtMy5TvsutYli3vv+wjIGW2D94XVJyyfWB/Fc8vj4kHn0SgEoFGRPpjHlQinW1PFL6IKDQiCiKi0MhY1BMDnmg0Kg6JauDcYhtCdAKe2Hhio6OIiYwgNspD02gPraM8xEV7iI+JJCYywlrVxpgjUp2CXllVqdg8rc4+qOo4YBzAgAEDDt3EPYweA8+GgWcfzVONMSasVWcSjkygbZntFGDrUexjjDGmBlWnoC8AuopIRxGJBkYAkyrsMwkYKY7BwN6a6D83xhhzaFV2uahqqYiMAabhnOd7Q1VXicgt/vhYYArOkMVUnGGLN9ZcysYYYypTrXHoqjoFp2iXfWxsmfsK3Bbc1IwxxhwJm8jaGGPChBV0Y4wJE1bQjTEmTFhBN8aYMCGqR3V9z7G/sMguIP0on54E7A5iOsESqnlB6OZmeR0Zy+vIhGNe7VU1ubKAawX9WIjIQlUd4HYeFYVqXhC6uVleR8byOjL1LS/rcjHGmDBhBd0YY8JEXS3o49xO4BBCNS8I3dwsryNjeR2ZepVXnexDN8YYE6iuttCNMcZUYAXdGGPCRJ0v6CJyr4ioiCS5nQuAiDwuIstFZKmIfC0ird3OCUBEnhWRtf7cPhORJm7nBCAiV4jIKhHxiYjrw8tEZLiIrBORVBF5wO18DhCRN0Rkp4isdDuXA0SkrYjMEJE1/v+Hd7qdE4CIxIrIfBFZ5s/rr27nVJaIeERkiYh8Eexj1+mCLiJtcRavznA7lzKeVdW+qtof+AJ4xOV8DvgG6K2qfYH1wJ9dzueAlcDlwGy3EymzIPr5QC/gKhHp5W5WB70FDHc7iQpKgXtUtScwGLgtRN6vIuBMVe0H9AeG+9dpCBV3Amtq4sB1uqAD/wbuo5Ll7tyiqvvKbMYTIrmp6teqWurfnIuzqpTrVHWNqq5zOw+/gwuiq2oxcGBBdNep6mxgj9t5lKWq21R1sf9+Lk6RauNuVs503qqa59+M8t9C4nsoIinAhcBrNXH8OlvQReQSYIuqLnM7l4pE5B8ishm4htBpoZd1E/CV20mEoEMtdm6qICIdgOOBeS6nAhzs1lgK7AS+UdWQyAt4HqcR6quJg1drgQu3iMi3QMtKQg8CfwHOrd2MHIfLS1U/V9UHgQdF5M/AGODRUMjLv8+DOD+V36+NnKqbV4io1mLnpjwRSQA+Be6q8AvVNarqBfr7zxV9JiK9VdXV8w8ichGwU1UXicjQmniNkC7oqnp2ZY+LSB+gI7BMRMDpPlgsIoNUdbtbeVXiA+BLaqmgV5WXiFwPXAScpbV4AcIRvF9us8XOj5CIROEU8/dVdaLb+VSkqjkiMhPn/IPbJ5RPAS4RkQuAWKCRiLynqtcG6wXqZJeLqq5Q1eaq2kFVO+B8EU+ojWJeFRHpWmbzEmCtW7mUJSLDgfuBS1Q13+18QlR1FkQ3fuK0pl4H1qjqv9zO5wARST4wiktEGgBnEwLfQ1X9s6qm+GvWCGB6MIs51NGCHuKeEpGVIrIcp0soJIZyAS8BDYFv/EMqx1b1hNogIpeJSCZwMvCliExzKxf/SeMDC6KvAT5W1VVu5VOWiHwI/AR0F5FMERnldk44Lc7rgDP9n6ml/tan21oBM/zfwQU4fehBHyIYiuzSf2OMCRPWQjfGmDBhBd0YY8KEFXRjjAkTVtCNMSZMWEE3xpgwYQXdGGPChBV0Y4wJE/8Pp3DRw+UPwvsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(grid, X.cdf(grid), 'C0')\n",
    "plt.plot(grid, Y.cdf(grid), 'C1')\n",
    "plt.title('Distribution CDFs')\n",
    "plt.legend(['x', 'y'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c44e1bfe-8ebe-41d1-b96d-90f970cb8311",
   "metadata": {},
   "source": [
    "Therefore, to test the alternative that $X$ is stochastically greater than $Y$, we pass `alternative='less'` into `ks_2samp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "97815b21-eb3e-401e-9cff-1da838141a70",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KstestResult(statistic=0.12484394506866417, pvalue=0.00022195733373729093)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\matth\\AppData\\Local\\Temp\\ipykernel_35644\\3756467023.py:1: RuntimeWarning: ks_2samp: Exact calculation unsuccessful. Switching to method=asymp.\n",
      "  res1 = stats.ks_2samp(x, y, alternative='less', method='exact')\n"
     ]
    }
   ],
   "source": [
    "res1 = stats.ks_2samp(x, y, alternative='less', method='exact')\n",
    "print(res1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33065824-d0c3-4ad6-9c08-296a160636fd",
   "metadata": {},
   "source": [
    "The $p$-value is tiny, confirming what we already know: the data are inconsistent with the null hypothesis, and we have evidence to reject it in favor of the alternative.\n",
    "\n",
    "The warning states that `ks_2samp` was unable to compute an exact $p$-value, and an asymptotic $p$-value is being returned instead. To determine whether the asymptotic $p$-value is accurate for these sample sizes, we can perform a permutation test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "649978b2-02ff-42a6-90d1-e04b3051c31c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def statistic(x, y):\n",
    "    return stats.ks_2samp(x, y, alternative='less').statistic\n",
    "\n",
    "# This would be extremely slow!\n",
    "# res2 = stats.permutation_test((x, y), statistic, alternative='greater')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69c52b5c-f7aa-4eaf-b2ad-9e53723aa9c7",
   "metadata": {},
   "source": [
    "The calculation above would be extremely slow to run. Unfortunately, `ks_2samp` does not accept an `axis` argument, so we can't speed it up using vectorization without truly implementing the statistic ourselves. However, lack of vectorization is not be the bottleneck here. Note that the call to `ks_2samp` is quite slow with the default parameters *even for 1D inputs*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4825ddb9-5bf9-445f-bb82-031d3b330fc8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "374 ms ± 2.57 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# No need for the warning; we know the exact calculation is unsuccessful\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "%timeit stats.ks_2samp(x, y, alternative='less', method='exact')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f728c880-bbde-4003-87a6-c4090925088d",
   "metadata": {},
   "source": [
    "By default, `permutation_test` needs to call `ks_2samp` 9999 times, which would take about an hour. We can speed this up dramatically by noting `permutation_test` only uses `ks_2samp` to compute the test statistic, so the `pvalue` attribute of the `ks_2samp` result object is not used at all. We can use `ks_2samp` to compute essentially the same value of the test statistic, but much faster, by specifying `method='asymp'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "abfe9b45-4e7b-4284-bd1b-2620eced89ff",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "167 µs ± 947 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# method='asymp' and method='exact' result in the same statistic value\n",
    "res1 = stats.ks_2samp(x, y, alternative='less', method='asymp')\n",
    "res2 = stats.ks_2samp(x, y, alternative='less', method='exact')\n",
    "np.testing.assert_allclose(res1.statistic, res2.statistic, atol=1e-15)\n",
    "\n",
    "# but method='asymp' is much faster\n",
    "%timeit stats.ks_2samp(x, y, alternative='less', method='asymp')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bafdc39b-6132-4121-b7b8-6853603b129a",
   "metadata": {},
   "source": [
    "Now we can run a randomized `permutation_test` in reasonable time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c0c64b0b-31c5-4586-b696-36f03642040b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002219573337372917 0.9999\n"
     ]
    }
   ],
   "source": [
    "def statistic(x, y):\n",
    "    return stats.ks_2samp(x, y, alternative='less', method='asymp').statistic\n",
    "\n",
    "res3 = stats.permutation_test((x, y), statistic, alternative='less')\n",
    "print(res1.pvalue, res3.pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9cb3241-fea3-4b3d-90af-a60a7852ffbd",
   "metadata": {},
   "source": [
    "This was much faster, but something is still wrong. Either the approximate $p$-value is wildly inaccurate, or we have set up our test incorrectly. The latter turns out to be the case: the value of `alternative` passed into `ks_2samp` changes *the definition of the test statistic*, but a *greater* statistic is always considered more extreme. Therefore, even if we wish to perform a test equivalent to `ks_2samp` with `alternative='less'`, we actually need to pass `alternative='greater'` into `permutation_test`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3fc268c6-db67-4bee-9c3e-b50ea2896e98",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002219573337372917 0.0005\n"
     ]
    }
   ],
   "source": [
    "# greater values of the statistic returned by `ks_2samp` are more extreme\n",
    "res4= stats.permutation_test((x, y), statistic, alternative='greater')\n",
    "print(res1.pvalue, res4.pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b06fd03-ca61-498c-9542-5dabc43bcddd",
   "metadata": {},
   "source": [
    "At last, `permutation_test` is invoked correctly. Indeed, the asymptotic $p$-value produced by `ks_2samp` appears to be reliable for these sample sizes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8013eef-9aba-475d-b83a-ffbbbe311e8d",
   "metadata": {},
   "source": [
    "### Other Tests\n",
    "As we can see, `permutation_test` with `permutation_type='independent'` is a versatile tool for comparing independent samples. Provided only data and a statistic, it can produce the null distribution and replicate the $p$-value of many such tests in SciPy, and it may be more accurate than these existing implementations, especially for small samples and when there are ties:\n",
    "\n",
    "- [`ttest_ind`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html)\n",
    "- [`cramervonmises_2samp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.cramervonmises_2samp.html)\n",
    "- [`ks_2samp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ks_2samp.html)\n",
    "- [`epps_singleton_2samp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.epps_singleton_2samp.html)\n",
    "- [`mannwhitneyu`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mannwhitneyu.html)\n",
    "- [`kruskal`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kruskal.html)\n",
    "- [`friedmanchisquare`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.friedmanchisquare.html)\n",
    "- [`brunnermunzel`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.brunnermunzel.html)\n",
    "- [`ansari`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ansari.html)\n",
    "- [`bartlett`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.bartlett.html)\n",
    "- [`levene`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.levene.html)\n",
    "- [`anderson_ksamp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.anderson_ksamp.html)\n",
    "- [`fligner`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.fligner.html)\n",
    "- [`median_test`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.median_test.html)\n",
    "- [`mood`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mood.html)\n",
    "\n",
    "In addition, `permutation_test` with `permutation_type='independent'` can be used to perform tests not yet implemented in SciPy.\n",
    "\n",
    "However, there are other types of permutation tests that do not assume that the samples are entirely independent. We continue the study of `permutation_test` with [Paired-Sample Tests](https://nbviewer.org/github/scipy/scipy-cookbook/blob/main/ipython/ResamplingAndMonteCarloMethods/resampling_tutorial_2b.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
}