{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Ordinary Least-Squares with Measurement Error\n", "====\n", "\n", "## Unit 12, Lecture 4\n", "\n", "*Numerical Methods and Statistics*\n", "\n", "----\n", "\n", "#### Prof. Andrew White, April 24 2018" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Goals:\n", "---\n", "\n", "0. Be able to plot with error bars\n", "1. Be able to choose the appropiate case for OLS with measurement error\n", "2. Are aware of attenuation error and its effect on independent variable measurement error" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from math import sqrt, pi, erf\n", "import seaborn\n", "seaborn.set_context(\"notebook\")\n", "seaborn.set_style(\"whitegrid\")\n", "import scipy.stats" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Plotting with Error Bars\n", "===\n", "\n", "Error bars are little lines that go up and down or left and right in graphs to indicate the standard error of a measurement. They may also indicate a confidence interval or standard deviation, however that will usually be specified in the figure caption. The code to make error bar plots is shown below with a constant y-error bar. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEBCAYAAABysL6vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEUdJREFUeJzt3XGM33V9x/FnKeU4dKbbjIEeIBrse8RD7H4iIQVlyGhIdBQVm7FFyaZidAsxUGIzGGKywezWmZkYGzMcCUnXWMplo27VtZsT6Yr+ckoP5C3Jih3XmpGZOqrn3XG9/fG7o/djlP5+v7vffX/99PlImtzvc9/ffd55c7z66ef7+36/S6anp5Eklee0qguQJHWHAS9JhTLgJalQBrwkFcqAl6RCGfCSVCgDXpIKZcBLUqEMeEkqlAEvSYUy4CWpUKcv9oT1er0PuBQ4BEwt9vySdJJaCpwDfLdWq4238oZFD3ga4f7tCuaVpBJcCTzayoFVBPwhgJUrV3LGGWdUMP3CGBkZYXBwsOoyeob9OMZeNLMfzTrtx8TEBD/60Y9gJkNbUUXATwGcccYZ9PX1VTD9wjnZ619o9uMYe9HMfjSbZz9a3tr2JKskFcqAl6RCGfCSVCgDXpIKZcBLUqEMeEkqlAEvSYWq4nPwknRKWrd5D0eOvMCO2uLM5wpekgplwEtSoQx4SSqUAS9Ji2BoeJThA4d58vlJVt+3m6Hh0a7PacBLUpcNDY+yYfs+JqaOAjB6eIwN2/d1PeQNeEnqso07k7HJ5ptAjk1OsXFndnVeA16Suuzg4bG2xheKAS9JXbZieX9b4wvFgJekLlu/JuhftrRprH/ZUtavia7O65WsktRla1cNAHDHtieYmDrKwPJ+1q+Jl8a7xYCXpEWwdtUAWx4/0LhVwW1XL8qcbtFIUqFcwUvSItl6y+XU6/VFm88VvCQVyoCXpEIZ8JJUqJb24CPic8AHgWngbzNzU0RcA2wC+oGtmXln98qUJLXrhCv4iHg3cDXwNuAdwB9HxCXA/cD1wEXApRFxXTcLlSS154QBn5nfAn4rM18E3kBj1b8ceCYz98+MPwjc2NVKJUltaWmLJjMnI+Ie4Hbga8AK4NCcQw4B57Yz8cjISDuH96TF/LjTycB+HGMvmtmPZovVj5Y/B5+Zd0fEXwD/CKyksR8/awlwtJ2JBwcH6evra+ctPaVer1OrLdKTc08C9uMYe9HMfjTrtB/j4+NtL4xb2YP/jYh4O0Bm/gLYDlwFnDPnsLOBg23NLEnqqlZW8G8G7omIK2is2q8HNgMbI+JCYD9wE42TrpKkHtHKSdavAzuAYaAOPJaZfw/cDDwEPAU8DWzrXpmSpHa1epL1s8BnXza2C7hk4UuSJC0Er2SVpEIZ8JJUKANekgplwEtSoQx4SSqUAS9JhTLgJalQBrwkFcqAl6RCGfCSVCgDXpIKZcBLUqEMeEkqlAEvSYUy4CWpUAa8JBXKgJekQhnwklQoA16SCmXAS1KhDHhJKpQBL0mFMuAlqVAGvCQVyoCXpEKd3spBEXE38KGZlzsy846I+CpwBfDzmfF7MvPhLtQoSerACQM+Iq4BrgVWAdPAP0fEDcA7gHdl5qHulihJ6kQrK/hDwG2ZOQEQET8Ezp/5c39EDAAP01jBH+1apZKktiyZnp5u+eCIeAvwHeBK4D7gk8DPgEeALZn5lRP9jHq9fgGwv5NiJUm8qVarPdvKgS3twQNExFuBHcD6zEzghjnf+yLwYeCEAT9rcHCQvr6+Vg/vOfV6nVqtVnUZPcN+HGMvmtmPZp32Y3x8nJGRkbbe09KnaCJiNbAL+ExmPhARF0fEB+YcsgSYbGtmSVJXtXKS9TxgCFiXmbtnhpcAX4iI3cAR4OPAA12rUpLUtla2aG4HzgQ2RcTs2JeBe2nsxy8DHsrMLV2pUJLUkRMGfGbeCtx6nG9/aWHLkSQtFK9klaRCGfBSgdZt3sO6zXuqLkMVM+AlqVAGvCQVyoCXpEIZ8FJhhoZHGT5wmL37f8rq+3YzNDxadUmqiAEvFWRoeJQN2/cxMdW479/o4TE2bN9nyJ+iDHipIBt3JmOTU01jY5NTbNyZFVWkKhnwUkEOHh5ra1xlM+ClgqxY3t/WuMpmwEsFWb8m6F+2tGmsf9lS1q+J47xDJWv5fvCSet/aVQMA3LHtCSamjjKwvJ/1a+KlcZ1aDHipMGtXDbDl8QMAbL3l8oqrUZXcopGkQhnwklQot2ikArk1I3AFL0nFMuAlqVAGvCQVyoCXpEIZ8JJUKANekgplwEtSoQx4SSqUAS9JhWrpStaIuBv40MzLHZl5R0RcA2wC+oGtmXlnl2qUJHXghCv4mSC/FlgFvB2oRcTvAvcD1wMXAZdGxHXdLFSS1J5WtmgOAbdl5kRmTgI/BFYCz2Tm/sx8EXgQuLGLdUqS2nTCLZrMfHL264h4C42tmi/SCP5Zh4Bz25l4ZGSkncN7Ur1er7qEnmI/jrEXzexHs8XqR8t3k4yItwI7gPXAizRW8bOWAEfbmXhwcJC+vr523tJT6vU6tVqt6jJ6hv04xl40sx/NOu3H+Ph42wvjlj5FExGrgV3AZzLzAeA54Jw5h5wNHGxrZklSV51wBR8R5wFDwLrM3D0zvLfxrbgQ2A/cROOkqySpR7SyRXM7cCawKeKlJ7N/GbgZeGjme18HtnWhPklSh1o5yXorcOtxvn3JwpYjSVooXskqSYUy4CWpUAa8JBXKgJekQhnwklQoA16SCmXAS1KhDHhJKpQBL0mFMuAlqVAGvCQVyoCXpEIZ8JJUKANekgplwEtSoQx4SSqUAS9JhTLgJalQBrx6wrrNe1i3eU/VZUhFMeAlqVAGvCQVyoCXpEIZ8Krc0PAowwcOs3f/T1l9326GhkerLkkqggGvSg0Nj7Jh+z4mpo4CMHp4jA3b9xny0gIw4FWpjTuTscmpprGxySk27syKKpLKcXqrB0bE64DHgPdm5rMR8VXgCuDnM4fck5kPd6FGFezg4bG2xiW1rqWAj4jLgK8AK+cMvwN4V2Ye6kZhOjWsWN7P6CuE+Yrl/RVUI5Wl1S2ajwGfAg4CRMRZwPnA/RHxRETcExFu96ht69cE/cuWNo31L1vK+jVRUUVSOZZMT0+3fHBEPAtcReMvhr8CPgn8DHgE2JKZXznRz6jX6xcA+9uuVD3rT//tfwD43FW/3tH7//3Hv+BL3/tfJo/C6886jd8bfC3veuNZC1miVJI31Wq1Z1s5sOU9+Lky8z+BG2ZfR8QXgQ/T2MZpyeDgIH19fZ1M3xPq9Tq1Wq3qMnrCa7+3hyNHXui4H7Ua/MfzjdsUbL3l8oUsrRL+bjSzH8067cf4+DgjIyNtvaejbZWIuDgiPjBnaAkw2cnPkiR1R0creBqB/oWI2A0cAT4OPLBgVUmS5q2jFXxmPgHcC3wHeAr4fmZuWcjCdHKYvQr1yecn53UV6tZbLi9ie0bqJW2t4DPzgjlffwn40kIXpJPH8a5CBVi7aqDK0iThlayaB69ClXqbAa+OeRWq1NsMeHXseFebehWq1BsMeHXMq1Cl3tbpxySll06k3rHtCSamjjKwvJ/1a8ITrFKPMOA1L2tXDbDl8QMcOfICO267uupyJM3hFo0kFcqAl6RCuUWjedt6y+XU6/Wqy5D0Mq7gJalQBrwkFcqAl6RCGfCSVCgDXpIKZcBLUqEMeEkqlAEvSYUy4CWpUAa8JBXKgJekQhnwklQoA16SCmXAS1KhDHhJKpQBL0mFaumBHxHxOuAx4L2Z+WxEXANsAvqBrZl5ZxdrlCR14IQr+Ii4DHgUWDnzuh+4H7geuAi4NCKu62aRkqT2tbJF8zHgU8DBmdfvBJ7JzP2Z+SLwIHBjl+qTJHXohFs0mflRgIiYHVoBHJpzyCHg3AWvTJI0L508dPs0YHrO6yXA0XZ/yMjISAdT9xYfNN3MfhxjL5rZj2aL1Y9OAv454Jw5r8/m2PZNywYHB+nr6+tg+t5Qr9ep1WpVl9Ez7Mcx9qKZ/WjWaT/Gx8fbXhh3EvB7gYiIC4H9wE00TrpKknpI25+Dz8xfAjcDDwFPAU8D2xa2LEnSfLW8gs/MC+Z8vQu4pBsFSZIWhleySlKhOtmDP+Wt27yHI0deYIfnjST1MFfwklQoA16SCmXAS1KhDPg2DQ2PMnzgME8+P8nq+3YzNDxadUmS9IoM+DYMDY+yYfs+JqYad2YYPTzGhu37DHlJPcmAb8PGncnY5FTT2NjkFBt3ZkUVSdLxGfBtOHh4rK1xSaqSAd+GFcv72xqXpCoZ8G1YvyboX7a0aax/2VLWr4njvEOSquOVrG1Yu2oAgDu2PcHE1FEGlvezfk28NC5JvcSAb9PaVQNsefxA41YFt11ddTmSdFxu0UhSoQx4SSqUWzQd2HrL5T5jUlLPcwUvSYUy4CWpUAa8JBXKgJekQhnwklQoA16SCnVSBfy6zXtYt3lP1WVI0knhpAp4SVLrDHhJKtS8rmSNiH8F3gBMzgzdkpl7512VJGneOg74iFgCrATemJkvLlxJr2z2YdcTU0dZfd9ub9MrSScwny2a2adcfCMifhARf7QQBb0SH3YtSe2bT8D/KrALuAF4D/CJiPjtBanqZXzYtSS1b8n09PSC/KCI+DRwfmZ++tWOq9frFwD72/nZH/zaT3ilKpcA2248u50fJUknuzfVarVnWzlwPnvwVwB9mblrZmgJx062ntDg4CB9fX0tHbvim7sZPTz2/8eX91Or1VqdckHV6/XK5u5F9uMYe9HMfjTrtB/j4+OMjIy09Z75bNEsBzZGxJkR8SvAR4CH5/HzjsuHXUtS+zoO+Mx8BNgBDAN14P7M7MplpmtXDXDv+y/mjKWNcgeW93Pv+y/2UzSS9Crm9Tn4zLwLuGuBanlVsw+7hsYTlSRJr84rWSWpUCfVM1lduUtS61zBS1KhDHhJKpQBL0mFMuAlqVAGvCQVyoCXpEIZ8JJUKANekgpVxYVOSwEmJiYqmHphjY+PV11CT7Efx9iLZvajWSf9mJOZS1/tuLkW7H7wrarX61cA317USSWpHFfWarVHWzmwihX8d4ErgUPA1AmOlSQ1LAXOoZGhLVn0FbwkaXF4klWSCmXAS1KhDHhJKpQBL0mFMuAlqVAGvCQVyoCXpEKdVM9k7QURcTfwoZmXOzLzjirr6RUR8ZfA6zPz5qprqVJEvA+4G3gN8I3MvLXikioTEb8PbJh5+U+ZeXuV9VQlIl4HPAa8NzOfjYhrgE1AP7A1M+/s1tyu4Nsw8x/mWmAV8HagFhE3VFtV9SLiPcBHqq6jahHxZuDLwFrgbcBvRsR11VZVjYg4C/gb4N3AJcCVM///nFIi4jLgUWDlzOt+4H7geuAi4NJu/o4Y8O05BNyWmROZOQn8EDi/4poqFRG/BvwZ8OdV19IDbqCxIntu5vdjHbC34pqqspRGvrwGWDbzZ6zSiqrxMeBTwMGZ1+8EnsnM/Zn5IvAgcGO3JneLpg2Z+eTs1xHxFhpbNaurq6gnbAb+BDiv6kJ6wIXARET8A42/+B8B7qq2pGpk5gsRcRfwNPAL4Fs0tilOKZn5UYCImB1aQWOhOOsQcG635ncF34GIeCvwTWB9Zj5TdT1ViYiPAv+VmbuqrqVHnA5cA/whcDlwGafo1lVEvA34A+CNNEJtCjgl9+Bf5jRg7g3AlgBHuzmZ2hARq4FdwGcy84Gq66nYOuDaiPg+8DngdyLiryuuqUo/Af4lM5/PzDHgYRr/JD8VrQF2ZeZ/Z+Y48HfAVZVW1Bueo3FHyFlnc2z7ZsG5RdOGiDgPGALWZebuquupWmb+9uzXEXEzcFVmfrq6iir3CPBARCwHXgCuo/H7cir6AfD5iHgNjS2a99HGbW4LtheIiLgQ2A/cROOka1e4gm/P7cCZwKaI+P7Mn09UXZR6Q2buBT5P41MTTwE/Br5aaVEVycxvAFuAOvAEjZOs91VaVA/IzF8CNwMP0fgdeRrY1q35vB+8JBXKFbwkFcqAl6RCGfCSVCgDXpIKZcBLUqEMeEkqlAEvSYUy4CWpUP8HXhG9lydeNpQAAAAASUVORK5CYII=\n", "text/plain": [ "