chapter_15/multiple-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiple Regression\n",
    "\n",
    "Chapter 15 of _Data Science from Scratch_ by Joel Grus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os.path\n",
    "import random\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from functools import partial\n",
    "from math import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "book_dir = '/Users/CBare/Documents/projects/data-science-from-scratch'\n",
    "sys.path.extend(os.path.join(book_dir, 'chapter_{:02d}'.format(i)) for i in [3,4,5,6,7,8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stats import mean, median, de_mean, standard_deviation, correlation\n",
    "from gradient_descent import minimize_stochastic\n",
    "from vector import dot, vector_add\n",
    "from normal import normal_cdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(x_i, beta):\n",
    "    \"\"\"assumes that the first element of each x_i is 1\"\"\"\n",
    "    return dot(x_i, beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predict Hours on Site"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal is to predict time spent on a site given 3 features, number of friends, work hours per day, and whether or not the person has a PhD.\n",
    "\n",
    "$$ minutes = \\alpha  + \\beta_1 friends + \\beta_2 work hours + \\beta_3 phd + \\epsilon $$\n",
    "\n",
    "Each person is represented as a vector of 4 elements, for example:\n",
    "```\n",
    "[ 1, # constant term\n",
    " 49, # number of friends\n",
    "  4, # work hours per day\n",
    "  0] # doesn't have PhD\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: The data below was snagged from Joel's source code here: https://github.com/joelgrus/data-science-from-scratch/blob/master/code-python3/multiple_regression.py#L95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [[1,49,4,0],[1,41,9,0],[1,40,8,0],[1,25,6,0],[1,21,1,0],[1,21,0,0],[1,19,3,0],\n",
    "     [1,19,0,0],[1,18,9,0],[1,18,8,0],[1,16,4,0],[1,15,3,0],[1,15,0,0],[1,15,2,0],\n",
    "     [1,15,7,0],[1,14,0,0],[1,14,1,0],[1,13,1,0],[1,13,7,0],[1,13,4,0],[1,13,2,0],\n",
    "     [1,12,5,0],[1,12,0,0],[1,11,9,0],[1,10,9,0],[1,10,1,0],[1,10,1,0],[1,10,7,0],\n",
    "     [1,10,9,0],[1,10,1,0],[1,10,6,0],[1,10,6,0],[1,10,8,0],[1,10,10,0],[1,10,6,0],\n",
    "     [1,10,0,0],[1,10,5,0],[1,10,3,0],[1,10,4,0],[1,9,9,0],[1,9,9,0],[1,9,0,0],\n",
    "     [1,9,0,0],[1,9,6,0],[1,9,10,0],[1,9,8,0],[1,9,5,0],[1,9,2,0],[1,9,9,0],[1,9,10,0],\n",
    "     [1,9,7,0],[1,9,2,0],[1,9,0,0],[1,9,4,0],[1,9,6,0],[1,9,4,0],[1,9,7,0],[1,8,3,0],\n",
    "     [1,8,2,0],[1,8,4,0],[1,8,9,0],[1,8,2,0],[1,8,3,0],[1,8,5,0],[1,8,8,0],[1,8,0,0],\n",
    "     [1,8,9,0],[1,8,10,0],[1,8,5,0],[1,8,5,0],[1,7,5,0],[1,7,5,0],[1,7,0,0],[1,7,2,0],\n",
    "     [1,7,8,0],[1,7,10,0],[1,7,5,0],[1,7,3,0],[1,7,3,0],[1,7,6,0],[1,7,7,0],[1,7,7,0],\n",
    "     [1,7,9,0],[1,7,3,0],[1,7,8,0],[1,6,4,0],[1,6,6,0],[1,6,4,0],[1,6,9,0],[1,6,0,0],\n",
    "     [1,6,1,0],[1,6,4,0],[1,6,1,0],[1,6,0,0],[1,6,7,0],[1,6,0,0],[1,6,8,0],[1,6,4,0],\n",
    "     [1,6,2,1],[1,6,1,1],[1,6,3,1],[1,6,6,1],[1,6,4,1],[1,6,4,1],[1,6,1,1],[1,6,3,1],\n",
    "     [1,6,4,1],[1,5,1,1],[1,5,9,1],[1,5,4,1],[1,5,6,1],[1,5,4,1],[1,5,4,1],[1,5,10,1],\n",
    "     [1,5,5,1],[1,5,2,1],[1,5,4,1],[1,5,4,1],[1,5,9,1],[1,5,3,1],[1,5,10,1],[1,5,2,1],\n",
    "     [1,5,2,1],[1,5,9,1],[1,4,8,1],[1,4,6,1],[1,4,0,1],[1,4,10,1],[1,4,5,1],[1,4,10,1],\n",
    "     [1,4,9,1],[1,4,1,1],[1,4,4,1],[1,4,4,1],[1,4,0,1],[1,4,3,1],[1,4,1,1],[1,4,3,1],\n",
    "     [1,4,2,1],[1,4,4,1],[1,4,4,1],[1,4,8,1],[1,4,2,1],[1,4,4,1],[1,3,2,1],[1,3,6,1],\n",
    "     [1,3,4,1],[1,3,7,1],[1,3,4,1],[1,3,1,1],[1,3,10,1],[1,3,3,1],[1,3,4,1],[1,3,7,1],\n",
    "     [1,3,5,1],[1,3,6,1],[1,3,1,1],[1,3,6,1],[1,3,10,1],[1,3,2,1],[1,3,4,1],[1,3,2,1],\n",
    "     [1,3,1,1],[1,3,5,1],[1,2,4,1],[1,2,2,1],[1,2,8,1],[1,2,3,1],[1,2,1,1],[1,2,9,1],\n",
    "     [1,2,10,1],[1,2,9,1],[1,2,4,1],[1,2,5,1],[1,2,0,1],[1,2,9,1],[1,2,9,1],[1,2,0,1],\n",
    "     [1,2,1,1],[1,2,1,1],[1,2,4,1],[1,1,0,1],[1,1,2,1],[1,1,2,1],[1,1,5,1],[1,1,3,1],\n",
    "     [1,1,10,1],[1,1,6,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],[1,1,4,1],[1,1,9,1],[1,1,9,1],\n",
    "     [1,1,4,1],[1,1,2,1],[1,1,9,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],[1,1,1,1],[1,1,1,1],\n",
    "     [1,1,5,1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = [68.77,51.25,52.08,38.36,44.54,57.13,51.4,41.42,31.22,34.76,\n",
    "     54.01,38.79,47.59,49.1,27.66,41.03,36.73,48.65,28.12,46.62,\n",
    "     35.57,32.98,35,26.07,23.77,39.73,40.57,31.65,31.21,36.32,20.45,\n",
    "     21.93,26.02,27.34,23.49,46.94,30.5,33.8,24.23,21.4,27.94,32.24,\n",
    "     40.57,25.07,19.42,22.39,18.42,46.96,23.72,26.41,26.97,36.76,\n",
    "     40.32,35.02,29.47,30.2,31,38.11,38.18,36.31,21.03,30.86,36.07,\n",
    "     28.66,29.08,37.28,15.28,24.17,22.31,30.17,25.53,19.85,35.37,\n",
    "     44.6,17.23,13.47,26.33,35.02,32.09,24.81,19.33,28.77,24.26,31.98,\n",
    "     25.73,24.86,16.28,34.51,15.23,39.72,40.8,26.06,35.76,34.76,16.13,\n",
    "     44.04,18.03,19.65,32.62,35.59,39.43,14.18,35.24,40.13,41.82,35.45,\n",
    "     36.07,43.67,24.61,20.9,21.9,18.79,27.61,27.21,26.61,29.77,20.59,\n",
    "     27.53,13.82,33.2,25,33.1,36.65,18.63,14.87,22.2,36.81,25.53,24.62,\n",
    "     26.25,18.21,28.08,19.42,29.79,32.8,35.99,28.32,27.79,35.88,29.06,\n",
    "     36.28,14.1,36.63,37.49,26.9,18.58,38.48,24.48,18.95,33.55,14.24,\n",
    "     29.04,32.51,25.63,22.22,19,32.73,15.16,13.9,27.2,32.01,29.27,33,\n",
    "     13.74,20.42,27.32,18.23,35.35,28.48,9.08,24.62,20.12,35.26,19.92,\n",
    "     31.02,16.49,12.16,30.7,31.22,34.65,13.13,27.51,33.2,31.57,14.1,\n",
    "     33.42,17.44,10.12,24.42,9.82,23.39,30.93,15.03,21.67,31.09,33.29,\n",
    "     22.61,26.89,23.48,8.38,27.81,32.35,23.84]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error(x_i, y_i, beta):\n",
    "    return y_i - predict(x_i, beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def squared_error(x_i, y_i, beta):\n",
    "    return error(x_i, y_i, beta) ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def squared_error_gradient(x_i, y_i, beta):\n",
    "    \"\"\"the gradient (with respect to beta) corresponding to the ith squared error term\"\"\"\n",
    "    return [-2 * x_ij * error(x_i, y_i, beta)\n",
    "            for x_ij in x_i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_beta(x, y):\n",
    "    beta_initial = [random.random() for x_i in x[0]]\n",
    "    return minimize_stochastic(squared_error,\n",
    "                               squared_error_gradient,\n",
    "                               x, y,\n",
    "                               beta_initial,\n",
    "                               0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[30.619881701311712, 0.9702056472470465, -1.8671913880379478, 0.9163711597955347]\n"
     ]
    }
   ],
   "source": [
    "random.seed(0)\n",
    "beta = estimate_beta(x, y)\n",
    "print(beta)\n",
    "# [30.63, 0.972, -1.868, 0.911]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using the model to make predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 6, 3, 1]\n",
      "39.43\n",
      "predicted= 31.755912580475684\n"
     ]
    }
   ],
   "source": [
    "i = 100\n",
    "print(x[i])\n",
    "print(y[i])\n",
    "print('predicted=', dot(beta, x[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Are friends and work hour correlated? Not really."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02529930456534995"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "correlation([x_i[1] for x_i in x], [x_i[2] for x_i in x])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PhD's have no friends. Sad."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.6107431572233505"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "correlation([x_i[1] for x_i in x], [x_i[3] for x_i in x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAF/FJREFUeJzt3X9wXeV95/H3p8JIcTBg1SolIcawgHvpNcGNJklBzaAkDbRNcbrJkjo0Q4rGrjOtaGromlqTBXcrN55mTDdOd128cnB2NjdAtgkuEwopURJEGho5EGPQQKjBdTAEuxib4NgR5rt/3GMjG/24tnXukfR8XjN3fO9zfn01c62PzvOc8xxFBGZmlq5fKLoAMzMrloPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgKbUCTNlfSIpJclXTfCOmslfTqn44ek8/LY93iRNCer86Sia7GpwUFgE81/BXojYkZEfG64FSJiSUT89zrXVVeSLpP0mqSfZqH4hKQ/rHHbmyUNZtu9LOlJSZ+XdGbeddvk5CCwieZs4LGRFkpqqGMtRdsREacApwLLgHWSLqxx29sjYgbQDPwe8MvAJoeBDcdBYBOGpG8C7cDns7+EL5B0m6T/Jenrkl4B2rO2vxqy3Qez7qSXJH1X0kVDlj0j6QZJmyXtkXS7pKYhy/9c0nOSdki69qh6flvS49lf1c9KumGYmhuz45aHtLVI+pmkX5I0S9Ld2TovSnpA0jH9v4uqrwG7gaFBcLWkf5e0S1LXCNsORsRjwEeBncD1x3JsS4ODwCaMiHgv8ADwJxFxSkQ8mS36GNANzAD6hm4jaT6wHvgj4BeBvwc2SmocstpVwBXAOcBFwCeyba8AbgB+EzgfeP9RJfUAf5T9ZV0GvjlMzQeAfwAWHnW8b0fEC1R/8f4YaAHOAJYDxzSvi6RfkPR7wOnAo0MWtQFzgfcB/01SaaR9RMRB4C7gN47l2JYGB4FNBndFxIMR8VpE7D9q2WLg7yPioYg4GBEbgAPAu4es87mI2BERLwL/CFyctV8FfCEitkTEK8DNR+17ELhQ0qkRsTsifjBCfV8Cfn/I549lbYf2cSZwdvbX+QNR+wRfb5H0ErALuAn4eEQ8MWT5ioj4WUT8EPgh8PYx9reDaleR2REcBDYZbB9l2dnA9VnXy0vZL863AW8Zss7zQ97vA07J3r/lqH1vO2rfHwZ+G9gm6duSfn2EGnqB6ZLeJWkO1aD5arbsb4CngPskbZV04yg/y9F2RMTpEdEcERdHxJePWj7SzzWStwIvHsPxLREOApsMRvsLejvQnf3CPPSaHhGVGvb7HNXQOGT2EQeN+H5ELAB+CfgacMewxVW7Xe6g2j20ELg7Il7Olr0cEddHxLnAlcBSSe+robZxlY1L/C7VrjezIzgIbLJbByzJ/hqXpDdL+h1JM2rY9g7gE5IulDSdavcLAJJOlnS1pNMiYhDYC7w2yr6+RHVA9mpe7xY6NJB9niQBe4CDY+xnXEk6KRs7qFC9cmh1vY5tk4eDwCa1iOgHFgGfp3pVzVNkg8E1bHsP8LdUB4Gf4o2DwR8HnpG0F1hC9Zf8SPt6CHiFanfTPUMWnQ/8M/BT4F+A/xkRvQCS7pG0vJZaj8NHJf2UavhsBP4DeEdE7MjpeDaJyQ+mMTNLm88IzMwS5yAwM0ucg8DMLHEOAjOzxE2KaWxnzZoVc+bMKboMM7NJZdOmTbsiomWs9SZFEMyZM4f+/v6iyzAzm1QkHX23/LDcNWRmlrhczwgkPQO8TPVuylcjolVSM3A7MAd4BrgqInbnWYeZmY2sHmcE7dmEWa3Z5xuB+yPifOD+7LOZmRWkiK6hBcCG7P0G4EMF1GBmZpm8gyCoTr+7SdLirO2MiHgue/881Yd1vIGkxZL6JfXv3Lkz5zLNzNKVdxC0RcSvAb8F/LGk9wxdmD2gY9jJjiLi1ohojYjWlpYxr36yGlQqFcrlMg0NDZTLZSqVWmZqNrOpLtfB4oh4Nvv3BUlfBd4J/ETSmRHxXPYg7RfyrMGqKpUKXV1d9PT00NbWRl9fHx0dHQAsXLhwjK3NbCrL7Ywgmxd+xqH3wAeALVSnxL0mW+0aqs9RtZx1d3fT09NDe3s706ZNo729nZ6eHrq7u4suzcwKlts01JLO5fXH9Z0EfCkiuiX9ItUHgsym+mjAq7JnyY6otbU1fEPZiWloaGD//v1MmzbtcNvg4CBNTU0cPHiwwMrMLC+SNg25YnNEuXUNRcRWhnmYdkT8B1D3R/WlrlQq0dfXR3t7++G2vr4+SqVSgVWZ2UTgO4sT0dXVRUdHB729vQwODtLb20tHRwddXV1Fl2ZmBZsUcw3ZiTs0INzZ2cnAwAClUonu7m4PFJvZ5HhUpccIzMyOXa1jBO4aMjNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIMgIZVKhXK5TENDA+VymUqlUnRJZjYBnFR0AVYflUqFrq4uenp6aGtro6+vj46ODgAWLlxYcHVmViRFRNE1jKm1tTX6+/uLLmNSK5fLrFmzhvb29sNtvb29dHZ2smXLlgIrM7O8SNoUEa1jrucgSENDQwP79+9n2rRph9sGBwdpamri4MGDBVZmZnmpNQhyHyOQ1CDpYUl3Z5/PkfSQpKck3S7p5LxrMCiVSqxYseKIMYIVK1ZQKpWKLs3MClaPweI/BQaGfF4F3BIR5wG7gY461JC89vZ2Vq1axbXXXsvLL7/Mtddey6pVq47oKjKzNOUaBJLOAn4H+N/ZZwHvBb6SrbIB+FCeNVhVb28vy5YtY/369cyYMYP169ezbNkyent7iy7NzAqW6xiBpK8Afw3MAG4APgF8LzsbQNLbgHsiojzMtouBxQCzZ89+x7Zt23KrMwUeIzBLT+FjBJI+CLwQEZuOZ/uIuDUiWiOitaWlZZyrS0+pVKKvr++Itr6+Po8RmFmuXUOXAldKegb4MtUuof8BnC7p0P0LZwHP5liDZbq6uujo6KC3t5fBwUF6e3vp6Oigq6ur6NLMrGC53VAWEX8B/AWApMuAGyLiakl3Ah+hGg7XAHflVYO97tBNY52dnQwMDFAqleju7vbNZGZWn/sIhgTBByWdSzUEmoGHgT+IiAOjbe/7CMzMjl2tYwR1mWIiIr4FfCt7vxV4Zz2Oa2ZmY/Okc2ZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeJyCwJJTZL+VdIPJT0maUXWfo6khyQ9Jel2SSfnVYOZmY0tzzOCA8B7I+LtwMXAFZLeDawCbomI84DdQEeONZiZ2RhyC4Ko+mn2cVr2CuC9wFey9g3Ah/KqwczMxpbrGIGkBkmPAC8A3wD+DXgpIl7NVvkx8NYRtl0sqV9S/86dO/Ms08wsabkGQUQcjIiLgbOAdwK/cgzb3hoRrRHR2tLSkluNZmapq8tVQxHxEtAL/DpwuqSTskVnAc/WowYzMxtenlcNtUg6PXv/JuA3gQGqgfCRbLVrgLvyqsHMzMZ20tirHLczgQ2SGqgGzh0Rcbekx4EvS/or4GGgJ8cazMxsDLkFQURsBuYP076V6niBmZlNAMfUNSRppqSL8irGzMzqb8wgkPQtSadKagZ+AKyTtDr/0szMrB5qOSM4LSL2Av8Z+GJEvAt4f75lmZlZvdQSBCdJOhO4Crg753rMzKzOagmCvwTuBZ6KiO9LOhf4Ub5lmZlZvYx51VBE3AncOeTzVuDDeRZlZmb1M2IQSFpDdZK4YUXEdblUZGZmdTVa11A/sAloAn6NanfQj6hOKe1nCJiZTREjnhFExAYASZ8E2g7NGCppLfBAfcozM7O81TJYPBM4dcjnU7I2MzObAmqZYuIzwMOSegEB7wFuzrMoMzOrn1quGvqCpHuAd2VNyyLi+XzLMjOzeql1rqEGYCfVZwxfIOk9+ZVkealUKpTLZRoaGiiXy1QqlaJLMrMJYMwzAkmrgI8CjwGvZc0BfCfHumycVSoVurq66Onpoa2tjb6+Pjo6OgBYuHBhwdWZWZEUMeKtAtUVpCeAiyLiQH1KeqPW1tbo7+8v6vBTQrlcZs2aNbS3tx9u6+3tpbOzky1bthRYmZnlRdKmiGgda71auoa2AtNOvCQr0sDAAG1tbUe0tbW1MTAwUFBFZjZR1HLV0D7gEUn3A4fPCnxn8eRSKpXo6+s74oygr6+PUqlUYFVmNhHUEgQbs5dNYl1dXXR0dLxhjKC7u7vo0sysYLVcProhe/j87Ih4og41WQ4ODQh3dnYyMDBAqVSiu7vbA8VmVtNg8e8CnwVOjohzJF0M/GVEXFmPAsGDxWZmx2M8B4tvpvqw+ZcAIuIR4NwTqs7MzCaMWoJgMCL2HNX22rBrmpnZpFPLYPFjkj4GNEg6H7gO+G6+ZZmZWb3UckbQCfwq1UtHK8Be4FN5FmVmZvVTy1VD+4Cu7GVmZlPMaI+q/NuI+JSkf2SYR1bW86ohMzPLz2hnBP8n+/ez9SjEzMyKMdqjKjdJagAWR8TVdazJzMzqaNTB4og4CJwtyQ+rNzObomqdffRBSZ+WtPTQK+/CbPx1dnbS1NSEJJqamujs7Cy6JDObAEYMAkmHxgiuBO7O1p0x5GWTSGdnJ2vXrmXlypW88sorrFy5krVr1zoMzGzkuYYkPQ68H/gn4LKjl0fEi7lWNoTnGjpxTU1NrFy5kqVLXz+ZW716NcuXL2f//v0FVmZmeal1rqHRguA64JPAOcCOoYuAiIi6zTfkIDhxknjllVeYPn364bZ9+/bx5je/mbEmHjSzyemEJ52LiM9FRAn4QkScO+R1Ti0hIOltknolPS7pMUl/mrU3S/qGpB9l/848pp/MjktjYyNr1649om3t2rU0NjYWVJGZTRRjDhZHxCePc9+vAtdHxIXAu4E/lnQhcCNwf0ScD9yffbacLVq0iGXLlrF69Wr27dvH6tWrWbZsGYsWLSq6NDMr2JjPIxi3A0l3AZ/PXpdFxHOSzgS+FRFzR9vWXUPjo7Ozk3Xr1nHgwAEaGxtZtGgRa9asKbosM8vJCY8RjHMxc4DvAGXg3yPi9KxdwO5Dn4/aZjGwGGD27Nnv2LZtW+51mplNJeP5YJoTLeQU4P8Bn4qIvUOXRTWFhk2iiLg1IlojorWlpSXvMpPg+wjMbDi5BoGkaVRD4P9GxD9kzT/JuoTI/n0hzxqsyvcRmNlIcguCrNunBxiIiNVDFm0ErsneXwPclVcN9rp169axatUqli5dyvTp01m6dCmrVq1i3bp1RZdmZgXLbYxAUhvwAPAorz/acjnwEHAHMBvYBlw11s1pHiw+cb6PwCw9tY4R1PKoyuMSEX1Ubz4bzvvyOq4Nr7GxkQ984AP09/cfvmqotbXV9xGYWf6DxTYxXHDBBTz44INcfvnl7Ny5k8svv5wHH3yQCy64oOjSzKxguZ0R2MTy5JNPcumll3LvvffS0tJCY2Mjl156Ke5yMzMHQSIOHDjAfffdN+wYgZmlzV1DifBcQ2Y2Ep8RJOLQXEMAS5YsYe3atSxbtowlS5YUXJmZFc1BkIhDcwotX76c66+/nsbGRpYsWeK5hsysfpPOnQjfR2BmduwmzFxDNnFUKhXK5TINDQ2Uy2UqlUrRJZnZBOCuoURUKhW6urro6emhra2Nvr4+Ojo6AFi4cGHB1ZlZkdw1lIhyucyaNWtob28/3Nbb20tnZydbtmwpsDIzy8uEeh7BiXIQnLiGhgb279/PtGnTDrcNDg7S1NTEwYMHC6zMzPLiMQI7QqlUYu7cuUg6/Jo7dy6lUqno0sysYA6CROzdu5enn36aSy65hB07dnDJJZfw9NNPs3fv3rE3NrMpzUGQiO3btzN//nz27NnDWWedxZ49e5g/fz7bt28vujQzK5ivGkrIfffdx6xZsw5/3rVrF34MqJn5jCAhhy4XHemzmaXJQZCIefPmsXHjRhYsWMCuXbtYsGABGzduZN68eUWXZmYFc9dQIjZv3sxFF13Exo0bD3cHzZs3j82bNxdcmZkVzUEwRUkjPSX0dY8++ugb1psM95WY2fhyEExRo/1Cl+Rf+GZ2mMcIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxOUWBJLWS3pB0pYhbc2SviHpR9m/M/M6vpmZ1SbPM4LbgCuOarsRuD8izgfuzz6bmVmBcguCiPgO8OJRzQuADdn7DcCH8jq+mZnVpt5jBGdExHPZ++eBM0ZaUdJiSf2S+nfu3Fmf6iah5uZmJB3TCzim9Zubmwv+Kc0sT4U9mCYiQtKIT0eJiFuBWwFaW1v9FJUR7N69O/eHzNTytDMzm7zqfUbwE0lnAmT/vlDn45uZ2VHqHQQbgWuy99cAd9X5+GZmdpTcuoYkVYDLgFmSfgzcBHwGuENSB7ANuCqv46cibjoVbj4t/2OY2ZSVWxBExMIRFr0vr2OmSCv21mWMIG7O9RBmViDfWWxmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4gqbYsLGT95TQMyc6dnCzaYyB8Ekdzz3EEjK/d4DM5s83DVkZpY4B4GZWeLcNTRFjTVuMNJydxlZPRzPuJa/m/lxEExR/k9jE9lI30+PXxXDXUNmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmeWiubkZScf0Ao5p/ebm5oJ/yqnBN5SZWS52796d+81hec+8mwqfEZiZJc5BYGaWOHcNmVku4qZT4ebT8j+GnTAHgZnlQiv21mWMIG7O9RBJcNeQmVniHARmZolz15CZ5SbvyztnzpyZ6/5T4SAws1wcz/iAH0xTDHcNmZklzkFgZpY4dw2ZWd2NNnYw0jJ3GeWnkDMCSVdIekLSU5JuLKIGMytORBzzy/JT9yCQ1AD8HfBbwIXAQkkX1rsOMzOrKuKM4J3AUxGxNSJ+DnwZWFBAHWZmRjFB8FZg+5DPP87ajiBpsaR+Sf07d+6sW3FmZqmZsFcNRcStEdEaEa0tLS1Fl2NmNmUVEQTPAm8b8vmsrM3MzApQRBB8Hzhf0jmSTgZ+H9hYQB1mZkYB9xFExKuS/gS4F2gA1kfEY/Wuw8zMqgq5oSwivg58vYhjm5nZkTQZbtSQtBPYVnQdU8gsYFfRRZgNw9/N8XV2RIx5tc2kCAIbX5L6I6K16DrMjubvZjEm7OWjZmZWHw4CM7PEOQjSdGvRBZiNwN/NAniMwMwscT4jMDNLnIPAzCxxDoKESLpZ0rOSHpG0RdKVWfttkj4yzPpzJP1M0sOSBiT9q6RP1L1wmxIkHRzy3btT0vTsO7ZlhPVvk/S0pB9KelLSFyWdVe+6U+AgSM8tEXEx8F+A9ZLG+g78W0TMj4gS1XmhPiXpD3Ov0qain0XExRFRBn4OLKlhmz+PiLcDc4GHgW9mc5TZOHIQTCLZX08DktZJekzSfZLelC27WNL3JG2W9FVJM0fbV0QMAK9SvZMT4D2Svitp63BnB9k2W4GlwHXj+GNZmh4AzsveNwz3nR4qqm4Bnqf6dEMbRw6Cyed84O8i4leBl4APZ+1fBJZFxEXAo8BNo+1E0ruA14BDT/05E2gDPgh8ZpRNfwD8ynFXb8mTdBLVX+aPZk0jfaeH4+9fDgqZdM5OyNMR8Uj2fhMwR9JpwOkR8e2sfQNw5wjb/5mkPwBeBj4aESEJ4GsR8RrwuKQzRjm+TvxHsES9SdKh7+4DQA/wFob5To+yD3//cuAgmHwODHl/EHjDafQYbomIz46x39H+s80HBo7xmGaQjREMbcj+CDmW7/R84P7xLy1t7hqaAiJiD7Bb0m9kTR8Hvj3KJsdF0hzgs8Ca8d632WhUdR3VLsx/KrqeqcZnBFPHNcBaSdOBrcB4XdnznyQ9DDRR7U76XETcNk77NhvL30j6NDAd+B7QHhE/L7imKcdTTJiZJc5dQ2ZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpa4/w/mz4CkgH65kwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "friends_no_phd = [x_i[1] for x_i in x if not x_i[3]]\n",
    "friends_phd = [x_i[1] for x_i in x if x_i[3]]\n",
    "\n",
    "plt.boxplot([friends_no_phd, friends_phd], labels=['no PhD', 'PhD'])\n",
    "plt.title('friends vs. PhD')\n",
    "plt.ylabel('friends')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PhD's and non-PhD's work about the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.031950793986044564"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "correlation([x_i[2] for x_i in x], [x_i[3] for x_i in x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean work hrs non-PhDs: 4.70\n",
      "mean work hrs     PhDs: 4.50\n"
     ]
    }
   ],
   "source": [
    "wrk_hrs_non_phds = mean([x_i[2] for x_i in x if not x_i[3]])\n",
    "wrk_hrs_phds = mean([x_i[2] for x_i in x if x_i[3]])\n",
    "\n",
    "print(f'mean work hrs non-PhDs: {wrk_hrs_non_phds:.2f}')\n",
    "print(f'mean work hrs     PhDs: {wrk_hrs_phds:.2f}'.format())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### People with more friends spend more time on the site"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XucXHV9+P/Xe247s5dkd5NNyOZOEsIlSNAARtRSRfESjPWuiGD5FvTbqm2liLZfq/22FZEfFNp+VaptAFFAUBPoBYFqVYxAIlE2N0lINslukr1lk73M7Nzevz/O2cnsZs5mZzOzszvzfj4eYWfOnJnzObPLeZ/P7f0RVcUYY0zl8pW6AMYYY0rLAoExxlQ4CwTGGFPhLBAYY0yFs0BgjDEVzgKBMcZUOAsEZUBEviAi3yp1OQpFRPpF5OxSl2MqEJEvich3Snj8vxWRLhE54vH6J0XkqPs7m+WxT1F+nyJyvYj8otCfW4kCpS6AOT0R6c96Wg0MASn3+U2q+veTX6riUdXa4ccisgE4pKp/VboSVSYRWQR8Flisqh05Xg8CdwKvVdXfeH1O9u/TTE0WCKaBURfG/cD/UtWnS1ciMx2JSEBVk3m8ZRHQnSsIuOYCYWB7gY5nSsSahspAdvOBiCwRERWRj4vIQRE5JiKfEJFLROS3ItIrIv806v1/KCI73X2fFJHFYxzrXSKy3f2cn4rIeVmv7ReRm93jHBeRh0Uk7PE5y0Xkf9z9ukTk4azX1H39RuAa4Ba3eeFx9/VmEXlMRDpFZJ+IfNrjGJeJyBER8Wdt+wMR+a37+FIR2SIiJ9zmjTvz+L4fEZH7RaTP/T7WjC5/1vMNIvK37uMrROSQiNwiIh0iclhE3i0i7xCR34lIj4h8YdQhw+532ScivxaRi7I+2/O7cMv5qIh8R0ROANfnOJeZ7nl0ikiriPyViPhE5ErgKaDZ/e43jHrfOcBu92mviPx31rn/sYi8DLw8+vsQkSoRuUNEDrjf+TdEJDLqu/ls1nfz8axjzhKRTe7v63lgWdZrIiJ3ue87ISIviciq0/82DQCqav+m0T9gP3DlqG1fAr7jPl4CKPANnLu1twIx4EfAHGA+0AH8nrv/emAPcB5ODfGvgF96HPscYAB4CxAEbnHfG8oq2/NAM9AI7AQ+4fFZ3wP+EudmJAy8Pus1BZa7jzcAf5v1mg/YCnwRCAFnA68AV3kcZy/wlqzn3wdudR9vBq51H9fiNHGM53fwJfc7fQfgB74C/CpX+UefA3AFkHTLHwT+COgEvgvUARcAUWBp1rESwPvc/W8G9rmPx/wust77bnffSI5zuR/Y6B57CfA74Iassh4a43tY4p5rYNS5P+X+/iM5fp93AZvc1+uAx4GvjPpu/sY9v3cAg0CD+/pDwCNADbAKaAN+4b52lftd1AOC8/c8r9T/v06Xf1YjKF//V1VjqvpjnIv391S1Q1XbgJ8DF7v7fQLnf8Sd6lTj/x5Y7VEr+CDw76r6lKomgDuACPC6rH3uUdV2Ve3B+Z98tUf5EsBioNkt53g7/S4BmlT1b1Q1rqqvAP8CfMhj/+8BHwYQkTqci8v3ssqwXERmq2q/qv5qnGUA5wL0H6qaAh4ALjrdG7IkgL9zv8OHgNnA3arap6rbgR2jPm+rqj7q7n8nTuB8LeP7Ljar6o9UNa2q0exCuDWlDwGfd4+9H/j/gGvzOJdcvqKqPTmOJ8CNwJ+5r/fh/L1llzcB/I2qJlT1P4B+YKVb1vcCX1TVAVVtAe4b9b464FxA3L/nw2d4HhXDAkH5Opr1OJrj+XC/w2Lgbreppxfowbmjmp/jM5uB1uEnqpoGDo7aN3t0yWDWcUa7xT3O827Tyh+e9oxOlrd5uLxumb+A016dy3eB94hIFfAe4NeqOnwON+DUcnaJyAsism6cZYBTzzMsIuPtc+t2Awg4vwvw/v2A8x0Dme/8EM7vYjzfxUG8zca5827N2tZK7t99PryO2YQz2GFrVnn/y90+rFtH9isM/w014dRYsz87+2/xv4F/Av4Z6BCRe0VkxhmeR8WwzmJzEOfu9MFx7NsOXDj8xL3DW4hTRc+Lqh7BaRZBRF4PPC0iP1PVPaN3zVHefaq6YpzH2SEircDbgY/gBIbh114GPiwiPpwg8aiIzFLVgXzPZ5RBnAvesLNwLt4TtXD4gVvWBTi/iySn/y7GSi/cxcma2Q532yIm8Psc5zG7cILcBW7NNB+dOOe7ENjlbls04qCq9wD3iMgcnCakvwD+T57HqUhWIzDfAD4vIhdApvPw/R77PgK8U0TeLM7Qwc/iDGX9Zb4HFZH3i8gC9+kxnItHOseuR3Havoc9D/SJyOdEJCIifhFZJSKXjHG47wKfAd6I00cwXIaPikiTe5fd627OVYZ8bQM+4pbtbcDvneHnvUZE3uPWOP4U5zv/FRP7LjLcWskjwN+JSJ3bHPjnQFHmLbjf878Ad7kXa0RkvohcNc6y/gD4kohUi8j5wHXDr4szGOIy9+9yAKcPpxC/y4pggaDCqeoPga8CD7kjS1pw7p5z7bsb+Cjwjzh3d1cDV6tqfAKHvgR4Tpw5EpuAz7ht3KN9GzjfbUr4kXtBWIfT97DPLce3gJljHOt7OBfj/1bVrqztbwO2u2W4G/jQcLu2O1LmDRM4L3CCztU4weUanI76M7ERp3/mGE77/XvcNvSJfBejfQrnwvkK8AucoPmvZ1jesXwOZ4DBr9y/t6eBleN875/gNBMdwemA/7es12bgBJljOE1G3cDXClPk8ieqtjCNMcZUMqsRGGNMhbNAYIwxFc4CgTHGVDgLBMYYU+GmxTyC2bNn65IlS0pdDGOMmVa2bt3apapNp9tvWgSCJUuWsGXLllIXwxhjphV3MuVpWdOQMcZUuKIFAhFZKSLbsv6dEJE/FZFGEXlKRF52fzYUqwzGGGNOr2iBQFV3q+pqVV0NvAYn/8oPgVuBZ9z8KM+4z40xxpTIZDUNvRnY62Z9XM/J9LH34eRKN8YYUyKTFQg+xMkc8HOz8oQfwSN9sIjcKM7qUVs6Ozsno4zGGFORij5qSERCwLuAz49+TVVVRHImO1LVe4F7AdasWWMJkYwxFaOlZR+PP7GZtvYu5jfP5up1a1m1amnRjjcZNYK34ywGMrzwxlERmQfg/vRaGNsYYypOS8s+vv7NTfT29jN3TgO9vf18/ZubaGnZV7RjTkYg+DAnm4XASTk8nEf8OpwUu8YYY4DHn9hMbU2EurpqRIS6umpqayI8/sTmoh2zqIFARGpwFjr/Qdbm24C3iMjLwJXuc2OMMUBbexe1tZER22prI7S1d3m848wVtY/AXfJv1qht3TijiIwxxowyv3k2vb391NWdXO20vz/K/ObZRTumzSw2xpgp5Op1a+kfiNLXN4iq0tc3SP9AlKvXrS3aMS0QGGPMFLJq1VI+edO7qK+v5WjHMerra/nkTe8q6qihaZF0zhhjKsmqVUuLeuEfzWoExhhT4SwQGGNMhbNAYIwxFc4CgTHGVDgLBMYYU+EsEBhjTIWzQGCMMRXOAoExxlQ4CwTGGFPhLBAYY0yFs0BgjDEVzgKBMcZUOAsExhhT4SwQGGNMhbNAYIwxFc4CgTHGVDgLBMYYU+EsEBhjTIWzQGCMMRXOAoExxlQ4CwTGGFPhihoIRKReRB4VkV0islNE1opIo4g8JSIvuz8bilkGY4wxYyt2jeBu4L9U9VzgImAncCvwjKquAJ5xnxtjjCmRogUCEZkJvBH4NoCqxlW1F1gP3Ofudh/w7mKVwRhjzOkVs0awFOgE/k1EXhSRb4lIDTBXVQ+7+xwB5uZ6s4jcKCJbRGRLZ2dnEYtpjDGVrZiBIAC8Gvi6ql4MDDCqGUhVFdBcb1bVe1V1jaquaWpqKmIxjTGmshUzEBwCDqnqc+7zR3ECw1ERmQfg/uwoYhmMMcacRtECgaoeAQ6KyEp305uBHcAm4Dp323XAxmKVwRhjzOkFivz5nwIeFJEQ8ArwcZzg84iI3AC0Ah8ochmMMcaMoaiBQFW3AWtyvPTmYh7XGGPM+BW7RmDMpGtp2cfjT2ymrb2L+c2zuXrdWlatWlrqYhkzZVmKCVNWWlr28fVvbqK3t5+5cxro7e3n69/cREvLvlIXzZgpywKBKSuPP7GZ2poIdXXViAh1ddXU1kR4/InNpS6aMVOWBQJTVtrau6itjYzYVlsboa29q0QlMmbqs0Bgysr85tn090dHbOvvjzK/eXaJSmTM1GeBwJSVq9etpX8gSl/fIKpKX98g/QNRrl63ttRFM2bKskBgysqqVUv55E3vor6+lqMdx6ivr+WTN73LRg0ZMwYbPmrKzqpVS+3Cb0werEZgjDEVzgKBMcZUOAsExhhT4ayPwBScpXgwZnqxGoEpKEvxYMz0Y4HAFJSleDBm+rFAYArKUjwYM/1YIDAFZSkejJl+LBCYgrIUD8ZMPzZqyBTUcIqH7FFD1370LTZqyBSMjUorPFHVUpfhtNasWaNbtmwpdTGMMSU2PCqttiZCbW2E/v4o/QNRyyflQUS2qmqu5YJHsBqBmdLs7s9kyx6VBmR+Pv7EZvu7OAMWCMyUlX33lz0nwe7+Kldbexdz5zSM2JY9Ks1uHCbGOovNlGVzEsxoY41Ks8mME2eBwExZNifBjDbWqDS7cZg4CwRmyrI5CWa0sRYeshuHiStqH4GI7Af6gBSQVNU1ItIIPAwsAfYDH1DVY8Ush5merl63lq9/cxPAiBEi1370LSUumSklr4WH5jfPpre3P9OBDHbjMF6TUSP4fVVdnTWE6VbgGVVdATzjPjfmFLbspMmHTWacuKLOI3BrBGtUtStr227gClU9LCLzgJ+q6sqxPsfmERhjxsNGDY003nkExQ4E+4BjgALfVNV7RaRXVevd1wU4Nvx81HtvBG4EWLRo0WtaW1uLVk4zMfY/nTFT21QJBPNVtU1E5gBPAZ8CNmVf+EXkmKo2eH4IViOYDPle1G2GpzFT33gDQVH7CFS1zf3ZAfwQuBQ46jYJ4f7sKGYZzOlNZPy1DdUzpnwULRCISI2I1A0/Bt4KtACbgOvc3a4DNharDGZ8JnJRt6F6xpSPYg4fnQv80OkGIAB8V1X/S0ReAB4RkRuAVuADRSyDGYfTTdvPxYbqGVM+ihYIVPUV4KIc27uBNxfruCZ/E7moX71uLbff8TA7dh4gnkgQCgZpaKi1Mf7GTEM2s9icwfhrBQFx/uM8N8ZMO5Z9dBrwGtFTqOGbE1lM5vEnNjPvrFmcs2JhZltf36ClAzZmGrIawRTnNaJn48ZnS5pp0TqLjSkfFgimOK8RPQ88+HTBhm9OZPioJYQzpnxYIJjivO68u3tOFOyOfCLDRy2vizHlwwLBFOd15z2rcUbB7sgn0sxjCeGMKR/WWTzFeaZivuZKfvz0llO3T2D45kTnBHilAzbGTC9WI5jivO6816+/vGB35NbMY0xlK2rSuUIpp6RzUzVj51QtlzFm4qZE9tFCKZdAYBk7J4cFNWMc4w0E1kdQJLkuRtmjc4DMT5uEVTjZwTZ7KKwFW2O8WSA4A2PN+M11Merrj7Ji2fwRn2GTsArLgq0x+bPO4gkaaxKW17j8wYEhm4RVZDbj2Zj8WSCYoLEmYXldjGpqwjY6p8hsxrMx+RtXIBCR20VkhogEReQZEekUkY8Wu3BT2Vh3nl4Xo3NXLrRJWEVmQ2GNyd94+wjeqqq3iMgfAPuB9wA/A75TrIJNdWNNwvKcBOZm9LQLf/FMJJOqMZVuvIFgeL93At9X1ePuymMV63QXe7sYlY4FW2PyM94+gidEZBfwGuAZEWkCYsUr1tRnuXaMMeVi3BPKRKQROK6qKXcx+jpVPVLU0rmm24QymzhmjJkKCjqhTESqgf8NLAJuBJqBlcATZ1LIcmVj2QvHZgkbU3zj7SP4N2Ar8Dr3eRvwfSwQ5NTW3sXcOQ0jtk3nseyluhjbLGFjJsd4+wiWqertQAJAVQdxlys3p5rfPJu29i5e3LaHX/yyhRe37clcRKebiaxeVigTWTDHGJO/8dYI4iISARRARJYBQ0Ur1TR3/nmLefKpLYSrgoQjVQwMROnuOcHbr7p0Qp+3ceOzPPDg03T3nGBW4wyuveZK1q+/fFLu1EvZzFVuNStjpqrx1gi+BPwXsFBEHgSeAT43njeKiF9EXhSRJ9znS0XkORHZIyIPi0hoIgWfynbsbKW5eRbRaJzDh7uJRuM0N89ix87WvD9r48Znueuex+jvH2TmzBr6+we5657H+Of/96NJuVMvZcoGmyVszOQYVyBQ1R/jTCK7HvgesEZVfzLOY3wG2Jn1/KvAXaq6HDgG3DDu0k4TO3cf4MiRHlLpND6fkEqnOXKkh527D+T9WQ88+DThqiDV1WF8IlRXhwlXBXnk0f8hlUqzZ287z27ezp697aRS6YI3m5TyYmyzhI2ZHONNMfGMqnar6r+r6hOq2iUiz4zjfQtwJqF9y30uwJuAR91d7gPePbGiT109PX309Q0yNJQgkUgxNJSgr2+Qnp6+vD+ru+cE4UjViG1Oc1OMffsOE48niESqiMcT7Nt3mF27DxbqNIDSXoxtroYxk2PMPgIRCQPVwGwRaeBkB/EMYL7nG0/6B+AWoM59PgvoVdWk+/yQ1+eIyI04Q1VZtGjROA41dfT1DZJKpvH5fYgIqko6laavbzDvz3IWqR+kujqc2RaLDuHzOZ8dDDq/wmAwQCKRZGCgsPP8Sj1L2mYJG1N8p+ssvgn4U5x5A7/O2n4C+Kex3igi64AOVd0qIlfkWzBVvRe4F5wJZfm+v5QSiRTBUNAJAOk0Pp8Pv99PIpHK+7OuveZK7rrnMcCpCcSiQ8SGEjQ1zSSZTBFPJAkG/CSSKRSorqka+wMnwOtibGP8jSkPYwYCVb0buFtEPqWq/5jnZ18OvEtE3gGEcWoRdwP1IhJwawULcOYkTAmFurBFwkEGBocIBvyIBFBVkqk0kXAw789av/5ygBGjhj5505Xs2NnK/tYjdHWdYGAwRk11mPnNs1iy+Ky8jzERNsbfmPIx3uGjx0XkY6M3qur9Xm9Q1c8DnwdwawQ3q+o1IvJ94H3AQ8B1wMZ8C10MhbywXXzxCl54YReptBMA/H4fVaEAF1+8YszjewWh9esvzwSEYcuWNfP1b25i+bLmEWksJqsjtdJnT1ttyJST8Q4fvSTr3xtwhpO+a4LH/Bzw5yKyB6fP4NsT/JyCKuTkpes/dhWLF59FQ30dM2ZU01Bfx+LFZ3H9x67KuX9Lyz5uv+Nhfrl5By/vaeOXm3dw+x0PjzkUtNQdqZW8ElgpJ9kZUwzjqhGo6qeyn4tIPc4d/bio6k+Bn7qPXwEmNrOqiAo5eWnVqqXccvMHx33HuOH+J2lv7yISqaK6OkwykaS9vYsN9z/JHbd/YszjlOoudKz1GMpdpdeGTPmZ6OL1A0BZ/cWPdWGbSDNAPhfp7dv3IyIcPzFAMpEiEPRTFQqyffv+MzmlohprPYZyZzOeTbkZ7zyCx0Vkk/vvCWA38MPiFm1yeY2XP/+8xUVvBojHE/T1R0m7/QnpVJq+/ijxeKJgxyi0UjdNlZLNeDblZrw1gjuyHieBVlU9VITylIzXePnHn9jM8eMDtLTsJzYUJ1wVYuHCpkzfQUFGGUXC9PXHUL+48w4gnVYiEWfuQKk7Jr2OX6lj/Cu5NmTK07gXpimlUi5M8/4PfpnWA0cJ+P34A35SyRTJVIq5cxpYtGhOQRafue7jX+WVfe3EYonMvINwOMjZS5v5i89+wHORGyhMIBqLLbKTW6mDszHjUeiFad6DkyNoDs7sYgFUVWecUSmngY7OXgCSqTTxRBKfz2lNO9rRy/nnLfbsMMznQtHUNJO9r7SP2KbqbPfqmNxw/5NEo0NFH8dvHaO5VWptyJSn8TYN3Q5crao7T7tnmYnHkyNmBKdS6cxjr+GTw8NBjx3rJ55I0Nrawc5dB7jl5g/mvHgc7x0gPpQgEPDj9wdJpdLEhxIc7x0gFovn7Jj85ebtvOrCs4t+gbaOUWPK33gDwdFKDAJj8fl89PdHc44y2nD/k7yyr52hoWSmqedY7wnP4aCH2jqZWV/L0FCcZCJFMOintjbCobZO3viGV+UczQTegWgivGowlTxM1JhKMd5AsEVEHgZ+RNaCNKr6g6KUqoi8Fnnxkkzmzg8kAv0DJy/I2R2Gf/zpexgYGELEaeJRTTEwkOKFLbtzXnABqkIB6rIu7PF4glQq7dkxecH5SzwDUb7GmlVtHaPGlL/xBoIZwCDw1qxtCkyrQDC8yEu4KjhikRfAMxj4fEIg4CedTpNOKz6f4PP58PnEMytn9tBCJ/O20yHf3x/l9jse4ljvAIl4ktYDHezc1cqCBU3s23cYESEQDJBMJInF4qxevZxVq5by1ivXnBK8hlNMQH4X6FyBaKx+gM/f+pGSZh81xhRfRY0aet8HvkxPz3GSqXRm4lbA76OxcSaPPvLXOd/zznd9ga7O4wSD/kxK6UQixeymmfz7pr/P+Z6L19zkWYZZjXWEI1WZjKGx6BBLl85jaCiR6VMIBYM0NNRyy80fBCjYqCGvEUB9/VFWLJvvBi2HqnK04xj/dM+nT/u9FoKNwjGm8AoyakhEblHV20XkHxm+rc2iqpNzlSiQI0d7SKXS+H2SmbgVTaQ4crTH8z033vBOvnbnw6TTiqYVEQiG/Nx4wzs93xMI+D2blMKRKkLuGgLDPw8d6uRrX70p54XwK7d9d8y79Xwull53/keP9hasmWkiLJOpMaV1uqah4Q7iLeQIBNON3+8jmUjhC/gBp9kmmUziD3lPsM6VBvp0/Qp+v5BMerw4ugbmPvcajljIUTten1VTE/bs75gMNkTVmNI63XoEj7sPdwBfAJZkvUcBzzTUU9GcpnpaDxwlmUxlJoelNc2cpvox35crDfRYAoEAyWSadDqNqtOx7PP5QJVYLJ6zL8BLIXMgeX3WuSsXZvoKStEPYENUjSmt8XYWfwf4C+AlIH2afaeEXKODXn3xCqqqghw82JlJF3H22fO44PwlQOHaqc+a25izL6K6JkJjQx3HjvUzOBgjFAzS3DzbMz01eKczeN3aC/JuThlrBJBlMjWmco2rs1hEfqGqr5+E8uSUb2fxxo3Pcsed3yeVSpFWxSeC3+/nQx+8gh07Wz07XrNH9ARDARrqa7jl5g95XiC9hqJmj07KXl7yzz79XpYta8472HiN9Bl98ezrG6S+vpbP3/qRMb+bfJq5JoOlsTCmOMbbWTzeQPBm4MPAM5RgHkG+gWDd+r+ko6OXUDCA+ARNK/FEkjlz6rnt7/5XzgvxzX/xDV7YuptUKk0qlcLv9+P3+7jkNSu5/rqrTnnP3r3tnhf74WBQzAvun3z6HubOachrpM9UvuDaqCFjCq+guYaAjwPnAkFONg1N2XkEnZ29iCjxeCJTIxCfs92rCeTF3+whGh1yJ4ApyWQaEXhh626isVNz+uxvPUq4Kkh1tZMhdPjnAw8+zfr1l7NsWTOvefWKzIVt2bLmgp7jRJpTpnKnrOXuMaZ0xhsILlHVlUUtSYElEs4wUTc7HsmEEgz6PfcfGIg5nch+f+YuO5lM0d8fzXnx7Oo6zrx5s0Z8RjhSRXfPiTGHQ0L+GUNz3S1PZMavdcoaY3IZ75rFvxSR84takgKaOaMGVJ3xriLOT1Vn+xhUyTQNpVJpt3aQO6eP3+8nFh0asT0WHWJW4wzP9Y833P8kt9/xEL/81Q727Gnnl7/awe13PDTmIjde6+MCeS8MYwuqGGNyGW8geC2wTUR2i8hvReQlEfltMQt2JubNm0U4EiKdVpLJFOm0Eo6ETrmDzxZyawuqmmkeAifFRK6L5wXnLyY2lGBwMEZalcHBGLGhBNdec6Xnwu7btu2hvb0bVSUSCaGqtLd3s+G+Jz3L5RVUhhfGyYfXKmzD+Y6MMZVpvE1DbytqKQqsqWkmra1H8YmPVDqN3+cjEPDT1DTTs1NyxoxaBqPxTBqJ4Z+NjXU5J1t95lPvYe/e9hEdwp+8yekQ3rGzNWf7/WA0TmNj3Skzi7fv2O95Ll7NObt2H2TnrtZT8haNNcrJaxW2qdA2b53FxpTOuAKBqrYWuyAFpZBMpRiKJzIjgBA4fnzAc52AWbNmMDgYJZ5IkU6l8fl9hIJ+5p01y/PiuWrV0pwjgbza7yPhoOfM4rHSQO9vPUJX1wkGBmPUVIeZPXsG3d0niMWGCEeqiERCJJKpTO3ijq+dmup62FTslLUUE8aU1nhrBNPKvtYjxGJxAER8qDujd8eO/QSDAZKpdOZi39c3wIb7n+TclQuJREJ0dh7PXHCbmmayeNHcvC+eXnfeG+5/km3b9pwys3jp0nmeF8Lzz1vMk09tyQxTHRiI0t1zgnQqRVW4ihMnBjKT1qqqQmPWLqaqqTyayZhKUJaBoKvrOAG/n1Do5OnF40kGo3FCqTSBQAB/wE86rURjcbZt28M/3PnHfP2bm1i+rHnEXfxE289zBY/rP3YVt3f0njKzeOaMGvx+X84LIcCC+bM5eLCTY739hKtCLFzYxJ69bcQTKQJ+XyaBXl/fIDXVVRMqbynZaCZjSmu8ncV5E5GwiDwvIr8Rke0i8mV3+1IReU5E9ojIwyISKvSx02kFcX5q1nNwcv4MDyv1u2sLDEbjmbv4fEbh5GvVqqXccvMHed3a81mxfD6vW3s+t9z8QWJDceLxBC9u28MvftnCi9v2EI8naGvvYtfug3R29jJjRjXNzbOZMaOazs5e0mn3vBAnmRFCOq1E3PkM04mNZjKmtIpZIxgC3qSq/SISBH4hIv8J/Dlwl6o+JCLfAG4Avl7IA8+d08CxY30kkimSySSBQICaqiCJRBJVJZ1OZzqDVZXqiBOLvJqAvNrvJ9LBmesY4aoQ236zx2nvD4cYiid4qWUfqy9azpEjxxARgm7HcjAYIJFI4vMJ1ZEq4vGkM//B56O2JsysxhkF+hYnj62CZkxpFa1GoI5+92nQ/afAm4BH3e33Ae8u9LGvveZKUuk01ZEQc+Y0UB0JkUqnOWfFAqpCARAhmUqDCFWhABdfvMJ68y8wAAAc5UlEQVTzs7zG8W/c+GzO7WPNCfAkWT9Ehp+CQHVNFfFEko7OXg4f6aajs5d4IklNdRXnnLOAuXMbmDmzhrlzGzjnnAWcu3Jh/scvscmojRljvBW1j0BE/MBWYDnwz8BeoFdVh7P1HwLme7z3RuBGgEWLFuV13FxrCHzyJmd5x+xRQ8MrgQ1n/8yVH2jHzlaSyRR79raPGLXzwINPs+zseQXp4IzF4ixePJc9e9ozWVGXL28mFos7qbP3Hxn+VgBIJVOcc85C/H5fwfo0Sm0qjmYyplIUNRCoagpYLSL1wA9x8hWN9733AveCk3Qu32N7rSFwy80fzNmc47WecU1NmGQyRTAYyDTb7G89ysBAjItedfaIzx5PB2eu5qRwOMSuXQeYMaOaxuAMkokkra1HnXUK1GkOqhu1vOXMGTWnJMObKnMCjDHTy6SMGlLVXhH5CbAWqBeRgFsrWAC0TUYZTueBB5/OmUSus/M4jQ21IyaBJRNJ/H5f3ss7eo2Xjw8lSSbTHOvtz8x7CPj9oBAbinPhqqUcONiZqZGsWNZMbChud9HGmIIoWiAQkSYg4QaBCPAW4KvAT4D3AQ8B1wEbi1WG0caauNTdc4KZM0fmIgpHqkh3n2BgcIjunj7S6TQ+n49wOMRZcxvyXt7Ra7z81t0vI9kdAzgDgTq7jnPuyoX09vZzcdYqZn19g8wZNdzyTNnMXmMqVzFrBPOA+9x+Ah/wiKo+ISI7gIdE5G+BF4FvF+PgXou5eE1cmtU4g/7+wUxNAJwkcsFggMGBGMpwDqIUgwOxTIdmPk0zXuPlB6NDzKiL0NBQl9k+OBhjYCDG1evWcvsdD7Nj54ER/RqFHFFjM3uNqWxFCwSq+lvg4hzbXwEuLdZxwbmw5UolgcCKZSP7pofb9a+95kruuucxgBELzfh8zjBTn98Zqw9KOqXs23+EvXvb2frrl+nuOcGRI8c4/7zFE1ozOBIJoUA8kcz0AyjOiCGHMw9CnP84zwtoojN7rRZhTHkoy5nFG+5/kvb2LiKRKqqrwyQTSdrbuwiHq2iLdOVMI7F+/eUcauvkkUf/h+6ePqqrq/jA+36P+7/zFKGqAOmUuovc+AhUCcd7+7njzkdIpdNoWmk/3MUddz4C4LkSmdd4+YsvWk5f/+CIfELzm2exZPFZPP7EZuadNYtzVpwcFtrXN5iZdVyIC/FEZvZaLcKY8lG0eQSltH37fsQnHD8xyNGjPRw/MYj4hL7+QXbtPsjAQJRwOMTAQJRduw9y/nmLaWnZx3PP72JGXQ2NjXXMqKvhued3Ac6deDgcojpSRTgcQnDmIQzFk4j4CAQDiPgYiif5l2//h2e5vMbLX3/dVQQCfpYva+bytRewfFkzgYCfq9etpa29K+es4127DxZsHsNEZvYWMj22Maa0yrJGEE8k6e+PEQj48fv9pFJpTpyIIqJc9KplJ++8ayIsXjyXHTtb+fmzL+WsRVRXVzEwEEOSKfwBP6lkilTaWa3T5xMSiQTptOJz01V0dvWOWTavkT5e/Q1es46r3A7rQsxjsNXOjKlsZRkIwuEQ/X1RBHWb1zWzYtmC+U0sXDAns6+q0tbexd697YTDoRGpHIazls6Z08DRo8eIxuIEAn7mzm3g6NEe4vEkfr9k0lXE40lCIe/lMMG7Xd1zKKg4S2b29vZnMqYG/D40Fs+5+M1ELsQTWadgImsmG2OmprIMBLNnzWRwIObkGkql8ft9VFdXjTn2f+/edrLGcDrcZG6NDbWICIl4kmAoQEN9Df39UY4fH0DV6UAGJ6jU1TlDUHNd8IG829U7O4/j8/kgmXaX3HQS5wnOXXmu/o5C5UAai+UHMqZ8lGUg8FpboLY24jn2v7v7BC9s3Z1Zs9jv9+P3+6itjeTsrG1v76G2JkxsKJGZX1BbE6Z53izPjtRIuCpnuoqxmnMGBmKEQgHq62sz24aXx9y1++Ap6xRcuGrppHTiTuXVzowx+SnLQDB8tzo6D89wTqFcF683vP5CfvX8jkxq53Q6DSjhcIh4IsmL2/ZkLt6LFs0hGPSzfPmSnHfkXsMxt/z6d5kmpFQqRX9/lGO9fUTdRXRyqa6poq9/8JShpalkipXnLDilv+MnP/1NwXIgnY7NbDamPJRlIDjd3Wqui9eOna2sumDpKRf2Q21dvPTSK0QiVUQiVcTjCV566RWWLp3H4GCMeCIJ6nRQDw46E8C+ce/jOTtSBwZipNNKKBggEAigaaW/P0ZPT5/nuZy3chGpZIpDh7oyCekWLJhNT09fzv6O373cNqEcSMaYylWWgQDyv1sdDhgL5jdltqkqe/Y4qZCc9nkdNZVLcPqj1d3B6WOY3zyb1gNHTwkqzkekGYonUFVEnP6FWHTIs1zDS1VGIiEaGuuIRYdoa+9mxbL5Ofs7nBnS+eVAMsZUtrINBPkuJuN18Q4GA8yb18jevYczd+TLls3j0KFOFi+ey8BAjEQySSgUoLq6isef2Oy5zrDPJ6TTMmL9ehEZc57wjp2tOZuA6uqqc/d3XHMlP356y6nbrRPXGOOhLAOBV2ftW69cw4+f3pLXIvHz5jZy4EAHdTOqaQzUkUimOHCgg0Qixb59hwmFgpkmo337DhONOu39uS7eLdv3o3pyhJEIqDrDXb0CVFt7V84moKMdx3jrlWtOWT9h/frLWbas2TpxjTHjJqqFzVtTDGvWrNEtW7aMe/+v3PZd9rceGZGyYfbsGXR3943oSAVnBNDwiJxc72lr7yYWHSISqSLgpqCORoeIxuLU1UZGJKkbHIxRW1vNWWc1MHdOg9v041BV/vPJF4jH44BTKxA3b9BZcxtZvHgutTWREXfxw/0co8fr9/UNkkqniUaHcr7HLvrGGAAR2aqqa063X1nWCHbuPsChQ52njM6Jx5NjdqR6db6uvmgZBw50ZALE8uXz2bL1d56J4rwmWwWDfiLh2hFrKQcDfvr6o55J35zsow+xY9eBEfMY5jQ1cPz4AC0t+zNNVgsXNhVldJAxpryVZSDo6e6jr89pP1dVksk0Q0MJQqHAmB2puS7esxpnEAoGTlkPYO6cBmbNqsuZKM5rslXT7Hr6+gaorj5Zu4jHE8QTydPMEj61U3rX7gN0dPYS8PsJhYIkEkl+9/IhhoYShf9CjTFlrSyTzvX1D5JKpTIjc1SVVCqFiNA/EKWvbxBVpa9vMLPO79Xr1uZ87dprrvTc7pUoziu53GtevYKlS+cRCgWJRocIhYIsXTqPuXMaPJO+OdlHG7ns0nN5w+UXctml5zLvrEY6Oo87mVADfgQIBPz4xEdH59i5jowxZrSyrBEkEikCAb87SziNzycEAn7SaR1zfoHXa16dr2N1ynoNX8010W2skT5ecxJSqTT+UIB0Op0JduITUql0kb9dY0y5KctAEAj4iMVS+Hw+/H6nwzaZTBEOh8acX+D1Wr7bvZxuotvoEUBjDWutrQ3j9znpsJOJFIGgn1DIR2PjzHGXxxhjoEwDQV1dNYODQ87wzDSIz0nUlt3+n0u+cw8mIlfwaGnZx4+f3sKys+dx0avOpr8/6jxf1uw5rPWKN17E//z8t84aCY0nV1S79porJ1QuY0zlKstAMHvWTE4c7yc2lERJI/iIhAPMnuV9tzyRuQdQmBXCxloqEnLPSQiHQ7z/vW88ZUU1r9XRjDHGS1kGgqammbS2HqUqJKTSafw+H36/P5PmIRevi/EDDz6dM4nbhvufzIzjP9Msn6db5MWZq3Ai81p1dZhduw9SWxvmtZeel+lX2LGzlZaWfTZ81BiTl7IcNYQ6/QT1DbWcNbeB+oZaAgHfmGu+t7V35RzC2d1zIuf27dv3F2ypxvnNs2lr7xqxHOVwLSMcDvHSS68QjydGJL3r7j5hS0UaYwqiLANBbCjOhauWUhUKEo3FqQoFuXDVUmJD3umevdbtHU7iNno7kElPnVlLOJGcUJbP889b7LmW8nDwGp30LhqNFWyFMmNMZSvLpqHhmb2jJ4HNmdPg2fHrueKWx9DOBQua+PWvf0cqraRSafr7o3R19XLJJefmXV6vxHI7drZmgtqBg52ZUUMrljWza/dBzxXKjDEmH2UZCLwu6q9bewG33/EQx3oHSMSTtB7oYOeuVm65+UNjDu3MNV/g7nt+QDQaxx9wVjJLpdLEh1Ic7x3Iu7xeieWGj5crqC1Y0JRzhbK3X3XpmX+BxpiKUrSkcyKyELgfmIvTsnGvqt4tIo3Aw8ASYD/wAVU9NtZn5Zt0DnIPBd1w35O8sGVX5i7e7/fh9wmXrDmXO772ibw+/+3vvJVUKk1sKJHJGxSuCuL3+/jPf78tr3J5JZarr6/NBLXRyeUiYWflstFJ8pYsPovP3/qRvM7FGFOepkLSuSTwWVX9tYjUAVtF5CngeuAZVb1NRG4FbgU+V+iD5xqv/+K2lxmKJwn4fQT8PlSVoXiKF7e9PKFjhEKBEe308XgiM7M3n8Xrh4eowqkzi71qKt+493HPWoQxxuSjaIFAVQ8Dh93HfSKyE5gPrAeucHe7D/gpRQgEuURjCUScyWXgLAoj6TTR2NiJ2nJd1C+4YAnbtu1BRDIJ5GKxOKtXL8978fodO1vzXgjeK8OprURmjMnXpPQRiMgS4GLgOWCuGyQAjuA0HeV6z43AjQCLFi0qSDmqIyH6+qOk005eHk0r6bRSVxsCYOPGZ09J87BsWbPnXXxHRy/HjvUzOBgjFAzS3Dyb6z921ZiL1wcCfoLBAJFwiKF4gv2tR4nG4nz+1o/kvPCfbqIb2EpkxpgzU/ThoyJSCzwG/Kmqnsh+TZ0OipydFKp6r6quUdU1TU1NuXbJ2+rVy4mEqwAn9xBAJFzF6tXL2bjxWe648/u0t3cRjQ7R3t7FHXd+n7vv+UHO8fo7drZy2aXncqJvgJ6ePk70DXDZpedmVhXLNbQzGo0jQCgYQEScn8DggPeaxdlBZfTxc2U4tclkxph8FbVGICJBnCDwoKr+wN18VETmqephEZkHdBSzDNmu/9hVmbv4eCJBKBikoaGW6z92Fbf+5bcYiicIBQP43drCUDzBS9v38c63Xzbic2prI/z6xZfp7DpOJByiocFZVP77j/2MBfObPJttqiNVqCqJRDLTnKSq1NSERxc1o629i2DAz4vb9mSakxYtbKKtvSvvpHfGGJNL0WoE4qzT+G1gp6remfXSJuA69/F1wMZilWG0VauW8t4/eANVVUHi8SRVVUHe+wdvYNWqpXR29hII+PD5BAE3dbWPhDtJbPSs347OXsJVQaqrw/hEqK4OE64K8sCDT3P1urUcPtLDc8/v4ufPvsRzz+/i8JEeVl+0LOd6BOeuXOhZ5nBViJda9jEUT2Sak15q2Ue4KjRZX5sxpswVs0ZwOXAt8JKIbHO3fQG4DXhERG4AWoEPFLEMI4yV5dPv9zN6KK0q+P3+nOP1h4YSNDTUjdg/HKmiu2e49UtBQJz/AMobXn8hP356yynrEQyPKMpJsn6IDD/l5ANjjDkzxRw19Au8L1dvLtZxh3mN1/fK8nnBeYv5zW/3AmQmiKWSKWbMrGHZ2fNOmfW7a/dBYtGhEYvXx6JDzGqc4a4qNotzVpy80+/rG5zQ6KBYLM6FF559yprJsZh3ugxjjMlHWc4s9hpp09cfZcWy+SP2Hc7P85lPv4e//psNHD3aSyKRJBgM0NzcSE1NJOd4/Z6ePjq7jgNOTWB4PYBP3nQlTz2zdcxsovnwTJfRVJ/3ZxljTC5lmXTOa6TN4MCQ59rAq1Yt5ctfvJ53vv0yLlmzkne+/TK+/MXrOW/lopzvefXFK/izT7+X2tpqjh8foLa2mj/79HtZv/5yzwR24XCIr39zE729/SMCVEvLPs9z8VpLeczmJGOMyUNZ1gi88vvX1IQ5fKSbHbsOkIgnCYYCNNTXZMbej7XO8PBnjJ71m2shGK9cR5FwlWfT1FjLZ+bbnGSMMfkoy0DgNXyzqWkmHR29oKDOfzhdr+tELsSrVi3lrVeuOWVy2kSbjGyYqDGmmMoyEIx1Rz7vrEbOWbEgs29f3+CYd+SQ/4W4pWUfj/3wZwzFE1SFggzFEzz2w58xp6mB/v6opYUwxkwpZdlHMHwXP3rWbWwoPimLuWy470na27tRVSKREKpKe3s3x08MWHu/MWbKKcsaAeS+i5+sRG3bd+zH5/Nx4sTgiBTVhw518rWv3pSzmclrwRxjjCm2sg0EuS6snquQuZ3FhboYJxJJ+vqjBAIB/H4/qVSa4ycGqKuN5AxQXsNdLXeQMWYylGXT0PCFdfQwTcAzUZvXe1pa9tHSso+v3PZd/uTT9/CV27475nBPcOYVqJ7shhacWcrhSFXO/b2Gu9pC9MaYyVCWNYKxZhB7pXv2es+G+58kGh3K6269sbGOwcEY8XiSZDKJ3++ntjZMY2Ndzv29hrvaIjPGmMlQloGgrb2LYDAwMmPnojljXli93rNr1wEWLZpzymIyY400Om/lIpLJFG2HukgkkoSCPs46q5HzVuZeV8EWmTHGlFJZBoJwOMS2bXuIRKqIRKqIxxO89NIrrM5K05DrPZs3byeeSJFOp+n2+ThypJtAMMD+1qM5F5Px6lM4/7zFPPnUFiKREA2NTorq9vZuPvi+K3Ie+3R9F8YYU0xlGQiGl7pRANWTK9/kXALH0Xaoi8Fo3ElDLUI6nWYwmiKQSBIJhwgFna8q5K4j0NPT59nBu2NnK+euXEhn5/FMorolS85ix85Wli1rzhk8bPawMaZUyjIQxIbiXLhqKQcOdmaac1YsayY25H0X33qwg1AoiKqSTqfx+/0EAkLCXTxm9GIyMbffIFc/xPBnL5h/cmU1VWXX7oMcauv07G+wC78xphTKctTQ/ObZhEJBLl69nNe/bhUXr15OKBQcM+lbKpXC7xdnsZlIFeGqIH6/IELOxWSCwYDn5DSvpHMDAzEbHWSMmXLKMhB4ZexE8bwQNzXVk0ymSaedpqR0Wkkm08yaNQO/38fyZc1cvvYCli9rxu/3ccH5SzwzmXodv7qmalJmNhtjTD7KMhBMJMXEH/3hO6gKBUmn0yQSSdLpNFWhIH/8ifU5P+v6667yTBfhdXyvlNY2OsgYU0pl2UcA3ikmWg8cPdmJWx2mqWkmixfNzaSTHp0xdHh7rvb7sTp4801pbYwxpSKj1+mditasWaNbtmw548/ZuPFZ7rrnscz6w8Orig0vKDMZLKeQMWayiMhWVV1zuv3KtkaQy46draw8Z8Ep6w/v2Nk6aYHARgcZY6aaigoEbe1dOdcfts5aY0wlK8vOYi9ewzqts9YYU8kqKhDYQvDGGHOqojUNici/AuuADlVd5W5rBB4GlgD7gQ+o6rFiHN+rU9ZSORhjzEhFGzUkIm8E+oH7swLB7UCPqt4mIrcCDar6udN9Vr6jhlpa9nH7HQ9xrHeARDxJMBSgob6GW27+kF30jTEVY7yjhorWNKSqPwN6Rm1eD9znPr4PeHcxju21ZvCG+54sxuGMMWZam+xRQ3NV9bD7+AgwtxgH2b5jP+FI1YiMocPbjTHGjFSy4aOqqiLi2S4lIjcCNwIsWpR7QZfTHCDnc5vQZYwxI032qKGjIjIPwP3Z4bWjqt6rqmtUdU1TU5PXbjldcMESYrG4k0IaZzH5WCzOggVNntlHjTGmUk12INgEXOc+vg7YWIyDXP+xq5hZX8ux3n7a27s41tvPzPpaZs6osTTQxhgzStECgYh8D9gMrBSRQyJyA3Ab8BYReRm40n1eFDXVYRrq66ivr6Whvo6a6jAdXb2WBtoYY0YpWh+Bqn7Y46U3F+uYwx5/YjPzzmrknBULMtv6+gbZ+8ph+vujtki8McZkKcuZxW3tXTnv/Gtqwjaz2BhjRinLQOCVU+jclQtzLhhjo4aMMZWsLLOPXr1urecCMJYG2hhjRirLGoHXUpEWAIwx5lRlWSMAWwDGGGPGqyxrBMYYY8bPAoExxlQ4CwTGGFPhLBAYY0yFs0BgjDEVrmgrlBWSiHQCrafZbTZQqUmDKvncobLP3869co3n/Ber6mnTN0+LQDAeIrJlPEuylaNKPneo7PO3c6/Mc4fCnr81DRljTIWzQGCMMRWunALBvaUuQAlV8rlDZZ+/nXvlKtj5l00fgTHGmIkppxqBMcaYCbBAYIwxFW7aBwIReZuI7BaRPSJya6nLU2wi8q8i0iEiLVnbGkXkKRF52f3ZUMoyFouILBSRn4jIDhHZLiKfcbeX/fmLSFhEnheR37jn/mV3+1IRec79+39YREKlLmuxiIhfRF4UkSfc55V07vtF5CUR2SYiW9xtBfu7n9aBQET8wD8DbwfOBz4sIueXtlRFtwF426httwLPqOoK4Bn3eTlKAp9V1fOB1wJ/7P6+K+H8h4A3qepFwGrgbSLyWuCrwF2quhw4BtxQwjIW22eAnVnPK+ncAX5fVVdnzR0o2N/9tA4EwKXAHlV9RVXjwEPA+hKXqahU9WdAz6jN64H73Mf3Ae+e1EJNElU9rKq/dh/34VwU5lMB56+Ofvdp0P2nwJuAR93tZXnuACKyAHgn8C33uVAh5z6Ggv3dT/dAMB84mPX8kLut0sxV1cPu4yPA3FIWZjKIyBLgYuA5KuT83aaRbUAH8BSwF+hV1aS7Szn//f8DcAuQdp/PonLOHZyg/2MR2SoiN7rbCvZ3X7YrlFUqVVURKesxwSJSCzwG/KmqnnBuDh3lfP6qmgJWi0g98EPg3BIXaVKIyDqgQ1W3isgVpS5PibxeVdtEZA7wlIjsyn7xTP/up3uNoA1YmPV8gbut0hwVkXkA7s+OEpenaEQkiBMEHlTVH7ibK+b8AVS1F/gJsBaoF5HhG7py/fu/HHiXiOzHaf59E3A3lXHuAKhqm/uzA+cm4FIK+Hc/3QPBC8AKd/RACPgQsKnEZSqFTcB17uPrgI0lLEvRuO3C3wZ2quqdWS+V/fmLSJNbE0BEIsBbcPpIfgK8z92tLM9dVT+vqgtUdQnO/+P/rarXUAHnDiAiNSJSN/wYeCvQQgH/7qf9zGIReQdO+6Ef+FdV/bsSF6moROR7wBU4KWiPAn8N/Ah4BFiEk677A6o6ukN52hOR1wM/B17iZFvxF3D6Ccr6/EXkVTgdgn6cG7hHVPVvRORsnLvkRuBF4KOqOlS6khaX2zR0s6quq5Rzd8/zh+7TAPBdVf07EZlFgf7up30gMMYYc2ame9OQMcaYM2SBwBhjKpwFAmOMqXAWCIwxpsJZIDDGmApngcBUNBH5tIjsFJEHR21fIyL3FOgY14vIPxXis4wpBksxYSrd/wauVNVDwxtEJKCqW4AtpSuWMZPHagSmYonIN4Czgf8UkeMi8oCIPAs8ICJXZOW9r3HXgXjezYe/3t1+vYj8QET+y80Jf3vWZ39cRH4nIs/jpEgY3v5+EWlx1xX42eSesTG5WY3AVCxV/YSIvA34feBPgKtxkntFRyU3+0uctAZ/6KZ5eF5EnnZfW42TBXUI2C0i/4izbsKXgdcAx3FSIbzo7v9F4Co3gVh9cc/QmPGxGoExJ21S1WiO7W8FbnVTQP8UCONM6wdnYZDjqhoDdgCLgcuAn6pqp7tOxsNZn/UssEFE/ggnXYQxJWc1AmNOGvDYLsB7VXX3iI0il+HUBIalOM3/U24t5DKcRVa2ishrVLX7DMpszBmzGoExp/ck8Ck3+ykicvFp9n8O+D0RmeWmzX7/8AsiskxVn1PVLwKdjEyjbkxJWI3AmNP7vzgZbn8rIj5gH7DOa2dVPSwiXwI2A73AtqyXvyYiK3BqGc8AvylWoY0ZL8s+aowxFc6ahowxpsJZIDDGmApngcAYYyqcBQJjjKlwFgiMMabCWSAwxpgKZ4HAGGMq3P8PsNuTO1r1SbAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter([x_i[1] for x_i in x], y, color='#30336699')\n",
    "plt.title('Time on site vs. number of friends')\n",
    "plt.xlabel('friends')\n",
    "plt.ylabel('minutes')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### People who work more spend less time on the site"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xt8nFd16P3fmvtI8l2yY8l3O9hxHGK3TsCkKZRcATuhlAY4EBwOp6GcFtIWTgqcc1rojZDmlCa0h8M9Jg0EyuV17AYSEsgVk+CQhMi3Jo4t2ZZ1tWRJo9Fc1/vH84wsGTuWbe15rJn1/XyUmXk0o71HVmY9z95rry2qijHGmOoVCroDxhhjgmWBwBhjqpwFAmOMqXIWCIwxpspZIDDGmCpngcAYY6qcBQLzG0TkUyLy1aD7MVFEZFBElgTdj3ITkZtE5MlxPvfTIvJvrvtkzk2RoDtgyk9EBkc9rAEyQMF//CFV/Yfy98odVa0r3ReRu4GDqvq/guuRMecWCwRV6LgPxv3Af1PVh4PrkZloInJO/r8tIhFVzQfdDzOWDQ2Z3zB6mEBEFomIisgHROSAiPSKyB+LyCUi8msR6RORfznu9f9VRHb5z31QRBa+SlvXicgO/+c8KiIXjPrefhH5uN/OURH5jogkTvJzlonIY/7zukXkO6O+p/73bwbeC9zqDxdt8b/fKCLfF5EuEdknIh89SRuvE5F2EQmPOvb7IvJr//6lIrJdRPpFpENE/mkcv+uEiKRFpN5//D9FJC8iU/3Hfysi/+zfnyYi3/T72SIi/0tEQv73bhKRp0Tk8yLSA3z6BG39o4g8KSLTTtKdmP/zB/x/k7WjXnuB/+/T53/vulHfe1RE/tuox2OGpPzf/5+IyEvAS+L5vIh0+r+rF0Vk1al+V8YdCwRmvF4HnA+8C/hn4H8CVwIXAjeIyBsBROR64FPAO4AG4Ang2yf6gSLyGv97f+Y/9wFgi4jERj3tBuBaYDHwWuCmk/Tvb4GHgBnAPOALxz9BVb8M3Avcrqp1qrrB/yDdArwANAFXAH8mItec4PVPAyngzaMO/xfgW/79O4E7VXUqsBT47kn6OvpnDgO/BN7oH3oj0AJcNurxY/79LwDTgCX+8fcDHxj1414HvALMAf6+dFBEQiLyFbzf39WqevQk3bkOuA+YDtwP/Iv/+ije7+ghYDbwEeBeEVl+qvc3ytv9/q0ErgZ+F3iN/35uAHpO42eZCWaBwIzX36rqsKo+hPdh+G1V7VTVQ3gf9mv85/0x8FlV3eUPAfwDsPokVwXvAv5DVX+iqjngDiAJvGHUc+5S1TZVPYL3YbT6JP3LAQuBRr+f45okBS4BGlT1b1Q1q6qvAF8B3n2S538beA+AiEwB3sqxQJcDlolIvaoOquovxtmHx4A3+sM5rwXu8h8n/P497l+FvBv4pKoOqOp+4P8AN476OW2q+gVVzatq2j8W9fs3E9igqkOv0o8nVfUBVS0A9wAX+8dfD9QBt/m/o58CW0u/h3H6rKoe8fuVA6YAKwDx/1YOn8bPMhPMAoEZr45R99MneFyad1gI3OkPIfQBRwDBO9s+XiPe2S8AqloEDhz33PZR94dGtXO8W/12nvGHLv7rKd/Rsf42lvrr9/lTeGfVJ/It4B0iEse76vmVqpbewwfxznJ3i8gvRWT9OPvwGPAm4LeAF4Gf4J3xvx54WVV7gHq8D/WWUa9rYezv6sAJfvYy4HrgM6qaPUU/jv9dJ/zg1Agc8P99Ttb2qYz0zQ8k/wL8K9ApIl8uDYWZYFggMBPtAF7m0fRRX0lV/fkJntuG90EMgIgIMB84dLqNqmq7qv6RqjYCHwL+r4gsO9FTT9Dffcf1d4qqvvUk7ezE+xB8C2OHhVDVl1T1PXjDJ58DvicitePo/s+B5cDvA4/5bSzAu9ooDQt1c+yqp2QBY39XJyolvAtv+OhHpzmUM1obML80H3GCtlN42Wcl553gZ4zpm6repaq/jTdU9Brgf5xh38wEsEBgJtr/Az4pIhfCyATnH57kud8F3iYiV/jj0B/DS2U9UdB4VSLyhyIyz3/Yi/fBUzzBUzvwxthLngEGROQvRSQpImERWSUil7xKc98CbsEb5/73UX14n4g0+GfOff7hE/VhDH+45lngTzj2wf9zvGG2x/znFPB+X38vIlP8oba/AE6Z+6+q38a7ynlYRJae6vkn8DTeFcKtIhIVkTcBG/DmEwCex7tKqvGD7wdf7YeJl2jwOv/fPAUMM47fk3HHAoGZUKr6Q7yz4ftEpB9oxjt7PtFz9wDvw5sE7cb7cNkwjiGME7kEeFq8NRL3A7f44/3H+xqw0h8G+v/8D9j1eHMP+/x+fBVvEvNkvo03dPNTVe0edfxaYIffhzuBd5fG6v0spctf5Wc+hjf088yox1OAx0c95yN4H5yvAE/iBaSvv8rPHKGqm4C/AX4qIovG85pRr83i/du8Be/383+B96vqbv8pnweyeEF2E96E/KuZijcP04t3ddUD/OPp9MlMLLGNaYwxprrZFYExxlQ5CwTGGFPlLBAYY0yVs0BgjDFV7pwsTHW8+vp6XbRoUdDdMMaYSeXZZ5/tVtWGUz1vUgSCRYsWsX379qC7YYwxk4qItJz6WTY0ZIwxVc9ZIBCR5SLy/KivfhH5MxGZKSI/EZGX/NsZrvpgjDHm1JwFAlXdo6qrVXU18Nt4S9R/CHwCeERVzwce8R8bY4wJSLmGhq4A9vpVGq/HW4aOf/v2MvXBGGPMCZQrELybYzXb54yqPd7OScr9isjN4u32tL2rq6scfTTGmKrkPGvI323qOuCTx39PVVVETljsyN9N6ssAa9eutYJI5oSam/exZes2DrV109RYz4b161i1anHQ3TJmUinHFcFb8DbvKG1k0iEicwH8284y9MFUoObmfXzxS/fT1zfInNkz6Osb5Itfup/m5n1Bd82YSaUcgeA9jN2z9n5go39/I7C5DH0wFWjL1m3U1SaZMqUGEWHKlBrqapNs2bot6K4ZM6k4DQT+7kxXAT8Ydfg24CoReQlv8/PbXPbBVK5Dbd3U1SXHHKurS3KorfskrzDGnIjTOQJVTQGzjjvWg5dFZMxZaWqsp69vkClTju2SODiYpqmxPsBeGTP52MpiM2ltWL+OwVSagYEhVJWBgSEGU2k2rF8XdNeMmVQsEJhJa9WqxXz4Q9cxfXodHZ29TJ9ex4c/dJ1lDRlzmiZF0TljTmbVqsX2wW/MWbIrAmOMqXIWCIwxpspZIDDGmCpngcAYY6qcBQJjjKlyFgiMMabKWSAwxpgqZ4HAGGOqnAUCY4ypchYIjDGmylkgMMaYKmeBwBhjqpwFAmOMqXIWCIwxpspZIDDGmCpngcAYY6qcBQJjjKlyFgiMMabKWSAwxpgqZ4HAGGOqnAUCY4ypck4DgYhMF5HvichuEdklIutEZKaI/EREXvJvZ7jsgzHGmFfn+orgTuDHqroCuBjYBXwCeERVzwce8R8bY4wJiLNAICLTgN8FvgagqllV7QOuBzb5T9sEvN1VH4wxxpyayyuCxUAX8A0ReU5EvioitcAcVT3sP6cdmHOiF4vIzSKyXUS2d3V1OeymMcZUN5eBIAL8FvBFVV0DpDhuGEhVFdATvVhVv6yqa1V1bUNDg8NuGmNMdXMZCA4CB1X1af/x9/ACQ4eIzAXwbzsd9sEYY8wpOAsEqtoOHBCR5f6hK4CdwP3ARv/YRmCzqz4YY4w5tYjjn/8R4F4RiQGvAB/ACz7fFZEPAi3ADY77YIwx5lU4DQSq+jyw9gTfusJlu8YYY8bP9RVBYJqb97Fl6zYOtXXT1FjPhvXrWLVqcdDdMsaYc05Flphobt7HF790P319g8yZPYO+vkG++KX7aW7eF3TXjDHmnFORgWDL1m3U1SaZMqUGEWHKlBrqapNs2bot6K4ZY8w5pyIDwaG2burqkmOO1dUlOdTWHVCPjDHm3FWRgaCpsZ7BwfSYY4ODaZoa6wPqkTHGnLsqMhBsWL+OwVSagYEhVJWBgSEGU2k2rF8XdNeMMeacU5GBYNWqxXz4Q9cxfXodHZ29TJ9ex4c/dJ1lDRljzAlUbProqlWL7YPfGGPGoSKvCIwxxoyfBQJjjKlyFgiMMabKVewcgakOVkrEmLNnVwRm0rJSIsZMDAsEZtKyUiLGTAwLBGbSslIixkwMCwRm0rJSIsZMDAsEZtKyUiLGTAzLGjKTVqmUyOisoRvfd5VlDVUYywxzT1Q16D6c0tq1a3X79u1Bd8MYU2alzLC62iR1dUkGB9MMptJWO2ycRORZVT3RdsFj2BWBOWt2xmZcGZ0ZBozcbtm6zf7GJpAFAnNWRp+xjc7ltzM2MxEOtXUzZ/aMMcfKmRlWLSc5Nllszorl8huXgswMq6YFixYIzFmxXH7jUpCZYdV0kmOBwJwVy+U3LgW5yVQ1neQ4nSMQkf3AAFAA8qq6VkRmAt8BFgH7gRtUtddlP4w7G9av44tfuh9gTFbHje+7KuCemUoR1CZTTY319PUNjkxQQ+We5JTjiuD3VHX1qBSmTwCPqOr5wCP+YzNJ2bagplJV04JFp+sI/CuCtaraPerYHuBNqnpYROYCj6rq8lf7ObaOwBgThMmeNTTedQSuA8E+oBdQ4Euq+mUR6VPV6f73BegtPT7utTcDNwMsWLDgt1taWpz1s1JM9j9aY8zEOlcCQZOqHhKR2cBPgI8A94/+4BeRXlWdcdIfwuS7IgjiA9lWYBpjjjfeQOB0jkBVD/m3ncAPgUuBDn9ICP+202Ufyi2o3ONqSnUzxkwsZ4FARGpFZErpPnA10AzcD2z0n7YR2OyqD0EI6gO5mlLdjDETy2X66Bzgh940ABHgW6r6YxH5JfBdEfkg0ALc4LAPZRfUkvhqSnUzxkwsZ4FAVV8BLj7B8R7gClftlgQ1cRrUB7Ll8xtjzlRFriwOskZIULnHls9vjDlTFRkItmzdRj5f4OW9bTy1bQcv720jny+UZeLUPpCNMZNNRZah3rWnlY6OXqLRCMlEjEw2x/6WDtLD2bK0H8SS+GotB21rJ4w5exV5RTCUyiBALBpBRLxb/3ilqsb00WoqE2yMSxUZCGprE6gquVweBe9WldraRNBdc6Ya00erMfgZ40JFBoIVy+ezePFcYrEo6XSGWCzK4sVzWbF8ftBdc6Yay0FXY/AzxoWKDAQb1q8jHA6xbGkjl627kGVLGwmHQxVZNbCkmiolllRj8DPGhYoMBNWYuVON77kag58xLlRk1hAEt5mFKZ9S8BudNXTj+66yf3djTpPT6qMTxaqPjq9Nqz5qjBltvNVHK/aKIKj88qDy+Udn0AAjt1u2brNAYIx5VRU5RxBkfnlQq5otg8YYc6YqMhAEmV++a08r+1s6yGRzY1Y179rT6rRdy6AxxpypigwEQZ4dB7Wq2TJojDFnalyBQERuF5GpIhIVkUdEpEtE3ue6c2cqyLPjoFY1V2P6qDFmYox3svhqVb1VRH4f2A+8A3gc+DdXHTsbQdbmX7F8PslkjK6uo6SGhqmtSdDUVM/CBXOct20ps8aYMzHeoaFSwHgb8O+qetRRfyZEkGfH1biq2RgzuY33imCriOwG0sCHRaQBGHbXrbMX1NnxqlWLufrKtdxz78P0HOln1syp3PjeK+1M3RhzzhrXFYGqfgJ4A7BWVXPAEHC9y45NVs3N+3jo4e0sXTKXN79pNUuXzOWhh7dbaWRjzDlrXFcEIlID/HdgAXAz0AgsB7a669rZCWpBmS3sMsZMNuMdGvoG8CzeVQHAIeDfOUcDQXPzPm6/4z56+1LksnlaWjvZtbuFWz/+bucfxofaupkze8aYY7awyxhzLhvvZPFSVb0dyAGo6hAgznp1lu7e9CBtbT2oKslkDFWlra2Huzc96LztpsZ6DrV189zzL/Pkz5t57vmXR65KjDHmXDTeK4KsiCQBBRCRpcA5u+/jjp37yReKdHb2UigUCYdD1NQk2LFzv/O2V16wkAd+9AyFQoGiKv1HU3R29vGWay513rYxxpyJ8QaCTwM/BuaLyL3AZcAHxvNCEQkD24FDqrpeRBYD9wGz8IabblTVCd1VPpVKM5gaJhQSRIRisUh/f4pioTCRzZzQE0+9SCgESggtFAmFQ4RC3vHrr7/MefvVxjavN+bsjTdr6CG8RWQ3Ad/Gyx762TjbuAXYNerx54DPq+oyoBf44Lh7O05FBVUoFqFQKFIs+o/LUHF7x479xOMxQqEQIhAKhYjHY+zYsd9941XGNq83ZmKMN2voEVW9AviPExx7tdfNw1uE9vfAX4iIAG8G/ov/lE14VxtfPP2un1yhUASgtNdC6bZ03KVcrsDAYJpoJEw4HKZYKHK0f4gpx9U+qiSWoWXM5PaqVwQikhCRmUC9iMwQkZn+1yKgaRw//5+BW4HSJ/AsoE9V8/7jgyf7OSJys4hsF5HtXV1d42jqmFgsQjTqfYXDoZH7sZj77Rfiiah3+YGOuY0nos7bDkKQZ+VWetuYiXGqoaEP4Y3jrwB+5d9/FtgM/MurvVBE1gOdqvrsmXRMVb+sqmtVdW1DQ8NpvXZ2w3RAiYRDJJNxIuEQoP5xt2bNmkpNbRxEyBeKIEJNbZxZs6Y6bzsIW7Zuo1Aojtl/oVAolqXkt5XeNmZivOopsqreCdwpIh9R1S+c5s++DLhORN4KJICpwJ3AdBGJ+FcF8/DWJEyo31pzPvF4lAMHuhjOZEnEYyxZMpcLVy6a6KZ+wwXLF5BMxOju7h8pOldfP5VFC89z3nYQdu85QHt7D7FYlGQyTjabY9++w6TTEzr/f0Ib1q/j9jvuY+fuVnLZPNFYhBnTa8tSXNCYSjLesZKjIvL+4w+q6jdP9gJV/STwSQAReRPwcVV9r4j8O/BOvMyhjXhXFxNqw/p1HDzURVNj/Zjqo+Uo/FaqfLpsaWPZ2w5irD6VGkZEiEa9P6VoNEIulyeVKlcpKvFG37z/cA4vb5kwlillJtp4F5RdMurrcrwJ3uvOsM2/xJs4fhlvzuBrZ/hzTirI6qNBtR3UWH1NbZxsLk9nVx+H23vo7Oojm8t7w2OObdm6jbnnzeR1l67g8ssu4nWXrmDueTPLMiwVFMuUMi6M64pAVT8y+rGITMc7ox8XVX0UeNS//wrgfHVVkLX5g2g7qAya2fXTadnf7j/yzsYL+QKz693Px1RjOQ/LlDIunGkaTQo4p//qqu3yObAPRfGGg6Ym40SiEfK5POl0piwjNE2N9bS0dozZBKihYVpZNgEKSjUGP+PeeLeq3CIi9/tfW4E9wA/ddu3MVePlc1AZNMPDWS66aAmxWJR0OkMsFuWii5YwPOx+snjlBQvZvecAqVSaRCJGKpVm954DrLxgofO2g2KZUsaF8V4R3DHqfh5oUdWDDvozIbZs3UY+X+DlvW1jMnfKdfkcxNXIhvXr+MtPfYXOzj4KhQLhcJjZs6fzuX/4I6ftNjXW09c3yJrVy0aODQwMlSVVd+euFpa/Zt6xDK3aJAsXzmHnrpaKLecR5DaspnKNt8TEY6O+njqXgwDArj2t7G/pIJPNkUzEyGRz7G/pYNeeVudtl0pg//wXO3n55TZ+/oud3H7Hfc6vRh57/AU62o9QLBYRvPpKHe1HeOzxF5y2u2H9OgZTaQYGhlBVBgaGypYldaitm5qaxJhjNTWJih4mCTIRwlSu8ZaYeAdejaDZeKO/AqiqnpOrpIZSGQSI+SmNMX/seijlvmBqqQR2IhknmYyRyxdGSmDf8Y9/7Kzd737vMSKRsauns9k83/3eY/zJf3+7s3aD3JozkYjxi6d3ksnkKRaL9IRCdHQe4fWvW+m87SAFmQhhKtN400dvB65T1WmqOlVVp5yrQQCgtjaBqpLL5VHwblWprU2c8rVna8fO/UgoRH9/io6OXvr7U0go5LwE9tBQBgnBcCbLUDrDcCaLhLzjLjU37+P7P3ycTDZHPBYlk83x/R8+Xp4SEwe7SaUy3lWQX2U2lcpw6GDlXhEY48J4A0GHqu469dPODSuWz2fx4rljJjAXL57LiuXznbedzebo6xsknc6QzXoZNH19g2SzOaftxmIRhodz5PNFCoUi+XyR4eGc8/pKd296kJbWTnp7BzjaP0hv7wAtrZ1l2QSo5UAn0WiYcMj7Mw6HQkSjYVoOdDpv25hKMt5AsF1EviMi7xGRd5S+nPbsLGxYv46hIe+DWFXJZvMMDWXKMm4dCofJ5wtevTnx6s7l8wVC4bDTdhvqpwG/WXG1dNyV5154mVQqQyaTI5crkMnkSKUyPPfCy07bBSgUCqgqhaIX/ArFove4DPtOGFNJxnu6OBUYAq4edUyBH0x4jyZIamiY3r4Bcrm8V3k07r7yKIAWlWg0jKr3YeyVXwijjjdDyGRy1NbGSaezFItKKCQkkzEyGbdXIoODaQqFAuIHvdLt8SmOLiSTcY4eTY08LhaVYlGZNq1yS34b48J4VxaPazeyc8Xd33yQ7q4+wNsYBqC7q4+7v/kgd9zubsIWIBoNU1OTYGjIG7sOhULU1MSJRt1eEQBMnz6F2Q2jJotzeQp5t2fHOiq+edtN6G8cd0VOsmjtZMfN5FRti0ODcKr9CG71b78gIncd/1WeLp6+5557iUw2D+qVokaVTDbPc8+95LztefMaGEoNjymBPZQaZt680yulfbouXLmI4XSGrD8xns3lGU5nnFdcDYVKn7riD0fJccfdSQ2euLDdyY6byacaF4cG4VRXBKUJ4u2UTvUmgfRwDhEZuRoQEaSopIfdDpMATJtWS01Ngny+QKFQJBwKUVMTZdq0Wqft3rTxGm6/o5fevhTpdJZoLEJj4yxu2niN03ZraxMUi8UxQ2EilCVDq6jeEFhpshigUCxSLMfliCkLq61UHqfaj2CLf3cn8Clg0ajXKHDSMtRBqknGGBhMUyiGCIXEHzsulmW7yOHhLGvWLKO1tXNkVfOCBbOdl1xYtWoxt3783WW/hF5z8TJ+8cwuMpkcxWKRcDhMLBZlzcXLTv3is1STjDGYGqaoRfylLV6acI37yqemPKy2UnmMN2vo34Bv4G1gv97/2uCqU2dr9eplJBMxRLxKmCKQTMRYvdr9h1NTY71XdG2UdDpTsbVgLv+di8jn8yNXA6pKPp/n8t+5yHnbl156gf/vHPLbD5FMxLj00guct23Kw2orlcd4U2m6VPV+pz2ZQDe9/xo6O/vo7R0km8sRi0aZMaOOm97vdpgEvEJoD/zoaS+Vsaj096fo7OzlLde4rbzd3LyPT//tJtrbe8lmc8RiUZ597j/59P/e6PSq4IknXyQSiVAo5EYmxyORCE88+aLzej9B/jvbBGZ5WG2l8hjvFcFfi8hXJ8s6glWrFvMHv3858XiUbDZPPB7lD37/8rL8j/rEky9SVCWXK5DN5cnlChRVeeLJF522e+cXfkBrayfFQpF4PEaxUKS1tZM7v+A2w/e5F14mlysQi0WprU0Si0XJ5QplWUfgDYe9izesW8n5y5p4w7qV3Prxd1XsJkDVyGorlcd4rwg+gLeBfRQo+sfO2XUEzc37eOjh7SxdMpeLX7uEwcG093hpo/M/oOdeeJlMJo9IKZ0SMpm88w/GHTtbEIR8oUg2lycUCiEIO3a2OG03nc4SEgj7WULhkFAUyrJncVBsArO8rLaSe+MNBJeo6nKnPZlAW7Zuo1AojilD3dAwrSz/ow4OpikWC4RDIS+RUr2JatcLrHK5vJelFA4dG6svuM+gqUnG6R9IkR/OjskamlrrfmK+dGZeV5scc2bu+ozRJjBNpRnv0NDPRWTSlHTcvecA+/YdJpvNkUzGyWZz7Nt3mN17Djhv27sKEPxLAn91k4xcHbji1RRS/DxOf0WXOq81tHjxeX5to8KY28WLz3PaLozdd+KpbTt4eW8b+XzB+Z7FNoFpKs14PyVeDzwvIvuADMfKUL/WWc/OQio17Jd28N5eNBohl8uTSrlfaFRXm6C/MIRyLK8+Eg5R5ziv/rzzZrJv32G/3s6x1bXnnTfTabuo98cQDoeOtatallUnu/a00tHRSzQaGbPvRNpxqm61TmDaBHnlGu8VwbXA+Xi1hjZwjqeP1tTGURizylb9466tXr1spK5RqfBbLB5xnrq6eOF51NQkCIfDiEA47JW6WLzQ7Zn5wUNdzJg5lZqaOPFYhJqaODNmTuXgoS6n7cLYfSdExLv1j7tUjROY1TpBvnnzU7zzhs/we1d+jHfe8Bk2b34q6C45Md5aQ25nHCfYBcsXkM8XOHSwm+FMlkQ8RtO8ei5YvsB525dfdhFPP72baCRMUUOERBBCXH6Z47x6gUg4DDEoFL3SzJFwuCybyMdjkTGL9bLZHIVC8VVeMTFqaxO0d/TQc2RgJHU1Ho8wZ47jqyCqbwIz6O1fg7B581N8/q7vk4h7lQEGB4f4/F3fB6i4rVDHe0Uwqay8YCGtrR3k8l72TC6fp7W1oyybmu/c1cLKlQuYO3cW06fXMXfuLFauXMDOXW5jaVfXUUIhvOkJvNtQyDvu0oUXLmJ4ODtmE6Dh4SwXXrjIabvgXQkMD+fGlN4eHs6N7EznUnPzPj5727f404/exWdv+1bFnxkHuf1rUO6592ES8Sg1NQlCItTUJEjEo9xz78NBd23CVWQgeOLJF4mEI36tIfUWOYUjznP5wcsomdfUwJrVy/idN6xizeplzGtqcJ5RkkoNE4tFmd0wnfPOm8nshunEYlHn8yI3vf8aGhvrAWFoaBgQGhvry7Ko63DHEe+Ka9RXSITDHUectluNwyRBDcMFqedIP4nk2OHkRDJOz5H+gHrkjrNTJxFJAI8Dcb+d76nqX4vIYuA+YBbwLHCjqk7o7N6Onfupm5Icc2aYzeWdbxcJXkZJX9/gSG45lCejpKY2zsDgENlcnmgkTC5fKMu8SGlRVxCTiEePpkolho4RxuxR4EI1DpPU1iYYGEiRy+WJ+HuAl2v716DMmjmVwcEhamqOvcfhdIZZM8/ZXXrPmMtr6AzwZlUdFJEo8KSI/Aj4C+DzqnqfiPw/4IPAFye89ePz58tUkTKojJIg50WCGi8vlbourZso7YcTIBnFAAAd7klEQVTgugR2UNlKQVqxfD7JZIyurqMjwa+pqZ6FC+YE3TVnbnzvlSNzAolknOF0huFMjg9/6MqAezbxnA0NqWfQfxj1vxR4M/A9//gm4O0T3XaQ49ZBZZSsvGAhbW09JJMx5s6dRTIZo62tpyzzIkGpq0uSyxVGtqYsFArkcgXqHFeZDXqYJIj5iQ3r1xEOh1i2tJHL1l3IsqWNhMOhsmz/GpTrr7+MP//oH1BXV8PRoynq6mr484/+QcVNFIPbKwJEJIw3/LMM+FdgL9Cnqnn/KQeBppO89mbgZoAFC07vrPam91/D7X4xsqGhYWLRaNnGrSGYM+Sdu1poapzFwYPd9PYNkojHmDevnp27Wpz/4QaVX55MxEnEI2Rz3iK2UEhIxCMkE26Hw4IcJglqNXXpBGf0v/ON77uqLP/OQa5fuP76yyryg/94TieLVbWgqquBecClePWKxvvaL6vqWlVd29Bwert7BVWMLEi79rTS3dPPlKk1NM6dxZSpNXT39DvP6gh64rRQVOKxKLW1CeKxKAXHe0ODN0yyePFcYrEo6XSGWCzK4sVzWbF8vvO2R9c5EhGmTKmhrjbpfDV1UIL++6oWZdnRXVX7RORnwDpguohE/KuCecAhF20GmecdxBnM6OEK/Nt8Lu98uCLoAmxT6mrI5Qvk897ZeTLpftHghvXruP2O75DN+gsWs3mGhjJlGSYJqs5RUFciQf99VQuXWUMNQM4PAkngKuBzwM+Ad+JlDm0ENrvqQxCC2hcgqOGKIAuwld5zTU185D1ns7kyZbKov2ajtHCjPMkIQWWlBVXI0Qr8lYfLK4K5wCZ/niAEfFdVt4rITuA+Efk74Dngaw77UHZ3fuEH7N/fTmlJ7/Bwlv3727nzCz/gK1/6mLN2g8rqaGqsp6W1Y0y7DQ3TypJNEtR73rJ1G3PPm8Vrzj82FDQwMFSWs9SgstJ27zlAe3sPsVh0TCFH1+XGgwp8JdVSX8ll1tCvVXWNqr5WVVep6t/4x19R1UtVdZmq/qGqOhm7CGrl569ffIViUfFrkCJAsaj8+sVXnLYbVFbHygsWsnvPAVKpNIlEjFQqze49B8qSrbRh/Tq6e47S0dnLkSP9dHT20t1z1Pl7PtTW/RuZSeU6Sw0qKy2VGiafL3K0f4iOjiMc7R8iny86X7C4Yf06BlNpBgaGUFUGBoYYTKXLMgxXTfMTZZkjKLegxjMB8vkCodDYstOhkJDPF5y2u2rVYlZesJDvfu8xhoYy1NTEueGdb3T+fnfuamH5a+bR3d3vnZXXJlm4cE5ZspX27m2jvf0IxaISDocpFou0tx9h7942p+876LPUoOa/BgaHiITDhMIh8rk8w8MZZsyY4rTNILOVgpyf2Lz5Ke6592F6jvQza+ZUbnzvlU7/f6rIQBDkys9oNEI2m0PCSmnZa7GoxGJRp+1u3vwU933nUQqFAuFwiFw2z33feZR5TQ1O/4BKJTXmz5s9ckxVy3J2fM+9D1NXmxyz8nNoaJh77n3Y6Xuu1jLU8ViETDZPIZsjHA4Rd7zXRUlQgS+o+Ykgit1VZK2hIAtkXbRqMaGQ+HvDKKreFcFFjv+Qv/L1B8hkc4RCIaJRr85SJpvjK19/wGm7QW7SElQtmFWrFnP1lWvZ+8phfvro8+x95TBXX7m2IseOSxQlky0QDodJJuOEw2Ey2QJapknyIAT1tx1EsbuKDARBrvy85SPvYNGi80gkYt4Cp0SMRYvO45aPvMNpu11dfUQiIW9YCi/4RCIhurr6nLYb5BjurJlTGU6P/TctRy2Y5uZ9fP+HT5DJ5IjFImQyOb7/wycqcuy4RBDq6hJEoxGKxSLRaIS6uoSXNVWhgvrb7jnSTy6X5+ChblpaOzh4qJtcLu/0BKciA0FtbQJVHVNiolwrP1etWsyn//dG3vaW13HpJSt421te5zx1FLyNaE5UXikcDjttN8hNWm5875WkUhna24/QdriH9vYjpFIZbnyv21owd3/zQdraugH1h6WUtrZu7v7mg07bDVJtbYJoJMy0qTXMmTOTaVNriEbCZVtNHUTiR1BXfpFwmCO9gxSLRUSEYrHIkd5Bb38RV206+8kBCrpAVhBjmhdesJAXfr0X8LaNLBSKFPIFLn7tUudtBzWGu3RpI7PnTKOjo49CLk8oEmb2nGksXdrotN0dO/aTSMTGbIWqquzYsd9pu0FasXw++Xyeg6OKGs6bV+98NXVz8z5uv+M+evtS5LJ5Wlo72bW7hVs//m7nf3PNzft46OHtLF0yl4tfu4TBwbT3eGmj07YldGw/qZGCiv5xVyryiqAaC2Td8tF3MH9BgzdRnMsTDoeYv6CBWz7qdkgqSFu2bmPZkiauvfoSNrxtHddefQnLljSVpdxCNpenu6ef9o4jdPf0k83lT/2iSWzlBQs5dFxRw0NlKGp496YHaWvrQVVJJmOoKm1tPdy9yf3VV1DlPHK5AlOmJlFV8vkCqsqUqV6BRVcq8oogyJSzoKxatZjP/NVNVbH4pSSorI558xp44YW9RCLHUin70hkuvtj91RcEs8hp564WViyff+wquzbJokXnOU8T3rFzP4lkfEzplNJx14L6+6qtSdDR2UsoFELEK68+NJT9jb5MpIoMBFB9e8pC9b3noFY1T5taSzIZo1BUCoUi4XCIZDTMtKm1TtuF4IZKSkFnXtOxApDlShMOam+RoNaLTJniBZvSeqRisUixqEyZ4q68esUGgmpZGl7NVl6wkAd/sp1EPEoiGSeVStNzpJ+3XHOp03aHM1mWLm1k797D5HJ5YtEIS5fOZTjjfmOa0lBJIhknmYyRyxdGhkru+Mc/dtZuUEH3wgsX8fzzLyMiI/WkhoezrF69zGm7UCoueB87d7eSy+aJxiLMmF7rfL1IX98gtbVx0unsSHn12to4fX2Dp37xGarIOYJqWhpezUp7MKTTWQ4f7iGdztLUOIudu1qctpuIx2ht7RxT8ru1tZNEPOa0XRg7VFJKjU4k486HSoIqJRLkntgAQ+kMvb0DHOntp7d3gKG0+xT0bDZHNpsnFPLSwUOhENmsV1DRlYoMBNVWs71a7d5zgO7uo0ydWkNjYz1Tp9bQ3X2U3XsOuG1YvJTkvr5B2jt66esbJJfLU7aU+gCGSkqlRGprkwwPZ6mtTbL8NfOcB90g9xa5+5sP0tc7yIzpdTQ1NTBjeh19vYPO04TDkTD5nDdJHAqFvEnjXIFwxNJHT4uVrq0OqdQwIjImjTOXyzsvhNbZ1Uc4EqaQzVMqPx2OhOl0vHgPvKGS7dv/k3w+RaFYJBwKEYmEWbv2NU7bDbKUSFB27NhPoVCko7PPL90SprYm7jxNuFgoEo6EgGN7cYcjIYqForM2K/KKIMiyB6Z8amrjKF4qp6p6t/5xl4ZSGWLRCLMbpjP3vFnMbphOLBopy8r1yy+7iEwmy3AmSzaTYziTJZPJcvllFzltN6j/p4Ic5h1KZzjan6JYLBIKhSgWixztTzkfHopGI8yYPoVkMk4s5m22NGP6lJETHhcqMhAEWfbAlM8FyxewaOEc4rEo6eEs8ViURQvncMHy09vj+nQFuXL9iSdfJBaPkkjEiMUj/m2UJ5580Wm7Qf0/FeQwb6lisLfvkI6M/LmuJHzhykUUi0WmTq1hzuwZTJ1aQ7FY5MKVi5y1WZGBIMiyB6Z8NqxfRyQSHrNwMBIJO/9wCnLP4h079zNlSs2Yq5EpU2qcTxYH9f9UkHs/iAihkFAoFskXihSKxd8oMe/CTRuvobFxFiJCOp1FRGhsnMVNG91NkFfkHAFUX059NQpq4WCQexYDgeXVB/H/VJA74E2bVksmkyUcFopFJRSSkeMueRPk7y5r+nvFBgJTHYIK+KmhNL19A+RyeaLRCLG42+J+Jd5k8R7yhRTFQpFQOEQkHGLt2uVlab/cglorAtA4d5ZX1Tcc9pID8gXyhQKNc2c5b7vcf9cVOTRkjEt3b3qQo30ppk+vo3HuLKZPr+NoX6os9W8uv+wihjM5Uqlh0sNZUqlhhjM555PFQQkqbRUgHo+yZPFc8gVvS858ociSxXOJx91uMhUEuyIw5jQFWf9m6wO/IJ/3igqW5PN5tj7wC+dbgwYhyLTVRCJGe3svsWiEcDhEOBSivb2XefMaTv3iScYCgTFnIqBx+h07W4hGImNSCXO5PDt2uj9DDkKQ+0MfPZpiMOWlzJby+cl4xyuNBQIzqQVRUyrI+jeFQoHIcStMQ+GQ85TGoAS5P/S+fYcBRrKESnsDlI5XEgsEZtJqbt7HX//N3XR09JHL5Xk+upftv9rDZ/7qJqfB4Kb3X8PtnX309g4yNDRMLBotW/2b+vppdHZ6K5hLC50AZs+e7rztIARZUn4onSUajRAdFXhz+QJDaffFBcvNAoGZtO686wccaO0iHAkTjUYoFIocaO3izrt+wFe+/DFn7Zbq3wRR3fat117K177xI0S8EgTFoqJa5K3Xus+iCUpQmWHJRJSBwTT5fH4kfVREmFLnrhx0UJxlDYnIfBH5mYjsFJEdInKLf3ymiPxERF7yb93ttmAq2o5dLV4QiIQJiRCNeGl+O8qQUbJ3bxvP/uolduxs4dlfvcTevW3O2wTo7x9ixYoF1NUmiERC1NUmWLFiAf39Q2Vpv5osWdJIsVD0p38EVa8O0JIlbrdCDYLLK4I88DFV/ZWITAGeFZGfADcBj6jqbSLyCeATwF867IepUIVCARFhOJMbGSYJhbyNPFzavPkpPn/X90nEo0ybVsvg4BCfv+v7AM4zdw61dXPB8gWsXHGs/HOlF38LyrRptdTV1ZDPF8YU+HO9oCwIzgKBqh4GDvv3B0RkF9AEXA+8yX/aJuBRLBCYMzBtWi1dXUcprfgvFr2zt4aGaU7bvefeh0nEo9TUeLWFSrf33Puw80AQZBZNtW32NDycZc2aZbS2do6sal6wYDbDw+7nCMr9uy7LgjIRWQSsAZ4G5vhBAqAdOOFacRG5WUS2i8j2rq6ucnTTTDKzZkwFvMzN0tfo4670HOknkRxb4TSRjNNzpN9pu+Bl0Rxu7+HpX+7myaeaefqXuznc3uO8vEU1bvbU1FhPLBphzepl/M4bVrFm9TJi0UhFVlx1HghEpA74PvBnqjrm/xRVVUoF3Y+jql9W1bWqurahofIWcJiz19nVS+i4v+BQyDvu0qyZUxk+rhTxcDrDrJluA1BJELtmVeNmT17QPcLTz+zmiade5OlndnO4/UhFVlx1GghEJIoXBO5V1R/4hztEZK7//blAp8s+mMqVGsqgCpFIeORL1Tvu0o3vvZJUKkN7+xHaDvfQ3n6EVCrDje+90mm7ENyuWUFWAQ2WgoB4/+Ek560TKojftcusIQG+BuxS1X8a9a37gY3+/Y3AZld9MMaFpUsbmT1nGqGwt41gKBxi9pxpLF3qPptkx479JBIxotEIgreJSSIRc75rVpCbPW3e/BTvvOEz/N6VH+OdN3yGzZufct4meGfmyVIpEfFKiSSTcedXQUH8rl1eEVwG3Ai8WUSe97/eCtwGXCUiLwFX+o+NOW21tYmRVbbqTxBEImHnG8Rs2bqNZUuauPbqS9jwtnVce/UlLFvSVL5hkuPr4Tuujw/BbUxTytAaHBwak6FVjmCwa08r+1s6yGRzJBMxMtkc+1s62LWn1Wm7QfyunQUCVX1SVUVVX6uqq/2vB1S1R1WvUNXzVfVKVT3iqg+msq25eBnJZJx4PEo0GiYej5JMxllzsdtSD0EOk1y4chHD6cyY7TmH0xmnu1dBcBvTjM7QColQU5MgEY9yz70PO20XvC1JBe9KQES8W/+4S6tWLebqK9ey95XD/PTR59n7ymGuvnKt7UdgzInctPEaOrt66e1LkcvmicYizJhe63QnJwg2hfOmjddw+x3ee06ns0RjEee7V5UEscK350g/iUSU7p5+8vk8kUiE2tryZGjV1iYYGEiRy+VHakqVY0vS5uZ9PPTwdpYumcvFr13C4GDae7y00dnv3/YjMJNWaSenN7x+JcuWNfKG16/k1o+/uyw7lAW1J3ZQ7zkotbUJenoGKBaKhMNhioUiPT0DZdkfOqgtSYPIGrIrAjOpBXGWGmQhtFL7lfrBf7y5c2bS2dHrreyVEIViEVSZO2em87a9LUnvI5vLg0I2l2doaNh5wD/U1s2c2WMr77geerRAYMwZCPLDuJpW+MYTUVatWsTevYcZzmRJxGMsXTGfeKJcu4QJKKj3H/wcUqeCGHq0QGDMJFJadVpXmxyz6rQcE7dBBKDSh+KSxcdScwcGhpg+vc5pu+AN0cw9byavOX/emLa3bN3m9H0HsQeDzREYM4kEtcI3qBITQc7HBJUdFkSGll0RGDOJBDF+DGMDEDBy6/rsOMj5mKbGenbs3M+BA10jw1Lz5zc4T9WF8g89WiAwZhJpaqxnf0s73d39IxUx6+unsmjheU7bDSoAQXDzMVOn1rBzV4u/CRAM5ofZuauFda9fWfa+uGZDQ8ZMIisvWMie/zxIKpUmkYiRSqXZ858HWXnBwlO/+CwEWWKiuXkfn73tW/zpR+/is7d9q2wVTx/48TOEwyG/sKESCkE4HOKBHz9TlvbLyQKBMZPIzl0trFg+n9raJMPDWWprk6xYPp+djndlC2qsPsjy193dR4nFoiQTcWprEiQTcWKxKN3dR523XW42NGTMGQgqhfNQWzfJ4/ZCSCbjZZvALPdY/Zat28jnC7y8t23MUJjruQlgZAFbeFSt89LCNtfK/fdlgcCY09TcvI/b7/gOvb2DZHM5Wlo62bW7lVs//i7nH06JRIznn3+ZZDJOMhknm83x4ouvsHq12/pKEMxY/a49rXR09BKNRsYUfkuXYZewC1cu5LnnXyKXL4w6qqxZfb7Tdr2/r/tGSqe0tHaya3eL0xXkNjRkzGm6+5sP0tbWDai/TaXS1tbtfE8AYKQcvgKoHquO775MfiCCKvwGsP6trycW8xaulfbBjsWirH/r6522e/emB2lt7aK3d4Cj/YP09g7Q2trF3Zvc/X3ZFYGZ1IIYohm9JwB4ewKoqvM9AQCGM1kuWrWY1gNdI0Ml5y9tZDjj/gw5CEEVfgNvPuaiVYt/I0Nr564Wp3tTP//CXtLDWaKRMJFIBC0q6eEsz7+w11mbFgjMpBXkKtsg9gSAYytt14waChoYGGL2camdlWLF8vkkkzG6uo6OfBg3NdWzcMEJtzqfUIfaupnX1MD8ebNHjqmq8/mYoXSGkEAo5P1NSUgICU63JLWhITNpBbXKNqg9ASDYlbZB2LB+HeFwiGVLG7ls3YUsW9pIOBwqy/sNKmU2mYxRVCgUvaG/QlEpqnfcFQsEZtIKqgTATRuvobFxFiJCOp1FRMq6J0AQG8QEJcj3G1TQXXPxMhKJKCKQz+cRgUQi6nTDJRsaMpNWUBvElPYECKoCaDWVoYbg3m9pp7B77n2YniP9zJo5lRvfe6XzvgSx4ZKU9no9l61du1a3b98edDfMOWb0HMHoKo2VfIZsyifIv6+JSoIQkWdVde2pnmdXBGbSCnqDGFPZgiq0B1Z0zpjTUm3DJKZ8giy0V242WWyMMScQZKG9crNAYIwxJ1BNqbrOhoZE5OvAeqBTVVf5x2YC3wEWAfuBG1S111UfjKlE1bRncZCqaQ7KWdaQiPwuMAh8c1QguB04oqq3icgngBmq+pen+lmWNWSMpxozpSzwnbnxZg05GxpS1ceBI8cdvh7Y5N/fBLzdVfvGVKKgVlMHJcj9CKpJuecI5qjqYf9+O+C+YIgxFSSo1dRBqbbAF5TA0kdVVUXkpONSInIzcDPAggULytYvY85lQa2mDkrQKZzVMixV7iuCDhGZC+Dfdp7siar6ZVVdq6prGxoaytZBY85l1ZTJAsHvlVwtw1LlDgT3Axv9+xuBzWVu35hJrVT/Zu8rh/npo8+z95XDXH3l2oo8S4VgA181DUs5CwQi8m1gG7BcRA6KyAeB24CrROQl4Er/sTFmnJqb9/HQw9tZumQub37TapYumctDD2+vyLNUCLb6aDXNxzibI1DV95zkW1e4atOYShdk/ZugBFVGpJrmY2xlsTGTSDWdpQatmuZjLBAYM4lUU/2boFXTJkBWfdSYSWTD+nV88Uv3A4xZWXzj+64KuGeVqVqq29oVgTGTSDWdpZrysSsCYyaZajlLNeVjVwTGGFPlLBAYY0yVs0BgjDFVzgKBMcZUOQsExhhT5ZztUDaRRKQLaDnDl9cD1bbs0t5zdai291xt7xfO/j0vVNVTlm+eFIHgbIjI9vFs1VZJ7D1Xh2p7z9X2fqF879mGhowxpspZIDDGmCpXDYHgy0F3IAD2nqtDtb3nanu/UKb3XPFzBMYYY15dNVwRGGOMeRUWCIwxpspVdCAQkWtFZI+IvCwinwi6P66JyHwR+ZmI7BSRHSJyS9B9KgcRCYvIcyKyNei+lIOITBeR74nIbhHZJSKVt2XWcUTkz/2/6WYR+baIJILu00QTka+LSKeINI86NlNEfiIiL/m3M1y0XbGBQETCwL8CbwFWAu8RkZXB9sq5PPAxVV0JvB74kyp4zwC3ALuC7kQZ3Qn8WFVXABdT4e9dRJqAjwJrVXUVEAbeHWyvnLgbuPa4Y58AHlHV84FH/McTrmIDAXAp8LKqvqKqWeA+4PqA++SUqh5W1V/59wfwPiCagu2VWyIyD3gb8NWg+1IOIjIN+F3gawCqmlXVvmB7VRYRICkiEaAGaAu4PxNOVR8Hjhx3+Hpgk39/E/B2F21XciBoAg6MenyQCv9QHE1EFgFrgKeD7Ylz/wzcChSD7kiZLAa6gG/4w2FfFZHaoDvlkqoeAu4AWoHDwFFVfSjYXpXNHFU97N9vB+a4aKSSA0HVEpE64PvAn6lqf9D9cUVE1gOdqvps0H0powjwW8AXVXUNkMLRcMG5wh8Xvx4vCDYCtSLyvmB7VX7q5fo7yfev5EBwCJg/6vE8/1hFE5EoXhC4V1V/EHR/HLsMuE5E9uMN/b1ZRP4t2C45dxA4qKqlK73v4QWGSnYlsE9Vu1Q1B/wAeEPAfSqXDhGZC+DfdrpopJIDwS+B80VksYjE8CaX7g+4T06JiOCNHe9S1X8Kuj+uqeonVXWeqi7C+/f9qapW9JmiqrYDB0RkuX/oCmBngF0qh1bg9SJS4/+NX0GFT5CPcj+w0b+/EdjsopGK3bxeVfMi8qfAg3hZBl9X1R0Bd8u1y4AbgRdF5Hn/2KdU9YEA+2Qm3keAe/0TnFeADwTcH6dU9WkR+R7wK7zMuOeowHITIvJt4E1AvYgcBP4auA34roh8EK8U/w1O2rYSE8YYU90qeWjIGGPMOFggMMaYKmeBwBhjqpwFAmOMqXIWCIwxpspZIDDmVYjIfhGpP8VzBsvVH2NcsEBgzEn4FWyDarti1/iYc48FAlNxROR/iMhH/fufF5Gf+vffLCL3+vffIyIv+vXtPzfqtYMi8n9E5AVg3ajjSRH5kYj80Una/HsReUFEfiEic/xji0TkpyLyaxF5REQW+MfvFpF3jm7Tv32TiDwhIvcDO0WkVkT+w/+5zSLyron+XRkDFghMZXoCuNy/vxao82swXQ48LiKNwOeANwOrgUtEpFTetxZ4WlUvVtUn/WN1wBbg26r6lRO0Vwv8QlUvBh4HSsHiC8AmVX0tcC9w1zj6/lvALar6Grza9G1+X1YBPx7n+zfmtFggMJXoWeC3RWQqkAG24QWEy/GCxCXAo34Rszzeh/Tv+q8t4BXtG20z8A1V/eZJ2ssCpd3RngUW+ffXAd/y798D/M44+v6Mqu7z778IXCUinxORy1X16Dheb8xps0BgKo5foXIfcBPwc7wP/98DlnHqYmXDqlo47thTwLV+wbMTyemxWi0FTl3DK4///56IhIDYqO+lRr2P/8S7QngR+DsR+atT/FxjzogFAlOpngA+jjdU8wTwx8Bz/gf2M8AbRaTenxB+D/DYq/ysvwJ68bY+PR0/59iWiu/1+wGwH/ht//51QPREL/aHsIZU9d+Af6Tyy02bgFggMJXqCWAusE1VO4Bh/xj+jk+fAH4GvAA8q6qnKu97C95WibefRh8+AnxARH6NVxX2Fv/4V/ACUWlCOnWS118EPONXkv1r4O9Oo21jxs2qjxpjTJWzKwJjjKlyFgiMMabKWSAwxpgqZ4HAGGOqnAUCY4ypchYIjDGmylkgMMaYKvf/A3cjTuBMNBEoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter([x_i[2] for x_i in x], y, color='#30336699')\n",
    "plt.title('Time on site vs. work hours')\n",
    "plt.xlabel('work hours')\n",
    "plt.ylabel('minutes')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32.33459183673472\n",
      "26.203428571428574\n"
     ]
    }
   ],
   "source": [
    "y_phd = [y_i for x_i, y_i in zip(x,y) if x_i[3]]\n",
    "y_not_phd = [y_i for x_i, y_i in zip(x,y) if not x_i[3]]\n",
    "print(mean(y_not_phd))\n",
    "print(mean(y_phd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFT1JREFUeJzt3X20XXV95/H3pyEWEHkIpJkAYqhQDDIaZi4+tNQOxDpYESijKKNOoJnJOA8oM62KZM3SmZqOMqwi1VntiqQlrmJ4UgSsy4FJg5rliE2EKhIdKoiABC7KkzxIgO/8cXbgmrkPJyH7nNy736+1su7e++y9f997OZzP2U+/X6oKSVJ3/cqwC5AkDZdBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQaGCSnJvkomHXsbMk+XmSXx92HTsiyRlJ1g+7Du0aDALtNM0H49Z/zyZ5Ysz8u6rqT6rqXw+7zp2lqvaqqtsBklyc5GPDrmmsJB9NsqX5+z+U5BtJXt/ntjckeTLJo0keSbIxyTlJfrXtujV4BoF2muaDca+q2gv4MfDWMcsuGXZ9HXVZ899jLrAe+EKS9Lntf6yqlwDzgT8E3gl8eTu21zRhEGhgmm+of91ML0hSSc5McleSB5O8N8kxSb7TfIP99Dbb/0GSTc26/yvJyyZp66Qk32v2c0OShWNe+1GSP2raeTjJZUl2n2A/hyX5arPeA0kuG/NaNa8vA94FfLD59n1t8/qBST6fZDTJHUneN0Ebr02yOcmsMct+P8l3munXJNnQfDO/L8mf9vP3HquqtgCrgX8E7D+mnfObv+cdSd48wbaPVdUNwEnA64G3bG/72rUZBBq21wKHA+8APgksB94IvBI4LcnvACQ5GTgXOJXet9uvA2vG22GS32heO7tZ98vAtUleNGa104ATgEOBVwFnTFDfHwPXAfsBBwOf2naFqloJXAKc1xz9vDXJrwDXAn8PHAQsBs5O8s/H2f5G4DHg+DGL/yXwuWb6QuDCqtobeDlw+QS1Tqg5pXMGcFdVPdAsfi3wA+AA4Dxg1WTf9qvqx8AG4Le3t33t2gwCDdsfV9WTVXUdvQ/DNVV1f1XdQ+/D/uhmvfcC/72qNlXV08CfAIsmOCp4B/A3VXV98034fGAP4DfHrPNnVfWTqvoZvQ/sRRPUtwV4GXBgU2e/F1iPAeZW1X+rqqeaawmfoXd6ZTxrgNMBkrwE+D2eD7otwGFJDqiqn1fVN/usAXph+hBwF/BPgd8f89qdVfWZqnqG3tHCfGDeFPv7CTBnO9rXNGAQaNjuGzP9xDjzezXTLwMubE71PAT8DAi9b9vbOhC4c+tMVT1L74Nw7Lqbx0w/PqadbX2waedbzammP5jyN3q+3gO31tvUfC4Tf9B+Dji1+eZ+KvDtqtr6OywFfgP4fpK/S3JinzUAXF5V+1bVr1XV8VW1ccxrz/0NqurxZnKiv8NWB9H722sG2W3YBUh9ugtY0edF558A/3jrTHO646XAPdvbaFVtBv5Ns59jgf+d5GtV9Q/brjpOvXdU1eF9tnNrkjuBN/PLp4WoqtuA05vTTacCVybZv6oe297f54VI8lJ6RxWfGGS7ap9HBJou/gL4cJJXAiTZJ8nbJ1j3cuAtSRYnmU3vjpdfAN/Y3kaTvD3Jwc3sg/Q+8J8dZ9X7gLHPFHwLeDTJh5LskWRWkqOSHDNJc58D3g+8AbhiTA3vTjK3ObJ5qFk8Xg2tSLJnc63manq/15cH1bYGwyDQtFBVV9H7JnppkkeAW+h9ex5v3R8A76Z3YfcB4K30bmV9ageaPga4McnPgWuA9299dmAbq4Ajm9NAX2zOu59I79rDHU0dFwH7TNLWGuB3gL8dc0EXehe1v9fUcCHwzqp6Ap57dqOti7efTvIovZD7JPB54IQmkDSDxIFpJKnbPCKQpI4zCCSp4wwCSeo4g0CSOm5aPEdwwAEH1IIFC4ZdhiRNKxs3bnygquZOtd60CIIFCxawYcOGYZchSdNK85DilDw1JEkd11oQJDkiyc1j/j2S5Owkc5Jcn+S25ud+bdUgSZpaa0FQVT+oqkVVtYhe/ySPA1cB5wBrmz5Y1jbzkqQhGdSpocXAD5veFE+m1+Utzc9TBlSDJGkcgwqCd/J83+rzqureZnozE3TLm2RZMyrThtHR0UHUKEmd1HoQNKNCncSY3hS3ql5HR+N2dlRVK6tqpKpG5s6d8u4n9WHNmjUcddRRzJo1i6OOOoo1a8Yd4EtSxwzi9tE30xtkY+uAI/clmV9V9yaZD9w/gBo6b82aNSxfvpxVq1Zx7LHHsn79epYuXQrA6aefPuTqJA3TIE4Nnc4vjy17DbCkmV5Cr49ztWzFihWsWrWK4447jtmzZ3PcccexatUqVqxYMezSJA1Zq91QJ3kx8GPg16vq4WbZ/vQGDjmE3nCCpzXjxk5oZGSkfKDshZk1axZPPvkks2fPfm7Zli1b2H333XnmmWeGWJmktiTZWFUjU63X6hFBVT1WVftvDYFm2U+ranFVHV5Vb5wqBLRzLFy4kPXrf3nc9fXr17Nw4cIhVSRpV+GTxR2xfPlyli5dyrp169iyZQvr1q1j6dKlLF++fNilSRqyadHXkF64rReEzzrrLDZt2sTChQtZsWKFF4olTY+hKr1GIEnbb5e4RiBJ2vUZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSx7UaBEn2TXJlku8n2ZTk9UnmJLk+yW3Nz/3arEGSNLm2jwguBL5SVa8AXg1sAs4B1lbV4cDaZl6SNCStBUGSfYA3AKsAquqpqnoIOBlY3ay2GjilrRokSVNr84jgUGAU+KskNyW5KMmLgXlVdW+zzmZg3ngbJ1mWZEOSDaOjoy2WKUnd1mYQ7Ab8E+DPq+po4DG2OQ1UVQXUeBtX1cqqGqmqkblz57ZYpiR1W5tBcDdwd1Xd2MxfSS8Y7ksyH6D5eX+LNUiSptBaEFTVZuCuJEc0ixYDtwLXAEuaZUuAq9uqQZI0td1a3v9ZwCVJXgTcDpxJL3wuT7IUuBM4reUaJEmTaDUIqupmYGSclxa32a4kqX8+WSxJHWcQSFLHGQSS1HFtXyzWkCTZoe16j3ZI6hKDYIaa7AM9iR/4kp7jqSFJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOa3XM4iQ/Ah4FngGerqqRJHOAy4AFwI+A06rqwTbrkCRNbBBHBMdV1aKqGmnmzwHWVtXhwNpmXpI0JMM4NXQysLqZXg2cMoQaJEmNtoOggOuSbEyyrFk2r6rubaY3A/PG2zDJsiQbkmwYHR1tuUxJ6q5WrxEAx1bVPUl+Dbg+yffHvlhVlaTG27CqVgIrAUZGRsZdR5L0wrV6RFBV9zQ/7weuAl4D3JdkPkDz8/42a5AkTa61IEjy4iQv2ToNvAm4BbgGWNKstgS4uq0aJElTa/PU0DzgqiRb2/lcVX0lyd8BlydZCtwJnNZiDZKkKbQWBFV1O/DqcZb/FFjcVruSpO3jk8WS1HEGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1XF9BkOS8JHsnmZ1kbZLRJO9uuzhJM1OS7f6n9vR7RPCmqnoEOBH4EXAY8IG2ipI0s1XVuP+mek3t6DcItg5g8xbgiqp6uKV6JEkD1u8IZV9K8n3gCeDfJZkLPNleWZKkQenriKCqzgF+Exipqi3A48DJbRYmSRqMfi8W7wn8e+DPm0UHAiNtFSVJGpx+rxH8FfAUvaMCgHuAj7VSkSRpoPoNgpdX1XnAFoCqehzwfi5JmgH6DYKnkuwBFECSlwO/aK0qSdLA9BsEHwW+Arw0ySXAWuBD/WyYZFaSm5J8qZk/NMmNSf4hyWVJXrQjhUuSdo5+7xq6DjgVOANYQ+/uoXV9tvF+YNOY+U8AF1TVYcCDwNK+q5Uk7XT93jW0tqp+WlV/U1VfqqoHkqztY7uD6T2EdlEzH+B44MpmldXAKTtWuiRpZ5j0gbIkuwN7Agck2Y/nLxDvDRzUx/4/CXwQeEkzvz/wUFU93czfPdF+kiwDlgEccsghfTQlSdoRUx0R/FtgI/AK4NvN9EbgauDTk22Y5ETg/qrauCOFVdXKqhqpqpG5c+fuyC4kSX2Y9Iigqi4ELkxyVlV9ajv3/VvASUl+D9id3lHEhcC+SXZrjgoOpvdMgiRpSPrta+jhJP9q24VV9dmJNqiqDwMfBkjyz4A/qqp3JbkCeBtwKbCE3tGFJGlI+g2CY8ZM7w4spneqaMIgmMSHgEuTfAy4CVi1A/uQJO0kfQVBVZ01dj7JvvS+0felqm4Abmimbwde03eFkqRW7ehQlY8Bh+7MQrRj5syZs0MjPW3P+nPmzBnybympTX0dESS5lqZ7CXrhcSRweVtFqX8PPvhg66M3OUygNLP1e43g/DHTTwN3VtXdLdQjSRqwfq8RfLXtQiRJw9FvFxOnJrktycNJHknyaJJH2i5OktS+fk8NnQe8tao2TbmmJGla6feuofsMAUmamfo9ItiQ5DLgi4wZkKaqvtBKVZKkgek3CPYGHgfeNGZZAQaBJE1z/d41dGbbhUiShmOq8Qg+WFXnJfkUzz9Q9pyqel9rlUmSBmKqI4KtF4g3ME4QSJKmv6nGI7i2mbwVOBdYMGabYsd6H5Uk7UL6vVj818AHgO8Cz7ZXjiRp0PoNgtGquqbVSiRJQ9FvEHwkyUXAWnyOQJJmlH6D4Ex6A9jP5vlTQz5HIEkzQN9DVVbVEa1WIkkain77GvpGkiNbrUTSjOLoedNHv0cErwNuTnIHvWsEAaqqXtVaZZKmNUfPmz76DYITWq1CkjQ0/fY1dGfbhUiShqPfawSSpBmq31ND2kXVR/aGj+7TfhuSZqzWgiDJ7sDXgF9t2rmyqj6S5FDgUmB/YCPwnqp6qq06Zrr810cGckGuPtpqE5KGqM1TQ78Ajq+qVwOLgBOSvA74BHBBVR0GPAgsbbEGSdIUWguC6vl5Mzu7+VfA8cCVzfLVwClt1SBJmlqrF4uTzEpyM3A/cD3wQ+Chqnq6WeVu4KAJtl2WZEOSDaOjo22WKUmd1moQVNUzVbUIOBh4Db3+ivrddmVVjVTVyNy5c1urUZK6biC3j1bVQ8A64PXAvkm2XqQ+GLhnEDVIksbXWhAkmZtk32Z6D+B36Q19uQ54W7PaEuDqtmqQJE2tzecI5gOrk8yiFziXV9WXktwKXJrkY8BNwKoWa5AkTaG1IKiq7wBHj7P8dnrXCyRJuwC7mJCkjjMIJKnjDAJJ6jg7nZPUCjtEnD4MAkmtsEPE6cNTQ5LUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nLePzgBJWt3/fvvt1+r+JQ2XQTDN7ch92klav79b0vThqSFJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjmstCJK8NMm6JLcm+V6S9zfL5yS5PsltzU97NJOkIWrziOBp4A+r6kjgdcB/SHIkcA6wtqoOB9Y285KkIWktCKrq3qr6djP9KLAJOAg4GVjdrLYaOKWtGiRJUxtIN9RJFgBHAzcC86rq3ualzcC8CbZZBiwDOOSQQ9ovUtJO51gZ00PrQZBkL+DzwNlV9cjYN0ZVVZJxO8avqpXASoCRkRE7z5emGcfKmD5avWsoyWx6IXBJVX2hWXxfkvnN6/OB+9usQZI0uTbvGgqwCthUVX865qVrgCXN9BLg6rZqkCRNrc1TQ78FvAf4bpKbm2XnAh8HLk+yFLgTOK3FGiRJU2gtCKpqPTDRlaLFbbUrSdo+PlksSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxxkEktRxAxm8XoM31aDhE73ueLFS9xgEM5Qf6JL65akhSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjqutSBI8pdJ7k9yy5hlc5Jcn+S25ud+bbUvSepPm0cEFwMnbLPsHGBtVR0OrG3mJUlD1FoQVNXXgJ9ts/hkYHUzvRo4pa32JUn9GfQ1gnlVdW8zvRmYN+D2JUnbGNrF4ur1gTBhPwhJliXZkGTD6OjoACuTpG4ZdBDcl2Q+QPPz/olWrKqVVTVSVSNz584dWIGS1DWDDoJrgCXN9BLg6gG3L0naRpu3j64B/g9wRJK7kywFPg78bpLbgDc285KkIWqtG+qqOn2Clxa31aYkafs5HoGkgZts4CQHTRo8g0DSwPmhvmuxryFJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgkqeMyHR7sSDIK3DnsOmaQA4AHhl2ENA7fmzvXy6pqyu6bp0UQaOdKsqGqRoZdh7Qt35vD4akhSeo4g0CSOs4g6KaVwy5AmoDvzSHwGoEkdZxHBJLUcQaBJHWcQdAhST6a5J4kNye5JclJzfKLk7xtnPUXJHkiyU1JNiX5VpIzBl64ZoQkz4x5712RZM/mPXbLBOtfnOSOJH+f5P8m+WySgwdddxcYBN1zQVUtAt4O/GWSqd4DP6yqo6tqIfBO4OwkZ7ZepWaiJ6pqUVUdBTwFvLePbT5QVa8GjgBuAv42yYvaLLKLDIJppPn2tCnJZ5J8L8l1SfZoXluU5JtJvpPkqiT7TbavqtoEPE3vSU6ANyT5RpLbxzs6aLa5HfjPwPt24q+lbvo6cFgzPWu89/RY1XMBsBl48yAL7QKDYPo5HPifVfVK4CHgXzTLPwt8qKpeBXwX+MhkO0nyWuBZYLRZNB84FjgR+Pgkm34beMUOV6/OS7IbvQ/z7zaLJnpPj8f3XwscvH76uaOqbm6mNwILkuwD7FtVX22WrwaumGD7/5Tk3cCjwDuqqpIAfLGqngVuTTJvkvbzwn8FddQeSba+d78OrAIOZJz39CT78P3XAoNg+vnFmOlngP/vMHoKF1TV+VPsd7L/2Y4GNm1nmxI01wjGLmi+hGzPe/poYO3OL63bPDU0A1TVw8CDSX67WfQe4KuTbLJDkiwAzgc+tbP3LU0mPe+jdwrzK8OuZ6bxiGDmWAL8RZI9gduBnXVnz8uT3ATsTu900p9V1cU7ad/SVP5Hkv8C7Al8Eziuqp4ack0zjl1MSFLHeWpIkjrOIJCkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4/4fkQgt8XNQCeAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.boxplot([y_not_phd, y_phd], labels=['no PhD', 'PhD'])\n",
    "plt.title('Time on site vs. PhD')\n",
    "plt.ylabel('minutes')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is kinda misleading and points out something interesting. The coefficient for PhD is about 0.92, implying that, all other things being equal, a user with a PhD is likely to spend an extra hour on the site compared to a similar user without a PhD. But, the mean minutes spend on the site for PhDs is lower that than of non-PhDs. Whaaa?\n",
    "\n",
    "The explaination is that all other things aren't equal. We saw above that PhDs are sadly friendless creatures and number of friends has a strong influence on minutes-on-site."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Goodness of Fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def total_sum_of_squares(y):\n",
    "    \"\"\"the total squared variation of y_i's from their mean\"\"\"\n",
    "    return sum(v ** 2 for v in de_mean(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multiple_r_squared(x, y, beta):\n",
    "    sum_of_squared_errors = sum(error(x_i, y_i, beta) ** 2\n",
    "    for x_i, y_i in zip(x, y))\n",
    "    return 1.0 - sum_of_squared_errors / total_sum_of_squares(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6800074955952597"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multiple_r_squared(x, y, beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrap\n",
    "We take a slight digression into bootstrapping in order to compute bootstrapped estimates of standard errors for our regression coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap_sample(data):\n",
    "    \"\"\"randomly samples len(data) elements with replacement\"\"\"\n",
    "    return [random.choice(data) for _ in data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap_statistic(data, stats_fn, num_samples):\n",
    "    \"\"\"evaluates stats_fn on num_samples bootstrap samples from data\"\"\"\n",
    "    return [stats_fn(bootstrap_sample(data)) for _ in range(num_samples)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "close_to_100 = [99.5 + random.random() for _ in range(101)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "far_from_100 = ([99.5 + random.random()] + [random.random() for _ in range(50)] +\n",
    "[200 + random.random() for _ in range(50)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100.06586584030461, 100.03937022319404, 0.3008288923397343)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "median(close_to_100), mean(close_to_100), standard_deviation(close_to_100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD8CAYAAACSCdTiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAADkdJREFUeJzt3X+MZfVZx/H3U7aIaCu/ppuRBQeFtCFGaJmQbarVgm0oaEGlBGLsErfuH7UJRE1djSYYa1yaaH+kxrIW4sZYCqLN0oJt1y1ETRrawS6wsNBdNpCyWdgBoS3/tC48/nG/wO1kd+fOuffOnX32/Upu7jnf8+v5zpn5zJnvufdOZCaSpJpeN+kCJEnjY8hLUmGGvCQVZshLUmGGvCQVZshLUmGGvCQVtmqQlSLiCeD7wEvAwcycjYhTgNuAGeAJ4KrMfH48ZUqSuljKlfy7MvP8zJxt8xuB7Zl5DrC9zUuSVpAY5B2v7Up+NjOf7Wt7DPiVzNwfEdPAvZn55iPt57TTTsuZmZnhKpakY8z999//bGZOddl2oOEaIIGvRkQCN2XmZmB1Zu5vy58GVi+2k5mZGebm5rrUKUnHrIh4suu2g4b8L2bmvoh4E7AtIh7tX5iZ2X4BHKq4DcAGgDPPPLNrnZKkDgYak8/Mfe35APAF4ELgmTZMQ3s+cJhtN2fmbGbOTk11+mtDktTRoiEfET8REW94ZRp4D7ATuBNY11ZbB2wdV5GSpG4GGa5ZDXwhIl5Z/3OZ+eWI+CZwe0SsB54ErhpfmZKkLhYN+czcC5x3iPbngIvHUZQkaTR8x6skFWbIS1JhhrwkFWbIS1Jhg74ZSiptZuNdEzv2E5sum9ixVZ9X8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUZ8pJUmCEvSYUNHPIRcVxEfCsivtTmz4qI+yJiT0TcFhHHj69MSVIXS7mSvw7Y1Td/I/DxzDwbeB5YP8rCJEnDGyjkI2INcBnw2TYfwEXAHW2VLcAV4yhQktTdoFfynwA+Arzc5k8FXsjMg23+KeD0Q20YERsiYi4i5ubn54cqVpK0NIuGfET8GnAgM+/vcoDM3JyZs5k5OzU11WUXkqSOVg2wzjuA90XEpcAJwBuBTwInRcSqdjW/Btg3vjIlSV0seiWfmX+SmWsycwa4GvhaZv42cA9wZVttHbB1bFVKkjoZ5nXyfwz8QUTsoTdGf/NoSpIkjcogwzWvysx7gXvb9F7gwtGXJEkaFd/xKkmFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVNiiIR8RJ0TENyLigYh4OCL+orWfFRH3RcSeiLgtIo4ff7mSpKUY5Er+B8BFmXkecD5wSUSsBW4EPp6ZZwPPA+vHV6YkqYtFQz57Xmyzr2+PBC4C7mjtW4ArxlKhJKmzgcbkI+K4iNgBHAC2AY8DL2TmwbbKU8Dph9l2Q0TMRcTc/Pz8KGqWJA1ooJDPzJcy83xgDXAh8JZBD5CZmzNzNjNnp6amOpYpSepiSa+uycwXgHuAtwMnRcSqtmgNsG/EtUmShjTIq2umIuKkNv3jwLuBXfTC/sq22jpg67iKlCR1s2rxVZgGtkTEcfR+KdyemV+KiEeAz0fER4FvATePsU5JUgeLhnxmPgi89RDte+mNz0uSVijf8SpJhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9JhRnyklSYIS9Jha2adAFaeWY23jWxYz+x6bKJHVuqaNEr+Yg4IyLuiYhHIuLhiLiutZ8SEdsiYnd7Pnn85UqSlmKQ4ZqDwB9m5rnAWuD3I+JcYCOwPTPPAba3eUnSCrJoyGfm/sz8nzb9fWAXcDpwObClrbYFuGJcRUqSulnSjdeImAHeCtwHrM7M/W3R08DqkVYmSRrawCEfET8J/CtwfWZ+r39ZZiaQh9luQ0TMRcTc/Pz8UMVKkpZmoJCPiNfTC/h/zsx/a83PRMR0Wz4NHDjUtpm5OTNnM3N2ampqFDVLkgY0yKtrArgZ2JWZf9u36E5gXZteB2wdfXmSpGEM8jr5dwC/AzwUETta258Cm4DbI2I98CRw1XhKlCR1tWjIZ+Z/A3GYxRePthxJ0ij5sQaSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFGfKSVJghL0mFLRryEXFLRByIiJ19badExLaI2N2eTx5vmZKkLga5kv9H4JIFbRuB7Zl5DrC9zUuSVphFQz4z/xP43wXNlwNb2vQW4IoR1yVJGoGuY/KrM3N/m34aWD2ieiRJIzT0jdfMTCAPtzwiNkTEXETMzc/PD3s4SdISdA35ZyJiGqA9Hzjcipm5OTNnM3N2amqq4+EkSV10Dfk7gXVteh2wdTTlSJJGaZCXUN4KfB14c0Q8FRHrgU3AuyNiN/CrbV6StMKsWmyFzLzmMIsuHnEtkqQR8x2vklSYIS9JhRnyklTYomPy0nKa2XjXpEtYdpPq8xObLpvIcbW8vJKXpMIMeUkqzJCXpMIMeUkqzJCXpMIMeUkqzJCXpMIMeUkqzJCXpMIMeUkqzJCXpMIMeUkqzJCXpMIMeUkqzJCXpMIMeUkqzJCXpML8z1CLmOR/KvI/90gallfyklSYIS9JhRnyklSYY/Ir2CTvB0gVTepnapL317ySl6TCDHlJKsyQl6TCjpoxecenpdHyPSDHBq/kJakwQ16SCjPkJamwo2ZMXlId3mNbPkNdyUfEJRHxWETsiYiNoypKkjQanUM+Io4D/g54L3AucE1EnDuqwiRJwxvmSv5CYE9m7s3MHwKfBy4fTVmSpFEYJuRPB77TN/9Ua5MkrRBjv/EaERuADW32xYh4bNzHHKHTgGcnXcQYVe8f1O+j/TsKxI1HXDxIH3+m67GHCfl9wBl982ta24/IzM3A5iGOMzERMZeZs5OuY1yq9w/q99H+Hf3G3cdhhmu+CZwTEWdFxPHA1cCdoylLkjQKna/kM/NgRHwY+ApwHHBLZj48ssokSUMbakw+M+8G7h5RLSvRUTnMtATV+wf1+2j/jn5j7WNk5jj3L0maID+7RpIKO2ZCPiKui4idEfFwRFzf2s6LiK9HxEMR8cWIeOOg27b2GyJiX0TsaI9Ll6s/7fi3RMSBiNjZ13ZKRGyLiN3t+eTWHhHxqfYRFA9GxNsOs88L2tdjT1s/jrTfQv2byLkcUx//KiK+ExEvLmj/sYi4rW1/X0TMjLNv7ZjL2b9rI2K+7xx+cLy9G33/IuLEiLgrIh5tebOpb1m385eZ5R/AzwM7gRPp3Yf4D+Bseq8Q+uW2zu8Cfznotm3ZDcAfTbBf7wTeBuzsa/sYsLFNbwRubNOXAv8OBLAWuO8w+/xGWx5t/fceab+F+jeRczmmPq4FpoEXF7R/CPhMm74auK1Y/64FPn00nz96OfOuNn088F9936Odzt+yfkNP6gG8H7i5b/7PgY8A3+W1+xJnAI8Mum2bnkgwLKhvZsE32GPAdJueBh5r0zcB1xxqvb62aeDRvvlrgJuOtN9C/ZvYuRxlHxfsd2EIfgV4e5teRe8NOFGof9eyzCE/zv61dT4J/N4w5+9YGa7ZCfxSRJwaESfS+416BvAwr33ezvv50Td3LbbtKz7c/vS6ZTmGMAawOjP3t+mngdVtepCPoTi9tR9qncPtd7mNq3+wcs7lMH08kle3z8yD9C5yTh2u1E7G1T+A32rn8I6IONTP83IYSf8i4iTg14HtC7dfyvk7JkI+M3cBNwJfBb4M7ABeojdE86GIuB94A/DDJWwL8PfAzwHnA/uBvxlrR5Yoe7/yR/7yqXHtd8J1rMhzuVK+1uMy4v59EZjJzF8AtgFbRrTfzrr2LyJWAbcCn8rMvcPUcEyEPEBm3pyZF2TmO4HngW9n5qOZ+Z7MvIDeF/TxQbdt7c9k5kuZ+TLwD/Q+mXPSnomIaYD2fKC1D/IxFPta+6HWOdx+l9tY+rfCzuUwfTySV7dvIfJTwHNDV7t0Y+lfZj6XmT9os58FLhhBrV2Mon+bgd2Z+Ym+tk7n75gJ+Yh4U3s+E/hN4HN9ba8D/gz4zKDbtvnpvtV+g97QzqTdCaxr0+uArX3tH2h3+NcC3+37kxKANv+9iFjbXnXygQXbH2q/y20s/Vth57JzH5ew3yuBr7UrzeU2lv4tOIfvA3aNotgOhupfRHyUXoBfv2BRt/O33DcpJvWgd5f6EeAB4OLWdh29q/JvA5t47SbsTwN3H2nb1v5PwEPAg+0EHPEmyhj6dCu9oYX/oze+t57eGN12YDe9VwKd0tYNev/k5fFW82zffnb0Tc/SC7jHgU/3fU0Oud9C/ZvIuRxTHz/W9vVye76htZ8A/Auwh96rjH62WP/+mt59tgeAe4C3HG39o3d1n/R+Qe1ojw8Oc/58x6skFXbMDNdI0rHIkJekwgx5SSrMkJekwgx5SSrMkJekwgx5SSrMkJekwv4f4FxNiWYZnoAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(bootstrap_statistic(close_to_100, median, 100))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100.37472916630605, 100.4997364160267, 99.91542590974795)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "median(far_from_100), mean(far_from_100), standard_deviation(far_from_100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAADTtJREFUeJzt3W2spPVZx/HvTxaqaatAOSEbFj1gsYY3AtkQTB9eQKsUKqA2DaTRNZIQkzaBVFO3kpia+AI0tmrS2KCQrgYLtQ+BtDUWkdqYWOoBludSFoQIWXa3D0gbTZX28sXcWw+bHWbO08xw9ftJTs7MPffZ+8p/Zr875z4zZ1NVSJJe+X5k3gNIkjaHQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1MS2WR7spJNOquXl5VkeUpJe8e65556vV9XSpP1mGvTl5WVWVlZmeUhJesVL8vQ0+3nKRZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpqY6TtFN2J59+fmctynrrt4LseVpLXyGbokNWHQJakJgy5JTbxizqFL0kZ1/1mcz9AlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUxNRBT3JMkvuSfHa4flqSu5PsS3JrkuO2bkxJ0iRreYZ+NfDoquvXAx+uqtcD3wKu3MzBJElrM1XQk+wALgb+arge4Hzgk8Mue4DLtmJASdJ0pn2G/qfA+4HvD9dfBzxfVS8O158BTtnk2SRJazAx6EneARysqnvWc4AkVyVZSbJy6NCh9fwRkqQpTPMM/Y3AJUmeAm5hdKrlz4Djk2wb9tkBPHu0L66qG6pqZ1XtXFpa2oSRJUlHMzHoVfWBqtpRVcvA5cA/VdW7gbuAdw677QJu27IpJUkTbeR16L8LvC/JPkbn1G/cnJEkSeuxbfIu/6+qvgh8cbj8JHDu5o8kSVoP3ykqSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSExODnuRHk3wlyf1JHk7yB8P205LcnWRfkluTHLf140qSxpnmGfp3gfOr6ueAs4ALk5wHXA98uKpeD3wLuHLrxpQkTTIx6DXyneHqscNHAecDnxy27wEu25IJJUlTmeocepJjkuwFDgJ3AE8Az1fVi8MuzwCnbM2IkqRpTBX0qvpeVZ0F7ADOBX522gMkuSrJSpKVQ4cOrXNMSdIka3qVS1U9D9wF/DxwfJJtw007gGfHfM0NVbWzqnYuLS1taFhJ0njTvMplKcnxw+UfA94GPMoo7O8cdtsF3LZVQ0qSJts2eRe2A3uSHMPoH4BPVNVnkzwC3JLkD4H7gBu3cE5J0gQTg15VDwBnH2X7k4zOp0uSFoDvFJWkJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJamJiUFPcmqSu5I8kuThJFcP209MckeSx4fPJ2z9uJKkcaZ5hv4i8NtVdSZwHvCeJGcCu4E7q+oM4M7huiRpTiYGvar2V9W9w+VvA48CpwCXAnuG3fYAl23VkJKkydZ0Dj3JMnA2cDdwclXtH256Djh5UyeTJK3J1EFP8hrgU8A1VfXC6tuqqoAa83VXJVlJsnLo0KENDStJGm+qoCc5llHMb66qTw+bDyTZPty+HTh4tK+tqhuqamdV7VxaWtqMmSVJRzHNq1wC3Ag8WlUfWnXT7cCu4fIu4LbNH0+SNK1tU+zzRuDXgAeT7B22/R5wHfCJJFcCTwPv2poRJUnTmBj0qvoXIGNuvmBzx5EkrZfvFJWkJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJamJiUFPclOSg0keWrXtxCR3JHl8+HzC1o4pSZpkmmfoHwMuPGLbbuDOqjoDuHO4Lkmao4lBr6ovAd88YvOlwJ7h8h7gsk2eS5K0Rus9h35yVe0fLj8HnLxJ80iS1mnDPxStqgJq3O1JrkqykmTl0KFDGz2cJGmM9Qb9QJLtAMPng+N2rKobqmpnVe1cWlpa5+EkSZOsN+i3A7uGy7uA2zZnHEnSek3zssWPA/8KvCHJM0muBK4D3pbkceCtw3VJ0hxtm7RDVV0x5qYLNnkWSdIG+E5RSWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNTHxjkTRLy7s/N5fjPnXdxXM5rrSZfIYuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqYkNBT3JhkseS7Euye7OGkiSt3bqDnuQY4CPA24EzgSuSnLlZg0mS1mYjz9DPBfZV1ZNV9T/ALcClmzOWJGmtNhL0U4D/WHX9mWGbJGkOtm31AZJcBVw1XP1OksfW8cecBHx986aaXq5/2ZvnNtcEzrVGuX5hZ3OutVnIuTbh8fVT0+y0kaA/C5y66vqOYdtLVNUNwA0bOA5JVqpq50b+jK3gXGuzqHPB4s7mXGvzwz7XRk65/BtwRpLTkhwHXA7cvjljSZLWat3P0KvqxSTvBf4BOAa4qaoe3rTJJElrsqFz6FX1eeDzmzTLy9nQKZst5Fxrs6hzweLO5lxr80M9V6pqFseRJG0x3/ovSU0sdNAX5VcLJDk1yV1JHknycJKrh+0fTPJskr3Dx0Vzmu+pJA8OM6wM205MckeSx4fPJ8x4pjesWpe9SV5Ics081izJTUkOJnlo1bajrk9G/nx4zD2Q5JwZz/XHSb46HPszSY4fti8n+e9V6/bRGc819n5L8oFhvR5L8otbNdfLzHbrqrmeSrJ32D7LNRvXiNk+zqpqIT8Y/aD1CeB04DjgfuDMOc2yHThnuPxa4GuMft3BB4HfWYC1ego46YhtfwTsHi7vBq6f8335HKPX0s58zYC3AOcAD01aH+Ai4O+BAOcBd894rl8Atg2Xr1811/Lq/eawXke934a/B/cDrwJOG/7OHjPL2Y64/U+A35/Dmo1rxEwfZ4v8DH1hfrVAVe2vqnuHy98GHmXx3xV7KbBnuLwHuGyOs1wAPFFVT8/j4FX1JeCbR2wetz6XAn9dI18Gjk+yfVZzVdUXqurF4eqXGb2/Y6bGrNc4lwK3VNV3q+rfgX2M/u7OfLYkAd4FfHyrjj/OyzRipo+zRQ76Qv5qgSTLwNnA3cOm9w7fMt0069MaqxTwhST3ZPTOXICTq2r/cPk54OT5jAaM3qOw+i/ZIqzZuPVZpMfdbzJ6FnfYaUnuS/LPSd48h3mOdr8t0nq9GThQVY+v2jbzNTuiETN9nC1y0BdOktcAnwKuqaoXgL8Afho4C9jP6Nu9eXhTVZ3D6DdfvifJW1bfWKPv8ebycqaM3nR2CfB3w6ZFWbMfmOf6jJPkWuBF4OZh037gJ6vqbOB9wN8m+fEZjrRw99tRXMFLnzjMfM2O0ogfmMXjbJGDPtWvFpiVJMcyuqNurqpPA1TVgar6XlV9H/hLtvBbzZdTVc8Onw8CnxnmOHD4W7jh88F5zMboH5l7q+rAMONCrBnj12fuj7skvwG8A3j3EAGGUxrfGC7fw+hc9c/MaqaXud/mvl4ASbYBvwLcenjbrNfsaI1gxo+zRQ76wvxqgeHc3I3Ao1X1oVXbV5/z+mXgoSO/dgazvTrJaw9fZvRDtYcYrdWuYbddwG2znm3wkmdNi7Bmg3Hrczvw68OrEM4D/nPVt8xbLsmFwPuBS6rqv1ZtX8ro/yAgyenAGcCTM5xr3P12O3B5klclOW2Y6yuzmmuVtwJfrapnDm+Y5ZqNawSzfpzN4ifA6/1g9JPgrzH6l/XaOc7xJkbfKj0A7B0+LgL+Bnhw2H47sH0Os53O6FUG9wMPH14n4HXAncDjwD8CJ85htlcD3wB+YtW2ma8Zo39Q9gP/y+hc5ZXj1ofRqw4+MjzmHgR2zniufYzOrR5+nH102PdXh/t3L3Av8Esznmvs/QZcO6zXY8DbZ31fDts/BvzWEfvOcs3GNWKmjzPfKSpJTSzyKRdJ0hoYdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJamJ/wPSnT3gPfbcHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(bootstrap_statistic(far_from_100, median, 100))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standard Errors of Regression Coefficients\n",
    "\n",
    "Now, we can bootstrap estimates for our betas, getting some idea of their distributions and estimating standard error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_sample_beta(sample):\n",
    "    \"\"\"sample is a list of pairs (x_i, y_i)\"\"\"\n",
    "    x_sample, y_sample = list(zip(*sample)) # magic unzipping trick\n",
    "    return estimate_beta(x_sample, y_sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(0)\n",
    "bootstrap_betas = bootstrap_statistic(list(zip(x, y)),\n",
    "                                      estimate_sample_beta, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.953551702104508,\n",
       " 0.06288763616183773,\n",
       " 0.11722269488203318,\n",
       " 0.8591786495949066]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bootstrap_standard_errors = [\n",
    "    standard_deviation([beta[i] for beta in bootstrap_betas])\n",
    "    for i in range(4)]\n",
    "bootstrap_standard_errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the book, we should get the following:\n",
    "```\n",
    "# [1.174,      # constant term, actual error = 1.19\n",
    "#  0.079,      # num_friends, actual error = 0.080\n",
    "#  0.131,      # unemployed, actual error = 0.127\n",
    "#  0.990]      # phd, actual error = 0.998\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_value(beta_hat_j, sigma_hat_j):\n",
    "    if beta_hat_j > 0:\n",
    "        # if the coefficient is positive, we need to compute twice the\n",
    "        # probability of seeing an even *larger* value\n",
    "        return 2 * (1 - normal_cdf(beta_hat_j / sigma_hat_j))\n",
    "    else:\n",
    "        # otherwise twice the probability of seeing a *smaller* value\n",
    "        return 2 * normal_cdf(beta_hat_j / sigma_hat_j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      beta     b.s.e.      p val\n",
      "     30.62       0.95       0.00\n",
      "      0.97       0.06       0.00\n",
      "     -1.87       0.12       0.00\n",
      "      0.92       0.86       0.29\n"
     ]
    }
   ],
   "source": [
    "print('{:>10} {:>10} {:>10}'.format('beta', 'b.s.e.', 'p val'))\n",
    "for i in range(4):\n",
    "    print('{:10.02f} {:10.02f} {:10.02f}'.format(beta[i],\n",
    "                              bootstrap_standard_errors[i],\n",
    "                              p_value(beta[i], bootstrap_standard_errors[i])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected p-values from the book.\n",
    "```\n",
    "# ~0   (constant term)\n",
    "# ~0   (num_friends)\n",
    "# ~0   (work_hours)\n",
    "# 0.36 (phd)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bootstrap distributions of betas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xt8lOWZ8PHflWSSSUJIgAAC4RAO4ZyAUI4GKViLtYtua3dtFUFXfdtPLd1t12pX19OufbXW3be1ulbXw7b1UA+tpavV1gOCimhAlDMEDBDOBBLIYZLJzP3+MQcmk5lkMpknM0+4vp+PnyTPPHmeO6NeuXLd93PdYoxBKaVU75KW7AEopZRKPA3uSinVC2lwV0qpXkiDu1JK9UIa3JVSqhfS4K6UUr1Qp8FdRJ4UkWMisiXK6yIivxCRShH5TETOT/wwlVJKdUUsmfvTwJIOXr8EGOf/50bgv7o/LKWUUt3RaXA3xqwBTnZwymXAr43Ph0CBiAxJ1ACVUkp1XUYCrjEMOBDydbX/2OHwE0XkRnzZPbm5uTMmTJgQ902PNBzhVPMpJvafyJGGI9Q21zKhf/zXa8fdBMd3QH4R1FVDwQjIGZC46yulVBw2bNhwwhgzsLPzEhHcY2aMeQx4DGDmzJmmoqIi7mv99OOf8vKul1l/1Xru++g+VlWu4oNvfZCoocKRLfDofFhyK7x+K/zNXTBjeeKur5RScRCRfbGcl4jVMgeB4SFfF/mPWcoYg4gAIAiGRPfICVxPEnxdpZSyXiKC+yrgGv+qmTlAnTGmXUnGChISeBMf3AM30eCulLKfTssyIvIcsBAoFJFq4E7AAWCMeRR4DfgKUAk0AtdaNdhQBhMM7iKC9d0ttXumUso+Og3uxphvdvK6Ab6bsBHFyBgTrJhYU5YJ0MxdqVBut5vq6mpcLleyh9KrOZ1OioqKcDgccX1/j06oJlowc7ciAAf+EgiUZbTvvVIAVFdXk5eXx6hRo4LzXiqxjDHU1NRQXV1NcXFxXNewbfuB0Ey9Z8oySikAl8vFgAEDNLBbSEQYMGBAt/46sm9wt3y1TLs7Wnx9pexDA7v1uvse2ze4h0yoIliXuet/xEopG7JtcIe2NXfL17lr2UepXu8nP/lJsoeQMLYN7j1WltHMXalzhgb3FCMiPVAS18xdqVTx61//mtLSUsrKyli2bBlVVVUsWrSI0tJSFi9ezP79+wFYsWIFK1euZN68eYwePZqXXnoJgMOHD7NgwQKmTZvGlClTWLt2LbfeeitNTU1MmzaNq666Kpk/XkLYdilkm9Uyus5dqaS4+09b2XbodEKvOWloX+78m8lRX9+6dSv//u//zgcffEBhYSEnT55k+fLlwX+efPJJVq5cySuvvAL4Avl7773Hjh07WLp0KVdccQXPPvssX/7yl7ntttvweDw0NjZSXl7OL3/5SzZt2pTQnydZbB3cLW0/oOvclUpJb7/9Nt/4xjcoLCwEoH///qxbt47f//73ACxbtowf/ehHwfMvv/xy0tLSmDRpEkePHgXgC1/4Atdddx1ut5vLL7+cadOm9fwPYjH7BvfQmrul69w1c1cqmo4y7FSRlZUV/DwQJxYsWMCaNWt49dVXWbFiBT/4wQ+45pprkjVES9i65m7tahmlVCpatGgRL774IjU1NQCcPHmSefPm8fzzzwPwzDPPUF5e3uE19u3bx+DBg7nhhhu4/vrr2bhxIwAOhwO3223tD9BD7Ju5hzUOs4yullEqpUyePJnbbruNCy+8kPT0dKZPn85DDz3EtddeywMPPMDAgQN56qmnOrzG6tWreeCBB3A4HPTp04df//rXANx4442UlpZy/vnn88wzz/TEj2MZ2wZ3oE3jMGhbquk+XeeuVKoKTJ6Gevvtt9ud9/TTT7f5ur6+Pur3A9x///3cf//9iRtoEtm2LBNaYw8EdEtKM5q5K6VsyL7BPbQs4//oNV5L76iUUnZh3+Ae9oQqWLUbk2buSin7sW9wjzShmsjYruvclVI2ZtvgDu0369DMXSmlfGwb3MMfYgILN8n2X10ppezCvsE9QrC15ClVXS2jVMr5xS9+wcSJE9s1+KqoqGDlypUJucfTTz/NTTfdlJBrJYNt17lHWi2T2Mxd17krlaoeeeQR3nzzTYqKioLHWltbmTlzJjNnzkziyFKHfTP3SGUZzdyV6vW+/e1vs3fvXi655BLy8/NZtmwZ8+fPZ9myZaxevZqvfvWrADQ0NHDdddcxa9Yspk+fzh//+EfAl5F/7WtfY8mSJYwbN65Nk7GnnnqKkpISZs2axfvvvx88/uKLLzJlyhTKyspYsGBBz/7AceoVmXtaj/yO0sxdqXb+fCsc2ZzYa543FS65L+rLjz76KK+//jrvvPMOv/zlL/nTn/7Ee++9R3Z2NqtXrw6ed++997Jo0SKefPJJamtrmTVrFhdddBEAmzZt4pNPPiErK4vx48fzve99j4yMDO688042bNhAfn4+X/ziF5k+fToA99xzD2+88QbDhg2jtrY2sT+vRWybuWNol7lb8xCTZu5KpbKlS5eSnZ3d7vhf/vIX7rvvPqZNm8bChQtxuVzBTTwWL15Mfn4+TqeTSZMmsW/fPtavX8/ChQsZOHAgmZmZ/P3f/33wWvPnz2fFihU8/vjjeDyeHvvZuqNXZO6hxxJ4Ax9d565UdB1k2D0lNzc34nFjDC+//DLjx49vc3z9+vVt2gCnp6fT2tra4T0effRR1q9fz6uvvsqMGTPYsGEDAwYM6P7gLWTbzD18J6bwY1bcUSllH1/+8pd56KGHgnNxn3zySYfnz549m3fffZeamhrcbjcvvvhi8LU9e/Ywe/Zs7rnnHgYOHMiBAwcsHXsi2Ddz1wlVpVQH/vVf/5V//Md/pLS0FK/XS3FxMf/7v/8b9fwhQ4Zw1113MXfuXAoKCtrsznTzzTeze/dujDEsXryYsrKynvgRusW+wT3CUkhrb6iZu1KpoqqqCoC77rqrzfGFCxeycOFCALKzs/nVr37V7ntXrFjBihUrgl+HBvxrr72Wa6+9tt33BLbwsxPblmWAdr1lEpu5h61zV0opG7FtcI+0MYe1/dw1c1dK2Yd9g3uPTahq5q6Ush/7BveQzD1NfD+GpZt1aM1dKWUj9g3uVk+ohvdzV0opG+kdwd3KpZBozV0pZT+2De6h7QfOHtJ17kqp2I0aNYoTJ050eE6fPn16aDSJFVNwF5ElIrJTRCpF5NYIr48QkXdE5BMR+UxEvpL4obbVc5l74IaauSvVmySzR0xn7Q4SodPgLiLpwMPAJcAk4JsiMinstNuBF4wx04ErgUcSPdBw1q+WCdTc09p+rZRKqgceeIBf/OIXAPzTP/0TixYtAuDtt98Obt7x3HPPMXXqVKZMmcItt9wS/N4+ffrwwx/+kLKyMtatWxc83tTUxCWXXMLjjz8e8Z633XYbZWVlzJkzh6NHjwK+B6kWLVpEaWkpixcvDjYlW7FiBS+99FKbewKsXr2a8vJyli5dyqRJk2hoaODSSy+lrKyMKVOm8Lvf/S5RbxEQ2xOqs4BKY8xeABF5HrgM2BZyjgH6+j/PBw4lcpCRtGk/YOlyRS3LKBXN/R/dz46TOxJ6zQn9J3DLrFuivl5eXs6DDz7IypUrqaiooLm5Gbfbzdq1a1mwYAGHDh3illtuYcOGDfTr14+LL76YV155hcsvv5yGhgZmz57Ngw8+GLxefX09V155Jddccw3XXHNNu/s1NDQwZ84c7r33Xn70ox/x+OOPc/vtt/O9732P5cuXs3z5cp588klWrlzJK6+80uHPtnHjRrZs2UJxcTEvv/wyQ4cO5dVXXwWgrq4uzncssljKMsOA0C451f5joe4CrhaRauA14HuRLiQiN4pIhYhUHD9+PI7hnmUwwT7uWpZR6twR6Mp4+vRpsrKymDt3LhUVFaxdu5by8nI+/vjjYOvejIwMrrrqKtasWQP4OkB+/etfb3O9yy67jGuvvTZiYAfIzMwMbgAyY8aMYOuDdevW8a1vfQuAZcuW8d5773U69lmzZlFcXAzA1KlT+etf/8ott9zC2rVryc/Pj+v9iCZRvWW+CTxtjHlQROYCvxGRKca0XXhujHkMeAxg5syZ3YqWkTJ3nVBVqmd1lGFbxeFwUFxczNNPP828efMoLS3lnXfeobKykokTJ7J79+6o3+t0OklPT29zbP78+bz++ut861vfardII3C/wPFY2gNnZGTg9fpCn9frpaWlJfhaaHvikpISNm7cyGuvvcbtt9/O4sWLueOOOzp/A2IUS+Z+EBge8nWR/1iofwBeADDGrAOcQGEiBhhNpAnVhD7E1C5T18xdqVRRXl7Oz372MxYsWEB5eTmPPvoo06dPR0SYNWsW7777LidOnMDj8fDcc89x4YUXRr3WPffcQ79+/fjud7/bpTHMmzeP559/HoBnnnmG8vJywLcCZ8OGDQCsWrUKt9sd8fsPHTpETk4OV199NTfffDMbN27s0v07E0tw/xgYJyLFIpKJb8J0Vdg5+4HFACIyEV9w717dpRPGmGA53NrM3f8WaWxXKmWUl5dz+PBh5s6dy+DBg3E6ncHgOmTIEO677z6++MUvUlZWxowZM7jssss6vN7Pf/5zmpqa2uyn2pmHHnqIp556itLSUn7zm9/w85//HIAbbriBd999NzhpG20zkc2bNzNr1iymTZvG3Xffze233x7zvWNijOn0H+ArwC5gD3Cb/9g9wFL/55OA94FPgU3AxZ1dc8aMGaY7rn/jenP1q1cbY4z5Y+UfzZSnp5j9dfu7dc029n1ozJ19jdn5hu/j6p8m7tpK2di2bduSPYRzRqT3GqgwMcTtmGruxpjX8E2Uhh67I+TzbcD8bv6e6RLTUzX3s3e08NpKKZVYtn1C1fI9VNHeMkop+7J3cO+JbfYCv0B0KaRSykZsG9y9xtuuK6Q1E6qBTzS4K6Xsw7bBPbTmHujnntDgHszUtSyjlLIf2wZ34OwTqlhYlhEtyyil7Me2wd1rvGeT6mD81X7uSp3LYmnh25Xz7My2wT3STkyWrJZRSikb0uDeGS3LKJVyqqqqmDBhAldddRUTJ07kiiuuoLGxEfA9OXr++eczdepUduzwdaysqanh4osvZvLkyVx//fXWNhlMEYlqHNbzQnZiCk6oWvYvTCdVlYrkyE9+QvP2xLb8zZo4gfP+5V86PW/nzp088cQTzJ8/n+uuu45HHvFtI1FYWMjGjRt55JFH+NnPfsZ///d/c/fdd3PBBRdwxx138Oqrr/LEE08kdMypyLaZu9d4g0Fdn1BV6twzfPhw5s/3PRh/9dVXB1vufu1rXwPatudds2YNV199NQCXXnop/fr16/kB9zDbZu5ezq5zt2RCNfRaIlqWUSqCWDJsq4S35w18nZWVBcTWnrc3s23mbozpocxd0LKMUqln//79wa3ynn32WS644IKo5y5YsIBnn30WgD//+c+cOnWqR8aYTPYN7iETqpY8xBThjkqp1DF+/HgefvhhJk6cyKlTp/jOd74T9dw777yTNWvWMHnyZH7/+98zYsSIHhxpcti2LBOxK2RCSydallEqlWVkZPDb3/62zbFAjR1g5syZrF69GoABAwbwl7/8pQdHl3y2zdy9hEyoWtk4TAJlGQ3uSin7sG1wN8bqlr8htO2vUill1KhRbNmyJdnDSGn2Du5WbtYR/leAlmWUCjoXHgJKtu6+x7YN7qFLIa19iEnLMkqFcjqd1NTUaIC3kDGGmpoanE5n3New9YRqjz3EpBOqSgUVFRVRXV3N8ePHkz2UXs3pdFJUVBT399s3uIfsxGRtV8iQGyilcDgcFBcXJ3sYqhP2LctYvhNT+LU0c1dK2Yetg3u6pANna+5e4038jUS0LKOUsh3bBvdIq2Wso2UZpZS92Da4W/4QkzYOU0rZmH2De0jL3wBtHKaUUj62De6hT6hqP3ellGrLtsE9NHO3fCcmLcsopWzGtsHdYCxuHBZ6LX1CVSllL7YN7tavc/cT0ZK7Usp2bBvcPcbTQxOqgYtr5q6Usg/bBvc2vWWsXgqpZRmllM3YNri3mVDFym329AlVpZT92Da4Wz+hGkozd6WUvdg2uHu8Hms369A9VJVSNhZTcBeRJSKyU0QqReTWKOf8nYhsE5GtIvJsYofZnpezjcMsbfkr+oSqUsp+Ou3nLiLpwMPAl4Bq4GMRWWWM2RZyzjjgx8B8Y8wpERlk1YADemwpZJBm7kop+4glc58FVBpj9hpjWoDngcvCzrkBeNgYcwrAGHMsscNsL+ITqlbtoaplGaWUzcQS3IcBB0K+rvYfC1UClIjI+yLyoYgsiXQhEblRRCpEpKI7W3QFyi+Bskwwc9c9VJVSCkjchGoGMA5YCHwTeFxECsJPMsY8ZoyZaYyZOXDgwLhv5jEe4OwqmeD1dQ9VpZQCYgvuB4HhIV8X+Y+FqgZWGWPcxpjPgV34gr0l2mXu2ltGKaXaiCW4fwyME5FiEckErgRWhZ3zCr6sHREpxFem2ZvAcbYRnrlb+hBTYJs9pZSykU6DuzGmFbgJeAPYDrxgjNkqIveIyFL/aW8ANSKyDXgHuNkYU2PVoAN7pVqaubfbH1szd6WUfXS6FBLAGPMa8FrYsTtCPjfAD/z/WC4Q3INPqFq6FFInVJVS9mPLJ1S9tA3ulj7EBP4JVWsurZRSVrBncPe2De7W1Nx1QlUpZV/2DO5hmXug5h4o1ySUaFdIpZT92DO4B2rutK25W0czd6WUvdg6uAe7QooFE6pt2g+gmbtSylZsHdzD2w9YUpbR1TJKKRuydXDvsc06tOaulLKZ3hHcLam562oZpZR92TK4BzJ0S1v+BgiauSulbMeWwT1aV8iE1tyNZu5KKfuyZXAPrnOnh9oPaOaulLIZWwb3YFkmrYcmVDVzV0rZjC2De/hDTMGau1X93DVzV0rZjK2De4/sxCS6zl0pZT+2Du490/IXzdyVUrZjz+BOT2zWoatllFL2Zcvg3m6du5Xb7OlqGaWUDdkyuAdr7oQ1DrN0tYxSStmHrYN7cCcmP8s269DMXSllM7YM7oEg3q7lrxUBWFfLKKVsyJbBPTxzt6Tmrv3clVI2Zu/gTk88oaqZu1LKfmwZ3ANBPFiWsWSdu9bclVL2ZcvgHljnHt7H3bInVCUNzdyVUnZiz+Ae2GYvzfcQkyW9ZcIfYrJkCz+llLKGLYN7sCwTvs5dH2JSSinApsE9am8Zq7pC6oSqUspmbBncg+vcsXJC1U80c1dK2Y89g3tYb5lAWUa32VNKKR9bBvfwfu66zZ5SSrVlz+BOWM09sGlHQuNv6Dr3NF0to1Qvsr+mkbv/tJVdR88keyiWsWVwb7daxuqau3aFVKrXMMbwf367gafer2LFkx/hcnuSPSRL2DK4tyvLWF1z17KMUr3G5oN1bD98mktLh3CozsWrnx1O9pAsYcvgHsjQ00KGL4h1NXedUFWq1/jrtqOkpwn3Xj6FYQXZ/HmLBveUEd5bBnz190Stc2+pPsiB+3/LwXUFtJ6q1cxdqV7kk/21TBySR0FOJl+aNJi1u0/Q3Nr7SjMxBXcRWSIiO0WkUkRu7eC8r4uIEZGZiRtie5E260hU5m5aW6n+7ndp+LSSMweyOfDj/4vxGjRzV8r+vF7DpwdqKSsqAGDemAE0t3r5rLouySNLvE6Du4ikAw8DlwCTgG+KyKQI5+UB3wfWJ3qQ4cJXy/gGkJgnVM+8+SbNO3cy9Lt/y5BZtbi27aJuS71m7kr1AgdONXKmuZWpw/IBmFXcH4D1e2uSOSxLxJK5zwIqjTF7jTEtwPPAZRHO+zfgfsCVwPFFFL5aBnz190Rk7rUvvYxj6FDyZk+m78gmnOPHUPPhKYy39/3ZptS5Zs/xegDGDe4DQEFOJmMG5rLpQG0yh2WJWIL7MOBAyNfV/mNBInI+MNwY82pHFxKRG0WkQkQqjh8/3uXBBoRvsxf4vLuZu+fMGRo+/JC8JUuQ9DREoP/fX0bLSTcN+yz/naWUsljlMV9wHzOwT/DYtOH92HSgzqLNfpKn2xOqIpIG/Afww87ONcY8ZoyZaYyZOXDgwLjvaVXNvXH9emhtJe+LC4PH8i6cR5ozjbotjd26tlIq+fYca6CwTyYFOZnBY2XD8zlR38yhut6VwMUS3A8Cw0O+LvIfC8gDpgCrRaQKmAOssnJSNVJZJhGZe8OH6xGnk+yysmCNPS0rk74T+lK/x4XX1bv+5St1rqmqaWDUgNw2xyYP9dXftx7sXZOqsQT3j4FxIlIsIpnAlcCqwIvGmDpjTKExZpQxZhTwIbDUGFNhyYgJWecelrkHJlrj1bRxI9llZUhmJmdXxwh5E/Pxug0N69Z16/pKqeSqPtXE8P45bY5NHJKHCGw9dDpJo7JGp8HdGNMK3AS8AWwHXjDGbBWRe0RkqdUDjCT4hGoCM3dvUxOunTvJnjat7Qsi5IzKJc0h1L+zOu7rK6WSq6XVy6G6Job3y25zPCczg9GFub0uuGfEcpIx5jXgtbBjd0Q5d2H3h9WxSDX37q6WcW3fAR4P2WWlvgMhvyjSMtLJGeGg4f33Mca0mchVStnDodomjIGisMwdYOKQvr1uxUyveUK1u+vcXdu2AeCcPLn9iyLkjnDgPngQd3V13PdQSiVP9akmAIb3ixzcq081cdrl7ulhWcaewZ0IE6rdXC3j2r6N9P79yRg0KHiXs4TcIt8fOY0fWzaVoJSy0KFaX3AfVpDd7rWJQ/IA2HWk97QAtmVwD+8KCd3vLdO8cxdZ40tCesP7ryUCkkZmv3TSCwpo3KDBXSk7OlTXhAgMzs9q91rJYF9w39mL+rvbMrgnuiuk8XhorqzEWVIS4VXfTkyCwVlWStOnn8Z1D6VUch2qbaKwTxZZGentXhtWkE1eVgY7NXNPrkg19+6slnFXV2NcLrLGjQu9S/hdyS4ro6VyD54zvec/AKXOFYfrXAyNUJIBX/woOS+PHRrckyvaE6rxrnNv3rMHgKyxY9u/KGf3UM2eMgUA1/btcd1HKZU8h2qbGJrvjPp6yeA+7D56pte0IbBlcI84odqNzD0Q3DPHjAm5SdsJVTA4J/maYQZW1iil7MEYw+E6F+d1ENzHDcrjVKOb4/XNPTgy69gzuJv2T6imdeNHadmzl4xBg0jPywu9i//j2cw9o7CQjEGDaNbMXSlbOdPcSmOLhyEdBPfx5/n+/999tL6nhmUpWwb3iOUXiX8P1ebP95I5enTkF/2rZQLBPmv8eFy7dsd1H6VUchz1NwU7Lz9yzR3OtgHe1UtWzNgyuEfK3ONdLWOMoWXv52SNLg5/IeQLAf8vDuf4EloqKzHu3vOwg1K93eFAcO8bPXMf2CeLghwHuzRzT55IjcPiXefuqanBe+YMmaPCgnuEsgxA1rhxGLeblgMHUErZw5EYgruIUDIoj92auSdPxMZhcWbuLVVVAGQWhwf3wIWFwIQqQOZo36Rrc2Vll++llEqOI6d9wT3SA0yhxg7uw65esmLGlsE9kevczwb3UeE3Oft5aObuL99ocFfKPg7XuRiQmxnxAaZQJYP6cNrVyrEz9l8xY8vgHphQTVTmLg4HjiFDwl4JLcukBWvuabm5ZAwdQsvez+MZulIqCY6edjG4g5JMQKANQW9YMWPL4E6w7UvbzD2e1TLNVVU4RoxA0qP8Rg8rywBkjSqm5XMN7krZxeE6V4fLIAPG+YN7b1gxY8vgHsjc202oxpG5u/ftI3PkyPYvtCvLnP0ys9gX3HtDXU6pc8HR0y4GxxDcffurOth9TIN7UniNt01JBvxlmS4GW+P10rL/AJmjRkV6NXjl8Mw9s7gYb2MjrcePd+l+Sqme53J7ONnQ0uFKmYDAipnesBzSlsHdGNMmawf/hGoXM/fWw4cxzc0dZ+7SdikkEPxl4N63r0v3U0r1vKP+lTKxlGXA9zBTb+gxY8/gTvut7oSu19xb/ME5YnDvKHMf5Tu/2b/SRimVugIPMA3p4OnUUCWD83rFihlbBnev8bbrJROeycciGNxHdS1zdwwZgjgcmrkrZQOBB5iGFMSeuYP9J1VtGdwjbVIdz2qZlqp9iNMZsrVeJIHgfvbakp6OY/jw4C8HpVTqiqX1QKjgrkw27+1uz+BOhJp7PGWZqioyR41C0jp7G9qWZQAyR4ygZd/+Lt1PKdXzDtc10deZQW5WRkznF/bJYkBupu3XutsyuEdaLRPPUshAcI+og7IM+Or0Lfv3237SRane7uCpJob1y+nS95QMzmOHlmV6ntd4I06odiXQGreblurqKJOp0NGEKoBj5AiMy0XrMV0OqVQqO1jbxLAY6+0B48/zNRDzeu2bvNkyuEP7zTm6uhSy5UA1eDzte8oEtMnc09pn7sNHAODer3V3pVKZL7jHtlImYPx5eTS2eDhY22TRqKxny+AeLXPvSs29pcrXPiArWjdIwp9QbXvtzJG+4N6yX+vuSqWq0y43Z1ytUTfGjiYwqbr98GkrhtUjek1w72rNPdAbJqaae6SyzNChkJGhk6pKpbCDp3yZ97B+XQvuE/xb7u2w8YoZWwZ3g4lYlvF6Y8/cm/fuJX3AANLz86PexX/liBOqkpGBY9hQWg5ocFcqVe0/2QjAiP5dm1DNzcpg5IAcdhzRzL1HRVzn3sWWv76t9aLsm+q7if/CkTN3gMwRI3Fr5q5Uytpf4wvuI/vndvl7J57Xl+2HNXPvUV687da5d6UsY4yheW8Hm2LD2Rq7pPknVNufkjlihC6HVCqF7T/ZSF9nBvk5ji5/76Shffn8RAP1za0WjMx6tgzuxphudYX0nDiBt66OrDFjOrqL/8JpESdUwTep6q2vx3PqVKxDV0r1oH0nGxnexZJMwOShfQH7TqraM7hHahzWhfYDgS3yssaN7eAmgWtFL8s4RvhXzGgbAqVS0ucn6iku7HpJBmDyUN983NaDdYkcUo+xZXD3mu6VZZp3+4P72I6Ce3jmHrnmDuDW5ZBKpRyX20P1qSbGDOwT1/cP7ptFYZ8sNh/UzL3HROwKSVrMZZnm3btJLyggvbAw+knBmnvgL4QIwb1oGKSlaeZrn83fAAATvUlEQVSuVAqqqmnAGBgzKL7gLiKUFuWz+WBtgkfWM2IK7iKyRER2ikiliNwa4fUfiMg2EflMRN4SkWjP9CdEpNUyCDGXZVy7dpJVUtL+Gm3v4r9uWsQnVAEkMxPH0KG0VGlwVyrVVB7zNf4aHWdZBmDqsHx2H6unwYaTqp0GdxFJBx4GLgEmAd8UkUlhp30CzDTGlAIvAT9N9EBDGdpPqKYRW1nGeL00764kq6SkkxPD17lH/sWROXKkZu5KpaCdR86QniaMjTNzB5g2ogBj4LNq+9XdY8ncZwGVxpi9xpgW4HngstATjDHvGGMa/V9+CBQldphtRa25x1CWadm3D9PYiHPihE5u4vFfOJ1oE6pwNrjrckilUsv2w2cYXZiL05Ee9zWmFRUA8MkB+62IiyW4DwMOhHxd7T8WzT8Af470gojcKCIVIlJxvBubS0fbQzWWskzzjh0AOCdO7OQmHoJZe5QJVfC1L/DW1+OpqYlp7EqpnrH98GnG+9sIxKtfbiajC3PZuK93BveYicjVwEzggUivG2MeM8bMNMbMHDhwYNz38dI+iAsS8Xg417Zt4HCQ2dFKGfCVYdL8v/EljaiZu7/xWIvup6pUyjhR38zB2iZKi6K1F4ndzFH92LDvlO3a/8YS3A8Cw0O+LvIfa0NELgJuA5YaYyzdWdZrvKRL2z+1Yi3LuLZuJWvcWNIyMzu5Sas/qAN0UHMvHgVAs78RmVIq+T6r9q1wKfOXVbpj5sj+nGp0s+e4vXZmiiW4fwyME5FiEckErgRWhZ4gItOBX+EL7McSP8y2ou2h2tmEqjGGpq3byJ48ufObeFohPStw8ahlGceQIUhmJi2fV8UydKVUD/hkfy1pAlMTkLnPKu4PwIefn+z2tXpSp8HdGNMK3AS8AWwHXjDGbBWRe0Rkqf+0B4A+wIsisklEVkW5XEJEnFAlrdOau/vAAbx1dTinTO38Jq0uSA/0o4g+oSrp6WSOGkXL3r0xjFwp1RPWf36SKcPyycmMbd/UjowckMPgvlms32uvebWYfnJjzGvAa2HH7gj5/KIEj6tD0RqHdRbcmz7bDEB2aQzB3d0Imf71sYF7GRPyUNNZmaNH49q+rfNrKqUs53J72HSgluVzE/O4jYgwZ/QA3q+sifyMTYqy5ROqkRqHxRLcXZs/Q7KyOm47EOBuBIe/4VDgX2aU0kzW6GLcB6rxtrR0fl2llKUqqk7R0upl7pgBCbvmvDEDOFHfHHwwyg5sGdyjrXOPJXN3Tp6MOGJo/9naDBmBmnvgXlFWzIweA16v1t2VSgFrK4/jSBdmFycyuPtalXywxz6lGXsG9whlmc4mVI3bjWvbNrKnxlCSAfC4w2ruRF0xkzXW1zq4ZU9lbNdWSlnm3Z3HmTGyH7lZ3a+3Bwzvn0NRv2w+2HMiYde0mi2De6S6V2eZu2vXLkxzM9llpbHdxOuGNH9wD/YO62Cte1oazZV7Yru2UsoSR+pc7Dhyhi+OH5Twa88bM4B1e2rw2GS9uy2De7SukB0G982+yVRnrJm71wNp/t/8nZRl0rKyyBwxgubdu2O7tlLKEqt3+lZiL7QguM8fW8hpV6ttNu+wZXCPus69g4eYmj7bTHq/fjiKYmx74209+4RqJ2UZ8G38ocFdqeR6e8cxhhVkUzI4/mZh0cwd7avhv19pj9KMLYN7tKWQHuOJ+j2uzZtxTp0S+zImryek/UDHq2UAssaV0LJ/P16XK7brK6USqqXVy/uVJ7hw/EBLlisO6utkzMBcPrTJendbBneP19MuuKdLetTM3dvYSPOePWTH8vBSgPGAhPaWgWhlGcDXQtjr1bq7UklSUXWShhaPJfX2gDmjB/Bx1SlaPbHtHZFMtgzuhshdIaNl7q4dO8DrxTllSuw3Cc3cYynLjPf1h2/euSP2eyilEubdXb4lkPMSuL493JzRA6hvbmXLodSvu9syuEeaUE2X9KgTqq4tWwFwxtJTJsB422fuHZRlMkeMQHJycO3YGfs9lFIJ8+6u48wc2T+hSyDDzR7t6zPz0eepX5qxb3CP9BBTlJa/rm3bSC8sxDG4C3+ueT2Q5r+HdJ65S3o6zpISbUOgVBIcO+NbAlle0sG+yAkwKM/J6MJc1u9N/SZivSu4e6MH90435wgXKXPvpOukc9JEmrfvwEQZh1LKGoEVLOVj498nIlazR/fn46qTKd/fvXcF9wiZu7e5mea9e3FOCt/2tRPG07afO3RYlgFwTpqEt6EB9/79XbuXUqpb1u4+Qb8cB5OH9rX8XrOK+3Pa1cqOI2csv1d32Da4t3tCNcpDTM2VldDa2vmeqeHa7cREh2UZOFvTb/LX+JVS1jPG8N7uE8wbW0hamvUdG2f5e9aket3dnsGdCBOqaem0elvbnRvcM3VCF4O7NyRzj2GdO0DW2LGI0xl8GlYpZb3dx+o5dqaZ8rHW1tsDhhVkM6wgm4+qUrvubs/g7vWSntZ2m71oq2Vc23cgOTk4Rozo2k0irpbpOHMXhwPnpEk0ffZZ1+6llIrbml3HASgvsb7eHjC7uD8ffX4ypq09k8WWwd1jIjzElJaOx3javdnNO3bgLClB0rr4oxpvl1bLBGSXluLauhWjvd2V6hFrd59g9MBchhVk99g9ZxX350R9S0rvq2rL4B5pQjWwYXZo9m6MwbVzJ1njx8dxk9CyTGyrZQCyp03DtLTg2qZLIpWymsvt4cO9NVzYg1k7+B5mAliXwksibRncPcYTDOYBGf4Ojm/ufzN4rPXQIbxnzuCcEEdwj6MsA5B9/nQAGjd+0vV7KqW6ZN3eGppbvT0e3EcOyGFIvpMPU3jzDtsG90hLIQH++d1/Dh5z7fQ9LZo1vouTqdB2KWTgozd6Y7IAx6BBOEaMoHHDhq7fUynVJW9tP0pOZnowk+4pIsK8MYV8sOdEyq53t2Vw93q9wUw9IEPaP3Ls2rEDRHxNvboq0lLIGMoyADkzZ9JUUaEPMyllIa/X8Nb2Y1wwthCnI73zb0iw8nGFnGp0szVF+8zYMrhHytwd6Wf3RQ3U3Zu37yBzxAjS++R2/SZeb/vMPYayDEDunNl46uqCyzCVUon3aXUth+tcfHnyeUm5/wXjfEsvAxuEpBpbBnev8UatuQM0e5oBcG3fTlZX2w4ERGr5G+Oyp5zZcwBoWLcuvnsrpTr12ubDZKQJF00cnJT7F/bJoqwon7d2aHBPmFbT2q4sE7oEsqm1CU9dHe7q6q63HQjwes4ugexi5u4YPIiscWOpf++9+O6tlOqQ12v406eHWTh+IPk5js6/wSIXTRzMpgO1HD2depv02DK4R9qsw+11Bz9vbm3GtX07QNcbhgWY0D1U/UE+hgnVgNzyBTRWbMDb0BDf/ZVSUb2/5wRHTru4fPqwpI7jkqm+ktBrmw8ndRyR2DK4e4233QRqH8fZPRNdHhdNgQ2xp3Shh3ubm4Rus+f/GGPmDtBnwQJwu6n/4IP47q+Uiur5jw5QkOPgS5OSU5IJGDsoj0lD+vKHTw4mdRyR2DK4t5rWdpn7paMvZc4QX6272dOMa/MWHMOHk9GvX3w3iVhzjz1zz5lxPmn5+dS/+VZ891dKRXTstIs3th7hGzOKyMro+VUy4a6YUcRn1XVsPVSX7KG0Ybvgboyh1dvaZnUM+CZUl09eDoDL3UTTpk1kl5bGd5PAEsZA5h742IWyjDgc5C1cyJl33tFWBEol0K/X7cNjDFfPGZnsoQDwtfOH4XSk8T8fVCV7KG3YLri3Gl/nx0jr2rPSswBoPnSQ1mPHgk+LdlkgQ29Xlok9uAPkXbIE7+nT1L//fnzjUEq1Ud/cym8+3MfFkwYzckAcS5wtUJCTyTdmDOcPnxzkYG1TsocTZL/g7m/rG565A2Rn+BsHbfDV23NmzozvJoHWwRKeuXftoaQ+8+aRXlBA3apV8Y1DKdXGMx/uo67JzXcWjk32UNr49sIxAPz8zV1JHslZtg3ukTJ3Z7oTgLSNW0nv35+scePiu0kguAdWy6TFl7lLZiZ9v/pV6t98i9ZTp+Ibi1IK8GXtv1qzlwUlA5k2vCDZw2ljWEE2y+eO4sUN1Ww5mBq1d9sF9xaPr36dmZ7Z7jVnhhPxGpwV28mdN6/rbX4DAsE98NdBIMh73JHP70DBN76Bcbupe/nl+MailALgsTV7OdnQwg+/FEc7kR7wvcXjGJCbxY9/v5lWT/Jbj9g2uAfq66GyM7KZeMCQcbqRvMWL4r+JJzxz9wf5CDs9dcY5voSc2bM5+Zvf4tWJVaXiUn2qkcfW7OGrpUMoS7GsPSA/28HdSyez+WAdj6zek+zh2C+4uzy+J8EiZe65jlwWbjZ4nA76LFwY/038v0CCmXt6/MEdYMCNN9B69Ci1z/8u/jEpdY4yxnDbH7aQJsKPvxLnQ4k95NLSISwtG8rP39rNR58nt9e77YJ7U6tvNjo4eRoi7fgp5m8zHLxgHGnZ3diVxd+bhsBfB4FfJK3NcV0ud948cufN5fhDD+E+ejT+cSl1Dnp87V7e3XWcW5ZM6NHdluJ1799OYWT/HL7z2w18fiJ5T6jHFNxFZImI7BSRShG5NcLrWSLyO//r60VkVKIHGtDg9r1ZoU+kAhivl6P33gsCn32puHs3cfv7RDh8E7Rk+D+2xtc/QkQ47447MB4P1StX4qlP3a25lEoVxhh+++E+fvLaDr4y9TyumZsa69o7k+d08N/LZ2KAKx9bx6cHapMyjk6Du4ikAw8DlwCTgG+KSHg3rn8AThljxgL/Cdyf6IEG1Lh8O5/0c/qePPW6fK0GDn7/+9S/+RZvLTmPPbndDJ4t/u/P9P8CycxtezwOmaNGMeyn9+Pauo3Pv/51av/wCu6DBzGt8ZV6lOqN3B4vh2qbeH3LEVY89TG3v7KFheMH8h9/Nw0J9HiygdED+/DcDXPISEvja//1AT966VPW7DrOsTMuPD20uUf79YTtzQIqjTF7AUTkeeAyIHST0MuAu/yfvwT8UkTEWLA1eFVdFWmSRlFeEa0nT7J73nwAJCuLQTffTN2EfXxS9TpnWs6Ql5kX301OH/J97DPI9zHX17eZM90rqeRddBEjnnyCI3ffw+Ef/zh4fNA//5AB11/frWsrZXff/s0GXt96JPh1/9xMbr90ItfOLyY9zT6BPWD8eXm8trKcB/+6kxcqDvBCRXXwtX+7fArLLH7CVjqLvyJyBbDEGHO9/+tlwGxjzE0h52zxn1Pt/3qP/5wTYde6EbjR/+V4YGeifpBuKAROdHpW8uk4E0vHmVg6zsTqaJwjjTGdbhobS+aeMMaYx4DHevKenRGRCmNMnI+y9hwdZ2LpOBNLx5lYiRhnLBOqB4HhIV8X+Y9FPEdEMoB8IHW3BVdKqV4uluD+MTBORIpFJBO4EghvlrIKWO7//ArgbSvq7UoppWLTaVnGGNMqIjcBbwDpwJPGmK0icg9QYYxZBTwB/EZEKoGT+H4B2EVKlYk6oONMLB1nYuk4E6vb4+x0QlUppZT92O4JVaWUUp3T4K6UUr3QOR3cO2urkCpEpEpENovIJhGpSPZ4AkTkSRE55n/OIXCsv4j8VUR2+z/GuYlt4kQZ510ictD/nm4Ska8kc4z+MQ0XkXdEZJuIbBWR7/uPp8x72sEYU/H9dIrIRyLyqX+sd/uPF/vbpFT626a070KY/DE+LSKfh7yf07p8cWPMOfkPvsnhPcBoIBP4FJiU7HFFGWsVUJjscUQY1wLgfGBLyLGfArf6P78VuD9Fx3kX8M/JHlvYOIcA5/s/zwN24Wv5kTLvaQdjTMX3U4A+/s8dwHpgDvACcKX/+KPAd1JwjE8DV3Tn2udy5h5sq2CMaQECbRVUjIwxa/Ctjgp1GfA//s//B7i8RwcVQZRxphxjzGFjzEb/52eA7cAwUug97WCMKcf4BBpCOfz/GGARvjYpkPz3M9oYu+1cDu7DgAMhX1eTov+R4vuX/RcR2eBv4ZDKBhtjDvs/PwIMTuZgOnGTiHzmL9skvXwUyt9ZdTq+TC4l39OwMUIKvp8iki4im4BjwF/x/bVea4wJdOxL+v/34WM0xgTez3v97+d/ikj73Yk6cS4Hdzu5wBhzPr7OnN8VkQXJHlAsjO9vzVRda/tfwBhgGnAYeDC5wzlLRPoALwP/aIw5HfpaqrynEcaYku+nMcZjjJmG78n6WcCEJA+pnfAxisgU4Mf4xvoFoD9wS1evey4H91jaKqQEY8xB/8djwB/w/Ueaqo6KyBAA/8djSR5PRMaYo/7/qbzA46TIeyoiDnxB8xljzO/9h1PqPY00xlR9PwOMMbXAO8BcoMDfJgVS6P/7kDEu8Ze/jDGmGXiKON7Pczm4x9JWIelEJFdE8gKfAxcDWzr+rqQKbUWxHPhjEscSVSBY+v0tKfCeiq9h+RPAdmPMf4S8lDLvabQxpuj7OVBECvyfZwNfwjdH8A6+NimQ/Pcz0hh3hPwyF3xzAl1+P8/pJ1T9y7X+H2fbKtyb5CG1IyKj8WXr4GsX8WyqjFNEngMW4mtPehS4E3gF32qEEcA+4O+MMUmdzIwyzoX4SggG32qk/xNS104KEbkAWAtsBrz+w/+Cr6adEu9pB2P8Jqn3fpbimzBNx5fIvmCMucf//9Tz+ModnwBX+zPkVBrj28BAfKtpNgHfDpl4je3a53JwV0qp3upcLssopVSvpcFdKaV6IQ3uSinVC2lwV0qpXkiDu1JK9UIa3JVSqhfS4K6UUr3Q/wfs7qVpGMkoXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(4):\n",
    "    sns.distplot([beta[i] for beta in bootstrap_betas], hist=False,\n",
    "                 label=['const', 'friends', 'work hours', 'phd'][i])\n",
    "plt.ylim(0,1)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regularization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ridge_penalty(beta, alpha):\n",
    "    return alpha * dot(beta[1:], beta[1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def squared_error_ridge(x_i, y_i, beta, alpha):\n",
    "    \"\"\"estimate error plus ridge penalty on beta\"\"\"\n",
    "    return error(x_i, y_i, beta) ** 2 + ridge_penalty(beta, alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ridge_penalty_gradient(beta, alpha):\n",
    "    \"\"\"gradient of just the ridge penalty\"\"\"\n",
    "    return [0] + [2 * alpha * beta_j for beta_j in beta[1:]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def squared_error_ridge_gradient(x_i, y_i, beta, alpha):\n",
    "    \"\"\"the gradient corresponding to the ith squared error term including the ridge penalty\"\"\"\n",
    "    return vector_add(squared_error_gradient(x_i, y_i, beta),\n",
    "                          ridge_penalty_gradient(beta, alpha))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_beta_ridge(x, y, alpha):\n",
    "    \"\"\"use gradient descent to fit a ridge regression with penalty alpha\"\"\"\n",
    "    beta_initial = [random.random() for x_i in x[0]]\n",
    "    return minimize_stochastic(partial(squared_error_ridge, alpha=alpha),\n",
    "                               partial(squared_error_ridge_gradient, alpha=alpha),\n",
    "                               x, y,\n",
    "                               beta_initial,\n",
    "                               0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[30.619881701311712,\n",
       " 0.9702056472470465,\n",
       " -1.8671913880379478,\n",
       " 0.9163711597955347]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random.seed(0)\n",
    "beta_0 = estimate_beta_ridge(x, y, alpha=0.0)\n",
    "beta_0\n",
    "# [30.6, 0.97, -1.87, 0.91]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.267438780018153"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dot(beta_0[1:], beta_0[1:]) # 5.26"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6800074955952597"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multiple_r_squared(x, y, beta_0) # 0.680"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### With increasing regularization, coefficients get smaller"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[30.55985204967343, 0.9730655363505671, -1.8624424625144256, 0.9317665551046306]\n",
      "5.2837373774215655\n",
      "0.680010213297079\n"
     ]
    }
   ],
   "source": [
    "beta_reg = estimate_beta_ridge(x, y, alpha=0.01)\n",
    "print(beta_reg)\n",
    "# [30.6, 0.97, -1.86, 0.89]\n",
    "print(dot(beta_reg[1:], beta_reg[1:])) # 5.19\n",
    "print(multiple_r_squared(x, y, beta_reg)) # 0.680"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[30.894860179735474, 0.9490275238632391, -1.8501720889216575, 0.5325129720515789]\n",
      "4.607360065077926\n",
      "0.6797276241305292\n"
     ]
    }
   ],
   "source": [
    "beta_reg = estimate_beta_ridge(x, y, alpha=0.1)\n",
    "print(beta_reg)\n",
    "print(dot(beta_reg[1:], beta_reg[1:]))\n",
    "print(multiple_r_squared(x, y, beta_reg))\n",
    "# [30.8, 0.95, -1.84, 0.54]\n",
    "# 4.60\n",
    "# 0.680"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[30.666778908554885, 0.908635996761392, -1.6938673046100265, 0.09370161190283018]\n",
      "3.7035858123105934\n",
      "0.6757061537631815\n"
     ]
    }
   ],
   "source": [
    "beta_reg = estimate_beta_ridge(x, y, alpha=1)\n",
    "print(beta_reg)\n",
    "print(dot(beta_reg[1:], beta_reg[1:]))\n",
    "print(multiple_r_squared(x, y, beta_reg))\n",
    "# [30.7, 0.90, -1.69, 0.085]\n",
    "# 3.69\n",
    "# 0.676"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[28.372861060795607, 0.7307660860322116, -0.9212163182015426, -0.018495551723207087]\n",
      "1.3830006628491893\n",
      "0.5752138470466858\n"
     ]
    }
   ],
   "source": [
    "beta_reg = estimate_beta_ridge(x, y, alpha=10)\n",
    "print(beta_reg)\n",
    "print(dot(beta_reg[1:], beta_reg[1:]))\n",
    "print(multiple_r_squared(x, y, beta_reg))\n",
    "# [28.3, 0.72, -0.91, -0.017]\n",
    "# 1.36\n",
    "# 0.573"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lasso penalty\n",
    "\n",
    "Lasso can force coefficients to zero, but is not amenable to gradient descent. :("
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lasso_penalty(beta, alpha):\n",
    "    return alpha * sum(abs(beta_i) for beta_i in beta[1:])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}