{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# HIDDEN\n",
    "from datascience import *\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plots\n",
    "plots.style.use('fivethirtyeight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Further Generalization: New Test Statistics\n",
    "\n",
    "Thus far, we have been looking at entire distributions, whether the variable is categorical, such as employment status, or quantitative, such as age. When we analyze quantitative variables, we can summarize the distributions by using measures such as the median and the average. These measures can be useful as test statistics – sometimes, working with them can be more efficient than working with the whole distribution. We will now generalize our methods of testing so that we can replace the total variation distance between distributions by other test statistics of our choice.\n",
    "\n",
    "As part of the survey we have been studying, the sampled people were asked to rate their satisfaction with their spouse or partner. One of the questions they answered was, \"Every relationship has its ups and downs. Taking all things together, how satisfied are you with your relationship with your spouse or partner?\"\n",
    "\n",
    "The question was in a multiple choice format, with the following possible answers:\n",
    "\n",
    "- 1: very satisfied\n",
    "- 2: somewhat satisfied\n",
    "- 3: neither satisfied nor dissatisfied\n",
    "- 4: somewhat dissatisfied\n",
    "- 5: very dissatisfied\n",
    "\n",
    "The answers given by the sampled people are in the column ``rel_rating`` of the table ``couples``. More than 63% of the sampled people gave their satisfaction the highest possible rating:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table border=\"1\" class=\"dataframe\">\n",
       "    <thead>\n",
       "        <tr>\n",
       "            <th>mar_status_code</th> <th>gender_code</th> <th>rel_rating_code</th> <th>age</th> <th>mar_status</th> <th>gender</th> <th>rel_rating</th> <th>count</th>\n",
       "        </tr>\n",
       "    </thead>\n",
       "    <tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>1          </td> <td>1              </td> <td>51  </td> <td>married   </td> <td>male  </td> <td>very satisfied    </td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>2          </td> <td>1              </td> <td>53  </td> <td>married   </td> <td>female</td> <td>very satisfied    </td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>1          </td> <td>1              </td> <td>57  </td> <td>married   </td> <td>male  </td> <td>very satisfied    </td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>2          </td> <td>1              </td> <td>57  </td> <td>married   </td> <td>female</td> <td>very satisfied    </td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>1          </td> <td>1              </td> <td>60  </td> <td>married   </td> <td>male  </td> <td>very satisfied    </td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>2          </td> <td>1              </td> <td>57  </td> <td>married   </td> <td>female</td> <td>very satisfied    </td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>1          </td> <td>1              </td> <td>62  </td> <td>married   </td> <td>male  </td> <td>very satisfied    </td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>2          </td> <td>1              </td> <td>59  </td> <td>married   </td> <td>female</td> <td>very satisfied    </td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>1          </td> <td>2              </td> <td>53  </td> <td>married   </td> <td>male  </td> <td>somewhat satisfied</td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "        <tr>\n",
       "            <td>1              </td> <td>2          </td> <td>2              </td> <td>61  </td> <td>married   </td> <td>female</td> <td>somewhat satisfied</td> <td>1    </td>\n",
       "        </tr>\n",
       "    </tbody>\n",
       "</table>\n",
       "<p>... (2058 rows omitted)</p"
      ],
      "text/plain": [
       "<Table(8 rows): | mar_status_code | gender_code | rel_rating_code | age | mar_status | gender | rel_rating | count |>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "columns = ['mar_status', 'gender', 'rel_rating', 'age']\n",
    "couples = Table.read_table('married_couples.csv').select(columns)\n",
    "\n",
    "def describe(column, descriptions):\n",
    "    \"\"\"Relabel a column of codes and add a column of descriptions\"\"\"\n",
    "    code = column + '_code'\n",
    "    couples.relabel(column, code)\n",
    "    couples[column] = np.choose(couples[code]-1, descriptions)\n",
    "    \n",
    "describe('mar_status', ['married', 'partner'])\n",
    "describe('gender', ['male', 'female'])\n",
    "describe('rel_rating', [\n",
    "    'very satisfied', \n",
    "    'somewhat satisfied', \n",
    "    'neither satisfied nor dissatisfied', \n",
    "    'somewhat dissatisfied', \n",
    "    'very dissatisfied', \n",
    "])\n",
    "couples['count'] = 1\n",
    "couples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6353965183752418"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "couples.where('rel_rating', 'very satisfied').num_rows/couples.num_rows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us examine whether this rating is related to age. Perhaps older people are more appreciative of their relationships, or more frustrated; or perhaps high satisfaction with the spouse or partner has nothing to do with age. As a first step, here are the histograms of the ages of sampled people who gave the highest rating, and those who did not. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEqCAYAAABnZEX7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9UVWW+P/D3ES6IGiMefhwQsRyPl44JaHKARgvIssuq\nkLUgMosrkTMKrKFRGsRpTOwKSsaFKz/yAvKNWeIKljigabVunqtgBEy1YGLkx2TI166eQxCSXIE8\n8v3Dxf564nhAYdrxnPdrLf7w2Z9nP/sDytt99tn7KPr6+kZAREQkmBlyHwAREdE/AgOOiIiExIAj\nIiIhMeCIiEhIDDgiIhISA46IiITEgCMiIiFNKOCKiorg4+MDlUqF4OBg1NXVWaxvaWlBWFgY3N3d\nodFokJmZabK9trYWTz75JBYtWgR3d3dotVocOHBgzH6qqqoQEBAANzc3BAYG4sSJE3fRGhERWbNx\nA66yshKpqalITk5GTU0NtFotoqKicOnSJbP1/f39iIiIgEqlgk6nQ0ZGBg4cOIDc3FypZs6cOdiy\nZQtOnTqF+vp6JCcnY+/evSgsLJRqGhoaEBcXh+joaNTW1iIqKgobN27EZ599NgVtExGR6BTjPcnk\n8ccfx7Jly5CdnS2NPfzwwwgPD8fOnTvH1BcXFyMtLQ0dHR2wt7cHAOzfvx+HDh3C3/72tzuu8+KL\nL8LBwUEKudjYWFy9ehWVlZVSzbp16+Ds7IyioqK765KIiKyOxTO44eFhNDU1ISQkxGQ8NDQU9fX1\nZuc0NDQgKChICrfR+suXL6Orq8vsnKamJjQ2NiI4OFgaa2xsvKt1iYiIbmdraWNPTw+MRiNcXV1N\nxp2dnWEwGMzOMRgM8PT0NBlzcXGRtnl5eUnjGo0GPT09+OGHH/Daa69hw4YNJvv58bouLi53XJeI\niOh2FgPuXigUignXfvDBBxgYGEBjYyP++Mc/wsXFBa+88spUHxIREVkhiwGnVCphY2Mz5qypu7sb\nbm5uZue4urqarR/ddrvRs7kHH3wQBoMB//Ef/yEF3J328+N9EBERmWPxGpydnR38/Pyg0+lMxnU6\nHQICAszO0Wq1qKurw9DQkEm9h4eHycuTP2Y0GnHz5k2T/ZhbNzAw0NIhExERAZjAbQIJCQkoKytD\naWkp2trakJKSAoPBgNjYWABAWloawsPDpfrIyEg4ODggPj4e58+fR3V1NXJychAfHy/VHDx4EB9+\n+CG++uorfPXVVygtLUVeXh6ef/55qWbz5s04e/YssrOz0d7ejqysLNTW1mLLli1T2f+019HRIfch\nyMZae7fWvgH2Tndn3GtwERER6O3txf79+6HX66HRaFBeXi69kUSv16Ozs1Oqd3R0xLFjx5CcnIyQ\nkBA4OTkhMTERCQkJUs3Nmzexa9cudHV1wcbGBosWLcKuXbuk0ARuncEVFxdjz549SE9Px6JFi1BS\nUoIVK1ZMYftERCSqce+Do5+3jo4OqNVquQ9DFtbau7X2DbB3a+39XvFZlEREJCQGHBERCYkBR0RE\nQmLAERGRkBhwREQkpCl/VBcR0T/MDDu0tF+UZWkX5S/gqpwry9p0bxhwRDRt9F69hqyialnW3r3t\nJXT3XJVlbQCws7GTbe3pigFHRDQBvX3fI7v4z7Ktv/WVZ2Vbe7riNTgiIhISA46IiITEgCMiIiEx\n4IiISEgMOCIiEhIDjoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiExIAjIiIhMeCIiEhIDDgiIhIS\nA46IiITEgCMiIiEx4IiISEgMOCIiEhIDjoiIhGQr9wEQ0V2aYYeW9ouyLe+i/AVclXNlW59oohhw\nRNNM79VryCqqlm39na++wICjaYEvURIRkZAmFHBFRUXw8fGBSqVCcHAw6urqLNa3tLQgLCwM7u7u\n0Gg0yMzMNNleXV2NiIgILF68GAsWLMCaNWtw6tQpk5rDhw/DycnJ5GvevHkYHh6+yxaJiMgajRtw\nlZWVSE1NRXJyMmpqaqDVahEVFYVLly6Zre/v70dERARUKhV0Oh0yMjJw4MAB5ObmSjWffPIJgoOD\nUVFRgZqaGjzxxBN48cUXxwTnrFmz0NHRgfb2drS3t6OtrQ12dnaTbJmIiKzBuNfg8vLysGHDBsTE\nxAAAMjMz8fHHH+PQoUPYuXPnmPqKigoMDg6ioKAA9vb28Pb2RkdHB/Lz85GYmAgA2Lt3r8mclJQU\nfPTRR3j//fcRFBQkjSsUCjg7O0+qQSIisk4Wz+CGh4fR1NSEkJAQk/HQ0FDU19ebndPQ0ICgoCDY\n29ub1F++fBldXV13XOv777+Hk5OTydj169exbNkyLF26FNHR0Whubh63ISIiImCcgOvp6YHRaISr\nq6vJuLOzMwwGg9k5BoNhTL2Li4u0zZzCwkJcuXIF0dHR0tiSJUuQl5eHI0eOoKioCDNnzsRTTz2F\nCxcujN8VERFZvSm/TUChUNxVfVVVFd544w2UlJTA09NTGvf394e/v7/054CAAKxevRoHDx7Evn37\npux4iYhITBYDTqlUwsbGZsyZV3d3N9zc3MzOcXV1NVs/uu12VVVV2LJlC9555x2sXbvW4oHOmDED\nvr6+Fs/gOjo6LO5DVNbaN2C9vQ8MDMi29rVr12T9vsvV+w3jDVm/74D1/X1Xq9WTmm8x4Ozs7ODn\n5wedTofw8HBpXKfTYd26dWbnaLVa7Nq1C0NDQ9J1OJ1OBw8PD3h5eUl1x44dQ3x8PN555x08++yz\n4x7oyMgIvvzyS/j6+t6xZrLfjOmoo6PDKvsGrLf3+s9bMHv2bNnWnzNnDtTqhbKsLWfvtja2sn7f\nAev8HTcZ494mkJCQgLKyMpSWlqKtrQ0pKSkwGAyIjY0FAKSlpZmEX2RkJBwcHBAfH4/z58+juroa\nOTk5iI+Pl2qOHj2KTZs2YdeuXQgMDIRer4der8d3330n1ezduxenT59GZ2cnmpubkZiYiNbWVrz8\n8stT2T8REQlq3GtwERER6O3txf79+6HX66HRaFBeXi5dL9Pr9ejs7JTqHR0dcezYMSQnJyMkJARO\nTk5ITExEQkKCVFNSUoKbN29i+/bt2L59uzS+atUqHD9+HMCt++mSkpJgMBjg6OgIX19fnDx5EsuX\nL5+q3omISGATepNJXFwc4uLizG7Lz88fM6bRaHDy5Mk77u/EiRPjrpmeno709PSJHB4REdEYfBYl\nEREJiQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdEREJiwBERkZAYcEREJCQGHBERCYkBR0REQmLA\nERGRkBhwREQkJAYcEREJiQFHRERCmtDH5RARjbK1sUFL+0WZVuf/yWniGHBEdFd6+75HdvGfZVk7\ncePTsqxL0xP/O0RERELiGRzRPTD09KG756pMq/P/pUQTwYCj6WuGnWzXgoaGf8C+/ApZ1ubLdEQT\nw4Cjaav36jVkFVXLsvarcetkWZeIJo6vdRARkZAYcEREJCQGHBERCYkBR0REQmLAERGRkBhwREQk\nJAYcEREJiQFHRERCYsAREZGQJhRwRUVF8PHxgUqlQnBwMOrq6izWt7S0ICwsDO7u7tBoNMjMzDTZ\nXl1djYiICCxevBgLFizAmjVrcOrUqTH7qaqqQkBAANzc3BAYGIgTJ07cRWtERGTNxg24yspKpKam\nIjk5GTU1NdBqtYiKisKlS5fM1vf39yMiIgIqlQo6nQ4ZGRk4cOAAcnNzpZpPPvkEwcHBqKioQE1N\nDZ544gm8+OKLJsHZ0NCAuLg4REdHo7a2FlFRUdi4cSM+++yzKWibiIhEN27A5eXlYcOGDYiJiYFa\nrUZmZibc3Nxw6NAhs/UVFRUYHBxEQUEBvL29ER4ejqSkJOTn50s1e/fuRVJSEpYvX477778fKSkp\n8PPzw/vvvy/VFBQU4NFHH8XWrVuhVquxbds2rFq1CgUFBVPQNhERic7iw5aHh4fR1NSE3/72tybj\noaGhqK+vNzunoaEBQUFBsLe3N6nfs2cPurq64OXlZXbe999/DycnJ+nPjY2N+M1vfjNm3cLCQssd\n0U9G3o+MAXgJmYgssRhwPT09MBqNcHV1NRl3dnaGwWAwO8dgMMDT09NkzMXFRdpmLuAKCwtx5coV\nREdHm+znx+u6uLjccV366XX3XMXu7DLZ1ufHxhCRJVP+cTkKheKu6quqqvDGG2+gpKRkTDDerY6O\njknNn67k6vvatWEMDAzIsvYouda/Ybwha+9yrm2tvcvdN2B9v+PUavWk5lsMOKVSCRsbmzFnTd3d\n3XBzczM7x9XV1Wz96LbbVVVVYcuWLXjnnXewdu3aCe3nx/u43WS/GdNRR0eHbH23tF/E7NmzZVl7\nlFzr29rYytq7nGtba+9y9w1Y5++4ybB4EcPOzg5+fn7Q6XQm4zqdDgEBAWbnaLVa1NXVYWhoyKTe\nw8PD5OXJY8eOYfPmzSgoKMCzzz5rdj/m1g0MDBy/KyIisnrjXqVPSEhAWVkZSktL0dbWhpSUFBgM\nBsTGxgIA0tLSEB4eLtVHRkbCwcEB8fHxOH/+PKqrq5GTk4P4+Hip5ujRo9i0aRN27dqFwMBA6PV6\n6PV6fPfdd1LN5s2bcfbsWWRnZ6O9vR1ZWVmora3Fli1bprJ/IiIS1LjX4CIiItDb24v9+/dDr9dD\no9GgvLxcul6m1+vR2dkp1Ts6OuLYsWNITk5GSEgInJyckJiYiISEBKmmpKQEN2/exPbt27F9+3Zp\nfNWqVTh+/DiAW2dwxcXF2LNnD9LT07Fo0SKUlJRgxYoVU9U7EREJbEJvMomLi0NcXJzZbbff3zZK\no9Hg5MmTd9zfRJ9IEh4ebnJ2SERENFG8kYiIiITEgCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIi\nITHgiIhISAw4IiISEgOOiIiExIAjIiIhMeCIiEhIDDgiIhLSlH+iNxERTb3Zs2ahpf2iLGu7KH8B\nV+VcWdaeDAYcEdE08F3/NeT+n4l9EstU2/nqC9My4PgSJRERCYkBR0REQmLAERGRkBhwREQkJAYc\nEREJiQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdEREJiwBERkZAYcEREJCQGHBERCYkBR0REQmLA\nERGRkCYUcEVFRfDx8YFKpUJwcDDq6uos1re0tCAsLAzu7u7QaDTIzMw02a7X6/HKK69Aq9VCqVQi\nPj5+zD4OHz4MJycnk6958+ZheHj4LtojIiJrNW7AVVZWIjU1FcnJyaipqYFWq0VUVBQuXbpktr6/\nvx8RERFQqVTQ6XTIyMjAgQMHkJubK9UMDQ1BqVTid7/7HVauXAmFQmF2X7NmzUJHRwfa29vR3t6O\ntrY22NnZ3WOrRERkTcYNuLy8PGzYsAExMTFQq9XIzMyEm5sbDh06ZLa+oqICg4ODKCgogLe3N8LD\nw5GUlIT8/HypxsvLC/v27cP69esxd+6dPyVWoVDA2dkZLi4u0hcREdFEWAy44eFhNDU1ISQkxGQ8\nNDQU9fX1Zuc0NDQgKCgI9vb2JvWXL19GV1fXXR3c9evXsWzZMixduhTR0dFobm6+q/lERGS9LAZc\nT08PjEYjXF1dTcadnZ1hMBjMzjEYDGPqR8+87jTHnCVLliAvLw9HjhxBUVERZs6ciaeeegoXLlyY\n8D6IiMh62U71Du90Pe1u+fv7w9/fX/pzQEAAVq9ejYMHD2Lfvn1m53R0dEzJ2tONXH1fuzaMgYEB\nWdYeJdf6N4w3ZO1dzrWttXe5+wbk6/3atWuy/J5Rq9WTmm8x4JRKJWxsbMaceXV3d8PNzc3sHFdX\nV7P1o9vu1YwZM+Dr62vxDG6y34zpqKOjQ7a+W9ovYvbs2bKsPUqu9W1tbGXtXc61rbV3ufsG5Ot9\nzpw5UKsXyrL2ZFh8idLOzg5+fn7Q6XQm4zqdDgEBAWbnaLVa1NXVYWhoyKTew8MDXl5e93ygIyMj\n+PLLL6FSqe55H0REZD3GfRdlQkICysrKUFpaira2NqSkpMBgMCA2NhYAkJaWhvDwcKk+MjISDg4O\niI+Px/nz51FdXY2cnJwx97o1NzejubkZ/f396O3tRXNzM1pbW6Xte/fuxenTp9HZ2Ynm5mYkJiai\ntbUVL7/88lT1TkREAhv3GlxERAR6e3uxf/9+6PV6aDQalJeXw9PTE8Ctm7Y7OzulekdHRxw7dgzJ\nyckICQmBk5MTEhMTkZCQYLLfxx57DMCta3YjIyP44IMP4OXlhaamJgC37qdLSkqCwWCAo6MjfH19\ncfLkSSxfvnyqeiciIoFN6E0mcXFxiIuLM7vt9vvbRmk0Gpw8edLiPr/77juL29PT05Genj6RwyMi\nIhqDz6IkIiIhMeCIiEhIDDgiIhISA46IiITEgCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIiITHg\niIhISAw4IiISEgOOiIiExIAjIiIhMeCIiEhIDDgiIhISA46IiITEgCMiIiEx4IiISEgMOCIiEhID\njoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiExIAjIiIh2cp9ADRJM+zQ0n5RlqWHhn+QZV0ioolg\nwE1zvVevIauoWpa1X41bJ8u6REQTwZcoiYhISAw4IiIS0oQCrqioCD4+PlCpVAgODkZdXZ3F+paW\nFoSFhcHd3R0ajQaZmZkm2/V6PV555RVotVoolUrEx8eb3U9VVRUCAgLg5uaGwMBAnDhxYoJtERGR\ntRs34CorK5Gamork5GTU1NRAq9UiKioKly5dMlvf39+PiIgIqFQq6HQ6ZGRk4MCBA8jNzZVqhoaG\noFQq8bvf/Q4rV66EQqEYs5+GhgbExcUhOjoatbW1iIqKwsaNG/HZZ59Nol0iIrIW4wZcXl4eNmzY\ngJiYGKjVamRmZsLNzQ2HDh0yW19RUYHBwUEUFBTA29sb4eHhSEpKQn5+vlTj5eWFffv2Yf369Zg7\nd67Z/RQUFODRRx/F1q1boVarsW3bNqxatQoFBQX32CoREVkTiwE3PDyMpqYmhISEmIyHhoaivr7e\n7JyGhgYEBQXB3t7epP7y5cvo6uqa8IE1Njbe1bpERES3sxhwPT09MBqNcHV1NRl3dnaGwWAwO8dg\nMIypd3FxkbZN1J32czf7ICIi6zXl98GZu572U+no6JBtbTkNDAzIsu4N4w3Z1h5lrb3Luba19i53\n34B8vV+7dk2W369qtXpS8y0GnFKphI2NzZizpu7ubri5uZmd4+rqarZ+dNtE3Wk/lvYx2W/GdFT/\neQtmz54ty9q2NrayrT3KWnuXc21r7V3uvgH5ep8zZw7U6oWyrD0ZFl+itLOzg5+fH3Q6ncm4TqdD\nQECA2TlarRZ1dXUYGhoyqffw8ICXl9eED0yr1ZpdNzAwcML7ICIi6zXuuygTEhJQVlaG0tJStLW1\nISUlBQaDAbGxsQCAtLQ0hIeHS/WRkZFwcHBAfHw8zp8/j+rqauTk5Iy51625uRnNzc3o7+9Hb28v\nmpub0draKm3fvHkzzp49i+zsbLS3tyMrKwu1tbXYsmXLVPVOREQCG/caXEREBHp7e7F//37o9Xpo\nNBqUl5fD09MTwK2btjs7O6V6R0dHHDt2DMnJyQgJCYGTkxMSExORkJBgst/HHnsMwK1rdiMjI/jg\ngw/g5eWFpqYmALfO4IqLi7Fnzx6kp6dj0aJFKCkpwYoVK6aqdyIiEtiE3mQSFxeHuLg4s9tuv79t\nlEajwcmTJy3u87vvvht33fDwcJOzQyIioonisyiJiEhIDDgiIhISA46IiITEgCMiIiEx4IiISEgM\nOCIiEhIDjoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiExIAjIiIhMeCIiEhIDDgiIhISA46IiITE\ngCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiExIAjIiIhMeCIiEhI\nDDgiIhISA46IiITEgCMiIiFNKOCKiorg4+MDlUqF4OBg1NXVWaxvaWlBWFgY3N3dodFokJmZOaam\ntrYWjz32GFQqFfz8/FBSUmKy/fDhw3BycjL5mjdvHoaHh++iPSIisla24xVUVlYiNTUVb7/9NoKC\nglBYWIioqCh8+umn8PT0HFPf39+PiIgIrFq1CjqdDm1tbUhMTMSsWbOQmJgIAOjs7MRzzz2Hl156\nCUVFRairq8O2bdugVCrx7LPPSvuaNWsWmpqaMDIyIo3Z2dlNRd9ERCS4cQMuLy8PGzZsQExMDAAg\nMzMTH3/8MQ4dOoSdO3eOqa+oqMDg4CAKCgpgb28Pb29vdHR0ID8/Xwq4kpISeHh4YN++fQAAtVqN\nv/zlL8jNzTUJOIVCAWdn5ylplIiIrIvFlyiHh4fR1NSEkJAQk/HQ0FDU19ebndPQ0ICgoCDY29ub\n1F++fBldXV1Sjbl9fvHFFzAajdLY9evXsWzZMixduhTR0dFobm6+u+6IiMhqWQy4np4eGI1GuLq6\nmow7OzvDYDCYnWMwGMbUu7i4SNsAoLu722zNjRs30NPTAwBYsmQJ8vLycOTIERQVFWHmzJl46qmn\ncOHChbtoj4iIrNW4L1HeLYVCMSX78ff3h7+/v/TngIAArF69GgcPHpRe2iQiIroTiwGnVCphY2Mz\n5mytu7sbbm5uZue4urqarR/dZqnG1tYWSqXS7H5nzJgBX19fi2dwHR0dltoR1sDAgCzr3jDekG3t\nUdbau5xrW2vvcvcNyNf7tWvXZPn9qlarJzXfYsDZ2dnBz88POp0O4eHh0rhOp8O6devMztFqtdi1\naxeGhoak63A6nQ4eHh7w8vKSak6cOGEyT6fTYcWKFbCxsTG735GREXz55Zfw9fW94/FO9psxHdV/\n3oLZs2fLsratja1sa4+y1t7lXNtae5e7b0C+3ufMmQO1eqEsa0/GuPfBJSQkoKysDKWlpWhra0NK\nSgoMBgNiY2MBAGlpaSbhFxkZCQcHB8THx+P8+fOorq5GTk4O4uPjpZrY2FhcvnwZqampaGtrQ2lp\nKY4cOSK9yxIA9u7di9OnT6OzsxPNzc1ITExEa2srXn755ansn4iIBDXuNbiIiAj09vZi//790Ov1\n0Gg0KC8vl+6B0+v16OzslOodHR1x7NgxJCcnIyQkBE5OTkhMTERCQoJUs3DhQpSXl2PHjh04dOgQ\n3N3dkZmZiWeeeUaq6e/vR1JSEgwGAxwdHeHr64uTJ09i+fLlU9g+ERGJakJvMomLi0NcXJzZbfn5\n+WPGNBoNTp48aXGfv/rVr3DmzJk7bk9PT0d6evpEDo+IiGgMPouSiIiExIAjIiIhMeCIiEhIDDgi\nIhISA46IiITEgCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiExIAj\nIiIhMeCIiEhIDDgiIhISA46IiITEgCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIiITHgiIhISAw4\nIiISEgOOiIiExIAjIiIhMeCIiEhIDDgiIhISA46IiIQ0oYArKiqCj48PVCoVgoODUVdXZ7G+paUF\nYWFhcHd3h0ajQWZm5pia2tpaPPbYY1CpVPDz80NJScmYmqqqKgQEBMDNzQ2BgYE4ceLEBNv66RiN\nRgwN/yDblwIKub8FREQ/S7bjFVRWViI1NRVvv/02goKCUFhYiKioKHz66afw9PQcU9/f34+IiAis\nWrUKOp0ObW1tSExMxKxZs5CYmAgA6OzsxHPPPYeXXnoJRUVFqKurw7Zt26BUKvHss88CABoaGhAX\nF4cdO3bgmWeeQXV1NTZu3IgPP/wQDz/88BR/G+7dN/pe5L8rT/AuXeIF7196yLI2EdHP3bgBl5eX\nhw0bNiAmJgYAkJmZiY8//hiHDh3Czp07x9RXVFRgcHAQBQUFsLe3h7e3Nzo6OpCfny8FXElJCTw8\nPLBv3z4AgFqtxl/+8hfk5uZKAVdQUIBHH30UW7duBQBs27YNNTU1KCgoQFFR0dR0PwVGRkbw9f+9\nIsvaSqf7GHBERHdg8SXK4eFhNDU1ISQkxGQ8NDQU9fX1Zuc0NDQgKCgI9vb2JvWXL19GV1eXVGNu\nn1988QWMRiMAoLGx8a7WJSIiup3FgOvp6YHRaISrq6vJuLOzMwwGg9k5BoNhTL2Li4u0DQC6u7vN\n1ty4cQM9PT0W93OndYmIiG437kuUd0uhsK43PSyc74r38lNlPYb38h+Ube2gh+VbGwBWa5fJtrac\nvcvZN2C9vVvz3/fpyOIZnFKphI2NzZizpu7ubri5uZmd4+rqarZ+dJulGltbWyiVSos1Pz6rIyIi\nMsdiwNnZ2cHPzw86nc5kXKfTISAgwOwcrVaLuro6DA0NmdR7eHjAy8tLqjG3zxUrVsDGxsZiTWBg\n4ARbIyIiazbufXAJCQkoKytDaWkp2trakJKSAoPBgNjYWABAWloawsPDpfrIyEg4ODggPj4e58+f\nR3V1NXJychAfHy/VxMbG4vLly0hNTUVbWxtKS0tx5MgR6V2WALB582acPXsW2dnZaG9vR1ZWFmpr\na7Fly5ap7J+IiASl6OvrGxmvqLi4GDk5OdDr9dBoNEhPT0dQUBAAID4+HufOnUNTU5NU/7e//Q3J\nycn4/PPP4eTkhNjYWPz+97832ee5c+ewY8cOtLa2wt3dHa+++io2btxoUlNVVYU9e/ags7MTixYt\nwuuvv46nn356CtomIiLRTSjgiIiIpptp8yzKc+fO4fnnn4dGo4GTkxPKysrG1GRkZODBBx+Eu7s7\nnn76abS2tspwpFMrKysLISEh8PLywuLFi/H888/j/PnzY+pE7L2wsBC/+tWv4OXlBS8vLzz55JP4\n6KOPTGpE7PvHsrKy4OTkhNdee81kXMTeMzIy4OTkZPLl7e09pka0vkdduXIFmzdvxuLFi6FSqRAY\nGIhz586Z1IjW/7Jly8b8zJ2cnBAdHQ3g1sM07rXnaRNw//u//4uHHnoIGRkZcHBwGHM7QnZ2NvLz\n85GZmYnTp0/DxcUFERERuHbtmkxHPDXOnTuHTZs24aOPPkJ1dTVsbW2xbt069PX1STWi9j5//nzs\n3r0bZ8/5Cxv+AAAJUklEQVSexX//93/j0UcfxYYNG/DXv/4VgLh9366xsRHvvvsuli5davJ3XuTe\nlyxZgvb2dunrk08+kbaJ3HdfXx/Wrl0LhUKBiooKNDQ0IDMzU7qPGBCz/zNnzpj8vM+cOQOFQoGI\niAgAQE5Ozj33PC1fovT09MRbb72F9evXA7iV8N7e3vjNb34jPdprcHAQarUab7755phre9PZwMAA\nvLy8UFZWhrVr11pV7wDwwAMPYNeuXYiJiRG+76tXryI4OBgHDhzA3r17pQeXi/wzz8jIwPHjx01C\nbZTIfQPA7t27UVdXh1OnTpndLnr/o/bv34/c3Fy0tbXBzs5uUj1PmzM4Sy5evAiDwYDQ0FBpbObM\nmXjkkUeEe7TX999/j5s3b2Lu3LkArKd3o9GIo0ePYmhoCI888ohV9P3qq69i3bp1WLVqFUZG/v//\nQ0XvvbOzEw8++CB8fX0RFxeHzs5OAOL3/f7772PFihWIjY2FWq3G6tWrUVhYKG0XvX/gVoj/6U9/\nwnPPPQd7e/tJ9zzlTzKRg16vBwCTU3ng1iPFrlyR50HI/yjbt2+Hj48PtFotAPF7b2lpwZNPPomh\noSE4ODigpKQEarVa+sstat/vvvsuOjs7pQeL3/7ypMg/c39/fxQUFECtVqO7uxtvvfUW1q5di08/\n/VTovoFbwV5cXIyEhARs3boVzc3NSElJAQBs2rRJ+P6BW/c6d3V14V//9V8BTP7vuhABZ4lIjw7b\nsWMHGhoacOrUqQn1JULvS5Yswblz53D16lVUVVUhLi4Ox48ftzhnuvfd0dGBN998Ex988IH04IOR\nkRGTs7g7me69r1mzxuTP/v7+8PX1RVlZGVauXHnHedO9bwC4efMmHn74Yfzxj38EcOvNFxcuXEBR\nURE2bdpkca4I/QO3/mP38MMPY+nSpePWTqRnIV6iHH1s2OgjwUaJ9Giv1NRUHDt2DNXV1Vi4cKE0\nLnrv//RP/4T7778fvr6+2LlzJ1auXInCwkKh+25oaEBPTw8CAwPh7OwMZ2dnfPLJJyguLoaLi4v0\nODsRe/+xWbNmwdvbG19//bXQP3MAUKlU+Od//meTMbVajUuXLgEQ/996d3c3Tp06JX00GzD5noUI\nuIULF8LNzQ2nT5+WxgYHB/Hpp5/e8ZFi00lKSooUbosXLzbZJnrvP2Y0GnHz5k3cf//9wvb99NNP\no66uDrW1taitrUVNTQ2WL1+OyMhI1NTU4Je//KWwvf/Y4OAg2tvb4ebmJvTPHAACAwPR3t5uMvb3\nv/9desSh6P/Wy8rKMHPmTERGRkpjk+3ZZvv27bv+EQc71QYGBtDa2gq9Xo8//elP0Gg0uO+++/DD\nDz/gF7/4BYxGI/793/8dixcvhtFoxB/+8AcYDAZkZ2fDzs5O7sO/Z8nJyXjvvfdQUlKC+fPnY2Bg\nAAMDA1AoFLCzs4NCoRC29127dsHe3h43b97EN998g4KCAlRUVGD37t144IEHhO175syZ0pmbs7Mz\nXFxcUF5ejgULFuCFF14Q+mf++uuvSz/zv//973jttdfw9ddfIzs7G46OjsL2DQALFizAvn37YGNj\nA5VKhTNnzuDf/u3fsHXrVqxYsULon/vIyAgSEhLw1FNP4ZlnnpHGJ9vztLkG9/nnn0uf9q1QKJCR\nkYGMjAy88MILyMvLQ1JSEq5fv47XXnsNfX19WLlyJSorKzF79myZj3xyiouLoVAoTJ73Cdx6s8no\nBWhRezcYDPj1r38Ng8EAR0dHPPTQQzh69Kj0Qbii9m2OQqEwueYgau+XL1/GK6+8gp6eHjg7O8Pf\n3x//9V//BU9PTwDi9g0Ay5cvx+HDh7F792689dZbWLBgAV5//XXExcVJNaL2X1NTg6+//lp6U9Xt\nJtPztLwPjoiIaDxCXIMjIiL6MQYcEREJiQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdEREJiwBER\nkZAYcEREJCQGHBERCYkBRySTCxcu4Ne//jV8fX3h7u4OPz8/bNu2DX19fWNq8/PzsWzZMqhUKjz+\n+OOor6/HsmXLEB8fb1LX2dmJTZs2YfHixXBzc8Pq1atx4sSJn6olop+VafOwZSLRXLlyBfPnz8ee\nPXswb948dHZ2IisrC3/961/x0UcfSXWlpaX4wx/+gJiYGKxbtw4XLlzApk2b0N/fb/IA5kuXLmHN\nmjVwdXVFRkYGnJ2dcfToUcTExODw4cP4l3/5FznaJJINH7ZM9DNx48YNNDY2IiwsDGfOnIGPjw9u\n3ryJZcuW4aGHHsJ7770n1R4/fhwxMTHSp2kAQGJiIj788EM0NjZi7ty5Um1ERAS+/fZb1NTU/OQ9\nEcmJL1ESyWR4eBhvv/02/P394e7uDhcXF4SFhQEAvvrqKwDAN998g//5n/8Z83FJYWFhsLU1fQHm\n448/xhNPPIH77rsPN27ckL5CQ0Px5Zdf4tq1az9NY0Q/E3yJkkgmaWlpKCwsREpKCrRaLe677z5c\nunQJL730EgYHBwEAer0eAODi4mIy18bGBkql0mSsu7sbR44cwZEjR8aspVAo0Nvbizlz5vyDuiH6\n+WHAEcmksrIS69evx7Zt26Sx/v5+kxo3NzcAt8LrdkajEd9++63JmFKpxCOPPIKkpCSz66lUqqk4\nbKJpgwFHJJPr16+PeZnx8OHDJn+eP38+5s+fjz//+c944YUXpPETJ07AaDSa1D7++ONobGyEt7c3\nZs6c+Y87cKJpggFHJJM1a9bgyJEj0Gg0eOCBB3D8+HE0Njaa1MyYMQO///3vkZSUhN/+9rcIDw9H\nZ2cnsrOz4ejoiBkz/v9l9B07duDxxx9HWFgYNm3ahAULFqCvrw/nz5/HxYsXkZub+1O3SCQrBhyR\nTDIzMzEyMoI333wTAPDkk0+iuLgYoaGhJnUxMTEYGBhAfn4+ysvLodFo8J//+Z9Yv349HB0dpTpP\nT0/odDrs3bsXb775Jr799lvMmzcPGo0G69ev/0l7I/o54G0CRNPQF198gdDQUBw8eBDPPfec3IdD\n9LPEgCP6mbt48SIKCwsRFBSE++67D+3t7cjKyoK9vT3q6up4vY3oDvgSJdHPnIODA1pbW/Hee++h\nr68Pc+fORUhICN544w2GG5EFPIMjIiIh8UkmREQkJAYcEREJiQFHRERCYsAREZGQGHBERCQkBhwR\nEQnp/wEW28GphvZQ/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1061eec18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEqCAYAAABnZEX7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1QFHeeP/D3CAeihhV5GlAxcR2PTBTQyABZTYCYmKOy\n4lRBWGOWk6BZBWrJKgmSZI2YFZS4Hpw8xOPpli2xAiUEJJpNXZwzYAhwMQUbVh72DOHnHs4QEFFO\nIAK/Pyz7nDAOKGw6fOf9quIPvv3p/vaHQd9093SPoq+vbwxERESCmSX3DhAREf09MOCIiEhIDDgi\nIhISA46IiITEgCMiIiEx4IiISEgMOCIiEtKkAi4vLw9eXl5QKpUIDAxEbW2t2frm5maEhITAzc0N\narUaaWlpRstramrw7LPPYunSpXBzc4NGo8HRo0fHbaeiogJ+fn5wdXWFv78/qqqq7qM1IiKyZBMG\nXFlZGZKSkpCQkIDq6mpoNBqEh4fj8uXLJuv7+/uh1WqhVCqh0+mQmpqKo0ePIjMzU6qZN28edu7c\niTNnzqCurg4JCQk4ePAgcnNzpZr6+npER0cjIiICNTU1CA8Px9atW/HFF19MQ9tERCQ6xURPMnn6\n6aexcuVKpKenS2OPP/44QkNDsXfv3nH1+fn5SE5ORnt7O2xtbQEAhw8fRkFBAf7yl7/cc56XXnoJ\ndnZ2UshFRUXh2rVrKCsrk2o2bdoEJycn5OXl3V+XRERkccwewQ0PD6OxsRFBQUFG48HBwairqzO5\nTn19PQICAqRwu1Pf1dWFzs5Ok+s0NjaioaEBgYGB0lhDQ8N9zUtERHQ3a3MLe3p6MDIyAhcXF6Nx\nJycnGAwGk+sYDAYsWrTIaMzZ2Vla5uHhIY2r1Wr09PTgu+++w2uvvYYtW7YYbef78zo7O99zXiIi\noruZDbgHoVAoJl370UcfYWBgAA0NDfjtb38LZ2dnbNu2bbp3iYiILJDZgHN0dISVldW4o6bu7m64\nurqaXMfFxcVk/Z1ld7tzNPfoo4/CYDDgX//1X6WAu9d2vr8NIiIiU8xeg7OxsYGPjw90Op3RuE6n\ng5+fn8l1NBoNamtrMTQ0ZFTv7u5udHry+0ZGRjA6Omq0HVPz+vv7m9tlIiIiAJM4RRkbG4tf/epX\nWL16Nfz8/FBQUACDwYCoqCgAQHJyMi5cuICKigoAQFhYGA4dOoSYmBgkJCSgvb0dGRkZSExMlLZ5\n7NgxPPzww1i2bBkA4Pz588jKysIrr7wi1ezYsQMhISFIT09HSEgIqqqqUFNTgz/96U/T+gOY6drb\n26FSqeTeDVlYau91F5pxJK9Stvn3vvoiHlu+RJa5LfU1Byy79wc1YcBptVr09vbi8OHD0Ov1UKvV\nKCkpkd5Iotfr0dHRIdXb29ujvLwcCQkJCAoKgoODA+Li4hAbGyvVjI6OYt++fejs7ISVlRWWLl2K\nffv2SaEJ3D6Cy8/Px4EDB5CSkoKlS5eisLAQq1evnsb2iYhIVBPeB0c/bpb8V52l9s4jOMt7zQHL\n7v1B8VmUREQkJAYcEREJiQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdEREJiwBERkZAYcEREJCQG\nHBERCYkBR0REQmLAERGRkBhwREQkJAYcEREJiQFHRERCYsAREZGQJvxEbyKiu1lbWaG57RtZ5rax\nspFlXpqZGHBEdF96+64jPf8DWebetW2jLPPSzMRTlEREJCQGHBERCYkBR0REQmLAERGRkBhwREQk\nJAYcEREJiQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdEREKaVMDl5eXBy8sLSqUSgYGBqK2tNVvf\n3NyMkJAQuLm5Qa1WIy0tzWh5ZWUltFotli1bhsWLF2P9+vU4c+aMUc3x48fh4OBg9LVgwQIMDw/f\nZ4tERGSJJgy4srIyJCUlISEhAdXV1dBoNAgPD8fly5dN1vf390Or1UKpVEKn0yE1NRVHjx5FZmam\nVPPZZ58hMDAQpaWlqK6uxjPPPIOXXnppXHDOmTMH7e3taGtrQ1tbG1pbW2Fjw6eJExHRxCb8NIGs\nrCxs2bIFkZGRAIC0tDR88sknKCgowN69e8fVl5aWYnBwEDk5ObC1tYWnpyfa29uRnZ2NuLg4AMDB\ngweN1klMTMTHH3+MDz/8EAEBAdK4QqGAk5PTlBokIiLLZPYIbnh4GI2NjQgKCjIaDw4ORl1dncl1\n6uvrERAQAFtbW6P6rq4udHZ23nOu69evw8HBwWjs5s2bWLlyJR577DFERESgqalpwoaIiIiACY7g\nenp6MDIyAhcXF6NxJycnGAwGk+sYDAYsWrTIaMzZ2Vla5uHhMW6d3NxcXLlyBREREdLY8uXLkZWV\nhRUrVuD69et477338Nxzz6GmpgZLly6dXHdEJJS5c+bI9mGrzo4/gYvjfFnmpgcz7R94qlAo7qu+\noqICb7/9NgoLC42C0dfXF76+vtL3fn5+WLduHY4dO4ZDhw5N2/4S0cxxtf8GMv+9Spa59776IgNu\nhjEbcI6OjrCyshp3tNbd3Q1XV1eT67i4uJisv7PsbhUVFdi5cyfee+89bNiwweyOzpo1C97e3rh0\n6dI9a9rb281uQ1SW2jdgub0PDAzINvetkVuyzi/X3Ddu3JD9903u+X9oKpVqSuubDTgbGxv4+PhA\np9MhNDRUGtfpdNi0aZPJdTQaDfbt24ehoSHpOpxOp4O7u7vR6cny8nLExMTgvffew8aNE38M/djY\nGL766it4e3vfs2aqP4yZqL293SL7Biy397oLzZg7d65s81tbWcs6v1xzz5s3DyrVElnmBiz3930q\nJrxNIDY2FsXFxSgqKkJraysSExNhMBgQFRUFAEhOTjYKv7CwMNjZ2SEmJgYXL15EZWUlMjIyEBMT\nI9WcPHkS27dvx759++Dv7w+9Xg+9Xo+rV69KNQcPHsTZs2fR0dGBpqYmxMXFoaWlBS+//PJ09k9E\nRIKa8BqcVqtFb28vDh8+DL1eD7VajZKSEul6mV6vR0dHh1Rvb2+P8vJyJCQkICgoCA4ODoiLi0Ns\nbKxUU1hYiNHRUezZswd79uyRxteuXYtTp04BuH0/XXx8PAwGA+zt7eHt7Y3Tp09j1apV09U7EREJ\nbFJvMomOjkZ0dLTJZdnZ2ePG1Go1Tp8+fc/tVVVNfJE4JSUFKSkpk9k9IiKicfgsSiIiEtK03yZA\nZAkMPX3o7rkm0+z8u5RoMhhwRA+gu+ca9qcXyzJ33NbnZZmXaKbhn4JERCQkBhwREQmJAUdEREJi\nwBERkZAYcEREJCQGHBERCYkBR0REQmLAERGRkBhwREQkJAYcEREJiQFHRERC4rMoiYgmwdrKCs1t\n38g2v42VjWxzz1QMOCKiSejtu470/A9km3/Xto2yzT1T8RQlEREJiQFHRERCYsAREZGQGHBERCQk\nBhwREQmJAUdEREJiwBERkZB4HxzNXLNsZLvxdmj4O1nmJaLJY8DRjNV77QaO5FXKMver0ZtkmZeI\nJo+nKImISEgMOCIiEhIDjoiIhMSAIyIiIU0q4PLy8uDl5QWlUonAwEDU1taarW9ubkZISAjc3Nyg\nVquRlpZmtLyyshJarRbLli3D4sWLsX79epw5c2bcdioqKuDn5wdXV1f4+/ujqqrqPlojIiJLNmHA\nlZWVISkpCQkJCaiuroZGo0F4eDguX75ssr6/vx9arRZKpRI6nQ6pqak4evQoMjMzpZrPPvsMgYGB\nKC0tRXV1NZ555hm89NJLRsFZX1+P6OhoREREoKamBuHh4di6dSu++OKLaWibiIhEN2HAZWVlYcuW\nLYiMjIRKpUJaWhpcXV1RUFBgsr60tBSDg4PIycmBp6cnQkNDER8fj+zsbKnm4MGDiI+Px6pVq/Dw\nww8jMTERPj4++PDDD6WanJwcPPnkk9i1axdUKhV2796NtWvXIicnZxraJiIi0ZkNuOHhYTQ2NiIo\nKMhoPDg4GHV1dSbXqa+vR0BAAGxtbY3qu7q60NnZec+5rl+/DgcHB+n7hoaG+5qXiIjobmYDrqen\nByMjI3BxcTEad3JygsFgMLmOwWAYV+/s7CwtMyU3NxdXrlxBRETEhNu51zaIiIjuNu1PMlEoFPdV\nX1FRgbfffhuFhYVYtGjRlOZub2+f0vozlaX2DQADAwOyzHtr5JZscwPy9Q1Ybu9y9w1Y3r91lUo1\npfXNBpyjoyOsrKzGHTV1d3fD1dXV5DouLi4m6+8su1tFRQV27tyJ9957Dxs2bJjUdr6/jbtN9Ycx\nE7W3t1tk3wBQd6EZc+fOlWVuaytr2eYGIOvcltq73H0Dlvl/3FSYPUVpY2MDHx8f6HQ6o3GdTgc/\nPz+T62g0GtTW1mJoaMio3t3dHR4eHtJYeXk5duzYgZycHGzcuNHkdkzN6+/vP3FXRERk8SZ8F2Vs\nbCyKi4tRVFSE1tZWJCYmwmAwICoqCgCQnJyM0NBQqT4sLAx2dnaIiYnBxYsXUVlZiYyMDMTExEg1\nJ0+exPbt27Fv3z74+/tDr9dDr9fj6tWrUs2OHTvw6aefIj09HW1tbThy5Ahqamqwc+fO6eyfiIgE\nNeE1OK1Wi97eXhw+fBh6vR5qtRolJSXS9TK9Xo+Ojg6p3t7eHuXl5UhISEBQUBAcHBwQFxeH2NhY\nqaawsBCjo6PYs2cP9uzZI42vXbsWp06dAnD7CC4/Px8HDhxASkoKli5disLCQqxevXq6eiciIoFN\n6k0m0dHRiI6ONrns7vvb7lCr1Th9+vQ9tzfZJ5KEhoYaHR0SERFNFp9FSUREQmLAERGRkBhwREQk\nJAYcEREJiQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdEREJiwBERkZAYcEREJCQGHBERCYkBR0RE\nQmLAERGRkCb1cTlEphh6+tDdc03GPeDfZ0R0bww4emDdPdewP71Ytvnjtj4v29xE9OPHP4GJiEhI\nDDgiIhISA46IiITEgCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiE\nxIAjIiIhMeCIiEhIkwq4vLw8eHl5QalUIjAwELW1tWbrm5ubERISAjc3N6jVaqSlpRkt1+v12LZt\nGzQaDRwdHRETEzNuG8ePH4eDg4PR14IFCzA8PHwf7RERkaWaMODKysqQlJSEhIQEVFdXQ6PRIDw8\nHJcvXzZZ39/fD61WC6VSCZ1Oh9TUVBw9ehSZmZlSzdDQEBwdHfGb3/wGa9asgUKhMLmtOXPmoL29\nHW1tbWhra0NraytsbGwesFUiIrIkE34eXFZWFrZs2YLIyEgAQFpaGj755BMUFBRg79694+pLS0sx\nODiInJwc2NrawtPTE+3t7cjOzkZcXBwAwMPDA4cOHQIAfPDBB/ecW6FQwMnJ6YEasxizbNDc9o0s\nUw8NfyfLvEREk2E24IaHh9HY2Ihf//rXRuPBwcGoq6szuU59fT0CAgJga2trVH/gwAF0dnbCw8Nj\n0jt38+ZNrFy5EqOjo1ixYgXefPNNeHl5TXp9S9B77QaO5FXKMver0ZtkmZeIaDLMnqLs6enByMgI\nXFxcjMadnJxgMBhMrmMwGMbVOzs7S8sma/ny5cjKysKJEyeQl5eH2bNn47nnnsOlS5cmvQ0iIrJc\nE56ivF/3up52v3x9feHr6yt97+fnh3Xr1uHYsWPS6c3va29vn5a5Z5qBgQFZ5r01cku2ue+w1N7l\nnNtSe5e7b8Dy/o9TqVRTWt9swDk6OsLKymrckVd3dzdcXV1NruPi4mKy/s6yBzVr1ix4e3ubPYKb\n6g9jJqq70Iy5c+fKMre1lbVsc99hqb3LObel9i5334Bl/h83FWZPUdrY2MDHxwc6nc5oXKfTwc/P\nz+Q6Go0GtbW1GBoaMqp3d3e/r+tv3zc2NoavvvoKSqXygbdBRESWY8LbBGJjY1FcXIyioiK0trYi\nMTERBoMBUVFRAIDk5GSEhoZK9WFhYbCzs0NMTAwuXryIyspKZGRkjLvXrampCU1NTejv70dvby+a\nmprQ0tIiLT948CDOnj2Ljo4ONDU1IS4uDi0tLXj55Zenq3ciIhLYhNfgtFotent7cfjwYej1eqjV\napSUlGDRokUAbt+03dHRIdXb29ujvLwcCQkJCAoKgoODA+Li4hAbG2u03aeeegrA7Wt2Y2Nj+Oij\nj+Dh4YHGxkYAt++ni4+Ph8FggL29Pby9vXH69GmsWrVqunonIpox5s6ZI9stQc6OP4GL43xZ5p6K\nSb3JJDo6GtHR0SaXZWdnjxtTq9U4ffq02W1evXrV7PKUlBSkpKRMZveIiIR3tf8GMv+9Spa59776\n4owMOD6LkoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiExIAjIiIhMeCIiEhIDDgiIhISA46IiITE\ngCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiExIAjIiIhMeCIiEhI\nDDgiIhISA46IiITEgCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiE\nNKmAy8vLg5eXF5RKJQIDA1FbW2u2vrm5GSEhIXBzc4NarUZaWprRcr1ej23btkGj0cDR0RExMTEm\nt1NRUQE/Pz+4urrC398fVVVVk2yLiIgs3YQBV1ZWhqSkJCQkJKC6uhoajQbh4eG4fPmyyfr+/n5o\ntVoolUrodDqkpqbi6NGjyMzMlGqGhobg6OiI3/zmN1izZg0UCsW47dTX1yM6OhoRERGoqalBeHg4\ntm7dii+++GIK7RIRkaWYMOCysrKwZcsWREZGQqVSIS0tDa6urigoKDBZX1paisHBQeTk5MDT0xOh\noaGIj49Hdna2VOPh4YFDhw5h8+bNmD9/vsnt5OTk4Mknn8SuXbugUqmwe/durF27Fjk5OQ/YKhER\nWRKzATc8PIzGxkYEBQUZjQcHB6Ours7kOvX19QgICICtra1RfVdXFzo7Oye9Yw0NDfc1LxER0d3M\nBlxPTw9GRkbg4uJiNO7k5ASDwWByHYPBMK7e2dlZWjZZ99rO/WyDiIgsl/V0b9DU9bQfSnt7u2xz\ny2lgYECWeW+N3JJt7jsstXc557bU3uXuG5Cv9xs3bsjy/6tKpZrS+mYDztHREVZWVuOOmrq7u+Hq\n6mpyHRcXF5P1d5ZN1r22Y24bU/1hzER1F5oxd+5cWea2trKWbe47LLV3Oee21N7l7huQr/d58+ZB\npVoiy9xTYfYUpY2NDXx8fKDT6YzGdTod/Pz8TK6j0WhQW1uLoaEho3p3d3d4eHhMesc0Go3Jef39\n/Se9DSIislwTvosyNjYWxcXFKCoqQmtrKxITE2EwGBAVFQUASE5ORmhoqFQfFhYGOzs7xMTE4OLF\ni6isrERGRsa4e92amprQ1NSE/v5+9Pb2oqmpCS0tLdLyHTt24NNPP0V6ejra2tpw5MgR1NTUYOfO\nndPVOxERCWzCa3BarRa9vb04fPgw9Ho91Go1SkpKsGjRIgC3b9ru6OiQ6u3t7VFeXo6EhAQEBQXB\nwcEBcXFxiI2NNdruU089BeD2NbuxsTF89NFH8PDwQGNjI4DbR3D5+fk4cOAAUlJSsHTpUhQWFmL1\n6tXT1TsREQlsUm8yiY6ORnR0tMlld9/fdodarcbp06fNbvPq1asTzhsaGmp0dEhERDRZfBYlEREJ\niQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdEREJiwBERkZAYcEREJCQGHBERCYkBR0REQmLAERGR\nkBhwREQkJAYcEREJiQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdEREJiwBERkZAYcEREJCQGHBER\nCYkBR0REQmLAERGRkBhwREQkJAYcEREJiQFHRERCYsAREZGQGHBERCSkSQVcXl4evLy8oFQqERgY\niNraWrP1zc3NCAkJgZubG9RqNdLS0sbV1NTU4KmnnoJSqYSPjw8KCwuNlh8/fhwODg5GXwsWLMDw\n8PB9tEdERJbKeqKCsrIyJCUl4fe//z0CAgKQm5uL8PBwfP7551i0aNG4+v7+fmi1WqxduxY6nQ6t\nra2Ii4vDnDlzEBcXBwDo6OjACy+8gF/+8pfIy8tDbW0tdu/eDUdHR2zcuFHa1pw5c9DY2IixsTFp\nzMbGZjr6JiIiwU0YcFlZWdiyZQsiIyMBAGlpafjkk09QUFCAvXv3jqsvLS3F4OAgcnJyYGtrC09P\nT7S3tyM7O1sKuMLCQri7u+PQoUMAAJVKhf/6r/9CZmamUcApFAo4OTlNS6NERGRZzJ6iHB4eRmNj\nI4KCgozGg4ODUVdXZ3Kd+vp6BAQEwNbW1qi+q6sLnZ2dUo2pbX755ZcYGRmRxm7evImVK1fiscce\nQ0REBJqamu6vOyIislhmA66npwcjIyNwcXExGndycoLBYDC5jsFgGFfv7OwsLQOA7u5ukzW3bt1C\nT08PAGD58uXIysrCiRMnkJeXh9mzZ+O5557DpUuX7qM9IiKyVBOeorxfCoViWrbj6+sLX19f6Xs/\nPz+sW7cOx44dk05tEhER3YvZgHN0dISVldW4o7Xu7m64urqaXMfFxcVk/Z1l5mqsra3h6Ohocruz\nZs2Ct7e32SO49vZ2c+0Ia2BgQJZ5b43ckm3uOyy1dznnttTe5e4bkK/3GzduyPL/q0qlmtL6ZgPO\nxsYGPj4+0Ol0CA0NlcZ1Oh02bdpkch2NRoN9+/ZhaGhIug6n0+ng7u4ODw8PqaaqqspoPZ1Oh9Wr\nV8PKysrkdsfGxvDVV1/B29v7nvs71R/GTFR3oRlz586VZW5rK2vZ5r7DUnuXc25L7V3uvgH5ep83\nbx5UqiWyzD0VE94HFxsbi+LiYhQVFaG1tRWJiYkwGAyIiooCACQnJxuFX1hYGOzs7BATE4OLFy+i\nsrISGRkZiImJkWqioqLQ1dWFpKQktLa2oqioCCdOnJDeZQkABw8exNmzZ9HR0YGmpibExcWhpaUF\nL7/88nT2T0REgprwGpxWq0Vvby8OHz4MvV4PtVqNkpIS6R44vV6Pjo4Oqd7e3h7l5eVISEhAUFAQ\nHBwcEBcXh9jYWKlmyZIlKCkpwRtvvIGCggK4ubkhLS0NP//5z6Wa/v5+xMfHw2AwwN7eHt7e3jh9\n+jRWrVo1je0TEZGoJvUmk+joaERHR5tclp2dPW5MrVbj9OnTZrf5s5/9DOfOnbvn8pSUFKSkpExm\n94iIiMbhsyiJiEhIDDgiIhISA46IiITEgCMiIiEx4IiISEgMOCIiEhIDjoiIhMSAIyIiITHgiIhI\nSAw4IiISEgOOiIiExIAjIiIhMeCIiEhIDDgiIhISA46IiITEgCMiIiEx4IiISEgMOCIiEhIDjoiI\nhMSAm+EUUMi9C0REP0rWcu/ATPf/urpR+P7Hsszt+dPFeGSxsyxzExH92DHgpmh0dAzNbZ2yzG03\n25YBR0R0DzxFSUREQmLAERGRkBhwREQkJAYcEREJiQFHRERCYsAREZGQGHBERCSkSQVcXl4evLy8\noFQqERgYiNraWrP1zc3NCAkJgZubG9RqNdLS0sbV1NTU4KmnnoJSqYSPjw8KCwvH1VRUVMDPzw+u\nrq7w9/dHVVXVJNsiIiJLN2HAlZWVISkpCQkJCaiuroZGo0F4eDguX75ssr6/vx9arRZKpRI6nQ6p\nqak4evQoMjMzpZqOjg688MIL8Pf3R3V1NXbt2oXXX38dlZWVUk19fT2io6MRERGBmpoahIeHY+vW\nrfjiiy+moW0iIhLdhAGXlZWFLVu2IDIyEiqVCmlpaXB1dUVBQYHJ+tLSUgwODiInJweenp4IDQ1F\nfHw8srOzpZrCwkK4u7vj0KFDUKlUiIyMxObNm41CMCcnB08++SR27doFlUqF3bt3Y+3atcjJyZmG\ntomISHRmA254eBiNjY0ICgoyGg8ODkZdXZ3Jderr6xEQEABbW1uj+q6uLnR2dko1prb55ZdfYmRk\nBADQ0NBwX/MSERHdzWzA9fT0YGRkBC4uLkbjTk5OMBgMJtcxGAzj6p2dnaVlANDd3W2y5tatW+jp\n6TG7nXvNS0REdLdpf9iyQmFZH9+yZKEL3s9OknUf3s9+VLa5Ax6Xb24AWKdZKdvccvYuZ9+A5fZu\nyb/vM5HZIzhHR0dYWVmNO2rq7u6Gq6uryXVcXFxM1t9ZZq7G2toajo6OZmu+f1RHRERkitmAs7Gx\ngY+PD3Q6ndG4TqeDn5+fyXU0Gg1qa2sxNDRkVO/u7g4PDw+pxtQ2V69eDSsrK7M1/v7+k2yNiIgs\n2YTvooyNjUVxcTGKiorQ2tqKxMREGAwGREVFAQCSk5MRGhoq1YeFhcHOzg4xMTG4ePEiKisrkZGR\ngZiYGKkmKioKXV1dSEpKQmtrK4qKinDixAnExcVJNTt27MCnn36K9PR0tLW14ciRI6ipqcHOnTun\ns38iIhKUoq+vb2yiovz8fGRkZECv10OtViMlJQUBAQEAgJiYGJw/fx6NjY1S/V/+8hckJCTgwoUL\ncHBwQFRUFF5//XWjbZ4/fx5vvPEGWlpa4ObmhldffRVbt241qqmoqMCBAwfQ0dGBpUuX4q233sLz\nzz8/DW0TEZHoJhVwREREM82MeRbl+fPn8Ytf/AJqtRoODg4oLi4eV5OamopHH30Ubm5ueP7559HS\n0iLDnk6vI0eOICgoCB4eHli2bBl+8Ytf4OLFi+PqROw9NzcXP/vZz+Dh4QEPDw88++yz+Pjjj41q\nROz7+44cOQIHBwe89tprRuMi9p6amgoHBwejL09Pz3E1ovV9x5UrV7Bjxw4sW7YMSqUS/v7+OH/+\nvFGNaP2vXLly3Gvu4OCAiIgIAMDY2NgD9zxjAu5///d/sWLFCqSmpsLOzm7c7Qjp6enIzs5GWloa\nzp49C2dnZ2i1Wty4cUOmPZ4e58+fx/bt2/Hxxx+jsrIS1tbW2LRpE/r6+qQaUXtfuHAh9u/fj08/\n/RT/+Z//iSeffBJbtmzBn//8ZwDi9n23hoYG/OEPf8Bjjz1m9Dsvcu/Lly9HW1ub9PXZZ59Jy0Tu\nu6+vDxvfMutEAAAJGElEQVQ2bIBCoUBpaSnq6+uRlpYm3UcMiNn/uXPnjF7vc+fOQaFQQKvVAgAy\nMjIeuOcZeYpy0aJFePfdd7F582YAtxPe09MTv/rVr7Br1y4AwODgIFQqFd55551x1/ZmsoGBAXh4\neKC4uBgbNmywqN4B4JFHHsG+ffsQGRkpfN/Xrl1DYGAgjh49ioMHD0oPLhf5NU9NTcWpU6eMQu0O\nkfsGgP3796O2thZnzpwxuVz0/u84fPgwMjMz0draChsbmyn1PGOO4Mz55ptvYDAYEBwcLI3Nnj0b\nTzzxhHCP9rp+/TpGR0cxf/58AJbT+8jICE6ePImhoSE88cQTFtH3q6++ik2bNmHt2rUYG/u/v0NF\n772jowOPPvoovL29ER0djY6ODgDi9/3hhx9i9erViIqKgkqlwrp165CbmystF71/4HaI//GPf8QL\nL7wAW1vbKfc87U8ykYNerwcAo0N54PYjxa5cuSLHLv3d7NmzB15eXtBoNADE7725uRnPPvsshoaG\nYGdnh8LCQqhUKumXW9S+//CHP6CjowN5eXkAjJ8QJPJr7uvri5ycHKhUKnR3d+Pdd9/Fhg0b8Pnn\nnwvdN3A72PPz8xEbG4tdu3ahqakJiYmJAIDt27cL3z9w+17nzs5O/PM//zOAqf+uCxFw5oj06LA3\n3ngD9fX1OHPmzKT6EqH35cuX4/z587h27RoqKioQHR2NU6dOmV1npvfd3t6Od955Bx999JH04IOx\nsTGjo7h7mem9r1+/3uh7X19feHt7o7i4GGvWrLnnejO9bwAYHR3F448/jt/+9rcAbr/54tKlS8jL\ny8P27dvNritC/8DtP+wef/xxPPbYYxPWTqZnIU5R3nls2J1Hgt0h0qO9kpKSUF5ejsrKSixZskQa\nF733f/iHf8DDDz8Mb29v7N27F2vWrEFubq7QfdfX16Onpwf+/v5wcnKCk5MTPvvsM+Tn58PZ2Vl6\nnJ2IvX/fnDlz4Onpia+//lro1xwAlEol/vEf/9FoTKVSSZ+9KXr/3d3dOHPmDCIjI6WxqfYsRMAt\nWbIErq6uOHv2rDQ2ODiIzz///J6PFJtJEhMTpXBbtmyZ0TLRe/++kZERjI6O4uGHHxa27+effx61\ntbWoqalBTU0NqqursWrVKoSFhaG6uho//elPhe39+wYHB9HW1gZXV1ehX3MA8Pf3R1tbm9HYX//6\nV+kRh6L/Wy8uLsbs2bMRFhYmjU21Z6s9e/bs+3vs7HQbGBhAS0sL9Ho9/vjHP0KtVuOhhx7Cd999\nh5/85CcYGRnBv/zLv2DZsmUYGRnBm2++CYPBgPT0dNjY2Mi9+w8sISEB77//PgoLC7Fw4UIMDAxg\nYGAACoUCNjY2UCgUwva+b98+2NraYnR0FH/729+Qk5OD0tJS7N+/H4888oiwfc+ePVs6cnNycoKz\nszNKSkqwePFivPjii0K/5m+99Zb0mv/1r3/Fa6+9hq+//hrp6emwt7cXtm8AWLx4MQ4dOgQrKyso\nlUqcO3cOv/vd77Br1y6sXr1a6Nd9bGwMsbGxeO655/Dzn/9cGp9qzzPmGtyFCxewceNGALebTk1N\nRWpqKl588UVkZWUhPj4eN2/exGuvvYa+vj6sWbMGZWVlmDt3rsx7PjX5+flQKBRGz/sEbr/Z5M4F\naFF7NxgMeOWVV2AwGGBvb48VK1bg5MmT0gfhitq3KQqFwuiag6i9d3V1Ydu2bejp6YGTkxN8fX3x\nH//xH1i0aBEAcfsGgFWrVuH48ePYv38/3n33XSxevBhvvfUWoqOjpRpR+6+ursbXX38tvanqblPp\neUbeB0dERDQRIa7BERERfR8DjoiIhMSAIyIiITHgiIhISAw4IiISEgOOiIiExIAjIiIhMeCIiEhI\nDDgiIhISA46IiITEgCOSyaVLl/DKK6/A29sbbm5u8PHxwe7du9HX1zeuNjs7GytXroRSqcTTTz+N\nuro6rFy5EjExMUZ1HR0d2L59O5YtWwZXV1esW7cOVVVVP1RLRD8qM+Zhy0SiuXLlChYuXIgDBw5g\nwYIF6OjowJEjR/DnP/8ZH3/8sVRXVFSEN998E5GRkdi0aRMuXbqE7du3o7+/3+gBzJcvX8b69evh\n4uKC1NRUODk54eTJk4iMjMTx48fxT//0T3K0SSQbPmyZ6Efi1q1baGhoQEhICM6dOwcvLy+Mjo5i\n5cqVWLFiBd5//32p9tSpU4iMjJQ+TQMA4uLi8Kc//QkNDQ2YP3++VKvVavHtt9+iurr6B++JSE48\nRUkkk+HhYfz+97+Hr68v3Nzc4OzsjJCQEADAf//3fwMA/va3v+F//ud/xn1cUkhICKytjU/AfPLJ\nJ3jmmWfw0EMP4datW9JXcHAwvvrqK9y4ceOHaYzoR4KnKIlkkpycjNzcXCQmJkKj0eChhx7C5cuX\n8ctf/hKDg4MAAL1eDwBwdnY2WtfKygqOjo5GY93d3Thx4gROnDgxbi6FQoHe3l7Mmzfv79QN0Y8P\nA45IJmVlZdi8eTN2794tjfX39xvVuLq6ArgdXncbGRnBt99+azTm6OiIJ554AvHx8SbnUyqV07Hb\nRDMGA45IJjdv3hx3mvH48eNG3y9cuBALFy7EBx98gBdffFEar6qqwsjIiFHt008/jYaGBnh6emL2\n7Nl/vx0nmiEYcEQyWb9+PU6cOAG1Wo1HHnkEp06dQkNDg1HNrFmz8PrrryM+Ph6//vWvERoaio6O\nDqSnp8Pe3h6zZv3fZfQ33ngDTz/9NEJCQrB9+3YsXrwYfX19uHjxIr755htkZmb+0C0SyYoBRyST\ntLQ0jI2N4Z133gEAPPvss8jPz0dwcLBRXWRkJAYGBpCdnY2SkhKo1Wr827/9GzZv3gx7e3upbtGi\nRdDpdDh48CDeeecdfPvtt1iwYAHUajU2b978g/ZG9GPA2wSIZqAvv/wSwcHBOHbsGF544QW5d4fo\nR4kBR/Qj98033yA3NxcBAQF46KGH0NbWhiNHjsDW1ha1tbW83kZ0DzxFSfQjZ2dnh5aWFrz//vvo\n6+vD/PnzERQUhLfffpvhRmQGj+CIiEhIfJIJEREJiQFHRERCYsAREZGQGHBERCQkBhwREQmJAUdE\nREL6/+XmTmCIkurqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1061eec50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compare ages of those who rate satisfaction highly and those who don't\n",
    "\n",
    "couples['high_rating'] = couples['rel_rating'] == 'very satisfied'\n",
    "\n",
    "for condition in [True, False]:\n",
    "    mrr = couples.where('high_rating', condition).select(['age'])\n",
    "    mrr.hist(bins=np.arange(15,71,5),normed=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two distributions are different, but they are also similar in many ways. This raises the question of how the two distributions compare in the U.S. population.\n",
    "\n",
    "To answer this, we start by setting up null and alternative hypotheses.\n",
    "\n",
    "**Null hypothesis.** Among married and cohabitating people in the United States, the distribution of ages of those who give the highest satisfaction rating is the same as the distribution of ages who do not. The difference in the sample is due to chance.\n",
    "\n",
    "**Alternative hypothesis.** Among married and cohabitating people in the United States, the distribution of ages of those who give the highest satisfaction rating is different from the distribution of ages who do not.\n",
    "\n",
    "Let us see if we can answer the question by just comparing the medians of the two distributions in the sample. As our test statistic, we will measure how far apart the two medians are.\n",
    "\n",
    "Formally, the **test statistic** will be the absolute difference between the medians of the ages of the two groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(couples.where('high_rating', True)['age'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(couples.where('high_rating', False)['age'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The median age of those who gave the highest rating is one year more than the median of those who did not. The observed value of our test statistic is therefore 1.\n",
    "\n",
    "To see how this statistic could vary under the null hypothesis, we will perform a permutation test. The test will involve repeatedly generating permuted samples of ages of the two groups, and calculating the test statistic for each of those samples. So it is useful to define a function that computes the test statistic – the absolute differeence between the medians – for every generated sample.\n",
    "\n",
    "Like the function ``tvd``, the function ``abs_diff_in_medians`` takes as its arguments the table of data, the name of a column of conditions that come in two categories, and the name of the column of values to be compared. It returns the absolute difference between the median values under the two conditions. By applying this function to the original sample, we can see that it returns the observed value of our test statistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def abs_diff_in_medians(t, conditions, values):\n",
    "    \"\"\"Compute the difference in the median of values for conditions 1 & 2.\"\"\"\n",
    "    a = t.where(conditions, True)[values]\n",
    "    b = t.where(conditions, False)[values]\n",
    "    return abs(np.median(a)-np.median(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "abs_diff_in_medians(couples, 'high_rating', 'age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run the permutation test, we will define a function ``permutation_stat`` that is similar to ``permutation_ tvd`` except that it takes a fourth argument. This argument is the name of a function that computes the test statistic. The code is the same as for ``permutation_tvd`` apart from the number of repetitions and the replacement of ``tvd`` by the generic function ``stat``.\n",
    "\n",
    "Because ``permutation_test`` takes a function as an argument, it is called a *higher order function*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def permutation_stat(original, conditions, values, stat):\n",
    "    repetitions=1600\n",
    "    test_stats = []\n",
    "\n",
    "    for i in np.arange(repetitions):\n",
    "        shuffled = original.sample()\n",
    "        combined = Table([original[conditions], shuffled[values]], \n",
    "                         [conditions, values])\n",
    "        test_stats.append(stat(combined, conditions, values))\n",
    "\n",
    "    observation = stat(original, conditions, values)\n",
    "    p_value = np.count_nonzero(test_stats >= observation) / repetitions\n",
    "    \n",
    "    print(\"Observed value of test statistic:\", observation)\n",
    "    print(\"Empirical P value:\", p_value)\n",
    "    t = Table([test_stats], ['Empirical dist. of test statistic under null'])\n",
    "    t.hist(normed=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us run the permutation test using ``abs_diff_in_medians`` to compute the test statistic. The empirical histogram of the test statistic shows that our observed value of 1 is typical under the null hypothesis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observed value of test statistic: 1.0\n",
      "Empirical P value: 0.7025\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEqCAYAAABAysQTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcU/f9P/BXQMALouEu0mCVWKRgkSqVqrNImVhbxT68\nADqs2toKdOs6HOA2V6w1FtTWIlAv2IuKfrWKyOqtKls7SmWOVh1zmulEbF1ARSxh4iX9/eGD/IyE\nhIREPiav5+PBQ3PyOee8P/kQXjmXnCO5du3aTyAiIhKMQ1cXQEREpA8DioiIhMSAIiIiITGgiIhI\nSAwoIiISEgOKiIiExIAiIiIhGQ2oVatWISoqCjKZDIGBgYiPj8epU6cMzlNTUwOpVNrm58iRIxYr\nnIiIbFs3Yw3Ky8vxyiuvIDw8HBqNBsuWLUNcXByOHj2Kvn37Gpx3165dCAkJ0T421p6IiKiV0YDa\nuXOnzuO1a9dCJpPh6NGjGD9+vMF5pVIpvLy8OlchERHZJZOPQf3444/QaDQd2hqaNWsW5HI5YmNj\nUVJSYlaBRERknySmXovvpZdewn/+8x/8+c9/hkQi0dvm6tWr2Lp1K0aOHAlHR0fs3bsXK1euREFB\nAaZPn26RwomIyLaZFFCLFi3C7t27sW/fPgQEBJi0ooULF+Lrr79GeXm5yUUSEZH96fAuvszMTBQX\nF2PPnj0mhxMADBs2DOfOnTN5PiIisk9GT5IAgPT0dJSUlKC0tBSBgYFmrejkyZPw9fU1a14iIrI/\nRreg0tLSsHXrVqxbtw5ubm5QqVRQqVRQq9XaNllZWZg8ebL2cVFRET777DOcPn0aSqUSubm5KCws\nxPz5863Ti4eQUqns6hK6jL323V77DbDvZB6jW1CFhYWQSCQ6AQQAGRkZSE9PBwCoVCqcP39e+5xE\nIsGKFStQW1sLR0dHBAYGIi8vD9OmTbNs9UREZLOMBlRDQ4PRheTn5+s8TkhIQEJCgvlVERGR3eO1\n+IiISEgMKCIiEhIDioiIhMSAIiIiITGgiIhISAwoIiISEgOKiIiExIAiIiIhMaCIiEhIDCgiIhIS\nA4qIiITEgCIiIiExoIiISEgMKCIiEhIDioiIhMSAIiIiITGgiIhISAwoIiISEgOKiIiExIAiIiIh\nMaCIiEhIDCgiIhISA4qIiITEgCIiIiExoIiISEgMKCIiEhIDioiIhMSAIiIiIXXr6gLo4VZ35Rrq\nrzSaNE9T001Un6mxUkV3eXn0gbdHX6uug4isiwFFnVJ/pRFL3i8yaR61Wo1evXpZqaK7Fr+RyIAi\neshxFx8REQmJAUVEREJiQBERkZAYUEREJCQGFBERCYkBRUREQmJAERGRkAwG1KpVqxAVFQWZTIbA\nwEDEx8fj1KlTRhdaXV2N5557Dv369UNwcDCys7MtVjAREdkHgwFVXl6OV155BQcPHsSePXvQrVs3\nxMXF4dq1a+3Oc/36dUyZMgW+vr4oKyuDQqFAbm4u1qxZY/HiiYjIdhm8ksTOnTt1Hq9duxYymQxH\njx7F+PHj9c6zY8cO3LhxAwUFBXBxcUFQUBCUSiXy8/ORmppqucqJiMimmXQM6scff4RGo0Hfvu1f\nQqayshKRkZFwcXHRThs3bhwuXbqECxcumF8pERHZFZMCKiMjA0OHDkVERES7berq6uDt7a0zzcvL\nS/scERFRR3T4YrGLFi1CZWUl9u3bB4lE0m47Q88RERF1VIcCKjMzE7t370ZpaSkCAgIMtvX29m6z\npVRfX699rj1KpbIjpdgUW+hzU9NNqNVqk+czZx5TNDU1Cfn6iljTg8K+2w+5XG6R5RgNqPT0dJSU\nlKC0tBSBgYFGFxgREYG33noLLS0t2uNQZWVl8PPzg0wma3c+S3XoYaFUKm2iz9Vnaky+dcaDuN2G\nq6sr5HLDH6YeNFsZc3Ow7/bZ984yeAwqLS0NW7duxbp16+Dm5gaVSgWVSqXz6TcrKwuTJ0/WPp46\ndSp69OiB5ORknDp1Cnv27MHq1auRnJxsvV4QEZHNMbgFVVhYCIlEohNAwN2TJdLT0wEAKpUK58+f\n1z7n5uaG4uJipKWlISoqClKpFKmpqUhJSbF89UREZLMMBlRDQ4PRBeTn57eZFhwcjL1795pfFRER\n2T1ei4+IiITEgCIiIiExoIiISEgMKCIiEhIDioiIhMSAIiIiITGgiIhISAwoIiISEgOKiIiExIAi\nIiIhMaCIiEhIDCgiIhISA4qIiITEgCIiIiExoIiISEgMKCIiEhIDioiIhMSAIiIiITGgiIhISAwo\nIiISEgOKiIiExIAiIiIhMaCIiEhIDCgiIhISA4qIiITEgCIiIiExoIiISEgMKCIiEhIDioiIhMSA\nIiIiITGgiIhISAwoIiISEgOKiIiExIAiIiIhMaCIiEhIDCgiIhISA4qIiIRkNKDKy8sRHx+P4OBg\nSKVSFBUVGWxfU1MDqVTa5ufIkSMWK5qIiGxfN2MNmpubERISgoSEBCxYsAASiaRDC961axdCQkK0\nj/v27Wt+lUREZHeMBlRMTAxiYmIAACkpKR1esFQqhZeXl/mVERGRXbPaMahZs2ZBLpcjNjYWJSUl\n1loNERHZKKNbUKbq3bs3li5dipEjR8LR0RF79+7F3LlzUVBQgOnTp1t6dUREZKMsHlDu7u46uwLD\nwsLQ0NCA1atXGwwopVJp6VKEZwt9bmq6CbVabfJ85sxjiqamJiFfXxFrelDYd/shl8stshyLB5Q+\nw4YNw+bNmw22sVSHHhZKpdIm+lx9pga9evUyaR61Wm3yPKZydXWFXB5g1XWYylbG3Bzsu332vbMe\nyPegTp48CV9f3wexKiIishFGt6DUajXOnj0LANBoNKitrcWJEyfg7u4Of39/ZGVloaqqSnsiRFFR\nEZydnREaGgoHBwfs378fhYWFyMrKsm5PiIjIphgNqKqqKkyaNAkAIJFIoFAooFAokJiYiLy8PKhU\nKpw/f17bXiKRYMWKFaitrYWjoyMCAwORl5eHadOmWa0TRERke4wG1JgxY9DQ0NDu8/n5+TqPExIS\nkJCQ0PnKiIjIrvFafEREJCQGFBERCYkBRUREQmJAERGRkBhQREQkJAYUEREJiQFFRERCYkAREZGQ\nGFBERCQkBhQREQmJAUVEREJiQBERkZAYUEREJCQGFBERCYkBRUREQmJAERGRkBhQREQkJAYUEREJ\niQFFRERCYkAREZGQGFBERCQkBhQREQmJAUVEREJiQBERkZAYUEREJCQGFBERCYkBRUREQmJAERGR\nkBhQREQkJAYUEREJiQFFRERCYkAREZGQGFBERCQkBhQREQmJAUVEREJiQBERkZAYUEREJCSjAVVe\nXo74+HgEBwdDKpWiqKjI6EKrq6vx3HPPoV+/fggODkZ2drZFiiUiIvthNKCam5sREhIChUKBHj16\nQCKRGGx//fp1TJkyBb6+vigrK4NCoUBubi7WrFljsaKJiMj2dTPWICYmBjExMQCAlJQUowvcsWMH\nbty4gYKCAri4uCAoKAhKpRL5+flITU3tfMVERGQXLH4MqrKyEpGRkXBxcdFOGzduHC5duoQLFy5Y\nenVERGSjLB5QdXV18Pb21pnm5eWlfY6IiKgjjO7iM5WxY1TtUSqVFq5EfLbQ56amm1Cr1SbPZ848\npmhqahLy9RWxpgeFfbcfcrncIsuxeEB5e3u32VKqr6/XPtceS3XoYaFUKm2iz9VnatCrVy+T5lGr\n1SbPYypXV1fI5QFWXYepbGXMzcG+22ffO8viu/giIiJQUVGBlpYW7bSysjL4+flBJpNZenVERGSj\njAaUWq3GiRMncOLECWg0GtTW1uLEiRO4ePEiACArKwuTJ0/Wtp86dSp69OiB5ORknDp1Cnv27MHq\n1auRnJxsvV4QEZHNMRpQVVVVGDt2LMaOHYsbN25AoVBg7NixUCgUAACVSoXz589r27u5uaG4uBiX\nLl1CVFQU0tPTkZqa2qFT1ImIiFoZPQY1ZswYNDQ0tPt8fn5+m2nBwcHYu3dv5yojIiK7xmvxERGR\nkBhQREQkJAYUEREJiQFFRERCYkAREZGQGFBERCQkBhQREQmJAUVEREJiQBERkZAYUEREJCQGFBER\nCYkBRUREQmJAERGRkBhQREQkJAYUEREJyej9oB6Uwm0HuroEAECP7s6YGB2BPr17dXUpRER2TZiA\nOvhlVVeXAADo07snJkSN6OoyiIjsHnfxERGRkBhQREQkJAYUEREJiQFFRERCYkAREZGQGFBERCQk\nBhQREQmJAUVEREJiQBERkZAYUEREJCQGFBERCYkBRUREQmJAERGRkBhQREQkJAYUEREJiQFFRERC\nYkAREZGQGFBERCQkBhQREQmpQwG1YcMGDB06FL6+vnjmmWdQUVHRbtuamhpIpdI2P0eOHLFY0URE\nZPu6GWuwa9cuZGZmYuXKlYiMjMT69esxbdo0fPPNN/D39zc4X0hIiPZx3759LVMxERHZBaNbUHl5\neZg5cyaSkpIgl8uRnZ0NHx8fbNy40eB8UqkUXl5e2h8nJyeLFU1ERLbPYEDdvHkTx48fR1RUlM70\ncePG4ejRowYXPGvWLMjlcsTGxqKkpKTzlRIRkV0xuIvvypUruHPnDry9vXWme3p6oq6uTu88vXv3\nxtKlSzFy5Eg4Ojpi7969mDt3LgoKCjB9+nTLVU5ERDbN6DEoU7m7uyMlJUX7OCwsDA0NDVi9ejUD\nioiIOsxgQHl4eMDR0bHN1lJ9fT18fHw6vJJhw4Zh8+bNBtuo1eoOL8+anByBK5cv43LdJauvS6lU\nWn0d1tbUdNOssbP2eDc1NQn5+opY04PCvtsPuVxukeUYDChnZ2eEhYWhrKwMkydP1k4vKytDXFxc\nh1dy8uRJ+Pr6GmzTq1evDi/Pmnr27AkPT09I+7hadT1KpdJig9iVqs/UmDx2arXa6uPt6uoKuTzA\nquswla2MuTnYd/vse2cZ3cWXkpKCV199FeHh4XjqqaewceNG1NXVYc6cOQCArKwsVFVVaU+EKCoq\ngrOzM0JDQ+Hg4ID9+/ejsLAQWVlZ1u0JERHZFKMBNWXKFFy9ehUrVqyASqVCcHAwtm/frv0OlEql\nwvnz57XtJRIJVqxYgdraWjg6OiIwMBB5eXmYNm2a1TpBRES2p0MnScybNw/z5s3T+1x+fr7O44SE\nBCQkJHS+MiIismu8Fh8REQmJAUVEREJiQBERkZAYUEREJCQGFBERCYkBRUREQmJAERGRkBhQREQk\nJAYUEREJiQFFRERCYkAREZGQGFBERCQkBhQREQmJAUVEREJiQBERkZAYUEREJCQGFBERCYkBRURE\nQmJAERGRkBhQREQkJAYUEREJiQFFRERCYkAREZGQunV1ASJqUv8PP6iuWHcdTTdRfabGpHm8PPrA\n26OvlSoiIhILA0qPKw3XocjbbtV1qNVq9OrVy6R5Fr+RyIAiIrvBXXxERCQkBhQREQmJAUVERELi\nMSgiO1d35RrqrzRadR08KYjMwYAisnP1Vxqx5P0iq66DJwWRObiLj4iIhMSAIiIiITGgiIhISAwo\nIiISEgOKiIiExIAiIiIhMaCIiEhIHQqoDRs2YOjQofD19cUzzzyDiooKg+2rq6vx3HPPoV+/fggO\nDkZ2drZFiiUiIvthNKB27dqFzMxMpKWl4auvvkJERASmTZuGixcv6m1//fp1TJkyBb6+vigrK4NC\noUBubi7WrFlj8eKJiMh2GQ2ovLw8zJw5E0lJSZDL5cjOzoaPjw82btyot/2OHTtw48YNFBQUICgo\nCJMnT8avfvUr5OfnW7x4IiKyXQYvdXTz5k0cP34cv/zlL3Wmjxs3DkePHtU7T2VlJSIjI+Hi4qLT\n/p133sGFCxcgk8ksUDbRQ8jB2eTr0T0ILTdvdXUJRHoZDKgrV67gzp078Pb21pnu6emJuro6vfPU\n1dXB399fZ5qXl5f2OQYU2aurjU1YtWFPV5fRxhvz4rq6BCK9LH6xWIlEYtZ8/5efaeFKOke0ekT1\n+OAAvlYd9FT44/i//Me7ugy9Ip8c0tUl2Cy5XN7VJTy0DB6D8vDwgKOjY5utpfr6evj4+Oidx9vb\nW2/71ueIiIg6wmBAOTs7IywsDGVlZTrTy8rK8NRTT+mdJyIiAhUVFWhpadFp7+fnx917RETUYUbP\n4ktJSUFRURE+/fRTnD59Gunp6airq8OcOXMAAFlZWZg8ebK2/dSpU9GjRw8kJyfj1KlT2LNnD1av\nXo3k5GTr9YKIiGyO0WNQU6ZMwdWrV7FixQqoVCoEBwdj+/bt2hMhVCoVzp8/r23v5uaG4uJipKWl\nISoqClKpFKmpqUhJSbFaJ4iIyPZIrl279lNXF0FERHQ/q16Lz54vkWRK32tqaiCVStv8HDly5AFW\n3Hnl5eWIj49HcHAwpFIpioqM30bcVsbc1L7bypivWrUKUVFRkMlkCAwMRHx8PE6dOmV0PlsYd3P6\nbgvjvn79eowaNQoymQwymQw///nPcfDgQYPzmDveFj/NvFXrJZJWrlyJyMhIrF+/HtOmTcM333zT\n5ntSwP+/RNLo0aNRVlaG06dPIzU1FT179kRqaqq1yrQKU/t+73whISHax3379n0Q5VpMc3MzQkJC\nkJCQgAULFhj9yoEtjbmpfW/1sI95eXk5XnnlFYSHh0Oj0WDZsmWIi4vD0aNH2+2LrYy7OX1v9TCP\ne//+/bFkyRIMGjQIGo0GRUVFmDlzJo4cOYLQ0NA27Tsz3lbbxRcdHY3Q0FC8//772mlPPvkkJk+e\njMWLF7dpX1hYiKysLCiVSu1VKFasWIGNGzfin//8pzVKtBpT+15TU6M9WzIsLOxBlmo1/v7+yMnJ\nQUJCQrttbGnM79WRvtvimAOAWq2GTCZDUVERxo8fr7eNrY57R/puq+P+6KOP4q233sLs2bPbPNeZ\n8bbKLr7WSyRFRUXpTDfnEkmXLl3ChQsXrFGmVZjT91azZs2CXC5HbGwsSkpKrFmmEGxlzDvD1sb8\nxx9/hEajMbhFYKvj3pG+t7KVcb9z5w527tyJlpYWPP3003rbdGa8rRJQ5l4i6f72914i6WFhTt97\n9+6NpUuX4pNPPsGOHTvws5/9DHPnzsX27dsfRMldxlbG3By2OuYZGRkYOnQoIiIi2m1jq+Pekb7b\nyrhXV1ejf//+8PHxwRtvvIGPPvqo3StmdGa8rXYMylTmXiLJFri7u+uchh8WFoaGhgasXr0a06dP\n78LKrItjbltjvmjRIlRWVmLfvn0Gx9YWx72jfbeVcR88eDDKy8vR2NiIkpISzJs3D6WlpRg2bFib\ntp0Zb6tsQdnzJZLM6bs+w4YNw7lz5yxdnlBsZcwt5WEe88zMTBQXF2PPnj0ICAgw2NbWxt2Uvuvz\nMI67k5MTBgwYgCeeeAKLFy/G8OHDsX79er1tOzPeVgkoe75Ekjl91+fkyZPw9fW1dHlCsZUxt5SH\ndczT09O1f6ADAwONtrelcTe17/o8rON+rzt37kCj0eh9rjPjbbXvQdnzJZJM7XtRURE+++wznD59\nGkqlErm5uSgsLMT8+fO7qgtmUavVOHHiBE6cOAGNRoPa2lqcOHFCe/dlWx5zU/tuK2OelpaGrVu3\nYt26dXBzc4NKpYJKpYJarda2sdVxN6fvtjDub731FioqKlBTU4Pq6mpkZWWhvLwcM2bMAGDZ8bba\nMSh7vkSSqX2XSCRYsWIFamtr4ejoiMDAQOTl5WHatGld1APzVFVVYdKkSQDu9kmhUEChUCAxMRF5\neXk2Peam9t1WxrywsBASiUTnDxJw94SB9PR0ALb7Xjen77Yw7nV1dZg/fz7q6urg5uaGkJAQ7Ny5\nU3vmsiXHm5c6IiIiIVn1UkdERETmYkAREZGQGFBERCQkBhQREQmJAUVEREJiQBERkZAYUEREJCS7\nCagtW7bovZOlVCrFgAEDHng9CoUCUqm0Q22/+uorSKVSlJeXW7yO1jt8duTut/pMnDgRzz//vPax\nObX+6U9/Ql5enlnrN1VlZSWio6PRv39/SKVS/OMf/+iSmmpqaqBQKHS+0NgZna33q6++wvLly/HT\nT7pfi2z9/di6dWuXLKurmfI+FY0l3ptdzW4CqtUnn3yCQ4cO6fzs3r37gdcxe/ZsHDp0qENtw8LC\ncOjQIQwdOtRq9Zh7xWGJRKIzrzm1fv7558jPzzdr/aZ6/fXXodFosG3bNhw6dAiDBg3qkpouXLiA\n7Oxs1NTUWGR5na33r3/9K9599902odKvXz8cOnSo3RvwWXtZInhYr75+/3vzYSTM7TYelKFDh3bJ\nFtP9/Pz84OfnZ7DNnTt3ANy9h8yTTz75IMoy2U8//aTzJhC5Vo1Gg3//+99IS0vDmDFjurocId0f\nKs7OzmaPpyWX1ZXu74e5bt26BScnJ4ssqyPuf28+jOxuC8qY1l2BR48eRVJSEh555BEMHjwY7733\nHgBg//79GDVqFPz8/DBu3Dh89913OvNPnDgREyZMwOeff47IyEj4+PggIiKizVaavl0HUqkUS5cu\nxXvvvYehQ4fC29sb//znP9vdNC8tLcX48ePh7+8PmUyG6Oho7Nu3T/v8unXrEBMTg0cffRQBAQGI\niYnBwYMHzX5tdu7ciREjRsDHxweRkZEoLS1t00ZfrYcPH8bPf/5zyGQy+Pv7Y8SIEcjOzgYALFiw\nANu2bcMPP/yg3eX6xBNPmFzb9evXsXDhQgQFBcHHxwcjRozQ2aLYsmULPDw8oNFokJ2dDalU2u5W\nnrGaLl++jF//+tcIDg7Wju8nn3yiswyVSoXXXnsNQ4YMgY+PD4KCgjBjxgxcvnwZX331lfa6fXFx\ncdp1GNr1Yupr2Nq3lpYWZGZm4umnn4a/vz8ee+wxxMfHQ6lUapetUCi0y/L09NQuA9C/C7iqqgpx\ncXEYOHAg+vXrh7CwMKSlpZm1LODuFldcXBxkMhn69++P0aNHY9OmTe2+FgAQGhqq92KjUqkUy5cv\n1+mbVCrFuXPnMH36dPj7+yM0NBTZ2dltguf48eOYMGECfH19ERwcjJycHL3hdPv2baxatUr7Xhgy\nZAh+//vf61ytu7WvhYWFWLx4sfb3srGxUW9/Wt83+/btw8KFCzFo0CAMGjQI8+fP15mnvdfwYdx9\n1xF2twV1+/Zt3L59W2eag4MDHBx0szo5ORkJCQmYN28edu/ejSVLlqCurg5/+ctfsHDhQvTs2RN/\n/OMfMWvWLHz77bfaT0YSiQTnzp1DRkYGMjMz4eXlhQ0bNmDu3Lnw8PDQ+eSu79NNUVERBgwYgHfe\neQe9evWCr68vrl271qbd2rVrkZGRgeeffx6pqalwdXXFd999h9raWm2bCxcuYNasWRgwYAA0Gg32\n7t2LGTNm4LPPPkN0dLRJr9uf//xnvPzyy4iNjcWyZctQX1+PzMxM3L59u907aQLA+fPnkZCQgLi4\nOGRkZMDJyQlnz57V7tpKT0/H1atXUVVVhW3btgG4+0nbFBqNBjNmzMCJEyewaNEiPP7449i/fz9+\n97vf4cqVK/jDH/6A2NhY7N+/H7GxsUhKSkJSUlK76zFU0/Xr1xEbG4uWlhZkZGQgICAAhw8fxptv\nvomWlhbtValfffVVfP/993j77bfRv39/1NXV4csvv0RzczPCwsKwYsUKpKWlITs7G+Hh4QDu3gTO\n0q9hS0sLmpqa8Oabb6Jfv364du0aNmzYgJiYGFRWVsLb2xuzZ8/GpUuXsGnTJhw4cACOjo5tamj9\nXW1qasKLL76IESNGoKCgAK6urqipqcHf/vY3ADBpWcDdXZNJSUmIjIzE6tWr4e7ujlOnTmmvAt8e\nQ7uv9E2fNWsWZs6ciZSUFOzbtw8KhQL9+/fHzJkzAdy9E/akSZPQr18/fPjhh3ByckJubi5qa2vb\nLG/+/Pk4cOAA3njjDUREROD06dN45513cOHCBXz66ac6bVeuXInw8HB88MEHuHPnjs5tz/XJyMhA\nbGwsCgsLcebMGfzxj3+Eo6MjCgoKjPbRFtldQI0YMaLNtPHjx2vf2K3i4+O1nwpHjx6N0tJSbNiw\nAX//+9+19zDRaDRITExEZWUlRo0aBeDuZnVdXR0OHTqk3Z0RHR2NkSNHYtmyZTpbOO3tOiguLtb5\nRf7Xv/6l8/z169fx9ttv44UXXtB5Q7ReTbjV0qVLtf/XaDQYM2YMzp49i8LCQpMDSqFQICgoSOcA\n9+DBgxETE2MwoI4fP45bt25h1apVcHV1BQCdkB4wYADc3d07tfvn4MGD+Oabb5Cfn4+EhAQAwDPP\nPIPm5masWbMGKSkp8PDwQJ8+fQDc3b1qaF2Gavrwww9x8eJFVFRU4NFHHwUAjB07Fo2NjXj33Xfx\n8ssvw8HBAceOHcPixYsxdepU7bz3XvW6NYwGDx5stN+deQ3d3NyQm5urfazRaBAVFYXHHnsMn332\nGZKTk+Hn54d+/foBAIYPH97mw9q9lEolGhsbkZWVheDgYADAqFGjkJiYCAAmLeunn35CRkYGnnji\nCfzpT3/STh87dqzB18Mcqamp2hrHjh2LL7/8Ejt37tQGVH5+Pm7cuIFdu3Zpd71HRUUhJCREZzlf\nf/01iouLsXbtWu0dcMeOHQupVIr58+fj5MmTCA0N1bb39vbG5s2bO1znqFGj8O677wK4+zusVCqx\nadOmNgFlL+xuF9+WLVtQVlam83PvLoFWMTEx2v87Ojpi4MCBkMvlOjfYav3D/MMPP+jM6+/vr/OH\nwsHBAZMmTUJVVZXR+qKjo41+yqqsrIRarcZLL71ksN13332HGTNmYPDgwfD09ISXlxfKyspw9uxZ\no3Xc686dO/j222+1u6VaDR8+3OgNx4YOHQonJyfMnTsXJSUl2jtpWtLXX38NBweHNrcsmDZtGm7e\nvKn9dG8Jhw8f1va7dWv89u3bGDduHK5evar9MDFs2DB88MEH+PDDD1FdXd2p4xidfQ2Li4sRHR2N\ngIAAeHh4oH///mhqajL59wAABg4ciD59+uCNN97A9u3bjW7pGKJUKnHx4kUkJSWZvYyOuv/EjCFD\nhujUXllDksJSAAAIOUlEQVRZieHDh+scF+7ZsydiY2N1xu7w4cNwdnbGCy+8oDP+rR8Ov/76a531\nTJw4sVN1BgcHo6WlxSrvm4eB3QVUcHAwwsLCdH70nTTRt29fncdOTk7aT+D3TgOAGzdu6EzXdxtj\nb29v3Lx5E5cvXzZYX0duC3/16lUAMHiSxcWLFzFp0iQ0NjYiJycHX3zxBcrKyvDss8+2qdeYK1eu\n4NatW3r75eXlZXDeRx99FDt37oRGo8Frr72Gxx57DDExMRbdV97Q0ACpVIpu3XR3CLS+lg0NDRZb\nV319PcrLy7WB3/rz0ksvQSKRaMfmo48+woQJE/DBBx9g9OjRCA4O1nvcoyM68xru27cPc+fORVBQ\nEAoLC3HkyBGUlZXB09PT5N8DAOjTpw9KS0vh6+uLtLQ0hIaG4umnn8aePXtMXlZHfo8t5f7jvc7O\nzjr9V6lU7b5v71VfX4+bN2/Cz89PZ/zlcjkkEkmb37WOvJ/vdf/fndZdteaMlS2wu118D4JKpWoz\nra6uDs7OzvD09DQ4b0f2LXt4eAC4u+UWFBSkt83hw4fx448/4qOPPtLucgGgc7fPjvLw8ICTkxPq\n6uraPFdXV4eAgACD848ZMwZjxozBrVu3UFFRAYVCgRkzZuDkyZMW+Y6JVCpFQ0MDbt++rRNSreNg\nye+xeHh4wMfHBwqFQu/zrbf99vT0RE5ODnJycnD27FkUFRVBoVDA09MTc+fONXm95r6Gu3btwqBB\ng3S+I3Xr1i1tOJgjNDQUn376KTQaDaqqqrBq1SrMmTMHf/3rXzFkyJAOL+fe32NTde/eHbdu3dKZ\n1pk++fr6tvu+vZe7uzu6d++us6v+/uXcy9LHirp37w4AFu27yOxuC+pB+P7773Hs2DHt4zt37qCk\npMRip9hGRETA1dUVH3/8cbttmpubAUDnD/a///1vHD161OT1OTo6Ijw8HCUlJTpbAMeOHdM5KcMY\nJycn/OxnP8Prr78OtVqtPcjv4uKC//3vfybX1Wr06NHQaDQoLi7Wmb5jxw64uLggIiLC5GW2V1N0\ndDROnz4Nf3//NlviYWFh2mNE9xo0aBD+8Ic/oG/fvjh16pR2+YDpn4xNfQ2bm5vbnKiwbds2aDSa\nNv1tbd9RDg4OGD58OBYtWgSNRoMzZ86YtKzAwEDIZLI2JxZ0xCOPPILq6mqdaQcOHDB5Oa0iIiJw\n7NgxfP/999pparUa+/fv1wmZ1j0QjY2Nesf//oCyNG9vb7i4uLTpe2fOzhWZ3W1BHT9+XO/+3PDw\ncL1nHJnD29sbc+bMQWZmJjw8PLBx40acO3dOe6p6Z/Xu3RuLFy/Gb3/7WyQlJWHq1Kno3bs3Tp48\nie7du2P+/PmIiopCt27d8NprryElJQX//e9/sXz5cjzyyCNt/jh1RGZmJqZMmYLExES89NJLuHz5\nMpYvXw4fHx+Du602btyIiooKxMTEwM/PD1euXMF7770HPz8/7aftoKAgfPLJJ9i4cSPCwsLg4uKC\nxx9/HMDdT9mJiYk6B/rvFxMTg5EjR+LNN9/E5cuXERQUhIMHD2LTpk34zW9+Y9YWVHs1JScno7i4\nGBMmTEBycjIGDRqE5uZmKJVKVFRUoKioCI2NjYiLi8P06dMhl8vh5OSEzz//HNeuXcO4ceMA3P3j\n3K1bN2zatAl9+vSBi4sL5HK53oDrzGsYExODvXv3YtGiRRg/fjy+/fZbrF+/Hn369NEZt9Yt8TVr\n1uDZZ5+Fo6Mjhg0b1qaW/fv34+OPP8bzzz8PmUyG5uZmrF27Fr1799Z+EOjosiQSCRQKBX7xi1/g\nhRdewNy5c+Hu7o4zZ87g8uXLyMzMbHd8XnzxRaSmpmr7dfLkyU5doSI5ORkbNmzAiy++qD1TMjc3\nFz179tTZbTd69GhMnToVs2fPRkpKCsLDw+Hg4IALFy7giy++QFZWVrtf/rYEiUSCKVOmYPPmzQgM\nDERgYCAOHDjQ7u5eS32Hq6vYTUC1fgrSd2KBRCLB2bNntX/I9G2Wm7KpPnDgQPzqV7/CkiVLcPbs\nWQQEBKCwsBCjR4/WWZ4py7y/7SuvvAJvb2/k5ubi1VdfhZOTEx577DEsXLgQwN0/EuvXr8eyZcuQ\nmJiIgQMHIisrC1988YVZx3/Gjh2L9evXY/ny5UhKSsKgQYOwfPlyFBQUtKnt3sehoaE4dOgQlixZ\ngvr6ekilUkRGRqKwsFD7STspKQnHjh3DkiVL0NjYCJlMhuPHjwO4e9aZsUCVSCTYvn07lixZgtWr\nV+Pq1asICAjAsmXLsGDBApP7aqgmNzc3HDhwANnZ2Xj//fdx6dIl9OnTB3K5XHsSSY8ePRAWFoZP\nP/0UtbW1cHBwgFwux4YNGzBhwgQAd3cV5eTk4P3338fzzz8PjUaD0tJS7dmg9+rMazh79mxcvHgR\nW7Zswccff4zw8HBs3boVs2bN0hmn2NhYvPzyyygsLNR+j0nfbqPAwED07NkTOTk5UKlUcHV1xZNP\nPondu3drdyV3dFkA8Nxzz6G4uBg5OTl4/fXXAdw95mZs3BITE/H9999j06ZN+Pjjj/H0009jy5Yt\nbYKwvffZ/dPd3d1RUlKCjIwMLFiwAB4eHpgzZw5u3bqFnJwcnXnXrVuHtWvXYvPmzVi5ciWcnZ0h\nk8nw7LPPGj0ma0hHT5tfvnw5NBqN9t8XX3wR2dnZiI+PN9hHQ+sQleTatWsPd8QKZuLEidBoNO3u\noyYioo7hMSgiIhISA8rCHrZNaCIiUXEXHxERCYlbUEREJCQGFBERCYkBRUREQmJAERGRkBhQREQk\nJAYUEREJ6f8BljfhkluTCKEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1084f56d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "permutation_stat(couples, 'high_rating', 'age', abs_diff_in_medians)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The absolute difference between the medians does not reveal a difference between the two underlying distributions of age. As you can see from the empirical histogram above, the test statistic is rather coarse – it has a small number of likely values. This suggests that it might be more informative to compare means instead of medians.\n",
    "\n",
    "The functions that we have developed make it easy for us to run the permutation test again, this time with a different test statistic.\n",
    "\n",
    "The **new test statistic** is the absolute difference between the means.\n",
    "\n",
    "As before, we will define a function that computes the test statistic based on a sample. Applying this function to the orignial sample shows that the observed means of the two groups differ by about 1.07 years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def abs_diff_in_means(t, conditions, values):\n",
    "    \"\"\"Compute the difference in the mean of values for conditions 1 & 2.\"\"\"\n",
    "    a = t.where(conditions, True)[values]\n",
    "    b = t.where(conditions, False)[values]\n",
    "    return abs(np.mean(a)-np.mean(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0743190048811186"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "abs_diff_in_means(couples, 'high_rating', 'age')"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}