{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "remove_input"
    ]
   },
   "outputs": [],
   "source": [
    "path_data = '../../../../data/'\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "from scipy import stats\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "remove_input"
    ]
   },
   "outputs": [],
   "source": [
    "def r_scatter(r):\n",
    "    plt.figure(figsize=(5,5))\n",
    "    \"Generate a scatter plot with a correlation approximately r\"\n",
    "    x = np.random.normal(0, 1, 1000)\n",
    "    z = np.random.normal(0, 1, 1000)\n",
    "    y = r*x + (np.sqrt(1-r**2))*z\n",
    "    plt.scatter(x, y)\n",
    "    plt.xlim(-4, 4)\n",
    "    plt.ylim(-4, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlation ###\n",
    "\n",
    "In this section we will develop a measure of how tightly clustered a scatter diagram is about a straight line. Formally, this is called measuring *linear association*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table `hybrid` contains data on hybrid passenger cars sold in the United States from 1997 to 2013. The data were adapted from the online data archive of [Prof. Larry Winner](http://www.stat.ufl.edu/%7Ewinner/) of the University of Florida. The columns:\n",
    "\n",
    "- `vehicle`: model of the car\n",
    "- `year`: year of manufacture\n",
    "- `msrp`: manufacturer's suggested retail price in 2013 dollars\n",
    "- `acceleration`: acceleration rate in km per hour per second\n",
    "- `mpg`: fuel econonmy in miles per gallon\n",
    "- `class`: the model's class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "hybrid = pd.read_csv(path_data + 'hybrid.csv')"
   ]
  },
  {
   "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>vehicle</th>\n",
       "      <th>year</th>\n",
       "      <th>msrp</th>\n",
       "      <th>acceleration</th>\n",
       "      <th>mpg</th>\n",
       "      <th>class</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Prius (1st Gen)</td>\n",
       "      <td>1997</td>\n",
       "      <td>24509.74</td>\n",
       "      <td>7.46</td>\n",
       "      <td>41.26</td>\n",
       "      <td>Compact</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Tino</td>\n",
       "      <td>2000</td>\n",
       "      <td>35354.97</td>\n",
       "      <td>8.20</td>\n",
       "      <td>54.10</td>\n",
       "      <td>Compact</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Prius (2nd Gen)</td>\n",
       "      <td>2000</td>\n",
       "      <td>26832.25</td>\n",
       "      <td>7.97</td>\n",
       "      <td>45.23</td>\n",
       "      <td>Compact</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Insight</td>\n",
       "      <td>2000</td>\n",
       "      <td>18936.41</td>\n",
       "      <td>9.52</td>\n",
       "      <td>53.00</td>\n",
       "      <td>Two Seater</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Civic (1st Gen)</td>\n",
       "      <td>2001</td>\n",
       "      <td>25833.38</td>\n",
       "      <td>7.04</td>\n",
       "      <td>47.04</td>\n",
       "      <td>Compact</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>148</th>\n",
       "      <td>S400</td>\n",
       "      <td>2013</td>\n",
       "      <td>92350.00</td>\n",
       "      <td>13.89</td>\n",
       "      <td>21.00</td>\n",
       "      <td>Large</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>149</th>\n",
       "      <td>Prius Plug-in</td>\n",
       "      <td>2013</td>\n",
       "      <td>32000.00</td>\n",
       "      <td>9.17</td>\n",
       "      <td>50.00</td>\n",
       "      <td>Midsize</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>C-Max Energi Plug-in</td>\n",
       "      <td>2013</td>\n",
       "      <td>32950.00</td>\n",
       "      <td>11.76</td>\n",
       "      <td>43.00</td>\n",
       "      <td>Midsize</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>151</th>\n",
       "      <td>Fusion Energi Plug-in</td>\n",
       "      <td>2013</td>\n",
       "      <td>38700.00</td>\n",
       "      <td>11.76</td>\n",
       "      <td>43.00</td>\n",
       "      <td>Midsize</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>152</th>\n",
       "      <td>Chevrolet Volt</td>\n",
       "      <td>2013</td>\n",
       "      <td>39145.00</td>\n",
       "      <td>11.11</td>\n",
       "      <td>37.00</td>\n",
       "      <td>Compact</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>153 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                   vehicle  year      msrp  acceleration    mpg       class\n",
       "0          Prius (1st Gen)  1997  24509.74          7.46  41.26     Compact\n",
       "1                     Tino  2000  35354.97          8.20  54.10     Compact\n",
       "2          Prius (2nd Gen)  2000  26832.25          7.97  45.23     Compact\n",
       "3                  Insight  2000  18936.41          9.52  53.00  Two Seater\n",
       "4          Civic (1st Gen)  2001  25833.38          7.04  47.04     Compact\n",
       "..                     ...   ...       ...           ...    ...         ...\n",
       "148                   S400  2013  92350.00         13.89  21.00       Large\n",
       "149          Prius Plug-in  2013  32000.00          9.17  50.00     Midsize\n",
       "150   C-Max Energi Plug-in  2013  32950.00         11.76  43.00     Midsize\n",
       "151  Fusion Energi Plug-in  2013  38700.00         11.76  43.00     Midsize\n",
       "152         Chevrolet Volt  2013  39145.00         11.11  37.00     Compact\n",
       "\n",
       "[153 rows x 6 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hybrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The graph below is a scatter plot of `msrp` *versus* `acceleration`. That means `msrp` is plotted on the vertical axis and `accelaration` on the horizontal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdQAAAEfCAYAAAD4CDcbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABFvElEQVR4nO3de1iUdd748fc4IBKIGCJICC2KqGRZlpLtesJUMjOV54Ha1R5S8ZDrYQUPmbi4rqc8bpEaurtlJxOtxXKtyxWPGbRtCr8yosfVkFUQchSU4zC/P3i4ZYYBZnBg7sHP67q8Lue+PzPz/c4M85nv8dbodDoDQgghhLgj7exdACGEEKItkIQqhBBC2IAkVCGEEMIGJKEKIYQQNiAJVQghhLABSahCCCGEDUhCFUIIIWxAEqoQQghhA5JQ7SwnJ8feRbgjjl5+cPw6OHr5wfHrIOW3PzXUwa4J9dSpU0RHR9OnTx88PT159913lXOVlZWsWLGCwYMH4+fnR0hICNOmTSM3N9foMcrLy4mPjycoKAg/Pz+io6PJy8szitHpdMTGxhIQEEBAQACxsbHodDqjmNzcXKKiovDz8yMoKIhFixZRUVFhFPPtt9/y1FNP4evrS58+fVi3bh0Gg2w0JYQQws4J9ebNm/Tt25e1a9fi6upqdO7WrVucPXuWuLg4jh07xnvvvUdeXh6RkZFUVVUpcUuXLuXAgQPs2rWLgwcPUlxcTFRUFHq9XomZNm0amZmZ7N27l5SUFDIzM5kxY4ZyXq/XExUVRUlJCQcPHmTXrl2kpqaybNkyJebGjRtMmDCBrl27cuTIEdauXctrr73G66+/3oKvkBBCCEfhZM8nHzVqFKNGjQJg9uzZRuc6derExx9/bHRs8+bNhIWFkZ2dTWhoKNevX2f37t0kJSUxfPhwAHbs2EG/fv04evQo4eHhZGdnc/jwYQ4dOsSgQYOUx4mIiCAnJ4fg4GCOHDnCuXPnyMrKwt/fH4DExETmzp3L8uXL8fDwYO/evZSWlrJt2zZcXV3p27cvP/zwA2+88QZz5sxBo9G08KslhBBCzRxqDLW4uBgAT09PAM6cOUNlZSUjRoxQYvz9/QkJCSE9PR2AjIwM3N3dlWQKEBYWhpubm1FMSEiIkkwBwsPDKS8v58yZM0rM448/btSSDg8P5/Lly1y8eLFF6iuEEMJxOExCraio4JVXXmHMmDHcd999ABQUFKDVavHy8jKK9fb2pqCgQInx8vIyakFqNBq6dOliFOPt7W30GF5eXmi12kZjam/XxgghhLh72bXL11JVVVXExsZy/fp13n///SbjDQZDvQTanBjT46YxtROSGuvutWTmmRpmp90JRy8/OH4dHL384Ph1kPLbX2vUITg4uMFzqk+oVVVVTJ06le+++45PPvmEe++9VznXtWtX9Ho9RUVFdOnSRTleWFjI4MGDlZjCwkKjBGowGCgqKlJamF27dlW6f2sVFRWh1+uNYkxbooWFhQD1Wq51NfbiA8o4rqNy9PKD49fB0csPjl8HNZX/wo1KYo/rKCzX08VFS/JQTwI7Ojd6HzWVv7nUUAdVd/lWVlYSExPDt99+y4EDB/Dx8TE6379/f5ydnUlLS1OO5eXlkZ2drYyZDhw4kJKSEjIyMpSYjIwMbt68aRSTnZ1ttNwmLS0NFxcX+vfvr8ScPn2asrIyo5hu3boRGBho87oLIURzxB7XkXG1gvM39GRcrWD6MZ29i3TXsGtCLSkpITMzk8zMTKqrq7l06RKZmZnk5uZSVVXFCy+8wD//+U927tyJRqMhPz+f/Px8SktLgZqZwJMnTyYhIYGjR49y9uxZZsyYQWhoKMOGDQMgJCSEkSNHsmDBAr766isyMjJYsGABo0ePVn7NjBgxgj59+jBz5kzOnj3L0aNHSUhIYMqUKXh4eAAQGRmJq6srs2fP5rvvviM1NZUtW7Ywe/ZsmeErhFCNwnJ9o7dFy7FrQv3mm28YMmQIQ4YMobS0lDVr1jBkyBBWr15NXl4eBw8e5PLlywwbNoyQkBDl3/79+5XHWL16NU8//TQxMTGMGTMGNzc3PvjgA7RarRKTnJzMAw88wMSJE5k0aRIPPPAAO3bsUM5rtVr27NnDPffcw5gxY4iJieHpp59m1apVSkynTp346KOPuHz5MsOHDyc+Pp6XXnqJOXPmtM6LJYQQFujiom30tmg5Gp1OJ1v92JEa+v3vhKOXHxy/Do5efnD8Oqip/BeLK5l+TMZQ7UH1k5KEEEJYLrCjM58/3fBESdFyVD0pSQghhHAUklCFEEIIG5CEKoQQQtiAJFQhhBDCBiShCiGEEDYgCVUIIYSwAUmoQgghhA1IQhVCCCFsQBKqEEIIYQOSUIUQQggbkIQqhBBC2IAkVCGEEMIGJKEKIYQQNiAJVQghhLABSahCCCGEDUhCFUIIIWxAEqoQQghhA5JQhRBCCBuQhCqEEELYgCRUIYQQwgYkoQohhBA2IAlVCCGEsAEnexdACCGEaGkXblQSe1xHYbmeLi5akod6EtjR2abPIS1UIYQQbV7scR0ZVys4f0NPxtUKph/T2fw5JKEKIYRo8wrL9Y3etgXp8hVCCKFoja5Re+jiouU8eqPbtiYtVCGEEIrW6Bq1h+Shngz0bk+Qh5aB3u1JHupp8+eQFqoQQghFa3SN2kNgR2c+f9q7RZ/Dri3UU6dOER0dTZ8+ffD09OTdd981Om8wGFizZg29e/fG19eXsWPHcu7cOaOY8vJy4uPjCQoKws/Pj+joaPLy8oxidDodsbGxBAQEEBAQQGxsLDqdzigmNzeXqKgo/Pz8CAoKYtGiRVRUVBjFfPvttzz11FP4+vrSp08f1q1bh8FgsN0LIoQQdmbaFdoSXaNtlV0T6s2bN+nbty9r167F1dW13vmtW7eSlJTEunXrOHLkCN7e3kyYMIHi4mIlZunSpRw4cIBdu3Zx8OBBiouLiYqKQq+//atq2rRpZGZmsnfvXlJSUsjMzGTGjBnKeb1eT1RUFCUlJRw8eJBdu3aRmprKsmXLlJgbN24wYcIEunbtypEjR1i7di2vvfYar7/+egu9OkII0fpao2u0rbJrl++oUaMYNWoUALNnzzY6ZzAY2LZtG/Pnz2f8+PEAbNu2jeDgYFJSUoiJieH69evs3r2bpKQkhg8fDsCOHTvo168fR48eJTw8nOzsbA4fPsyhQ4cYNGgQAJs3byYiIoKcnByCg4M5cuQI586dIysrC39/fwASExOZO3cuy5cvx8PDg71791JaWsq2bdtwdXWlb9++/PDDD7zxxhvMmTMHjUbTWi+bEEK0mNboGm2rVDsp6eLFi+Tn5zNixAjlmKurK4MHDyY9PR2AM2fOUFlZaRTj7+9PSEiIEpORkYG7u7uSTAHCwsJwc3MzigkJCVGSKUB4eDjl5eWcOXNGiXn88ceNWtLh4eFcvnyZixcv2v4FEEII4VBUOykpPz8fAG9v419K3t7eXL58GYCCggK0Wi1eXl71YgoKCpQYLy8voxakRqOhS5cuRjGmz+Pl5YVWqzWK8fPzq/c8tefuv/9+s/XIyclpsq6WxKiZo5cfHL8Ojl5+cPw6SPntrzXqEBwc3OA51SbUWqZdqQaDocnuVdMYc/GWxJgeN1eWxu4Ljb/4gNLt7Kgcvfzg+HVw9PKD49dBym9/aqiDart8fXx8AJQWYq3CwkKlZdi1a1f0ej1FRUWNxhQWFhrNxjUYDBQVFRnFmD5PUVERer2+0ZjCwkKgfitaCCHE3Ue1CTUwMBAfHx/S0tKUY2VlZZw+fVoZD+3fvz/Ozs5GMXl5eWRnZysxAwcOpKSkhIyMDCUmIyODmzdvGsVkZ2cbLbdJS0vDxcWF/v37KzGnT5+mrKzMKKZbt24EBgba/gUQQgjhUOyaUEtKSsjMzCQzM5Pq6mouXbpEZmYmubm5aDQaZs2axZYtW0hNTeW7775j9uzZuLm5ERkZCUCnTp2YPHkyCQkJHD16lLNnzzJjxgxCQ0MZNmwYACEhIYwcOZIFCxbw1VdfkZGRwYIFCxg9erTSPTBixAj69OnDzJkzOXv2LEePHiUhIYEpU6bg4eEBQGRkJK6ursyePZvvvvuO1NRUtmzZwuzZs2WGrxDC4Vy4UcmoT67yyL4rvHjWhYvFlfYuksOz6xjqN998w7hx45Tba9asYc2aNTz33HNs27aNefPmUVpaSnx8PDqdjgEDBrB//346duyo3Gf16tVotVpiYmIoKytjyJAhbN++Ha329mLk5ORkFi9ezMSJEwGIiIhg/fr1ynmtVsuePXuIi4tjzJgxdOjQgcjISFatWqXEdOrUiY8++oi4uDiGDx+Op6cnL730EnPmzGnJl0gIIVpE7RaDNbRMP6aT5TJ3SKPT6WSrHztSw0D6nXD08oPj18HRym9u8/WKKxccqg6mHO09AHhk3xXO37i9AU6Qh5Z/TfJt8n5q3TxfDe+BasdQhRBtU1vdfN3RNHeLQXn/Gqb6ZTNCiLalrW6+rhaWtiCTh3oy/VhNnJuhwuItBuX9a5gkVCFEq2qN61LezeqOjZ5H3+DYaN0tBnNycizutpX3r2HS5SuEaFWy+XrLaukWpLx/DZMWqhCiVZnbfD3nip0K0wa1dAtSNs9vmLRQhRCiDZEWpP1IC1UIIdoQaUHajyRUIUSrMjcLVYi2QLp8hRCtStYxirZKEqoQolXJOkbRVklCFUK0qubu0COE2klCFUK0KpmFKtoqmZQkhGhVsg5VtFXSQhVCCCFsQBKqEEIIYQPS5SuEMEut170UQq2khSqEMEvWiwphHUmoQgizZL2oENaRhCqEMEvWiwphHUmoQgizZL2oENaRSUlCCLPkqiVCWEdaqEIIIYQNSEIVQgghbEC6fIUQQihk/XHzSQtVCCGEQtYfN58kVCGEEApZf9x8klCFEEIoZP1x80lCFUIIoZD1x80nk5KEEEIoZP1x86m6harX61m1ahUPPvggPj4+PPjgg6xatYqqqiolxmAwsGbNGnr37o2vry9jx47l3LlzRo9TXl5OfHw8QUFB+Pn5ER0dTV5enlGMTqcjNjaWgIAAAgICiI2NRafTGcXk5uYSFRWFn58fQUFBLFq0iIqKiharvxBCCMeh6oS6ZcsWdu7cybp168jIyGDt2rUkJyezadMmJWbr1q0kJSWxbt06jhw5gre3NxMmTKC4uFiJWbp0KQcOHGDXrl0cPHiQ4uJioqKi0OtvD7ZPmzaNzMxM9u7dS0pKCpmZmcyYMUM5r9friYqKoqSkhIMHD7Jr1y5SU1NZtmxZ67wYQgghVE3VXb4ZGRmMGTOGiIgIAAIDA4mIiODrr78Galqn27ZtY/78+YwfPx6Abdu2ERwcTEpKCjExMVy/fp3du3eTlJTE8OHDAdixYwf9+vXj6NGjhIeHk52dzeHDhzl06BCDBg0CYPPmzURERJCTk0NwcDBHjhzh3LlzZGVl4e/vD0BiYiJz585l+fLleHh4tPbLI4QQQkVU3UINCwvj5MmT/PDDDwB8//33nDhxgieffBKAixcvkp+fz4gRI5T7uLq6MnjwYNLT0wE4c+YMlZWVRjH+/v6EhIQoMRkZGbi7uyvJtPa53dzcjGJCQkKUZAoQHh5OeXk5Z86caZkXQAghhMNQdQt1/vz5lJSUMGjQILRaLVVVVcTFxTFt2jQA8vPzAfD2Nh5A9/b25vLlywAUFBSg1Wrx8vKqF1NQUKDEeHl5odFolPMajYYuXboYxZg+j5eXF1qtVokxJycnp8l6WhKjZi1d/kulGhJ+aI+uUoOns4E/9KrgPleDTZ9D3gP7c/Q6SPntrzXqEBwc3OA5VSfU/fv388EHH7Bz50569+5NVlYWS5YsISAggClTpihxdRMh1HQFmx4zZRpjLt6SmMaOQ+MvPqB0KTuq1ij/S59cJau4ZvJXbhms/qmTTWchyntgf45eBym//amhDqru8k1ISGDOnDlMmjSJ0NBQoqOjeemll9i8eTMAPj4+APVaiIWFhUprsmvXruj1eoqKihqNKSwsxGC43eoxGAwUFRUZxZg+T1FREXq9vl7LVdiW7NwihHAEqk6ot27dQqs13qVDq9VSXV0N1ExS8vHxIS0tTTlfVlbG6dOnlfHQ/v374+zsbBSTl5dHdna2EjNw4EBKSkrIyMhQYjIyMrh586ZRTHZ2ttFym7S0NFxcXOjfv79tKy6MyM4tQghHoOou3zFjxrBlyxYCAwPp3bs3mZmZJCUlER0dDdR0tc6aNYuNGzcSHBxMz5492bBhA25ubkRGRgLQqVMnJk+eTEJCAt7e3nTu3Jlly5YRGhrKsGHDAAgJCWHkyJEsWLCArVu3YjAYWLBgAaNHj1a6EEaMGEGfPn2YOXMmq1at4tq1ayQkJDBlyhSZ4dvCkod6Mv2Y8dUvhBBCbVSdUNevX88f//hHFi5cSGFhIT4+PrzwwgssWrRIiZk3bx6lpaXEx8ej0+kYMGAA+/fvp2PHjkrM6tWr0Wq1xMTEUFZWxpAhQ9i+fbtR6zc5OZnFixczceJEACIiIli/fr1yXqvVsmfPHuLi4hgzZgwdOnQgMjKSVatWtcIrcXeTnVuEEI5Ao9PpbDtdUlhFDQPpd8LRyw+OXwdHLz84fh2k/PanhjqoegxVCCGEcBSSUIUQQggbUPUYqhB3qws3Kok9bjwRK7Cjs72LJYRohLRQhVCh2OM6Mq5WcP6GnoyrFUw/prN3kYQQTZCEKoQKyWYWQjieZnX5VlVV8d577/H555+Tm5sLQPfu3Rk1ahTPPfcczs7SNSXEnejiouU8eqPbwjzpHhdqYXULNT8/n6FDhzJv3jxOnjwJ1GzTd/LkSebNm8fQoUOVTeuFEM2TPNSTgd7tCfLQMtC7vWxm0QjpHhdqYXULddGiReTk5PDaa6/x3HPPKZsj6PV63nvvPRYuXMiiRYt46623bF5YIe4WspmF5aR7XKiF1Qn1H//4BzNmzOA3v/mN0XGtVsvkyZP5/vvvefvtt21WQCGEaIx0jwu1sLrL18XFhe7duzd4PjAwEBcXlzsqlBBCWEq6x4VaWN1CnThxIvv27SMmJqbe5KOKigr27dvHhAkTbFZAIYRojHSPC7WwOqE+88wzfPHFFwwfPpwXX3yRoKAgNBoNP/74I3/5y18AGD9+PF9//bXR/QYMGGCbEgshhANpy7OQ23LdmqNZCbXWwoUL0Wg0AEYX564bYzAY0Gg0/Pzzz3dSTiGEcEi1s5ABzqNn+jFdm2lRt+W6NYfVCfX1119XkqgQQojGteVZyG25bs1hdUL99a9/3RLlEEKINqktz0Juy3VrDqtm+ZaWlnLvvfeycePGliqPEEK0KW1tFvKFG5WM+uQqj+y7Qrm+mgc7O7WZut0pq1qorq6ueHt707Fjx5YqjxBCtCltbRZy3XFTgIHe7Tn+bNup352weh3qhAkT+Oijj6iurm6J8gghhFAxGTdtmNVjqGPHjuX48eOMGTOGKVOmcP/99+Pq6lovTpbJCCFE2yPjpg27o2UzX331Vb0Zv7JMRggh2q7koZ5MP2a89lTUsDqhJiUltUQ5hBBCOIC2NiZsS1Yn1Oeff74lyiGEEEI4tGZdYNycjIwMdDodTzzxBG5ubrZ6WCGEEColWw8as3qW7/r16+ttfh8VFcWYMWOIiopi4MCB/PTTTzYroBBCtKS66ypHfXKVi8WV9i6Sw5CLuxuzOqF+/PHH9O3bV7l98OBBPv/8c+bNm8euXbuoqKhg/fr1Ni2kEEK0FEkKzSdLaIxZ3eV76dIlgoODldsHDhygR48erFixAoCcnBzeeecd25VQCCFakCSF5pMlNMaaNYaq199+AY8dO8a4ceOU235+fly9evXOSybaPBl/EWpgq6RwN36eZQmNMasTas+ePfn000+ZOnUqhw8f5sqVK4wcOVI5n5eXh6enpy3LKNooufSTY2srCcRWSaE5n2dHfw1lCY0xqxPqb3/7W6ZOnUpgYCC3bt2iV69eDB8+XDl/7Ngx+vXrZ9NCirZJutrUrakv+7byg8hWSaE5n+e28hqKGs3ay3f//v08//zz/O53vyM1NRUnp5q8fO3aNby8vJg8ebLNCnjlyhVmzpxJjx498PHxYdCgQZw8eVI5bzAYWLNmDb1798bX15exY8dy7tw5o8coLy8nPj6eoKAg/Pz8iI6OJi8vzyhGp9MRGxtLQEAAAQEBxMbGotPpjGJyc3OJiorCz8+PoKAgFi1aREVFBaJ5TLvW7vbxF7VparKO/CAy1pzPc0u8hjJr2X6aNYY6bNgwhg0bVu94586dbTohSafTMXr0aMLCwvjwww/x8vLi4sWLeHvf/gW3detWkpKSSEpKIjg4WFnW89VXXylXxVm6dCkHDx5k165ddO7cmWXLlhEVFcWxY8fQams+9NOmTePSpUvs3bsXjUbD3LlzmTFjBnv27AFqxo2joqLo3LkzBw8e5Nq1a8yaNQuDwcCrr75qszrfTWT8xX4s6Wps6steJqQYa87nuSVeQ2n12k+zEurBgwfZvXs3Fy5cQKfTYTAYjM5rNJp6rcTm+NOf/oSvry87duxQjt1///3K/w0GA9u2bWP+/PmMHz8egG3bthEcHExKSgoxMTFcv36d3bt3k5SUpHRN79ixg379+nH06FHCw8PJzs7m8OHDHDp0iEGDBgGwefNmIiIiyMnJITg4mCNHjnDu3DmysrLw9/cHIDExkblz57J8+XI8PDzuuL53Gxl/sR9LvnSb+rKXH0TGmvN5bonXUHoO7MfqhLpu3TrWrVtHp06deOCBBwgKCmqJcgHw6aefEh4eTkxMDCdOnMDX15cpU6Ywffp0NBoNFy9eJD8/nxEjRij3cXV1ZfDgwaSnpxMTE8OZM2eorKw0ivH39yckJIT09HTCw8PJyMjA3d1dSaYAYWFhuLm5kZ6eTnBwMBkZGYSEhCjJFCA8PJzy8nLOnDnDkCFDWux1EMLWLPnSberLXn4Q3TmTtohNSM+B/VidUJOTkxk6dCgffPABLi4uLVEmxYULF9i1axezZ89m/vz5ZGVlsXjxYgBiY2PJz88HMOoCrr19+fJlAAoKCtBqtXh5edWLKSgoUGK8vLyMrpyj0Wjo0qWLUYzp83h5eaHVapUYc3JycpqspyUxaubo5QfHr4O15XczuADaOrcrzD5GUsjt/1dc0ZFzpbklbNrd9h4AvHjWhazimvfhPHomf3aZXQ+VA3CpVEPCD+3RVWrwdDbwh14V3OfadAZ+OUBDQll7rv3f/V4OuE5Ojq5Fyq82rVGHuvswmLI6oVZWVvLMM8+0eDIFqK6u5uGHH1Y2jXjooYc4f/48O3fuJDY2Volr6BJyjTGNMRdvSUxjx6HxFx9QupQdlaOXHxy/Ds0p/zu+lUatz4QB7rz0dYndlm/cje8BwM3MK1CnNVmiaU9wcAAAL31ylazimm753DJY/VMni3oEgoHjD1pXDkd//UEddbB6lu+IESP45ptvWqIs9fj4+BASEmJ0rFevXly6dEk5D9RrIRYWFiqtya5du6LX6ykqKmo0prCw0Ggs2GAwUFRUZBRj+jxFRUXo9fp6LVch1K62u/Zfk3z5/GlvVn5dItvv2UFjM4NlLNTxWJ1QX331Vb755hvWrl1Lbm5uvQlJthQWFsaPP/5odOzHH3+ke/fuAAQGBuLj40NaWppyvqysjNOnTyvjof3798fZ2dkoJi8vj+zsbCVm4MCBlJSUkJGRocRkZGRw8+ZNo5js7Gyj5TZpaWm4uLjQv39/21ZciFYmX972kTzUk4He7Qny0DLQu73ROLUsK3M8Vnf5dunShUmTJrFy5coGN8HXaDT1WoTNMXv2bEaNGsWGDRuYOHEimZmZvPnmmyxfvlx5nlmzZrFx40aCg4Pp2bMnGzZswM3NjcjISAA6derE5MmTSUhIwNvbW1k2Exoaqiz9CQkJYeTIkSxYsICtW7diMBhYsGABo0ePVroQRowYQZ8+fZg5cyarVq3i2rVrJCQkMGXKFJnhK1TH2h14ZCKLfTQ2sUtmUTseqxPq73//e/70pz8RGBjIgAEDWjSZPPLII7z77rusXLmSV199FX9/f15++WWmTZumxMybN4/S0lLi4+PR6XQMGDCA/fv3K2tQAVavXo1WqyUmJoaysjKGDBnC9u3blTWoUDPZavHixUycOBGAiIgIox8MWq2WPXv2EBcXx5gxY+jQoQORkZGsWrWqxeovRHNZuxZRvrzVR2ZROx6NTqezqs82KCiIwYMHyxVlbEQNA+l3wtHLD45fB3Plf2TfFc7fuN3iDPLQ8q9Jvq1dNIvZ+j1o7T1y2+JnyNGooQ5Wj6FWV1cTHh7eEmURQtjI3T7+Jtc4FfZgdUKNiIgw2ktXCKE+jU12sZQj7wkrk6wcnyN+/qweQ124cCEvvvgi8+bNY/Lkyfj7+xuNRdaSpSRC2I8txt8ceU9YmWTl+Bzx82d1Qn3ssccAyMrKYvfu3Q3G/fzzz80vlRDCLHNjgy3FkVt5MslK3WxxcQY1sjqhLlq0qMldiIQQLcP0V/sTH1/lnYc0tMRUDEdu5ckMWXWzxcUZ1MjqhLp06dKWKIcQwgKmv9JLqgwk/NDe6q3mLCGtPNFSbHFxBjVq1uXbhBD2YfqrHeBaZcv0GFnbymvtpSrCcVnS+nTEXgarZ/kKIewneagn7k7GCdTTueW2/7SGLFURlrLFLHQ1khaqEA4ksKMzp571NuoKezngur2LBTjmJBJhH47Y+rSEJFQhHIzpl5El17psDfaaRCJdzUItpMtXCGET9urGk65moRbSQhVC2IS9uvHU0tUsLWUhLVQhhEMzt2+xPbatk5aykIQqhHBo5rqa7ZHc1NJSFvYjXb5CCIdmrqvZHsnNEXf2EbYlLVQhRJtjj8vXtdW1lcJy0kIVQrQ59ti2rq2urRSWk4QqhGhzJLkJe5CEKoQwu+TDYMCqZSCybETc7SShCiHMXk4LsOoCz454QWghbEkSqhDColmxTc2UlWUj4m4ns3yFsCN7bEBgjrlZsdbOlDU9X3Cr2m71EcIeJKEKYUdq2V3H3JIPa5eBmF5arqTKILsFibuKdPkKYUdq6SZtaFasNWOggR2d6XpPO0pu3K6DdPuKu4m0UIWwI3tsQNCS2lp9hLCGtFCFsIG6S0bcnTRggBK9ocnlI/bYgKAlWVIfc8trhGgLJKEKYQN1l4zU1dTykba2AYEl9TG3vCYppDVKJ0TLki5fIWygsbFCGUc0ppZxYyFszaES6saNG/H09CQ+Pl45ZjAYWLNmDb1798bX15exY8dy7tw5o/uVl5cTHx9PUFAQfn5+REdHk5eXZxSj0+mIjY0lICCAgIAAYmNj0el0RjG5ublERUXh5+dHUFAQixYtoqKifqtEOIbaJSsT/9nhjpesNDZWaM9xRLUsy6lLxllFW+UwCfWrr77irbfeIjQ01Oj41q1bSUpKYt26dRw5cgRvb28mTJhAcXGxErN06VIOHDjArl27OHjwIMXFxURFRaHX3/5lPG3aNDIzM9m7dy8pKSlkZmYyY8YM5bxerycqKoqSkhIOHjzIrl27SE1NZdmyZS1feQekxi9yU7Vdj7ll7e54yUrdJSYP3evEg52dVHHVEbUsy6lLrsoi2iqHGEO9fv0606dP57XXXmP9+vXKcYPBwLZt25g/fz7jx48HYNu2bQQHB5OSkkJMTAzXr19n9+7dJCUlMXz4cAB27NhBv379OHr0KOHh4WRnZ3P48GEOHTrEoEGDANi8eTMRERHk5OQQHBzMkSNHOHfuHFlZWfj7+wOQmJjI3LlzWb58OR4eHq38qqibI2xDZ8uuR7WOhaqxe9VgsHcJhGgZDtFCrU2YQ4cONTp+8eJF8vPzGTFihHLM1dWVwYMHk56eDsCZM2eorKw0ivH39yckJESJycjIwN3dXUmmAGFhYbi5uRnFhISEKMkUIDw8nPLycs6cOWPzOjs6NX6Rm7obuh7VWEc1tpqFsAXVJ9S33nqL8+fPm+1azc/PB8Db27hl4O3tTUFBAQAFBQVotVq8vLwajfHy8kKjub3Li0ajoUuXLkYxps/j5eWFVqtVYsRtavwiN1Xb9di9Q3WLdz2evFyK/+7/0OWvefjv/g+nLpeajbN1V3nCAHfcnTQ4acDdScOKAe539Hi24Ag/toRoDlV3+ebk5LBy5Ur+/ve/0759+wbj6iZCqOkKNj1myjTGXLwlMY0dh5o6NMWSGDUzV/6XAzQklLXnWqUGT2cDLwdcJydH1/qFa8Lt5RplVFzRkXOlZZ7nv75wpbS65nNSUmXgvz4v4tjg+kn1xbMuZBXX/Pg4j57Jn10msVcFCT+0R/d/r+UfelVwn6txv2lDn6GlZ10oqdIqz7v01FV2PVRuy6pZzc3gAmjr3K4ZGmiLfweOxNHLD61Th+Dg4AbPqTqhZmRkUFRUxOOPP64c0+v1fPHFF/z5z3/myy+/BGpaj3W7YgsLC5XWZNeuXdHr9RQVFdGlSxejmMGDBysxhYWFRgnUYDBQVFRk9Di13b+1ioqK0Ov19VqudTX24gPKGK2jaqj8wcDxB1u/PM3RGu9B5SnjWeUVBo3Z57z+zWWgWrmtMziz4t8uZBVXAZBbBiv+3ZFj432UmMbKfzPzCnC7BViiaU9wcMAd1OTOveNbWW/zh4orF9rk34GjcPTygzrqoOou37Fjx/LFF19w4sQJ5d/DDz/MpEmTOHHiBD179sTHx4e0tDTlPmVlZZw+fVoZD+3fvz/Ozs5GMXl5eWRnZysxAwcOpKSkhIyMDCUmIyODmzdvGsVkZ2cbLbdJS0vDxcWF/v37t+TL4JAcYZZva+qg1TR6u9bP5cYtz5/LDGTrqoyOfW9yuzFq7HqvncD1r0m+fP60t1yEXLQZqm6henp64unpaXTsnnvuoXPnzvTt2xeAWbNmsXHjRoKDg+nZsycbNmzAzc2NyMhIADp16sTkyZNJSEjA29ubzp07s2zZMkJDQxk2bBgAISEhjBw5kgULFrB161YMBgMLFixg9OjRyi+eESNG0KdPH2bOnMmqVau4du0aCQkJTJkyRWb4muEIs3xb056RnYk6fI0yvYEOWg17RnY2G+fpAiVVxrevlpkENT6aYcRRtja8VKrhpU+uGpVTEq1wNKpOqJaYN28epaWlxMfHo9PpGDBgAPv376djx45KzOrVq9FqtcTExFBWVsaQIUPYvn07Wu3tX+vJycksXryYiRMnAhAREWG0REer1bJnzx7i4uIYM2YMHTp0IDIyklWrVrVeZR1IW5h4Ym7P2eZ+yT/RzZWT452Ux0v8uoTkoU71Hs/vHmcu3awwun1v+2oyr93OsiEelv/ZqnU5j6mEH9qTVSw/wIRj0+h0OlkVZkdq6Pe/Ew2Vf9QnV432th3o3V61X5CtVQdLHu9icf3xRaDesbqJ2NE/QwD93s8lt+z2CFSQh5Z/TfK1Y4ms4+jvgaOXH9RRB4dvoQp1ak5Xo6UtQlu2HBtj61a2JY9ni+uSOpLa9zK/3LgfWw1jvUJYSxKqsJk7TXSWjru21vhsFxct5+vMkL3TL3lbP15bcPu9rEmoLlp46F7ZjlA4JkmowmbuNNGZttjOFlXwyL4r9ZJza43P2npCT1OPt/d/i5lx/AbV1Ey/Tx7iwaQeHVutRW4Ppu/dfW7aNtsaF22fJFRhM3ea6ExbcOXVcP6Gvl5ybq2Wnq0n9DT1eLXJFGpWok4/foNJPTq26RnT0moXbYmq16EKx3Knax7rXoXE9K51k3NbvVpJdQO328KM6Ya05vaPQrQ0aaEKm7nTLtK6LTjTGbF1k7OjLAWxVjuMk2rtr9223IqrfS9rZmh2t3dxhLgjklCFzdgy0TnKhgSWsmQcNHmIB9NNxlChZoP76DqbQqhhg3shRH2SUIUqmEs4bakVask46KQeHZnUo2O9+678uoSSqprl4iVVBhK/LuHzp11bvtBCCKtIQhWqYEnCac3ZrrZ+rv/cqjR7++TlUqPW556RnXmim3GybMtjqEK0JTIpSaiCJUmjNS9Mbevn0pWbvx19+BolVQaqDDWtz6jD1+rdV40b3Ash6pOEKlTBkqRh65Za7RVxJv6zQ70r4tj6ue7toDF7u0xvvPOn6W24s1nNctUfIVqPdPkKVWhsElJt92veTeOkZklLrbGu29vdzO3ILato0bWuvq5O/FRSYXQbai7jVjs+CtC+Xc0M5//cqkRXXpN4PZ3b0cg17BvVltewCqE2klCFKjQ2Q7huUgBwaQcPeVnWUmssoTTWCm3pXZISBrgz6pOruLQzUPJ/Me0Anw4Y1bWkxMBPdRbTWJsUZfxViNYjCVWoXr3t6dwt356usYTSWCu0pXdJMl1nCzVrUP9jeu1TM6xJipa2tK2dhNWWt0MUorkkoQrVM00K7loNoyy8GHVjCaW21Xi5uIxuHTu06lrXBpOiBRdTbKz72TTRrRjgTuLXJU22tK3tGpauZCHqk4QqVM+0u7RcX23xl3ljXbf23KXHNNHX6u3phIu2ndkx1OIqQ5Pdz6aJrmbNatOJztquYTV1JUtrWaiFJFQVki8IY6bdpY/su2J0vrEvc7VuU1ib6I0SZ/t2GAw19fG7x5lPI5p+3y/cqOTFsy7czKy5Ks+V0iqj81duVVnUmrd2EpZpfMGtai4WV9rlcyqtZaEWklBVSL4gGtfQl78j/RAxl+jrjqvWvu9vDvFstE6xx3VkFWuBmqvyuDsZTwf+udzATzeb/ixZOwkreagnT3x81WgHJ0s/p+bepzuhptayuLtJQlUhR/iCqP1SvFzcgW7ZV1s1eTX05d8aP0RaMmmbe9+bqpPpfTxdoG/n9kr5/nOrkpIqGoyvZW1LPrCjM13vaUfJjduPZ+nn1FydkkIsfup62vLFA4RjkYSqQvb+grAkaTS2hrOlNfTl3xo/RFoyaZt735uqk+l9/O6pP5v40s3bs4kLblWbvWi7JUw/Fx1NWsNdXLT1YhIGuLPSZFKUrd+ntnYhBeG4JKGqkK2+IJrbmrIkaaixFd0aP0Rast5133d3rYZyfTV5JcaPX3CrmpOXS5Uk1dFJQ8g9eiqd2pv9rNR9zIJb1ZRUGSgxc9F2S5h+Lh7s7MRA7/ZGn6/px4xjardWrL09/ZjO5u+TWsfJxd1HEqoK2eoLormtKUuShr1b0ea0RkulJevd2PVga5VUGYySFMCDHSF9kq/ZxzTUWYZTaTBek2PtjwHT+BK9gePPNv5Dy3QrxcJyPR+N8qr3PlVc0VlVFiHUSBJqG9bc1pQlSaM113Ba2tKuTUi1V3AZsK+gwSu4NFdL9h4YDCjHTLdZrMs0SV2rbHhfQtNdpuqyZj0vWPa5MI0x3Vqxi4vW7A/GHOOJ20I4JEmobVhzd8mxZDOA1lzDaW1Lu24LrvYKLpcmW5ZQm7qcWkv2HgANJr+6TJOUp3PDu0GY/ohyaVez05S163nBsh8TpjGWbiwhRFsgCbUNs7Q11dzNAFqLpS3t2h8GdZMNmL+CS0PuJBlbw5I61Sa/jk4aDIaaLlZzSerlgOsNPo/pj6qHvNor760163nBsh8T5mLkYujibiEJtQ2ztDWlxglGdZnbRMDcTNWGujc7aC2/VIsll1OzhYZ6DxpKfqbqJqmcHF2Dz9PYjyo1joML4cgkoQrVf7E2NlP1iY+v0vWedsqaS1PuTjXdtpYy7U61JBk3ZzZ1Q4nO1pOqGvtRJctNhLAtSahC9V+sdZPCI/uuGG0mUDe5mu4SNNC74RZeQ/aM7EyUyRhqU0y7zOsm+aYmUJlqza52WW4ihG1JQhUO9cXa0KbyULMfbl9X5zv6YXCfmxN9O99+DH/3pv9E6i0nuYO1nkIIx9XO3gVozKZNmxg+fDjdu3enR48eREVF8d133xnFGAwG1qxZQ+/evfH19WXs2LGcO3fOKKa8vJz4+HiCgoLw8/MjOjqavLw8oxidTkdsbCwBAQEEBAQQGxuLTqczisnNzSUqKgo/Pz+CgoJYtGgRFRVNz8oUtpMwwB13Jw1OmvofXl9XJz5/2pt/TfLl86e9m7UlYG1r8/wNPRlXK5TZt41prItcbePRQoiWo+qEevLkSaZOncpnn31GamoqTk5OPPvss1y7dk2J2bp1K0lJSaxbt44jR47g7e3NhAkTKC4uVmKWLl3KgQMH2LVrFwcPHqS4uJioqCj0+ttfdtOmTSMzM5O9e/eSkpJCZmYmM2bMUM7r9XqioqIoKSnh4MGD7Nq1i9TUVJYtW9Y6L4YAYOXXJZRUGagy1FyQ291JQ5CHloHe7W3SVd2cCVrJQz0Z6N2eIA9tvW5ntY1HCyFajqq7fPfv3290e8eOHQQEBPDll18SERGBwWBg27ZtzJ8/n/HjxwOwbds2goODSUlJISYmhuvXr7N7926SkpIYPny48jj9+vXj6NGjhIeHk52dzeHDhzl06BCDBg0CYPPmzURERPzfOstgjhw5wrlz58jKysLf3x+AxMRE5s6dy/Lly/Hw8GjFV+buZZrgut7Tjn81sEuQOU1d6aQ5E7TqdplfLK5U9Xi0EKLlqLqFaqqkpITq6mo8PT0BuHjxIvn5+YwYMUKJcXV1ZfDgwaSnpwNw5swZKisrjWL8/f0JCQlRYjIyMnB3d1eSKUBYWBhubm5GMSEhIUoyBQgPD6e8vJwzZ860VJVFHRduVFJwq9romLUtwKa6dOu2NpvT6q1NrnfS7SyEcEyqbqGaWrJkCf369WPgwIEA5OfnA+DtbTzpw9vbm8uXLwNQUFCAVqvFy8urXkxBQYES4+XlhUZzu7tOo9HQpUsXoxjT5/Hy8kKr1Sox5uTk5DRZL0ti1Ky55b9UqiHhh/boKjV4Ohv4Q68K7nNteN3ni2ddKKm6nUDbYyCzqByvv1yifTuI9innr5ddMAAaYGXPcrq4wO/OuVBRDe3bgbvWQN3fkZeLy+rVoe6lxCqu6BxiW7za8lv7mqrJ3fp3oBaOXn5onToEBwc3eM5hEurLL7/Ml19+yaFDh9BqjVsldRMh1ExUMj1myjTGXLwlMY0dh8ZffEDpUnZUd1L+lz65SlZxzaSu3DJY/VOnRmfE3sy8AnW6YyvRYPi/BmtpNfzlcgflnAFY8WMH7nHSUFptUGLKq43fq24dOwBlqn4PmlrnWvc9sPY1bWmWrtG9m/8O1MDRyw/qqINDdPkuXbqUffv2kZqayv33368c9/HxAajXQiwsLFRak127dkWv11NUVNRoTGFhIYY6V+MwGAwUFRUZxZg+T1FREXq9vl7L9W534UYloz65yiP7rjDqk6tcLK6/4QJYPwHItHu3qXZXNfV3OtJouKMuXXuwZuax2na9as6saSEcleoT6uLFi0lJSSE1NZVevXoZnQsMDMTHx4e0tDTlWFlZGadPn1bGQ/v374+zs7NRTF5eHtnZ2UrMwIEDKSkpISMjQ4nJyMjg5s2bRjHZ2dlGy23S0tJwcXGhf//+Nq+3I7P0S9Q0QTY1Hmo6vnmPBcOnpjsduWo1DjfGaU2StPY1bWlqS/BCtCRVd/nGxcWxZ88e3nnnHTw9PZUxUzc3N9zd3dFoNMyaNYuNGzcSHBxMz5492bBhA25ubkRGRgLQqVMnJk+eTEJCAt7e3nTu3Jlly5YRGhrKsGHDAAgJCWHkyJEsWLCArVu3YjAYWLBgAaNHj1a6EEaMGEGfPn2YOXMmq1at4tq1ayQkJDBlyhSZ4WvC0i9Ra3doMt2A4tTlUqJMrg1al5OmeTsfqU1jM48v3KjkxbMu3My8osqru6h9W0shbEnVCXXnzp0AypKYWosXL2bp0qUAzJs3j9LSUuLj49HpdAwYMID9+/fTsWNHJX716tVotVpiYmIoKytjyJAhbN++3WgsNjk5mcWLFzNx4kQAIiIiWL9+vXJeq9WyZ88e4uLiGDNmDB06dCAyMpJVq1a1WP0dlbtJq9D0dq073aHpiW6uXJrsqixVOVNUQUWdScChnZ2UGEfW2A+P2OM6soq1gF6VVwpS+7aWQtiSRqfTOcYUwDZKDQPptZqzybu58g/9Wz5nf65Sbj90rxPHxvu0SJnrMrcG1JIuXTW9B9Z6ZN8VztfZ2zjIQ2vVuly1sPV70JzP8p1w5M8QOH75QR11UHULVbQuay/k3ZBiky5Y09stxZH2JLYV6VI1z1afZSGsofpJSaL12GoCidomxrRlyUM9ebCj3qFmLbcGmQwl7EFaqHcJS7rAbNXakXGz1hPY0ZldD5UTHBxg76KoirTchT1IQr1LWNIFZqtE2FTXa2uPb4m7j/yoE/YgCfUuYUkXWGuNQcr4lmhpd+N4urA/GUO9S6hpXFPGt4QQbZEk1LvEnV5FxZbUlNyFEMJWpMv3LqGmLjAZ3xJCtEWSUIXFmro4t6XUlNyFEMJWpMtXWEyuHCKEEA2ThCosJpOJhBCiYZJQhcVkMpEQQjRMEqqwmJpmCgshhNrIpCRhMXOTiXKu2KkwQgihMtJCFUIIIWxAEqoQQghhA5JQhRBCCBuQhCqEEELYgCRUIYQQwgY0Op3OYO9CCCGEEI5OWqhCCCGEDUhCFUIIIWxAEqoQQghhA5JQhRBCCBuQhCqEEELYgCRUO7ly5QozZ86kR48e+Pj4MGjQIE6ePGnvYllEr9ezatUqHnzwQXx8fHjwwQdZtWoVVVVV9i6aWadOnSI6Opo+ffrg6enJu+++a3TeYDCwZs0aevfuja+vL2PHjuXcuXN2Kq15jdWhsrKSFStWMHjwYPz8/AgJCWHatGnk5ubascTGmnoP6po3bx6enp689tprrVjCpllShx9//JHf/OY3BAQE0K1bN4YMGUJ2drYdSltfU+UvKSkhPj6evn374uvry6OPPkpSUpKdSlvfpk2bGD58ON27d6dHjx5ERUXx3XffGcXY+29ZEqod6HQ6Ro8ejcFg4MMPPyQ9PZ3169fj7e3d9J1VYMuWLezcuZN169aRkZHB2rVrSU5OZtOmTfYumlk3b96kb9++rF27FldX13rnt27dSlJSEuvWrePIkSN4e3szYcIEiouL7VBa8xqrw61btzh79ixxcXEcO3aM9957j7y8PCIjI1XzI6ep96DW3/72N/71r3/RrVu3ViydZZqqw4ULFxg9ejSBgYGkpqZy+vRpXnnlFdzc3OxQ2vqaKv+yZcv4/PPP2b59O+np6SxcuJDExEQ++OADO5S2vpMnTzJ16lQ+++wzUlNTcXJy4tlnn+XatWtKjL3/lmUdqh2sXLmSU6dO8dlnn9m7KM0SFRVF586d2b59u3Js5syZXLt2jT179tixZE277777WL9+Pb/+9a+Bml+0vXv3Zvr06cTFxQFQWlpKcHAwf/jDH4iJibFncc0yrYM533//PWFhYZw6dYrQ0NBWLF3TGir/Tz/9xOjRo/n444+JjIwkNjaW3/72t3YqZePM1WHatGloNBqSk5PtWDLLmCv/448/zrhx43j55ZeVY0899RShoaG8+uqr9ihmo0pKSggICODdd98lIiJCFX/L0kK1g08//ZQBAwYQExNDz549+eUvf8mbb76JweAYv23CwsI4efIkP/zwA1Dz5X3ixAmefPJJO5fMehcvXiQ/P58RI0Yox1xdXRk8eDDp6el2LNmdqf1F7unpad+CWKiqqopp06YRFxdHSEiIvYtjterqag4dOkRISAiTJk2iR48eDB8+nP3799u7aBYLCwvj0KFDXLp0CYD09HT+3//7f4SHh9u5ZOaVlJRQXV2tfMbV8Lcs10O1gwsXLrBr1y5mz57N/PnzycrKYvHixQDExsbauXRNmz9/PiUlJQwaNAitVktVVRVxcXFMmzbN3kWzWn5+PkC97nZvb28uX75sjyLdsYqKCl555RXGjBnDfffdZ+/iWGTNmjV07tyZqVOn2rsozXL16lVKSkrYtGkTL7/8MitWrOD48eNMnz6de+65hzFjxti7iE1at24dCxYs4IEHHsDJqSY1rF+/XrVlX7JkCf369WPgwIGAOv6WJaHaQXV1NQ8//DArVqwA4KGHHuL8+fPs3LnTIRLq/v37+eCDD9i5cye9e/cmKyuLJUuWEBAQwJQpU+xdvGbRaDRGtw0GQ71jjqCqqorY2FiuX7/O+++/b+/iWOTkyZO89957nDhxwt5Fabbq6mqgpot0zpw5ADz44IOcOXOGnTt3qjYp1bVjxw7S09N5//336d69O1988QXLly8nICCAkSNH2rt4Rl5++WW+/PJLDh06hFarNTpnz79lSah24OPjU69bq1evXkpXi9olJCQwZ84cJk2aBEBoaCi5ubls3rzZ4RKqj48PAAUFBfj7+yvHCwsLHWaSWK2qqiqmTp3Kd999xyeffMK9995r7yJZ5MSJE1y5csXob0Kv17NixQq2bdtWbyanGnl5eeHk5GT279oRun1LS0tZuXIlf/3rX4mIiADggQceICsri9dee01VCXXp0qXs37+fAwcOcP/99yvH1fC3LGOodhAWFsaPP/5odOzHH3+ke/fudiqRdW7dulXvV6FWq1V+pTuSwMBAfHx8SEtLU46VlZVx+vRpBg0aZMeSWaeyspKYmBi+/fZbDhw4oHy5OIJp06Zx6tQpTpw4ofzr1q0bs2fP5m9/+5u9i2eR9u3b88gjj5CTk2N03FH+risrK6msrFT93/XixYtJSUkhNTWVXr16GZ1Tw9+ytFDtYPbs2YwaNYoNGzYwceJEMjMzefPNN1m+fLm9i2aRMWPGsGXLFgIDA+nduzeZmZkkJSURHR1t76KZVVJSwvnz54GarrlLly6RmZlJ586d6d69O7NmzWLjxo0EBwfTs2dPNmzYgJubG5GRkXYu+W2N1aFbt2688MILfPPNN7z//vtoNBplPMnDw6PRZSqtpan3wLQF4eTkhI+PD8HBwfYorllN1WHu3LnExMQwePBghgwZwokTJ9i/f3+ja25bU1Plf+KJJ0hMTMTNzY3u3btz6tQpPvjgAxITE+1c8hpxcXHs2bOHd955B09PT+Uz7ubmhru7OxqNxu5/y7Jsxk4+++wzVq5cyY8//oi/vz/Tp09nxowZDjFuV1xczB//+Ec++eQTCgsL8fHxYdKkSSxatIgOHTrYu3j1nDhxgnHjxtU7/txzz7Ft2zYMBgNr167lr3/9KzqdjgEDBrBhwwb69u1rh9Ka11gdlixZwkMPPWT2fklJSY0ur2ktTb0Hpvr166e6ZTOW1OHdd99l06ZN5OXlERQUxO9+9zvV/DBrqvz5+fkkJiaSlpbGtWvX6N69O1OmTGHOnDmq+F5qaMb64sWLWbp0KYDd/5YloQohhBA2IGOoQgghhA1IQhVCCCFsQBKqEEIIYQOSUIUQQggbkIQqhBBC2IAkVCGEEMIGJKEKIepZs2aNaq9Uo+ayibubJFQhhOrk5eWxZs0aMjMz7V0UISwmCVUIoTr/+c9/WLduHVlZWfXOxcfHc+XKFTuUSojGSUIVQrS4W7du2eyxnJycVLnFpRCSUIWwk59++omFCxfy2GOP0a1bNwICAoiKiuLcuXP1YisqKnj11Vd57LHH6Nq1K8HBwTz33HNGsQaDgeTkZH75y1/i6+tLUFAQzz77LF988YXRY+3bt4/w8HCj5/z+++8tKnNaWhpPP/00/v7++Pn58fTTT5Oenm4UUzvG+f333zNz5kx+8YtfEBYWZnGdT5w4wZNPPgnASy+9hKenJ56enqxZs8bo8U29/fbbDB48GB8fH3r27MmMGTPqXVh61qxZ+Pj4UFBQQExMDN27dycwMJB58+ZRVlZm0WsgREPkajNC2Mk333zDqVOnGDduHAEBAVy+fJm//OUvPPXUU3z55ZfKJdiqq6t57rnn+Mc//sEzzzzD9OnTKS0t5cSJE5w5c4Y+ffoAMG/ePN5++22GDRvG888/j8FgICMjg9OnTzN48GAAtmzZwu9//3vGjRtHdHQ0N2/eZOfOnYwePZpjx44ZXV/SVEpKCrGxsfzqV79i2bJlVFdX8+677/LMM8/w6aef8uijjxrFx8TEEBAQwLJly6ioqLC4ziEhISxZsoS1a9fyP//zPzz++ONAzXV3G7J582YSExMZPHgwK1eu5NKlSyQnJ3P69GmOHz9ulICrq6uZMGECoaGhJCYm8s9//pO33noLLy8vEhISrH4fhaglm+MLYSe3bt3innvuMTp2/vx5Hn/8ceLj44mLiwNqrmDy0ksv8corryjHahkMBjQajXIlkRdeeIGtW7eajcnNzeXhhx9m4cKFytU5AK5cucLAgQN55plneP3114GaVuC6devQ6XQA3Lx5k9DQUCIiIoyuDnPr1i3CwsK4//77SU1NNbrvuHHj2L17d7Pq/NVXX/Hkk0+avVqOadmKioro27cvjz32GB9//DFOTjXthE8//ZRf//rXxMXF8corrwA1LdT333/f6BjA888/T3p6Ov/7v/9b730SwlLS5SuEndRNLLdu3eLnn3+mU6dO9OjRgzNnzijnUlNT6dSpk9lLmdVeVqs2mdVNEqYxBw4coKqqikmTJlFUVKT8c3Z25tFHH+X48eMNljUtLQ2dTsd///d/G923tLSUYcOGcfr0aSorK43uM3Xq1GbX2RpHjx6lvLyc2bNnK8kUYOzYsQQHB/PZZ5/Vu49p2Z544gmKioooLi5uVhmEAOnyFcJuysrKWL16NR9++GG9WateXl7K///973/Ts2dPXFxcGnysf//733h7e9e7UHddta2vgQMHmj1v2nI0d98JEyY0GHP9+nW6dOmi3DbXfWxpna3x008/AdCrV69653r16sXJkyeNjjk7O9OtWzejY7VdwteuXaNjx47NKocQklCFsJMlS5bw9ttvExsbS1hYGB4eHrRr146lS5dSXV2txNV22TbGkpjax0xJSTFqydVq167hDqva+77xxhv4+fmZjfHw8DC67erqWi/G0jrbisFQf0SrsXqaixfCUpJQhbCT/fv3Ex0dzdq1a42O63Q67r33XuV2UFAQ6enpVFRU0L59e7OPFRQUxD/+8Q+uXr3aYCv1F7/4BQD+/v707t3bqrLW3rdLly4MGzbMqvvWZWmdm/pxUFdAQAAAP/zwAz179jQ6l5OTo5wXoqXJGKoQdqLVauu1iFJSUuot9XjmmWfQ6XQkJSXVe4za+z/zzDMArF69utEYJycn1qxZY7Y1WFhY2GBZw8PD6dSpExs2bKC8vNyq+9ZlaZ1ru59rJx41ZtiwYbi4uLB9+3b0er1y/O9//zs5OTmMHj3aorIJcaekhSqEnURERPDBBx/QsWNH+vbtS1ZWFvv376839hgdHc2HH35IYmIiZ8+e5YknnqCsrIyTJ08yYcIEoqOj+dWvfsXzzz/PX/7yFy5cuMCoUaOAmtmyoaGhLFy4kPvvv5/ExESWLVvGyJEjGTduHJ07dyY3N5fPP/+cRx99lM2bN5sta8eOHdm6dStTp07ll7/8Jf/1X/+Fj48PeXl5nDhxAjc3N1JSUmxW5x49euDh4cGf//xn3N3dcXd3p0+fPvTt27feY3p5ebFkyRISExMZP34848aNIy8vjzfffJOAgADmzJlj2RsixB2ShCqEnaxduxZnZ2c++ugj3nnnHfr378++fftYvny5UZxWq2XPnj1s3LiRlJQUPv30Uzp37syjjz5K//79lbjXX3+d0NBQdu/ezYoVK3B3d+ehhx7iiSeeUGJeeuklevbsyWuvvcamTZuoqqqiW7duhIWFMXny5EbL++yzz9KtWzc2bdrEG2+8QWlpKT4+Pjz66KNMmTLFpnV2cXFhx44d/OEPfyAuLo7KykoWL15sNqECLFiwAC8vL7Zv387y5ctxd3dn/PjxrFixQjbSF61G1qEKIYQQNiBjqEIIIYQNSEIVQgghbEASqhBCCGEDklCFEEIIG5CEKoQQQtiAJFQhhBDCBiShCiGEEDYgCVUIIYSwAUmoQgghhA1IQhVCCCFs4P8DRAnkiCHfIC4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hybrid.plot.scatter('acceleration', 'msrp')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the positive association. The scatter of points is sloping upwards, indicating that cars with greater acceleration tended to cost more, on average; conversely, the cars that cost more tended to have greater acceleration on average. \n",
    "\n",
    "The scatter diagram of MSRP versus mileage shows a negative association. Hybrid cars with higher mileage tended to cost less, on average. This seems surprising till you consider that cars that accelerate fast tend to be less fuel efficient and have lower mileage. As the previous scatter plot showed, those were also the cars that tended to cost more. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hybrid.plot.scatter('mpg', 'msrp')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Along with the negative association, the scatter diagram of price versus efficiency shows a non-linear relation between the two variables. The points appear to be clustered around a curve, not around a straight line. \n",
    "\n",
    "If we restrict the data just to the SUV class, however, the association between price and efficiency is still negative but the relation appears to be more linear. The relation between the price and acceleration of SUV's also shows a linear trend, but with a positive slope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suv = hybrid[hybrid['class'] == 'SUV']\n",
    "suv.plot.scatter('mpg', 'msrp')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suv.plot.scatter('acceleration', 'msrp')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will have noticed that we can derive useful information from the general orientation and shape of a scatter diagram even without paying attention to the units in which the variables were measured.\n",
    "\n",
    "Indeed, we could plot all the variables in standard units and the plots would look the same. This gives us a way to compare the degree of linearity in two scatter diagrams.\n",
    "\n",
    "Recall that in an earlier section we defined the function `standard_units` to convert an array of numbers to standard units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standard_units(any_numbers):\n",
    "    \"Convert any array of numbers to standard units.\"\n",
    "    return (any_numbers - np.mean(any_numbers))/np.std(any_numbers)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this function to re-draw the two scatter diagrams for SUVs, with all the variables measured in standard units.  \n",
    "\n",
    "(*employing Pandas plotting function instead of matplotlib*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SUVs = pd.DataFrame(\n",
    "    {'mpg (standard units)':standard_units(suv['mpg']), \n",
    "    'msrp (standard units)':standard_units(suv['msrp'])}\n",
    ").plot.scatter(0, 1)\n",
    "\n",
    "plt.xlim(-3, 3)\n",
    "\n",
    "plt.ylim(-3, 3)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(\n",
    "    {'acceleration (standard units)':standard_units(suv['acceleration']), \n",
    "    'msrp (standard units)':standard_units(suv['msrp'])}).plot.scatter(0, 1)\n",
    "\n",
    "plt.xlim(-3, 3)\n",
    "\n",
    "plt.ylim(-3, 3)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The associations that we see in these figures are the same as those we saw before. Also, because the two scatter diagrams are now drawn on exactly the same scale, we can see that the linear relation in the second diagram is a little more fuzzy than in the first.\n",
    "\n",
    "We will now define a measure that uses standard units to quantify the kinds of association that we have seen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The correlation coefficient ###\n",
    "\n",
    "The *correlation coefficient* measures the strength of the linear relationship between two variables. Graphically, it measures how clustered the scatter diagram is around a straight line.\n",
    "\n",
    "The term *correlation coefficient* isn't easy to say, so it is usually shortened to *correlation* and denoted by $r$.\n",
    "\n",
    "Here are some mathematical facts about $r$ that we will just observe by simulation.\n",
    "\n",
    "- The correlation coefficient $r$ is a number between $-1$ and 1.\n",
    "- $r$ measures the extent to which the scatter plot clusters around a straight line.\n",
    "- $r = 1$ if the scatter diagram is a perfect straight line sloping upwards, and $r = -1$ if the scatter diagram is a perfect straight line sloping downwards."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function ``r_scatter`` takes a value of $r$ as its argument and simulates a scatter plot with a correlation very close to $r$. Because of randomness in the simulation, the correlation is not expected to be exactly equal to $r$.\n",
    "\n",
    "Call ``r_scatter`` a few times, with different values of $r$ as the argument, and see how the scatter plot changes. \n",
    "\n",
    "When $r=1$ the scatter plot is perfectly linear and slopes upward. When $r=-1$, the scatter plot is perfectly linear and slopes downward. When $r=0$, the scatter plot is a formless cloud around the horizontal axis, and the variables are said to be *uncorrelated*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_scatter(0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_scatter(0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_scatter(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_scatter(-0.55)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating $r$ ###\n",
    "\n",
    "The formula for $r$ is not apparent from our observations so far. It has a mathematical basis that is outside the scope of this class. However, as you will see, the calculation is straightforward and helps us understand several of the properties of $r$.\n",
    "\n",
    "**Formula for $r$**:\n",
    "\n",
    "**$r$ is the average of the products of the two variables, when both variables are measured in standard units.**\n",
    "\n",
    "Here are the steps in the calculation. We will apply the steps to a simple table of values of $x$ and $y$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(1, 7, 1)\n",
    "y = np.array([2, 3, 1, 5, 2, 7])\n",
    "t = pd.DataFrame({'x':x,\n",
    "        'y':y})\n",
    "t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the scatter diagram, we expect that $r$ will be positive but not equal to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t.plot.scatter(0, 1, s=30, color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1.** Convert each variable to standard units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_su = t.copy()\n",
    "\n",
    "t_su['x (standard units)'] = standard_units(x)\n",
    "\n",
    "t_su['y (standard units)'] = standard_units(y)\n",
    "\n",
    "t_su"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 2.** Multiply each pair of standard units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_product = t_su.copy()\n",
    "\n",
    "t_product['product of standard units'] = t_su.iloc[:,2] * t_su.iloc[:,3]\n",
    "\n",
    "t_product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3.** $r$ is the average of the products computed in Step 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# r is the average of the products of standard units\n",
    "\n",
    "r = np.mean(t_product.iloc[:,4])\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, $r$ is positive but not equal to 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Properties of $r$\n",
    "\n",
    "The calculation shows that:\n",
    "\n",
    "- $r$ is a pure number. It has no units. This is because $r$ is based on standard units.\n",
    "- $r$ is unaffected by changing the units on either axis. This too is because $r$ is based on standard units.\n",
    "- $r$ is unaffected by switching the axes. Algebraically, this is because the product of standard units does not depend on which variable is called $x$ and which $y$. Geometrically, switching axes reflects the scatter plot about the line $y=x$, but does not change the amount of clustering nor the sign of the association."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t.plot.scatter('y', 'x', s=30, color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `correlation` function ###\n",
    "We are going to be calculating correlations repeatedly, so it will help to define a function that computes it by performing all the steps described above. Let's define a function ``correlation`` that takes a table and the labels of two columns in the table. The function returns $r$, the mean of the products of those column values in standard units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def correlation(t, x, y):\n",
    "    return np.mean(standard_units(t[x])*standard_units(t[y]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's call the function on the ``x`` and ``y`` columns of ``t``. The function returns the same answer to the correlation between $x$ and $y$ as we got by direct application of the formula for $r$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation(t, 'x', 'y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we noticed, the order in which the variables are specified doesn't matter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation(t, 'y', 'x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calling ``correlation`` on columns of the table ``suv`` gives us the correlation between price and mileage as well as the correlation between price and acceleration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation(suv, 'mpg', 'msrp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation(suv, 'acceleration', 'msrp')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These values confirm what we had observed: \n",
    "\n",
    "- There is a negative association between price and efficiency, whereas the association between price and acceleration is positive.\n",
    "- The linear relation between price and acceleration is a little weaker (correlation about 0.5) than between price and mileage (correlation about -0.67). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Correlation is a simple and powerful concept, but it is sometimes misused. Before using $r$, it is important to be aware of what correlation does and does not measure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Association is not Causation ###\n",
    "\n",
    "Correlation only measures association. Correlation does not imply causation. Though the correlation between the weight and the math ability of children in a school district may be positive, that does not mean that doing math makes children heavier or that putting on weight improves the children's math skills. Age is a confounding variable: older children are both heavier and better at math than younger children, on average."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlation Measures *Linear* Association ###\n",
    "Correlation measures only one kind of association – linear. Variables that have strong non-linear association might have very low correlation. Here is an example of variables that have a perfect quadratic relation $y = x^2$ but have correlation equal to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_x = np.arange(-4, 4.1, 0.5)\n",
    "nonlinear = pd.DataFrame(\n",
    "        {'x':new_x,\n",
    "        'y':new_x**2}\n",
    "    )\n",
    "nonlinear.plot.scatter('x', 'y', s=30, color='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation(nonlinear, 'x', 'y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlation is Affected by Outliers ###\n",
    "Outliers can have a big effect on correlation. Here is an example where a scatter plot for which $r$ is equal to 1 is turned into a plot for which $r$ is equal to 0, by the addition of just one outlying point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line = pd.DataFrame(\n",
    "        {'x':np.array([1, 2, 3, 4]),\n",
    "        'y':np.array([1, 2, 3, 4])}\n",
    "    )\n",
    "line.plot.scatter('x', 'y', s=30, color='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation(line, 'x', 'y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outlier = pd.DataFrame(\n",
    "        {'x':np.array([1, 2, 3, 4, 5]),\n",
    "        'y':np.array([1, 2, 3, 4, 0])}\n",
    "    )\n",
    "outlier.plot.scatter('x', 'y', s=30, color='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation(outlier, 'x', 'y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ecological Correlations Should be Interpreted with Care ###\n",
    "Correlations based on aggregated data can be misleading. As an example, here are data on the Critical Reading and Math SAT scores in 2014. There is one point for each of the 50 states and one for Washington, D.C. The column ``Participation Rate`` contains the percent of high school seniors who took the test. The next three columns show the average score in the state on each portion of the test, and the final column is the average of the total scores on the test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sat2014 = pd.read_csv(path_data + 'sat2014.csv').sort_values(by=['State'])\n",
    "sat2014"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scatter diagram of Math scores versus Critical Reading scores is very tightly clustered around a straight line; the correlation is close to 0.985. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sat2014.plot.scatter('Critical Reading', 'Math')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation(sat2014, 'Critical Reading', 'Math')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's an extremely high correlation. But it's important to note that this does not reflect the strength of the relation between the Math and Critical Reading scores of *students*. \n",
    "\n",
    "The data consist of average scores in each state. But states don't take tests – students do. The data in the table have been created by lumping all the students in each state into a single point at the average values of the two variables in that state. But not all students in the state will be at that point, as students vary in their performance. If you plot a point for each student instead of just one for each state, there will be a cloud of points around each point in the figure above. The overall picture will be more fuzzy. The correlation between the Math and Critical Reading scores of the students will be *lower* than the value calculated based on state averages.\n",
    "\n",
    "Correlations based on aggregates and averages are called *ecological correlations* and are frequently reported. As we have just seen, they must be interpreted with care."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Serious or tongue-in-cheek?\n",
    "\n",
    "In 2012, a [paper](http://www.biostat.jhsph.edu/courses/bio621/misc/Chocolate%20consumption%20cognitive%20function%20and%20nobel%20laurates%20%28NEJM%29.pdf) in the respected New England Journal of Medicine examined the relation between chocolate consumption and Nobel Prizes in a group of countries. The [Scientific American](http://blogs.scientificamerican.com/the-curious-wavefunction/chocolate-consumption-and-nobel-prizes-a-bizarre-juxtaposition-if-there-ever-was-one/) responded seriously whereas\n",
    "[others](http://www.reuters.com/article/2012/10/10/us-eat-chocolate-win-the-nobel-prize-idUSBRE8991MS20121010#vFdfFkbPVlilSjsB.97) were more relaxed. You are welcome to make your own decision! The following graph, provided in the paper, should motivate you to go and take a look."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove_input"
    ]
   },
   "outputs": [],
   "source": [
    "\n",
    "from IPython.display import Image\n",
    "Image(\"chocoNobel.png\")"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}