DS_Lecture3.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis Tests\n",
    "\n",
    "**Overview**\n",
    "\n",
    "1. Hypothesis Tests\n",
    "    * 1-sample\n",
    "    * 2-sample\n",
    "2. Permutation Tests"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WarmUp/Review Exercises\n",
    "\n",
    "Load the *Auto.csv* data (into a dataframe *cars*) for the following tasks and create a scatter plot `mpg ~ weight`.\n",
    "Take a look at row with index $25$ (i.e. row 26); we will refer to that data point as $x_{25}$ from now on.\n",
    "\n",
    "* Compute the standard score for \n",
    "    * $x_{25}$[\"mpg\"]\n",
    "    * $x_{25}$[\"weight\"]\n",
    "* Compute the product of these standard scores (call it $p_{25}$).\n",
    "* If you repeated this process for all rows of the *cars* dataframe and averaged all products $p_{i}$, what would the resulting number tell you ? (What is it called?)\n",
    "* Take a bootstrap sample from cars.mpg and compute the mean"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------\n",
    "\n",
    "**Importing Standard Libraries**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "#pd.options.mode.chained_assignment = None # disable chained assignment warning\n",
    "import seaborn as sns"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "\n",
    "**Loading our Functions**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'%.3f'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%run ./ourFunctions.py\n",
    "%precision 3"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------\n",
    "\n",
    "**Read in the Cars DF**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(392, 10)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#cars = pd.read_csv('../data/Auto.csv')\n",
    "cars = pd.read_csv(\"https://raw.githubusercontent.com/markusloecher/DataScience-HWR/main/data/Auto.csv\")\n",
    "#d\n",
    "cars.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mpg</th>\n",
       "      <th>cylinders</th>\n",
       "      <th>displacement</th>\n",
       "      <th>horsepower</th>\n",
       "      <th>weight</th>\n",
       "      <th>acceleration</th>\n",
       "      <th>year</th>\n",
       "      <th>origin</th>\n",
       "      <th>name</th>\n",
       "      <th>Manufacturer</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>18.0</td>\n",
       "      <td>8</td>\n",
       "      <td>307.0</td>\n",
       "      <td>130</td>\n",
       "      <td>3504</td>\n",
       "      <td>12.0</td>\n",
       "      <td>70</td>\n",
       "      <td>1</td>\n",
       "      <td>chevrolet chevelle malibu</td>\n",
       "      <td>chevrolet</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>15.0</td>\n",
       "      <td>8</td>\n",
       "      <td>350.0</td>\n",
       "      <td>165</td>\n",
       "      <td>3693</td>\n",
       "      <td>11.5</td>\n",
       "      <td>70</td>\n",
       "      <td>1</td>\n",
       "      <td>buick skylark 320</td>\n",
       "      <td>buick</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>18.0</td>\n",
       "      <td>8</td>\n",
       "      <td>318.0</td>\n",
       "      <td>150</td>\n",
       "      <td>3436</td>\n",
       "      <td>11.0</td>\n",
       "      <td>70</td>\n",
       "      <td>1</td>\n",
       "      <td>plymouth satellite</td>\n",
       "      <td>plymouth</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>16.0</td>\n",
       "      <td>8</td>\n",
       "      <td>304.0</td>\n",
       "      <td>150</td>\n",
       "      <td>3433</td>\n",
       "      <td>12.0</td>\n",
       "      <td>70</td>\n",
       "      <td>1</td>\n",
       "      <td>amc rebel sst</td>\n",
       "      <td>amc</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>17.0</td>\n",
       "      <td>8</td>\n",
       "      <td>302.0</td>\n",
       "      <td>140</td>\n",
       "      <td>3449</td>\n",
       "      <td>10.5</td>\n",
       "      <td>70</td>\n",
       "      <td>1</td>\n",
       "      <td>ford torino</td>\n",
       "      <td>ford</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    mpg  cylinders  displacement  horsepower  weight  acceleration  year  \\\n",
       "0  18.0          8         307.0         130    3504          12.0    70   \n",
       "1  15.0          8         350.0         165    3693          11.5    70   \n",
       "2  18.0          8         318.0         150    3436          11.0    70   \n",
       "3  16.0          8         304.0         150    3433          12.0    70   \n",
       "4  17.0          8         302.0         140    3449          10.5    70   \n",
       "\n",
       "   origin                       name Manufacturer  \n",
       "0       1  chevrolet chevelle malibu    chevrolet  \n",
       "1       1          buick skylark 320        buick  \n",
       "2       1         plymouth satellite     plymouth  \n",
       "3       1              amc rebel sst          amc  \n",
       "4       1                ford torino         ford  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cars.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mpg                  10.0\n",
      "cylinders               8\n",
      "displacement        360.0\n",
      "horsepower            215\n",
      "weight               4615\n",
      "acceleration         14.0\n",
      "year                   70\n",
      "origin                  1\n",
      "name            ford f250\n",
      "Manufacturer         ford\n",
      "Name: 25, dtype: object\n"
     ]
    }
   ],
   "source": [
    "print(cars.iloc[25,])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzuElEQVR4nO2df5RdZXnvP8+cnCRnwDITHTUMiYNUSUUgkSh0xXoJbUUF6ShFS8XS1oq3S9ctSGMHayVYVk2bUu1d7XIttN5i8dKA4IjELvCScFUEbOIkhAgpVX455JJpYSiQIZnMPPePs/dkn332z3P2Pmefc57PWlkzs8/e+332nsmz3/08z/t9RFUxDMMweoe+dhtgGIZhtBZz/IZhGD2GOX7DMIwewxy/YRhGj2GO3zAMo8dY1G4DkvCqV71KR0ZG2m2GYRhGR7Fz587/UNUh//aOcPwjIyPs2LGj3WYYhmF0FCLyRNB2C/UYhmH0GOb4DcMwegxz/IZhGD2GOX7DMIwewxy/YRhGj9ERVT3dzvjEJJvv3MfT0zMcP1Bhw7knM7pmuN1mGYbRpZjjbzPjE5NcddseZmbnAJicnuGq2/YAmPM3DCMXLNTTZjbfuW/B6bvMzM6x+c59bbLIMIxuxxx/m3l6eibVdsMwjGYxx99mjh+opNpuGIbRLOb428yGc0+mUi7VbKuUS2w49+Q2WWQYRrdjyd024yZwrarHMIxWYY6/AIyuGTZHbxhGy8g91CMiJRGZEJE7nJ83isikiOxy/r0nbxsMwzCMo7Rixv9HwMPAL3i2fUFV/7oFYxuGYRg+cp3xi8gJwHnAV/IcxzAMw0hO3jP+LwKfAl7h2/4JEfkdYAdwpao+l7MdmWDSCoZhdAO5zfhF5HzggKru9H30JeAkYDWwH7gu5PjLRGSHiOyYmprKy8zEuNIKk9MzKEelFcYnJps657pN2zhxbCvrNm1r6lyGYRhJyTPUsw64QEQeB/4ZOEdEblTVZ1R1TlXngS8Dbws6WFWvV9W1qrp2aKiuZWTLyVpaIY8HiWEYRhJyc/yqepWqnqCqI8BvAdtU9RIRWe7Z7X3AQ3nZkCVZSyuEPUiuvHm3vQEYhpEr7ajj/ysRWQ0o8DjwsTbYkJrjBypMBjj5RqUVwh4Yc6qAqXQahpEfLZFsUNV7VPV85/sPq+qpqnqaql6gqvtbYUOzZC2tkOSBYSqdhmHkgWn1JGR0zTCff/+pDA9UEGB4oMLn339qw7PxoAdJEGlDSZYwNgwjDpNsSEGW0gp+jZ4+kYUwj5c0oSRr6mIYRhLM8edEkpp/74PE77QhfSgpqvLIHL9hGC7m+HOgkZl3Fiqd1tTFMIwkmOPPgUZn3s2GkrKuPDIMozux5G4OtGvmbU1dDMNIgjn+HGhXO8WsK48Mw+hOLNSTAxvOPbnpRG2jWFMXwzDiMMefA9ZO0TCMImOOPyds5m0YRlExx18wTPPfMIy8McdfIGzlrWEYrcCqegpE1pr/hmEYQZjjLxC28tYwjFZgoZ4C4Mb16yXaqtjKW8MwssQcf5sJEmfzYitvDcPIGnP8bSYoru8ybFU9hmHkgDn+NhMWvxfg3rFzWmuMYRg9Qe7JXREpiciEiNzh/LxMRL4rIo86XwfztqHItEvXxzCM3qUVVT1/BDzs+XkMuFtV3wDc7fzcs5iipmEYrSZXxy8iJwDnAV/xbP4N4Abn+xuA0TxtKDqmqGkYRqvJO8b/ReBTwCs8216jqvsBVHW/iLw66EARuQy4DGDlypU5m9l6TJrBMIx2kduMX0TOBw6o6s5GjlfV61V1raquHRoayti69uKWcE5Oz6AclWYYn5hst2mGYfQAec741wEXiMh7gKXAL4jIjcAzIrLcme0vBw7kaENi0szAm52tW1N0wzDaSW4zflW9SlVPUNUR4LeAbap6CXA7cKmz26XAt/KyISlpZuBZzNZNmsEwjHbSDq2eTcCvi8ijwK87P7eVNOJoWQipWQmnYRjtpCWOX1XvUdXzne//U1V/VVXf4Hx9thU2RJFmBp7FbN1KOA3DaCemzkm6GXgWs3Ur4TQMo52YZAPpmqMn2fcz43u46YGnmNOjept+3R1/a8bxiUnWbdpm5Z0psJJYw2gMc/yka44et+9nxvdw4/1P1h0X1U3LOm+lx+6ZYTSOqIapwBeHtWvX6o4dO9ptRiJOuuo7NTN9P4P9ZfoXL6p5aGy+cx+TATmC4YGKCbWFsG7TNrtnhhGDiOxU1bX+7Tbjz5gopw/w3MFZnjs4CxydpYbJMlt5ZzhWEmsYjWPJ3YwpiaTaf2Z2LvSYgf5yFiZ1JVYSaxiNY44/Yy4+c0XqY+ZUKZfqnf+LLx8xGYcQrCTWMBrHHH/GXDt6KutOWpbqmOGBCscsro+6zc5rqoVhSXEriE4c28q6Tds68uFiJbGG0TgW48+Y8YlJfvzk84n3d2epV2zZFfh51jHrbqqG8ZfEGoaRDJvxZ0xUD10/JZGFWWqrYtZZSE4YhtHZmOPPmDQz9HnVhRlrq2LWnVoN0w3hKcMoChbqCcG7+rYkwsVnruDa0VNjjzt+oBJYXx62r0uaRWTNEGZfkathuik8ZRhFwBZwBRC2+rZS7uPz7z8t0tn4nVQYlXKpJhnplR84rlJGBKYPznL8QIX1q4bY/shU6gdCkKQBECg5UeTEqC3WMozGCFvAZaGeAG564KnA7TOz87Ha+6NrhrnwjGH8xZnlPmGwvxxYgeLX+J+eqS7ycvX+b7z/ydT6/2F9A4COq4bp1PCUYRQVC/UEELX61t8pK2hWvf2RKfxnmJ1X+hcvYuKz76w7Z5qEcJANQUQlce8dO6fQjt5PJ4anDKPI2Iw/gLjVt+5MM2xWHRbjz3LmGndMN82SbbGWYWSLOf4A4lbfujPNsFl12IMjy5LNuGO6SdLAFmsZRrbkFuoRkaXA94AlzjjfUNWrRWQj8FFgytn106r6nazHb0Sr3XvMMYtLvHS4PvzinWmGzZ7nVKmUS4n0/SFY4z+KJLPdND0GXIqsb2+LtQwjO/Kc8R8CzlHV04HVwLtE5Cznsy+o6mrnXy5OP21DdP8xLx2eo1IucclZK0NnmmGzZ3e/pDNU/4x2oFKuSQRH2RBG2llyFk3kDcPoDFpSziki/cAPgD8E3g28qKp/nfT4tOWcjZT/NXJMUOlmXqWRec/GrWTSMLqPtujxi0gJ2An8IvD3qvqAiLwb+ISI/A6wA7hSVZ8LOPYy4DKAlStXphq3kcRmI8ekXXSV1Hn791u/aohbd07muoCpm5LBhmFEk6vjV9U5YLWIDADfFJE3A18C/hxQ5+t1wO8HHHs9cD1UZ/xpxm2k/K/RksGkseekq0+D9gtaTJakpDMNVjJpGL1DS6p6VHUauAd4l6o+o6pzqjoPfBl4W9bjNVL+l2XJ4PjEJKuvuYuRsa2MjG1lzefu4ppv700kjpampn9yeiaxbk2c1o2VTBpG75BnVc8QMKuq0yJSAX4N+EsRWa6q+53d3gc8lPXYjejeZKWVMz4xyYZbdjM7f/QlxW21GIQ/lJI2tJIk7JPkbaNVWkGGYbSf3JK7InIacANQovpmcbOqfk5E/olqlY8CjwMf8zwIAumkZuthSdIw/MnT1dfcxfRM+IMi6XmS2GSJW8Poblqe3FXVB4E1Ads/nNeY7WZ8YjKV0xeOhmvckMpLh480NHbWiWuXItf2G4bRGKbVkxFuOCWKgUqZY5YsYnJ6BoEFPR839LK03MfsXGNvYHkkrk0O2TC6E5NsyIi4pGy5JGy84BTuHTuH4YFKnYjbzOxcZC4girwS19atyzC6E5vxh5C2EUtU2GSwv8zV7z1lYZacZW38cI6J6zA73fCUhX8MozMxxx+AvxHLnOrCz2HOP6rzVv/iRYn2HaiUOXRkPlE5Z9oVwo1o3YTZ6eYmIL/wj+UWDCM/LNQTQFgjlrDtEBxOcfHr3oSFXjZecEqNvk6Yyqe3SXueBNnpzU24ZB3+Md0gw8gXm/EHENaIJapBizecEjRL9q+0XVruW5jZD1TKbLzglLqa+hPHtgaO5W3SDuEtFpudMQeFiNL2GmiEqNyCzfoNo3nM8QdQEgl08nENWtxwyoljW+tmxVB1jkHCboeOzLPjiWfrHHVUNY7r7IMqhDbcshuEhQqhZsIx/hBR2JqALKUdTDfIMPLFQj0BhDVimVNNJJFwXKUc/IEQKt3w9YC+uutXDQWGhNavGqrp9BXU5tFfFppVOCYo/FPuEw4ePhIqB5GWbmoiYxhFxBx/ANeOnsolZ60MnOEHxZu9Ojirr7mLFw4FL8JSDZdvCIqbb39kKlBTf/sjU6l69LpkMWMO6h2AUNMcvtl4vOkGGUa+tESPv1naKdkQJ3cQFLrJCgEe23Re3fawUFIcWUo0eENNeYxlVT2G0Txt0ePvBuLizWnUNMMIqpSB6JBHGmkIyHbGnORh1+zbhbVaNIz8sFBPDHHx5mYd3GB/mQ+dtTJVaCOqdNR73maak0fJOCd52Fk83jCKS0/M+NOEDeK6X0GtU25k9u2lf/Eirh09lbWvW5bYRn/pqP+NoVIu1awUTkuQRs8VW3Zx+ZZdDCe4XjcBbat7DaOYdH2MP01f3LB9LzxjmO2PTAU6sWZj/GFx/LBrCXo4ZB0Pj5OWDgtNQfXtIuxh2YpFZ4ZhHKVnY/xpFgOF7bv9kanQRKV7jitv3h25wCuMpCGROKXMLB1qXPhKqXf+Xse+btO2wi7AsqSxYfRAjD/NYqBGFw6NrhlmvgGnL8D6VUOR+7ix9su37GqZUmaSh5FCaA6hqAuwTArCMKp0veNPsxiomYVDjSQzFbh152So4/E6qjDycKZJksduueZjm87j3rFzambNRV2AZTLThlElN8cvIktF5EcisltE9orINc72ZSLyXRF51Pk6mJcNkG4xUDMLh5I4yyCiHE+7qme8i7Sg+mbipZWN67OkqG8ihtFq8ozxHwLOUdUXRaQM/EBE/gV4P3C3qm4SkTFgDPiTvIxIo0XfTMPxukobqa7UhehkKFRDDuMTkwvn8PYCiCKJMx2fmOSab+9dWDHsF4SLuh5vArsdjeuzJq4TmcX/jV6hJVU9ItIP/AD4Q+BrwNmqul9ElgP3qGqk9ypCs/W0JaFB1UFLFvWFNlIv9wmbLzqdHU88W9MLIIwkDVjGJybZ8I3ddbo97lh5ObWiOtCoCi8gcfWXYXQKbanqEZESsBP4ReDvVfUBEXmNqu4HcJz/q/O0oVnGJybZePveGocdp3YZFkteWu4Lnf3PzmvdOEGkcUab79wX2MN3dl5zq7Apcp/eqDeRIlciGUbW5Or4VXUOWC0iA8A3ReTNSY8VkcuAywBWrlyZj4ExRNXoRzmFsJjxtCNkFkac008yy09iR9xnzVB0Lf2w0leL/xu9REuqelR1GrgHeBfwjBPiwfl6IOSY61V1raquHRqKLnnMi7jkaphTiKpqGW4wGVsSqaueiSMq8ZtXhU2nOtCiViIZRh7kWdUz5Mz0EZEK8GvAI8DtwKXObpcC38rLhmaJc1Z9IgulmB/68n2MjG1lZGwrk9Mz9PlKYVwZg4OHgyWb43j9UH+odk4YG849mXKpXlq63CeRSeEonZ44OtWBFrUSyTDyIFGoR0TeH7D5eWCPqgbO2IHlwA1OnL8PuFlV7xCR+4CbReQjwJPARQ3Y3RLidHjmVLnqtj38/fZHefTASzWfzSssWdTH4SPzoZo/LgOVMiLhWv1AzfmTxs3dz9JU9UTF6CG+UmfDuScHJkmL7kCLWolkGHmQqKpHRLYCvwxsdzadDdwPvBH4nKr+U14GQvuqerLQ2n/c0eHJQ9c/S319lzA7B/vLvDw7n1jzyByoYbSfZqt65oFfUtVnnJO9BvgScCbwPSBXx98u/LPAZgpf42Lfcc3a05yzGcLOGfQ24l185nf0WT+QDMPIjqQx/hHX6TscAN6oqs8C0aUoHc7omuEFaYJGE7OQLPbtjpV0nDzi5mnP6YaCTP/GMDqHpDP+74vIHcAtzs+/CXxPRI4BpvMwrIiExa9PGFxaF+MHWHfSsthj3dh30tW6LuWS8NKhI4yMbaUkwpxq6nJPSN5/IGrxWVD55sbb93ZUeMfCU0YvkXTG/3HgfwGrgTXADcDHVfUlVV2fk22Fw99o3FWl/O4nz65x8lB1+l//6C/HHju6ZpjPjO/hxvufjHX6g/1lxPmKHq37d49LO9sOUqu8deckF54xXGfnxgtOoewvVYpgemY2Unyu0aqhPDDVTqPXSCzZICKvpRrTnwf+VVX/X56GeSmCZEOenHTVdxLN9N1kblyjlKRJ37iEs581n7srsvIoyXnSNMZpFWnvg2F0CmHJ3UQzfhH5A+BHwPuohnnuF5Hfz9bE3iVpeMdNvMYldZMmfZMstvLOztM4/bDzF1EauVMXnRlGoySN8W8A1qjqfwKIyCuBHwJfzcuwXsKN0cfhJl7j1hccVyknGjfsPN6Fac2UswYlioviZL0x/b6Q+1/0RWeG0ShJHf/PgRc8P78APJW9OcUkz8Tf+MQkSxYJB2fjHf/k9AzrNm2LXAwG8NLhIzUyz2EEJZzh6MK0JYv6muolHNRw/bhKOTBJHOdks/wd+MNNQU6/ExadGUajJF3A9TXgVI7KK1xANfTzbwCq+jd5GQjtjfHnGZNudIGYtwF82Mw/aXx6fGKy4X7BcfiVSMt9wjwwN59OJjrr30FYTL8kwryqVfUYXUNTMX7gp8A41cTuPNUHwNPAK5x/XUueMekkHbaC8DaAD6uzSRo6aaRf8PBApVpZFIP/rLPzWuf0AY5duijSyWb9Owi7N/Oqga0ki0TRKqKMziRpqOc7wKeBEc8xqqqn5WFUkcgzJt3MOdxj47pKJSHsHIP9ZV48dKRG079cEtavGkrULCYp007SOCyck/XvIIt71g6K3OvA6CySzvhvpJrIfT9wvvPvvXkZVSTyVJts5hzusVmoSoad47zTltdP2xW2Pri/IZvDOH6gEllLn/XvoFOVOItYEWV0Jkkd/5SqfltVH1PVJ9x/uVpWEPJ0EkkatK87aVnk+FELw5ISdo7tj0wx6wvNzM5r6rJOl3Kf1MlEu9cS5dSy/h1kcc/aQVEqoozOJ2mo52oR+QpwN9Um6gCo6m25WFUgRtcMs+OJZxfkFEoiXHhGcBen8YlJPn3bgxycnQdABD505kquHT019NxwVOCsf3GJg4fnUKqJxovPXMHa1y1j79N7F5ziMYtLzKty+ZZdXL5l18K5yhGP8KAQinfcMGG1KzznT0J/uQ9Fahy4m+Adjhh3dM1w6FhPT8/kIpkc1omryHRqiMooHkkd/+8Bq4Ay1eQuVP8/d73jH5+Y5NadkwtVL3Oq3LpzkrWvW1bjOMYnJvnkzbvwTpBVWYiFRzn/NNUsLx0OTgY7z5q6uG9QXHjDN3aDsjCbD4sVx60X8FIpl/gLp2l5nIMOut44p9aJjjprOrXXgVE8kjr+01U12HN1OUl7yG6+cx8BBSsA3PTAU6GOP+3YSfDaF3SOoAbsQdcUVufvR6DmLagRfX5zavFYsxgjK5I6/vtF5E2q+pNcrSkgSeOqUXHWRmvkm4ndurPnNOfw7+s6lMtjQj4KbH9kKvCzpJUo/n4EJZGaxGW3OLdmF6LZm4+RBUmTu28HdonIPhF5UET2iMiDeRpWFJJWlETFWUuSXNUy6TnjcMdMW9bpZ3TNcKL+AP6HhltvfvmWXYkrUUbXDC8kchtVHC0ypgJqFIWkM/53pT2xiKwAvga8lmpe4HpV/VsR2Qh8FHCniJ9W1e+kPX+rSBqC2HDuyXUxfpeLz1wRO07QTDBpqCUI13GuXzXE1+9/snYFbUlqYvxB1+S1Z2lU5thBBEbGtjrn6uPIvAaGlFzC3kTCQmtX3rybK7bs6ujwRtKwoWHkTaIZv7eEM0U55xHgSlX9JeAs4OMi8ibnsy+o6mrnX2GdPiQv/RtdM8zffGA1/R4nKQKXnBVe1eMSNhMEFsZOy7BTG3/rzskapy/AB9+6gs0XnR56TX57ZtzMcQTeB97M7Hyk04fwN5GwB8KcasfPkq0c0ygKSWf8qVHV/cB+5/sXRORhoCOnNUnjqo3GX6Nmua52TBSVcinwjSTovG48/trR8Lr1JEllV9cmTNkyzt6wpG2SSqJOnSVbOaZRFJLG+JtCREaodu56wNn0CSdX8FURGQw55jIR2SEiO6amghOH3UKSWW4UYW8kjc4wk8xAXV2btDo/g/3lyMVSSRa1wVGl0k7SrOnUFcNG95HbjN9FRI4FbgUuV9X/EpEvAX9OdfL558B1QF1TF1W9HrgequqcedvZTtLUywcR9qbR6AwziT1JewP46V8cLcjmL1lEqush/AhHH4idollj5ZhGUUjcerGhk4uUgTuAO4Okm503gTtU9c1R5+n21oufGd9Tl4BNyhtefQwHD88HOpLxiUk2fGN3ncjaB9+6gu2PTIWu5B3oL/Piy0fq5BpcvJLI4xOTXHnL7kDVzSAEeGzTeXXbw1YXb7hld6gdfgb7y0x89p2J9jWMXiBMljm3Gb+ICPAPwMNepy8iy534P1RbOT6Ulw2dQFACNowli/o4dORoovUNrz6Gnz/3cnSNvO/Es3Nao6w5OT3Dhlt2gxxd2PXcwVnC+qoP9pe5+r2n1MxS+4CkdUdBbxv+B597HUsW9SV2+q7dSRrQ5EGezXoMI2vyDPWsAz4M7BGRXc62TwMXi8hqqi7pceBjOdpQeJKuzhXgLy88rcaZrNu0LbI8cPOd+xI5zqB9wg7zh2qSjgHB8ezxicnAt52Z2bmGyljbkfQ1uWSj08izqucHENgnpNDlm60maSmfUu9E4pK3eZQJplmx7GWgUmbjBafUXcPmO/elDnEN9pdDFULbURpp9flGp9GSqh4jnKSlfEG1/HGrivMoE0yzYtnL9MwsG2/fW1d9E+WoB/vLgVUwV7/3FAZCGsq3ozTS6vONTiP3qh4jmiSrc90QiRtHdrVs5lTr+tp6wylpVv6WS1K36KrUJzVJ23JJeOnQEU4c28rxAxXWrxri4OEjia91ema2LgQSVRWkWp05u9fq1e85//TldQ3n21Ua2Un1+ZaLMMBm/G0naGXwJWetrKvLBxZW08JRSQblaDzNvwLXe+4oBvvLfPCtK+ricn3OZ+J8RavO211bcOP9T6ZuyuLX6Qmr2++jOpb3Wr36PbfunOTCM4YL0UylU+rzTSvIcMm1nDMrur2cMwnrNm2LrJcfHqjUNVLxE1Q26pZmum8SYeeNGz8N/pJO/yz0pUNHFpx+FEmuuVV0wkw67HdYpPtoZEvLyzmNbMlite21o6ey9nXLUnfASnr+pPhDIP4FaCc6Ym9xhNnUCU64HVguwnAxx98hxK2QTRpPbnSVb7Ori13iQiDjE5OJ9X+CrrkdpZWdUs7ZSbkII18sxt8hRGnYZBFPjotTJ9XQ8R8flK+IazWZxOmHXXNUaWVetGPMRuiUXISRPzbj7xCCOlTNqS40MY9ypu4xfXJ0YZa/rt79es239y4kbJcsOjovCNKZWb9qqEb6wf9znF3esYAa+7yIQH+5tNBvOGxNADQezkgaHgrar1NCKKYVZLiY4+8g0so++0MQXqc6PTNblWqgNhzxskd7319+mVXbvyANIb99XlRrm8x7ZSv8NBLOSBqqCdtvIGRBWRFDKNa60QAL9XQ1cXIQs/NaE46IC1m47RSTSCGH7Ts+McmVN9c7/TREhVEaCWeEXffG2/cm2k8VC6EYHYU5/i4mSajBu09UyCJNDXjYvp8Z35M4hp/Gbi9JO6YlOdf0zGzN9YXt9/zMbOoxDaOdWKini0mjqx+1//EDlVR6NGH73vTAU4mcvtvdK6qmPyqMkjacEXWfvNcXdX8shGJ0Ejbj7yDShFogvhKn3CesXzW0cM7nXjoUuN/6VUOhjnFyeqbOnqiOYnGUS8J1Hzidxzadx71j57DxglPqrsFtwuKOmfa++IkKyXivxapijG7BVu52CP7EItQ2RIk6zl8947LupGX8+MnnY7V8BvvLTB+cjVXRjFsFXIqpzw/S+nevwT2nX5uo3Cc1vQSAhX3iKp68rPncXYH3yL+q1RaHGZ1E2Mpdc/wdQjPL7cOOjXPEjeA626CH1IVnDAcKqyWNhzciG+HvFhbmtNM8WIPE8oYHKoy8ssL9P3tuQVDu4jNXcO3oqansNYwsMcmGDqeZWvFmQi9peXp6JrJePEwyIum50+KtAIoq2Uxa4+5/QHiF47wPpTk92unMnL9RNGzG3yG0a8ZfKZdYsqgvkWhaUntc0oZNVl9zV2I7vAjhidkge6PsSvvWURLhp59/T+BnFjYy8qYdPXdXAF8DXgvMA9er6t+KyDJgCzBCtfXiB1T1ubzs6BbCwiducjbKeQQdK8BZrx+si/GX+4Rjly5i+uBsZNPzPqDk0/BPmugcn5hk4+17a5y4fwbud4rrVw3xUoD2f5AdftyS0iDcUlV3rOMqZV46fGThfH670r51hD1Y4xaN2UPByJM8Qz1HgCtV9cci8gpgp4h8F/hd4G5V3SQiY8AY8Cc52tEVhEkmeGPmYStOR9cMs+OJZ2skmRX48ZPPc+EZw7EyC+MTk3VNNEsl4YNvXZFYosF7rrDmMFFhmaC+vAAIfPCtKxKXivrpX1yqGSvojcJbtppWrK4kwV3r4xbLdYLom9G55Nlzdz+w3/n+BRF5GBgGfgM429ntBuAezPEnwl8rHtds3cv2R6YCG5pvf2QqNjSz+c59dTPq2TlNdGzQuaKqiJ6engncJ8ylzyvcunOy4XzFwcNziXr+ujP99auGFmL3Sbj4zBWR5wvabj18jbxpSR2/iIwAa4AHgNc4DwX34fDqkGMuE5EdIrJjamqqFWZ2HGkSvnkkhxtJtsYdc/xAJfV53faMjZD0ceEuGNv+SPDfYkmEdSctW7CjJMIlZ60MTexG9UuOu9/NrlswjNwdv4gcC9wKXK6q/5X0OFW9XlXXquraoaGh/AzsYOKarTe6b5bHpjnGzRGE7RPl2udUU8tGJ8WbuwhzyvOqXLR2Ja89bikCvPa4pax93bK6/Vyn7a5J8OIuTOsLeYgdP1Cx9olGJuTq+EWkTNXpf11Vb3M2PyMiy53PlwMH8rShmwlbmfv08zOM+GaDzaw6DTrW23g9zawzzOZjFlerh67YsotnQ1YQV8p99IV4f1cfx9XLGewvVxd3RZD0HWFmdo7Lt+xi3aZtDPSXQ0+24ZbdkQ7Z67Shtl+yd2FaUNjK/V11iva/UWxyc/wiIsA/AA+r6t94ProduNT5/lLgW3nZ0O24gmQDlVpn5PoNr/NpRLzMP47XqfobryeddYY1l5/3nG9mNlh2+eDsPCURyqVal+06xdE1w9w7dg6PbTqPic++k80XnR7ZxD5tVmByeoYXXz5SNz5U7/msT1fa75DDchclkUhbBirlhd9Vp2j/G8Umtzp+EXk78H1gD9VyToBPU43z3wysBJ4ELlLVZ6POZXX80WTRiD2L8RodJ21t/EClzDFLFjVd6njSVd9pKCk8UCnzwstHEh3rbSx/4tjW1A8bqL2v1jDdSEPL6/hV9QeEv03/al7j9iJZNGLPYrxGx0l73PMzs+y6+p0NjeWl0Uqg51MsIvPmK8IatsQtpJt01hqMrhkOXc/RiFCcrRXoXUydswuIS7Jm3Qkqy2RvI8dldT3DTdibxAavQx6fmOTFl+sXoJVLVU2fuHxDFiE7L5Yk7m3M8XcBeTdiTzJeM+OkaeSe5fVsOPfk2ARw2PhhNrtn8zvkzXfuq8sBAByzeBHXjp4aGwLy5gu8uYx7x84JXHAXV+5pSeLexkTaCkIzr92NNmJvlEaadocpWnqP88omiMD0wVkq5T5mZucXHKNbYXP5ll0M9pc577TlkauHo8YF6oKR5ZLwtpFB7vvZs3U9gIOavPslr91Dnp6e4fItu9h8577IfgbTM7OcOLY1kW5SVEgsTLo6bNWvJYl7GxNpKwCNau13ClEyDVHXOT4xWacRFIdfhjlq3KXlvlAN/rBYut/WRqSiGyUsgRt1nWHHWpK4NzBZ5gLT7Uv0o2Qaoq4zLDwShfd8ceOGfZZGNqFVM2Q3xBT0ZhgngwFHO5aF6Tx5xwgj7VupJY+Lizn+AtDtr92NVh01WyXU6PFJZBO8++Y94/eGpoLE2+KcPhxdFewed+vOyUQCfS5xaqLN7m+0FkvuFoCsq2SKRqNVR81WCTVyvECkbIR/e5rEdCO4oZewN5gkOkX+dpXuca7IXliS2EvaZLAlj4uNOf4C0O1NvButOmqm6iZu3DAUFurl45q8w9GVyINhUg5N4O23cOLY1tA3izidorBgmbs+II7xicnIfgZZbHfHMfG51mCOvwBkVZtdVLzXB0c16uOuc3TNMJsvOr1OksJlsL9cJ8PgPV+YpEUUro1+m4MqZbzOf+Kz7+SSs1YmHselJBIqKeH2KHZr7aNs9v79DFTKgbISQcTV7rshmzDSvq2Fbbd1Ba3FqnqMriesgsUfAgmrMEpTARPVHrLcJzXJ6rjKrSQVQ0HnWPO5uwKrlcIoiTCvGhjnj7IhriIrTaVaVHvQ6z5wetdMglpNWFWPzfiNricsvKCQ6C0rTdgiSs7BLxoX91YXFRYJO8f4xGQqpw/VUFHYLDvKhri3tTRvsWHjzKnazD8HrKrH6HrSNFpPc7wbtvCWLfaFLMQaHqjUdVDLw+6o5GmSRWL+ktUoG+KuJc31RlVHdVNpc1Ewx290PY0Im3md+YCj7e8P07h19d5zR2npJyHpuGHHRLn1s14/yA9/+mysPIR39h1374Jq9SF8Fbn7mX//lw7V6xiF2RR03WnXCfT6GgOL8Rs9QZr/6EHx6XJJOGbxIp6fma05PkkcPqoFY6PjRh0TRH+5D0US1fz73yjC7l2gvX0CQl2PZvda8PUtiNo/zqZGV7t3+0p5L7Zy1+hp0oQdgmrQZ+eUY5YsqpODTrJILKxPbzPjRh3jp1IusXhRX2DSOSjB7X+jCLt3gfZGrLQOcu5JVmYH2dTMavduXymfBEvuGoaPNMncJIvEkq4gbqT2PWkCOCzprFCzBmHJonqXEFZfn/fK8qikcDOrvbt9pXwSbMZvGD7ikrlegmLgSY5rdty4Y/yhETfm7megUuZlT6vL6ZnZGmmFKOmF4yrl0NLVZolLvDdyr7I4tlvIzfGLyFeB84EDqvpmZ9tG4KOA++77aVX9Tl42GIaXpHH+sISmu4rWn5hcsqgv1PGXSxKY2A2yJSqR6t3fK1vdv7h+xa47pv+YUp8w542vlwQRIsMeYWGRK2/ezdJyuoBBUIw/jPWrhmp+9t+vRkTmXLLsYtap5Nlz9x3Ai8DXfI7/RVX96zTnsuSu0SxpE3pJHE0SRxaU2I2yBYKrXpKKsQH0Cfz2mSvr7PXjrxjy4vYKbrRPsEtQVc/G2/fGvinEyWtXyqVUInN+eqWqpx09d78nIiN5nd8w0pA2oedPaK7btC0w8RpHUGI3ypYgsbSgsaOYV7jpgadia/Zn5zW0tt8rdNeo+mhYuGbznftiHX+Stw5XZK4R0q6p6Dbakdz9hIg8KCJfFZHBsJ1E5DIR2SEiO6amklVFGEYYzSb0smwkn6WwWRhJG8kHCbw1K3TnnsMrMNdIUjhOXruXkrFZ02rH/yXgJGA1sB+4LmxHVb1eVdeq6tqhoaGw3QwjEc1KP2fZSD4rYbMo4qSaXfwCb2FCd0nP547tF5jzykGkvefdLlveDlrq+FX1GVWdU9V54MvA21o5vtG7NCt9HXR8uSSRstFh509rS9pZd7kkXHzmithj3DHjmrePrhnmug+cnuj6K+US133gdLY/MhUazkpyPXFvHb2WjM2alpZzishyVd3v/Pg+4KFWjm/0BlGJuywa2gfJE/irbaLOH3Su9auG2HznPq7Ysqtudawb43bj8QOecQb6y7w8O8eMU5I52F/m6vdWG8Kvfd2yujHSJkO997JS7qtd8KXK4kUlZg8fdfDuOoCwvMDk9AzXfHsvF54xzNYH9y8IylVXFrNwHYeOzNU0q/dXTrmVRTueeDbRqmjvtQRJSfRavD/Pqp6bgLOBVwHPAFc7P6+m+rfzOPAxz4MgFKvqMZLSicvxo6pWgkoWW3UtSeUg/JRLEpv4LvUJfSQr7YwjiSRG1LUU/e+jGVouy6yqF6vqclUtq+oJqvoPqvphVT1VVU9T1QuSOH3DSEMntvwLs/mmB55q67UkkYMIIkm109y8ZuL0oVrBFEfUtRT97yMPTLLB6Co6sQIkSos+zf5ZU+R75iVJBVPctXTKtWaFOX6jq+jECpAw28IqaVp1LUW+Z16SVBzFXUunXGtWmFaP0VU0oh/f7thumM1hMf64awladRyU1I3bL2i1chLSSDNkwZwq6zZti/xdxmkqTU7PsG7TNtavGqpJOg9Uymy84JS2/I3k+bdqevxG15FGP74oib2kTjzuWoIeFn7S7PeWlcclat7iEtRw5bhKmcNH5jjoEYMLok+qq45rbejj8JzWaAxF2Rv1uxyfmOSab+9N3Zqy3Cdsvqi1fX+z+lsNS+6a4zd6hjRN04tOVHPyJDHvPPb76effE7lPkqY1jY7vEve7bMSGJOfNmqz+Vq3ZutHzdGLiN4y0CeEi7Jen9ETSMbKU3siTvP9WzfEbPUMnJn7DSJsQLsJ+eUpPJB0jS+mNPMn7b9Ucv9EzdNPS/7BrSSrVkOV+ABefuaIhm13KfVJNCjc4vnuOuN9lI6JzSc7r71L2oS/fx+uv2srIWPXfL/3ZvyyI1CUh779Vq+oxeoZmZRuKRNS1JJVqaHS/kVdWuP9nzzGnVVnni89ckUg2wWtzkGxC3PXExuYTvBz475tfZqORqp6gLmV+W2dm5/nkll01NqSx06p6DMPoSZIkZtuRqE+TMG61fZbcNQyjo2mmkXqepBmzKIUE5vgNw+gIkjZSbzVpxixKIYE5fsMoKP6EYZrkYCcTdt1xidm45Gde9zNNwthdIdzu36Uldw2jgAQlDK+6bQ+QLDnYqSS57rT9D5Ket1GCErEjr6xw38+erVuJnPXYjWLJXcMoIN20yjgNeV13O+9nO8e25K5hdBDdtMo4DXlddzvvZxF/l+b4DaOAdNMq4zTkdd3tvJ9F/F3m5vhF5KsickBEHvJsWyYi3xWRR52vg3mNbxidTDetMk5DXtfdzvtZxN9lnjP+fwTe5ds2Btytqm8A7nZ+NgzDx+iaYT7//lMZHqggVOPBRZCPzpvRNcNceMbwgkZPSYQLzxhu6rr9DeuhtfeziL/LXJO7IjIC3KGqb3Z+3gecrar7RWQ5cI+qxj72LLlrGL1B1j0TityDoRUUJbn7GrfBuvP11S0e3zCMAhPWeL7RZuhZn69bKGxyV0QuE5EdIrJjamqq3eYYhtECsq6AKWJFTRFoteN/xgnx4Hw9ELajql6vqmtVde3Q0FDLDDQMo31kXQFTxIqaItBqx387cKnz/aXAt1o8vmEYBSbrCpg05+sliYzcJBtE5CbgbOBVIvJz4GpgE3CziHwEeBK4KK/xDcPoPLLWoU96vl6TyDDJBsMwep5ulcgoSlWPYRhG4ei1JLA5fsMwep5eSwKb4zcMo+cpoqxCnpgev2EYPU/ezc2Lhjl+wzAMqs6/Wx29Hwv1GIZh9Bjm+A3DMHoMc/yGYRg9hjl+wzCMHsMcv2EYRo/REZINIjIFPOHb/CrgP9pgTqN0kr1maz6YrflgtobzOlWtkzfuCMcfhIjsCNKgKCqdZK/Zmg9maz6YremxUI9hGEaPYY7fMAyjx+hkx399uw1ISSfZa7bmg9maD2ZrSjo2xm8YhmE0RifP+A3DMIwGMMdvGIbRYxTK8YvIV0XkgIg85Nm2UUQmRWSX8+89ns+uEpF/F5F9InKuZ/sZIrLH+ex/iojkYOsKEdkuIg+LyF4R+SNn+zIR+a6IPOp8HWy3vRG2Fu7eishSEfmRiOx2bL3G2V7E+xpma+Huq2eckohMiMgdzs+Fu68Rthb5vj7ujLNLRHY42wp7b1HVwvwD3gG8BXjIs20j8McB+74J2A0sAU4EfgqUnM9+BPwyIMC/AO/OwdblwFuc718B/Jtj018BY872MeAv221vhK2Fu7fOeY91vi8DDwBnFfS+htlauPvqseGTwP8G7nB+Ltx9jbC1yPf1ceBVvm2FvbeFmvGr6veAZxPu/hvAP6vqIVV9DPh34G0ishz4BVW9T6t38mvAaA627lfVHzvfvwA8DAw7dt3g7HaDZ+y22RthaxjttFVV9UXnx7LzTynmfQ2zNYy2/s2KyAnAecBXfDYV6r5G2BpGW22Nsatw9xYKFuqJ4BMi8qBUQ0Hu69Iw8JRnn58724ad7/3bc0NERoA1VGd8r1HV/VB1uMCri2Svz1Yo4L11XvF3AQeA76pqYe9riK1QwPsKfBH4FDDv2VbI+xpiKxTzvkL1gX+XiOwUkcucbUW9tx3h+L8EnASsBvYD1znbg2JfGrE9F0TkWOBW4HJV/a+oXQO2tdTeAFsLeW9VdU5VVwMnUJ0JvTli9yLaWrj7KiLnAwdUdWfSQwK2tdvWwt1XD+tU9S3Au4GPi8g7IvZtu72Fd/yq+ozzn2se+DLwNuejnwMrPLueADztbD8hYHvmiEiZqiP9uqre5mx+xnllw/l6oAj2Btla5Hvr2DcN3AO8i4Le1yBbC3pf1wEXiMjjwD8D54jIjRTzvgbaWtD7CoCqPu18PQB807GtiPd2weBC/QNGqE3uLvd8fwXV2BjAKdQmSH7G0QTJv1JNsrkJkvfkYKdQjcF90bd9M7UJnb9qt70Rthbu3gJDwIDzfQX4PnB+Qe9rmK2Fu68+u8/maMK0cPc1wtZC3lfgGOAVnu9/SHWyUth7m8svq4kbeBPVV7hZqk+/jwD/BOwBHgRu9/3y/5RqRnwfnuw3sBZ4yPns73BWKGds69upvoY9COxy/r0HeCVwN/Co83VZu+2NsLVw9xY4DZhwbHoI+KyzvYj3NczWwt1Xn91nc9SZFu6+RthayPsKvJ6qI98N7AX+tOj31iQbDMMweozCx/gNwzCMbDHHbxiG0WOY4zcMw+gxzPEbhmH0GOb4DcMwegxz/IbRACLyFRF5U8w+/ygivxmwfUREfjs/6wwjGnP8htEAqvoHqvqTBg8fAczxG23DHL/R04jIp0Tkfzjff0FEtjnf/6qI3Cgi7xSR+0TkxyJyi6N3hIjcIyJrne8/IiL/5mz7soj8nWeId4jID0XkZ57Z/ybgVxzt9itaeLmGAZjjN4zvAb/ifL8WONbRNXo71VWinwF+TasCXDuoasQvICLHA39GdZn9rwOrfOdf7pzrfKoOH6rL97+vqqtV9QuZX5FhxLCo3QYYRpvZCZwhIq8ADgE/pvoA+BWqsgBvAu51GiEtBu7zHf824P+q6rMAInIL8EbP5+NaFRX7iYi8Js8LMYykmOM3ehpVnXVUIH+PqrjWg8B6qvK/j1HV2L844hRxrfEOpdjXMFqChXoMoxru+WPn6/eB/05VyO5+YJ2I/CKAiPSLyBt9x/4I+G8iMigii4ALE4z3AtUWmIbRFszxG0bV2S8H7lPVZ4CXqcbgp4DfBW4SkQepPghqYviqOgn8BdWOZv8H+AnwfMx4DwJHpNqk3ZK7RssxdU7DaBIROVZVX3Rm/N8Evqqq32y3XYYRhs34DaN5Njp9dx+imhcYb6s1hhGDzfgNwzB6DJvxG4Zh9Bjm+A3DMHoMc/yGYRg9hjl+wzCMHsMcv2EYRo/x/wGG3W8y8kbW5wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(\"weight\", \"mpg\",data=cars)\n",
    "plt.xlabel(\"weight\")\n",
    "plt.ylabel(\"mpg\")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "code_folding": []
   },
   "source": [
    "\n",
    "## Hypothesis Tests\n",
    "\n",
    "We have learned about the **bootstrap** as a slick way of resampling your data to obtain sampling distributions of various measures of interest, \n",
    "Without having to learn highly specific distributions (such as the $\\chi^2$, Poisson, Binomial or F-distribution) the bootstrap enables us to \n",
    "\n",
    "\n",
    "1. get **confidence intervals**\n",
    "2. perform **one-sample tests**\n",
    "3. perform **two-sample tests**\n",
    "\n",
    "Imagine the EPA requiring the average mpg for 4-cylinder cars to be at least $\\mu_0 = 30$ and needs to decide -based on this sample only- whether the manufacturers need to implement some improvements.\n",
    "In statistical jargon: can the EPA **reject the claim** that the **true mean** is at least 30 ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "29.28391959798995"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cars4=cars[cars[\"cylinders\"]==4]\n",
    "empirical_mean = np.mean(cars4.mpg)\n",
    "empirical_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.656280603443601"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "empirical_sd = np.std(cars4.mpg)\n",
    "empirical_sd"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute a confidence interval of the true mpg of 4-cyl cars via the bootstrap !\n",
    "(Is there just \"THE CI\" ??)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg_bs = draw_bs_reps(cars4.mpg, func = np.mean, size=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([28.498, 30.093])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.percentile(mpg_bs, [2.5, 97.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.959963984540054"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "norm.ppf(0.975)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------\n",
    "\n",
    "**There, we have it, the empirical mean is in fact below 30.**\n",
    "Is this sufficient evidence for the EPA to demand improvements ? \n",
    "Does $\\bar{x} < \\mu_0$ \"prove\" that $\\mu < \\mu_0$ ??\n",
    "\n",
    "Think about an experiment which tries to establish whether a pair of dice in a casino is biased (\"gezinkt\" in German). You toss the dice $50$ times and observe an average number of pips of $\\bar{x} = 3.8$ which is clearly \"way above\" the required value of $\\mu_0 = 3.5$ for a clean pair of dice. Should you go to the authorities and file a complaint ? Since no one really knows the true value of $\\mu$ for these dice, how would a court decide whether they are unbiased ?\n",
    "\n",
    "1. **Innocent until proven guilty**: Always assume the opposite of what you want to illustrate ! That way you can gather evidence against this **Null hypothesis**, e.g. $H_0: \\mu = 3.5 (\\equiv \\mu_0)$. \n",
    "2. Define a **test statistic** which is pertinent to the question you are seeking to answer (computed entirely from your sample). In this case, your test statistics could be e.g. the (scaled?) difference between the observed sample mean and the claim: $\\bar{x} - \\mu_0$.\n",
    "3. We define the **p-value** as the probability of observing a test statistic equally or more extreme than the one observed, given that $H_0$ is true.\n",
    "4. If the observed value of your test statistic seems very unlikely under $H_0$ we would favor the **alternative hypothesis** which is usually complementary, e.g. $H_A: \\mu \\neq 3.5$ \n",
    "\n",
    "\n",
    "Notice that we used vague expressions such as *unlikely* in the procedure above. In a legal context, we would want to make statements such as \"beyond a reasonable doubt\". Well, in statistics we ry to quantifiy that notion by limiting the type-I error probability $\\alpha$  to be at most e.g. $0.01$ or $0.05$. So none our decisions are ever \"proofs\" or free of error.\n",
    "\n",
    "It also means that you, the person conducting the hypothesis test, need to specify a value of $\\alpha$, which clearly introduces a certain amount of subjectivity. Later on, we will discuss the inherent tradeoffs between type-I and type-II errors: only then will one fully understand the non-arbitrary choice of $\\alpha$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parametric Tests\n",
    "\n",
    "1. From the sample, compute $\\bar{x}$ and $\\hat{\\sigma}$.\n",
    "2. Compute the test statistic \n",
    "$$t = \\frac{\\bar{x} - \\mu_0}{\\hat{\\sigma}/\\sqrt{n}}$$\n",
    "3. Reject $H_0$ \n",
    "    * Two-sided test: if $|t| > t_{n-1, 1-\\alpha/2}$ \n",
    "    * Left-sided test: if $t   < t_{n-1, \\alpha}$ \n",
    "    * Right-sided test: if $t   > t_{n-1, 1-\\alpha}$\n",
    "    * If $n>50$ one may replace the t distribution with the normal\n",
    "4. Define the p-value as twice the tail area of a t distribution \n",
    "5. Alternatively one rejects the Null if p-val $< \\alpha$.\n",
    "\n",
    "Back to the casino example, assume the following empirical data:\n",
    "\n",
    "$\\bar{x} = 3.8, n = 100, \\hat{\\sigma} = 1.7$\n",
    "\n",
    "$H_0: \\mu = 3.5$, $H_A: \\mu \\neq 3.5$, $\\alpha = 0.01$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "two-sided cv:  2.5758293035489004\n",
      "one-sided cv:  2.3263478740408408\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "# that is the critical value that we compare our test statistic with\n",
    "print('two-sided cv: ',norm.ppf(1-0.01/2)) # returns the norm score\n",
    "print('one-sided cv: ',norm.ppf(1-0.01))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Two sided test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The test statistic t = 2.496 is not larger than the critical value 2.576 so we fail to reject the Null\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.840131867824506"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We care about both sides since we do not know how the casino is using the biased dice.\n",
    "a=0.01\n",
    "n=200\n",
    "\n",
    "critValue = norm.ppf(1-a/2)\n",
    "z = np.abs((3.8-3.5)/(1.7/np.sqrt(n))) \n",
    "\n",
    "if (z < critValue):\n",
    "    s = \"not\"\n",
    "    s2 = \"fail to\"\n",
    "else:\n",
    "    s = \"\"   \n",
    "    s2 = \"\"\n",
    "    \n",
    "print(\"The test statistic t =\", np.round(z,3), \"is\", s ,\"larger than the critical value\",np.round(critValue,3), \"so we\",s2,\"reject the Null\" ) \n",
    "norm.cdf(1-a/2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### One sided test!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The test statistic t = 2.496 is  larger than the critical value 2.326 so we  reject the Null\n"
     ]
    }
   ],
   "source": [
    "# We care only about one side since we do know that the casino makes \n",
    "# money off upward biased dice.\n",
    "\n",
    "critValue = norm.ppf(1-a)\n",
    "z = (3.8-3.5)/(1.7/np.sqrt(n))\n",
    "if (z < critValue):\n",
    "    s = \"not\"\n",
    "    s2 = \"fail to\"\n",
    "else:\n",
    "    s = \"\"   \n",
    "    s2 = \"\"\n",
    "    \n",
    "print(\"The test statistic t =\", np.round(z,3), \"is\", s ,\"larger than the critical value\",np.round(critValue,3), \"so we\",s2,\"reject the Null\" ) \n",
    "#norm.cdf(1-a)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### p-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The p value  0.006 is less than alpha 0.01 so we reject the Null\n"
     ]
    }
   ],
   "source": [
    "#p-value is the right tail probability of values equal or larger than your test statistic\n",
    "pVal = 1-norm.cdf(z)\n",
    "print(\"The p value \", np.round(pVal,3), \"is less than alpha\", a, \"so we reject the Null\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-------\n",
    "\n",
    "### Non parametric Tests\n",
    "\n",
    "Parametric Tests require many assumptions, that may not always be true, and are also a bit abstract. It often helps (i) robustness and (ii) the understanding of the testing process to use simulations instead. \n",
    "We now learn about two different ways of such simulations, each following the basic \"pipeline\":\n",
    "\n",
    "1. Clearly define the Null hypothesis.\n",
    "2. Define the test statistic.\n",
    "3. Generate many sets of simulated data (bootstrap replicates) **assuming the null hypothesis is true**\n",
    "4. Compute the test statistic for each simulated data set\n",
    "5. The p-value is the fraction of your simulated data sets for which the test statistic is at least as extreme as for the real data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Bootstrap Hypothesis Tests\n",
    "\n",
    "Besides computing confidence intervals we can also use the bootstrap to perform hypothesis test.\n",
    "We need to \"fake\" the process of drawing new samples again and again under the Null hypothesis, which might appear impossible since our one sample likely will have a sample mean not equal to $\\mu_0$.\n",
    "\n",
    "The trick is to \"shift\" our data such that $\\bar{x} = \\mu_0$!\n",
    "As a test statistic we use the difference of the mean of the bootstrap value minus $\\mu_0$. Notice that there is no need to scale this difference!\n",
    "\n",
    "We illustrate this with our casino example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5, 6, 4, 5, 2, 2, 3, 2, 5, 6, 4, 1, 3, 4, 5, 5, 6, 5, 5, 4])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.choice(6,20)+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The sample mean is 3.65 , so we shift our data by -0.15\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(123)\n",
    "n=400\n",
    "#x = np.random.choice(6,n,p=[0.14,0.14,0.15,0.18,0.19,0.20])+1 #add a bias\n",
    "x = np.random.choice(6,n)+1 \n",
    "#cheating:\n",
    "x = x + 0.2 #add a bias\n",
    "a=0.01\n",
    "\n",
    "mu0=3.5\n",
    "xBar = np.round(np.mean(x),2)\n",
    "obsDiff = np.round(mu0-xBar,2)\n",
    "print(\"The sample mean is\", xBar, \", so we shift our data by\", obsDiff)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAARkklEQVR4nO3df4xld1nH8ffHLa0oSrd0XMvuxi24xCwJLGQsNf5CK/0Z3ZIQbI2yIU3WH23U6D+LmFRRkmJUIkmtrnbj1iClgthNW6nL2gRJLO0Uy9JtrR1gSXdduiOFgmKqi49/zHf1uszs3JmduXeG7/uV3NxznvM99zxnd/K5Z84990yqCklSH75p3A1IkkbH0Jekjhj6ktQRQ1+SOmLoS1JHzhl3A2dy4YUX1pYtW8bdhiStKY888si/VtXEXMtWdehv2bKFqampcbchSWtKks/Nt2zB0ztJvjnJQ0k+meRwkt9s9YuTfDzJdJL3Jzm31c9r89Nt+ZaB13pbqz+Z5Ipl2DdJ0iIMc07/eeBHq+rVwHbgyiSXAu8C3l1V3w18Ebihjb8B+GKrv7uNI8k24DrglcCVwB8mWbeM+yJJWsCCoV+z/q3NvqA9CvhR4AOtvg+4tk3vaPO05ZclSavfWVXPV9VngWngkuXYCUnScIa6eifJuiSPAieAA8CngS9V1ck25CiwsU1vBJ4GaMufA14yWJ9jncFt7UoylWRqZmZm0TskSZrfUKFfVV+rqu3AJmaPzr9npRqqqj1VNVlVkxMTc374LElaokVdp19VXwIeAL4POD/Jqat/NgHH2vQxYDNAW/5i4AuD9TnWkSSNwDBX70wkOb9NvxB4A/AEs+H/pjZsJ3B3m97f5mnL/65mb+W5H7iuXd1zMbAVeGiZ9kOSNIRhrtO/CNjXrrT5JuCuqronyePAnUl+G/hH4PY2/nbgz5NMA88ye8UOVXU4yV3A48BJ4Maq+try7o4k6Uyymu+nPzk5WX45S5IWJ8kjVTU517JV/Y1caTXbsvvesWz3yC3XjGW7+sbgDdckqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0JakjfiNXa9q4vhUrrVUe6UtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktSRBUM/yeYkDyR5PMnhJL/U6r+R5FiSR9vj6oF13pZkOsmTSa4YqF/ZatNJdq/MLkmS5jPMH1E5CfxqVX0iybcBjyQ50Ja9u6p+d3Bwkm3AdcArgZcCH0nyirb4VuANwFHg4ST7q+rx5dgRSdLCFgz9qjoOHG/TX0nyBLDxDKvsAO6squeBzyaZBi5py6ar6jMASe5sYw19SRqRRZ3TT7IFeA3w8Va6KcmhJHuTrG+1jcDTA6sdbbX56qdvY1eSqSRTMzMzi2lPkrSAoUM/yYuADwK/XFVfBm4DXg5sZ/Y3gd9bjoaqak9VTVbV5MTExHK8pCSpGeoPoyd5AbOB/96q+iuAqnpmYPmfAPe02WPA5oHVN7UaZ6hLkkZgmKt3AtwOPFFVvz9Qv2hg2BuBx9r0fuC6JOcluRjYCjwEPAxsTXJxknOZ/bB3//LshiRpGMMc6X8/8DPAp5I82mq/BlyfZDtQwBHgZwGq6nCSu5j9gPYkcGNVfQ0gyU3A/cA6YG9VHV62PZEkLWiYq3c+BmSORfedYZ13Au+co37fmdaTJK0sv5ErSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHVkwdBPsjnJA0keT3I4yS+1+gVJDiR5qj2vb/UkeU+S6SSHkrx24LV2tvFPJdm5crslSZrLMEf6J4FfraptwKXAjUm2AbuBg1W1FTjY5gGuAra2xy7gNph9kwBuBl4HXALcfOqNQpI0GguGflUdr6pPtOmvAE8AG4EdwL42bB9wbZveAdxRsx4Ezk9yEXAFcKCqnq2qLwIHgCuXc2ckSWe2qHP6SbYArwE+DmyoquNt0eeBDW16I/D0wGpHW22++unb2JVkKsnUzMzMYtqTJC1g6NBP8iLgg8AvV9WXB5dVVQG1HA1V1Z6qmqyqyYmJieV4SUlSM1ToJ3kBs4H/3qr6q1Z+pp22oT2faPVjwOaB1Te12nx1SdKIDHP1ToDbgSeq6vcHFu0HTl2BsxO4e6D+lnYVz6XAc+000P3A5UnWtw9wL281SdKInDPEmO8Hfgb4VJJHW+3XgFuAu5LcAHwOeHNbdh9wNTANfBV4K0BVPZvkt4CH27h3VNWzy7ETkqThLBj6VfUxIPMsvmyO8QXcOM9r7QX2LqZBSdLy8Ru5ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjw9xPX9IqsmX3vWPb9pFbrhnbtrU8PNKXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6smDoJ9mb5ESSxwZqv5HkWJJH2+PqgWVvSzKd5MkkVwzUr2y16SS7l39XJEkLGeZI/8+AK+eov7uqtrfHfQBJtgHXAa9s6/xhknVJ1gG3AlcB24Dr21hJ0ggteJfNqvpoki1Dvt4O4M6qeh74bJJp4JK2bLqqPgOQ5M429vHFtyxJWqqzOad/U5JD7fTP+lbbCDw9MOZoq81X/zpJdiWZSjI1MzNzFu1Jkk631NC/DXg5sB04DvzecjVUVXuqarKqJicmJpbrZSVJLPGPqFTVM6emk/wJcE+bPQZsHhi6qdU4Q12SNCJLOtJPctHA7BuBU1f27AeuS3JekouBrcBDwMPA1iQXJzmX2Q979y+9bUnSUix4pJ/kfcDrgQuTHAVuBl6fZDtQwBHgZwGq6nCSu5j9gPYkcGNVfa29zk3A/cA6YG9VHV7unZEkndkwV+9cP0f59jOMfyfwzjnq9wH3Lao7SdKy8hu5ktQRQ1+SOmLoS1JHDH1J6siSrtOXTrdl973jbkHSEDzSl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHVkw9JPsTXIiyWMDtQuSHEjyVHte3+pJ8p4k00kOJXntwDo72/inkuxcmd2RJJ3JMEf6fwZceVptN3CwqrYCB9s8wFXA1vbYBdwGs28SwM3A64BLgJtPvVFIkkZnwdCvqo8Cz55W3gHsa9P7gGsH6nfUrAeB85NcBFwBHKiqZ6vqi8ABvv6NRJK0wpZ6Tn9DVR1v058HNrTpjcDTA+OOttp89a+TZFeSqSRTMzMzS2xPkjSXs/4gt6oKqGXo5dTr7amqyaqanJiYWK6XlSSx9NB/pp22oT2faPVjwOaBcZtabb66JGmElhr6+4FTV+DsBO4eqL+lXcVzKfBcOw10P3B5kvXtA9zLW02SNELnLDQgyfuA1wMXJjnK7FU4twB3JbkB+Bzw5jb8PuBqYBr4KvBWgKp6NslvAQ+3ce+oqtM/HJYkrbAFQ7+qrp9n0WVzjC3gxnleZy+wd1HdSZKWld/IlaSOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakj54y7AUlrx5bd945lu0duuWYs2/1G5JG+JHXE0JekjpxV6Cc5kuRTSR5NMtVqFyQ5kOSp9ry+1ZPkPUmmkxxK8trl2AFJ0vCW40j/R6pqe1VNtvndwMGq2gocbPMAVwFb22MXcNsybFuStAgrcXpnB7CvTe8Drh2o31GzHgTOT3LRCmxfkjSPsw39Av42ySNJdrXahqo63qY/D2xo0xuBpwfWPdpq/0+SXUmmkkzNzMycZXuSpEFne8nmD1TVsSTfARxI8k+DC6uqktRiXrCq9gB7ACYnJxe1riTpzM7qSL+qjrXnE8CHgEuAZ06dtmnPJ9rwY8DmgdU3tZokaUSWHPpJvjXJt52aBi4HHgP2AzvbsJ3A3W16P/CWdhXPpcBzA6eBJEkjcDandzYAH0py6nX+oqo+nORh4K4kNwCfA97cxt8HXA1MA18F3noW25YkLcGSQ7+qPgO8eo76F4DL5qgXcONStydJOnvee+cbzLjujSJpbfA2DJLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JH/Bu5kla9cf7t5yO3XDO2ba8Ej/QlqSOGviR1xNCXpI4Y+pLUET/IXQHj/NBJks7EI31J6sjIQz/JlUmeTDKdZPeoty9JPRtp6CdZB9wKXAVsA65Psm2UPUhSz0Z9Tv8SYLqqPgOQ5E5gB/D4SmzMc+uSzta4cmSlvhQ26tDfCDw9MH8UeN3ggCS7gF1t9t+SPDmi3hZyIfCv425iCdZi32uxZ1ibfa/FnqGDvvOus9rOd823YNVdvVNVe4A94+7jdEmmqmpy3H0s1lrsey32DGuz77XYM9j32Rj1B7nHgM0D85taTZI0AqMO/YeBrUkuTnIucB2wf8Q9SFK3Rnp6p6pOJrkJuB9YB+ytqsOj7OEsrLpTTkNai32vxZ5hbfa9FnsG+16yVNW4e5AkjYjfyJWkjhj6ktQRQ38eSS5IciDJU+15/RxjvivJJ5I8muRwkp8bR6+n9TRM39uT/EPr+VCSnxxHrwP9LNhzG/fhJF9Kcs+oezytjzPeSiTJeUne35Z/PMmWMbR5ek8L9fxD7Wf5ZJI3jaPHuQzR968kebz9HB9MMu/16aMyRM8/l+RTLTc+NvK7ElSVjzkewO8Au9v0buBdc4w5FzivTb8IOAK8dA30/Qpga5t+KXAcOH8199yWXQb8OHDPGHtdB3waeFn7//8ksO20Mb8A/FGbvg54/5h/JobpeQvwKuAO4E3j7HeRff8I8C1t+ufXyL/1tw9M/wTw4VH26JH+/HYA+9r0PuDa0wdU1X9W1fNt9jxWx29Ow/T9z1X1VJv+F+AEMDGqBuewYM8AVXUQ+MqIeprP/95KpKr+Ezh1K5FBg/vzAeCyJBlhj6dbsOeqOlJVh4D/HkeD8xim7weq6qtt9kFmv/szTsP0/OWB2W8FRno1zWoIqdVqQ1Udb9OfBzbMNSjJ5iSHmL29xLtaiI7TUH2fkuQSZo9IPr3SjZ3Bonoes7luJbJxvjFVdRJ4DnjJSLqb2zA9r0aL7fsG4G9WtKOFDdVzkhuTfJrZ33J/cUS9AavwNgyjlOQjwHfOsejtgzNVVUnmfDeuqqeBVyV5KfDXST5QVc8sf7f/Zzn6bq9zEfDnwM6qWtEjvOXqWZpLkp8GJoEfHncvw6iqW4Fbk/wU8OvAzlFtu+vQr6ofm29ZkmeSXFRVx1s4nljgtf4lyWPADzL7K/2KWY6+k3w7cC/w9qp6cIVa/V/L+W89ZsPcSuTUmKNJzgFeDHxhNO3Naa3e/mSovpP8GLMHDz88cLp1XBb7b30ncNuKdnQaT+/Mbz//9+67E7j79AFJNiV5YZteD/wAMO67gg7T97nAh4A7qmpF36CGtGDPq8gwtxIZ3J83AX9X7VO7MVmrtz9ZsO8krwH+GPiJqloNBwvD9Lx1YPYa4KkR9ufVO/M9mD0He7D9h3wEuKDVJ4E/bdNvAA4x+wn9IWDXGun7p4H/Ah4deGxfzT23+b8HZoD/YPZc6RVj6vdq4J+Z/Rzk7a32DmaDB+Cbgb8EpoGHgJetgp+LhXr+3vZv+u/M/lZyeNw9D9n3R4BnBn6O96+Bnv8AONz6fQB45Sj78zYMktQRT+9IUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktSR/wEhsyyDaLetQwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Alternatively, if we had not shifted our data\n",
    "#TS = bs_mean_dice-xBar\n",
    "\n",
    "#The test statistic\n",
    "bs_mean_dice = draw_bs_reps(x+obsDiff,np.mean,10000)\n",
    "TS = bs_mean_dice-mu0\n",
    "plt.hist(TS);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The two sided p value of 0.012 is not smaller than alpha= 0.01 so we fail to reject the Null.\n"
     ]
    }
   ],
   "source": [
    "#two sided p value\n",
    "pVal = np.mean( abs(TS) > abs(obsDiff))\n",
    "print(\"The two sided p value of\", np.round(pVal,3), \"is not smaller than alpha=\",a,\"so we fail to reject the Null.\" ) "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------\n",
    "\n",
    "#### Tasks\n",
    "- Test $H_0: \\mu \\geq 30, H_A: \\mu < 30$ for the mean mpg of 4-cylinder cars\n",
    "    * using bootstrap replicates\n",
    "    * via standard testing theory.\n",
    "- Compute the corresponding **p-values**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu0=30\n",
    "empirical_mean = np.mean(cars4.mpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#assume the claim is true !!\n",
    "shift = mu0-empirical_mean\n",
    "cars4Shifted= cars4.mpg+shift\n",
    "bs_mean_mpg = draw_bs_reps(cars4Shifted,np.mean,10000)\n",
    "\n",
    "#The test statistic\n",
    "TS = bs_mean_mpg-mu0\n",
    "pVal = np.mean( TS < -shift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARS0lEQVR4nO3db8yddX3H8fdHUEY2yWAUrG1j0XSLhUQcdzo2k4VNIx0+KD5gqQ+kD0iqBLOZ+KRsyf48aILL5hKWQdJFQkk2WRN1NApTbDTGBMEbA5aCjCqddG3oPd0mPmFSv3tw/7qclNP7nPtPz7lvfu9XcnJd53t+13X9fufSD6e/c53rTlUhSerDm6bdAUnS5Bj6ktQRQ1+SOmLoS1JHDH1J6siF0+7AKJdffnlt3rx52t2QpDXlySef/M+qWnd2fdWH/ubNm5mdnZ12NyRpTUny78PqTu9IUkcMfUnqiKEvSR0ZGfpJfinJE0meTnIkyV+2+mVJHk3yQlteOrDNnUmOJnk+yY0D9euSHG6v3Z0k52dYkqRhxvmk/yrw+1X1HuBaYHuS64E9wKGq2gIcas9JshXYCVwNbAfuSXJB29e9wG5gS3tsX7mhSJJGGRn6Ne9n7emb26OAHcD+Vt8P3NzWdwAPVtWrVfUicBTYlmQ9cElVPVbzd3l7YGAbSdIEjDWnn+SCJE8Bp4BHq+px4MqqOgnQlle05huAlwY2P95qG9r62fVhx9udZDbJ7Nzc3CKGI0layFihX1Wnq+paYCPzn9qvWaD5sHn6WqA+7Hj7qmqmqmbWrXvdbwskSUu0qKt3quq/gW8wPxf/cpuyoS1PtWbHgU0Dm20ETrT6xiF1SdKEjPxFbpJ1wM+r6r+TXAx8APg0cBDYBdzVlg+1TQ4C/5TkM8Dbmf/C9omqOp3klfYl8OPArcDfrfSA1J/Ne748leMeu+tDUzmutBzj3IZhPbC/XYHzJuBAVX0pyWPAgSS3AT8CbgGoqiNJDgDPAq8Bd1TV6bav24H7gYuBR9pDkjQhI0O/qr4HvHdI/cfA+8+xzV5g75D6LLDQ9wGSpPPIX+RKUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1ZJxbK0saYlr38Qfv5a+l85O+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjoyMvSTbEry9STPJTmS5I9b/S+S/EeSp9rjpoFt7kxyNMnzSW4cqF+X5HB77e4kOT/DkiQNM85dNl8DPlVV303yVuDJJI+21/62qv56sHGSrcBO4Grg7cDXkvx6VZ0G7gV2A98GHga2A4+szFAkSaOM/KRfVSer6rtt/RXgOWDDApvsAB6sqler6kXgKLAtyXrgkqp6rKoKeAC4ebkDkCSNb1Fz+kk2A+8FHm+lTyT5XpL7klzaahuAlwY2O95qG9r62fVhx9mdZDbJ7Nzc3GK6KElawNihn+RXgM8Dn6yqnzI/VfMu4FrgJPA3Z5oO2bwWqL++WLWvqmaqambdunXjdlGSNMJYoZ/kzcwH/j9W1RcAqurlqjpdVb8A/gHY1pofBzYNbL4RONHqG4fUJUkTMs7VOwE+CzxXVZ8ZqK8faPZh4Jm2fhDYmeSiJFcBW4Anquok8EqS69s+bwUeWqFxSJLGMM7VO+8DPgocTvJUq/0J8JEk1zI/RXMM+BhAVR1JcgB4lvkrf+5oV+4A3A7cD1zM/FU7XrkjSRM0MvSr6lsMn49/eIFt9gJ7h9RngWsW00FJ0srxF7mS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6MjL0k2xK8vUkzyU5kuSPW/2yJI8meaEtLx3Y5s4kR5M8n+TGgfp1SQ631+5OkvMzLEnSMON80n8N+FRVvRu4HrgjyVZgD3CoqrYAh9pz2ms7gauB7cA9SS5o+7oX2A1saY/tKzgWSdIII0O/qk5W1Xfb+ivAc8AGYAewvzXbD9zc1ncAD1bVq1X1InAU2JZkPXBJVT1WVQU8MLCNJGkCFjWnn2Qz8F7gceDKqjoJ8/9hAK5ozTYALw1sdrzVNrT1s+vDjrM7yWyS2bm5ucV0UZK0gLFDP8mvAJ8HPllVP12o6ZBaLVB/fbFqX1XNVNXMunXrxu2iJGmEsUI/yZuZD/x/rKovtPLLbcqGtjzV6seBTQObbwROtPrGIXVJ0oSMc/VOgM8Cz1XVZwZeOgjsauu7gIcG6juTXJTkKua/sH2iTQG9kuT6ts9bB7aRJE3AhWO0eR/wUeBwkqda7U+Au4ADSW4DfgTcAlBVR5IcAJ5l/sqfO6rqdNvuduB+4GLgkfaQJE3IyNCvqm8xfD4e4P3n2GYvsHdIfRa4ZjEdlCStHH+RK0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqyDh/OUsaafOeL0+7C5LG4Cd9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqyMjQT3JfklNJnhmo/UWS/0jyVHvcNPDanUmOJnk+yY0D9euSHG6v3Z0kKz8cSdJCxvmkfz+wfUj9b6vq2vZ4GCDJVmAncHXb5p4kF7T29wK7gS3tMWyfkqTzaGToV9U3gZ+Mub8dwINV9WpVvQgcBbYlWQ9cUlWPVVUBDwA3L7HPkqQlWs6c/ieSfK9N/1zaahuAlwbaHG+1DW397PpQSXYnmU0yOzc3t4wuSpIGLTX07wXeBVwLnAT+ptWHzdPXAvWhqmpfVc1U1cy6deuW2EVJ0tmWFPpV9XJVna6qXwD/AGxrLx0HNg003QicaPWNQ+qSpAlaUui3OfozPgycubLnILAzyUVJrmL+C9snquok8EqS69tVO7cCDy2j35KkJRj5R1SSfA64Abg8yXHgz4EbklzL/BTNMeBjAFV1JMkB4FngNeCOqjrddnU781cCXQw80h6SlmBaf7Tm2F0fmspxtXJGhn5VfWRI+bMLtN8L7B1SnwWuWVTvJEkryl/kSlJHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6sjI0E9yX5JTSZ4ZqF2W5NEkL7TlpQOv3ZnkaJLnk9w4UL8uyeH22t1JsvLDkSQtZJxP+vcD28+q7QEOVdUW4FB7TpKtwE7g6rbNPUkuaNvcC+wGtrTH2fuUJJ1nI0O/qr4J/OSs8g5gf1vfD9w8UH+wql6tqheBo8C2JOuBS6rqsaoq4IGBbSRJE7LUOf0rq+okQFte0eobgJcG2h1vtQ1t/ez6UEl2J5lNMjs3N7fELkqSzrbSX+QOm6evBepDVdW+qpqpqpl169atWOckqXdLDf2X25QNbXmq1Y8DmwbabQROtPrGIXVJ0gQtNfQPArva+i7goYH6ziQXJbmK+S9sn2hTQK8kub5dtXPrwDaSpAm5cFSDJJ8DbgAuT3Ic+HPgLuBAktuAHwG3AFTVkSQHgGeB14A7qup029XtzF8JdDHwSHtIkiZoZOhX1UfO8dL7z9F+L7B3SH0WuGZRvZMkrSh/kStJHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0JakjFy5n4yTHgFeA08BrVTWT5DLgn4HNwDHgD6vqv1r7O4HbWvs/qqqvLOf4er3Ne7487S5IWsVW4pP+71XVtVU1057vAQ5V1RbgUHtOkq3ATuBqYDtwT5ILVuD4kqQxnY/pnR3A/ra+H7h5oP5gVb1aVS8CR4Ft5+H4kqRzWG7oF/DVJE8m2d1qV1bVSYC2vKLVNwAvDWx7vNVeJ8nuJLNJZufm5pbZRUnSGcua0wfeV1UnklwBPJrk+wu0zZBaDWtYVfuAfQAzMzND20iSFm9Zn/Sr6kRbngK+yPx0zctJ1gO05anW/DiwaWDzjcCJ5RxfkrQ4Sw79JL+c5K1n1oEPAs8AB4Fdrdku4KG2fhDYmeSiJFcBW4Anlnp8SdLiLWd650rgi0nO7Oefqupfk3wHOJDkNuBHwC0AVXUkyQHgWeA14I6qOr2s3kuaqGldEnzsrg9N5bhvREsO/ar6IfCeIfUfA+8/xzZ7gb1LPaYkaXn8Ra4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SerIhdPuwBvR5j1fnnYXJGkoQ1/SqjfND1LH7vrQ1I59Pji9I0kdMfQlqSMTD/0k25M8n+Rokj2TPr4k9WyioZ/kAuDvgT8AtgIfSbJ1kn2QpJ5N+ovcbcDRqvohQJIHgR3As+fjYF5FI2m5ppUj5+sL5EmH/gbgpYHnx4HfOrtRkt3A7vb0Z0meX6HjXw785wrta1rW+hjWev9h7Y9hrfcfOhhDPr3s/b9jWHHSoZ8htXpdoWofsG/FD57MVtXMSu93ktb6GNZ6/2Htj2Gt9x8cw3JM+ovc48CmgecbgRMT7oMkdWvSof8dYEuSq5K8BdgJHJxwHySpWxOd3qmq15J8AvgKcAFwX1UdmWAXVnzKaArW+hjWev9h7Y9hrfcfHMOSpep1U+qSpDcof5ErSR0x9CWpI2/o0E9yS5IjSX6R5JyXRiU5luRwkqeSzE6yj6MsYgyr8vYWSS5L8miSF9ry0nO0W1XnYNT7mXl3t9e/l+Q3p9HPhYwxhhuS/E97z59K8mfT6Oe5JLkvyakkz5zj9bVwDkaNYfLnoKresA/g3cBvAN8AZhZodwy4fNr9XeoYmP9S/AfAO4G3AE8DW6fd99a3vwL2tPU9wKdX+zkY5/0EbgIeYf63J9cDj0+730sYww3Al6bd1wXG8LvAbwLPnOP1VX0OxhzDxM/BG/qTflU9V1Ur9WveqRhzDP9/e4uq+l/gzO0tVoMdwP62vh+4eXpdGds47+cO4IGa923gV5Osn3RHF7Ca/zcxlqr6JvCTBZqs9nMwzhgm7g0d+otQwFeTPNluAbHWDLu9xYYp9eVsV1bVSYC2vOIc7VbTORjn/VzN7zmM37/fTvJ0kkeSXD2Zrq2Y1X4OxjXRc7Dm/3JWkq8Bbxvy0p9W1UNj7uZ9VXUiyRXAo0m+3/4LPRErMIaxbm9xvizU/0XsZqrn4CzjvJ9Tfc/HME7/vgu8o6p+luQm4F+ALee7YytotZ+DcUz8HKz50K+qD6zAPk605akkX2T+n8YTC5wVGMNUb2+xUP+TvJxkfVWdbP/0PnWOfUz1HJxlnPdztd9SZGT/quqnA+sPJ7knyeVVtVZuZLbaz8FI0zgH3U/vJPnlJG89sw58EBj6Tfsqtppvb3EQ2NXWdwGv+5fLKjwH47yfB4Fb2xUk1wP/c2Yaa5UYOYYkb0uStr6N+Tz48cR7unSr/RyMNJVzMO1vt8/nA/gw858GXgVeBr7S6m8HHm7r72T+yoangSPMT6lMve+LGUN7fhPwb8xfsbFqxgD8GnAIeKEtL1sL52DY+wl8HPh4Ww/zfxDoB8BhFrg6bBWP4RPt/X4a+DbwO9Pu81n9/xxwEvh5+//AbWvwHIwaw8TPgbdhkKSOdD+9I0k9MfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSR/4PS5R8nr0P72AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#assume the claim is true !!\n",
    "shift = mu0-empirical_mean\n",
    "cars4Shifted= cars4.mpg+shift\n",
    "bs_mean_mpg = draw_bs_reps(cars4Shifted,np.mean,10000)\n",
    "\n",
    "#The test statistic\n",
    "TS = bs_mean_mpg-mu0\n",
    "pVal = np.mean( TS < -shift)\n",
    "\n",
    "plt.hist(TS);\n",
    "\n",
    "#np.percentile(bs_mean_mpg,5)\n",
    "#p value is simply the left tail beyond xBar\n",
    "#np.mean(bs_mean_mpg < empirical_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shift = 0.7160804020100429\n",
      "pValue = 0.0375\n"
     ]
    }
   ],
   "source": [
    "print('shift =',shift)\n",
    "print(\"pValue =\", pVal)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------------------------------------\n",
    "\n",
    "### From one sample to 2 samples\n",
    "\n",
    "From Auto to birth weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "preg = pd.read_csv('data/pregNSFG.csv.gz', compression='gzip')\n",
    "\n",
    "#only look at live births\n",
    "firsts = preg[(preg.outcome == 1) & (preg.birthord == 1)]\n",
    "#live[live.babysex == 1].babysex = \"male\"\n",
    "\n",
    "#we reduce the sample size further by conditioning on \n",
    "#the mother's age at the end of pregnancy\n",
    "firsts = firsts[(firsts.agepreg < 30) & (firsts.prglngth >= 30)]\n",
    "bwt = firsts[[\"babysex\",\"totalwgt_lb\"]] \n",
    "\n",
    "bwt.babysex.replace([1.0],\"male\",inplace=True)\n",
    "bwt.babysex.replace([2.0],\"female\",inplace=True)\n",
    "bwt = bwt.dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape: (3758, 2)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>babysex</th>\n",
       "      <th>totalwgt_lb</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>male</td>\n",
       "      <td>9.1250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>male</td>\n",
       "      <td>8.5625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>male</td>\n",
       "      <td>7.5625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>male</td>\n",
       "      <td>7.8125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>female</td>\n",
       "      <td>7.0000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   babysex  totalwgt_lb\n",
       "2     male       9.1250\n",
       "5     male       8.5625\n",
       "8     male       7.5625\n",
       "10    male       7.8125\n",
       "11  female       7.0000"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('shape:',bwt.shape)\n",
    "bwt.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "babysex\n",
       "female    7.103830\n",
       "male      7.378682\n",
       "Name: totalwgt_lb, dtype: float64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grouped = bwt[\"totalwgt_lb\"].groupby(bwt[\"babysex\"])\n",
    "grouped.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzd0lEQVR4nO3de3ycdZ3o8c93ZjKZ3O9N0vQWaOkFSguWgoJIV7kpCuqy4uKqi2yXBdzj7tmzix7Z4zkc9sh5eVh11a0oHMQbKwerdbdy8wbKrUVqaUtb0jZt06bN/X6Zmcz3/PHMJJN0kkyamTyT5Pt+kVdmfs/ze+Y7pJNvftdHVBVjjDFmLI/bARhjjMlMliCMMcYkZAnCGGNMQpYgjDHGJGQJwhhjTEI+twNIpfLycl22bJnbYRhjzKzx2muvtahqRaJjcypBLFu2jJ07d7odhjHGzBoicnS8Y9bFZIwxJiFLEMYYYxKyBGGMMSahOTUGYYwxAKFQiIaGBgYGBtwOJWMEAgEWLVpEVlZW0nUsQRhj5pyGhgYKCgpYtmwZIuJ2OK5TVVpbW2loaKC2tjbpemntYhKR60TkgIjUicg9CY7fKCK7RWSXiOwUkSvijtWLyBuxY+mM0xgztwwMDFBWVmbJIUpEKCsrm3KLKm0tCBHxAl8HrgYagB0isk1V98Wd9gtgm6qqiFwI/AhYFXd8k6q2pCtGY8zcZclhtLP5/5HOFsRGoE5VD6tqEHgcuDH+BFXt0ZH9xvMA23vczG6D3fDyFmh5y+1IjJm2dCaIGuB43POGaNkoIvJBEdkP/AdwW9whBZ4RkddEZPN4LyIim6PdUzubm5tTFLoxU6OqPLHzOAce+2t46h/g0Rsg2Od2WMZFX/3qV1m9ejW33nprWq7/hS98gS996UtpuXZMOhNEovbMGS0EVd2qqquAm4D74g5drqoXA9cDd4nIlYleRFUfUtUNqrqhoiLhanFj0u57rxzji//vBc5p+CltubXQcwoObHc7LOOib3zjG2zfvp3vf//7body1tKZIBqAxXHPFwEnxztZVZ8HzhWR8ujzk9HvTcBWnC4rYzLOQGiIrzz3FrdVHSFLhri7bzOaVwH7/93t0IxL7rjjDg4fPswHPvAB7r//fm677TYuueQSLrroIn76058C8Oijj3LTTTfx/ve/n9raWr72ta/x4IMPctFFF3HZZZfR1tYGwLe+9S0uueQS1q1bx4c//GH6+s5smR46dIjrrruOt73tbbzzne9k//79KXkf6ZzmugNYISK1wAngFuBP408QkeXAoegg9cWAH2gVkTzAo6rd0cfXAP8jjbEac9a27TpJS88gNy+pI9hXwktdi2lcdgULj7wAqmCDpa767z/by76TXSm95pqFhfy3958/7vEtW7bw1FNP8atf/YoHH3yQP/qjP+KRRx6ho6ODjRs38p73vAeAPXv28PrrrzMwMMDy5ct54IEHeP311/mbv/kbHnvsMT7zmc/woQ99iL/4i78A4POf/zwPP/wwn/70p0e93ubNm9myZQsrVqzglVde4c477+SXv/zltN9n2hKEqoZF5G7gacALPKKqe0XkjujxLcCHgY+LSAjoBz4STRaVwNboqLsP+IGqPpWuWI2Zjp/tPsmyslwquvcRWXop/jd9vB5ZzsK+rdBxFEqWuR2icdEzzzzDtm3bhscLBgYGOHbsGACbNm2ioKCAgoICioqKeP/73w/A2rVr2b17N+Akkc9//vN0dHTQ09PDtddeO+r6PT09vPjii9x8883DZYODgymJPa0L5VR1O7B9TNmWuMcPAA8kqHcYWJfO2IxJhfbeIC8eauXOK2qQHW/hXf1+1vcU82zXIt4H0LDTEoTLJvpLfyaoKk8++SQrV64cVf7KK6+QnZ09/Nzj8Qw/93g8hMNhAD75yU/yk5/8hHXr1vHoo4/y61//etR1IpEIxcXF7Nq1K+Wx215MxkzDs/tOMxRRbqzpBh2CyvO5aEkJzzSXoh4fNO2b/CJmTrv22mv5l3/5F2Iz+l9//fUp1e/u7qa6uppQKJRwwLuwsJDa2lqeeOIJwElIf/jDH6YfOJYgjJmW7XsaWVyaw7mReqegci0rq/LpG/ISKlwKLQddjc+479577yUUCnHhhRdywQUXcO+9906p/n333cell17K1VdfzapVqxKe8/3vf5+HH36YdevWcf755w8PhE+XjKxTm/02bNigdsMgM1M6+oJs+J/P8akravms5zHY+X/hcyfY09jDDf/yW16ufZiq8Am46xW3Q5133nzzTVavXu12GBkn0f8XEXlNVTckOt9aEMacpWf2niYcUW64cCG01kHZcvB4Wb4gH4/AcU8NtB6CobDboRpzVixBGHOWfrb7JEvLcrmgphDaj0LJUgACWV6WlOayL1QFkZAzk8mYWcgShDFnoS06e+l9a6udLQM6jkHx0uHj51Tks6t/gfOk+YArMRozXZYgjDkLT+89xVBEed+F1dDTBOH+4RYEwOKSHHZ2FztP2utdidGY6bIEYcxZ+PfdJ6ktz2NNdaHTeoBRLYhFJbkcH8hBs/Kg8/g4VzEms1mCMGaKWnoGeelQKzdcWO3ssR8bYyiJTxA5gDCYv3AkgRgzy1iCMGaKntpziojidC/BSBdS8ZLhcxaX5gLQlW0Jwkzdr3/9a2644Qa3w7AEYcxU/cfuRs6tyGNlZYFT0HEU8irAnzd8jtOCgBZvpSUIM2tZgjBmCjr7Qrxa38Z1F1SN3MKx/eio8QeAopws8rN9nNByGOiAgdTuJmoyX319PatWreL222/nggsu4NZbb+W5557j8ssvZ8WKFbz66qu8+uqrvOMd7+Ciiy7iHe94BwcOnDnjrbe3N+F24TMhrZv1GTPX/PpgE0MR5d2rK0cKO47CwotHnSciVBZmcyxS7hR0HoeAu5vGzVs/vwdOvZHaa1atheu/OOlpdXV1PPHEEzz00ENccskl/OAHP+C3v/0t27Zt45/+6Z947LHHeP755/H5fDz33HN87nOf48knnxx1jfvvvz/hduF5eXnjvGrqWIIwZgp+c6CZsjw/6xcVOwWRIehsgPM/eMa5lYUB3uovdZ50HINKSxDzTW1tLWvXrgXg/PPP593vfjciwtq1a6mvr6ezs5NPfOITvPXWW4gIoVDojGuMt134TGwlYgnCmCn4/bF2NiwrweOJdi91nYBI+IwuJnASxN62IueJjUO4J4m/9NNlsu287733XjZt2sTWrVupr6/nqquuOuMa420XPhNsDMKYJLX1Bqlv7eOiJSUjhe1nTnGNWVCQzYHuAOrLsQRhEurs7KSmpgZwbkGayHS3C58OSxDGJGnX8XYA1i8uHilMsEguZkFhgGBYiRRUQ9e4t2M389jf//3f89nPfpbLL7+coaGhhOdMd7vw6bAuJmOStOdEFyKwtqZopLDjKCBQtPiM8ysLne6EwZxKci1BzDvLli1jz549w8/jWwjxxw4eHLlnyH333QfAVVddNdzdlJOTwze/+c30B5yAtSCMSdKh5h4WFuWQlx33d1X7USisAZ//jPMrCwMAdPsXQLclCDP7pDVBiMh1InJAROpE5J4Ex28Ukd0isktEdorIFcnWNWamHW7u5ZyKMVMLO44mHH8AqCxwEkS7rxy6GiESSXeIxqRU2hKEiHiBrwPXA2uAj4rImjGn/QJYp6rrgduAb0+hrjEzRlU50tLLuRX5ow8kWCQXsyDaxdREmXNfiL7WdIdp4sylu2Wmwtn8/0hnC2IjUKeqh1U1CDwO3Bh/gqr26EjUeYAmW9eYmdTcPUjPYHh0CyI8CN2N47YgAlleinKyODFU7BR0nUh/oAaAQCBAa2urJYkoVaW1tZVAIDCleukcpK4B4vc5bgAuHXuSiHwQ+F/AAuB9U6lrzEw53NILQG15XILoOA7oqE36xlpQkE19sNh50t0IrE9XiCbOokWLaGhooLm52e1QMkYgEGDRokVTqpPOBCEJys5I56q6FdgqIlcC9wHvSbYugIhsBjYDLFky/gfVmOk40d4PwOKS3JHCjnrn+zhdTABl+X7qg9FZT9aCmDFZWVnU1ta6Hcasl84upgYgfu7fImDcqRyq+jxwroiUT6Wuqj6kqhtUdUNFRcX0ozYmgcZOJ0FUFcU10SdYJBdTlp/Nof4cEK8zUG3MLJLOBLEDWCEitSLiB24BtsWfICLLJbolpohcDPiB1mTqGjOTGjsHKM3zE8jyjhR2HANPFhRUj1uvPM9Pc+8QFFTZYjkz66Sti0lVwyJyN/A04AUeUdW9InJH9PgW4MPAx0UkBPQDH4kOWiesm65YjZlMY+cA1UVjBvg6jkLxYvB4E1fCaUF09oeILFyIx9ZCmFkmrSupVXU7sH1M2Za4xw8ADyRb1xi3NHYOUFM8JkG010PJsgnrleU7C+iCuZUE2g9OeK4xmcZWUhuThMbO/tHjDwBtR6Bk4oHQsjxnLURv9gIbgzCzjiUIYybRHxyioy9EdVFOXGG7c6e4SVoQ5dEWRKevAoLddmc5M6tYgjBmEk3dA8DI3kqA070EUDpJCyLfaUG0esucgm5rRZjZwxKEMZNo6QkCI+MJgNO9BEmPQZzW6J3lbC2EmUUsQRgzidaeQQDK80buDjbcgpgkQRRk+/B7PRwfiiUIa0GY2cMShDGTaOtN0IJoPwJ5FZBdMGFdEaEs38+xYPQ8WwthZhFLEMZMojWaIErz4hNE/aSth5iyfD+n+wVyy+y+EGZWsQRhzCRae4LkZ/tGr6Juq590imtMWV62001VsNBaEGZWsQRhzCRaewdHtx7CQehqmFILoqUnCIV2b2ozu1iCMGYSrT3B0eMPncdBI5NOcY0pz8+mtXcQza+C7lNpitKY1LMEYcwkWnuDwyuiAWeAGpJvQeT5GQhFCOVVQW8zDIVSH6QxaWAJwphJtPYMUpaXaA1EkmMQ0cVy3VkVgELP6RRHaEx6WIIwZgKqSlvvmC6m9nrwBSC/MqlrxOq2+6KrqW0thJklLEEYM4Gu/jDhiCae4upJ7uMTW2DXTHSxnE11NbOEJQhjJtDSG11FnT9mFXWS3Usw0oJojJQ4BTZQbWYJSxDGTKA9ukiuJNaCUI1u870s6WvEWh8ngznOHehsqquZJSxBGDOBzn5nxlFRTpZT0NsCod4pJYhAlpeCbB8tvWHn9qS2o6uZJSxBGDOBWIIojiWI2BTXJNdAxDiL5Qade1NbgjCzhCUIYyZwRgtiilNcY8rzs2kdXk1tCcLMDpYgjJlAR5+TIApHtSAESpZO6Tpl+X5ae6P7MdkgtZkl0pogROQ6ETkgInUick+C47eKyO7o14sisi7uWL2IvCEiu0RkZzrjNGY8nf0hCrJ9eD3iFLQdgcIa8GVPXHGMslgLoqDKufXoYHcaojUmtXzpurCIeIGvA1cDDcAOEdmmqvviTjsCvEtV20XkeuAh4NK445tUtSVdMRozma7+0EjrAaa0zXe88jw/bX1BhvKr8ILTzVQx8b0kjHFbOlsQG4E6VT2sqkHgceDG+BNU9UVVbY8+fRlYlMZ4jJmyzv7QyPgDOF1MpcumfJ3ygmxUodu/wCmwgWozC6QzQdQAx+OeN0TLxvMp4OdxzxV4RkReE5HN41USkc0islNEdjY3N08rYGPGGpUggr3OPkpTHKAGhjf7a/XEVlNbgjCZL50JQhKUacITRTbhJIh/iCu+XFUvBq4H7hKRKxPVVdWHVHWDqm6oqKiYbszGjNLZH6I4NzZAXe98n+IUVxhZTd2k0dXUtljOzALpTBANwOK454uAMz4VInIh8G3gRlVtjZWr6sno9yZgK06XlTEzalQL4iynuAKUxxLEoA+yC20mk5kV0pkgdgArRKRWRPzALcC2+BNEZAnwY+DPVPVgXHmeiBTEHgPXAHvSGKsxCXXEJ4gp3gciXmwvp5aeYHQ1tbUgTOZL2ywmVQ2LyN3A04AXeERV94rIHdHjW4B/BMqAb4gIQFhVNwCVwNZomQ/4gao+la5YjUlkIDREMByJWwNRD4EiyC2d8rUKA1n4POLcm7qw2loQZlZIW4IAUNXtwPYxZVviHt8O3J6g3mFg3dhyY2ZSwlXUZ9G9BODxCKV5/uhaiGo48kKqwjQmbWwltTHjOCNBtB85qwHqmLL87Oh+TNXQcwoikVSEaUzaWIIwZhyjEkQkAh3HoXhqW2zEK8/309IbhMKFEAlDn60BNZnNEoQx4+jsi0sQPachEoLixZPUGp+zYV90R1ewqa4m41mCMGYcw1t952ZBZ4NTWHT2CaIsz9nyWwuqnQJbLGcynCUIY8YxqoupM7opQNHZ7wZTlp/NQChCf3Z0QaclCJPhLEEYM46OaIIoCMS3IM4+QcQWy7VQAuKx+0KYjGcJwphxdPWHKAhEt/rubHBWQAeKzvp6w4vl+ocgb4G1IEzGswRhzDhGbbPR2TCt1gOM7MfU0m23HjWzgyUIY8YxOkEcT0GCiO7oGpvqal1MJsNZgjBmHKMSRNeJ6SeIvPgWRLW1IEzGswRhzDiGt/oO9kFf67QTRCDLS2HAR3NsNXV/G4QGUhStMalnCcKYcQy3ILpOOAXTWAMRU1kY4HTXgLNhH1grwmQ0SxDGJKCqdMbuR52CNRAxToKIW01tu7qaDGYJwpgEBkIRguFIdJHc9NdAxCwozKapawAKFjoFdl8Ik8EsQRiTwKhV1LE9k2JbZExDZWGApu5BIvmxLiZrQZjMZQnCmARGJYie05BbBt6saV+3siCbcERpi+SAL2Ab9pmMZgnCmARGJ4gmyK9MyXUrCwMAnLaprmYWsARhTAJntCDyF6TkuguiCaKpa9BZLGddTCaDWYIwJoHhrb5z/NEEkaoWhLOa+nTXgDOTybqYTAZLa4IQketE5ICI1InIPQmO3yoiu6NfL4rIumTrGpNOwy2IgA96miGvIiXXrSiIJYi4LibVlFzbmFRLKkGIyJMi8j4RSTqhiIgX+DpwPbAG+KiIrBlz2hHgXap6IXAf8NAU6hqTNp19QUSgwNMP4f6UtSCyfV5K8/yc7h5wEkR4AAY6UnJtY1It2V/4/wr8KfCWiHxRRFYlUWcjUKeqh1U1CDwO3Bh/gqq+qKrt0acvA4uSrWtMOnX2hyjI9uHpbXYKUpQgABYURNdCxFZT26Z9JkMllSBU9TlVvRW4GKgHno12Cf25iIw3968GOB73vCFaNp5PAT+fal0R2SwiO0VkZ3Nz8+RvxpgkdPaHKMqNDlBDygapAaqKYqupY4vlLEGYzDSVLqMy4JPA7cDrwFdwEsaz41VJUJaws1VENuEkiH+Yal1VfUhVN6jqhoqK1PQTGzO8D9NwgkhdC6KyIDAySA2WIEzG8iVzkoj8GFgFfBd4v6rG/kX/m4jsHKdaAxC/u9ki4IwpGyJyIfBt4HpVbZ1KXWPSpaM/REmuH9LQxVRZmE1LzyDhvErnA2hdTCZDJZUggG+r6vb4AhHJVtVBVd0wTp0dwAoRqQVOALfgjGPEX2MJ8GPgz1T14FTqGpNOnf0hFhbnOC0I8UJOScquvaAwQEShZcBDVU6ptSBMxkq2i+l/Jih7aaIKqhoG7gaeBt4EfqSqe0XkDhG5I3raPwJlwDdEZFesNTJe3SRjNWbaOvtCFMcvkvOkbkb48GrqrgFbTW0y2oQtCBGpwhkczhGRixgZGygEcie7eLTVsX1M2Za4x7fjjGkkVdeYmaCqdMRuFtTSlNIBaoDqIidBNHb2s66w2hbLmYw1WRfTtTgD04uAB+PKu4HPpSkmY1zVMxhmKKJx22ykbvwB4hNEdKD61J6UXt+YVJkwQajqd4DviMiHVfXJGYrJGFeN3majCarWpvT6pXl+/D4Ppzqj94XobYKhMHiTHRI0ZmZM1sX0MVX9HrBMRP527HFVfTBBNWNmtY6+2EZ9XmcWU4pbECJCdVGAk50DsKAaNOIkicKFKX0dY6Zrsj9Z8qLf89MdiDGZItaCKPf0QiQMeakdgwCoKgxwqrN/5CZEXY2WIEzGmayL6ZvR7/99ZsIxxn2xFkQJ0V1gUjxIDbCwOIdXj7SNJIjuk8DbUv46xkxHspv1/W8RKRSRLBH5hYi0iMjH0h2cMW4YHoMIxxJEaruYILbdxoDdetRktGQnd1+jql3ADTirnM8D/kvaojLGRR39QQDyw21OQRoSxMKiAOGI0kIBeHw21dVkpGQTRGxDvvcCP1TVtjTFY4zrOvtCZPs8+AdanII0dDFVFeUA0NgZhPwqWyxnMlKyCeJnIrIf2AD8QkQqgIH0hWWMezr64jbq8wUguyDlr3HGWghLECYDJbvd9z3A24ENqhoCerH7M5g5qjO2ironuopaEm0uPD3xq6kprLYN+0xGmsrKnNU46yHi6zyW4niMcV1HfzDl96Ie64zFcoefT8vrGDMdyW73/V3gXGAXMBQtVixBmDmooy/E4tJcpwVRek5aXiO2WK6xcwAWVcFgJwR7wZ83eWVjZkiyLYgNwBpVu7u6mfu6YjcLOnUallyWttepKgw4XUxrogvkuhqhfHnaXs+YqUp2kHoPUJXOQIzJFB39IUqzBfra0tbFBM5iOWeQOrYWwsYhTGZJtgVRDuwTkVeBwVihqn4gLVEZ45JgOEJfcIgqXw+gkJe+29iOLJY71/lLzRKEyTDJJogvpDMIYzJFe5+zSK7K2+kUpLEFUV0UIDSktHrLqABLECbjJJUgVPU3IrIUWKGqz4lILuBNb2jGzLyWHqeBvMAzEwnCWSx3ss9HhT/fttswGSfZvZj+Avh/wDejRTXAT9IUkzGuae1xWhBldDgFaVhFHWOL5UymS3aQ+i7gcqALQFXfAtL3yTHGJW29ToIoHErfTq4xsQQxvO23tSBMhkk2QQyqajD2JLpYbtIpryJynYgcEJE6EbknwfFVIvKSiAyKyN+NOVYvIm+IyC4R2ZlknMZMS6yLKT/UBtlFkJWTtteKLZYbbkHYhn0mwyQ7SP0bEfkckCMiVwN3Aj+bqIKIeIGvA1fj7AC7Q0S2qeq+uNPagL8GbhrnMptUtSXJGI2ZttbeID6P4B9ohvz0zWCCMYvlyqItCNW0bO1hzNlItgVxD9AMvAH8JbAd+PwkdTYCdap6ONr6eJwx+zepapOq7gBCU4ramDRp6wlSmudHeprSOkAdM7xYrqAahgahvz3tr2lMspLdrC+CMyh9p6r+sap+K4lV1TXA8bjnDdGyZCnwjIi8JiKbxztJRDaLyE4R2dnc3DyFyxtzptbeQcrys6P7MKV/mG24BVEQXYdq4xAmg0yYIMTxBRFpAfYDB0SkWUT+MYlrJ2onT2WrjstV9WLgeuAuEbky0Umq+pCqblDVDRUV6e0SMHNfS0+Q8nx/dCfX9LcgqotzoovlYgnCZjKZzDFZC+IzOLOXLlHVMlUtBS4FLheRv5mkbgOwOO75IiDpUThVPRn93gRsxemyMiat2nqDVOYoDHbNWAsiNKS0eUudAmtBmAwyWYL4OPBRVT0SK1DVw8DHoscmsgNYISK1IuIHbgG2JROUiOSJSEHsMXANzn5QxqRVa88gi/29zpO8mUgQ0TvLDRU7BdaCMBlksllMWYlmEalqs4hkJaoQd05YRO4GnsZZdf2Iqu4VkTuix7eISBWwEygEIiLyGWANzt5PW8WZzeEDfqCqT03trRkzNf3BIXqDQ9Rk9TsFM9HFFF0LcaIH1gaKrQVhMspkCSJ4lscAUNXtODOe4su2xD0+hdP1NFYXsG6y6xuTSq290W02pMMpmIEupqr4xXKFC60FYTLKZAlinYh0JSgXIJCGeIxxTVP3zO3DFFOW58fv9dDYZdttmMwzYYJQVduQz8wbpzoHACjVDkAgrzztrykiVBUFaOyI3hei+UDaX9OYZCW7UM6YOS+WIArD7ZBbBt4Jh9lSprooEL03dZUzBhGJzMjrGjMZSxDGRJ3uGsDv85A92Dwj4w8x1UUBTsZWU+sQ9NnuMiYzWIIwJqqxc4CqwgDSPTOrqGNGFstFxzxsHMJkCEsQxkSd6hpwZhX1nB65T/QMiC2W6/RFdwKwqa4mQ1iCMCbqVOcAVQXZzi/oGZjBFFNVGJ3qqsVOgbUgTIawBGEMEIkop7oGOCdvECKhkc3zZsDCYmc19fFgPiDWgjAZwxKEMThrIILhCOfm9jgFM9mCiC6WO9kVhrwKu3GQyRiWIIwBjrY6+y8t83c7BTM4BnHmYjlrQZjMYAnCGOBoWx8AVbFV1AUz14KILZZz1kJUQ7e1IExmsARhDHCstQ+vRyiNtDoF+TM3BgGMrKYuqoHOEzP62saMxxKEMcCxtj4WFgfw9jZBdiH4c2f09RcWBWjs6oeixdDfBoM9M/r6xiRiCcIYoL61lyWludAzs1NcY6qKcjjVOUCkaIlT0Hl84grGzABLEGbeCw9FOHCqm5WVhdB9ekanuMYsLHYWy3VkR5NThyUI4z5LEGbeO9LSy2A4wvkLC51Fai4kiJroWogTGt3io+PojMdgzFiWIMy8t6/RueXJmuqC6DYbLiSIEidBHB3MA68fOo7NeAzGjGUJwsx7+0524fd6WJ4/COEBZ6B4hsVaEA0dg87r2xiEyQCWIMy8t/dkF+dV5ZPV3eAUFCW6C256FQSyKMrJ4kR7PxQvthaEyQhpTRAicp2IHBCROhG5J8HxVSLykogMisjfTaWuMamgquxr7GJNdSF0upcgABaV5NDQ3gfFS2yQ2mSEtCUIEfECXweuB9YAHxWRNWNOawP+GvjSWdQ1ZtpOdw3S1hvk/IVFcQli5ruYwOlmamjvh6Il0NsEwT5X4jAmJp0tiI1AnaoeVtUg8DhwY/wJqtqkqjuA0FTrGpMKe086W2usWRhtQWTlQk6JK7EsKsnlREc/WlrrFLQddiUOY2LSmSBqgPh2ckO0LKV1RWSziOwUkZ3Nzc1nFaiZv/addGYwra4udAaGixaDiCux1JTk0BccoitvmVPQ+pYrcRgTk84EkehTpqmuq6oPqeoGVd1QUVGRdHDGgDNAvawsl/xsn9OCcGn8AZwxCIDjEv1bqKXOtViMgfQmiAYgvjN3EZDsNpXTqWtM0vY1djndS+AsTit2Z/wBRqa6Hu8BChdZC8K4Lp0JYgewQkRqRcQP3AJsm4G6xiSlayDEsbY+Z4C6vx36WqFsuWvxLC5xNghsaO+H8hXQctC1WIwB8KXrwqoaFpG7gacBL/CIqu4VkTuix7eISBWwEygEIiLyGWCNqnYlqpuuWM389ObJ2ArqQmiNDgiXnutaPIU5PvKzfZzoiCaIXT8EVdfGRIxJW4IAUNXtwPYxZVviHp/C6T5Kqq4xqRTbYuP8hYVQf8gpdLEFISIjayEuuACC3c5MpjL3kpaZ32wltZm39p7sojzfT0VBNrTWgXigZJmrMQ2vhVi43ilo3OVmOGaeswRh5q19J7tYXV2IiEDrIWeKq8/vakyLSnKc7TYqVjub9p3c5Wo8Zn6zBGHmpfBQhLrmHmf9AzgDwi52L8UsLs2lezBM2yBQeT6cfN3tkMw8ZgnCzEv1rb0EwxFWVhbAUAia90PVBW6HRW15HuDco4KatzkJYmjsRgPGzAxLEGZe2n+qG4CVVQVO62EoCJVrXY4KlkUTRH1LL9S+C4I90LDD5ajMfGUJwsxLB0514/UIyxfkw6k3nMIq9xPE4pJcvB5xWhC1VzoD54d+5XZYZp6yBGHmpQOnullWlksgy+skCF8gI8Yg/D4Pi0pyONLaCznFULMB3nrG7bDMPGUJwsxLB053s6oqOkB9chcsWAPetC4LSlpteR5HmnudJ6tvcKa6th5yNSYzP1mCMPNOXzDMsbY+zqssgHAQTuyEJZe5HdawZWV51Lf2oqpwwR8DAm884XZYZh6yBGHmnYOne1CNDlCf2u3chzqDEsQ5FXn0BYdo6h6EohpYdgXs/pGz7YYxM8gShJl3DpxytthYVVUAx15yChdnToIYNdUVYO3N0HbI1kSYGWcJwsw7+091k5PlZUlpLhx5AUrPgYJKt8MaFksQh5p7nII1H3BWVVs3k5lhliDMvHPwdDfnVebjiQSh/gVY/h63QxqlpjiH/GwfB6JrNcgpgRXXwJ4nITLkbnBmXrEEYeadA6e6nfGHYy9BqA/OfbfbIY0iIqyqKmB/Y/dI4dqboec0HHnevcDMvGMJwswrLT2DtPQEWVlVCHXPOV03y65wO6wzrKou4M1TXc5MJoDzrgN/Aez7iatxmfnFEoSZV2LdNquqCqDuF7Dk7ZCd73JUZ1pVVUj3QJiTnQNOQVYAznmXE7PNZjIzxBKEmVdiezCtyu2Gpn0ZN/4Qs7q6AID90ZsaAU6snceh+YBLUZn5xhKEmVcOnOqiLM9P2akXnILlmTX+EHNepZMg9p4ckyDA6RozZgZYgjDzyoHTPc4A9cGnoXCRs8VGBioIZLFiQT6vH2sfKSxeDBWrLEGYGWMJwswbkYjy1ulu1izww6FfwsrrQMTtsMb1tqUl/P5YB5FI3JjDOZuc2VfhQfcCM/NGWhOEiFwnIgdEpE5E7klwXETkq9Hju0Xk4rhj9SLyhojsEpGd6YzTzA/H2/voCw5xpW+fM7115fVuhzShi5eW0Nkf4nBLz0hh7ZXO1iB2jwgzA9KWIETEC3wduB5YA3xURMa2568HVkS/NgP/Oub4JlVdr6ob0hWnmT9iA9Sru18Efz4se6fLEU1sw9ISAHbWx3UzLX2Hc48IWw9hZkA6WxAbgTpVPayqQeBx4MYx59wIPKaOl4FiEalOY0xmHnOmuCplJ34J524CX7bbIU2otjyP0jw/rxxpGynMKYbq9ZYgzIxIZ4KoAY7HPW+IliV7jgLPiMhrIrJ5vBcRkc0islNEdjY3N6cgbDNX7T/VxXuKT+HpboTzMrt7CZwV1VeuKOc3B5sZih+HqL3S6WIK9roXnJkX0pkgEo3+jV3hM9E5l6vqxTjdUHeJyJWJXkRVH1LVDaq6oaKi4uyjNXPeH4538qHc3YDAede6HU5SNq1aQFtvkN0NHSOFtVdCJDyyE60xaZLOBNEALI57vgg4mew5qhr73gRsxemyMuastPQMcqKjj3cM/AaWXg555W6HlJR3nVeBR+BX+5tGCpdcBp4s62YyaZfOBLEDWCEitSLiB24Bto05Zxvw8ehspsuATlVtFJE8ESkAEJE84BpgTxpjNXPc7oYOzpd6ivvq4cKb3Q4nacW5fi5ZVsp/vNE4si+TPw8WXWIJwqRd2hKEqoaBu4GngTeBH6nqXhG5Q0TuiJ62HTgM1AHfAu6MllcCvxWRPwCvAv+hqk+lK1Yz9+063slN3hdRTxas/oDb4UzJTRfVcKi5lzdOdI4U1l4JjX+A/vbxKxozTWldB6Gq21X1PFU9V1Xvj5ZtUdUt0ceqqndFj69V1Z3R8sOqui76dX6srjFna8+xJv4463fIimsgt9TtcKbkvWur8fs8/Pj3J0YKa68EjcDRF90LzMx5tpLazHnBcISKo9sp0Q7YeLvb4UxZUU4W71m9gJ/94SShoYhTuGgD+HKsm8mklSUIM+ftOtrKn7ONnoJzna0qZqGb37aY1t4gz+w97RT4sp3BaksQJo0sQZg5r/Xl77HKcxzPps9m9N5LE7nyvApqinP43stHRwprr3S2LO9pGr+iMdNgCcLMbT3NvL3uQeqyziN3/YfdjuaseT3Cn166hJcOt1LXFL0Vae27nO/1L7gXmJnTLEGYuSsSoeeJvyIn0sfejf8LPLP7n/tHLllMllf43svHnILqdZBdCId/425gZs6a3Z8YYyby7L3kH32WB4Zu5R2XZfbGfMkoz8/m+guqefL3DfQFw+D1OXtKHdgOQyG3wzNzkCUIMze98hC89DW+F7mWzrW3UVGQ2RvzJetjly2leyDMT3dFNyW48BbobXbub2FMilmCMHPP/u3w1D/wRv7l3Df0cf5q03K3I0qZS5aVsLq6kId/e8S5kdDy90BuGez6gduhmTnIEoSZW5oPoE/ezun81fxJy6f4q03nsXxBvttRpYyIsPnKWuqaevjVgSbw+WHtnzjdTF2Nbodn5hhLEGbuGOwm8vitdEf8fKD5r7h63Tn89R+tcDuqlLvhwoUsLArw0POHnYLL7oDIELz0NXcDM3OOJQgzN6jS9+RdaOsh7ui/k49f+3a+/JH1eDyzc93DRLK8Hm67opZXjrSx63gHlCyDC/8Ednwb2o64HZ6ZQyxBmDmh8dkvk3vwp3yVW7jzz2/jrk3L52RyiLll4xKKc7N48NmDTsG7/xE8Pvj534OOve2KMWfHEoSZ9V7/3dOU/+4+nvds5L1/+QBXrJgd93qYjvxsH3ddtZznDzbzYl0LFC6ETf8V3nrGaUkYkwKWIMysNRRRvvfsK1Q/85c0exew6o7vsbK60O2wZsyfvX0pC4sCfPGp/c6MpkvvgPOug6fugWMvux2emQN8bgdgzGR6B8N8+4Uj/Gz3Sdp7g+RmeynPz6a9o5OvDPxXir39DH1iK3kLKt0OdUYFsrz83bUr+dsf/YEfvHqMj122FD74TfjWJvjhLfDJ7VC5xu0wzSxmLQiT0eqaenjvV1/gn587SHVRgOsuqOJtS0pY4mnh2577udBzhOyPPELe0ovcDtUVH7yohsuXl/HFn+/nREc/5BTDx34M3mz47geh9ZDbIZpZzBKEyVgvHmrhQ9/4Hb2DYf5t82V891OXcv81C/ly8RN8pfl2lofrkJsfRVa9z+1QXSMifPFDFxJR5c7v/56B0BCU1sLHfwKREDx8DZz4vdthmlnKEoTJSD/aeZyPP/wqlYUBtt55OZcuCsDzX4KvrodX/tWZ1vnp1+D8m9wO1XWLS3P554+s5w/HO/jUd3bQ2ReCBavhtqfBnwv/93p49Vs2u8lMmegc+kezYcMG3blzp9thmGlo6h7g/zx9kH/beZx3rijn6zevpHDv9+F3X4Ge07Dyvc6UzgWr3Q414zyx8zif2/oGhYEsbrqohkuWlbK2eJCFv/5bpO45qF4Pb78bzrsGAkVuh2syhIi8pqobEh5LZ4IQkeuArwBe4Nuq+sUxxyV6/L1AH/BJVf19MnUTsQSR+XoGw7y4v4FDb7xMoHUv1f0HWRo6gkfD9GqA5qFcOjSf2poqNhR14TnyGwj2ODfHuepzsPTtbr+FjPZGQydffu4gL9S1EAw7tyctCvi4s+RVPtr/QwoHove1LloMxUugsMaZIlu0aPixlp2LZBe4+C7MTHIlQYiIFzgIXA00ADuAj6rqvrhz3gt8GidBXAp8RVUvTaZuIpYgUi8SUXqDYfoGBunr6aK/twuCPeTJIDl+H4GcPHJy88nKzoWsHOc+yV4fRCLQ10pf2wma6vfRUr+bUOOblPXVcQ4n8Ynzy6tX8mnIPpeQN4cc7adQeyjx9OILdkNhNSy9HNbfCksudfn/xOwyEBpi/6lu9pzoZO/JLl4/1s7BU51cLAe5Ou8Qa7NPUTHURHGomeKhFrIIj6rfQCUN/nPoKlpJZMEa8havo3rpcpZUlOL3jfRMqyq9wSHae4N09PQT7u9EB7opkH4KvYPk5+WSk1eMJ1Dg3LsiKyfxXf0iQ84fAsFeGOxxHnuzIKfE+crKnbV3A8x0EyWIdE5z3QjUqerhaBCPAzcC8b/kbwQeUydLvSwixSJSDSxLom7qfPNKCA0AI8myPximuWcQgTP6biXuvPjHsfqS8Nz4OiPHYpdOdM3RHwdNcN7oY+MdHz9GPeM8ja/j/IefMJWS/P0GwngBxUeEXJwf5hIVTnsW0FOynFOLP0jVqo34Fq4jr3gJK+2Dn3KBLC/rFxezfnHxcNnJjn5+deBCXjjYwr939DMYHiIn30u+38Mify813nZqpIXKwXoKOg+yuO8tqppfwdscgb3ONfrVT5PkMxQdvowoeBmilD4Wy+CkcYXx0Ecug5KNlzBZGsZPiGyCE9YL4mOAACHxE8RHWLJQkvt3k/x56VOWn01Bdhp/3eaUwm0/T/ll05kgaoDjcc8bcFoJk51Tk2RdAERkM7AZYMmSJWcXaflKGArGLgjAYF+Ipkhv7FXG/raO/qNL8A9PZPi4JCgbS0WGzxybhs587fHLRv2qH/UyE9WRuD/KRq4zHLkIWV4Pviwf+AvwBvLxZBeg/jz6NJuBUJjwYD9Dg30MhfqIBPvRYD9+HcTn9TAQqEDzqyheeC5LV65nYcXcX+GcyRYW53DrpUu59dKlyVcK9tFzYg9th3fR3Xycvs4WGOzCoxFEwOcR/FlZdOYU4c0phOwiyC6ghxw6w1n0Dwww1N+NDnShwW6yh3oJDPWSFRlgSHyExUeILAY9OQx4chn05DhfkoOXMHlDXeQOdZMX6cYf6cenIbI0hFfDY/7wGU9yv/aTu9bZyy/OoyDfn74XSNOYUjoTRKK0PfanMN45ydR1ClUfAh4Cp4tpKgEO+/C3zigqBhK2uYyZT/y55NduJL92o9uRGBekM0E0AIvjni8CTiZ5jj+JusYYY9IonesgdgArRKRWRPzALcC2MedsAz4ujsuATlVtTLKuMcaYNEpbC0JVwyJyN/A0zlTVR1R1r4jcET2+BdiOM4OpDmea659PVDddsRpjjDmTLZQzxph5bKJprrbVhjHGmIQsQRhjjEnIEoQxxpiELEEYY4xJaE4NUotIM3DU7TjGKAda3A4iRebKe7H3kXnmynuZje9jqapWJDowpxJEJhKRnePNEJht5sp7sfeReebKe5kr7yPGupiMMcYkZAnCGGNMQpYg0u8htwNIobnyXux9ZJ658l7myvsAbAzCGGPMOKwFYYwxJiFLEMYYYxKyBDEDROQLInJCRHZFv97rdkxTISLXicgBEakTkXvcjmc6RKReRN6I/hxmzc6OIvKIiDSJyJ64slIReVZE3op+L3EzxmSM8z5m3edDRBaLyK9E5E0R2Ssi/ylaPut+JhOxBDFz/llV10e/trsdTLJExAt8HbgeWAN8VETWuBvVtG2K/hxm03z1R4HrxpTdA/xCVVcAv4g+z3SPcub7gNn3+QgD/1lVVwOXAXdFPxez8WcyLksQZjIbgTpVPayqQeBx4EaXY5p3VPV5oG1M8Y3Ad6KPvwPcNJMxnY1x3seso6qNqvr76ONu4E2ghln4M5mIJYiZc7eI7I42sWdTs7MGOB73vCFaNlsp8IyIvCYim90OZpoqo3dgJPp9gcvxTMds/XwgIsuAi4BXmFs/E0sQqSIiz4nIngRfNwL/CpwLrAcagf/jZqxTJAnKZvPc6MtV9WKcLrO7RORKtwMys/fzISL5wJPAZ1S1y+14Ui1ttxydb1T1PcmcJyLfAv49zeGkUgOwOO75IuCkS7FMm6qejH5vEpGtOF1oz7sb1Vk7LSLVqtooItVAk9sBnQ1VPR17PJs+HyKShZMcvq+qP44Wz4mfSYy1IGZA9B9KzAeBPeOdm4F2ACtEpFZE/MAtwDaXYzorIpInIgWxx8A1zK6fxVjbgE9EH38C+KmLsZy12fj5EBEBHgbeVNUH4w7NiZ9JjK2kngEi8l2c5rMC9cBfxvopZ4PotMMvA17gEVW9392Izo6InANsjT71AT+YLe9FRH4IXIWznfRp4L8BPwF+BCwBjgE3q2pGDwCP8z6uYpZ9PkTkCuAF4A0gEi3+HM44xKz6mUzEEoQxxpiErIvJGGNMQpYgjDHGJGQJwhhjTEKWIIwxxiRkCcIYY0xCliCMMcYkZAnCGGNMQv8fjhSJ43nOLiwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tmp=grouped.plot(kind='density', legend=True);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------------\n",
    "\n",
    "#### A two-sample bootstrap hypothesis test for difference of means\n",
    "\n",
    "A one sample test compares a data set to one fixed number !\n",
    "\n",
    "We now want to compare two sets of data, both of which are samples! In particular test the hypothesis that male and female babies have the same biological weight (but not necessarily the same distribution). \n",
    "\n",
    "$H_0: \\mu_m = \\mu_f, H_A: \\mu_m \\neq \\mu_f$\n",
    "\n",
    "To do the two-sample bootstrap test, we shift both arrays to have the same mean, since we are simulating the hypothesis that their means are, in fact, equal (**equal to what value ??**). We then draw bootstrap samples out of the shifted arrays and compute the difference in means. This constitutes a bootstrap replicate, and we generate many of them. The p-value is the fraction of replicates with a difference in means greater than or equal to what was observed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOWklEQVR4nO3dUYxc1X3H8e+vEHggVIV6TRwwNa1ctfBQEm0RElJFRRsQPBgeiMxDYilITiSQEil5MMlDkJAlUpVUqlSQjEBxqxRqKUFYhbYhViQUqQEMImBDKE5wYbFlOyUt5IUW59+HvSiTZdYzuzPjmT18P9Jo7pw5597/2Wv/9vrOnetUFZKktvzWtAuQJI2f4S5JDTLcJalBhrskNchwl6QGGe6S1KCB4Z5kY5IfJHk5ycEkX+za70zyZpLnu8f1PWPuSHIoyStJrp3kBCRJH5RB17kn2QBsqKrnkpwLPAvcCHwa+GVV/fWS/pcCDwFXAB8Hvg/8YVWdHH/5kqR+Bh65V9XRqnquW34HeBm48BRDtgAPV9W7VfUacIjFoJcknSZnrqRzkk3AJ4CngKuA25N8FtgPfLmqfsFi8P+oZ9gCp/5lwLp162rTpk0rKUWSPvSeffbZn1fVXL/3hg73JB8FvgN8qareTnIfcBdQ3fM9wOeA9Bn+gXM/SbYD2wEuvvhi9u/fP2wpkiQgyX8u995QV8sk+QiLwf7tqvouQFUdq6qTVfUr4H5+feplAdjYM/wi4MjSdVbVrqqar6r5ubm+v3gkSas0zNUyAR4AXq6qb/a0b+jpdhNwoFveC2xNcnaSS4DNwNPjK1mSNMgwp2WuAj4DvJjk+a7tq8AtSS5n8ZTLYeDzAFV1MMke4CXgPeA2r5SRpNNrYLhX1Q/pfx798VOM2QnsHKEuSdII/IaqJDXIcJekBhnuktQgw12SGmS4S1KDVnT7AenDaNOOx6ay3cN33zCV7aoNhrvWhGkFrLRWeVpGkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDRoY7kk2JvlBkpeTHEzyxa79/CRPJHm1ez6vZ8wdSQ4leSXJtZOcgCTpg4Y5cn8P+HJV/TFwJXBbkkuBHcC+qtoM7Ote0723FbgMuA64N8kZkyhektTfwHCvqqNV9Vy3/A7wMnAhsAXY3XXbDdzYLW8BHq6qd6vqNeAQcMWY65YkncKKzrkn2QR8AngKuKCqjsLiLwBgfdftQuCNnmELXZsk6TQZOtyTfBT4DvClqnr7VF37tFWf9W1Psj/J/hMnTgxbhiRpCEOFe5KPsBjs366q73bNx5Js6N7fABzv2heAjT3DLwKOLF1nVe2qqvmqmp+bm1tt/ZKkPoa5WibAA8DLVfXNnrf2Atu65W3Aoz3tW5OcneQSYDPw9PhKliQNcuYQfa4CPgO8mOT5ru2rwN3AniS3Aq8DNwNU1cEke4CXWLzS5raqOjnuwiVJyxsY7lX1Q/qfRwe4ZpkxO4GdI9QlSRqB31CVpAYZ7pLUIMNdkhpkuEtSgwx3SWrQMJdCSpqCTTsem9q2D999w9S2rfHwyF2SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNGhjuSR5McjzJgZ62O5O8meT57nF9z3t3JDmU5JUk106qcEnS8oY5cv8WcF2f9r+pqsu7x+MASS4FtgKXdWPuTXLGuIqVJA1nYLhX1ZPAW0OubwvwcFW9W1WvAYeAK0aoT5K0CqOcc789yQvdaZvzurYLgTd6+ix0bZKk02i14X4f8AfA5cBR4J6uPX36Vr8VJNmeZH+S/SdOnFhlGZKkflYV7lV1rKpOVtWvgPv59amXBWBjT9eLgCPLrGNXVc1X1fzc3NxqypAkLWNV4Z5kQ8/Lm4D3r6TZC2xNcnaSS4DNwNOjlShJWqkzB3VI8hBwNbAuyQLwdeDqJJezeMrlMPB5gKo6mGQP8BLwHnBbVZ2cSOWSpGUNDPequqVP8wOn6L8T2DlKUZKk0fgNVUlqkOEuSQ0y3CWpQYa7JDVo4AeqUq9NOx6bdgmShuCRuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwaGe5IHkxxPcqCn7fwkTyR5tXs+r+e9O5IcSvJKkmsnVbgkaXnDHLl/C7huSdsOYF9VbQb2da9JcimwFbisG3NvkjPGVq0kaSgDw72qngTeWtK8BdjdLe8Gbuxpf7iq3q2q14BDwBXjKVWSNKzVnnO/oKqOAnTP67v2C4E3evotdG2SpNNo3B+opk9b9e2YbE+yP8n+EydOjLkMSfpwW224H0uyAaB7Pt61LwAbe/pdBBzpt4Kq2lVV81U1Pzc3t8oyJEn9rDbc9wLbuuVtwKM97VuTnJ3kEmAz8PRoJUqSVurMQR2SPARcDaxLsgB8Hbgb2JPkVuB14GaAqjqYZA/wEvAecFtVnZxQ7ZKkZQwM96q6ZZm3rlmm/05g5yhFSZJG4zdUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSggTcOk/Ths2nHY1PZ7uG7b5jKdlvkkbskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapC3/F2DpnU7Vklrh0fuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUEjXQqZ5DDwDnASeK+q5pOcD/wTsAk4DHy6qn4xWpmSpJUYx5H7n1fV5VU1373eAeyrqs3Avu61JOk0msRpmS3A7m55N3DjBLYhSTqFUcO9gO8leTbJ9q7tgqo6CtA9rx9xG5KkFRr19gNXVdWRJOuBJ5L8ZNiB3S+D7QAXX3zxiGVIknqNdOReVUe65+PAI8AVwLEkGwC65+PLjN1VVfNVNT83NzdKGZKkJVYd7knOSXLu+8vAp4ADwF5gW9dtG/DoqEVKklZmlNMyFwCPJHl/Pf9YVf+a5BlgT5JbgdeBm0cvU5K0EqsO96r6GfAnfdr/C7hmlKIkSaPxG6qS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAaN+h9kS9LYbNrx2FS2e/juG6ay3UnyyF2SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAZ5nfsIpnVNriQN4pG7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yG+oSvrQm+a3zSf1v0BN7Mg9yXVJXklyKMmOSW1HkvRBEzlyT3IG8HfAXwILwDNJ9lbVS5PYnvd4kaTfNKkj9yuAQ1X1s6r6X+BhYMuEtiVJWmJS4X4h8EbP64WuTZJ0GkzqA9X0aavf6JBsB7Z3L3+Z5JVl1rUO+PkYa5slzm1tcm5r00zOLd8YafjvLffGpMJ9AdjY8/oi4Ehvh6raBewatKIk+6tqfrzlzQbntjY5t7Wp5bn1M6nTMs8Am5NckuQsYCuwd0LbkiQtMZEj96p6L8ntwL8BZwAPVtXBSWxLkvRBE/sSU1U9Djw+hlUNPHWzhjm3tcm5rU0tz+0DUlWDe0mS1hTvLSNJDZqJcE9yfpInkrzaPZ+3TL8HkxxPcmBJ+51J3kzyfPe4/vRUPtgY5jbU+GlYwdz63opi1vbboFtmZNHfdu+/kOSTw46dthHndjjJi90+2n96Kx9siLn9UZJ/T/Jukq+sZOyaVlVTfwB/BezolncA31im358BnwQOLGm/E/jKtOcxobkNNX5W58biB+o/BX4fOAv4MXDprO23U9XZ0+d64F9Y/B7HlcBTw45dq3Pr3jsMrJv2PEaY23rgT4GdvX/eZn2/jfqYiSN3Fm9NsLtb3g3c2K9TVT0JvHWaahqXUec21PgpGaa2tXIrimHq3AL8fS36EfA7STYMOXaaRpnbrBs4t6o6XlXPAP+30rFr2ayE+wVVdRSge16/inXc3v1z8sFZOnXB6HMbx89mUoapbdCtKGZlvw1zy4zl+sz67TZGmRssfrv8e0me7b5ZPktG+dnP+n4byWm7n3uS7wMf6/PW18aw+vuAu1j8Q3gXcA/wuTGsdygTnttUjWFup7oVxVT32xIDb5lxij7DjJ2mUeYGcFVVHUmyHngiyU+6f2nOglF+9rO+30Zy2sK9qv5iufeSHEuyoaqOdv8UPL7CdR/rWdf9wD+vvtKVm+TcgFHHj2QMc1v2VhTT3m9LDLxlxin6nDXE2GkaZW5U1fvPx5M8wuLpjFkJ92HmNomxM29WTsvsBbZ1y9uAR1cyeMm5wZuAA8v1nYKR5jaG8ZM0TG3L3opixvbbMLfM2At8truy5Ergf7rTUbN+u41Vzy3JOUnOBUhyDvApZuvv1yg/+1nfb6OZ9ie6tfip9e8C+4BXu+fzu/aPA4/39HsIOMriByMLwK1d+z8ALwIvsLhzNkx7TmOcW9/xs/BYwdyuB/6DxSsTvtbTPlP7rV+dwBeAL3TLYfE/oflpV/f8oDnOymO1c2PxSpIfd4+Da3RuH+v+Tr0N/He3/NtrYb+N8vAbqpLUoFk5LSNJGiPDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBv0/ysjuNl5cpNQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "meanNull = np.mean(bwt.totalwgt_lb)# pooled mean\n",
    "w_m = bwt[bwt[\"babysex\"]==\"male\"].totalwgt_lb\n",
    "w_f = bwt[bwt[\"babysex\"]==\"female\"].totalwgt_lb\n",
    "empirical_diff_means = np.mean(w_m)-np.mean(w_f) \n",
    "#shift:\n",
    "\n",
    "w_m_shifted = w_m - np.mean(w_m) + meanNull\n",
    "w_f_shifted = w_f - np.mean(w_f) + meanNull\n",
    "\n",
    "\n",
    "# Compute 10,000 bootstrap replicates from shifted arrays\n",
    "M=1000\n",
    "bs_replicates_m = draw_bs_reps(w_m_shifted, np.mean, M)\n",
    "bs_replicates_f = draw_bs_reps(w_f_shifted, np.mean, M)\n",
    "\n",
    "# Get replicates of difference of means: bs_replicates\n",
    "bs_replicates = bs_replicates_m - bs_replicates_f\n",
    "\n",
    "tmp=plt.hist(bs_replicates)\n",
    "# Compute and print p-value: p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The one sided p value of 0.0 is much smaller than alpha= 0.01 , so we fail to reject the Null\n",
      "-> The observed difference of  0.275 is exremely unlikely to have occurred by chance alone\n"
     ]
    }
   ],
   "source": [
    "#debugging:\n",
    "if 0:\n",
    "    meanNull\n",
    "    #plt.hist(w_f_shifted)\n",
    "    #plt.hist(bs_replicates_m)\n",
    "    #bs_replicates_m\n",
    "    w_m_bs = np.random.choice(w_m, size=len(w_m))\n",
    "    np.nanmean(w_m_bs)\n",
    "    np.argwhere(np.isnan(w_m_bs))\n",
    "    np.argwhere(np.isnan(w_m))\n",
    "\n",
    "#p-value (one-sided):\n",
    "pVal = np.mean(bs_replicates> empirical_diff_means)\n",
    "#cutoff right tail\n",
    "#np.percentile(bs_replicates, 95)\n",
    "print(\"The one sided p value of\", np.round(pVal,3), \"is much smaller than alpha=\",a,\", so we fail to reject the Null\" ) \n",
    "\n",
    "print(\"-> The observed difference of \", np.round(empirical_diff_means,5), \"is exremely unlikely to have occurred by chance alone\")"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3.9.0 64-bit",
   "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.9.0"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "949777d72b0d2535278d3dc13498b2535136f6dfe0678499012e853ee9abcab1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}