{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Ordinary Least-Squares in 1-D and N-D\n", "====\n", "\n", "## Unit 12, Lecture 1\n", "\n", "*Numerical Methods and Statistics*\n", "\n", "----\n", "\n", "### Reading\n", "\n", "Bulmer: Chapter 12 (Regression & Correlation)\n", "\n", "----\n", "\n", "\n", "#### Prof. Andrew White, April 12 2018" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Goals:\n", "---\n", "\n", "1. Understand the nomenclature of the regression equation and terms\n", "2. Be able to recognize which terms represent noise and which are model terms\n", "3. Perform linear regression in 1 or more dimensions\n", "4. Be able to construct probability distributions representing parameter uncertainty\n", "5. Compute hypothesis tests and confidence intervals" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "%matplotlib inline\n", "import random\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from math import sqrt, pi, erf\n", "import scipy.stats" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Least Squares Linear Regression\n", "====\n", "\n", "Let's define the problem setting for linear regression." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We have some data, $x$ and $y$, and we want to fit a line to the data. Our model will be:\n", "\n", "$$\\hat{y} = \\hat{\\alpha} + \\hat{\\beta}x$$\n", "\n", "where a $\\hat{}$ indicates our best estimate of something. We make an assumption that the process that generates our data looks like:\n", "\n", "$$y = \\alpha + \\beta x + \\epsilon$$\n", "\n", "where $\\epsilon$ is coming from a normal distribution" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Nomenclature\n", "-----\n", "\n", "* $\\alpha$ - an intercept\n", "* $\\hat{\\alpha}$ - our estimated best intercept\n", "* $\\beta$ - a slope\n", "* $\\hat{\\beta}$ - our estimated best slope\n", "* $\\hat{y}$ - our estimated y's\n", "* $\\epsilon$ - the variable which is creating normally distributed noise.\n", "* $x,y$ - the data\n", "* Residual - the difference between $\\hat{y}$ and $y$ at a particular point ($\\hat{y}_i - y_i$)\n", "* Sum of squared residues (SSR) - $\\sum_i (\\hat{y}_i - y_i)^2$\n", "* Sum of squared error (SSE) - the same as SSR\n", "* Total sum of squares (TSS) - $\\sum_i (\\bar{y} - y_i)^2$, which is a kind of variance\n", "* $S_{\\textrm{something}}$ - The standard error of something" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Ordinary Least Squares Regression (OLS-1D)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "One way to view this problem is as an optimization. Can we write down an equation that results in a single value that we minimize? Our ultimate goal is to make $\\hat{y}$ as close as possible to $y$. Mathematically, that means we want $\\sum_i(y_i - \\hat{y}_i)^2$ to be small. \n", "\n", "That is the objective function to minimize, but what are the dimensions to change? Those are $\\hat{\\alpha}$ and $\\hat{\\beta}$. So we need to write down a function which takes in $\\hat{\\alpha}$ and $\\hat{\\beta}$ and returns how good the fit is:" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "$$f(\\hat{\\alpha}, \\hat{\\beta}) = \\sum_i \\left(y_i - \\hat{\\alpha} - \\hat{\\beta} x_i\\right)^2$$\n", "\n", "We can minimize this equation using any of our minimization techniques or we can do it analytically. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Using calculus you can show that the minimum to $f(\\alpha, \\beta)$ is:\n", "\n", "$$\\hat{\\beta} = \\frac{\\sum_i(x_i - \\bar{x})(y_i - \\bar{y})}{\\sum_i(x_i - \\bar{x})^2}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "With a little bit of algebra, you can show this is\n", "\n", "$$\\hat{\\beta} = \\frac{\\sigma_{xy}}{\\sigma_x^2}$$\n", "\n", "where $\\sigma_{xy}$ is the sample covariance of $x$ and $y$ and $\\sigma_x^2$ is the sample variance of $x$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "To find the intercept, you can just take the average of the residuals (not their squares!) given the model so far:\n", "\n", "$$\\hat{\\alpha} = \\frac{1}{N}\\sum_i (y_i - \\hat{\\beta}x_i)$$\n", "\n", "Let's see this in action" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEIFJREFUeJzt3W9oXfd9x/HPp4rL1GRDCVaCrSSzN4LWMNOoiJBNULJmnZJs1K6h0MBC2Arug2RLR9Cw+yRlY9jgLt0elILbeDEsSymJqpgl1A12WOiTUDkKtTPPJGRp4msvVghawxDEdr57oKtWdvTnXp17/v3O+wXi3nt0pPO92P743O/vnN/PESEAQP19ouwCAAC9QaADQCIIdABIBIEOAIkg0AEgEQQ6ACSCQAeARBDoAJAIAh0AEnFVkQfbuHFjbNmypchDAkDtHT9+/L2IGFxrv0IDfcuWLZqeni7ykABQe7Z/0cl+tFwAIBEEOgAkgkAHgEQQ6ACQCAIdABJR6FUuANA0UzMt7T9yWmfn5rV5oF8T48PaMTKUy7EIdADIydRMS3smT2j+wiVJUmtuXnsmT0hSLqFOywUAcrL/yOlfhfmi+QuXtP/I6VyOt2ag277J9ou2T9l+zfbD7e3ftN2y/Wr7695cKgSAmjo7N9/V9qw6ablclPRIRLxi+zclHbf9Qvt7346Ib+VSGQDU3OaBfrWWCe/NA/25HG/NM/SIOBcRr7SffyDplKR8OvoAkJCJ8WH1b+i7bFv/hj5NjA/ncryueui2t0gakfRye9NDtn9u+6Dta1f4mV22p21Pz87OZioWAOpkx8iQ9u7cpqGBflnS0EC/9u7clttVLo6Izna0r5H0H5L+ISImbd8g6T1JIenvJW2KiL9c7XeMjo4Gk3MBQHdsH4+I0bX26+gM3fYGSc9IejIiJiUpIt6NiEsR8ZGk70m6PUvBAIBsOrnKxZIel3QqIh5bsn3Tkt2+JOlk78sDAHSqk6tcxiTdL+mE7Vfb274h6T7bt2mh5fKWpK/lUiEAoCNrBnpE/FSSl/nW870vBwCwXtwpCgCJYC4XACsqcmIpZEegA1hW0RNLITsCHcCyVptYqkmBXqdPKQQ6gGUVPbFUFdXtUwqDogCWtdIEUnlNLFVFRU9/mxWBDmBZRU8sVUV1+5RCoANYVtETS1VR3T6l0EMHsKIdI0ONCvArTYwPX9ZDl6r9KYVAB4AVLP5nxlUuAJCAOn1KoYcOAIngDB1Asup0U1AvEOgAklS3m4J6gZYLgCTV7aagXuAMHUBlZWmZ1O2moF7gDB1AJS22TFpz8wr9umUyNdPq6OfrdlNQLxDoACopa8ukiVMX0HIBUElZWyZ1uymoFwh0AJW0eaBfrWXCu5uWSZ1uCuoFWi4AKqmJLZOsOEMHUElNbJlkRaADqKymtUyyouUCAIkg0AEgEQQ6ACSCQAeARKwZ6LZvsv2i7VO2X7P9cHv7dbZfsP16+/Ha/MsFAKykkzP0i5IeiYhPS7pD0oO2b5W0W9LRiLhF0tH2awBASdYM9Ig4FxGvtJ9/IOmUpCFJ2yUdau92SNKOvIoEAKytq+vQbW+RNCLpZUk3RMQ5aSH0bV/f8+qAhmvaijvIpuNAt32NpGckfT0ifmm705/bJWmXJN18883rqRFopCauuINsOrrKxfYGLYT5kxEx2d78ru1N7e9vknR+uZ+NiAMRMRoRo4ODg72oGWiEJq64g2w6ucrFkh6XdCoiHlvyrcOSHmg/f0DSs70vD2iuJq64g2w6OUMfk3S/pM/bfrX9da+kfZK+YPt1SV9ovwbQI01ccQfZrNlDj4ifSlqpYX5Xb8sBsGhifPiyHrrE9LFYHbMtAhXF9LHoFoEOVBjTx6IbzOUCAIkg0AEgEQQ6ACSCQAeARBDoAJAIAh0AEkGgA0AiCHQASAQ3FgEJYz71ZiHQgUQxn3rz0HIBEsV86s1DoAOJYj715iHQgUQxn3rzEOhAoibGh9W/oe+ybUXPpz4109LYvmPauvs5je07pqmZVmHHbiIGRYFElT2fOoOyxSPQgYSVOZ/6aoOyBHo+CHQgR02+DpxB2eLRQwdysthyaM3NK/TrlkNT+sgMyhaPQAdy0vTrwKswKNs0tFyAnDS95VD2oGwTEehATjYP9Ku1THg3qeXAItfFouUC5ISWA4rGGTqQE1oOKBqBDuSIlgOKRMsFABLBGTqS1uQbe9A8a56h2z5o+7ztk0u2fdN2y/ar7a978y0T6F7Tb+xB83TScnlC0t3LbP92RNzW/nq+t2UB2TX9xh40z5qBHhEvSXq/gFqAnmr6jT1oniyDog/Z/nm7JXNtzyoCeoS5RNA06w3070r6XUm3STon6R9X2tH2LtvTtqdnZ2fXeTige9zYg6ZZV6BHxLsRcSkiPpL0PUm3r7LvgYgYjYjRwcHB9dYJdG3HyJD27tymoYF+WdLQQL/27tzGVS5I1rouW7S9KSLOtV9+SdLJ1fYHysKNPWiSNQPd9lOS7pS00fYZSY9KutP2bZJC0luSvpZjjQCADqwZ6BFx3zKbH8+hFgBABtz6DwCJINABIBEEOgAkgkAHgEQQ6ACQCKbPBVbB9LuoEwIdWMHi9LuLMzYuTr8riVBHJdFyAVbA9LuoGwIdWAHT76JuCHRgBUy/i7oh0IEVMP0u6oZBUWAFiwOfXOWCuiDQgVUw/S7qhJYLACSCQAeARBDoAJAIAh0AEkGgA0AiCHQASASBDgCJINABIBEEOgAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEgE0+ei0qZmWsxHDnRozTN02wdtn7d9csm262y/YPv19uO1+ZaJJpqaaWnP5Am15uYVklpz89ozeUJTM62ySwMqqZOWyxOS7r5i225JRyPiFklH26+Bj5maaWls3zFt3f2cxvYd6yqM9x85rfkLly7bNn/hkvYfOd3rMoEkrBnoEfGSpPev2Lxd0qH280OSdvS4LiQg6xn22bn5rrYDTbfeQdEbIuKcJLUfr19pR9u7bE/bnp6dnV3n4VBHWc+wNw/0d7UdaLrcr3KJiAMRMRoRo4ODg3kfDhWS9Qx7YnxY/Rv6LtvWv6FPE+PDmWsDUrTeQH/X9iZJaj+e711JSEXWM+wdI0Pau3Obhgb6ZUlDA/3au3MbV7kAK1jvZYuHJT0gaV/78dmeVYRkTIwPa8/kicvaLt2eYe8YGSLAgQ6tGei2n5J0p6SNts9IelQLQf5D21+V9LakL+dZJOppMYi5jhwohiOisIONjo7G9PR0YccDgBTYPh4Ro2vtx63/AJAIAh0AEkGgA0AiCHQASASBDgCJINABIBEEOgAkgkAHgESwYhFWxYpBQH0Q6BVXZqAuzme+OBfL4nzmkgh1oIJouVRY2UuwsWIQUC8EeoWVHaisGATUC4FeYWUHKisGAfVCoFdY2YHKikFAvRDoFVZ2oLJiEFAvXOVSYVVYIIIVg4D6INArjkAF0CkCPXHcGAQ0B4GeMG4MApqFQdGElX0dO4BiEegJK/s6dgDFItATVvZ17ACKRaAnrOzr2AEUi0HRhFXhOnYAxSHQE8d17EBz0HIBgEQQ6ACQCAIdABKRqYdu+y1JH0i6JOliRIz2oigAQPd6MSj6RxHxXg9+DwAgA1ouAJCIrIEekn5i+7jtXcvtYHuX7Wnb07OzsxkPBwBYSdZAH4uIz0q6R9KDtj935Q4RcSAiRiNidHBwMOPhAAAryRToEXG2/Xhe0o8k3d6LogAA3Vv3oKjtqyV9IiI+aD//E0l/17PKEsECEwCKkuUqlxsk/cj24u/5t4j4cU+qSgQLTAAo0roDPSLelPSZHtaSnNUWmCDQAfQaly3miAUmABSJQM8RC0wAKFLygT4109LYvmPauvs5je07pqmZVmHHZoEJAEVKej70sgclWWACQJGSDvQqDEqywASAoiTdcmFQEkCTJB3oDEoCaJKkA70Xg5JlDqoCQDeS7qFnHZQse1AVALqRdKBL2QYlqzCoCgCdSrrlkhWDqgDqhEBfBYOqAOqEQF8Fd3oCqJPke+hZcKcngDoh0NfAnZ4A6oKWCwAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEgEgQ4AiSDQASARBDoAJIJAB4BEEOgAkAgCHQASkSnQbd9t+7TtN2zv7lVRAIDurTvQbfdJ+o6keyTdKuk+27f2qjAAQHeynKHfLumNiHgzIj6U9ANJ23tTFgCgW1kCfUjSO0ten2lvAwCUIEuge5lt8bGd7F22p21Pz87OZjgcAGA1WQL9jKSblry+UdLZK3eKiAMRMRoRo4ODgxkOBwBYTZZA/5mkW2xvtf1JSV+RdLg3ZQEAurXuNUUj4qLthyQdkdQn6WBEvNazygAAXcm0SHREPC/p+R7VAgDIgDtFASARBDoAJIJAB4BEZOqhF2FqpqX9R07r7Ny8Ng/0a2J8WDtGuH8JAK5U6UCfmmlpz+QJzV+4JElqzc1rz+QJSSLUAeAKlW657D9y+ldhvmj+wiXtP3K6pIoAoLoqHehn5+a72g4ATVbpQN880N/VdgBoskoH+sT4sPo39F22rX9DnybGh0uqCACqq9KDoosDn1zlAgBrq3SgSwuhToADwNoq3XIBAHSOQAeARBDoAJAIAh0AEkGgA0AiHPGxdZ3zO5g9K+kX6/zxjZLe62E5dcB7bgbeczNkec+/HRFrLspcaKBnYXs6IkbLrqNIvOdm4D03QxHvmZYLACSCQAeARNQp0A+UXUAJeM/NwHtuhtzfc2166ACA1dXpDB0AsIpaBLrtu22ftv2G7d1l15M32zfZftH2Kduv2X647JqKYLvP9oztfy+7liLYHrD9tO3/av9Z/0HZNeXN9t+0/06ftP2U7d8ou6Zes33Q9nnbJ5dsu872C7Zfbz9em8exKx/otvskfUfSPZJulXSf7VvLrSp3FyU9EhGflnSHpAcb8J4l6WFJp8ouokD/LOnHEfF7kj6jxN+77SFJfy1pNCJ+X1KfpK+UW1UunpB09xXbdks6GhG3SDraft1zlQ90SbdLeiMi3oyIDyX9QNL2kmvKVUSci4hX2s8/0MI/9KTnELZ9o6Q/lfT9smspgu3fkvQ5SY9LUkR8GBFz5VZViKsk9du+StKnJJ0tuZ6ei4iXJL1/xebtkg61nx+StCOPY9ch0IckvbPk9RklHm5L2d4iaUTSy+VWkrt/kvS3kj4qu5CC/I6kWUn/0m4zfd/21WUXlaeIaEn6lqS3JZ2T9L8R8ZNyqyrMDRFxTlo4YZN0fR4HqUOge5ltjbg0x/Y1kp6R9PWI+GXZ9eTF9p9JOh8Rx8uupUBXSfqspO9GxIik/1NOH8Orot033i5pq6TNkq62/eflVpWWOgT6GUk3LXl9oxL8mHYl2xu0EOZPRsRk2fXkbEzSF22/pYWW2udt/2u5JeXujKQzEbH4yetpLQR8yv5Y0n9HxGxEXJA0KekPS66pKO/a3iRJ7cfzeRykDoH+M0m32N5q+5NaGEQ5XHJNubJtLfRWT0XEY2XXk7eI2BMRN0bEFi38+R6LiKTP3CLifyS9Y3txxfO7JP1niSUV4W1Jd9j+VPvv+F1KfCB4icOSHmg/f0DSs3kcpPJrikbERdsPSTqihVHxgxHxWsll5W1M0v2STth+tb3tGxHxfIk1off+StKT7ROVNyX9Rcn15CoiXrb9tKRXtHAl14wSvGPU9lOS7pS00fYZSY9K2ifph7a/qoX/2L6cy7G5UxQA0lCHlgsAoAMEOgAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEgEgQ4Aifh/89OQK7JVGsgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Make some data -> this is problem setup\n", "# Do NOT copy this because it only generates random data\n", "# This does not perform regression\n", "\n", "x = np.linspace(0,10, 20)\n", "y = 1 + x * 2.5 + scipy.stats.norm.rvs(scale=2, size=20)\n", "plt.plot(x,y, 'o')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "alpha_hat = 0.79 (1)\n", "beta_hat = 2.5 (2.5)\n" ] } ], "source": [ "cov = np.cov(x,y, ddof=2)\n", "\n", "#recall that the diagonal is variances, so we use that directly\n", "beta_hat = cov[0,1] / cov[0,0]\n", "alpha_hat = np.mean( y - beta_hat * x)\n", "\n", "print(f'alpha_hat = {alpha_hat:.2} ({1})')\n", "print(f'beta_hat = {beta_hat:.2} ({2.5})')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAH3RJREFUeJzt3Xl4VOX5xvHvkxAgsiOLEEAQKYq4INGiKKIgi7tYFbQutS1d1NZWWd1wBUWpC1qLimLdq4hUgbALiKBBrIiIIiAQtiCyL9ne3x9n8MeSkEkyM2fOzP25Lq7MnJyT8wzLzTvvvIs55xARkeBL8bsAERGJDAW6iEiCUKCLiCQIBbqISIJQoIuIJAgFuohIglCgi4gkCAW6iEiCUKCLiCSISrG8Wb169Vzz5s1jeUsRkcBbsGDBJudc/dLOi2mgN2/enOzs7FjeUkQk8Mzsh3DOU5eLiEiCUKCLiCQIBbqISIJQoIuIJAgFuohIglCgi4gkCAW6iEiCUKCLiETTjlyYOAD2bI36rRToIiLRUFgA8/8FT7eHz16AH+ZG/ZYKdBGRSPthLow6Byb2Z2PNNvSp9A9avFREx2HTGbcwJ2q3jenUfxGRhLZ9PUy+Gxa9DbWaMv+0J7lxXkN25xcBkLNlN4PGLgLgsnYZEb99qS10M2tqZjPMbImZLTazv4aODzGzHDP7IvTrgohXJyISBIX5MPdpr3vl63HQqR/c/Cl/X9Ts5zDfZ3d+IcOzlkaljHBa6AXA7c65z82sBrDAzKaEvvcP59xjUalMRCQIls+ECf1h01Jo1R16DIUjWwKwdsvuYi8p6XhFlRrozrl1wLrQ4+1mtgSI/HsFEZEg2boGsu70WuR1mkOft6B1jwNOaVw7nZxiwrtx7fSolFSmD0XNrDnQDpgfOnSLmX1pZqPNrE6EaxMRiT8Fe2HWYzDyNPg2C869C/48/5AwB+jXvTXpaakHHEtPS6Vf99ZRKS3sQDez6sC7wG3OuW3AP4GWwCl4LfjHS7iur5llm1l2bm5uBEoWEfHJd1Pg2Q4w/QE4tgvc8imc0w/SqhZ7+mXtMhja60QyaqdjQEbtdIb2OjEqH4gCmHOu9JPM0oAPgCzn3Ihivt8c+MA51/ZwPyczM9NpgwsRCZzNKyBrMCydAEe2gp6PeIEeI2a2wDmXWdp5pfahm5kBLwJL9g9zM2sU6l8HuBz4qrzFiojEpfzdMOcJ+PgJsFToeh90+DNUqux3ZcUKZ5RLR+A6YJGZfRE6NhjoY2anAA5YCfwhKhWKiMSac15rfNJA2LIK2l4B3R6Emo39ruywwhnlMgewYr41IfLliIj47MfvYWJ/WDYV6h8PN3wALc72u6qwaKaoiAhA3k5v9MonI6FSVeg+FE7/PaSm+V1Z2BToIlKicQtzGJ61lLVbdtO4djr9ureO2ggN3zjnjSXPuhO25cDJ10DXIVCjod+VlZkCXUSKNW5hDoPGLmJ3fiEQ/XVIfLHxG5jYD1bMgqNOhF+NhmYd/K6q3BToIlKs4VlLfw7zffatQxL4QN+zDT56BOY/B5WrwQWPQeZNkJJ6yKlBepeiQBeRYsV6HZKYcA6+fBum3A07NsKp10OXe6BavWJPD9q7FK2HLiLFKmm9kWitQxJ16xfBSz3hvb5QMwN+Nw0uearEMIfDv0uJRwp0ESlWrNchiZrdW7zVEP/VCXKXwsVPeWHepH2plwbtXYq6XESkWPu6FILSf3yIoiL44jWYOgR2b/b6yM+9E46oG/aPiPVqiRWlQBeREl3WLiM4Ab6/nM9hQj/IyYamv4QLxkKjk8v8Y/p1b31AHzrE97sUBbqIJI6dP8L0+2HBGKhWHy57Dk7uDVbcZPfSBe1digJdRIKvqBAWvOwta7tnG3T4E3QeCFVrVfhHB+ldigJdRIJt9acw4Q5Y9z9ofjb0fBQatvG7Kl8o0EUkmHZs9D7w/OI1qNHYm+V5Qq8DuleCNCkoEhToIhIshQXw2fMw42FvvfKOt0GnflCl+gGnBW1SUCQo0EUkOFbO8UavbPwaWp7nda/Ua1XsqQm9dEEJFOgiErf2dZkUbFnLA0e8Sbei2VCrGVz9Khx30WFHrwRtUlAkKNBFJC6NW5jDPWMX0rvoA/5S5T3SCgt5xl1Bs7MGc/Hxx5Z6fdAmBUWCpv6LSFyaOeEt3rM7GJz2BvOKjuf8vEcZnncFw6atCuv6hFm6oAzUQheR+LJlNWQN5on88aykIb/J68eMonY/fzvcLpOgTQqKBAW6iMSH/D3wydMw63EA/lXpGkbs6MZeKh9wWlm6TII0KSgSFOgi4r9vs2DiAPhpBRx/CXR/mIYrUkgZuwgCso5KPFCgi4h/Ni+HSYPg20lwZCu47j1vOCJwWaiXJZm6TCpKgS4isZe3C+b8Az5+ElLT4Pz74Zd/gkoHdq8kW5dJRSnQRSR2nINvPoBJg2HrKjjxSi/Mazb2u7KEoEAXkdjY9B1M7A/fT4cGbeDGD6H5WX5XlVAU6CISXXt3wKzh8MkzkJYOPR6B034HqYqfSNPvqIhEh3Pw1bsw+W7YvhZOuRa6DoHqDfyuLGGVGuhm1hR4BTgKKAJGOeeeNLO6wFtAc2AlcJVz7qfolSoigbHha697ZeVsb+u3q8ZA09P9rirhhTP1vwC43Tl3PNABuNnM2gADgWnOuVbAtNBzEUlme7Z6wxCfOwvWL4ILR8DvZyjMY6TUFrpzbh2wLvR4u5ktATKAS4HOodPGADOBAVGpUkTiW1ERfPkWTLkHduZC+xvhvLuh2pF+V5ZUytSHbmbNgXbAfKBhKOxxzq0zM3WMiURYIHbcWfeltwXc6vmQkQnXvAUZp/pdVVIKO9DNrDrwLnCbc26bhbmLtpn1BfoCNGvWrDw1iiSluN9xZ9dmmPEQZI+G9Lpw6TNw8jWQokVc/RLW77yZpeGF+WvOubGhwxvMrFHo+42AjcVd65wb5ZzLdM5l1q9fPxI1iySFw+2446uiIljwMjzd3gvz034Pt2ZDu18rzH0WzigXA14EljjnRuz3rfHADcCw0Nf3o1KhSJKKyx131izwulfWfg7NzoQLHoWjTvSvHjlAOF0uHYHrgEVm9kXo2GC8IH/bzH4LrAKujE6JIskprnbc2bkJpg6Bhf+G6g2h1/PetP0wu14lNsIZ5TIHKOlPrUtkyxGRffp1b31AHzr4sHxsYQEseAmmPwB5O+GMW+CcAVC1ZuxqkLBppqhInPJ9x51V8+DDO2DDImjRCXoOhwbHxebeUi4KdJE45svysdvXw5R74cs3oWYTuHIMtLlU3SsBoEAXEU9hPnw6CmYMhcK9cPbt3q/K1fyuTMKkQBcRWDELJvSH3CVw7PnQ8xE4sqXfVUkZKdBFktnWHJh8FyweC7WPht5vQOue6l4JKAW6SDIq2OutTz5rOLgi6DwIOv7VW69cAkuBLpJslk2FiQPgx2Vw3EXQ/SGo09zvqiQCFOgiyeKnHyBrsLenZ92WcO270Kqr31VJBCnQRRJd/h74+EmYMwIsBbrcC2fcDJWq+F2ZRJgCXSRROce8Sa/R9NP7yXAbmJbSkYKu99P9zEy/K5MoUaCLJKIfv2f927fRYcMsvi3KoE/BnXxSdALpEzcxND0nPpbflYhToIskkrydMHsEzH2KGoWpPJB/LWMKu1MQ+qe+b/ldBXpiUqCLJALnYMl4mDQYtq2Bk3pz7qdns5E6h5zq6/K7ElUKdJGgy10KE/vD8pnQsC1c8QIcfQZp306HeFl+V2JCgS4SVHu3w0ePwLx/euutXPAYtP8NpHr/rONh+d1A7ImaQBToIkHjHCx6x5uyv2M9tLvOG4pY/cAtHv1efjfu90RNQAp0kSDZsBgm9IMfPobG7aD3a9Ck5GGIviy/G3K4PVEV6NGhQBeJooh1OezeAjOHwqfPQ9VacPGTXss8JTXyRUdIXO6JmuAU6CJREpEuh6Ii+N8bMPVeb1/PzJvgvLvgiLrRKjti4mpP1CSR4ncBIonqcF0OYVm7EEZ3g/f/DHVaQN+ZcNGIQIQ5eB/Kpqcd+A4i5nuiJhm10EWipNxdDrs2e5syZ78E1erDZc/BSVdDSrDaX35/KJuMFOgiUVLmLoeiQvh8DEy7H/Zsgw5/gs4DvT7zgPLzQ9lkFKz/8kUCpExdDqs/g+fPgw/+Bg1OgD/Ohh5DAx3mEntqoYtESVhdDjtyYeoQ+OJVqNEIrngR2l6hLeCkXBToIlFUYpdDYQFkvwjTH4L8Xd72b536Q5XqsS9SEoYCXSTWVn7sTQ7auBiOORd6Pgr1f+F3VZIAFOiS0OJqLZFt62DKPbDobajVFK76Nxx/sbpXJGJKDXQzGw1cBGx0zrUNHRsC/B7IDZ022Dk3IVpFipRH3KwlUpjvLaD10SPe40794ay/QeUjYleDJIVwRrm8DPQo5vg/nHOnhH4pzCXuVHhiTyR8PwP+2RGm3A3Nz4Kb58F5dyrMJSpKbaE752aZWfPolyISWb6uJbJlNUy+E75+H+o0hz5vQevi2kUikVORPvRbzOx6IBu43Tn3U4RqEokIX9YSKdgLc5+G2Y97y9yeexeceSukVY3ePUVCyjux6J9AS+AUYB3weEknmllfM8s2s+zc3NySThOJuJivJfLdFHi2gzdt/9iucMuncE4/hbnETLla6M65Dfsem9nzwAeHOXcUMAogMzPTled+IuURs7VENq+ArMGwdAIc2Qp+PRaO7RLZe4iEoVyBbmaNnHPrQk8vB76KXEkikRPVtUTyd8Ocf8CcJyClEnS9Dzr8GSpVjs79REoRzrDFN4DOQD0zWwPcC3Q2s1MAB6wE/hDFGkXii3PwzYeQNQi2rIK2v4JuD0DNxn5XJkkunFEufYo5/GIUahGJf5uWwaQBsGwqNGgDN37oDUcUiQOaKSoSjr07YPZjMHckpKVDj2Fw2u8gNc3vykR+pkAXORznYPF7MPku2JYDJ18D598H1Rv4XZnIIRToIiXZuMRbRGvlbDjqJPjVS9Dsl35XJVIiBbrIwfZs89Zdmf8cVK4OF46A9jdCSmqpl4r4SYEuso9z8OVbMPlu2JkL7W+A8+6Bakf6XZlIWBToIgDrvvS6V1bPg4z2cM2bkNE+tPzu9PhYflekFAp0SW67f/J2Dcp+EdLrwCUj4ZRrISUlfpbfFQmTAl2SU1GRt4/n1CFeqJ/2Ozh3sBfqIYdbfleBLvFIgS7JJ2cBfHgHrP0cmp0BFwyHo0485DRfl98VKQcFuiSPnT/CtPvg81e8ceSXj4KTripxCzhflt8VqYDyLp8rEhxFhfDZC/D0qfDFa3DGzXBLNpx89WH384z58rsiFaQWuiS2VfNhwu2wfhG06AQ9h0OD48K6NGbL74pEiAJdEtP2DTD1XvjfG1Azw5vlecLlh22RFyeqy++KRJgCXRJLYT58+jzMHOqtV37W3+Hs26FKdb8rE4k6BbokjhWzYWJ/2Pi1twVcj0eg3rF+VyUSMwp0Cb6tOTDlbvjqXajdDHq/Dq0vKHP3ikjQKdAluAryYN4z8NFwKCqAcwbCWbd565WLJCEFugTTsmle98qPy6D1hdD9Iajbwu+qRHylQJdg2bIKsgbDkv9C3WPg2neg1fl+VyUSFxToEgz5e2DuUzD7cbAU6HIPnHELVKrid2UicUOBLvFv6SRvY+afVnpjybs9CLWa+F2VSNxRoEv82rwcJg6E77KgXmu4/n04prPfVYnELQW6xJ+8XTBnBHz8JKRW9lrkv/wjpKb5XZlIXFOgS/xwDpaMh6w7YetqOOlqOP9+qHGU35WJBIICXeJD7rfeMMTlM6BhW+g1Co4+0++qRAJFgS7+2rsdPnoU5j0LadWg56OQ+VtI1V9NkbLSvxrxh3PeVP3Jd8H2ddDu19BlCFSvf8Bp3ibNWr5WJBylBrqZjQYuAjY659qGjtUF3gKaAyuBq5xzP0WvTEkoGxbDhH7ww8fQ6BS4+lVoknnIadqkWaRswmmhvwyMBF7Z79hAYJpzbpiZDQw9HxD58iTo9m9h/6JWEc9mZNFyxetQtSZc9AScej2kpBZ7rTZpFimbUgPdOTfLzJofdPhSoHPo8RhgJgp0Oci+Fvae/Hx6pcxh4J7XOfL77axocRUtrhoGR9Q97PXapFmkbMrbh97QObcOwDm3zswaRLAmSRDDs5ZyTMEy7q/8Mu1TvuPzomO5MX8AW9a34eNSwhy0SbNIWUV9k2gz62tm2WaWnZubG+3bSbzYtZk/7RjJfyvfxdG2gX75fbkibwiLXYuwW9japFmkbMrbQt9gZo1CrfNGwMaSTnTOjQJGAWRmZrpy3k+CoqgQPn8Fpt1P70pbeLmgO08UXME2qv18SrgtbG3SLFI25Q308cANwLDQ1/cjVpEE15psmHAHrF0IR3dkVos7GD4tj938/webZW1ha5NmkfCFM2zxDbwPQOuZ2RrgXrwgf9vMfgusAq6MZpES53bkwrQhsPBVqH4U9HoBTvwV55kxtJbGkYvEijkXu16QzMxMl52dHbP7SZQVFkD2aJjxIOTthA5/hnP6Q5UaflcmklDMbIFz7tDJGgfRTFEpnx/mepODNnzlLWnbczjU/4XfVYkkNQW6lM329TDlHvjyLajZBK56BY6/BMz8rkwk6SnQJTyF+TD/OZg5DArz4Ow74Oy/Q+VqpV8rIjGhQJfSLf/I617ZtBRadYMew+DIln5XJSIHUaBLybau8Tab+Hoc1D4a+rwJrXv6XZWIlECBLocq2AufjIRZj4ErgnPvhDP/AmlV/a5MRA5DgS4H+m6qt3PQ5u/huIug+8NQ52i/qxKRMCjQxfPTSpg0GJZ+CHVbwq/fhWO7+l2ViJSBAj3Z5e+GOU/Ax0+ApULXId4EoUpVAO0YJBIkCvRk5RwsnQCTBsKWVXBCL+j2INT6/7DWjkEiwaJAj3NRaSH/+D1MHADLpkD94+D68XDMOYecph2DRIJFgR7HIt5CztsJsx+HuU9DahXvA8/T+0JqWrGna8cgkWCJ+gYXUn6HayGXiXOw+D0YeZoX6Cf0glsXwBk3lxjmUPK65doxSCQ+KdDjWERayLlL4ZVL4T83Qnpd+M0k6PUvqNGw1Eu1Y5BIsKjLJY5VaE/NPdvgo0e89VcqV/NWQ8y8CVLD/yPXjkEiwaJAj2P9urc+oA8dwmghOweL/gOT74Yd66Hddd5QxGr1ylWDdgwSCQ4Fehwrcwt5/VfeIlqr5kLjdtD7dWjSPoYVi4ifFOhxLqwW8u4tMONh+Ox5qFobLn4S2l0PKSmaGCSSRBToQVZUBP97HabcC7s3e33k594JR9QFNDFIJNko0IMq53OveyUnG5qcDheOhUYnH3CKJgaJJBcFetDs2gzT7oMFY6BafbjsOTjpakg5dASqJgaJJBcFelAUFcKCl2H6A96QxA5/gs4DoWqtEi+p0LBHEQkcTSwKgtWfwvPnwod/hwYnwB/nQI+hhw1z0MQgkWSjFno827ERpg6BL16DGo3giheh7RVgFtblmhgkklwU6PGosAA+e8Ebipi/CzreBp36QZXqZf5RmhgkkjwU6PFm5cfe6JWNi6HledDzUajXyu+qRCQAFOjxYts6mHK3N22/VlO4+lVvT88wu1dERCoU6Ga2EtgOFAIFzrnMSBSVVAryvAW0PnoECvOhU384629Q+Qi/KxORgIlEC/1c59ymCPyc5PP9DJjYHzZ9C7/oCT0ehrrH+F2ViASUulz8sGU1ZA2GJeOhTgu45m34RXe/qxKRgKtooDtgspk54F/OuVERqClxFeyFuU/BrMe95+feBWfeCmlV/a1LRBJCRQO9o3NurZk1AKaY2TfOuVn7n2BmfYG+AM2aNavg7QLs28kwaQBsXg7HXwLdH4LaSfz7ISIRV6GZos65taGvG4H3gNOLOWeUcy7TOZdZv379itwumDavgDf6wOtXgqXCde/B1f9WmItIxJU70M2smpnV2PcY6AZ8FanCAi9vF8x4mMKRp7Nr6XSG5ffhnB0PMW6bpt2LSHRUpMulIfCeeeOkKwGvO+cmRaSqIHMOvvkQJg2CrauYWNSRB/L6sIG6sLVA65GLSNSUO9Cdc8uBk0s9MZlsWuYNQ/x+GjRow82VH+DDbS0POEXrkYtItGi1xUjYu8NbROvZDrDmM+jxCPxhNhMOCvN9tB65iESDxqFXhHOweCxk3QXb18Ip10LXIVC9AaD1yEUkttRCL6+NS2DMxfDOTVCtHtw0GS579ucwB61HLiKxlfAt9Ijver9nG8wc5q2/UqUGXDgC2t8IKamHnKr1yEUklhI60CO6671z8OVbMPlu2JkL7W+A8+6Bakce9jKtRy4isZLQgR6xXe/XfemtUb56HmRkwjVvQcapEa5WRKRiEjrQK7zr/e6fYPpDkP0ipNeBS0Z6H3ym6KMHEYk/CZ1MJY0mKXWUSVERLBgDT7fHffYi/0npycmbh9Jxcgbj/rcuCpWKiFRcQgd6uUaZ5CyAF7rAf//CpqpHc3nhUPrtvJatVP+5D37cwpwoVy4iUnYJHeiXtctgaK8TyaidjgEZtdMZ2uvE4vvPd/4I42+F57vAthzo9TyX7ryLL/KbHnDavj54EZF4k9B96BDGKJOiQsgeDdMfhLwdcMbNcM4AqFqTta9/WOwlmukpIvEo4QP9sFbNgwl3wPpF0KIT9BwODY77+dua6SkiQZLQXS4l2r4B3vsjjO4OuzbDlS/D9eMPCHPQTE8RCZbkaqEX5sOno2DGUCjYA2f9HTrdAZWrFXu6ZnqKSJAkT6CvmAUT+kPuEji2q7ciYr1jS71MMz1FJCgSP9C35sDku7xVEWs3g96vQ+sLwNuYQ0QkYSRuoBfkwbxn4KPh4Aqh8yDo+FdI0weaIpKYEjPQl02FiQPgx2XQ+kLo8TDUae53VSIiUZVYgf7TD5A1GL75AOq2hGvfgVbn+12ViEhMJEag5++Bj5+EOSPAUqDLPXDGLVCpit+ViYjETPADfelEmDQQfloJJ1wO3R6EWk38rkpEJOaCG+g/fg+TBsF3WVCvNVz/PhzT2e+qRER8E7xAz9sJs0fA3KcgtQp0ewh++QdITfO7MhERXwUn0J2DJeNh0mDYtgZOuhrOvx9qHOV3ZSIicSEYgZ67FCb2h+UzoWFbuOJ5OPpMv6sSEYkrwQj0T0bC2oXeaoiZN0FqMMoWEYmlYCRj1/vgvHugen2/KxERiVsVWj7XzHqY2VIzW2ZmAyNV1CGOqKswFxEpRbkD3cxSgWeAnkAboI+ZtYlUYSIiUjYVaaGfDixzzi13zuUBbwKXRqYsEREpq4oEegawer/na0LHRETEBxUJ9OIWFHeHnGTW18yyzSw7Nze3ArcTEZHDqUigrwGa7ve8CbD24JOcc6Occ5nOucz69fXBpohItFQk0D8DWplZCzOrDPQGxkemLBERKatyj0N3zhWY2S1AFpAKjHbOLY5YZSIiUiYVmljknJsATIhQLSIiUgEVmlgkIiLxQ4EuIpIgFOgiIgki7hfnGrcwh+FZS1m7ZTeNa6fTr3trLmun+UsiIgeL60AftzCHQWMXsTu/EICcLbsZNHYRgEJdROQgcd3lMjxr6c9hvs/u/EKGZy31qSIRkfgV14G+dsvuMh0XEUlmcR3ojWunl+m4iEgyi+tA79e9NelpqQccS09LpV/31j5VJCISv+L6Q9F9H3xqlIuISOniOtDBC3UFuIhI6eK6y0VERMKnQBcRSRAKdBGRBKFAFxFJEAp0EZEEYc4dsq9z9G5mlgv8UM7L6wGbIlhOEOg1Jwe95uRQkdd8tHOu1E2ZYxroFWFm2c65TL/riCW95uSg15wcYvGa1eUiIpIgFOgiIgkiSIE+yu8CfKDXnBz0mpND1F9zYPrQRUTk8ILUQhcRkcMIRKCbWQ8zW2pmy8xsoN/1RJuZNTWzGWa2xMwWm9lf/a4pFsws1cwWmtkHftcSC2ZW28zeMbNvQn/WZ/hdU7SZ2d9Cf6e/MrM3zKyq3zVFmpmNNrONZvbVfsfqmtkUM/su9LVONO4d94FuZqnAM0BPoA3Qx8za+FtV1BUAtzvnjgc6ADcnwWsG+CuwxO8iYuhJYJJz7jjgZBL8tZtZBvAXINM51xZIBXr7W1VUvAz0OOjYQGCac64VMC30POLiPtCB04Flzrnlzrk84E3gUp9riirn3Drn3Oehx9vx/qEn9BrCZtYEuBB4we9aYsHMagKdgBcBnHN5zrkt/lYVE5WAdDOrBBwBrPW5nohzzs0CNh90+FJgTOjxGOCyaNw7CIGeAaze7/kaEjzc9mdmzYF2wHx/K4m6J4D+QJHfhcTIMUAu8FKom+kFM6vmd1HR5JzLAR4DVgHrgK3Oucn+VhUzDZ1z68BrsAENonGTIAS6FXMsKYbmmFl14F3gNufcNr/riRYzuwjY6Jxb4HctMVQJOBX4p3OuHbCTKL0NjxehfuNLgRZAY6Camf3a36oSSxACfQ3QdL/nTUjAt2kHM7M0vDB/zTk31u96oqwjcImZrcTrUjvPzF71t6SoWwOscc7te+f1Dl7AJ7KuwArnXK5zLh8YC5zpc02xssHMGgGEvm6Mxk2CEOifAa3MrIWZVcb7EGW8zzVFlZkZXt/qEufcCL/riTbn3CDnXBPnXHO8P9/pzrmEbrk559YDq81s347nXYCvfSwpFlYBHczsiNDf8S4k+AfB+xkP3BB6fAPwfjRuEvd7ijrnCszsFiAL71Px0c65xT6XFW0dgeuARWb2RejYYOfcBB9rksi7FXgt1FBZDvzG53qiyjk338zeAT7HG8m1kAScMWpmbwCdgXpmtga4FxgGvG1mv8X7j+3KqNxbM0VFRBJDELpcREQkDAp0EZEEoUAXEUkQCnQRkQShQBcRSRAKdBGRBKFAFxFJEAp0EZEE8X8vgUjs7Y2mDQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(x,y, 'o')\n", "plt.plot(x, alpha_hat + beta_hat * x)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "----\n", "Notice that we didn't get exactly the correct answer. The points were generated with a slope of 2.5 and an intercept of 1, whereas our fit was a little bit off" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Hypothesis Test for Normality\n", "====\n", "\n", "Another important hypothesis is determining if a sample is normally distributed. We said in class many times that we require things to be normal for many methods. Here's how we actually test that." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Shapiro–Wilk Test\n", "====\n", "\n", "**Data Type:** single group of samples\n", "\n", "**Compares:** If the samples came from an unknown parent normal distribution (are they normally distributed)\n", "\n", "**Null Hypothesis:** The samples are from the unknown parent normal distribution\n", "\n", "**Conditions:** None\n", "\n", "**Python:** `scipy.stats.shapiro`\n", "\n", "**Notes:** There are many other tests for normality. This one is not the simplests, but is the most effective" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "(0.9825882911682129, 0.9762558937072754)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = [12.4, 12.6, 11.8, 11.5, 11.9, 12.2, 12.0, 12.1, 11.8]\n", "\n", "scipy.stats.shapiro(data)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The $p$-value is quite high, so we don't reject the null hypothesis" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "(0.9547258019447327, 0.0017220161389559507)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = np.linspace(0,1, 100)\n", "\n", "scipy.stats.shapiro(data)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The $p$-value is 0.002, so we correctly reject the null hypothesis" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Linear Regression Assumption Check\n", "=====\n", "\n", "One of our assumptions is that the noise is normally distributed. Recall our model:\n", "\n", "$$y = \\alpha + \\beta x + \\epsilon$$\n", "$$\\epsilon = y - \\alpha - \\beta x \\approx y - \\hat{\\alpha} - \\hat{\\beta} x = y - \\hat{y}$$\n", "\n", "Where the $\\approx$ is because we are using our estimates for $\\alpha$ and $\\beta$. We can now check our assumption by histogramming the residuals, which should be the same as looking at the $\\epsilon$ distribution." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADWdJREFUeJzt3X+I34V9x/HXazHij1oczXdbZ7ydgyIrro1yiF2gbNGW2EjKxgqRWbqtcP+0nY6OLk6YlDHI6Og62Ng4rOtAa2mtocP0hylVpNCmvcRok56WzqUatctJceoKdbGv/XHftJfz+73v5/Q+38/3nT4fcJhv7uN9X0TvyTef+3y/XycRAKCOX+p6AABgbQg3ABRDuAGgGMINAMUQbgAohnADQDGEGwCKIdwAUAzhBoBizmrji27atCnT09NtfGkAOCMdPHjw2SS9Jse2Eu7p6WnNz8+38aUB4Ixk+wdNj+VUCQAUQ7gBoBjCDQDFEG4AKIZwA0AxjcJt+89tH7V9xPZdts9pexgAYLCR4bZ9kaQ/kzST5DJJGyTtansYAGCwpqdKzpJ0ru2zJJ0n6en2JgEAVjMy3EmekvT3kp6Q9Iyk/0lyX9vDAACDjXzmpO1flvRuSZdIek7S52zfkOSOFcfNSpqVpKmpqRamYr1N797X2X0f27Ojs/sGqmtyquQaSf+VZDHJ/0m6R9LvrDwoyVySmSQzvV6jp9sDAF6FJuF+QtJVts+zbUlXS1podxYAYJgm57gPSLpb0iFJ3+n/O3Mt7wIADNHo1QGT3Crp1pa3AAAa4JmTAFAM4QaAYgg3ABRDuAGgGMINAMUQbgAohnADQDGEGwCKIdwAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACiGcANAMYQbAIoh3ABQzMhw277U9uFlH8/bvmkc4wAArzTyrcuSPCZpiyTZ3iDpKUl7W94FABhiradKrpb0n0l+0MYYAMBoaw33Lkl3tTEEANBM43DbPlvSTkmfG/L5WdvztucXFxfXax8AYIW1POK+VtKhJP896JNJ5pLMJJnp9Xrrsw4A8AprCff14jQJAHSuUbhtnyfpHZLuaXcOAGCUkZcDSlKSH0t6Q8tbAAAN8MxJACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMUQbgAopulbl11o+27bj9pesP22tocBAAZr9NZlkv5R0peT/KHtsyWd1+ImAMAqRobb9uslvV3SH0tSkpckvdTuLADAME1OlfympEVJ/2b7Idu32T5/5UG2Z23P255fXFxc96EAgCVNwn2WpCsk/UuSyyX9r6TdKw9KMpdkJslMr9db55kAgFOahPu4pONJDvRv362lkAMAOjAy3El+KOlJ25f2f+tqSd9tdRUAYKimV5V8SNKd/StKHpf0J+1NAgCsplG4kxyWNNPyFgBAAzxzEgCKIdwAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAimn0Dji2j0l6QdLLkk4m4d1wAKAjTd9zUpJ+L8mzrS0BADTCqRIAKKZpuCPpPtsHbc8OOsD2rO152/OLi4vrtxAAcJqm4d6a5ApJ10r6gO23rzwgyVySmSQzvV5vXUcCAH6uUbiTPN3/5wlJeyVd2eYoAMBwI8Nt+3zbF5z6taR3SjrS9jAAwGBNrir5VUl7bZ86/tNJvtzqKgDAUCPDneRxSW8dwxYAQANcDggAxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMU0DrftDbYfsn1vm4MAAKtbyyPuGyUttDUEANBMo3Db3ixph6Tb2p0DABilybu8S9InJH1E0gXDDrA9K2lWkqampl77MpzRpnfv6+R+j+3Z0cn9Autp5CNu29dJOpHk4GrHJZlLMpNkptfrrdtAAMDpmpwq2Sppp+1jkj4jaZvtO1pdBQAYamS4k9ycZHOSaUm7JH0tyQ2tLwMADMR13ABQTNMfTkqSkjwg6YFWlgAAGuERNwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMUQbgAohnADQDGEGwCKIdwAUAzhBoBimrxZ8Dm2v2X7YdtHbX90HMMAAIM1eQecn0jaluRF2xslfd32l5J8s+VtAIABRoY7SSS92L+5sf+RNkcBAIZrdI7b9gbbhyWdkLQ/yYF2ZwEAhmn0ZsFJXpa0xfaFkvbavizJkeXH2J6VNCtJU1NT6z4UqGx6977O7vvYnh2d3TfasaarSpI8p6V3ed8+4HNzSWaSzPR6vXWaBwBYqclVJb3+I23ZPlfSNZIebXsYAGCwJqdK3ijp321v0FLoP5vk3nZnAQCGaXJVySOSLh/DFgBAAzxzEgCKIdwAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAimnynpMX277f9oLto7ZvHMcwAMBgTd5z8qSkDyc5ZPsCSQdt70/y3Za3AQAGGPmIO8kzSQ71f/2CpAVJF7U9DAAw2JrOcdue1tIbBx9oYwwAYLTG4bb9Okmfl3RTkucHfH7W9rzt+cXFxfXcCABYplG4bW/UUrTvTHLPoGOSzCWZSTLT6/XWcyMAYJkmV5VY0iclLST5ePuTAACrafKIe6uk90raZvtw/+NdLe8CAAwx8nLAJF+X5DFsAQA0wDMnAaAYwg0AxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgmCbvOXm77RO2j4xjEABgdU0ecX9K0vaWdwAAGhoZ7iQPSvrRGLYAABrgHDcAFDPyXd6bsj0raVaSpqamXvXXmd69b70mlXFsz46uJ+AM9ov4PdWVcX0vr9sj7iRzSWaSzPR6vfX6sgCAFThVAgDFNLkc8C5J35B0qe3jtt/f/iwAwDAjz3EnuX4cQwAAzXCqBACKIdwAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAimkUbtvbbT9m+/u2d7c9CgAwXJP3nNwg6Z8lXSvpzZKut/3mtocBAAZr8oj7SknfT/J4kpckfUbSu9udBQAYpkm4L5L05LLbx/u/BwDowMh3eZfkAb+XVxxkz0qa7d980fZjyz69SdKza5/XuonY5b877eZEbBrgjNi14s+6LWfEn9UYTeKuV7XpNf7/9RtND2wS7uOSLl52e7Okp1celGRO0tygL2B7PslM01HjMom7JnGTxK61mMRNErvWYhI3LdfkVMm3Jb3J9iW2z5a0S9J/tDsLADDMyEfcSU7a/qCkr0jaIOn2JEdbXwYAGKjJqRIl+aKkL76G+xl4CmUCTOKuSdwksWstJnGTxK61mMRNP+PkFT9nBABMMJ7yDgDFjD3ctv/CdmxvGvd9D9jyN7YfsX3Y9n22f73rTZJk+2O2H+1v22v7wq43SZLt99g+avuntjv9ifskvgyD7dttn7B9pOsty9m+2Pb9thf6//1unIBN59j+lu2H+5s+2vWm5WxvsP2Q7Xu73jLIWMNt+2JJ75D0xDjvdxUfS/KWJFsk3Svpr7se1Ldf0mVJ3iLpe5Ju7njPKUck/YGkB7scMcEvw/ApSdu7HjHASUkfTvJbkq6S9IEJ+PP6iaRtSd4qaYuk7bav6njTcjdKWuh6xDDjfsT9D5I+ogFP4OlCkueX3Txfk7PrviQn+ze/qaVr5zuXZCHJY6OPbN1EvgxDkgcl/ajrHSsleSbJof6vX9BSkDp99nOWvNi/ubH/MRHff7Y3S9oh6bautwwztnDb3inpqSQPj+s+m7D9t7aflPRHmpxH3Mv9qaQvdT1iwvAyDK+S7WlJl0s60O2Sn52OOCzphKT9STrf1PcJLT3A/GnXQ4ZpdDlgU7a/KunXBnzqFkl/Jemd63l/Tay2KckXktwi6RbbN0v6oKRbJ2FX/5hbtPTX3DvHsanprgnQ6GUYcDrbr5P0eUk3rfjbZieSvCxpS/9nOHttX5ak058P2L5O0okkB23/bpdbVrOu4U5yzaDft/3bki6R9LBtaemv/odsX5nkh+u5oemmAT4taZ/GFO5Ru2y/T9J1kq7OGK/ZXMOfV5cavQwDfs72Ri1F+84k93S9Z7kkz9l+QEs/H+j6B7tbJe20/S5J50h6ve07ktzQ8a7TjOVUSZLvJPmVJNNJprX0jXdF29Eexfablt3cKenRrrYsZ3u7pL+UtDPJj7veM4F4GYY18NKjpU9KWkjy8a73SJLt3qmrpWyfK+kaTcD3X5Kbk2zud2qXpK9NWrQlruPeY/uI7Ue0dBqn88uk+v5J0gWS9vcvVfzXrgdJku3ft31c0tsk7bP9lS529H9we+plGBYkfXYSXobB9l2SviHpUtvHbb+/6019WyW9V9K2/v9Ph/uPKLv0Rkn397/3vq2lc9wTeendJOKZkwBQzC/6I24AKIdwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMX8P4m6RMe5lwB/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.hist(y - beta_hat * x - alpha_hat)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "----\n", "\n", "At this point, it's unclear if they are normally distributed. Luckily we just learned the Shapiro–Wilk Test!" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "(0.9317496418952942, 0.16685707867145538)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scipy.stats.shapiro(y - beta_hat * x - alpha_hat)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "----\n", "It looks like the residuals may indeed be normally distributed. So, our assumption was valid." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Linear Regression Goodness of Fit\n", "====\n", "\n", "There are a few ways of measuring goodness of fit. One way is to just compute the SSR, the sum of the squared residuals. However, this has the negative that the units of $y$ appear in the goodness of fit. Here is what people typically use, the coefficient of determination:\n", "\n", "$$R^2 = 1 - \\frac{\\textrm{SSR}}{\\textrm{TSS}} = 1 - \\frac{\\sum_i \\left(\\hat{y}_i - y\\right)^2}{\\sum_i \\left(\\bar{y} - y\\right)^2}$$\n", "\n", "This equation has the property that it's unitless, it's $1$ when the fit is perfect, and $0$ when the fit is awful. In the case of linear regression, $R$ is the same as the correlation coefficient." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.9454803725748993 0.9723581503617374\n", "[[1. 0.97235815]\n", " [0.97235815 1. ]]\n" ] } ], "source": [ "ssr = np.sum((y - alpha_hat - beta_hat * x)**2)\n", "tss = np.sum((np.mean(y) - y)**2)\n", "\n", "rsq = 1 - ssr / tss\n", "\n", "print(rsq, sqrt(rsq))\n", "print(np.corrcoef(x,y))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Linear Regression Error Analysis\n", "=====\n", "\n", "Knowing goodness of fit is not quite useful yet though. We'd like to do hypothesis tests, build confidence intervals, etc. We need to know how the $\\hat{\\beta}$ are distributed!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Standard Error\n", "---" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "As you may have noticed, sometimes it gets confusing that the difference between the sample mean and the true mean, $\\mu - \\bar{x}$, are distributed according to a normal or $T$ distribution if the samples are distributed according to a normal distribution. Thus we end up with two distributions: one with parameters $\\mu$, $\\sigma$ that describes the samples and one with parameters $\\mu$, $\\sigma / \\sqrt{N}$ that describes the difference between sample mean and true mean. Often this $\\sigma / \\sqrt{N}$ term is called **Standard Error** to distinguish it from the standard deviation of the population (true/hidden) distribution. Thus you can write:\n", "\n", "$$T = \\frac{\\bar{x} - \\mu}{\\sigma_x / \\sqrt{N}} = \\frac{\\bar{x} - \\mu}{S} $$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We're to be referring to **Standard Error** often now. You then can use **Standard Error** in confidence intervals or hypothesis tests. The same rules as previously apply: $N < 25$ requires a $t$-distribution and above is normal." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "One other important thing to remember is that the denominator in the standard error should be the square root of the degrees of freedom. Usually this is written as $N - D$, where $D$ is the deducted degrees of freedom. In the case of linear regression, we have $D$ being the number of coefficients we're fitting. That $N - D$ is also the degrees of freedom in the $t$-distribution. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Error Equation for OLS-1D\n", "---\n", "\n", "You may read the derivation in Bulmer on page 226. The variance in our estimated values ($S^2$) for slope and intercept are:\n", "\n", "$$S^2_{\\epsilon} =\\frac{\\sigma^2_{\\epsilon}}{N - D} = \\frac{1}{N - D}\\sum_i \\left(\\hat{y}_i - y_i\\right)^2$$\n", "\n", "$$S^2_{\\alpha} = S^2_{\\epsilon} \\left[ \\frac{1}{N - D} + \\frac{\\bar{x}^2}{\\sum_i\\left(x_i - \\bar{x}\\right)^2}\\right]$$\n", "\n", "$$S^2_{\\beta} = \\frac{S^2_{\\epsilon}}{\\sum_i \\left(x_i - \\bar{x}\\right)^2}$$\n", "\n", "$D$ here is the number of fit coefficients. 2 in our case." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Example: Compute the standard error of the intercept and create a probability distribution for the true intercept\n", "----\n", "\n", "We'll continue using the data from above" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The standard error for the intercept is 0.8439285616547768\n" ] } ], "source": [ "df = len(x) - 2\n", "s2_epsilon = np.sum((y - alpha_hat - beta_hat * x) ** 2) / df\n", "s2_alpha = s2_epsilon * (1. / df + np.mean(x) ** 2 / (np.sum((np.mean(x) - x) ** 2)))\n", "\n", "print('The standard error for the intercept is', np.sqrt(s2_alpha))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's just visualize now what the distribution for where the true intercept looks like. We know it is distributed according to:\n", "\n", "$$P(\\alpha) = T(\\mu=\\hat{\\alpha}, \\sigma=S_\\alpha, df=18)$$\n", "\n", "The degrees of freedom here is " ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VFX+BvD3m15IhYRUSCCdJi2CosAiCohiXUXXFcWCytp+qGBdl1WxrLuKKOu6uuqq2BWRIiIIqAihk0YKhISQkALpfc7vj4xuhEAmYWbOlPfzPHkymbmZ++Yyebm5c+65opQCERE5FhfdAYiIyPxY7kREDojlTkTkgFjuREQOiOVOROSAWO5ERA6I5U5E5IBY7kREDojlTkTkgNx0rbhPnz4qJiZG1+qJyCi/rA4AMCDEt2dPkJ3d/jkx0UyJ6HS2b99erpQK6Wo5beUeExODtLQ0XasnIqNr/vkTAODD28f27AkmTGj/vGGDWfLQ6YlIgSnL8bAMEZEDYrkTETkgljsRkQNiuRMROSCWOxGRA2K5ExE5IJY7EZED0jbOnYj0aGptQ+aRGpRUNaKsphGHjzXAxQVYtvUQArzdkRTuj5jePhAR3VHpDLDciZxAQUUdlu8qxk/5FdhecAxNrYaTlpn/2d5fbwf6uOOs6EBMTumL6UMjEODtbs24ZAYsdyIH1WZQWJtRivd+LsCmnHKIAMlh/rj+7P5IjQ1CVJAPQv08MfeDnTAohZevHY7KumakF1dh56Hj2HqgEo98vg9PfpWBC1P64uZxsRjRL0j3j0UmYrkTOaAN2UfxzMosZJfWIDzAC/dPTsA1o6PR19/rpGUFgKsIIgK9ERHojcGRAbhmdD8opbD3cBU+3V6EL3cXY8WeI7gguS/mXZSApDB/6/9Q1C0sdyIHkl9WiyeWp2NTTjn69/bBK9cNx5RBYXBz7f7YCRHB0KhADI0KxENTk/DWDwex9Ps8TH1pE64dHY2HpyXDz4uHa2wVy53IASil8N7Ph/DXrzPg4eqCRy9Oxh/HxsDDzTwD4nw83HDXxDj84ez+eGV9Dv69+QA27i/Hc1cNxblmWQOZG8udyM6V1zbhoU/2YF3WUZwX3wcvXD2s08Mv5hDg445HLk7B1CHhmPfxblz/xs/YUFGP/sE+4Nga28JyJ7JjWSXVmP2fNJTVNuGJS1Jw49gYuLhYvmZH9AvCyrvPw6JVWSh5vwH1za0Ir2tGkK+HxddNpuFJTER2an3WUVz56o9oNRjw6ZxzcNO5sVYp9l94ubviz5cOwoCQXqhpbMWMJT8gu6TGauun02O5E9mh/24pwOy3tyGmjy++vGschkQFaMsS6ueJlAh/NLa04eqlP2LHoWPastD/sNyJ7My/Nx/Ao1/sw8TEUHw8ZyzCAixzfL07/Dzd8Nmd5yDI1wN/eONn/JRXoTuS02O5E9mRf36fh4UrMjB1cBiW3jASPh6287ZZVJAPPr59LCIDvTHrra1Yn3VUdySnxnInshNLv8/DM6uyMH1oOF6eORzuPRi7bmmh/l748PaxiO/bC7f/dzt+zCvXHclp2d6rg4hO8tG2QixalYVLhkXgH9ecZZPF/otgXw+8e/PZ6B/sg1vfTsOeouO6Izkl232FEBEA4NuMUiz4fC/OTwjB364e1qOzTa0tyNcD784+G0G+Hpj11jbkHq3VHcnp2P6rhMiJbS+oxF3v78DgCH+8dv0Is51xag1hAV54d/bZcBHgj//+GUdrGnVHcir280ohcjKFlfW49Z3tiAj0xpuzRsPX03bePDVVbB9fvDUrFZX1zbj93e1obGnTHclpsNyJbFBdUytufScNrW0GvDlrNHr38tQdqceGRAXgxd+fhZ2HjmPBZ3uhlNIdySmw3IlsjFIK8z7ejf2lNVh83QjE9vHVHemMTRsSjvsnJ+DznYfx2vd5uuM4BZY7kY155btcrNpXggVTkzE+IUR3HLP50+/icMmwCDy/Jhubcsp0x3F4LHciG7I5pxwvfrsflw+PxC3nxeqOY1YigmevHIL40F64d9kulFTxDVZLYrkT2Yij1Y2498OdiAvphacuH+yQF6j28XDDq9ePQENLG+a+vwMtbSdfy5XMg+VOZAPaDAr3LNuF2qZWLLl+hE1NK2BucaF+eOaKIUgrOIbn12TrjuOwWO5ENuDldTn4Kb8CC2cMRkJfP91xLG7GWZG4/ux+eH1jPr7fz+PvlsByJ9Is7WAlFn+XgytGROLqUdG641jNY9NTkNC3F+Z9vBsVtU264zgcljuRRrVNrbj/o92IDPLGX2YM1h3HqrzcXfGPa4ajqr4F8zn+3exMKncRmSIi2SKSKyLzT7PcaBFpE5GrzBeRyHEt/CoDRcfq8fffn4VedngG6plKifDHAxclYm1GKT7cVqg7jkPpstxFxBXAEgBTAaQAmCkiKadY7lkAa8wdksgRfZNegg/TCjFn/ECMignWHUeb2eNicW5cbzz5VQYKKup0x3EYpuy5pwLIVUrlK6WaASwDMKOT5f4E4FMAnKGfqAsVtU1Y8NlepIT7494LEnTH0crFRfDC1cPg5iJ48JM9MBh4eMYcTCn3SAAd/14qMt73KxGJBHA5gKXmi0bkuJ78KgPVjS148ZphdjXTo6WEB3jjsekp+PlAJf77c4HuOA7BlFdVZ2dSnPhf6z8APKSUOu2UbyJym4ikiUhaWRmHP5FzWptRiuW7izF3YjySwvx1x7EZV4+KwvkJIVi0KguFlfW649g9U8q9CEDH8VlRAIpPWGYUgGUichDAVQBeFZHLTnwipdTrSqlRSqlRISGOM2cGkamqGlrwyOd7kRTmhzsmDNQdx6aICJ65YghchIdnzMGUct8GIF5EYkXEA8C1AJZ3XEApFauUilFKxQD4BMCdSqkvzJ6WyM499XUGKuqa8fxVPBzTmchAbzxycTJ+yq/AR2kcPXMmunx1KaVaAcxF+yiYTAAfKaXSRWSOiMyxdEAiR/FjXjk+SivCbecPwJCoAN1xbNa1o6ORGhuMp1dmoqyGJzf1lEm7DkqplUqpBKXUQKXUU8b7liqlTnoDVSk1Syn1ibmDEtmzxpY2PPr5PvTv7YN7JsXrjmPTRARPXz4EjS0GLFyRoTuO3eLfhURW8NqGPOSX1+Gvlw2Gl7ur7jg2Ly60F+6YMBDLdxdjQzZHV/cEy53IwvLKavHahjxcOiwC58VzIIGp7pw4EANCfPHYl/vQ0Mxrr3YXy53IgpRSeOTzvfByd8Gj05N1x7Ernm6uePryISisbMCS9bm649gdljuRBX25qxhb8ivx0NQkhPp56Y5jd8YM6I0rhkfi9Y35yC+r1R3HrrDciSykprEFT63MxLCoAMwc3U93HLs1f1oSPN1c8MTydM4c2Q0sdyILeXldDsprm/CXGYPh4uJ4l8yzllA/L9x/YQI25ZRj9b4S3XHsBsudyAJySmvw1g8Hcc2oaAyLDtQdx+7dMKY/ksP98ZcVGahvbtUdxy6w3InMTCmFx79Mh6+nGx6ckqQ7jkNwc3XBwhmDcKSqkW+umojlTmRmK/eW4Kf8Csy7KBHBvh664ziMUTHBuHx4JP616QAOVXBisa6w3InMqLGlDU+vzERSmB+uS+WbqOb20JQkuLkInlrJM1e7wnInMqN/bczH4eMNeOKSQXDlm6hmFxbghbsmxmFNeil+yC3XHcemsdyJzORIVQNe3ZCHqYPDMHZgb91xHNbscbGIDvbGk1+lo7XNoDuOzWK5E5nJs6uy0KYUHp7GM1EtycvdFY9MS8H+0lq8v/WQ7jg2i+VOZAY7Dh3DF7uKcdt5AxAd7KM7jsO7aFBfnDOwN/6+dj+q6lt0x7FJLHeiM6SUwsIVGQjx8+TVlaxERPDIxck43tCCV9bn6I5jk1juRGdoxZ4j2HnoOB64MBG+nm664ziNQREBuHpkFP7z40EcLK/THcfmsNyJzkBjSxsWrcpCcrg/rhwZpTuO05l3YSLcXV2waFWW7ig2h+VOdAbe+uEgDh9vwKMXJ3Poowah/l64Y/xArE4vwZb8Ct1xbArLnaiHymub8Or6XExKCsW5cX10x3Fat54/AOEBXnh6ZSYMBs4a+QuWO1EPvfRtDupb2rCAQx+18nJ3xbwLE7GnqApf7SnWHcdmsNyJeiCvrH2M9czUaMSF9tIdx+ldPjwSKeH+eG51NhpbeEk+gOVO1CPPrc6Cl5sL7pmUoDsKAXBxaR8aefh4A9756aDuODaB5U7UTdsOVmJNeinmjB+IED9P3XHI6Ny4PpiQGILF3+XiWF2z7jjasdyJukEphadXZqKvvyduOW+A7jh0ggVTk1HX1IrF33HOd5Y7UTes2leCnYeO4/7JCfD2cNUdh06QGOaHq0dG490tB1FY6dxzvrPciUzU0mbAc6uzkNC3F64aGa07Dp3CfZMT4OoieOGbbN1RtGK5E5lo2dZDOFhRj/lTk3jCkg0LC/DC7HGx+HJXMfYWVemOow3LncgEtU2teGldDs6ODcbExFDdcagLt48fiCAfdyxanQmlnPPEJpY7kQn+tTEf5bXNWDAtGSLca7d1/l7u+NPv4vFDbgU25jjnFZtY7kRdOFrTiH9tyse0IWE4KzpQdxwy0fVj+iE62BuLVmU55bQELHeiLixel4vmVgMeuChJdxTqBk+39mkJMo9UY/lu55uWgOVOdBoHyuvwwdZDmJnaD7F9fHXHoW66ZGgEBkX444VvstHU6lzTErDciU7jhW+y4eHmgj9NitMdhXrAxUUwf2oSio414L0tznW9VZY70SnsKTqOr/ccwS3jYhHq56U7DvXQefEhODeuNxZ/l4PqRue53qpJ5S4iU0QkW0RyRWR+J4/PEJE9IrJLRNJEZJz5oxJZj1IKi1ZlIdjXA7eez2kG7N1DU5JwrL4F/9qYrzuK1XRZ7iLiCmAJgKkAUgDMFJGUExZbB2CYUuosADcDeMPcQYmsaVNOOX7Mq8DciXHw83LXHYfO0NCoQEwfGo43Nh3A0ZpG3XGswpQ991QAuUqpfKVUM4BlAGZ0XEApVav+d6aALwDnG3dEDsNgUHh2dRaigrxx/Zh+uuOQmcy7MBEtbQYsXucck4qZUu6RAAo7fF1kvO83RORyEckC8DXa995PIiK3GQ/bpJWVlfUkL5HFrdh7BOnF1fi/CxPg6cbJwRxFTB9fXJsajQ+2HsLB8jrdcSzOlHLv7HS8k/bMlVKfK6WSAFwGYGFnT6SUel0pNUopNSokJKR7SYmsoLnVgL99k42kMD/MGHbSPgzZubsnxcPd1cUpJhUzpdyLAHScAi8KwCnPCFBKbQQwUER4xWCyOx9uO4SCino8NCUJLpwczOGE+rVPKrZizxGHn1TMlHLfBiBeRGJFxAPAtQCWd1xAROLEOOGGiIwA4AGgwtxhiSyprqkVL63LRWpsMCYk8i9LR3Xb+AEI8nHHc2uydEexqC7LXSnVCmAugDUAMgF8pJRKF5E5IjLHuNiVAPaJyC60j6y5RjnrVGxkt/69+QDKa5swf2oSJwdzYP5e7rhrYhw25ZTjh1zHnVTMpHHuSqmVSqkEpdRApdRTxvuWKqWWGm8/q5QapJQ6Syk1Vim12ZKhicytorYJr2/Mx0WD+mJEvyDdccjC/jCmPyIDHXtSMZ6hSgRgyfo81De34oGLEnVHISvwcnfFfZMTsPdwFVbuO6I7jkWw3MnpFVbW479bCnD1yGjEhfrpjkNWcvnwSCT29cMLa7LR0mbQHcfsWO7k9P6+dj9EgHsnx+uOQlbk6iJ4cEoiDlbU48NthV1/g51huZNTyzxSjc93Hcasc2MQHuCtOw5Z2e+SQjE6JggvrctBfXOr7jhmxXInp/b8mmz4ebrhzvGc0tcZibRPCVxW04Q3Nx/QHcesWO7ktLbkV+C7rKO4Y0IcAnw4OZizGtk/GBem9MXS7/NRWdesO47ZsNzJKf0ypW+YvxduOjdGdxzS7MEpiahvbsUr3znOpGIsd3JKq/eVYFfhcdw/OQFe7pwczNnFhfrh96Oi8e6WgyisrNcdxyxY7uR0WtoMeH5NNuJDe+GKEZwcjNrde0ECXETwNweZVIzlTk7nw22FyC+vw0NTkuDmyl8BahcW4IWbzo3Fl7uLkV5s/5OK8ZVNTqV9crAcjI4JwqTkUN1xyMbcMWEgArzdsWiV/U8qxnInp/LGpgMoq2nC/KnJnByMThLg7Y65xknFNufY96RiLHdyGmU1TXh9Yx6mDArDyP6cHIw6d8PY9knFnlmVadeTirHcyWm8vC4Hja0GPDiFk4PRqXm6uWLeRQlIL67GV3tOeV0im8dyJ6eQX1aLD7YewnWp/TAgpJfuOGTjZgyLREq4P55fk42m1jbdcXqE5U5O4fk12fB0c8Hdkzg5GHXNxUXw8LRkFB1rwDs/FuiO0yMsd3J42wsqsWpfCW47fyBC/Dx1xyE7MS6+D85PCMHi73JwvN7+piVguZNDU0rhqa8zEerniVvPj9Udh+zMw9OSUNtkn9MSsNzJoa3aV4Idh47j/y5MgI+Hm+44ZGeSwvxx1cgovPNTAQ5V2Ne0BCx3cljNrQYsWpWFxL5+uGpktO44ZKfun5wIFxfguTX2dWITy50c1rtbCnCosh4LpiXB1YUnLFHPhAV44dbzBmDFniPYeeiY7jgmY7mTQ6qqb8Hi73JwXnwfjE8I0R2H7Nzt4weiTy9P/PXrTChlHyc2sdzJIb38XQ6qGlqwgNMMkBn08nTDvAsTsL3gGFbuLdEdxyQsd3I4B8rr8M5PB3HNqGikRPjrjkMO4upR0UgK88Oi1Zl2cWITy50czjMrM+Hh6oL7L0zQHYUciKuL4NGLU1BY2YC3fzyoO06XWO7kUH7Kq8A3GaW4c2IcQv28dMchBzMuvg8mJoZg8bpcVNQ26Y5zWix3chgGg8Jfv85ARIAXZo/jCUtkGY9cnIz6lja8uHa/7iinxXInh/HJ9iKkF1fjoalJvC4qWUxcqB9uGNMfH2w9hKySat1xTonlTg6hprEFz63Jwsj+Qbh0WITuOOTg7r0gHv7e7li4IsNmh0ay3MkhvLI+F+W1zXh8egqHPpLFBfp44L4LEvBDbgXWZpTqjtMpljvZvYPldXhz8wFcNTIKw6IDdcchJ3Hd2f0QF9oLT620zaGRLHeye08Zhz4+eBGvsETW4+7qgsemp6Cgoh5vbj6oO85JWO5k1zbuL8PaX4Y++nPoI1nX+IQQXJDcF4u/y0FpdaPuOL/Bcie71dxqwJ+/SkdMbx/cch6HPpIej09PQatB4ZmVmbqj/IZJ5S4iU0QkW0RyRWR+J49fLyJ7jB8/isgw80cl+q23fjiA/LI6PHHJIHi6cegj6dGvtw9uO28AvthVjLSDlbrj/KrLchcRVwBLAEwFkAJgpoiknLDYAQDjlVJDASwE8Lq5gxJ1VFrdiJfX5WBSUigmJoXqjkNO7s6JAxEe4IXHv0xHm8E2hkaasueeCiBXKZWvlGoGsAzAjI4LKKV+VEr9MtHxFgBR5o1J9FuLVmWhpU3h8UtO3M8gsj4fDzc8PC0ZGUeq8f7WQ7rjADCt3CMBFHb4ush436nMBrCqswdE5DYRSRORtLKyMtNTEnWw9UAlPt95GLePH4D+vX11xyECAEwfGo6xA3rj+dVZKLeBeWdMKffOzgjp9O8OEZmI9nJ/qLPHlVKvK6VGKaVGhYTwAgrUfS1tBjz2xT5EBnrjjgkDdcch+pWIYOFlg9DQ0oZFq/Rfks+Uci8C0PEClFEAik9cSESGAngDwAylVIV54hH91n9+OIjs0hr8+dJBvOA12Zy4UD/MHjcAn2wv0v7mqinlvg1AvIjEiogHgGsBLO+4gIj0A/AZgBuUUrY9VRrZrSNVDfj7t/txQXIoJqf01R2HqFN3T4pDRIAXHv1iH1rbDNpydFnuSqlWAHMBrAGQCeAjpVS6iMwRkTnGxR4H0BvAqyKyS0TSLJaYnNbCFRloMyg8cckg3VGITsnHww2PX5KCrJIa/EfjRT1M+rtWKbUSwMoT7lva4fYtAG4xbzSi/9mQfRQr95Zg3oUJiA720R2H6LQuGhSG3yWF4sW1+zF1SDgiA72tnoFnqJLNq29uxaNf7MPAEF/cev4A3XGIuiQiePLSQVAK+PPydC0ZWO5k8176NgdFxxrwzBVDeSYq2Y3oYB/cNzkeazNKsSa9xOrrZ7mTTUsvrsIbmw/g2tHRSI0N1h2HqFtuOjcWSWF+eOLLdNQ2tVp13Sx3slltBoUFn+1FkI87FkxN1h2HqNvcXV3wzBVDUFrTiOdXW3fsO8udbNZ/fjyIPUVVePySQQjwcdcdh6hHhvcLwo1jY/DOlgKrjn1nuZNNOlRRjxfWZGNiYgguGRquOw7RGXngokREBHjjoU/3oLHFOldtYrmTzVFKYf5ne+DqInjq8iG8JirZPV9PNzx9xRDkldVhyfpcq6yT5U42Z9m2QvyYV4EF05IQoWF8MJEljE8IwRUjIvHahjxkFFdbfH0sd7IpJVWNePrrTIwZEIyZo/vpjkNkVo9dnIJAH3d8uqPI4uvizEtkM345HNNiMGDRFUPh4sLDMeRYgnw98MVd51rljFXuuZPN+HBbITZkl2HB1GTE9OE87eSYooJ8rPI+EsudbEJhZT0WrsjA2AG9ccOY/rrjENk9ljtpZzAoPPDJbogInruKh2OIzIHlTtq9/dNBbMmvxGPTkznjI5GZsNxJq+ySGjyzKguTkkLx+1HRXX8DEZmE5U7aNLW24Z5lO+Hv5YZnrxrKk5WIzIhDIUmb51dnI6ukBm/NGo0+vTx1xyFyKNxzJy0255Tjjc0HcMOY/piYFKo7DpHDYbmT1ZXVNOG+j3YhLrQXHp7GqXyJLIGHZciqDAaF+z/aheqGFrw7OxXeHryyEpElcM+drOq17/OwKaccT1wyCElh/rrjEDksljtZTdrBSry4dj8uHhqOmakc9khkSSx3sory2ibMfX8nIgO98cwVnKOdyNJ4zJ0srrXNgLs/2Ilj9c349I5z4O/FS+YRWRrLnSzuxbX78WNeBZ67aigGRwbojkPkFHhYhixqbUYpXt2Qh5mp0ZxegMiKWO5kMblHa3Dfh7swJDIAT1wySHccIqfCcieLqKpvwS1vp8HL3QX/vGEkvNw5np3ImnjMncyutc2AuR/swOHjDfjg1jG8yDWRBix3MrtFq7KwKaccz1wxBKNignXHIXJKPCxDZvXfLQV4Y/MB3Di2P2am9tMdh8hpsdzJbDZkH8UTy9MxMTEEj01P0R2HyKmx3MksskqqMff9nUjs64fF142AmytfWkQ68TeQzljx8Qbc9NY29PJ0w5uzRqOXJ9/KIdLNpHIXkSkiki0iuSIyv5PHk0TkJxFpEpF55o9JtupYXTP++OZW1Da24s1ZoxEW4KU7EhHBhNEyIuIKYAmAyQCKAGwTkeVKqYwOi1UCuBvAZRZJSTapvrkVN7+9DYcq6/HOzalIieAUvkS2wpQ991QAuUqpfKVUM4BlAGZ0XEApdVQptQ1AiwUykg1qbjXgzvd2YHfhcbx87VkYM6C37khE1IEp5R4JoLDD10XG+7pNRG4TkTQRSSsrK+vJU5ANaG0z4J5lO7Ehuwx/vWwIpgwO1x2JiE5gSrl3NvG26snKlFKvK6VGKaVGhYSE9OQpSLM2g8K8j3dj1b4SPHpxMq47m2PZiWyRKeVeBKDjdH5RAIotE4dsmcGg8Mjne/HFrmI8cFEibjlvgO5IRHQKppT7NgDxIhIrIh4ArgWw3LKxyNYYDAoLPtuLZdsKMXdiHO6aGKc7EhGdRpejZZRSrSIyF8AaAK4A3lRKpYvIHOPjS0UkDEAaAH8ABhG5F0CKUqragtnJStoMCg98shuf7TiMu38Xh/smJ+iORERdMOlsE6XUSgArT7hvaYfbJWg/XEMOpqXNgHkf78aXu4px/+QE3D0pXnckIjIBTyWkU2psacNd7+3AuqyjeHBKIu6cwEMxRPaC5U6dqmpowS1vb0NawTEsvGwwbhjTX3ckIuoGljudpLS6ETe+uRV5ZbVYPHM4pg+N0B2JiLqJ5U6/kXmkGjf/ZxuqG1rw7xtH4/wEno9AZI9Y7vSr7/eX4a73dqCXpxs+mjMWgyICdEcioh5iuROUUnjnpwL8ZUUGEvr64c1ZoxAewOueEtkzlruTa2xpw2Nf7MPH24twQXJf/OPaszgfO5ED4G+xEztS1YA7/rsDuwqP4+5J8bh3UjxcXDqbSoiI7A3L3UltyD6K+z7cheZWA5b+YQRndiRyMCx3J9PaZsA/vs3BK+tzkRTmhyXXj8DAkF66YxGRmbHcncihinrc99EubC84hmtGRePPlw6Ct4er7lhEZAEsdyeglMIn24vw5+XpcHERvHTtWZhxVo+ut0JEdoLl7uBKqxvx6Bf7sDajFGfHBuNvvx+GqCAf3bGIyMJY7g7ql731hSsy0NRqwMPTkjB73AC4cjQMkVNguTugvLJaPP7lPvyQW4HUmGAsunIIBvBNUyKnwnJ3IA3NbXhlfQ5e35gPL3dXLLxsMK5P7cex60ROiOXuAAwGhS92HcZzq7NRUt2IK0ZEYsHUZIT4eeqORkSasNzt3M/5FXhqZSb2FFVhSGQAXp45HKmxwbpjEZFmLHc7tbeoCs9/k42N+8vQ198Tf7t6GC4fHslDMEQEgOVud/YdrsLi73KwJr0UgT7uWDA1CX8cG8OTkYjoN1judkAphbSCY1iyPhcbssvg5+WGeybFY/Z5sfD3ctcdj4hsEMvdhrW2GbBqXwne2HwAuwuPI9jXAw9clIgbxvZnqRPRabHcbdDRmkZ8uLUQH2w9hOKqRsT09sHCGYNw5cgo+Hjwn4yIusamsBFtBoWNOWX4OK0Q36SXotWgMC6uD56cMRiTkkL5RikRdQvLXSOlFLJKavDlrmJ8vrMIpdVNCPJxx6xzYnD9mP6I7eOrOyIR2SmWuwa5R2uwam8Jlu8uRs7RWri6CCYkhODJS6Pwu6S+8HBz0R2RiOwcy90K2gwKuwqP47usUqzeV4K8sjoAwOiYICycMQjThoSjdy+eTUpE5sNyt5DS6kb8kFuOTTl7dgUNAAAHV0lEQVTl+H5/GSrrmuHqIhgzIBizzonB5JQwhAV46Y5JRA6K5W4mZTVN2HqgEj8fqMCW/ArsL60FAAT7emBCQggmJoXi/PgQBPhwCCMRWR7LvQeaWw3YX1qDnYXHsbPgGHYcOoaDFfUAAB8PV4zsH4QrR0RhXHwfJIf5c6QLEVkdy70LdU2tyCqpQeaRamQcqUb64SpkHqlBc5sBANCnlydG9AvEzNR+SI0NxuDIALi78g1RItKL5Y72IYmVdc04UF6HvLJa5JfVIedoLbJLanD4eMOvy/l7uSElwh83nRuDIVEBGBoZiOhgb4hwz5yIbIvTlHttUyuKjzeg+HgDio41oPBYPYoqG1BQWYeC8nrUNLX+uqyHqwsGhPhiZP8gzEyNRkJfP6RE+CMykEVORPbBpHIXkSkAXgLgCuANpdSiEx4X4+PTANQDmKWU2mHmrCdpMyhUNbSgorYJ5bXNKK9tQllNE8qMn0urG1FS1YiS6kbUNLb+5ns9XF0QGeSNfsE+GNkvCP17+yK2jy8GhvRCZJA3rzVKRHaty3IXEVcASwBMBlAEYJuILFdKZXRYbCqAeOPH2QBeM342u/XZR7FwRQaO1TXjeEMLlDp5GTcXQYifJ0L9vTAgxBfnDOyN8EBvRAR6IyLAC5FB3ujr58U3OonIYZmy554KIFcplQ8AIrIMwAwAHct9BoB3lFIKwBYRCRSRcKXUEXMHDvR2R3K4P4J9PBDk444gXw/07uWJPsbPIX6eCPR2Z3ETkVMzpdwjARR2+LoIJ++Vd7ZMJACzl/vwfkFYcl2QuZ+WiMihmDJmr7Nd4BMPhpiyDETkNhFJE5G0srIyU/IREVEPmFLuRQCiO3wdBaC4B8tAKfW6UmqUUmpUSEhId7MSEZGJTCn3bQDiRSRWRDwAXAtg+QnLLAfwR2k3BkCVJY63ExGRabo85q6UahWRuQDWoH0o5JtKqXQRmWN8fCmAlWgfBpmL9qGQN1kuMhERdcWkce5KqZVoL/CO9y3tcFsBuMu80YiIqKc4CQoRkQNiuRMROSCWOxGRAxLV2fn71lixSBmAgh5+ex8A5WaMYy62mguw3WzM1T3M1T2OmKu/UqrLseTayv1MiEiaUmqU7hwnstVcgO1mY67uYa7uceZcPCxDROSAWO5ERA7IXsv9dd0BTsFWcwG2m425uoe5usdpc9nlMXciIjo9e91zJyKi07CLcheR50UkS0T2iMjnIhJ4iuWmiEi2iOSKyHwr5LpaRNJFxCAip3znW0QOisheEdklImk2lMva2ytYRNaKSI7xc6cT81tre3X18xsnwnvZ+PgeERlhqSzdzDVBRKqM22eXiDxupVxvishREdl3isd1ba+ucunaXtEisl5EMo2/j/d0sozltplSyuY/AFwIwM14+1kAz3ayjCuAPAADAHgA2A0gxcK5kgEkAtgAYNRpljsIoI8Vt1eXuTRtr+cAzDfent/Zv6O1tpcpPz/aJ8NbhfbrFYwB8LMV/u1MyTUBwAprvZ46rPd8ACMA7DvF41bfXibm0rW9wgGMMN72A7Dfmq8xu9hzV0p9o5T65QrXW9A+X/yJfr0coFKqGcAvlwO0ZK5MpVS2JdfREybmsvr2Mj7/28bbbwO4zMLrOx1Tfv5fLx+plNoCIFBEwm0glxZKqY0AKk+ziI7tZUouLZRSR5RSO4y3awBkov0KdR1ZbJvZRbmf4Ga0/093olNd6s8WKADfiMh2EblNdxgjHdurrzLO82/8HHqK5ayxvUz5+XVsI1PXOVZEdovIKhEZZOFMprLl30Gt20tEYgAMB/DzCQ9ZbJuZNOWvNYjItwDCOnnoEaXUl8ZlHgHQCuC9zp6ik/vOeCiQKblMcK5SqlhEQgGsFZEs496GzlxW317deBqzb69OmO3ykWZmyjp3oP0U9FoRmQbgCwDxFs5lCh3byxRat5eI9ALwKYB7lVLVJz7cybeYZZvZTLkrpS443eMiciOA6QAmKePBqhOYdKk/c+cy8TmKjZ+PisjnaP/T+4zKygy5rL69RKRURMKVUkeMf3oePcVzmH17dcJsl4+0dq6OBaGUWikir4pIH6WU7jlUdGyvLuncXiLijvZif08p9Vkni1hsm9nFYRkRmQLgIQCXKqXqT7GYKZcDtDoR8RURv19uo/3N4U7f1bcyHdtrOYAbjbdvBHDSXxhW3F62evnILnOJSJiIiPF2Ktp/jyssnMsUNnm5TV3by7jOfwPIVEq9eIrFLLfNrP0Ock8+0H75vkIAu4wfS433RwBY2WG5aWh/RzoP7YcnLJ3rcrT/z9sEoBTAmhNzoX3Uw27jR7qt5NK0vXoDWAcgx/g5WOf26uznBzAHwBzjbQGwxPj4XpxmRJSVc801bpvdaB9gcI6Vcn0A4AiAFuPra7aNbK+ucunaXuPQfohlT4fummatbcYzVImIHJBdHJYhIqLuYbkTETkgljsRkQNiuRMROSCWOxGRA2K5ExE5IJY7EZEDYrkTETmg/we9FehBh3NMwQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "alpha_grid = np.linspace(-2, 2, 100)\n", "P_alpha = scipy.stats.t.pdf(alpha_grid, loc=alpha_hat, scale=np.sqrt(s2_alpha), df=len(x) - 2)\n", "\n", "plt.plot(alpha_grid, P_alpha)\n", "plt.axvline(1, color='red')\n", "plt.axvline(alpha_hat)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Example: Construct a 95% confidence interval for the slope\n", "---\n", "\n", "Once we compute $S_\\beta$, we can use the formulas we've seen before:\n", "\n", "$$P(\\beta = \\hat{\\beta} \\pm y) = 0.95$$\n", "\n", "$$T = \\frac{y}{S_\\beta}$$\n", "\n", "$$y = TS_\\beta$$\n", "\n", "The confidence interval will then be:\n", "\n", "$$\\beta = \\hat{\\beta} \\pm TS_\\beta$$ with 95% confidence" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0202139147447661 2.10092204024096\n", "beta = 2.511951711962952 +/- 0.29869995143838984 with 95% confidence\n" ] } ], "source": [ "s2_beta = s2_epsilon / np.sum((x - np.mean(x))**2)\n", "T = scipy.stats.t.ppf(0.975, len(x) - 2)\n", "\n", "print(s2_beta, T)\n", "\n", "print('beta = ', beta_hat, '+/-', T * np.sqrt(s2_beta), ' with 95% confidence')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Notice that just like in confidence intervals, once $N$ becomes large we can replace the $t$-distribution with a normal distribution.\n", "\n", "The next interesting consequence of the $t$-distribution is that we can construct a hypothesis test. For example, we could test if the intercept should be $0$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Example: A hypothesis test for the existence of an intercept\n", "---\n", "\n", "Our null hypothesis is that the intercept is $0$. Then, we can compute how big an interval would have to be constructed around $0$ to just capture what we calculated for $\\alpha$. The distribution for that interval is:\n", "\n", "$$P(\\alpha) = T(\\mu=0, \\sigma=S_\\alpha, df=N - D = 19)$$\n", "\n", "We take $D = 1$ here because we work under the assumption of the null hpyothesis (i.e., only slope needs to be fit). For the actual regression analysis, you should still use $D = 2$. Only when you do the hypothesis test itself do you use the $D = 1$\n", "\n", "To convert to a standard $t$-distribution (variance is 1), so we have:\n", "\n", "$$T = \\frac{\\hat{\\alpha}}{S_\\alpha}$$\n", "\n", "$$\\int_{-T}^T p(T)\\,dT = 1 - p$$" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "alpha = 0.7923280760210064 T = 0.93885680852821 p-value = 0.3595872204956516\n" ] } ], "source": [ "df = len(x) - 2\n", "s2_epsilon = np.sum((y - alpha_hat - beta_hat * x) ** 2) / df\n", "s2_alpha = s2_epsilon * (1. / df + np.mean(x) ** 2 / (np.sum((np.mean(x) - x) ** 2)))\n", "#ensure our T-value is positive, so our integral doesn't get flipped\n", "T = abs(alpha_hat / sqrt(s2_alpha))\n", "p = 1 - (scipy.stats.t.cdf(T, len(x) - 1) - scipy.stats.t.cdf(-T, len(x) - 1))\n", "print('alpha = ', alpha_hat, ' T = ', T, ' p-value = ', p)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "---\n", "\n", "*Depends on random date above!*\n", "\n", "When I ran this, the $p$-value is $\\approx 0.04$ which says the evidence is weak but we can reject the null hypothesis. There is an intercept, since the $p$-value is below our threshold of 5%." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Extending Least-Squares to Multiple Dimensions in Domain - OLS-ND\n", "====\n", "\n", "Our governing equation is now:\n", "\n", "\n", "$$y = {\\mathbf X\\beta} + \\epsilon$$\n", "\n", "where $\\mathbf X$ is an $N\\times D$ matrix, where $N$ is the number of data points and $D$ is the number of dimensions. $\\beta$ is then a $D$ length column vector. We want to find $\\hat{\\beta}$ and get a model that looks like:\n", "\n", "$$\\hat{y} = {\\mathbf X\\hat{\\beta}}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "This can be done with optimization just like last time, but we can more easily do it with matrix algebra. The equation for $\\hat{\\beta}$ is:\n", "\n", "\n", "$$\\hat{\\beta} = (\\mathbf{X}^T \\mathbf{X})^{-1}\\mathbf{X}^Ty$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Example\n", "---\n", "\n", "Let's try to fit the equation:\n", "\n", "$$ y = 3 + 3 x_1 + x_2 + \\epsilon$$\n", "\n", "with \n", "\n", "$$\\hat{y} = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2$$" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "array([ 8.59319365, 8.2434237 , 12.21993175, 0.08357724, 10.2652057 ,\n", " 2.56500889, 2.10507327, 0.20001633, -0.80413993, 5.18654211,\n", " 9.00227094, 3.56647433, -0.94086123, 5.75923108, 4.54752716])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#NOTE: THIS IS NOT PART OF REGRESSION!!!!\n", "#DO NOT COPY PASTE THIS CODE INTO HW/EXAM\n", "#generate data\n", "#I add some noise to the x coordinate to just spread the points out a little.\n", "x1 = np.linspace(0,1,15)+ scipy.stats.norm.rvs(size=15)\n", "x2 = np.linspace(0,1,15) + scipy.stats.norm.rvs(size=len(x1))\n", "y = 3 * x1 - 2 * x2 + 3 + scipy.stats.norm.rvs(size=len(x1))\n", "y" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 1.48872779, -0.36739723],\n", " [ 1. , 0.64261885, -1.06126457],\n", " [ 1. , 1.16473294, -1.89915635],\n", " [ 1. , -1.21189489, -0.06291139],\n", " [ 1. , 1.14115123, -1.92180539],\n", " [ 1. , -0.38402235, 0.01482152],\n", " [ 1. , -0.12382456, 0.55549077],\n", " [ 1. , 0.28917314, 1.90087508],\n", " [ 1. , -0.22601835, 1.51870358],\n", " [ 1. , 1.51075808, 0.74723702],\n", " [ 1. , 1.55939315, -0.0256333 ],\n", " [ 1. , -0.59703574, -0.80557387],\n", " [ 1. , -0.19451134, 1.31367909],\n", " [ 1. , 0.35150002, -0.63603463],\n", " [ 1. , 1.10545322, 1.71878621]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy.linalg as linalg\n", "\n", "x_mat = np.column_stack( (np.ones(len(x1)), x1, x2) )\n", "x_mat" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Now we have our $X$ matrix set-up. Now we need to evaluate the matrix equation for $\\hat{\\beta}$ above:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "#dot -> matrix multiplication\n", "#transpose -> take a transpose\n", "#linalg.inv -> compute a matrix inverse\n", "beta_hat = linalg.inv(x_mat.transpose() @ x_mat) @ x_mat.transpose() @ y" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Since it is tedius to type that whole equation out, you can instead use a shortcut:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n", "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n", " \"\"\"Entry point for launching an IPython kernel.\n" ] } ], "source": [ "beta_hat, *_ = linalg.lstsq(x_mat, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `*_` symbol means put the rest of the return value into the `_` variable, which recall is how we indicate that we're making a variable which we will not use." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's now see how well the regression did!" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "y_hat = x_mat @ beta_hat" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The first plot will be a $\\hat{y}$ vs $y$ plot. This is called a **parity** plot. If $y$ and $\\hat{y}$ are the same, you would see a $y=x$ line. How far the deviation is from that line is how bad the fit is." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEMCAYAAADeYiHoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHL5JREFUeJzt3Xt0VfWd9/H314AaEE1FUBLEoMMTbcE2mHGstDNOvQQVCwtdfbTUOrWV9axxrO24ItLW8qzWFms6bbWtjtQqduDxUodBBy94rXcdoqGiYhRvQEAIoxmqRAnJ9/ljn0B2zg65nXP2Pud8Xmu5YP+yyf4uxfPJd//277fN3REREelpn7gLEBGRZFJAiIhIJAWEiIhEUkCIiEgkBYSIiERSQIiISCQFhIiIRFJAiIhIJAWEiIhEGhZ3AUNxyCGHeGVlZdxliIjklRdeeGGbu4/p67y8DojKykoaGhriLkNEJK+Y2bv9OU+3mEREJJICQkREIikgREQkkgJCREQiKSBERCSSAkJERCIpIEREJJICQkQkn7S1wYIFsGFD1i8VS0CY2c1mttXMXu42Vm9mr5nZS2b2H2ZWFkdtIiKJdcstMGIE/OhHsHJl1i8XVwexGJjeY+whYLK7Hwu8DszPdVEiIonU2gpmcOGFwfGcOfCtb2X9srEEhLs/AbzfY+xBd9+VOnwOGJ/zwkREkuaaa+BTn9pz/OabsGRJTi6d1L2YLgTuiLsIEZHYbN4M5eV7juvqgrDIocRNUpvZ94FdwNJevj7XzBrMrKGlpSW3xYmI5MI//3M4HN57L+fhAAkLCDO7AJgBzHF3jzrH3Re5e42714wZ0+dutSIi+WPdumCu4Ze/DI5//nNwh0MPjaWcxNxiMrPpwDzg79x9R9z1iIjk1Hnnwe237zlubYWDDoqvHuJ7zPU24Fmgysw2mtk3gd8Ao4CHzGy1mf1rHLWJiORUY2PQNXSFw+LFQdcQczhATB2Eu58XMfz7nBciIhKXzk446SR48sngePRo2LgR9t8/1rK6S9QchIhIUXjsMSgp2RMOK1bAtm2JCgdI0ByEiEjBa2+HY44J1jIATJkS3GIqKYm3rl4oIESkaC1vbKZ+ZRObWtsoLyulrraKWdUV2bnYsmVw9tl7jp96CqZNG/C3yWXNCggRKUrLG5uZv2wNbe0dADS3tjF/2RqAzH7g7tgBhxwSbLIHUFsL998fTEwPUM5qTtEchIgUpfqVTbs/aLu0tXdQv7IpcxdZtAhGjtwTDmvWwAMPDCocIEc1d6MOQkSK0qbWtgGND8j77wdPJXX5xjfg5puH/G2zWnMEdRAiUpTKy0oHNN5vP/lJOBzefjsj4QBZrLkXCggRKUp1tVWUDg8/PVQ6vIS62qrBfcPm5uDW0Q9+EBx/73vBgrfKyqEV2k3Ga+6DbjGJSFHqmtTNyBNBl1wCv/nNnuOtWyELe8VltOZ+sF72xMsLNTU13tDQEHcZIlKsmprg6KP3HP/qV3DppfHV009m9oK71/R1njoIEZGBcodzzgnWNnTZvh1GjYqvpizQHISIyEA0NMA+++wJh6VLg8AosHAAdRAiIv3T2QknngjPPx8cjxsXPKG0337x1pVF6iBERPry8MPBfkld4XD//bBpU0GHA6iDEBHp3c6d8Fd/BRs2BMfHHReEREI318s0dRAiIlHuvDPoELrC4dlng/mHIgkHUAchIhL24YdQVgYdqT2PzjoL7r570Psn5TN1ECIiXa6/PngaqSscXn0V7rmnKMMB1EGIiMB//3ewJXeXuXPhxhvjqych1EGISHG7/PJwOKxfr3BIUUCISHF65png1lF9fXC8YEGw4O3ww+OtK0FiucVkZjcDM4Ct7j45NXYwcAdQCbwDfMXdP4ijPhEpcD3nFLZsgbFj46klweLqIBYD03uMXQE84u6TgEdSxyIimbNsWTgcPvOZoGtQOESKpYNw9yfMrLLH8EzgpNTvbwX+BMzLWVEiUrjcg/2TulPX0KckzUEc6u6bAVK/6r+ciAzdtdeGw+Hss9U19FPePeZqZnOBuQATJkyIuRoRSaz2dth33/DYhx/CyJHx1JOHktRBbDGzcQCpX7dGneTui9y9xt1rxmThjU0iUgC++91wOMybF3QNCocBSVIHcQ9wAXB16te74y1HRPLOX/4CBx4YHmtvh2FJ+qjLH7F0EGZ2G/AsUGVmG83smwTBcKqZvQGcmjoWEemfM88Mh8MNNwRdg8Jh0OJ6ium8Xr50ck4LEZH8t3kzlJeHxzo7i3b/pExK0hyEiMjAHHVUOBz+8z+DrkHhkBHqvUQk/7z6arDIrTv3eGopYOogRCS/mIXD4b/+S+GQJQoIEckPjz8evnU0YkQQDH/91/HVVOB0i0lEkq/nnMK6dcH8g2SVOggRSa7bbw+HQ01N0DUoHHJCHYSIJE/U5notLeEX+0jWqYMQkWSprw+Hw5w5QWAoHHJOHYSIZNXyxmbqVzaxqbWN8rJS6mqrmFVdkX5iW1sw8dzdjh1QWpqbQiWNOggRyZrljc3MX7aG5tY2HGhubWP+sjUsb2wOnzhhQjgcfvjDoGtQOMRKHYSIZE39yiba2jtCY23tHdSvbAq6iC1b4LDDwn9o1y4oKclhldIbdRAikjWbWtt6HzcLh8PFFwddg8IhMdRBiEjWlJeV0twjJCa1vMtDN18cPlGb6yWSOggRyZq62ipKh+/pCN752YxwOPz619pcL8HUQYhI1nQ9rfTcL37P1UsWhL+o/ZMSTwEhIlk1a+p4ZnUfuPdeOOOMuMqRAVBAiEh2XHNN8C7o7nrpGvq9VkJySgEhIpnXc07huefgb/4m8tSutRJdj8N2rZUAFBIx0yS1iGTO+eenh4N7r+EAe18rIfFSByEiQ9fRAcN6fJy8+26wQroPe10rIbFSByEiQzN3bno4uPcrHCBYKzGQccmdRAWEmX3XzF4xs5fN7DYz2z/umkSkFzt2BLeTfve7PWPbtw/48dWeayUASoeXUFdblYkqZQgSExBmVgF8G6hx98lACXBuvFWJSKQTT4SRI/ccX3RREAyjRg34W82qrmDh7ClUlJViQEVZKQtnT9EEdQIkbQ5iGFBqZu3ACGBTzPWISHdbt8Khh4bHOjrSX+4zQLOqKxQICZSYDsLdm4GfA+uBzcD/uPuD8VYlIruNHBkOh/p6cGf5nzcz7epHmXjFvUy7+tH0rbwlbyWmgzCzTwEzgYlAK/BHM/uauy/pcd5cYC7AhH5OgonIELz2GhxzTHgsNc+gNQyFLTEdBHAK8La7t7h7O7AMOLHnSe6+yN1r3L1mzJgxOS9SpKiYhcPhj38MTUJrDUNhS0wHQXBr6QQzGwG0AScDDfGWJFKkHn8cTjopPBbxdJLWMBS2xHQQ7v48cBfwIrCGoLZFsRYlUozMwuHwzDO9PrqqNQyFLTEBAeDuC9z9aHef7O7nu/sncdckUjSWLo3eJuPzn+/1j2gNQ2FL0i0mEYmDe/pjquvWwVFH9flHuyaitRNrYVJAiBSzq66CK6/cczxmTLDWYQC0hqFwKSBEitGuXTB8eHhs2zYYPTqeeiSREjUHISI58PWvh8PhlFOC20wKB+lBHYRIsfjww/S9ktraYH/tiSnR1EGIFIOpU8PhcMklQdegcJC9UAchUsg2b4by8vBYBjbXk+KgvyUihcosHA7XXRf9SKtIL9RBiBSal1+GKVPCYwN8iY8IqIMQKSxm4XC4+26FgwyaOgiRQvDww3DqqeExBYMMkQJCJN/13D9p1SqoqYmnFikousUkkq8WL47eXE/hIBmiDkIk30Q9ifTOO3DEEbGUI4VLHYRIPyxvbE7Ge5evvDIcDkccEQSGwkGyQB2ESB8S8d7l9nbYd9/w2AcfQFlZbq4vRUkdhEgfYn/v8le+Eg6Hs84KugaFg2SZOgiRPsT23uXt2+Ggg8Jjn3yS3kmIZIk6CJE+xPLe5aOPDofD5ZcHXYPCQXJIHYRIH+pqq0JzEJDF9y5v3AiHHx4e6+xMf5xVJAfUQYj0YVZ1BQtnT6GirBQDKspKWTh7SuYnqM3C4XDjjUHXoHCQmCSugzCzMuAmYDLgwIXu/my8VUmxy+p7l1evhurq8Ji2yZAESFxAANcCD7j7OWa2LzAi7oJEBmN5YzP1K5vY1NpGeVkpdbVV6SHTszu4/36YPj13RYrsRaICwswOBP4W+AcAd98J7IyzJpHB6HPtxP33wxlnhP+QugZJmKTNQRwJtAC3mFmjmd1kZiPjLkpkoPa6dsIsHA6rVyscJJGSFhDDgKnADe5eDXwEXNH9BDOba2YNZtbQ0tISR40ifYpaIzGn8T6enn9yeNAdPvvZHFUlMjCJusUEbAQ2uvvzqeO76BEQ7r4IWARQU1OjH7skkcrLSmnuCgl33rnmrPAJGzbA+PG5L0xkABLVQbj7e8AGM+t6wPxk4NUYSxIZlLraKkqHl3DH0nmhcNh+5P8KugaFg+SBpHUQAJcAS1NPML0FfCPmekQGbNYxo5l11emhsRVPvsaML2RhcZ1IlgwoIMzsYeAyd/9zlurB3VcDeuOJ5K9hw6Cj2wT16NGwbRsz4qtIZFD2eovJzD5tZku6DV0O/NLMbjGzcdktTSTPbNkSPKHUPRw++QS2bYuvJpEh6GsO4hHgB10H7v6iu38JWAE8YGYLzCyLO5aJ5AkzOOywPcdnnKHN9STv9RUQpwE/6T5gZgY0ATcQzBe8YWbnZ6c8kYRbsyZ9NXRnJ9x7bzz1iGTQXgPC3de4+5yuYzN7CmgGfglUEKx4Pgk43swWZa9MkQQyg2OP3XP8ve9pcz0pKAN9iun/AK+4py37vMTM1maoJpFku+8+OPPM8JhWQksBGtA6CHd/OSIcupzZy7hI4TALh8Mf/qBwkIKVsYVy7v5Wpr6XSOJce236rSN3OF/Tb1K4krhQTiRZegbDk0/CF74QTy0iOZSorTZEEuWii6K7BoWDFAl1ECI9dXZCSUl47M034cgj46lHJCbqIES6mzo1PRzcFQ5SlNRBiAB89BEccEB4rLUVDjoonnpEEkAdhIhZOByOOCLoGhQOUuTUQUjxam5Ofy9De3uwG6uIqIOQImUWDodzzgm6BoWDyG76v0GKy4svwnHHhce0ElokkjoIKR5m4XD48Y8VDiJ7oQ5CCt9DD8Fpp4XHFAwifVIHIYXNLBwO992ncBDpJwWEFKZFi6K3yTj99HjqEclDusUkhadnMPz5z+EX+4hIvySugzCzEjNrNLMVcdcieeayy6K7BoWDyKAksYO4FFgLHBh3IZInOjrS1y9s2gTjxsVTj0iBSFQHYWbjCd5Md1PctUieqK0Nh8PYsUHXoHAQGbKkdRC/Ai4HRsVdiCTchx/CqFHpYyNHxlOPSAFKTAdhZjOAre7+Qh/nzTWzBjNraGlpyVF1kigHHxwOh9NPD7oGhYNIRiWpg5gGfNnMzgD2Bw40syXu/rXuJ7n7ImARQE1NjR5oLyZRm+vt2pX+/gYRyYjEdBDuPt/dx7t7JXAu8GjPcJAi1nNzvbq6oGtQOIhkTZI6CJF0q1dDdXV4TCuhRXIiMR1Ed+7+J3efEXcdEjOzcDj87ncKB5EcSmRASHF79ro/RC94+9a34ilIpEjpFpMkixmf73b41f99FY2TjmNhYzOzqitiK0ukGKmDkGS47rq0rqFy3gqeqfwcbe0d1K9siqkwkeKlDkLi5Q77hH9OOeWb17PukAmhsU2tbbmsSkRQByFxuvjitHCYtvCRtHAAKC8rzVVVIpKiDkJyb9cuGD48PLZlC4wdS11jM/OXraGtvWP3l0qHl1BXWzXkyy5vbKZ+ZRObWtsoLyulrrZK8xoie6EOQnLri18Mh0NlZXCbaexYAGZVV7Bw9hQqykoxoKKslIWzpwz5g3x5KniaW9twoLm1jfnL1rC8sXlI31ekkKmDkNzYvh0OOig8tmMHlKbfOppVXZHxn+zrVzaFuhJg9+S3ugiRaOogJPv22y8cDrNnB11DRDhkS2+T3Jr8FumdOgjJnnffDW4hddfRkTYxnQvlZaU0R4SBJr9FeqcOQrLDLBwOV14Z+UhrrtTVVlE6PLyxX6Ymv0UKlToIyaxVq+D448NjCdg/qWueQU8xifSfAkIyp+f+Sf/2b/C15OzYno3Jb5FCpoCQoVu2DM4+OzyWgK5BRIZGASFD07NreOKJYK2DiOQ9TVLL4NTXR2/JrXAQKRjqIGRgop5Eev11mDQpnnpEJGvUQUj/XXhh5OZ6yz8cEVNBIpJN6iCkbzt3Bquhu/nct/8fraUHQmpPI0BPCIkUGHUQsndTp4bC4a1DK6mctyIIhxS90EekMKmDkGgffAAHHxwe+/hjTl7wcOTp2tNIpPAkqoMws8PN7DEzW2tmr5jZpXHXVJTMwuEwZ04wOb3ffr3uXaQ9jUQKT6ICAtgFXObuxwAnABeb2adjrql4vPVW+qOrnZ2wZMnuQ+1pJFI8EhUQ7r7Z3V9M/f4vwFpAM5+5YAZHHbX78MbTLmT5ixvTAiNbL/QRkeRJ7ByEmVUC1cDz8VZS4J55BqZNCw1VzlsBQGkvTydpTyOR4pCoDqKLmR0A/DvwHXff3uNrc82swcwaWlpa4imwUJiFwuEfZ16xOxxATyeJFLvEBYSZDScIh6Xuvqzn1919kbvXuHvNmDFjcl9gIbjttrRbRxPnreC+o7+QdqqeThIpXokKCDMz4PfAWnf/Rdz1FCQz+OpX9xw/+yy46+kkEUmTqIAApgHnA18ys9Wpf86Iu6iC8OMfR2+ud8IJgJ5OEpF0iZqkdvenAOvzROm/qM313noLJk4MDemNayLSU6ICQjLs3HPhjjvCY3t5kY+eThKR7hQQhejjj6G0x9zBBx9AWVk89YhIXkraHIQMVVVVOByOPz7oGhQOIjJA6iAKxbZt0POx3507YfjweOoRkbynDqIQmIXD4aKLgq5B4SAiQ6AOIp81NcHRR4fHOjvTH2cVERkEdRD5yiwcDv/yL0HXoHAQkQxRB5Fv/vQn+Pu/D4/t5dFVEZHBUgeRT8zC4bB8ucJBRLJGAZEPFi+O3iZj5sxYyhGR4qBbTEnXMxgaGuC44+KpRUSKijqIpPr+96O7BoWDiOSIOoik6eyEkvCuqqxfD4cfHk89IlK01EEkycyZ4XA44ICga1A4iEgM1EEkwY4dMHJkeGz7dhg1Kp56RERQBxG/8ePD4XDSSUHXoHAQkZipg4jLe+/BuHHhsfZ2GKb/JCKSDOog4mAWDodLLgm6BoWDiCSIPpFy6ZVXYPLk8JhWQotIQqmDyKLljc1Mu/pRJl5xb9A1dA+H669XOIhIoqmDyJLljc3MX7aGmtdX8fSdPwx/UcEgInkgUQFhZtOBa4ES4CZ3vzrmkgatfmUTa686PTT2D+f8X9447os8HVNNIiIDkZiAMLMS4LfAqcBGYJWZ3ePur2byOssbm6lf2cSm1jbKy0qpq61iVnVFJi8BN93E0/MvCg1VzlsBgLW2ZfZaIiJZkpiAAI4H1rn7WwBmdjswE8hYQHTd9mlr7wCgubWN+cvWAGQuJHrsnzT9G7/mtbETdx+Xl5Vm5joiIlmWpEnqCmBDt+ONqbGMqV/ZtDscurS1d1C/smno3/zWW0Ph8EnJMCrnrQiFQ+nwEupqq4Z+LRGRHEhSBxH1rsy02VwzmwvMBZgwYcKALrCpl9s7vY33S8Tmesdeejvb9z8gNFaRrdtZIiJZkqQOYiPQfVe68cCmnie5+yJ3r3H3mjFjxgzoAr3d3hn0bZ8f/SgUDrcfexqV81akhYMBT1/xJYWDiOSVJAXEKmCSmU00s32Bc4F7MnmButoqSoeHf9of1G2fjz8ObictWBAa+/V58yJP17yDiOSjxASEu+8C/glYCawF7nT3VzJ5jVnVFSycPYWKslKM4LbPwtlTBvaT/QUXQGm3D/yf/jRY17DffpkLIBGRBDDP40VbNTU13tDQkJuLvf8+jB4dHuvogH3CGZuTx2hFRIbAzF5w95q+zkvSJHUs+vWBPm0aPPPMnuMlS2DOnMjvN6u6QoEgIgWhqAOiz3URb78NRx4Z/kN53HGJiAxEYuYg4rDXdREHHhgOh0ceUTiISFEp6g4iav3D5PfWseJn3wkPKhhEpAgVdUCUl5XS3C0k3vnZjPAJL70EU6bkuCoRkWQo6ltMXY+lnvTmqlA47Dh0XNA1KBxEpIgVdQcx63PlnPDtr3PYU4/uHnvggVVMr+3z6S8RkYJXvAHxwgtQU8NhXcenngoPPsj0OGsSEUmQ4gyId96BmlSXMGYMrF8P++8fa0kiIklTnHMQo0bBySfDihWwdavCQUQkQnF2EKNHw8MPx12FiEiiFWcHISIifVJAiIhIJAWEiIhEUkCIiEgkBYSIiERSQIiISCQFhIiIRFJAiIhIpLx+J7WZtQDvxl1HD4cA2+IuYhDytW7I39pVd27la92Q+dqPcPcxfZ2U1wGRRGbW0J+XgSdNvtYN+Vu76s6tfK0b4qtdt5hERCSSAkJERCIpIDJvUdwFDFK+1g35W7vqzq18rRtiql1zECIiEkkdhIiIRFJAZJCZTTezJjNbZ2ZXxF1Pf5jZ4Wb2mJmtNbNXzOzSuGsaCDMrMbNGM1sRdy39ZWZlZnaXmb2W+vf++bhr6i8z+27q78nLZnabmSXybVtmdrOZbTWzl7uNHWxmD5nZG6lfPxVnjVF6qbs+9XflJTP7DzMry1U9CogMMbMS4LfA6cCngfPM7NPxVtUvu4DL3P0Y4ATg4jypu8ulwNq4ixiga4EH3P1o4LPkSf1mVgF8G6hx98lACXBuvFX1ajGkvWL+CuARd58EPJI6TprFpNf9EDDZ3Y8FXgfm56oYBUTmHA+sc/e33H0ncDswM+aa+uTum939xdTv/0LwYVURb1X9Y2bjgTOBm+Kupb/M7EDgb4HfA7j7TndvjbeqARkGlJrZMGAEsCnmeiK5+xPA+z2GZwK3pn5/KzArp0X1Q1Td7v6gu+9KHT4HjM9VPQqIzKkANnQ73kiefNB2MbNKoBp4Pt5K+u1XwOVAZ9yFDMCRQAtwS+rW2E1mNjLuovrD3ZuBnwPrgc3A/7j7g/FWNSCHuvtmCH4wAsbGXM9gXAjcn6uLKSAyxyLG8uYRMTM7APh34Dvuvj3uevpiZjOAre7+Qty1DNAwYCpwg7tXAx+RzFsdaVL37GcCE4FyYKSZfS3eqoqHmX2f4Jbw0lxdUwGRORuBw7sdjyeh7XdPZjacIByWuvuyuOvpp2nAl83sHYLbeV8ysyXxltQvG4GN7t7Vpd1FEBj54BTgbXdvcfd2YBlwYsw1DcQWMxsHkPp1a8z19JuZXQDMAOZ4DtcmKCAyZxUwycwmmtm+BJN398RcU5/MzAjuh69191/EXU9/uft8dx/v7pUE/64fdffE/zTr7u8BG8ysKjV0MvBqjCUNxHrgBDMbkfp7czJ5MsGecg9wQer3FwB3x1hLv5nZdGAe8GV335HLaysgMiQ1ifRPwEqC/2nudPdX4q2qX6YB5xP8BL469c8ZcRdV4C4BlprZS8DngJ/GXE+/pLqeu4AXgTUEnx+JXJ1sZrcBzwJVZrbRzL4JXA2camZvAKemjhOll7p/A4wCHkr9//mvOatHK6lFRCSKOggREYmkgBARkUgKCBERiaSAEBGRSAoIERGJpIAQEZFICggREYmkgBDJEDObYmZPdzueamaPxlmTyFBooZxIhpjZPgT7b1W4e4eZPUbwro0XYy5NZFCGxV2ASKFw904zewX4jJlNAtYrHCSfKSBEMus5gv2t/pH0N4OJ5BUFhEhmPUfw2sjfpl6wI5K3NAchkkGpW0uPA5Pc/aO46xEZCj3FJJJZlwLzFQ5SCBQQIhlgZkeZ2WtAqbvfGnc9IpmgW0wiIhJJHYSIiERSQIiISCQFhIiIRFJAiIhIJAWEiIhEUkCIiEgkBYSIiERSQIiISKT/D/j9siUN4BrsAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(y, y_hat, 'o')\n", "plt.plot(y, y, 'r')\n", "plt.xlabel('$y$')\n", "plt.ylabel('$\\hat{y}$')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Of course we can also look at the histogram of residuals" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAD95JREFUeJzt3X+oX3d9x/HnyzSiTKEbudosyW0cC2NWUMslthRGELe1aVk2qBBhVsrg0lJBQRiZg4r/1f0ho8Y1BC22zCmC2oU2meuYxfaPVJMsra1Rlkm3XhLWWjE1tCjR9/64x3H59nvzPffec3NzP30+4HDPj8895/3hNC9OP/f8SFUhSWrLG9a6AEnS8Ax3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoOuWKsDb9q0qbZv375Wh5ekden48eM/qaqpSe3WLNy3b9/OsWPH1urwkrQuJfnvPu0clpGkBhnuktQgw12SGmS4S1KDDHdJalDvcE+yIcl/JHl4zLYkuTfJ6SRPJ7l22DIlSUuxlCv3jwGnFtl2E7Cjm2aB+1ZYlyRpBXqFe5KtwM3AFxZpsgd4sOYdBa5MsnmgGiVJS9T3yv3vgb8Gfr3I9i3A8wuW57p1kqQ1MPEJ1SS3AC9U1fEkuxZrNmbda768nWSW+WEbpqenl1CmdGlt3/fImhz3uXtuXpPjqj19rtxvAP4syXPAV4H3J/nHkTZzwLYFy1uBM6M7qqqDVTVTVTNTUxNfjSBJWqaJ4V5Vf1NVW6tqO7AX+Peq+suRZoeA27q7Zq4DzlXV2eHLlST1sewXhyW5A6CqDgCHgd3AaeAV4PZBqpMkLcuSwr2qHgMe6+YPLFhfwF1DFiZJWj6fUJWkBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGTQz3JG9K8t0kTyV5Nsmnx7TZleRckpPddPfqlCtJ6qPPZ/Z+Aby/qs4n2Qg8keRIVR0dafd4Vd0yfImSpKWaGO7d91HPd4sbu6lWsyhJ0sr0GnNPsiHJSeAF4NGqenJMs+u7oZsjSa4ZtEpJ0pL0Cveq+lVVvQfYCuxM8q6RJieAq6vq3cDngIfG7SfJbJJjSY69+OKLK6lbknQRS7pbpqp+BjwG3Diy/uWqOt/NHwY2Jtk05vcPVtVMVc1MTU0tv2pJ0kX1uVtmKsmV3fybgQ8APxxpc1WSdPM7u/2+NHy5kqQ++twtsxl4IMkG5kP7a1X1cJI7AKrqAHArcGeSC8CrwN7uD7GSpDXQ526Zp4H3jll/YMH8fmD/sKVJkpbLJ1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQX2+ofqmJN9N8lSSZ5N8ekybJLk3yekkTye5dnXKlST10ecbqr8A3l9V55NsBJ5IcqSqji5ocxOwo5veB9zX/ZQkrYGJV+4173y3uLGbRj9+vQd4sGt7FLgyyeZhS5Uk9dXnyp0kG4DjwO8Dn6+qJ0eabAGeX7A81607O7KfWWAWYHp6epkl61Lbvu+RNTnuc/fcvCbHlVrQ6w+qVfWrqnoPsBXYmeRdI00y7tfG7OdgVc1U1czU1NTSq5Uk9bKku2Wq6mfAY8CNI5vmgG0LlrcCZ1ZUmSRp2frcLTOV5Mpu/s3AB4AfjjQ7BNzW3TVzHXCuqs4iSVoTfcbcNwMPdOPubwC+VlUPJ7kDoKoOAIeB3cBp4BXg9lWqV5LUw8Rwr6qngfeOWX9gwXwBdw1bmiRpuXxCVZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhrU5xuq25J8O8mpJM8m+diYNruSnEtyspvuXp1yJUl99PmG6gXgE1V1IslbgeNJHq2qH4y0e7yqbhm+REnSUk28cq+qs1V1opv/OXAK2LLahUmSlm9JY+5JtjP/sewnx2y+PslTSY4kuWaR359NcizJsRdffHHJxUqS+ukd7kneAnwd+HhVvTyy+QRwdVW9G/gc8NC4fVTVwaqaqaqZqamp5dYsSZqgV7gn2ch8sH+5qr4xur2qXq6q8938YWBjkk2DVipJ6q3P3TIBvgicqqrPLtLmqq4dSXZ2+31pyEIlSf31uVvmBuDDwPeTnOzWfRKYBqiqA8CtwJ1JLgCvAnurqlahXklSDxPDvaqeADKhzX5g/1BFSZJWxidUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUF9vqG6Lcm3k5xK8mySj41pkyT3Jjmd5Okk165OuZKkPvp8Q/UC8ImqOpHkrcDxJI9W1Q8WtLkJ2NFN7wPu635KktbAxCv3qjpbVSe6+Z8Dp4AtI832AA/WvKPAlUk2D16tJKmXPlfu/y/JduC9wJMjm7YAzy9YnuvWnR35/VlgFmB6enpplV4mtu97ZE2O+9w9N6/JcSWtT73/oJrkLcDXgY9X1cujm8f8Sr1mRdXBqpqpqpmpqamlVSpJ6q1XuCfZyHywf7mqvjGmyRywbcHyVuDMysuTJC1Hn7tlAnwROFVVn12k2SHgtu6umeuAc1V1dpG2kqRV1mfM/Qbgw8D3k5zs1n0SmAaoqgPAYWA3cBp4Bbh9+FIlSX1NDPeqeoLxY+oL2xRw11BFSZJWxidUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN6vOZvfuTvJDkmUW270pyLsnJbrp7+DIlSUvR5zN7XwL2Aw9epM3jVXXLIBVJklZs4pV7VX0H+OklqEWSNJChxtyvT/JUkiNJrhlon5KkZeozLDPJCeDqqjqfZDfwELBjXMMks8AswPT09ACHliSNs+Ir96p6uarOd/OHgY1JNi3S9mBVzVTVzNTU1EoPLUlaxIrDPclVSdLN7+z2+dJK9ytJWr6JwzJJvgLsAjYlmQM+BWwEqKoDwK3AnUkuAK8Ce6uqVq1iSdJEE8O9qj40Yft+5m+VlCRdJnxCVZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho0MdyT3J/khSTPLLI9Se5NcjrJ00muHb5MSdJS9Lly/xJw40W23wTs6KZZ4L6VlyVJWomJ4V5V3wF+epEme4AHa95R4Mokm4cqUJK0dEOMuW8Bnl+wPNetkyStkSsG2EfGrKuxDZNZ5odumJ6eXvYBt+97ZNm/q/Xj9XieX499fj167p6bV/0YQ1y5zwHbFixvBc6Ma1hVB6tqpqpmpqamBji0JGmcIcL9EHBbd9fMdcC5qjo7wH4lScs0cVgmyVeAXcCmJHPAp4CNAFV1ADgM7AZOA68At69WsZKkfiaGe1V9aML2Au4arCJJ0or5hKokNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1qFe4J7kxyY+SnE6yb8z2XUnOJTnZTXcPX6okqa8+31DdAHwe+GNgDvhekkNV9YORpo9X1S2rUKMkaYn6XLnvBE5X1Y+r6pfAV4E9q1uWJGkl+oT7FuD5Bctz3bpR1yd5KsmRJNcMUp0kaVkmDssAGbOuRpZPAFdX1fkku4GHgB2v2VEyC8wCTE9PL7FUSVJffa7c54BtC5a3AmcWNqiql6vqfDd/GNiYZNPojqrqYFXNVNXM1NTUCsqWJF1Mn3D/HrAjyTuSvBHYCxxa2CDJVUnSze/s9vvS0MVKkvqZOCxTVReSfBT4FrABuL+qnk1yR7f9AHArcGeSC8CrwN6qGh26kSRdIn3G3H8z1HJ4ZN2BBfP7gf3DliZJWi6fUJWkBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QG9Qr3JDcm+VGS00n2jdmeJPd2259Ocu3wpUqS+poY7kk2AJ8HbgLeCXwoyTtHmt0E7OimWeC+geuUJC1Bnyv3ncDpqvpxVf0S+CqwZ6TNHuDBmncUuDLJ5oFrlST11CfctwDPL1ie69YttY0k6RK5okebjFlXy2hDklnmh20Azif5UY/jrwebgJ+s5gHymdXc+0Wtet/WkH1bn9Z93y7y77lP367uc4w+4T4HbFuwvBU4s4w2VNVB4GCfwtaTJMeqamat61gN9m19sm/r05B96zMs8z1gR5J3JHkjsBc4NNLmEHBbd9fMdcC5qjo7RIGSpKWbeOVeVReSfBT4FrABuL+qnk1yR7f9AHAY2A2cBl4Bbl+9kiVJk/QZlqGqDjMf4AvXHVgwX8Bdw5a2rjQ31LSAfVuf7Nv6NFjfMp/LkqSW+PoBSWqQ4b4MST6Y5Nkkv06y6F+2J7224XKU5HeSPJrkP7ufv71Iu+eSfD/JySTHLnWdS9Hy6zN69G1XknPdeTqZ5O61qHOpktyf5IUkzyyyfT2fs0l9G+acVZXTEifgD4E/AB4DZhZpswH4L+D3gDcCTwHvXOvae/Tt74B93fw+4DOLtHsO2LTW9fboz8TzwPzNAEeYf17jOuDJta57wL7tAh5e61qX0bc/Aq4Fnllk+7o8Zz37Nsg588p9GarqVFVNegCrz2sbLkd7gAe6+QeAP1/DWobQ8usz1ut/YxNV1XeAn16kyXo9Z336NgjDffWs11cyvL26ZxS6n29bpF0B/5rkePfk8eWq5ddn9K37+iRPJTmS5JpLU9qqW6/nrK8Vn7Net0K+HiX5N+CqMZv+tqr+uc8uxqy7LG5NuljflrCbG6rqTJK3AY8m+WF3RXK5Gez1GZehPnWfAK6uqvNJdgMPMf/21vVuvZ6zPgY5Z4b7IqrqAyvcRa9XMqyFi/Utyf8m2VxVZ7v/zX1hkX2c6X6+kOSbzA8RXI7hPtjrMy5DE+uuqpcXzB9O8g9JNlXVun43C+v3nE001DlzWGb19Hltw+XoEPCRbv4jwGv+LyXJbyV562/mgT8Bxv7l/zLQ8uszJvYtyVVJ0s3vZP7f/EuXvNLhrddzNtFQ58wr92VI8hfA54Ap4JEkJ6vqT5P8LvCFqtpdi7y2YQ3L7use4GtJ/gr4H+CDAAv7Brwd+Gb3398VwD9V1b+sUb0Xtdh5aOH1GT37ditwZ5ILwKvA3upuybicJfkK83eNbEoyB3wK2Ajr+5xBr74Ncs58QlWSGuSwjCQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalB/wdFF523SuSl2QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.hist(y - y_hat)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Error Analysis For Multidimensional Least Squares - OLS-ND\n", "====" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Just like before, if we can find the standard error for the noise/residual, $\\epsilon$, then we can find everything else. The equation for that is:\n", "\n", "$$S^2_{\\epsilon} =\\frac{\\sigma^2_{\\epsilon}}{N - D} = \\frac{1}{N - D}\\sum_i \\left(\\hat{y}_i - y_i\\right)^2$$\n", "\n", "where $D$ is the dimension of $\\beta$. Knowing that, the standard error for $\\hat{\\beta}$ is \n", "\n", "$$S^2_{\\beta} = S^2_{\\epsilon} \\left(\\mathbf{X}^T\\mathbf{X}\\right)^{-1}$$\n", "\n", "where the diagonal elements are the standard errors. The off-diagonal elements are unrelated.\n", "\n", "And again, $P(\\beta - \\hat{\\beta}) \\propto T(0, \\sigma = S_\\beta, df = N - D)$\n", "\n", "The $N - D$ term is called the degrees of freedom. $D$ is the number of fit coefficients." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Transforming Problems into Multiple Domain Dimensions\n", "====\n", "\n", "One of the most common uses for this form of least-squares is if your problem is non-linear and you want to linearize it. For example, let's say I am in 1 dimension and the model I'd like to have is:\n", "\n", "$$y = \\beta_2 x^2 + \\beta_1 x + \\beta_0 + \\epsilon$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "I can package that into a vector, like this: $[x^2, x, 1]$ and create an $\\mathbf{X}$ matrix by creating a vector for each $x_i$. If my $x$ values are 1,2,3, that would like:\n", "\n", "$$\\left[ \\begin{array}{lcr}\n", "1^2 & 1 & 1\\\\\n", "2^2 & 2 & 1\\\\\n", "3^2 & 3 & 1\\\\\n", "\\end{array}\\right]$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Example: Linearizing a Polynomial for OLS-ND\n", "----" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#NOTE THIS IS NOT PART OF REGRESSION!\n", "#make some data to regress\n", "x = np.linspace(-3, 3, 25)\n", "y = 2 * x ** 2 - 3 * x + 4 + scipy.stats.norm.rvs(size=len(x), loc=0, scale=1.5)\n", "#END" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VGX6xvHvk0aCIEVCl6KyCFICIqJRRECwoALqKmLfBV37qqzYFjusuOqiriwqomtjf4JYQJAiqAgKCAoIEVFQgkozSEkk5f39cYZIICHJZCZnZnJ/rmuuZM6cmfMMCXfOvOct5pxDRESiX5zfBYiISGgo0EVEYoQCXUQkRijQRURihAJdRCRGKNBFRGKEAl1EJEYo0EVEYoQCXUQkRiRU5sHq1avnWrRoUZmHFBGJekuWLNninEstbb9KDfQWLVqwePHiyjykiEjUM7P1ZdlPTS4iIjFCgS4iEiMU6CIiMaJS29BFxD+5ubls2LCBnJwcv0uREiQnJ9O0aVMSExODer4CXaSK2LBhAzVr1qRFixaYmd/lyH6cc2zdupUNGzbQsmXLoF4j4gN9ytJMRs/IYGNWNo1rpzCsb2v6d2rid1kiUScnJ0dhHsHMjMMOO4zNmzcH/RoRHehTlmZyx+TlZOfmA5CZlc0dk5cDKNRFgqAwj2wV/flE9EXR0TMyCsN8r+zcfEbPyPCpIhGRyBXRgb4xK7tc20Uketx77708+uijJT4+ZcoUvvrqq0qsKPpFdKA3rp1Sru0iEjpTlmaSPmoOLYdPJX3UHKYszazc4yvQyy2iA31Y39akJMYX2ZaSGM+wvq19qkikath7/SozKxvH79evKhrqDz30EK1bt6Z3795kZHhNp88++yzHHXccHTt25LzzzmP37t188sknvP322wwbNoy0tDTWrl1b7H5SVEQHev9OTRg5sD1NaqdgQJPaKYwc2F4XREXCLBzXr5YsWcLrr7/O0qVLmTx5MosWLQJg4MCBLFq0iC+++II2bdrw/PPPc+KJJ3LOOecwevRoli1bxpFHHlnsflJURPdyAS/UFeAilSsc168++ugjBgwYQPXq1QE455xzAFixYgV33303WVlZ7Ny5k759+xb7/LLuV5VF9Bm6iPgjXNeviuuWd8UVV/DUU0+xfPlyRowYUeJI1rLuV5Up0EXkAOG4ftW9e3fefPNNsrOz2bFjB++88w4AO3bsoFGjRuTm5vLKK68U7l+zZk127NhReL+k/eR3Ed/kIiKVb28zZyhHaXfu3JkLL7yQtLQ0mjdvzsknnwzAAw88wPHHH0/z5s1p3759YYhfdNFFDBkyhDFjxvDGG2+UuJ/8zpxzlXawLl26OC1wIeKPVatW0aZNG7/LkFIU93MysyXOuS6lPVdNLiIiMSI6Av23HbB+gd9ViIhEtOgI9HdvgVcvhN3b/K5ERCRiRUegn3wL7NkB8x7xuxIRkYgVHYFevw10uhQWPQdb1/pdjYhIRCo10M0s2cw+M7MvzGylmd0X2F7XzGaa2ZrA1zphrfTUOyE+CWbfF9bDiIhEq7Kcof8G9HTOdQTSgNPNrBswHJjtnGsFzA7cD5+aDSH9JvjqLfj+07AeSkSiy9y5c/nkk08q/Drr1q2jXbt2pe738MMPF7l/4oknVvjYoVBqoDvPzsDdxMDNAecCLwa2vwj0D0uF+zrxeqjREN6/Gyqx/7yIRLZQBXpZ7R/olXnsgylTG7qZxZvZMmATMNM59ynQwDn3I0Dga/3wlRmQdAj0vAs2fOadqYtI1Hj55Zfp2rUraWlpXH311eTn57N+/XpatWrFli1bKCgo4OSTT+b9998HoH///hx77LEcc8wxjBs3rvB1pk+fTufOnenYsSO9evVi3bp1jB07lscff5y0tDQ++uijIsedN28eaWlppKWl0alTJ3bs2IFzjmHDhtGuXTvat2/PxIkTD6h3woQJXH/99YX3+/Xrx9y5cxk+fDjZ2dmkpaUxePBgAGrUqAFQ4uvOnTuXHj16cP7553P00UczePBgwjGos0xD/51z+UCamdUG3jSz0j+TBJjZUGAoQLNmzYIqsoi0wbDwGZg1AlqfCQlJB+yihaVFSvHecPhpeWhfs2F7OGNUsQ+tWrWKiRMnMn/+fBITE7n22mt55ZVXuOyyy7j99tu55pprOP7442nbti19+vQBYPz48dStW5fs7GyOO+44zjvvPAoKChgyZAgffvghLVu2ZNu2bdStW5drrrmGGjVqcNtttx1w7EcffZSnn36a9PR0du7cSXJyMpMnT2bZsmV88cUXbNmyheOOO47u3buX6W2OGjWKp556imXLlh3w2MFed+nSpaxcuZLGjRuTnp7O/PnzOemkk8r6r1sm5erl4pzLAuYCpwM/m1kjgMDXTSU8Z5xzrotzrktqamoFywXi4qHPA/DLOq/Xy37CNTG/iARv9uzZLFmyhOOOO460tDRmz57Nt99+C8Cf//xnduzYwdixY4ssSTdmzBg6duxIt27d+OGHH1izZg0LFy6ke/futGzZEoC6deuWeuz09HRuueUWxowZQ1ZWFgkJCXz88ccMGjSI+Ph4GjRowCmnnFI4P3tFHOx1u3btStOmTYmLiyMtLY1169ZV+Hj7K/UM3cxSgVznXJaZpQC9gX8AbwOXA6MCXyuvDeSo3nBkT5j3D0gbBCm/d7A52MT8OksXCSjhTDpcnHNcfvnljBw58oDHdu/ezYYNGwDYuXMnNWvWZO7cucyaNYsFCxZQvXp1evToQU5ODs65YqfgPZjhw4dz1llnMW3aNLp168asWbPK1NyRkJBAQUFB4f2yTNd7sNetVq1a4ffx8fHk5eWV+nrlVZYz9EbAB2b2JbAIrw39XbwgP83M1gCnBe5XntMegJzt8GHRRWa1sLRI5OnVqxdvvPEGmzZ5H+S3bdvG+vXrAbj99tsZPHgw999/P0OGDAFg+/bt1KlTh+rVq7N69WoWLlwIwAknnMC8efP47rvvCl8HDpxqd19r166lffv23H777XTp0oXVq1fTvXt3Jk6cSH5+Pps3b+bDDz+ka9euRZ7XokULli1bRkFBAT/88AOfffZZ4WOJiYnk5uYecKyyvG44lXqG7pz7EuhUzPatQK9wFFUmDdt57emfjYOuQ6BOC8CbgD+zmPDWwtIi/mnbti0PPvggffr0oaCggMTERJ5++mnWrVvHokWLmD9/PvHx8UyaNIkXXniBiy++mLFjx9KhQwdat25Nt27dAEhNTWXcuHEMHDiQgoIC6tevz8yZMzn77LM5//zzeeutt3jyyScLp+YFeOKJJ/jggw+Ij4+nbdu2nHHGGSQlJbFgwQI6duyImfHII4/QsGHDIs0g6enptGzZkvbt29OuXTs6d+5c+NjQoUPp0KEDnTt3LjI3+4ABA4p93dWrV4f/H5lonz73140wpjMcfSacPx74vQ1932aXlMR4rUUqVZ6mz40OVXf63EMbw4k3wIpJsMH7Q6GFpUWkqor+FYvSb4QlE7zBRle+B2ZaWFpEqqToPkMHqFYTTr0Dvl8Aq6f6XY1IRKvMJlYpv4r+fKI/0AE6XQb1WsPMv0P+gVeeRQSSk5PZunWrQj1COefYunUrycnJQb9G9De5AMQneIONXv0jLH4Bjh/qd0UiEadp06Zs2LCBzZs3+12KlCA5OZmmTZsG/fzYCHSAVn2gZXeYOxI6XgjJtfyuSCSiJCYmFo6wlNgUG00uAGbQ50HI/gU+eszvakREKl3sBDpAo47Q4UJv8q6sH/yuRkSkUsVWoAP0vNs7W5/zgN+ViIhUqtgL9NqHQ7dr4cuJsHGp39WIiFSa2At0gJP+CtXrwfQ7tLKRiFQZsRnoyYdC7xHeYKMv/+d3NSIilSI2Ax0g7RJo3Blm3gM5v/pdjYhI2MVuoMfFwVmPws5N3kIYIiIxLnYDHaDJsdD5Mvh0LGyqnPmIRUT8EtuBDtBrBCTVgPeG6QKpiMS02A/0Qw7z+qZ/9yF8NcXvakREwib2Ax2gy1XQsD3MuAv27PK7GhGRsKgagR4XD2c+Cr9mHrCotIhIrKgagQ7QrBt0HASfPAlbvvG7GhGRkKs6gQ7Q+z5ITIHpt+sCqYjEnKoV6DUbQI874JtZkPGe39WIiIRUqYFuZoeb2QdmtsrMVprZTYHt95pZppktC9zODH+5IdB1CKS28c7Sc7P9rkZEJGTKcoaeB9zqnGsDdAOuM7O2gcced86lBW7TwlZlKMUnwpmjIet7mP8vv6sREQmZUgPdOfejc+7zwPc7gFVAk3AXFlYtT4Z258HHj8Mv6/yuRkQkJMrVhm5mLYBOwKeBTdeb2ZdmNt7M6pTwnKFmttjMFkfU4rSnPQAWD9Pv9LsSEZGQKHOgm1kNYBJws3PuV+AZ4EggDfgR+Gdxz3POjXPOdXHOdUlNTQ1BySFSqwmcMgwypsKamX5XIyJSYWUKdDNLxAvzV5xzkwGccz875/KdcwXAs0DX8JUZJt2ug8OOgvduh7zf/K5GRKRCytLLxYDngVXOucf22d5on90GACtCX16YJSTBGY/AtrWw4Gm/qxERqZCEMuyTDlwKLDezZYFtdwKDzCwNcMA64OqwVBhuR/WCo/vBh6Ohwx+hVlO/KxIRCUqpge6c+xiwYh6Kjm6KZdH3YXi6K7x/N1wwwe9qRESCUrVGipakTnM4+VZY+SasneN3NSIiQVGg73Xijd4F0ndu1hS7IhKVFOh7JSbD2WMgaz188LDf1YiIlJsCfV8t0uHYK2HhvyFzid/ViIiUiwJ9f6fdBzUawNs3Qn6u39WIiJSZAn1/ybXgrH/Czys0eZeIRBUFenGOPgva9od5j8CWNX5XIyJSJgr0kpzxiLe60ds3QkGB39WIiJRKgV6Smg2g70Pw/Sfw+QS/qxERKZUC/WDSBkPLU2DmCPh1o9/ViIgclAL9YMzg7Ce83i5Tb9XC0iIS0RTopal7BJx6J2RMg6+m+F2NiEiJFOhl0e1aaJQG04bB7m1+VyMiUiwFelnEJ8A5T3ph/v49flcjIlIsBXpZNeoA6TfBspdh7Qd+VyMicgAFenmc8jeoeyS8cxPs2e13NSIiRSjQyyMxBc7ZOyPjQ35XIyJShAK9vFqcBMdeEZiR8XO/qxERKaRAD8Zp98Mh9eHtGzQjo4hEDAV6MPadkfGTMX5XIyICKNCDNiWnE3PiTuC3WQ9zycMvMGVppt8liUgVl+B3AdFoytJM7pi8nOq5lzG92kruynmMCyfXA6B/pyY+VyciVVWpZ+hmdriZfWBmq8xspZndFNhe18xmmtmawNc64S83MoyekUF2bj5bqcXtuUNpE/c917rXGT0jw+/SRKQKK0uTSx5wq3OuDdANuM7M2gLDgdnOuVbA7MD9KmFjVnbh93MKOvNqXk+Gxk/l8O3q9SIi/ik10J1zPzrnPg98vwNYBTQBzgVeDOz2ItA/XEVGmsa1U4rcfzDvEta7+jxRbSzkbPepKhGp6sp1UdTMWgCdgE+BBs65H8ELfaB+qIuLVMP6tiYlMb7w/m6SGe5uoL5tg2l/87EyEanKynxR1MxqAJOAm51zv5pZWZ83FBgK0KxZs2BqjDh7L3yOnpHBxqxsGtdOYVDfgcRl7YR5o6D16XDMAJ+rFJGqxlwZFm0ws0TgXWCGc+6xwLYMoIdz7kczawTMdc61PtjrdOnSxS1evDgEZUeo/Fx4vg9s+xauXQCHNva7IhGJAWa2xDnXpbT9ytLLxYDngVV7wzzgbeDywPeXA28FU2ikmLI0k/RRc2g5fCrpo+YE1688PhEGPgv5e2DKtVpcWkQqVVna0NOBS4GeZrYscDsTGAWcZmZrgNMC96PS3n7lmVnZOCAzK5s7Ji8PLtTrHQV9HoRvP4BFz4a8VhGRkpTahu6c+xgoqcG8V2jL8cfefuX7ys7NZ/SMjOAGCnW5Cr6eATP/7i0yXf/oEFUqIlIyDf2naL/ysmwvlZm3wlHSITB5COTtqUB1IiJlo0DnwH7lpW0vk5oN4Owx8NOXMHdk8K8jIlJGCnQO7FcOkJIYz7C+B+20U7o2/aDTJTD/CVi/oGKvJSJSCk3ORfH9yof1bR2aibZOHwXrPoY3r+bdE/+PkXMyQ38MERHK2A89VGK+H3pJvl+IG38Gkwu6c+ueoYWbUxLjGTmwvUJdRA4qZP3QJQSadeOl+AGcFzeXvnGfFW7e25NGRCQUFOiV5KFd57K8oAUjE58jlV8Ktwfdk0ZEZD8K9EqSWrsmN+deR3V+45+JY4nDG0VaoZ40IiL7UKBXkmF9W7MxoRn35l1O9/jl3BD/Zmh60oiIBKiXSyUp7EkzPZFJuzO4KXEyXU84nXRdEBWREFEvFz/s2QXP9YadP8PVH0EthbqIlEy9XCJZ0iHwx5cg7zf4vyu8aXdFRCpIge6Xeq28+V42fAYzR/hdjYjEAAW6n9oNhK5Xw8Kn4auonk5eRCKAAt1vfR6EJl1gynWwda3f1YhIFFOg+y0hCS6YAPEJ8L/LIFcDjUQkOAr0SFD7cBj4HPy8Eqbd5nc1IhKlFOiRolVv6D4Mlr4Mn//X72pEJAop0CNJj+HeknXTboOflvtdjYhEGQV6JImLh/Oeh5Q6Xnt6zna/KxKRKKJAjzQ1UuH8F+CX9fDW9VCJI3lFJLop0CNR8xPgtPtg1duw8Bm/qxGRKKFAj1QnXA9H94OZ98D3n/pdjYhEgVID3czGm9kmM1uxz7Z7zSzTzJYFbmeGt8wqyAzOfRpqHe7N97Jzs98ViUiEK8sZ+gTg9GK2P+6cSwvcpoW2LAEgpbY3iVf2Npg4GHJz/K5IRCJYqYHunPsQ2FYJtUhxGnWAAf+BHz6Fd27URVIRKVFF2tCvN7MvA00ydUJWkRzomP7Q8x74ciJ8+Kjf1YhIhAo20J8BjgTSgB+Bf5a0o5kNNbPFZrZ482a1Awft5Fuhw0XwwYOwYrLf1YhIBApqCTrn3M97vzezZ4F3D7LvOGAceCsWBXO8qmzK0kxGz8hgY1Y2zWsN5I26X1Nvyl+gdjNoWuoCJiJShQR1hm5mjfa5OwBYUdK+ErwpSzO5Y/JyMrOyccC67Xmcvfkv7EpKhdcGQdb3fpcoIhGkLN0WXwMWAK3NbIOZ/Ql4xMyWm9mXwKnAX8NcZ5U0ekYG2bn5Rbb9mHsIQ/KGecvXvXoR/LbDp+pEJNKU2uTinBtUzObnw1CL7GdjVvFzoy/4tR4MmQAvnw9v/AkGvebNAyMiVZpGikawxrVTSt5+ZE84czSsmQHv313JlYlIJFKgR7BhfVuTklj0zDslMZ5hfVt7d477Exz/F1j4b1ikD00iVV1QvVykcvTv1ASgsJdL49opDOvbunA7AH0fgm3fwrRhULeld+YuIlWSuUocedilSxe3ePHiSjtelfHbDni+L2zfAH+eCamt/a5IRELIzJY450rtp6wml1hQrSZc/DokVINXLoBdW/yuSER8oECPFbWbeb1ddv4MEy/xujWKSJWiQI8lTbtA/2fg+wXw9g1QUOB3RSLy2w5452bYHf45DhXosabdQDj1bm8ir/fv0uyMIn7KzfFGdX/+EmR+HvbDKdBjzJSlmaR/0pkX8vrCwn+z+vU7/S5JpGrKz/UWp1n3kffJuVXvsB9SgR5DCud+2Z7D/XmX8r+8Uzg649+s+L8H/S5NpGopyIc3r4Gv32NZh7+T/l4qLYdPJX3UHKYszQzbYdUPPYbsO/eLI47heUOobjn0WzkaWjaBLlf6XKFIFeAcvPtXWPEGK9reyqClx5Cd603jkZmVzR2TlwMUHU8SIjpDjyH7z/1SQBx/zb2OD/I7er9gy9/wqTKRKsI5byqOz1+Ek2/j6m9POmCCvezcfEbPyAjL4RXoMaS4uV9ySeD+6ndA83SYPBQy3vOhMpEqYt4jsOAp6Ho19Ly7xAn2StpeUQr0GFLS3C83nd7B66PeqCP873L4dp5PFYrEsAVPw9yHIW0wnD4KzA4+wV4YKNBjSP9OTRg5sD1NaqdgQJPaKYwc2N5rq0s+FC6ZBIcd6XWj+mGR3+WKxI4lL8KMO6HNOXD2GIjzorXUCfZCTHO5VDU7foIXzoDdW+GKqdCwfZFl7oqdAExESrZikrcuwVG94KJXvSk49hGK/19lnctFgV4VZX0P40+H/D3MOn4CN8zcWeTCTUpi/O9n9iJSsq9nwOsXQ9Ou3ifgpOphOYwm55KS1W4Gl70FztF+zmXUzf2pyMPhvAovEjO++xAmXgoN28PFE8MW5uWhQK+q6rWCy6aQ7LJ5OelhUskq8nC4rsKLRLMpSzNJHzWHAXf8i90vXsCv1Q+HSyZ716gigAK9KmvYnmHV7qG+ZfHfpJHUYmfhQ+G6Ci8SrfaOxK65PYMXkv7BpoJa9Mu6jSkZkXPyo0Cv4s4841yuLxhGS/uJV5Ie5jC2h/UqvEi0Gj0jgyPy1vJy0sNkU41Lcu/k+9xDI6p5UkP/qzjvwufF3DHNeGjPP3gz5QFW9X6JvrogKlJEk+1LeS5pNDuoziV77mSDSwUiq3lSZ+hC/05NeOyu20i56m2aJe2i76dXwJZv/C5LJHKsmclL1Uax2dXm/N/u5TvXqPChSGqeLDXQzWy8mW0ysxX7bKtrZjPNbE3ga53wlimVovkJcMW7kJcD4/vCj1/4XZFI2Oy9wFnqLIgrJsFrF5FT6ygudffzI4cVPhRpzZNlOUOfAJy+37bhwGznXCtgduC+xIJGHeCqGZCQDBP6wfoFflckEnKFU01nZeP4fRbEA0J9yQRv0FDTrtT+y3T+NjC9+JHYEaJMA4vMrAXwrnOuXeB+BtDDOfejmTUC5jrnSv0zpYFFUWT7Bnipv/f1wpcrZXJ+kcqSPmoOmcW0fTepncL84T29Ox8/AbNGQKs+cMGLvvYzD/fAogbOuR8BAl/rH6SQoWa22MwWb968OcjDSaWr1RSufM/rr/7aRbBist8ViYTMQWdBdA5m3euFebvz4MJXImLQUFmE/aKoc26cc66Lc65LampquA8noVQj1WtTb9oF3rjK+/gZUOb2R5EIVNKFzCa1qsHUW+Djx+HYK2Hgs5CQVMnVBS/YQP850NRC4Oum0JUkESW5ljcS7qje8M5NMP9fZW9/FIlQxc2CWDPR8Vq98bB4PJz0V+j3OMTFl/AKkSnYQH8buDzw/eXAW6EpRyJSUnVvFrljBsLMv7Pj3bvJzs0rsovmf5Fosv9U00fUimNm4/9weOZU6H2vdzPztcZglDqwyMxeA3oA9cxsAzACGAX8z8z+BHwPXBDOIiUCJCTBec9B8qFcumQClrCde/KuxO1zThBJAyxEStO/UxOvh0rOdnj1Ivh+gXdW3uUqv0sLWqmB7pwbVMJDvUJci0S6uHjo9wT/XbadS3mTmpbNrbnXkBf4NYqkARYiZbJrC7w8EH5e6Z2wtD/f74oqRCNFpXzMqNnvIf5ZcDHnxn/Cy0kjqcuvETfAQqRUP62AZ0+FzRlw0WtRH+agQJcg9O/UhCMH3M19iTeTZt8wLfke/t0zPqIGWIgc1Mo34fnTID8XrpgGf+jjd0UhoRWLpGI2LoXXB8PubXDuUyWe5WiZO4kIBQXwwUPw0aPeKkMX/hdqNvS7qlJpxSKpHI07wdC50DgNJv0JZv4dCvKL7KJujhIRcrbD64O8MO98mTfGIgrCvDwU6FJxNerDZW9Dlz/B/H/Bq3+E7F8KHx49I6PImqWgbo4SWqUOdNuyBp7tBd/MgjMfhbPHHLCYcyxQoEtoJCRBv8eg3xPw7Tx4tidsWg2UMsxapIJK/QT49fve72P2Nm8t3a5DorKPeVko0CW0ulwJl78Dv+2E53rD6mkldmdUN0cJhRI/AU5fDR895n1irNPCaxpscZIfJVYaBbqEXvMTvP889Y6C1wfxXIs5VE8sekakbo4SKsV90kshh+G7H4HZ90G7gd6U0LWb+VBd5VKgS3jUauLN1tjhItqsfpJZTcdzVC0idh5piV77f9JrapuZlHQfZ8V/Cr3vg/Oej5rZEitKa4pK+CSmwICx0KgDjd+/m1mpmXD1K1D3CL8rkxgyrG9r7pi8nOzcfLrFfcXTif8ikXwWdhvLiSdd5Hd5lUqBLuFlBidcB/XbwP9dCeN6wFmPxcSoPCmbcI9B6N+pCVaQx6b3RnJl7kR+iGvMNz3HcdrJ6SE7RrTQwCKpPNu+g8lDYMMiOGaAF+zV6/pdlYTR3h4o+160TEmMD22T2+av4c2rYePn0P4C7/cq+dDQvHaE0MAiiTx1W8KV06HnPbDqXfh3N/h6ht9VSRiFdQxCQQEs+Df852T4ZZ23TFxgRtCqSoEulSs+AbrfBkM/gOr1vC5lb10POb/6XZmEQdjGIPyyHl46B2bcAUf0gGsXwjH9K/aaMUBt6OKPhu29UJ870htd+t086P9MkX7Cmv8l+jWunVLsYsxBj0FwDpa+DNPv8O6f8xR0uiRmBwqVl87QxT8J1byVYa6cDhYPE/rBjLsgN0fzv8SI4pZ6C3oMwo6f4bVB8Pb10Kgj/GU+dL5UYb4PnaGL/5odD9d87E3steApWDOTd379M9m5Rc/G97a96iw9euz9WZX3k9b+n84eb7+OrisegNzd0HckHH8NxOl8dH/q5SKR5ZvZ8Nb15P36E0/l9+epvP6FKyKBNzDpu1Fn+VefhN2+PWMOZSf3Jb7IgPj5/FK7HXUGj4fUqjfCWL1cJDod1Quu/YRZ8Sdzc8JkJieNoJVtKHxY87/Evr09Y06J+4IZ1YbTL24hj+Wez7nZI6pkmJeHAl0iT0odcs5+hhvzb6Gpbea9pOHcn/ACjRN3af6XKiB5+zc8lziaF5P+wQ6XwoA99zEmfyA/bM/1u7SIpzZ0iUheG+ufuXR6By7c9QoXJ8xmUMICEncNg9xrIDG52OepZ0wU27UF5o5kRrXx7HbVGJk7iAn5ffmNJECfzspCbegSHTZneBdNv54OtZpB7xHQ7rwiPRwqZVSihF5uDnw6Fj76J+zZxbfN/8gla3uyMfeQwl2q+s+xUtrQzWydmS03s2VmpqSW8EltDRdP9BYoSK7lLXf3/Gnww2eFu2hlpCjjHKyYBE8fB7MmecMyAAAJ+UlEQVRGQLMT4NoFHHHFWP42MJ0mtVM0O2c5haLJ5VTn3JYQvI5I6Y7oAVfPgy9eg9kPeKF+zADofa9WRoomP3wGM+705vVp0N77Q31Ej8KH+3dqogAPgtrQJfrExXujA9v2h0+e9Eaarp7KQ4ecwahd/fiVQ4rsrrbXCPLLOph1L6x8E2o09EZ6pl3s/Uylwiray8UB75vZEjMbGoqCRMqsWg049Q640Ztlb1D+28yrdguXxc8gkTxAKyNFjOxf4P174KnjIGM6nHI73LDEG+mpMA+ZCl0UNbPGzrmNZlYfmAnc4Jz7cL99hgJDAZo1a3bs+vXrK1KvSMl+/ILNk4aRuuVTNrnavJV4Bk16X8uZ3TqE9DCx1JMm7O9lyzfeBc9lr3qjPNMuhp53w6GNQ3eMKqCsF0VD1svFzO4FdjrnHi1pH/VykbBzDtbOhoXPwDezIL4adLgAjv8LNGxX4ZcPtidNJP4RCFuvIOfg2w+8n8Ga9yE+Cdqd7y10EoKfQVUU9kA3s0OAOOfcjsD3M4H7nXPTS3qOAl0q1eYM+PQ/3gXU3N3Q4mTo9hf4w+lBf8xPHzWn2NkDm9ROYf7wnsU+J1K7UwbzXg5qz274cqJ3Rr55NRxSH477E3S5CmrUD0HFVVdZA70iF0UbAG+a1w84AXj1YGEuUulSW0O/x6DXPfD5S/DpOHj9YqjTwpvcKW1wuRdDCKYnzcG6U/oZ6MG8l2I/aRwBLHoWlkzw2sobdoD+Y6HdQG9GTak0QQe6c+5boGMIaxEJj5Q6kH4TdLsOVr8DC8fC9OEw5yGvt8zxQ8u8cHUw83tHanfK8r6X/T9p1N/+JUlvPkJB/GfE4eDoft4noGYnaEpbn6jbolQd8Qlen/VjBkDm517TwKLnvK+t+njb/9D3oOuc7rvC/F6l9aQJ+SIPJShvO31538voGRkclvcTp8Uv5tz4T0iLW8uvrjqvx/Xj4usfgDrNQ/p+pPwU6FI1NekMA8fBaffDoue9VXDWzPAW2mh+Ihx9FrQ+84CQCmZ+72D+CJTX/mfPexcE2bfm/ZXpvTgHPy2H1VMZt/s1jqnm9VJbXXA49+RewaT87mT/lszFCvOIoLlcRMALro1LYfVU77Z5lbe9YXtofZYX8A3bB92UEO5eLiG9wJmfB99/Evi3mAbbvweML+xo3tnTmZkFx7LeNazYMaRcKuOiqEjsMPPO2pt09i6ibl0LGdO8UJv3D5g3ypsU7OgzvXBvdqLXhFNG4R7KXuF2+j27vMVFVk/1Pqlk/wIJyXDEqXDK3+APp/Pdmj28Mnk52fnh+6QhFaNAFynOYUfCiTd4t52bvVkeV0+FxS94be4pdbxukI06eL06GraHmo18uxhYrnb6ggLY9i38vNxrTtm4DNbPh7wcSK4Nrc/w/mgd2ROSfp9GoX8n72uk9aeX36nJRaQ89uyCtXO8cP9+Ifzy3e+PVT8MGrTzwn1vyNdrBfGJYS+rpL7uj5xzFGc3yvKCe+/t55WQu8vbKS4B6rWGlt0DnzxOKNcnD6kclT5StCwU6BJzcn6FTV8FwvLLQGB+Bfm/eY/HJ0H9Nl64N2gPhzbypv9NrgXVDvXOiJMPDS70836DnO2Ft09WrmX64gyq787k2Gob6FZ9IzV3rQNX4O1f7dDAH5t9bqlHq694FFCgi/glPw+2rtkn5Fd4X3dvLfk5iYf8HvT73hJTYM/OIsFdeMvLKfn1ajU7MLxrN1P/8Cili6IifolP8M7K67eBDn/0tjkHOzfB7i3Fh3POdsjJ+v37nT/BlgxvOH21mr8H/KFNign+2kXv12zgtfFLlaNAF6kMZl7Q1mxQpt0P6OZ4si4+SukU6CIRJphBQiJQ8QUuRCTEtDaqBEuBLhJhInUyL4l8CnSRCFPSpF1aG1VKo0CXKm3K0kzSR82h5fCppI+aw5SlmX6XxLC+rUlJLLoAh4bYS1nooqhUWZF68TGYGR1FQIEuVVikriQE4Z/MS2KTmlykytLFR4k1CnSpsnTxUWKNAl2qLF18lFijNnSpsnTxUWKNAl2qNF18lFiiJhcRkRhRoUA3s9PNLMPMvjGz4aEqSkREyi/oQDezeOBp4AygLTDIzNqGqjARESmfipyhdwW+cc5965zbA7wOnBuaskREpLwqEuhNgB/2ub8hsK0IMxtqZovNbPHmzZsrcDgRETmYivRyKW5xwgMWKHXOjQPGAZjZZjNbH+Tx6gFbgnxupNF7iTyx8j5A7yVSVeS9NC/LThUJ9A3A4fvcbwpsPNgTnHOpwR7MzBaXZZHUaKD3Enli5X2A3kukqoz3UpEml0VAKzNraWZJwEXA26EpS0REyivoM3TnXJ6ZXQ/MAOKB8c65lSGrTEREyqVCI0Wdc9OAaSGqpTTjKuk4lUHvJfLEyvsAvZdIFfb3Ys4dcB1TRESikIb+i4jEiKgKdDN7wMy+NLNlZva+mTX2u6ZgmdloM1sdeD9vmlltv2sKhpldYGYrzazAzKKyN0KsTGFhZuPNbJOZrfC7loows8PN7AMzWxX43brJ75qCZWbJZvaZmX0ReC/3hfV40dTkYmaHOud+DXx/I9DWOXeNz2UFxcz6AHMCF5f/AeCcu93nssrNzNoABcB/gNucc4t9LqlcAlNYfA2chtcVdxEwyDn3la+FBcHMugM7gZecc+38ridYZtYIaOSc+9zMagJLgP5R+jMx4BDn3E4zSwQ+Bm5yzi0Mx/Gi6gx9b5gHHEIxA5mihXPufedcXuDuQrx+/FHHObfKOZfhdx0VEDNTWDjnPgS2+V1HRTnnfnTOfR74fgewimJGoUcD59kZuJsYuIUtt6Iq0AHM7CEz+wEYDPzd73pC5CrgPb+LqKLKNIWF+MPMWgCdgE/9rSR4ZhZvZsuATcBM51zY3kvEBbqZzTKzFcXczgVwzt3lnDsceAW43t9qD6609xLY5y4gD+/9RKSyvI8oVqYpLKTymVkNYBJw836fzqOKcy7fOZeG9ym8q5mFrTks4lYscs71LuOurwJTgRFhLKdCSnsvZnY50A/o5SL4YkY5fibRqNxTWEj4BdqbJwGvOOcm+11PKDjnssxsLnA6EJYL1xF3hn4wZtZqn7vnAKv9qqWizOx04HbgHOfcbr/rqcI0hUWECVxIfB5Y5Zx7zO96KsLMUvf2YDOzFKA3YcytaOvlMglojderYj1wjXMu09+qgmNm3wDVgK2BTQujsceOmQ0AngRSgSxgmXOur79VlY+ZnQk8we9TWDzkc0lBMbPXgB54s/r9DIxwzj3va1FBMLOTgI+A5Xj/1wHuDIxMjypm1gF4Ee93Kw74n3Pu/rAdL5oCXUREShZVTS4iIlIyBbqISIxQoIuIxAgFuohIjFCgi4jECAW6iEiMUKCLiMQIBbqISIz4f+z7PqheWEItAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(x,y, 'o', label='data')\n", "plt.plot(x,2 * x ** 2 - 3 * x + 4, '-', label='exact solution')\n", "plt.legend(loc='upper right')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "array([[ 9. , -3. , 1. ],\n", " [ 7.5625, -2.75 , 1. ],\n", " [ 6.25 , -2.5 , 1. ],\n", " [ 5.0625, -2.25 , 1. ],\n", " [ 4. , -2. , 1. ],\n", " [ 3.0625, -1.75 , 1. ],\n", " [ 2.25 , -1.5 , 1. ],\n", " [ 1.5625, -1.25 , 1. ],\n", " [ 1. , -1. , 1. ],\n", " [ 0.5625, -0.75 , 1. ],\n", " [ 0.25 , -0.5 , 1. ],\n", " [ 0.0625, -0.25 , 1. ],\n", " [ 0. , 0. , 1. ],\n", " [ 0.0625, 0.25 , 1. ],\n", " [ 0.25 , 0.5 , 1. ],\n", " [ 0.5625, 0.75 , 1. ],\n", " [ 1. , 1. , 1. ],\n", " [ 1.5625, 1.25 , 1. ],\n", " [ 2.25 , 1.5 , 1. ],\n", " [ 3.0625, 1.75 , 1. ],\n", " [ 4. , 2. , 1. ],\n", " [ 5.0625, 2.25 , 1. ],\n", " [ 6.25 , 2.5 , 1. ],\n", " [ 7.5625, 2.75 , 1. ],\n", " [ 9. , 3. , 1. ]])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x_mat = np.column_stack( (x**2, x, np.ones(len(x))) )\n", "x_mat" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "scrolled": true, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1.98114133 -2.87140368 4.227188 ]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n", "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n", " \"\"\"Entry point for launching an IPython kernel.\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XlcVNX/x/HXYRhWWQVcQEXLFBVFcs2tssRvi5lWrpUtZvUrs9KvVraYmZa22fI1Tb+WmtlXTc191zJ3cVfccgMXQEFA9jm/P0CKBEFguDPweT4e8xjmzp17Pxfy3Z1zzzlXaa0RQghh/xyMLkAIIUTZkEAXQogKQgJdCCEqCAl0IYSoICTQhRCigpBAF0KICkICXQghKggJdCGEqCAk0IUQooJwLM+d+fn56eDg4PLcpRBC2L2dO3fGaa39i1qvXAM9ODiYHTt2lOcuhRDC7imlThVnPWlyEUKICkICXQghKggJdCGEqCDKtQ1dCGGczMxMzp49S1pamtGliEK4uLgQFBSE2Wwu0ecl0IWoJM6ePYuHhwfBwcEopYwuR/yD1pr4+HjOnj1L3bp1S7QNmw/0BZHRjF8RRUxCKjW9XRkW0YDuzQONLksIu5OWliZhbsOUUlStWpXY2NgSb8OmA31BZDRvzN9HamY2ANEJqbwxfx+AhLoQJSBhbttK+/ex6Yui41dE5YX5NamZ2YxfEWVQRUIIYbtsOtBjElJvarkQwn689957TJgwodD3FyxYwMGDB8uxIvtn04Fe09v1ppYLIcrOgsho2o1bS90RS2g3bi0LIqPLd/8S6DfNpgN9WEQDXM2mfMtczSaGRTQwqCIhKodr16+iE1LR/HX9qrShPmbMGBo0aMA999xDVFRO0+mUKVNo2bIlzZo1o2fPnly9epU//viDRYsWMWzYMMLCwjh+/HiB64n8bDrQuzcPZGyPUAK9XVFAoLcrY3uEygVRIazMGtevdu7cyU8//URkZCTz589n+/btAPTo0YPt27ezZ88eQkJCmDp1KnfccQfdunVj/Pjx7N69m1tuuaXA9UR+Nt3LBXJCXQJciPJljetXv/32Gw8//DBubm4AdOvWDYD9+/czcuRIEhISSE5OJiIiosDPF3e9ysymz9CFEMaw1vWrgrrlDRgwgK+++op9+/bx7rvvFjqStbjrVWYS6EKI61jj+lXHjh355ZdfSE1NJSkpiV9//RWApKQkatSoQWZmJrNmzcpb38PDg6SkpLzXha0n/mLzTS5CiPJ3rZmzLEdph4eH06tXL8LCwqhTpw4dOnQAYPTo0bRu3Zo6deoQGhqaF+K9e/dm4MCBTJw4kblz5xa6nviL0lqX285atGih5QYXQhjj0KFDhISEGF2GKEJBfyel1E6tdYuiPitNLkIIUUHYR6CnJ8GpzUZXIYQQNs0+An3xa/BjL7h6yehKhBDCZtlFoGe2G0yUToMNHxtdihBC2Cy7CPR3j87muaAgkndMhfjjRpcjhBA2qchAV0q5KKW2KaX2KKUOKKVG5S73VUqtUkodzX32sVaR/Rr147LO4jtvL1gzylq7EUIIu1acM/R04G6tdTMgDOiqlGoDjADWaK3rA2tyX1tF46qNefCWB5nh6U70kcVwequ1diWEsEPr16/njz/+KPV2Tp48SZMmTYpc78MPP8z3+o477ij1vstCkYGucyTnvjTnPjTwEPB97vLvge5WqTDX4OaDcTCZ+dy/OqwcCeXYf14IYdvKKtCL65+BXp77vpFitaErpUxKqd3ARWCV1norUE1rfQ4g9znAemVCNfdqPNXkaZa7mNgduxsOLrTm7oQQZWzmzJm0atWKsLAwBg0aRHZ2NqdOnaJ+/frExcVhsVjo0KEDK1euBKB79+7cfvvtNG7cmMmTJ+dtZ/ny5YSHh9OsWTM6d+7MyZMnmTRpEp999hlhYWH89ttv+fa7YcMGwsLCCAsLo3nz5iQlJaG1ZtiwYTRp0oTQ0FDmzJlzXb3Tp0/npZdeynv9wAMPsH79ekaMGEFqaiphYWH069cPgCpVqgAUut3169dz55138sgjj9CwYUP69euHNQZ1Fmvov9Y6GwhTSnkDvyiliv5Okksp9RzwHEDt2rVLVOQ1AxoPYO6RuYyv5sDM1e+gGtwHjk7XrSc3lhaiCMtGwPl9ZbvN6qHwr3EFvnXo0CHmzJnDpk2bMJvNvPjii8yaNYsnnniC4cOH8/zzz9O6dWsaNWpEly5dAJg2bRq+vr6kpqbSsmVLevbsicViYeDAgWzcuJG6dety6dIlfH19ef7556lSpQpDhw69bt8TJkzg66+/pl27diQnJ+Pi4sL8+fPZvXs3e/bsIS4ujpYtW9KxY8diHea4ceP46quv2L1793Xv3Wi7kZGRHDhwgJo1a9KuXTs2bdpE+/bti/vbLZab6uWitU4A1gNdgQtKqRoAuc8XC/nMZK11C611C39//1IV62Z2Y3D4YPaaLCzPiIXt3123jrUm5hdClNyaNWvYuXMnLVu2JCwsjDVr1nDixAkAnn32WZKSkpg0aVK+W9JNnDiRZs2a0aZNG86cOcPRo0fZsmULHTt2pG7dugD4+voWue927drx2muvMXHiRBISEnB0dOT333+nT58+mEwmqlWrRqdOnfLmZy+NG223VatWBAUF4eDgQFhYGCdPniz1/v6pyDN0pZQ/kKm1TlBKuQL3AB8Bi4AngXG5z+XSBtLtlm78eOhHPuMod234CJewPuD6VwebG03ML2fpQuQq5EzaWrTWPPnkk4wdO/a6965evcrZs2cBSE5OxsPDg/Xr17N69Wo2b96Mm5sbd955J2lpaWitC5yC90ZGjBjB/fffz9KlS2nTpg2rV68uVnOHo6MjFosl73Vxpuu90XadnZ3zfjaZTGRlZRW5vZtVnDP0GsA6pdReYDs5beiLyQnye5VSR4F7c19bnYNyYFjLYZwji5nOFtiY/yazcmNpIWxP586dmTt3Lhcv5nyRv3TpEqdOnQJg+PDh9OvXj/fff5+BAwcCkJiYiI+PD25ubhw+fJgtW7YA0LZtWzZs2MCff/6Ztx24fqrdvzt+/DihoaEMHz6cFi1acPjwYTp27MicOXPIzs4mNjaWjRs30qpVq3yfCw4OZvfu3VgsFs6cOcO2bdvy3jObzWRmZl63r+Js15qKPEPXWu8FmhewPB7obI2iitKyekvurnU337GB7ju+w6/VQPAJBnIm4I8uILzlxtJCGKdRo0Z88MEHdOnSBYvFgtls5uuvv+bkyZNs376dTZs2YTKZmDdvHv/973/p27cvkyZNomnTpjRo0IA2bdoA4O/vz+TJk+nRowcWi4WAgABWrVrFgw8+yCOPPMLChQv58ssv86bmBfj8889Zt24dJpOJRo0a8a9//QsnJyc2b95Ms2bNUErx8ccfU7169XzNIO3ataNu3bqEhobSpEkTwsPD89577rnnaNq0KeHh4fnmZn/44YcL3O7hw4et/0vGjqfPPXXlFN0XPET3pGTerdYRHpkG/NWG/vdmF1ezSe5FKio9mT7XPlTK6XPreNahd8M+zK/ixpGoRXA2538UcmNpIURlZdd3LHq+2fMsOr6QCQEBfLvyLdRTy0EpubG0EKJSstszdAAvZy9eaPYim51M/B67Gw4vMbokIYQwjF0HOkCvBr2o41GbCf7VyFr1NmRff+VZCCEqA7sPdLPJzOsthnLCpJmbeRF2/NfokoQQwhB2H+gAd9a6k1bVW/KNnx9XNoyFtESjSxJCiHJXIQJdKcXQFsNIQPOdswV++9TokoQQ/3BtAquyMn36dGJiYsp0m/auQgQ6QEjVEB669SFmentxZvtkSDhjdElCCCuylUC3xhD+kqowgQ7wcvOXcTQ58blPFVg72uhyhBCFGD9+PC1btqRp06a8++67ecsLmjI3OzubAQMG5E1J+9lnnzF37lx27NhBv379CAsLIzU1/+jwiRMn0qhRI5o2bUrv3r0BiI+Pp0uXLjRv3pxBgwZRp04d4uLirrupxYQJE3jvvfcAmDJlCi1btqRZs2b07NmTq1evAjBgwABee+017rrrLoYPH05KSgpPP/00LVu2pHnz5ixcmDO11YEDB/KmDG7atClHjx612u8U7Lwf+j8FuAXwVOgzfLP7GyKjfqF5zAtQ87pZC4So9D7a9hGHL5XtcPSGvg0Z3mp4keutXLmSo0ePsm3bNrTWdOvWjY0bN9KxY8cCp8w9efIk0dHR7N+/H4CEhAS8vb356quvmDBhAi1aXD+Acty4cfz55584OzuTkJAAwKhRo2jfvj3vvPMOS5YsyTfHemF69OiRN7/MyJEjmTp1Ki+//DIAR44cYfXq1ZhMJt58803uvvtupk2bRkJCAq1ateKee+5h0qRJvPLKK/Tr14+MjAyys7NvtLtSq1Bn6JAzZ3qAqz/j/f2xLB8hdzYSwsasXLmSlStX0rx5c8LDwzl8+HDemWtBU+bWq1ePEydO8PLLL7N8+XI8PT2L3EfTpk3p168fM2fOxNEx57x148aN9O/fH4D7778fH5+ib4O8f/9+OnToQGhoKLNmzeLAgQN57z366KOYTKa8Yxo3bhxhYWF5M0OePn2atm3b8uGHH/LRRx9x6tQpXF2tO6dUhTpDB3B1dGXI7a/y5u9vsvTiXh7Y+zM062V0WULYlOKcSVuL1po33niDQYMG5Vte2JS5Pj4+7NmzhxUrVvD111/z888/M23atBvuY8mSJWzcuJFFixYxevTovCAuaOrdG02TO2DAABYsWECzZs2YPn0669evz3vP3d093zHNmzePBg0a5Nt2SEgIrVu3ZsmSJURERPDdd99x9913F/1LKqEKd4YOcH+9+2lctTGf+geQvPptSLtidElCiFwRERFMmzaN5OScWxVHR0dz8eLFQqfMvXZ7up49ezJ69Gh27doFFD5l7rXpbu+66y4+/vhjEhISSE5OpmPHjnkzIy5btozLly8DUK1aNS5evEh8fDzp6eksXrw4b1tJSUnUqFGDzMzMfLMqFnRMX375Zd586JGRkQCcOHGCevXqMXjwYLp168bevXtL++u7oQp3hg45c6aPbDOSvkv68o05g39v+AgixhhdlhAC6NKlC4cOHaJt27ZATnfGmTNn0rVr1wKnzI2Ojuapp57KO4u+dpOMAQMG8Pzzz+Pq6srmzZvzmjOys7Pp378/iYmJaK159dVX8fb25t1336VPnz6Eh4fTqVOnvFtims1m3nnnHVq3bk3dunVp2LBhXq2jR4+mdevW1KlTh9DQ0ELnXH/77bcZMmQITZs2RWtNcHAwixcvZs6cOcycOROz2Uz16tV55513rPNLzWW30+cWx/ub32f+kbn8HHOB257ZAAENi/6QEBWUTJ+bX3BwMDt27MDPz8/oUvKplNPnFscr4a/g6eTJGD9f9LKhcoFUCFGhVehA93L24tUWr7PLyZFFsbvg4AKjSxJC2IiTJ0/a3Nl5aVXoQAd46NaHaObflE/9/Ehc+RZkpBhdkhCGKc8mVnHzSvv3qfCBnnOB9G0SFHzpePW6m0oLUVm4uLgQHx8voW6jtNbEx8fj4uJS4m1UyF4u/9TQtyF9Qvry46FZPLzjPzQO6wd+txpdlhDlKigoiLNnzxIbG2t0KaIQLi4uBAUFlfjzFbqXy98lZSTRbf4D1Ei6yEyXEBz6z4MCBhkIIYStkV4u/+Dh5MHrrYaxz8mReRe3QNQyo0sSQogyVWSgK6VqKaXWKaUOKaUOKKVeyV3+nlIqWim1O/dxn/XLLZ37695Pi4Db+aJqVS6vGA6ZqUV/SAgh7ERxztCzgNe11iFAG+D/lFKNct/7TGsdlvtYarUqy4hSirfajCTFwYHPTcmw6QujSxJCiDJTZKBrrc9prXfl/pwEHAICrV2Ytdzqcyv9Gz3BfI8q7N72JVw+aXRJQghRJm6qDV0pFQw0B7bmLnpJKbVXKTVNKVXgXJRKqeeUUjuUUjts5er6C81eIMDFjzE+nmQte8PocoQQokwUO9CVUlWAecAQrfUV4D/ALUAYcA74pKDPaa0na61baK1b+Pv7l0HJpedmduPfrUdw2MmROec3wtFVRpckhBClVqxAV0qZyQnzWVrr+QBa6wta62yttQWYArSyXpllr0udLrSt3pqvfH2JW/5vyEo3uiQhhCiV4vRyUcBU4JDW+tO/La/xt9UeBvaXfXnWo5TizTYjSXcw8YlDImz+2uiShBCiVIpzht4OeBy4+x9dFD9WSu1TSu0F7gJetWah1hDsFcyAJk+zuIo727d8BolnjS5JCCFKrNKMFC1MalYq3ec/gNuVGH72bIX5se+NLkkIIfKRkaLF5Oroyog2IzlmdmTW2TVwfK3RJQkhRIlU+kAHuKv2XXSq2Z5vfH04v/gVmWJXCGGXJNBzjWjzFtpk5gOnVPRauf+oEML+SKDnCvII4uXwIWxwc2XJ/ukQvdPokoQQ4qZIoP9Nv5B+NK3amHFVfYlb9BJkZxpdkhBCFJsE+t+YHEyMbv8hV02OjLVckMm7hBB2RQL9H+p51+P5sBdZWcWdNdu+gLijRpckhBDFIoFegKeaPEVDr1v4wNeTxEUvgcVidElCCFEkCfQCmB3MvN9hLJdNJj6+egR2TTe6JCGEKJIEeiFCqobwdOizLPKowu8bR8OVGKNLEkKIG5JAv4Hnmz1PvSpBjPJyJXnxECjHaRKEEOJmSaDfgJPJiVEdxnLB0cTn8dvg4AKjSxJCiEJJoBchLCCM/g37McfTg+2rhsPVS0aXJIQQBZJAL4aXb3+FINcA3qtiInXFm0aXI4QQBZJALwZXR1dGdRjLabMjX59ZBsfXGV2SEEJcRwK9mFrVaMWjtz7MDE9P9i4dDBlXjS5JCCHykUC/Ca+1/Df+Lj6845JBxtr3jS5HCCHykUC/CVWcqvBO+w847uTE5KjZEL3L6JKEECKPBPpN6hjUkQfrRDDVy4OoX1+UGRmFEDZDAr0Ehrd9Gy8nD952uEzWps+NLkcIIQAJ9BLxcvbi3ppDOOTsxNQdX9H/w/+yIDLa6LKEEJWcBHoJLIiMZuYaL9SVhnzr60F/y6e8M3+XhLoQwlBFBrpSqpZSap1S6pBS6oBS6pXc5b5KqVVKqaO5zz7WL9c2jF8RRWpmNknnH4FsV76tlsFz/Mj4FVFGlyaEqMSKc4aeBbyutQ4B2gD/p5RqBIwA1mit6wNrcl9XCjEJqQDo7CokxvTlmJMTV/w3UytRer0IIYxTZKBrrc9prXfl/pwEHAICgYeA73NX+x7obq0ibU1Nb9e8n7NTbsMS34bZXh708Z4MaYkGViaEqMxuqg1dKRUMNAe2AtW01ucgJ/SBgLIuzlYNi2iAq9mU9zol9n6c03351N+R+MVDDKxMCFGZFTvQlVJVgHnAEK31lZv43HNKqR1KqR2xsbElqdHmdG8eyNgeoQR6u6KAQC9PBoZ+SLLJzDuxv6H3zze6RCFEJaR0MW7aoJQyA4uBFVrrT3OXRQF3aq3PKaVqAOu11g1utJ0WLVroHTt2lEHZtmnWgRmM2/ExbyWm0nvARvCsaXRJQogKQCm1U2vdoqj1itPLRQFTgUPXwjzXIuDJ3J+fBBaWpFBbsSAymnbj1lJ3xBLajVtboi6IfRv1p71/OBM8XTj+y7Nyc2khRLkqTpNLO+Bx4G6l1O7cx33AOOBepdRR4N7c13ZpQWQ0b8zfR3RCKhqITkjljfn7bjrUlVKMvusT3B3dGJ5+nIytk6xTsBBCFKA4vVx+11orrXVTrXVY7mOp1jpea91Za10/99lub+VzrV/536VmZpeoX7mfqx+jO40nytmJiTs+gYuHy6pMIYS4IRkpyl/9you7vCgda3Wid72H+N7Tjc0LBkBWRimqE0KI4pFAJ3+/8uIsL47X246knksAb5kSubzmvRJvRwghiksCnev7lQO4mk0Mi7hhp50bcnF04aN7vybB0cy7J35Gn/yjtGUKIcQNORpdgC3o3jwQyGlLj0lIpaa3K8MiGuQtL6mGvg15JewlJuz+knlLB+HSYi5j10aX6T6EEOKaYvVDLysVvR96QSzawqBfe7M7/gDPnq3N2JSX8t5zNZsY2yNUQl0IcUNl1g9dlI6DcmDMPV+BNrO22nHuddiS915Je9IIIURBJNDLQYBbAEnnenHQ2ZlbAmbjz+W890rak0YIIf5JAr2cVDW1xiWhMT96uzDI8yscyBlFWpqeNEII8XcS6OVkWEQDkuP74pbpzqxqV3jGaU6pe9IIIcTfSaCXk+7NAxn78O2Q8DKJDo4cCdzKt+0uywVRIUSZkUAvR92bB7JlWH8+aDeKXS4ubDkxChLlPqRCiLIhgW6A+257mL7B9zHD3czyub0gO9PokoQQFYAEukGGtv+AMPdavONwmePLXjO6HCFEBSCBbhCzycwn903HzeTCkHOrSN43x+iShBB2TgLdQAFuAYzv/BVnzGbe2fQOOu6Y0SUJIeyYBLrBWga25dUmz7LK1Ynvf+kNmTLQSAhRMhLoNuCJ2wdzb9VmfGZOY/vCZ40uRwhhpyTQbYBSitER31LH7MXQK5Fc2PqN0SUJIeyQBLqNcDe78/l900kzOfL6nolkxkQaXZIQws5IoNuQej71eb/12+xxNjN+8ZOQlmh0SUIIOyKBbmMiGj7KE0H3MttZs3h+XyjH+eqFEPZNAt0GDbnrI8JdqjMq/SRRGz4wuhwhhJ2QQLdBZgcznzz4Ix4OTrx27EeuHF9ndElCCDtQZKArpaYppS4qpfb/bdl7SqlopdTu3Md91i2z8vFz82fCXV8Q4+jIyDUvYUm6YHRJQggbV5wz9OlA1wKWf6a1Dst9LC3bsgRAeO2OvN7wcdY5O/DN3O6QmWZ0SUIIG1ZkoGutNwKXyqEWUYB+rYfxsF843zoks3B+b7lIKoQoVGna0F9SSu3NbZLxKbOKRD5KKd7+13e0dqnBe6nH2LZiqNElCSFsVEkD/T/ALUAYcA74pLAVlVLPKaV2KKV2xMbGlnB3lZvZwcyn3f9HHZMbQ2KWc2L7t0aXJISwQSUKdK31Ba11ttbaAkwBWt1g3cla6xZa6xb+/v4lrbPSWhAZTbtxa2n27u/En3sVR2Xixb1fEHd8jdGlCSFsTIkCXSlV428vHwb2F7auKLkFkdG8MX8f0QmpaODMJU/Sop8mzuTA4HWDSY07YnSJQggbUpxui7OBzUADpdRZpdQzwMdKqX1Kqb3AXcCrVq6zUhq/IorUzOx8yy4m16fape7sd1S8ubAXFpkeQAiRqzi9XPporWtorc1a6yCt9VSt9eNa61CtdVOtdTet9bnyKLayiUkoeG70gxfbMrTew6x2zOKzn7uBJbvA9YQQlYuMFLVhNb1dC13+eIf36e3TjOn6Ej8veLycKxNC2CIJdBs2LKIBrmZTvmWuZhPDIhqglGL4A9Pp6BTAmCt7+W3tSIOqFELYCgl0G9a9eSBje4QS6O2KAgK9XRnbI5TuzQMBcHRwZHyPBTRQLgw99QtRe2YYW7AQwlBKl+PIwxYtWugdO3aU2/4qiwuXT9B3QXfQ2cy6ZwrVa99hdElCiFwWbWH+0fl0v7U7jg6OJdqGUmqn1rpFUevJGXoFUM2nHt/cPZFkpXhp1SBSEk4bXZIQAtBaM3brWEZtHsXq06utvj8J9AqiQZ07+ST8dY6ZNMMW9CArI8XokoSo9L6M/JKfon5iQOMBRNSJsPr+JNArkPbNnuLNOg/ym0pn3NyH0NnSnVEIo0zdN5Up+6bwyC0P8VrMaVTqZavvUwK9gnnsrrE85dmYOZkXmDj/EbTFYnRJQlQ6cw7P4fNdn/OvOhGMPBaJipwB0busvl8J9ApmQWQ0cw8/S2iiN99dPcaYHx81uiQhKpVfj//KmK1juDOwI2Oiz2A6+Tt0/w/Uv8fq+5ZAr0Cuzf0Sk5jO5pihNL1ShTnZRxj3Y3+jSxOiUlhzeg1vb3qbltVaMOHyVcxHl7O76Tu0W+ZP3RFLaDduLQsio622fwn0CuTvc79oHNkcPZzQJBdmZe7hxxWDDa5OiIptc8xmhm0YRuOqjZiY7oLzgfnsb/Q6fSIb502wF52Qyhvz91kt1CXQK5B/zv1iwcz2s28Qmmxm7Pl1zFv3pkGVCVGx7b64m1fWvUKwZzDfOATiHjkLOgxl0In2102wl5qZzfgVUVapQwK9Ailo7pdMnIm5Mop22oVRpxbx6+9jDKhMiIrr8KXDvLj6RQLcAphcpSleWydDq0Fw98hCJ9grbHlpSaBXIIXN/TKkazifP7KUltqJkcdms2LrZwZVKETF8mfinwxaNQh3J3emVG2P32+fQVg/6DoOlLrhBHvWIIFegdxo7heXKv582fNXmllMjDg0lXW7JhtdrhB2LSY5hoErBwIwpeZ91Fg7FkK6wYMTwSEnWm80wZ41yFwulUzypeMM/KU7USbNly1G0K5JfxZERjN+RRQxCanU9HZlWESDvAnAhBDXi0uN48llT3I5/TL/vaUfDZa+Bbd2ht4/gqNzvnXL4t9XcedykUCvhBIvHuCZX3tx0gTP1niNL36vme/CjavZlG9WRyHEXxLTExmwfADRydFMafgszZa8AUGtoP88cHKzyj5lci5RKK+AxkzuOp2gbM3UmE/xdNyZ731rXoUXwp4lpicyaNUgTl05xcSQZ2m2dCRUD4W+c6wW5jdDAr2S8g1swZR7/oNflgVqzaaqy6F871vrKrwQ9ir2aiw9FvTnQFwUfqfuoumit7niVgv6zwcXT6PLAyTQKzX/Oh0JuDQAL4sFc+3peDofz3vPWlfhhbBHZ5PO0nNhPy5cjcHvzH38nDmbixYvHkgYyoIo2zn5kUCv5B7u0g+/6Eeooi241ZmCj+sBq16FF8LeHE84nnMBNC2R6qfv539ZP5CKM/0z3+R0pqdNNU9KoFdy3ZsH8li3QdSPfQzf7CxMtX/guQ5n5YKoEMCBuAMMWD4ACxbqnryX+ZbJpGOmb8ZbnNX+gG01T0qgC7o3D+Q/I95jxp2fc0tWNt/HfM6i3VOMLksIQ20/v51nVj6Du9mdH0Ke42f1LbHam0fS3+NPXSNvPVvyaLJfAAAWOklEQVRqniwy0JVS05RSF5VS+/+2zFcptUopdTT32ce6ZYryULV+BNO6TqdFhoW39kxk+uYPjS5JCKtZEBlNu3FrC5wFcePZjbyw+gWquVXj+7q9qPXLS6R53crj+n3OUTVvPVtrnizOGfp0oOs/lo0A1mit6wNrcl+LCsC9Vmu+eWgeXdI1nxyZzadrX6c8xyoIUR6uTTVd0CyIy/5cxitrX+EW71uYXiOCaotehaBWeL+wnH/3aFfgSGxbUayBRUqpYGCx1rpJ7uso4E6t9TmlVA1gvda6yP9NycAi+5GdcIqx/+vGHCcLD/m34r2u35b4juVC2Jp249YSXUDbt3+NSNK9fya8WjhfuTWiytoxUL8LPPq9of3MrT2wqJrW+hxA7nPADQp5Tim1Qym1IzY2toS7E+XN5F2Ht3qv4sVMFxbGbmPIol6kZtnOxR8hSqOgC5lm3w2kec+hfWB7Jplq5YR5k57Qa5ZNDBoqDqtfFNVaT9Zat9Bat/D397f27kQZUh4BvNB/NSMtXmxMiGLQL91JTE8Ebtz+KISty38hU+PkvwKXasswpzbni1QzLpsmwu1PQY8p4OhkWJ03q6SBfiG3qYXc54tlV5KwKS5e9Oq3kvGmQPalRDNg/oP8sG1voe2PQtiDv2ZBtOBcbSHOfusgsRXzTZmYd06H9q/CA5+Bg6moTdmUkgb6IuDJ3J+fBBaWTTnCJjm5EdHnV/7j0oCYtHi+3fcUaepcvlVk/hdhT7o3D2RU9/r4BP+Mk+8WXJM7sdJ8geDopXDPezkPpYwtsgSK021xNrAZaKCUOquUegYYB9yrlDoK3Jv7WlRkjk60eexnpnm3wlGl4l/nC0wup/OtYksDLIS4kZjkGP4XM4Js1z283uxFtroeo8aF33POytu/anR5JVZkoGut+2ita2itzVrrIK31VK11vNa6s9a6fu7zpfIoVhjMwUTj7tPocaEpPjod7zr/walK3vAEmxpgIURhdpzfQe/FvYlOiuab9uMYsHU26ux26PkdtHja6PJKRUaKipujFHW6fkL76LbckpmOc62ZePktwdWsbGqAhRAFmXN4DgNXDsTL2YsfW4+i/a/DITYKes+G0EeMLq/UpGOxuGk5AylGk7n0P9Tzms0S/99oViueuxvLdAHCNmVmZzJ221j+d+R/dAjswEf+HfD46Qlw8YIBSyHodqNLLBNyxyJRKjp6F3MWPM5H7ooaLr58HvEdt/ncdt16cps7YZT41HheW/8auy7u4pnGT/NywhVMv3+ac5ehXjPAo7rRJRZJ7lgkyoUKDKf3k+uYZvEnNSWO/r8+xrLjS/Ktc6Nh1kJY06H4Q/Re0puD8Qf5qM27DDmyJSfMw5+AAYvtIsxvhpyhi7KRlUHs0ld5/dxKIl1ceOK2x3i19Rs4OjgWOsw60NuVTSPuNqBYUdEU9A3QxXsfb296Gy9nLyaGD6XR0pFw+U/oOg5aPmtX3RKLe4YubeiibDg64d/ta6Zu+47xWz/khyM/cyjuAOPv+brQ7ozSzVGUhWvfAK/d6Dw6IYW31n+Mg+9awgPC+aRWN/x+HggmMzyxEILbG1yx9UiTiyhT5lbP8ma3WXx4JYO98fvpteAhAvwvFLiudHMUZWH8iqi8MMchDdegH3DwXYs5pS3fuTfBb+4z4BMMz62v0GEOEujCGuq05cEn1jIj0xvHlHgy/L7Aveq2fKvY2jzSwn5d+6bn4Hwet+CvMVU5Qvb5+3jv3EnMaz+AJj3g6RXgXdvYQsuBBLqwDq9AQgas4ieftrS6moJDwHwCas9DqUybnEda2K8a3s6YfTfgFvwlypRKlTOPMDdlGQ+YtsE9o6DnVLuZLbG0JNCF9Zhd8e4xlW+avcLAhCukum+ndevvmfFCsIS5KBNnk87ie+tUXKotIyu5ISEn72NJ1rcEqVi2tJkE7YfY1cXP0pJAF9alFKY7XmbwA//l80tXOXv5CI8u7MH0/dPJtmQbXZ0oB9aYallrzS9Hf6Hnop7EZfxJj8BXeS3Ri58cPuWKgxdb75nLHV17l0H19kW6LYryc+lP4uY/zfsZp1nn7kZT30aM7jiWel71jK5MWMk/e6BAzvWT0jS5xaXGMWrzKNafWU/L6i35IORpai5/G2J2QeijcP+n4OJZVodgE2RgkbA9vnXxe2oVXzR+no9iL3Mq7gCPLuzJf/f/V87WK6h8PVBylWaq5TWn19BzUU/+iP6DYbcP5TuP26n5Q0+4fDLnNnE9v6twYX4zJNBF+TI5ojoN477+y1iQVoX2yVf4dOenPLGkLycSTxhdnShjZTUGITkjmZG/j2TIuiFUc6vGnI6f8cSO/+Gw8k2odye8uAUady99wXZOAl0Yo3oofgM38Hn9/nwUG8+puIPXna3Lbe7sX2FjDW5mDML289vpuagnv574lYGhA5lV4z5undUXYnZDt6+gz0/gUa2sSrZr0oYujHd6K3ELnuMDx2TWuLvRtGoTOlUdzGdLE8u07VWUv9K0oadnpzNx10RmHJxBLY9ajAkfStgf38KRZVCnPXT/BnzqWPsQbEJx29Al0IVtSE9Gr3yb5Yd/4kM/PxKVibSLXci41IG/f5GU+V/sT0lm2hy/cSEzj3yJxXwBc0o7vqjdnA4HxkHmVej8LrR+HhwqTwODBLqwT8fWELfoJUa7ZLLW3RVLahCpMY9iycj5Sq2AP8fdb2yNwmpOJJxg6NoPOJq0HUtGVRzOd+H9jN952LSJy95N8Ok3Dfwr3whj6eUi7NOtnfF74Q+6xoXw8cU4PJ3O4l7vc5yrz0c5XpH5XyqohLQEPtz6IT0W9eBo4n7SLtxH+MnOrMyawgMOW/g08xEeSn23Uob5zZDZFoXtcfUh/cFJrJ3/HbOuTmGOjyM/eW/DyWsXYdUfISmjJR5OHkZXKcpAZnYmsw/PZtLeSaRkpvDobY+yZrEHb6lfuMccyRFLIAMzX2O/rodKzDS6XJsngS5sUk4b67O8vLwpvS7OYuGV9Xzt68vyiz9x3/xlDGr6PI81eAwnk1O+z8mdkeyD1pr1Z9bzyc5POHXlFO1qtmNo42e4NfInRjh+wlXtzNjMPkzPjiCdnL+xfDsrmrShC/sQGwWr3uHAyTV8FlCDrWYIrBLI4OaD6Vq3Kw7KwSqjEkXZi7oUxfjt49l6fiv1vOoxtPlgOpzZB799AhkpnKjzGP2P301MpnveZyr737FcLooqpU4CSUA2kFXUDiXQRamdWI9e8RZ/XDnGZ9VqEqWyCPEN4dXbX2XoD6lyZyQbFpcax1eRX/HLsV/wcPLgxWYv8miWI+Y1oyHhNNSPgC6jwb+BfNP6h/IM9BZa67jirC+BLsqEJRv2zMayZjRLuMKXATU5pzPISq5P+sV/YUmvmW916RljrLjUOOZEzWHGwRmkZ6XTJ6QPg6q2xGvtGDi7HaqFQsQHOSM+RYHkFnSi4nIwQfP+ODTqzoN/fEmXTV/wUxVnvvI8gWPdiWReaU5G3F1YMgIAaXs1ytHLR5lxcAaLTywm05JJ59qdGXLrowRvmQJLxkCV6jkjPcP65vxNRamV9gz9T+AyoIFvtdaTb7S+nKELq7gSA2s/IHHvbL728mOOpxsWBwtZKbegrrRjdJde9Ayv+HersQVaazbFbOKHAz+w+dxmXEwuPHTrQ/Sv143g3f+DrZNAmaDdYLhjMDhXMbpku1BeTS41tdYxSqkAYBXwstZ64z/WeQ54DqB27dq3nzp1qsT7E+KGzu0hdt4wHC5tZ4aHH3M8PUl2zKCaWzUea/AYPev3pKpr1VLvpiK175bVsaRlpbH4xGJmHpzJ8cTj+Lv60zekL49UDcc78kfY/WPOKM+wvnD3SPCsWfRGRZ5yHymqlHoPSNZaTyhsHTlDF1anNRxfA1v+Q9ax1Wyo4slPNYLZkpWAo4MjEcER9G7Qm2b+zVAluJNNSXvS2OL/BMqiV9C19vE5h+dwOf0yIb4hPN7ocbpanDFvmwJHV4LJCZo8Am3/D6o3sdbhVGhWD3SllDvgoLVOyv15FfC+1np5YZ+RQBflKjYKtn4Le2ZzggzmBDZgoSmDFEs6Ib4h9GnYh651u+LqWPw29nbj1t50Txpb7U5ZkmO55sjlI8w4OIMlJ5aQZcmiU61OPFH/MVqcP4ra9i3EHgb3AGj5DLR4GqoEWOswKoXyCPR6wC+5Lx2BH7XWY270GQl0YYjUy7DrB9g6mZSkaBYH1OYnbx+OZVzC08mTh299mMcaPEZtz6Lb2euOWEJB/2Ju1JOmNMFpTTdzLFprjicc58ut81l/di0WczRoM638Ini7eXeCDy2DndNzftfVm0KbF6FJD3B0Lo9DqfCs3stFa30CaFbSzwtRblx9oN0r0Ob/cD/8K722TOKxqC3sqOLD7KDqzDw0g+8Pfs9tPrfRMagjnYI6EeoXiqmAnhc1vV0LDOcb9aQpq5s8lLWijkVrzYH4A6w+tZo1p9dw8spJtFZkZ9Ym69L9hCT60PPEOmrv/BbQ0PABaPMC1G5bqW7MbEtkpKionKJ35fS42D+fC0qzNDiMja4uRF6NJltn4+3sTYfADnSs1ZE7at6Bp1PObc1K0nxSXmfoN9tOX/CxKAZ1cSDNeTdrTq/hfMp5TMpEq+qtiDwchPmiLxH6MA+Z/iDM4ThXtBuLHe+l70ujK83c5EaQ6XOFKI6k87B9KkTOhKQYEk2ObK7VlA1evvyWeo7EzCRMykR4tXA6BXWiY1BHdh83M2HlkVIGZ9m2oZfmYu3HKw5wMXM/nn6HcfI4SEp2Ik4OTtwReAf31O7MnU4BeB3fwIF1s2nskNNL7bClFrOyOzMvuyOpuMjALSuTQBfiZmgNMZFweEnOI/YQ2cC+miFs8K/DBksSR5PPAFDLoxadgjrRqnor6vvUp2aVmjioG89Ebe1eLsX9FpBpyeREwgkOxh/k0KVDHIw/SNSlKNKy03BzdKNTUCc617qLDhYTbkfXwOGlkHgaUOxRDfk1I5xVlts5pasXug9R9iTQhSiN+OMQtTQn3E9vATTnfGqzMagxG5xgW+Ix0rPTAXB1dOVW71vzHvV96lPfpz5VXaqWqGtkSRR8gTMLk8sFPu5bNS/Aoy5FkWHJAMDN0Y2Gvg1pVLURbfzDaJOSjPORlXB0Rc7FTUcXqHcXNLwfbuvKgqMZNtlbpzKQQBeirCTHwpHlOeF+fC1kp5Pq6sOR2uEc8/TjmKMjR7OTOJp0mktpl/I+5u3snS/kb/G+haouVXE3u+NmdsPV0bXIM/vCZGZnkpiRSEJaAgnpCbw4+zcupyegTFdR5kuYXKJxcDmPUjnh62H2IKRqCCG+ITTyDSHEVIU6yXE4XDiQc7PlU5sgKw1cvKHBv3JC/Ja7wck9335tsT99ZSCBLoQ1ZKTkhPq1M/fLf/71nltV4quFcNwniKNuHhxTWRxLj+NYwnGSM5ML3Jyro2tOwDu64W52z3t9LfTNDmauZFwhMT2RhPSEvOeUzJRCS9TZrmSnBuKQGUTfJi14vLYfQYkXcLiwH87vgwsH4NrnHRzBrwHU7ZgT4rXbgkmmeLI1EuhClIe0K3DxYE5Qnt+bG5gHIbc5BpMTOqAh5wNu45hHAInOLlx1MHFVKVLQpGDhqs7ialYaKVkpXM28SkpmCqlZqaRkppBhycDTyRNvZ2+8nL3wdvbO+9nL0R1vBzNemPDWiuMn49m8OxqPlPO0cI6mjVsMHiknQVtyanH2hOqh+R/+DaWvuB2QQBfCKNlZEH/0byG/P+f5anzhnzG7g4vX9Q+zK2QkQ1ri9Y+stMK351X7+vD2ri39w+2UTJ8rhFFMjhAQkvNo+ljOMq0h+SJcjSs4nNMSIS3hr5+Tz0NcFGRcBWePvwLeM7CA4PfO/9qjWs5gKlHpSKALUR6Uyglaj2rFWv26i48d5OKjKJoEuhA25p+DhKITUnlj/j4ACXVxQyXrMyWEsJrxK6Ly9fUGSM3MZvyKKIMqEvZCAl0IG2Ork3kJ2yeBLoSNKWzmRrk3qiiKBLqo1BZERtNu3FrqjlhCu3FrWRAZbXRJDItogKs5/9S9rmYTwyIaGFSRsBdyUVRUWrZ68fHavmWIvbhZEuii0rrRxUejw7N780DDaxD2R5pcRKUlFx9FRSOBLiotufgoKhoJdFFpycVHUdFIG7qotOTio6hoJNBFpSYXH0VFIk0uQghRQZQq0JVSXZVSUUqpY0qpEWVVlBBCiJtX4kBXSpmAr4F/AY2APkqpRmVVmBBCiJtTmjP0VsAxrfUJrXUG8BPwUNmUJYQQ4maVJtADgTN/e302d1k+SqnnlFI7lFI7YmNjS7E7IYQQN1KaXi4F3ZzwuhuUaq0nA5MBlFKxSqlTJdyfHxBXws/aGjkW21NRjgPkWGxVaY6lTnFWKk2gnwVq/e11EBBzow9orf1LujOl1I7i3CTVHsix2J6Kchwgx2KryuNYStPksh2or5Sqq5RyAnoDi8qmLCGEEDerxGfoWusspdRLwArABEzTWh8os8qEEELclFKNFNVaLwWWllEtRZlcTvspD3IstqeiHAfIsdgqqx+L0vq665hCCCHskAz9F0KICsKuAl0pNVoptVcptVsptVIpVdPomkpKKTVeKXU493h+UUp5G11TSSilHlVKHVBKWZRSdtkboaJMYaGUmqaUuqiU2m90LaWhlKqllFqnlDqU+9/WK0bXVFJKKRel1Dal1J7cYxll1f3ZU5OLUspTa30l9+fBQCOt9fMGl1UiSqkuwNrci8sfAWithxtc1k1TSoUAFuBbYKjWeofBJd2U3CksjgD3ktMVdzvQR2t90NDCSkAp1RFIBn7QWjcxup6SUkrVAGporXcppTyAnUB3O/2bKMBda52slDIDvwOvaK23WGN/dnWGfi3Mc7lTwEAme6G1Xqm1zsp9uYWcfvx2R2t9SGsdZXQdpVBhprDQWm8ELhldR2lprc9prXfl/pwEHKKAUej2QOdIzn1pzn1YLbfsKtABlFJjlFJngH7AO0bXU0aeBpYZXUQlVawpLIQxlFLBQHNgq7GVlJxSyqSU2g1cBFZpra12LDYX6Eqp1Uqp/QU8HgLQWr+lta4FzAJeMrbaGyvqWHLXeQvIIud4bFJxjsOOFWsKC1H+lFJVgHnAkH98O7crWutsrXUYOd/CWymlrNYcZnN3LNJa31PMVX8ElgDvWrGcUinqWJRSTwIPAJ21DV/MuIm/iT266SkshPXltjfPA2ZprecbXU9Z0FonKKXWA10Bq1y4trkz9BtRStX/28tuwGGjaiktpVRXYDjQTWt91eh6KjGZwsLG5F5InAoc0lp/anQ9paGU8r/Wg00p5QrcgxVzy956ucwDGpDTq+IU8LzWOtrYqkpGKXUMcAbicxdtscceO0qph4EvAX8gAdittY4wtqqbo5S6D/icv6awGGNwSSWilJoN3EnOrH4XgHe11lMNLaoElFLtgd+AfeT8Wwd4M3dkul1RSjUFvifnvy0H4Get9ftW2589BboQQojC2VWTixBCiMJJoAshRAUhgS6EEBWEBLoQQlQQEuhCCFFBSKALIUQFIYEuhBAVhAS6EEJUEP8PsXE+OapvEM0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "beta,*_ = linalg.lstsq(x_mat, y) \n", "print(beta)\n", "plt.plot(x,y, 'o', label='data')\n", "plt.plot(x,2 * x ** 2 - 3 * x + 4, '-', label='exact solution')\n", "plt.plot(x,x_mat.dot(beta), label='least squares')\n", "plt.legend(loc='upper right')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Example Error Analysis \n", "---" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.2894065474411858\n" ] } ], "source": [ "yhat = x_mat @ beta\n", "resids = yhat - y\n", "SSR = np.sum(resids**2)\n", "se2_epsilon = SSR / (len(x) - len(beta))\n", "print(se2_epsilon)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0.01088978 0. -0.03539179]\n", " [ 0. 0.02817731 0. ]\n", " [-0.03539179 0. 0.20659959]]\n" ] } ], "source": [ "se2_beta = se2_epsilon * linalg.inv(x_mat.transpose() @ x_mat)\n", "print(se2_beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Now that we have the standard error matrix, we can create confidence intervals for the $\\beta$ values." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "beta_0 is 1.98 +/- 0.22 with 95% confidence\n", "beta_1 is -2.87 +/- 0.35 with 95% confidence\n", "beta_2 is 4.23 +/- 0.94 with 95% confidence\n" ] } ], "source": [ "for i in range(len(beta)):\n", " #get our T-value for the confidence interval\n", " T = scipy.stats.t.ppf(0.975, len(x) - len(beta)) \n", " # Get the width of the confidence interval using our previously computed standard error\n", " cwidth = T * np.sqrt(se2_beta[i,i]) \n", " # print the result, using 2 - i to match our numbering above\n", " print(f'beta_{i} is {beta[i]:.2f} +/- {cwidth:.2f} with 95% confidence') " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Goodness of Fit\n", "====\n", "\n", "The same equations above apply here. For the example of the polynomial:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "(0.967408962984229, 0.9835695008408043)" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "TSS = np.sum( (np.mean(y) - y)**2)\n", "R2 = 1 - SSR / TSS\n", "R2, np.sqrt(R2)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The Statistics of Regression\n", "===\n", "\n", "You can justify a regression and its coefficients by:\n", "\n", "1. A Spearmann correlation test\n", "2. A Shapiro-Wilks normality test on residuals\n", "3. A Goodness of fit" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "You can justify a particular **model** by:\n", "\n", "1. Plotting\n", "2. Hypothesis tests on individual coefficients\n", "\n", "You CANNOT use a goodness of fit to compare models, since each time you add a new parameter you get a better fit" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Conducting a good regression\n", "===\n", "\n", "You should do the following when you want to do a good job:\n", "\n", "1. Justify correlation\n", "2. Compute your model coefficients, their standard errors, confidence intervals, and p-values for existing\n", "3. Plot your regressed model (use $y$ vs $\\hat{y}$ for N-dimensions)\n", "4. Histogram and compute Shapiro-Wilks normality test on residuals" ] } ], "metadata": { "celltoolbar": "Slideshow", "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.1" } }, "nbformat": 4, "nbformat_minor": 1 }