notebooks/ols_prediction_error.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Verify the Error Distribution for Ordinary Least Squares\n",
    "*This notebook verifies the equation for ordinary least squares prediction error variance*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.linear_model import LinearRegression\n",
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set Up Simulation Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = 2*np.random.random_sample((25, 3))-1\n",
    "x_prime = 2*np.random.random_sample(3)-1\n",
    "b = np.array([-2, 1, 0.5])\n",
    "sigma = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate a Prediction Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_prediction_error(X, b, x_prime, sigma):\n",
    "    y = np.dot(X, b) + np.random.normal(scale=sigma, size=X.shape[0])\n",
    "    model = LinearRegression(fit_intercept=False)\n",
    "    model.fit(X, y)\n",
    "    y_pred = model.predict([x_prime])[0]\n",
    "    y_prime = np.dot(x_prime, b) + np.random.normal(scale=sigma)\n",
    "    return y_prime - y_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Simulation Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100000\n",
    "errors = np.array([generate_prediction_error(X, b, x_prime, sigma) for _ in range(N)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Running Mean and Variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "means = []\n",
    "variances = []\n",
    "total = 0\n",
    "total2 = 0\n",
    "for i, e in enumerate(errors):\n",
    "    total += e\n",
    "    total2 += e*e\n",
    "    means.append(total / (i + 1))\n",
    "    variances.append(total2 / (i + 1))\n",
    "means = np.array(means)\n",
    "variances = np.array(variances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Expected Variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.linalg.inv(np.dot(X.T, X))\n",
    "expected_error_variance = sigma**2*(1 + np.dot(x_prime, np.dot(A, x_prime)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the Running Stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD4CAYAAAD//dEpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAf50lEQVR4nO3de5xV5X3v8c+XGRjkIneQgAoq0ZA0NTpR09ysCF6SqD01Dakviz3wsrZNbWPbU6iJnhrTgyc9uZ22SYjRkh6TaGxaabRSwJCYNl4GRYtBBfE2iDBcRZHLML/zx34G1gxrz6xh72GG4ft+vfZrr/WsZ631LNYw33nWVRGBmZlZZ/r1dAPMzOzo4MAwM7NCHBhmZlaIA8PMzApxYJiZWSG1Pd2Aaho9enRMmjSpp5thZnZUWbFixeaIGNNZvT4VGJMmTaKhoaGnm2FmdlSR9HKRej4kZWZmhTgwzMysEAeGmZkV4sAwM7NCHBhmZlaIA8PMzApxYJiZWSEOjOSRdVtYu2lnTzfDzKzX6lM37lVi5oJHAHhp/sd6uCVmZr2TexhmZlaIA8PMzApxYJiZWSEODDMzK8SBAURETzfBzKzXc2C04/AwM8vnwACyGfHztZt7riFmZr2YA6Odvc0tPd0EM7NeyYFhZmaFODAAn7UwM+ucA8PMzApxYND2yqiXt+zqwZaYmfVeVQkMSRdLek7SWklzc6bXSbo7TX9U0qRUPl3SCkn/lb4vyMyzPC1zZfqMrUZbO/Pa9rePxGrMzI46FT+tVlIN8HfAdKAReFzSooj4ZababGBbRJwmaSZwG/ApYDPwiYh4TdJ7gMXAhMx8V0VEQ6Vt7Ez2HMb7J4/s7tWZmR2VqtHDOAdYGxHrImIv8APg8nZ1LgcWpuF7gWmSFBFPRsRrqfwZYKCkuiq06bD5vj0zs3zVCIwJwKuZ8Uba9hLa1ImIZmAHMKpdnd8EnoyIPZmyO9PhqM9LUt7KJV0rqUFSQ1NTUyXbkTgxzMzyVCMw8n6Rt/+t22EdSe+mdJjq9zLTr4qIXwE+nD5X5608IhZERH1E1I8ZM6ZLDT+4jPxhMzM7qBqB0QicmBmfCLxWro6kWmAYsDWNTwT+GfidiHihdYaIWJ++dwLfo3Toq9s5L8zM8lUjMB4HpkiaLGkAMBNY1K7OImBWGr4SeCgiQtJw4H5gXkT8R2tlSbWSRqfh/sDHgVVVaGuuyMSEexhmZvkqDox0TuIzlK5wWg3cExHPSLpF0mWp2neAUZLWAjcArZfefgY4Dfh8u8tn64DFkp4GVgLrgW9X2tZC2+M+hplZroovqwWIiAeAB9qV3ZQZ3g18Mme+W4Fbyyz27Gq0rQifwzAz65zv9DYzs0IcGMC6prcODLuDYWaWz4EB7Hh734Fhv3HPzCyfAwPIvyXQzMyyHBjtuINhZpbPgUHb29B9Wa2ZWT4HBpB9TJV7GGZm+RwY7TgwzMzyOTDa6Wpe7N63nxUvb+2WtpiZ9SZVudO7L9n85p7OK2Wc8fkHAfjQaaP5f3PO7Y4mmZn1Cg4M2l5W+8TL2wrNM/sfHmfZs5sOjP987WZe3vIWJ48aXO3mmZn1Cj4k1c6FU8cVqpcNi1Yf/dJynt+4s9pNMjPrFRwY5L/d6XDN+MrPqrg0M7Pew4HRDf5z7WZe3brrwPibe5p7sDVmZtXhcxjtFOlt7G1uaTO+7E8/yrT/89MD4799+6MArLxpOmfesuRA+Qt/fSk1/fwcEjM7OrmHAZwwbGCX6v/zk40Hhp+79WJOHTOEl+Z/7JB62bAAOPUvH6Bp5x72NO8vtJ6b71vFpLn3c9ejL/PdX7zErr3N/PrfLG8TWJd+7WEmzb2fSXPv54a7V/LlJc/T0hJEBC0tpYuEn3hl24E6k+be36VtNTNrpb70dNb6+vpoaGjo8nzrt7/NB+c/BMCHp4zmH2d3fHls9pduNii68su44XMXUlfbj6ED+x8yLSKYPO+BnLmq74bp7+T6aVMOjK9av4NH1m3h104dzaVffxiApTd8lHcMH8igAe6QmvVFklZERH1n9fwboJ2H12wuXPeG6e9sM/7H06bwtWVrDqn30vyPHRIm9bcuBWDJZz/Cn9/7NP/0+79GTT+xe9/+A/d2HAlfXvI8X17yfId1LvzyT9uMN3zuQkYPqevOZplV1Zt7mvnY1x/m5S27uPq8k3no2U38t7Mm8NiLWzn3lFGcOmYwe/a18MbufYwaMoD3TxrJxBGDiIg2jw46UvY2t7Dj7X3sad7P6zt2M3zQACaOOI6B/WuOeFuyqtbDkHQx8DWgBrg9Iua3m14HfJfSq1e3AJ+KiJfStHnAbGA/cH1ELC6yzPaq0cMAWP5n57N09UaaW4LZH5rMb3/7EeZecgZnnzySHW/v41f/6t8BWH3LxRw3oO0OnDT3fqaOP57BdTU8/tI2nrp5BsOO639gWkceu3Ea53xx2YHx6VPH8fdXnUWNxFt7m3l7334+98+r+Pdfbmwz3+pbLmZg/37sSYeq8gLnqZtmMGxQqR1v7mnmPTcvLvrPk+vZL1zc4z+8Vl5E8Pa+/W16hZvf3MOy1Rv5zxe28LsfnMzjL27liw+sBuDkUYO49Yr38MFTR9Ovn3hrTzN1tf2orelHRLBp5x5GDBrAgNr8o9jd9Yv1qVe3s/y5Js6ZPJKTRw1iYP8atr61h1Xr3+DlLbv41s9e4L0Th3H6uKHs2rufH65o5MZL34UEt96/uurtafXudxzP6eOG8qMn1+dOnzRqEOOOH8ijLx58CsRF7x7Hz9ds5q29+5kw/DjWb3/7kPmmjB3Cmk1vdrk9p44ZzOc+PpVfP31sl+eF4j2MqgSGpBrgeWA60Ag8Dnw6In6ZqfMHwHsj4jpJM4HfiIhPSZoKfB84B3gHsBRo/dO9w2W2d7iB0bhtFx+67Se5074280z++AcrAbh+2hQu+9V3HPiLO++8RRHfWP4Ctz34bKf1Dnf5XdG8v4V+Ej99vonnN+7k2o+ccsh//FXrd7B91z5OGjmIj3wp/9+p1dxLzuC6j56aO+31Hbv5wv2l3Xf/0xsKzVNNe5tb2LRzN89v3MnufS38wV1P8LM//3W27drLNXc+xrZdpRdpDR/Un3FDB3LPdR9g7aadfOpbj9Cczgf90QWn8X8fWnvIsme+/0Ted9Jwfqv+xG77i/Ttvft5101HrvdZibNOGs4Tr2zPnfaXl57BYy9uZenqQ+9l6k5D62q57vxT+Y33TeDO/3iR+1a+xm+8bwLDBw3g9BOG8PqOPdz+83Vt3sCZ9bH3jm/zc9uZM04YyrOv59+XdeaJw1n56qH/PuedMpJH1h0MmXeOG8Ke5hZe3lK66nLc8XUMHlDLus2HtvGPLjiNP51xeuH2ZR3pwPgA8D8j4qI0Pg8gIv5Xps7iVOcXkmqB14ExwNxs3dZ6abYOl9le/dCh0XD22V1u/57mFp58pdgd3q1OGDaQSRXe1f3Iui0AHDeghrf3tj0Rfu4po6p6f0i1BPBoand3OWfyKJpbWti3v4XBOedN9u5v4bXtuzl51KAD/0YtETz2Yt94ptepY4Yc6Lk2bnub7bv2Vm3Zktq8VXL4oAG8vXd/hxdiDOxfQ/+afuzcva9snSPp+OP6s29/C6OH1PHqtrcPPDF0zNA6QIwc3J+62hoG1PajtpuvSty1dz81/UT/GtFPoiX1tFrXGpSuvGz97i6VLl8//ekRPYcxAXg1M94ItD9zfKBORDRL2gGMSuWPtJt3QhrubJlIuha4FuC9dUfuuPq2XfuYNKqyZZx3ysEFNLcEDS9tPaS8txGl9u3c08wz63d0yzoee7FYIL2+49AufSXGDB1I/xrxWs6hgtFD6tjXEuzYtZfamn68c9xQWiI4fmB/nnxlG/v2t+Qs8fC80NTxIYm6/jUMHVjLuOMHMnhALf108KGZveWPjKAU4jXZVwdQOixW00/U9hOD62rbTN+Xertv7N7H8EEDCm3LhOHHVb3tXTGo3SHpfu16l2r33V2O1H6vVmDktbd916VcnXLleQdLD+kORcQCYAGUDkmxfHmHDc3TtG0XM8sckirnges/DO84vsvrKqcWOK9qS+t+Qynf3q8seZ6vLVvDVeeexLOv72TB1Wez4Gfr+NMZp+ceA9/61l7O+sISpk8dx5J252e66pTRg1m3+S0G1PRj7/4W3jX+eFZveIM1X7yE/jXFriI/qYvrLNenbf1Lfk9zC2s3vcm7xh/f5j6c/S3Bk69sY+nqTdzx8xf52swz+f27nii1YeQgXkk3fy757Ec4ZcyQDu/h6S1B0UqUTjy2LxvTwTyt1wuO6JYWWYcKHkatVmA0AidmxicCr5Wp05gOSQ0DtnYyb2fL7DFTqxgWfc1np7+Tz7a7gmzepe8qW3/k4AEdnq9paQmWP7+Jpp17GDW47pDnfW17ay8jBg+orNHdoPVcxsD+NbxnwrBDptf0E/WTRlI/aSRzLzkDODLnrcwOV7UC43FgiqTJwHpgJvDb7eosAmYBvwCuBB6KiJC0CPiepC9TOuk9BXiM0h8knS3TjgH9+okLzij/UMjeGBZmfVFV7vSOiGbgM8BiYDVwT0Q8I+kWSZelat8BRklaC9zAwZPdzwD3AL8EHgT+MCL2l1tmNdp7aPs7r9PwuQu7Y9VmZkeNqt24FxEPAA+0K7spM7wb+GSZeb8IfLHIMntK/34Hs/XOa97fgy0xM+sZfpYUsPGN3R1OnzF13IGb3gDeP3lkdzfJzKzXcWAAd/zHix1Obz0h2WpInZ+oYmbHHv/m68S6v76UfulyxqdumsFbe/1uCzM7NjkwOnDWScMPhAXAsEH92xyaMjM7lviQVAd64imVZma9lQOD8pfVThk75Mg2xMysF3Ng5DguPbb7I+/s6EEGZmbHFp/D4NDHqJx/+hj+4uIzmDS6sqfRmpn1Je5h5Ljq3JMdFmZm7TgwcnxoyuieboKZWa/jwKDtSe/+Nb4yyswsjwOjHfW6NwuYmfUODox2vnn1WT3dBDOzXsmB0U5H710wMzuWOTDMzKwQB0bGcD8nysysLAdGxtihdT3dBDOzXsuBYWZmhVQUGJJGSloiaU36HlGm3qxUZ42kWalskKT7JT0r6RlJ8zP1r5HUJGll+syppJ1mZla5SnsYc4FlETEFWJbG25A0ErgZOBc4B7g5Eyx/ExFnAO8DPijpksysd0fEmelze4XtNDOzClUaGJcDC9PwQuCKnDoXAUsiYmtEbAOWABdHxK6I+AlAROwFngAmVtiew1Lu8eZmZnZQpYExLiI2AKTvsTl1JgCvZsYbU9kBkoYDn6DUS2n1m5KelnSvpBMrbGchvsvbzKy8Th9vLmkpcELOpBsLriPvt/CBv+kl1QLfB74eEetS8b8C34+IPZKuo9R7uaBM+64FrgU46aSTCjYpX+CuhplZOZ0GRkRcWG6apI2SxkfEBknjgU051RqB8zPjE4HlmfEFwJqI+GpmnVsy078N3NZB+xakZVBfX+/f+GZm3aTSQ1KLgFlpeBZwX06dxcAMSSPSye4ZqQxJtwLDgD/JzpDCp9VlwOoK21mID0mZmZVXaWDMB6ZLWgNMT+NIqpd0O0BEbAW+ADyePrdExFZJEykd1poKPNHu8tnr06W2TwHXA9dU2E4zM6tQRa9oTYeOpuWUNwBzMuN3AHe0q9NI/vkNImIeMK+StpmZWXX5Tm8zMyvEgZHhq6TMzMpzYJiZWSEODA72LHyVlJlZeQ6MDDkvzMzKcmCYmVkhDowMP4TQzKw8BwY+d2FmVoQDA19Oa2ZWhAMjwye9zczKc2CYmVkhDgwzMyvEgZHhq6TMzMpzYOCrpMzMinBg4KukzMyKcGCYmVkhDgx8SMrMrAgHRoYPTZmZlefAwDfsmZkVUXFgSBopaYmkNel7RJl6s1KdNZJmZcqXS3pO0sr0GZvK6yTdLWmtpEclTaq0rZ3xZbVmZuVVo4cxF1gWEVOAZWm8DUkjgZuBc4FzgJvbBctVEXFm+mxKZbOBbRFxGvAV4LYqtDVXaw/DeWFmVl41AuNyYGEaXghckVPnImBJRGyNiG3AEuDiLiz3XmCa5INHZmY9pRqBMS4iNgCk77E5dSYAr2bGG1NZqzvT4ajPZ0LhwDwR0QzsAEa1X7CkayU1SGpoamqqfGvMzCxXbZFKkpYCJ+RMurHgevJ6Bq1HgK6KiPWShgL/BFwNfLeTeQ4WRCwAFgDU19f7qJKZWTcpFBgRcWG5aZI2ShofERskjQc25VRrBM7PjE8Elqdlr0/fOyV9j9I5ju+meU4EGiXVAsOArUXaa2Zm1VeNQ1KLgNarnmYB9+XUWQzMkDQineyeASyWVCtpNICk/sDHgVU5y70SeCiie65j8o17ZmadK9TD6MR84B5Js4FXgE8CSKoHrouIORGxVdIXgMfTPLekssGUgqM/UAMsBb6d6nwH+EdJayn1LGZWoa0d6qY8MjPrEyoOjIjYAkzLKW8A5mTG7wDuaFfnLeDsMsvdTQofMzPreb7T28zMCnFgwIHrsXxAysysPAcG+dfvmplZWw6MLHcxzMzKcmAAfuKImVnnHBj4clozsyIcGFnuaJiZleXAyHJHw8ysLAcGB89hOC/MzMpzYOAjUWZmRTgwMnzy28ysPAeGmZkV4sAwM7NCHBiA79szM+ucA8PMzApxYGT4lLeZWXkODHxZrZlZEQ4MYMTgAQBMHHFcD7fEzKz3qigwJI2UtETSmvQ9oky9WanOGkmzUtlQSSszn82SvpqmXSOpKTNtTt5yq+WDp44G4Ibpp3fnaszMjmqV9jDmAssiYgqwLI23IWkkcDNwLnAOcLOkERGxMyLObP0ALwM/ysx6d2b67RW2s5ABNe5wmZmVU+lvyMuBhWl4IXBFTp2LgCURsTUitgFLgIuzFSRNAcYCD1fYHjMz6yaVBsa4iNgAkL7H5tSZALyaGW9MZVmfptSjyF6o9JuSnpZ0r6QTyzVA0rWSGiQ1NDU1Hd5WmJlZpzoNDElLJa3K+VxecB15FyG1v4J1JvD9zPi/ApMi4r3AUg72Yg5dUMSCiKiPiPoxY8YUbJKZmXVVbWcVIuLCctMkbZQ0PiI2SBoPbMqp1gicnxmfCCzPLONXgdqIWJFZ55ZM/W8Dt3XWTjMz616VHpJaBMxKw7OA+3LqLAZmSBqRrqKakcpafZq2vQtS+LS6DFhdYTs75Bv2zMw612kPoxPzgXskzQZeAT4JIKkeuC4i5kTEVklfAB5P89wSEVszy/gt4NJ2y71e0mVAM7AVuKbCdhbiZ0qZmZVXUWCkQ0fTcsobgDmZ8TuAO8os45ScsnnAvEraZmZm1eUbD8zMrBAHhpmZFeLAMDOzQhwYZmZWiAPDzMwKcWCYmVkhDgwzMyvEgQG0feahmZnlcWCYmVkhDgwzMyvEgWFmZoU4MMzMrBAHhpmZFeLAMDOzQhwYZmZWiAPDzMwKcWCYmVkhDgz8Tm8zsyIqDgxJIyUtkbQmfY8oU+9BSdsl/bhd+WRJj6b575Y0IJXXpfG1afqkStva+bZ09xrMzI5e1ehhzAWWRcQUYFkaz/Ml4Oqc8tuAr6T5twGzU/lsYFtEnAZ8JdUzM7MeUo3AuBxYmIYXAlfkVYqIZcDObJkkARcA9+bMn13uvcC0VN/MzHpANQJjXERsAEjfY7sw7yhge0Q0p/FGYEIangC8mpbbDOxI9duQdK2kBkkNTU1Nh7kJZmbWmdoilSQtBU7ImXRjhevP6zFEgWkHCyIWAAsA6uvrff7azKybFAqMiLiw3DRJGyWNj4gNksYDm7qw/s3AcEm1qRcxEXgtTWsETgQaJdUCw4CtXVi2mZlVUTUOSS0CZqXhWcB9RWeM0puLfgJcmTN/drlXAg+F33RkZtZjqhEY84HpktYA09M4kuol3d5aSdLDwA8pnbxulHRRmvQXwA2S1lI6R/GdVP4dYFQqv4HyV1+ZmdkRUOiQVEciYgswLae8AZiTGf9wmfnXAefklO8GPllp+8zMrDp8pzfgA11mZp1zYGQo98IsMzMDB4aZmRXkwDAzs0IcGGZmVogDw8zMCnFgmJlZIQ4MMzMrxIFhZmaFODDMzKwQBwbgt3qbmXXOgZHh9/mZmZXnwDAzs0IcGGZmVogDw8zMCnFgmJlZIQ4MMzMrxIFhZmaFVBQYkkZKWiJpTfoeUabeg5K2S/pxu/K7JD0naZWkOyT1T+XnS9ohaWX63FRJO83MrHKV9jDmAssiYgqwLI3n+RJwdU75XcAZwK8Ax5F5BzjwcEScmT63VNhOMzOrUKWBcTmwMA0vBK7IqxQRy4CdOeUPRAI8BkyssD2Hxe/0NjPrXKWBMS4iNgCk77GHs5B0KOpq4MFM8QckPSXp3yS9u4N5r5XUIKmhqanpcFafWVZFs5uZ9Wm1nVWQtBQ4IWfSjVVsx98DP4uIh9P4E8DJEfGmpEuBfwGm5M0YEQuABQD19fXuK5iZdZNOAyMiLiw3TdJGSeMjYoOk8cCmrjZA0s3AGOD3Mut8IzP8gKS/lzQ6IjZ3dflmZlYdlR6SWgTMSsOzgPu6MrOkOcBFwKcjoiVTfoJUOkAk6ZzUzi0VtrWs/ekkRktLJxXNzI5hlQbGfGC6pDXA9DSOpHpJt7dWkvQw8ENgmqRGSRelSd8ExgG/aHf57JXAKklPAV8HZqYT493iwVWvA7C26c3uWoWZ2VGv00NSHYmILcC0nPIGMpfIRsSHy8yfu/6I+Fvgbytpm5mZVZfv9DYzs0IcGGZmVogDw8zMCnFg4Dd6m5kV4cDI8I3eZmblOTDMzKwQB4aZmRXiwDAzs0IcGGZmVogDw8zMCnFgmJlZIQ4MMzMrxIEB1NWU/hlq+vlODDOzcip6Wm1fcdMnpjL2+IHMmDqup5tiZtZrOTCA4YMGMPeSM3q6GWZmvZoPSZmZWSEODDMzK8SBYWZmhVQUGJJGSloiaU36HlGm3oOStkv6cbvyf5D0Ynqf90pJZ6ZySfq6pLWSnpZ0ViXtNDOzylXaw5gLLIuIKcCyNJ7nS8DVZab9eUScmT4rU9klwJT0uRb4RoXtNDOzClUaGJcDC9PwQuCKvEoRsQzY2cXlfjdKHgGGSxpfUUvNzKwilQbGuIjYAJC+xx7GMr6YDjt9RVJdKpsAvJqp05jKDiHpWkkNkhqampoOY/VmZlZEp4EhaamkVTmfy6uw/nnAGcD7gZHAX7SuNqdu7ptUI2JBRNRHRP2YMWOq0CQzM8vT6Y17EXFhuWmSNkoaHxEb0iGjTV1ZeWvvBNgj6U7gz9J4I3BipupE4LXOlrdixYrNkl7uShsyRgObD3Peo5W3+djgbT42VLLNJxepVOmd3ouAWcD89H1fV2bOhI0onf9YlVnuZyT9ADgX2JEJl7Ii4rC7GJIaIqL+cOc/Gnmbjw3e5mPDkdjmSgNjPnCPpNnAK8AnASTVA9dFxJw0/jClQ09DJDUCsyNiMXCXpDGUDkGtBK5Ly30AuBRYC+wCfrfCdpqZWYUqCoyI2AJMyylvAOZkxj9cZv4LypQH8IeVtM3MzKrLd3oftKCnG9ADvM3HBm/zsaHbt1mlP+bNzMw65h6GmZkV4sAwM7NCHBiApIslPZcedljueVi9kqQTJf1E0mpJz0j641Se+2DIjh7sKGlWqr9G0qxM+dmS/ivN8/V0GXSPk1Qj6cnWh1pKmizp0dT+uyUNSOV1aXxtmj4ps4x5qfw5SRdlynvdz4Sk4ZLulfRs2t8f6Ov7WdJn08/1KknflzSwr+1nSXdI2iRpVaas2/druXV0KCKO6Q9QA7wAnAIMAJ4CpvZ0u7rQ/vHAWWl4KPA8MBX438DcVD4XuC0NXwr8G6VLmc8DHk3lI4F16XtEGh6Rpj0GfCDN82/AJT293aldNwDfA36cxu8BZqbhbwK/n4b/APhmGp4J3J2Gp6b9XQdMTj8HNb31Z4LS89rmpOEBwPC+vJ8pPQ7oReC4zP69pq/tZ+AjwFnAqkxZt+/XcuvosK09/Z+gpz/pH3JxZnweMK+n21XB9twHTAeeA8ansvHAc2n4W8CnM/WfS9M/DXwrU/6tVDYeeDZT3qZeD27nREpPSL4A+HH6z7AZqG2/X4HFwAfScG2qp/b7urVeb/yZAI5PvzzVrrzP7mcOPlNuZNpvPwYu6ov7GZhE28Do9v1abh0dfXxIqgsPOuztUhf8fcCjlH8wZLnt7ai8Mae8p30V+B9ASxofBWyPiOY0nm3ngW1L03ek+l39t+hJpwBNwJ3pMNztkgbTh/dzRKwH/obSTcEbKO23FfTt/dzqSOzXLj881oHRhQcd9maShgD/BPxJRLzRUdWcsjiM8h4j6ePApohYkS3OqRqdTDtqtpnSX8xnAd+IiPcBb1H+/TPQB7Y5HVO/nNJhpHcAgym9K6e9vrSfO9Oj2+jAOMwHHfYmkvpTCou7IuJHqXij0jtE1PbBkOW2t6PyiTnlPemDwGWSXgJ+QOmw1FcpvTel9ekF2XYe2LY0fRiwla7/W/SkRqAxIh5N4/dSCpC+vJ8vBF6MiKaI2Af8CPg1+vZ+bnUk9mu5dZTlwIDHgSnpyosBlE6WLerhNhWWrnj4DrA6Ir6cmdT6YEho+2DIRcDvpKstzuPggx0XAzMkjUh/2c2gdHx3A7BT0nlpXb9DFx8yWW0RMS8iJkbEJEr766GIuAr4CXBlqtZ+m1v/La5M9SOVz0xX10ym9IbHx+iFPxMR8TrwqqTTU9E04Jf04f1M6VDUeZIGpTa1bnOf3c8ZR2K/lltHeT15Uqu3fChdefA8pSsmbuzp9nSx7R+i1MV8mtIDHFem7RlF6aTwmvQ9MtUX8HdpW/8LqM8s679TeuDjWuB3M+X1lJ4k/ALwt7Q78drD238+B6+SOoXSL4K1wA+BulQ+MI2vTdNPycx/Y9qu58hcFdQbfyaAM4GGtK//hdLVMH16PwN/BTyb2vWPlK506lP7Gfg+pXM0+yj1CGYfif1abh0dffxoEDMzK8SHpMzMrBAHhpmZFeLAMDOzQhwYZmZWiAPDzMwKcWCYmVkhDgwzMyvk/wMb7FvqYtOuLAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x1a15352110>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAdgUlEQVR4nO3de3wV9Z3/8dcnCUkkcieWu1xMtVQQMSJeqoI3lLbYrr/fom7rum2pXdi61dbLz8vPaq0tbq29YFfbanfbIlrbbalFqVq00lolCgpBUgIoIAgBuYqQ22f/OJMwSU5yTi5wDjPv5+PBg5nvfOeczzDhnTnfmTNj7o6IiMRDTqYLEBGRw0ehLyISIwp9EZEYUeiLiMSIQl9EJEbyMl1Ac/379/fhw4dnugwRkSPKq6++us3di1P1y7rQHz58OGVlZZkuQ0TkiGJmb6fTT8M7IiIxotAXEYkRhb6ISIykFfpmNsXMKsys0sxuSrL8GjNbbmbLzGyxmY0O2ruZ2X8Fy940s5u7egNERCR9KUPfzHKBOcDFwGjg8oZQD5nr7mPcfRwwG7gvaP8/QIG7jwFOAb5oZsO7qHYREWmndI70JwCV7r7W3auBecC0cAd33x2aLQIa7uLmQJGZ5QFHAdVAuK+IiBxG6VyyORjYEJrfCJzWvJOZzQSuA/KByUHzEyR+QWwGugNfcff3OlOwiIh0XDpH+pakrcX9mN19jruPAm4Ebg2aJwB1wCBgBHC9mY1s8QZmM8yszMzKqqqq0i6+uWdXbmHL7v0dXl9EJOrSCf2NwNDQ/BBgUxv95wGXBtNXAE+7e427bwX+ApQ2X8HdH3L3UncvLS5O+YWyVn3+v8v49AN/7fD6IiJRl07oLwFKzGyEmeUD04H54Q5mVhKanQqsDqbXA5MtoQiYCKzqfNmte2fnB4fy5UVEjmgpx/TdvdbMZgELgVzgYXcvN7M7gTJ3nw/MMrPzgRpgB3BVsPoc4BFgBYlhokfc/Y1DsB0iIpKGtO694+4LgAXN2m4PTV/bynp7SVy2KSIiWUDfyBURiRGFvohIjCj0RURiRKEvIhIjCn0RkRhR6IuIxIhCX0QkRhT6IiIxotAXEYkRhb6ISIwo9EVEYkShLyISIwp9EZEYUeiLiMSIQl9EJEYU+iIiMaLQFxGJEYW+iEiMKPRFRGJEoS8iEiMKfRGRGFHoi4jEiEJfRCRGFPoiIjGi0BcRiZHIhH5tXX2mSxARyXqRCf3yTbszXYKISNaLTOiLiEhqCn0RkRhR6IuIxIhCX0QkRhT6IiIxklbom9kUM6sws0ozuynJ8mvMbLmZLTOzxWY2OrRsrJm9ZGblQZ/CrtwAERFJX8rQN7NcYA5wMTAauDwc6oG57j7G3ccBs4H7gnXzgF8A17j7R4FzgZquKz9c56F4VRGRaEnnSH8CUOnua929GpgHTAt3cPfwRfJFgAfTFwJvuPvrQb/t7l7X+bJFRKQj0gn9wcCG0PzGoK0JM5tpZmtIHOl/OWj+MOBmttDMXjOzG5K9gZnNMLMyMyurqqpq3xaIiEja0gn9ZAMn3qLBfY67jwJuBG4NmvOAs4Arg78/ZWbnJVn3IXcvdffS4uLitIsXEZH2SSf0NwJDQ/NDgE1t9J8HXBpa9wV33+bu+4AFwPiOFCoiIp2XTugvAUrMbISZ5QPTgfnhDmZWEpqdCqwOphcCY82se3BS9xxgZefLFhGRjshL1cHda81sFokAzwUedvdyM7sTKHP3+cAsMzufxJU5O4CrgnV3mNl9JH5xOLDA3f9wKDbEko5CiYhIWMrQB3D3BSSGZsJtt4emr21j3V+QuGzzkPKWpxlERKQZfSNXRCRGFPoiIjGi0BcRiRGFvohIjCj0RURiJDKhr0s2RURSi0zoi4hIagp9EZEYUeiLiMSIQl9EJEYU+iIiMRKZ0NfjEkVEUotM6IuISGoKfRGRGFHoi4jEiEJfRCRGFPoiIjGi0BcRiRGFvohIjCj0RURiRKEvIhIjCn0RkRhR6IuIxIhCX0QkRhT6IiIxotAXEYkRhb6ISIwo9EVEYkShLyISIwp9EZEYUeiLiMRIWqFvZlPMrMLMKs3spiTLrzGz5Wa2zMwWm9noZsuHmdleM/tqVxXesoZD9coiItGRMvTNLBeYA1wMjAYubx7qwFx3H+Pu44DZwH3Nln8XeKoL6hURkU5I50h/AlDp7mvdvRqYB0wLd3D33aHZIsAbZszsUmAtUN75ckVEpDPSCf3BwIbQ/MagrQkzm2lma0gc6X85aCsCbgS+3tYbmNkMMyszs7Kqqqp0axcRkXZKJ/STjZZ7iwb3Oe4+ikTI3xo0fx34rrvvbesN3P0hdy9199Li4uI0ShIRkY7IS6PPRmBoaH4IsKmN/vOAHwXTpwGXmdlsoDdQb2b73f2HHSlWREQ6J53QXwKUmNkI4B1gOnBFuIOZlbj76mB2KrAawN0/FupzB7D3UAW+Jf1AIiIiYSlD391rzWwWsBDIBR5293IzuxMoc/f5wCwzOx+oAXYAVx3KokVEpGPSOdLH3RcAC5q13R6avjaN17ijvcWJiEjX0jdyRURiRKEvIhIjCn0RkRhR6IuIxIhCX0QkRhT6IiIxotAXEYkRhb6ISIwo9EVEYkShLyISI5EJfT0uUUQktciEvoiIpKbQFxGJEYW+iEiMKPRFRGJEoS8iEiMKfRGRGIlM6OuSTRGR1CIT+iIikppCX0QkRhT6IiIxotAXEYkRhb6ISIxEJvQNXb4jIpJKZEJfRERSi0zoO57pEkREsl5kQl9ERFKLZOivenc32/ceyHQZIiJZJ5KhP+X+F7no/hczXYaISNaJTOg3v3pnm470RURaiEzoi4hIammFvplNMbMKM6s0s5uSLL/GzJab2TIzW2xmo4P2C8zs1WDZq2Y2uas3QERE0pcy9M0sF5gDXAyMBi5vCPWQue4+xt3HAbOB+4L2bcAn3H0McBXw8y6rXERE2i2dI/0JQKW7r3X3amAeMC3cwd13h2aLIHHRvLsvdfdNQXs5UGhmBZ0vW0REOiIvjT6DgQ2h+Y3Aac07mdlM4DogH0g2jPMPwFJ3b3GG1cxmADMAhg0blkZJIiLSEekc6Se7qU2Lr7+6+xx3HwXcCNza5AXMPgp8G/hisjdw94fcvdTdS4uLi9MoSUREOiKd0N8IDA3NDwE2tdIXEsM/lzbMmNkQ4H+Az7r7mo4UmQ49LlFEJLV0Qn8JUGJmI8wsH5gOzA93MLOS0OxUYHXQ3hv4A3Czu/+la0oWEZGOShn67l4LzAIWAm8Cj7t7uZndaWafDLrNMrNyM1tGYlz/qoZ24DjgtuByzmVmdkzXb4aIiKQjnRO5uPsCYEGztttD09e2st43gG90pkAREek6+kauiEiMRCb0XbfTFxFJKTKhv2d/TaZLEBHJepEJ/WSqa+szXYKISFaJTOgnu07/z3+vOvyFiIhksciEfrIvDr+2fkcG6hARyV6RCf1kR/pjh/Q6/IWIiGSxyIR+QV7LTdm444MMVCIikr0iE/rNH5cI8MfyLRmoREQke0Un9JMM79TW6+odEZGwyIR+si9nTTlxwOEvREQki0Um9JP55oJVmS5BRCSrRDr0RUSkKYW+iEiMRDr0Z04alekSRESySmRCP9nVO93z03pcgIhIbEQm9JPRl7NERJqKdOg/+sr6TJcgIpJVIhP6eoiKiEhqkQl9ERFJTaEvIhIjCn0RkRiJTOgnu2RzVHHR4S9ERCSLRSb0k1lT9X6mSxARySqR/fZSbo7xkYE9Ml2GiEhWieyR/uiBPVnxzm4+8YPFmS5FRCRrRCb0a+uaXqifm5MY5F/+zq5MlCMikpUiE/qbdjW95UKyE7siInEXmdBvnvFL1+/MSB0iItksMqEvIiKpKfRFRGIkMqGv+62JiKSWVuib2RQzqzCzSjO7Kcnya8xsuZktM7PFZjY6tOzmYL0KM7uoK4tP19Y9+zPxtiIiWSdl6JtZLjAHuBgYDVweDvXAXHcf4+7jgNnAfcG6o4HpwEeBKcADwet1ubYu1nnohbWH4i1FRI446RzpTwAq3X2tu1cD84Bp4Q7uvjs0W8TB0ZZpwDx3P+Du64DK4PUOq9p6Df6IiEB6t2EYDGwIzW8ETmveycxmAtcB+cDk0Lp/a7bu4CTrzgBmAAwbNiydultoK9br9YQVEREgvSP9ZCMnLVLU3ee4+yjgRuDWdq77kLuXuntpcXFxGiW1T01dfZe/pojIkSid0N8IDA3NDwE2tdF/HnBpB9ftMuvuuaRx+tFXNrTRU0QkPtIJ/SVAiZmNMLN8Eidm54c7mFlJaHYqsDqYng9MN7MCMxsBlACvdL7slpp/pDDdh0FEpIWUY/ruXmtms4CFQC7wsLuXm9mdQJm7zwdmmdn5QA2wA7gqWLfczB4HVgK1wEx3rzsUG6JRexGR1NK6n767LwAWNGu7PTR9bRvr3g3c3dECO6NfUT7b36/OxFuLiGSlyHwjNzyYc0yPAgC+csGHG9tcV/CIiEQn9I/KP/idr+svTIT95BOOaWwr37S7xToiInETmdDvWditcbrhJO7AXoWNbQvL3z3sNYmIZJvse0ZuRQWce267Vxt1oJZ5wVOyRj5zNPQowIB5a7cnOswFRvbrsjJFRI5EkTnSDzucV2tW1+qLXyJy5Mi+I/3jj4fnn2/3ams27GT6nL8kpr95CQTPyJ1+0x8a+7z1raldUiLArn01nHTnHxvnV901hcJu7b+X3La9B+h/dEGX1SUimVe5dS+DehfSPb/tiN2+9wB5OTn88pW32f1BLTPOHknfovyOvWmaR7vZF/qdNHZIr8aHogOMGdyrxcPR397+Ph/qWdihkD5QW8eXfvEaf1q1tUn7Cbc9TeXdF5OXm/6Hp+GhX0iXjBnAA1ee0u56RKR19fXO+9W1HF2Q1+QLmy+t2c7lP07cFqzXUd34+NiB/PLl9Qzv151JJxzDkrfe4+ySYh54fk3jOoN6FXJsvyJe37iTfdV1XHteCd97bnXj8o8O6sm+6jrWbXs/ZV2XjBnAguUtzzP+be12fjvzzM5sckqRC/3mLj15cJPQX1u1l8nfeYHTR/bj0RkT2/Va+2vqOOG2p1tdftwtT7X4NPG7Ze8wflgfhvbtztY9+/mHH/2VF2+Y3CTwARYsf5f6eicnR98kFmmv+a9v4quPv051B+6zteuDGn758noA3tq+j0f+8hYAK95pesXfpl372bTr4LM5woEP7btCsCHwR/Yvolf3bixdv5OJI/sya1JJijU7L/KhX9jt4JF3TV09k7/zAgAvrd3O/yzdyNQxg8jPa/vo/BtPruQni9clXTb7srHc8MQbjfPujplxoLaOhxe/xbefXgUkhpYm3P0cQIvAbzDy/x38/tvyOy5kzB1/5Oefm8DHSop5cXUV185bxmu3XZDGVotkt5q6ehYs38zmXfv51lOrGtsL8nI4UFvPPZ8ewydOGsTS9TvYsvsAS9a9x2NlG5jy0QE8Xf4uw/p2Z/17+zpVw0cG9uSK04ZxbN/uDOpdyF/XbOexJRso37SbHgV5HNu/O2eXFPPsm1v4+5a95OYYdfXOT68qZfOu/eyrruX3r2/mZ1efSr82hmj319SxsPxdxg3tTb0nhnROHNyrQyMNXcGy7UtLpaWlXlZW1u71Xt+wk2lz/sLYIb2YP+usxvbd+2sYe0di7P3Hny3lC//d9LXPGNWPuV9o/Yi/4ZNBc3OuGM/UsQMBWL99H2ffu6hx2ZP/dhYf/8HitOpe881LeHPz7rT7TxjRl8e/eHqTtu8/t5r7nvk7kPq8RW1dPU+teJd6d6aNa3GX61hp+AVdX++Uvb2DjwzsQb1Dbo6xctNu7nnqTZau39livQeuHM+Ywb345cvr6X90PkP6dGfKiQMaf9Ye+swpTD7hGHZ9UNNmGGSz2rp6HNixr5rrH3+dF1dvA2Dq2IH84Y3Nba77pXNHMaBnId3zc7nz9yvZc6D2MFQME0f25ZpzRnHGqP7k5hgf1NRRlJ/bOKxTU1dPt3YMvx5pzOxVdy9N2S8qob9sw04uTRL6AOPveob32rgdQ1tBmeyo/Lczz2Tc0N4p+6USPvl70Xf/TMWWPWmve/WZw/mXM0fwsdmLWiz70/XnMKJ/ESNuPvjJYe4XTuOKH7/c5ms+e93ZvLZ+JycO6sXoQT3TruVIUFfvTLj72YzeluOcDxfzwt+rGN6vO29tTxyl/vlrkxjWr3uTfvX1ztY9B8jNMYp7FFBX73znjxWcMLAnX350KT0K89izP70g/dnVpzJ6YE96de/G9r3VDOp9VOO30+ct2cDNv1netRvZTkcX5HFMzwLuvWwsdz35Jvm5OazYtIt91U1v0TV6YE9Wbj44fDKgZyElHzqan3+uxaM9YkuhH5IqkNfdc0mrd+U8/tanOBBclnn3p07kUycPTnpG3t2bhGyDP3z5LKZ+/+BR/Pemj+PjYwc1Odnc4NW3d5Bj8KkH/grQ+HEyU9r6d8kmDfv+UPjWp8dwVkl/nl7xLr9d9k6Lcd64eODK8Tz6yvrGI/6F/342xw/o0bjc3Xlt/U56HZXHvQsrWFi+BYCi/FyOH9CD719+MgN7HUVtfT0FeZkZ1oi6dEM/8mP6AJOOL2ZRRVWryx94fg0zJx2XdFlD4N/xidFcedqxrb6GmXFMjwK27jkAJD4Gz7liPJAIzwO19SnH8E45tg8Ai2+cxFnfXsSvrjmd8cP6sOuDGgq75VCQl8uU+//MqnebfiI4eVhvfvOlM3i+ooqrf7akybLvTR/HtfOWNc6PKi7iuevPZee+aq78yct8ZuKx3NTK0d6ImxdQ8Y0pKf+THqit4/5nV/Oj4EqH1n5ZrNv2Pi+t2c4/njoUA3JyjJ37qllUsZUxg3tx22/LuWHK8WzbW82+6lp2vF9NTo7x08XreOFrk1i/fR+TvvM8j82YyG+WvsPc4ORbe3xv+jg+MXYQ1cFH/WS/fPfX1FGQl4M7jSfWP/+xkXz+YyNbfd2Gg6dk211TV0/JLU8BiStAhvcv4q9rtre7doD+R+czflgfBvQq5F/PPY7FlduYOmYgB2oTV42cMKAn7+7eT3GPAk78/wvTes25XziNM0b1b9KW7KKCS8YMbPU1zKzx5/fBz7SeO7k5CvxMi8yRfmtj+gAb3tuXdBgkrGGIZ1HFVq5+ZAn3fHpMk4++6V7jv+uDGo7qlpvy5HBn7dlfw4Hael5e+17juQVIBPDxtyauMDppaG9+147Lv55duYVTju1DTo5x0tf/2GSZGay7Z2q7hrGW33EhW/ccYGT/Im789Rs8XrYx7XU74rnrz+G80PmXr110PJ87a0TGTpi1R8P5BUj80gnXHF7WWdW19XTLtSPiE5y0T+yGd9oKfUg9xHPvZWM5dXhfzv2P51ss69O9G0tvv7DdNR3JauvqOS44Os12K75+EUcXxOJDq0irYje8095fXf82+Th+8KfKxvmvhS67bG7HvpoOVnXkysvN4a1vTWXH+9WcfNczTZZ9ZuKx/PxvbzfOP/LPpzIpuKPpMyu3tLhCKuytb01tHAo5UFtPQV6OjjpFDqPIhH6D1uJj0VfPZVLoKH7mpOPoWdiNuxe8eVjqOlL1KcpvHNr6xwdf4orThjFt3GDuuvTEpP0vGP2hxv4HauuoqfMWR+ENIX8kDLuIRE10L1ptZkT/oibzhd1y+cLZrZ+Yu+vSE7n54hMOdVlHlMe+eHq7ru0vyMvVsItIlonc/8i2hnlW3TWFE257mk+eNKixbdak4/jhosom/Zbccj7FwdO3vnjOqENRpohIRsTmSB8SR/ev3HIe3/m/JzW2XXt+CUP7HsU/TRwGwGkj+tKvo3e5ExHJcpE70k91SvCYHoVN5rvl5vDiDZMB+MalYw5RVSIi2SFWR/oiInGn0BcRiZHIhH62fclMRCQbRSb0G+mLPiIirYpe6OuIX0SkVZEJfX2VX0QktciEvsb0RURSi0zoN9IRv4hIq6IX+iIi0iqFvohIjKQV+mY2xcwqzKzSzG5Ksvw6M1tpZm+Y2XNmdmxo2WwzKzezN83s+3aIzrg2PPau4BA/sUpE5EiWMiHNLBeYA1wMjAYuN7PRzbotBUrdfSzwBDA7WPcM4ExgLHAicCpwTpdVHzJmcC++fF4J359+8qF4eRGRSEjnsHgCUOnua929GpgHTAt3cPdF7r4vmP0bMKRhEVAI5AMFQDdgS1cU3pyZcd0FH2ZAr8LUnUVEYiqd0B8MbAjNbwzaWvM54CkAd38JWARsDv4sdPcWj6oysxlmVmZmZVVVVenWLiIi7ZRO6Ccbg096UbyZ/RNQCtwbzB8HfITEkf9gYLKZnd3ixdwfcvdSdy8tLi5Ot3YREWmndEJ/IzA0ND8E2NS8k5mdD9wCfNLdDwTNnwL+5u573X0viU8AEztXsoiIdFQ6ob8EKDGzEWaWD0wH5oc7mNnJwIMkAn9raNF64BwzyzOzbiRO4upJ5CIiGZIy9N29FpgFLCQR2I+7e7mZ3Wlmnwy63QscDfzKzJaZWcMvhSeANcBy4HXgdXf/fVdvhIiIpMey7Z41paWlXlZWlukyRESOKGb2qruXpuqnbzKJiMSIQl9EJEaybnjHzKqAtzu4en9gWxeWcyTQNseDtjkeOrPNx7p7ymvesy70O8PMytIZ04oSbXM8aJvj4XBss4Z3RERiRKEvIhIjUQv9hzJdQAZom+NB2xwPh3ybIzWmLyIibYvakb6IiLRBoS8iEiORCf1Uj3TMZmY21MwWBY+ULDeza4P2vmb2jJmtDv7uE7Rb8OjJyuARleNDr3VV0H+1mV0Vaj/FzJYH6xyyx1a2l5nlmtlSM3symB9hZi8H9T8W3OQPMysI5iuD5cNDr3Fz0F5hZheF2rPuZ8LMepvZE2a2Ktjfp0d9P5vZV4Kf6xVm9qiZFUZtP5vZw2a21cxWhNoO+X5t7T3a5O5H/B8gl8SN3UaSeErX68DoTNfVjvoHAuOD6R7A30k8mnI2cFPQfhPw7WD6EhK3qTYSt6p+OWjvC6wN/u4TTPcJlr0CnB6s8xRwcaa3O6jrOmAu8GQw/zgwPZj+T+BLwfS/Av8ZTE8HHgumRwf7uwAYEfwc5GbrzwTwX8Dng+l8oHeU9zOJ52isA44K7d9/jtp+Bs4GxgMrQm2HfL+29h5t1prp/wRd9A9+OomncjXM3wzcnOm6OrE9vwMuACqAgUHbQKAimH4QuDzUvyJYfjnwYKj9waBtILAq1N6kXwa3cwjwHDAZeDL4gd4G5DXfryTu8np6MJ0X9LPm+7qhXzb+TAA9gwC0Zu2R3c8cfPJe32C/PQlcFMX9DAynaegf8v3a2nu09ScqwzvtfaRj1go+zp4MvAx8yN03AwR/HxN0a21722rfmKQ90+4HbgDqg/l+wE5P3M4bmtbZuG3B8l1B//b+W2TSSKAKeCQY0vqJmRUR4f3s7u8A/0Hi2RqbSey3V4n2fm5wOPZra+/RqqiEftqPdMxmZnY08Gvg3919d1tdk7R5B9ozxsw+Dmx191fDzUm6eoplR8w2kzhyHQ/8yN1PBt4n8ZG8NUf8NgdjzNNIDMkMAoqAi5N0jdJ+TiWj2xiV0E/rkY7ZzBJPFvs18Et3/03QvMXMBgbLBwINTyVrbXvbah+SpD2TzgQ+aWZvAfNIDPHcD/Q2s7ygT7jOxm0LlvcC3qP9/xaZtBHY6O4vB/NPkPglEOX9fD6wzt2r3L0G+A1wBtHezw0Ox35t7T1aFZXQT/lIx2wWnIn/KfCmu98XWjQfaDiDfxWJsf6G9s8GVwFMBHYFH+0WAheaWZ/gCOtCEuOdm4E9ZjYxeK/Phl4rI9z9Zncf4u7DSeyvP7n7lcAi4LKgW/Ntbvi3uCzo70H79OCqjxFACYmTXln3M+Hu7wIbzOz4oOk8YCUR3s8khnUmmln3oKaGbY7sfg45HPu1tfdoXSZP8nTxSZRLSFz1sga4JdP1tLP2s0h8XHsDWBb8uYTEWOZzwOrg775BfwPmcPBRlKWh1/oXoDL4c3WovRRYEazzQ5qdTMzw9p/Lwat3RpL4z1wJ/AooCNoLg/nKYPnI0Pq3BNtVQehqlWz8mQDGAWXBvv4tias0Ir2fga8Dq4K6fk7iCpxI7WfgURLnLGpIHJl/7nDs19beo60/ug2DiEiMRGV4R0RE0qDQFxGJEYW+iEiMKPRFRGJEoS8iEiMKfRGRGFHoi4jEyP8Cv6ayMm0R5CAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "start = 25\n",
    "x = np.array(range(len(errors)))[start:] + 1\n",
    "plt.plot(x, means[start:])\n",
    "plt.axhline(y=0, c='r')\n",
    "plt.show()\n",
    "plt.plot(x, variances[start:])\n",
    "plt.axhline(y=expected_error_variance, c='r')"
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}