{ "metadata": { "name": "", "signature": "sha256:3b3e54d686bec6c68a86e5ebc6212835523c11ef8591db89181071b79423a143" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import pylab as pl\n", "\n", "import warnings\n", "warnings.filterwarnings(action = \"ignore\", category=DeprecationWarning)\n", "\n", "pl.xkcd()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 28, "text": [ "" ] } ], "prompt_number": 28 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## t-test for comparison of sample means\n", "- [source](http://dataiap.github.io/dataiap/day3/hypothesis_testing.html)\n", "- steps: \n", " 1. get a hunch from plotting\n", " 2. use statistics to quantify the evidence backing the hunch" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "from scipy import stats\n", "\n", "## data comparing averages of heights from two twons\n", "town1_hts = [5, 6, 7, 6, 7.1, 6, 4]\n", "town2_hts = [5.5, 6.5, 7, 6, 7.1, 6]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "## plotting - use boxplot to show the difference between two samples\n", "fig, ax = pl.subplots(1, 1)\n", "_ = ax.boxplot([town1_hts, town2_hts])\n", "_ = ax.set_xticklabels([\"town1\", \"town2\"])\n", "_ = ax.set_ylim((3.5, 7.5))" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEECAYAAAA1X7/VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8FNXdP/DPJrvZ+32TTQBFQUAFKVWsoGIjF6VFUcQ+\nVlsrLTw8opSH1lpRfpVgWy/VVqH2QX3ZYtVar3hDRdsCglgFERAFGhElQO57v1+y5/fH5AxZs9mZ\nwO4kWb7v12tfSzaTndlh9jNnzpyLijHGQAghZMAq6+sNIIQQcnwoyAkhZICjICeEkAGOgpwQQgY4\nCnJCCBngKMgJIWSAU0stEIlEMGPGDFx22WWoqqpCeXk5AGD06NEYN25c1rI+nw9r167FhRdeCI1G\ngx07duC0007DGWecUZytJ4QQApVUO/Ldu3dj7Nix3V5/6623MH369KzX9u7dizPPPFP8+fTTT8eW\nLVvgcDgKtLmEEEK+TrJEfuDAAVx88cWYNm0a1Go13njjDZxzzjndQpwbMmQIxo8fjwkTJmDBggWw\nWCwF32hCCCFHyapaqaurw0UXXYR9+/bhjTfewD333NPj8rW1tVi9ejUYY9BoNAXdWEIIId1JVq1w\njDFceOGFuOuuuzBlypScy+zduxff+c53cNVVV+Hf//43hgwZgscffxxWq7WgG91XMpkMotEoPB4P\n2traEIlEEIvFEA6H4fF4EAgEEI/HkUwmkUgkEI/HkUqlEI1GEQqFEIvFkE6nkclkkMlkxPdVqVRQ\nq9UoLy9HeXk5tFotzGYzLBYL9Ho9TCYTrFYrTCYTzGYzzGYzjEYjHA4Hqqur4XA4UFZWOvetI5EI\n/H4/AoEAgsEgvF4vAoEAotEootEo4vE4wuEwAoEAIpEIwuEwotEoEomEuH8ZY/j6oc33s0ajgUaj\ngVqthl6vh8FggNFohMlkgsVigdVqFZ+1Wi1sNhuqqqpgsVhgMBigVkuWf/q1ZDIJj8eDYDCISCSC\nQCCAtrY2+Hw+RKNRBINBhEIhJBIJJJNJxONxxGIxJBIJpFIppFIpdHR0ZB3DgLB/y8rKoFarUVFR\nAa1WC61WC41GA51OB6PRCLPZnHUc833rdDphsVhgsVhKpgDIGEMwGEQgEEAoFEIgEIDH44HH40E4\nHEY8HhePZ54ZiUQC0WhU/D3f1wCwcePGHtclO8jfeustzJo1C+FwuMcDORAIYMeOHaitrcWrr76K\nK6+8EgsXLsQf//hH8YNdffXVcLvd0Ol0cLvdqKmpgcPhgMFggNVqhdVqFb9URqOxKAGVTqcRiUQQ\nDAYRDofR1tYGr9cLv98Pr9crHtStra3iTvd4PGhsbEQqlZJ8f5VKBa1WC51OB41GI34evV4PtVqN\nsrIy8cEYQyaTQUdHBzo6OpBOp5FIJBAKhRAKhcSAkqLX6+F0OlFZWQmbzYaamhrU1NSIr7lcLjgc\nDjidTlitVtjtdhgMBqhUqkLsUlEmk0EymUQsFkMgEEBLSwsaGxvR0tIi7l+fz4dwOIxwOCwe6F3D\nW84+BiCGgslkgsFgQEVFhbh/eajwz8f3czqdFsMonU4jFoshGo2KJwQ5Kioq4HK5MHjwYLhcLlit\nVjgcDthsNthsNvGka7fbYbPZYLfbxeWOJ6R4MITDYUQiEfEY4SHBAyMSiYiBHAgExP3e3NwMj8eD\neDwuuS69Xg+tVouKigrodDrxZ34SLC8vF/evSqUST558H3cNplQqhXg8Lm6XnP3btaBitVrFk63N\nZoPD4YDD4YDVahX3ucFggF6vF/e/Xq8vSHYwxsSMaGlpEQsVfD97PB60t7fD7/ejra0N7e3tWSdI\nOcdyeXl51onPYDDAZDKJ+cEbmEyePBl1dXU530N2kH/ve9/Dp59+ir179/a4zM6dO/H222/j5z//\nOerr6zFmzBicccYZ2LNnDwBg2rRpiMfj2Lx5M6xWK4LBoOR69Xo9dDqdeBbnH06tVmd9WflBxB+8\nxJBIJBCLxZBMJsWDX2rnqtVqOBwOuFwuuFwumM1mOBwODBo0CE6nE3a7HW63G0ajUSwx8y+yVqst\neIktnU6LX2D+xY1EImhvb0dzc7MYjO3t7fB4PPB6vWJ4xmKxvJ+T71P+OXgY8kAEskOQn2w6OjqQ\nSqWyvrD8CiWfiooKOJ1OMYB5yddgMGQFocPhEEto/MvLS886nQ4Gg0E8wAuFf2n5lUAgEEAikYDP\n50NbW5t4Yo1EImhtbcWRI0fELzYvCPDSU0/4SZ1/afnxzb+wjDExCGOxmBh+fP9KvT8gFCT0er24\n7/j+rKqqgtvthsViEfev0WiE1WqFy+WC0+kUQ6RYVx3pdFq8Qg2FQuK+5VcI/DiPRqPi8R0MBsV9\n4ff7ZeWGSqWCwWDoFo78+8kLUgDEY7trCZlfaUsdz+Xl5eJJ2uVyoaqqSixg8NdsNhvMZrNYgHK5\nXGJQy73C65qducj+39q4cSNOO+20rNeCwSDuvPNOXHvttTjvvPNwzz334Pnnn8f06dPFy1qdTicu\nr9Vq4fF4AAil93A4jObmZni9XkSjUfELxEtHvOTBdyrfsbw01fXymV828/8gXmLgXxZ+ljebzeLZ\n3WKxwGw2w+l0wul0wmazwel0wmw296uqCn5i+Xrrn0ceeQSnnHIKxo0bh29/+9s5/zYWi6G1tRVe\nrxft7e3w+XwIBALw+Xzw+XziJTS/nEsmk0in0+L+BY7uWx6evBpIo9GIJTa+f/nJQK/Xw2q1orKy\nEoMHD4bb7Ybdbodery/4VUChqFQq8ZL/WPATAQ8jXt3m9Xrh8XjEKw5+2cwLGbwKg+/viooKsZTJ\nCwtarVYsMPCqNb6tPCT4lUB/2Me8Wkqr1eLQoUPi62q1WjxBH6t0Og2v14tgMAi/3w+/3y9Wu/F9\nzDOE79uuP3etfuPbpNPp4HK5sgoKfN9bLBZUVlbC7XaLV7Q8R8xmc8EKFKtWrcKKFSsQi8UQCoXE\nqhWfz5eVnbnICvJYLIb29vaspoUA8NRTT2HFihXYuXMnNm7ciNraWqTTaYwdOxarVq0CACxatEhc\n3mAwZF1amUymbicHIt+CBQvEf/d0YaXX6zF06FAMHTpUqc06YR3viaCUyCk1Hyu1Wo2qqipUVVUV\nbR19YcGCBVnf6a6+np1fJyvItVotZs+ejUmTJmW9fs011+Ddd9/F/PnzAQD/8z//g0AggKuuugod\nHR146aWXcNVVV4nLW61WBAIBOaskhBDSSSo7ZQV5WVkZXnzxxW6vu1wuPP/881nLLVmypMf3sVgs\nRT1TE0JIKZLKTkUrgg0GA2KxWI/VAIQQQrqTyk5Fg9xoNIIxJnknmBBCyFFS2alokPObQHLb6xJC\nCJHOTkW7qPEenoFAoOTuOPeFZcuW9fUmEJITHZuFJZWdfRbk5Pj11MuLkL5Gx2ZhSWVnn1SthEIh\nJVdLCCEDmlR2UpATQkg/R0FOCCEDXL8KcqPRCEAYppQQQog8UtmpaJCbTCYAVCInhJDekMpORYPc\nYrGgvLwcXq9XydUSQsiAJpWdiga5SqWC0WikDkGEENILUtmp+KDbWq1W1ow3hBBCjsqXnYoHOR/8\nhRBCiHz5slPxIDcajdRqhRBCeilfdioe5Hq9XtYErOSocDhMQ/+SfomOTeXky05Fx1oBhPkIk8mk\n0qsdsBhjcLvdAISJPBwOB9xuN2pqasTJcvl8jUajEQaDQZzpfMqUKf1q7lFSWujYVFa+7JQM8k8+\n+QQ/+tGPMH36dNjtdnGi0enTp2PMmDHdlt++fTteeuklpNNpzJ8/v9ucnGq1WtZM4ETAL6Wi0Sga\nGhrQ0NAg+28zmUyxNosQOjYVli87JYO8vr4eu3btwq5du8TXzGZz1lyc3AcffIDa2lps3boVXq8X\no0ePxrZt2zB27Fhxmb6e3XugMZlMCIfDCIVCaG9vh8fjQUtLC5qamuDz+cSZw0OhECKRCGKxGFKp\nFNLpNO1rUlR0bCor3z6TDPIvvvgCN954IyZOnIjy8nLcd999WLp0KYYNG9Zt2QceeAAOhwNjx44F\nYwwdHR146qmncP/994vLdHR0iKV6Io9KpYLFYoHFYsm53wnpK3RsKidfdkpWUpnNZixduhQ/+tGP\noNPp8K1vfQvXXHNNzmW3bdsGrVYLQPgP1mq12LJlS9YyqVQKFRUVvf0MhBByQsuXnZIl8ptuugmA\nMKD5woULsWPHjh6XbWtrQ01NjfhzWVkZ2traspaJx+Ni2JPj09gInHYaUFkJHDzY11tDCCmmfNkp\n+7bxqlWrUFlZierq6h6XSafT3V77egP2WCwGvV6PpqYmqFSqbg+aWUQ+rRaIxYBgsK+3hBByvOrq\n6nJm4qBBgwAczc5cZAf5X//6V5x66ql5l7HZbFl3ozOZDOx2e9YyiUQCWq2WmiAWAD8504gHpL+h\npuWFw0vhPDtzkRXkHo8H+/btE+eN4/bt24cZM2Zg+/btAIBx48YhHo8DENqYJhKJrBYrfGN0Op24\nHDl2vBkuteQi/c2ttwIaDbBiRV9vycCn0+kAHM3OXGR1CGpsbASAbhXtDz/8MN58803U1NTg8ccf\nx8KFC3HVVVfhyy+/hNfrRSaTEevYuXA4DJPJhFGjRlGPsOPEb2BTs3zS33g8QDoNGAx9vSUDR11d\nXd6qZZ6ducgK8uHDh+P666/H97///azXly5dCrVajUWLFgEAZs6ciZdeegm//e1voVar8fbbb+OC\nCy4Ql0+lUohEIt1K9uTYUImc9FdNTcJzl7YP5DhIZaesIDcYDHjyySe7vV5TU4OHHnoo67UrrrgC\nV1xxRc73CQQCAIS6dHL81GogGj0a6IT0Fy0twnNnD35ynKSyU9GxVtrb2wEATqdTydWWLJUK6OEm\nNiEF9eKLwIcfAiYTYDYDDgdQVQU4nUB1tVDy7lrz+sgjQHMzMGpU321zKZHKTkWD3OfzAQAcDoeS\nqyWEHIdDh4BrrpGuwjMahWCvrBT+bbcD69YBLpdQV242AxaLcDKwWIRlDAbh3zfcIFxZbtggFFBI\nNqnsVDTI/X4/AFAdOSEDyFNPCSE+YQIwdSoQCgHt7UBbm3BTs7lZeEQiwqMXY2d1EwoJwU6ySWVn\nn5TIqWqFkIFj927h+cYbhZJzLpkMEA4L4e71Cv/mIe/1CgEfDgud18JhIBAQXotGhfDmPZNp9I7c\npLJT0SDnE4f21ISGENL/dBYG4XL1vExZmVCStliA4cN7vw6NRmiuSDfuc5PKTkV3G7/zaqFrp4JI\nJIQSDDUCIsXUWRgs6nHGR/eggVFzk8pORYPc7/ejvLycSuQFkkwCqZTwIKRYQiHhWYnyF5XIc5PK\nTsWD3Gaz0aDyBcJHOaAmiKSY+Fg+SgxaStGQm1R2KhrkXq+Xmh4WEC8p0QUOKSY+gKkS3e2pl3Ju\nUtmpeB05NT0snM77HxTkpKj4cWY0Fm8dvKBJQZ6bVHYqGuTt7e1UIi8gfslLVSukmJSoWqFxg/KT\nyk5Fg/zIkSMYPHiwkqssaeeeK9zw3LChr7eElDLeokRdxMbK/L1zzE1DIJ2dipfIq6qqlFxlydNo\nqGqFFBcfJrmYTQP5VeXXJhQjnaSyU7EgT6VSSCaTMBazoo0QUjTFbBpoNgvPvD6eHCUnOxULcj53\np4FGmidkQCpm00AeC5FI8dYxUMnJTsWCPNg5Q7CZn3oJIQNKMSf00miEZ+rc1p2c7FQsyFtbWwEA\nlZWVSq2y5NFMeUQJvCRezCkF+VSUNJVvd3KyU1aQx+NxrFy5EnPnzsUdd9yBaDSac7ktW7Zg2bJl\nePXVV/Hiiy9i8eLF4oDofBhGu93eqw9BekYT3BIlKDE3LG/ayJs6kqPkZKdkg6J0Oo3vfve7sNls\neOyxx3DBBRfAZrPhl7/8ZbdlGxoacNddd4k/L1myRBx2MdH5P6RVop/vCYImuCVKUKuF4yydLl5b\ncj58bTJZnPcfyORkp2SQP/3009iwYQPeffddpFIp+Hy+vINeXXLJJZg8eTKmTZuGs88+W3ydj95F\nPTsLhya4JUpQokTOq1aoRN6dnOyUDPI1a9YAEMbDffvtt7F161accsopPS4/efJkXHzxxdDr9chk\nMijrbLPEK+xpCNvCqawETjlFmDORkGJRorMOjwU+9jk5Sk52Sgb5/v37AQAbN27EOeecg9NPPx2v\nv/46pk2blnP5F154ASeddBLuvfdeOJ1OrF27FkajkUrkeWQywNVXA1u2CJeuFRXCs9kszHtoswkP\nu10Y3J9PfLtwIfCb31CJnBSXEkHOJ77xeou3joGqICXyjs7rqeuvvx7Dhw9HIpHAnXfemTPIL730\nUlx66aVwOBxobm7GLbfcgpUrV+L2229HMBiESqWCyWTCRRddhM2bN+dc37Jly1BXVyfn85WMNWuA\nl18+vvfINfFtTc3RiW//+U/g88+BTZuEUjwhcvGq2WLWX1OHIKCurg7Lly/Pem3Tpk1Z2dkTySCv\nrq5GfX09rFYr1J2n5i+++CLnsqtXr8aBAwfw8MMPo7rzev/jjz8GIFTNGI1GqFSqHlu9nKjuuUd4\n/v3vgdmzhXrCREKY39DvF2ZoCQaFm5sez/FNfLt+PfCTnxT/M5HSwYO8mE0DqYt+biaTKSs7eyIZ\n5BMnTsSmTZsQCoXEdoxDhw4FALz33ntYuHAh/u///g/nn38+/vSnP8Hr9WLFihXiHHODBg0CAEQi\nEbGLafhEPu3m8NlnwvPcucCx1DzJmfj2d78Tls037yIhuSjRooT3PqdoyGYwGLKysyeSQX7TTTfh\nkUcewdq1azF+/HioVCrcfvvtAIAXX3wRu3btwpo1a3D++edj7ty5CIfDaGtrwyuvvIKqqir87Gc/\nAwBEo1Gxi+m+ffuO9/OVFF73eKzD0cqZ+HbrVmDjxuKOKU1KE7+iL2b3eV7A6Oz7ckKqq6vLWa3c\nNTt7Itkh6OSTT8bWrVvh8Xiwbt06bN68GVdddRUA4Fe/+hUefPBB3HHHHQCAO+64A6NGjcLy5ctx\n8cUXY/fu3WILl1gsBj0NnJ2TEqPL8S8hBTnpLR7kxSwt8yDnEz2To+Rkp6wRhkeOHInf8WvzLpxO\nJxYvXiz+rFKpMGfOHMyZM6fbsul0WqxjJ7kVc3Q5JafrIqWFV60Us40377To8RRvHQOVnOxUdBhb\nDR8Zh+RUzNHl+PyeNGYZ6S1eGCzmzc7OW2lobCzeOgYqOdlJQd6PFHMQLApycqyU6KzDq1aoHXl3\n/SrIqWqlZ0qMLsdLUzTUDektPjFNS0vx1sHv3UQiNG/n1/WrqhXGmNhdn2RTYiwL/t50LiW9xaeK\nLGa1R1nZ0XryzgFTSSc52alosuZr0H4iU6ILtBItY0hpcruF52I3DeTDbdMNz+6kspOKyP2AEiVy\nfpKgICe9xUvKxR7QyuEQninIe0/RIGc0pU1OSpTIOardIr1VXS30OC52wPIbnlS10p1UdipWY6pS\nqcQBuEg2JYOcardIb40ZA/gPhxErMwIo3gHkdApXjNQpKJuc7FQ0yKlEnpsSo8txmQyVyknvqMsZ\n4HZDDxwdR9ntFobXdDqFXmZWq/AwGoWuoCaT8G+dTvi32Sz8W6MRSi4qlXAgqlRCu1vG8MSfGZ54\nnNEd+a+Rk52K7bGysjJkqF1RTkqMLse/LxTkpLfKYp3jO0SjwvCaUkNsHi8q8GWRk52KBblarUZa\nibqDAUiJ0eXKyoSbqXQuJb3FjCbEWsMwdISECmyPR2hU3tQk1INEo8Kd0FBIaAgeDguPSEQYG4L/\nHIsBqZRQh5grrFWqoyUOqgMUyclOxYJcq9UiSTOr5qTE6HJqtRDk6fTREwchcqhUQHWNCiqVBYcO\nWWAZNqwo6+FXjMhQ66qu5GSnYhfZBoOBJpTogRKjy9HA/eR4GI3C5CbFPEZfe00ocMyaVbx1DERy\nslOxIOczXZDulBhdjqbSIseDT3hSzBYlNLlEbnKyU7EgNxqNFOQ9UGJ0OT58bTGrb0jpstmE5855\ngIuCH6N01ZhNTnYqFuRWqxWhUIjakuegxOhyfPC0VKp46yCliwd5MUvk/KqRj9RJBHKyU9GqFQBU\nT56DEqPL6XTCczFL/aR08W76xSyRU9VKbnKyU1arlcOHD+Pee+9FNBqFw+HA3XffjYocTR9SqRTu\nu+8+HDx4EJFIBIsXL8a3vvUtAICuM0ni8TjMNCh2FiVGl+Nt1YtZD09KlxJXjUo0wx2I5GSnZIm8\nubkZEydOxKBBg/Dzn/8cf/7zn/HCCy/kXPauu+7Ck08+iVWrVmH27Nm46KKLcPjwYQAQJw+lEnl3\nSowuR18Scjx41UowWLx1UMuq3ORkp2SQ33fffTh8+DD++7//GyqVCjU1NTjppJO6LccYw6OPPoqR\nI0dCrVbjjDPOQCKRwGuvvZa1MTH6X+pGidHleNUKlcjJsVCi1QpdNeYmJzslq1Zef/11VFRU4L33\n3sPnn3+OV155BSNHjuy2XGtrK9ra2sQqF/788ccfAxAq7AHAS3M5daPE6HJKXBqT0qVEHTkfOoJ6\nH2eTk515S+SMMRw8eBDJZBKxWAyVlZX4xje+gQ8++KDbsnwlfAB0PqNFe+eYlI7OwYb9lCTd8NHl\n3t1YvDEmXC7hi0Ijy5FjwXuIF7Nqjvfap4Zt2eRkp2SJnI+6NXv2bNTX1yMej+PBBx/Ec889J+vv\n+GAvls4iYTAYxKOPPoobb7wx598tW7YMdXV1UptVUpQYXe7hlQwPr2Cd3xYaXY7IxxiwaJHw79Wr\ngR07ijL4IXiknH56333WvlRXV4fly5dnvXbzzTfjpptuAiBkZ0/yfqNVKhWGDBmCxsZGaLValHcO\ngMBvYHZV1dmGjrd15M/8dX55EAwGxVAnAhpdjvRnX+9EtnNncde3ZElx338g8fl8WdnZE8mi2aRJ\nk/D0008jGo2KI3CNGDECAPDoo4/iN7/5DdauXYtvfOMbGDFihHhnlT/z5odd63lOPvnkY/1cJYlG\nlyP9mckkVHcEg4DXW9zDc9Ik4PvfV/4z9leBQEDe/UUm4dNPP2Umk4n98Y9/ZL/97W+ZzWZje/fu\nZYwxNmvWLAaAPf7444wxxv72t78xu93OPvvsM7ZkyRI2cuRIFgqFxPfSaDRsyZIlUqs8IZnNjFks\njAUCxVtHR4fwIKSYfD7GNm5k7NNP+3pLSodUdko2Pxw9ejR27dqFSCQCm82GTz/9FKd3VmI99thj\nePXVVzFnzhwAwHXXXYc33ngDL7/8MoYOHYqtW7eKvZIAoVRONztzK/bocmvXCkODXn55cd6fEG7z\nZqC2FvjFL/p6S0qHVHbKuus1bNgw3Hbbbd1ed7lcmDlzZtZrEydOxMSJE3O+j81moyDvgdUKNDcL\nl6qDBhX+/XnzMZqhnBRbc7PwXFPTt9tRSqSyU9FJv4xGI/Xs7EGxR5fjJ4diDgNACCDUjWs0QosW\nUhhS2aloOzSr1YpAMXsUDGDFHl3O5RKeqT8WKbbFi4UHNY4qHKnsVLREbrFY8jahOZEVu+ccH1ku\nEqGec0QZ1DCqcKSyU9Egp1mCelbsLvRlZUdPFp2dbQkhA4RUdlKQ9xNKjC5XWSk80w1PQgaWfhXk\ndrsdHo9H7L5PjlJidLnOIRsoyAkZYKSyU9Egdzgc4gBcJJsSo8vxG55UtUKKKRSiKQULTSo7FQ1y\nW2f9gY+G4OuGTzxbzNaZvEROLVdIMU2eLExksm1bX29J6ZDKTkWD3N05FU5LMSenHKCUmOVeiVI/\nIbzQyCczIcdPKjsVDfJBnb1Sjhw5ouRqBwQlJkfmJfK2tuKtgxB+T67L6BzkOEllZ5+UyNsoSbrh\nB30xG/VQpyCihFBIeKY51gtHKjsVDXI+HCP17uyOtyMvZvND/sXiXzRCioFfVVLVSuFIZaeiQa7v\nnCabWq10p0SJnJ8sKMhJMWk0ws1Ojaavt6R0SGWnomOt8AmZk8Wc+G+AUiLIeTd96pNFiokGOC08\nqexUtEReVlYGvV5PvTtzUKL5IV8HXRARMrBIZaeiQQ4IwzFGitnGboDSaoXnYl6sUB05IQNXvuxU\nPMi1Wi1VreSg7qzkKmaPOKpaIWTgypedknXkzc3NmDVrFi666CLodDrs3LkTCxcuxLRp07otu337\ndtx6662YOHEiVCoVtm7dioceeghnnnlm1sbEi9lYeoDiQ34Wcxiazmq2opb6CSHFkS87JYM8mUzi\ngw8+wAcffIDy8nL84he/wNSpU3MuG41GsWHDBmzYsAE6nQ733XdfVogDQqU9lci7K+u8NirmWOGd\nN76pjpwUVTwuHMd6PY1JXkj5slNWq5WJEyfi//2//4dx48aJPYx6csUVV2DBggUYP348nDnmetJo\nNEjRiDrdlClQycXr4ROJ4q+LnLiMRiHIU6mjVYbk+OXLTlnxYbFY4PF4sHbtWnz55Zd5l7Xb7Th0\n6BBef/11NPNZWLsoLy9HR0eHnNWeUJQouShR6ieEH8t0nBVWvuyUDPKysjI0NDTg9NNPx1dffYUz\nzzwT69ev73HZ/fv3Y8KECVi/fj3OPPNMfPbZZ1nLqNVqpNPpY/gYpU2JG5C8/p3Oo6RQwuFwtzGy\n//IX4P77qVql0PJlp4rJmOUhlUpBo9Fg9+7dGDt2LCZOnIj333+/23KMMXR0dECtVmPt2rW4/PLL\nce211+KZZ54Rl5kwYQIsFgvWrl0LLb/W72LZsmWoq6vrxccb+BgTOgTxNuTjxgFuN1BTI8xEbjAI\nE09YrcJlq8kkPIxGoRu0ySQ0LdTphN50arXwJSorE54ZEx5PPAHMnQuMGQPs3t2nH5mUAMYYTJ09\n2VwuFxwOB9xuN2pqauB0OmEwGGC1WmG1WmE0GmEwGKDRaKDRaDBlyhSUKVGfOIDU1dVh+fLl3V7n\nNzl5dr7zzjvdlpGswVq1ahXuvfdebNy4Uezvf+DAgZzLLl26FC+//DI+/PBDWDr7g399WX5WUVPl\nmejrTUN37izu+pYsKe77kxMDb9McjUbR0NCAhoYG2X+boXoX2XgpPF+JXDJN33jjDfj9fmQyGfg7\n+96OGTME1bQzAAAfkUlEQVQGALB+/XpcdtlluPvuu7F48WLs2rULwWAQjDF4OucTGz16dNb7lZWV\nIZPJ0Nm4C5NJqFoJhYTZezweoKUFaGoSpn6LRoVuz6GQEPrhsPCIRIQWKPznWEy4wZRO527GqFIB\nkyYB3/++8p+RlB4+j2QoFEJ7ezs8Hg9aWlrQ1NQEn8+HaDQKv9+PUCiESCSCWCyGVCqFdDoNFdW7\nyMZPejw7c5EM8vvvvx8/+9nPsHfvXjz11FOoqanBgw8+CADYv38/YrEY9uzZAwD4xS9+gYcffhjv\nvfce/vjHP2LUqFG46667st6v638gzd15lEolDGplsQDDhkkvHwwC//mPMGnziBHF3z5CclGpVLBY\nLLBYLBgm58AlPaqrq8tbrZzv5CerjvzIkSNYt24dnE4npk6dKtaLMcawfft2jBkzBrrOMSvr6+ux\nadMmDB48GJMnT+5WD/7tb38bKpUKGzdulPHRSE9efx2YOROYPh14662+3hpCSLHly05ZFdWDBw/G\n3Llzu72uUqkwfvz4rNdGjhyJkSNH9vhemUyG6scLgLfsrKnp2+0ghCgjX3YqXlFNNzoLIxIRWqjk\n6HNFCClB+bJT8USlIC+MxYuFB91mIOTEkC87FS+RJ5NJcZB0cvzo5j8hJ4Z82al4kIfDYfFmKSGE\nEHnyZafiQR6LxWDgU9UQQgiRJV92Kh7k8Xg8Z9d8QgghPcuXnYoGOWMMoVAIZj7nGDlmoVBxZxMi\nhPQfUtmpaJBHo1Gk02nY7XYlV1uSJk8WZvzZtq2vt4QQUmxS2alokPt8PgCgIC8APstPZ4daQkgJ\nk8pORYM8EAgAgDiKIjl2fPxyagBESOmTyk5Fg5yPiEgl8uMXCgnPdLuBkNInlZ2KBrnX6wUgDEJP\njg+fTJuqVggpfVLZqWiQt7e3AwAcDoeSqy1JGo1ws1Oj6estIYQUm1R2KjroSVtbGwCgqqpKydWW\npM45PgghJwCp7FS8jlyv10Ov1yu5WkIIGdCkslPRIG9qaoLb7VZylYQQMuBJZadkkOeaQEhqUqGe\nft/U1IRBgwZJrZIQQkgXUtkpWUe+efNmLFiwAOeddx4AYPv27XjhhRdyzgIUi8WwYMECJBIJtLW1\n4ZJLLsGtt94qzjXn9/sxePDgY/0shBByQpLKTskgT6fT2LNnD/bs2QO3242HH364x6nc7r77bvzj\nH//AoUOHsH//fowaNQrjxo3DJZdcAkCosB87duwxfhTSVTwOZDKAXk9jkhNS6qSyU1Yd+fXXX49P\nPvkEDQ0NuPrqq3tc7tlnn4Xb7UZZWZl4GfDqq6+Kv/d4PNSGvECMRuHR0dHXW0IIKTap7JQV5Mlk\nEmvWrMEDDzyA+vr6nMuk02kcOHBAnIqIP//nP/8BIEwcGovFaFKJAuGl8Eymb7eDEFJccrJTMsi1\nWi0ymQwWLVqEQ4cO4Rvf+AY++uijbstFo1FkMhmxPrysTHhr3rWU90yi7vmFUV4uPFOJnJDSJic7\nJYP8/PPPxzPPPAO73Y7vfve7iMfjeOCBB7otx0vgvMVKprOoyF9vbm4GAFRXV+ORRx6BSqXK+air\nq+vFRzxxUZATUlrq6uq65eGiRYuysrMnkkE+b9481NTUwO/3w2g0AgAOHz4s/j7WOZ6qwWCAzWZD\nR2ey8GdeV951rADe3ZQcOz5RSCLRt9tBCCkem80ma4wqySD/6KOPMHjwYGi1WrS0tAAAzj33XADA\ngw8+CIPBgL/+9a8AgMmTJyMYDAI4OuxibW0tACASiQAAjEYjwnwMVnLMrFbAYKAgJ6SUGY3GrOzs\nEZPw4Ycfsu985zvsscceY+PGjWOTJk1iPp+PMcbYXXfdxQCwP/3pT4wxxvbu3cvcbjd7+OGH2fz5\n89l5553HIpEIY4yx5557jgFgn376qdQqCSGEdJKTnSrGJLppAmhtbcW2bdswZMgQjB07VryhyRjD\n4cOHcdJJJ4nLhsNhvP/++zCbzTj33HPFOvI///nPmDdvHg4ePIiTTz75mM9QhBByIpGTnbJGP6yq\nqsKMGTO6va5SqbJCHABMJpPYAagrXp1CzQ8JIUQ+Odmp2KBZvM7cYrEotUpCCBnw5GSnYkEeDoeh\n1WrFqhZCCCHS5GSnYkEei8VgMBiUWl3Jy2SEeTs7W38SQkqUnOxULMiDwSDMNFNwwSxYAFgswBNP\n9PWWEEKKSU52KhbkoVCIgryAnE7huXMGKEJIiZKTnYoFeSKRgJZ3RyTHjXfy8vn6djsIIcUlJzsV\nC/KOjg660VlAVqvwTEFOSGmTk52KBnk5H+mJHDdetdI5uCQhpETJyU7FgjydTlOQF5DTCVRUAHSR\nQ0hpk5OdisUAY0wco5wcvwsuEKZ7o2neCCltcrJT0fKcilKnYOicSMiJQyo7FY2DDM1LRgghvSaV\nnYoFeVlZGQU5IYT0kpzspCAnhJB+rF8FuVqtRjqdVmp1J4SODuDIEWHMFUJIaZKTnYoFuV6vF+f3\nJIVxxRXAkCHAu+/29ZYQQopFTnb2Osij0WiPv4vH4z3+bDabaa7OAuO9O6lTECGlS0529irI3377\nbYwcObLH3y9duhSjRo3CzJkzcdFFF2Hy5Mni7wwGgziJKCkMt1t4poGzCCldcrJTdpC3t7djzpw5\nkpXu9fX1eP3112G327FmzRrxdY1Gg1QqJXd1RAY+8xOdHwkpXXKyU1aHIMYYbrnlFuj1+m7VJ1/3\n4osvYvr06TAajVmv63Q6yb8lvaPTCc+0WwkpXXKyU1aQ/+Uvf8HFF1+ML7/8Evv378+77KZNm/Du\nu+/CYrFg4cKFqK6uBgBotVokEgmZm07ksNuF4Wx5oBNCSo+c7JSsWtm/fz/ee+893HDDDZIrHDRo\nEC688EKsWLECL7/8Ms455xxx4tCKigowxqgJYgEtWCDUjy9b1tdbQggpFjnZmTfIU6kUlixZggce\neEDWOClz587FZZddBpVKhdraWjQ2NuKll14CIFweAMIg6TNmzIBKpcr5qKur68VHJISQ0lBXV9ct\nD996662s7OxJ3iDfunUrPv/8c/zgBz/A9OnTsWPHDng8HsyePRsAkEwmcejQIQDCBKGnnnoqLr/8\ncgDCWQQA2jqbVFg728r5/X7Y7fbj+byEEHJCCAaDWdnZk7xBfsEFF2DXrl1Yt24d1q1bB6PRCKfT\nKZayr7vuOpx88sn417/+BbVajY6ODrHJ4ZEjRwAAF110EQDA2TkTgtfrhcPhOM6PRwghpa+9vT0r\nO3siq/nh7t27cd1118Hv98Pn82HevHkAhHBWq9Uwm83QaDRYuXIldu7ciTvuuAP/+Mc/sGLFCkyc\nOBEAYOpsKxcOh7Fy5UowxnI+qGqFEHIiqqur65aHN998c1Z29kTFGGNSK0ilUlkN0lUqFaxWKxhj\nSKVSYjUKIFSlHDlyBCNGjMhqgvjuu++itrYW//rXv7I6CpFjl04f7dXJOwcRQkqLnOyU1fxQo9HA\nZrN1e12lUmWFOABUVlaisrKy27I81Kl3Z+HU1wOjRwOjRgH79vX11hBCikFOdio2aBYP95aWFqVW\nWfIMBuE5z/A3hJABTk52KhbkvMLe5/MptcqSZzYLzzQWGSGlS052KhbkBoMBGo0GHhqqr2B4j04a\nHZiQ0iUnOxWdIWjQoEFobGxUapUlT6MRnmksMkJKl5zslHWzs1CcTmfetpCkd8rLgepqQK3o/yIh\nRGlS2aloBJhMJoRoXrKCKS8Hmpr6eisIIcUmlZ2KVa0Awlmlvb1dyVUSQsiAJ5WdigZ5VVUVBTkh\nhPSSVHYqGuR2ux1er1dyliFCCCFHSWWnokFeWVmJdDpN9eSEENILUtmpaJDzUQ+pLXlhZDJAYyPQ\n3NzXW0IIKSap7FS01QofryXfuLpEvkgEGDxY6KpPQ9gQUrqkspNK5AMYn4+V5uwkpLRJZafideT5\nNob0Dh9jhY+5QggpTVLZqWiQ8wHS6WZnYQSDwjMFOSGlTSo7FQ1yi8UCQJiHjhw/PnwtH86WEFKa\npLKzV0He0dGBw4cP510mGo1i9+7dOccFoMklCkulAk46iWYHIqTUSWVnr4L8nnvuwXXXXdfj79et\nW4exY8di8+bNmDJlClatWpX1e7VaDY1GgyjNhFAQEyYADQ3Aa6/19ZYQQopJKjtlB/nWrVvzToyc\nTqfx4x//GFdeeSVuuukmzJs3D//7v/+LQCCQtZzZbKY6ckII6aV82SkryCORCO6++25YrdYel9mx\nYweam5vFZex2O1KpFNavX5+1nNFopKoVQgjppXzZKSvIlyxZgmXLlkGr1fa4zFdffQUAKC8vF964\nTHjrgwcPZi1XUVGBZDIpZ7WkUzgcBmOsrzeDENKH8mWnZM/O1157DSeddBK++c1v5l2O192oVCoA\nR4P86z2R9Ho9YrEYfD6f2Mi9q2XLluWtwjnRMMbg7ryb6XK54HA44Ha7UVNTA6fTCYPBAKvVCqvV\nCqPRKE4LpdFoMGXKFPH/gRDSv9XV1WH58uXdXnc4HPB4PGJ25pI3yJubm/G3v/0NzzzzjORG6Dq7\nF/KSIx+lS/e1bodarRaJRCJv6Z4cxS+lotEoGhoa0NDQIPtvaZRJQga+RCIB4Gh25pI3yLds2YLW\n1lZcdtll6OjoQHNzM2KxGBYtWoSVK1eiubkZbW1tOOuss3DqqacCEG56dn0+5ZRTst6zrKwMmUyG\nSooymUwmhMNhhEIhtLe3w+PxoKWlBU1NTfD5fIhGo/D7/QiFQohEIojFYkilUkin0+LVESFk4OIF\nMp6dueQN8tmzZ2P27NkAhHpas9mMs846CytXrgQAXHjhhfjiiy/w6aef4uyzz8aQIUPQ2toKAGht\nbYVOp8PUqVOz3rO8vBwdHR3Q6XRU7yuTSqWCxWKBxWLBsGHD+npzCCFFUFdXl7damWdnLrKKxW++\n+SZmz54NvV6Pzz//HIsWLQIAnHfeeRg8eDCqqqqgVqvx9NNP45133sGvfvUrPPbYY3jiiSfgcrmy\nV5jnrEIIISS3fNmpYgUuFnd0dMDj8cDpdIotWLqaNGkS1Go1NmzYUMjVEkJIScuXnQUfj7y8vBxV\nVVU9/r6srIyqVAghpJfyZafidxzpBhwhhPRevuxUPMipNE4IIb2XLzsVD/JMJkOlckII6aV82dkn\nJXIKckII6Z182al4kKfTaajVis75TAghA16+7FQ8yFOpFDQajdKrJYSQAS1fdioe5PF4HHq9XunV\nEkLIgJYvOxUP8lAoJE4kSgghRJ582al4kIfDYQpyQgjppXzZSVUrhBAyAPSbqpVMJoNYLCbOCE0I\nIUSaVHYqGuR8kgSqWiGEEPmkslPRIA+Hw3k3hhBCSHdS2dknJXKqWiGEEPmkslPRLpaBQAAAYLFY\nlFxtyeo6mwhNWE36Ezo2C0sqOws+sUQ+69evx5QpU7BhwwbU1tYqtdqS1XXcBRpVkvQndGwWllR2\nyqpaicVi2LBhAzZv3ox4PN7jci0tLTh8+LD48969e8VLAkCYCR4ADAaD3O0nhJATnlR2SgZ5a2sr\nJk2ahIMHD+LOO+/EsGHDsGfPnpzLrl+/HieddBKqq6vhdDoxa9asrDOz1+sFADgcjl5/EEIIOVFJ\nZadkkD/77LPYvn07QqEQZs2ahaamJjz55JM9Lm+xWGA0GvFf//Vf2Lx5c9YZpKmpCQBQU1PTqw9B\nCCEnMqnslLzZefHFF2PmzJmYNGkS1qxZAwA466yzelz+zjvvxC233JLzdy0tLTAajdRqpUB++tOf\noqmpCe3t7X29KYRkCQaDSCQSSCQSfb0pJeG2227Dr3/962NvtXLWWWfh+eefx8svv4xnnnkGt956\nK6699toel3/ttddw6NAh1NfX44c//CGuu+468XderxdOp/MYPkb/kMlkEI1G4fF40NbWhkgkglgs\nhnA4DI/Hg0AggHg8jmQyiUQigXg8jlQqhWg0ilAohFgshnQ6jUwmg0wmI76vSqWCWq1GeXk5ysvL\nodVqYTabYbFYoNfrYTKZYLVaYTKZYDabYTabYTQaMX/+fFRXV5dcVVUkEoHf70cgEEAwGITX60Ug\nEEA0GkU0GkU8Hkc4HEYgEEAkEkE4HEY0GkUikRD3L2Os2002vp81Gg00Gg3UajX0ej0MBgOMRiNM\nJhMsFgusVqv4rNVqYbPZUFVVBYvFAoPBMODH008mk/B4PAgGg4hEIggEAmhra4PP50M0GkUwGEQo\nFEIikUAymUQ8HkcsFkMikUAqlUIqlUJHR0fWMQwI+7esrAxqtRoVFRXQarXQarXQaDTQ6XQwGo0w\nm81ZxzHft06nExaLBRaLpWSGuWaMIRgMIhAIIBQKIRAIwOPxwOPxIBwOIx6Pi8czz4xEIoFoNCr+\nnu9rAHmzU9YRqVKpYLfbMWLECPz73//GV199hWHDhnVbbvjw4bjlllswc+ZM/OxnP8MPfvAD2O12\nfOc73wEgnKX5xtxxxx2wWq2oqamBw+GAwWCA1WqF1WoVv1RGoxFlZYVv6p5OpxGJRBAMBhEOh9HW\n1gav1wu/3w+v1yse1K2treJO93g8aGxsRCqVknx/lUoFrVYLnU4HjUYjfh69Xg+1Wo2ysjLxwRhD\nR0cHEokEOjo6kE6nkUgkEAqFEAqFxICSotfr4XQ6UVlZCZvNhpqaGtTU1IivuVwuOBwOOJ1OWK1W\n2O12GAyGgs/WlMlkkEwmEYvFEAgE0NLSgsbGRrS0tIj71+fzIRwOIxwOiwd61/CWs48BiKFgMplg\nMBhQUVEh7l8eKvzz8f3MvxypVArpdBqxWAzRaFQ8IchRUVEBl8uFwYMHw+VywWq1wuFwwGazwWaz\niSddu90Om80Gu90uLnc8IcWDIRwOIxKJiMcIDwkeGJFIRAzkQCAg7vfm5mZ4PJ68DRY4vV4PrVaL\niooK6HQ68Wd+EiwvLxf3r0qlEk+emUwG6XQ6K5hSqRTi8bi4XXL2r9FohMPhQHV1NaxWq3iytdls\ncDgccDgcsFqt4j43GAzQ6/Xi/tfr9QXJDsaYmBEtLS1ioYLvZ4/Hg/b2dvj9frS1taG9vT3rBCnn\nWC4vL8868RkMBphMJjE/ysvLoVKp8hbYJIM8GAziyJEjuPTSSxEOh3H11Vdjzpw52LRpU7dl7Xa7\n+MU5++yzAQAvvPCCGOTPPPMMJk+ejNraWnz00UdZLVp6otfrodPpxLM4/3BqtTrry8oPIv7gJYZE\nIoFYLIZkMike/FI7V61Ww+FwwOVyweVyoaamBqNHj8agQYPgdDpht9vhdrthNBrFEjP/Imu12oKX\n2NLptPgF5l/cSCSC9vZ2NDc3i8HY3t4Oj8cDr9eL999/Hy0tLYjFYnk/J9+n/HPwMOSBCCDrC8pP\nNh0dHUilUllfWH6Fkk9FRQWcTqcYwBaLBUOHDoXBYMgKQofDIZbQ+JeXl551Oh0MBgPKy8sLup/5\nl5ZfCQQCASQSCfh8PrS1tYkn1kgkgtbWVhw5cgRtbW3Yv3+/WBDgpaee8JM6/9Ly45t/YRljYhDG\nYjEx/Pj+lXp/QChI6PV6cd9ZrVZUVVVhzJgxcLvdsFgs4v41Go2wWq1wuVxwOp1iiBTrqiOdTotX\nqKFQSNy3/AqBH+fRaFQ8vpubm8V94ff7EQwGZe0Dg8HQLRz595MXpACIx3bXEjK/0pY6nsvLy8WT\ntMvlwvDhw8UCBn/NZrPBbDaLBSiXyyUGdW+u8JYvX97j7yTfYdKkSfjkk0+wZcsWsXtofX09AKHb\n6BtvvIEZM2bAZDLhqquuwn/+85+sgO66I3Q6Hd5//33x53A4jObmZni9XkSjUfELxEtHvOTBdyrf\nsbw01fXymV828/8gXmLgXxZ+ljebzeLZ3WKxwGw2w+l0wul0wmazwel0wmw2F+VK4FjxE8uxVKHE\nYjG0trbC6/Wivb0dPp8PgUAAPp8PPp9PvITml3PJZBLpdFrcv8DRfcvDk1cDaTQascTG9y8/Gej1\nelitVlRWVmLw4MFwu92w2+3Q6/X9ds5WlUolXvIfC34i4GHEq9u8Xi88Ho94xcEvm3khg1dh8P1d\nUVEhljJ5YUGr1YoFBl61xreVhwS/EujP+1itVosn6GOVTqfh9XoRDAbh9/vh9/vFaje+j3mG8H3b\n9eeu1W98m3Q6HVwuV1ZBge97i8WCyspKuN1u8YqW54jZbC5YgWLVqlVYsWIFYrEYQqFQVtUKYwzL\nli3r8W8lg7yiogJutxvV1dXYu3cvAGDq1KkAhB5bv//973Hbbbfh3nvvxfDhw3HppZdCo9Fg3759\nAIBrrrmmx/c2mUw47bTTevVhSe/o9XoMHToUQ4cO7etNKXnHeyIg8qjValRVVaGqqqqvN6WgFixY\ngAULFhzT30r27Pzqq6+wePFilJWV4bPPPsOFF16IBx98EBaLBW+//TbmzZuH1atXY+rUqWhoaMBP\nf/pTlJWVYf/+/fjpT3+K+fPnH9OGEUIIkUd2F/1MJiO7uoEx1m8v7chRHR0dSCaTNNEH6Vf4TV2j\n0TjgWwgpRXZFcG/qjCnE5dm9ezdWrFih+Hq3bduGVatW4cILL8Svf/1rxddPBoa+uJr+6KOPcNll\nl+GGG25AZWUlli5d2q2ZI+mu/9zROwF9+OGHaGxsVHy9GzZsQHNzMz744AMa0IjkdPDgQWzfvl3x\n9c6bNw979uzB008/jR//+Me4++67sXr1asW3Y6Ch65Y+0tDQgPvvvx+1tbU4cOAABg8eDK1WCwDw\n+/3YvXs3KisrMWrUKKhUKrGVCSC0nw6FQgAArVYLo9EojsVQVVWFxsZGhMNhjB49Gn6/H42NjRgz\nZozYhvmXv/wlPv74Y9x111198MlJfxeJRHDbbbchkUjgwIEDYtNaQOhMtHPnTqhUKowbNw4ajQbx\neFwskGi1WjDGkEwmAQhdynn38srKSrS2tiIcDmPQoEEwGAzYt28fTjvtNFitVgDCEB+NjY3IZDIY\nP348AGDTpk2YO3eu0rthQKESeR9Zvnw56uvrsW3bNvzhD39Aa2srAODJJ5/EmDFj8MUXX2Dp0qWY\nPn06IpEIPvzwQ8ycORPDhw/H2rVrMX/+fAwfPhyrV6/Gtm3bMHXqVNx+++34+OOPcfnll+Occ87B\nQw89hKeeego//vGPMWnSJCp9E1mefPJJrFu3Di0tLfjDH/6Abdu2AQA++eQTjBgxAu+88w7efPNN\njBo1Cnv27EFTUxNuvvlmDB8+HHV1dXjssccwfPhw3HTTTTh8+DDmzJmD66+/Hrt378ZPfvITjBs3\nDr/5zW9w991348EHH8TQoUPR0NAAQLhaPHz4MCwWC7766isA+YcEIZ0Y6RMHDx5kANgvf/lL8bXG\nxkZWUVHBFi1axBhjbP/+/QwAu/POOxljjD377LMMAHv++efZ9u3bGQD2yiuvsI6ODjZp0iSWyWQY\nY4w9+uijDAC75557GGOM/fznP2cA2KFDh8R18b9fsmSJUh+ZDCAjRoxgZ599dtZrZ599NjvjjDMY\nY4xlMhl26qmnsvPPP58xxlhzczMrKytj8+fPZ5lMhjkcDnbzzTczxhibM2cOq6+vZ4wxduDAAQaA\nTZkyhTHG2FtvvcUAsNWrV2etKxaLsVNOOYXV1tayWCxWzI9aEqhE3o9s2rQJyWQSbrcbAMTnf/7z\nnwCAyy+/HEajEX//+9+xfv16mM1m/P3vf8eWLVtwwQUXdLvJzC9NeasUfrlLSG95PB58/PHHqK6u\nBiA0aHC73Xj//fcRiUTgdrsxZcoUvPjii/joo4+QSqXwwgsvIBKJ4NChQxgxYkTW+51zzjkAej42\nb7/9dowfPx5vvPEGDbwlAwV5P9DQ0IDPPvtM/Jnfpf96d2yDwYArrrgCb775JtasWYO6ujq8/vrr\nePzxx7MGJ+P6U+9UMjB1dHTgnXfeEQsJXY9J/m/+u2uvvRZerxfz58/HypUr0drailtvvRUzZszo\n9r75js1nn30WsVgMzz33HL766issXbq0kB+pJNE3vY/wG5vpdBo7duzA5s2bUVtbC71ej/379wMA\nvvjiCwDAd7/7XfHvrr32WiQSCYwePRo/+MEPEI/H8fHHH1M9IikorVYrjj/y2GOPweFw4Fvf+hYO\nHDggjk3y5Zdfora2VpxzYNasWaioqIDf78cNN9yAb37zm3j00Ufz9u7+ul27duHGG2/EF198galT\np+LKK6/E8OHDi/UxS0df1+2cqDKZDLvyyivZiBEj2OzZs9mBAwcYY4y98sorbMiQIezWW29l5557\nLvvhD3/I4vG4+HeJRII5HA72r3/9izHG2LRp09hvf/tb8fdbtmxhEyZMYADYtGnT2BNPPMHOOuss\nBoDNmjWLeb1edu+997IZM2YwAGz06NHsRz/6EWttbVV2B5B+7aGHHmIGg4HNmzePPffcc4wxxurr\n69m4cePYtddey773ve+x8ePHi8ctN2vWLHb77bczxhj73e9+xy6++GLxdz6fj82aNUs87p544gl2\nySWXMABswoQJbPPmzezMM89kALIer732mnIffIBSdPJl0t2RI0fE4Tm5VCqFQ4cOwW63w263d/ub\ntrY2uFwuqFQqBAIBcVAwAIjH42LTREAYqKzrsKUOhyPnCJAOh6PgowmSgc3n8yEWi2HQoEHia4wx\nHDlyBCqVCoMGDep2XyYQCIhDOCeTSUSjUfHYzmQy8Hg84rJfPzbNZjPC4XC31lU2m61kxigvFgpy\nQggZ4KiOnBBCBjgKckIIGeAoyAkhZICjICeEkAGOgpwQQgY4CnJCCBngKMgJIWSAoyAnhJAB7v8D\nwxDSaiWzjIYAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 29 }, { "cell_type": "markdown", "metadata": {}, "source": [ "** Observations: **\n", "1. since most statistic tests are based on a normal distribution assumption (e.g. t-test), we need to notice whether the population (or sample) is roughly normally distributed \n", "2. generally, if the boxes of each distribution overlap, and you haven't taken something on the order of a buttload (metric units) of measurements, you should doubt the differences in distribution averages, like in the case here\n", "3. however, the conclusions of the first two should consider the influcence of outliers, if any" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "** Run a statistical test**\n", "- In statistics, what we are asking is whether differences we observed are reliable indicators of some trend or just happend by lucky chance.\n", "- Statisticans usually do this by comparing the sample statistics (effect size - the average difference in this case) with the standard error and use a model, e.g. a t-test to estiamte the probability of this happening by pure chance\n", "- Different flavors of t-test - paired or unpaired:\n", " - Paired v.s. Unpaired\n", " 1. datasets are paired (also called dependent), e.g. when you measure the same set of itmes twice, usually before and after some change. \n", " 2. datasets are unpaired: e.g., random samples from differet/same populations\n", " - Equal v.s. Unequal size: whether the sizes of the sets are the same\n", " - Equal v.s. Unequal variances: whether the variance of the two samples under comparison are equal or not. If we made no assumption about the sizes of the datasets\n", " - In this case, we are running an *unpaired, unequal size, unequal variance* test. That's **Welch's T-Test**" ] }, { "cell_type": "code", "collapsed": false, "input": [ "## scipy implemented paired/unpaired\n", "## scipy.stats.ttest_rel implements paired test\n", "## scipy.stats.ttest_ind implements unpaired test\n", "%pdef stats.ttest_rel\n", "%pdef stats.ttest_ind\n", "print \"Welch's T-test p-value\", stats.ttest_ind(town1_hts, town2_hts, equal_var=False, )[1]" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " \u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mttest_rel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", " \u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mttest_ind\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mequal_var\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", " Welch's T-test p-value 0.347028503558\n" ] } ], "prompt_number": 32 }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Check the assumption for t-test**\n", "- The big assumption for most t-test is that the data we used came from a normal distribution\n", "- SO, actually the thing we should do before the t-test is to check the normality of the data, which is implemneted as Shapiro-Wilk test" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pdef stats.shapiro\n", "print 'Town1 Shapiro-Wilks p-value', stats.shapiro(town1_hts)[1]\n", "print 'Town2 Shapiro-Wilks p-value', stats.shapiro(town2_hts)[1]" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " \u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshapiro\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0ma\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreta\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", " Town1 Shapiro-Wilks p-value 0.380458295345\n", "Town2 Shapiro-Wilks p-value 0.562481462955\n" ] } ], "prompt_number": 34 }, { "cell_type": "markdown", "metadata": {}, "source": [ "**But what if the data sets are NOT normally distributed**\n", "- simple answer: IT does NOT matter. t-tests are resilient to breaking of the normality assumption\n", "- less simple answer: use the **nonparametric** equivalents that don't make normality assumptions. But be aware, nonparametric methods generally need more data than parametric version.\n", "- the NONparametric version of Welch's t-test is Mann-Whiteney U test" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pdef stats.mannwhitneyu\n", "print 'Mann-Whitney U p-value', stats.mannwhitneyu(town1_hts, town2_hts, )[1]" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " \u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmannwhitneyu\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0muse_continuity\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", " Mann-Whitney U p-value 0.253597522173\n" ] } ], "prompt_number": 35 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }