0018 Pooled Variance t-Test Recipe.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:d1280bf8eb4763d45038d4e9c899ab70037f7c87d4fac517beb3822ec0d89095"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import pandas as pd\n",
      "import scipy\n",
      "import matplotlib.pyplot as plt\n",
      "import seaborn as sns\n",
      "from scipy import stats\n",
      "import statsmodels.api as sm\n",
      "import math\n",
      "from __future__ import division\n",
      "\n",
      "# turn of data table rendering\n",
      "pd.set_option('display.notebook_repr_html', False)\n",
      "\n",
      "sns.set_palette(['#00A99D', '#F5CA0C', '#B6129F', '#76620C', '#095C57'])\n",
      "np.version.full_version"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 1,
       "text": [
        "'1.9.1'"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Data\n",
      "Here we create two samples of different sample size. We want to know if the sample means are significantly different."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sample_a = np.array([4.9, 6.2, 1.0, -4.8, 2.2, 5.4])\n",
      "sample_b = np.array([2.8, 6.9, 8.4])\n",
      "\n",
      "sns.violinplot([sample_a, sample_b], names=['sample a', 'sample b'])\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFVCAYAAADc5IdQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XeUE/UaxvHvJNned1k60hZC7yy9S5GOiKKigoAgWFCw\nSxFpgih2LKCCIiKIimKlKb13GHpZYHvvm2TuH4teC50kk/J+zvHcTTbJPHphn52ZX1E0TUMIIYQQ\n+jHoHUAIIYTwdlLGQgghhM6kjIUQQgidSRkLIYQQOpMyFkIIIXQmZSyEEELozHQjbzKbzT7AfKAi\n4AdMAQ4BnwA2YD8wWlVVmTclhBBCXMWNnhnfCySpqtoW6Aa8A8wGnr/4nAL0sU9EIYQQwrPdaBl/\nBUz422cUAY1UVf394nM/ArfeZDYhhBDCK9zQZWpVVXMAzGZzCMXF/CLw6t9ekg2E3XQ6IYQQwgvc\nUBkDmM3mCsDXwDuqqn5hNptn/u3bIUD6ld5vsVg1k8l4o4cXQggh3JFyqSdvdABXKeAXYJSqqmsu\nPr3LbDa3U1V1HXAbsOpKn5GWlnsjhxZCCCHcVnR0yCWfV25kowiz2fwGMABQ//b048CbgC9wEBh+\npdHUSUlZMtJaCCGEV4mODrnkmfENlbE9SBkLIYTwNpcrY1n0QwghhNCZlLEQQgihMyljIYQQQmdS\nxkIIIYTOpIyFEEIInUkZCyGEEDqTMhZCCCF0JmUshBBC6EzKWAghhNCZlLEQQgihMyljIYQQQmdS\nxkIIIYTOpIyFEEIInUkZCyGEEDqTMhZCCCF0JmUshBBC6EzKWAghhNCZlLEQQgihMyljIYQQQmdS\nxkIIIYTOpIyFEEIInUkZCyGEEDqTMhZCCCF0JmUshBBC6EzKWAghhNCZlLEQQgihMyljIYQQQmdS\nxkIIIYTOpIyFEEIInUkZCyGEEDqTMhZCCCF0JmUshBBC6EzKWAghhNCZlLEQQgihMyljIYQQQmcm\nvQMIIYS4Nqp6CIDKlavi6+urcxphTzdVxmazuRkwQ1XVDmazuSGwAjh68dvvqaq65GYDCiGEgLNn\nzzBmzCgMBhg9+km6d++ldyRhRzdcxmaz+WlgEJB98anGwGuqqr5mj2BCCCH+LzEx/pJfC89wM/eM\njwG3A8rFx42BHmazeZ3ZbP7IbDYH33Q6IYQQAKSkpAAQEqiQnJyocxphbzdcxqqqfg1Y/vbUFmCc\nqqrtgBPAxJvMJoQQ4qKkpOICLhFuIzEhTuc0wt7sOYBruaqqGRe//gZ480ovjogIxGQy2vHwQgjh\nuVJSLhAZphAdaePEhQtER4foHUnYkT3L+Cez2fyYqqrbgE7A9iu9OC0t146HFs6SnJzE6DEPExgY\nxPtvfygjOoVwkuPHVEpGWikVqbFlXwZxccn4+fnpHUtcp8v9EmWPecbaxf8dCbxuNpvXAC2AKXb4\nbOFiTp06SWZKCvFnz5CamqJ3HCG8gs1m48zZc5QvCeVLgabBmTOn9I4l7OimzoxVVT0FtLz49R6g\ntR0yCReWnZ3919c5OTk6JhHCe1y4cJ78/CLKl9KoUKr4/OfYsaNUq2bWOZmwF1mBS1yXnJz/l3F2\ndpaOSYTwHn8u9lGlvEbJSAgKUDh8+IDOqYQ9SRmL65KZmfG3rzN1TCKE99i7dyeB/grlokFRIKaC\nlX17t+kdS9iRlLG4LhkZ6XBxFHxGRrrOaYTwfJqmsXPHJmpUtmK4+BO7VhWNC/EpJCTI4h+eQspY\nXJek1BRsAQEApKWl6pxGCM938uRxkpIzqFdN++u5etWLv968eYNesYSdSRmL65KQlEhRgD9agD/J\nKcl6xxHC461btxpFgYbm/5dx6SgoVxLWrvlJx2TCnqSMxXVJSkrEGhiAJTCA8/EX9I4jhEezWq38\n9uv31I2xEfqvBYZb1rdyWD1GXNwZfcIJu5IyFtfMarWSlZaGNTgIS1CglLEQDrZ58wZS07Jo21j7\nz/da1dcwGmHFiuU6JBP2JmUsrllSUiKazYY1OAhrcBDpKcnYbDa9YwnhkTRNY8mXnxId8c9L1H8K\nC4FmdWz8/NP3/5jlINyTlLG4ZmfOnAbAEhaCJTQEm8UiozmFcJCtWzdz5OgJbmv1/1HU/9a9tY3C\nIgtLlnzu3HDC7qSMxTX7s4yLwkIoCg/9x3NCCPspKipi3kdvUipKoU2j/54V/6lcSWhZ38Z3337N\nhQvnnZhQ2JuUsbhmx44fRQsMQPPzwxJWXMYnThzXOZUQnmfZsi85GxfPwK4Wrra5Xf9ONoxGG2+/\nNRNNu3xxC9cmZSyu2f7DBymIDAdA8/XBGhrCAVmSTwi7OnnyOIs+/4TGNW00uMS94n+LCIX+HS3s\n3LWHn3763vEBhUNIGYtrkpOTTUpCPIVREX89VxAVwSH1kPw2LoSd5OfnMW3qiwT627i/57UPjuwY\nq1GrisbcuW9y+vQpxwUUDiNlLK7Jvn17QdMoLFnir+cKS5YgNyODc+fidEwmhGew2WzMmjmZc+fj\nGX675T/ziq/EYIDht1vx97EyaeI4GV3thqSMxTXZvXsnGI0URkf99VxB6ZIA7NmzS69YQniMjz9+\nn42bNnNXFyu1q17/1abwEHhkYBHJySlMfulZCgoKHJBSOIqUsbgqTdP4Y/MGCkpGgfH/o0msIUHY\ngoPYIOvjCnFTFi9eyNKlS+jQ1EaXFjd+26faLTC0n5WDhw4zdcrzFBUV2TGlcCQpY3FVZ86cJjUh\nnrwK5f75DUUht3wZ9uzeSV5enj7hhHBjmqbxxRcL+PTT+bSsb2NQdxuKcnOf2byuxv09rWzbvpMp\nU56XM2Q3IWUsrmr9+t8ByK9Q9j/fy7+lHDaLha1bNzk7lhBuzWaz8cH7b7Fgwce0rG/jwT62yy7u\ncb3aN7lYyNu288LzY8jJybbPBwuHkTIWV2Sz2fj+5+8pLBWNLTDgP98vjC6BLSiQH376QYd0Qrin\n3NxcXp78HN98u5xbm9kY2tf29ztAdtGhqcaI/lYOHz7MmMeHyUBLFydlLK5o3749pCclkRNT+dIv\nMCjkVKnIvj27SEpKdG44IdzQ+fPnePKJh9iydSv33Gblntvsd0b8b83qaoy730paaiJjHh/O9u1b\nHXMgcdOkjMUVffX1EjQ/P/JuKXfZ1+TGVALg22+/dlIqIdzT2rWreGT0gyQnnefJQVY6N9du+h7x\n1dSorDHhoSLCg/MZP/4Z5s9/H4vF4tiDiusmZSwu68yZ0+zYuplscxWutCafNTiIvFvK8d0P38q9\nKSEuITs7m1dnTeGVV6ZQLrqQl0YWUSfGeYvllIyEF4dZaNfYxldfLWbckyNlH2QXI2UsLuvzLxaC\nyUiOOeaqr82qY6YoP1/OjoX4l+3btzJyxL2sWbuKXm1tPDPEQlS483P4+cLg3jZGDrASF3eCUaMe\nZNmyJVitVueHEf+h6LWUYVJSlqyh6MKOHTvCo4+OIKt2dbIa1bum90Ss3UhoYioLPllEeHjE1d8g\nhAdLTU3lgw/eZN26dZSNVhjar4gql7/b41TpWfDpCiO7VYXq1Srz6GPPEhNTXe9YXiE6OuSSNybk\nzFj8h6ZpvDP3bTQ/P7Lr1Lzm92U2rEtRYQGfLJjvwHRCuDar1cp3333N8OH3sGH9Onq3szFphOsU\nMRSv1vXY3VYe6m8l/sIpHntsJO+9O4esrCy9o3ktOTMW/7Fq1S+8+up00mMbkmuuel3vDd26m+Aj\nx3lt9tvUrFnLQQmFcE3btm3hg/fnEHcunlpVNO7rYaV0iau/T0+5ebBstYE12wwEBwUw6L5hdO/e\nG5PJpHc0j3S5M2MpY/EPaWmpDBl+P9mB/iR1bc/1DvVUCoso9f2vlA6L4MP35uPr6+uYoEK4kGPH\njjJ/3jvs2r2HklEKd3W20LCG40dK29OZeFj8k5FDJxXKlyvF4CGjadmyNYo7/Uu4ASljcVWapvHi\nhGfZuWsHiT06YQkLvaHP8TsXT9Tq9fS7fQAPDR9l55RCuI64uLMs+PR9/li/gaAAhV5tLXSK1XDX\nk0pNg92qwpJfjMSnQLWYygx5cDQNGzbWO5rHkDIWV/XVV18wf/4HZDSpT07Najf1WWFbdhJ05AQT\nJ06hefNWdkoohGuIizvL4i8+Zc3a1fiYoEtzK91a2rjEInVuyWqFDXsUvl1rIjVDo07tGgy67yHq\n1WsgZ8o3ScpYXNH+/Xt5+pknyC1fhrS2za/78vR/WK1E/7yWoNx83nv7Q8qWdaHRK0LcoLNnz7Bo\n0cesW7cOH5NG+8Y2urexEXYdew+7k6IiWLNd4ccNJtKzNGrVrM69g4bTsGFjKeUbJGUsLisu7iyP\njnmYHKOBxNs6ovn62OVzjVk5lPxxFSUjo3h7znuEhobZ5XOFcLbDhw+y5MuFbN6yGR8TdIq10bWl\n55bwvxUVwe87FVauN5GaqRFTtSJ33jWEli1bY7T3otoeTspYXFJ6ehqjHhtJalYmid3aYw2x708X\n38RkSvz2B1VjqvPazDkyoEu4DZvNxo4d2/hqyafs23+IoACFjk2t3NrcRmiQ3un0UWSBjbsVftxo\nIiFFo2zpEvQfcD+dOnXBz89P73huQcpY/EdWVhZjnxnD2bOnSbq1LUXRUQ45jv+ps0T+sYXGTZsz\ncfxkfHzsc+YthCMUFhayevWvLFu6kLhzCUSEKnRtYaFtY40A6RsAbDbYcUhh5Xojp85DaEggPXvd\nQa9efWXBn6uQMhb/8FcRnz5FSvsWFJQr49DjBarHCN+6m8axzZn4ohSycD2pqSn88MO3fL9iGZlZ\nudxSGrq2tBJb231HRzuapoF6SuHnTQq7VQM+JiPtO3SiX787qVz5+tYo8BZSxuIvmZkZjHv2yeIi\nbteCgvKOLeI//b2QJ7zwklyyFi5BVQ/zzTdf8scfv2Oz2ahXzUbXlho1KrnXPGG9XUiGXzYZ2LjH\nQGER1K1Tk7797qFZsxZyX/lvpIwFAPHxF3j6+bEkJSaS0ra504r4T38WcvWatZny0nRCQkKcenwh\noPhS9Pr16/ju2y9RjxzH30+hTUMrnWJtlHLM3RqvkZ0Lf+xUWLXVREqGRnSJcHr2upNu3brLIE4c\nVMZms7kZMENV1Q5mszkG+ASwAfuB0aqqXvbDpYyd79ixIzzzwlPk5OeT0qElhSX1WafP/+RZIjdu\no1TZssyaNpvo6JK65BDeJykpiZUrv+PHlV+TkZlL6SiFjrEWWjfQCPDXO51nsVphl6qwaouBw6cU\nfHyMtG/Xkd597vDqTSnsXsZms/lpYBCQrapqS7PZ/B3wqqqqv5vN5veAn1VV/eZy75cydq7169fx\nyqxpFPqYSO7YGkv4ja2uZS++FxIp8fsmggICeXniVGrWrK1rHuG5NE1jz55drPjuKzZv2YKmadSv\nbqNTrEatKhoG2S7H4eISYPVWAxsuXsI2V69C7z4Dad26ndfdrnJEGd8O7AUWqqrawmw2x6mqWv7i\n93oDXVRVfeRy75cydg6r1cqnC+bz1ZJFFJWIJKVdC2wuskyQKS2DEus2YczN45FRj9O9ey+9IwkP\nkpOTw6pVv7Diuy+JO5dAcGDxpegOTW1Ey4BfXeTmwfrdCmu2FS+3GRoSQNdufejRow+lSpXWO55T\nOOoydSXgi4tlfE5V1XIXn+8IDFFV9b7LvVfK2PEyMzOY+srL7N25g5yYymTENgAXG0ihFBQSuX4L\nfucT6Nz1NkY//LjMVxQ35eTJE3z//desXvUz+QUWKpeDjrFWmtXWkEH8rsFmg0MnFVZvVdilGgCF\n2KaN6dV7AA0bNsHgwZcrLlfG9hywb/vb1yFA+pVeHBERiMnkWsXgSbZv387Tzz9PZkY66c0aklut\nys0vcekAmp8vKR1aE7J7P7/+/CMHDx/ktZkzqVpVpkWIa2exWFizZg2Lv1jInr0H8TFBszo2OsTa\nXGofYVHMYIDaVTVqV9VISbexdruBdTt3sGXrdsqVjebOuwbRu3dvrxrgac8z4++A2aqqrjObzXOB\nVaqqfnW598qZsWNYrVYWfvYJX375OdaQYFLaxGKJdI9rcn7n4onctB2TxcqokY9y2209Zf1bcUWp\nqan8+OMKVv6wjNS0LKIjFDo0tdCmoUZwoN7pxPUoshQvJLJqi5FjZ8HP10THjp3p1bu/R81ZduRl\n6kUXB3BVAz4EfIGDwHAZTe1cp0+fYvqsaZw+fpScqpXIbNoAzce9Visw5OUTsWErfhcSadAklnFj\nniIqysV3ZxdO9/e5wVarjboxGp2aFf+vB1/h9BqnzhcP+Nq8z0CRxbPmLMs8Yw9mtVpZuuxLFiz8\nGKvRSFpsA/IrVdA71o3TNIIOHSVs9wF8fX15bPQYOnbsLGfJXs5qtbJ+/Tq+Wb6Iw2rx3ODWDYrn\nBpeW39c8UnZu8QYVq7YWb+VYMjqC3n3uolu3HgQFuecuHVLGHurUqZO8Mns6p44dJe+WcmTENsTm\nIRMmjZlZRG7cjk9SCo2aNmPMo2OJjo7WO5ZwspycHH766Qe+/WYRSckZlIxUuLWZzA32Jn/OWf51\ns4EjpxUC/H3o2q03ffve4XajsKWMPUx+fj6fL1rAsmVfYvPxIa1pA/IrlXfJQVo3xaYRdPgoYXsO\nYDKaGDp4OL169XX7S1Xi6lJTU/j66yX8uPIbcvMKMVfU6NrSRv3qcinam506Dz9vMrBtvwENhdat\nWnPnXfdTtWqM3tGuiZSxB9m+fSuz33iV9OQkcqtWIrNRXWz+nj0dyJiVTfjWXfidT6BC5So89cTT\nVKtm1juWcIDz58+xdOkifv31J6xWG7G1i/cOriyjosXfpGTAqi0G1m43kleg0bhxQwYOHEydOvX0\njnZFUsYeICEhnnfmvsW2zRuxhoWS1qwhhaW86LKtpuF/Oo6I7XtR8vPp1q0HQwYPk/VuPcT58+f4\n7LN5rF27FqNBo3VDG7e1slEyUu9kwpXl5sHqbQq/bDaRlaNRu5aZ+x8YSb16DfSOdklSxm4sPz+f\nJUsWseSrL7ACmXXMZNeq7nILeDiLUlhIyJ6DBKvH8QsMZNjgYdx2Wy+5dO2mEhLiWfT5x/y26leM\nBo2OTYvPhCP0XbFVuJmCwuLBXivXm0jP0qhfrw73PzCCWrXq6B3tH6SM3ZCmafzxx1reef8dMlNT\nyK1UofiSdJBMoITi5TTDt+3GNyGJsrdU5PHRY1z2t2HxX9nZ2SxevIBvv12Ggo32TWx0b20j3HvW\neRAOUFgEa7cr/PCHicwcjRbNmzF02COUK1de72iAlLHbOXz4IG+++yYnj6pYIsJJb1rfuy5JXytN\nw//MOcJ37MWQk0uT2BaMfGiUy/zFE/9ltVr54Yfv+Gzhh2Tn5NGqgY1+HWxEyt0GYUcFhfDLpuIz\n5SKrQq+efbh30IMEB+s7JUrK2E0kJMTz4fz32fD7WrQAfzLq1ya3aiUweNgoaXuzWAk+dJTQA4dR\nrDZ69uzLfYMeICRErnW6kmPHjvLGnGkcO36KGpU1Bna1UtG5W2oLL5ORBcvXGPh9p4HwsGAeHjWW\n1q3b6bZugZSxi8vKymLxl5/xzbdfY9M0smpVI7u2GU1Wtr8uhrx8QvYcIOjYKXz9/Rl093306dPf\n67ZpczVFRUUsXDCPZV8vITgQ7ulmIbaO5nEz8YTrOnkOPllh4swFiI1twpgxzxER4fzRgVLGLqqw\nsJAVK5azcNECCnJzya1SkcwGteW+8E0ypWUQunMf/ufjCY2MYviQ4XTocKsM8tLB+fPnmDF9PEeP\nnaRNIxt3dbER5Bq7eAovY7XCb1sUlq0yEhQUxFNPT6JRoyZOzSBl7GKsVitr167ig/kfkJmaQkHZ\nUmQ0qoslIlzvaB7F90IiYbv24ZOSRtlbKjJy2MM0aRIrS2s6yfbtW5k2dTwKRQzpY6FJLa/+ay9c\nRFwCvPeVD+eTNO6//0EGDhzktJ8JUsYuQtM0tm7dzNx5c4k/ewZLZATpjepQWKaU3tE8l6bhfyqO\n8D0HMGRlU71WbUYOG0XNmrX0TubRfvzxe95++zXKl9R47G4LUfJ7pnAhBYXwyQoDm/ca6NqlG488\nOhaTyfEb60gZu4ADB/Yx96P3OHb4ELbQYNLr1ya/ogcuYemqrDYCj54gbP9hlLx8mjZrybAHH+KW\nWyrqnczj/Pjj97z55mzqxmg8fKeVAM9eIE64KU2Dr1cb+P53Ax06dGDcuBcxOHitVSljHZ08eYIP\n57/Pru1b0QICyKhXg9yYysgCu/pQiiwEHTpC6MGjYLHQsVMXBt//INHRJfWO5hE2b97A5MnjqRNj\n47G7rZjkNr1wcd+tU1i+2sgd/e9k6LCHHXosKWMdJCTE88mn81i7dhWaj4msWtXJqVHN7fYY9lSG\n/AKC9x0i+MgJDAYDvXv14+6B98rymjchNTWFh4bfS3R4Ac8OseAng9iFG9A0WPiDgTXbDEyZ8gqN\nG8c67FiXK2M5NXOA9PR03n3vLR4cOoi1v68hu2Y14vt2I7tuTSliF2Lz9yOzaQMS+nQlu2I5vvlm\nKYMeGMgXXywkPz9P73hu6Z13ZpObV8CI/v8v4hkf//PUWB7LY1d7rChwd1cbZaNhzuvTyM/Px9mk\njO0oPz+PRYsWcP/gu1mxYjlZlSsQ36crmY3rofnJTTNXZQ0OIr1lUxJ73kpmiQgWLJjPoAfuZuXK\n77BYLHrHcxtxcWfZuHETkaEapUvonUaI6+PjA/d2t5KcksG6daucfny5TG0HFouFn39eyfwF88nN\nzCCvQlmyGtbBEiarP7kj38RkQnfuwzcphRKlyzBi6AhatWor06GuYv789/l62WJmP2khTNaXFm5I\n02DCeyaCwqow5415DjnG5S5TyzXTm6BpGps3b+DdD94lOf4ChSVLkNGiPUUl5bTAnRWWLEFy1/b4\nxV3Auns/U6dOonI1M4+MfMTldoBxJYcP7aZiWaSIhdtSFKgbY+PXLaewWCxOmer0JynjG6Sqh3nn\n/bc4eugg1rBQ0tu3pKB8GZmm5CkUhYIKZUkoV4bA46c4sfcgY8c+SrOWrXlo6EjKlpWd7v/t3Lk4\nala06R1DiJtSJlrDYrGRmJjg1L/nUsbXKSkpkQ8+msv639eg+fuR0ayhTFPyZAaF3GqVyatUgaBD\nR9iydTNbt2yiV88+3DdoiO47wLgSk9GITne9hLAb7eLvk848KwYp42tWUFDAV0sXs/jLz7HabGTV\nqVG8kYOvbOTgDTQfE9n1apEbU5mQPQf57tuv+WXVLwwfMpyuXXvImtdARGQkyenpescQ4qakZCgY\nDAqhoc4d8yMDuK5C0zTWr1/H23PfJjM1hbyK5clsVBdrcJDe0YSOTKnphG/fjW9CMmVvqcjjo8dQ\nr14DvWPp6sMP32X58q+IqaD950LRs0Osl3zPv6eZyOvl9Xq/fto8I5pPRd56+9NLvvZmyTzjG3Du\nXBzjnn2CadNeIlWzktylHWltm0sRCyyR4SR3bkdq2+bEpSTxzDNPMP2VKaSnp+kdTTctWrRB0yAz\nR+8kQtyYC0lw9IxCi5adnH5sOTO+hMLCQpYsWcQXX36O1aCQ0aA2udWqgkEGZ4n/UiwWgvcfJuTA\nEXz9/Bgx7GG6devh8DVuXY2maYwedT95WXG8PMqCk2+5CXHT5n1jYPM+EwsXLiU8PMIhx5Az42u0\nf/9eho4YzOeff0p2+TIk9OpCrjlGilhclmYykdWgDok9byU7NJi33nqNR8aM4syZ03pHcypFUXhg\n8CjiU4r3jBXCnZyIg/W7DPTp3d9hRXwlcmZ8UX5+Ph9/8iHffbccW3AQabENKSgr2xqK66RpBJw8\nQ/j2vRitVoY8MJR+/QZ41QCvCePHsmfPLiaOKKJstN5phLi6gkKY/IGJvKIQPvzoC4KCHHcrUjaK\nuIJDhw4wbeZUkuMvkGOuSmbDurKGtLgphrx8wrbsJODseapUN/PcUy9QvnwFvWM5RUpKMqMefoCw\nwFxeGCabRQjX98l3BtbtMDBlykwaN27q0GPJZepLsNlsfLH4M54c+yiJ2Zkk39qWjNiGUsTiptkC\n/Elr14K01rEcP32Sh0cPY9WqX/SO5RRRUSUYO248cYnF9+Bssg6IcGG/bVFYt8PAgAEDHV7EV+K1\nZ8bp6elMnTGZ/Xt2kVupAhnNGsmcYeEQhtw8ItdvwTchmY63duXR0WPw9/fXO5bDLV26mHnz3qd7\naxsDOksjC9ezW1V48wsjsU2bMn7CdKfcTpLL1H+jqod4cdLzZGdlkd6kPrnVKssylsKxbDZC9hwk\nZP9hylS4hWmTZ1C6dBm9UzmUpmm8/darrPxxJXd1sdKtlcvcmRIC9RTMXmiiYsVKzJz1LgEBAU45\nrlymvmjDhj8Y+9TjZFosJHZrT271KlLEwvEMBrIa1iGlU2vOJ8Yz+rERqOohvVM5lKIojBr9JG1a\nt+bLX4ys3ip/z4RrOH4W5izyoVSpUkyZOsdpRXwlXlXGy5d/xZSpE8kPCyGxW3sskc4fvi68W0HZ\n0iR27UC2pjH2qcfZsOEPvSM5lNFo5KmnJxAb24SFPxhZs00KWejrRBzMXuhDeHgU06a/RVhYmN6R\nAC8q488XLeCDD94lr0JZkju3xRbg+ffshGuyXvxlMD8shClTJ7J27Wq9IzmUj48PL7wwlaZNGrHg\neyOrZA6y0MnRM/DqAh9CwyOZOesdoqNdZ+6dV5Tx0mVf8tnCj8mtcgtpbZqjydJAQme2AH+SO7el\nsGQUM2dN9fgzZF9fX14cP53mzWL5bKWRnzdKIQvnOnwKXlvoQ0RkCWbNeo/o6JJ6R/oHjy/jH39c\nwbyP5pJXsTzpLZrISlrCZWgmEykdWlEYFcHU6S+xa9cOvSM5lK+vLy+8OJVWLVuy+Gcj361VZMtF\n4RT7jiq8ttBEiehSzHr1PZc6I/6T3cvYbDbvNJvNay7+M8/en389VPUQb775GvllS5HWOhYMBqJ+\nWfuP18hjeazn48g1G0ju2ApLSDCTp00iKSkJT2YymXju+cl07NCR5WuMfPWrQQpZONSOQ8XTl8qX\nr8CsV99ieZ94AAAgAElEQVQjMjJK70iXZNfrtWaz2R9AVdUO9vzcG5GVlcnElyegGQyktW7Gf/Z0\nE8JFaL6+JLdtjvHHVUx6eTxvvPa20zc2dyaj0cjYcS/g5+fPjz+tpKAI7r3NJn9Fhd1t2qvw0XIj\n1apW5eWprxMSEqJ3pMuy6zxjs9ncDPgUOE1x0T+vquqWS73W0fOMX5k1jbVrV5HUtT1FJSIdeSgh\n7ML/1Fki/9jC4MHDuOuue/WO43CapvHRh+/y9fKltGloY3BvKWRhP7/vUPhkhZE6tWsx6aVZBAYG\n6h0JcNKiH2azuQ7QTFXVeWazuRrwI1BdVdX/LL9jsVg1k8kxq50cPnyYQffdR1atamQ1queQYwjh\nCBFrNhCanMbKFSuIiPD8qXeapvH++3P56KN5NKtrY3g/G160p4ZwkFVbFD5baaR5sya8OnuOq614\n55RFP44AnwOoqnoUSAEuucxQWlouSUlZDvln+qxX0Xx9ya5Tw87/ekI4VlajuhQVFDDnjXcc9vfD\nlf5JTs6mf/9BDB48jC37DIx51YjF+v//HjM+/mczy2N5fLXHP2/8s4hjee75qWRlFen+5/zv/1yO\nvct4CDAbwGw2lwVCgQt2PsYVJSTEc2DPbrJqxqD5ynYxwr1YwkLJrVyBX377icLCQr3jOM1dd93L\n8OEPk52r8N4SAxaL3omEO0rJgMU/G2nVsiUvvDgVXzfqAHtfpjYBHwMVLz71tKqqmy/1WkfdM16y\nZBEff/whCf1uwxrsuD0phXAUv/PxRK1az4QJL9OiRWu94zjVt98uY+7ct2lUw8bDA2x48Dg2YWc/\nrldY8quRNm3a8MwzE112D/HL3TO26x91VVUtwH32/MzrtW7jHxRFRUoRC7dVULokmp8vGzdv9Loy\n7tOnP4qi8N57bzF3KYwcYMNBQ0uEB/l5o3sU8ZV43NjF8+fiKIxwjbVGhbghBgNFYaGcPHNK7yS6\n6N37dkaMeIQdhwx8+LUBq/Xq7xHe67ctCot/NtKmdWu3LWKw85mx3nJzc8nPzsYaUknvKELcFEtI\nEBcunNc7hm769u2PxVLEvHnvYzLC0L4y7Un817odCp+vNNK8eTOeduMiBg8rY4ulCADNjf8PEQJA\nM5qwFBXpHUNXd9wxkIKCAj777BP8fOG+HjbZ7VT8ZdNehU9XGGncqCHPP/+y2y+U497p/yUoKBgA\npdC7f4gJ96cUFhIQJOMe7rnnfvLzclm6bAn+vjCgsxSygJ2Hi1fWqlO7FuMnTMfHx0fvSDfNo8rY\naDTiGxiIMTdP7yhC3BRTXh7hYZ6/6MfVKIrCg0NHkpuXw8qVPxDgp9GrnSxm7c0OHFd4b4mRmKpV\nmPTSLPz8/PSOZBcedxemZs3a+Ccl6x1DiBtnteKbnEajug30TuISFEVh9Ogn6dC+A1+vNvKb7Ifs\ntY6dhbcWmyhXtixTpr7uMktc2oPHlXFsoyYYM7Iw5OTqHUWIG+KbmAxWKw0aNNQ7isswGAyMHfcC\nzZvF8vlKI+t3SyF7mzPx8PpnPkRGRjFtxluEhITqHcmu7Lrox/Vw1KIf58+fY+iw+7AE+l9yrnFK\nl/aXfN+/t7aT18vr9Xp9xO+bCUtMYfHnywgICLjk67xVYWEhEyeMY8/efYy600qTWnLJ2htcSIYZ\n833w8Qtl9mtzKVWqtN6RbtjlFv3wuDPjsmXL0bBxU4x5+chGqcLdGHNyCThzjh639ZIivgRfX18m\nTnqFGuZqzF1qZM8ROUP2dImpMOtTHxRjIDNeedOti/hKPO7MGGDPnl08++yTZDSqS05ts6MOI4Td\nhW3aQfCJ03zy8SJKliyldxyXlZ2dzbPPPMKZM6d5/B4rtavKL96eKCUdZnzsQ77Fj5kz36Zy5ap6\nR7ppXnNmDFC/fkMaNoklbO8huXcs3IZPUgpBx07St29/KeKrCA4OZtr0NyhbthxvLDJy8IScIXua\nlAx45RMfcgt9mTp1jkcU8ZV4ZBkDPDZ6DEZFIWzrbrlcLVyf1UrEll2ERERw36AheqdxC6GhYcx4\n5W3KlinLG4uMHDguhewpktNh5ic+5OT7MG3a61Sv7vlXOD22jEuXLsPgB4YSEHeeoMPH9I4jxBWF\nbd+LKS2dJx8bJ/eKr0N4eDgzZr5NmTLFZ8hyD9n9JaYWX5rOzvdl6rTXMZtr6h3JKTy2jAFu7zeA\nRk2bEbZzHz5JKXrHEeKS/E+eJejIcfrdPoDmzVvqHcfthIdHMHPWO1SsWJG3FhvZul8K2V3FJcD0\n+T4UWvx55ZU3qVGjlt6RnMajy9hgMPDsUy8QHhVFiXWbMWbl6B1JiH/wSUohcvN2qppr8OCQh/SO\n47b+vGRdw2xm7lIja7ZJIbub42eLz4gVYzAzZ71DTEx1vSM5lUeXMUBISAgzXp6Jv6JQYvV6DPkF\nekcSAgBTRiYl1mwkKqoEUyZNc/uF7vUWFBTMlKlzaNK4EQu+N/L1KoMMF3ETOw8rzPzUREhoJK+9\n/j6VKlXWO5LTeXwZA1SsWImpk2fgm5tH1OoNspGE0J0xO4cSqzcQ5OvLrOmvER4u61Dbg7+/PxMm\nvkKXzl1Y8buBj5YbKLLonUpcyW9bFN5ebKRixcq8PudDSpcuo3ckXXjkPOPL2bx5I5OnTKAwPIzk\nW1uj+fo6O4IQGLNyiP7td/ytNmbPnON1l+OcQdM0Fi1awGeffUK1W+CRgRZCZRMsl2K1wqIfDaze\nZqBZbFOefe4l/P09f/Di5eYZe1UZQ3EhvzxlAgXhoSR3aoPmJ4UsnMeYlU30r78TgMKrM16TInaw\ndetWM3v2NMKCbDw6sIhbvPOky+Vk5cDcpcXzw/v3v5MhQx7C6CX70EsZ/83WrZt56eXxFAUHktSx\nNbYgz9n5Q7guU2oa0as3EGAwMnvmHKpUidE7kldQ1cO8PPkZsrIyGdzbSot6ciNZT6cvwFuLfcjM\nNvDIo2Pp0uU2vSM5lZTxv+zZs4vxk56nwGAguWMrLBFhesYRHs7vfAJRv28iLCSMmdNnc8stFfWO\n5FXS0lKZOuV5DhxU6RhrY2BXGz4yXs6pNA1+36nw+UojoaFhjJ8wA7O5ht6xnE7K+BJOnjzO08+N\nJSsvl5S2zSksI0sQCvsLOH6KiM07KFO+ArOmzSYqqoTekbySxWJh/vy5LF++jMrl4OEBFqJl3JxT\nFBTCwh8MbNhtoEH9ujzz7EteO2hRyvgyEhLieebFp0g4d4702AbkVvfs9U+FE9k0QnbvI+TAEWrW\nrc/LE6cQFBSsdyqvt2HDH7z22lQ0ayGDe1uIreMSP4o81pkL8N5SEwkpcPfd93HPPQ94zf3hS5Ey\nvoKcnBxenjaJPTu3k22OIbNJPTB4xawv4SBKkYWI9Vvwj7tAt9t6MnrU4zKP2IXEx19gxrTxqEeP\n06aRjXu62fD30zuVZ7HZYNVWha9+NRIcHMLTz0yiQYNGesfSnZTxVVitVj78aC7ffrOUwtLRpLZp\njk3+doobYMzMosS6TZgysxk54hF69+6ndyRxCRaLhYUL5/HVV4uJjlAY0b+IKuX1TuUZ0rNg3jdG\n9h9TaNq0MU8++SLh4eF6x3IJUsbX6LfffmbOG69S5OdLSrvmFEVF6h1JuBG/uPNEbdiGv68fE194\nSc4E3MDevbuZNXMSqWkZ9Gpro2dbGybvvYp607YdUFjwvYlCi5Hhw0fTo0cfFEWWJ/2TlPF1OHpU\nZfxLL5CRnk560wbkVvO+pdnEdbJphOw7SMjeQ5SvVJkpk6ZRqlRpvVOJa5Sdnc2777zGmrVrqFwO\nhvezUCZa71TuJTcPPltpYNNeAzFVK/H0My9RocItesdyOZcrY7kxegnVqpmZ+/ZH1KpVh/DNOyj1\n1QqwWP/6ftQva//xenns5Y9/Wk3U6vWE7D1E+46deXvOe1LEbiY4OJinn5nA889PIjkjgIlzTfy8\nUcFm0zuZe9h3VOHFd33Yst/Ivffez+tzPpQivk5SxpcRHh7OzOmzuWvgIIz5BZT8aQ3GzGy9YwkX\n45OUgm9KOgFJKTz22JM8Pe45/PxkrIG7atOmHXPf/4xGjZqy+GcjMz4uHgUsLi0vHz75zsBrnxkJ\nDi3F66+/x6BBQ2Sw4g2Qy9TXYNu2LUx95WUKigpJbd6Y/IoyysPraRpBh44StmsfEVElmDxhiixt\n6UE0TWPVql94773XsRQVMOBWKx1jNZlk8TcHjit8/J2J1Azof/sA7rt/KL6y3v9VyT3jm5SQEM+k\nKRM4dexo8fSnxnXBi+fKeTOlsJDwjdsJOHuexrHNeWbc84SEhOgdSzhAcnISc+ZMZ8eOXZgraTzY\nx0pJLx/TmVcAS34xsHa7gXJlSzJ23ERq1qyldyy3IWVsB0VFRXw0by7fffs1RVGRpLZphjVEtoLx\nJj4paUT9vgVjbi7Dho6gX78BMlLUw2maxq+//sT7c+dgsRQyoLOVjk298yz54AmF+d+aSM3QuL1f\n8dmw3Ja5PlLGdrRhwx/MnD2dAquVtJZNyK9QVu9IwtE0jUD1OOE79hIaHs6kFydTs2ZtvVMJJ0pK\nSuKNOdPYsXM3NSprDO1jpYSXrOiYXwBLfjWwZpuBsmWiGTtuArVq1dE7lluSMrazCxfOM/Hl8Zw9\neYLsmtXIbFRXVu3yUEphEeGbdxBwOo4GTWJ5/ukXCAkJ1TuW0IGmafzyy0ren/smmq2IO7tYaN9E\nw5MvjqinYN43PiSna/TpczuDBz8kZ8M3QcrYAQoLC3n/w3dZ+f23FEVHkdKmmWzH6GFMqemU+GMz\nxuxcBj8wjDvuuAuD/NLl9RITE3j9tans3rOPujEaQ/pYifCw38+KimDZagO/bDJQqmQUY8dNoE6d\nenrHcntSxg60bt1qXn19JkVAauumFJSVOaaeIODYSSK27iY4JISJL7wkP4jEP9hsNr7//lvmzXsX\nH6OV+3t6zqYTZy7AB1+bOJcI3bv3YNiw0QQEBOgdyyM4pYzNZrMBeBeoBxQAw1RVPX6p13pSGQPE\nxZ1h/EsvEh93lsx6tciuVxOPvnblySxWwrbuIuj4KWrVrc/45yd67XZv4uri4s4ya+Ykjhw9Qcv6\nNgZ1txHgr3eqG2Ozwc+bFJatMhIaEsKTY1+kSZNYvWN5FGeV8e1AT1VVHzSbzc2A51RV7Xup13pa\nGQPk5+fzxpuzWbvmN/LLlSatVSyan8y7cyfG7Byift+MKSWNgQMHMWjQYK/e7k1cG4vFwhdfLOCL\nLz4jKgxGDiiiqpstR5CWCR9+beTQSYUWzZvx+JjnCAsL0zuWx3FWGc8GtqiquuTi4zhVVS/5R9IT\nyxiKB3j88MO3vDv3bawBASS3a44lUnYrcQd+5+OJXL8NP4OB558ZT7NmLfSOJNzMwYP7mTF9Aqmp\nadzeyUq3lu4xBWrvUYWPlpsoKDIycuTjdOvWQ6bsOYizyvhDYJmqqj9dfHwaqKyq6n9WeLVYrJrJ\ng7dG2bdvH4+PHUtmZgapLZqQX6mC3pHE5WgawQePELpzH+UrV+Lt1+dQvrybndYIl5GVlcWUlyex\navU66lXTGH67lWAXHddptcLyNQZ++MNA1SoVmPHKa1SuLBvjOJjTzow3q6r61cXHZ1VVvWQLeeqZ\n8d+lpqYy/qUXOHHkMFl1apDVoLbcR3Y1Fivhm7YTeOosLVq15elxz+LvLwNVxM358wrZ+++/RViw\nxqg7i6hSTu9U/5SRDXO/MnL4lELXrt0YNeoJWc7SCZy1a9MGoDuA2WxuDuy18+e7lcjISF6f9Qa3\ndrmNkP2HiVy7EaXIoncscZEhN4/oX9YSeDqOBx4YxvgXJkkRC7tQFIWePfvy6qvvYDCFMWO+iY17\nXOcX8VPnYfIHPhw/Z+LJJ59hzJhnpIh1Zu8zY4X/j6YGGKKq6pFLvdYbzoz/pGkaK1Ys5725b2OJ\nCCO5QytsgfJDX0+mtHRKrNmIb5GFF56bSPPmLfWOJDxUeno606a+wL79B+nW0saAzjZd7yNv2a8w\nb7mRsLAwJk6aKRucOJnMM3YB27ZtYfLUiRQajSR3aCkDu3Tidz6eqN+3EBwUxIwps6haNUbvSMLD\nWSwW5s59gx9++J7GNW0Mv92GsydaaBqsXK+w9DcjtWpWZ/yEGTJlTwdSxi7ixIljPPviM2TlZJHc\nviWFpaL1juRVAk6cJmLTDspVuIXpL88kOlr++wvn0DSNb75ZyocfvkuV8vDEvRaCnHSBzGaDz1ca\nWL3NQNu2bRk79gW5LK0TKWMXkpSUyLhnnyQxIZ7UNs1kowknCTp0lLDte6hRpx5TJk0jKEh23BLO\nt379Ol55ZTKlo2yMu89CmIN337RYYd5yA5v3Gejf/04efHCELOmqI2cN4BLXIDq6JG+9/g4VK1Uh\nct0mAk6c1juSZ9M0gvccJGz7HmKbt+SVqbOkiIVuWrdux+TJM0lO92HGJyYyshx3LIsV3l9aXMSD\nBw9j2LCHpYhdlPy/opPQ0DBen/UGNWvXJWLDNilkR9E0QvYcJHTvQTp06syEFyfL5Tmhu4YNG/Py\nlFdJyzQx81MTmdn2P4bNBh8sM7D9oIHhwx/mrrvutf9BhN1IGesoICCAaS+/Qq269YsL+bgUsr2F\n7DlAyL5DdOrcjXFPPitLWwqXUadOPSa/PIvkDBOvf24iv8B+n61p8NkPBrYdMPDggw9x++132u/D\nhUNIGevM39+fqZNnULteAyI2bcf/zDm9I3mMoAMqIfsOc2uX23hyzFNyeU64nHr1GvDccy9x+oLC\nu18ZsVrt87k//KGwZruBAXfcxYABd9vnQ4VDyU8nF+Dv78+Ul6ZROaYakeu34hufqHcktxdw7BRh\nO/fRonU7nnh8nBSxcFnNm7fkkUeeYN9RhaW/3fyf052Hi3ddateuHUMeHGGHhMIZ5CeUi/D3D+CV\nqbMoVboMJdZuwpSWoXckt+V3Lp6IzTuoU78hzz39ghSxcHndu/eiZ49e/LTRwNb9N75SV0IKfPi1\niZiqlXjiiedkswc3Ij+lXEhISCizpr9GUGAgJdZuQimw400kL2HKyCRq/RbK3VKRyROn4uPjo3ck\nIa7JiJGPYa5WhU9XmEhOu/73Wywwd6kJk8mP8RNm4OfnZ/+QwmGkjF1MdHQ0UyZNw5SfT9S6zcVD\nIsU1UQoLiVq7iUA/f6ZNnkFAgCw5KtyHyWTi2eengOLDvG+NXO8SEN//YeDUeRjzxHOULFnKMSGF\nw0gZu6AaNWrx5Jin8E1IImTPAb3juAdNI3zTDkw5uUyeOFV+GAm3VLp0GYYNf4TDJxV+33ntl5jP\nJRaXcfv27WnVqq0DEwpHkTJ2UZ06dbm425OK74UEveO4vMCjJwg4c46hQ4ZTu3ZdveMIccO6detJ\n7Vo1WPabidy8q79e02DRj0YC/P0ZOXKM4wMKh5AydmGjH36U6DJlidqwHaWgUO84LsuYkUX49r3U\nbdCIfv0G6B1HiJuiKAoPjxpLdp7Git+v/iN671GFgycUBt03jLCwMCckFI4gZezC/P0DmPD8JAwF\nBYTu9OqtoS9P04jYvAM/Pz+efUpGTgvPULVqDB06dGT1NgMZV1idS9Ng+WojpUtF0aNHH+cFFHYn\nP7lcXExMNW7vN4CgY6dk/vElBB49iW9iMqNGPEJkZKTecYSwm3vuGUKRReGXTZf/Mb3/mMLpCzDw\n7gcxmUxOTCfsTcrYDQwaNJiI6JJEbNsto6v/RikoJGz3fmrUrkvnzt30jiOEXZUrV54WzZuxboeR\nwqJLv+bXLQbCw4Lp2LGzc8MJu5MydgP+/v6MHjEaY3omgcdO6R3HZYTsO4RSWMRjox+XxQ2ER+rT\n9y5y8jS2H/zvn+/kdNh3VKFHz/4yn94DSBm7iZYt2xBToyZhew+iWCx6x9GdMTuHYPU4HTt1oXLl\nqnrHEcIh6tatT6mSEWzY/d8f1Zv2FBe0XBXyDFLGbkJRFEYMfRglL5/Aoyf1jqO74AMqBkVh8P1D\n9Y4ihMMoikL7Dt05dFIhO/ef39t20EitmtUoVaq0PuGEXUkZu5E6depSvVZtQg8ewW7bu7ghQ24e\nQcdOcWvnbkRHR+sdRwiHatWqDZoGu4/8/1J1chqcjYdWrW/VMZmwJyljN/PAvYNRcvMIOHVW7yi6\nCVKPgaZxt2yWLrxATEx1wsOCOHDs/2W8/3jx102bNtcrlrAzKWM307BhY6LLliVYPc51L17rCaxW\ngo+domlsC0qXLqN3GiEcTlEU6jdowqGT/1+v+tBJhciIYMqXr6BvOGE3UsZuRlEU7uw3AJ+UNHyS\nU/WO43QBp+JQ8gu4vW9/vaMI4TR16jQgI1sjOb348bGzRmrXaSizCDyIlLEb6tixC0ZfXwKPn9Y7\nitMFHT9FZMlS1K/fUO8oQjiN2VwTgJPnFDKzITVDo0aNOjqnEvYkZeyGAgMDad2qDYGnz4LFewZy\nGbOy8U1Iome3HnJGILxKxYqVMRoVzsYrnE0o/rNftWo1nVMJe5IydlNdO9+GUliE/7kLekdxmoCT\nxYPWOnXqonMSIZzL19eXMqVLcD65eLtEgIoVK+maSdiXlLGbqlevAYGhoV41qjrodBwxNWrJXsXC\nK5UvX5n4ZCMJKQpBgX6EhYXrHUnYkZSxmzIajXRq34mAc/Eol1u41oOY0jMwpmfQpaPMqxTeqXSZ\nW0hO10hKg1KlouVWjYeRMnZjbdt0AKsVv3PxekdxOP/T50BRaNWqnd5RhNBFqVKlKCyCpHQDpUqV\n0zuOsDMpYzdWq1bt4kvVZ+L0juJwQWfOEWOuKdskCq8VFVW82lx6FpSIljn2nkbK2I0ZDAbat2lP\nwPkEj948wpiZhTE9g07tOugdRQjdREQU/yKaXwAREVE6pxH2JmXs5lq3agsWC37nE/SO4jABZ84B\n0KpVW52TCKGf8PCwv74OCwu7wiuFOzLpHUDcnLp16+MXGIT/mXPk3+KZ95ECzpynQpWqREeX1DuK\nELopU6Ycdw64i+zsTJo3b6V3HGFnUsZuzmQy0bpVG35bt4p0qxWMRr0j2ZUxOweflFQ697lD7yhC\n6MpoNDLkwZF6xxAOIpepPUC7Nu1RCovwu+B5l6r9L16ibt1aRlELITyX3c6MzWazAsQBRy4+tUlV\n1eft9fni8ho0aIRvQAABp+MoKF9W7zh2FXg6jnKVKlOmjGf9ewkhxN/Z8zJ1VWCHqqq97fiZ4hr4\n+PjQpnU7Vq39jXSLFUyecanamJWNT3Iq3foM0DuKEEI4lD0vUzcGypnN5tVms/kHs9lc3Y6fLa6i\nU4dbociCf9x5vaPYzZ9rUbdtK1OahBCe7YbOjM1m81BgzL+eHgVMU1V1mdlsbgV8BsTeZD5xjerV\na0BweDj5J06TX8kDNhzXNIJOnqZ6rTqyFrUQwuPdUBmrqjoPmPf358xmcwBgufj9DWaz+Yo3+SIi\nAjF5yOVUV9G/Tx8+XbAAQ24etsAAvePcFN+kFIyZ2dxz5x1ER4foHUcIIRzKnveMJwCpwCyz2Vwf\nOHOlF6el5drx0AKgdetOfPrppwQeP0123Rp6x7kpgUdPYvLzo169WJKSsvSOI4QQdnG5kwt73jOe\nAbQ1m81rgFeBwXb8bHENypevgLl2XYKPnQRN0zvODVMKCgk8HcetHTsTEODeZ/hCCHEt7HZmrKpq\nBtDLXp8nbsztvfsxffpk/M7FU1DePReTDzx+CqxWevXso3cUIYRwCln0w8O0bNmG4PAIgg8d1TvK\njbHZCFGPE1OjFlWqxOidRgghnELK2MOYTCbu6Nsfv/hETGkZese5bv5nz2PIzmHgHQP1jiKEEE4j\nZeyBunfvhcnXl+ADqt5Rro+mEXLgCBElS9K8eUu90wghhNNIGXugkJBQevboTeCpsxizsvWOc838\nLiTik5LKfQPvw+hhG14IIcSVSBl7qDv634XBYCB4v5ucHWsaIfsPERweQadOXfROI4QQTiVl7KGi\nokrQrVtPgk6cdouzY9/4RHwTkhk08F58fX31jiOEEE4lZezB7rl7EEaDgZB9h/SOcmWaRtieg4RE\nRHDbbTI7TgjhfaSMPVhUVAl69upL4IkzmNIz9Y5zWX5xF/BJSmHwoMFyViyE8EpSxh7u7rsG4ePn\nR8iu/XpHuTSbjfDd+ylRugxdunTXO40QQuhCytjDhYWFcc9d9xIQdx7fhCS94/xH4PHTGNMzGTls\nJCaTPZdKF0II9yFl7AX69r2D4PAIwnbsdak1q5WiIsL2HKRKdTMtW7bRO44QQuhGytgL+Pv7M3LY\nSHxS0gg4ecXNtJwq+ICKkpfHow8/hqIoescRQgjdSBl7iQ4dbqVClaqE7TqAUmTROw7G7BxCDh6l\nVdv21KhRS+84QgihKyljL2EwGHh81OMYcnNdYpnM0J37MBoMjBj2sN5RhBBCd1LGXqR27bq0atue\nkINHMGbn6JbDNyGJgNNxDLzzHqKjS+qWQwghXIWUsZcZMWwURqOR0B179Qlg0wjftofQyCjukJ2Z\nhBACkDL2OtHR0cVTnc6cw/dCotOPH3jsBKa0dB4Z+Qj+/v5OP74QQrgiKWMvdMcdAwkvUYLwHXvA\nZnPacZWCQsJ2H6R6rdq0bt3OaccVQghXJ2XshXx9fXlk5KOY0jIIPHrSaccN2XsQpbCQx0ePkalM\nQgjxN1LGXqplyzbUqF2XsD0HUQoKHX48U3omwepxunbrQZUqMQ4/nhBCuBMpYy+lKAqPjR6DUlhI\nyN6DDj9e2I49+Pr7M+SBYQ4/lhBCuBspYy9WuXIVunTtTvCRExgzshx2HL9z8fidT+CBQYMJCwtz\n2HGEEMJdSRl7uSEPDMXHx5fQnfsccwCbjfCde4ksWYpevfo55hhCCOHmpIy9XHh4BPcMdNyuToHH\nT2FMz+Th4Q/j4+Nj988XQghPIGUs6Nv3DkIiIgjbuc+uuzopRRbC9hyiSnUzrVq1tdvnCiGEp5Ey\nFlxg3vAAAAmPSURBVPj7+zP0gaH4JKfif/a83T43SD2GkpfHw8NHy1QmIYS4AiljAcCtt3ajROky\nhO05ALabPztWCgsJOXiE+o2aUKdOXTskFEIIzyVlLAAwGo0MGzwMY3om/mfibvrzgg4fQykoZNiQ\nh+yQTgghPJuUsfhL69btiC5blrB9h2/q3rFSWETIoWM0jm1OTEw1OyYUQgjPJGUs/mI0Gnnwvgcx\npmfc1L3jIPU4SmEhg+8bYsd0QgjhuaSMxT+0adOeiOhoQg4eubGzY6uVEPUYdRo0JCamuv0DCiGE\nB5IyFv9gNBq558578ElKwTcx+brfH3jiDEpePvfcea8D0gkhhGeSMhb/ceut3fALDCLo8LHre6Om\nEXL4GGUrVqJBg0aOCSeEEB5Iylj8h7+/Pz1u60HA2fMYcnKv+X2+ickY0zO4s98AmVcshBDXQcpY\nXFKvnn0BCDp27fsdBx45gV9gIO3adXRULCGE8EhSxuKSSpcuQ+269Qk6ceaaBnIpBYUEnj1P505d\n+F979xsjR13Hcfy9t7d3e9frHdfrpVdaoC1Jv6004QARiiQUaAoohhBMiEoikUThgQaN0WhUNEae\naDTEKERThQcYQ0SjjzRgpAgRlGBISPAXTYoU8eiV/r0/beG6Pti9djna6+7d2rnOvV/JZmd+M7Pz\nnQe//ezs7vymXC6fgQolKT/mHMYRcWtEPFo3f2VEPBcRz0TEN1pTnrJ080030zY2TsfI7tOu2/Xq\nLpia4oZtHzoDlUlSvswpjCPiAeB+oP6HwQeBj6WUrgauiIjhFtSnDG3efDWlcpnunbtOu273ztdY\nsfo8B/mQpDmY65nxs8A91MI4InqBzpTS9A+MfwC2zr88Zamzs5OrNn+Qrl1vwNSxU67XNj5Bx+hb\n3Hj9tjNYnSTlR/tsCyPiLuDeGc13ppQei4gtdW29wMG6+UPAupZUqExdt+V6dvzpj3SOvMmRVStP\nuk7Xa/8BqgOGSJKaN2sYp5S2A9sbeJ2DwNK6+V5g/2wb9Pd3095ebOCllaWtW6/h2/d3Un595JRh\nXH79DYbOW83w8MYzXJ0k5cOsYdyolNLBiDgaEeuAncA24JuzbbNvX+PXrypbw8OX8teXX+JAZRhm\nXD9cOPo2nbv3cPVttzM6eiijCiXp7DA4uPSk7fO5tKlSe0y7G3gUeB54MaX0t3m8thaQq664irbx\nCdoPvjdsO0d2w7EKH7j8ygwqk6R8mPOZcUppB7Cjbv55YHMritLCMj20ZcfIbt7p633Xso6RUYql\nEhs2vC+L0iQpFxz0Q6c1NLSS3oEBOkdG37Os681RNmy8iFKplEFlkpQPhrFOq1AoMLzpYsp79r5r\nNK7C0bcp7j/ApRdfkmF1knT2M4zVkE0XbaIwMUmx7sYRHXveAmDjxouyKkuScsEwVkPWr69etlR6\na9/xtunp9esjk5okKS8MYzVk7dp1FNraKO07cfl4ae9++pYPsmRJT4aVSdLZzzBWQzo6OhhceS6l\nvQdOtO0/SDgWtSTNm2Gshl24Zh0d09caT01RPDTGhWsvzLYoScoBw1gNW7dmLW1j4/DOFO0Hx6BS\n4fzzL8i6LEk66xnGatiqVauhUqF9bJz2Q2Mn2iRJ82IYq2ErV54LQHFsjOLYeK1tVZYlSVIuGMZq\n2NBQ9a5N7WMTtI+NUyp30dPjP6klab4MYzWsr+8c2opF2iYmaBufpH/5QNYlSVIuGMZqWKFQoOec\nforjk7RPTrJicEXWJUlSLhjGasrAwADFycO0Tx5mxfLBrMuRpFwwjNWU5f0DtB8+QuHwYQaWLcu6\nHEnKBcNYTVm+bIC2w0fgWIW+vv6sy5GkXDCM1ZTe3l4KR44cn5YkzZ9hrKYsXdpbN700w0okKT8M\nYzWl/rri7m6vMZakVjCM1ZTu7iXHp5csWTLLmpKkRhnGakp9AHd3d2dYiSTlh2GsppTL5ePTXV1d\nGVYiSflhGKsp9WFcLhvGktQKhrGa0tl5IoxLpVKGlUhSfhjGakpHRycAbcUihUIh42okKR8MYzWl\np6eHtmKRbm+dKEktU6hUKpnseHT0UDY71ryNjIzQ0dHBMsemlqSmDA4uPelXiu1nuhCd/YaGhrIu\nQZJyxa+pJUnKmGEsSVLGDGNJkjJmGEuSlDHDWJKkjBnGkiRlzDCWJCljhrEkSRmb86AfEXEr8NGU\n0ifq5r8L7Kqtcl9K6en5lyhJUr7NKYwj4gFgG/D3uuZLgS+llH7disIkSVos5vo19bPAPUD9GJuX\nAZ+KiKcj4nsRUZx3dZIkLQKznhlHxF3AvTOa70wpPRYRW2a0PwH8JqX0akQ8BNwN/KhllUqSlFOz\nhnFKaTuwvcHX+llK6UBt+rfAbbOtfKo7V0iStNi05N/UEVEAXoqIVbWmrcALrXhtSZLybj5hXKk9\nSClVgLuAxyPiKaAT+Om8q5MkaREoVCqVrGuQJGlRc9APSZIyZhhLkpQxw1iSpIwZxpIkZcwwVstE\nxFMREVnXIS02jfa9iNgSET88EzWpOYaxWun45W6SzqhG+579c4Ga812bdHaJiPXAz4G3qX4I+zjw\nBvATYDWwEvhdSunrEfEwcBS4gOo1478EPgKcD9xSe/4CUAZWAA+mlB6q21cf1ZHbltWaPpdSerlu\nedvJ9vt/OXApYwup79VcEhFPAH3Aj1NKD7f4kDUHnhkvHluB52rP91HtiOcBf0kp3QhcQXU8cah+\net6ZUroBeAVYk1L6MPA41TeGCrAcuAnYDHwxIgZr2xaArwJPppSuAz4DPDijllPtV8qjhdT3CsAx\nqnfduwb4SkQsb/kRq2meGS8e24EvA78HDlDttPuAyyPiWuAg1U/i016sPe+n+qZAbf1ybXpHSmkK\nmIiIl4F1ddtuAq6NiNtr8/0zapltv1LeLKS+VwGeqY2aOBkRr1A9C98zv0PUfHlmvHjcAvw5pbQV\n+BXVN4dPAvtTSncA3we6T7HtyW7q8X6AiOgGNgL/rFv2D+AHKaVrgTuAR2Zse2eD+5XyYCH1vQLV\nDwGFiOgBNgD/mtNRqaU8M148XgAeiYijVD+EfZ7qb1O/iIjLgH8DL0TEubX16//ocbLp3trvTv3A\nt1JKe2t/5qwA3wG2R8SngV6qX83Ve/Ik+12ZUvpvi45VWkgWUt+b/qPXE7XlX6u7254y5NjUalrt\nXta3pZQ+m3Ut0mJi38svv6bWXHgJk5QN+15OeWYsSVLGPDOWJCljhrEkSRkzjCVJyphhLElSxgxj\nSZIy9j8/u9FtuiuBhQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3729dd8>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Find the mean for both samples\n",
      "xbar_a = sample_a.mean()\n",
      "xbar_b = sample_b.mean()\n",
      "\n",
      "xbar_a, xbar_b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "(2.4833333333333338, 6.0333333333333341)"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Hypothesis\n",
      "Is there a statistical significant difference between the two sample means at an alpha level of .05?\n",
      "\n",
      "|Hypothesis||$\\alpha = .05$|\n",
      "|---:|:---:|---|:---|\n",
      "|$H_0:$|$\\mu_{a} = \\mu_{a}$|There is no significant difference between a and b|\n",
      "|$H_a:$|$\\mu_{b} \\neq \\mu_{b}$|There is a significant difference between a and b|"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Sum of Squares\n",
      "First we need to find the sum of squares from both samples. We do this by summing all the differences from the mean in the sample, then square them. The following formula applies:\n",
      "\n",
      "$$SS=\\sum(x_i-\\bar{x})^2$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Find the sum of squares for both samples\n",
      "ss_a = np.sum((sample_a - xbar_a) ** 2)\n",
      "ss_b = np.sum((sample_b - xbar_b) ** 2)\n",
      "\n",
      "ss_a, ss_b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "(83.48833333333333, 16.806666666666668)"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Pooled Variance\n",
      "Since the sample size between a and b is different, we can use the pooled variance to account for these difference in calculations. The following formula applies:\n",
      "\n",
      "$$S_p^2=\\frac{SS_a+SS_b}{df_a + df_b}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the Degrees of Freedom for both samples\n",
      "n_a = len(sample_a)\n",
      "n_b = len(sample_b)\n",
      "dof_a = n_a - 1\n",
      "dof_b = n_b - 1\n",
      "\n",
      "dof_a, dof_b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "(5, 2)"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the Pooled Variance\n",
      "Sp2 = (ss_a + ss_b) / (dof_a + dof_b)\n",
      "Sp2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "14.327857142857143"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Standard Error\n",
      "Now we can calculate the Standard Error by using the Pooled Variance instead of the individual variance. The following formula applies:\n",
      "\n",
      "$$sem=\\sqrt{\\frac{S_p^2}{n_a}+\\frac{S_p^2}{n_b}}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the Standard Error of the Mean\n",
      "sem = np.sqrt(Sp2 / n_a + Sp2 / n_b)\n",
      "sem"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "2.6765516194216339"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## The t-Statistic\n",
      "The t-statistic for the mean difference between these samples is calculated as follows:\n",
      "\n",
      "$$t=\\frac{(\\bar{x}_a-\\bar{x}_b)-(\\mu_a-\\mu_b)}{sem}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the t-statistic\n",
      "observed_diff = xbar_a - xbar_b\n",
      "expected_diff = 0  # we expect no difference between population means\n",
      "\n",
      "t = (observed_diff - expected_diff) / sem\n",
      "t"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "-1.3263334711127697"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## The t-critical value\n",
      "To test our hypothesis, we need to get the t-critical value."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alpha = .05\n",
      "t_critical = np.abs(stats.t.ppf(alpha / 2, dof_a + dof_b))  # two sided\n",
      "t_critical"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "2.3646242510103002"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Hypothesis testing"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# We compare our t-statistic against the t-criticalk value to test the \n",
      "# acceptance of the alternative hypothesis\n",
      "-t_critical > t or t > t_critical"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "False"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since the t-statistic lies between the t-critical values of alpha, we fail to reject the null hypothesis. The mean difference from both samples is probably due to chance and therefore not statistical significant."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Using scipy.stats or statsmodels"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# A quicker way to get these results is to call scipy stats\n",
      "# ttest_ind on both samples\n",
      "t_statistic, p = stats.ttest_ind(sample_a, sample_b)\n",
      "t_statistic, p"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "(-1.3263334711127697, 0.22636455334070857)"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Another alternative is using Statsmodels\n",
      "t_statistic, p, dof = sm.stats.ttest_ind(sample_a, sample_b, \n",
      "                   alternative='two-sided', usevar='pooled')\n",
      "t_statistic, p, dof"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "(-1.3263334711127697, 0.22636455334070857, 7.0)"
       ]
      }
     ],
     "prompt_number": 12
    }
   ],
   "metadata": {}
  }
 ]
}