Linear_Regression_and_Confidence_Intervals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Regression and Confidence Intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#simulate a random set of (x,y) points with \n",
    "import numpy as np\n",
    "\n",
    "#will simulate n many points\n",
    "n = 45\n",
    "# true values are slope = 4, intercept = 3, with random error added as well, drawn from N(0,sigma2)\n",
    "b = 3\n",
    "m = 4\n",
    "sigma2 = 1.5\n",
    "flat = 2\n",
    "X = flat* np.random.rand(n, 1) #rand is flat, between [0,flat]\n",
    "y = m * X + b + (np.sqrt(sigma2)*np.random.randn(n, 1)) #randn is from a standard normal dist N(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEVCAYAAAAGrllxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFPFJREFUeJzt3X+Qb3dd3/HnK/fGhBLjkGbNoLK5RUURi1F3tClFMk0qDdWBEdQoVJxSrgMTUVEKncklNz8wIyoMUxC9Y0KQWjF2YqogzFRqtK2X4mYY2kZCRggEjMGbH425Ibn59e4f391272b37ve793zPOd9zno+Znb37/X52P5899+x5nfP5fM7npKqQJI3TKV03QJLUHUNAkkbMEJCkETMEJGnEDAFJGjFDQJJGzBCQpBEzBCRpxAwBSRqxvV03YCdnn3127du3r+tmSNJCueWWW+6pqqWdyvU+BPbt28fq6mrXzZCkhZLkC9OUsztIkkbMEJCkETMEJGnEDAFJGjFDQJJGzBCQpBEzBCRpxAwBSRoxQ0CSRqzREEhyaZLVJMeSXL9NmcuTVJKLmqxbkjS7ppeNuAu4Gngx8LTNbyb5RuAVwN80XK8kaRcavRKoqhur6ibg3m2KvBt4M/Bok/VKknantTGBJD8MPFpVf9RWnZKkE2tlFdEkZwC/CHz/lOX3A/sBlpeX59gySRq3tq4ErgA+UFV3TFO4qg5V1UpVrSwt7bgctiRpl9oKgQuBNyS5O8ndwLOAG5K8uaX6JUlbaLQ7KMnetZ+5B9iT5HTgcSYhcOqGon8BvBH4SJP1S5Jm0/SVwGXAw8BbgFet/fuyqrq3qu5e/wCeAO6vqqMN1y9JmkGjVwJVdRA4OEW5fU3WK0naHZeNkKQRMwQkacQMAUkaMUNAkkbMEJCkETMEJGnEDAFJGjFDQJJGzBCQpBEzBCRpxAwBSRoxQ0CSRswQkKQRMwQkacQMAUkaMUNAkkbMEJCkETMEJGnEDAFJGrFGQyDJpUlWkxxLcv2G1/9Rkv+c5L4kR5L8XpJnNlm3JGl2TV8J3AVcDVy36fVnAIeAfcC5wIPA+xquW5I0o71N/rCquhEgyQrwDRte/8jGckneDfxpk3VLkmbX1ZjA9wG3dlS3JGlNo1cC00jyfOCtwEtPUGY/sB9geXm5pZZJ0vi0eiWQ5JuAjwA/U1X/dbtyVXWoqlaqamVpaam9BkrSyLQWAknOBf4YuKqqPtBWvZKk7TXaHZRk79rP3APsSXI68DhwDvBfgPdU1a83WackafeaHhO4DLh8w9evAq4ACng2cHmS//d+VZ3RcP2SpBk0PUX0IHBwm7evaLIuSZq3w4fh5pvhggvg/PO7bs18tD47SJIWweHDcOGF8Oij8FVfBR/72DCDwLWDJGkLN988CYAnnph8vvnmrls0H4aApF46fBiuuWbyuQsXXDC5AtizZ/L5ggu6ace82R0kqXf60BVz/vmTeh0TkKSWbdUV08VB+Pzzh3vwX2d3kKTemUdXTNfdS33llYCk3mm6K6YP3Ut9ZQhIalwT8+ub7IqZZ/fSot9LYAhIalQfz7rXu5fW29TUTJ8+/q6zckxAUqP6ML9+c///evfSVVc1e6Duw+96srwSkNSoeZ11T2u7s/N5zPTp+ndtgiEgqVFdz6+fR///dv3+Xf+uTTAEJDWuy/n1uz073+5Av1O///q/17uCFi0IDAFJg7Kbs/MTHeh3urJY9MFhQ0DS4Mx6JXKiA/1OVxZ9ubt5twwBSaN3ogP9TlcWiz44nKrqug0ntLKyUqurq103Q9LAncxNX328YSzJLVW1smM5Q0CShmfaEPBmMUm95sJv89VoCCS5NMlqkmNJrt/03oVJbkvylSR/kuTcJuuWNDzrM28OHJh8Ngia1/SVwF3A1cB1G19McjZwI3AAOAtYBX634bolDcwQlmXou0ZDoKpurKqbgHs3vfVDwK1V9XtV9QhwEPiOJN/aZP2ShqVvj3gcYtdUW1NEnwd8av2LqnooyWfXXr+tpTZIWjB9WpZh0W8K205bIXAGcGTTaw8AX71V4ST7gf0Ay8vL822ZpF7ryyMeF/2msO20NTvoKHDmptfOBB7cqnBVHaqqlapaWVpamnvjJGknfeuaakpbVwK3Aq9e/yLJ04FvXHtdknqvT11TTWo0BJLsXfuZe4A9SU4HHgd+H/jlJC8HPgy8FfifVeV4gKSF0ZeuqSY13R10GfAw8BbgVWv/vqyqjgAvB94G3A98L3BJw3VLkmbU6JVAVR1kMv1zq/f+GHBKqCT1iMtGSNKIGQKSNGKGgCSNmCEgaWEMcdmGrvlkMUkLYajLNnTNKwFpBIZwBu2KovPhlYA0cH08g97N4xgX/Vm+fWUISAPXt4XPdhtKQ122oWuGgDRwfTuDPplQGuKyDV0zBKSB69sZdN9CaewMAWkE+nQG3bdQGjtDQFLr+hRKY+cUUUkaMUNAkkZsqhBI8utJKsnXbfHetyR5NMm7mm+eJGmepr0SWL/P8Hu2eO+dwN+xzXMEJEn9NW0IfHzt83EhkORfABcDb62q+5tsmCRp/qYKgar6DHAfG0IgyanAO4D/DfzGXFonSZqrWaaIfhx4QZJUVQE/AzwHuKiqnphL67RwdrMmjKTuzBoCLwG+Jcl9wAHgpqr62LQ/IMk+4NeA84FjwH8EfraqHp+hHeqpPi5UJunEZpkiunFw+BeB04Cfn7G+XwP+FngmcB7wIuD1M/4M9ZRL/Q7PEJag1onNciXwP4AngdcA/wT45ar63Iz1/QPg3VX1CHB3ko8Cz5vxZ6inXBNmWLyyG4eprwSq6kHgL4HvY3I2/7Zd1Pcu4JIkfy/J1zOZWfTRXfwc9dD6mjBXXeUBYwi8shuHWdcO+gTw7cC/XQuFWf0p8Fom9xXsAd4P3LS5UJL9wH6A5eXlXVSjrrSxJszQBp/7+vt4ZTcOmUz0maLgZErobaxNFa1pv/H/f/8pwOeZTCf9FeAM4DrgM1X1b7b7vpWVlVpdXZ2lKrWki4PX0Loo5vH7NPn/0teA0s6S3FJVKzuVm+VK4BeY9Om/ctYAWHMW8CwmYwLHgGNJ3gdcDWwbAuqnLg7Ghw/DwYNw7Bg8+WQ7T8ma90Gw6ad+Nf3/4mqfw3fCEEhyFvBi4PnAm4B3VNXHT/Q926mqe5LcAbwuyfqVwKuBT+3m56lbbT+ycP3gth4Ap5wy/y6KNoKu6S6Xvj1KUv2308Dwi4H/APwrJmsEvfkk6/sh4J8DR4C/Ah4Hfu4kf6Y6sH7w2rOn2YPxdlMS1w9u6wFw0UXNHZR3qnOeA6NND6bP6/9FwzX1mEBXHBPor6a7Sk505j2vs/Iu6pw3+/EF8xkTkI7TdH/xiboy5vVIwi7qnDf78TULQ0C9sVP/+DwObl3UKfWJIaDe6OLMe1HP9qWmOCYgSQM07ZiAzxiWpBEzBCRpxAwBSRoxQ0BqkOvva9E4O0hqyDQ3l3kjl/rGEJAastO6PYt6B7KGze4gqQGHD8Odd07W7Nlu3R4f0qI+8kpAOkkbz/D37oXXvhZ+4ieeepbvQ1rUR4aAOjWEPvKNZ/gAy8tb/y7enaw+MgTUmaH0kc9yhu9aROobQ0CdGcoDUDzD1yIzBNSZIfWRe4avRWUIqDOeQUvdMwTUKc+gpW55n4AkjVjrIZDkkiSfTvJQks8meWHbbZAkTbTaHZTknwG/BPwo8AngmW3WL0k6XttjAlcAV1bVx9e+/uuW61fDhnCzlzRmrYVAkj3ACvAHSf4KOB24CXhTVT3cVjvUnKHc7DVEhrOm1eaYwDnAqcArgBcC5wHfCVy2uWCS/UlWk6weOXKkxSY2a+hry7sgWj+th/OBA5PPQ93/1Iw2Q2D9bP/fVdXfVNU9wDuAl2wuWFWHqmqlqlaWlpZabGJzxvCHuH6z13arZk5r6GHZNsNZs2itO6iq7k/yJaDaqrNLQ1kS4USauNnLLqXmDelObM1f2wPD7wN+OslHgceAnwU+1HIbWjGWP8STvdlrDGHZNu/E1izaDoGrgLOB24FHgBuAt7Xchlb4hzidsYRl27wTW9NKVb97Z1ZWVmp1dbXrZmiOnMkiNS/JLVW1slM51w5S5/p81mpAaegMAc1kTAfF7Qatx7QNNHyGgKY2tpk82021HNM20PC5iqimNrb551vdBzG2baDh80pgAXXVHTG2mTzbzfAa0zbQ8Dk7aMF03SVjf7jbQIvB2UED1fXNVX2eydMWt4GGxDGBBdPUej2LxvWFpPnwSmDBjPFO5K67wHbLbiMtAkNgAW3VHTHkA07XXWC7sajBpfExBAZg6AecRZyVtIjBpXEyBAZg6AecRewCW8Tg0jgZAg3psjtmDAecRZuRs4jBpXEyBBrQdXeMB5x+WrTg0jgZAg3oQ3eMBxxJu+F9Ag0Y69x9SYvPK4EG2B0jaVEZAg3Z2B0z5Dn7kobFEGhY14PEkjSLTsYEknxzkkeS/Psu6p8n15uXtEi6Ghh+D/AXHdU9Vw4SS1okrXcHJbkE+D/AnwPf1Hb98+YgsaRF0moIJDkTuBK4EHhNm3W3aTdz9h1MltSFtq8ErgKuraovJtm2UJL9wH6A5eXllprWHQeTJXWltTGBJOcBFwHv3KlsVR2qqpWqWllaWpp/4zrmYLKkrrR5JXABsA+4c+0q4AxgT5Jvq6rvarEdvTOGBeAk9VObIXAI+OCGr3+BSSi8rsU29JKDyZK60loIVNVXgK+sf53kKPBIVR1pqw0na56Dty4AJ6kLnd0xXFUHu6p7Nxy8lTREriI6JQdvJQ2RITAl7wSWNEQuIDclB28lDZEhMAMHbyUNjd1BkjRihoAkjZghIEkjZghI0ogZApI0YoaAJI2YISBJI2YIzOjwYbjmmslnSVp03iw2AxeRkzQ0XgnMwEXkJA2NITADF5GTNDR2B83AReQkDY0hMCMXkZM0JHYHSdKIGQKSNGKthUCS05Jcm+QLSR5M8skkF7dVvyTpqdq8EtgLfBF4EfA1wAHghiT7WmyDJGmD1kKgqh6qqoNV9fmqerKqPgTcAXx3W21ogncMSxqSzmYHJTkHeA5wa1dtmJV3DEsamk4GhpOcCvw28P6qum2L9/cnWU2yeuTIkfYbuA3vGJY0NK2HQJJTgA8AjwKXblWmqg5V1UpVrSwtLbXavhPxjmFJQ9Nqd1CSANcC5wAvqarH2qz/ZHnHsKShaXtM4L3Ac4GLqurhlutuhHcMSxqSNu8TOBf4KeA84O4kR9c+XtlWGyRJx2vtSqCqvgCkrfokSTtz2QhJGjFDQJJGzBCQpBEzBCRpxAwBSRoxQ0CSRswQkKQRMwQkacQMAUkaMUNAkkbMEJCkETMEJGnEDAFJGjFDQJJGzBCQpBEzBCRpxAwBSRoxQ0CSRswQkKQRazUEkpyV5PeTPJTkC0l+vM36JUnHa+1B82veAzwKnAOcB3w4yaeq6taW2yFJosUrgSRPB14OHKiqo1X134A/AP5lW22QJB2vze6g5wBPVNXtG177FPC8FtsgSdqgze6gM4AHNr32APDVmwsm2Q/sX/vyaJLPzLltG50N3NNifYvC7bI1t8tTuU221vZ2OXeaQm2GwFHgzE2vnQk8uLlgVR0CDrXRqM2SrFbVShd195nbZWtul6dym2ytr9ulze6g24G9Sb55w2vfATgoLEkdaS0Equoh4EbgyiRPT/IC4KXAB9pqgyTpeG3fLPZ64GnA3wK/A7yuh9NDO+mGWgBul625XZ7KbbK1Xm6XVFXXbZAkdcRlIyRpxAwBSRqx0YXAtOsXZeKXkty79vH2JGm7vW2ZYbscTPJYkqMbPp7ddnvbkOTSJKtJjiW5foeyP5fk7iQPJLkuyWktNbN1026XJD+Z5IlN+8oF7bW0PUlOS3Lt2t/Og0k+meTiE5Tvzf4yuhDg+PWLXgm8N8lWdy3vB17GZBrr84EfAH6qrUZ2YNrtAvC7VXXGho/PtdbKdt0FXA1cd6JCSV4MvAW4ENgHPBu4Yt6N69BU22XN4U37ys3zbVpn9gJfBF4EfA1wALghyb7NBfu2v4wqBGZcv+jVwK9W1Zeq6q+BXwV+srXGtsh1nbZWVTdW1U3AvTsUfTVwbVXdWlX3A1cx0H0FZtouo1FVD1XVwar6fFU9WVUfAu4AvnuL4r3aX0YVAsy2ftHz1t7bqdwQzLqu0w8muS/JrUleN//m9d5W+8o5Sf5+R+3pk+9Mck+S25McSNL2ysWdSHIOk7+rrabA92p/GVsITL1+0RZlHwDOGOi4wCzb5QbgucAS8FrgrUl+bL7N672t9hXYevuNyZ8B3w58LZMrzR8D3tRpi1qQ5FTgt4H3V9VtWxTp1f4ythCYev2iLcqeCRytYd5YMcu6Tn9ZVXdV1RNV9efAu4BXtNDGPttqX4Gt96vRqKrPVdUda90j/wu4koHvK0lOYbIKwqPApdsU69X+MrYQmGX9olvX3tup3BCczLpOBQzx6mgWW+0rX64q+8yPN+h9Za2X4FomkyteXlWPbVO0V/vLqEJgxvWLfgt4Y5KvT/J1wM8D17fW2BbNsl2SvDTJM9am0H4P8AbgP7Xb4nYk2ZvkdGAPsCfJ6dv0af8W8Jok35bkGcBlDHRfgem3S5KL1/rGSfKtTGbMDHJfWfNeJl2lP1hVD5+gXL/2l6oa1QdwFnAT8BBwJ/Dja6+/kEl3z3q5AG8H7lv7eDtry2wM8WOG7fI7TGaFHAVuA97QddvnuE0OMjl73fhxEFhe+/2XN5R9I/Bl4O+A9wGndd3+rrcL8Ctr2+Qh4HNMuoNO7br9c9om565th0fWtsH6xyv7vr+4dpAkjdiouoMkScczBCRpxAwBSRoxQ0CSRswQkKQRMwQkacQMAUkaMUNAkkbMEJCkETMEJGnEDAFpSkmeluRLSe7c/EzYJL+59jzdS7pqn7QbhoA0pZqsDHk58Czg9euvJ7kGeA3w01X1wY6aJ+2KC8hJM0iyh8njAL+WyQPC/zXwTuDyqrqyy7ZJu2EISDNK8gPAHwIfA/4p8O6qekO3rZJ2xxCQdiHJLcB3AR9k8uwF/5C0kBwTkGaU5EeA89a+fNAA0CLzSkCaQZLvZ9IV9IfAY8APA/+wqj7dacOkXTIEpCkl+V4m4wCfAC4GvgH4NPBHVfWyLtsm7ZbdQdIUkjwX+DBwO/CyqjpWVZ8FrgVemuQFnTZQ2iWvBKQdJFkG/jvwKPCPq+rLG957JvBZ4JNVZRBo4RgCkjRidgdJ0ogZApI0YoaAJI2YISBJI2YISNKIGQKSNGKGgCSNmCEgSSNmCEjSiBkCkjRi/xci8r4chBLqFgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# To plot\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['axes.labelsize'] = 14\n",
    "plt.rcParams['xtick.labelsize'] = 12\n",
    "plt.rcParams['ytick.labelsize'] = 12\n",
    "\n",
    "#plot the data\n",
    "plt.plot(X, y, \"b.\")\n",
    "plt.xlabel(\"$x$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([-0.2, 2.2, 0, 15])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate the Variance of Y\n",
    "\n",
    "Variance is defined as $V[Y] = E[(Y-E[Y])^2]$, where $E$ is the expectation value. \n",
    "\n",
    "It can also be written as $V[Y] = E[Y^2] - E[Y]^2$.\n",
    "\n",
    "In a loose sense, the variance gives an idea of how much uncertainty there is around a random variable. The variance is a function of the probability distribution, but because in a real experiment we wouldn't know what the probability distribution of Y is (or maybe it's not easy to formulate) we can calculate the sample variance instead as an estimate:\n",
    "\n",
    "$$ \\hat{V}[Y] = \\frac{1}{N-1}\\sum_i(y_i-\\bar{y})^2)  \\qquad \\bar{y} = \\frac{1}{N}\\sum_i y_i $$\n",
    "\n",
    "It is $N-1$ because there is one less degree of freedom because $\\bar{y}$ is derived from the sample measurements as well. To ge the real variance of Y requires doing error propagation on the known formula\n",
    "\n",
    "$$ V[mX + b + \\epsilon] = m^2 V[X] + V[\\epsilon] + 2m\\,Cov[X, \\epsilon] $$\n",
    "\n",
    "b is a constant, so has variance 0, and the errors between $\\epsilon$ and $X$ are uncorrelated, and $\\epsilon$ comes from a normal distribution of $N(0,\\sigma^2)$ and $X$ comes a flat uniform distribution from a to b. So \n",
    "\n",
    "$$ V[Y] = m^2 \\frac{1}{12}(b-a)^2 + \\sigma^2 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample variance: [6.40667819] \n",
      "Population variance: 6.833333333333333\n"
     ]
    }
   ],
   "source": [
    "# The real variance of Y\n",
    "y_var = (m**2)*(1./12.)*(flat**2)+sigma2\n",
    "# Our sample variance of Y\n",
    "y_samp_mean = sum(y)/len(y)\n",
    "y_samp_var = np.dot(np.transpose(y-y_samp_mean),y-y_samp_mean)/(len(y)-1.)\n",
    "\n",
    "print('Sample variance: {} \\nPopulation variance: {}'.format(y_samp_var[0],y_var))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use Normal equation to get the best fit line\n",
    "\n",
    "The Population Mean Square Error is the average square difference between the random variable and the regression function. For any regression function f(X) it is \n",
    "\n",
    "$$ MSE[f(X)] = E[(Y-f(X))^2] = E[ E[(Y-f(x))^2|X=x] ] $$\n",
    "\n",
    "The population MSE is closely related to the variance of Y as defined above. Instead of calculating the average of the fluctuations around $E[Y]$ we are instead calculating the averge fluctuations around each point $f(X)$ for each X (and then averaging those together weighted by $X$'s probability distribution). There is an optimal $f(X)$ to minimize the MSE, and that is when $f(x) = E[Y|X=x]$.\n",
    "\n",
    "We can find the optimal $f(X)$ which gives this minimum MSE by using the Normal formula. The Normal formula gives consistent estimators for the parameters of the regression function $f(X)$ when $f(X)$ is a line. It essentially gives us the minimum of the **Sample** Mean Square Error (not the minimum of the population MSE). The Sample Mean Square Error is defined as:\n",
    "\n",
    "$$ \\hat{MSE} = \\frac{1}{n} \\sum_i{(y_i - b_{pred} - m_{pred}x_i )^2} = \\frac{1}{n}*(y-X\\hat{\\theta})^T(y-X\\hat{\\theta})$$\n",
    "\n",
    "And the Normal Formula for getting $b_{pred}$ and $m_{pred}$ is\n",
    "\n",
    "$$ \\hat{\\theta} = [b_{pred},m_{pred}] = \\left( X^T * X \\right)^{-1}*X^T*y$$\n",
    "\n",
    "X is the \"design matrix\" which encodes our data points for the $N$ events in our sample. You have to add a column of 1's to design matrix so that you can dot it with the $\\theta$ vector of parameters, i.e. $[1 \\; x_i] \\cdot \\begin{bmatrix} b_p \\\\ m_p \\end{bmatrix} = b_p + x_i m_p$.  The normal equation is a linear equation where $ L = \\left( X^T * X \\right)^{-1}*X^T $ is the linear operator.\n",
    "\n",
    "$$ [b_{pred},m_{pred}] = L*y $$\n",
    "\n",
    "If you wanted to recover the exact b and m, in the high N limit, \n",
    "\n",
    "$$ \\left( X^T * X \\right)^{-1}*X^T*y \\rightarrow \\frac{Cov[X,Y]}{V[X]} = \\begin{bmatrix}b \\\\m \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "population intercept: 3, population slope: 4\n",
      "best fit intercept: 3.626, best fit slope: 3.723\n",
      "gradient intercept: 3.626, gradient slope: 3.723\n"
     ]
    }
   ],
   "source": [
    "#Need to add a column of constant 1 to fit the intercept\n",
    "#Normal equation says inv(X^T*X)*X^T*y gives the best fit m and b.\n",
    "X_b = np.c_[np.ones((n, 1)), X]  # add x0 = 1 to each instance\n",
    "theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)\n",
    "print('population intercept: {}, population slope: {}'.format(b,m))\n",
    "print('best fit intercept: {:.3f}, best fit slope: {:.3f}'.format(theta_best[0][0],theta_best[1][0]))\n",
    "\n",
    "#Rather than using the normal equation, let's also use gradient descent\n",
    "from sklearn.linear_model import LinearRegression\n",
    "lin_reg = LinearRegression()\n",
    "lin_reg.fit(X, y)\n",
    "print('gradient intercept: {:.3f}, gradient slope: {:.3f}'.format(lin_reg.intercept_[0], lin_reg.coef_[0][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the best fit line\n",
    "\n",
    "Now that the Normal equation has calculated the $b_{pred}$ and $m_{pred}$ that minimizes the sample MSE, we can plot the regression line over the points in our sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Draw a line between two points, x=0 and x=2\n",
    "X_new = np.array([[0], [2]])\n",
    "X_new_b = np.c_[np.ones((2, 1)), X_new]  # add x0 = 1 to each instance\n",
    "y_predict = X_new_b.dot(theta_best)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEVCAYAAAAGrllxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH/5JREFUeJzt3XmYFNW9//H3lwERRXABwW0g7sZdJyIuiKIiSFxRUVExE8k1MXv8aZ64oHL1ahavicbIYwsICi5RrhoxXlTUxFEz6HUhccctrCpBdpiZ7++PM4QRema6Z6qrurs+r+fpZ5zu6j6Hsqc+VadOfcvcHRERSacOSXdARESSoxAQEUkxhYCISIopBEREUkwhICKSYgoBEZEUUwiIiKSYQkBEJMUUAiIiKdYx6Q60pkePHt63b9+kuyEiUlJmzZr1mbv3bG25og+Bvn37Ultbm3Q3RERKipl9lMtyGg4SEUkxhYCISIopBEREUkwhICKSYgoBEZEUUwiIiKSYQkBEJMUUAiIiKaYQEBFJsUhDwMwuMbNaM1ttZhOaWeZqM3MzOzbKtkVEJH9Rl42YC4wFBgNdNnzRzHYBhgPzIm5XRETaINIjAXd/yN2nAZ83s8itwGXAmijbFRGRtontnICZnQGscffH42pTRERaFksVUTPrClwPHJ/j8qOB0QCVlZUF7JmISLrFdSRwDTDJ3efksrC7j3P3Knev6tmz1XLYIiLSRnGFwCDgB2Y238zmAzsB95vZZTG1LyIiWUQ6HGRmHRs/swKoMLNNgTpCCHRqsujfgJ8A06NsX0RE8hP1kcAVwErgcmBk439f4e6fu/v8dQ+gHljs7ssibl9ERPIQ6ZGAu48BxuSwXN8o2xURkbZR2QgRkRRTCIiIpJhCQEQkxRQCIiIpphAQEUkxhYCISIopBEREUkwhICKSYgoBEZEUUwiIiKSYQkBEJMUUAiIiKaYQEBFJMYWAiEiKKQRERFJMISAikmIKARGRFFMIiIikmEJARCTFIg0BM7vEzGrNbLWZTWjy/KFm9r9m9oWZLTKzB8xsuyjbFhGR/EV9JDAXGAvctcHzWwHjgL5AH2ApMD7itkVEJE8do/wwd38IwMyqgB2bPD+96XJmdivwbJRti4hI/pI6JzAAmJ1Q2yIi0ijSI4FcmNl+wFXAyS0sMxoYDVBZWRlTz0RE0ifWIwEz2xWYDvzQ3Z9vbjl3H+fuVe5e1bNnz/g6KCKSMrGFgJn1AWYA17n7pLjaFRGR5kU6HGRmHRs/swKoMLNNgTqgF/A0cJu7/yHKNkVEpO2iPidwBXB1k99HAtcADuwMXG1m/37d3btG3L6IiOQh6imiY4Axzbx8TZRtiYgUWk0NzJwJAwdC//5J96YwYp8dJCJSCmpqYNAgWLMGNtkEnnqqPINAtYNERLKYOTMEQH19+DlzZtI9KgyFgIgUpZoauOGG8DMJAweGI4CKivBz4MBk+lFoGg4SkaJTDEMx/fuHdnVOQEQkZtmGYpLYCPfvX74b/3U0HCQiRacQQzFJDy8VKx0JiEjRiXooJtHhpbo6ePxxWLUKzjwzpkZzpxAQkchFMb8+yqGYQg4vNftvffddyGRg4kSYPx/69VMIiEj5K4aTuhtaN7y0rk9RzfTZ8N/6zGPL6ffJg2Hj//zzYTxr6FCorg4/i5BCQEQiVQwndTfcOy/UTJ+ZM2HNaueghr/x7VUZDhg6BVYvhV13DScgzj8ftt8+msYKRCEgIpEq1F53rpo7Eol8ps/nnzNiwWSGeYZ9eYMV3oUlx5zBtpdXw5FHglmEjRWOQkBEIpX0/PpCHIn8+8hiQAP9l88Iwz3TpvG1NWtYttc3mL77H9jmeyM45LjuEfwL4qUQEJHIJTm/vq1HIs2d4K2pgVFHf8S5a8azA+PBP4att4aLL4bqarruuy9b1sBTM6G+a/LnP/KlEBCRstKWI5GsQ0gHrYZp0+h9RYZ/rJ4BwAyOY9aIX3LqhJOhc+fm31tCQaAQEJGyk++RSNMhpD1Wv07Hn2bg7cnwxRds36uSGzpeRabhQuZ37sNTPwA6Z39vklc3t5VCQERSb1DVEuZ1mMJ59Rm+0VBLQ+0mcOopUF1N50GDOOblCjrMzH5kkfSJ8PYyd0+6Dy2qqqry2trapLshIuXGHZ57LpzkffBBWLmShb33ZflZ1XztypGwzTY5f1Qx3nzGzGa5e1Vry+lIQETSZe7ccBXvXXfBe+9Bt25hPn91NdtWVbVpamcpF5pTCIhIUYtkL3vtWvjTn8Je//TpYQB/wAC48koYPhw22yzCHpeWSEPAzC4BRgH7AlPcfVST1wYBtwGVwEvAKHf/KMr2RaS8tHvmzdtvhw3/3XfDggWw3XZw6aXwrW/BbrsVrN+lJOojgbnAWGAw0GXdk2bWA3gI+DbwKHAdcB9waMTti0gZadPMm+XL4f77w8b/r38N9XuGDQv1e4YMgY4aAGkq0rXh7g8BmFkVsGOTl04DZrv7A42vjwE+M7M93f2tKPsgIuUj55k37vDSS2HDP3UqLFsGu+8ON94Yxvt7946kP8V4Ari94orEvYHX1v3i7svN7P3G5xUCIpJVqxd+LVoEkyaFk7yzZ4ex/TPPDHv9hx8eaf2eUr8orDlxhUBXYNEGzy0Btsi2sJmNBkYDVFZWFrZnIlLUNpp5U18PTz4Z9vofeSSc9O3XD8aNg7POCrN9CqDULwprTlwhsAzY8P9MN2BptoXdfRwwDsJ1AoXtmoiUhDlzYPz48Pj0U+jRAy65JJzk3Wefgjdf6heFNSeuEJgNXLDuFzPbHNil8XkRkexWrYKHHw57/U89FYZ3Bg+Gm2+Gk04KW+OYJF0dtVCiniLasfEzK4AKM9sUqAMeBn5pZqcDfwKuAl7XSWERyer//i9s+O+5BxYvhr594dprYdQo2GmnxLpVyheFNSfqI4ErgKub/D4SuMbdxzQGwK3AZMJ1AiMibltEStm//gX33hs2/q+8Eqp0nnpqOMl7zDHQoUPSPSxLUU8RHQOMaea1GcCeUbYnIiWuoQGefTZs+P/4xzD8s//+8Nvfwrnnhrr9UlC6akJE4vfPf8KECWFq5wcfQPfucOGFYa//oINK5taM5UAhICLxWLMGHnss7PU/8UQ4Chg4EK65Bk4/Hbp0afUjJHoKAREprH/8Y339nkWLYPvt4fLLw9TOXXZJuneppxAQkegtWwb33Rc2/jU1oV7PN78ZhnsGD25z/Z5yLNuQNIWAiETDHV58Ee68MwTA8uWw557wy1/CeedBr17t+vhyLduQNIWASAoUdA964cJQvyeTCUM/m28eyjdUV4fGIjrJW65lG5KmEBApcwXZg66vhz//eX39nrq68KF33hkKuG2RtSzYV/qUbyiVa9mGpCkERMpcpHvQH3wQpnVOmBCmefbsCT/8YTjJ+/Wv5/QRbQ2lci3bkDSFgEiZa/ce9MqV8NBDYa//mWfClbsnnBAu6Bo2LO/6Pe0JpXIs25A0hYBImWvzHvQrr4QN/733hpIOO+8MY8fCBRfAjju2/v5maFinuCgERFIg5z3oxYtD0bZMJhRx23TTcCFXdTUcdVQk9Xs0rFNcFAIiadfQEIZ5Mpkw7LN6NRx4INx2G5x9Nmy1VeRNalineCgERNLqk0/CCd7x48MNW7bcEr797bDXf+CBSfdOYqIQEEmTNWvClM5MJkzxdA9lmseODWWbVb8ndXIa4DOzP5iZm9n2WV7bw8zWmNkt0XdPRCIxezb85Cewww5wxhnw5pvwi1/A+++HAfpzzlEApFSuRwI1wHeAQ4BpG7x2M/AlzdxHQEQSsnQpTJ0a9vpfegk6dQq3ZKyuhuOPh4qKpHsoRSDXEHix8edXQsDMTgSGAN9z98UR901E8uUOf/1r2PDffz+sWBEu4vr1r0P9np49k+6hFJmcQsDd3zazLwghAICZdQJ+A7wJ3FGY7olIThYsgIkTw9W8b78NXbuGIZ7qaujXTzdpkWblc2L4ReBwMzN3d+CHwO7Ase5eX5DeSclRqd8Y1dWFm7NkMuFmLXV1cPjhcNllYdy/a9ekeyglIN8QGArs0XhUcCUwzd2fyvUDzKwv8HugP7AaeBD4kbvX5dEPKVIq9RuT995bX79n3jzYdlv48Y9D/Z49dRtvyU8+IVDT+PMQYADQGfhpnu39HlgIbAdsCfwv8F3gt3l+jhQhlfotoBUrwo3YM5lwY/YOHWDo0DDcc+KJ4aRvAejIrvzlEwIvAQ1ANXAE8Et3/yDP9r4G3Oruq4D5ZvYEsHeenyFFSjVhIuYOs2atr9/z5ZfhdozXXx/q92y/0YztSOnILh1yDgF3X2pmfyccBcwH/rMN7d0CjDCzmcBWhJlFV7bhc6QIqSZMRD7/fH39ntdfD/P3hw8Pe/0DBsR2kldHdumQ7xXDLwP7AD9396VtaO9Z4CLCdQUVwEQ2vu4AMxsNjAaorKxsQzOSlDhqwpTbEEVNDcx8uoGTuz7F12sy8PDDYatbVQW33x7q93TvHnu/dGSXDhYm+uSwYJgS+hbwBXCI5/rG9e/vAHxImE76K6ArcBfwtrv/v+beV1VV5bW1tfk0JTFJYmNcbkMUsx7+mOlnjmdk3Xj68hF1W2xFx1Ejw17//vu36TOj/P9SboGbJmY2y92rWl3Q3XN6AD8nnBM4NNf3bPD+HoAD3Zs8dwrwZkvvO/jgg12KzwsvuHfp4l5REX6+8EI8bR5/vHuHDu4Q2r7++sK3ef31Ef/7Vq1yv+8+9+OP9wYzd/AnOdbPtil+4zUr2/XRSfx/keIE1HoO2+YWh4PMbGtgMLAfcCnwG3d/saX3tBA2n5nZHOBiM1t3JHAB8FpbPk+SFfd48bojgNWrQ+XjDh0KP0QR+VHHG2+Ecf7Jk8O4/0478emFV3LcPRfyXl1fNtkEvn9c+/qscXzJV2sF5AYD9wLfItQIuqyd7Z0GnAAsAt4D6oAft/MzJQHrxosrKqLdGNfUwA03hJ9Nrdu4rQuAY4+NbiiotTabblDz9uWXcMcdcMghsN9+8Pvfh6qdTzwBc+awU+Yaxj/Tl+uui+bfU6j/L1LGcjlcSPKh4aDiFfVQSUtDGYUa5ihImw0N7s8+637++eGN4L7PPu433+y+aFE0HW9BQYawpOQQxXCQSEuingnU0lBGoaafRtrmvHnr6/e8+y5ssUUo2lZdDd/4RmxTO3XXLsmHQkCKRmtTEguxcWt3m3V18PjjcOed4Wd9PRx5ZKjVP3w4bL55tB0WiZhCQIpGEhebtbnNd94Je/wTJ8L8+dC7N/zsZ6F+z+67F7DHItHK+TqBpOg6ASkay5fDgw+GGT7PPx/Ovp54YhjuGTKkYPV7RNoi1+sEdCQg0hJ3+NvfwoZ/ypRwt67ddoP/+i84/3zYbrukeyjSLgoBkWw++yzM589kwv14u3SBM88Me/1HHKGbtEjZUAiIrFNfDzNmhA3/tGmwdm2Y33/HHTBiBHTrlnQPRSKnEBD58EMYPz48PvkEttkGvve9cJJ3333z+ijV2pFSoxCQdFq1KuztZzJhehDAccfBr34FJ58MnTvn/ZG5lJlQSEixUQhIurz22vr6PYsXQ58+cPXVMGpU+O92aK1uT7lVQJXyoBCQ8vevf4WZPZlMuFPXJpvAqaeGk7yDBoViRO1UUwMffxxmjUL2C89U3E2KkUJAypN7uBdvJhPm9q9aFQq43XILnHtuGPePSNM9/I4d4aKLwuzRDTfwukmLFCOFgCQq8jHyuXNhwoRwNe/774cZPaNGhb3+gw8uyNTOpnv4AJWV2f8tuv2mFCOFgCQmsjHytWvhscfCXv/06aHe9FFHhbH+00+HzTaLvO9N5bOHr+JuUmwUApKYdo+Rv/VW2PDffTcsXBiu3r3ssjC1c9ddC9TrjWkPX0qZQkAS06Yx8mXL4IEHQtXOF14Ig/DDhoXhnhNOCL8nQHv4UqoUApKYnPeg3eGll8Je/9SpIQj22ANuuimcge3VK8Zei5QXhYAkqsU96EWLYNKksPH/+9/D2P5ZZ4W9/sMOU/0ekQgoBKS41NfDk0+GDf8jj4STvv36wbhxIQBUv0ckUrGHgJmNAK4GKoH5wCh3fz7ufkiRmTMnTOucMAE+/RR69IBLLgl7/XvvnXTvRMpWrCFgZscBNwJnAS8DKsaeZqtWwUMPhb3+p58OwzuDB8PNN8NJJ4WzxSJSUHEfCVwDXOvuLzb+/s+Y25eItelir1dfDRv+e+4JJR369oVrrw0Xde20U8H6KiIbiy0EzKwCqAIeMbP3gE2BacCl7r4yrn5IdPK62GvxYrj33rDxf/XVUKXztNPCcM/RR0dSv0fWU7VSyVWcRwK9gE7AcOBIYC3wP8AVwC+aLmhmo4HRAJWVlTF2MVrl/ofY6sVeDQ3hyUwmDPusWgUHHAC/+x2ccw5svXUyHS9zqlYq+YgzBNbt7f/O3ecBmNlvyBIC7j4OGAfhRvMx9jEyafhDbPZir08/DSd4x4+HDz6A7t3DVbzV1XDQQRt9TrmHZdxUrVTyEVsIuPtiM/sUKMmNer7S8IfY9GKvow9fw6FzH4WhGfjzn8NRwNFHh7H+004L9+jNIg1hGTdVK5V8xH1ieDzwfTN7gjAc9CPgsZj7EIu0/CH27/53+n+WgeGTwsVdO+wAP/85XHgh7LJLq+9PQ1jGTbWMJB9xh8B1QA/gHWAVcD/wnzH3IRZl/Ye4dCncd18Y63/xxVCv56STwnDP4MHr76ySg7SEZdxUy0hyZe7FPTpTVVXltbW1SXdD3MPYTSYTAmD5cthrr7DhP+882HbbNn+0zgmIRM/MZrl7VWvLqWyEtGzhwlCqOZMJpZs33xxGjAgb/0MPjaR+TzHvtSqgpNwpBGRjdXXh5G4mA48+Gn4/7DDIZHipz5k8/XJXBgL9y7x+W3MnrRUMUk4UArLe+++vr98zdy707Ak/+lGY3rnXXqmbyZPtpDWkax1I+VMIpN3KlfDHP4a9/pkzw5W7J5wQLugaNuwr9XvSNpMn20nrtK0DKX8KgRLU7uEId3jllbDhv/deWLIEdt4Zxo4N9Xt22CHr29I2k6e5GV5pWgdS/jQ7qMS0a0jmiy9C0bZMBl57DTbdNNyIvbo63Jg9h/o9Gg/XOpDSoNlBZSrv4YiGhlCmOZOBhx+G1atD6Ybbbgv1e7bcMq/2i3kmT1y0DqScKARKTM5DMp98Emr3jB8PH34IW20FF10U9voPOCC+DkdEe98ihaEQKDEtXom8enW4JWMmE27R6B7Gjq6/Hk49NQz/lKBSnZWk4JJSoBAoQRsNR7z5JnPHZtjy0clstuIz2HFHuOKKUL/na19LrJ9RKcUZOaUaXJI+CoFS9eWXMHVq2Ot/+WV60IlHOJlJm1Rz+ZTj6H9E7vV7il0pzkoqxeCSdFIIlBJ3+Mtfwob/gQdgxQrYe29mDP0NI58YyYKGnlTUw6HPQ/8jku5sdEqxGF8pBpekk0IgIgUd/50/HyZODFfzvvMOdO0K554bTvIecgibv2h8+QxUlPEGp9Rm5JRicEk6KQQiUJDx37o6mD497PU/9lgYVzjiiFCr/4wzQiG3RtrgFKdSCy5JJ4VABCId/3333bDHP3EizJsHvXrBT38a6vfssUezb9MGR0TaQiEQgXaP/65YAQ8+GPb6n3suXLl74olhuGfoUOjUqQC9FhFRCESiTcMx7lBbGzb8U6aE2T677hrm9F9wAWy/fYF7LSKiEIhM0+GYFk8Sf/45TJ4cNv5vvBFuwD58eNjrHzAgkpu0iIjkSiEQsawnifs1wIwZYcM/bVp4saoKbr8dzj4bundPutsiklKJhICZ7Qa8ATzo7iOT6EOhND1J3Hv1R9RdMR7eGw8ffwxbbw3/8R9hr3+//ZLuqohIYkcCtwF/S6jtgjr6sNWc3WEa59dnGNQwA3sGOPZYuOkmOOUU6Nw56S6KiPxb7CFgZiOAfwEvALvG3X7BvP46ZDIcOnkyh679giVbVvLP065ip6suhD59ku6diEhWsYaAmXUDrgUGAdVxtl0QS5aEmT2ZTJjps8kmYW+/uprugwbRvSL3+j2qOCkiSYj7SOA6IOPun1gLs2DMbDQwGqCysjKmruXIHZ5/fn39npUrYd994b//G0aOhG22yfsjVXFSRJISWwiY2QHAscCBrS3r7uOAcRBuL1ngruVm3rz19XvefRe6dYPzzw8neauq2jW1UxUnRSQpcR4JDAT6Ah83HgV0BSrM7OvuflCM/cjd2rXw+ONhr//xx8NWesCAUKt/+HDYbLNImlHFSRFJSpwhMA6Y2uT3nxFC4eIY+5Cbt99eX79nwQLo3RsuvTTU79ltt8ibUwE4EUlKbCHg7iuAFet+N7NlwCp3XxRXH1o1eTLccUeo2V9RAcOGheGeIUOgY8dw8vbBwmyoVQBORJKQ2BXD7j4mqbabdd99sHAh3HhjGO/v3fvfL+nkrYiUI5WNaGrSpFDCIctJXp28FZFy1CHpDhSVLbdsdpbPupO3FRU6eSsi5UNHAjnSyVsRKUcKgTzo5K2IlBsNB4mIpJhCQEQkxRQCIiIpphAQEUkxhYCISIopBEREUkwhICKSYgqBPNXUwA03hJ8iIqVOF4vlQUXkRKTc6EggD9mKyImIlDKFQB5URE5Eyo2Gg/KgInIiUm4UAnlSETkRKScaDhIRSTGFgIhIisUWAmbW2cwyZvaRmS01s1fNbEhc7YuIyMbiPBLoCHwCHAV0B64E7jezvjH2QUREmogtBNx9ubuPcfcP3b3B3R8D5gAHx9WHKOiKYREpJ4nNDjKzXsDuwOyk+pAvXTEsIuUmkRPDZtYJuAeY6O5vZXl9tJnVmlntokWL4u9gM3TFsIiUm9hDwMw6AJOANcAl2ZZx93HuXuXuVT179oy1fy3RFcMiUm5iHQ4yMwMyQC9gqLuvjbP99tIVwyJSbuI+J3A7sBdwrLuvjLntSOiKYREpJ3FeJ9AH+A5wADDfzJY1Ps6Nqw8iIvJVsR0JuPtHgMXVnoiItE5lI0REUkwhICKSYgoBEZEUUwiIiKSYQkBEJMUUAiIiKaYQEBFJMYWAiEiKKQRERFJMISAikmIKARGRFFMIiIikmEJARCTFFAIiIimmEBARSTGFgIhIiikERERSTCEgIpJiCgERkRSLNQTMbGsze9jMlpvZR2Z2Tpzti4jIV8V2o/lGtwFrgF7AAcCfzOw1d58dcz9ERIQYjwTMbHPgdOBKd1/m7n8BHgHOi6sPIiLyVXEOB+0O1Lv7O02eew3YO8Y+iIhIE3EOB3UFlmzw3BJgiw0XNLPRwOjGX5eZ2dsF7ltTPYDPYmyvVGi9ZKf1sjGtk+ziXi99clkozhBYBnTb4LluwNINF3T3ccC4ODq1ITOrdfeqJNouZlov2Wm9bEzrJLtiXS9xDge9A3Q0s92aPLc/oJPCIiIJiS0E3H058BBwrZltbmaHAycDk+Lqg4iIfFXcF4t9F+gCLASmABcX4fTQRIahSoDWS3ZaLxvTOsmuKNeLuXvSfRARkYSobISISIopBEREUix1IZBr/SILbjSzzxsfN5mZxd3fuOSxXsaY2VozW9bksXPc/Y2DmV1iZrVmttrMJrSy7I/NbL6ZLTGzu8ysc0zdjF2u68XMRplZ/QbflYHx9TQ+ZtbZzDKNfztLzexVMxvSwvJF831JXQjw1fpF5wK3m1m2q5ZHA6cQprHuBwwDvhNXJxOQ63oBuM/duzZ5fBBbL+M1FxgL3NXSQmY2GLgcGAT0BXYGril05xKU03ppVLPBd2VmYbuWmI7AJ8BRQHfgSuB+M+u74YLF9n1JVQjkWb/oAuDX7v6pu/8T+DUwKrbOxkh1nbJz94fcfRrweSuLXgBk3H22uy8GrqNMvyuQ13pJDXdf7u5j3P1Dd29w98eAOcDBWRYvqu9LqkKA/OoX7d34WmvLlYN86zp908y+MLPZZnZx4btX9LJ9V3qZ2TYJ9aeYHGhmn5nZO2Z2pZnFXbk4EWbWi/B3lW0KfFF9X9IWAjnXL8qy7BKga5meF8hnvdwP7AX0BC4CrjKzswvbvaKX7bsC2ddfmjwH7ANsSzjSPBu4NNEexcDMOgH3ABPd/a0sixTV9yVtIZBz/aIsy3YDlnl5XliRT12nv7v7XHevd/cXgFuA4TH0sZhl+65A9u9Varj7B+4+p3F45A3gWsr8u2JmHQhVENYAlzSzWFF9X9IWAvnUL5rd+Fpry5WD9tR1cqAcj47yke27ssDdNWb+VWX9XWkcJcgQJlec7u5rm1m0qL4vqQqBPOsX3Q38xMx2MLPtgZ8CE2LrbIzyWS9mdrKZbdU4hfYQ4AfA/8Tb43iYWUcz2xSoACrMbNNmxrTvBqrN7OtmthVwBWX6XYHc14uZDWkcG8fM9iTMmCnL70qj2wlDpd9095UtLFdc3xd3T9UD2BqYBiwHPgbOaXz+SMJwz7rlDLgJ+KLxcRONZTbK8ZHHeplCmBWyDHgL+EHSfS/gOhlD2Htt+hgDVDb++yubLPsTYAHwJTAe6Jx0/5NeL8CvGtfJcuADwnBQp6T7X6B10qdxPaxqXAfrHucW+/dFtYNERFIsVcNBIiLyVQoBEZEUUwiIiKSYQkBEJMUUAiIiKaYQEBFJMYWAiEiKKQRERFJMISAikmIKARGRFFMIiOTIzLqY2adm9vGG94Q1szsb76c7Iqn+ibSFQkAkRx4qQ14N7AR8d93zZnYDUA18392nJtQ9kTZRATmRPJhZBeF2gNsSbhD+beBm4Gp3vzbJvom0hUJAJE9mNgx4FHgKOAa41d1/kGyvRNpGISDSBmY2CzgImEq494L+kKQk6ZyASJ7M7EzggMZflyoApJTpSEAkD2Z2PGEo6FFgLXAGsK+7/yPRjom0kUJAJEdm1o9wHuBlYAiwI/AP4HF3PyXJvom0lYaDRHJgZnsBfwLeAU5x99Xu/j6QAU42s8MT7aBIG+lIQKQVZlYJ/BVYAxzm7guavLYd8D7wqrsrCKTkKARERFJMw0EiIimmEBARSTGFgIhIiikERERSTCEgIpJiCgERkRRTCIiIpJhCQEQkxRQCIiIpphAQEUmx/w+aKBPtX5M7lQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#now plot the line\n",
    "plt.plot(X, y, \"b.\")\n",
    "plt.plot(X_new, y_predict, \"r-\")\n",
    "plt.xlabel(\"$x$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([-0.2, 2.2, 0, 15])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate the Confidence interval on the best-fit parameter\n",
    "\n",
    "So we have our regression estimators, $m_{pred}$ and $b_{pred}$.\n",
    "\n",
    "Each estimator is a random variable. They were calculated from our N random points ${(x_i,y_i)}$, and so would change if we collected a new sample of points. That means these estimators have their own probability distributions associated with them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distributions of $m_p$ and $b_p$\n",
    "\n",
    "We want to know the expectation value and variance of $m_{p}$ and $b_{p}$, our best fit values. The estimates of these two parameters are often correlated, which means that we model their probability distribution as a 2D distribution. By the central limit theorem, if our sample size N is big enough, the shapes of the distributions should be a multivariate Gaussian of 2 dimensions, the means of each parameter should approach their real parameters, and their variances and covariances should also approach $\\frac{\\Sigma}{N}$ (where $\\Sigma$ is for the full covariance matrix).\n",
    "\n",
    "The expectation values will be the real parameters \n",
    "\n",
    "$E(m_p) = m $ and $E(b_p) = b $\n",
    "\n",
    "The variance can be gotten from the normal equation using standard error propagation.\n",
    "\n",
    "$$ V[\\begin{bmatrix} b_p \\\\m_p \\end{bmatrix}] = V[\\left( X^T * X \\right)^{-1}*X^T*Y] = V[L*Y] $$\n",
    "\n",
    "We could try to do the exact propagation of errors using the math fact that $V[AB] = Cov[A^2,B^2] + E[A^2]E[B^2] - (E[AB])^2$, and the fact that Y can be written as $Y = mX + b + \\epsilon$... but we are going to skip all of that derivation here and also assume that $\\epsilon$ is uncorrelated with $X$ and that  $V[\\epsilon]$ is a constant $\\sigma^2$.\n",
    "\n",
    "$$ V[L*Y] = \\sigma^2*L*L^T = (X^T*\\frac{1}{\\sigma^2}*X)^{-1}  = \\sigma^2(X^T*X)^{-1} $$ \n",
    "\n",
    "This is also know as the inverse of the Hessian, where $H = X^T*\\frac{1}{\\sigma^2}*X$. \n",
    "\n",
    "So now let's calculate the sample Variance for the estimators. If $X = \\begin{bmatrix}1 & x_1\\\\1 & x_2 \\\\ .. & ..\\end{bmatrix} $, then $X^T*X = \\begin{bmatrix} n & \\sum_i{x_i}\\\\ \\sum_i{x_i} & \\sum_i{x_i^2} \\end{bmatrix}$ \n",
    "\n",
    "and using the math fact that $ \\begin{bmatrix}a & b \\\\ c & d \\end{bmatrix}^{-1} = \\frac{1}{ad-bc} \\begin{bmatrix}d & -b \\\\ -c & a  \\end{bmatrix} $\n",
    "\n",
    "$ \\hat{V}[\\begin{bmatrix} b_p \\\\m_p \\end{bmatrix}] = \\sigma^2 (X^T*X)^{-1} = \\frac{\\sigma^2}{n \\hat{V}(x)} \\begin{bmatrix}\\bar{x^2} & -\\bar{x}\\\\ -\\bar{x} & 1 \\end{bmatrix} $, where we have used the sample Variance of x $\\hat{V}(x) = \\bar{x^2} - \\bar{x}^2$\n",
    "\n",
    "So finally, \n",
    "\n",
    "$\\hat{\\sigma}_{b_p} = \\sqrt{\\hat{V}_{00}} = \\sigma\\sqrt{\\frac{\\bar{x^2}}{n*\\hat{V}[x]}}$ \n",
    "\n",
    "and \n",
    "\n",
    "$\\hat{\\sigma}_{m_p} = \\sqrt{\\hat{V}_{11}} = \\sigma\\frac{1}{\\sqrt{n*\\hat{V}[x]}}$\n",
    "\n",
    "We are not done just yet since we don't know what $\\sigma^2$ is exactly, so we have to come up with an estimator for that too. In many toy problems, $\\sigma^2$ is just a given, but here we will calculate it.\n",
    "\n",
    "The estimator for $\\sigma^2$ is actually going to be very close to our sample MSE($f(X)$), but not exactly that. When we plug in the $b_p$ and $m_p$ predictions into the regression function, since they were both calculated from the n measurements in our sample, we have to reduce the degrees of freedom by two for an unbiased measurement. \n",
    "\n",
    "$$ \\hat{MSE} = \\frac{1}{n-2} \\sum_i{(y_i - b_{pred} - m_{pred}x_i )^2} $$\n",
    "\n",
    "\n",
    "So our estimate $\\hat{\\sigma}^2 = \\hat{MSE}_{n-2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Mean Square Error, and  get estimate of $\\sigma^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample MSE: 1.313\n",
      "sample MSE w/ N-2 (predicted sigma^2): 1.374\n",
      "MSE from sklearn: 1.313\n",
      "\n",
      "population sigma^2: 1.5\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "resids= X_b.dot(theta_best) - y\n",
    "sample_MSE = ((1.0/n) * sum(np.power(resids,2)))\n",
    "print('sample MSE: {:.3f}'.format(sample_MSE[0]))\n",
    "# To get an estimator for sigma^2, we should take the MSE but use 1/(m-2)\n",
    "# n-2 comes from the fact that since there are two parameters being esimated, m and b, \n",
    "# there are two linear constraints that reduce the degrees of freedom\n",
    "sample_MSE_n2 = (1.0/(n-2.0)) * sum(np.power(resids,2))\n",
    "print('sample MSE w/ N-2 (predicted sigma^2): {:.3f}'.format(sample_MSE_n2[0]))\n",
    "\n",
    "#get the score from gradient descent\n",
    "from sklearn.metrics import mean_squared_error\n",
    "reg_score = mean_squared_error(X_b.dot(theta_best), y)\n",
    "print('MSE from sklearn: {:.3f}'.format(reg_score))\n",
    "print('')\n",
    "print('population sigma^2: {}'.format(sigma2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll use $\\hat{\\sigma}^2$ to calculate the standard error of $b_p$ and $m_p$, by just subsituting in $\\hat{\\sigma}^2$ for $\\sigma^2$ in our variance matrix.\n",
    "\n",
    "With total best-fit V being\n",
    "\n",
    "$\\hat{\\sigma}_{b_p} = \\sqrt{\\hat{V_{00}}} = \\hat{\\sigma}\\sqrt{\\frac{\\bar{x^2}}{n*\\hat{V}[x]}}$  \n",
    "\n",
    "and \n",
    "\n",
    "$\\hat{\\sigma}_{m_p} = \\sqrt{\\hat{V_{11}}} = \\hat{\\sigma}\\frac{1}{\\sqrt{n*\\hat{V}[x]}}$\n",
    "\n",
    "\n",
    "\n",
    "And we'll use this to create our confidence intervals.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2 = sample_MSE_n2\n",
    "s = math.sqrt(sample_MSE_n2)\n",
    "#calculate the estimate for the Variance of our parameters b and m\n",
    "V_best = np.linalg.inv(X_b.T.dot(X_b))*s2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating the confidence Intervals\n",
    "\n",
    "We need to first set our significance level. We know the mean and variance of our best fit parameters, now we just need to know what distribution (e.g. Gaussian, Student's t, etc) to apply these to.\n",
    "\n",
    "When the sample size is large enough, by the central limit theorem, one can use a standard Gaussian distribution. You can get the Z-score from the confidence level by using the percentile function of the gaussian. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set the confidence level at 68%, and 95%\n",
    "alpha1 = 0.32\n",
    "alpha2 = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For 68% level, Z-score is: 0.99\n",
      "For 95% level, Z-score is: 1.96\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "# The Z-score for the given confidence level, taken from a Gaussian\n",
    "sig1 = round(norm.ppf(1.0-(alpha1/2.0)),2) #percentile function\n",
    "sig2 = round(norm.ppf(1.0-(alpha2/2.0)),2) #percentile function\n",
    "print('For 68% level, Z-score is:',sig1)\n",
    "print('For 95% level, Z-score is:',sig2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the 95% confidence interval for a Gaussian with a known and $\\sigma^2$ would be\n",
    "\n",
    "$$\\left( m_{pred} - 1.96*\\sigma_{m_p},m_{pred} + 1.96*\\sigma_{m_p} \\right)$$\n",
    "\n",
    "Because we estimated $\\sigma^2$ from the sample, that means $\\hat{\\sigma}^2$ has its own probability distribution with its own variance, which will have the effect of making the tails of our distribution for our estimators bigger. To take this uncertainity from $\\hat{\\sigma}^2$ into account, we will use a Student's t distribution instead. Using the student's t distribution will make our confidence interval a bit wider.\n",
    "\n",
    "We can calculate the 95% confidence points by using n-2 as our degrees of freedom for student's t (since our calculation of $\\hat{\\sigma}^2$ from the sample MSE had n-2 degrees of freedom).\n",
    "\n",
    "$$\\left( m_{pred} - t(.975,n-2)*\\hat{\\sigma}_{m_p},m_{pred} + t(.975,n-2)*\\hat{\\sigma}_{m_p} \\right)$$\n",
    "\n",
    "Student's t is distributed as $t(v) = \\frac{N(0,1)}{\\sqrt{\\chi^2(v)/v}}$, where $v$ is the degree of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For 68% level, T-score is: 1.01\n",
      "For 95% level, T-score is: 2.02\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import t\n",
    "t1 = t.ppf(1.0-(alpha1/2.0),n-2) #percentile function\n",
    "t2 = t.ppf(1.0-(alpha2/2.0),n-2) #percentile function\n",
    "print('For 68% level, T-score is:',round(t1,2))\n",
    "print('For 95% level, T-score is:',round(t2,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### intercept confidence interval calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "population intercept: 3\n",
      "central: [3.6261219]\n",
      "confidence interval using Students t\n",
      "lower 68%: [3.28186624] upper 68% [3.97037755] In interval? [False]\n",
      "lower 95%: [2.93606735] upper 95% [4.31617644] In interval? [ True]\n"
     ]
    }
   ],
   "source": [
    "#calculate the confidence interval on m, using student's t and estimated sigma_m_p\n",
    "SE_b_p = math.sqrt(V_best[0][0])\n",
    "b_pred = theta_best[0]\n",
    "interval_b1 = np.array([b_pred - t1*SE_b_p,b_pred + t1*SE_b_p])\n",
    "interval_b2 = np.array([b_pred - t2*SE_b_p,b_pred + t2*SE_b_p])\n",
    "print('population intercept:',b)\n",
    "print('central:',b_pred)\n",
    "print('confidence interval using Students t')\n",
    "print('lower 68%:',interval_b1[0],'upper 68%',interval_b1[1],'In interval?',(b >=interval_b1[0] and b <=interval_b1[1]))\n",
    "print('lower 95%:',interval_b2[0],'upper 95%',interval_b2[1],'In interval?',(b >=interval_b2[0] and b <=interval_b2[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### slope confidence interval calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "population slope: 4\n",
      "central: [3.72264689]\n",
      "confidence interval using Students t\n",
      "lower 68%: [3.42849815] upper 68% [4.01679562] In interval? [ True]\n",
      "lower 95%: [3.1330308] upper 95% [4.31226297] In interval? [ True]\n"
     ]
    }
   ],
   "source": [
    "#calculate the confidence interval on m, using student's t and estimated sigma_m_p\n",
    "SE_m_p = math.sqrt(V_best[1][1])\n",
    "m_pred = theta_best[1]\n",
    "interval_m1 = np.array([m_pred - t1*SE_m_p,m_pred + t1*SE_m_p])\n",
    "interval_m2 = np.array([m_pred - t2*SE_m_p,m_pred + t2*SE_m_p])\n",
    "print('population slope:',m)\n",
    "print('central:',m_pred)\n",
    "print('confidence interval using Students t')\n",
    "print('lower 68%:',interval_m1[0],'upper 68%',interval_m1[1],'In interval?',(m >=interval_m1[0] and m <=interval_m1[1]))\n",
    "print('lower 95%:',interval_m2[0],'upper 95%',interval_m2[1],'In interval?',(m >=interval_m2[0] and m <=interval_m2[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate confidence interval on regression line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we've found the confidence interval on the fit value of the slope and intercept.\n",
    "\n",
    "But what if we want to find the confidence interval on the predicted y, at every point x?\n",
    "\n",
    "$y_{p} = b_{p} + m_{p}x$\n",
    "\n",
    "So to find the expectation value and Variance of $y_{p}$ is just a matter of doing error propagation.\n",
    "\n",
    "$E[aA+bB] = aE[A] + bE[B]$ \n",
    "\n",
    "and \n",
    "\n",
    "$V[aA+bB] = a^2*V[A]+b^2*V[B]+2ab*Cov[A,B]$\n",
    "\n",
    "so\n",
    "\n",
    "$E[y_{p}] = E[b_{p}] + xE[m_{p}] = b + mx $\n",
    "\n",
    "and \n",
    "\n",
    "$V[E[Y|X=x]] = V[y_{p}|X=x] = V[b_{p}] + x^2 V[m_{p}] + 2x Cov[b_{p}, m_{p}] = V_{00} + x^2V_{11} + 2xV_{01}$\n",
    "\n",
    "$V[y_{p}|X=x] = \\frac{\\sigma^2}{n*V[x]}\\left(\\bar{x^2} + x^2 - 2x\\bar{x}\\right) $\n",
    "\n",
    "Again, since we don't know $\\sigma^2$, we'll use $\\hat{\\sigma}^2$ and get the standard error for $y_{p}$\n",
    "\n",
    "$\\hat{\\sigma}_{y_p|x} = \\sqrt{\\hat{V_{00}}+x^2\\hat{V_{11}}+2x\\hat{V_{01}}} $\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SE2_y_p(x):\n",
    "    return V_best[0][0]+(x**2)*V_best[1][1]+2*x*V_best[0][1]\n",
    "\n",
    "X_new = np.linspace(0, 2, 201).reshape(201, 1)\n",
    "X_new_b = np.c_[np.ones((201, 1)), X_new]  # add x0 = 1 to each instance\n",
    "y_predict = X_new_b.dot(theta_best)\n",
    "y_mu_upper = np.array(list(map(lambda x,y: y + t2*math.sqrt(SE2_y_p(x)), X_new, y_predict)))\n",
    "y_mu_lower = np.array(list(map(lambda x,y: y - t2*math.sqrt(SE2_y_p(x)), X_new, y_predict)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEVCAYAAAARjMm4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XlclVX+wPHPYRcB910RFcRdU3JfcEVxB0vTNJ0mK61s2pvJSZvKqWmvX9M0M20uZQruilsulUtquaGi4L6LKyD7Pb8/zhXJAe4FLxfE7/v14qU897nnOffhcr+c5XuO0lojhBBCFMSlpCsghBCi9JNgIYQQwiYJFkIIIWySYCGEEMImCRZCCCFskmAhhBDCJgkWQgghbJJgIYQQwiYJFkIIIWxyK+kK2FK1alUdEBBQ0tUQQog7yo4dOxK11tUcVV6pDxYBAQFs3769pKshhBB3FKXUMUeWJ91QQgghbJJgIYQQwiYJFkIIIWySYCGEEMImCRZCCCFskmAhhBDCJgkWQgghbJJgIYQQwiYJFkIIIWxyaLBQSj2hlNqulEpXSn2VzzmvKqW0UqqPI68thBCi+Dh6uY/TwOtAGFDu1geVUo2AEcAZB19XCCFEMXJoy0JrHa21XghczOeUT4AXgQxHXlcIIUTxctqYhVLqPiBDa73cWdcUQgjhGE5ZdVYp5QO8CfSz8/yJwEQAf3//YqyZEEIIezirZTEdmKm1PmLPyVrrz7XWIVrrkGrVHLYcuxBCiCJyVrDoDTyllDqrlDoL1AO+V0q96KTrCyGEuA0O7YZSSrlZy3QFXJVSXkAWJli45zp1G/AMsMKR1xdCCFE8HN2yeAVIBV4CHrT+/xWt9UWt9dkbX0A2cFlrnezg6wshhCgGDm1ZaK2nAdPsOC/AkdcVQghRvGS5DyGEEDZJsBBCCGGTBAshhBA2SbAQQghhkwQLIYQQNkmwEEIIYZMECyGEEDZJsBBCCGGTBAshhBA2SbAQQghhkwQLIYQQNkmwEEIIYZMECyGEEDZJsBBCCGGTBAshhBA2SbAQQghhkwQLIYQQNkmwEEIIYZMECyGEEDY5NFgopZ5QSm1XSqUrpb7KdbyjUmq1UuqSUuqCUmqeUqqWI68thBCi+Di6ZXEaeB344pbjlYDPgQCgPpAEfOngawshhCgmbo4sTGsdDaCUCgHq5jq+Ivd5SqlPgA2OvLYQQojiU1JjFt2B2BK6thBCiEJyaMvCHkqpVsBfgaEFnDMRmAjg7+/vpJoJIYTIj1NbFkqpQGAFMEVr/WN+52mtP9dah2itQ6pVq+a8CgohhMiT01oWSqn6wBrgb1rrmc66rhBC3E0s2kLCpQSHl+vQYKGUcrOW6Qq4KqW8gCygBvAD8H9a688ceU0hhBBGwqUEen7dk+uZ1x1etqO7oV4BUoGXgAet/38F+CPQEHhVKZV848vB1xZCiLvK+ZTzfLjlQz7d9ikAARUD6F6/O5+Ef+LwaymttcMLdaSQkBC9ffv2kq6GEEI43ObNsH49hIZCp072PSc1M5XFcYuZuXsmMfExZOtsBjUexJIHlvzuPKXUDq11iKPq6vTZUEIIIUyg6N0bMjLAwwPWrs0/YFi0BRdlOoImL5/Mlzu/pK5fXZ7r/BxjW42lefXmxV5fCRZCCFEC1q83gSI72/y7fv3/BosDiQeYuWsms/bMYtnoZbSo3oIpHaYwpuUYQgNCcXVxdVp9JVgIIcq0onT1OENoqGlR3GhZhIaa48kZyXy18yu+2fUN205vw0W50LdhXzKyMwBoXbN1idRXgoUQoswqTFePs3XqZOqzfj106pZG9eBTQCO01ryw+gUaV2nMu/3e5YEWD1DLt+TXXZVgIYQos+zp6ikpWmuy6/zM0RYzefvn72m4ryE7Ju7A19OXg08epK5fXduFOJEECyFEmZVfV48jFaWb65td3zBt/TSOXDmCt7s3kU0jGdtqLFprlFKlLlCABAshRBmWu6unOMYs7O3munj9It/Hfk9ks0iql68OQGDlQKaHTmd40+H4ePg4rlLXrsGyZY4rz0qChRCiVHH0gHSnTsXX9VRQN1d6VjrLDy3nm93fsOzgMjItmXi5eTHhngmMaz2Oca3HFepaBd6XxERYvBiiomDNGlMZB5NgIYQoNUrzgHRe8uvmupZ+jYYfNuRi6kVqlK/BE+2fYFzrcbSuUbSZTHnel/qnYcECiI6GDRtMxAoIgCefhIgI6NLFUS8TkGAhhChFSvOAdF5/2d/o5opal8BV/9nMvnyeTnyCn6cfUzpM4d4699KnYR/cXG7vo/bGfamXfYQRadHUGxUNxzeZB5s0gZdegshIaNMGlLqta+VHgoUQotRwxoB0UeT1l33j1mYcYua+mWzO3IxKUPTRfci2ZOPq4srUHlMdc/H9+xlzNJoBlija8BtoSPa6B/72NxMgmjZ1zHVskGAhhCg1intAuqhyWjwqjfQsF9av92Bd1r/4yw9/oUX1FrzV5y0eaPEA9SrUK/I1clouPTSdyu004w/R0bB/P/5ApRadWFvvHSr9YThtRzR01EuzmwQLIUSpUpwD0kVh0Rb8Wm2EIbMgeB5uMf8hNPQ+GrZ8mIFBA2lVoxXKRtePrUH7zT9b+HOvLQzOjKYm0aCPgIuLecLkyTB8OL61a+O9GVauh/TNzr9HEiyEECIP6VnpTFs/jdl7ZnPi2gm82/rQREfy5N8bWT+oa1DDp4bNcvIdtM/Kgo0bISqKFjMXsC7jDBm4s4a+7It4hYH/GgJVq9oux0kkWAghhNXppNPsPb+Xfo364eHqweKDi2lZoyVv932bIcFD8Hb3LnSZuQftSU/nxGdr6PSfaFi0CC5eBG9vMtoPYPzPESzOHkiaZwXWPgdUzb+ckhj8l2AhhLirJaUnseDAAmbtnsXaI2vx8/Tj3HPn8HD14LdHf8PD1eO2yu/VIYXdrjEMsUQx0LIUv2+SwM8PBg82A9RhYVTx9ubRzRC8Pv+uqpIe/JdgIYS4a3218ysmLZtEalYqDSo24C/d/sKYlmNyAkSRA8XVq7B0KURF0SEmhm8zUrnuXZWk3iPxezzC9Cd5/L5sW2M19gz+a63Zd2EfcRfjilbvAkiwEELcFbTWbD+9nVm7ZzGm1Rja12lPi+oteKj1Q4xtPZZOdTvZHKgu0IULpmspOtpkUWdmQu3a8PDDEBGBd7dueLvd3kduXgHlStoV0rPSqeFTg80nN9Pliy5F6i6zRYKFEKJMO3L5CLP3zGbW7lnEXYzDDU+8kpvT/r72hNQOIaT2bew8euqUyaKOijKD1RYLNGgAU6aYLqb27c2sJgeyaAs7Tu9gZcJKYuJj2HJyC1M6TOHdsHe5t/a9/HvwvwlrFIb/X/wdel2H7sGtlHoCGA+0BL7VWo/P9Vhv4P8Af2ArMF5rfcxWmbIHtxCisLIsWbi5uJFlyaLmOzW5mHqRtpVC2TPnQbL3ROKpKxZ9NtHhwzdzILZsMceaNTPBISICWrd2eBZ1amYq5dzLobWm2afNOJB4AICQ2iH0b9Sf4U2H07ZW2989p7TvwX0aeB0IA8rdOKiUqgpEA38ElgB/A+YCHR18fSHEXSotK41lB5cxc/dM4i7GsW/SPtxc3JgdMZum1Zoy+1N/dm0HSzZkuBZiNpHWsG+fCQ5RUbBrlzneti288YYJEE2aOPS1ZGZnsunEJmLiY1iZsJLkjGQOPnkQpRRPtn+Sil4V6duwL9XKV3PodQvi0GChtY4GUEqFALkXZI8AYrXW86yPTwMSlVJNtNYHHFkHIcTdZc+5PXy09SPm7ZvH1fSr1PSpyegWo0nNSsXb3ZuwwDCgkLOJtIZff73ZgoiLM62Fzp3hvfdg+HCzaF8x+GjrR7y8+hWuZyfhqtzo6t+F+5vfn9NamnTvpGK5ri3OGrNoDuy68Y3WOkUplWA9LsFCCFEo+y7so3K5ytT0qUn8pXi+3fstkc0iebDlg/Rq0AtXF9f/eY7N2UQWi8l8uxEgjh0DV1fo2dOMQQwbBrUct73p9czrrD+6Pqf1sGDkAppVa0bGuUakbx+Dy8H+uJ/uyYwVfqUio91ZwcIHuHDLsauAb14nK6UmAhMB/P0dO0gjhLgznU46zdy9c5m1Zxa/nvmV6aHT+WuPvzKo8SDOPXeO8h7lbZbxP7OJMjPN8t7R0Wag+uxZ0+zo1w+mTTO5EFWqOPR1xF+KZ9KySWw8tpH07HTKuZUjNCCU9Kx0U6XYgbB0IJZsyCxMd1kxc1awSAb8bjnmByTldbLW+nPgczAD3MVbNSFEaWbRFgbNGURMfAwaTUjtED4I+4CRLUYC4O7qjruru/0FpqXB6tUmQCxeDJcugbc3hIebQerwcJM05wBX0q6w5vAaYuJj6FCnA4+0e4Sq3lU5m3yWyfdOpn9gf7rV74aXm1fOc0o6+S4/zgoWscBDN75RSpUHGlmPCyFEjhsD1b+e+ZU3er+Bi3IhsHIgf+3xVx5o8QDBVYMLX2hyMqxYYQLE0qXm+woVYMgQM0AdFgblytkux05v/fQWiw8uZuvJrWTrbPw8/ajnZ1akrehVkd2P7873uaV15V1HT511wwSgVzED3I8AWUAlIB74A7AMmA700FrbnA0lU2eFKPuyLdmsO7qOOXvmELU/imvp16jpU5O4J+Lw8yziX/mXL+dkUbNypWlRVKtmxh4iI81YhMftLeUBcCHlAisTVpJwKYFXQ18FoO/MvlxJu0L/Rv0JCwyjQ50OhWv9OICjp846OlhMwwSK3KZrracppfoAnwD1uZlncdRWmRIshCibtNZk62zcXNz4945/M3HpRHw9fIlsFsnoFqPp2aBn4XeYO38eFi40LYi1a83KrnXqmNZDZCR07WoGrW9T7PlYvo/9nhXxK9h+ejsaTS2fWhyechgvNy8yszOdHhxuVaqDRXGQYCFE2XLo4iHm7JnDnL1zeK7TczzS7hEuXr/IuqPrGBg0kHLuhewOOnHiZhb1Tz+ZWU2NGpngEBkJISG3nUV9LvkcKxNWMjBoIFW8q/DBlg94dtWzdKjTgQGBAxgQNIC2tdriohybrX07SntSnhBC/A+tNR//8jGzds9i2+ltKBShAaHU8asDQBXvKoxoNsL+AuPjbybJ/fKLOdaiBUydaloRLVveVhZ1liWLrSe3siJ+BSviV/DrmV8B+DbyW0a1GMX4NuMZ13oclctVLvI17jQSLIQQxeJq2lV2nNlBrwa9UEoxe89ssixZ/KPvPxjVYhR1/eraLuQGrWHvXhMgoqNht3WAOCQEZswwAaJx49uq75mkM6RkphBYOZDjV4/T9cuuuCpXOtXrxBu93mBA4ABa12wNmEHqu40ECyGEw6RnpbP80HLm7J3DkrglAJx77hwVvCqwZuwafD3zTK3Km9awffvNJLlDh0xroWtXeP99EyBuIw8rMzuTzSc3s+KQaT3sOreL+5vfz9wRc2lYqSFLHlhCV/+ud2VgyIsECyGEQyw6sIjxi8ZzJe0K1ctXZ2K7iYxuOTpnNpNdgSI7GzZtuhkgTpwANzczc+nZZ2HoUKhZs8h1vJJ2JefD/56PuxN7dUvOkhozes9gUONBOefm/r+QYCGEKAKtNTvP7mT2ntkMCBxA74a9Ca4azODGgxnTcgy9G/a2fyZTZiasW2eCw8KFcO4ceHqa3Ie//c1kUVcu2thARnYGPx//OWfs4fjV4yQ+n8j2X9w59PWzuGS64H6yNzNWVCg1+QyllQQLIQRglkWylQgWlxjHd3u/49u935q9IVzcqFG+Br0b9qZJ1SZ8M/wb+y6WmmqyqKOiTBb1lStQvjwMHGi6l8LDwbcQXVa5aK1RSjFz10wmLZ9EckYy7i7udPXvyrhW48jIzmD9eney94wodUtqlGYSLIQQbN5sdvq8scRE7r0erqVfw8/TD601YbPCOH71OD0CevCnjn9iRLMRVPG2c+2kpCRYvty0IJYtg5QUqFjRZFFHRkLfvnZnUecObG3vTeen4z/ltB4+7P8hfRr2IbhqMGNajmFA4AB6Nej1u26w0rqkRmkmwUIIwfr15oMzO9v8u2TdWba7zOO72O+IvxTPqWdO4ebixqyIWTSs1JDavrXtK/jSJViyxLQgVq2C9HSoXh0efNC0IHr2BPfCJa/dCGzpbhcg9mE8m/xAanYKHq4edK/fHXcXU177Ou1pX6d9nmWU1iU1SjMJFkKInL+002v8hA6dxluZ67DEWGhVoxVPd3ia9Kx03Dzc6Orf1XZhZ8/ezKJet85kUderB489ZloQnTsXOos6MzuTn47/xPJDy9n7S3UyMp7HklYZyp+mlX6IVx4YQGhAKD4ePnaXmdd+1iJ/EiyEuIslpSexKG4R7QLbsXZtU/77QyarPI4zIeQVRrYYSbNqzewr6Pjxm0lyP/9spr0GBcFzz5kWREhIkZLk5sXOY27sXFYfXs219Gu4u7jTv9Y4axeSKx4zt/P+Wuh0eykWwg4SLIS4y6RmprL80HK+i/2OpQeXkpaVxtTuU3mt52t07BgKxKHs+WA/ePBmgLixJE+rVvDqq6YF0bx5oQJEtiWb7ae3s/HYRp7r/BxKKZYcXMKWk1sY2Xwk4UHh9G7QG19PXzY3ly4kZ5O1oYS4i1i0hfof1OfktZNUL1+d+5vdz6gWo+hUr5PtdY20hj17buZA7N1rjrdvb4LD8OGmNVEIV9KuEBMfw7JDy4iJjyHxeiIuyoWEpxIIqBhAUnoSPh4+9gUv8TuyNpQQwi7Zlmw2HtvId3u/Y3/ifjaM34CLcmFaj2nUr1if0IBQ27kQWpu1l24ssxEfb1oL3brBhx+aAFGvnt110lqz5/weqnlXo5ZvLdYcXsMDUQ9QpVwVBgQNIDwwnLDAsJw1lwqV8S2KlQQLIcqY/Rf2868d/+L72O85k3wGb3dvhgYPJTUrFW93bx5u+3DBBWRnm9Vbo6LMaq4nT5os6t694fnnTRZ1jRp21yc5I5m1h9ey/NBylscv5+S1k7zV5y1e6PICYY3C2PzwZu6tfW+e+2aL0sOuYKGU+gx4FKijtT59y2PBwB7gn1rrKY6vohCiIFprtp3eRv0K9anhU4Pfzv7GP7f/k/CgcB5o8QADgwba3p86IwN++OFmFvWFC+DlZbKo33wTBg2CSpXsrlNSehK+nr6kZaVR691aJGck4+vhS79G/ZgeOp3woHDAtBw61rW5B5ooBextWWzGBIv2wMJbHnsfuAZMc1y1hBAF0Vqz48wOvo/9nu9jv+fY1WM5f60PbzKc8KBw2wvgpaaaHeSiokwuxNWr4ONjAkNEBAwYYL63Q1pWGhuPbWTZwWUsj19OjfI1+OkPP+Hl5sXbfd6mSdUmdPHvgofr7e9MJ0qGvcFii/Xf3wULpdRAYAAwWWt92cF1E0LkISM7gxaftuDQpUO4ubjl/LU+tMlQAMq5l8t/A6Fr10wWdVSU+ff6ddNiGD7cDFL36WNaFIUwff103t70Ntczr+Pl5kXPgJ4MCR6S8/jj9z5e5NcqSg+7goXWOk4pdQkTLABQSrkD7wF7gX8VT/WEuLtprfnt7G/Mi53H2ZSzfDn0SzxcPbi/+f00qtSIoU2G2t6A5+JFs/5SdLTJos7IMCu3PvSQaUH06GFXFrVFW/j1zK8siVtiZi89GENV76o0qNSA8a3HEx4UTs8GPfF293bQqxelSWEGuLcAXZRSSpv5tlOAxkAfrXV2sdROlAr2LDAnHGv/hf3M3D2T72O/J+FyAq7KlbDAMLIsWbi5uPF6r9cLLuDMGTP2EBVlfnjZ2VC/PkyebFoQHTvanUWdcCmBv//0d5YeWsrZ5LO4KBc61e3EueRzVPWuyrjW4xjXetztv2hRqhU2WIQDwdZWxlRgodZ6rb0FKKUCgE+BTkA6MB94WmudVYh6CCcqaIE54Thaa3af203DSg3x9fQlJj6Gt39+m94Ne/NS15cY3mS47QX7jh69OcV10yYz7bVxY3jhBRMg2ra1K0nuxNUTLD24lOCqwfRq0AuA7/d9T//A/gwKGsSAoAFU9a7qgFct7iSFCRabrf+2B7oDnsCzhbzep8B5oBZQEVgNTAI+KmQ5wkluXWBOlnJ2HK01e8/vNYPU+77n4MWDfDPsG8a2HsuEeyYwtvVY2x/KBw7czKL+1ewTTevWMH266WJq1syuALH15FaWHFzC0oNL2XVuFwCTQibRq0EvGlVuROLzibi7FtxVJS3Qsq0wwWIrYAEeBroC/9BaHy7k9RoAn2it04CzSqkYoHkhyxBOJEs5F48raVfo+J+OxF2Mw0W5EBoQyjMdn2FA0ACggD2etYZdu25mUe/bZ4537Ahvv20CRKNGNq+fnJFM7PlYOtTtAMAfFv+BA4kH6FKvC2/3eZvBwYMJrhKcc749gUJaoGWb3cFCa52klNqHaVWcBd4owvU+BEYppdYDlTAzqaYWoRzhJLKU8+2zaAtbT24lan8UAO/0e4eKXhXp6t+VKR2mENE0gho+BSS5WSywdevNLqbDh8HFBbp3h8cfh2HDoG5dm/U4duVYTuth3dF1eLp6kvhCIh6uHsyOmI1/BX/bg+X5kBZo2VfYDO5fgBbAy1rrpCJcbwPwCCYvwxX4mv/N20ApNRGYCOB/GxuyC8coiaWcy0KXxi+nfmHmrpksOLCAU0mncHdxZ3jT4TmP/2fIf/J83ubNsGFtFoMr/kjzuGiTRX3qlJmx1KcPvPyyyaKuVq3A62dbzLwTVxdX3t30Ls+tfg6AxlUa88S9TzA4eDCuygxyt6nZ5rZeq7RA7wJaa7u+AHcgAdiGdQHCwnwBLsBx4C+Y8Y4qwCLg7YKe165dOy2K36ZNWr/5pvm3pG3apHW5clq7upp/S0Od7JGRlaFXxa/SGVkZWmutX1j1gvZ63UsP+26Ynrlrpr6cerngAtLS9L53l+svXR/W56mqNegsz3JaDx+u9axZWl+28XytdUpGil64f6GesHCCrvZ2Nb0mYY3WWutdZ3fpdze9q+MS44rtZ12a3kNCa2C7LuTndEFfhfmwfxkzZtGxSBeCqoAGKuQ6NgzYW9DzJFgUv9Ly4Xzjw+axx0xdwPz75pslW5+C7kdqZqpefGCxfmjBQ7rS3ytppqFXxa/SWmudmJKok9KTCr5ISorWUVFajxmjtZ+f1qCv4qtnMVqPcJmv/zEt2a66nk8+r4d8O0SXe72cZhq6wowK+oH5D+jtp7b/z2sqDT9rUfwcHSwK7IZSSlUGwoBWwPPAe1rrLQU9p4AWTKJS6gjwuFLqHcAHeAjYVZTyhOOUhv7m3AOkbm43UwBKqkvDngHbA4kHaP/v9iRlJFHBswJDgocQ2TQyZze5fKe6Xr1q9qCOioIVK8yyG1WqwIgR7G8aQaepfUjO9MTDA57pl3cRcYlxLIpbRHn38kxuP5nK5Spz8tpJ/tj2jwwNHmq2F81jULo0/KzFncnWmEUYMAcz3fV94KXbvF4E8AHwIpANrAP+dJtlitvkrP7mgsYhcn+IATzyCPj7F9+Yha0xkVs/VGPWXeOozzKi9kcRXCWYN3q/QVDlICa0mcCAoAH0atCr4HWPEhNh0SIzQL1mjSm0Vi2YMMHkQHTvDm5uNAVWdMm7bttPb2f+vvksPLCQuItxAAxrMozJ7Sfj6uLKjok7bL5uGVsQRebIZkpxfEk3lHMUd3+zre4PZ3aP2HOtG+eoVnO0y9hw7T7dQzMNXeudWnr6+un2XejkSa0/+UTrnj21dnEx/WoBAVo/+6zWP/+sdXZ2gU+/nnFdr4xfqS0Wi9Za67HRY7Xba266zzd99MdbP9bHrhwr7EvPeW0ytlD24cxuKHH3KO4ZT7a6P5w5RbeguiRcSmD14dU82vFR1q5VTF4TwznP/YxsPZnIppG2d5Q7cuRmDsRmax5rkyZmBlNkJLRpU2CSXOL1RJYeXMqiuEWsSljF9czrxE6KpVm1ZrzR6w0+GvCR7dVkbSiJ2W3izifBQjiFPd0fzvoQy10Xdw9Nrba/MfWHBSyMW8je82ar0F4NetGpU2N+bPcp3u7eBW/ruX//zQDx22/m2D33wOuvmyS5pk0LrI9FW3BRLqw9vJZ+s/ph0Rbq+NZhfOvxDG0ylMDKgQDUq2D/jnRCOJrswS2cprTkTmRZstj4czpbfyqParaQl3cOx0W50M2/G8OaDGNYk2EEVAzIvwCtTVC4sczGgQPmeOfOJjhERECDBgU83SzzEb0/mugD0YxpOYYXurzAtfRrvLPpHYYGD6Vtrbay77S4LbIHt7hjlWT3x/XM66xKWMWCAwtYenApz3d+npdffomk9N7U8P+CQY0HUa18AUluFgts2XKzBXH0qJmy1aMHPPmkyaKuXbvAOmiteeWHV5gbO5eEywkoFF38u+QEJj9PP17r+ZrjXrQQDiTBQpRpWmtGRY1iSdwSUrNSqehVkcGNB+ds5enr6cuEeybk/eSsLNiwwQSHBQvMst/u7tC3L0ydCkOGQNX8F/rLzM5k47GN7Dq3i2c6PYNSip3ndhJYOZAXurzA0OChBS/zIUQpIsFClCnHrhxjUdwiDl8+zAf9P0Aphbe7Nw/f8zDDmgzLN/8gR3q6mdoaFWWmul66BN7eZovRiAgYOBAqVMj36amZqaw+vJro/dEsObiES6mX8PHwYWK7ifh4+LDkgSUFD5ALUUpJsBB3vP0X9vPt3m9ZHLc4Z3ntltVbkp6VjqebJ18O/bLgAlJSTHJcdDQsXQpJSeDnB4MHmxlMYWEmYOTjatpV3F3d8Xb35t+//pspMVNyWjARTSPo16hfzu5xEijEnUqChbjjpGam8sORH+hUrxOVy1Vm9eHVvPHjGznLaw9rMoygKkEFF3LligkMUVEQEwNpaaZLaeRI04Lo3dtMmcrH+ZTzLI5bTPT+aNYcXsN/hvyHca3HcX/z+9EXmpK0O5TeNd3p1MTBL16IEiKzocQd4VzyOZYdWsaSg0ty8g++HvY141qP42raVTItmbY3CrpwwXQtRUWZpI7MTDMoHRFhWhBdu5q1RgqQkpFC+Jxwfjr+ExZtoWGlhgxvMpwJbSbQvHpzu/d1KC0zw0TZJbOhxF1Ba01KZgo+Hj49SpaPAAAgAElEQVScTjpN3ffqotHU86vHhDYTGBI8hB71ewBQwSv/MQROnjSD09HRsHGjmdXUsCE8/bQJEu3bm70h8nH0ylGi9kWRlJHEtNBplPcoT/Xy1Xml2ytENI2gVY1Wv5vias/aS7JRkLgTSbAQpUZmdiY/Hv+RxXGLWXJwCW1rtWXeffOo7VubjwZ8RFf/rrSu0dp2/kFCws0ciK1bzbFmzeAvfzEtiFatCsyiPnL5CPP2zWPevnlsP21atd3rd+dV/SpKKebdNy/P523eDMeP214EURbzE3ciCRaiVPjz2j/z6bZPuZp+FU9XT3o37M3gxoNzHn+i/RP5P1lrs73ojRyIXdaFjNu1gzfeMC2IJgUPHsRfiqd+hfq4u7rzrx3/4q2f3+Le2vfyVp+3iGwaSaPKBW9VeuuquY88AuPG5R0EZDE/cSeSYCGKRUF98ocvH2ZJ3BIzxXRkNB6uHlT0qkhE0wiGBA+hT8M++Hj4FHwBrWHHjpstiIMHTWuhc2d47z0YPhwCAgosIi4xjvn75jNv3zx2ndvFygdX0q9RP57q8BSPhTxWcBb3LW5dNdffP//WgmxVK+5EEiyEw+XVJ189OIHPd3zO0kNL2XdhHwDNqjXj+NXjOUlqNmVnm8JvtCBu9Pn07GnGIIYNM8t+23Dy2knCZ4ez5/weADrV7cR7/d6jVY1WANT2LTgTOy+FbS3IYn7iTiPBQjjc+vWQ7pqIpXkM6YktWL++DT3qnuP9Le/TvX53/njPHxkcPDhngbwCZWaaLOqoKFi4EM6eNZ/G/frB9OkmF6JKPpsMYQbK913Yx7x98/D18OXZzs9Sy6cWDSo14OF7HiayWSR1/ere9muW1oIo6yRYCIfQWrPn/B6WHlzKd97LsDyzBVwsuGx5kdDQNrSv04HEFxLx8/SzXVhaGqxebVoPixffzKIODzcD1OHhJmmuAHvP72Xu3rnM3z+fA4kHUChGtRgFgKuLK4tGLXLEy/4daS2IskyChSiy65nXOXL5CM2rN8eiLfT+pjeJ1xNpV6sdfwicitfxgYye2s76AepacKBITobly02AWLbMfF+hgll/KSLCZFGXK1dgfQ5dPERg5UCUUryz6R1m7p5Jj/o9eKr9UwxvOpyaPjUd+vqFuJtIUp4olGNXjrHs0DKWHlzKuqPrqOVTi4SnElBKseHoBhpXaUwtX9vjBgBcvgxLlpgAsXKlaVFUq2bGHiIjzVhEAVnUYGYxzd07l7mxc9lzfg+7HttFqxqtOHblGF5uXrJQn7hrSVKecKosSxauyhWlFH9e+2dm/DQDgEaVGvFou0cZGDQw59weAT1sF3ju3M0s6h9+MCu71q0LEyeaFkTXrjcTFQoQlxjHmOgx7Dhj9p3uUq8LH/X/iDq+dQCoX7F+EV6tECI/Tg8WSqlRwKuAP3AWGK+1/tHZ9RD5u5R6iZj4GJYeXEpMfAw/TviR5tWbE9YojKreVRkYNJDGVRrbvznPiROm9RAdDT/+aKa9NmoEzzxjWhAhIQVmUQOcunaKefvmUb18dUa3HE1dv7p4u3vzbr93ua/ZfbKLnBDFzKnBQinVF3gLGAn8AtjZXyGc4eDFgzy8+GE2ndiERVuo5l2NIcFDclZK7RHQw77WA8ChQzdzILZtM8datIC//tW0IFq2LDCLGsxiffP3zWdu7Fx+PPYjGs2YlmMY3XI05T3Ks3HCxtt5uUKIQnB2y2I68JrWeov1+1NOvr6wSslI4YcjP7AifgVtarZhYruJHN9Xk6Mn0nmowV94tOdA7q1zr/1LamsNe/feDBB7TA4DISEwY4YJEI0b2ywmKT0JX09fAB5a+BAx8TE0q9aM6aHTGdliJI2r2C5DCOF4TgsWSilXIARYrJSKB7yAhcDzWutUZ9Xjbvfptk9ZeGAhG45tICM7g/Lu5ankVYnNm2FImB8ZGb/wnQc8shZcbKUfaG1aDTe6mA4dMq2Frl3h/fdNgPD3t1mnK2lXWHhgIXNj57LuyDqOTDlCLd9avN7zdf7R9x+0qN7CMS9e5JBVb0VhObNlUQNwB0YA3YBMYBHwCvCX3CcqpSYCEwH87fiwKaqy/gtzPfM6646sY9+FfTzf5XkAFhxYwKlrp3ji3icIDwqnq39XPN08mTHDzsXtsrPh559N62HBAjMe4eZmZi49+ywMHQo17ZuieiDxAC+ueZGY+BgysjOoX6E+UzpMyXm8Xe12DrgL4lay6q0oCmcGixuth4+11mcAlFLvkUew0Fp/DnwOZupscVSmrP7C3NhWdPmh5aw/up707HR8PXyZdO8kynuUZ9GoRTm7tuVW4HIVmZmwbt3NLOrz58HT0+Q+/O1vJou6cuU865M7IN9zbxrLDy2ncrnKhAaE4uvhy69nfmVSyCRGthhJhzod7B80F0Umq96KonBasNBaX1ZKnQRKRWJHWfmFSc1MZf3R9bSv054q3lVYHLeYKTFTCK4SzKR7JzEgcADd63fH080TIM9AAXksV9EmFRatuplFfeUKlC9v9qCOjDR7Uvv6Fli3zZuhV58sMmqvQ+2Yg9c90aRkXWNk85GEBoRSx68Ox58+LgHCyWTVW1EUzh7g/hJ4UikVg+mGehpY6uQ6AHf2L0z8pXhWHFrBivgVrDu6jrSstJxd40a3HM3AxgNpWKlhocvt1CKJTseXw/tRJps6JQUqVTJdSxER0LevzSzq3Navh7QR/aHhWkjzI9gSwd8fHE3PBj1zzpFA4XyyjpUoCqdmcCul3IEPgdFAGvA98ILWOi2/5xRnBvedMmaRlpXGpdRL1PatzYmrJ/D/wIzjBFUOIjwonAGBA+gR0AMvN6/CF37pkmk5REfDqlWQng7Vq5slviMjzc1xd7erqNjzsczZM4clB5ew+eHN7N5RntBHF5FlycbjWDg/rPIq1fdZiLLE0RncstxHKXX48uGc1sMPR34gPCic+ffPB+CrnV/R1b+rfau25uXsWTP2EBVlxiKys6FevZt7UXfubFcWNZi9sb/c+SVz9sxhz/k9uCpX+jTsw2eDPiOgYsAdE5CFKGskWJRR2ZZsXF3MB/TwucNZeGAhYJbVCA8KZ1iTYfRq0KvoFzh27OYU159/NtNeg4JMcIiIMPkQdnYJnU85T1pWGv4V/PntzG+0/bwtnet15oEWD3Bfs/vKzHpMEujEnUzWhipDjl45mtN6+OXULxz/03E8XD0YGjyU0PqhhAeFE1QlqOgXOHjQtB6iosyucmD2n371VXY2imTF8eaE9lR0utd2UdfSr7HwwELm7JnDmsNrGNd6HF8M/YI2Ndtw7Olj+FcovinOJaGg2XISRMTdSIJFCYiJj+GZlc+wP3E/AA0qNmBEsxEkZyRTuVxlxrcZX7SCtYbdu29mUcfGmuPt28Nbb5lxiKCgQk8bnrxsMl/s/IK0rDQCKgbwQpcXGNNyDGAGqMtaoID8Z8uV1SnXQtgiwaKYJVxKICY+hpUJK3m649P0atCLil4VqetXl4ntJjIgcEDhFuW7lcVisqhvbDWakGAW5evWDT780ASIer9fZK+gacNZlizWHVnHskPLeC/sPVyUCzV8avDHe/7I6Jaj6Vi3410xgym/2XJlZcq1EIUlwaIYJKUn8ee1fyYmIYb4S/EAVKIBHbwu0asBdKzbkVVjVxX9AtnZZvXWG2MQp06ZLOreveHFF81U1+rV8336rR+EPXpotpzcyrd7vmVu7FzOpZzD18OXx0MeJ7hqMH/t8dei1/UOld/00jt5yrUQt0MGuG+T1prYC7HExMdQzq0ck9tPxqItNPiwAS2qtyDYtT+fPtOfzHOBeHqoondbZGSY/R+iosx+EBcugJcX9O9vBqgHDTI5EXbavBl+WJdNr56uJFdfTb9Z/fB09WRg44GMbjGa8KBwyrnbn1NxN5ExC3EnkAHuUmLFoRVE7Y8iJj6GU0lm8dyBQQOZ3H4yLsqFI1OO4KJcmDEDss6BpSjdFtevmx3koqPNjnJXr4KPjwkMEREmi9rHp1D1PpN0hu/2fsfsPbMJ6xBGp05vkJkdytfDvmZo8FAqeFUo3I24C8le2+JuJMHCDhZtYcfpHWw4toFnOz2LUoq5sXNZeGAhfRv1JaxRGGGNwn63Ac+Npb0L3W1x7ZrZgzoqClasMAGjUqWbSXJ9+pgWRSHN3TuXL3Z+wZrDa7BoCyG1QwiuGgyAu6s741qPK3SZjiZ/sQtRekk3VD4uXr/I8kPLiUmIYVXCKhKvJwJw6MlDBFYO5FLqJfw8/XBzsR1vbX4IXrxosqijomD1ahNZatY0ASIiAnr0sDuL+obM7Ew2ndiUs1nRyPkj2XZqG2NajmFMqzE0qdqkUOUVt7I8y0iCoCgJ0g1VTDKzM9lycgsBFQOoV6EeG49tZNzCcVQvX50BgQMIaxRG30Z9qV7eDBxXLpf3Kqt5ybPb4swZWLCAK19G4/frelws2VC/PkyebFoQHTvanUV9g9aaLSe3MHvPbObGziXxeiIHJh8guGownw/6HD9Pv1I7k6mszjIqy0FQ3F3u6mBx7MoxViasJCY+hrVH1nIt/Rpv9nqTl7u9TN9GfdkxcQdtaraxf7c4W44evZkDsXkzaM05Fcw/eZGlnhG8M6ctnToX7cP81zO/cv+8+0m4nICXmxdDgofwYMsHaVCpAUCpH4soq7OMymoQFHefuypYpGWlcSbpDA0qNSA1M5XGnzQmIzuDen71GNV8FGGBYfRu0BsAHw8f2tZqe/sXPXDgZg7Er7+aY23awPTpfH4xkkkfNyXbonDNgvUboFNn+4o9l3yO7/Z+Rx2/OoxoNoKGlRoSWDmQV7q/QkTTCPw8/W6/7k5UVldCLatBUNx9yvSYhdaaQ5cO8c/VMSw/GMMxtZ62tduw6eFNAMyLnUeL6i1oUrWJ47pntIadO2+2IPabLG06djTdS8OHQ6NGQOG7KJIzkll4YCGz98xmdcJqsnU2D7V+iK+GfeWYuotiIWMWoiTIQoI2XM+8nrPBz/iF4/l619fmgcRgXI+G8Y+J4fxpcJhjK2mxwNatN1sQR46YLOru3U2AGDYM6ua9obWtDxKLtuR0g/Wb2Y/Vh1fjX8HfDFS3HEPz6s0d+1qEEGWCDHDfQmvN7nO7iYmPISYhhk0nNnHiTyeoXr469zW7jyt7O7LkgzAsFxuAK6R1BgY74MJZWbBxowkOCxbA6dNmxlKfPvDnP5ss6mrVbBaT1+C31pptp7cxe/dsovZHsfOxnVT1rsrU7lOZ2n0qXfy7OG4cRQgh7HBHB4t1R9YxOno0Z5PPAtCmZhue6fgMFm0BYGDjgVTuDavehAxXB/QZp6ebvqLoaJNFnZhodo7r39+0IAYOhIoVi1z8+ZTzfLb9M2btnsWhS4fwdPVkUONBXEu/RlXvqnSr3+02Ki+EEEV3RwSLbEs2205vM62H+BgeC3mM8W3G07BSQ0IDQunfqD/9GvWjlm+t/3nubQ+cpqSYLOqoKFi61CTN+frC4MEmB6J/f7M3dRFdSLnA1fSrBFYOJDkjmekbptO9fnde7PIikc0iqehV9OAjhBCOUurHLCo3qqyZCJfTLuOiXOhQpwPPdHqGEc1GFLlMmwOOV6+awBAdbbKoU1OhShXTtRQZaUalPT2LfP3rmddZHLeYWbtnsTJhJYMaD2LByAWAmeVUVjYPEkKUnLtuzCJbZxPZJJL+gf3p07BPoZLh8pLvDKQLF25mUa9ZA5mZUKsW/OEPpgXRvbtZ2fU2vfLDK3y49UOSM5Kp61eXZzo+w4OtHsx5XAKFEKI0KpFgoZQKAvYA87XWDxZ0blDlIL4Y+oXDrp07Sapq+imSZiyA5GjYsMHMagoIgKeeMi2IDh3MrKbbsPvcbr7b+x2v9ngVTzdPKnlVYmTzkTzY6kG61+8uA9VCiDtCSbUs/g/YVhIXDgs6zDUVzVCi6GjZAkuApk3h5ZdNgGjTxu69qPNzJukMc/bM4Zvd37D73G7cXNwY1HgQnet15tnOzzrmhQghhBM5PVgopUYBV4BNQKBTLrpvX04ORNudO2kLnK19D8cHvY7/0xEmWDhI7PlYWn3WCou20L5Oez4e8DEjm4+kWnnb02iFEKK0cmqwUEr5Aa8BvYGHi+1CWsNvv91MkjtwwBzv3BneeQciIqjZoEGRis49ON6+Qzbrjq5j5u6Z1Chfg7f7vk2zas2Y0XsGQ4OH5iwBLoQQdzpntyz+BvxXa32ioOU1lFITgYkA/v7+9pVssZhP8htbjR49alZt7dEDnnzSZFHXrn1blb8xOJ5eIRaXH7+mYvfZJKafpoJnBSa2m3ij7rzQ5YXbuo4QQpQ2TgsWSqk2QB/gHlvnaq0/Bz4Hs9xHvidmZpos6qgok0V99qzJou7bF6ZOhSFDoGpVh9T/bPJZ1q2rQUaGwhLyEZZ7vqBqxgA+HfEBg4MH4+VW+A2JhBDiTuHMlkUoEAAct7YqfABXpVQzrbX9y7ump5sNgm5kUV+6BN7eZovRiAiTRV3BMctxp2SksPDAQmbunsnqw6v5172b8PDoQPpPU/HY9DpfLK1GJ1maSQhxF3BmsPgc+C7X989hgsfjNp+ZnAwxMaYFsWwZJCWBn5/Joo6MhLAwEzAc5HzKeZ5f/TzR+6NJzkimfoX6vNz1ZcLa1bZmg9eVFUSFEHcVpwULrfV14PqN75VSyUCa1vpCgU9MSDAL8qWlmS6lkSNNC6J3b5NVh2OWgN57fi/nU87Tq0Ev/Dz92HB0AyObj2Rsq7F0q98tJx+iXl673gkhRBlX6pf7CPHw0NsffdS0ILp2/Z8s6tvZtvJs8lnm7JnDzN0z2Xl2J82rNWfvpL3A75cGF0KIO81dt9wHrVrBxx/n+3BRt618bcNrTN8wHYu2EFI7hI/6f8SoFqNyHpdAIYQQN5X+YGGDPdtWWrSF9UfXM3P3TF7t8SoBFQNoX6c9L3V5ibGtx9KkahNnV1sIIe4od3ywKGgJ8tjzsczcPZPZe2Zz8tpJ/Dz9iGwaSUDFAPoH9qd/YP+SqrYQQtxR7vhgAb/fbS7bko2riyuXUi/R+rPWAPQP7M87fd9hSPAQyrmXK8GaCiHEnalMBIvUzFQWxS3im13fkGnJZPXY1VQuV5n598+nc73OVC9fvaSrKIQQd7Q7OlhsO7WNf27/J/P3zScpI4l6fvUY13ocWmuUUgxrMqykqyiEEGXCHRcs9l/YT70K9fDx8GHTiU3M3zef+5rdx9jWY2V/CCGEKCalP88iJEQv37Ccb/d8y8zdM9lxZgdfDPmCCfdMICUjBaUU3u6Oy94WQoiy4K7Ls4i/FE/td2uTrbNpW6st74e9T3hQOADlPcqXcO2EEOLuUOqDBcBznZ9jbKuxNK8uq/YJIURJKPXBIrByIH/v8/eSroYQQtzVysxo8ObNMGOG+VcIIYRjlfqWhT1uZzFBIYQQtpWJlkVeiwkKIYRwnDIRLG4sJujqmv9igkIIIYquTHRDFbSYoBBCiNtXJoIF/H4xQSGEEI5VJrqhhBBCFC8JFkIIIWxyWrBQSnkqpf6rlDqmlEpSSv2mlBrgrOsLIYQoOme2LNyAE0APoAIwFfheKRXgxDoIIYQoAqcFC611itZ6mtb6qNbaorVeChwB2jmifMngFkKI4lNis6GUUjWAxkDs7ZYlGdxCCFG8SmSAWynlDswGvtZaH8jj8YlKqe1Kqe0XLlywWZ5kcAshRPFyerBQSrkAM4EM4Im8ztFaf661DtFah1SrVs1mmZLBLYQQxcup3VBKKQX8F6gBhGutMx1RrmRwCyFE8XL2mMU/gaZAH611qiMLlgxuIYQoPs7Ms6gPPAq0Ac4qpZKtX2OcVQchhBBF47SWhdb6GKCcdT0hhBCOI8t9CCGEsEmChRBCCJskWAghhLBJgoUQQgibJFgIIYSwSYKFEEIImyRYCCGEsEmChRBCCJskWAghhLBJgoUQQgibJFgIIYSwSYKFEEIImyRYCCGEsEmChRBCCJskWAghhLBJgoUQQgibJFgIIYSwSYKFEEIImyRYCCGEsMmpwUIpVVkptUAplaKUOqaUGu3M6wshhCgaNydf7/+ADKAG0AZYppTapbWOdXI9hBBCFILTWhZKqfJAJDBVa52stf4JWAyMdVYdhBBCFI0zu6EaA9la64O5ju0CmjuxDkIIIYrAmd1QPsDVW45dBXxvPVEpNRGYaP02XSm1t5jr5ghVgcSSroQdpJ6OdSfU806oI0g9HS3YkYU5M1gkA363HPMDkm49UWv9OfA5gFJqu9Y6pPird3ukno4l9XScO6GOIPV0NKXUdkeW58xuqIOAm1IqKNex1oAMbgshRCnntGChtU4BooHXlFLllVJdgKHATGfVQQghRNE4OylvElAOOA98Czxux7TZz4u9Vo4h9XQsqafj3Al1BKmnozm0nkpr7cjyhBBClEGy3IcQQgibJFgIIYSwqUSChb1rRCnjLaXURevX20oplevxNkqpHUqp69Z/25RQPZ9XSu1VSiUppY4opZ6/5fGjSqlUpVSy9WtVCdRxmlIqM1cdkpVSDXM9Xlru5Ypb6pihlNqT6/Fiu5fW8p9QSm1XSqUrpb6yce6flFJnlVJXlVJfKKU8cz0WoJRaZ72fB5RSfZxdR6XUQ9af5TWl1Enr749brsfXK6XSct3LOEfVsZD1HK+Uyr7l5x6a6/Fiu5eFrOdnt9QxXSmVlOvx4r6fnkqp/1p/f5KUUr8ppQYUcL5j359aa6d/YQa352IS9bpikvOa53Heo0AcUBeoA+wDHrM+5gEcA/4EeAJPWb/3KIF6vgC0xeStBFvrMSrX40eBPiV8L6cBs/Ipo9Tcyzyetx74qzPupbX8CGAY8E/gqwLOCwPOYVYgqGSt599zPb4ZeA8zoSMSuAJUc3IdHwe6WX++dYAdwEu33Ns/loJ7OR74qYDHi+1eFqaeeTzvK+ALJ97P8tbf4wDMH/qDMHlqAc54fxbLi7LjBWcAjXMdm5n7heQ6vgmYmOv7h4Et1v/3A05hHaS3HjsO9Hd2PfN47kfAx7m+L5YPuELey2nkHyxK5b20/lJkAw2K+17mce3XbXzAzQHezPV9b+Cs9f+NgXTAN9fjP2L9Q8dZdczj/GeAJbm+L9YPt0Lcy/HkEyycdS8Lez+t7+kkoIez7+ct9dgNROZx3OHvz5LohirMGlHNrY/ldV5zYLe2vlKr3fmUU9z1zKGUUpi/5m6dEjxbKXVBKbVKKdW6hOo4WCl1SSkVq5R6PNfxUnkvgXHAj1rrI7ccL457WVh5vTdrKKWqWB87rLVOuuXxkl4HrTv/+76coZRKVEr9nLvrpwTcY63HQaXU1FzdZaX1XkYCF4CNtxx32v1UStXA/G7llX7g8PdnSQQLu9eIyuPcq4CP9QO5MOUUdz1zm4a5r1/mOjYG81dyfWAdsFIpVdHJdfweaApUAx4B/qqUeqAI5RR3PXMbh2nq51Zc97Kw8npvgnlNxX0/C00pNQEIAd7JdfhFoCGmi+pzYIlSqlEJVG8j0AKojvkQfgC4Me5X6u6l1UPAN7f8geW0+6mUcgdmA19rrQ/kcYrD358lESzsXiMqj3P9gGTrD6gw5RR3PQEzUIb5gBuotU6/cVxr/bPWOlVrfV1rPQPTP9jNmXXUWu/TWp/WWmdrrTcBHwIjCltOcdfzBqVUV6AmMD/38WK8l4WV13sTzGsq7vtZKEqpYcDfgQFa65wF8LTWW7XWSVrrdK3118DPQLiz66e1Pqy1PqK1tmit9wCv4bz3ZqEppeoBPYBvch931v1USrlgunEzgCfyOc3h78+SCBaFWSMq1vpYXufFAq2srYwbWuVTTnHXE6XUH4CXgN5a65M2ytaAsnGOw+tYQB1K1b20egiI1lon2yjbUfeysPJ6b57TWl+0PtZQKeV7y+NOXwdNKdUf+Dcw2PpBXJCSupe3uvW9WSruZS7jgE1a68M2znP4/bT+jv4Xs4FcpNY6M59THf/+dOZgTK7BlO8ws2PKA13IfwbPY8B+TLOutvXF3DobagpmBs8TOH4Gj731HAOcBZrm8Zi/9bkegBemeX0BqOLkOg7FzIpQQHvMgPZDpe1eWs8th2kx9HLmvbRew81a9gzMX29egFse5/W3/sybWe/rD/x+tskWTJePFzAcx86GsreOvYCLQPc8HquImTHjZS1vDJACBJfAvRwA1LD+vwmwF3jVGfeyMPXMdX4c8Adn30/rdT6z3g8fG+c5/P3psBdRyBdcGVhovZnHgdHW490w3Uw3zlPA28Al69fb/H7Gzj2Y6YCpwK/APSVUzyNAJqZ5d+PrM+tjzTGDxSnWX9y1QEgJ1PFb6/WTgQPAU7eUUyrupfXYA5hgpW45Xqz30nqNaZi/CHN/TcMEqmTAP9e5z2CmJ17DjFF55nosADM7JhXz4eKwGVz21hEzppN1y/tyhfWxasA2TNfDFcyHR9+SuJeYD61z1p/rYUw3lLsz7mURfuadrPX0vaUMZ9zP+ta6pd3yMx3jjPenrA0lhBDCJlnuQwghhE0SLIQQQtgkwUIIIYRNEiyEEELYJMFCCCGETRIshBBC2CTBQgghhE0SLIQQQtgkwUIIIYRNEiyEEELYJMFCCDsppcpZ97I+nns/Y+tj/7HuIz2qpOonRHGSYCGEnbTWqcCrQD1g0o3jSqkZmC1/n9Raf1dC1ROiWMlCgkIUglLKFbMFZXXMrmh/BN7HLKn9WknWTYjiJMFCiEJSSg0ClmCWSO8FfKK1fqpkayVE8ZJgIUQRKKV2AG0xmzqN1vKLJMo4GbMQopCUUvcDbazfJkmgEHcDaVkIUQhKqX6YLqglmN0R7wNaaq33l2jFhChmEiyEsJNSqgNmnOIXzL7RdTF7xC/XWg8ryboJUdykG0oIOyilmgLLgIPAMK11utY6AfgvMFQp1aVEKyhEMZOWhRA2KKcyoJMAAABmSURBVKX8gZ+BDKCz1vpcrsdqAQnAb1prCRiizJJgIYQQwibphhJCCGGTBAshhBA2SbAQQghhkwQLIYQQNkmwEEIIYZMECyGEEDZJsBBCCGGTBAshhBA2SbAQQghhkwQLIYQQNv0/Hrqf/W0pQZYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#now plot the confidence bands\n",
    "plt.plot(X, y, \"b.\")\n",
    "plt.plot(X_new, y_predict, \"r-\")\n",
    "plt.plot(X_new, y_mu_upper, \"g--\")\n",
    "plt.plot(X_new, y_mu_lower, \"g--\")\n",
    "plt.xlabel(\"$x$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([0, 2, 0, 15])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X: [1.] y pred [7.34876878] y down: [6.99632288] y up: [7.70121469]\n"
     ]
    }
   ],
   "source": [
    "print('X:',X_new[100],'y pred',y_predict[100],'y down:',y_mu_lower[100],'y up:',y_mu_upper[100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The red line is a very good fit of the true slope (4) and offset (3), so the green bars are very small. \n",
    "\n",
    "The green bars are a prediction of the mean response at $x$, not the 95% of where a new prediction will land.\n",
    "\n",
    "The claim is: 95% of the time, CI covers the most optimal regression $\\mu = mx + b$ at that point $x$.\n",
    "\n",
    "## Confidence Interval of an individual observations, y\n",
    "\n",
    "If you want to draw the 95% bands of where an individual prediction will land, rather than the mean of the predictions at a point x, that is a different interval to calculate.\n",
    "\n",
    "For that, you want to predict, $ y_p = b_p + m_p x + e $\n",
    "\n",
    "What's different is now you want to take into account the random noise fluctuations, $e$, as well. $e$ is a gaussian centered at 0 with variance $\\sigma^2$ (but we don't know the true value of $\\sigma$).\n",
    "\n",
    "the $y_p$ random variable now has an extra term. But this term is independet of the other two. We can make a new covariance matrix, which is the original 2x2 matrix, but with a third row which is just $\\sigma^2$ on the diagonal and no covariant terms.\n",
    "\n",
    "$V[Y|X=x] = V[(b_p+m_px) + e] = V[b_p+m_px] + Var[e] + 2Cov[b_p+m_px, e] = \\frac{\\sigma^2}{n*\\hat{V}[x]}\\left(\\bar{x^2} + x^2 - 2x\\bar{x}\\right) + \\sigma^2 + 0$\n",
    "\n",
    "and\n",
    "\n",
    "$ E(y) = E(b) + xE(m) + E(e) = b + mx $\n",
    "\n",
    "So since we don't know $\\sigma^2$ exactly, we'll use $\\hat{\\sigma}^2$ and get the standard error for $y(x)$\n",
    "\n",
    "$\\hat{\\sigma}_{y|x} = \\sqrt{(\\hat{V_{00}}+ x^2\\hat{V_{11}}+2x\\hat{V_{01}}) + \\hat{\\sigma}^2} $\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SE2_new_y(x):\n",
    "    return SE2_y_p(x) + s2\n",
    "\n",
    "y_new_upper = np.array(list(map(lambda x,y: y + t2*math.sqrt(SE2_new_y(x)), X_new, y_predict)))\n",
    "y_new_lower = np.array(list(map(lambda x,y: y - t2*math.sqrt(SE2_new_y(x)), X_new, y_predict)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEVCAYAAAARjMm4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xdc1df5wPHPYYOg4hYVURyAcePeO26G2qRNY5NGE7P3TtQmaZo0TTrSNk1/bWaTVL2Ae0czHLFqokaGW+MeKKLse8/vjwNIDHoB7wKe9+vFS7n3y/f78AXuc894zlFaa4QQQojr8XJ3AEIIITyfJAshhBB2SbIQQghhlyQLIYQQdkmyEEIIYZckCyGEEHZJshBCCGGXJAshhBB2SbIQQghhl4+7A7CnUaNGOiIiwt1hCCFEtbJt27azWuvGjjqfxyeLiIgItm7d6u4whBCiWlFKHXbk+aQbSgghhF2SLIQQQtglyUIIIYRdkiyEEELYJclCCCGEXZIshBBC2CXJQgghhF2SLIQQQtglyUIIIYRdDk0WSqn7lVJblVL5Sqn3r3HMbKWUVkqNdOS1hRBCOI+jl/s4DrwMjAECr35SKRUJTAFOOPi6QgghnMihLQutdZLWOgU4d41D3gaeAgoceV0hhBDO5bIxC6XUVKBAa73MVdcUQgjhGC5ZdVYpFQz8FhhdweNnAjMBwsPDnRiZEEKIinBVy2Iu8JHW+mBFDtZav6u1jtVaxzZu7LDl2IUQQlSRq5LFCOBBpdRJpdRJoBUwTyn1lIuuL4QQ4gY4tBtKKeVTfE5vwFspFQAUYZKFb5lD/wc8Cix35PWFEEI4h6NbFs8DucDTwG3F/39ea31Oa32y5AOwAue11pccfH0hhBBO4NCWhdZ6DjCnAsdFOPK6QgghrrDarA4/p8fvwS2EEKJiLuRdYHrKdPZn7nf4uWVtKCGEqKay8rL4aMdHvLP1HQDq+dfjXM45RrUd5fBrKa21w0/qSLGxsXrr1q3uDkMIITxCVl4WizIWMT91Piv3r6TAWkDP5j3ZOvPHr5NKqW1a61hHXVe6oYQQwsNl5WVR178uSimeXvM072x7h1Z1W3Ffr/uYGjOVPi37OD0GSRZCCOGBLuRduNKC2LeSr+74ij4t+/BQ34e4vevt9GnZBy/lupEESRZCCOFBjl08xj1L72HlvpUU2goJrxfOg30epHEds5pFVKMot8QlyUIIIdzofO55FmYsxN/bn1s730rDoIYcunCIB/s8yNSYqfRu0RullLvDlGQhhBCuVpIg5qfOZ/X+1RTaChkdOZpbO99KgE8Au2btcneIPyHJQgghXCA7P5sQ/xAA7lp8F0lpSbSu15qH+jzE1E5T6RXWy80RXp9MnRVCCCfJzM1kYbppQaw5sIY9D+whon4E209sx2qzEhsW67QuJpk6K4QQHi7jbAYPr3yYNQfWUGQrIqJ+BA/3fRhfL7Oeao/mPdwcYeVJshBCiBuUmZtJSnoKzYKbMa79OBoENmBf5j4e6fsIU2OmOrUF4SqSLIQQogpKEsS83fNYe3AtRbYift7554xrP47GdRqz5/491T5BlCXJQgghKiinMIcg3yAAJn06iQ0/bKBN/TY82vdRpnaaSs/mPUuPrUmJAiRZCCHEdZ3LOUdyejLzU+ez4cgGjj92nLr+dXl1xKsE+QbRo3mPGpcYyiPJQgghyrH9xHaeWfsMaw+sxaqtRIZG8kDvByiwFgAwqPUgN0foWpIshBACOJtzluS0ZGIaxzAgfACBPoHsz9zPE/2fYFqnaXRr1q1WtCCuRZKFEKLWOnP5TGkX07qD67BqKw/1eYgB4QOIbhzN3gf21uoEUZYkCyFErVJ2kLr/v/uzL3Mf7Rq048kBTzKt0zS6Nu1aeqwkiiscmiyUUvcDvwI6A59qrX9V/Hhf4CWgJ2AF1gMPaq1POPL6QghRnuPZx0lKS2JB6gLSzqZx7NFj+Hj58Jexf6FZcDO6Nu0qicEOR7csjgMvA2OAwDKPhwLvAiuBIuBt4D3gZgdfXwghSq07uI4X17/IhiMb0Gg6Ne7ErNhZ5BXlEewXzM3t5CWoohyaLLTWSQBKqVigZZnHl5c9Tin1NvCFI68thBCHLhzCkmphZNuRdG1mupMu5l/kN8N+Q2J0ItGNo90cYfXlrjGLwcBuN11bCFGD7M/cz4LUBSxIW8DW42bR0ddHvk7XZl0ZGjGUHffscHOENYPLk4VSqgvwIjD5OsfMBGYChIeHuygyIUR1kZWXRb2AehRaC+n5bk+y8rPoFdaL10e+TmJMIm1D2wIyQO1ILk0WSql2wHLgIa31V9c6Tmv9LmaMg9jYWM9eQ10I4RJpZ9JYkLqA+anzKbQVknZfGr7evnyS+AmdGneidf3W7g6xRnNZslBKtQbWAC9prT9y1XWFENXbvN3zmPvFXFLPpKJQDAwfyJSYKVhtVry9vBnXfpy7Q6wVHD111qf4nN6At1IqADP7qSnwOfBXrfU7jrymEKLm0Frz3cnvWJC6gBk9ZxBRPwKAJnWacG/svcRHxxMWEubeIGspR7csngdml/n8NmAuoIG2wGylVOnzWutgB19fCFHNaK3Zenxr6SD1gfMH8FbedG7amYj6EUzrNI1pnaa5O8xaT7ZVFUK4nE3byMzNpFFQI87lnKPpG01RSjGizQimxkxlctRkGgU1cneYTrdpE6xfD0OHQr9+jj23bKsqhKiWbNrGxh82siB1AZY0C+0atGPd9HU0DGrI4lsX06dlHxoENnB3mC6zaROMGAEFBeDnB2vXOj5hOJIkCyGE0/35mz/zu69/x4lLJ/D39mdMuzFMi7nStTS2/Vg3Ruce69ebRGG1mn/Xr5dkIYSoRYpsRXxx6AsWpC7g1ZGvUj+gPn7efvRr1Y8p0VMY32E8df3ruiweZ3b13IihQ02LoqRlMXSouyO6PhmzEELcsEJrIesOrWP+7vkkpydzLvccQb5BLPv5MoZEDHFbXJ7e1SNjFkKIGq/AWsD53PM0DW7KgfMHGPPxGIL9gpnYYSJTYqZwc7ubS5cCdxdP7+rp18+z4rkeSRZCiArLK8pj1f5VLEhdwKKMRYxtP5ZPEz+lY6OOrLptFYNaDyLAJ8DdYZZyRVePp3ZzOZokCyFEhTy1+in+tvVvXCq4RGhAKAnRCdx6062lz4+KHOXG6MrXr5/penLWi7lHdnNdvAhLlzr8tJIshBA/cangEsv2LmP5vuW8O+FdfL19CQ0M5ZZOtzAlZgrD2wzH19vXKdd29Dt1Z3b1uLKb67r35exZWLQILBZYs8YE42CSLIQQgNn3YWH6QixpFlbuX0leUR6Ngxqz//x+ohpF8fTAp50eg0e+U78OV81oKve+tD4OycmQlARffGEyVkQEPPAAJCTAgAEOjUGShRC12NmcsxRaC2ke0pxvT3zL7Sm30yKkBTN6zCAxOpGB4QPx9vJ2WTyePCBd3jt7Z3dzlSi5L62sB5mSl0SrW5LgyEbzZFQUPP00JCZCt27gpGXZJVkIUcucyD5BcnoyljQLXxz6ggd6P8BbN7/FwPCBbP71Znq16IWX8nJLbJ5ae3C9Fo/TZzSlpfGLQ0mMtVnoxreg4VJAd3jpJZMgol2z+58kCyFqkUmfTmLJniVodGnX0s86/QwAby9v+rTs49b4XPVOvbJc0eIpbbkM0fQL/M6MPyQlQVoa4UDoTf1Y2+oNQu+Mp8eUto69eAVIshCihtp7bi+WNAvbTmxj3pR5KKXo2bwnvcJ6kRiTSEzjGHeHWC5PrD240RaPvUH7TRtsPDt8MxMLk2hGEuiD4OVlvuC++yA+npCwMII2wcr1kL/J9fdIkoUQNcjB8wf5cMeHWNIs7Dq9C4BeYb24kHeB0MBQZg+dbecMojw30uK5ZhdWURF8+SVYLNz0UTLrCk5QgC9rGEVqwvOM/8ckaNTI/nlcRJKFENWY1ppvT35LWEgYzYKbseXYFuZ+MZcB4QN4a8xbJEQnEF5P9rF3hKq2eMp2YZGfzw/vrKHf/yXBwoVw7hwEBVHQeyy/2pDAIut48vzrsfZxoNG1z+OOwX9JFkJUMzZt45uj32BJs2BJs3DowiF+N+J3PDXwKSZ1nMSxR4/RPKS5u8MUxYb3ucxO7xVMslkYb1tC3Q+zoW5dmDjRDFCPGUPDoCDu3gQd11+75eLuwX9ZSFCIaqTAWkDHtzty6MIhfL18GRU5isToRCZ3nEzDoIbuDk+UyMqCJUvMIPWKFZCbS05QI7JHxNF0VoLpT/Lzq/RpK1OwKAsJClFLFFoL+fzg51jSLGTmZrJg2gL8vP24o9sdtA1ty4QOE6gfUN/dYYoSZ86YrqWkJFNFXVgIYWHw619DQgJBgwYR5HNjL7nuHPyXZCGEh9lwZAP/3P5PFmYs5ELeBYL9gpnccTJWmxVvL29eHPKiu0Osthy+6N+xY6aK2mIxg9U2G7RpAw89ZLqYevc2s5pqAIcmC6XU/cCvgM7Ap1rrX5V5bgTwVyAc+Ab4ldb6sCOvL0R1dLngMsv3LWdEmxGEBoay9fhWFmYsZFLHSSRGJzI6crRHreRaXTlsNtGBA1dqIDZvNo/FxMBzz5llNrp2dVoVtTs5umVxHHgZGAMEljyolGoEJAF3AYuBl4D/An0dfH0hqoWsvCyW7FmCJc3Cin0ryC3K5cO4D/ll118yo+cMZvWahZ935fu0xbVVeTaR1pCaapKDxQI7dpjHe/SAV14xCSIqyomRewaHJgutdRKAUioWaFnmqQRgt9Z6fvHzc4CzSqkorXW6I2MQwlPZtA0v5cXJSydp/cfWFFgLaB7cnDu730lidCKDWg8CcPuGQTVVpWYTaQ3bt19pQWRkmNZC//7w5psQH28W7XMST9wjw1VjFp2AHSWfaK0vK6X2Fz8uyULUWCcvnSQ5zazD1KROEz5J/IRmwc14adhLDAwfSN+Wfd22DlNtY7ewzmYzr9IlCeLwYfD2hmHDzBhEXBw0d/6UZHcX312Lq5JFMHDmqseygJDyDlZKzQRmAoSHS0GRqH4+2fUJ72x9h6+PfI1G06FhB4ZGDC19/skBT7ovuFrsJ7OJCgvN8t5JSWag+uRJ8wo9ejTMmWNqIRq6dkqyu4vvrsVVyeISUPeqx+oC2eUdrLV+F3gXTJ2Fc0MT4sbtz9xPSnoKD/R5AD9vP3ad2sWFvAvMHjKbxJhEOjXuhKqBg57VUl4erF5tEsSiRZCZCUFBMG6cmcE0bpwpmnMTdxffXYurksVuYHrJJ0qpOkBk8eNCVDtaa1LPpJKUloQlzcKOU6aXtU/LPgwMH8hLw1/i1ZGvujlKUerSJVi+3CSIJUvM5/XqwaRJZoB6zBgIDLR/Hhfw1JV3HT111qf4nN6At1IqACgCkoHfK6USgaXAi8BOGdwW1YnWmpzCHOr41eF/x/9Hn//rg0LRr1U//jD6DyREJxBRPwIAHy8pYXK78+evVFGvXGlaFI0bw623mhbEsGFVqqJ2BU9cedfRv9HPA2WXtbwNmKu1nlOcKN4GPsbUWdzi4GsL4XBWm5Wvj3xNUloSyenJxEXF8eexfyY2LJZ/Tvwn49uPl3WYPMnp05CSYloQa9ealV1btIAZM0yCGDjQDFqLSnP01Nk5wJxrPLcGqPmTkUWN8dTqp3jvu/c4k3OGAJ8AxkSOYVjEMAC8lBd39bjLzREKAH744UoV9ddfm1lNkZHw6KMmQcTG1pgqaneStrIQmCrqFftW8OXhL/njzX9EKUVeUR4j244kITqBm9vdTLBfsLvDFCX27btSJLdli3nsppvghRfMGETnzjWyitqdZNVZUWtdyLvA4ozFJKUnsWLfCvKK8mgY2JAd9+ygRd0W7g5PlKU1fP+9SRBJSbBzp3k8Nta0HhISoEMH98boYWTVWSFuwMlLJ/Hx8qFRUCPWHljL7Sm30yKkBTN6zCA+Kp5BrQfJ4LSn0Bq2br1SJLd3r2ktDBwIb71lEoTUYbmM/FWIGu/QhUMkpyWTlJ7EhiMbeHn4yzw76FnGth/LN3d9Q2xYrFRRewqrFTZuvJIgfvgBfHzMzKXHHoPJk6FZM4dcyhOX1PBkkixEjWW1WRn43kA2HzUrg3Zp2sUUyUUnAmYNpt4terszRAGminrdOpMcUlLg1Cnw9ze1Dy+9ZKqoGzRw6CU9dUkNTybJQtQIWmu2n9hOUloSR7OP8kHcB3h7eTOg1QASoxOJj4onskGku8P0aC59p52ba6qoLRZTRX3hAtSpA+PHm+6lceMgpNzVgBzCU5fU8GSSLES1tuPkDt7/7n2S0pM4knUEb+XN8DbDKbQW4uvtyxuj33B3iNWCS95pZ2fDsmWmBbF0KVy+DPXrmyrqxEQYNarCVdQ3mtg8dUkNTybJQlQrBdYC1h1cR+8WvQkNDOWrI1/x961/Z3TkaOYOncvEDhNlL+oqcNo77cxMWLzYtCBWrYL8fGjSBG67zbQghg0DX99KndIRic1Tl9TwZJIshMfLKcxh5b6VJKUnsThjMVn5Wfx70r+5o/sdTO86neldpxPi77wui9rAoe+0T568UkW9bp2pom7VCu65x7Qg+ve/oSpqRyU2T1xSw5NJshAeSWuNUoqzOWdp/cfW5BTm0CCwAQnRCSREJzCy7UgASRIOcsPvtI8cuVIkt2GDmfbavj08/rhpQcTGOqxITrqQ3EOK8oTHOH35NAvTF5KUnkRoQCifJH4CwCtfvkK/Vv0Y3Hqw1EB4kj17riSIkr/RLl1MckhMhE6dnFZFLdNe7ZOiPFHjfLrrU97ZZjYKsmkbbUPbclvn20qff27wc26MTpTSGnbtulID8f335vHeveG118xWo+3buyQU6UJyPUkWwuUyzmaQkp7Cw30fxt/Hn91ndpOZm8nzg54nMSaRzk06y0ZBnkJrs/ZSyTIb+/aZ1sKgQfCnP5kE0aqVu6MULiDJQjid1prvTn5HUloSSelJpJ5JBaB/q/4Maj2IOUPn8PLwl90cpShltZrVWy0Ws5rr0aOminrECHjiCVNF3bSpu6MU5cgtzGVB6gL2nNvj8HNXKFkopd4B7gZaaK2PX/VcR2AX8Het9UMOj1BUS1ablUsFl6gXUK90oyAv5cWQ1kOYFTuLuKg4WtZtCchGQR6hoAA+//xKFfWZMxAQYKqof/tbmDABQkPdHaUox+ELhzmTc4bYsFis2sqMxTNoUqeJw69T0b/STZhk0RtIueq5t4CLXGMfC1F75Bfls+bAGpLTk1mUsYhpnabx9ri3iQ2L5YO4DxjXfhyNghq5O0xRIjfX7CBnsZhaiKwsCA42iSEhAcaONZ8Lj6K1ZtfpXaSkp5CSnsK3J7+lV1gvtszYQrBfMLtm7SKyQSTejzp2k6eKJovNxf/+KFkopcYDY4H7tNbnHRqZqFbuW3ofH+78kEsFlwjxC2F8h/GMaz8OMBsF3d71djdHKAC4eNFUUVss5t+cHNNiiI83M5hGjjQtCuFRbNpWutjlHQvv4IMdH6BQ9G/Vn9+P+j2TO04uPbZ9Q+dMMqhQstBaZyilMjHJAgCllC/wJvA98A+nRCc80qlLp1iUsYgNP2zgvcnvoZQi2C+YW2+6lfioeIa3GY6/j7+7wxQlzp0z6y8lJZkq6oICs3Lr9OmmBTFkSKWrqIXz5RTmsHr/alIyUli6Zyk7Z+2kWXAzbrnpFgaFD2JChwk0DXbd2FFlOos3AwOUUkqb4oyHgA7ASK211SnRCY+waRMsXHeMnDb/5du8ZDYc2YBG0za0Lacun6JZcDNeG/Wau8MUZZ04YcYeLBZTkGC1QuvWcN99pgXRt6/sRe2h0s+m88zaZ1i5byW5RbnUD6jPhA4TyCnMAeDmdje7Ja7KJotxQMfiVsYLQIrWem1FT6CUigD+BvQD8oEFwMNa66JKxCFcoKRf9FBqI24ZH0Z+22+wTX2MdsFdeHHIi8RHxdOlaReZ4upJDh26MsV140Yz7bVDB3jySZMgevSQrUY90MHzB1mYsZCODTsytv1Y6vjWYdvxbfy6+6+Ji4pjcOvB+Hq7v+VXmWSxqfjf3sBgwB94rJLX+xtwGmgO1AdWA/cCf67keYQTWG1WNh3dRHJaMikZKRw4f4CR3i9RUPA8toyxeP1lP3c+1pZnhro7UlEqPf1KFfX27eaxrl1h7lzTxRQT47IEIVXVFffdye9KB6h3nNoBwH297mNs+7G0qteKww8f9rg3YpVJFt8ANuDXwEDg91rrA5W8Xhvgba11HnBSKbUC6FTJcwgHKlmDqchWROSfIzmSdQQ/bz9GtBnB0wOeJix7Ehv8oKAgEL+ctrIOj7tpDTt2XKmiTjU1K/TtC6+/bhJEpOv37ZDNhK6vyFZE+tl0bmpyE2AGqXec3MHA8IG8MeoNJkdNpl2DdqXHe1qigEokC611tlIqFdOqOAm8UoXr/Qm4RSm1HgjFzKR6oQrnETcgOz+b5fuWk5yeTGZuJitvW4mPlw/3xt5L6/qtGdd+HHX965YeL0s5u5nNBt98c6WL6cAB8PKCwYNh1iyIi4OWLd0aomwm9FOXCy6zav8qUjJSWLJnCXlFeZx94iyBvoH8e9K/aVG3hVPqIZylstVQW4CbgGe01tlVuN4XwAxMXYY38AE/rdtAKTUTmAkQLhuyO8zKfSv5y5a/sObAGvKt+TQOakxcVBxWmxVvL2+eGvhUuV/njnV4anOXxqZN8MXaIibW/4pOGUmmivrYMTNjaeRIeOYZU0XduLG7Qy0lK8H+2Ec7PmLmkpnkFeURGhDKhA4TiIuKw9vLTCro3ry7myOsAq11hT4AX2A/8D+KV6utzAfgBRwBnsOMdzQEFgKvX+/revbsqUXVHDx/UL+16S19+tJprbXWf93yV936rdb64eUP6y8PfamLrEWlx27cqPVvf2v+dbeNG7UODNTa29v86wkxuURenk79wzL9nvev9WkaaQ26yD9Q6/h4rT/+WOvz5x1yGWf9rD3pd8iV9mfu129ufFMPfm+wXrF3hdZa6x0nd+gHlz2oPz/wuS4oKnBLXMBWXcnX6et9VKZl8ThmzOEXxYFUVgOgFWbMIh/IV0q9B7wMPFmF84mraK35/vT3JKcnk5yezHcnvwMgLCSMaZ2mMaPHDGbFzvpJf6in9DeXtCaOHPGMLg2XtG5ycmDFCtO9tHgx0Rcv0oIQFjORFK8E+jxzM4/PruOwyznzZ12bVoLNKczhd1//jpT0FHad3gVAl6ZdKLAWlP7/T2P/5M4QHe66yUIp1QAYA3QBngDe1Fpvvt7XXIvW+qxS6iAwSyn1BhAMTAd2VOV8wrBpG+dzz9MwqCHHso/R5Z0uKBT9WvXj96N+T1xUXOnA2bWm33lCf3PZFzEfnyslAO7q0nBqAs3KMntQWyywfLlZdqNhQ5gyhbToBPq9MJJLhf74+cGjox10zWKe8LOujgqthXx15CvO5ZxjaqepBPgE8N537xEZGslbY95icsfJtAlt4+4wncpey2IM8AlmuutbwNM3eL0E4I/AU4AVWAc8coPnrHVK9qFOTk9mYcZC+rToQ8otKbSs25L/Tvkvg1sPpllwswqfz1X9zdd7p172RQxgxgwID3feu3p7rQaHv6iePQsLF5oWxJo15qTNm8Mdd5gaiMGDwceHaGD5AOe1aGRsoeIuFVxi5b6VpQPUF/Iu0L5Be6Z2moqX8mLvA3sJ8KlFS6M4sk/LGR8yZvFjz699Xtd7tZ5mDrrOK3X0lHlT9Pzd82/4vM7ub7Y3DuHKcYqKXMsh8Rw9qvXbb2s9bJjWXl5ag9YREVo/9pjWGzZobbXe8PdSFbV1bKEiTl06pW02m9Za6xmLZmjmoBu81kBPT56uk9OS9eWCy26OsOJw45iFcLGzOWdZnLGYZfuW8VH8RwT4BBAaGEpCdALxUfGMbDuSQN9Ah1zL2f3N9t6p3/Ae0A6M5YbiOXjwSg3EpuI61qgoM4MpMRG6dXN7FXVtGluoiP2Z+0lJTyE5PZmNP2xk+93b6dasGw/0foCfd/45A8MHyjL6yOZHHufM5TN8susTktOT+erIV9i0jfB64Rw4f4CYxjE82u9Rd4dYJRXp/nDVi1hFu2IqHE9a2pUE8e235rHu3eHll02RXHS0gyIXjpR+Np1p86eVDlB3bdqV2UNmly6j37lpZ3eG53EkWbiZ1prUM6n4+/jTrkE7Dl44yMMrH6ZT4048O/BZ4qPj6d6su0dWdFaGK1sOTo9Fa5MUSpbZSE83j/fvD2+8YRJEm5o92FndWG1WNvywgeS0ZKIaRXF37N2E1wunaXBT7ux+Z60YoL5RSldpFqzrxMbG6q1bt7o7DIeyaRubj24uXRtmb+Ze7ul5D3+f8Hds2sb+zP1OW5NeVJHNBps3X2lBHDpkpmwNGWK6l+LiICzM3VGKq6zct5L5qfNZlLGIMzln8PP24/5e9/OHMX9wd2hOp5TaprWOddT5pGXhIrp4DSaAHv/owY5TO/D18mV4m+E80vcRJkeZzUu8lJckCk9RVARffGGSQ3KyWfbb1xdGjYIXXoBJk6CR7PznSbLysth0dFPpMt5/3vJnvjr8FeM7jCc+Kp6x7cYS4h/i5iirJ0kWTnQ+9zzL9i4jJSOF1DOpfD/re5RS3NfrPkL8Qxjbbiz1Auq5O0xRVn6+mdpqsZiprpmZEBRkthhNSIDx46Ge/Mw8yalLp1iYsZDk9GTWHlhLoa2QY48eIywkjHcnvEujoEayGZcDSLJwgs8Pfs6rX7/K+kPrKbIV0Sy4GZM7TuZy4WWC/YKZ0XOGu0MUZV2+bIrjkpJgyRLIzoa6dWHiRNPFNGaMSRjCY5S01C2pFqbOn1q6GdeDfR4kPiqepnXMDnIt6rZwc6Q1hySLG6S1ZveZ3aSkpzAlZgpRjaLIzs/mh6wfeLzf48RFxdGrRa/S/XOFh7hwwSQGi8Ust5GXZ7qUfvYz04IYMcJMlaqC2rwIorNordl5amfpUjYP9H6Au3rcRf9W/Zk9ZDbx0fF0btK52k8E8WSSLKqgZJOgkgHq/ef3A9CkThOiGkUxqeOk0jEI4UHOnDFdSxaLmQ5VWGgGpe+6y7QgBg40a43cgIouEyIJpWJw5kENAAAgAElEQVSsNitPrn6S5PRkDl44iEIxIHxA6fTW5iHNmT10tpujrB0kWVRQbmEux7OPE9kgktyiXEZ+OBKNZnib4TzR/wkmdpxIWIiZDSPvbjzI0aNmcDopCb780sxqatsWHn7YtCB69zZ7QzhIRQr+PGXhRk+UX5TP2oNrOXThEPf2uhdvL282Ht1IdONonh30LJM6TqpWe0DUJJIsriMzN5Ole5aSkpHCyn0riWkcw5YZWwj2C2bVL1fRrVm3H20SJDzE/v1XaiC++cY8FhMDzz1nWhBdujilinrTJrNirr1FEGUxvx+7mH+R5XvNZlzL9i4juyCbpnWaMrPnTHy8fNhw5wbpxvUAkiyu4dm1z/L6htexaithIWHc3vV24qLiSp8f3HqwG6MTP6K12V60pAZiR/FCxj17wiuvmBZEVJRTQ7h61dwZM+D228tPArKYn5nBVD+gPv4+/vxh4x/4zZe/oUmdJtxy0y3ER8UzvM3w0iU2JFF4hlqfLHTxHhAp6SkszFjIkp8voVlwM3o278mTA54kLiqO2LBY+YWtJKf3yWsN27ZdaUHs2WNaC/37w5tvQnw8REQ44cLlu3rV3PDwa3/fnlTN7koHzx8sHaDecGQDi29dzPgO47mz+52MihxFv5b9SneSE56n1iaLoxeP8uamN1mYsZAD5w+U7gFx+vJpmgU3IzEmkcSYRHeHWS05rU/eajUnL2lBlPT5DBtmxiDi4syy325Q2dZCbVrM7+Slk4z5eAw7T+0ErqzBFNM4BoDW9VvTun5rd4YoKqDWJIvcwlxWH1hNaEAog1oPQmvN37f+nRFtRvD0gKeZ2HFipfaAENfm0D75wkJTRW2xQEoKnDxpXo1Hj4a5c00tRMOGDoy+ampra+FqJTMFk9OSqR9QnxeGvECTOk2IDI1ketfpxEXF0Ta0rbvDFFVQo5PFuZxzLNmzhIUZC1m5fyU5hTlM6zSNQa0H0apeK849eY4gXym2crQb7pPPy4PVq03rYdGiK1XU48aZAepx40zRnIepTa2Fq607uI5Pdn3Coj2LOH35NH7efvyi8y8AM+aQ9LMkN0coblSNSxbncs7RMMi80xz98Wi2n9hOi5AW3NHtDuKi4hjSekjpsZIonKNK77IvXYJly0yCWLrUfF6vnll/KSHBVFEHOmbvDnHjsvOzWX1gNfFR8Sil+PT7T/nv7v8yrv04swZT+7EyU7CGqfarzmqt2XFqBwvTF5KSkcLec3s588QZAn0D+fzg59T1r0vP5j2l9sETnT8PixebBLFypWlRNG5sxh4SE81YRBWrqIXjnb58mkUZi0hOT2bNgTUUWAvYOmMrPcN6cjbnLCF+IbIGkweRVWfLWLZ3GfcuvZfDWYdLKzvnDJ1Dka0IgOFthrs5QvETp05dqaL+/HOzsmvLljBzpmlBDBx4pVBBuJ1N2/BSXnx5+EuGvj8UjaZN/Tbc3+t+4qPj6dasG0BpRbWouVyeLJRStwCzgXDgJPArrfVX9r4upzCHVftXkZKewu1db2d4m+E0C25G56adeWHwC0zoMIGmwU2dHb6oih9+MK2HpCT46isz7TUyEh591LQgYmMdWkUtbkzamTSS0pJISk9icsfJvDjkRXqF9WL2kNnERcXRpWkXaanXQi5NFkqpUcBrwM+ALYDdeY7ncs4R91kcq/avIrcol/oB9UsL4no078HiWxc7NWZRRXv3XqmB+N//zGM33QQvvmhaEJ07u30vavFjL33xEp98/wnpZ83Of31b9qVNfbN7XKBvoKzBVMu5umUxF/iN1npz8efH7H3B0eyjWE9YuavHXcRFxTEofBC+3r7OjbKWuqFCOq3h+++vJIhdZl9jYmPh1VdNgujQwcERi6qy2qxs/GEjW45t4bH+jwHw7clvaRHSggd6P8DkjpNleW/xIy4b4FZKeQO5wIvAXUAAkAI8obXOvdbXdeneRe/YvkOavU5WpUI6rU2roaSLae9e01oYONAkh4QEU8osPEKBtYB1B9eRlJbE/O9TOF9wGj+vAE48fowGgQ1KxydEzVCdB7ibAr7AFGAQUAgsBJ4Hnit7oFJqJjATIDw83GmJQpaJvqLChXRWK2zYYFoPyclmPMLHx8xceuwxmDwZmklxo6fIKcwBzDTx9797n7uX3E2QdzD5349HpSbgdWQsGSND6NdP1mASdmitXfIBhAIamF7msUTg2+t9Xc+ePbUzbNyodWCg1t7e5t+NG51ymWrjuvejoEDrlSu1njlT6yZNtAat/f21njRJ6/ff1/rcueue97e/lfvrSudzz+uPdnyk4z+L14EvB+p/bP2H1lrrM5fP6CUZS/TcV3K1t7f5MXp7m5+PqHmArdqBr+Eua1lorc8rpY4WJwy3k2Wif+wnhXTdcmHhqitV1BcuQJ06Zg/qxESzJ3XI9Te+l30bXKvAWsDkzyaz5sAaimxFhIWEcWf3O4kNMz0RjYIaMb7DeBoMg9/V8lVvReW5eoD7PeABpdQKTDfUw8ASF8cAyDLR5el3Uzb9jiyDtyymmvryZQgNNV1LCQkwalSlqqglITvX4QuHSU5P5szlM7wy4hX8vP0I8Qvhkb6PkBCdQO8WvcvtWpJ1rERVuDpZvAQ0AvYAecA84BUXxwDIH0ypzEzTckhKglWrID8fmjSB224zLYihQ8G3arPPJCE73r7MfczbPY+ktCS2ndgGQO8WvXlJv4SX8mLe1HkVOk9tXsdKVE21X+5DVMHJk2YFV4sF1q0zb/1btTKth8REsyeEg6qoZRLBjdFa8+3Jb4lpHEOATwBz189lzhdz6NuyLwlRCcRHx9OuQTt3hyk8kKNnQ0myqC0OH74yxXXDBjPttX17kxwSEkw9hExP/hF3JbqSGoiktCSS05M5nHWYRbcsYmLHiZy+fJpCa6HUQAi7qvPUWeFqe/aY1oPFYnaVA7P/9OzZfBeZyPIjnRg6TNGvl3vD9ETXG5x3ZhI5knWEXv/sVbrM9+jI0cweMpsB4QMAaFKniWMvKEQFSbKoSbSGnTuvVFHv3m0e790bXnvNbDXavr3MUqqAaw3OO/Le5RTmsHLfSpLSk2ge3JzXR71Oq7qtiI+KZ1jEMFnmW3gUSRbVnc1mqqhLthrdv98syjdoEPzpTyZBtGr1oy+RWUr2XWtw3hH3bmH6Qj7c+SHL9y4ntyiXBoENuLPbnQAopXhnwjsO/E6EcAxJFi7g8G4Lq9Ws3loyBnHsmKmiHjECnnrKTHVtcu3uCpmlZN+1ZstV5d6dvnya5XuX88uuv8RLebFi3wo2H93Mnd3vJCE6gcGtB+PjJX+KwrPJALeTOazboqDA7P9gsZj9IM6cgYAAuPlmM0A9YYKpiahEXDJLqWoqcu+OZB0hOS2ZpPQkvj7yNTZtK90oKDs/mzp+dWR5DeFUMsBdzdxQt0VOjtlBLinJ7CiXlQXBwSYxJCSYKurg4CrFJfPsq+5a967IVoSPlw/rD61n2AfDAOjcpDPPD3qehOgEujTtAkCI//Ur34XwRJIsnKzS3RYXL5o9qC0WWL7cJIzQUDP2kJgII0eaFkUNVN1aO7p4S19LqgVLmoWpMVOZO2wufVv25fWRrxMXFUf7hu3dHaYQDiHJwskqVCl+7pyporZYYPVqk1maNYPp000LYsiQKldRVxfVbYbW7HWz+XjXxxw4fwAv5cXg1oOJaRwDQIBPAE8MeKL02OqWBIUojyQLFyi32+LECUhO5sJ7SdTdvh4vmxVat4b77jMtiL59a9Ve1J48Q6ukSG7z0c2lSSD9XDodGnbgmYHPMLnjZBrXaVzu11a3JCjEtUiycKVDh67UQGzaBFpzSnXk7zzFEv8E3vikB/36184qak+boVVoLWT9ofVY0iykpKdw6vIpAnwCuKP7HTQKasRniZ9VaJ8VT06CQlSGJAtnS0+/UgOxfbt5rFs3mDuXd88lcu9forHaFN5FsP4L6NffveG6iycs7JhflI9VWwnyDeLDHR9y1+K7qONbh3Htx5EYnci49uNKB6cruiGXpyVBIaqqViQLl/YZaw3ffXelBZGWZh7v2xd+/3szUB0ZCUDnTeD3rryQlHDHDK2cwhyW712OJc3Ckj1LeG3ka8zqNYu4qDgaBjVkTOQYAn0rviz71TwhCQrhCDU+Wbikz9hmg2++udKCOHjQVFEPHgz33gtxcdCy5U++TF5I3KfIVsStlltZumcpuUW5NAxsyLRO0+jRvAcADYMaEhcV55BryTRlURPU+GThtD7joiL48kuTHJKT4fhxM2Np5Eh49llTRd24/EHPsuSFxDUyczNZmL6QoxeP8sKQF/Dx8qHIVsQd3e4gMSZRqqiFsKPG/3U4tM84P980BZKSTBX12bNm57ibbzYzmMaPh/r1HRS5uFGnL58mOS0ZS5qFzw9+jlVbad+gPc8MegYfLx+Sf5bs7hCFqDZqfLK44a6ey5dNFbXFAkuWmKK5kBCYONHUQNx8s9mbWniEoxeP0iioEQE+Afxj6z94cf2LtG/Qnif6P0FiTCI9m/es8OC0EOKKWrk2lN0B76wskxiSkkwVdW4uNGxoupYSE80giL+/Q2MSVbc/cz+WNFNFveXYFpKmJREfHc/x7OOcyznHTU1ukgQhah1ZG+oGXXPA+8yZK1XUa9ZAYSE0bw533mlaEIMHm5Vdhcc4c/kMoz4axY5TOwDo2bwnvx3+W3qG9QQgLCSMsJAwd4YoRI3hllc/pVR7YBewQGt9myuvXXbAu1H+MbJfTYZLSfDFF2ZWU0QEPPigaUH06WNmNQm3K7sOk5+3Hy8MeYFGQY1o16Adt3e9nYToBCLqR7g7TCFqLHe9Vf4r8D93XHhM+wNcVElMxkJf22ZYDERHwzPPmATRrZvsRe1Btp/Yzmfff4YlzVK6DlNidCJgCuMWTFvg5giFqB1cniyUUrcAF4CNQDuXXDQ1tbQGosd339EDOBnWnSMTXib84QSTLIRHsNqsbDq6if6t+uOlvPjX9n/xz+3/ZETbEXbXYRJCOI9LB7iVUnWBrcAI4NdAO3vdUFUa4NYavv32SpFcerp5vH9/M/6QkABt2lTlW5AVRJ2gvHWYNt65kX6t+nEi+wSBvoHUD5ApyUJURnUf4H4J+JfW+ofrzU5RSs0EZgKEh4dX7Mw2m3klL9lq9NAhs2rrkCHwwAOmijrsxgY7ZQVRx9t5aifDPhhGZm4mdXzrML7DeBKiEujctDMAzUOauzlCIQS4MFkopboBI4Hu9o7VWr8LvAumZXHNAwsLTRW1xWKqqE+eNFXUo0bBCy/ApEnQqJGjvgVZQfQGXS64zIp9K7CkWejatCtPDXyKjg07EtcxjkkdJzE6cvQNrcMkhHAeV7YshgIRwJHiVkUw4K2UitFa96jwWfLzzQZBJVXUmZkQFGS2GE1IMFXU9eo5I35ZQbSKktKS+GTXJyzbu6x0HaaoRlEA+Pv486/J/3JzhEIIe1yZLN4FPivz+eOY5DHL7ldeugQrVpgWxNKlkJ0NdeuaKurERBgzxiQMJ5OF/yomKy+LLw9/ycSOEwH4z67/sPGHjbIOkxDVmNsquJVSc6jIAHdoqN6alwd5eaZLKS7OtCBGjDBv75FBZ09wPvc8izIWsSBtAav2r6LAWsCBBw/QJrQN53LOUT+gPt5etWfnPyHcrboPcJfSWs+p0IGXL8Pdd5sWxMCBP6milkFn91uxbwUTP51Ika2I8Hrh3N/rfqbETKF1/daAWe5bCFG9eX5fQJcu8Je/XPNpGXR2rZKVXBekLWBK9BTujr2bXmG9eKTvI0yNmUpsWKyswyREDeT5ycIOGXR2Pq0172x9h3mp8/jy8JfYtI12Ddrh6+0LmJbD66Ned3OUQghnqvbJQgadnePoxaNsPb6VuKg4lFJ8vOtjMnMzeW7Qc0yJmULnJp2lBSFELVIrlygX5Tt84TCWNAsLUhew6egmfL18OfvkWer61+Vi/kXq+td1d4hCiAqqMQPcwjNorVFK8X/b/48Zi2cA0L1Zd14Z/gqJ0YmlCUIShRC1mySLWmjPuT0sSF3AgtQFPDfoORJjEhkWMYzXRr5GYnQikQ0i3R2iEMLDSLKoJQqsBbz29WvMT53PrtO7AOjbsi8BPgEARDaI5MkBT7ozRCGEB5NkUUNprfn+9PfsP7+fuKg4fL18+WjnRzQNbsofx/yRhOgEWtVr5e4whRDVhCSLGkRrzXcnvzNdTGkL2HNuD42DGjOxw0S8vbzZcc8OWahPCFElkiyquZLZbEopXlz3Ii9/9TLeypthbYbxaN9HiYuKK11mQxKFEKKqJFlUQzZtY8uxLaWD1P9J+A8DwgcwJWYKEfUjmBw1mUZBjluaXQghakyyqA2LCZ7PPc/cL+ZiSbNw9OJRfL18GR05unQF167NutK1WVc3RymEqIlqRLKoqYsJWm1WvjryFdn52UzsOJE6fnX47+7/0qdFH14d8SoTOkyQ7UaFEC5RI5JFTVpMsMhWxPpD61mQuoDk9GROXz5N16ZdmdhxIn7efhx5+EjpmkxCCOEqNSJZVPfFBItsRaVdSXctuosPdnxQuh/1lOgpjGs/rvRYSRRCCHeoEcmiOi4mWGgtZO3BtczbPY+FGQvZctcWIhtEcnfPu5nccTJj2o0hyNf5u/8JIURF1IhkASZBVIckcfTiUV5c9yIp6SmczztPXf+6TOo4Cau2AtCvVTX4JoQQtU6NSRaeqsBawNoDa/H19mVk25EE+QaxKGMR4zuMZ1rMNEZHjsbfx9/dYQohxHVJsnCCkgQxL3UeKekpXMi7wJjIMYxsO5IGgQ04+fjJ0jEKIYSoDlz2iqWU8gf+BowEGgD7gGe11stdFYMzWW3W0krpuM/iWL5vOXX96zK542SmdZrGqLajSo+VRCGEqG5c+arlA/wADAGOAOOAeUqpzlrrQy6Mw2EKrAWs3r+a+anzWb5vOen3pRMaGMojfR/h3l73MqrtKOliEkLUCC5LFlrry8CcMg8tUUodBHoCh270/K6s4N57bi+vfPUKKekpZOVnUc+/HnFRcVwquERoYCijIkfZP4kQQlQjbusPUUo1BToAu2/0XM6u4M4vymf1gdU0DmpMn5Z9UEqxMGMhcVFxTOs0jZFtR+Ln7ee4CwohhIdxS7JQSvkC/wE+0Fqnl/P8TGAmQHh4uN3zOaOCO78on1X7VzE/dT6LMhaRlZ/FL7v8kj4t+9CuQTtOP35aCuSEELWGy5OFUsoL+AgoAO4v7xit9bvAuwCxsbHa3jkdVcFt0za8lBcAA/49gG0ntlE/oD7x0fFMi5nGiLYjSo+VRCGEqE1cmiyUUgr4F9AUGKe1LnTEeW+kgjuvKK+0BfH1ka/JuD8DP28/nh30LIE+gYxoO0K6mIQQtZ6rWxZ/B6KBkVrrXEeeuLIV3DtP7eT3G3/PooxFXMy/SGhAKHFRcVzMv0ijoEYkRCc4MjwhhKjWXFln0Rq4G8gHTppGBgB3a63/4+zr5xXlsXLfSto3bE9M4xiy8rJYumcpU6KnMLXTVIa3GS4tCCGEuAZXTp09DCi7BzpQSYKYlzqPxRmLyS7I5rF+j/HG6DcYED6AU4+fkrEHIYSogBpXSqy1RimF1pqot6M4nHWYBoENmNZpGtM6TWNYxDAAvJQXXt5ebo5WCCGqhxqRLHILc1mxbwXzU+eTfjadbTO3oZRi7tC5NA9pzrCIYdKCEEKIG1Ctk8WWY1v44+Y/snjPYi4VXKJhYEMSohPILcolyDeI6d2muztEIYSoEapVssgtzGX5vuXEhsUSXi+cI1lHWLV/FbfedCtTY6YyNGKotCCEEMIJPD5Z2LQNS6qF+anzWbJnCZcLL/PayNd4csCTTO44mckdJ0uCEEIIJ1Na2y2Qdiuflj7aOsNqah+iEpjaybQgZJlvIYS4NqXUNq11rKPO5/GvuGEhYbz3y/cYEjFEEoQQQriJx7/6NqnT5EdrMgkhhHA9KTQQQghhlyQLIYQQdkmyEEIIYZckCyGEEHZJshBCCGGXJAshhBB2SbIQQghhlyQLIYQQdkmyEEIIYZckCyGEEHZJshBCCGGXS5OFUqqBUipZKXVZKXVYKfVzV15fCCFE1bh6IcG/AgVAU6AbsFQptUNrvdvFcQghhKgEl7UslFJ1gETgBa31Ja3118Ai4JeuikEIIUTVuLIbqgNg1VrvKfPYDqCTC2MQQghRBa7shgoGsq56LAsIufpApdRMYGbxp/lKqe+dHJsjNALOujuICpA4Has6xFkdYgSJ09E6OvJkrkwWl4C6Vz1WF8i++kCt9bvAuwBKqa2O3BrQWSROx5I4Hac6xAgSp6MppbY68nyu7IbaA/gopdqXeawrIIPbQgjh4VyWLLTWl4Ek4DdKqTpKqQHAZOAjV8UghBCialxdlHcvEAicBj4FZlVg2uy7To/KMSROx5I4Hac6xAgSp6M5NE6ltXbk+YQQQtRAstyHEEIIuyRZCCGEsMstyaKia0Qp4zWl1Lnij9eVUqrM892UUtuUUjnF/3ZzU5xPKKW+V0plK6UOKqWeuOr5Q0qpXKXUpeKPVW6IcY5SqrBMDJeUUm3LPO8p93L5VTEWKKV2lXneafey+Pz3K6W2KqXylVLv2zn2EaXUSaVUllLq30op/zLPRSil1hXfz3Sl1EhXx6iUml78s7yolDpa/PfjU+b59UqpvDL3MsNRMVYyzl8ppaxX/dyHlnneafeyknG+c1WM+Uqp7DLPO/t++iul/lX895OtlPpWKTX2Osc79vdTa+3yD8zg9n8xhXoDMcV5nco57m4gA2gJtABSgXuKn/MDDgOPAP7Ag8Wf+7khzieBHpi6lY7FcdxS5vlDwEg338s5wMfXOIfH3Mtyvm498KIr7mXx+ROAOODvwPvXOW4McAqzAkFocZy/K/P8JuBNzISOROAC0NjFMc4CBhX/fFsA24Cnr7q3d3nAvfwV8PV1nnfavaxMnOV83fvAv114P+sU/x1HYN7oT8DUqUW44vfTKd9UBb7hAqBDmcc+KvuNlHl8IzCzzOe/BjYX/380cIziQfrix44AN7s6znK+9s/AX8p87pQXuEreyzlcO1l45L0s/qOwAm2cfS/LufbLdl7gPgF+W+bzEcDJ4v93APKBkDLPf0XxGx1XxVjO8Y8Ci8t87tQXt0rcy19xjWThqntZ2ftZ/DudDQxx9f28Ko6dQGI5jzv899Md3VCVWSOqU/Fz5R3XCdipi7/TYjuvcR5nx1lKKaUw7+aunhL8H6XUGaXUKqVUVzfFOFEplamU2q2UmlXmcY+8l8DtwFda64NXPe6Me1lZ5f1uNlVKNSx+7oDWOvuq5929Dtpgfvp7+apS6qxSakPZrh836F4cxx6l1Atluss89V4mAmeAL6963GX3UynVFPO3VV75gcN/P92RLCq8RlQ5x2YBwcUvyJU5j7PjLGsO5r6+V+axX2DeJbcG1gErlVL1XRzjPCAaaAzMAF5USt1ahfM4O86ybsc09cty1r2srPJ+N8F8T86+n5WmlLoDiAXeKPPwU0BbTBfVu8BipVSkG8L7ErgJaIJ5Eb4VKBn387h7WWw68OFVb7Bcdj+VUr7Af4APtNbp5Rzi8N9PdySLCq8RVc6xdYFLxT+gypzH2XECZqAM8wI3XmudX/K41nqD1jpXa52jtX4V0z84yJUxaq1TtdbHtdZWrfVG4E/AlMqex9lxllBKDQSaAQvKPu7Ee1lZ5f1ugvmenH0/K0UpFQf8DhirtS5dAE9r/Y3WOltrna+1/gDYAIxzdXxa6wNa64Naa5vWehfwG1z3u1lpSqlWwBDgw7KPu+p+KqW8MN24BcD91zjM4b+f7kgWlVkjanfxc+UdtxvoUtzKKNHlGudxdpwope4EngZGaK2P2jm3BpSdYxwe43Vi8Kh7WWw6kKS1vmTn3I66l5VV3u/mKa31ueLn2iqlQq563uXroCmlbgb+CUwsfiG+Hnfdy6td/bvpEfeyjNuBjVrrA3aOc/j9LP4b/RdmA7lErXXhNQ51/O+nKwdjygymfIaZHVMHGMC1Z/DcA6RhmnVhxd/M1bOhHsLM4Lkfx8/gqWicvwBOAtHlPBde/LV+QACmeX0GaOjiGCdjZkUooDdmQHu6p93L4mMDMS2G4a68l8XX8Ck+96uYd28BgE85x91c/DOPKb6vn/Pj2SabMV0+AUA8jp0NVdEYhwPngMHlPFcfM2MmoPh8vwAuAx3dcC/HAk2L/x8FfA/MdsW9rEycZY7PAO509f0svs47xfcj2M5xDv/9dNg3UclvuAGQUnwzjwA/L358EKabqeQ4BbwOZBZ/vM6PZ+x0x0wHzAW2A93dFOdBoBDTvCv5eKf4uU6YweLLxX+4a4FYN8T4afH1LwHpwINXnccj7mXxY7dikpW66nGn3svia8zBvCMs+zEHk6guAeFljn0UMz3xImaMyr/McxGY2TG5mBcXh83gqmiMmDGdoqt+L5cXP9cY+B+m6+EC5sVjlDvuJeZF61Txz/UAphvK1xX3sgo/837FcYZcdQ5X3M/WxbHlXfUz/YUrfj9lbSghhBB2yXIfQggh7JJkIYQQwi5JFkIIIeySZCGEEMIuSRZCCCHskmQhhBDCLkkWQggh7JJkIYQQwi5JFkIIIeySZCGEEMIuSRZCVJBSKrB4L+sjZfczLn7u/4r3kb7FXfEJ4UySLISoIK11LjAbaAXcW/K4UupVzJa/D2itP3NTeEI4lSwkKEQlKKW8MVtQNsHsinYX8BZmSe3fuDM2IZxJkoUQlaSUmgAsxiyRPhx4W2v9oHujEsK5JFkIUQVKqW1AD8ymTj/X8ockajgZsxCikpRS04BuxZ9mS6IQtYG0LISoBKXUaEwX1GLM7ohTgc5a6zS3BiaEk0myEKKClFJ9MOMUWzD7RrfE7BG/TGsd587YhHA26YYSogKUUntHfTEAAACFSURBVNHAUmAPEKe1ztda7wf+BUxWSg1wa4BCOJm0LISwQykVDmwACoD+WutTZZ5rDuwHvtVaS8IQNZYkCyGEEHZJN5QQQgi7JFkIIYSwS5KFEEIIuyRZCCGEsEuShRBCCLskWQghhLBLkoUQQgi7JFkIIYSwS5KFEEIIuyRZCCGEsOv/AeVEfhrtj7W7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#now plot the confidence intervals\n",
    "plt.plot(X, y, \"b.\")\n",
    "plt.plot(X_new, y_predict, \"r-\")\n",
    "plt.plot(X_new, y_new_upper, \"g--\")\n",
    "plt.plot(X_new, y_new_lower, \"g--\")\n",
    "plt.xlabel(\"$x$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([0, 2, 0, 15])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X: [1.] y pred [7.34876878] y down: [4.95849009] y up: [9.73904747]\n"
     ]
    }
   ],
   "source": [
    "print('X:',X_new[100],'y pred',y_predict[100],'y down:',y_new_lower[100],'y up:',y_new_upper[100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coverage of the confidence interval\n",
    "\n",
    "The claim here is that 95% of the time, the predicted interval will cover a new y that gets dropped in the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.93"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_coarse = np.round(X,2)\n",
    "\n",
    "def in_ci(x,y):\n",
    "    bi = int(x*100)\n",
    "    if y >= y_new_lower[bi] and y <= y_new_upper[bi]:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "y_test = np.array(list(map(lambda x,y: in_ci(x,y), X_coarse, y)))\n",
    "coverage = sum(y_test)/len(y_test)\n",
    "round(coverage,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}