week05a-corr-regress-least-squares.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Correlation and linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation\n",
    "\n",
    "Correlation $\\neq$ Causation\n",
    "\n",
    "* If two variables are correlated, variability in one variable *might* cause variability in the other. Or there may be another underlying variable that causes variations in both.\n",
    "* Statistically significant correlations may be *spurious* - see examples at https://tylervigen.com/spurious-correlations\n",
    "\n",
    "![images/correlation-causation-cage.png](images/correlation-causation-cage.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The caclulation of correlation is based on covariance.\n",
    "\n",
    "The __variance__ $s_x^2$ has to be positive.\n",
    "\n",
    "$$ s_x^2 = \\frac{1}{N-1}\\sum_{i=1}^N (x_i - \\bar{x})(x_i - \\bar{x})$$\n",
    "\n",
    "The __covariance__ $s_{xy}^2$ can be positive or negative.  If $x_i$ tends to be above the mean $\\bar{x}$ when $y_i$ tends to be above the mean, then the covariance is positive\n",
    "\n",
    "$$ s_{xy}^2 = \\frac{1}{N-1}\\sum_{i=1}^N (x_i - \\bar{x})(y_i - \\bar{y})$$\n",
    "\n",
    "The _Pearson correlation coefficient_ is the normalized covariance between two variables.\n",
    "\n",
    "$$ r_{xy} = \\frac{s_{xy}^2}{s_x s_y}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Varying degrees of correlation\n",
    "\n",
    "There may be a pattern in the data, but a linear model may not be appropriate.\n",
    "\n",
    "![images/Correlation_examples.png](images/Correlation_examples.png)\n",
    "\n",
    "[source](http://wikipedia.org/wiki/Correlation_and_dependence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parametric statistic (like Pearson's correlation are sensitive to outliers. For example, two radnomly generated sets of numbers might have a very low correlation coefficent, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'r = 0.035')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAARdUlEQVR4nO3dcWxd5X3G8eeZMeVqULkTXtsY0mQScscIYMlCVP6nzdAcYFpNtEqgrlO3avlnXalWuUvUVVs3ackUqaoqMW1Zi9qpFaxTg2GE1aMKWUTFoA6GhCz1xBiI3LDFtHUhnQuJ+e0PX1PbsZ17c997z3l9vx/pqvY5h/e8Ob15cs7vfc85jggBAPL1C0V3AADQHIIcADJHkANA5ghyAMgcQQ4AmSPIASBzBDkAZI4gR0ezfaPtI7b/r/a/N66x7Tts32v7Ndv/Y/uPF6270vb3bP/Q9oztJ2wPLVr/cdtzts8s+nywxX88dAiCHFmxfUnCti6V9KCkb0h6l6SvS3qwtnwlfy7pGknvk/QhSZ+1va227oyk35fUW2vrryX987L+PhERly/6HEr1Z0FnI8hRerZftP0nto9K+mnCMP+gpEskfSki3oiIL0uypK2rbP+7kv4yIn4cESck/b2kj0tSRPwsIqYi4q1aG3OaD/RfStRXYFUEOXJxl6TbJfVExLnlK20frZU0Vvr8zSpt/pqko7H0ORVHa8uXt/8uSRskPbto8bPLt639Y/MzSQ9J+kpEnF60esD2q7b/0/bnU15doLPxRUIuvhwRL6+2MiKuv4g2L5f0k2XLfiLpilW2XVi/6rYRcb3tyyTdIWlxieawpOskvaT58P9HSeck7b6IfgNLcEaOXKwa4k04I+mdy5a9U9Lrq2y7sH7NbWtllvsk7bR9Q23ZCxHx3xHxVkQck/QXkn672T8AIBHkyMeaj+m0fXzZjJDFn79d5T87Lul621607Pra8qU7j/ixpFck3bBo8Q0rbbtIt6RfWWVdaL6WDjSN0grWhYg4r65dh0OaH5T8VC3s/6C2/OAq2/+DpD+1PSHp3bXtf0+SbN+s+b9PT0nqkvSp2jZP1tbfKunpiPhf2++X9HlJ/3QRfQbOwxk5OlZEvClpRPOzUWY0P31wpLZctj9qe/EZ959J+i/N17n/TdLeiPhObd07JN0j6YeSqpJuk3R7RJyqrf91SUdt/1TSI5L2S/qrFv7x0EHMiyUAIG+ckQNA5ghyAMgcQQ4AmSPIASBzhUw/vPLKK2PTpk1F7BoAsnXkyJFXI6J3+fJCgnzTpk2amJgoYtcAkC3bL620nNIKAGSOIAeAzBHkAJA5ghwAMkeQA0DmePohgIaNTVa1d3xKp2ZmtaGnotHhfo0M9BXdrY5FkANoyNhkVbv2H9Ps2TlJUnVmVrv2H5MkwrwglFYANGTv+NTbIb5g9uyc9o5PFdQjJDkjt/2i5l95NSfpXEQMpmgXQPmcmpltaDlaL2Vp5UMR8WrC9gCU0IaeiqorhPaGnkoBvYFEaQVAg0aH+1Xp7lqyrNLdpdHhfo1NVjW056A27zygoT0HNTZZLaiXnSXVGXlI+lfbIenvImLf8g1s75C0Q5I2btyYaLcA2m1hQHP5rBVJDIIWJMmr3mxviIhTtn9Z0qOS/igiDq+2/eDgYPDQLGB9GdpzcMWSS19PRd/bubWAHq0/to+sNAaZpLSy8ILZiDgt6QFJN6VoF0A+GAQtTtNBbvsXbV+x8LOk35D0XLPtAsjLaoOdDIK2Xooz8ndLetz2s5KeknQgIr6ToF0AGVlrEBSt1fRgZ0S8IOmGBH0BkLHVBkEZ6Gw9btEHkMzIQB/BXQDmkQNA5ghyAMgcQQ4AmSPIASBzBDkAZI4gB4DMMf0QDeEVX0D5EOSoG6/4AsqJ0grqxiu+gHIiyFE3nm4HlBNBjrrxdDugnAhy1I2n2wHlxGAn6sbT7YByIsjREJ5uB5QPpRUAyBxBDgCZI8gBIHMEOQBkjiAHgMwR5ACQuWRBbrvL9qTth1O1CQC4sJRn5HdLOpGwPQBAHZIEue2rJN0u6Ssp2gMA1C/VGfmXJH1W0lurbWB7h+0J2xPT09OJdgsAaDrIbf+mpNMRcWSt7SJiX0QMRsRgb29vs7sFANSkOCMfkvRbtl+UdL+krba/kaBdAEAdmg7yiNgVEVdFxCZJd0o6GBG/03TPAAB1YR45AGQu6WNsI+KQpEMp2wQArI0zcgDIHEEOAJkjyAEgcwQ5AGSOIAeAzBHkAJA5ghwAMkeQA0DmCHIAyBxBDgCZI8gBIHMEOQBkjiAHgMwR5ACQOYIcADJHkANA5ghyAMgcQQ4AmSPIASBzBDkAZI4gB4DMNR3kti+z/ZTtZ20ft/2FFB0DANTnkgRtvCFpa0Scsd0t6XHb/xIR/56gbQDABTQd5BERks7Ufu2ufaLZdgEA9UlSI7fdZfsZSaclPRoRT6ZoFwBwYUmCPCLmIuJGSVdJusn2dcu3sb3D9oTtienp6RS7BQAo8ayViJiRdEjSthXW7YuIwYgY7O3tTblbAOhoKWat9Nruqf1ckXSLpB802y4AoD4pZq28V9LXbXdp/h+Gb0XEwwnaBTra2GRVe8endGpmVht6Khod7tfIQF/R3UIJpZi1clTSQIK+AKgZm6xq1/5jmj07J0mqzsxq1/5jkkSY4zzc2QmU0N7xqbdDfMHs2TntHZ8qqEcosxSlFQCJnZqZbWh5p6HstBRBDpTQhp6KqiuE9oaeSgG9qV87Apay0/korQAlNDrcr0p315Jlle4ujQ73F9SjC1sI2OrMrEI/D9ixyWrS/VB2Oh9BDpTQyECfdm/for6eiiypr6ei3du3lPqMs10BS9npfJRWgJIaGegrdXAv166AzbXs1EqckQNIYrUgTR2wOZadWo0gB5BEuwI2x7JTq1FaAZDEQpC2Y1pgbmWnViPIASRDwBaD0goAZI4gB4DMEeQAkDmCHAAyR5ADQOYIcgDIHNMPgTbjEaxIjSAH2ohHsKIVKK0AbcQjWNEKBDnQRjyCFa2QTWmFuiLWAx7BilbI4oy8XW8eAVqNR7CiFZoOcttX237M9gnbx23fnaJji1FXxHrBI1jRCilKK+ckfSYinrZ9haQjth+NiP9I0LYk6opYX3hCIFJrOsgj4hVJr9R+ft32CUl9kpIFOXVFAGVX5Dhe0hq57U2SBiQ9ucK6HbYnbE9MT0831C51xXIZm6xqaM9Bbd55QEN7DjJWgY5X9DhesiC3fbmkb0v6dES8tnx9ROyLiMGIGOzt7W2obeqK5VH0FxYoo6LH8ZJMP7TdrfkQ/2ZE7E/R5nLUFcthrS8s//+gUxU9jpdi1oolfVXSiYj4YvNdQpkV/YUFymi18bp2jeOlKK0MSfqYpK22n6l9bkvQLkqo6C8sUEZFj+OlmLXyuCQn6AsyMDrcv+ShTxIDz8BCWbGoWSvZ3KKPcij6CwuUVZHjeAQ5GsbAM1AuWTxrBQCwOs7I0TCeRAmUC0GOhvCGG6B8KK2gIUXfwQbgfAQ5GsINQUD5EORoCDcEAeVDkKMhRd/BBuB8DHaiIdwQBJQPQY6GcUMQUC6UVgAgc5yRA+sUN251DoIcWIe4cauzUFoB1iFu3OosBDmwDnHjVmchyIF1iBu3OktHBvnYZFVDew5q884DGtpzkDfAY93hxq3O0nGDnQwCoRNw41Zn6bggX2sQqJO/5ExVW3+4catzdFyQl30QqIhA5SoFyFuSGrnte22ftv1civZaqcyDQAuBWp2ZVejngdrqGj5T1YC8pRrs/JqkbYnaaqkyDwIVFahlv0oBsLYkQR4RhyX9KEVbrTYy0Kfd27eor6ciS+rrqWj39i2lKCEUFahlvkoBcGFtq5Hb3iFphyRt3LixXbtdUVkHgTb0VFRdIbRbHaijw/1LauRSea5SAFxY2+aRR8S+iBiMiMHe3t527TYrRZV9ynyVAuDCOm7WSpkVOfe3rFcpAC6MIC8ZAhVAo1JNP7xP0hOS+m2ftP2JFO0CAC4syRl5RNyVoh0AQOM68qFZALCeEOQAkDmCHAAyR5ADQOYIcgDIHEEOAJkjyAEgcwQ5AGSOW/RRCF4tB6RDkKPteLUckBalFbQdr5YD0iLI0Xa8Wg5IiyBH2/FqOSAtghxtV+YXYKM1xiarGtpzUJt3HtDQnoMam6wW3aV1hcFOtF2Rb0JC+zG43XoEOQrBm5A6x1qD23wH0qC0AqClGNxuPYIcQEsxuN16BDmAlmJwu/WokQNoKQa3W48gB9ByDG63FqUVAMhckiC3vc32lO3nbe9M0SYAoD5NB7ntLkn3SLpV0rWS7rJ9bbPtAgDqk+KM/CZJz0fECxHxpqT7JX04QbsAgDqkCPI+SS8v+v1kbdkStnfYnrA9MT09nWC3AAApTZB7hWVx3oKIfRExGBGDvb29CXYLAJDSBPlJSVcv+v0qSacStAsAqEOKIP++pGtsb7Z9qaQ7JT2UoF0AQB2aviEoIs7Z/qSkcUldku6NiONN9wwAUJckd3ZGxCOSHknRFgCgMdzZCQCZI8gBIHMEOQBkjiAHgMwR5ACQOYIcADJHkANA5ghyAMgcQQ4AmeOdnQDQBmOT1Za9gJogB4AWG5usatf+Y5o9OydJqs7Matf+Y5KUJMwprQBAi+0dn3o7xBfMnp3T3vGpJO0T5ADQYqdmZhta3iiCHABabENPpaHljSLIAaDFRof7VenuWrKs0t2l0eH+JO0z2AkALbYwoMmsFQDI2MhAX7LgXo7SCgBkjiAHgMwR5ACQOYIcADLXVJDb/ojt47bfsj2YqlMAgPo1e0b+nKTtkg4n6AsA4CI0Nf0wIk5Iku00vQEANKxtNXLbO2xP2J6Ynp5u124BYN274Bm57e9Kes8Kqz4XEQ/Wu6OI2CdpnyQNDg5G3T0EAKzpgkEeEbe0oyMAgIvD9EMAyFyz0w/vsH1S0gckHbA9nqZbAIB6NTtr5QFJDyTqCwDgIlBaAYDMEeQAkDmeR46sjU1WW/awfiAXBDmyNTZZ1a79x95+O3l1Zla79h+TJMIcHYUgR7b2jk+9HeILZs/Oae/4FEHeZlwZFYsgR7ZOzcw2tBytwZVR8RjsRLY29FQaWo7WWOvKCO1BkCNbo8P9qnR3LVlW6e7S6HB/QT3qTFwZFY8gR7ZGBvq0e/sW9fVUZEl9PRXt3r6Fy/k248qoeNTIkbWRgT6Cu2Cjw/1LauQSV0btRpADaMrCP6TMWikOQQ6gaVwZFYsaOQBkjiAHgMwR5ACQOYIcADJHkANA5hzR/hfa256W9FLbd1yfKyW9WnQnSopjszaOz+o4Nmur9/i8LyJ6ly8sJMjLzPZERAwW3Y8y4tisjeOzOo7N2po9PpRWACBzBDkAZI4gP9++ojtQYhybtXF8VsexWVtTx4caOQBkjjNyAMgcQQ4AmSPIl7H9EdvHbb9lm+lSNba32Z6y/bztnUX3p0xs32v7tO3niu5L2di+2vZjtk/U/l7dXXSfysT2Zbafsv1s7fh84WLaIcjP95yk7ZIOF92RsrDdJekeSbdKulbSXbavLbZXpfI1SduK7kRJnZP0mYj4VUk3S/pDvjtLvCFpa0TcIOlGSdts39xoIwT5MhFxIiJ4a+xSN0l6PiJeiIg3Jd0v6cMF96k0IuKwpB8V3Y8yiohXIuLp2s+vSzohiQeX18S8M7Vfu2ufhmegEOSoR5+klxf9flL8ZUSDbG+SNCDpyWJ7Ui62u2w/I+m0pEcjouHj05FvCLL9XUnvWWHV5yLiwXb3JwNeYRnzVlE325dL+rakT0fEa0X3p0wiYk7SjbZ7JD1g+7qIaGi8pSODPCJuKboPmTkp6epFv18l6VRBfUFmbHdrPsS/GRH7i+5PWUXEjO1Dmh9vaSjIKa2gHt+XdI3tzbYvlXSnpIcK7hMyYNuSvirpRER8sej+lI3t3tqZuGxXJN0i6QeNtkOQL2P7DtsnJX1A0gHb40X3qWgRcU7SJyWNa36w6lsRcbzYXpWH7fskPSGp3/ZJ258ouk8lMiTpY5K22n6m9rmt6E6VyHslPWb7qOZPmB6NiIcbbYRb9AEgc5yRA0DmCHIAyBxBDgCZI8gBIHMEOQBkjiAHgMwR5ACQuf8HBvDWfkrG08AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "np.random.seed(437)\n",
    "x = np.random.normal(1,1,20)\n",
    "y = np.random.normal(2,2,20)\n",
    "\n",
    "r = stats.pearsonr(x,y)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x,y,'o')\n",
    "plt.title('r = '+str(round(r[0],3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, adding an outlier to this data set greatly inflates the correlation coefficient (to nearly 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'r = 0.96')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAASjElEQVR4nO3df6zdd33f8eerjhGG0Nk0N1Fyk9RZlVn8CNjTVZYt00ShYJeyxe3GaMZYpKG5lUCFinlymCqaiq6RssI0qWNzR4YlMmgk0iQrqGnkgGirDrjBaX7MTUNZSOx48W2DSxAeOOa9P+73ZjfX9/iee+85997POc+HdHTO93O+53zf/ih++Zv398dJVSFJas+PrHcBkqSVMcAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwCUiyM8mDSb7XPe88z7qTSe5J8lySY0l+ccH7m5J8NMkzSZ5PciTJ1uH/KTRuDHA1KckFA/yulwH3AJ8GtgGHgHu68cV8GvjfwCXAzwD/LslPznv/FuDvAX8X+FHgPcD/HVS90px4JaZakeRJ4BPAu4EdwCur6oUBfO/bgP8GXF7dX4gkTwH7qur3F6x7IfA8cHFVzXRjB4EtVfWeJNuAp4E3VtVfrLY26XzcA1drbmR2r3frYuGd5OEkp3o8/lOP73wd8HC9dG/m4W78nE0seJ57/fru9TXAC8A/SfJ/kvx5kvct488n9W1g/xsqrZH/WFVP93qzqt6wgu+8EPjrBWN/Dbxqke9/PskfA7+SZD/wWuAfAzPdKpcDfwP4W8BVwNXA4SR/XlX3r6A2qSf3wNWanuG9Ct9ltlc9348y2ypZzLuZDeenmW3p3AEc69473T3/WlWdrqqHgc8Cbx9oxRIGuNpz3oM2SR5L8t0ej//c42OPAW9IMr8t8oZu/NwCqr5VVe+oqomq+jvAjwFf7d5+uJ86pUGwhaKRUlWL9a2X8iXgLPBLXcj/q278gcVWTvIaZve4vw/8U+BtwGu67f9Fkj8E/m2SXwL+JvAuZnv30kC5B66xV1U/APYC/wI4BfxLYG83TpJ3J5m/N74b+CbwbeAXgT1zZ6R0bgR+HPgr4PPAr1TV4aH/QTR2PI1QkhrlHrgkNcoAl6RGGeCS1CgDXJIataanEV500UW1ffv2tdykJDXvwQcf/Muqmlg4vqYBvn37dqanp9dyk5LUvCTfWmzcFookNcoAl6RGGeCS1CgDXJIaZYBLUqO8G6EkDdHdR45z232P88yp01y2dQv7d+9g767JgXy3AS5JQ3L3kePcfNcjnD5zFoDjp05z812PAAwkxG2hSNKQ3Hbf4y+G95zTZ85y232PD+T7DXBJGpJnTp1e1vhyGeCSNCSXbd2yrPHlMsAlaUj2797Bls2bXjK2ZfMm9u/eMZDv9yCmJA3J3IFKz0KRpAbt3TU5sMBeyBaKJDXKAJekRhngktQoA1ySGrVkgCd5eZKvJvnTJI8luaUbf3WS+5M80T1vG365kqQ5/eyBfx94c1W9EdgJ7ElyHXAAOFxVVwOHu2VJ0hpZMsBr1ne7xc3do4AbgEPd+CFg71AqlCQtqq8eeJJNSR4CTgL3V9VXgEuq6gRA93xxj8/uSzKdZHpmZmZQdUvS2OsrwKvqbFXtBC4Hrk3y+n43UFUHq2qqqqYmJiZWWqckaYFlnYVSVaeALwF7gGeTXArQPZ8ceHWSpJ76OQtlIsnW7vUW4KeAPwPuBW7qVrsJuGdYRUqSztXPvVAuBQ4l2cRs4N9ZVb+X5E+AO5O8F3gKeOcQ65QkLbBkgFfVw8CuRcb/CnjLMIqSJC3NKzElqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSoJQM8yRVJvpjkaJLHknygG//VJMeTPNQ93j78ciVJc5b8VXrgBeBDVfX1JK8CHkxyf/fex6vq3w+vPElSL0sGeFWdAE50r59PchSYHHZhkqTzW1YPPMl2YBfwlW7o/UkeTnJ7km0Drk2SdB59B3iSC4HPAR+squ8AnwB+AtjJ7B76b/b43L4k00mmZ2ZmBlCyJAn6DPAkm5kN7zuq6i6Aqnq2qs5W1Q+B3wauXeyzVXWwqqaqampiYmJQdUvS2OvnLJQAnwSOVtXH5o1fOm+1nwUeHXx5kqRe+jkL5XrgPcAjSR7qxj4M3JhkJ1DAk8AvDKVCSdKi+jkL5Y+ALPLWFwZfjiSpX16JKUmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSoJQM8yRVJvpjkaJLHknygG391kvuTPNE9bxt+uZKkOf3sgb8AfKiqXgNcB7wvyWuBA8DhqroaONwtS5LWyJIBXlUnqurr3evngaPAJHADcKhb7RCwd1hFSpLOtaweeJLtwC7gK8AlVXUCZkMeuLjHZ/YlmU4yPTMzs7pqJUkv6jvAk1wIfA74YFV9p9/PVdXBqpqqqqmJiYmV1ChJWkRfAZ5kM7PhfUdV3dUNP5vk0u79S4GTwylRkrSYfs5CCfBJ4GhVfWzeW/cCN3WvbwLuGXx5kqReLuhjneuB9wCPJHmoG/swcCtwZ5L3Ak8B7xxOiZKkxSwZ4FX1R0B6vP2WwZYjSeqXV2JKUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNWrJAE9ye5KTSR6dN/arSY4neah7vH24ZUqSFupnD/xTwJ5Fxj9eVTu7xxcGW5YkaSlLBnhVfRl4bg1qkSQtw2p64O9P8nDXYtnWa6Uk+5JMJ5memZlZxeYkSfOtNMA/AfwEsBM4AfxmrxWr6mBVTVXV1MTExAo3J0laaEUBXlXPVtXZqvoh8NvAtYMtS5K0lBUFeJJL5y3+LPBor3UlScNxwVIrJPkM8CbgoiTHgI8Ab0qyEyjgSeAXhlijJGkRSwZ4Vd24yPAnh1CLJGkZvBJTkhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNWjLAk9ye5GSSR+eNvTrJ/Ume6J63DbdMSdJC/eyBfwrYs2DsAHC4qq4GDnfLkqQ1tGSAV9WXgecWDN8AHOpeHwL2DrguSdISVtoDv6SqTgB0zxf3WjHJviTTSaZnZmZWuDlJ0kJDP4hZVQeraqqqpiYmJoa9OUkaGysN8GeTXArQPZ8cXEmSpH6sNMDvBW7qXt8E3DOYciRJ/ernNMLPAH8C7EhyLMl7gVuBtyZ5AnhrtyxJWkMXLLVCVd3Y4623DLgWSdIyeCWmJDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIateSPGo+iu48c57b7HueZU6e5bOsW9u/ewd5dk+tdliQty6oCPMmTwPPAWeCFqpoaRFHDdPeR49x81yOcPnMWgOOnTvPLv/MQ0996jo/uvWadq5Ok/g2ihfKTVbWzhfAGuO2+x18M7zkF3PE/n+LuI8fXpyhJWoGx64E/c+r0ouPFbLhLUitW2wMv4A+SFPBfqurgwhWS7AP2AVx55ZWr3NziltPTvmzrFo73CPFe4S5JG9Fq98Cvr6q/Dfw08L4k/2DhClV1sKqmqmpqYmJilZs711xP+/ip0xSzPe2b73qkZztk/+4dpMd3XbZ1y8Drk6RhWVWAV9Uz3fNJ4HeBawdR1HIs1tM+feZsz3bI3l2TvPu6K88J8S2bN7F/944hVSlJg7fiAE/yyiSvmnsNvA14dFCF9atX2+N87ZCP7r2Gj79rJ5NbtxBgcusWfuPnrvFUQklNWU0P/BLgd5PMfc9/r6rfH0hVy9Crp71UO2TvrkkDW1LTVrwHXlXfrKo3do/XVdWvD7Kwfu3fvYMtmze9ZMx2iKRx0PyVmHN70V5ZKWncNB/gYDtE0ngauwt5JGlUGOCS1KiRaKH0wzsQSho1qao129jU1FRNT0+v2fbmLLwDIcDmTeGVL7uAU6fPsCnhbBWTBrukDSjJg4vdMHAsWiiLXa155mxx6vQZAM52/4gtdRm+JG0kYxHgvW5etZjzXYYvSRvJhu+Br7Z3ffeR44TZ2yb2y7sSSmrBht4D7/dOg3cfOc71tz7AVQc+z/W3PvCS92+77/FlhTfMhv3C75GkjWZDB3g/dxpcKuRXujdtP1zSRrehA7yfOw0uFfK9bmq17RWbmeze25TF7xBuP1zSRrahe+D93Gmw1wHKufH9u3eccwrhls2b+Mg/fN1LeulXHfj8oq0W++GSNqoNvQfez50Ge+09z43v3TXJb/zcNS/Z257bs57fHum1p+6v9EjaqDZ0gPcK31v+x2PsvOUPuOrA5188h3uh+eN7d02++I9Br3O+vS2tpNZs6ACHxcP32987w6nTZ857dsnkgj3npXrl8/+x8Fd6JLVgQ/fA5ywWvucTOGfPuZ8Dot6WVlJLNvweOCz/QGLBOUFsj1vSqGkiwFcSsgsvxLHHLWnUNBHgi4XvUhYepFzY4966ZTMv3/wj/PLvPORVl5KatKoAT7InyeNJvpHkwKCKWmhh+G57xWa2btn8YhBve8XmRT+38EKcvbsm+eMDb+bj79rJ91/4Id/+3pnzXqIvSRvZig9iJtkE/BbwVuAY8LUk91bV/xpUcfMtdYBxORfinO+MFA9iSmrFavbArwW+UVXfrKofAJ8FbhhMWcu3nIOU/ZyRIkkb3WoCfBJ4et7ysW5sXSznIKVnpEgaBasJ8MWuYT+ni5FkX5LpJNMzMzOr2Nz5LedCnPOF/fluTStJG8lqLuQ5Blwxb/ly4JmFK1XVQeAgzP4m5iq2t6R+L8SZW2fhD0UAL7nx1dzBzfmfkaSNYjUB/jXg6iRXAceBnwf+2UCqWgOLhf31tz7gwU1JzVhxgFfVC0neD9wHbAJur6rHBlbZOvDgpqSWrOpeKFX1BeALA6pl3fVz/3FJ2iiauJnVcvTzI8i91un14w9ebi9pIxqpAJ/7fczzHYTsZ52l/gGQpI1gpAK8nyssl1rHW8pKakUTN7PqVz8HIT1QKWlUjFSA93OFpVdhShoVIxXg/VxO733BJY2KkeqB93MQ0gOVkkZFqsevug/D1NRUTU9Pr9n2JGkUJHmwqqYWjo9UC0WSxokBLkmNMsAlqVEGuCQ1ygCXpEat6VkoSWaAbw1xExcBfznE72+F8+AczHEeZrU+Dz9eVRMLB9c0wIctyfRip9qMG+fBOZjjPMwa1XmwhSJJjTLAJalRoxbgB9e7gA3CeXAO5jgPs0ZyHkaqBy5J42TU9sAlaWwY4JLUqJEI8CR7kjye5BtJDqx3PWslye1JTiZ5dN7Yq5Pcn+SJ7nnbeta4FpJckeSLSY4meSzJB7rxsZmLJC9P8tUkf9rNwS3d+NjMwXxJNiU5kuT3uuWRnIfmAzzJJuC3gJ8GXgvcmOS161vVmvkUsGfB2AHgcFVdDRzulkfdC8CHquo1wHXA+7r/BsZpLr4PvLmq3gjsBPYkuY7xmoP5PgAcnbc8kvPQfIAD1wLfqKpvVtUPgM8CN6xzTWuiqr4MPLdg+AbgUPf6ELB3TYtaB1V1oqq+3r1+ntm/uJOM0VzUrO92i5u7RzFGczAnyeXAzwD/dd7wSM7DKAT4JPD0vOVj3di4uqSqTsBssAEXr3M9ayrJdmAX8BXGbC66tsFDwEng/qoauzno/Afg3wA/nDc2kvMwCgGeRcY8N3IMJbkQ+Bzwwar6znrXs9aq6mxV7QQuB65N8vr1rmmtJXkHcLKqHlzvWtbCKAT4MeCKecuXA8+sUy0bwbNJLgXonk+ucz1rIslmZsP7jqq6qxsey7moqlPAl5g9PjJuc3A98I+SPMlsO/XNST7NiM7DKAT414Crk1yV5GXAzwP3rnNN6+le4Kbu9U3APetYy5pIEuCTwNGq+ti8t8ZmLpJMJNnavd4C/BTwZ4zRHABU1c1VdXlVbWc2Cx6oqn/OiM7DSFyJmeTtzPa9NgG3V9Wvr3NJayLJZ4A3MXurzGeBjwB3A3cCVwJPAe+sqoUHOkdKkr8P/CHwCP+/7/lhZvvgYzEXSd7A7MG5TczumN1ZVb+W5McYkzlYKMmbgH9dVe8Y1XkYiQCXpHE0Ci0USRpLBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElq1P8DwxyWU0ui5hwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xo = np.append(x,45)\n",
    "yo = np.append(y,32)\n",
    "\n",
    "ro = stats.pearsonr(xo,yo)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(xo,yo,'o')\n",
    "plt.title('r = '+str(round(ro[0],3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing significance of linear correlation\n",
    "\n",
    "The following test statistic, will follow a t distribution for randomly generated sets of data.\n",
    "\n",
    "$$ t = |r|\\frac{\\sqrt{N-2}}{\\sqrt{1-r^2}}$$ \n",
    "\n",
    "To test for significant at confidence level $1-\\alpha$, calculate a t statistic and compare with a critical t value $t_{\\nu,1-\\alpha/2}$, where $\\nu = N-2$ is the degrees of freedom. The degrees of freedom is $N-2$ because you can always draw a line through two data points. Three data points are needed to obatin one independent piece of information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression (Type I) ###\n",
    "### Least-Squares Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Type I regression minimizes the variance along one axis__\n",
    "\n",
    "- x is the independent variable\n",
    "- y is the dependent variable\n",
    "- Goal: minimize sum of squared errors, SSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 7)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAD8CAYAAAC8TPVwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3DU9b3/8ecbCHKNUoxXVCRqExCMGFHrOVSsVdEeaj04R3tg7K/2RMULUFoq1Wp1WqunUmm9M2DlBxZHEZRB8R6lMCIkmnIRDj9T9CBewMELoEIgn98f72Q2C0g2sLvf/e6+HjMZsp/9kn3X6osPn6uFEBARkXhoF3UBIiKSOoW2iEiMKLRFRGJEoS0iEiMKbRGRGFFoi4jESKuhbWbfNrO6Fl9fmNmYbBQnIiLJrC3rtM2sPbAeOC2E8F7GqhIRkT1q6/DI94B6BbaISDQ6tPH5S4GZe3rDzKqAKoCuXbueUlZWtp+liYjEQGMjrF8PGzYkt3fqBMccA926pfRjamtrPwkhlLT2XMrDI2bWEfgA6BdC+Hhvz1ZWVoaampqUfq6ISGzNnQvXXgvr1iXaOnWC3/wGfvEL6Ngx5R9lZrUhhMrWnmtLT3so8GZrgS0ikvfWr4frroM5c5LbzzkHHngAjjsuYx/dljHty/iGoRERkYKwcyfccw+UlycHdkkJzJgBL7yQ0cCGFHvaZtYF+D5wZUarERHJVXV1UFUFS5cmt19xBfz3f8O3vpWVMlIK7RDCl0DPDNciIpJ7tm6FW26BSZO8p92srAweeggGD85qOW1dPSIiUjieeQZGjYL//d9EW8eOcNNNMH48HHBA1ktSaIuI7OqDD2DMGHjiieT2IUPgwQfhhBOiqQudPSIiktDYCPff7xONLQO7Z0945BF4+eVIAxvU0xYRccuWwZVXwuLFye0/+Qn88Y9w8MGRlLUr9bRFpLB9+SXccAOcckpyYB9/PLzyCvz1rzkT2KCetogUsuee84nGtWsTbUVFMGGCf3XqFF1t30ChLSKF56OPYOxYeOyx5PbBg30ZXw6fnaThEREpHI2NMHmyTzS2DOwePWDqVKiuzunABvW0RaRQrFzpE42LFiW3jxgBEyfCIYdEU1cbqactIvntq6/gxhuhoiI5sEtL4cUXYfr02AQ2qKctIvnspZfgqqugvj7R1qED/OpXHuSdO0dX2z5SaItI/tmwAcaN85P3WjrzTJ9o7NcvmrrSQMMjIpI/QvAJxbKy5MA+6CCfgFywINaBDeppi0i+WLXKh0IWLEhuv+wy+NOf4LDDoqkrzdTTFpF4+/prPzr1pJOSA/vYY2H+fPjb3/ImsEE9bRGJs+pq712vWZNo69DB72f8zW+gS5foassQhbaIxM8nn3gwT5uW3H766T523b9/NHVlgYZHRCQ+QvCgLitLDuziYj9SddGivA5sUE9bROJizRofCqmuTm6/5BL485/h8MOjqSvL1NMWkdy2bRvcdpv3oFsG9jHHwLx58PjjBRPYoNAWKXj19X46aXExtGvnv44albyJMDILFvj281tuge3bva19ex/PXrkSLrww2voioNAWKWDz58OAATBlCmze7EPGmzf76wED/P1IbNoEP/sZfPe7sHp1ov3UU6Gmxm+S6do1ouKilVJom9lBZjbLzFab2SozOyPThYlIZtXXw/DhfnFLQ0Pyew0N3j58eJZ73CH4TsayMt/Z2Kx7d7jnHnj9de95F7BUe9p/Bp4LIZQBJwGrMleSiGTDxIm7h/WuGhrg7ruzUw/vvAPnngsjR8LGjYn2iy+Gt9+Ga6/1oZEC12pom1kxMBiYChBC2B5C+CzThYlIZs2YkVpoT5+e4UK2b4ff/x5OPNFP5WvWqxc8/TQ8+aR/L0BqPe0+wEbgr2b2lplNMbPdBpPMrMrMasysZmPLPyVFJCdt2ZLe5/bJwoVw8slw002+SgR8NnTsWO9dDxuWwQ+Pp1RCuwMwEHgghHAysBW4YdeHQgiTQwiVIYTKkpKSNJcpIunWrVt6n2uTTz/1W2T+9V89nJsNHAhLlvgBT927Z+CD4y+V0H4feD+E8EbT61l4iItIjI0Y4ReP701RkQ8xp00IfjdjeblvN2/WtasPnr/xBpxySho/MP+0GtohhI+AdWb27aam7wFv7+W3iEgMjBuXWmiPHZumD1y7FoYO9aNSP/440T5smB+rOmaMH/Yke5Xq6pHrgEfNbBlQAdyeuZJEJBtKS2HWLD8Ib9fwLiry9lmz/Ln90tAAd97plw88/3yi/cgjYfZseOopOOqo/fyQwpFSaIcQ6prGqweEEC4KIXya6cJEJPOGDoVly6CqKrGarrjYXy9b5u/vl8WLfbjjhhv8gl0AM7juOh/L/tGP/LWkTH8XESlwpaVw772wYoW/fvXVNPzQzz+HCRPgwQd9HLtZRYXf0ThoUBo+pDBpG7uIpE8I8MQTvqPxgQcSgd2lC9x1FyxdqsDeT+ppi0h6vPuu71p85pnk9gsugPvug969o6gq76inLSL7Z8cO70X365cc2Icd5semzpunwE4j9bRFZN8tWeKbZOrqEm1mcPXVcPvtcOCB0dWWpxTaItJ2X3zhW8/vvTd5orF/f59oPEMHgWaKhkdEJHUhwJw50LevH5XaHNidO8Mdd0BtrQI7w9TTFpHUrFvnE41z5ya3n3eeX6rbp080dRUY9bRFZO927IBJk/y8kJaBfeihMHOmX2+jwM4a9bRF5JvV1vr2yDffTG6vqvLhkB49oqmrgKmnLSK727LFT4oaNCg5sPv29TOwH3pIgR0R9bRFJNncuT52vW5dou2AA+Dmm/0W9I4do6tNFNoi0mTbNnjn/8EPf5jcfs45viX9uOOiqUuSKLRFCt3Onb76Y2mFf9+spMQvJvjxj3USXw7RmLZIIaur83XV11+fHNhXXOEXE/znfyqwc4x62iKFaOtW+O1vvSfdFNYV1EGPb8FTr8HgwdHWJ99IoS1SaJ59FkaNgvfeS7R17MikmzbD+Kt90lFylkJbpFB8+CGMHu3nXbc0ZIhfVnDCCdHUJW2iMW2RfNfY6Ks/ysqSA7tnT3jkEXj5ZQV2jKinLZLPli/33YuLFye3X365n4F98MHR1CX7TD1tyXv19T6EW1wM7dr5r6NGeXve+vJLv0x34MDkwD7+eHjlFe9hK7BjSaEteW3+fBgwAKZMgc2b/STRzZv99YAB/n7eef55OPFEuPNOP+wJoKjIdzQuW+Zj2BJbKQ2PmNm7wGZgJ7AjhFCZyaJE0qG+HoYP907nrhoa/Gv4cM+x0tLs15d2H30EP/+5n7zX0uDBPtFYXh5NXZJWbelpDwkhVCiwJS4mTvRg3puGBl+qHGuNjTB5sodyy8Du0QOmToXqagV2HtHwiOStGTNSC+3p07NTT0asXOk96SuvhM8+S7SPGAGrV8NPf+oD+ZI3Uv1/MwAvmFmtmVXt6QEzqzKzGjOr2bhxY/oqFNlHW7ak97mc8tVXcOONUFEBixYl2ktL4cUX/U+iQw6Jrj7JmFRD+8wQwkBgKHCNme22xzWEMDmEUBlCqCwpKUlrkSL7olu39D6XM156yS/Qvf32xERjhw4e4suX+6l8krdSCu0QwgdNv24A5gCDMlmUSDqMGOGLJvamqAhGjsxOPfttwwYv9vvfT16veOaZfvDT737nF+xKXms1tM2sq5l1b/4eOBdYkenCRPbXuHGphfbYsdmpZ5+FAA8/7JOJM2Yk2g86yCcgFyyAfv2iq0+yKpWe9qHAQjP7B7AEeCaE8FxmyxLZf6WlMGsWdOmye3gXFXn7rFk5vtxv9Wo46yy44grGbPoNY2ha6nLZZX506n/9lyYaC0yr67RDCP8ETspCLSJpN3Sor8O++25fqrxzp++IHDnSe9g5G9hffw1/+IN/NS2BqaMCOnWCOfPh/PMjLlCiorNHJO+VlsK998KKpkG9V1+NtJzWVVfDVVfBmjWJtg4d4PCjofcxcH776GqTyOnvVSK54pNP4Cc/gbPPTg7s00/3G9H79IF2CuxCp9AWiVoIMG2aH506bVqivbjY725ctMiX+Img4RGRaK1Z40Mh1dXJ7ZdcAn/+Mxx+eDR1Sc5ST1skCtu2wW23eQ+6ZWAfcwzMmwePP67Alj1ST1sk2xYs8LNCVq9OtLVv78tZfvtb6No1stIk9ym0RbJl0yYYP95P3mvp1FN9k0xFRTR1SaxoeEQk00LwnYxlZcmB3b073HMPvP66AltSpp62SCa98w5cfbUf8tTSxRf7RGOvXtHUJbGlnrZIJmzfDr//vV/71TKwe/WCp5+GJ59UYMs+UU9bJN0WLvSJxrffTrS1awejR8Ott/qwiMg+UmiLpMunn/oN6JMnJ7cPHOhtp5wSTV2SVzQ8IrK/QoDHHvOjU1sGdteuflLVG28osCVt1NMW2R9r1/pE4/PPJ7cPG+anVB11VDR1Sd5ST1tkXzQ0wJ13+uUDLQP7yCNh9mx46ikFtmSEetoibbV4MVRV+X2Mzczg2mv9yq/i4uhqk7yn0BZJ1eefw4QJfptCCIn2igp46CEYpKtTJfM0PCLSmhDgiSd8R+MDDyQCu0sXuOsuWLpUgS1Zo562yN68+64PezzzTHL7BRfAffdB795RVCUFTD1tkT3ZscN70f36JQf2YYf5sanz5imwJRLqaYvsaskS39FYV5doM/OlfbffDgceGF1tUvDU0xZp9sUXcP31fidjy8Du39+v/LrvvkgCu74eRo3y3fGvveaLU0aN8nYpPCmHtpm1N7O3zGxeJgsSyboQYM4c6NvXj0ptnmjs3BnuuANqa+GMMyIpbf58GDAApkyBnTu9bfNmfz1ggL8vhaUtPe3RwKpMFSISiXXr4KKL/KjU9esT7eedBytWwK9+BUVFkZRWXw/Dh8OXX/penpYaGrx9+HD1uAtNSqFtZr2AC4EpmS1HJEt27IBJk/y8kLlzE+2HHgozZ3oXtk+f6OoDJk7cPax31dDgx5tI4Ui1pz0JGA80ZrAWkeyorYXTTvM7GbduTbRXVcGqVXDppT7xGLEZM1IL7enTs1OP5IZWV4+Y2Q+ADSGEWjM7ay/PVQFVAEcffXTaChRJm507/YCnQYOgsUX/o29fP53vzDOjq20PtmxJ73OSH1LpaZ8JDDOzd4HHgLPNbMauD4UQJocQKkMIlSUlJWkuU2Q/zZ0LS5fA+vcTgX3AAX67zFtv5VxgA3Trlt7nJD+0GtohhAkhhF4hhN7ApcArIYQRGa9MJB3Wr4d//3f44Q9h27ZE+znn+ETjr38NHTtGV99ejBjR+hxoURGMHJmdeiQ3aJ225KedO335Xnm5H5UKVFBHRef/8cHiF16A446LuMi9GzcutdAeOzY79UhuaNOOyBDCq8CrGalEJF3q6nxScenSpOZJV6zwM7B79oyosLYpLYVZs3xZX0ND8qRkUZF/zZrlz0nhUE9b8sfWrfDLX0JlZXJgl5X5VsIpU2IT2M2GDoVly/zPoOJivx+4uNhfL1vm70thsdDyXOA0qaysDDU1NWn/uSLf6NlnfW/3e+8l2jp2hJtugvHjfdJRJIeZWW0IobK153RglMTbhx/C6NF+3nVLQ4b4ZQUnnBBNXSIZouERiafGRr+QoKwsObB79oRHHoGXX1ZgS15ST1viZ/lyH9RdvDi5/fLL/Qzsgw+Opi6RLFBPW+Ljyy/hhhtg4MDkwD7+eHjlFe9hK7Alz6mnLfHw/PN+CcHatYm2oiK/aHfCBOjUKbraRLJIoS257aOP4Oc/95P3Who82Ccay8ujqUskIhoekdzU2OiHOJWXJwd2jx4wdSpUVyuwpSCppy25Z+VKv6Nx0aLk9hEj/JDpQw6Jpi6RHKCetuSOr76CG2+EiorkwC4t9bNCpk9XYEvBU09bcsNLL8FVVyXfndWhg1/3deONfl+jiCi0JWIbNvhxdjNmMAa/N2sSY+E73/Ex7X79Ii5QJLcotCUaIcBf/+oHPG3aBEAdFdC+A9z/EPzsZ346kogkUWhL9q1e7RONCxYkt5cc4mdcV/1LNHWJxIC6MpI9X38Nt9wCAwYkB/axx/rt53375uwtMiK5QqEt2VFdDSedBLfdljjNv317n2hcsQLOPz/a+kRiQsMjklmffAK/+AVMm5bcfvrp8NBD3usWkZSppy2ZEYIHdVlZcmAXF8P99/s6bAW2SJuppy3pt2aNr7murk5uv+QSmDQJjjgimrpE8oB62pI+27b5mHX//smBfcwxMG8ePP64AltkPym080B9vV+P2PLi11GjkjcXZtyCBb79/JZbYPt2b2vf3sezV66ECy/MYjEi+avV0DazTma2xMz+YWYrzezWbBQmqZk/34eGp0yBzZt9KHnzZn89YIC/n1GbNvlGmO9+19dfNzv1VKipgT/+Ebp2zXARIoUjlZ72NuDsEMJJQAVwvpmdntmyJBX19TB8uF/o0ryKrllDg7cPH56hHncI8OijPtE4dWqivXt3+Mtf4PXXvectImnVamgHt6XpZVHTV8hoVZKSiRN3D+tdNTTA3Xen+YPfeQfOPdePSt24MdF+8cXw9ttw3XU+NCIiaZfSmLaZtTezOmAD8GII4Y3MliWpmDEjtdCePj1NH7h9O9x+u080vvRSor1XL3j6aXjySf9eRDImpdAOIewMIVQAvYBBZnbirs+YWZWZ1ZhZzcaWvS/JmC1bWn+mLc/t1aJFfqHujTf6dnTwWc+xY713PWxYGj5ERFrTptUjIYTPgFeB3fYchxAmhxAqQwiVJSUlaSpP9qZbt/Q+t0effuqHO/3Lv/gqkGYDB8KSJfCnP/k4tohkRSqrR0rM7KCm7zsD5wCr9/67JBtGjPALyfemqAhGjtyHHx4CPPaY38M4eXKivWtXHyR/4w045ZR9+MEisj9S6WkfDlSb2TJgKT6mPS+zZUkqxo1LLbTHjm3jD167Fi64AC67DD7+ONE+bBisWgVjxvitMiKSda3+lxdCWAacnIVapI1KS2HWLF/W19CQPClZVORfs2b5cylpaPDhjltv9fsamx15JNxzD1x0EZil9X+DiLSNdkTG3NChsGwZVFUlVtkVF/vrZcv8/ZQsXuzDHTfckAhsM1++9/bb8KMfKbBFcoD+jpsHSkvh3nv9WGqAV19tw2/+/HOYMAEefNDHsZtVVPjRqYMGpbNUEdlP6mkXqhB87KS8HB54IBHYXbrAXXfB0qUKbJEcpJ52IXrvPbjmGnjmmeT2Cy6A++6D3r0jKUtEWqeediHZscP3vvftmxzYhx3mx6bOm6fAFslx6mkXiqVLfXayri7RZgZXX+1b0w88MLraRCRlCu1898UXcNNNPlPZcqKxf3+faDzjjOhqE5E20/BIPpszx4dC7rknEdidO8Mdd0BtrQJbJIbU085H69b5+uqnn05uP+88v1S3T59o6hKR/abQzichwPr1voxv69ZE+6GH+oW6//Ef2iAjEnMK7XxRWwtvAls2Ay0Cu6rKh0N69IiqMhFJI41px92WLX4i1KBBTYHdpG9fWLjQJxsV2CJ5Qz3tOJs7F6691sewgQrq/ACS237vt6B37BhxgSKSbgrtOFq/Hq6/HmbPTmqedM4zviX9uOMiKkxEMk3DI3Gyc6evty4vTw7skhK/MPKFFxTYInlOPe24qKvza7+WLEluv+IKuPNO6NkzmrpEJKvU0851W7fCL38JlZXJgV1WBq+9BlOmKLBFCoh62rns2Wdh1Cg/la9Zx46+LX38eDjggOhqE5FIKLRz0YcfwujR8MQTye1DhvhlBSecEE1dIhI5DY/kksZGX/1RVpYc2D17wiOPwMsvK7BFCpx62rli+XLfvbh4cXL75Zf7TTIHHxxNXSKSU9TTjtqXX/plugMHJgf28cfDK694DzvPA7u+3ofuFy70udXiYn9dXx91ZSK5R6EdpeefhxNP9CV7O3Z4W1ER3HyzX6U+ZEi09WXB/PkwYIAvgtm509s2b/bXAwb4+yKS0Gpom9lRZlZtZqvMbKWZjc5GYXnto4/gxz+G88+HtWsT7YMHwz/+AbfeCp06RVdfltTXw/Dh/peNhobk9xoavH34cPW4RVpKpae9AxgXQigHTgeuMbO+mS0rTzU2wuTJvqNx5sxEe48eMHUqVFf7ewVi4sTdw3pXDQ1w993ZqUckDloN7RDChyGEN5u+3wysAo7MdGF5Z+VK70lfeSV89lmifcQIWL0afvpTaFdYo1UzZqQW2tOnZ6cekThoU0qYWW/gZOCNPbxXZWY1ZlazcePG9FSXD776Cm68ESoqYNGiRHtpqZ8VMn06HHJIdPVFaMuW9D4nUghSDm0z6wY8CYwJIXyx6/shhMkhhMoQQmVJSUk6a4yvl17yC3Rvvz0x0dihg4f48uXw/e9HW1/EunVL73MihSCl0DazIjywHw0hzG7t+YK3YQOMHOmh3HIW7Tvf8YOffvc7v2C3wI0Y4Ytl9qaoyP9RiohLZfWIAVOBVSGEP2W+pBgLAR5+2CcTZ8xItB90kN8g8/e/Q79+0dWXY8aNSy20x47NTj0icZBKT/tMYCRwtpnVNX1dkOG64mf1ajjrLD8qddOmRPull8KqVb7bscAmGltTWgqzZkGXLruHd1GRt8+a5c+JiEtl9cjCEIKFEAaEECqavp7NRnGx8PXXcMstvhNkwYJE+7HH+s6QmTPhsMOiqy/HDR3q+4iqqnwnZLt2/mtVlbcPHRp1hSK5xUIIaf+hlZWVoaamJu0/N+dUV8NVV8GaNYm29u39fsabb/auoohICsysNoRQ2dpzOjBqX3zyiQfztGnJ7aef7mPXAwZEU5eI5D0NsrZFCB7UZWXJgV1cDPff7+uwFdgikkHqaadqzRofCqmuTm6/5BKYNAmOOCKaukSkoKin3Zpt2+C223yTTMvAPuYYmDcPHn9cgS0iWaOe9t4sWOBnhaxenWhr394XDv/2t9C1a2SliUhhUmjvyaZNfnHu1KnJ7aee6qf0VVREU5eIFDwNj7QUAjz6qE80tgzs7t3hL3+B119XYItIpNTTbvbOO3D11X7IU0s/+pEHdq9e0dQlItKCetrbt/spfP37Jwd2r17w1FMwe7YCW0RyRmH3tBct8onGlSsTbe3awfXX+4qR7t2jq01EZA8KM7Q//dRvQJ88Obl94EBvO+WUaOoSEWlFYQ2PhACPPeZHp7YM7K5d/SLCN95QYItITiucnvbatTBqFDz3XHL7v/0b3HsvHH10NHWJiLRB/ve0Gxrgzjv98oGWgX3EET7J+PTTCmwRiY387mkvXuwHMy9fnmgzg2uv9Su/ioujq01EZB/kZ2h//jlMmAAPPujj2M1OOsnHsgcNiq42EZH9kF/DIyH4/VTl5fDAA4nA7tIF7roLamoU2CISa/nT037vPbjmGnjmmeT2Cy6A++6D3r0jKUtEJJ3i39PesQMmToS+fZMD+7DD/NjUefMU2CKSN+Ld01661Cca6+oSbWZ+WcEf/gAHHhhdbSIiGRDPnvYXX/hW89NOSw7s/v19a/r99yuwRSQvtRraZvawmW0wsxXZKKhVc+b4UMg99yQmGjt3hjvugNpaOOOMaOsTEcmgVHrajwDnZ7iO1q1bBxddBBdfDOvXJ9rPOw9WrIBf/QqKiqKrT0QkC1oN7RDCAmBTFmrZs507/eLc8nLfvdjs0ENh5kyYPx/69ImsPBGRbErbRKSZVQFVAEena1t4ba0fnVpbm9xeVeXDIT16pOdzRERiIm0TkSGEySGEyhBCZUlJyf79sC1b/PLcQYOSA7tvX1i4EB56SIEtIgUp95b8zZ3rZ4OsW5doO+AAuPlm+MUvoGPH6GoTEYlY7oT2+vW+jG/27OT2c87xLenHHRdNXSIiOSSVJX8zgdeBb5vZ+2Z2RVor2LnTz7MuL08O7JISmDEDXnhBgS0i0qTVnnYI4bKMfXpdnU80LlmS3H7FFX4Gds+eGftoEZE4imZH5Nat8MtfQmVlcmCXlcFrr8GUKQpsEZE9yP6Y9rPP+rVf772XaOvYEW66CcaP90lHERHZo+yF9ocfwujR8MQTye1DhvhlBSeckLVSRETiKvPDI42NvvqjrCw5sHv2hEcegZdfVmCLiKQosz3t5ct99+Lixcntl1/uN8kcfHBGP15EJN9kJrQbG+GGG/xygh07Eu3HH++7GYcMycjHiojku8yE9sqV8NZbiddFRX7R7oQJ0KlTRj5SRKQQZCa0t29PfD94sE80lpdn5KNERApJ5iYie/SAqVOhulqBLSKSJpnpaX/rW7BqFRxySEZ+vIhIocpMT/vYYxXYIiIZEM+LfUVECpRCW0QkRhTaIiIxotAWEYkRhbaISIwotEVEYkShLSISIwptEZEYUWiLiMSIQltEJEYU2iIiMZJSaJvZ+Wb2P2b2jpndkOmiRERkz1oNbTNrD9wHDAX6ApeZWd9MFyYiIrtLpac9CHgnhPDPEMJ24DHgh5ktS0RE9iSV87SPBNa1eP0+cNquD5lZFVDV9HKbma3Y//IicTDwSdRF7AfVHy3VH6041//tVB5KJbRtD21ht4YQJgOTAcysJoRQmUoBuSbOtYPqj5rqj1ac6zezmlSeS2V45H3gqBavewEf7EtRIiKyf1IJ7aXA8WZ2rJl1BC4F5ma2LBER2ZNWh0dCCDvM7FrgeaA98HAIYWUrv21yOoqLSJxrB9UfNdUfrTjXn1LtFsJuw9MiIpKjtCNSRCRGFNoiIjGS1tCO83Z3M3vYzDbEdX25mR1lZtVmtsrMVprZ6Khragsz62RmS8zsH0313xp1TW1lZu3N7C0zmxd1LW1lZu+a2XIzq0t16VkuMbODzGyWma1u+m/gjKhrSpWZfbvpn3vz1xdmNuYbn0/XmHbTdvc1wPfxZYJLgctCCG+n5QMyzMwGA1uA/xtCODHqetrKzA4HDg8hvGlm3YFa4KIY/fM3oGsIYYuZFQELgdEhhMURl5YyM/s5UAkUhxB+EHU9bWFm7wKVIYRYbkwxs2nA30MIU5pWuXUJIXwWdV1t1ZSj64HTQgjv7emZdPa0Y73dPYSwANgUdR37KoTwYQjhzabvNwOr8N2ssRDclqaXRU1fsZklN0THi1MAAAIDSURBVLNewIXAlKhrKTRmVgwMBqYChBC2xzGwm3wPqP+mwIb0hvaetrvHJjTyiZn1Bk4G3oi2krZpGl6oAzYAL4YQ4lT/JGA80Bh1IfsoAC+YWW3TkRRx0gfYCPy1aXhqipl1jbqofXQpMHNvD6QztFPa7i6ZZWbdgCeBMSGEL6Kupy1CCDtDCBX4rttBZhaLYSoz+wGwIYRQG3Ut++HMEMJA/DTPa5qGC+OiAzAQeCCEcDKwFYjVnBpA07DOMOCJvT2XztDWdveINY0FPwk8GkKYHXU9+6rpr7avAudHXEqqzgSGNY0LPwacbWYzoi2pbUIIHzT9ugGYgw93xsX7wPst/mY2Cw/xuBkKvBlC+HhvD6UztLXdPUJNE3lTgVUhhD9FXU9bmVmJmR3U9H1n4BxgdbRVpSaEMCGE0CuE0Bv/9/6VEMKIiMtKmZl1bZq8pmlY4VwgNquoQggfAevMrPmUvO8BsZiA38VltDI0Aqmd8peSfdzunjPMbCZwFnCwmb0P3BJCmBptVW1yJjASWN40Lgzw6xDCsxHW1BaHA9OaZs/bAY+HEGK3dC6mDgXm+J/7dAD+FkJ4LtqS2uw64NGmDuM/gf8TcT1tYmZd8JV3V7b6rLaxi4jEh3ZEiojEiEJbRCRGFNoiIjGi0BYRiRGFtohIjCi0RURiRKEtIhIj/x8QKBFIxh98AAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.array([2,3,3.5,4.5,5])\n",
    "y = np.array([3,2,5,4,6])\n",
    "\n",
    "slope,intercept,r,p,se = stats.linregress(x,y)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x,y,'bo',ms=10)\n",
    "plt.plot(np.array([0,7]),slope*np.array([0,7])+intercept,'r-',lw=3)\n",
    "for xp,yp in zip(x,y):\n",
    "    plt.plot([xp,xp],[yp,slope*xp+intercept],'b')\n",
    "plt.xlim([0,7])\n",
    "plt.ylim([0,7])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Type I regression assumes that the x variable is exactly known (error free)\n",
    "\n",
    "```python\n",
    "from scipy import stats\n",
    "help(stats.linregress)\n",
    "```\n",
    "\n",
    "Examples of (nearly) error free measurements:\n",
    "- Time\n",
    "- Chemical standard (when compared to an error-prone sensor)\n",
    "- Distance\n",
    "\n",
    "Type 1 regression minimizes the sum of square error.\n",
    "\n",
    "Be careful when choosing a subset!\n",
    "\n",
    "![images/temperature_trends_short_vs_long.png](images/temperature_trends_short_vs_long.png)\n",
    "\n",
    "_source_: Walsh, J., et al. (11 January 2013), “Figure 6: Short-term Variations Versus Long-term Trend, in: D. Is the global temperature still increasing? Isn’t there recent evidence that it is actually 1 cooling?, in: Appendix I: NCA Climate Science - Addressing Commonly Asked Questions from A to Z”, in Federal Advisory Committee Draft Climate Assessment. A report by the National Climate Assessment Development Advisory Committee (NCADAC)[1], Washington, DC, USA: U.S. Global Change Research Program, p.1065.\n",
    "\n",
    "### Linear Regression Assumptions\n",
    "\n",
    "* Validity of linear model\n",
    "\n",
    "* Constant variance: same variance regardless of x value (homoscdastic)\n",
    "\n",
    "* Independence of errors (errors are uncorrelated)\n",
    "\n",
    "Anscombe's quartet shows data sets with the same mean and standard deviation in both variables, as well as the same regression line.\n",
    "\n",
    "The code for this example comes from:\n",
    "http://matplotlib.org/examples/pylab_examples/anscombe.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(3, 12, 'IV')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAfK0lEQVR4nO3deZhU1bnv8e8rg4gSAcFuBREDSEQckrTGKDkOKBKDoFdjIg6IAzHXKMSBSRyiImpC5KDGAVEkBqdcQjRHgh4RTXDsBoMoMqogkaZRjBgZu9f94+0G0vRY067a9fs8Tz127S56L4vFr1e9e621LYSAiIjEy25RN0BERFJP4S4iEkMKdxGRGFK4i4jEkMJdRCSGmmbyZO3atQudO3fO5Cklj5SUlKwLIbSP4tzq25JOifTtjIZ7586dKS4uzuQpJY+Y2cdRnVt9W9Ipkb6tsoyISAwp3EVEYih24W5mmFnUzRBJqdr6tfq71CZ24S4iIgp3EZFYUriLiMSQwl1EJIYU7iIiMaRwFxGJIYW7iEgMKdxFRGJI4S4iEkP1hruZPWJma81s4U7Hfm1mH5jZAjP7k5m1Tm8zRUSkMRoycp8C9K127EWgZwjhcGAJMCrF7RIRkSTUG+4hhFeBz6sdeyGEsK3y6RtAxzS0TUREEpSKmvvFwMzavmlmQ8ys2MyKy8rKUnA6keygvi3ZLKlwN7PrgW3AH2p7TQjhoRBCUQihqH379N8kJ4RACCHt5xHJZN+urV+rv0ttEr4Tk5kNAvoBvYN6l4hIVkko3M2sLzACOD6E8HVqmyQiIslqyFTIJ4DXge5m9omZXQLcC7QCXjSzd8zsgTS3U0REGqHekXsI4dwaDk9OQ1tERCRFtEI1RgoLC7ffdm3nR2FhYdRNE5EMU7jHSGlpaaOOi0h8KdxFRGJI4S4iEkMKdxGRGFK4i4jEkMI9RgoKChp1XETiK+HtByT7rFmzJuomiEiW0MhdAM2RlzyRR9tgKdwF2DEXvj/wLnBEteMiOe3LL2HMGDjxRKioiLo1GaFwFwCOBf4G/BnoCVwVbXNEUmPzZvjv/4YuXWDsWHjlFZg7N+pWZYTCPd8tWgRnnMFcoBdQBlwJXB5tq0SSU1EB06bBIYfAsGGwbh306gWvvQY/+EHUrcsIXVDNV//8J9x8M0yeDBUV/BsYD/wG2BBty0SS87//CyNGwLx5/rxHD7jjDujXD8yibVsGaeSeb/71L7j+eujaFSZN8s5++eV0BW5CwS7RSMkF/fnzoU8fOOUUD/YOHXzw8o9/wOmn51Wwg8I9f2zeDBMmeO3x9tth40Y46yx47z24/36C5shLmtUV4PVteldn+H/4IZx3HnznO/Dii7D33j5SX7IELr4YmuZngSI//6/zSUUFPPGEzxT46CM/9oMfwF13wTHHbH+Z5shLuiWza2lNr9kHGFlaCt27w9at0Lw5XHkljBoF++yTbHNznsI9zl54wWuP77zjzw891Ec0P/pR3n1ElfhoCQwDhgN7A2zbBhdcALfeCgceGGnbsonKMhmWkcVCJSVedzz1VA/2jh3hkUe89phnF5Ukc9Ldt5sAlwJLgbF4sM8Er7VPnapgr0bhnmFpvaHGihUwcCAUFfmMgb33hjvv9Nrj4MHQpEny5xCpRdr6dggMwBfXTQL2B94GTgJOAwpPPVWrq2ugskwclJXBbbfB/fd77XH33XfUHtu2jbp1IvUqKCio8ZdAvzZtoFcvZlQ+XwaMBv4IVG0koDuQ1Uwj91z27397qHfpAhMneu3xwgth8WL49a8V7JJV6tq1dM2aNYQQdjzee4/Qvz/PrV8Pr73GOjOuAHoAz7Aj2DWbq3YaueeYwsJCPist5WLgZmC/yuMvNW9O77ffhsMPj65xInVo0Iys1avhppvg0Ud9pteee8I119Du2mu5r1Ur7qvhj5iuIdVI4Z5LQuDY0lLGAd0rD70FjADmbNlCULBLrvriC5+eO2GCr8Fo0gR+/nO48UbI89p5olSWybD6bqhR24yDAW3bwrHHMh0P9qXAj4HvAXMy03TJc/XNhknoZjGbN8Pdd3tpcdw4D/azz4b334ff/S4lwZ6v21knNXI3s1/is5MCfjF7cAhhUyoaFlf1fTStfhGoBzAO6L9+PbzxBqXAr/BZA9vS1EaRmtR34bJRC+GqNvYaMwY+/tiPHX+8z+763vca1a7aLsZW/VLJ1wuuCY/czawDvjNsUQihJz4N9aepali+6wA8DCzA91j/CuCmm+gK3I+CXXJUCDBrlm8VcMEFHuw9e8Jf/gIvv9zoYAd2vRhb+cj3VdfJlmWaAnuYWVN84dg/k29SftsbH6kvBS4BKoD7gK4AN9/sIS+Si0pK4OSToW9fX1DXsaNfOH3nHa2aToOEwz2EsBrfIXYl8CnwrxDCC9VfZ2ZDzKzYzIrLysoSb2ncbdoE48ezAhgJ7AE8jZdlfgFUfYBM5ibY+Vp7TBf17QZavhzOPdcX182eDa1b+8XTJUvgoou0uC5davo405AH0AaYDbQHmgEzgPPr+jPf/e53g1SzbVsIjz0WQqdOIfiH1vAyhKP8OsZ/PJJV/eel8mdnA6A4JNifk33kQ99udP8pLQ3hyitDaNbM+/buu4dw3XUhfPZZdrc7CyXSt5Mpy5wMfBhCKAshbAWm43drk4YIAWbO9NrjoEGwciX07Ml5rVtzIr68emdarCFRa/Cnxq++8k28unSBe+7xxXUXXeQj9bvuyvjiumQ+7eayZMJ9JXCMmbU0X0XQG1iUmmZFL60ljLffht694bTTYMECOOAAmDIF3nmHP6xfr4tDEpm6+n29Fy63boUHHvAbwdx4o4f8aad5ff3RR6FTp0j+n/L1gmvCUyFDCG+a2R+BefjkjfnAQ6lqWNTSMn1q2TK/C9LTT/vzNm1g9Gj4xS+gRYvEf65IiiTU70OA6dO9Ly9Z4seOPtqnNZ5wQuobKQ2S1Dz3EMJN+N3ZZCfV7yzTHrgR+Bl+cYLdd4ehQ2HkSA94kVz16qswfDi8+aY/79bN7/R11lma/RIxrVBNg6pg3xMP9eX4jJcm4FvvLl3qo5oMB3u+1h4lDRYu9PuSHn+8B3tBga8ofe89X2GqYI+c9pZJg6bAZfhHmqrYfA4YBSx85JGomhX7GqNkwKpVvrHXY4/5KtO99oLrroOrr/avJWso3FMpBPjjH3kf6FZ56A18Y69Xo2uVSNJag9+yceJEX5PRtKlv7HXDDT5ql6yjskwtGl3CmDPHbzh9zjl0A5YAZwHfR8EuuaN6/94duAZYYebTGDdtgnPOgUWL4N57FexZTOFeiwZPn3r3XV86feKJ8NZbUFDA5cCh+MR/kVyyvd9v20aYMoVNBxzAb4A2Iezo40895dMdJasp3BO1cqUvzDjiCHj+ea833nILLFvGjIKCGjf20oVLyXoheH/+9re9f69a5TeAmTkTXnoJjjoq6hZKA6nm3liff+77Tt9zj+9FXVV7HDMG9t0X0IVLyV7Vp+lWKSgoYM2zz3pdfc4cP9ipk9/GceBA7f+SgxTuDbVxowf6uHF+1xiAn/wExo71ZdYiOaCmYO8KjC0t3bHdbps2vtjuiiu0uC6HKdzrU14OU6f6cupPPvFjJ53k89SLiqJtm0gS9sXXYQyhcnFdixY7Fte1bh1p2yR5CvfaVNUeR470BRvg9fU774Q+fbRIQ3LWXvgMmGsrvy4HJgOXLF3qe6xLLOTtBdU6NwZ7803fE6NfPw/2Aw+E3/8e5s2DU09VsEtu2rKF/wssA27Gg/1Z4HD8XpkK9njJ25F7TbXHbsDtpaU+Xx18a9IxY/yCqWqPkqtCgGeegdGjua/y0OvAcODvETZL0itvw31nBfhWAZdR+Ya0aAHDhvnMAdUeJZe9/LJv7FVcDMCyJk0YXl7On6q9TNN04ydvyzLgH0t/hW/s9XPA8JtSs2yZz4pRsEuuWrAAfvhDv/hfXAyFhfDAA3TdtInpebi3eT7Kz3DfsoUr8FC/Ed+98c/AYfjonQ4domubSDI+/hguvBCOPBL++ldo1crvirRsGfzsZ74uQ/JCfoV7RYUvne7Rg3vxqWCvAb2AM4jRbaQk/3z2GVx7LRx8sF/8b9oUrrrKb049ZgzsuWfULZQMy59f47Nne+2xpASApU2acF15OX+u9jLVHiWnbNzoOzWOGwf/+pcfO/dcX1n6zW9G2zaJVPxH7v/4B/Tt6/csLSnx2uODD9Jt0yZmJFh7TOv9VUUaorwcHnnE73w0cqQHe+/eXl+fNk3BLnkQ7ldfDbNmee3xttu89jhkSFK1x7TcX1WkMT780Pvx6tVeX581C158Eb77XQ0+BMiHsswdd8Af/uB7ZbRvH3VrRFKja1e/I9I3v+llmN12jNM0+BDIh3A/6ihtUyrxdMMNUbdAslj8yzIiInkotuGuuqOI5LPYhns6646Nvr+qiEiGJR3uZtbEzOab2V9S0aBc0OD7q4pEQIMPgdRcUB2KL+78Rgp+logkSYMMgSRH7mbWEfgRlfttiYhIdki2LDMB3xa6orYXmNkQMys2s+KysrIkTyeSPdS3JZslHO5m1g9YG0Ioqet1IYSHQghFIYSi9hlcRKS6o6RbVH1bpCGSqbkfB/Q3s9OAFsA3zOzxEML5qWlaclR3FJF8lvDIPYQwKoTQMYTQGfgpMDtbgl1EJN/Fdp67SFxpgZ40REr2lgkhzAHmpOJniUjdtDGYNIRG7iIiMaRwFxGJIYW7iEgMKdxFRGJI4S6SY7RATxoi/ndiEokZLdCThtDIXUQkhhTuIiIxpHAXEYkhhbuISAxZCCFzJzPbACzO2AlzXztgXdSNyCHdQwitojix+najqW83TqP7dqZnyywOIRRl+Jw5y8yK9X41nJkVR3h69e1GUN9unET6tsoyIiIxpHAXEYmhTIf7Qxk+X67T+9U4Ub5f+rtqHL1fjdPo9yujF1RFRCQzVJYREYkhhbuISAwp3EVEYkjhLiISQwp3EZEYUriLiMSQwl1EJIYU7iIiMaRwFxGJIYW7iEgMKdxFRGJI4S4iEkMKdxGRGFK4i4jEUEZvs9euXbvQuXPnTJ5S8khJScm6EEL7KM6tvi3plEjfzmi4d+7cmeLiKG9zKXFmZh9HdW71bUmnRPq2yjIiIjEUu3A3M8ws6eMi2aR6Px04cCBmxv3331/vnz3llFMwM2bMmJHOJkqWiV24i+SDIUOGADBp0qQ6X/fRRx/x0ksvsd9++9GvX79MNE2yhMJdJAedcMIJHHzwwcyfP5958+bV+rqHH36YEAKDBw+madOMXmKTiCncRXLUZZddBtQ+ei8vL2fKlCmYGZdeemkmmyZZQOGeRwoLC7fXbnd+FBYWRt00ScCgQYNo3rw506ZN4+uvv97l+88//zyrV6/m5JNP5qCDDoqghRIlhXseKS0tbdRxyW7t27fnjDPO4Msvv+Tpp5/e5fsPP/wwsKM+L/lF4S6Sw6qCuyrIq3z66ac8//zzFBQUMGDAgCiaJvVI9ydphbtIDjvppJPo0qULc+fOZdGiRduPP/roo2zbto2LLrqIZs2aRdhCqU26P0kr3EVy2M4XS6tG7yEEJk+erAupOaAlMBr4G6kPY4W7SI4bPHgwzZo1Y+rUqWzZsoXZs2ezYsUKTjzxRLp27Rp186Qm27ZxKbAUGAv0Avqk+BT1hruZPWJma81s4U7Hfm1mH5jZAjP7k5m1TnG7JA0KCgoadVxyQ0FBAf3792fdunXMmDFj+9RIXUjNQiHAjBlw2GFMAvYHioHewF9TfKqGjNynAH2rHXsR6BlCOBxYAoxKcbskDdasWUMIYZfHmjVrom6aJKlqzvv48eOZMWMG7dq148wzz4y4VfIf5s6FXr3gzDPhgw9YDvwEOBqYnYbT1RvuIYRXgc+rHXshhLCt8ukbQMc0tE1EGqhPnz4cdNBBvPXWW2zevHn7HHjJAu+/DwMGeLC/9hq0awcTJ3LCvvvyNBCqvTxVn6RTUXO/GJhZ2zfNbIiZFZtZcVlZWQpOJ5IdsqlvmxmXXHLJ9udVI3mJ0OrVcOmlcNhh8Oyz0LIl3HADLF8OV17JqtLStH6SthCq/96o4UVmnYG/hBB6Vjt+PVAE/J/QgB9UVFQUtOe1pIuZlYQQiqI4t/q2bPfFF3DXXTBhAmzcCE2awGWXwU03QYJz2BPp2wnvJGRmg4B+QO+GBLuISKxt3gy/+x3cdht8XlnJPussuP12OPjgjDcnoXA3s77ACOD4EMKum1qIiOSLigqYNg3GjIGPK2+Y9F//BXfeCcccE1mzGjIV8gngdaC7mX1iZpcA9wKtgBfN7B0zeyDN7RQRyS4hwKxZ8J3vwAUXeLAfeig89xzMmRNpsEMDRu4hhHNrODw5DW0REckNJSUwfDjMrpzE2LEj3HILXHih19izgHbvFxFpqOXLvfzy5JP+vHVrGDUKrrwS9tgj2rZVo3AXEanP2rV+ofSBB2DrVth9dw/0UaOgbduoW1cjhbuISG2++gruvtunNn71FZjBoEFegunUKerW1UnhLiJS3datMHky3HwzVG3Be9ppMG4cHH54pE1rKIW7iEiVEGD6dBg9GpYs8WNHHeUj9xNOiLRpjaVwFxEBePVVnwHz5pv+vGtXX4B09tlejskx2s9dRPLbwoVw+ulw/PEe7PvuC/fd5xt+/fjHORnsoJG7iOSrVat8v5fHHvNVpnvuCdddB9dcA3vtFXXrkqZwF5H8sn493HEHTJwImzZB06bw85/7jo0xunGNyjKyq9deg23b6n+dSC7ZtAl+8xvo0sUvkG7aBOec4+WXe++NVbCDwl12tnAh9OsHxx0HU6ZE3RqR1Cgvh6lToXt3L7usX+8zX956C556Crp1i7qFaaGyjHjt8cYbvfYYgtcbN22KulUiyQkB/vpXGDEC3n3Xjx12mO/W2Ldvzl4obSiFez5bv94XZUyc6HtRV9Uex4zxGQMiuertt31a45w5/rxTJ7j1VjjvvKzZ2CvdFO75aNMmuOcen8P7xRd+7Cc/gbFjvR4pkquWLfMFSM8848/btIHrr4crroAWLaJtW4Yp3PNJeTk8/rjPCli1yo+ddJJ/TC2K5O50IqlRWur7vTz0kE8GaNEChg71kkybNlG3LhIK93wQAsycCSNH7qg9HnGEh3qfPrGvPUqMbdgAv/2tz4L56ivYbTe4+GLfE+aAA6JuXaQU7nH31ltee3zlFX9+4IG+denAgf4PQSQXbd0KkybBr37l2/GCrzK9/Xbo2TPatmUJ/euOq6VLfen0977nwd62rY9wPvgAzj9/l2AvLCzEzHZ5FCZ4t3aJn6zoIyF4Pb1HD6+jr13rt7N75RV49lkF+040co+b0lIfzUyatKP2OGyY1x5bt67jj5U26rjkn8j7yJw5/in07bf9+cEH+2yvM89UabEGCve42LDB647jx8O//+0j80su8aDv0CHq1okkbsECv140c6Y/Lyz0mvrFF0OzZpE2LZsp3HPdli0+Sr/llh21x/79fUTTo0e0bRNJxsqVvrhu6lQvx7Rq5SP3X/7SN/mSOincc1VV7XH0aL9pL8D3v+97ZvTqFW3bRJLx+ec+OLnnHl9c16zZjsV17dtH3bqcoXDPRS+/7COY4mJ/3r2773I3YIBqj5K7Nm70QB83bsfiunPP9ZWlWlzXaJotk0sWLIAf/tAXHhUXw377wYMP+oZfZ5xB4X77JTyboaCWHfFqOy6SMuXl8OijfoF0xAgP9t69vY9Pm6ZgT1BSI3cz+yVwKRCAd4HBIQTtOJVqH3/sq0off3xH7XHECJ8Fs1PtMZnZDGvWrElZc0UaJAT4n//xi6XvvefHjjzSF9edcoo+hSYp4ZG7mXUArgKKQgg9gSbAT1PVMAE++wyuvdZHNL//vW/sNXQorFjh+2XoopLkqjfe8G13Tz/dg71zZx+8lJRo1XSKJFuWaQrsYWZNgZbAP5NvkrBxo49eunTxqY1btviK0sWLKXzySax9ey02kty0eLHfcPr73/cbUu+zD9x9ty+uO+88rZpOoYTLMiGE1Wb2G2AlsBF4IYTwQvXXmdkQYAhAp06dEj1dfigv9z3Vb7wRVq/2Y6ec4kH/7W8DWbCQRLZT326ETz/16bqTJnk/32MPn9I4fDjsvXfUrYulZMoybYABwEHA/sCeZnZ+9deFEB4KIRSFEIraaxpTzUKA556Dww/3hUerV3uYv/CCPyqDXbKL+nYDfPmlD1a6doUHHvC+fumlvj3G2LEK9jRK5jPQycCHIYSyEMJWYDpwbGqalUdefx2OP94XHr3/Pit3242BwG7z52N9+mjGi+SmLVt8WmPXrj6V8euvfaruwoU+eteq6bRLJtxXAseYWUszM6A3sCg1zcoDixfDWWfBscfC3/7mtccJE+hWUcET+PSjnTV0xksIYZeHZsJIxlRUwJNPwiGHwFVXQVmZ9/G//x1mzPDjkhEJh3sI4U3gj8A8fBrkbsBDKWpXfH36KVx+ORx6KEyf7rXH66/3VaZDh7Il6vaJJOqll+Doo33h0YoV8K1veaD//e9+03XJqKQuTYcQbgohfCuE0DOEcEEIYXOqGhY7X34JN9zA1/vvDw8+yLbych4E9t+4ERs7lsLu3Rv0Y1R6kWxzBPgNp08+2acy7r+/l17efVerpiOk7QfSbcsWv5B0662wbh0t8YsTo4HFO72sobNdVGKRbHEgcCtwHsCsWfCNb+xYXNeyZbSNE4V72lRUwFNPecnlww/92HHHcezcubwebctEktIWuB64Atgd2AI0HzbM+3q7dpG2TXbQioE0+HGbNpQ0aeILjz78kPeB/kDh0qX1BrvKLpLtpgJX48H+ONAdfCGSgj2raOSeSvPnw8iRPFO5o91q4EbgMaAcduy3XgeVXSTb3YYHx0jgnYjbIrVTuCegsLDwP2rknfHaY9UKri+AO4CJ+NJdkTh5A+gbdSOkXirLJKAq2PcBfgt8gAf7ZoCrr6YLcCc1B7vKLiKSCQr3BOwBjAKWA78EmuF1yO4A48fzeR1/VguNRCQTVJapRfXSC/iexkNbtWIZvpkOwEy89rggs80TEamTRu61qB7sA/BluOM3bGB/oBg4CTiNXYNdpRcRiZpG7vU4FrgLqFo8vRxfgPQMu+7/UkUlFomjgoKCGhfbadCSnfJ25F5YWFjn/Ua/BfwJmIsH+1rgF8AhwCsFBTUGuzq5iGSLvB2517bcv0lpKVx2GQvxGvtXwPjKx4bK12hkLvlIN4rJLXk7cq/uG8BYYCnAww8TgN8BXYGb2RHsIiK5ILbhXl/ZpUpzYBiwAq+ltwQ4+2z+a599uAKoPiZR6UVEckFsyzL1fYQ0YCC+lLpz5fdeAYYDbz7zDK+lvYUiIukT23Cv06xZzAOOrHz6Lj5X/fnoWiQiklKxLcvU5DvAiwB9+3IksAq4CA/5qmBX2UVE4iAvwv0gYBpQgt/Vm9at4a67OODrr5kSAuXaBkBEYib2ZZk78P1fmgOb8J0ahy9fDm3bRtouEZF0iu3Ivaq8EvDfYI8CBwO/LShQsItI7MV25L69vPLFF7BqFYMPO4zB0TZJJKdp+4HcEttw3651a3+ISFJ0PSq3xLYsIyKSzxTuIiIxlHS4m1kTM5tvZn9JRYNERCR5qRi5DwUWpeDniIhIiiQV7mbWEfgR8HBqmiMiIqmQ7Mh9Ar7XVkVtLzCzIWZWbGbFZWVlSZ5OJHuob0s2SzjczawfsDaEUFLX60IID4UQikIIRe3bt0/0dCJZR31bslkyI/fjgP5m9hHwJHCSmT2eklaJiEhSEg73EMKoEELHEEJn4KfA7BDC+SlrmYiIJEzz3EVEYigl2w+EEOYAc1Lxs0REJHkauYuIxJDCXUQkhhTuIiIxpHAXEYkhhbuISAwp3EVEYkjhLiISQwp3EZEYUriLiMSQwl1EJIYU7iIiMWQhhMydzGwDsDhjJ8x97YB1UTcih3QPIbSK4sTq242mvt04je7bKdk4rBEWhxCKMnzOnGVmxXq/Gs7MiiM8vfp2I6hvN04ifVtlGRGRGFK4i4jEUKbD/aEMny/X6f1qnCjfL/1dNY7er8Zp9PuV0QuqIiKSGSrLiIjEkMJdRCSGMhLuZtbXzBab2TIzG5mJc+YaM3vEzNaa2cKdjrU1sxfNbGnlf9tE2cZsYmYHmNnLZrbIzN4zs6GVxzP6nqlv1039unFS2a/THu5m1gS4D/gh0AM418x6pPu8OWgK0LfasZHASyGEbsBLlc/FbQOuCSEcAhwDXFHZrzL2nqlvN8gU1K8bI2X9OhMj96OBZSGEFSGELcCTwIAMnDenhBBeBT6vdngA8Fjl148BZ2S0UVkshPBpCGFe5dcbgEVABzL7nqlv10P9unFS2a8zEe4dgFU7Pf+k8pjUryCE8Cn4Xzqwb8TtyUpm1hn4NvAmmX3P1LcTo37dAMn260yEu9VwTPMvJSXMbC/g/wHDQghfZvr0NRxT35akpaJfZyLcPwEO2Ol5R+CfGThvHJSa2X4Alf9dG3F7soqZNcP/AfwhhDC98nAm3zP17cSoX9chVf06E+H+NtDNzA4ys+bAT4FnM3DeOHgWGFT59SDgzxG2JauYmQGTgUUhhN/u9K1Mvmfq24lRv65FKvt1RlaomtlpwASgCfBICGFs2k+aY8zsCeAEfCvUUuAmYAbwNNAJWAn8OIRQ/eJUXjKzXsDfgHeBisrDo/H6ZMbeM/XtuqlfN04q+7W2HxARiSGtUBURiSGFu4hIDCncRURiSOEuIhJDCncRkRhSuIuIxJDCXUQkhv4/6p3jfNTuvdQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5])\n",
    "y1 = np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])\n",
    "y2 = np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74])\n",
    "y3 = np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73])\n",
    "x4 = np.array([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8])\n",
    "y4 = np.array([6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89])\n",
    "\n",
    "def fit(x):\n",
    "    return 3 + 0.5*x\n",
    "\n",
    "xfit = np.array([np.amin(x), np.amax(x)])\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(221)\n",
    "plt.plot(x, y1, 'ks', xfit, fit(xfit), 'r-', lw=2)\n",
    "plt.axis([2, 20, 2, 14])\n",
    "plt.setp(plt.gca(), xticklabels=[], yticks=(4, 8, 12), xticks=(0, 10, 20))\n",
    "plt.text(3, 12, 'I', fontsize=20)\n",
    "\n",
    "plt.subplot(222)\n",
    "plt.plot(x, y2, 'ks', xfit, fit(xfit), 'r-', lw=2)\n",
    "plt.axis([2, 20, 2, 14])\n",
    "plt.setp(plt.gca(), xticklabels=[], yticks=(4, 8, 12), yticklabels=[], xticks=(0, 10, 20))\n",
    "plt.text(3, 12, 'II', fontsize=20)\n",
    "\n",
    "plt.subplot(223)\n",
    "plt.plot(x, y3, 'ks', xfit, fit(xfit), 'r-', lw=2)\n",
    "plt.axis([2, 20, 2, 14])\n",
    "plt.text(3, 12, 'III', fontsize=20)\n",
    "plt.setp(plt.gca(), yticks=(4, 8, 12), xticks=(0, 10, 20))\n",
    "\n",
    "plt.subplot(224)\n",
    "xfit = np.array([np.amin(x4), np.amax(x4)])\n",
    "plt.plot(x4, y4, 'ks', xfit, fit(xfit), 'r-', lw=2)\n",
    "plt.axis([2, 20, 2, 14])\n",
    "plt.setp(plt.gca(), yticklabels=[], yticks=(4, 8, 12), xticks=(0, 10, 20))\n",
    "plt.text(3, 12, 'IV', fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minimizing the sum of squared errors, SSE\n",
    "\n",
    "The goal of any \"least squares\" regression is to minimize the sum of squared errors.\n",
    "\n",
    "$$SSE  = \\sum_{i=1} ^N (y_i - \\hat{y_i})^2$$ \n",
    "\n",
    "Where:\n",
    "\n",
    "* $y_i$ is the y value of the data point\n",
    "* $\\hat{y_i}$ is  the y value predicted by the linear model (i.e. the y position of the line at the same x position as the data)\n",
    "* $(y_i - \\hat{y_i})^2$ is the squared error\n",
    "\n",
    "On an $x$-$y$ plot, the SSE can be thought of the sum of the vertical distances between the data values $y_i$ and the corresponding model values $\\hat{y}_i$.\n",
    "\n",
    "### Regression slope\n",
    "\n",
    "The regression slope depends on the covariance $s_{xy}^2$ and the variance $s_x^2$.\n",
    "\n",
    "$$ \\hat{a}_2 = \\frac{s_{xy}^2}{s_x^2} = \\frac{ \\sum_{i=1} ^N (x_i - \\bar{x})(y_i - \\bar{y})}{\\sum_{i=1} ^N (x_i - \\bar{x})^2}$$\n",
    "\n",
    "The regression slope can also be throught of as scaling the correlation coefficient by the standard deviations $s_y$ and $s_x$.\n",
    "\n",
    "$$ \\hat{a}_2 = \\frac{s_{xy}^2}{s_x^2} = r_{xy}\\frac{s_y}{s_x}$$\n",
    "\n",
    "Either way, the slope depends on both the spread of the data and the degree to which $x$ and $y$ covary. If the correlation between $x$ and $y$ is zero, the least squares regression line is flat. The regression line can have a steep slope if there is a strong correlation between $x$ and $y$, and the standard deviation of $y$ is large relative to the standard deviation of $x$.\n",
    "\n",
    "### Regression Intercept\n",
    "\n",
    "The least squares regression line is guaranteed to go through the mean values of $x$ and $y$. This property can be used to find the intercept of the regression line.\n",
    "\n",
    "$$ \\hat{a}_1 = \\bar{y} - \\hat{a}_2\\bar{x} $$\n",
    "\n",
    "### Standard error of the estimate\n",
    "\n",
    "The standard error of the estimate is a measure of the goodness of fit of the model. It depends on the sum of squared errors, which is minimized in least squares regression, as well as the degrees of freedom $\\nu = N - 2$.\n",
    "\n",
    "$$s_e = \\left( \\frac{SSE}{N-2} \\right)^{1/2}$$\n",
    "\n",
    "Note that this is different from the standard error of the data, $SE_x = s_x/\\sqrt{N}$, which was used earlier to describe how well the mean of a variable $x$ is known.\n",
    "\n",
    "### Standard error of the regression slope\n",
    "\n",
    "The standard error of the regression slope (sometimes called the standard error of the gradient) describes how well $a_2$ is known. It depends on the standard error of the estimate, and the variance in $x$.\n",
    "\n",
    "$$s_{a_2} = \\sqrt{\\frac{s_e^2}{\\sum_{i=1}^N (x-\\bar{x})^2}}$$\n",
    "\n",
    "Note that the `stderr` result given by the `stats.linregress` function is Python refers to the standard error of the slope, $s_{a_2}$.\n",
    "\n",
    "### Confidence intervals for the regression slope\n",
    "\n",
    "The standard error of the slope $s_{a_2}$ can be used to calculate confidence intervals for the regression slope. \n",
    "\n",
    "$$\\hat{a}_2 - t_{\\nu,1 - \\alpha /2} s_{a_2} < a_2 < \\hat{a}_2 + t_{\\nu,1- \\alpha /2} s_{a_2}$$\n",
    "\n",
    "Based on the data, these intervals describe a range in which the true slope $a_2$ is thought to lie. This depends on the estimated slope $\\hat{a}_2$, as well as the degrees of freedom $\\nu$ and significance level $\\alpha$. This is similar to how confidence intervals are calculated for an estimate of the mean, $\\bar{x}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Review: relationship between correlation and regression ###\n",
    "\n",
    "The concepts of correlation and regression are closely related. The correlation coefficient $r$ is given by the relationship\n",
    "\n",
    "$$ r = \\frac{s_{xy}}{s_x s_y}$$\n",
    "\n",
    "where $s_{xy}$ is the covariance between the $x$ variable and the $y$ variable, $s_x$ is the standard deviation of the $x$ variable and $s_y$ is the standard deviation of the $y$ variable.\n",
    "\n",
    "In the linear model\n",
    "\n",
    "$$ \\hat{y} = \\hat{a}_1 + \\hat{a}_2 x,$$\n",
    "\n",
    "the ordinary least squares regression slope $\\hat{a}_2$ is given by the relationship\n",
    "\n",
    "$$ \\hat{a}_2 = \\frac{s_{xy}}{s_x s_x}$$\n",
    "\n",
    "This means that the relationship betweeen the correlation coefficient $r$ and the least squares regression slope $\\hat{a}_2$ is \n",
    "\n",
    "$$ r = \\hat{a}_2 \\frac{s_y}{s_x}.$$\n",
    "\n",
    "This has several important implications:\n",
    "\n",
    "* $r$ and $a_2$ always have the same sign ($r$ is negative when the regression line has negative slope)\n",
    "\n",
    "* If $r = 0$, then  $\\hat{a}_2 = 0$ and $s_{xy} = 0$\n",
    "    \n",
    "There is also a close relationship between testing the null hypothesis of zero correlation and the confidence intervals on the regression slope. If the null hypothesis can be rejected at a given confidence level, then the confidence intervals of the slope will not include zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear regression (Type II) ###\n",
    "\n",
    "Type 2 regression is the case where there are potentially errors or uncertainty in both the x and y variables. In a Type 1 regression, one variable needs to be selected as the independent variable ($x$) and another needs to be the dependent variable ($y$). The least-squares method in the Type I approach minimizes the sum of squared errors $\\sum(\\hat{y}-y)^2$. Errors in the $x$ variable are not considered in the Type 1 approach.\n",
    "\n",
    "![images/type_2_reg.png](images/type_2_reg.png)\n",
    "\n",
    "Image: Emery and Thomson\n",
    "\n",
    "Reference for Geometric Mean Function Regression (GMFR, a.k.a. neutral regression):\n",
    "\n",
    "Ricker, W. E. Computation and uses of central trend lines\n",
    "Can. J. Zool., 1984, 62, 1897-1905 \n",
    "\n",
    "__Calculating Geometric Mean Function Regression__\n",
    "\n",
    "- $\\hat{a}_{2yx}$ : slope of regression of y on x\n",
    "\n",
    "- $\\hat{a}_{2xy}$ : slope of regression of y on x\n",
    "\n",
    "- Geometric mean: $\\hat{a}_{2GM} = \\sqrt{\\frac{\\hat{a}_{2yx}} {\\hat{a}_{2xy}}}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simplest Type 2 regression - the geometric mean\n",
    "\n",
    "The difference between the two Type 1 regression lines above can be used to illustrate one potential solution, the geometric mean regression. This is the line that bisects the two Type 1 regression lines. This can be thought of as a \"compromise\" between the two lines. An important assumption in this approach is that the uncertainty in each variable is assumed to be proportional to the sample standard deviation. If the uncertainty is dominated by natural variability, this might be a reasonable assumption. In cases, where the uncertainty is dominated by measurement error, there might be better ways of quantifying the uncertainty based on the analytical techniques. In this case, more sophisticated iterative approaches can be taken to determine the appropriate Type 2 regression model.\n",
    "\n",
    "Detailed discussion can be found in Ricker, W. E. (1984) Computation and uses of central trend lines, *Can. J. Zool.*, 62, 1897-1905.\n",
    "\n",
    "Glover, Jenkins and Doney provide an algorithm for determining the geometric mean regression slope:\n",
    "\n",
    "* Determine the slope of the regression of $y$ on $x$, $slope_{yx}$\n",
    "\n",
    "* Determine the slope of the regression of $x$ on $y$, $slope_{xy}$\n",
    "\n",
    "* Calculate the geometric mean regression slope, $slope_{GMR} = \\sqrt{\\frac{slope_{yx}}{slope_{xy}}}$\n",
    "\n",
    "* Determine the intercept by noting that the regression line goes through the mean of $x$ and the mean of $y$"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}