diff --git a/examples/mann_whitney.ipynb b/examples/mann_whitney.ipynb new file mode 100644 index 000000000..f4993bb9d --- /dev/null +++ b/examples/mann_whitney.ipynb @@ -0,0 +1,407 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Mann Whitney U test\n", + "\n", + "Allen Downey\n", + "\n", + "[MIT License](https://en.wikipedia.org/wiki/MIT_License)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "import numpy as np\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "sns.set(style='white')\n", + "\n", + "np.random.seed(18)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "[I saw a tweet recently](https://twitter.com/MaartenvSmeden/status/1100382052367691776) that make me realize that many people are confused about the [Mann-Whitney U test](https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test). A common misconception is that it tests for a difference in medians. Apparently that is true under a specific assumption, but not true in general.\n", + "\n", + "I have never been confused about the Mann-Whitney U test, because I believe \"There is Only One Test\". As I explain in [this article](https://allendowney.blogspot.com/2016/06/there-is-still-only-one-test.html), all hypothesis tests fit into a general framework; the only differences are:\n", + "\n", + "1. The test statistic, and\n", + "2. A model of the null hypothesis.\n", + "\n", + "If you have these two pieces, it is easy to compute the sampling distribution of the test statistic under the null hypothesis, and from that you can compute a p-value.\n", + "\n", + "To demonstrate, I'll use this framework to implement a version of the Mann-Whitney U test." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Example data\n", + "\n", + "First I'll generate samples from two gamma distributions with slighly different parameters." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "n1 = 90\n", + "n2 = 110" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(2.62509416921364, 1.257132162349856)" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "group1 = np.random.gamma(5, 0.5, size=n1)\n", + "group1.mean(), group1.std()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(2.4520533570340324, 1.0190315744898875)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "group2 = np.random.gamma(4.9, 0.48, size=n2)\n", + "group2.mean(), group2.std()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's what the distributions look like for the two groups." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEXCAYAAACDChKsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8lNW9+PHPzGQmy2TfSEJYAgknIIQ1LIKi4gK4C2qrrdrb1vZWbxdvF++19WfttbXtta2trd7Wtra1qK37BiqCyiL7EiBwQlizr2SbJJPM8vtjJjSEJJOETCaZfN+vV15knuecZ76ThPnOWZ5zDG63GyGEEKK/jIEOQAghxMgkCUQIIcSASAIRQggxIJJAhBBCDIgkECGEEAMiCUQIIcSAhAQ6ABFYSik3cBBwdjl1k9b6ZC/13gfu0FpXK6XeBb6ttc4fhHhygS9qrb/az3pPAdVa60e6HL8HeBI4AbgBA2DzxvupUuoR4D6gxFvFDOQBD2utj3qv8REwAajvfG2t9aw+xmYCXgWmAr/WWj/Vn9fW6ToPA/u11m8opZ4DDmqt/3cg1/Je70uARWv9u4FeQ4xukkAEwOVa6+p+1rmq4xut9cpBjOUiIH0QrwewSWt9XccDpdT1wKtKqXHeQy9pre/vdP7zwAal1EVa6wbv4e9orV8e4POPBa4BrFrrrom6P64ALjhJd7IEz4cHIQZEEojokVIqEvgzkAW4gN3AV4A/eotsVEqtBDYBq4FI4CfAaUDh+aT/OPB17+NXtNbfUkoZgV8CC4EoPK2CL3nrPQrEKKX+rLX+gvfN/vuABWjmXy2HaOBZYCZQBjiAzX18aR8CKUBsdye11n/zJpE7gGd8/Ixu8cbnwtOK+47W+pNO56OAdXhaNruVUquANODnQATQBnxfa73O21r6ImAF6rXWl3e6zn3APODnSqmOJHSxUmorMAZPIrhDa21TSk3F0+pKAEx4Wj1/6hL3zcANwFVKqVDgv4Ax3vr/B2RrrZd6yx71ljUCT3mv6wae0Fr/tZufSS7wNJ7f2TE8rbcHvKefxPN3EQnkAnfj+ftwAhXA/Vrrgq4trM6PlVIngRfwfIiJ9cbxdE9/r1prV9cYxeCQMRABnkSwr9PXa97jNwNR3q6aXO+xSVrrL3i/v1xrXdTlWrnA4946DXjemK4F5gD3KaXSgAV43kQXaa2nAX8BHvRe62E8LYYvKKWygB8DK7XWs4F78bQcrMAPgRYgG7gVT4LySSll8F7noI9W135gRqfHP+/yM+podf0c+JrWeh7wA+CyzhfRWjcCK4EW78+kDngZ+IbWOgfPG+jzSqkMb5WLgMs6Jw/vdX4L7MKToDp+P2OBK4EpeFpttyilQrzXf1BrPRdYCnxbKbWwy/VeA94Efqm1fgzYAXQ852WeH5WKVEpNA9qBo97yv/HGvQL4sVJqUefrep//VeAH3nK/Bjp39U0HPus9txj4Lp6/o5nAGuB17+/Il3g8f2uXAY8qpWbQw99rH64lBkhaIAJ67sLajOdN4iPgA+BXWutCH9c6obXe6/3+GJ5P0m1AtVKqAYj3tiC+D3xFKTUZz5tAYzfXugpIBT5U6mx+cAGZeN44v6m1dgNVnZJedy5RSu3D86k5FDgCrPLxOtx4WjwdeurCehF4TSn1Dp6f0c98XHcBUKi13g6gtT6klNqC52fgBvI6dZv58rrWuhlAKXUQSMaTTCYDf+r0MwsHZgPbernWa8AKpdQxPONBB/EknxzgFe91w7TWr3rjLlVKvQIsBz7tdJ0Z3vNrvf9u9MbWoUhrfcr7/XI83YdV3rLPKaWeBCb24bX/1vu7L1ZKrQOuxpO4+vv3Ki6AJBDRI631CaVUJp43tyuA9Uqpe7XWb/VSzd7lcXvXAkqpa/F0ZTwBvIHnDf1z3VzLBHyotb69U91xQKn3YedPqo5eYjpnDKSPcvF0h/RKa/2QUupPeJLdPcB/AvN7qWLCkyg6M+Lp4moDmvoRY+efbccEAROepH32U79SagxdJgB04zXgE6AAz5vvGTxvyvOBr/qIuzMH5/5e4NwJGp1fnwnPa+7M4L1mx+vpYOnmeTrH4Rzg36u4ANKFJXqklPp3PG+i72utvwe8h6crCjxvCl3fPPrqKuAtrfXTeLplbsLzZgKeN4aO634IXK2UyvbGsxLPDKlwYC3wRaWUUSkVB9w4wFjOo5T6Ip6uj3/4KBfi7Y+P0Fo/A3wNyPGOKfTkUyBbKTXfe42LgEuBj/oQWuefTU800KKU+pz3+uPwtCbm9nY9rXUxUI0nWbyP53e9CkjQWu/Hk+TbvWM+eLsiV+FJNp0dBuxKqeXecvPxtEq6W7V1HfAZpVSSt+wXgBqgEKjCM+bT8VxLu9S9y3tuPJ5Et9bH36vwA0kgAs4fA+no4/8rnjf2fKXUbiAGT582wD+Bj5VS0wfwfM8AlymlDgB78HR1ZXgH17cBk5RSr3qnBd8LvKiU2g/8CLhBa90EPILnE/gR4C3gwMBeOgC3e1/zXu/zXINnHKK1t0paawfwTWCNUmoPnp/Jv2mtu7bCOtepxjNm8xvv618DfEFrXdCHON8EfqKUuruX67fhSaZfUkrl4UkGP9Bab+mm+Frgq0qp//I+fg1IAvZqrU/gGWN6zXvddjyJ/hve664HHtVab+zy/A48ieURpdRePC2ycs7tDuwo+wGeyRQblFKH8IwHXecd9P4NkKqU0niSwoYu1TO8f5PrgK9rrTW9/70KPzDIcu5CiMGklPo58L9a6wpvC2g/nskXdYN0/ZPAaq31rsG4nhg4GQMRQgy2U3gmPrTjnaI9WMlDDC/SAhFCCDEgMgYihBBiQCSBCCGEGJCgGgPxTp/MxbO0xYWsOSSEEKOJCc9Nuzt7m0XYVVAlEDzJY1OggxBCiBHqEvq+plzQJZAygL///e+kpKQEOhYhhBgRysvLufPOO8H7HtpXwZZAnAApKSmkpw/2iuBCCBH0+tX1L4PoQgghBsSvLRCl1B149kow41kZ87ddzt+MZ1luE7ATuFdr3eZdquFxPPsDALyjtX7In7EKIYToH78lEKXUWOAxPIu42YGtSqmNHdueevd0eAqY413y4EU8q5n+Hs8iag9orV/wV3xCiJHL5XJRXFyMzWYLdCgjitlsJjk5mejo6EG5nj9bIFcCG7TWtQBKqZfx7Fr3KIB357OJWut2pVQEnr0Mznjr5gJZSqn/xrOOzn9orc+c9wxCiFGpuroag8GAUgqjUXri+8LtdtPS0kJJSQnAoCQRf/7k0zh3RL+MLntde5PHCqAISMSzcmhH2R/h2cymCE9LRQghAKirq2PMmDGSPPrBYDAQERHB2LFjqaysHJRr+rMFYuTcPQAMeHaTO4d357IEpdSP8eyjfIfW+uaO80qpn+FZ7lsIIQBwOp2YzQPdjmZ0Cw8Pp739vH3eBsSfCaQYz00pHVL4105yKKXigXla645Wx9+Bl5RSMXj2VPil97iB3nebE6OIy95Ca/ER7KVHcbW1YDCasKROJnzCdEzhUYEOTwwhg6EvW6eLrgbz5+bPBLIez6YySYANzyYz93Y6bwCeV0rN01qfxrPJzmY8W15+Vym11btv9P14N7URo1dbVRH129+iKX8z7nY7YMBgtuB2OsDlBFMIMbkriV28GlOYNdDhilHG4XDwhz/8gTfffBODwYDT6eTmm2/mK1/5ypAnOq01DzzwAO+8847fn8tvCURrXaKUegjYiGc/42e11juUUu8CD2utdyml7gXeVkq5gXzgq1prp1LqNuBppVQ4nj2a7/JXnGJ4czTVceajNTTmbcQQYiZy2hIiL1pCaFoWxtBw3E4H9tJCGvatp37bWzQd3EzKZx4idMzEQIcuRpEf/vCHVFdX89JLLxEdHU1TUxP33XcfUVFRHXd4D4nXX3+dJ554Ysi69/x6H4jWeg2eLTs7H1vZ6fvXgde7qbcJ2ct4VHO73TQd/JiaD/6My95KdO5K4havwhRx7swRgymEsHHZhI3LJmbucspf/imlf/0+Kbd+j/CJMwIUvRhNysvLefPNN/nkk0/OzmyKjIzk4YcfprCwEIAHH3yQuro6Tp06xXe+8x3i4+N57LHHsNvtxMXF8eijjzJhwgQ+//nPc//997NgwQKKi4u566672LBhAw8++CChoaEcOHAAm83Gv//7v3PTTTedE0djYyMffvghv/jFL/je9743JK892JYyEUHA2Wqj+t2nsR3+lNB0RdK1X8OS6HtpmtC0TMbe8zhlL/6I8n/+lLFfeLxP9cTItmHXaT7Ycdov175q/niumDe+1zJ5eXlMnjyZmJiYc45PnjyZyZMnn30cGxvLM888Q1tbG8uXL+dXv/oVOTk5rF27lgceeIBXXnml1+cpKiripZdeoqamhltuuYXFixeTlJR09nxUVBS/+c1vKC4uHsArHRiZAyeGFXvFSUqe/TY2vYP4yz9H2l3/068kEBKdQOpnvo/RHErFPx/H2dLkx2iF8Og8zrFu3TpuvPFGrr/+elatWnX2eE5ODgAnT54kOjr67OMVK1Zw+vRpGhsbe32OW265BbPZTEpKCnPmzGH37t1+eCX9Iy0QMWw0Hd5K1VtPYQyzknbX/xA2dsqArhMSnciYVd+h9Pn/R/W63zPm5gcGOVIxnFwxz3crwZ+mT5/OsWPHaGpqIjIykuXLl7N8+fKzXVAdwsLCAM9d9F253W6cTufZ78EzMN+ZyWQ6+73L5SIkJPBv39ICEQHndruo/egFKl99AkvyRMb+288GnDw6hI3LJm7JKmz5W2g+GvhPaiJ4paWlccMNN/C9732PhoYGwPPm/9FHH3V7o+OkSZOoq6sjLy8PgHfffZe0tDRiY2OJi4s7O26yfv36c+qtXbsWt9tNSUkJeXl5zJ0718+vzLfApzAxqrmd7VS+8Wtsh7cSNXMZicu/jCFkcGaQxF58M035W6ha93vGjf8VxtDwQbmuEF098sgj/PnPf+auu+7C6XRis9lYsGABf/jDH84ra7FY+OUvf8mPfvQjWlpaiImJ4Ze/9Nz29qUvfYkHH3yQV155hWXLlp1Tr7W1lVWrVtHW1sajjz5KXFzckLy23hg6mkvBQCk1ETjx4Ycfyn4gI4Cr3U7lq0/QXLib+Cs+T8zCGwd9znxr8RFK//IQsZfcRvyltw/qtUXgHD58mKlTpwY6jCHz4IMPMn/+fG655ZZBuV7Xn19xcXFHwsrQWp/s63WkC0sEhMveQvlLP6a5cA+JK75C7KKb/HLDVVh6NtbshdRvfxNnc8OgX1+I0UwSiBhyLnszZS88SuvpfJJu+A+i51zt1+eLW/pZ3O1t1G2VBQ3EyPT4448PWutjMEkCEUPK7Wyn4uWfYS8tZMwt3yZqxlK/P6clMZ3IGUtp2LUWR6PsCiDEYJEEIoaM2+2m6q3f0nLyAEnXfQ1r9oIhe+64xatwOx007F47ZM8pRLCTBCKGTP32N2k6tIm4pZ8lKufyIX1uc3wqEVNyadjzHq621iF9biGClSQQMSRaTh2kdsPzWLMXErt4le8KfhC78EZcLU005n0UkOcXIthIAhF+52xpovL1JzHHp5B03f0B28chNF0RmpZFw863cbvPvxtYCNE/kkCE39V88GectjqSb/xmQG/mMxgMROeupL22jNZThwIWhwg+DoeDp59+mhUrVrBy5UquueYannnmGYbyPjubzcY3vvENrr/+eq6//vqRvR+IEAC2o7toOvARsUtWE5o62XcFP7NmL6TmvT/SsPcDWe5dDJrhsB/I73//e9LS0njyySepqanhxhtvZMGCBSQmJvrtOSWBCL9xOdqoef+PmJPGEbdkdaDDAcAYYiEy5zIadq3DaavHZI3xXUkMa415H9G4f4Nfrh018wqici7rtcxw2Q9k/vz5ZGRkAJCQkEBsbCzV1dWSQMTI1LDjbRx1laTc8TAG09DskNYX0bOupGHH2zTmbSR20U2+KwjRi+GyH8jixYvPfv/uu+/S1tZGZmbmIL3K7kkCEX7haDrDmS2vEJGVS0TGzECHcw5L0jhC07NpzNvol/W3xNCKyrnMZyvB37ruB/L000/jcrmwWCxnE0Nv+4E8/PDDA9oPZPny5eeVW7t2LT/+8Y959tln/b7kuwyiC7+o+/R13O1tJFw5PLezj5qxlPbqYtrKTwQ6FDHCdd4PBGD58uW88cYbPP3005w586+VD4ZiP5C//e1v/PSnP+WPf/wj2dnZF/jKfJMEIgado6mOxj3vEzljKeb4tECH0y3r1IvBFELjwY8DHYoY4YbLfiDr16/nueee44UXXkAp5Y+Xeh7pwhKDrn77m7idDuIWD+3ib61tDlpaHThdbtxuiIowE2oxddtFZQqPJCJzLrZDm0lYdhcGo6mbKwrRN8NhP5Bf//rX2O12vvrVr5499j//8z/MmOG/2YZ+3Q9EKXUH8H3ADPxKa/3bLudvBn4ImICdwL1a6zal1HjgeSAZ0MCdWmufm1vLfiCB52xp5PRvvopVzSf5xm/45zmcLgpO15F/ooYTpQ2cKm+gqq4FW0v7eWUtIUZSEq1kpMaQOS6WeVOTSU+OAsCmt1Px8s9Iuf0hIjLn+CVW4R+yH8iFGaz9QPzWAlFKjQUeA+YCdmCrUmqj1jrfe94KPAXM0VpXKKVeBO4Bfg/8Dvid1vpFpdQPgB8A3/NXrGLwNO7fgLu9ddBnN7U7nOw6XMnHe4rZW1BJc6unfzgxNpyMtGimT0ogPiYMa7gZk9EIuGlqbqeuyU5JVROHTtTw8d5i/vgmpCVaWZY7nqtzp2MMi6Qpf4skECEGwJ9dWFcCG7TWtQBKqZeB1cCjAFprm1Jqota6XSkVgae1cUYpZQYuBTregZ4DPkYSyLDndjlp2L2OsPHTsCRPGJRr1tS38Nam47y37RRNLe3ERoZyyayxzJqSxIzJicREhvb5WpW1zezML2frgTL+tvYwL36g+Wb6ZMbpHbid7cNqqrEQnT3++OOBDqFb/kwgaUBZp8dlwPzOBbzJYwWe7qoS4H0gEWjQWjs61ZP+qBGg+dheHHWVxF/x+Qu+Vn2Tnb+/d4QPtp/C5XKzcEYqVy+YwKysJEymgc39SI6P4Nolk7h2ySROlzfw1uYTrN1dxJcj97P+jfe44saVA762GHput1umYA/AYA5b+DOBePoR/sUAnDd/TWu9FkhQSv0YeBr4Tpd6dFdPDD8Nu9ZiiozHOmW+78I9cDpdvLX5OC+8r2ltc3LNggncfFkmqYnWQYwUxqdEc9/qmRQtHo/tz5up2vcJ/10aw3/eOZfkuIhBfS4x+EwmE+3t7VgslkCHMuK0tLRgNg9Oa9ufH7eKgdROj1OA0o4HSql4pVTnvUz/DuQAlUCMUqpjWkxq53pieHI0VNNyfD9Rs5ZhMA3sc0l5jY3/+t0W/vjmIaZOjOepb1/O11bPHPTk0dm41Djipy8iN7KMU6V1fOOJj9h1uMJvzycGR2xsLBUVFd3eUyG653a7aW5upqSkhOTk5EG5pj9bIOuBR5RSSYANWAXc2+m8AXheKTVPa30auBXY7O3W2gTcDqwB7gJkG7lhrungJ4B7wHcEb95fwq9f2ovRYOA/75zLZXOGrtfSOnURTQc+4n9XJ/Lzj+z86E/b+Y9bZ3Hl/PFDFoPon8TERIqLi9FaBzqUEcVsNjNmzJiza3ZdKL8lEK11iVLqIWAjYAGe1VrvUEq9Czystd6llLoXeFsp5QbygY4JzF8D/qKU+j5wGvisv+IUF87tdtN44GNC07Mxx6X0u+4/Pizg+bVHmDoxnm9/bui7kCIyZmIIjSC0dC8/+dqX+clfdvLkS3upb7Kz6oqsIY1F9I3RaGT8eEnwgebXGwm11mvwtCI6H1vZ6fvXgde7qXcKuMyfsYnB01Z+nPbqYhJXfKVf9ZwuN7/95z4+2HGay+am8/XbZmEOGfob+gwhZqxZ87AVbCdxxb08/MWF/OqFPTz3Tj6hFhPXLZk05DEJMRLIlBNxwRoPfAymEM/yIH3kcrl56h+e5HH7VVN44LNzApI8OlizF+FqaaLl1CHMIUYeuGMOCy5K4fevH2Dz/pKAxSXEcCYJRFwQt9uF7fCnREyegyk8so913Dz9ah7rd57ms1crPrd8asCnY4ZPnoXBEobt8FYATCYj3/n8PKZOjOeJv+/h0PGagMYnxHAkCURcEHvJUZxNtVinLupznRc/KGDdpydZfUUWn716aBZ988UYYiEiax62gh24XZ5VUUPNJn7wbwtIjgvnZ3/bSV2jPcBRCjG8SAIRF8Smt4ExhIjMub4LA58eKGXNe0e4fG46d60MfMujs8jsRbiaG2gtOvyvYxEWHrw7l6bmdp74+26crqHb41qI4U4SiBgwt9uN7cg2wjNmYArzfa/GqbIGfrFmD1PGx3L/rbOGVfIACJ80C0OIBZvefs7xjLQYvnpLDvuOVvGPD2TaqBAdJIGIAWurOImjrhKrWuizrL3dyeN/3UlEWAj/fc98LObht3y60RJGeMZMmvWO85Z7uHL+eC6bm85L6ws4UVofoAiFGF4kgYgBs+ntYDBinZLrs+xf3smnuLKJb312Dgkx4UMQ3cBY1XwcDdXn7VRoMBi496YZRFkt/PqlvTidcge0EJJAxIC1HNtD6NgpmKwxvZbbV1DJW5uOc92SDGZNGZwlFPwlImseGIzndWMBREVY+MrNMygsrueNT44FIDohhhdJIGJAHE112MuOETF5dq/lmlvbefLFvYxNiuTua6cNUXQDZ4qIJmz8VGwFO7o9vzgnjQUXpfD3dUeoqG0e4uiEGF4kgYgBaTm+F8DnRkwvfVBAdX0r3/zMbMIsI2MHZeuU+bRXnaa9tuy8cwaDga/ekgMGA399Nz8A0QkxfEgCEQPSfGwvpsg4LGMyeixTVNHIG58c46r548meGD+E0V2YCOVZjr67bizw7IJ489LJfLK3hILTZ4YyNCGGFUkgot/cLictx/cRMXl2j1Nx3W43v3/9AGEWE3etHP5dV52ZY5KxjMnosRsL4JbLM4mNDOVPbx0a1A16hBhJJIGIfrOXFOBqtRE+uefuq20Hy9lXUMWdy6cSG9X3bWeHC6tagL24AEdT9y2MiDAzdyzP5tDxGrYdPL+rS4jRQBKI6LfmE3lgMBI+cUa3550uN39be5j05EhWXjxxaIMbJFa1AHDTXLCzxzJXzx/P2KRI1rynpRUiRiVJIKLfWk8eIDRlUo+LJ27aV0JRRSN3XJM9YvcYNyeNIyQuBZvuuRvLZDJy25VTOFnWwM582cVQjD4j83+3CBhXWwutJQWEZ/TQ+nC6WPPeESamRrM4J22Ioxs8BoMBq5pPy8kDuFptPZZbOnssY+IjeGm9tELE6CMJRPRL6+l8cDkJ66H7asOuIsqqbdy5PBujcXitddVf1ikLwOWg+djeHsuYTEZWX5FFwek69hVUDWF0QgSeJBDRLy0nD2AwmQlLzz7vnNPp4qX1BWSOi2XBRf3b2nY4Ch2bhcka2+N03g7LcseREBPGS+sLhigyIYYHSSCiX1pOHCB0XDZG8/kzq7YeKKOitpnblmUNu5V2B8JgNBGRNY/mY3twOdp6LGcOMXHT0kwOHa+hsKhuCCMUIrAkgYg+czY30lZ5stvZV263m9c+KiQ10cr8i1IDEJ1/WNUC3G2ttJ480Gu5q+aPJ8xi4q3Nx4coMiECTxKI6LPW4iMAhI+/6Lxz+SdqOVpUx01LJ2Ma4WMfnYVPnIHBEt7rbCwAa7iZy+eNY9O+EuqbZOdCMTr4dXEipdQdwPcBM/ArrfVvu5y/EfghYABOAF/QWp9RSt0NPA50zI18R2v9kD9jFb61Fh3GYDITmjr5vHOvfVRIVISFK+aNC0Bk/mMIMROROYfmoztxu+7FYOx5H5PrFmewdutJ3tt2ituunDKEUQoRGH5rgSilxgKPAUuAWcC9Sqlpnc5HA08D12qtZwJ5wCPe0/OAB7TWs7xfkjyGgdaiI4SmZWIIMZ9zvKSqie2Hyrl2ccaIWTCxP6xqAU5bPfaS3gfJx6dEMzMrkbVbT8h+IWJU8GcX1pXABq11rdbaBrwMrO503gzcp7Uu8T7OA8Z7v88F7lZKHVBKPa+UivNjnKIPXO127GXHCBt3/uyrtVtPYjIaRuxd575ETJ4NphCfs7EArl8yier6VrYfKh+CyIQILH8mkDSg8yJBZUB6xwOtdY3W+jUApVQ48CDweqeyPwJygCLgKT/GKfrAXlYILgehXabvtrU72bDrNAtnpBIXHRag6PzLGBpB+MQZ2LrZ6raredNSSIgJ44Mdp4coOiECx58JxAh0/t9mAM5r1yulYoB3gP1a678AaK1v1lpv0Vq7gZ8BK/wYp+iD1iLPAHpYujrn+Ja8Uhqb21mxcGIAoho61inzcdRV0FZ5qtdyJqOBK+aNY8+RCmrqW4YoOiECw58JpBjoPJ8zBSjtXEAplQpswtN99SXvsRil1Lc6FTMADj/GKfqgtegw5qRxmMKjzjm+dutJ0hKtzMhMDFBkQ8OqFni2uj281WfZZbnjcblh4+7iIYhMiMDxZwJZDyxTSiUppSKAVcC6jpNKKRPwFvAPrfU3va0NgCbgu0qpBd7H9wOv+TFO4YPb5cRerAlLn3rO8VNlDRw+Wcs1CyeO+GVLfDFZYwifOIOm/C0+u7HGJkUyLSOe9TtOy/pYIqj5LYF4B8cfAjYC+4A1WusdSql3lVLzgBuAOcBqpdQ+79ezWmsncBvwtFLqMDAX+K6/4hS+tVUV4bI3nzeAvm7bSUJMRpblBtfU3Z5Ypy3GcaYce9kxn2WvzB1PSVUT+pTsWCiCl1/nXGqt1wBruhxb6f12Fz0kMK31JjzJRQwDrUWHAQgb968WSLvDxcd7SlgwPYWYyJG3YdRAWNUCqtf+Hlv+ZsLSMnstu3hmGv/3+gE+2HF6RG3nK0R/yJ3owqfWosOYouIJiUk6e2zPkQoam9uC7sbB3pjCI4mYPMvbjdX7fR4RYWYunpHKlrxS2h1yT4gITpJARK/cbjetRYcJGzf1nAUSN+4uJibSwhyVHMDohl7ktCU4G2vPzkqAwkOWAAAgAElEQVTrzaWz07G1tLO3oHIIIhNi6EkCEb1yNFThbKw9Z/n2puY2th8qZ+nsdEJG6I6DAxUxZR6GEAu2/C0+y87MSiIqwsymvSU+ywoxEo2u//2i387e/9Fp/GPz/lIcTheXzx093VcdjJZwIrLm0nR4K26Xs9ey5hAjF+eksf1QGa1tMhNdBB9JIKJXrUWHMVjCsSSPP3tsw64ixo2JYnJ6TAAjC5zIaZfgam6g5eRBn2UvmTmWFruT3YelG0sEH0kgolf2kqOEjc06uwptZW0zh0/Wcvnc9KDYNGogwjNne5Z4z9/ss+z0zERio0LZtE+6sUTwkQQieuRqt9NWeYrQ1H9NWd2837OYwCWzxgYqrIAzhliwqgU0HdmGq733vT9MRgNLctLYmV9Oc2v7EEUoxNCQBCJ61FZxEtwuQtM6J5ASMsfFkpJgDVxgw0BUzmW47c00+9hoCjz3hLQ5XOw+It1YIrhIAhE9spcVApxtgVTUNnO0qI5LZqYFMqxhIWzCRYTEJNOYt8Fn2akZCcREWth2sMxnWSFGEkkgokf20kJMkXGERCcAsGW/px9/8czR233VwWAwEpVzOS0nDuCor+q1rMloYP60FHYdrpCbCkVQkQQiemQvKzxn/GPT/lKyxsUyJj4igFENH5E5lwFuGvM2+iy7cEYqza0ODhRW+z0uIYaKJBDRLVerjfaa0rPjH+U1NgqL6lgirY+zzLHJhGfk0LjvQ5/3hMzKSiLMYpJuLBFUJIGIbtnLjwMQmjoZgK15ntlXi2X84xxRc67G0VBNc+GeXstZzCbmZo9h+6EyXC5Z4l0EB0kgolv20nMH0LcdLGdyeox0X3VhzcrFFBlHw573fJZdOCOV2gY7BUWyxLsIDpJARLfsZYWExI7BFBHFmYZWjpyqZeH0VN8VRxmDKYSoWVfScmwf7XUVvZadN3UMJqOBbQekG0sEB0kgolv20sKz4x878stxu5EE0oPo2VeBwUDDrnW9losMN3PRpAR2He490QgxUkgCEedx2upxNFSf032VkhDBhJQoHzVHp5DoBKxTF9Gwbz0ue3OvZedNHcOp8kYqz/ReToiRQBKIOM/Z8Y+0yTS3trOvoIqF01NH7dpXfRGz4Abc9mYa9n3Ya7l5U8cAsFtaISIISAIR52ktKwSDkdCUSezRlTicLum+8iEsLZOwcVNp2PlOr1N605MjSUmIYKckEBEEJIGI89hLj2JOTMdoCWfbgXKirRbZ17sPYhbeiKO+iqZeNpsyGAzMyx7D/qPV2Nt7v3dEiOFOEog4h9vtxl52jNDUTJwuN7uPVJA7zTN7SPQuImsuluTx1G15pddWyLxpY2hrd3LwmNyVLka2EH9eXCl1B/B9wAz8Smv92y7nbwR+CBiAE8AXtNZnlFLjgeeBZEADd2qtm/wZq/Bw1Ffham4gLG0yR4vO0NTSzlw1JtBhjQgGg5HYJbdS+eoT2I5sI3La4m7LTZ+ciMVsYld+BXOz5WcrRi6/tUCUUmOBx4AlwCzgXqXUtE7no4GngWu11jOBPOAR7+nfAb/TWmcDu4Af+CtOca7OK/DuPVKJwQAzpyQFOKqRw5q9EHNiOmc2/xO3u/uFE0PNJmZmJbLzcAVut9yVLkYuf3ZhXQls0FrXaq1twMvA6k7nzcB9WuuOrdrygPFKKTNwqbc8wHPArX6MU3RiLy0EUwiWMRPYrSuZMi6OaKsl0GGNGAaDkbhLbqO9qoimg5t6LDdXJVNR20xZjW0IoxNicPkzgaQBnW+5LQPSOx5orWu01q8BKKXCgQeB14FEoEFr7eiunvAve1khockTabK7OXr6DLNVcqBDGnGsUxdhSZnEmY9fwO3ofhfCWd6f676C3peCF2I482cCMQKd2+cG4Lw2vVIqBngH2K+1/ks39eiunhh8brcLe9lxQtMy2VdQhcsNc7MlgfSXwWAk/orP4aivon5393enpyVaSYoLlwQiRrReE4hSas4FXLsY6HzzQApQ2uX6qcAmPN1XX/IergRilFIm7+PUrvWEf7TXlOJuayE0dTJ7dSXWcDNZ42IDHdaIFJExk/BJM6nb/DLO5obzzhsMBmZlJZF3tAqnrM4rRihfLZBnO75RSn2/n9deDyxTSiUppSKAVcDZj2PeBPEW8A+t9Te11m4ArXU7nqRyu7foXcDafj63GICOO9AtqZns0ZXMykrCZJKZ3gOVcOU9uOzN1H70QrfnZ09JxtbqoFBW5xUjlK93h86T/2/pz4W9g+MPARuBfcAarfUOpdS7Sql5wA3AHGC1Umqf96sjYX0Nz6ytfOASPFOBhZ/ZywoxWMIoa4+mpr6VOdJ9dUEsSeOJzl1J494PsJcdP+98TlYiIOMgYuTydR9I1zGMftFarwHWdDm20vvtLnpIYFrrU8Bl/X0+cWHspYWEpkxm+1HPDW5zZAD9gsVfchu2Q5uoXvsMaff8BIPRdPZcTGQok8bGsLegituvUgGMUoiB6U//hHTUBjG3s522ipOEpk1mz5FKxo2JIjE2PNBhjXjGMCsJV38Re9kx6re/dd752VOS0KdqabE7uqktxPDmqwWSrpT6dTffA6C1/rp/whJDra2yCLezHWPSJA69X8O1izMCHVLQsE69mIj8LZz5+EUipuRiSfjXvvKzpiTxysZCDh6rJndaSgCjFKL/fLVAfgvUeL86f9/xJYKEvfQoACfa4ml3uOT+j0FkMBhIvObLGMyhVL39u3PuUJ+WkYAlxCjjIGJE6rUForX+4VAFIgLLXlaIMTyKXUVOLGYT0yclBDqkoBISFUfCVV+g6q3f0LBrHTG5nqFAi9nEtEkJ7JUEIkYgn4spKqVygW8BM4Bm4ACehREP+jk2MYTsZYWEpmayR1cxfXICFrPJdyXRL5EzltKUv5najc8TkTkHc5yny2r2lCT+/HY+NfUtJMTIuJMYOXzdSLgMeBNP0vgensUOTwPvK6WW+j06MSRcba20VRXjiJtASVWTzL7yE4PBQNLKr2Iwmqh889dnl3yfNUWWNREjk68WyH8B12it8zodW6uUWgf8BFjmt8jEkGmrOAFuFyfa4wGHJBA/ColOJGH5l6l640nqtr5G3JLVTEyNJibSwr6CKpbljg90iEL0ma9B9OQuyQMArfUOQNa4CBKt3jvQd1WGkxQXTnpyZIAjCm6RF12Cddpizmz6B/bSQoxGAzMzk9h3tEqWdxcjiq8E0tuem7JFXZCwlxViikpg27EW5qhkDAb51fqTwWAgcfm9mKwxVL75JK52O7OmJFHXaOdUeWOgwxOiz3wlEPk4NArYSwtpixlPi126r4aKKTyS5Ov/g/aaUmo3/K3TOEhlgCMTou98jYFkK6XO68LC0/qY5Id4xBBztjThOFNOUdgMT1dKluw+OFTCM3KInn8dDTveJiVzLmOTItlbUMVNSzMDHZoQfeIrgazAM323DAjDs0S7CCL2smMA7K21kj0hDmu4OcARjS7xl99Jy4n9VL39W3InfYl1e6tpd7gwh8gqyGL48/VXOgl4FPgu8CQQqrX+uOPL79EJv+vYA317eZh0XwWAMcRC8g3fwNncyKKmD2htc3BUlncXI4SvBPJ1YLrWegFwPZ5tZ0UQsZcW0h6RRIvbIsuXBEhoSgbxSz9DWNk+5oceJ6+wOtAhCdEnPtvJWutS77+fAtJBHmTspYWUG5KJtlrITJeZ2YESs/AGwsZNZbV1B8eOFAY6HCH6pL+zsGTN6SDiaKzF2VTLwYZoZk1JwmiU6buBYjCaSLrh65iMRubXvkNLa1ugQxLCp/6O1Mm03iDSsYXtEVuMjH8MA+bYZFpn387kkApOvP9SoMMRwidfs7BylFINnR5HeB8bALfWOtp/oQl/s5cV4sZIiSNexj+GiclLV/Depx+Tc/At2hZeiiV5QqBDEqJHvlogk/FM4+346ng83fuvGMHsZYXUmhJIT4snPjos0OEIwBpuYW/8cuxuM1XvPnPO3iFCDDe+9gM5NVSBiKHldrtpLT3G0eZU5syX1sdwoqaM45VN8/hcyWYadr9PzLzlgQ5JiG7J3UqjlKOuAndrEycdCczJlgQynMzMSmKnPYO2pGxqNz6Po0E2/xTDk18TiFLqDqVUvlLqqFLqvl7K/VUpdU+nx3crpcqUUvu8X4/5M87RqGMAvdyQzNSJsvvgcKImxGEJMbEr9mpwOal+/4+BDkmIbvnckXCglFJjgceAuYAd2KqU2qi1zu9UJg34Pzz7imzoVH0e8IDW+gV/xTfatZYW0o6J5IwsWTZjmLGYTUzNiGf76Tauu+RWajf+HZvejlUtCHRoQpzDn+8cVwIbtNa1Wmsb8DKwukuZO4E3gH90OZ4L3K2UOqCUel4pFefHOEelxtOa4vZ4ZmenBDoU0Y2czCROljXgnnYNluTxVL//J1zt9kCHJcQ5/JlA0vAswtihDEjvXEBr/XOt9bPd1C0DfgTkAEXAU/4KcjRyu5w4K09y2pnAbBn/GJZmZiUCcPBkHQnXfAlnQzX1n74R4KiEOJc/E4iRc288NAB9mpOotb5Za71Fa+0GfoZnVWAxSNqrSzC62mgIG0taouw+OBxlpscSHhpCXmE14eMvwjr1Yuo+fQ1HveybLoYPfyaQYiC10+MUoNRXJaVUjFLqW50OGZAlVAaVrbgAgNgMFeBIRE9MJiPTJyeQd9STMBKW3QVAzYa/BTIsIc7hzwSyHlimlEpSSkUAq4B1fajXBHxXqbMjhvcDr/kpxlGpQh+k2WVmas7UQIciepGTmURptY2qMy2ExCQRs+gmbPlbaDmd77uyEEPAbwlEa10CPARsBPYBa7TWO5RS7yql5vVSzwncBjytlDqMZxbXd/0V52jUVnqUImcSM7Jk/GM46xgHySv0tEJiF92EKTqRmvf/hNvlDGRoQgB+nMYLoLVeA6zpcmxlN+Xu6fJ4EzDHn7GNVq62ViJaK7BFzSc81K+/fnGBJqREE221kFdYzbLc8RjNoSQsu4vK135BY95HRM9aFugQxSgnNwCMMlWF+RhxEz1Ruq+GO6PRwIzMRPKOVuF2e+ajWKdeTOjYKZz55EWZ1isCThLIKHP64H4AMmfNDnAkoi9mZiZSXd9KWbUNAIPBQPzln8PZWEvDrrUBjk6MdpJARpnW4gJq3NFMzBgb6FBEH8zM8mwCur/TNrfhEy4iInMudVtfxdnSGKjQhJAEMpo4nC6imotpjhovuw+OEKmJVhJjwth/9Nz7P+IvvxNXazN1W2WCoggcSSCjyPGCE0QbmrGOl/s/RgqDwUBOVhIHCqtxuf51X64leQKRM5bSsPNdHA3VvVxBCP+RBDKKnMjbC0DGjFkBjkT0R05mIg22Nk6VN5xzPG7p7bhxU/uxbH8rAkMSyChiO61xYCI+IyvQoYh+yMn0jIPkFZ7b0jDHJBMzdzlNBz6iraYkEKGJUU4SyCjR1NJOdHMxzdY0DCZzoMMR/ZAUF05aopW8o+d3VcVefAuGEDN1m/4ZgMjEaCcJZJTYr8tJN9UQkS7jHyNRTlYSB49X43Seux6pyRpD9NzlNB3aTFt1cYCiE6OVJJBRojDvIBaDkzHZ0wMdihiAnMxEmlsdFBbXnXcuduGNGMyhnNnUdVsdIfxLEsgo4Ha7aTx1BICIcdICGYlyMjvWxTq/G8tkjSEmdwW2/K20VZ0e6tDEKCYJZBQoqmgksb0UhyWKkOikQIcjBiAmMpSJqdHdjoMAxCy4EYNFWiFiaEkCGQV2Ha5kQkg1oWOzMBjkBsKRKicrkfwTNbQ7zl+J1xQRRcy8ldgOf0pb5akARCdGI0kgo8CBQycZY2ogeoIsoDiSzcxMos3h4sjJM92ej1l4AwZLuLRCxJCRBBLkbC3ttJV4xj/CxmUHOBpxIaZPTsBoNLC/sPttbU3hUcTMvxbbkW3Yy08McXRiNJIEEuT2FlQy0VSB2xhCaFpmoMMRFyAizExWemyP4yAAMfOvxxgaIa0QMSQkgQS5nfkVZIVWEZaWiTHEEuhwxAXKyUqk4PQZWuyObs+bwiOJmX89zQU7sJcfH+LoxGgjCSSIuVxuDhwpId1YQ/h4Gf8IBjmZiThdbg4dr+mxTMz8azGGRlC35ZUhjEyMRpJAglhhcR2xrSUYcRE2ThJIMJiWkYDFbGKvruyxjDHMSvS8FdiObKOtqmgIoxOjjSSQILYzv4LJ5krAQGi6DKAHA4vZxPTJCew+0nMCAYiZfx0Gcyh1W18dosjEaCQJJIjtOlzO9MhaLMkTMIVZAx2OGCRzVTIlVU1U1Db3WMYUEU30nKtpOrSZ9jPlQxidGE38mkCUUncopfKVUkeVUvf1Uu6vSql7Oj0er5T6RCl1RCn1hlIq0p9xBqPahlaOF58hzV1BmIx/BJXZKhmAPb10YwHELLgBjEbqPn19KMISo5DfEohSaizwGLAEmAXcq5Sa1qVMmlLqLWB1l+q/A36ntc4GdgE/8FecwWr34QrSTbWYXG0y/hFk0pMjSY4LZ8+Ril7LhUTFEzXzChr3b8TR0POguxAD5c8WyJXABq11rdbaBrzM+YniTuAN4OykdaWUGbjUWx7gOeBWP8YZlHYermBGVC2AJJAgYzAYmK2S2X+0GkeX5d27il10M7hd1G17Y4iiE6OJPxNIGlDW6XEZkN65gNb651rrZ7vUSwQatNaOnuqJ3rU7XOwrqCIn6gwhsWMIiYoPdEhikM3NTqbF7uDIydpey5ljk4mccSmNez/AaasfoujEaOHPBGIE3J0eG4DePy51X48+1hNe+SdqaLG3k9xeIuMfQSonMwmj0eBzHAQ8uxa6He3U73h7CCITo4k/E0gxkNrpcQpQ2od6lUCMUsrkfZzax3rCa9fhCtLMjRjbmqT7KkhZw81MnRjPrsO9j4MAWBLGYp26kPpda3G2NA1BdGK08GcCWQ8sU0olKaUigFXAOl+VtNbtwCbgdu+hu4C1fosyCO3Mr2BJqg2AsHHTfJQWI9X8aWM4UdpA5Zmep/N2iF28GndbCw275L+SGDx+SyBa6xLgIWAjsA9Yo7XeoZR6Vyk1z0f1r+GZtZUPXAJ8319xBpviykZKqpq4KKIGkzUGc3yq70piRMqdlgJ4PjD4EjpmIhGZc6nf+TauthZ/hyZGiRB/XlxrvQZY0+XYym7K3dPl8SngMn/GFqw+PVAGuIlrOk7YxOmygVQQS0+OJDXRyo78cq5dnOGzfOyS1ZQ+91807PmA2IU3DEGEItjJnehBZtvBMnLHunE31xE+cUagwxF+ZDAYmD8thbyj1T2uzttZ2NgphE2cQf22N3A52oYgQhHsJIEEkeq6FgpO13HpmAYASSCjwPyLxuBwuthX4Hs2FkDc4lU4bXU07tvg58jEaCAJJIhsO+i57WaCu4SQmGTMcSkBjkj427SMBKxhIew45HscBCBswnRCxyrqP30Nt9N3q0WI3kgCCSKfHihjfHIEhoojhE+cHuhwxBAIMRmZmz2GXYcrcLm63j51PoPBQNySVTgaqmk6+MkQRCiCmSSQIFHfZOfg8RqumAyuVhvhE3MCHZIYIvMvSqGuyc5hH3eldwifPAfLmAzObHkFt8vp5+hEMJMEEiR25pfjcrmZaa0CIExaIKNG7rQxmEOMbM3r2/22nlbIahxnyrEd3urn6EQwkwQSJDbtK2VMfAQR1UewjMkgJDIu0CGJIRIRZmaOSmZLXmmfurEAItR8zInpnlaIW1YKEgMjCSQI1DfZ2Xe0istzEmktPkLE5FmBDkkMsSUz06ipb0WfOtOn8gaDkdjFq2ivKqK5YKefoxPBShJIENjq/eS5KLEeXE7CJ80OdEhiiOVOSyHEZGRzXkmf60ROW0xIXApnNr+C2923losQnUkCCQKf7CshPTmSqDNHMFjCCUufEuiQxBCzhnu6sbbu73s3lsFoIvbim2krP0bL8X1+jlAEI0kgI1xNfQuHjtdw6cw0Wo7vI3ziDAwmc6DDEgGweGYa1fWtFJzuWzcWQNSMpZiiEzmz+WVphYh+kwQywm3eX4rbDYsnGnHUVxExScY/RqsFF3m6sT7Z1/duLIPJTOyim7AXH6H1dL4foxPBSBLICPfJ3mIy0qKJPnMYgIjMOQGOSASKNdzMgukpfLynmHZH32dWRc28ApM1lrotL/suLEQnkkBGsKKKRgpO13H53HHYCnZ6pu/GJAU6LBFAV+aOp8HW1qeNpjoYzaHELLyBlhN5tJYU+DE6EWwkgYxgH+48jdFo4NKp0diLNRFTcgMdkgiw2VOSiIsK5cOdp/tVL3rO1RjDI6nbLK0Q0XeSQEYop8vNxt3FzM1OxlJ+EHBjzZIEMtqZTEYunzuOXYcrqGu097me0RJOzPzraS7cjb38hB8jFMFEEsgIta+gktqGVpbljsdWsBNTVAKWFN+bCongd0XuOJwuNx/vLe5Xveh5KzCERshYiOgzSSAj1Ic7i4iKMDMvK5aWE/uxZs2T3QcFABNSoskaF8v6Haf7NTXXFGYlJncltiPbpBUi+kQSyAjU1NzGtoNlLJ2dTvup/bjb7VinLgp0WGIYuWbhBE6WNXDoeE2/6sUsuAFjmJUzn7zkp8hEMJEEMgJt2FVEu8PFlfPHY8vfgskaQ9j4aYEOSwwjl80dR1SEhTc3He9XPVOYlZgFN9B8dCetJUf9FJ0IFpJARhiXy827W0+gJsSRkRxG89HdWNVCDEZToEMTw0io2cTyRRPYdrCM8hpbv+rG5F6LMTyKM5+86KfoRLDwawJRSt2hlMpXSh1VSt3XzflZSqldSqkCpdSzSqkQ7/G7lVJlSql93q/H/BnnSJJXWEVJlY1rF2fQXLgbt6MN67SLAx2WGIauXZyB0WDg7c39G88whoYTu+gmWo7vo7XoiJ+iE8HAbwlEKTUWeAxYAswC7lVKde1neR64X2s9BTAAX/Yenwc8oLWe5f16yF9xjjTvbDlBtNXC4pw0bIe3YrLGEjZuaqDDEsNQQkw4S2aO5f3tp2hube9X3eh5KzBZY6n9+AU/RSeCgT9bIFcCG7TWtVprG/AysLrjpFJqAhCutd7mPfQccKv3+1zgbqXUAaXU80op2R0JqDrTwo5D5Vy9YAImRwu2o7uwTrtYuq9Ej25cOokWu4N3tvSzFWIOJfbim2k9dZCWE3l+ik6MdP5MIGlAWafHZUB6H8+XAT8CcoAi4Cn/hTlyrP30BG5g+aKJ2PK3gNNB1IzLAx2WGMayxsWRO20Mr2wspKm5rV91o+ZcTUhMEjUf/lV2LRTdCvHjtY1A50noBsDVl/Na65s7DiqlfgYc81+YI0Nzazvvbj3JwumpjImPoOSNjViSx8vNg8Knz6+Yyjd+8RGvbCzk7mv7PlvPGGIh/rI7qXzjVzQd/ISoGZcNalztDicnyxo4XlJPSZWNpuY2bK3tmE0mIsJDSIgOY0JqNBlpMSTHhct9TsOQPxNIMXBJp8cpQGmX86ldzyulYoB/01r/0nvcADj8GOeI8N62U9ha2ll9RRZt1cXYS48Sv+xu+U8lfMpIi+HSWem8uek4118yifjosD7XtV60mNAdb1G7cQ3W7EUYzaEXFEtrm4NPD5Tx6YEy9uhK7G1OACwhRiIjLFjDQ2h3uLC1OGjs1GJKTbAyd2oyF+ekMX1SgvzdDxP+TCDrgUeUUkmADVgF3NtxUmt9SinVqpRarLXeAnweWAs0Ad9VSm3VWm8H7gde82Ocw167w8nrHx8jJzORKePjqNnwNhiMRE6/xHdlIYA7l2ezeX8Ja947wv239n3PGIPBSPyyuyl7/mHqd7xD3OJbBvT8FbXNvL35OB/sOI2tpZ346DCumDeOmZlJTBobw5j4CIzGc5NCc2s7pysaKSyqY4+u5P3tp3l78wlSEiK4av4EVlw8kagIy4DiEYPDbwlEa12ilHoI2AhYgGe11juUUu8CD2utdwF3An9QSkUDe4Bfa62dSqnbgKeVUuFAAXCXv+IcCTbuLqa2oZVvfmY2Lkcbjfs3EJE1j5BImVsg+iY10cp1SybxxifHuHT2WHIy+77sf/iEi4iYkkvd1leJnrUMkzWmz3Wr61r4x/oC3t9+CoCLc9JYefFEpmUknJcwuooIM5M9IZ7sCfFct2QSrW0Oth0o44Mdp/nb2sP888MCrl44gVWXZ/WrVSUGjyGYtrFUSk0ETnz44Yekp6f7Kj4iOJ0u7vv5BsJCQ/jlN5fSdPBjqt78DSl3PExExsxAhydGkNY2B19/4iOcLjdPfftywkP7/vmxraaE4v/7JlGzryRpxVd8lre3O3llw1Fe2XAUl9vN1QsmcOuyKSTGhl/ISzjrZFkDr2w8yid7SwgxGbnx0knccnkWkeGynfNAFBcXs2zZMoAMrfXJvtaTO9GHuQ27iiipsnHbsikYDAYadq3DnJBG+MScQIcmRpgwSwjfuH02VWea+cs7/du+1pIwlui519C4dz32sp6XR3G73Xx6oJSv/fRDXnhfs3B6Kv/34JX8+6qZg5Y8ACamRvOfd8zlme8tY9H0VP754VG+8pP1/V5AUlwYSSDDmL3dyZr3jqDGx7FoRir20kLspUeJnrtCBhHFgFw0KYHrL5nEO1tOsHF3Ub/qxi39LKaIaKrefQa3y3ne+aKKRv7f7z/lx8/tJDw0hB9/bTHf+fw8kuMjBiv886QmWvn25+byq28tJS3RypMv7eXB327mVFmD355T/IskkGHsnc0nqK5v5e5rp2EwGKjb9gYGSzhRM5YGOjQxgt1z7TRyMhN58sW97NGVfa5nCrOScPW/8f/bu/PoqMo0j+PfW1v2fSUBggnwQtMxEAQMMewoboiDgooLM+PWDm4t2G492M4Zt0ZtsRUX3HrsVltFEHADZF9FgbDIGxKWAAkhIQGyV5FU/1FFD+0RSUpSt0Kezzk5VN2q4v6KE+q57633fa7zUFeOYVUAAA56SURBVBHHv/vyn9tr6py8MW8rd89YSkFxFbePy+TF3w4jMyO+LeL/pIzO0TwzJY97JvRlf1kN9zy/jLfmb6e+scNP4GxTUkACVE29i4+WFJDdK5HM7vE4j5RQ+8Naoi4YgyU4zOx4oh2z26w8MnkgXZIiePrdDezY0/KW72G9BxOS3o/KZX+j8WgFX6zZw+1PLWH+yt2MGtiVVx8axZV56Vit/v9osVgMRg9K49WHRjJqQFc+XVbIlBlLyS8s93uWjkIKSIB6/6ud1Da4uOUyz8KvY2vnYtjsRA64wuRk4lwQFmLn8dsuJDo8mEdeWc3c5UUt+u7AMAxiL7mV5qYmls16hlc+ySetUwR/un8YU67tS3TEL1sncjZEhjm4e0Jfnv6vi7BaDB6dtYbX5uTTIKORs04KSAAq3H+UBat2c2lON9JTozhxrJzqrcuJyBqBLTza7HjiHBEXFcLz9w3hgt5JvPnZNh5/Yx1bCytOW0iO1TTy5dq93D97B/OrM+nRvJvpFzt48je5pKe2fGqvv/RJj2PmA8MYm5fOgtV7uOe5Za2+wJb4eW25kFD4oKnZzcsfbyYqPIibvKOPyhUfYBgG0TnjTE4nzjXhoQ4e/feBfLZyNx8u0jwyazVdkiJQXWNIjgsFA44cbWBv6XF27qvE7YbzUiLJGjcJ++ZK4nf8nabcnIA9sAl22LhtXCYXZnbixQ828fArqxibl8FNl/UmyC5NSH8pKSABZuHq3RQeOMa0G/sTHmKnsWwvNfnLibpwLLaoli/+EqKlDMPgqiEZjMnpxorvD7D0uwN8r8uoPN4IeE4JJceFct1oxcA+yWSkRmEYBs5u93LwzQepWPgKSRMeDuiZgZkZ8bw0dTjvLNjOvBVFbPzhEPddn02vtFizo7VrUkACyP6yat5d+APZKpG8vqkAVH7zHpbgMKIH+9ZCQoiWCrJbGT0ojdGD0gDPwkPDME57pO5I6ErsiBs5suhtjm/8nKgBl/szbquFBNn4zfgsBmem8OLfN/G7l1Zy9bDu3HBJLxwyGvGJfAcSIJyuJp79v42EBFm5Z2JfDMOgdud66ndvIvqi8VhDws2OKDqYYIftjKd5IgdcRmiPARxZ/BcaDhb4Kdkvk9UzgT9PHc7oQWl8srSQ+15Yjt5XaXasdkkKSIB4a/529pYe577rsomLCqG5oZaKr2bjSDov4I/sRMdlGBYSrpyCLTKWsk9mcKKmyuxILRIabGfKtX35w2051De4mPbSSmbP2yYztVpJCkgA+GZjMQtX72Hc0Awu6J0EwJFv3qOp9ijxl90pVxwUAc0aEk7S+Gk0N9RQ9tEzNLsazY7UYtm9Enn5wRGMyenGvBVFTJmxlC0Fsm6kpaSAmGxLQTkzP9xMVo94bvbOuqrduZ7qTV8TNfAKglO6m5xQiDMLSk4n8ar7aCwp5PC8F3+y1UmgCg22c9f4LJ66KxerxeCx19Yw88NNHKtpP4XQLFJATLSn5BhPvruBzonhPHzLQOw2C66jhylf+DJBnTKIHX6D2RGFaLEwNZC40ZOp0+spX/hKu7sM7q8z4pk5dTjjh3dnycb93PG0Z4V9U1P7eh/+JAXEJLv2V/HorDWEBNmYfmsOYSF2mhtqKfvoadxuN4lX/xbDKq2pRfsSNfAKYoZMpCZ/GRWfv9auRiLgmYk2+Yo+vPTAMHp0jub1uVu59/llbC2sMDtaQJICYoJtRRWe4hFs46m7LiIhJgR3k4uyT/6Is+IASf82FXtMstkxhfBJ9EXXEp17DdWbF1M25zmaTzjP/KIA0zU5kifuyOGRyQOodzbxyKzVPPnOBvZKl99/IetA/MjtdvPlun28/ulWOsWH8j93DPbMuHI1cnjOc9Tv3UrClXcTmi4XihLtl2EYxA67HmtoBEcWvU3pe9NJGj8NW0T7WrRnGAY5mSlk90ri02WFfLqskLVbS8nNSuH60Yq0TpFmRzSdFBA/qWtw8eqcfJZ+d4BslcgDk/oTGeagqb6aso//SEPxDuLH3E7E+cPMjirEWRE18AqskXGUf/ZnDr45jYSx97TLg6Mgu5XrRisuzz2PecuL+GzlblZvKSE3K4UJI3sGZB8wf5EC0sbcbjer80t4Y+42qqobmDSmFxNG9sRiMWgsKaRsjmfufOK4ewnvk2d2XCHOqvBeOThiUymbM4ND7z9BRPbFxA6b1C4XxkaEOrjx0t6MHZLBvBVFzF9ZxOotJfRJj+Py3PO48NfJ2G0da8q9XBO9jbjdbvILK/hwUQFbiypIT43irvHno9JiaXY2ULXqI46tn48tPIbE8dNkuq44pzW7Gqla/j7H1i/AEhJOTN4EIvqNwmJzmB3NZzV1ThZ/W8yCVXsoq6wjItTO0H6dyc1KoXe3WFOuieIrX6+JLgXkLKutd7E6v4Sv1u2loPgosZFBTBjZkzGDz8M40Uj15sUcXTuXppoqIrJGEDvyZqwhEaZkFcLfGsv2cmTRWzTs2441LJrIAZcRcf4IbBExZkfzWVOzmy27ylmyoZi120pxnWgmItRO/95JDOqTTLZKJDQ4sGdUBmQBUUrdADwG2IE/aa1f/tHjfYHZQCSwArhTa31CKdUVeA9IBDQwSWtd04L9dcPPBcTtdlNaUcvmXeV8v/Mwm/RhnCeaSU0I56oh6YzITqG5VFPzw1pqdqzC3VhHcFofYofeQHCXXn7JKEQgcbvdNOzbxtE1c6jfkw+GheCuvQnNyCa0ezb2+C4B3dn359Q1uNhUUM6G7Yf4dkcZ1XVOLBaDbsmR9EyLoWeXaHqmxdA5MQKrJXDeY8AVEKVUKrAK6A80AmuA67XWO055zjbgVq31OqXUm8BGrfUspdQC4D2t9QdKqd8D4Vrr37Vgn91ogwLS3Oymus5JVXUjVccbKKus42B5DcWHqtm1v4rqOhfgplushdw0C/1TIZajOEuLaDhYgNvZgGFzENZ7MJH9RkvhEMLLVVlC9Zal1BVuxHm4GABreCyOpG4EJaXhSEjDFpOMLSIWa1gUhrX9fG3b1Oxm595KvteHKSiuYldxFbUNnl5bDruVTnGhdIoPIzkujJT4MBJiQokMcxAdHkRkuINgh//eq68FpC0TjgK+0VpXAiilPgauAZ7w3k8DQrTW67zPfwf4g1JqNjAEGHfK9uXAGQsIYAU4dOhQq8Mer3Xy0ZICjtc5cbqacTpP0Ohqpr7RxfE6F+7mfy20QVa4LnIjedZqgkIbsTU1QE0zbIeK7VBhWLDHpRKU0o/grr8iuEsvnDYHFQAHDrQ6nxDnrB5DocdQrDVV1O/bTuPBAlwHinHlfwvuUxciGhjBYViCQrDYgjDsDk9B8Y5WDKuDmLxrsEUlmvM+fkKUA4ZnRjA8MwK3uwtllXXsKTnG/rIaDlfVUbTnCOu+q8P1E6vdg+xWgoM8HZGDgqyeP+1WrFYDi+H5MQwDi8Vz22o1uOTCNFITWj9B4ZTPzFbNAmjLApIClJ5yvxQYeIbHOwPxwHGt9YkfbW+JTgCTJk3yJW+rPX7GZ+wElrR5DiGE1wufmJ3AVH/95X9FJ6CopU9uywJiAU49bDeA5hY8/uPt/Oh1P+dbIA9P0WlfPRSEEMI8VjzF49vWvKgtC8gBPB/mJyUDJT96vNNPPH4YiFJKWbXWTd7nnPq609JaN+L53kUIIUTrtHjkcVJbTlReDIxUSiUopUKB8cCXJx/UWu8DGpRSud5NNwFfaK1dwEpgonf7zcAXbZhTCCGED9qsgGitDwKPAkuBzcDftNYblFKfK6Uu8D5tEvCCUmonEA7M9G6/C7hdKbUDzyjmsbbKKYQQwjfn1EJCIYQQ/tN+1toLIYQIKFJAhBBC+EQKiBBCCJ9IARFCCOGT9tNY5hx2pqaTHY1SajowwXt3odb6QTPzBAql1AwgXms92ewsZlNKXQlMB8KAr7XW95ocyVRKqRuBh713v9BaT/XHfmUEYjJv08n/BS4C+uKZvvwrc1OZRyk1CrgY6Ifn36O/Uupqc1OZTyk1ErjF7ByBQCmVDryKp1/e+UC2UupSc1OZx7vObiYwFMgC8rz/j9qcFBDz/bPppNa6FjjZdLKjKgUe0Fo7vYtKfwC6mpzJVEqpWDwHGU+anSVAXA18qLU+4P0dmQisNzmTmax4PsvD8JzFsAP1/tixnMIy35maTnYoWuvtJ28rpXrgOZWVe/pXdAiv4VmU28XsIAGiO+BUSn2G5+BiAfB7cyOZR2td7b3sxU6gDk/38jX+2LeMQMx3pqaTHZJSqg+wCJimtd5ldh6zKKVuBfZrraWt8/+z4Rm5/yeQAwyiA5/eU0qdD/wHkIbngLQJkO9AOojTNZXssLz90ZYAD2mt3zU7j8kmAhcrpTbjuZbOWKXUCyZnMtshYLHWulxrXQ98SgcetQOXAEu01oe9DWXfAYb5Y8dyCst8i4HHlVIJQC2eppO3mxvJPEqpLsBcYKLW+huz85hNaz365G2l1GRgmNb6fvMSBYQFwLtKqWigGrgUz+9MR7UFeFYpFYbnFNaVtLItu69kBGKy0zWdNDeVqaYCwcDzSqnN3p87zQ4lAofWej3wLJ5LN+wA9gFvmxrKRFrrr4H3ge+AfDxfoj/tj31LM0UhhBA+kRGIEEIIn0gBEUII4RMpIEIIIXwiBUQIIYRPpIAIIYTwiRQQIYQQPpECIoQQwieyEl2INqaUugX4bzyttt3ARuAprfVfTA0mxC8kCwmF8AOl1F+BY0AQ0KS17rDtasS5Q0YgQvjHnXh6FtUD/U3OIsRZId+BCOEfSXh6fEXjabktRLsnp7CEaGNKKTueC/y8hueg7VYg13s1PSHaLRmBCNH2ngTKtNaztdavAxV4LlErRLsmIxAhhBA+kRGIEEIIn0gBEUII4RMpIEIIIXwiBUQIIYRPpIAIIYTwiRQQIYQQPpECIoQQwidSQIQQQvjkH1ln7a+814kAAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "sns.kdeplot(group1, label='Group 1')\n", + "sns.kdeplot(group2, label='Group 2')\n", + "\n", + "plt.xlabel('x')\n", + "plt.ylabel('PDF')\n", + "plt.title('Estimated PDFs for the two groups');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The test statistic\n", + "\n", + "The test statistic for the Mann-Whitney U test is the probability of superiority; that is, if I choose a random observation from each group, what is the probability that the observation from Group 1 exceeds the observation from Group 2.\n", + "\n", + "The following function takes two arrays representing the observed samples, and computes the probablity of superiority for Group 1 by comparing all pairs. This is not the most efficient implementation, but it is easy to compute using one of NumPy's `outer` methods (one of my favorite under-used features).\n", + "\n", + "In this example, because I sampled from a continuous distribution, the probability of a tie is small so I will ignore it." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "def test_stat(data):\n", + " \"\"\"Compute probability of superiority\n", + " \n", + " data: tuple of arrays\n", + " returns: float probability\n", + " \"\"\"\n", + " group1, group2 = data\n", + " array = np.greater.outer(group1, group2)\n", + " return np.mean(array)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's the test statistic for the observed data." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.5253535353535354" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "data = group1, group2\n", + "actual = test_stat(data)\n", + "actual" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The test statistic is far enough from 0.5 to be suspicious, but it's not clear whether I would expect that to happen by chance. That's the question we'll answer, but first we need a null hypothesis." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The null hypothesis\n", + "\n", + "There are several ways to model the null hypothesis. A simple choice is to assume that both samples were drawn from the same distribution.\n", + "\n", + "We can model that hypothesis by [permutation](https://en.wikipedia.org/wiki/Resampling_(statistics)#Permutation_tests); that is, we can combine the observations into a pool and then shuffle the pool.\n", + "\n", + "The following function implements this model. Each time it runs, it shuffles the pool, then splits it into two groups with the same sample sizes as the observed data." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "n1 = len(group1)\n", + "pool = np.hstack(data)\n", + " \n", + "def run_model(): \n", + " np.random.shuffle(pool)\n", + " return pool[:n1], pool[n1:]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's an example where we generate one batch of permuted data and compute the test statistic." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "fake = run_model()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.45494949494949494" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "test_stat(fake)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The sampling distribution\n", + "\n", + "To estimate the sampling distribution, we can run `run_model` and `test_stat` 1000 times and save the results. " + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "def sampling_dist(iters=1000):\n", + " \"\"\"Samples the distribution of the test statistic under the null hypothesis.\n", + "\n", + " iters: number of iterations\n", + "\n", + " returns: array\n", + " \"\"\"\n", + " return [test_stat(run_model()) for _ in range(iters)]" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "sample = sampling_dist()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's what the distribution of the test statistic looks like under the null hypothesis. The gray line shows the observed test statistic." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEXCAYAAACpuuMDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8HHd5+PHPanVZh2VZknVYtuXzkS/5ju0ktnOnSUggEI4SSgIFCi0/KC2ltEAp0BYKtPTHnR+UBkiAciSQkJMYO3Hi+L4t++tTtiXLkizLOqxbu78/ZrSaVSTt6tpZ7T7v18sv72hnZp+dnZ1nv8d8vx6/349SSinVK8HtAJRSSkUXTQxKKaWCaGJQSikVRBODUkqpIJoYlFJKBdHEoJRSKkjcJgYRKRGRHhE5YP87JCLbReQdjnW+KCLvDbGffxKRNw/yXGB7EfGLSO4wY1wjIt+3H68WkV8PZ/uREBGviPxORE6IyEeHiOcmETkyyteaLSK/GcX2WSLyx2FuE3gP4a4XzrF347MaIIYWESkZ5T5+ICKr7MdbReSBMQlueDE8KiKftB8P+L0ZyfdphLGMy/EQkftE5Jtjsa/xkOh2AC5rM8Ys710QkVnAZhHpMcb8xhjzT2Hs4xagfKAnwtx+KIuBYntfe4BIfEmnA3cC6caYnsHiGSOzABnF9tnAdcPcJtz3MNxj78ZnNR5uBx5xO4goMi7HwxjzFPDUWO93rHji9QY3+5fVEWNMRr+/vxv4mDFmnYg8aq/zdRH5AnA/0AnUAw8DbwX+HagD/gZ4MzAVmAv8Hsh3bO8H/h+wBquk9lljzO9F5GHgAWPMm+zXfxjrovIR4DUgC3gC+DHwbWPMEhHJAr4DLAf8wHPAPxpjukWkHfgKcAdQCHzVGPO9Ad7/BuBrQJr9nj5rv95OrIv1YeBtxpjT9vozBojnUWAHUAqkAh80xmwTkWT7uGwCvMB++5g2OV7fCxisRPSKMeZOEbne3i4d6AG+YB+jAuAnQO8vxGeMMZ8TkS3ARjvWVc5EJiI3Av9pv74f+DKwq997+HPgG8A6IBPwAB8Azg9x7MPZr3P9DOBbwA1AN/Bb4DPGmMAXT0Ru6l2//7KI/DNQYn+Ws4Aq4D3GmGr7M/yWHcdu4D3AImNMhYjca3+myUAr8EljzOv2/tYDRcBBY8x7HHH8K/B3wFngvfZncQFYABQAL9mfsW+wz4p+BjsfBzvvjTFv6ve98wN5xpjL/fbrB75rf3Y5wNeMMd8RkT8AvzTG/MBe77P28weBt2N993qP40PGmIsiUgx8zz7OHuDHxpivjcXxGOLcdb7ft9qflc/e9u+MMa/0P5aRFLdVSUM4CCx1/sG+KP41sMYYsxp4EVhrjPkOsAfrg3zSXj3NGLPYGPP3A+z7jDFmJdYX+McikjdYEMaYC8A/AduMMe/r9/Q3sZLTUmA1sAz4pP1cCnDZGHM9VoL5hoik9ns/OcCvgY8bY8qAh4DHsE7eu7FLUr1JYYh4ioFv2KWuR4B/tv/+aayL4CpjzDLgItbFwfn+erAuwqftpJAN/A/wZ/YxejPwPRGZCXzQcew2APPt5Pg+R6z9SzdfAP7TGLMKeD9wywDvYS3WBXK9MWYR1gX90yGOfTj7dfoiVtJciJXIb8BKmMOxAXi7MaYUuAZ82E6+vwL+1hizAtgCTAIQkfnAvwF32899CHhCRNLt/c0CVjiTAoAx5jNYn9WDxpid9p8zgevt+O8CbgjxWfUX8nwchTP253A/8B8ikoT1g+mDACKSgJX8e6sON2Gd84uAvVjfI4DHgS3GmKVYn897RORdY3Q8Bjt3nb4G/KV9bfkccNPYHJ6R08TwRn6sX1hOVVgJY5+IfB04YIz57SDbvzrEvr8PYIw5glX9tH6EMd6F9YvSb4zpsPd7l+P539n/78P6Yqb3234tcKr3ZDfGHMX6xXvTMOM47fjCHACm2Y/fhPXl2C8iB4C3AItC7Gs91i/K39rbPIv1WZQBzwNvE5Fngb/Aung3htjfL4HviMjjwCrgH/uvYIx5HeuX2l/Yn+sDQEb/9Ya7335uA/7bGNNjjOk0xmwyxmwNsU1/Wx2lrf1YpdKlQJcxZrP9Xn4ONNvr3I51LDfbx/JxrF+j8+zndxhjusN87f+1Y28FTmJ9xkN9VgMJdT6O1M/s/w/Y+50MPA3ki8gyrCrRs8YYY6/3ojHmhP34B8CddrK8ASuhYJ9XjxL8fXIa7vEI59z9BfCkiPwQq3r0qyM5GGNJE8MbrcGqmggwxviwfm08jPVL/RsiMtiH1zLEvp2/ahOALqwTyOP4e3IYMSbY2zmXkxzLbXbcves49w991SD995nE8HQ5Hjvfhxfrl9lyuzRxHaHr3L3Asd5t7O3WAS8YY3YDs7Gq4kqAXb0NgoMxxjyCdfH8A9YF4tAAJad7gGfsxd9hJdj+x2rY++2nG8exFpEZdonNKdQ50DbIuv1j7b3Ye4HNAxzL3s4CQ52j/Q30GQ/6WQ2yj4HOx5Gc9wPG5tyvXXJ8BKs09376SgvQd3zAOt977P/7H8ehvgvDOh7hnLt2yeRGrNqHhwFXq5FAE0MQEVmAVZT7j35/X4b1pTpmjPkyVr30GvvpbsK/oD5s728l1q+3nVjtE0tEJNUuCjsvoIPt+wXgoyLiEZEUrKqCP4QZA8DrQKmIXGfHsxirrn5riO3Cfa+98SXbxfkfYNXFD7W/HVjF7I12TMuxfpFNF5GvAJ+zS2kfB44CS+ztvSLyhou5iGzHqi55FOv4TMGqF3a+5u3A08Zqg9mDVbLxDvVew9yv00vAQyKSYH9Wv+aNVUl1wEwRmWa/l3cNsJ/+DgEeEbnbjus+rF+bAJuBO0Sk1H7ubnv9SWHsN5zPeNDPKoz99xrqvB+tH2JVL60CnnT8/VYR6Y3xw1iffTPW+/krsHq6YbUn9H6fRnU8hjh3sddNFJEKrCro7wN/CZTZ54pr4j0xTJK+7qr7sIqQ/2CMeca5kjHmIFYVwh4R2YP1S+Rv7KefAr4sIg+F8XpzRGQ/1on7LmPMFaz2ipeB41i/FPY41t9hb/NEv/18DKsIe9j+Z4B/DfM9YzfkvR34logcxiqSv89RzB7MYPH09yWgAqvaoxzrV9XfDrBeOdAuIruAy8DbgK+JyEHgp1h1thXAfwHLxeoeuwerMfAXQDVWw+/RAX6Ffwr4on28t2I1Blb0ew/fB26yj8E+4DQw205mg73XcPbr9AWsxv2D9vF41hgTtI4xphzrV+4eez9nBzhWQYwxXViJ7Et29cVbgVrH/j4E/MI+ll8C7jPGhFNSeAJ4TETuGOK16xj8swrXUOf9qBhjau39/dw+Tr0qgZ+KyDGsX+9/bf/9QaykcRjrfHoC61oAoz8eg527vdt223H8zL4G/Qp4v11F7Jq47ZWklIpNYt3fsBvYaHcMCOr15GZsE0W8lxiUUjFERD4IHMPqvnrB7XgmKi0xKKWUCqIlBqWUUkHGdUgMEZkMbAfeZKy7MW/Dumt0ElZ/4M+O5+srpZQavnGrShKRtVjdFEuxbiGvweo9swnrtvJngP8yxjwX5v5SsLqIVhN8P4BSSqnBebFuwNsdbm+n8SwxfBCrb/BP7eXrgJPGmLMAIvIYVpfJsBIDVlLYNtZBKqVUnNjA0CMzBIxbYjDGfABAJDB4ZhHWr/1e1QxvpM5qgMcff5yCgoKxCFGpmHbkSN+o6EuWLBliTRXLLl26xIMPPgjB198hRXLY7f7DOHiwxm8JVw9AQUEBxcVjOfKzUrGppqYm8Fi/M4phVMFHsldSJVY9V68CrJELlVJKRZFIlhh2AiIi87BuC3838KMIvr5SSqkwRKzEYIxpxxpE7jdYY+QcxxpQTCmlVBQZ9xKDMabE8Xgz1qQySimlopTe+ayUUipIJNsYlFK2zq4eztdYE66lJHnJy55EarJ+HVV00DNRqQhq7+jm1YNVbN1XybW2vqkCkhK9bFpZzK2rZ5Caol9L5S49A5WKkKq6Fh554hDNrZ1veK6ru4eXdp1jx+Fq3nn7ApbMzXUhQqUs2sagVARU1jbz3V8fDEoKWRkpzJk+hZysvhk3W9o6+dHTR9lnat0IUylASwxKjbuquha+95tDtLZbVUeTUpK4b+McVi/MJ9GbgM/nZ5+p5fevnqGxpQO/389jzx3D5/OzemG+y9GreKQlBqXGUXtnN//z9NGgpPCRt5WxbkkhiV7r65eQ4GH1wnz+5t2ryM9JB8Dv9/OzF45z4nyDa7Gr+KWJQalx9LuXT1Pf2AZASnIiH3lbGTPyMwdcd3J6Mh99YBlFeRmAlRx+8mw5V5tdnRdexSFNDEqNkyOnL7PjSN+Alg/cMn/QpNArIy2ZD71lKZlpyQBca+vi0WeO0t0znPEmlRodTQxKjYPW9i7+9w8nAsvLF0xjVem0sLbNykjhoXsW4fF4ADhX3cSzr50dlziVGogmBqXGwYs7z9HSZvVAmpyRwttvnR+40IdjbvEU7tswJ7C8ZV8lZy82jnmcSg1EE4NSY6y2oZVXD/SNKP/WTfNIS00a9n42rSxm/sxsa8Hv5+cvGrq6dVZbNf40MSg1xp7edoYen9UmMGf6FMrmj+xmNY/Hw7tuF1LsoTLqGlp5bnvFWIWp1KA0MSg1hk5duMqR05cDy2/eOGdYVUj9TZ2cGlSltHVfJVV1LaOKUalQNDEoNUb8fj/Pbu9rJF69sICZBZNHvd/1SwuZNyM78Bq/3nwSn88fYiulRk4Tg1Jj5OSFq4EGYm9CAndfXzIm+/V4PDxwy3y8CdbXtaK6kd3ll8Zk30oNRBODUmPA7/fz/OsVgeW1SwrInpw6ZvvPn5rGzauKA8tPbTsTNDqrUmNJE4NSY6B/aeG262aO+WvctnYWUzKtZNPa3sULO86N+WsoBZoYlBq1AUsLmWNXWuiVkuTlLZvmBpZfO1RF7ZXWMX8dpTQxKDVKZy82jXtpoVfZvFzmFk8BwOfz89S2M+P2Wip+aWJQapQ27zkfeLx6Uf64lBZ6eTwe3ryxr9Rw9MxlHYFVjTlNDEqNQvXla5SfqbcWPB5uWTVj3F9zRn4mqxcWBJaf3nYGv1+7r6qxo4lBqVHYsvdC4PGSOTlMm5oWkde954bZgfkcKmubOey4qU6p0dLEoNQIXW3uYO/xmsDyrWvGr22hvymZKdy4fHpg+bntFXrTmxozmhiUGqFtB6oCF+M506dQUjj6u5yH49bVM0hO8gJwqf4aB07WRfT1VezSxKDUCHR29QRNwnOT4+azSMlIS2bjir7Xff71Cnq01KDGgCYGpUZg7/GawDzOUydPYvHsHFfiuHlVMakpfaOvHtRSgxoDmhiUGia/38/L+6sCyxuWF5GQMPIRVEcjLTUpqNSwZc8F7aGkRk0Tg1LDdPLCVWrqrwGQnOTlusUFIbYYXxuWFQX1UDp54aqr8aiJTxODUsP0iqO0cN2ighHNzjaWMtKSWbukMLC8eff5IdZWKjRNDEoNQ31jG0fP1geWNzi6jLrp5lXFgQmBTpxvoLK22eWI1ESmiUGpYdh2oArsOvzSkqkRu6EtlJysSSxfkBdYfnlf1RBrKzU0TQxKham9s5sdR/omyNm4PPJdVIeyaWVfPAdO1NLe5XMxGjWRaWJQKky7y2vo6OwGIC87DZmV7XJEwWbmZ1I8LROA7h4fJ6raXY5ITVSuJAYReY+IHLX/fd2NGJQaDp/Pb1Uj2TYsn+5aF9XBeDweri/ra4QuP9+qXVfViEQ8MYhIGvBNYBOwDNggIrdFOg6lhuNU5VXqGqxJcVKSE1mzKN/liAa2QqaRkmzd8NbY2sPFKzr9pxo+N0oMXvt104Ek+1+bC3EoFTbn8BdrFuaTal98o01qciKrF/YlrWMX9Kulhi/iicEY0wx8DjgOVAIVwPZIx6FUuFraujh0qm9Y6/VLC4dY233O6qSzNR20d2ojtBoeN6qSyoD3A7OAIqAH+GSk41AqXHvKL9HTY11cZxZMpigvw+WIhlaUm8HMAmukV5/fz5lL2githseNqqQ7gc3GmFpjTAfwKHCTC3EoFZLf7+d1RzXSuiXRXVrotbq0rzrpZLUmBjU8biSGg8BtIpIuIh7gXmC3C3EoFdLZi03UXulrdF4heSG2iA7LJY/ePlOXGrqob9S2BhU+N9oYXgR+DuwFDmE1Pn8l0nEoFY7XD/eVFlZIXtQ2OveXmZZMcW5KYHnf8VoXo1ETjStnuTHm34F/d+O1lQpXa3tX0PwG6ydINVKv+UWpXLjcAcCe4zXcdt3MwHhKSg1F73xWahB7j9fS1d0DQFFeBjPyM12OaHhKpqWQ5LUSQe2VViprW1yOSE0UmhiUGoDf7w+6d2HdksIJ92s7KdFDSX5fdZLO7qbCpYlBqQFcqGnmYp31Czsp0cuq0mkuRzQycwpSA48PnbqsQ2SosGhiUGoAzi6qy+bnuT4Zz0gV5yQHqpPqGlqpsXtYKTUUTQxK9dPe2c1+42h0jvI7nYeS6PUwM89ZnXR5iLWVsmhiUKqf/aYuMLx2/tR0ZhdNdjmi0XG2Mxw6pe0MKjRNDEr1E9zoXDDhGp37m5mXjNdrfdUv1rXozW4qJE0MSjlcrGvh/KUmALzeBFYvKnA5otFLTkwImlTIOSCgUgPRxKCUg/NO57J5uWRMmpiNzv2Vzc0NPD6siUGFoIlBKVtXdw97HUNHTJQB88KxeG4u2FViFdVNXGvTCXzU4DQxKGU7ePIybR3WBTMnaxLziqe4HNHYyZiURIk9FLff7+dYxRWXI1LRTBODUrb+jc7RNqfzaC2aPTXwuPxsvYuRqGiniUEprLGETldeBcDj8bAmBhqd+1s0Jyfw+HhFAz0+vQtaDUwTg1IElxYWz8khKyNliLUnpqLc9MD7auvoouJio8sRqWiliUHFve4eH7vLawLLE/lO56F4PB4Wze4rNRzV6iQ1CE0MKu4dOX2ZlrZOAKZkplI6a2qILSau4HYGbYBWA9PEoOKe896FtYtjr9HZacHMbBLtu6Br6q/pXdBqQJoYVFyrb2zjxPkGa8HjYe3i2Gt0dkpO8jJvRl83XHOuwcVoVLTSxKDi2o4jlwKPS2dlkz05dYi1Y4Ozquy4JgY1AE0MKm71+PzsKu9LDLHa6NxfaUlfYjhxvoGeHp+L0ahopIlBxa1jZ+tpaukAICMtmcWOHjuxbFr2pEDJqKOzm4rqJpcjUtFGE4OKW85Z2tYuLggMTR3rPB6PViepIcXHN0Gpfq42d3DM0V0zlgbMC4dzGO7j57TbqgqmiUHFpZ1HL+H3W0NCzJ+RTe6USS5HFFkLZmYHJiCqrGmmubXT5YhUNNHEoOKOz+dn51HHgHlx0ujsNCklkZLCvilLtduqctLEoOLOifMNNDS1A5CWmsRSxyQ28WRhibOdQauTVB9NDCruOO90XrMon6TE+PwaONsZzLkGfDraqrLF5zdCxa3m1k6OnOmb2jLeGp2diqdlkm5PXdrS2snFyy0uR6SihSYGFVd2l18K/DKeXZRFQU66yxG5JyHBgzi7rVZoO4OyaGJQccPv9wcNgRGPjc79lWq3VTUATQwqbpyubKSuoRWA1JREls/Pczki9zlLDGcvNtLe2e1iNCpaaGJQccN5p/Oq0nySk7wuRhMdJqcnU5SXAVjdeE9euOpyRCoaaGJQcaG1vYtDJ+sCy/EyYF44gobHqNDqJKWJQcWJ3eU1dNujiM7Iz2S6/StZ9bufoaIhcEe4il+JbryoiNwLfB5IB140xnzcjThUfPD7/UHVSPHcRXUgJUWTSUlOpKOzmytNbdRdbWNadprbYSkXRbzEICJzgO8DbwHKgJUiclek41Dxo6K6iZr6a4A1g9nK0mkuRxRdEr0JzCvWWd1UHzeqku4H/tcYU2mM6QLeCex0IQ4VJ5xdVFfINFKTXSkoR7XSEu22qvq48Q2ZB3SKyFPATOD3wOdciEPFgfaObvab2sDyeq1GGpCzAfrk+at0dfvidqgQ5U6JIRG4DfhzYD2wFnjIhThUHNhraunq7gGgMDeDmQWZLkcUnXKnTCInyxp6vKu7h7MXG12OSLnJjcRwCXjJGFNnjGkDngSucyEOFQd2OAbMW7+kMDAHgXojHW1V9XIjMfweuFNEpoiIF7gL2OtCHCrGVdY2U1nbDFgNrKsWaqPzUIJGW9Vxk+JaxBODMWYn8FXgVaAcOAf8T6TjULFv19G+Rudl8/NIS01yMZroN2/GFLwJ1iXh4uUWGls6XI5IucWV7hnGmB8BP3LjtVV86O7xsfd4X6PzdYsLXIxmYkhNTmT29CxOXbBKC+Zcgx63OKXdDlRMOnb2Cq3tXQBMyUwN6qevBqejrSrQxKBi1K7yvmqk1QvzSUjQRudwlDoaoHVWt/iliUHFnJbWTsrP1geW1yzKdzGaiaUoN53MtGTAGnjwgt14r+KLJgYVc/aZ2sAv3ZLCLB33Zxg8Hp3VTWliUDFo19GawGMtLQxf8Gir2s4QjzQxqJhy8XILVXV99y4sX6CztA3XglnZYN8IeO5SU6ARX8WPIRODiKyMVCBKjYXd5X2lhaXzcvXehRHImJTEjGnWfBV+v5+T53VWt3gTqsTww94HIvLZcY5FqVHp8fnZc8xRjbRQ++CPVKkOjxHXQiUGZx+/t45nIEqN1vGKK7S0dgIwOT0laIgHNTzOdobyiis6q1ucCZUYnGeDdgRXUW23496FVQun6b0LozCzYHKgGq6ppYPK2haXI1KRNJzGZ/3JoKJWW0c3R8703btw3SKtRhoNb4KHhbP7Sg1HHcdWxb5QYyUVi8g3B3gMgDHmY+MTllLDc/RMPT09PgCKp2VSkJPuckQT3+LZOey122yOnqnnT9aXuBuQiphQieE7gzxWKqocPFkXeLxsvnZRHQulJVNJSPDg8/mprG3manMHUzJT3A5LRcCQicEY84VIBaLUSLV3dgfdiFU2P9fFaGLHpJRE5hZP4eR56+7n8rP1XF9W5HJUKhJCDrstImuATwBLgVbgMPBfxpgj4xybUmEpP3uFbrsaqSg3Q4fAGEOLZ+cEEsPRM5oY4kWoG9xuBZ7CSgZ/D/wzcB54UUQ2jXt0SoXhkKMaSUsLY2vJ3JzA4xPnG+jo6nExGhUpoUoM/wDcaYw55PjbcyLyPPBl4NZxi0ypMHR29XDMUY2k7QtjKydrEvk56dTUX6O7x4c5d4WyeXqMY12o7qrT+iUFAIwxuwCd+US57vi5K3Tav2KnTU0jf6pWI421pXP7SmGHTl12MRIVKaESw1DlRr17SLnu0Mm+C1XZvDw8Hj0tx1rZvL7EUH6mrz1Hxa7h3PmsVFTp6vYF3dS2TNsXxkXxtAyyJ6cC0NbRxalKHVQv1oVqYygVkTdUJWGVFuaMQzxKhe3E+QY6OrsBqy58el6GyxHFJo/HQ9m8XF7eVwlYpbRSx2Q+KvaESgx3YXVTrQZSgcpxj0ipMDnru8vm52o10jhaOrcvMRw+fZkHbpmvY1HFsFCJYQ7wReAkMBd4tzHmxXGPSqkQenp8HDndlxiWaU+ZcTW7KIuMtGRaWjtpae2korqJOdOz3A5LjZNQbQwfA5YYY9YC9wKfHv+QlArtVOXVwMxiUzJTmVmQ6XJEsS0hwRPUO8k5BImKPSFHVzXGXLT/fx3Qn2UqKhwM6o2k1UiR4Lx58MDJOnw+7ZsSq4bbK6l7vAJRKlw+n5/Dp4MTgxp/82dkkz6pb46GMxcbXY5IjZfhzMcA2n1VRYEzFxsDM7VlpiUzu0jruiPBm+AJuuv5wAmtTopVoRqfy0SkybGcZi97AL8xZvL4habUwJxjIy2dl6u9YyJoheTx+uGLgNXOcP9N8/Dq8Y85oRLD3IhEoVSYfD5/UDdVHRspsuZOn0JmWjLNdu+k05VXWTBT59aONaHmYzgXqUCUCse5S000tnQAkJaaxNxiHbIrkhISPCybn8erB6sA2H+iVhNDDBpuG4NSrnKWFpbOy9VqDBesLJ0WeHzo5GUdOykGaWJQE4bf7w8aNE9vanPHrILJTMm0xk5qbe8KGvZcxQZNDGrCqKxt4UpTGwCpKYnMn6nVSG5ISPCw2lFq2HusxsVo1HjQxKAmDGdvpMVzckj06unrllUL8wOPj56pD9yFrmKDa98sEfm6iDzq1uuricXv93PQ0b6wXHsjuaogJ53iadYwJN09vqA70dXE50pisOeSfsiN11YTU3X9NeoaWgFISU5EdNhn1zlLDXuPa3VSLIl4YhCRqcC/Av8W6ddWE5ez0XnR7KkkJWo1kttWlU4LjFF1uvIq9Y1tLkekxoob365HgM8ADS68tpqgDjjaF7Q3UnTITEtGZvXdw7C7XEsNsSKiiUFEPgBcMMZsjuTrqomt9korNfXXAEhK9FI6W6uRosV1iwoCj3cevaQjrsaISJcY3gncISIHsCYAuk9EvhHhGNQE4xz7v7RkKilJXhejUU5L5uYGRly92tzOifNaERALQo2VNKaMMbf3PhaRh4GbjDGfiGQMauJx9nhZNl+H2I4mSYkJrF5YwMv7LgDw+uFqSku0RDfRaQueimr1jW1U1TUD4PUmsHh2jssRqf7WLemrTjpy5jLN9pDoauJyLTEYYx41xjzs1uuricFZWpBZ2aSmRLSQq8JQkJMemBPD5/NrI3QM0BKDimoHtTfShLBuSWHg8Y4j1fj92gg9kWliUFHranMH5y9Z80QlJHhYMlerkaLVsgV5pCRbpbm6hlbOVOm0nxOZJgYVtQ6d6istzJ+RTVpqkovRqKGkJHlZ5RhYb8eRSy5Go0ZLE4OKWs72hbJ52hsp2jmrkw6erNOB9SYwTQwqKjW2dHDmolUd4fF4WKqJIerNyM9kep41sF5Xdw/7TK3LEamR0sSgotLBk3VgN2DOmZ5FZlqyyxGpcKxb2td1dcdhrU6aqDQxqKi03/S1L6yUaUOsqaLJSplGUqJ1Z3pVXTPnqptcjkiNhCYGFXWuNLVTUd1XjVSmcy9MGGmpSaxY0Pd5bTtQ5WI0aqQ0Maioc+BEX920zMomY5L2RppIblw+PfD4wIk6mq7pndATjSYGFXX2n+irRlqh1UgTzowAcOuwAAAXTUlEQVT8zMCd0D0+H9sPXXQ5IjVcmhhUVKltaKWypm9spKVztTfSROQsNWw/dJHuHp+L0ajh0sSgosp+RxfHhSVTmaRjI01Iy+blMjkjBYDm1s6goU1U9NPEoKKG3x88AJtWI01cXm8CN5QVBZa1EXpi0cSgosbZi02BeYNTUxK1GmmCW7+0EK/XusScq24KjHulop8mBhU1dpf33RC1YsE0khL19JzIMtOSWe7ouvqKlhomDP3mqajQ1d0T1BtpzaJ8F6NRY2Vjv66rOonPxKCJQUWFw6fq6ejsBiAvO42SwskuR6TGwsyCycyyP8ueHh/bD1W7HJEKhyYGFRV2OaqRVi/Mx+PxuBiNGksblxcHHm8/dJGubu26Gu00MSjX1Te2Yc43WAseD6sXajVSLCmbn8vkdKvratO1DvYd16k/o50mBuW67YeqAyOpls7KZurkVJcjUmMp0ZvAxhV9bQ1b9lXi8+nUn9FME4NyVVe3j11H+6qRnH3fVexYv7QwMPVnTf01jlVccTkiNRRNDMpVh09dpqXN6qkyJTOVRbN1XudYlJaaFDTD25a9F1yMRoWiiUG56jXHAGvrlhSSkKCNzrFq08rpgc/3dOVVveEtimliUK6pvnyNM1VXAWvehXVLCkJsoSay7MzUoGFO/rhHSw3RShODcs3WfZWBx0vn5pJlD7qmYtfNq2YEHh86dTkwBIqKLpoYlCuarnWy19Ft8aZVxUOsrWLF9LwMZNZUwBo0ceveyhBbKDdoYlCu2Hagih57jP6ZBZP1Tuc4csvqvlLDzqOXaGnrcjEaNRBNDCriOrp6gmb1unlVsd7pHEfmz5hCUV4GYI2R9dpBHVwv2mhiUBG3u/wSre3Wr8SpkyexdF5eiC1ULPF4PEGlhlf2V9Fuj5OlooMmBhVRPT0+tuzpq1feuGI6Xu2iGneWL5hGTtYkAFrbu3Re6CijiUFF1J5jNVxpsnqipKUmsVa7qMYlb4KHW9f0lRq27K2ks6vHxYiUkyYGFTE9Pj9/2HU+sHzzqmJSk3VO53i1ZlEBUzKtcbFaWjt5/bAOyR0tNDGoiNlvagP91tNSk7hx2fQQW6hYluhNCGpr+OPeCzokd5TQxKAiwufz8+LOc4HljSumk5qipYV4t25JAZlpyQA0tXQEzcuh3KOJQUXEnuM11DW0ApCSnMiG5VpaUJCU6A0qNWzefT5wf4tyjyuJQUQ+LyJH7X9fdSMGFTld3T6ef70isHzzqmLSUpNci0dFl/VlRaRPss6HhqZ29hzTiXzcFvHEICK3AXcAK4DlwCoRuT/ScajI2XGkmoamdgDSJyWxaaUOf6H6pCR5uclxTry0+zw9OpGPq9woMVQDf2uM6TTGdAHHgJkuxKEioKOrJ6ht4bY1M7UnknqDG5dND5QiL19tY7+pdTmi+BbxxGCMOWqM2QEgIvOBdwDPRjoOFRmv7KukpdWaiCcrI4UbtCeSGkBqSnC70ws7KrStwUWuNT6LyGLgD8DfGWNOuhWHGj/NrZ1sdoy5f+e6WSQlan8HNbCNK4JLDbvKta3BLW41Pt8AbAY+bYz5sRsxqPH3/OsVdNhj4OTnpHPd4sKhN1BxLS01iZsdw6+/uPOc3tfgEjcan2cAvwXebYz5RaRfX0XGpfprQXey3rdhjo6JpELasKKYDPu+hqvN7ezQu6Fd4UaJ4ZNAKvCfInLA/vdhF+JQ4+jpV8/g91s9S+bPyGZhyVSXI1ITQUqSl1vX9PVFeXHXOR151QUR7x5ijPk48PFIv66KnKNn6ik/Ux9YvnfDHJ1vQYXthrIitu69QGNLBy2tnWzdW8mfrC9xO6y4oi2Bakx1dft4cuupwPLaJYXMyM90MSI10SQlJnD3DbMDy1vsJKEiRxODGlNb914IDJQ3KSWJNzm+4EqFa3VpPkW51ixvnV09vLDjXIgt1FjSxKDGTENTe9Cw2nddXxJoSFRqOBISPNy7YU5geceRai7VX3MxoviiiUGNmd+9cpqubmuylaLcDK4vK3I5IjWRyaxs5s/MBsDv9/PE1lOBDg1qfGliUGPixPkGDp6sCyy/7Zb52j1VjYrH4+Etm+YGOi6cPN/AoVOXXY4qPmhiUKPW0+PjiS19Dc6rFuYzZ3qWixGpWFGUm8ENjpLn7145rVOARoAmBjVqL++vouaKVf+bkpzIvTfOCbGFUuG76/qSoGG5N+8+H2ILNVqaGNSo1Da08tz2s4HlO9fNIisjxcWIVKxJS03iHkfvts27L1B9WRuix5MmBjViPp+fX750gm57FMyivAw26sxsahysXVxISaFVPdnj8/GLPxh8OmfDuNHEoEbs9cPVnK68ClgNhX96h+D16imlxl5Cgod33r4gcH6dv9TEKweqXI4qdum3WI1IbUMrT207HVi+dc1MiqfpHc5q/BTkpHP7dbMCy8++dlbvbRgnmhjUsHV1+/jJM8cCvUPyp6Zzx9pZIbZSavRuXTODQvuO6K7uHn767DEdmnscaGJQw/bMa2eoqmsGwOtN4M/uXqgT8KiISPQm8J67Skm0q5QuXm7hqVdOh9hKDZd+m9WwHDpVx8v7KgPL922Yw/S8DBcjUvGmKDeDt2yaF1h+9WCVzhE9xjQxqLBV1bXw2HPHA8uL5uQEzdOrVKRcX1bI0nl5geWfvXCciuomFyOKLZoYVFiarnXyw98dCYyFlJM1iXffUarzLChXeDwe3nX7AvKy0wDo7vHx308dCYzsq0ZHE4MKqbW9i0eePMTV5nbAurv5A29eErgbVSk3pKUm8cE3LyUt1ToPW1o7eeTJw1xt1rkbRksTgxpSa3sX3/31IS7WtQDWL7WH71lEQU66y5EpBXnZk3j/vYsD9zfUNbTy7V8d4EpTu8uRTWyaGNSgGls6+N5vDgV6IOHx8K7bhVKdv1lFkbnFU3jo7kV4E6zLWX1jG9/65QGq7B8zavg0MagBVdY2842f76Oytjnwt3fdtoDrFhe4GJVSA1s6L5f3OUoOV5vb+b+/2M+eYzUuRzYxaWJQQfx+P7vKL/HN/z0QmGfX4/HwztuFtUsKXY5OqcEtnpPDB+5bQkpyImDdAPf488d47PljtLR2uhzdxKKJQQW0tHXx6O/L+fkLxwO9j1JTEvnQ/UtZp0lBTQClJVP5xJ+uZNrUtMDf9h6r4cs/3s3rhy/S06N3SYcj0e0AlPt8Pj+7yy/x9KtnuNbWFfh7XnYa7793sTY0qwklf2oan/jTlfxq80n2Hbeqklrbu/jlSyfYvPsCd66bxUqZpgM+DkETQ5w7V93Eb18+TUV1Y9Df1y8t4s2b5pKS5HUpMqVGLjU5kT+7ayGrS/P51R9P0GD3UqpvbONnLxznmdfOsnHFdNYtKQx0d1V9NDHEqfrGNn7/6lkOnAgeSmBKZioP3DKfxXNyXIpMqbGzcPZU/v69a9i2v5IteytpbbdKxI0tHTy97Qwv7jzP2sUFbFwxnZysSS5HGz00McSZ1vYuXtx5jlcPXKTH11ff6k1I4KZVxdy+dpaWElRMSUnyctt1s7hh2XRe2V/FqwerAo3RHZ3dvLK/km0Hqiibl8tNq2ZQUjjZ5Yjdp4khTnR1+3j1YBV/2Hmeto6uoOeWzc/jnhvmkJetv5hU7JqUksid62Zxy+oZ7D1ew9Z9ldTY8zn4/X4Onqzj4Mk6SgqzuGlVMUvn5pKQEJ9DvmhiiHE+n599ppbntldwpSl4HJmSwizu2ziH2UVZLkWnVOQlJSawbkkhaxcXcLyigS37LnDyfEPg+YrqRh79fSM5WZPYtKKYNYvzSU2Or0tlfL3bOOL3+zlyup5ntp8N/CrqlTtlEm+6cQ5l83J1EDwVtzweDwtnT2Xh7KlU1bXw8r5K9h2vDVSx1je28cTWkzy/o4KbVhazYcX0uEkQ8fEu44jf7+fkhas889pZzl8KHoY4LTWJP1lXwvqywsBEJ0opmJ6XwbvvLOWeG2bz6oEqth+uDjRUt7Z38ez2s2zdV8nNq2Zw4/KimE8Qsf3u4kiPz8/hU3X8cc8FLtQ0Bz2XkpzIppXF3LyymNQU/ciVGkxWRgr33DiH29bOYtfRS7y8rzIwlHdrexfPvHaGrfsucMvqGdywbHrMdtTQq8QE5vf7qbnSyp5jNewur6HpWvBww4neBG5YVsRta2aSkZbsUpRKTTwpSV42LJ/O9UsL2XOshhd3ng+00V1r6+LpbWfYsreSW1bPYP3SwpgrQcTWu4kDza2dnL3YyJmqRo6eqefy1TdOTJLoTWDNogJuXzuT7MxUF6JUKjZ4vQmsXVLI6oX57Cqv4cWd5wLzkrS0dvLUK6d5Ycc51i0p5PqyQqZlp4XY48SgicFFPp+fjq4eOjp7aO/sDjx2/q3pWieNLZ3UN7ZRc6V1yMHAMtOSWbe0kA3Lp5OpJQSlxozXm8D6pYWsWZTPrqOXeHHnucAgkx2d3by87wIv77tASWEWK0unsWRODtmTJ+6PMlcSg4i8G/gskAT8lzHmO27EMRa6un20tndxrb2L1rZu6//2Lto6urnW1k1rRxdt7d202v+cCaB3oLrRSE7yUjprKmsW5bNwdg7eOO13rVQkJHoTuL6siDWLCthVbrVB1DW0Bp6vqG6korqRJ7acpCg3gznTsygpmkzxtExys1InzPhMEU8MIjId+FdgFdABbBeRLcaY8vF6za5uH5W1zfj94PP78fv9fY99fX/z+f10d/vp7O6hs6uHzi4fXd09dNj/Wxf7Llrbu2lpsxJAZ9foL+7DkZTopXhaBrOLsphXPIV5M6aQlDgxTjalYkVSYgI3lBWxfkkh5lwDrx26yLGKenw+f2Cdi5dbuHi5hVcPVgHW6AJTs1KZnJ5MxqQkMtOSmZyeTGpKIsmJCSQleklKTCApMYGEBA8JHg94wJvgISXJS2FuesS6l7tRYrgN+KMx5gqAiPwaeAD4YojtvACXLl0a1ou1d/bwvV8fpKUtCsdj93hITvKSnJhASpKX5GTrcXKSl+QkLylJXtInJZKZlkJWegq5U1LIykhxnByt1FxqHfIlVPyqq6sLPK6srHQxktiWmQR/smoKGxenc/TMFU5caOBcddOAQ3w3Xx356yyak8vbbp437O0c18ywu1C5kRiKgGrHcjVwXRjbFQI8+OCD4xGTUkpFta3Ad/9lVLsoBE6Hs6IbiSEB8DuWPUA4s2fsBjZgJZLI1t8opdTE5cVKCrvD3cCNxFCJdYHvVQBcDLWRMaYDeHW8glJKqRgWVkmhlxuJ4SXgn0UkD7gGvA34kAtxKKWUGkDEu7MYY6qAzwBbgAPAz4wxuyIdh1JKqYF5/H5/6LWUUkrFDe0Ar5RSKogmBqWUUkE0MSillAqiiUEppVQQHV11mEINACgi9wNfwLqpZDfwIWNMp4g8BHwFqLFXfcYY85nIRR7+4IUicg/wbWPMbHt5CvA4MAeoA95hjBne2CSjMIq4NwFPABfsVfYbY94XgZCdMYU6Xz4PvB/onXT4B8aY74jIcuCHwGTgFeDDxpjuCRD3gH+PUNjhxC3AI0A2cAl4lzGmQURmAo8B0wADPGiMaYlU3KOMfcyvLVpiGAbHAIA3AsuBD4nIIsfz6cC3gduNMYuBVOBh++nVwN8YY5bb/yKdFIaM3bFePvB1rDvSe/0LsM0YsxD4AfB/xz/iQDyjiXs18HXHMY90Uggn9tVYX/DeGHsvBo8BHzXGLMB6Tx+cIHEP9nfX4xYRD/AU8BVjzDJgP/Bp++nvAt81xpQCe4DPRSruMYh9zK8tmhiGJzAAoDHmGtA7ACAA9t9KjDE1IpKG9euj95fTGuAhETksIo+JSHY0xe7wQ6wSj9M9WCUGgJ8Dd4lI0rhFGmw0ca8B7hCRQyLylIjMGOdY+wsn9tXAP9oxfltEUkVkFjDJGLPDXudR4O0Ri3qEcYf4ezTEvRK4Zox53l7+N+A79rm80V4fIn+8YYSx24/H/NqiiWF4BhoAsNi5gjGmS0Tuwqq+yAVedKz7JaDMfu7b4x5tsJCxi8jHgH3ADoIFtrWrM5qAvHGLdJDXtg0n7qvAt4wxZcCzwC/GMc6BDBm7iGRg/fL7O6wv/hSsX6oh3/M4G1HcQ7yfSAl13OYBl0Tkv0VkH/A9oAXre9rkqKqL9PGGkcfeu+6YXlu0jWF4whoA0BjzHJAjIv+G9QG+2xhzf+/zIvJVhjl2yRgYMnYRWYI1PMmtvPFL0X8Q+HAHPhwLI47bGPNhx+Pvi8hXRCTLGNM4viEHDBm7XYd9d++yiPwH8COsJDaSgSbHyojitqswBno/kao2DfX9TARuAjYaY/aIyJeA/7Tj63+nbySPN4w89ofH49qiJYbhqcQe/tsWNACgiEwVkTsczz8OlIlIloh8wvF3DxCxhkTbkLFjFZ0LsepXnwWKRGSb/VyVvT4ikghkAvXjHbBtRHGLSIKIfEZE+o9BH8njHup8mSki73c87wG6Qm0XASOKe4j3Eymhjtsl4KQxZo+9/HOsIf9rgSzHuVJIZI83jDD28bq2aGIYnpeAW0Ukz25DeBvwvON5D/CY3cMBrIvWq1hFvk+JyFr77x8FnoxQzL2GjN0Y83ljzAJjzHKsX30XjTG9o+A+C7zXfvxOrIboSH3hRxS3McYH3G+vj4i8F9hp199GSqjzpQ34qojMthsX/wp40hhzDmgXkRvs9f4MeC7a4x7i79ES93YgT0SW2cv3Anvtc3kb1rkN1rkeyeMNI4ydcbq2aGIYhsEGABSRZ0VktTGmHmuk2N+LyEFAgL83xvQA7wC+JyLHsKY1/VQ0xR5i888B60TkKPCXWF/4iBhl3A8Bf23H/T7gA+MbbbAwzpc64C+Ap7G6SHqA/7A3fxD4hogcBzKAb0Z73CHeTzTE3Yb1Y+EH9jlxC/C39uZ/idUTqBxrWoDPRiru0cQ+XtcWHURPKaVUEC0xKKWUCqKJQSmlVBBNDEoppYJoYlBKKRVEE4NSSqkgeuezmvBE5JtYY90ALALOYvWpB1hvd/Ubzv48wB+AB4wxV8PcZi7wZWPMO8Jdzx676WeO+0VGvb5SY0ETg5rwjDEf630sIhVYQybvGXSD0LxYQ2wMx2xgwXDWM8ZcwOozP5brKzVqeh+Diil2YnjAmRhEZDHWUOHZWBf9bxhjfiwimVgjac7FGpdmF/AR4CfAe4AjwJ3GmIv99vUDIAXrBq5HgP8ByoHpwBZjzN0i8jmsu1NTgXTgE1h3sgbWAz4G7DHGTAlnv/3WTwK+hnW3dzfWnbsfjeAd6SqGaRuDimn2BfRXWHeJrsIaiOwf7LumHwBS7OE0rsOaIKUE6y5pgA3OpGD7FPCEva832fvrBj4MGDspzAE2YQ14VgZ8HviCMabTud5w99tv/f+DNZpmGbAEyCHyQ0WrGKVVSSrWLcSaee7HItL7txRgBfBH4Esi8kessWq+bow5aw8UOJgngR+JyHp7m48ZY3yOfWOMOWMPJvceEZkHXI81rMVQQu63n9uAnxhj2u3lgeaoUGpEtMSgYp0XqHfMbrUcWA/81BhzGmuc+69izR3wRxHp/8s8iDHmt1h1/r/GmpTmiIg4R8XELo28hjUK7Qv2/vsPXT7s/fbTjWOYZhHJD7G+UmHTxKBiXTngE5F3AYg1O9pRYJmI/B/g/wEvGGM+BWzGKkn0YF103zBLnYj8EnirMebnWO0R17BKJN2O9W8CdhhjvoFV9/8WrARFv/WGu1+nl4AHRSRZRBLs96FVSWpMaGJQMc0Y0wHcB3xERA5hNQB/2hizE6vheRJwVET22o+/Y4zxA78BXhWRhf12+QXgYXv03B3AL40xr2E1VPtE5HWseTgK7ZE6jwCNWEMmp/dbb7j7dfoucAhr5rrDwHn6pnpUalS0V5JSSqkgWmJQSikVRBODUkqpIJoYlFJKBdHEoJRSKogmBqWUUkE0MSillAqiiUEppVQQTQxKKaWC/H8Jd/DrZQ/GfwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.axvline(actual, color='0.8', linewidth=3)\n", + "sns.kdeplot(sample, color='C0', linewidth=3, alpha=0.8)\n", + "\n", + "plt.xlabel('Test statistic')\n", + "plt.ylabel('PDF')\n", + "plt.title('Distribution of the test statistic under the null hypothesis');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The p-value\n", + "\n", + "Finally, we can compute the p-value, which is the probability that the test statistic exceeds the observed value." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.25" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.mean(sample > actual)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example, the p-value is fairly large, which means that the observed value (or higher) could plausibly occur even if the samples were drawn from the same distribution.\n", + "\n", + "In other words, with these sample sizes, this test is not able to rule out the possibility that the observed difference between the groups is due to chance." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}