Introduction to Linear Modeling with Python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   time  distance\n",
      "0  0.00  0.132007\n",
      "1  0.01  0.000000\n",
      "2  0.02  0.000000\n",
      "3  0.03  0.000000\n",
      "4  0.04  0.468150\n"
     ]
    },
    {
     "ename": "KeyError",
     "evalue": "-1",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-14-45b12d1d8519>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     10\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     11\u001b[0m \u001b[1;31m# Compute the total change in distance and change in time\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 12\u001b[1;33m \u001b[0mtotal_distance\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdistances\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mdistances\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     13\u001b[0m \u001b[0mtotal_time\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtimes\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mtimes\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     14\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files\\Anaconda3\\lib\\site-packages\\pandas\\core\\series.py\u001b[0m in \u001b[0;36m__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m    866\u001b[0m         \u001b[0mkey\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mapply_if_callable\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    867\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 868\u001b[1;33m             \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_value\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    869\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    870\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mis_scalar\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files\\Anaconda3\\lib\\site-packages\\pandas\\core\\indexes\\base.py\u001b[0m in \u001b[0;36mget_value\u001b[1;34m(self, series, key)\u001b[0m\n\u001b[0;32m   4373\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   4374\u001b[0m             return self._engine.get_value(s, k,\n\u001b[1;32m-> 4375\u001b[1;33m                                           tz=getattr(series.dtype, 'tz', None))\n\u001b[0m\u001b[0;32m   4376\u001b[0m         \u001b[1;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0me1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   4377\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m \u001b[1;32mand\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mholds_integer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mis_boolean\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_value\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_value\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.Int64HashTable.get_item\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.Int64HashTable.get_item\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mKeyError\u001b[0m: -1"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import sys, os\n",
    "\n",
    "os.chdir('C:/Users/T0230575/source/repos/DataCamp__/Python/Introduction to Linear Modeling with Python/')\n",
    "df = pd.read_csv('hiking_data.csv')\n",
    "print(df.head())\n",
    "\n",
    "# Compute the total change in distance and change in time\n",
    "total_distance = distances[-1] - distances[0]\n",
    "total_time = times[-1] - times[0]\n",
    "\n",
    "# Estimate the slope of the data from the ratio of the changes\n",
    "average_speed = total_distance / total_time\n",
    "\n",
    "# Predict the distance traveled for a time not measured\n",
    "elapse_time = 2.5\n",
    "distance_traveled = average_speed * elapse_time\n",
    "print(\"The distance traveled is {}\".format(distance_traveled))\n",
    "\n",
    "\n",
    "\n",
    "# Select a time not measured.\n",
    "time = 8\n",
    "\n",
    "# Use the model to compute a predicted distance for that time.\n",
    "distance = model(time)\n",
    "\n",
    "# Inspect the value of the predicted distance traveled.\n",
    "print(distance)\n",
    "\n",
    "# Determine if you will make it without refueling.\n",
    "answer = (distance <= 400)\n",
    "print(answer)\n",
    "\n",
    "\n",
    "\n",
    "# Complete the function to model the efficiency.\n",
    "def efficiency_model(miles, gallons):\n",
    "   return np.mean( miles / gallons )\n",
    "\n",
    "# Use the function to estimate the efficiency for each car.\n",
    "car1['mpg'] = efficiency_model(car1['miles'] , car1['gallons'] )\n",
    "car2['mpg'] = efficiency_model(car2['miles'] , car2['gallons'] )\n",
    "\n",
    "# Finish the logic statement to compare the car efficiencies.\n",
    "if car1['mpg'] > car2['mpg'] :\n",
    "    print('car1 is the best')\n",
    "elif car1['mpg'] < car2['mpg'] :\n",
    "    print('car2 is the best')\n",
    "else:\n",
    "    print('the cars have the same efficiency')\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   time  distance\n",
      "0  0.00  0.132007\n",
      "1  0.01  0.000000\n",
      "2  0.02  0.000000\n",
      "3  0.03  0.000000\n",
      "4  0.04  0.468150\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnW2MHdd53//PXt6NdHdtRLxkU9X23k0To40RoH5ZuHZdGInlFglb1O2HBA2WDCEZZbS0W9pwUajhh+QLgdZN3PJDSYG1xVDcCweC7SRGS6Rx7ARGvrhdOaosRxUcpFxGiSqRV3HFFzFL7p5+OHtyZ2fPmTln5sydO3P/P2Cwu3fn5dy75H+eeV5FKQVCCCHNZ67uBRBCCIkDBZ0QQloCBZ0QQloCBZ0QQloCBZ0QQloCBZ0QQloCBZ0QQloCBZ0QQloCBZ0QQlrCgbwdROQdAJ4G8NcB7AC4oJQ6KyK/DOBfALi+u+svKqWuZJ3r0KFDanl5udSCCSFk1nj22WdvKKUO5+2XK+gA7gP4jFLq2yLyFgDPisjXdn/3H5VSv+K7qOXlZWxsbPjuTgghBICIbPrslyvoSqlXALyy+/1NEXkRwNvKLY8QQkhsgnzoIrIM4D0AvrX70idF5HkReUpEHoq8NkIIIQF4C7qILAL4MoBPKaXeAHAewI8AeDe0Bf+rjuNOiMiGiGxcv37dtgshhJAIeAm6iHShxXyolPoKACilXlVKbSuldgD8FwDvtx2rlLqglFpRSq0cPpzr0yeEEFKQXEEXEQHwBQAvKqU+l3j94cRu/wzAC/GXRwghxBcfC/1DAI4B+IiIPLe7HQHwWRH5jog8D+AnAXy6yoUSQkgjGA6B5WVgbk5/HQ4ndmmfLJc/ACCWX2XmnBNCyMwxHAInTgB37uifNzf1zwCwulr55VkpSghpL5O2lk+fHou54c4d/foE8CksIoSQ5lGHtXztWtjrkaGFTghpJ3VYy0tLYa9HhoJOCGkndVjLZ84Avd7e13o9/foEoKATQtpJHdby6ipw4QIwGAAi+uuFCxMJiAIUdEJIW5m0tWwCsMeO6Z8vXwauXp2YmAMUdEJIW5mktWwCsJubgFLjAOwEc9ABQJRSE7vYysqKYvtcQkjrWF7WIp5mMNBWeklE5Fml1ErefrTQCSGkLDWnKxoo6IQQUpaa0xUNFHRCCClLzemKBgo6IYSUpeZ0RQNL/wkhJAarqxMX8DS00AkhpCw1tsxNQgudEELKUHPL3CS00Akh00edFm/otWtumZuEFjohZLqo0+Itcu0pyUEHaKETQqaNOi3ekGsbS95VbX/w4MSfMmihE0KmizotXt9rpy35NHNzwGikN2BiTxm00Akh00XRqksf33fePnnXNscfPeoW834f2NnZ//oknjKUUhPb3ve+9ylCCMlkfV2pXk8p7czQW6+nXy9zTNl9bL8L3UQKfSQANpSHxtJCJ4RMF0WqLn183z77uK4NAMePu61yXyru7cL2uYSQ5jM3Zw9OiozdHz772Mjzl/sioodeFPChs30uIWR28PG7u/Y5eDD73DbLPhQR4PHHK0+7pKATQpqPT7fDM2eAbnf/sTdvaivcFTAtm13T72vL/Ny5cufxgC4XQkg7GA61NX3tmrbGz5zZbxEfOjROJUyysKDdMUlLvNfT/vNTp+zHhNDpANvb2idvW1cOvi4XCjohZHZw+dFdLCwAW1vAvXvx1mBuFAGiTh86IWS2OXkSOHBA+68PHNA/h2aZ3L4dV8yBSvPRKeiEkPZx8iRw/rx2cwD66/nzwI/+aL3rMlRU9UpBJ4S0D5M7nub3f18HKcsyGJQ7T0X56BR0Qkj7MJa57fWzZ/dnxIRy6xbwl39Z7NgKZ41S0Akh7aPTcb+erAYtymikRT2Ufr/SWaO5gi4i7xCR3xORF0XkuyJyavf1gyLyNRH53u7XhypZISGEZGHLHzedDdMkOx5evQqsr5e31kN4881KT+9jod8H8Bml1I8B+ACAT4jIuwA8AeDrSql3Avj67s+EEBKXrA6JRrw3N3U64uYm8OijwDPP2M915YoOmJrznT6te7SUsdZDqLrjok8Hr+QG4LcA/AMALwF4ePe1hwG8lHcsuy0SMsOsrys1GOiOg4NBdvfE5DHz83s7Fs7Pj48dDMp1P0yeL0Y3Rd8tEHh2WwwqLBKRZQDfBPDjAK4ppX4w8bu/UEplul1YWETIjGJrcOVTYOOq7Oz3gRs3wguFXJjzmWrTzc3y53TR6QD37wcdEr2wSEQWAXwZwKeUUm8EHHdCRDZEZOP69eu+hxFCJsGkhjEXHSvnKrk3r8dK/xuN9PsHqvetuzJwIuAl6CLShRbzoVLqK7svvyoiD+/+/mEAr9mOVUpdUEqtKKVWDh8+HGPNhJAY2PzPJ05UI+pVjZWzNeUqSvL9p/uix6RCf71PlosA+AKAF5VSn0v86qsAju9+fxzat04IaQpVD2NOWv9zDqnJs7BdxTvm9bTw9vvA/HzhJePOHd2My5z76lXdKTEWFeagA8gPigL4+wAUgOcBPLe7HQHQh85u+d7u14N552JQlJApQiTqmLQ9+AQY88bKmfN0u3uP63bzx9GZ4GunUyxoac4fI1Bq1uAbCLYAz6Bo7g4xNwo6IVOEK0NkMKju3J1OWJaLUnq/fn98jn5//7FJEe/39ea6Yfls5jOIkUVjW28gvoJ+oDrbnxAy1Zw5Y888ieEScPnGd3ayx725SBbkjEZ7C4TSGTRle5cD4/XHaKI1GgGPPaa/58QiQkglFBnG7IvPSDhfsnz9w2Gc4c1pjN/f5fsPZWur2oKiXSjohMwyJvC3s6O/lhVzEwjd3NyfHeJj/dvSKF1W8uYmcOxYNWmA29vaYRLz3BW1zE1CQSeExCGZBgloQTSi7mP9nzypBTqdRpk1xFlFKCqaFBW1zE1CQSeE5ONTgGRzjSilxTzP+h8OgSef3C/Q5nyTbKBVFbduVVe4tQsFnRCSjW8BkqtcPu1qsN0cTp92W9ujkfaTxy7wsRFj+IULE8ytUNQp6ISQbHwKkIZDt+AmXQ2um0NW75S5Od0lsWr3ymBQbPhFt6uHSftQcbfFoOZcZWFzLkIaiKsBlsg4BdEEQm37XL48dre49ut0Ku1x4sXcnH4//T5w964eEO2LiP8NJ/m5eZ8+cnMuQkjL8G3M5ZOC6MrgUGqv79y1X91iDoxFdjTSYu6aemQjxDCuMDhKQSdkFglpzGVrgJVOQXSJVLoRVdZ+vm6LSRHjJlMkdbMEFHRC2kTedB/zO1sxjsu/61OAlCX6yeveumVvnnXjRrH3O830+3tTLiueJwqApf+EtIZ0Cbyxug3J37msT5dLZHU1W4jM706f1udYWhpboumy/G5Xi/rW1vh4m7/a+LR9MRkqMUr/y9LtAjdv7n2PFc8TBRgUJaQ9uAKOxu3hM4XH5IzbMOmFScFOi3x6n1u3igtsv69F0Kesv9vVTw9JAZ00nY6+AWW976zPNwMGRQlpAyEThbJK5H3E3OYiMdfM87kPh3pc3NGje/cpYy2//vpeV48rSNnpAG99a71iDowbj129qtduo+ryf5+WjLE2ts8lJABbL+6sHuJFWr2m29m6rplsX5tuM1vlcOVkm92s/u1ZrXJ9e6K73mNoy92sv0XB1sRgP3RCGk6oKKyvh/UAF9l/cyh6U6hCzJM3FNPv3PV5uH7n2xc966bls83P7/0sQ2/GOfgKOl0uhNSNy60SOodzdVVLhy9K7feBF3EJxMohd1WamuybI0fsvx+N3C6lmzfz0yFN1o7LTZLH4iLw1FN7P8sqWxNnwKAoIXWSzkwBtC/7wgUtYq4gpyuw5gqM2rCdx3V8vw+88QZw757fuUNYWPCryqyqmnQw8P/MXMcXCHSGwKAoIU0gq0/KkSPuwhSXVW/LBzdpgrbzpHHlk589W274sotu17/EvgoxFykn5sBE+px74+OXibXRh05Iiiz/rivQ2O9rn63LP5ucr2mCimtr+33drrmetuOVCvMpLyyMz1FmtmcTthgzWHMAg6KENICsYcqxhCUrCyUvUJcU99D1JG8w3W49Ylv1zaREoDMEX0Gny4WQOnG5OIq4F5KP/nll/oY7d3TuuC3HPZ17HorJU19dBS5ejDef05dOB/jIR+L1Ue/3gbW1iQc6g/BR/VgbLXRCLNhcHEXSB42FXjQvPG1tFllD1pomaaX3etrNFDM/fgKuFRegy4WQKcLll87aP1SM+v3iNwObaMVwV4joc8W4OYRsi4vxz2neSw34CjpdLoRUTUirWkMyjxnY7zawTckxI87KZG0k3TYx+nabboOhmSBl3SS3bhU/1tViQCm9rkOHKp8NWhQKOiFV4zPCzcbqqs5vVkpP/Un6bi9e1MKS5s6dsMEMaZLtXl2FPIC/4N68qcUv9OagVNj+McmLX4xGwGOPTaWoU9AJqZrQik8bRtwvX9Y/HzvmtsS3t+2BVhPQy8IIMKDneLrwFdytLXdOvYvBIH+dZfH5PLJujOZ9AWEN1KrGxy8Ta6MPncwksRo1+frVk4FVm88+z59t1hUz5S+dN+/aTGC2ioZfJhU0/XkUbfpleuFE7NniAgyKEjIlxPpP7xNY9DlvnlDXFcjsdHRmirkRLSwoNTfnf7xIftGVjaxOklmfQdbvI2fE+Ao6XS6EVE2sRk0+Lpq7d+155Um3QF4+uPF323LkY+V0p+n1dED30qVx8Pj2bf+JRZ2OdkfduKEbZZnpRQDw4IPu44ZD3aMmzfy8fv9nzugAtOv3MdxpMfFR/VgbLXRCShBqMYtoi3dtzd990u2OW84OBnstZvNzbMt8YaF8L/Lk+w15Ispqu2tYX9+7PpMemnV8TRY6BZ2QGBTJM0/vbxOOpKDa3AkxNuPW8HVXxBDftBDHPFeWCyVNlv/c9+/eJB86gKcAvAbghcRrvwzgzwA8t7sd8bkYBZ20ktD/1Lb9u12//i3dbphf2WfLGyCRtFbN+vPWMK0NudI33CwL2/cmHXozL0BMQf8wgPdaBP1f+1wguVHQSSsJfewuG2yswlLPm+yTnsaTJ+hra/Et+eRWdkpSVjaNq23AhBpx2fAV9NygqFLqmwAKjvIgZAYIDYyVDZiNRvEHIo9G2ZN9jh7VhUyLi/r7rGBlvw+cO6cDlOvr8XPKBwMdPHUFaPv9/cHcNKawyxWwvnKlWDFYzZTJcvmkiDwvIk+JyEPRVkRI03BVQaZfN5km+il3+sgrlx+N8odRiOhhGIbVVZ0NElK9mpeFc+2aPu/jj7tF/fjxsUhnnces8epVfZO6elX/PG3ZK54UFfTzAH4EwLsBvALgV107isgJEdkQkY3r168XvBwhU4yrBW5yIlCyn4uNbrdcyf40IKJFNpmOad53zGlD5kZ57pxOVUymKAL6xnPpkv78d3bcTwiuG25WameM/jZV4uOXAbCMhA/d93fpjT500lryAmN5BSrGnxs74GnbqrhGp7P/Pa+vl/d1pzdTnenz2Wa1E077w30qUxvgQy8k6AAeTnz/aQC/7nMeCjqpnVgZCaHn8U2Pm1R2SMyy+vl5u5jH7n9u8syLfLZFb7idTqXZK75EE3QAX4R2q9wD8DKAjwO4DOA7AJ4H8NWkwGdtFHRSK7FyhoucxzcTJiQDxghZqBUcu0Co2w3PVZ+bC8uCsT0BhH62tr9j3ni9GnugJ4lqocfaKOikVmJV9YVYc8n87rRw2G4CeY/+xl2SvEaIOFcxySf9Ga6v5+8/P++/bp+5p6E32JBGZ1MABZ2QNGWrAvPOEyKc/f7eEnvjQzfi7+oMqJR9P9c1bNWnVeWHh+bX5+0f4uoIdYHFanQ2ISjohKTJsqxD/uP6CleW2KYLg7pdv7L70LayVRb3lN2M+LrEvkry3Cw1+8zT+Ao6uy2S2cGWXgjolLq8kXA+57Gd10W6MOjevf2v2QpZbNOPXIjoFL5pRSldfJQmnfLposxgCVf64WCwNx+9YVDQyexgqgJt+d6hVYDJlqx5hTBl2NzcK1S+hS0iWjCrotPR1yibO58uVBLRRUFZYjoc6qrVo0fD5rQm8akdaCI+ZnysjS4XMhWU8aW7GmtV0QUx6XoxXRd99q96MEUyF7yKNMssd0teK+AQV80EmmrFAp4uF1FV3sVTrKysqI2NjYldjxAry8v2is3BQD9qFznWVCsWdXF0u9o6dfVo8bG4e73x4AzXOstiqkHPndM/V3EdEXuvmOFQz1LN+hxcxzYcEXlWKbWStx9dLmT2KPO47XJ5lPFVDwbAxYt60o4Ll4gZl0d6CpKvn9+Q1fNkbm7cuOry5bGYF7mODy7/9unT+Te1aS/NrxgKOpk9yoyEyxKMoqJ+7drYfx/SmVBE+407HW0lHz8OfPSj2mo+dkz7+fv98XvMIksod3a0kLsChVkj3gD99OHbdTHrxpr3JNAGH3hZfPwysTb60EmtuHyma2vjFMNOR6lHHnH7VotMo+90/NMHFxb2++NdPuOFhfzzJVMfy/jW00MuQj8Ln2pWW0FWcvRdlu88ORauhYB56IQkKBPMtDVy8hXCrEEKrs0219NWCenbZCuvSZXvzSY9XzRmjnu60jS9ziwxt/V3aRkUdEKSlM38KNJzxWbd5/UOcV3PtzrUtuU1qfK92VTZOCy5xtC/1QzgK+j0oZPZoOxggvTxPr7azU3g1CmdMz03p/3keT26XdczgyJ6vfDe4qa/9/Ky/jk9zCEZU8hCqbDrhrC0NC4UCsmaaXoP+cgwbZHMBmXT65IpjcOhFuci5zOphYAujHHR7++vooyRItjt6oyarABwpxMv9W9uLv9cvZ4O6F665F8Fm2SCGlYXTFsk7aRouXeZ9Lpk9sRwCDz6aHFhTc6yzKowvXt3/2u+1zRWqy0V8d49/dSQRYiYp6cFpXn66f2fe7e7N/vGNcMziSutMva80qbj45eJtdGHToKwZTqU6Wee7jxo2+bmsrNcYgQCjb84b7+0/z3Eh53nE3dl8IT2Y886JhmMzavIzKv+LPu3bzhgUJQ0mpBMh5Byb5f4iOxvZ2sjRgDQrDdPPPv9sRDGHuWWFsUi6ZiA+28VKrY+XRcbVKofGwo6aTahlqIvvlauS5Dy1pF3XpOTbm4gVQl1yFZ0HTHFNtY0qZbiK+j0oZPpJCQrxbfcezj074xo676Y569XKv+829u6olQp/XVnJ7vsfhKMRuGZM+mqzNXV/dkzIZSp3iV/BQWdTCcukU6LX1a5dzKAeugQ8NhjYcKVDkKGtNf1RSm9vtj9UJKYzyxWil9VYlv2pkAo6GRKcTXQevzxcWZDpzO2pNPW83Co+5yYftmjkb2TYaeTnalhcsir6l4I6JvM8ePV5VQrpT+zS5fK3zjMZ3XsWPhQCVI9Pn6ZWBt96CQIl1/Wx9/q64M3mRo+/u8qKyV9K0iLbibOkP5Ms7J2io7JI9EBg6KktfikyYUG9kKEsSrBLdqeYGEh/9jk3NSkqLsCovPz/r1bqp7/SbwFnS4X0jxcAdNr18auFh+S/nffAhWl8otpirC0VLz46c4d7XNWSrepdc1NPXZM+9OPHdvrirLFFba2dLFP0qf9+uv265dtq0CiQUEnzcMVMF1ayh+ivLhoz6LwFdN+H3jzTf+1Jq+3tgbMz+/fp9cDjhwZrz3Ul37w4Pj7rLmpSu39mocRahNcdh0340MlpgkKOmkeWROH8qzFH/iBscUJjLNgTp/2C0zevu3fb6TfB27e1Nc7cwZ45plxYNakTw4G4z4mJugamkJ48+be4OTqavg5bBw8uHcYs40yQyWKtnEgbnz8MrE2+tBJJiHFKa59fUrzXZWRRSolXVsyWLi+rgOKVfjebX7s0DYBts2nV3yZak0WEgUBBkVJZcQswU72AkmLUOh/8BDhdFVGlqncNOtPfyZle7H7XjvW9USUWlzM36/MvwPXGjudmSztz4OCTqohpmXl0z8kRp+WosIcss3N7Z2cEzLMIobQJj+nkBTMMjewMimMPmukxf5XUNBJNfg0USp7rrTwJMmyCmOIpunsV+TYso2u8rZHHrE/gczP+z0RJK1fk5IYe42+/w58r82USKWUoqCTinBZViENsvLO5foPnfd0ULbRVXKWZ9E2uYNBNUKZl6ceUnhV1Q0n5N9ByNg7QkEnFTFpC31xcSywrqHIocVBVW4i1RYfZf0+Kdpra+MbXKez1x1UpU8/5N9B8mnLdTOmha6UUt6CzrRFEkZWyqAPyVS1W7fsedlJbt3S/7VNZ0IbpqCo7q6FgM7JriovW6ns39+5o1MgT57UaZAmdXF7G3jySf355PWkcaVt+ny2oSmMyWZctj4zZVIiZ5U8xQfwFIDXALyQeO0ggK8B+N7u14d87h600FtC0ewG22N20s1R1GVSlZsjdOv1tCUcY6pRlVvWoBCXu8b1vlyZPZP8dzUDIJbLBcCHAbw3JeifBfDE7vdPAPj3PhejoM84ee6aIq4K42aosrGVazM9VLJG5E3rlpUi6hLWmO42EkQ0QdfnwnJK0F8C8PDu9w8DeMnnPBT0hhJjGk2eBV3Ewp6byxebqrbkk4X5TMpY5nmdD21b2ZtH6N80ZkCcBFG1oH8/9fu/yDj2BIANABtLS0sTefMkIq7ZnskgW+jxsbYDB/IzNw4cmIzAl705FHnKWF8v56YKhRZ6bfgKeuVBUaXUBaXUilJq5fDhw1VfjsTG1uxKKR1k8+m9kdcsK4t+H1hYcP/+/n3g1Cm9Dtd1RIp3R6xi4MTCAnDgwN7X7t3T/VJ8x+MB4319hlaETHnKomxAnFROUUF/VUQeBoDdr6/FWxKZKlzNrpTSGRV5ol60taoIcOOGznLJYjTSo+VcmRv37o3PF8rOTvxWuXfv6huRjZCGWjs74zbBFy7Y19nr6Xa6ly/7TXnKg3M/px8fMx77XS7/AXuDop/1OQ996A0kzzedV55d1Lfd71c/xacNW7rwyuUXZzOsRgNPl4vofd2IyBcB/ASAQwBeBfBLAH4TwDMAlgBcA/AzSilH9/sxKysramNjo8h9h0wS48K4dk23UB2NsvcfDMbtaG3neuwx+zzPLObnw4+ZRUTc+flJXPnnWX87MjWIyLNKqZW8/Q7k7aCU+jnHrx4JXhWZfszEH+OPHo20vzZLNPLcKjlGwz4WF/NdLb70esCDD+bflKYJM0QjGRMQsX+OvkVMWVOeSGtgpSjZiy24aHzJriChUu4BBadPj/3YvsQS88VF7eM9e7b8tPtJIaLXm/ZVP/54uYBk1pQn0h58/DKxNvrQG0BWrnFeCqLNJxvqA8/aP9mbxPdcybTGss27iqw55Li8dNAy9QD0oTcasDkXKURW61XTGjav65/P+UK3tPgUCRpWGWB1DenIG7gxyVJ3ltY3Fl9Bp8uF7MU1LHl7e5wmd/WqOw0w7ZN15S6vrYWlBBZNj0uup4x7IStHvNPR71OpcYqgcZVcvOh+nyYgaWacVp3+l2yGNYnrkcnjo/qxNlroDSHLPWEs3qyqwbQlaIYpuNLpQqxsw8KCv+WcNebOZ+v388vys9wXRdwdtKZJAtDlQkqR17fDNr+z27U3qMoSr6JC6dP3xLwHm68a8POpdzr5w5KzbjyG0Jsc/d0kAQV9Voll2eX5ydfW9gvd/LxbaF1i52sh+/jvQzZXq9gyW5lpPUnBZs8UkoKCPotMeoBzqNglBTk04yTEUg4R39g3CR+KthFmV8OZxVfQGRRtE7YcctO3I5Rk344YKKUbUJlqxZC+JWb/mJWjJkBqAoXr62HHl2l4lVfkw5xxUhAKepuIXQ1oxK5t2MQ35KbX6+lCn6INr/IEm10NSUEo6G2i6ZbdJGaCujoE5s3Z7Pf3dhg8d24svOZpY3NTp3bmiXqeYLOrISkIBb1NVGXZxW4ha6PT8WsyVQV5A6a3t3UbgXT+dlEXl49gM2ecFICC3ibKWnbDoe7JIqKHMJgp8T/7s9UMe0hirNxYPnsXNiv69Gnt48/C5rYq4+KiYJMKyG2fGxO2z51i0l0Wk/R6ulIyVtMsG4OBfpI4dWoynRGTbWPn5vIF3dZm9tAh+1r7fT2cg5BI+LbPpYVONFmj4u7cqVbMAS2Ajz22XyBDxrKFkLSiDx7M3jfptjJPMXNzwOu5IwAImSgU9DaRFBvTztb2mo26+2Lfvm1PS9zZCQ+WiugBGVmYQPFwCLzxhnu/pNvKPMVsbo6zw21Q6ElN5A64IA0h7TLZ3AQefVSLmxFK4z8G9vtsl5ayMz1s9HrFB0CHoJR7wEOabhd461vz3TbG4nb1a7e5TXwHXjclq4i0DlrobcEmNvfu7bd679wBjh3TAimi/cDDYXgmTKfjHk5cBT5iblIL88S83x/f0FxPJjYr2+cphvnipEYo6G0hxGWSFMfRSPuuAeCBB/zPsbOjRfHsWf9jqub738+vJu319q45JHfftW+nw3xxMhVQ0NtCmcf8rS3gF34BuHu3+utVadHntRPo9/cLbkjuvmvfS5eYfkimAgp6W3ANpvDl9m3/fbtdfb3hEDh+3P+4wUD7pZXSvVPyApexefPN/a+F5O6zgpNMOcxDbwPDofahb27qx//QxlehzM8DH/+4tkxDgqLr63vF7+RJ4Mkn/fzjsWCOOGkgzEOfFZKpdIAW86p7omxtAefPh4m5rdL0ypVqxDyrqnU0CkvnJKRBUNCbSFKMjh/fL6yxRDKmS2R7W681KaBV5L4bn3ZWC4FTp/bmk/s21SJkyqGgN410cUuV7pW3vCVub5Xt7b0CmlehGUoy6JmVOjgaxesbT8gUQUGfRrLcAb7FLTFK5l9/HThypBoXzp07OqumTCA3zeLi2Ee/uhqeUVN3tSwhJaGgTxtpCzztDvAVnRgi3O1WG7S8fRv44AfjPQWkP5uzZ+1phi6hZ4UnaTgU9GnCpAG63AHDob/lHcMVs7XlL+YiY6EMabX7jW9o94hJZbRZ7AsLfv78tCC70gxdQs8KT9J0fAaPxto4JDoDn6HMMYc2V7H1++P3EjIEOjlc2QxtFtFfzYDr5Ov9/v6B0aHDsF3XIWQKgeeQaOahTwvLy9nNsebm6pvoE8L6ursDIGEEAAAK0klEQVSvehah/w5N7v21a9oyP3OGBT6ktTAPvWnk+cabIOaAf9A2SZFpSNM88Yc57qQmSgm6iFwVke+IyHMiQtO7DG0JyBXJFNnebo/o5QW1CamQGBb6Tyql3u3zOEAyKNuLZVoomi6ZJXpNsniLDo4mJAJ0uUwLJiOjqpFrk6Jodo1L9Jpm8ZYZHE1IScqqhwLwOyLyrIiciLGg1pK2Mk+eHP986JDejh5tjq+8Cmyi1zSLN6S/OiGRKSvoH1JKvRfATwP4hIh8OL2DiJwQkQ0R2bh+/XrJyzUUm5V5/vz459FoMpPupwVXENQmek2zeEP6qxMSmVKCrpT6892vrwH4DQDvt+xzQSm1opRaOXz4cJnLNZcimR9NoYiLaHvbLnpHjuz3lTfN4mXPdFIjhQVdRBZE5C3mewD/EMALsRbWKqbVmozBgw+Gd2U0IpcUvePHdZfEtK/8yJHmWbzTnFJJWk0ZC/2HAPyBiPwvAP8DwH9TSv12nGW1iJBy/SZy+7YW4JBGWKYIKCl6V67YfeUXLmixp8VLSC6sFK0S4zuP4W6ZxCSiMgwGWqjz3u/CAnDr1v7X5+bc1aIi+nfmGhRzMmOwUnQacPnOOx1gbc2/y6AZ2rC+Hm9txtoNPcbFtWt7/ccu7t2zpxxm+cSN0E97yiIhNUMLvUqyrE7DwoIWua2t7H0eeCBuJoyvRQ1oIV9ayu41Mxho10mSQ4fsa7btG/I0YzuekBZDC30S5FUw+mRi3L6dLeZmn5hiboKKPhb1YDD2c7v2E7EHKV9/3b6/LUhs1uLT16XNQWZCSkBBL0peBeNwaPcV10WnYw8qmuCkrRd5OpvElmMtAjz+uN2v7bqhHTxovxGurmrXUl4LhGlNWSSkZijoRTl1KnsQxYkT01UslBdQNRZyMlvlwQft+yQzTi5fBs6ds5/TdgPodoGbN903wvRTQ9pvP+0pi4TUiU/T9FhbawZcrK9nD2yYm6t/2ETeoAzbQAfbkI3QwRG2cyYHSfT7+UMuso7nIAoyg4ADLiKTHKgg0vyeK7bAomvIRswgpCtQ3IbPlJCKYFA0Jml/+bQLz8JC/pBoE1hMBnZdWSwxg5BNK+UnpEFQ0H1oUi8WEZ3imPfktbSkOz4eOza+UWXtGws2ryKkMijoaWypiE1JkzMZJ650QYNphPXkk/nCH1ts2byKkMqgDz2Jrbil19PZHtOUsWIjWRafNXDa7Hf6dHahkCkmYqk9IbVDH3oRXMMU7t6tZz2+mKClEV6XW2N9fbxf1lNHspgIaM74N0JmHAq6YTh0W6y3b092LSHYqjR93Bouv3jyfGXGvzVpDighbcEntzHWNrV56Lb86yZtsd6ziFJra+N9BoOwnPGsc5fNZydkhoFnHjot9OEQ+Pmfb04Wi41Dh8ItYJ+qz6Lj35o2B5SQljDbQVEj5tOWVz4YADduhLl65ueBp56KG8AsWmjE4iFCosKgqA+nTk2nwGxuhj8xbG3Ft4CL5oyzeIiQWphdQT95crpTEYs8OW1uxg1CFs0ZZ/EQIbVwoO4F1MLJk8D583Wvohh5o+iS2ShAeRfM6mr4Ocz+pvcN89kJmQiz6UPvdKbT1ZJHr6cHJn/+83rKUR6c7ENIK6APfTjU2R8iejOZICdPNkvM04Mpzp0DLl7c27fcRVNaFhBCotBOC304BB591M+KrZKFhXJFSb2en896Em1vCSG1MbsWuklFrFvM+339VJBFt6v3Mxb42lqxplUMQhJC0LagqLHM63ap9HrA2bO6Na2LZDOtsjAISQhB010uySlCBw/Wm4YoojNMfLoe0hVCCAnA1+XSXAs93eq27pzypaX9In3mjL0dL10hhJAKaK4PfdqmCNkySjjMgRAyQZor6NOWkucqa19d1Za76S8+CTFn61pCZpJmCvpwWKw0PpReT2eeJHO+FxZ0I6z0ftPiRinTw5wQ0miaIehJi3NxETh6tLpr2Qp5btwYd/a+dUt3NZxWNwpb1xIys0x/lottzmdV+BbyTDNsXUtI62hPYdGkgp9zc80Xc4CtawmZYUoJuoj8lIi8JCJ/LCJPxFrUHrIm08difh54+unmiznAqlFCZpjCgi4iHQD/GcBPA3gXgJ8TkXfFWhiAOIG8bjf794NB/Ek/dcJUSUJmljKFRe8H8MdKqT8BABH5dQAfA/BHMRYGQE8UKkOnozsTnj49WxWbRXqYE0IaTxlBfxuAP038/DKAv1tuOSnKVn9eujQWNlZsEkJaThkfulhe25deISInRGRDRDauX79e4nKBPPLIWMzphiCEzABlBP1lAO9I/Px2AH+e3kkpdUEptaKUWjl8+HDYFXyGOAA6qNnp6O87HV0M9Lu/u3efOio2CSFkgpQR9P8J4J0i8sMiMg/gnwP4apxl7XL2rLaos+h2dVDz/n2df33/vi4GIoSQGaOwoCul7gP4JID/DuBFAM8opb4ba2EAtBV9+bIut09iRH4w0EFPWtuEENKASlFCCJlx2lMpSgghxAsKOiGEtAQKOiGEtAQKOiGEtAQKOiGEtISJZrmIyHUARdsnHgJwI+JymgDf82zA9zwblHnPA6VUbmXmRAW9DCKy4ZO20yb4nmcDvufZYBLvmS4XQghpCRR0QghpCU0S9At1L6AG+J5nA77n2aDy99wYHzohhJBsmmShE0IIyaARgj6RYdRThIi8Q0R+T0ReFJHvikjJWXzNQEQ6IvKHIvJf617LJBCRHxSRL4nI/979W3+w7jVVjYh8evff9Asi8kUReaDuNcVGRJ4SkddE5IXEawdF5Gsi8r3drw9Vce2pF/SJDKOePu4D+IxS6scAfADAJ2bgPQPAKehWzLPCWQC/rZT62wD+Dlr+3kXkbQD+FYAVpdSPA+hAz1FoG78G4KdSrz0B4OtKqXcC+Pruz9GZekFHYhi1UmoLgBlG3VqUUq8opb69+/1N6P/ob6t3VdUiIm8H8I8AfL7utUwCEXkrgA8D+AIAKKW2lFLfr3dVE+EAgAdF5ACAHixTzpqOUuqbAF5PvfwxAJd2v78E4J9Wce0mCLptGHWrxS2JiCwDeA+Ab9W7ksr5TwD+DYCduhcyIf4mgOsALu66mT4vIgt5BzUZpdSfAfgVANcAvALg/ymlfqfeVU2MH1JKvQJogw3AX6viIk0QdK9h1G1ERBYBfBnAp5RSb9S9nqoQkX8M4DWl1LN1r2WCHADwXgDnlVLvAXAbFT2GTwu7fuOPAfhhAH8DwIKIHK13Ve2iCYLuNYy6bYhIF1rMh0qpr9S9nor5EIB/IiJXoV1qHxGR9XqXVDkvA3hZKWWevL4ELfBt5qMA/o9S6rpS6h6ArwD4ezWvaVK8KiIPA8Du19equEgTBL36YdRThogItG/1RaXU5+peT9Uopf6tUurtSqll6L/vN5RSrbbclFL/F8Cfisjf2n3pEQB/VOOSJsE1AB8Qkd7uv/FH0PJAcIKvAji++/1xAL9VxUUOVHHSmCil7ouIGUbdAfBU9GHU08eHABwD8B0ReW73tV9USl2pcU0kPv8SwHDXUPkTAI/WvJ5KUUp9S0S+BODb0Jlcf4gWVoyKyBcB/ASAQyLyMoBfAvDvADwjIh+HvrH9TCXXZqUoIYS0gya4XAghhHhAQSeEkJZAQSeEkJZAQSeEkJZAQSeEkJZAQSeEkJZAQSeEkJZAQSeEkJbw/wEm8jGuhUm0xAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "ename": "NameError",
     "evalue": "name 'model' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-58cd5554ab9a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     25\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     26\u001b[0m \u001b[1;31m# Pass times and measured distances into model\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 27\u001b[1;33m \u001b[0mmodel_distances\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtimes\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmeasured_distances\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     28\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     29\u001b[0m \u001b[1;31m# Create figure and axis objects and call axis.plot() twice to plot data and model distances versus times\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'model' is not defined"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import sys, os\n",
    "\n",
    "os.chdir('C:/Users/T0230575/source/repos/DataCamp__/Python/Introduction to Linear Modeling with Python/')\n",
    "df = pd.read_csv('hiking_data.csv')\n",
    "\n",
    "print(df.head())\n",
    "\n",
    "times=df['time']\n",
    "distances=df['distance']\n",
    "\n",
    "# Create figure and axis objects using subplots()\n",
    "fig, axis = plt.subplots()\n",
    "\n",
    "# Plot line using the axis.plot() method\n",
    "line = axis.plot(times , distances , linestyle=\" \", marker=\"o\", color=\"red\")\n",
    "\n",
    "# Use the plt.show() method to display the figure\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "# Pass times and measured distances into model\n",
    "model_distances = model(times, measured_distances)\n",
    "\n",
    "# Create figure and axis objects and call axis.plot() twice to plot data and model distances versus times\n",
    "fig, axis = plt.subplots()\n",
    "axis.plot(times, measured_distances, linestyle=\" \", marker=\"o\", color=\"black\", label=\"Measured\")\n",
    "axis.plot(times, model_distances, linestyle=\"-\", marker=None, color=\"red\", label=\"Modeled\")\n",
    "\n",
    "# Add grid lines and a legend to your plot, and then show to display\n",
    "axis.grid(True)\n",
    "axis.legend(loc=\"best\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Visually Estimating the Slope & Intercept\n",
    "Building linear models is an automated way of doing something we can roughly do \"manually\" \n",
    "with data visualization and a lot of trial-and-error. The visual method is not the most efficient\n",
    "or precise method, but it does illustrate the concepts very well, so let's try it!\n",
    "\n",
    "Given some measured data, your goal is to guess values for slope and intercept, pass them into the\n",
    "model, and adjust your guess until the resulting model fits the data. Use the provided data xd, yd,\n",
    "and the provided function model() to create model predictions. Compare the predictions and data using\n",
    "the provided plot_data_and_model().\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "# Look at the plot data and guess initial trial values\n",
    "trial_slope = 1\n",
    "trial_intercept = 2\n",
    "\n",
    "# input thoses guesses into the model function to compute the model values.\n",
    "xm, ym = model(trial_intercept, trial_slope)\n",
    "\n",
    "# Compare your your model to the data with the plot function\n",
    "fig = plot_data_and_model(xd, yd, xm, ym)\n",
    "plt.show()\n",
    "\n",
    "# Repeat the steps above until your slope and intercept guess makes the model line up with the data.\n",
    "final_slope = 1\n",
    "final_intercept = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-15-94ab47993230>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     10\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     11\u001b[0m \u001b[1;31m# Compute the deviations by subtracting the mean offset\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 12\u001b[1;33m \u001b[0mdx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mx\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     13\u001b[0m \u001b[0mdy\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0my\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     14\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Mean, Deviation, & Standard Deviation\n",
    "The mean describes the center of the data. \n",
    "The standard deviation describes the spread of the data. But to compare two variables, \n",
    "it is convenient to normalize both. In this exercise, you are provided with two arrays of data,\n",
    "which are highly correlated, and you will compute and visualize the normalized deviations of each array.\n",
    "\n",
    "Re-centering / reducing\n",
    "\"\"\"\n",
    "\n",
    "# Compute the deviations by subtracting the mean offset\n",
    "dx = x - np.mean(x)\n",
    "dy = y - np.mean(y)\n",
    "\n",
    "# Normalize the data by dividing the deviations by the standard deviation\n",
    "zx = dx / np.std(x)\n",
    "zy = dy / np.std(y)\n",
    "\n",
    "# Plot comparisons of the raw data and the normalized data\n",
    "fig = plot_cdfs(dx, dy, zx, zy)\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Covariance vs Correlation\n",
    "Covariance is a measure of whether two variables change (\"vary\") together. \n",
    "It is calculated by computing the products, point-by-point, of the deviations seen in the\n",
    "previous exercise, dx[n]*dy[n], and then finding the average of all those products.\n",
    "\n",
    "Correlation is in essence the normalized covariance. In this exercise, you are provided with\n",
    "two arrays of data, which are highly correlated, and you will visualize and compute both the\n",
    "covariance and the correlation.\n",
    "\"\"\"\n",
    "\n",
    "# Compute the covariance from the deviations.\n",
    "dx = x - np.mean(x)\n",
    "dy = y - np.mean(y)\n",
    "covariance = np.mean(dx * dy)\n",
    "print(\"Covariance: \", covariance)\n",
    "\n",
    "# Compute the correlation from the normalized deviations.\n",
    "zx = dx / np.std(x)\n",
    "zy = dy / np.std(y)\n",
    "correlation = np.mean(zx * zy)\n",
    "print(\"Correlation: \", correlation)\n",
    "\n",
    "# Plot the normalized deviations for visual inspection. \n",
    "fig = plot_normalized_deviations(zx, zy)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that you've plotted the product of the normalized deviations, and labeled the plot with the correlation,\n",
    "a single value that is the mean of that product. The product is always positive and the mean is typical of\n",
    "how the two vary together.\n",
    "\"\"\"\n",
    "\n",
    "# Complete the function that will compute correlation.\n",
    "def correlation(x,y):\n",
    "    x_dev = x - np.mean(x)\n",
    "    y_dev = y - np.mean(y)\n",
    "    x_norm = x_dev / np.std(x)\n",
    "    y_norm = y_dev / np.std(y)\n",
    "    return np.mean(x_norm * y_norm)\n",
    "\n",
    "# Compute and store the correlation for each data set in the list.\n",
    "for name, data in data_sets.items():\n",
    "    data['correlation'] = correlation(data['x'], data['y'])\n",
    "    print('data set {} has correlation {:.2f}'.format(name, data['correlation']))\n",
    "\n",
    "# Assign the data set with the best correlation.\n",
    "best_data = data_sets['A']\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Define the general model as a function\n",
    "def model(x, a0=3, a1=2, a2=0):\n",
    "    return a0 + (a1*x) + (a2*x*x)\n",
    "\n",
    "# Generate array x, then predict y values for specific, non-default a0 and a1\n",
    "x = np.linspace(-10, 10, 21)\n",
    "y = model(x)\n",
    "\n",
    "# Plot the results, y versus x\n",
    "fig = plot_prediction(x, y)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Complete the plotting function definition\n",
    "def plot_data_with_model(xd, yd, ym):\n",
    "    fig = plot_data(xd, yd)  # plot measured data\n",
    "    fig.axes[0].plot(xd, ym, color='red')  # over-plot modeled data\n",
    "    plt.show()\n",
    "    return fig\n",
    "\n",
    "# Select new model parameters a0, a1, and generate modeled `ym` from them.\n",
    "a0 = 130\n",
    "a1 = 25\n",
    "ym = model(xd, a0, a1)\n",
    "\n",
    "# Plot the resulting model to see whether it fits the data\n",
    "fig = plot_data_with_model(xd, yd, ym)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Complete the function to convert C to F\n",
    "def convert_scale(temps_C):\n",
    "    (freeze_C, boil_C) = (0, 100)\n",
    "    (freeze_F, boil_F) = (32, 212)\n",
    "    change_in_C = boil_C - freeze_C\n",
    "    change_in_F = boil_F - freeze_F\n",
    "    slope = change_in_F / change_in_C\n",
    "    intercept = freeze_F - freeze_C\n",
    "    temps_F = intercept + (slope * temps_C)\n",
    "    return temps_F\n",
    "\n",
    "# Use the convert function to compute values of F and plot them\n",
    "temps_C = np.linspace(0, 100, 101)\n",
    "temps_F = convert_scale(temps_C)\n",
    "fig = plot_temperatures(temps_C, temps_F)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Generally we might use the average velocity as the slope in our model. \n",
    "But notice that there is some random variation in the instantaneous velocity values when \n",
    "plotted as a time series. The range of values v_max - v_min is one measure of the scale of\n",
    "that variation, and the standard deviation of velocity values is another measure. \n",
    "We see the implications of this variation in a model parameter in the next chapter of this course \n",
    "when discussing inference.\n",
    "\"\"\"\n",
    "\n",
    "# Compute an array of velocities as the slope between each point\n",
    "diff_distances = np.diff(distances)\n",
    "diff_times = np.diff(times)\n",
    "velocities = diff_distances / diff_times\n",
    "\n",
    "# Chracterize the center and spread of the velocities\n",
    "v_avg = np.mean(velocities)\n",
    "v_max = np.max(velocities)\n",
    "v_min = np.min(velocities)\n",
    "v_range = v_max - v_min\n",
    "\n",
    "# Plot the distribution of velocities\n",
    "fig = plot_velocity_timeseries(times[1:], velocities)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Import ols from statsmodels, and fit a model to the data\n",
    "from statsmodels.formula.api import ols\n",
    "model_fit = ols(formula=\"masses ~ volumes\", data=df)\n",
    "model_fit = model_fit.fit()\n",
    "\n",
    "# Extract the model parameter values, and assign them to a0, a1\n",
    "a0 = model_fit.params['Intercept']\n",
    "a1 = model_fit.params['volumes']\n",
    "\n",
    "# Print model parameter values with meaningful names, and compare to summary()\n",
    "print( \"container_mass   = {:0.4f}\".format(a0) )\n",
    "print( \"solution_density = {:0.4f}\".format(a1) )\n",
    "print( model_fit.summary() )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data\n",
    "x_data, y_data = load_data()\n",
    "\n",
    "# Model the data with specified values for parameters a0, a1\n",
    "y_model = model(x_data, a0=150, a1=25)\n",
    "\n",
    "# Compute the RSS value for this parameterization of the model\n",
    "rss = np.sum(np.square(y_model - y_data))\n",
    "print(\"RSS = {}\".format(rss))\n",
    "\n",
    "\"\"\"\n",
    "The value we compute for RSS is not meaningful by itself, but later it becomes meaningful in context\n",
    "when we compare it to other values of RSS computed for other parameterizations of the model.\n",
    "More on that next! Some notes about code style; notice you could have done the RSS calculation\n",
    "in a single line of python code, but writing functions than can be re-used is good practice.\n",
    "Notice also that we could have defined a parameter dictionary dict(a0=150, a1=25) and passed it into the model\n",
    "as model(x, **parameters) which would make it easier to pass around all the parameters together\n",
    "if we needed them for other functions\n",
    "\"\"\"\n",
    "\n",
    "# Complete function to load data, build model, compute RSS, and plot\n",
    "def compute_rss_and_plot_fit(a0, a1):\n",
    "    xd, yd = load_data()\n",
    "    ym = model(xd, a0, a1)\n",
    "    residuals = ym - yd\n",
    "    rss = np.sum(np.square(residuals))\n",
    "    summary = \"Parameters a0={}, a1={} yield RSS={:0.2f}\".format(a0, a1, rss)\n",
    "    fig = plot_data_with_model(xd, yd, ym, summary)\n",
    "    return rss, summary\n",
    "\n",
    "# Chose model parameter values and pass them into RSS function\n",
    "rss, summary = compute_rss_and_plot_fit(a0=150, a1=25)\n",
    "print(summary)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Visualizing the RSS Minima\n",
    "In this exercise you will compute and visualize how RSS varies for different values of model parameters.\n",
    "Start by holding the intercept constant, but vary the slope: and for each slope value,\n",
    "you'll compute the model values, and the resulting RSS. Once you have an array of RSS values,\n",
    "you will determine minimal RSS value, in code, and from that minimum, determine the slope that\n",
    "resulted in that minimal RSS.\n",
    "\"\"\"\n",
    "\n",
    "# Loop over all trial values in a1_array, computing rss for each\n",
    "a1_array = np.linspace(15, 35, 101)\n",
    "for a1_trial in a1_array:\n",
    "    y_model = model(x_data, a0=150, a1=a1_trial)\n",
    "    rss_value = compute_rss(y_data, y_model)\n",
    "    rss_list.append(rss_value)\n",
    "\n",
    "# Find the minimum RSS and the a1 value from whence it came\n",
    "rss_array = np.array(rss_list)\n",
    "best_rss = np.min(rss_array) \n",
    "best_a1 = a1_array[np.where(rss_array==best_rss)]\n",
    "print('The minimum RSS = {}, came from a1 = {}'.format(best_rss, best_a1))\n",
    "\n",
    "# Plot your rss and a1 values to confirm answer\n",
    "fig = plot_rss_vs_a1(a1_array, rss_array)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "The best slope is the one out of an array of slopes than yielded the minimum RSS value out of an array of RSS values.\n",
    "Python tip: notice that we started with rss_list to make it easy to .append() but then later converted to numpy.array()\n",
    "to gain access to all the numpy methods.\n",
    "\"\"\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Least-Squares with `numpy`\n",
    "The formulae below are the result of working through the calculus discussed in the introduction.\n",
    "In this exercise, we'll trust that the calculus correct, and implement these formulae in code using numpy. \n",
    "\"\"\"\n",
    "\n",
    "# prepare the means and deviations of the two variables\n",
    "x_mean = np.mean(x)\n",
    "y_mean = np.mean(y)\n",
    "x_dev = x - x_mean\n",
    "y_dev = y - y_mean\n",
    "\n",
    "# Complete least-squares formulae to find the optimal a0, a1\n",
    "a1 = np.sum(x_dev*y_dev) / np.sum( np.square(x-x_mean) )\n",
    "a0 = y_mean - (a1 * x_mean)\n",
    "\n",
    "# Use the those optimal model parameters a0, a1 to build a model\n",
    "y_model = model(x, a0, a1)\n",
    "\n",
    "# plot to verify that the resulting y_model best fits the data y\n",
    "fig, rss = compute_rss_and_plot_fit(a0, a1)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that the optimal slope a1, according to least-squares, is a ratio of the covariance to the variance.\n",
    "Also, note that the values of the parameters obtained here are NOT exactly the ones used to generate the\n",
    "pre-loaded data (a1=25 and a0=150), but they are close to those. Least-squares does not guarantee zero error;\n",
    "there is no perfect solution, but in this case, least-squares is the best we can do.\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "Optimization with Scipy\n",
    "It is possible to write a numpy implementation of the analytic solution to find the minimal RSS value.\n",
    "But for more complex models, finding analytic formulae is not possible, and so we turn to other methods.\n",
    "\n",
    "In this exercise you will use scipy.optimize to employ a more general approach to solve the same optimization problem.\n",
    "In so doing, you will see additional return values from the method that tell answer us \"how good is best\".\n",
    "Here we will use the same measured data and parameters as seen in the last exercise for ease of comparison\n",
    "of the new scipy approach.\n",
    "\"\"\"\n",
    "\n",
    "# Define a model function needed as input to scipy\n",
    "def model_func(x, a0, a1):\n",
    "    return a0 + (a1*x)\n",
    "\n",
    "# Load the measured data you want to model\n",
    "x_data, y_data  = load_data()\n",
    "\n",
    "# call curve_fit, passing in the model function and data; then unpack the results\n",
    "param_opt, param_cov = optimize.curve_fit(model_func, x_data, y_data)\n",
    "a0 = param_opt[0]  # a0 is the intercept in y = a0 + a1*x\n",
    "a1 = param_opt[1]  # a1 is the slope     in y = a0 + a1*x\n",
    "\n",
    "# test that these parameters result in a model that fits the data\n",
    "fig, rss = compute_rss_and_plot_fit(a0, a1)\n",
    "\n",
    "\"\"\"\n",
    "Notice that we passed the function object itself, model_func into curve_fit, rather than passing in the model data.\n",
    "The model function object was the input, because the optimization wants to know what form in general it's solve for;\n",
    "had we passed in a model_func with more terms like an a2*x**2 term, we would have seen different results for the\n",
    "parameters output\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'ols' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-16-f1a3a442f511>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     11\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     12\u001b[0m \u001b[1;31m# Pass data and `formula` into ols(), use and `.fit()` the model to the data\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 13\u001b[1;33m \u001b[0mmodel_fit\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mols\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mformula\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m\"y_column ~ x_column\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     14\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     15\u001b[0m \u001b[1;31m# Use .predict(df) to get y_model values, then over-plot y_data with y_model\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'ols' is not defined"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Least-Squares with `statsmodels`\n",
    "Several python libraries provide convenient abstracted interfaces so that you need not always be so\n",
    "explicit in handling the machinery of optimization of the model.\n",
    "As an example, in this exercise, you will use the statsmodels library in a more high-level,\n",
    "generalized work-flow for building a model using least-squares optimization (minimization of RSS).\n",
    "To help get you started, we've pre-loaded the data from x_data, y_data = load_data() and stored it\n",
    "in a pandas DataFrame with column names x_column and y_column using \n",
    "df = pd.DataFrame(dict(x_column=x_data, y_column=y_data))\n",
    "\"\"\"\n",
    "\n",
    "# Pass data and `formula` into ols(), use and `.fit()` the model to the data\n",
    "model_fit = ols(formula=\"y_column ~ x_column\", data=df).fit()\n",
    "\n",
    "# Use .predict(df) to get y_model values, then over-plot y_data with y_model\n",
    "y_model = model_fit.predict(df)\n",
    "fig = plot_data_with_model(x_data, y_data, y_model)\n",
    "\n",
    "# Extract the a0, a1 values from model_fit.params\n",
    "a0 = model_fit.params['Intercept']\n",
    "a1 = model_fit.params['x_column']\n",
    "\n",
    "# Visually verify that these parameters a0, a1 give the minimum RSS\n",
    "fig, rss = compute_rss_and_plot_fit(a0, a1)\n",
    "\n",
    "\"\"\"\n",
    "Note that the params container always uses 'Intercept' for the a0 key, but all higher order terms\n",
    "will have keys that match the column name from the pandas DataFrame that you passed into ols().\n",
    "Python style tip: notice that we used 'method chaining', which looks like this: ols().fit().\n",
    "Since the object returned by ols() has a method .fit(), this can be a conveniently shorter way\n",
    "to express calculations in python rather than using multiple lines.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Expected 2D array, got scalar array instead:\narray=50.7.\nReshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-2def854f655b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     23\u001b[0m \u001b[1;31m# Use the fitted model to make a prediction for the found femur\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     24\u001b[0m \u001b[0mfossil_leg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m50.7\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 25\u001b[1;33m \u001b[0mfossil_height\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpredict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfossil_leg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     26\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Predicted fossil height = {:0.2f} cm\"\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfossil_height\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     27\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\base.py\u001b[0m in \u001b[0;36mpredict\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m    211\u001b[0m             \u001b[0mReturns\u001b[0m \u001b[0mpredicted\u001b[0m \u001b[0mvalues\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    212\u001b[0m         \"\"\"\n\u001b[1;32m--> 213\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_decision_function\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    214\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    215\u001b[0m     \u001b[0m_preprocess_data\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mstaticmethod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0m_preprocess_data\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Program Files\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\base.py\u001b[0m in \u001b[0;36m_decision_function\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m    194\u001b[0m         \u001b[0mcheck_is_fitted\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"coef_\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    195\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 196\u001b[1;33m         \u001b[0mX\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcheck_array\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maccept_sparse\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'csr'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'csc'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'coo'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    197\u001b[0m         return safe_sparse_dot(X, self.coef_.T,\n\u001b[0;32m    198\u001b[0m                                dense_output=True) + self.intercept_\n",
      "\u001b[1;32mC:\\Program Files\\Anaconda3\\lib\\site-packages\\sklearn\\utils\\validation.py\u001b[0m in \u001b[0;36mcheck_array\u001b[1;34m(array, accept_sparse, accept_large_sparse, dtype, order, copy, force_all_finite, ensure_2d, allow_nd, ensure_min_samples, ensure_min_features, warn_on_dtype, estimator)\u001b[0m\n\u001b[0;32m    543\u001b[0m                     \u001b[1;34m\"Reshape your data either using array.reshape(-1, 1) if \"\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    544\u001b[0m                     \u001b[1;34m\"your data has a single feature or array.reshape(1, -1) \"\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 545\u001b[1;33m                     \"if it contains a single sample.\".format(array))\n\u001b[0m\u001b[0;32m    546\u001b[0m             \u001b[1;31m# If input is 1D raise error\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    547\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0marray\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: Expected 2D array, got scalar array instead:\narray=50.7.\nReshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample."
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import sys, os\n",
    "from numpy import reshape\n",
    "\n",
    "os.chdir('C:/Users/T0230575/source/repos/DataCamp__/Introduction to Modeling with Python/')\n",
    "df = pd.read_csv('femur_data.csv')\n",
    "\n",
    "heights = df['height']\n",
    "legs = df['length']\n",
    "\n",
    "# import the sklearn class LinearRegression and initialize the model\n",
    "from sklearn.linear_model import LinearRegression\n",
    "model = LinearRegression(fit_intercept=False)\n",
    "\n",
    "# Prepare the measured data arrays and fit the model to them\n",
    "legs = legs.values.reshape(len(heights),1)\n",
    "heights = heights.values.reshape(len(heights),1)\n",
    "model.fit(legs, heights)\n",
    "\n",
    "# Use the fitted model to make a prediction for the found femur\n",
    "fossil_leg = 50.7\n",
    "fossil_height = model.predict(fossil_leg)\n",
    "print(\"Predicted fossil height = {:0.2f} cm\".format(fossil_height[0,0]))\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that we used the pre-loaded data to fit or \"train\" the model, and then applied that model\n",
    "to make a prediction about newly collected data that was not part of the data used to fit the model.\n",
    "Also notice that model.predict() returns the answer as an array of shape = (1,1), so we had to index\n",
    "into it with the [0,0] syntax when printing. This is an artifact of our overly simplified use of sklearn here:\n",
    "the details of this are beyond the scope of the current course, but relate to the number of samples and features\n",
    "that one might use in a more sophisticated, generalized model.\n",
    "\"\"\"\n",
    "\n",
    "# Import LinearRegression class, build a model, fit to the data\n",
    "from sklearn.linear_model import LinearRegression\n",
    "model = LinearRegression(fit_intercept=True)\n",
    "model.fit(years, levels)\n",
    "\n",
    "# Use model to make a prediction for one year, 2100\n",
    "future_year = 2100\n",
    "future_level = model.predict(future_year)\n",
    "print(\"Prediction: year = {}, level = {:.02f}\".format(future_year, future_level[0,0]))\n",
    "\n",
    "# Use model to predict for many years, and over-plot with measured data\n",
    "years_forecast = np.linspace(1970, 2100, 131).reshape(-1, 1)\n",
    "levels_forecast = model.predict(years_forecast)\n",
    "help(plot_data_and_forecast)\n",
    "fig = plot_data_and_forecast(years, levels, years_forecast, levels_forecast)\n",
    "\n",
    "\"\"\"\n",
    "Note that with scikit-learn, although we could extract a0 = model.intercept_[0] and a1 = model.coef_[0,0],\n",
    "we do not need to do that in order to make predictions, we just call model.predict(). With more complex models,\n",
    "these parameters may not have easy physical interpretations. Notice also that although our model is linear,\n",
    "the actual data appears to have an up-turn that might be better modeled by adding a quadratic or even exponential\n",
    "term to our model. The linear model forecast may be underestimating the rate of increase in sea level.\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\n",
    "# Fit the model, based on the form of the formula\n",
    "model_fit = ols(formula=\"velocities ~ distances\", data=df).fit()\n",
    "\n",
    "# Extract the model parameters and associated \"errors\" or uncertainties\n",
    "a0 = model_fit.params['Intercept']\n",
    "a1 = model_fit.params['distances']\n",
    "e0 = model_fit.bse['Intercept']\n",
    "e1 = model_fit.bse['distances']\n",
    "\n",
    "# Print the results\n",
    "print('For slope a1={:.02f}, the uncertainty in a1 is {:.02f}'.format(a1, e1))\n",
    "print('For intercept a0={:.02f}, the uncertainty in a0 is {:.02f}'.format(a0, e0))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# build and fit a model to the df_monthly data\n",
    "model_fit = ols('Close ~ DayCount', data=df_monthly).fit()\n",
    "\n",
    "# Use the model FIT to the MONTHLY data to make a predictions for both monthly and daily data\n",
    "df_monthly['Model'] = model_fit.predict(df_monthly.DayCount)\n",
    "df_daily['Model'] = model_fit.predict(df_daily.DayCount)\n",
    "\n",
    "# Plot the monthly and daily data and model, compare the RSS values seen on the figures\n",
    "fig_monthly = plot_model_with_data(df_monthly)\n",
    "fig_daily = plot_model_with_data(df_daily)\n",
    "\n",
    "\"\"\"\n",
    "Notice the monthly data looked linear, but the daily data clearly has additional, \n",
    "nonlinear trends. Under-sampled data often misses real-world features in the data on smaller time\n",
    "or spatial scales. Using the model from the under-sampled data to make interpolations to the daily\n",
    "data can result is large residuals. Notice that the RSS value for the daily plot is more than 30 times\n",
    "worse than the monthly plot\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Extrapolation: Going Over the Edge\n",
    "In this exercise, we consider the perils of extrapolation. Shown here is the profile of a hiking trail on a mountain.\n",
    "One portion of the trail, marked in black, looks linear, and was used to build a model. But we see that the best fit\n",
    "line, shown in red, does not fit outside the original \"domain\", as it extends into this new outside data, marked in blue.\n",
    "If we want use the model to make predictions for the altitude, but still be accurate to within some tolerance, what are\n",
    "the smallest and largest values of independent variable x that we can allow ourselves to apply the model to?\"\n",
    "\"\"\"\n",
    "\n",
    "# Compute the residuals, \"data - model\", and determine where [residuals < tolerance]\n",
    "residuals = np.abs(y_data - y_model)\n",
    "tolerance = 100\n",
    "x_good = x_data[residuals < tolerance]\n",
    "\n",
    "# Find the min and max of the \"good\" values, and plot y_data, y_model, and the tolerance range\n",
    "print('Minimum good x value = {}'.format(np.min(x_good)))\n",
    "print('Maximum good x value = {}'.format(np.max(x_good)))\n",
    "fig = plot_data_model_tolerance(x_data, y_data, y_model, tolerance)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice the range of good values, which extends a little out into the new data, is marked in green on the plot.\n",
    "By comparing the residuals to a tolerance threshold, we can quantify how far out out extrapolation can go before\n",
    "the difference between model and data gets too large.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'model_fit_and_predict' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-4-b72873f18486>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     15\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     16\u001b[0m \u001b[1;31m# Build the model and compute the residuals \"model - data\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 17\u001b[1;33m \u001b[0my_model\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmodel_fit_and_predict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx_data\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my_data\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     18\u001b[0m \u001b[0mresiduals\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0my_model\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0my_data\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     19\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'model_fit_and_predict' is not defined"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels\n",
    "import os, sys\n",
    "\n",
    "os.chdir('C:/Users/T0230575/source/repos/DataCamp__/Python/Introduction to Linear Modeling with Python/')\n",
    "df = pd.read_csv('hiking_data.csv')\n",
    "\n",
    "\"\"\"\n",
    "Quantifying Goodness of fit: \n",
    "RSS (Building)\n",
    "RMSE (Evaluating) - most common metric\n",
    "R-squared (Evaluating)\n",
    "\"\"\"\n",
    "\n",
    "# Build the model and compute the residuals \"model - data\"\n",
    "y_model = model_fit_and_predict(x_data, y_data)\n",
    "residuals = y_model - y_data\n",
    "\n",
    "# Compute the RSS, MSE, and RMSE and print the results\n",
    "RSS = np.sum(np.square(residuals))\n",
    "MSE = RSS/len(residuals)\n",
    "RMSE = np.sqrt(MSE)\n",
    "print('RMSE = {:0.2f}, MSE = {:0.2f}, RSS = {:0.2f}'.format(RMSE, MSE, RSS))\n",
    "\n",
    "\"\"\"\n",
    "Notice that instead of computing RSS and normalizing with division by len(residuals) to get the MSE, \n",
    "you could have just applied np.mean(np.square()) to the residuals. Another useful point to help you remember;\n",
    "you can think of the MSE like a variance, but instead of differencing the data from its mean, you difference\n",
    "the data and the model. Similarly, think of RMSE as a standard deviation.\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "R-Squared\n",
    "In this exercise you'll compute another measure of goodness, R-squared. R-squared is the ratio of the variance\n",
    "of the residuals divided by the variance of the data we are modeling, and in so doing, is a measure of how much\n",
    "of the variance in your data is \"explained\" by your model, as expressed in the spread of the residuals.\n",
    "Here we have pre-loaded the data x_data,y_data and the model predictions y_model for the best fit model;\n",
    "you're goal is to compute the R-squared measure to quantify how much this linear model accounts for variation\n",
    "in the data.\n",
    "\"\"\"\n",
    "\n",
    "# Compute the residuals and the deviations\n",
    "residuals = y_model - y_data\n",
    "deviations = np.mean(y_data) - y_data\n",
    "\n",
    "# Compute the variance of the residuals and deviations\n",
    "var_residuals = np.mean(np.square(residuals))\n",
    "var_deviations = np.mean(np.square(deviations))\n",
    "\n",
    "# Compute r_squared as 1 - the ratio of RSS/Variance\n",
    "r_squared = 1 - (var_residuals / var_deviations)\n",
    "print('R-squared is {:0.2f}'.format(r_squared))\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that R-squared varies from 0 to 1, where a value of 1 means that the model and the data are perfectly\n",
    "correlated and all variation in the data is predicted by the model. A value of zero would mean none of the\n",
    "variation in the data is predicted by the model. Here, the data points are close to the line, so R-squared is\n",
    "closer to 1.0\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Uncertainty in predictions\n",
    "\n",
    "Variation Around the Trend\n",
    "The data need not be perfectly linear, and there may be some random variation or \"spread\" in the measurements,\n",
    "and that does translate into variation of the model parameters. This variation is in the parameter is quantified\n",
    "by \"standard error\", and interpreted as \"uncertainty\" in the estimate of the model parameter.\n",
    "In this exercise, you will use ols from statsmodels to build a model and extract the standard error for each parameter\n",
    "of that model.\n",
    "\"\"\"\n",
    "\n",
    "# Store x_data and y_data, as times and distances, in df, and use ols() to fit a model to it.\n",
    "df = pd.DataFrame(dict(times=x_data, distances=y_data))\n",
    "model_fit = ols(formula=\"distances ~ times\", data=df).fit()\n",
    "\n",
    "# Extact the model parameters and their uncertainties\n",
    "a0 = model_fit.params['Intercept']\n",
    "e0 = model_fit.bse['Intercept']\n",
    "a1 = model_fit.params['times']\n",
    "e1 = model_fit.bse['times']\n",
    "\n",
    "# Print the results with more meaningful names\n",
    "print('Estimate    of the intercept = {:0.2f}'.format(a0))\n",
    "print('Uncertainty of the intercept = {:0.2f}'.format(e0))\n",
    "print('Estimate    of the slope = {:0.2f}'.format(a1))\n",
    "print('Uncertainty of the slope = {:0.2f}'.format(e1))\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "The size of the parameters standard error only makes sense in comparison to the parameter value itself.\n",
    "In fact the units are the same! So a1 and e1 both have units of velocity (meters/second), and a0 and e0\n",
    "both have units of distance (meters).\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "Variation in Two Parts\n",
    "Given two data sets of distance-versus-time data, one with very small velocity and one with large velocity.\n",
    "Notice that both may have the same standard error of slope, but different R-squared for the model overall,\n",
    "depending on the size of the slope (\"effect size\") as compared to the standard error (\"uncertainty\").\n",
    "If we plot both data sets as scatter plots on the same axes, the contrast is clear. Variation due to the slope\n",
    "is different than variation due to the random scatter about the trend line. In this exercise, your goal is to\n",
    "compute the standard error and R-squared for two data sets and compare.\n",
    "\"\"\"\n",
    "\n",
    "# Build and fit two models, for columns distances1 and distances2 in df\n",
    "model_1 = ols(formula=\"distances1 ~ times\", data=df).fit()\n",
    "model_2 = ols(formula=\"distances2 ~ times\", data=df).fit()\n",
    "\n",
    "# Extract R-squared for each model, and the standard error for each slope\n",
    "se_1 = model_1.bse['times']\n",
    "se_2 = model_2.bse['times']\n",
    "rsquared_1 = model_1.rsquared\n",
    "rsquared_2 = model_2.rsquared\n",
    "\n",
    "# Print the results\n",
    "print('Model 1: SE = {:0.3f}, R-squared = {:0.3f}'.format(se_1, rsquared_1))\n",
    "print('Model 2: SE = {:0.3f}, R-squared = {:0.3f}'.format(se_2, rsquared_2))\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that the standard error is the same for both models, but the r-squared changes. The uncertainty in the\n",
    "estimates of the model parameters is indepedent from R-squred because that uncertainty is being driven not by\n",
    "the linear trend, but by the inherent randomness in the data. This serves as a transition into looking at statistical\n",
    "inference in linear models.\n",
    "\"\"\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Population mean 0.5, stdev 0.28\n",
      "    Sample mean 0.560, stdev 0.296\n",
      "Means:  center=  0.51, spread=  0.01\n",
      "Stdevs: center=  0.28, spread=  0.00\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADXJJREFUeJzt3X+sJeVdx/H3p6yoIKQQTrUC17ttkISSJuiNqTZWU4pdtUKN/QNizaKYmyZtqVrTQqpZof8YNdYmbdKsLUIihRisEZtaWatImgBxlx/lxxahgHQpytJtUo21dOvXP+7Z9HJ7d885M3N/Pff9Sk7uzJxnznyfTO4nT2bOPCdVhSRp63vZRhcgSRqGgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqxI71PNhZZ51V8/Pz63lISdryDhw48EJVjSa1W9dAn5+fZ//+/et5SEna8pL8+zTtvOQiSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNaDbQk42uQJLWV7OBLknbjYEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJasTEQE9yQ5Lnkzy8ynu/m6SSnLU25UmSpjXNCP1GYNfKjUnOBS4Bnhm4JklSBxMDvaruAo6s8taHgPcBNXRRkqTZdbqGnuRS4NmqenDgeiRJHe2YdYckpwAfAH5uyvaLwCLA3NzcrIfrxJkWJW1HXUborwZ2Ag8meRo4B7gvyQ+t1riq9lbVQlUtjEaj7pVKkk5o5hF6VT0EvOLY+jjUF6rqhQHrkiTNaJqvLd4C3A2cn+RQkqvWvixJ0qwmjtCr6ooJ788PVo0kqTOfFJWkRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpCnlupDrNu/sfwa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1YprfFL0hyfNJHl627Y+TfDHJF5L8TZKXr22ZkqRJphmh3wjsWrFtH3BhVb0W+Dfg2oHrkiTNaGKgV9VdwJEV2+6oqqPj1XuAc9agNknSDIa4hv4bwN8P8DmSpB56BXqSDwBHgZtP0GYxyf4k+w8fPtzncJKkE+gc6El2A28BfrWq6njtqmpvVS1U1cJoNOp6OEnSBDu67JRkF/B+4Geq6n+GLUmS1MU0X1u8BbgbOD/JoSRXAR8BTgP2JXkgycfWuE5J0gQTR+hVdcUqmz+xBrVIknrwSVFJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEdP8SPQNSZ5P8vCybWcm2Zfk8fHfM9a2TEnSJNOM0G8Edq3Ydg3wuao6D/jceF2StIEmBnpV3QUcWbH5MuCm8fJNwFsHrkuSNKOu19B/sKqeAxj/fcVwJUmSutix1gdIsggsAszNza3xsY6/vWpNDy1pC8t13x0etWe20Fj+GbPuO5SuI/T/TPJKgPHf54/XsKr2VtVCVS2MRqOOh5MkTdI10G8Hdo+XdwN/O0w5kqSupvna4i3A3cD5SQ4luQr4Q+CSJI8Dl4zXJUkbaOI19Kq64jhvXTxwLZKkHnxSVJIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS2parsuqsym2yECXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RG9Ar0JL+d5JEkDye5Jcn3DVWYJGk2nQM9ydnA1cBCVV0InARcPlRhkqTZ9L3ksgP4/iQ7gFOAr/QvSZLUxY6uO1bVs0n+BHgG+AZwR1XdsbJdkkVgEWBubq7r4U4ogarVt0tSF8sn9Ko9qwTMJtTnkssZwGXATuCHgVOTvH1lu6raW1ULVbUwGo26VypJOqE+l1zeBDxVVYer6lvAp4CfGqYsSdKs+gT6M8DrkpySJMDFwMFhypIkzapzoFfVvcBtwH3AQ+PP2jtQXZKkGXW+KQpQVXuAPQPVIknqwSdFJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1olegJ3l5ktuSfDHJwSQ/OVRhkqTZ9PpNUeDDwGer6m1JTgZOGaAmSVIHnQM9yenAG4ArAarqReDFYcqSJM2qzyWXVwGHgb9Icn+Sjyc5daC6JEkz6nPJZQfwY8C7q+reJB8GrgF+f3mjJIvAIsDc3FyPw71UMthHSdqmcl23IFm+X+2pocrprc8I/RBwqKruHa/fxlLAv0RV7a2qhapaGI1GPQ4nSTqRzoFeVf8BfDnJ+eNNFwOPDlKVJGlmfb/l8m7g5vE3XJ4Efr1/SZKkLnoFelU9ACwMVIskqQefFJWkRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEb0ffR/03D2RUkboeuMjWvBEbokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDWid6AnOSnJ/Uk+PURBkqRuhhihvwc4OMDnSJJ66BXoSc4BfhH4+DDlSJK66jtC/zPgfcD/DVCLJKmHzrMtJnkL8HxVHUjysydotwgsAszNzXU93CASqNrQEiRtkD6zIs6677H2tWd9A6fPCP31wKVJngZuBd6Y5C9XNqqqvVW1UFULo9Gox+EkSSfSOdCr6tqqOqeq5oHLgX+qqrcPVpkkaSZ+D12SGjHILxZV1Z3AnUN8liSpG0foktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYM8mDRekn3uXV67Stpa+kzEddW5ghdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiM6B3qSc5P8c5KDSR5J8p4hC5MkzabPXC5HgfdW1X1JTgMOJNlXVY8OVJskaQadR+hV9VxV3Tde/i/gIHD2UIVJkmYzyGyLSeaBi4B7V3lvEVgEmJubG+Jwg0igaqOrkNSy5bM+1p61D5zeN0WT/ADw18BvVdXXV75fVXuraqGqFkajUd/DSZKOo1egJ/kelsL85qr61DAlSZK66PMtlwCfAA5W1Z8OV5IkqYs+I/TXA78GvDHJA+PXLwxUlyRpRp1vilbV54Ht+TtPkrQJ+aSoJDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqRGDzLa4HrJGjzAd+9wqZ2CUNqtpZy1c3m47coQuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmN6BXoSXYleSzJE0muGaooSdLsOgd6kpOAjwI/D1wAXJHkgqEKkyTNps8I/SeAJ6rqyap6EbgVuGyYsiRJs+oT6GcDX162fmi8TZK0AfrMtrjatGbfNQ1akkVgcbz630ke63HM3pbP2rja8hSzOp4FvDBsVZvWduor2N8tIX/QaUbFDe9rx7qP+ZFpGvUJ9EPAucvWzwG+srJRVe0F9vY4zqaSZH9VLWx0HethO/UV7G/Ltktf+1xy+VfgvCQ7k5wMXA7cPkxZkqRZdR6hV9XRJO8C/gE4Cbihqh4ZrDJJ0kx6/WJRVX0G+MxAtWwVzVw+msJ26ivY35Zti76m/M01SWqCj/5LUiMM9LFJ0xgkeUeSh5I8kOTzy5+KTXLteL/Hkrx5fSvvpmt/k8wn+cZ4+wNJPrb+1c9u2mkqkrwtSSVZWLZtS53frn1t9dwmuTLJ4WX9+s1l7+1O8vj4tXt9K18DVbXtXyzd1P0S8CrgZOBB4IIVbU5ftnwp8Nnx8gXj9t8L7Bx/zkkb3ac17O888PBG92Ho/o7bnQbcBdwDLGzF89uzr02eW+BK4COr7Hsm8OT47xnj5TM2uk99Xo7Ql0ycxqCqvr5s9VS+8xDVZcCtVfXNqnoKeGL8eZtZn/5uRdNOU/FB4I+A/122baud3z593Yr6TEHyZmBfVR2pqq8B+4Bda1TnujDQl0w1jUGSdyb5Ekv/CFfPsu8m06e/ADuT3J/kX5L89NqWOoiJ/U1yEXBuVX161n03mT59hQbP7divJPlCktuSHHsgcqud24kM9CVTTWNQVR+tqlcD7wd+b5Z9N5k+/X0OmKuqi4DfAT6Z5PQ1q3QYJ+xvkpcBHwLeO+u+m1CfvjZ3bsf+DpivqtcC/wjcNMO+W4qBvmSqaQyWuRV4a8d9N4PO/R1fevjqePkAS9cvf3SN6hzKpP6eBlwI3JnkaeB1wO3jm4Vb7fx27muj55aq+mpVfXO8+ufAj0+775az0RfxN8OLpQesnmTpptexGyuvWdHmvGXLvwTsHy+/hpfeNHuSTXzTbID+jo71j6UbUc8CZ250n/r2d0X7O/nOjcItdX579rXJcwu8ctnyLwP3jJfPBJ5i6YboGePlTd3fSa9eT4q2oo4zjUGS61kKstuBdyV5E/At4GvA7vG+jyT5K+BR4Cjwzqr69oZ0ZEp9+gu8Abg+yVHg28A7qurI+vdielP293j7bqnz26evtHtur05yKUvn7whL33qhqo4k+SBL81IBXL/Z+zuJT4pKUiO8hi5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqxP8DhoDF0Wur4E0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "ename": "ValueError",
     "evalue": "too many values to unpack (expected 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-19-c248fe5eced0>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     99\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    100\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 101\u001b[1;33m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmeans\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbins\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mmean_bins\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'green'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    102\u001b[0m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstdevs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbins\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mstd_bins\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'red'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    103\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: too many values to unpack (expected 2)"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADkNJREFUeJzt23+s3Xddx/HnizXMoLJf3S/XNV2yRizyQ3ZSNCIgY6MzQifMZGikidMmhkUGwVCCOrfxx4aYGeLENBtZ2R9suECoIdiUIWKEzZ3ChBUYLUPT6xZWbJ1ZiMzq2z/Ot+Z+bs7dvb3ne+/hXp6P5OSc7+d87j3vb7b0ec/33JuqQpKkk5437QEkST9cDIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDXWTXuApVi/fn1t2rRp2mNI0qqxfv169u3bt6+qti20d1WGYdOmTQyHw2mPIUmrSpL1i9nnpSRJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElq9BKGJNuSPJbkcJJdY54/Pcl93fMPJdk05/mNSZ5J8u4+5pEkLd3EYUhyGnAHcBWwBXhrki1ztl0HHK+qS4HbgdvmPH878JlJZ5EkTa6PdwxbgcNV9XhVPQvcC2yfs2c7sKd7fD9weZIAJLkaeBw42MMskqQJ9RGGi4Ajs45nurWxe6rqBPA0cE6SHwfeA9zUwxySpB70EYaMWatF7rkJuL2qnlnwRZKdSYZJhkePHl3CmJKkxVjXw/eYAS6edbwBeGKePTNJ1gFnAMeAVwLXJPkAcCbwv0n+q6r+Yu6LVNVuYDfAYDCYGx5JUk/6CMPDwOYklwD/BlwL/MacPXuBHcCXgGuAz1VVAb90ckOSPwGeGRcFSdLKmTgMVXUiyfXAPuA04CNVdTDJzcCwqvYCdwH3JDnM6J3CtZO+riRpeWT0g/vqMhgMajgcTnsMSVpVkhyoqsFC+/zLZ0lSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEmNXsKQZFuSx5IcTrJrzPOnJ7mve/6hJJu69SuSHEjyte7+dX3MI0lauonDkOQ04A7gKmAL8NYkW+Zsuw44XlWXArcDt3Xr3wPeWFUvAXYA90w6jyRpMn28Y9gKHK6qx6vqWeBeYPucPduBPd3j+4HLk6SqvlJVT3TrB4EfS3J6DzNJkpaojzBcBByZdTzTrY3dU1UngKeBc+bseQvwlar6QQ8zSZKWaF0P3yNj1upU9iR5MaPLS1fO+yLJTmAnwMaNG099SknSovTxjmEGuHjW8Qbgifn2JFkHnAEc6443AJ8E3lZV357vRapqd1UNqmpw7rnn9jC2JGmcPsLwMLA5ySVJng9cC+yds2cvow+XAa4BPldVleRM4NPAe6vqH3uYRZI0oYnD0H1mcD2wD/gG8PGqOpjk5iRv6rbdBZyT5DDwLuDkr7ReD1wK/FGSR7rbeZPOJElaulTN/Tjgh99gMKjhcDjtMSRpVUlyoKoGC+3zL58lSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqRGL2FIsi3JY0kOJ9k15vnTk9zXPf9Qkk2znntvt/5Ykjf0MY8kaekmDkOS04A7gKuALcBbk2yZs+064HhVXQrcDtzWfe0W4FrgxcA24C+77ydJmpI+3jFsBQ5X1eNV9SxwL7B9zp7twJ7u8f3A5UnSrd9bVT+oqu8Ah7vvJ0makj7CcBFwZNbxTLc2dk9VnQCeBs5Z5NdKklZQH2HImLVa5J7FfO3oGyQ7kwyTDI8ePXqKI0qSFquPMMwAF8863gA8Md+eJOuAM4Bji/xaAKpqd1UNqmpw7rnn9jC2JGmcPsLwMLA5ySVJns/ow+S9c/bsBXZ0j68BPldV1a1f2/3W0iXAZuCfephJkrRE6yb9BlV1Isn1wD7gNOAjVXUwyc3AsKr2AncB9yQ5zOidwrXd1x5M8nHg68AJ4O1V9T+TziRJWrqMfnBfXQaDQQ2Hw2mPIUmrSpIDVTVYaJ9/+SxJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpMVEYkpydZH+SQ939WfPs29HtOZRkR7f2giSfTvLNJAeT3DrJLJKkfkz6jmEX8EBVbQYe6I4bSc4GbgReCWwFbpwVkA9W1YuAnwN+MclVE84jSZrQpGHYDuzpHu8Brh6z5w3A/qo6VlXHgf3Atqr6flX9HUBVPQt8Gdgw4TySpAlNGobzq+pJgO7+vDF7LgKOzDqe6db+X5IzgTcyetchSZqidQttSPJZ4IIxT71vka+RMWs16/uvAz4GfKiqHn+OOXYCOwE2bty4yJeWJJ2qBcNQVa+f77kk301yYVU9meRC4Kkx22aA18463gB8ftbxbuBQVf35AnPs7vYyGAzqufZKkpZu0ktJe4Ed3eMdwKfG7NkHXJnkrO5D5yu7NZK8HzgDuGHCOSRJPZk0DLcCVyQ5BFzRHZNkkOROgKo6BtwCPNzdbq6qY0k2MLoctQX4cpJHkvzOhPNIkiaUqtV3VWYwGNRwOJz2GJK0qiQ5UFWDhfb5l8+SpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVJjojAkOTvJ/iSHuvuz5tm3o9tzKMmOMc/vTfLoJLNIkvox6TuGXcADVbUZeKA7biQ5G7gReCWwFbhxdkCSvBl4ZsI5JEk9mTQM24E93eM9wNVj9rwB2F9Vx6rqOLAf2AaQ5CeAdwHvn3AOSVJPJg3D+VX1JEB3f96YPRcBR2Ydz3RrALcAfwZ8f8I5JEk9WbfQhiSfBS4Y89T7FvkaGbNWSV4OXFpV70yyaRFz7AR2AmzcuHGRLy1JOlULhqGqXj/fc0m+m+TCqnoyyYXAU2O2zQCvnXW8Afg88AvAZUn+pZvjvCSfr6rXMkZV7QZ2AwwGg1pobknS0kx6KWkvcPK3jHYAnxqzZx9wZZKzug+drwT2VdWHq+qnqmoT8CrgW/NFQZK0ciYNw63AFUkOAVd0xyQZJLkToKqOMfos4eHudnO3Jkn6IZSq1XdVZjAY1HA4nPYYkrSqJDlQVYOF9vmXz5KkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkRqpq2jOcsiRHgX+d9hynaD3wvWkPscI85x8NnvPq8D2Aqtq20MZVGYbVKMmwqgbTnmMlec4/GjzntcdLSZKkhmGQJDUMw8rZPe0BpsBz/tHgOa8xfsYgSWr4jkGS1DAMyyTJO5I8muRgkhu6tZcneTDJI0mGSbZOe86+zHO+L0vypSRfS/I3SV447TknleQjSZ5K8uistbOT7E9yqLs/q1tPkg8lOZzkq0leMb3Jl+4Uz/lF3X/zHyR59/SmXrpTPN/f7P7bfjXJF5O8bHqT98cwLIMkPwv8LrAVeBnwq0k2Ax8AbqqqlwN/3B2ves9xvncCu6rqJcAngT+Y3pS9uRuY+3vgu4AHqmoz8EB3DHAVsLm77QQ+vEIz9u1uFn/Ox4DfBz64YtP1724Wf77fAV5TVS8FbmGNfPZgGJbHzwAPVtX3q+oE8PfArwEFnPyp+QzgiSnN17f5zvengS90e/YDb5nSfL2pqi8w+sdvtu3Anu7xHuDqWesfrZEHgTOTXLgyk/bnVM65qp6qqoeB/165Cft1iuf7xao63q0/CGxYkSGXmWFYHo8Cr05yTpIXAL8CXAzcAPxpkiOMfqJ67xRn7NN85/so8KZuz693a2vR+VX1JEB3f163fhFwZNa+mW5tLZjvnNeqxZzvdcBnVnSqZWIYlkFVfQO4jdFPyX8L/DNwAvg94J1VdTHwTuCuqQ3Zo+c4398G3p7kAPCTwLNTG3I6MmbNXwNcg5L8MqMwvGfas/TBMCyTqrqrql5RVa9m9Lb0ELAD+ES35a8ZXZNfE8adb1V9s6qurKrLgI8B357ulMvmuycvEXX3T3XrM7Tvkjawdi4fznfOa9W855vkpYw+T9teVf8+pfl6ZRiWSZLzuvuNwJsZ/cP4BPCabsvrGMViTRh3vrPWngf8IfBX05twWe1lFH26+0/NWn9b99tJPw88ffJyxBow3zmvVWPPt/v//RPAb1XVt6Y0W/+qytsy3IB/AL7O6LLK5d3aq4AD3dpDwGXTnnOZz/cdwLe62610f1C5mm+MAv8kow9XZxhdPjiH0W+qHOruz+72BriD0TulrwGDac+/Aud8QbfnP4H/6B6/cNrnsIzneydwHHikuw2nPX8fN//yWZLU8FKSJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1/g+Zh5BB9uaTKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Inferential Statistics Concepts\n",
    "Sampling with Python\n",
    "\"\"\"\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels\n",
    "import os, sys\n",
    "\n",
    "os.chdir('C:/Users/T0230575/source/repos/DataCamp__/Python/Introduction to Linear Modeling with Python/')\n",
    "population = pd.read_csv('population.csv', header=None)\n",
    "\n",
    "population = np.random.random_sample(size=300)\n",
    "population = population.astype('float')\n",
    "\n",
    "\"\"\"\n",
    "Sample Statistics versus Population\n",
    "In this exercise you will work with a preloaded population. You will construct a sample by drawing points at random from the population. You will compute the mean standard deviation of the sample taken from that population to test whether the sample is representative of the population. Your goal is to see where the sample statistics are the same or very close to the population statistics.\n",
    "\"\"\"\n",
    "\n",
    "# Compute the population statistics\n",
    "print(\"Population mean {:.1f}, stdev {:.2f}\".format( population.mean(), population.std() ))\n",
    "\n",
    "# Set random seed for reproducibility\n",
    "np.random.seed(42)\n",
    "\n",
    "# Construct a sample by randomly sampling 31 points from the population\n",
    "sample = np.random.choice(population, size=31)\n",
    "\n",
    "# Compare sample statistics to the population statistics\n",
    "print(\"    Sample mean {:.3f}, stdev {:.3f}\".format( sample.mean(), sample.std() ))\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Variation in Sample Statistics\n",
    "If we create one sample of size=1000 by drawing that many points from a population. Then compute a sample statistic,\n",
    "such as the mean, a single value that summarizes the sample itself.\n",
    "If you repeat that sampling process num_samples=100 times, you get 100 samples. Computing the sample statistic,\n",
    "like the mean, for each of the different samples, will result in a distribution of values of the mean. The goal\n",
    "then is to compute the mean of the means and standard deviation of the means.\n",
    "\"\"\"\n",
    "\n",
    "num_samples = 100\n",
    "num_pts = 1000\n",
    "\n",
    "# Initialize two arrays of zeros to be used as containers\n",
    "means = np.zeros(num_samples)\n",
    "stdevs = np.zeros(num_samples)\n",
    "\n",
    "# For each iteration, compute and store the sample mean and sample stdev\n",
    "for ns in range(num_samples):\n",
    "    sample = np.random.choice(population, num_pts)\n",
    "    means[ns] = sample.mean()\n",
    "    stdevs[ns] = sample.std()\n",
    "\n",
    "# Compute and print the mean() and std() for the sample statistic distributions\n",
    "print(\"Means:  center={:>6.2f}, spread={:>6.2f}\".format(means.mean(), means.std()))\n",
    "print(\"Stdevs: center={:>6.2f}, spread={:>6.2f}\".format(stdevs.mean(), stdevs.std()))\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "If we only took one sample, instead of 100, there could be only a single mean and the standard deviation\n",
    "of that single value is zero. But each sample is different because of the randomness of the draws.\n",
    "The mean of the means is our estimate for the population mean, the stdev of the means is our measure\n",
    "of the uncertainty in our estimate of the population mean. This is the same concept as the standard error\n",
    "of the slope seen in linear regression.\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "Visualizing Variation of a Statistic\n",
    "Previously, you have computed the variation of sample statistics. Now you'll visualize that variation.\n",
    "We'll start with a preloaded population and a predefined function get_sample_statistics() to draw the samples,\n",
    "and return the sample statistics arrays.\n",
    "Here we will use a predefined plot_hist() function that wraps the matplotlib method axis.hist(), which both bins\n",
    "and plots the array passed in. In this way you can see how the sample statistics have a distribution of values,\n",
    "not just a single value.\n",
    "\"\"\"\n",
    "\n",
    "def get_sample_statistics(population, num_samples, num_pts):\n",
    "    for i in range(num_samples):\n",
    "        sample = np.random.choice(population, num_pts)\n",
    "        means[i] = sample.mean()\n",
    "        stdevs[i] = sample.std()    \n",
    "    return means, stdevs\n",
    "\n",
    "# Generate sample distribution and associated statistics\n",
    "means, stdevs = get_sample_statistics(population, num_samples=100, num_pts=1000)\n",
    "\n",
    "# Define the binning for the histograms\n",
    "mean_bins = np.linspace(97.5, 102.5, 51)\n",
    "std_bins = np.linspace(7.5, 12.5, 51)\n",
    "\n",
    "# Plot the distribution of means, and the distribution of stdevs\n",
    "plt.hist(means, bins = 20, color='green')\n",
    "plt.hist(stdevs, bins = 20, color='blue')\n",
    "plt.show()\n",
    "\n",
    "fig, ax = plt.hist(means, bins=mean_bins, color='green')\n",
    "fig, ax = plt.hist(stdevs, bins=std_bins, color='red')\n",
    "\n",
    "\"\"\"\n",
    "Notice you have to page through the plots to see both. Can you see the center and spread in the title\n",
    "and the plots? If you have not before, compute those values using e.g. means.mean() and means.std()\n",
    "to see that they match the figure.\n",
    "\"\"\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\"\"\"\n",
    "Estimation of Population Parameters\n",
    "Imagine a constellation (\"population\") of satellites orbiting for a full year, and the distance traveled in each\n",
    "hour is measured in kilometers. There is variation in the distances measured from hour-to-hour, due to unknown\n",
    "complications of orbital dynamics. Assume we cannot measure all the data for the year, but we wish to build a\n",
    "population model for the variations in orbital distance per hour (speed) based on a sample of measurements.\n",
    "\n",
    "In this exercise, you will assume that the population of hourly distances are best modeled by a gaussian, \n",
    "and further assume that the parameters of that population model can be estimated from the sample statistics.\n",
    "Start with the preloaded sample_distances that was taken from a population of cars.\n",
    "\"\"\"\n",
    "\n",
    "os.chdir('C:/Users/T0230575/source/repos/DataCamp__/Python/Introduction to Linear Modeling with Python/')\n",
    "sample_distances = pd.read_csv('distances.csv', header=None)\n",
    "print(sample_distances.head())\n",
    "\n",
    "# Compute the mean and standard deviation of the sample_distances\n",
    "sample_mean = np.mean(sample_distances)\n",
    "sample_stdev = np.std(sample_distances)\n",
    "\n",
    "def gaussian_model(x, mu, sigma):\n",
    "    coeff_part = 1/(np.sqrt(2 * np.pi * sigma**2))\n",
    "    exp_part = np.exp( - (x - mu)**2 / (2 * sigma**2) )\n",
    "    return coeff_part*exp_part\n",
    "\n",
    "# Use the sample mean and stdev as estimates of the population model parameters mu and sigma\n",
    "population_model = gaussian_model(sample_distances, mu=sample_mean, sigma=sample_stdev)\n",
    "\n",
    "# Plot the model and data to see how they compare\n",
    "fig = plot_data_and_model(sample_distances, population_model)\n",
    "\n",
    "\"\"\"\n",
    "Notice in the plot that the data and the model do not line up exactly. This is to be expected because the sample\n",
    "is just a subset of the population, and any model built from it cannot be a prefect representation of the population.\n",
    "Also notice the vertical axis: it shows the normalize data bin counts, and the probability density of the model. \n",
    "Think of that as probability-per-bin, so that if summed all the bins, the total would be 1.0.\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "# Compute sample mean and stdev, for use as model parameter guesses\n",
    "mu_guess = np.mean(sample_distances)\n",
    "sigma_guess = np.std(sample_distances)\n",
    "\n",
    "# For each sample distance, compute the probability modeled by the parameter guesses\n",
    "probs = np.zeros(len(sample_distances))\n",
    "for n, distance in enumerate(sample_distances):\n",
    "    probs[n] = gaussian_model(distance, mu=mu_guess, sigma=sigma_guess)\n",
    "\n",
    "# Compute the log-likelihood as the sum of the log(probs)\n",
    "loglikelihood = np.sum(np.log(probs))\n",
    "print('For guesses mu={:0.2f} and sigma={:0.2f}, the loglikelihood={:0.2f}'.format(mu_guess, sigma_guess, loglikelihood))\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Maximizing Likelihood, Part 1\n",
    "Previously, we chose the sample mean as an estimate of the population model paramter mu. \n",
    "But how do we know that the sample mean is the best estimator? This is tricky, so let's do it in two parts.\n",
    "In Part 1, you will use a computational approach to compute the log-likelihood of a given estimate. Then, \n",
    "in Part 2, we will see that when you compute the log-likelihood for many possible guess values of the estimate,\n",
    "one guess will result in the maximum likelihood.\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    " Although the likelihood (the product of the probabilities) is easier to interpret,\n",
    " the loglikelihood has better numerical properties. Products of small and large numbers can cause numerical\n",
    " artifacts, but sum of the logs usually doesnt suffer those same artifacts, and the \"sum(log(things))\"\n",
    " is closely related to the \"product(things)\"\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Maximizing Likelihood, Part 2\n",
    "In Part 1, you computed a single log-likelihood for a single mu. In this Part 2,\n",
    "you will apply the predefined function compute_loglikelihood() to compute an array of log-likelihood\n",
    "values, one for each element in an array of possible mu values.\n",
    "The goal then is to determine which single mu guess leads to the single maximum value of the loglikelihood array.\n",
    "\"\"\"\n",
    "\n",
    "# Create an array of mu guesses, centered on sample_mean, spread out +/- by sample_stdev\n",
    "low_guess = sample_mean - 2*sample_stdev\n",
    "high_guess = sample_mean + 2*sample_stdev\n",
    "mu_guesses = np.linspace(low_guess, high_guess, 101)\n",
    "\n",
    "# Compute the loglikelihood for each model created from each guess value\n",
    "loglikelihoods = np.zeros(len(mu_guesses))\n",
    "for n, mu_guess in enumerate(mu_guesses):\n",
    "    loglikelihoods[n] = compute_loglikelihood(sample_distances, mu=mu_guess, sigma=sample_stdev)\n",
    "\n",
    "# Find the best guess by using logical indexing, the print and plot the result\n",
    "best_mu = mu_guesses[loglikelihoods==np.max(loglikelihoods)]\n",
    "print('Maximum loglikelihood found for best mu guess={}'.format(best_mu))\n",
    "fig = plot_loglikelihoods(mu_guesses, loglikelihoods)\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    " Notice that the guess for mu that gave the maximum likelihood is precisely the same value as the sample.mean().\n",
    " The sample_mean is thus said to be the \"Maximum Likelihood Estimator\" of the population mean mu. We call that\n",
    " value of mu the \"Maximum Likelihood Estimator\" of the population mu because, of all the mu values tested, it\n",
    " results in a model population with the greatest likelihood of producing the sample data we have.\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Bootstrap and Standard Error\n",
    "Imagine a National Park where park rangers hike each day as part of maintaining the park trails. \n",
    "They don't always take the same path, but they do record their final distance and time. We'd like to build\n",
    "a statistical model of the variations in daily distance traveled from a limited sample of data from one ranger.\n",
    "\n",
    "Your goal is to use bootstrap resampling, computing one mean for each resample, to create a distribution of means,\n",
    "and then compute standard error as a way to quantify the \"uncertainty\" in the sample statistic as an estimator for\n",
    "the population statistic.\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "# Use the sample_data as a model for the population\n",
    "population_model = sample_data\n",
    "\n",
    "# Resample the population_model 100 times, computing the mean each sample\n",
    "for nr in range(num_resamples):\n",
    "    bootstrap_sample = np.random.choice(population_model, size=resample_size, replace=True)\n",
    "    bootstrap_means[nr] = np.mean(bootstrap_sample)\n",
    "\n",
    "# Compute and print the mean, stdev of the resample distribution of means\n",
    "distribution_mean = np.mean(bootstrap_means)\n",
    "standard_error = np.std(bootstrap_means)\n",
    "print('Bootstrap Distribution: center={:0.1f}, spread={:0.1f}'.format(distribution_mean, standard_error))\n",
    "\n",
    "# Plot the bootstrap resample distribution of means\n",
    "fig = plot_data_hist(bootstrap_means)\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that standard_error is just one measure of spread of the distribution of bootstrap resample means. \n",
    "You could have computed the confidence_interval using np.percentile(bootstrap_means, 0.95) and\n",
    "np.percentile(bootstrap_means, 0.05) to find the range distance values containing the inner 90%\n",
    "of the distribution of means.\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Estimating Speed and Confidence\n",
    "Let's continue looking at the National Park hiking data. Notice that some distances are negative because\n",
    "they walked in the opposite direction from the trail head; the data are messy so let's just focus on the overall trend.\n",
    "\n",
    "In this exercise, you goal is to use boot-strap resampling to find the distribution of speed values for a linear model,\n",
    "and then from that distribution, compute the best estimate for the speed and the 90th percent confidence interval of\n",
    "that estimate. The speed here is the slope parameter from the linear regression model to fit distance as a function\n",
    "of time.\n",
    "\"\"\"\n",
    "\n",
    "# Resample each preloaded population, and compute speed distribution\n",
    "population_inds = np.arange(0, 99, dtype=int)\n",
    "for nr in range(num_resamples):\n",
    "    sample_inds = np.random.choice(population_inds, size=100, replace=True)\n",
    "    sample_inds.sort()\n",
    "    sample_distances = distances[sample_inds]\n",
    "    sample_times = times[sample_inds]\n",
    "    a0, a1 = least_squares(sample_times, sample_distances)\n",
    "    resample_speeds[nr] = a1\n",
    "\n",
    "# Compute effect size and confidence interval, and print\n",
    "speed_estimate = np.mean(resample_speeds)\n",
    "ci_90 = np.percentile(resample_speeds, [5, 95])\n",
    "print('Speed Estimate = {:0.2f}, 90% Confidence Interval: {:0.2f}, {:0.2f} '.format(speed_estimate, ci_90[0], ci_90[1]))\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that the speed estimate (the mean) falls inside the confidence interval (the 5th and 95th percentiles).\n",
    "Moreover, notice if you computed the standard error, it would also fit inside the confidence interval. \n",
    "Think of the standard error here as the 'one sigma' confidence interval. Note that this should be very similar\n",
    "to the summary output of a statsmodels ols() linear regression model, but here you can compute arbitrary percentiles\n",
    "because you have the entire speeds distribution.\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "Visualize the Bootstrap\n",
    "Continuing where we left off earlier in this lesson, let's visualize the bootstrap distribution of speeds\n",
    "estimated using bootstrap resampling, where we computed a least-squares fit to the slope for every sample\n",
    "to test the variation or uncertainty in our slope estimation.\n",
    "To get you started, we've preloaded a function compute_resample_speeds(distances, times) to do the computation\n",
    "of generate the speed sample distribution.\n",
    "\"\"\"\n",
    "\n",
    "# Create the bootstrap distribution of speeds\n",
    "resample_speeds = compute_resample_speeds(distances, times)\n",
    "speed_estimate = np.mean(resample_speeds)\n",
    "percentiles = np.percentile(resample_speeds, [5, 95])\n",
    "\n",
    "# Plot the histogram with the estimate and confidence interval\n",
    "fig, axis = plt.subplots()\n",
    "hist_bin_edges = np.linspace(0.0, 4.0, 21)\n",
    "axis.hist(resample_speeds, hist_bin_edges, color='green', alpha=0.35, rwidth=0.8)\n",
    "axis.axvline(speed_estimate, label='Estimate', color='black')\n",
    "axis.axvline(percentiles[0], label=' 5th', color='blue')\n",
    "axis.axvline(percentiles[1], label='95th', color='blue')\n",
    "axis.legend()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that vertical lines marking the 5th (left) and 95th (right) percentiles mark the extent of the confidence\n",
    "interval, while the speed estimate (center line) is the mean of the distribution and falls between them. \n",
    "Note the speed estimate is the mean, not the median, which would be 50% percentile.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Test Statistics and Effect Size\n",
    "How can we explore linear relationships with bootstrap resampling? Back to the trail! For each hike plotted as one point,\n",
    "we can see that there is a linear relationship between total distance traveled and time elapsed. It we treat the distance\n",
    "traveled as an \"effect\" of time elapsed, then we can explore the underlying connection between linear regression and\n",
    "statistical inference.\n",
    "\n",
    "In this exercise, you will separate the data into two populations, or \"categories\": early times and late times. \n",
    "Then you will look at the differences between the total distance traveled within each population. This difference\n",
    "will serve as a \"test statistic\", and it's distribution will test the effect of separating distances by times.\n",
    "\"\"\"\n",
    "\n",
    "# Create two poulations, sample_distances for early and late sample_times.\n",
    "# Then resample with replacement, taking 500 random draws from each population.\n",
    "group_duration_short = sample_distances[sample_times < 5]\n",
    "group_duration_long = sample_distances[sample_times > 5]\n",
    "resample_short = np.random.choice(group_duration_short, size=500, replace=True)\n",
    "resample_long = np.random.choice(group_duration_long, size=500, replace=True)\n",
    "\n",
    "# Difference the resamples to compute a test statistic distribution, then compute its mean and stdev\n",
    "test_statistic = resample_long - resample_short\n",
    "effect_size = np.mean(test_statistic)\n",
    "standard_error = np.std(test_statistic)\n",
    "\n",
    "# Print and plot the results\n",
    "print('Test Statistic: mean={:0.2f}, stdev={:0.2f}'.format(effect_size, standard_error))\n",
    "fig = plot_test_statistic(test_statistic)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice again, the test statistic is the difference between a distance drawn from short duration trips and\n",
    "one drawn from long duration trips. The distribution of difference values is built up from differencing each\n",
    "point in the early time range with one from the late time range. The mean of the test statistic is not zero\n",
    "and tells us that there is on average a difference in distance traveled when comparing short and long duration\n",
    "trips. Again, we call this the 'effect size'. The time increase had an effect on distance traveled. :\n",
    "The standard error of the test statistic distribution is not zero, so there is some spread in that distribution,\n",
    "or put another way, uncertainty in the size of the effect.\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Null Hypothesis\n",
    "In this exercise, we formulate the null hypothesis as\n",
    "\n",
    "short and long time durations have no effect on total distance traveled.\n",
    "We interpret the \"zero effect size\" to mean that if we shuffled samples between short and long times,\n",
    "so that two new samples each have a mix of short and long duration trips, and then compute the test statistic,\n",
    "on average it will be zero.\n",
    "In this exercise, your goal is to perform the shuffling and resampling. Start with the predefined\n",
    "group_duration_short and group_duration_long which are the un-shuffled time duration groups.\n",
    "\"\"\"\n",
    "\n",
    "# Shuffle the time-ordered distances, then slice the result into two populations.\n",
    "shuffle_bucket = np.concatenate((group_duration_short, group_duration_long))\n",
    "np.random.shuffle(shuffle_bucket)\n",
    "slice_index = len(shuffle_bucket)//2\n",
    "shuffled_half1 = shuffle_bucket[0:slice_index]\n",
    "shuffled_half2 = shuffle_bucket[slice_index:]\n",
    "\n",
    "# Create new samples from each shuffled population, and compute the test statistic\n",
    "resample_half1 = np.random.choice(shuffled_half1, size=500, replace=True)\n",
    "resample_half2 = np.random.choice(shuffled_half2, size=500, replace=True)\n",
    "test_statistic = resample_half2 - resample_half1\n",
    "\n",
    "# Compute and print the effect size\n",
    "effect_size = np.mean(test_statistic)\n",
    "print('Test Statistic, after shuffling, mean = {}'.format(effect_size))\n",
    "\n",
    "\"\"\"\n",
    "Notice that your effect size is not exactly zero because there is noise in the data. \n",
    "But the effect size is much closer to zero than before shuffling. Notice that if you rerun your code,\n",
    "which will generate a new shuffle, you will get slightly different results each time for the effect size,\n",
    "but np.abs(test_statistic) should be less than about 1.0, due to the noise, as opposed to the slope, \n",
    "which was about 2.0\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "Visualizing Test Statistics\n",
    "In this exercise, you will approach the null hypothesis by comparing the distribution of a test statistic arrived\n",
    "at from two different ways.\n",
    "\n",
    "First, you will examine two \"populations\", grouped by early and late times, and computing the test statistic distribution.\n",
    "Second, shuffle the two populations, so the data is no longer time ordered, and each has a mix of early and late times,\n",
    "and then recompute the test statistic distribution.\n",
    "\n",
    "To get you started, we've pre-loaded the two time duration groups, group_duration_short and group_duration_long, and\n",
    "two functions, shuffle_and_split() and plot_test_statistic().\n",
    "\"\"\"\n",
    "\n",
    "# From the unshuffled groups, compute the test statistic distribution\n",
    "resample_short = np.random.choice(group_duration_short, size=500, replace=True)\n",
    "resample_long = np.random.choice(group_duration_long, size=500, replace=True)\n",
    "test_statistic_unshuffled = resample_long - resample_short\n",
    "\n",
    "# Shuffle two populations, cut in half, and recompute the test statistic\n",
    "shuffled_half1, shuffled_half2 = shuffle_and_split(group_duration_short, group_duration_long)\n",
    "resample_half1 = np.random.choice(shuffled_half1, size=500, replace=True)\n",
    "resample_half2 = np.random.choice(shuffled_half2, size=500, replace=True)\n",
    "test_statistic_shuffled = resample_half2 - resample_half1\n",
    "\n",
    "# Plot both the unshuffled and shuffled results and compare\n",
    "fig = plot_test_statistic(test_statistic_unshuffled, label='Unshuffled')\n",
    "fig = plot_test_statistic(test_statistic_shuffled, label='Shuffled')\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Notice that after you shuffle, the effect size went almost to zero and the spread increased, as measured\n",
    "by the standard deviation of the sample statistic, aka the 'standard error'. So shuffling did indeed have an effect.\n",
    "The null hypothesis is disproven. Time ordering does in fact have a non-zero effect on distance traveled.\n",
    "Distance is correlated to time.\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "Visualizing the P-Value\n",
    "In this exercise, you will visualize the p-value, the chance that the effect (or \"speed\") we estimated,\n",
    "was the result of random variation in the sample. Your goal is to visualize this as the fraction of points\n",
    "in the shuffled test statistic distribution that fall to the right of the mean of the test statistic (\"effect size\")\n",
    "computed from the unshuffled samples.\n",
    "\n",
    "To get you started, we've preloaded the group_duration_short and group_duration_long and functions\n",
    "compute_test_statistic(), shuffle_and_split(), and plot_test_statistic_effect()\n",
    "\"\"\"\n",
    "\n",
    "# Compute the test stat distribution and effect size for two population groups\n",
    "test_statistic_unshuffled = compute_test_statistic(group_duration_short, group_duration_long)\n",
    "effect_size = np.mean(test_statistic_unshuffled)\n",
    "\n",
    "# Randomize the two populations, and recompute the test stat distribution\n",
    "shuffled_half1, shuffled_half2 = shuffle_and_split(group_duration_short, group_duration_long)\n",
    "test_statistic_shuffled = compute_test_statistic(shuffled_half1, shuffled_half2)\n",
    "\n",
    "# Compute the p-value as the proportion of shuffled test stat values >= the effect size\n",
    "condition = test_statistic_shuffled >= effect_size\n",
    "p_value = len(test_statistic_shuffled[condition]) / len(test_statistic_shuffled)\n",
    "\n",
    "# Print p-value and overplot the shuffled and unshuffled test statistic distributions\n",
    "print(\"The p-value is = {}\".format(p_value))\n",
    "fig = plot_test_stats_and_pvalue(test_statistic_unshuffled, test_statistic_shuffled)\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Note that the entire point of this is compute a p-value to quantify the chance that our estimate for speed\n",
    "could have been obtained by random chance. On the plot, the unshuffle test stats are the distribution of\n",
    "speed values estimated from time-ordered distances. The shuffled test stats are the distribution of speed \n",
    "values computed from randomizing unordered distances. Values of the shuffled stats to the right of the mean\n",
    "non-shuffled effect size line are those that both (1) could have both occured randomly and (2) are at least\n",
    "as big as the estimate you want to use for speed.\n",
    "\"\"\""
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}