HW1_Overfitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mathematical exercises\n",
    "\n",
    "Solutions to the following exercises may be either typed into the submitted Jupyter notebook, or handwritten solutions may be scanned and submitted as a separate PDF file.\n",
    "\n",
    "1. Exercise 1.6 from the PRML textbook (3 points).\n",
    "2. Exercise 1.12 from the PRML textbook (3 points)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coding exercises\n",
    "\n",
    "In class, we generated some fake data and fit polynomials to the data. We noticed that higher-order polynomials did a good job of fitting the data, but a poor job of fitting new data that was generated in the same way. This is an example of **overfitting**.\n",
    "\n",
    "In this assignment, we will use **cross validation** to combat overfitting by selecting a model that optimally describes data that was *not* used to fit the model.\n",
    "\n",
    "Reminder: Once your notebook is complete, you should restart the kernel, run all of the cells in order, and make sure that no errors occur. Notebooks that run successfully without errors will be awarded one extra credit point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, as we did in class, let's make some fake data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEICAYAAABBBrPDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVNUlEQVR4nO3de4xcd3nG8edhY2AJlA14Feyxg902sgq4YNgGqCsUhYudNI1dc2mCSiACWa0IhBa5tUEiLapkt1ZRW0BQy0lJWiBEiVncktYFHBQuDfU6m2BsY2oiUnsSmiXBgcBCbOftH3M2Xq9nzs56z5zfmdnvR1p5zsWz7yjKPP5djyNCAAC08rTUBQAAqo2gAADkIigAALkICgBALoICAJCLoAAA5EoWFLYX277T9gHb+21f1+Qe2/4H24dtf9v2y1PUCgBz2TkJf/cJSe+PiHtsP0fSXttfiogDk+65VNKF2c8rJX0i+7Ol+fPnx5IlSzpUMgD0pr179/4oIgabXUsWFBHxkKSHstc/tX1QUk3S5KBYI+nmaKwKvNv2gO0F2d9tasmSJRoZGelk6QDQc2w/0OpaJcYobC+RtELSt6Zcqkk6Mun4aHZu6t9fb3vE9sjY2FjH6gSAuSh5UNh+tqTbJb0vIn5yNu8REdsiYigihgYHm7acAABnKWlQ2J6nRkh8OiJ2NLmlLmnxpONF2TkAQElSznqypBskHYyIj7S4baekq7PZT6+S9Fje+AQAoHgpZz2tlPQ2Sfts35ud+4CkCyQpIj4p6Q5Jl0k6LOnnkq4pv0wAmNtSznr6uiRPc09Ienc5FaUxPFrX1l2H9OCxcS0c6NeGVcu0dsUZ4/UAkEzKFsWcNzxa16Yd+zR+/KQkqX5sXJt27JMkwgJAZSSf9TSXbd116KmQmDB+/KS27jqUqCIAOBNBkdCDx8ZndB4AUiAoElo40D+j8wCQAkGR0IZVy9Q/r++0c/3z+rRh1bJEFQHAmRjMTmhiwJpZTwCqjKBIbO2KGsEAoNLoegIA5CIoAAC5CAoAQC6CAgCQi6AAAOQiKAAAuQgKAEAuggIAkIugAADkIigAALkICgBALoICAJCLoAAA5CIoAAC5CAoAQC6CAgCQi6AAAOQiKAAAuQgKAECupEFh+0bbD9v+TovrF9t+zPa92c+Hyq4RAOa6cxL//k9J+pikm3Pu+VpEXF5OOQCAqZK2KCLiLkmPpqwBAJCvG8YoXm37Ptv/bvvFzW6wvd72iO2RsbGxsusDgJ5W9aC4R9ILI+Klkj4qabjZTRGxLSKGImJocHCwzPoAoOdVOigi4icR8Xj2+g5J82zPT1wWAMwplQ4K2y+w7ez1RWrU+0jaqgBgbkk668n2ZyVdLGm+7aOSrpc0T5Ii4pOS3iTpj22fkDQu6cqIiETlAsCclDQoIuKqaa5/TI3pswCARCrd9QQASI+gAADkIigAALlSb+GBNg2P1rV11yE9eGxcCwf6tWHVMq1dUUtdFoA5gKDoAsOjdW3asU/jx09KkurHxrVpxz5JIiwAdBxBkcBMWwdbdx16KiQmjB8/qa27DhEUADqOoCjZ2bQOHjw2PqPzAFAkBrNLltc6aGXhQP+MzgNAkQiKkp1N62DDqmXqn9d32rn+eX3asGpZobUBQDMERcnOpnWwdkVNm9ctV22gX5ZUG+jX5nXLGZ8AUArGKEq2YdWy08YopPZaB2tX1AgGAEkQFCWb+LJnTQSAbkFQJEDrAEA3YYwCAJCLoAAA5CIoAAC5CAoAQC6CAgCQi1lPKBTboQO9h6BAYdgOHehNdD2hMGez4SGA6iMoUBi2Qwd6E0GBwrAdOtCbCAoUhu3Qgd7EYDYKw4aHQG9KGhS2b5R0uaSHI+IlTa5b0t9LukzSzyW9IyLuKbdKzAQbHgK9J3XX06ckrc65fqmkC7Of9ZI+UUJNAIBJkgZFRNwl6dGcW9ZIujka7pY0YHtBOdUBAKT0LYrp1CQdmXR8NDt3GtvrbY/YHhkbGyutOACYC3piMDsitknaJklDQ0ORuBx0IbYeAVqrelDUJS2edLwoOwcUhq1HgHxV73raKelqN7xK0mMR8VDqotBb2HoEyJd6euxnJV0sab7to5KulzRPkiLik5LuUGNq7GE1psdek6ZS9DK2HgHyJQ2KiLhqmush6d0lldM16E8v1sKBftWbhAJbjwANVe96whQT/en1Y+MKnepPHx5l6OZssfUIkI+g6DL0pxdv7YqaNq9brtpAvyypNtCvzeuW00oDMlWf9YQp6E/vDLYeAVojKCqs2VgE/ekAykZQVFSruf1vfEVNt++tn9b9VHZ/OoPpwNzCGEWbhkfrWrllt5Zu/KJWbtnd8cHjVmMRd353LGl/OoPpwNxDi6INKVbu5o1FpOxPzxtMp1UB9CZaFG1IMdOoqo8VZTAdmHsIijak+HKs6tz+qgYYgM4hKNqQ4suxqnP7qxpgADqHMYo2bFi17LQxCqmcL8cqzu3nudjA3ENQtIEvx9NVMcAAdA5B0Sa+HAHMVYxRAAByERQAgFwEBQAgF0EBAMhFUAAAcjHrqQuwWyuAlAiKikuxISEATEbXU8Xx6FMAqREUFcdurQBSIygqjt1aAaRGUFQcu7UCSI3B7IpjQ0IAqSUNCturJf29pD5J2yNiy5Tr75C0VdLEA5k/FhHbSy2yAtiQEEBKyYLCdp+kj0t6vaSjkvbY3hkRB6bc+rmIuLb0AgEAktKOUVwk6XBE3B8RT0i6RdKahPUAAJpI2fVUk3Rk0vFRSa9sct8bbb9G0vck/UlEHJl6g+31ktZL0gUXXNCBUlFlrFwHOqvqs57+VdKSiPhNSV+SdFOzmyJiW0QMRcTQ4OBgqQUirYmV6/Vj4wqdWrk+PFqf9u8CaE/KoKhLWjzpeJFODVpLkiLikYj4ZXa4XdIrSqoNXYKV60DnTRsUtt9j+7wO/O49ki60vdT20yVdKWnnlN+9YNLhFZIOdqAOdDFWrgOd106L4nw1ZiTdanu1bRfxiyPihKRrJe1SIwBujYj9tj9s+4rstvfa3m/7PknvlfSOIn43egcr14HOc0RMf1MjHN4g6RpJQ5JulXRDRHy/s+XN3NDQUIyMjKQuAyWZuruu1Fi5vnndcga0gRmwvTcihppda2uMIhpp8sPs54Sk8yTdZvtvCqsSOAtrV9S0ed1y1Qb6ZUm1gX5CAijYtC0K29dJulrSj9QYUB6OiOO2nybpfyLi1zpfZvtoUQDAzOW1KNpZR/E8Sesi4oHJJyPiSduXF1EgAKC6pg2KiLg+5xqzkACgx7F7bIbVvQDQHEEhnksNAHmqvoVHKVjdCwCtERRidS8A5CEoxOpeAMhDUIjnUgNAHgazxXOpASAPQZHhudQA0BxdTwCAXAQFACAXQQEAyEVQAAByERQAgFzMemqBTQIBoIGgaIJNAgHgFLqemmCTQAA4hRZFE2wSWF10CQLlo0XRBJsEVtNEl2D92LhCp7oEh0frqUsDehpB0QSbBFYTXYJAGnQ9NcEmgdVElyCQBkHRApsEVs/CgX7Vm4QCXYJAZ9H1hK7Rbpfg8GhdK7fs1tKNX9TKLbsZwwBmKWlQ2F5t+5Dtw7Y3Nrn+DNufy65/y/aSBGWiItauqGnzuuWqDfTLkmoD/dq8bvlpLT8GvIHiJet6st0n6eOSXi/pqKQ9tndGxIFJt71T0o8j4tdtXynpryX9QfnVoiqm6xLMG/CmKxE4OylbFBdJOhwR90fEE5JukbRmyj1rJN2Uvb5N0mttu8Qa0WUY8AaKlzIoapKOTDo+mp1rek9EnJD0mKTnT30j2+ttj9geGRsb61C56AasgQGK1xOD2RGxLSKGImJocHAwdTlIiDUwQPFSTo+tS1o86XhRdq7ZPUdtnyPpuZIeKac8dCPWwADFSxkUeyRdaHupGoFwpaS3Trlnp6S3S/ovSW+StDsiotQq0XVYAwMUK1lQRMQJ29dK2iWpT9KNEbHf9ocljUTETkk3SPpn24clPapGmAAASpR0ZXZE3CHpjinnPjTp9S8kvbnsugAAp/TEYDYAoHMICgBALoICAJCLoAAA5CIoAAC5CAoAQC6CAgCQi6AAAOQiKAAAuQgKAEAuggIAkCvpXk8Apjc8WmfbdCRFUAAVNjxa16Yd+556Dnj92Lg27dgnSYQFSkPXE1BhW3cdeiokJowfP6mtuw4lqghzEUEBVNiDx8ZndB7oBIICqLCFA/0zOg90AkEBVNiGVcvUP6/vtHP98/q0YdWyRBVhLmIwG6iwiQFrZj0hJYICqLi1K2oEA5Ki6wkAkIsWBTBHsHAPZ4ugAOYAFu5hNuh6AuYAFu5hNmhRAD1oajdTnYV7mAWCAugxzbqZLCma3MvCPbQjSVDYfp6kz0laIukHkt4SET9uct9JSfuyw/+NiCvKqhHoVs26mUI6IyxYuNc7Oj1RIdUYxUZJX4mICyV9JTtuZjwiXpb9EBJAG1p1J4Wk2kC/nP25ed1yBrJ7wEQLsn5sXKFTExWGR+uF/Y5UXU9rJF2cvb5J0lcl/XmiWoCe0mpMojbQr29svCRBReikvIkKRf1DIFWL4vyIeCh7/UNJ57e475m2R2zfbXttqzezvT67b2RsbKzoWoGuwv5Qc0sZOwx3rEVh+8uSXtDk0gcnH0RE2G42ziZJL4yIuu1flbTb9r6I+P7UmyJim6RtkjQ0NNTqvYA5gf2h5pZWLcgiJyp0LCgi4nWtrtn+P9sLIuIh2wskPdziPerZn/fb/qqkFZLOCAoAp2N/qLljw6plp81yk4pvQabqetop6e3Z67dL+sLUG2yfZ/sZ2ev5klZKOlBahQDQBdauqGnzuuUdnaiQajB7i6Rbbb9T0gOS3iJJtock/VFEvEvSb0j6R9tPqhFoWyKCoACAKTrdgkwSFBHxiKTXNjk/Iuld2etvSlpecmkAgCnY6wkAkIugAADkIigAALkICgBALnaPBdCzeKpfMQgKAD3pbJ/q18lw6dbgousJQM8ZHq3r/bfeN+On+nVyJ9YydnntFFoUAHrKxBfyyWi+7VveZnlF7sQ6tfXws1+eyH3vKrc2CAoAPaXZl/1keZvlFbUTa7Nur7zfebbdZGWh6wlAT8n7Up9us7xWITLTnVinC6up753XkqkCggJAT2n1pd5nT7tZXlHP8mi3BTLx3mU8U2I2CAoAPaXVl/3fvuWl03bjFLUTa6uwOu9Z85q+d1EtmU5hjAJAT5ntg5uK2Im11TMirv+9Fzd97zKeKTEbBAWAnpP6wU0zDauqP5XQ0WIKWbcaGhqKkZGR1GUAQKE6PX3W9t6IGGp2jRYFAFRc6umzDGYDQMWlnj5LUABAxbVasFfW9FmCAgAqbHi0Lre4Vtb0WYICACps665DajblyFJp02cZzAaAWej0bKRW3Uuh8vaBokUBAGepjK3DW3Uv1UpctU2LAgDaNLX18PMn8rcOL0IVVm0TFADQhpluHV6UKqzaJigAoA0z3Tq8SKm3JEkyRmH7zbb3237SdtMl49l9q20fsn3Y9sYyawSAydptJViN1sbKLbu74jGn7Ug1mP0dSesk3dXqBtt9kj4u6VJJL5J0le0XlVMeAJyuVSthoH/eUwPLlp6aytpNz8SeTpKgiIiDETHd2vOLJB2OiPsj4glJt0ha0/nqAOBMrZ5z8RdXvFjf2HiJagP9Z6x3qNJT6majytNja5KOTDo+mp07g+31tkdsj4yNjZVSHIC5ZbqHGlX9KXWz0bHBbNtflvSCJpc+GBFfKPJ3RcQ2SdukxjbjRb43AEzIG1ReONDfdCZUVZ5SNxsdC4qIeN0s36IuafGk40XZOQConCqsd+iUKk+P3SPpQttL1QiIKyW9NW1JANBcFdY7dEqSoLD9+5I+KmlQ0hdt3xsRq2wvlLQ9Ii6LiBO2r5W0S1KfpBsjYn+KegGgHanXO3RKkqCIiM9L+nyT8w9KumzS8R2S7iixNADAFFWe9QQAqACCAgCQi6AAAOQiKAAAuRzRW+vTbI9JekDSfEk/SlxO0fhM3aMXPxefqTuc7Wd6YUQMNrvQc0ExwfZIRLTcmbYb8Zm6Ry9+Lj5Td+jEZ6LrCQCQi6AAAOTq5aDYlrqADuAzdY9e/Fx8pu5Q+Gfq2TEKAEAxerlFAQAoAEEBAMjVk0Fhe7XtQ7YP296Yup7Zsn2j7Ydtfyd1LUWxvdj2nbYP2N5v+7rUNc2W7Wfa/m/b92Wf6S9T11QU2322R23/W+paimD7B7b32b7X9kjqeopie8D2bba/a/ug7VcX8r69NkZhu0/S9yS9Xo3Hp+6RdFVEHEha2CzYfo2kxyXdHBEvSV1PEWwvkLQgIu6x/RxJeyWt7fL/TpZ0bkQ8bnuepK9Lui4i7k5c2qzZ/lNJQ5J+JSIuT13PbNn+gaShiOipxXa2b5L0tYjYbvvpkp4VEcdm+7692KK4SNLhiLg/Ip6QdIukNYlrmpWIuEvSo6nrKFJEPBQR92SvfyrpoFo8E71bRMPj2eG87Kfr/yVme5Gk35W0PXUtaM32cyW9RtINkhQRTxQRElJvBkVN0pFJx0fV5V9Avc72EkkrJH0rcSmzlnXR3CvpYUlfioiu/0yS/k7Sn0l6MnEdRQpJ/2l7r+31qYspyFJJY5L+Kesm3G773CLeuBeDAl3E9rMl3S7pfRHxk9T1zFZEnIyIl6nxjPeLbHd1V6HtyyU9HBF7U9dSsN+JiJdLulTSu7Pu3W53jqSXS/pERKyQ9DNJhYzR9mJQ1CUtnnS8KDuHisn68W+X9OmI2JG6niJlTf47Ja1OXMpsrZR0Rdanf4ukS2z/S9qSZi8i6tmfD6vxtM2L0lZUiKOSjk5qxd6mRnDMWi8GxR5JF9pemg3mXClpZ+KaMEU28HuDpIMR8ZHU9RTB9qDtgex1vxoTKr6btKhZiohNEbEoIpao8f/S7oj4w8RlzYrtc7MJFMq6Zt4gqetnFEbEDyUdsb0sO/VaSYVMDknyzOxOiogTtq+VtEtSn6QbI2J/4rJmxfZnJV0sab7to5Kuj4gb0lY1ayslvU3SvqxPX5I+kD0nvVstkHRTNvPuaZJujYiemE7aY86X9PnGv1V0jqTPRMR/pC2pMO+R9OnsH8n3S7qmiDftuemxAIBi9WLXEwCgQAQFACAXQQEAyEVQAAByERQAgFwEBQAgF0EBAMhFUAAdZvu3bH87e17FudmzKrp6DyjMLSy4A0pg+68kPVNSvxr78WxOXBLQNoICKEG2pcIeSb+Q9NsRcTJxSUDb6HoCyvF8Sc+W9Bw1WhZA16BFAZTA9k41tuleqsYjYK9NXBLQtp7bPRaoGttXSzoeEZ/Jdpb9pu1LImJ36tqAdtCiAADkYowCAJCLoAAA5CIoAAC5CAoAQC6CAgCQi6AAAOQiKAAAuf4fgr9ytbQJ0b4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 30  # number of data points\n",
    "\n",
    "# Simulated data:\n",
    "x_data = 2*np.pi*np.random.rand(n)\n",
    "noise_amp = 0.3\n",
    "noise = noise_amp*np.random.randn(n)\n",
    "y_data = np.sin(x_data) + noise\n",
    "\n",
    "plt.plot(x_data, y_data, 'o')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In class, we used the function `np.polyfit()` to fit the data with a polynomial:\n",
    "$$\n",
    "\\hat{y} = \\beta_0 + \\beta_1 x + \\ldots + \\beta_{p-1} x^{p-1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Polynomial coefficients: \n",
      " [ 1.68902328e+01 -1.76281580e+02  7.16715959e+02 -1.48216770e+03\n",
      "  1.80647827e+03 -1.40422882e+03  7.27238033e+02 -2.56333095e+02\n",
      "  6.16864841e+01 -9.96968893e+00  1.03447185e+00 -6.22407800e-02\n",
      "  1.65033550e-03]\n"
     ]
    }
   ],
   "source": [
    "# A polynomial fit:\n",
    "p = 12\n",
    "fit = np.polyfit(x_data, y_data, p)\n",
    "beta_list = fit[::-1]  # flip the list so that it starts with beta0\n",
    "print('Polynomial coefficients: \\n', beta_list)\n",
    "\n",
    "\n",
    "# The following function gives us the y value of our fitting curve given the x value:\n",
    "def yhat(x, beta_list):\n",
    "    '''\n",
    "    Given the x value and the list of polynomial coefficients, return the y value.\n",
    "    '''\n",
    "    y = 0\n",
    "    for p in range(len(beta_list)):\n",
    "        y = y + beta_list[p]*x**p\n",
    "        \n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwZklEQVR4nO3deXxU9bnH8c+TZAhhkbAqBAVUFqWIYNQqrlQFcQF3rApupe0tapfLLWhvq92kF5fa1qoobnVfECkiuGBVtCphURYBAbUSQMIuEMj2u3/8ZjBAMpkkM3Nmku/79cprMmfOnPOMyzz5bc/PnHOIiIhUJyPoAEREJLUpUYiISFRKFCIiEpUShYiIRKVEISIiUSlRiIhIVIElCjM72MzeMrMlZrbYzG6q4hwzs7+Y2Qoz+8TM+gcRq4hIY5YV4L3LgF845+aZWUtgrpm97pxbUumcs4Hu4Z/jgfvCjyIikiSBtSicc2udc/PCv38DfArk7XPaUOBx530A5JpZxySHKiLSqAXZotjDzLoC/YAP93kpD/iq0vPV4WNrq7jGKGAUQPPmzY/p1atXfIMs2wXrP4XW3SAnN77XFhGJZsNnYAZtD0/YLebOnbvBOde+qtcCTxRm1gJ4Efipc25bXa/jnJsITATIz893BQUFcYowbONK+Gt/uOD30Hd4fK8tIhLNgwOhaS5cNTlhtzCzL6t7LdBZT2YWwieJJ51zVf0TKAQOrvS8c/hY8mVl+8fykkBuLyKNWHkJZDYJ7PZBznoyYBLwqXPurmpOmwqMCM9++i6w1Tm3X7dTUkT+JZXtDuT2ItKIlZdCZiiw2wfZ9TQAuApYaGYLwsduBg4BcM7dD0wHhgArgJ3ANckPMyySKMpLA7n9lPmFTJi5jDVbiumUm8OYQT0Z1m/fsX8RaZACblEEliicc7MBq+EcB/wkORHVYE+iSH6LYsr8QsZNXkhxaTkAhVuKGTd5IYCShUhjUF4WaItCK7NjFRmjKEv+GMWEmcv2JImI4tJyJsxclvRYRCQA5SVKFGkhIxMsI5DB7DVbimt1XEQamMY6mJ2WMrMD6XrqlJtTq+Mi0sCUlypRpI2sJoEMZo8Z1JOcUOZex3JCmYwZ1DPpsYhIAALuegp8wV1ayWwSyPTYyIC1Zj2JNFIVpZChRJEeMrMDW3A3rF+eEoNIY1RRDq5CXU9pIzOkldkiklyR7xzNekoTWdlamS0iybUnUahFkR4yQ4GtzBaRRqq8zD8qUaSJgKbHikgjpq6nNJOVHcjKbBFpxFIgUWjWU21khqBkR2KuveU/MO8f4Mqh/who3TUx9xGR9BLp7m6MRQHTUlZT2Lkp/tddtwgeHQK7t/tdrD64Hy6eBD3Pjv+9RCS9pECLQl1PtZHV1G+JGk8lO+CpyyDUHG4ogJs+hvY94PlrYM2C+N5LRNJPRfAtCiWK2khEonjvHti2Gi5+GNocCq06w/efh5zWMOW/NMtKpLFLga4nJYraCDWF0jgmit3fwL/vhSPOhy4nfHu8RXs45w5YvxjmTIrf/UQk/US6njKCGylQoqiNrKbxXXC38AUo2Q4n3rD/a73Oga4nw+y745ucRCS9aMFdmslqCmVx3ANi3mPQoTd0Prbq10/9JWxf588TkcZpT9dTIx3MNrOHzWy9mS2q5vXTzGyrmS0I//w62THuJaupz+4VFfW/1pavYM18OOpSP9OpKt1Ohi4D/DhGZHWmiDQukV6MyC6bAQi6RfEoMLiGc951zh0d/vltEmKqXqipf4zHgPay6f6x17nRzzvhJ7CtEJbPqP89RST9RKpBZDbSROGcewdIwMKEBMmKY6JY+gq06wntDo9+XvdBcEBnmPNQ/e8pIuknUg0iS2MU0ZxgZh+b2atm1jvQSOKVKEqL4T//hu5n1nxuZhbkXw2r3oINK+p3XxFJP429RRGDeUAX51xf4K/AlOpONLNRZlZgZgVFRUWJiSZeieKrj/xYR9eTYzu/3wiwDFjwZP3uKyLpZ0+LQomiSs65bc657eHfpwMhM2tXzbkTnXP5zrn89u3bJyagyBhFfaerfjHbf/FXXjsRTcsD4bDvwSfPxWcgXUTSx54WhbqeqmRmB5n5KUFmdhw+3o2BBZSV4x/r26L4YjZ0PBqator9PX2H+xXcX86u371FJL3smfXUNLAQAi0KaGZPA6cB7cxsNfAbIATgnLsfuBj4sZmVAcXAcOecCyjcb5t+9UkU5WWwdgH0H1m79/U6B5q0hI+fgW6n1P3+IpJeIomisZYZd85dXsPrfwP+lqRwahaKQ4tiwzIo3Ql5/Wt/795DYfEUGHIHNGlW9xhEJH2U7/YD2dWtt0qClO56SjmRFkV9xigK5/nHTrVMFABHDfclP5a+Uvf7i0h6KSsJdCAbtB9F7cRjjGLNPMg+wFeKra0uA+CAPFj4PBx1yX4vT5lfyISZy1izpZhOuTmMGdSTYf3y6h6riASvfHegA9mgFkXtxGOMonAedDoaMurwjz4jA75zEax8E3bsPaY/ZX4h4yYvpHBLMQ4o3FLMuMkLmTK/sO6xikjwUqBFoURRG/UdoyjbDV8vjrnbacr8QgaMn0W3sa8wYPws/6V/1KVQUQaLJ+917oSZyyguLd/rWHFpORNmLqtbrCKSGtSiSDP1HaP4epHfrSqGgexqWwhrWkP7I3z3UyVrtlRd1ba64yKSJsp2Bzo1FpQoaqe+YxRfL/GPB/Wp8dRqWwivLffjE199CJu/2PNap9ycKq9T3XERSRPlJYHWeQIlitqp7xjF+k99ssntWuOpUVsIfcID2ZVaFWMG9SQnlLnXuTmhTMYM6lm3WEUkNZTtCrTOEyhR1I5Z/fbNXr8E2veMaSA7agsh9xA45AT45HkIrz8c1i+P2y/sQ15uDgbk5eZw+4V9NOtJJN2lwGC2psfWVlZ23cco1n8Khw2M6dQxg3oybvLCvbqf9moh9LkEXvk5rPsEOvYFfLJQYhBpYMp3Qyg30BDUoqitrJy6tSh2bvLbmnY4IqbTa2wh9L7Ab7b+yXO1j0VE0odaFGkoK7tuiaJoqX/scGTMb4naQmjWBg4/Exa9CGf+FjIyqz5PRNJb+e7AE4VaFLUVqmOLYn14xlOMLYqYHHUJfLPWV6MVkYapbLcGs9NOXcco1n8K2a3ggE7xi6XH2dCkBSxU95NIg6XpsWkoKwfK6rCIrWgZtO8R3wqQTZrBEefBkn/WfzOlBKpyhbmIxEbTY9NQKKduX8qbVkHbw+MfT59LYPdW+Oy1+F87DlSDSqSeUmAwW4mitpo0h5IdtXtPaTFsK6xbxdiadDsVmndI2e4n1aASqSfVekpDTZpDaS0TRaTURiISRWaWryi7fCYUb4n/9etJNahE6qGiwhcBVYsizdSlRbFxpX9MRKIAP/upvASWvJyY69eDalCJ1EN5ZL9sJYr0UpdEsWmVf0xUoujUH9r1hHmPJeb69aAaVCL1sGe/7EacKMzsYTNbb2aLqnndzOwvZrbCzD4xszrsHxpnoeZ+z+uKitjfs2klNGsLObmJickMjr0eCufC6rmJuUcdqQaVSD2Ul/jHgKfHBr0y+1Hgb8Dj1bx+NtA9/HM8cF/4MThNmvvH0p2Q3SK292xalbjWRETf4fDmbTDnQeh8TGLvVUuqQSVSR2pRgHPuHWBTlFOGAo877wMg18w6Jie6ajRp5h9r0/20cRW0OSwx8UQ0PcAni0Uvwo4Nib2XiCRHmcYoYpEHfFXp+erwsf2Y2SgzKzCzgqKiosRF1CTcioh15lNpMWxbnfgWBcCxP/BN1YJHEn+vBkgLAyXlRAazNT02PpxzE51z+c65/Pbt2yfuRpGup1hbFImcGruvDr2g+1nwwd9h9/bE368B0cJASUl7WhTaCjWaQuDgSs87h48FJxTpetoZ2/mRGU9tk5AoAE4ZA8WbYK5aFbWhhYGSkjSYHZOpwGgzewY/iL3VObc20IgiXU8lMf7Fnuipsfs6+Djodiq73r6b894+lBVb/ZqFMYN6akA5Ci0MlJSUIoPZgSYKM3saOA1oZ2argd8AIQDn3P3AdGAIsALYCVwTTKSV1HYwe+NKyGkDOa0TF9M+3u78I079/HLOL3uOO7l0TzcKEL9ksWsrLHsVvvoQNnwGu7+B7JbQugt0GQA9z07qZ66vTrk5FFaRFLQwUAIV2dIgFGzXU6CJwjl3eQ2vO+AnSQonNpWnx8YiGVNj93HznKaMKT+RUZmv8Gz56ax27fd0o9Q7UayZD7P/DMum+2Zxdiu/D3jz9r6VtXQ6zH/CV9ntexmcOhYOCHaiWixq3HpWJAil4T9esoL9gyXVu55ST126nrqcmLh48AOxE2YuY82W4j1/Gf+Jyzkrey5/zHqIkaW/xJFRv26UTZ/DzJt9gshu5Rf4feciyDtm79LpFRWwdoFfJb7gKb9V65m/9efHs8R6nEUSaOV/juquk8CpRZGmYhzMnjK/kHtmLOTNXat5eIfRrlthQr50IrN1In8JF24pxoC1tOUPZVfwh9DDXFXxOo+XD6pbN0ppMcy+27ciMkMw8Fdw3Cho2qrKWL79oh3Kb04fzllf3AHT/xtWzoILH4x9kWIAtDBQUk6KtChSfdZT6gnVPEYR+fLO3PYlGeb4ZGfbhE21rGq2jgMMeLL8e7xV3pdfZT3BKaGltetGcQ4+nQb3Hgdv/8lvkDR6jp9VVU2S2Hd66U0ztzLlyHtg0O2wfAY8cjZsC3YugkhaSZEWhRJFbWVk+GQRpesp8uXd1b4G4Et3YMKmWlbXneSAvNxm/LT0JxRmdGRS9t0M67Autotu+AyeuBCevcLXtho5DS6eFHUb12qnl762HE74L7j8WT+w//Ag2Kq1CSIxUYsijTVpHnUwO/Ll3dX8F/MX7qC9jsdTdd1Jebk5vDd2IB+Pv4xuN71KqEVreGSIHzdwruqLbf4CXvkF/P0EWF0Ag/8EP5oN3U6uMY4ap5f2OAtG/hOKN8Nj56llIRKLsl2AqYRHWqqh1Hjky7urrWOza8FWWux1PJ5iKuOdewhcPws69YMpP4YHT4cP7ocVb8CKN+HDB+CJi+Av/WHuY9DvCrhhLnz3R35jpBjEtO9E52Pgihdg+9e+xbJrW60/r0ijUlrsV2UHPBFEiaIuQtETReTLu4t9zZfuQCBxUy1jLuPdoj1c/Qqc9xc/ED/jlz45PHEhvPo/vrvpxBvgp5/AefdAiw61iiPmfScOOR4uewI2LIfnr4bystp/aJHGomxX4K0J0Kynuslu4ReYVSPyJX34y1/zYXkP8hI81TLm2ToZmXDMSP+z5SvYtsZvs9imG7TsWK+/Wmo1vfSw0+Gcu+CfN8KrY/zvKTx1ViQwpcUQCn7RpxJFXTRtBdvXRz1lWJ928PIGhg28jmGnD0xSYLWQe7D/iaNaTS89ZqRfY/Len+GgoyA/+EX3IimnbFfgBQFBXU910zQXdm2Jfs7mLwAHbRO8D0U6+95v4PAzfNdXYWrtzCeSElKkRaFEURdNW0HxlujnJLsYYDrKyPCL8FocBM+OgB0bg45IJLWoRZHGcnJh97bo+2ZvXOkflSiia9YGLnscdhTBi9dBRXnN7xFpLEp3qUWRtprmgquAkuoHtNm0yp/XrE1SQkrr3dk69YMhE2DVW/Cv24OORiR1lBWrRZG2cnL9Y7Tup00rk9aaaBC7sx0zEvpdCe9MgOUzg45GJDWoRZHGIrWOog1ob1qVtIHsBrM725A74KA+MHkUbP4y6Ghi4xxsWOEXLy6ZCp+/Czs3BR2VNBQp0qLQ9Ni6aJrrH3dtrfr1st2wdTX0jbrdRtw0mN3ZQjlw6ePwwGnw3Ai4dmbgxdCqtfYTmPsofDrVj6/sxaBzPhz7A1+KPcbV7SL7Kd2VEv8P6L/guqip62nzF34Mo+3hSQmnQe3O1uZQuOA+eOb7MGMsnPfnoCPa2/qlMOt3sHSa/0uv1znQ7VS/eVOoGezc4OtkLXweXhoF7/8Vzv8L5PUPOnJJR2XFgRcEBCWKuqmp62njCv+YpK6nBrc7W69zYMBN8N49cPDxcHRyWmZRlZX4fTnemeATxGnj4PgffftHQ2WHDYST/9u3NmaMhUlnwbl3Q/+rkh62pLkUaVEEOkZhZoPNbJmZrTCzsVW8frWZFZnZgvDP9UHEuZ9I11N1LYpIomiTnEQRc72ndDLw19DlJJj2M/h6cbCxFC2HBwfCv/4IvYfBTQvgtLFVJ4mIjAx/7o/fh64nwdTR8O6dyYlXGgbn1KIws0zgXuBMYDUwx8ymOueW7HPqs8650UkPMJrslmCZ0VsUzdtH/yKJswa3O1tmFlz8MDxwCjx7FYx6q8oNkxJuyVRfcTerKQx/GnoNqd37m7XxFXOn/Aje/K2/zgmptQ28pKiy3f6xkbcojgNWOOdWOedKgGeAoQHGEzsz/6VV3WD2xpVJG59o0FoeCJc8Alu+DFeaLY16elzXklSUwxu3wnNX+fGHH75T+yQRkZkFw+6HI4fCzFtg6fS6xyWNR1lqbFoEwSaKPOCrSs9Xh4/t6yIz+8TMXjCzaqvYmdkoMysws4Kion1noSRAtDIeGz5Tjad46XIinPtnv+f2Kz+vdtOluK4l2bHRl1+ffTcccw1c8yq0qmdrLTMLLngAOh0Nk3/gp9SKRFOaGtugQuqvo/gn0NU5dxTwOvBYdSc65yY65/Kdc/nt27dPfGTN2voZLvvatRV2rFeLIp76XwUn/wLmPQ7v3lHlKXFbS7JmPkw8Fb78N5z/Vz/rKl77AYRy4LInITPkk0UNLSRp5CL7ZTfyFkUhULmF0Dl8bA/n3EbnXLijjoeAY5IUW81aHgTfVLEHdaTGU9vuyY2noTv9V9DnUpj1e/j3vfu9HJe1JAuegkmDfKvl2leh/4i6Rlu9Vnl+BtSaeRrclujK1KIAmAN0N7NuZtYEGA5MrXyCmXWs9PR84NMkxhfdAZ3gmyr2fd6TKNSiiKuMDBh2X7if/2Z4/297vRzTVqzVKSvxe4VP+bHfge+Hb0NeAv8m6X0BHDUc3v4/WLcwcfeR9FaqMQqcc2XAaGAmPgE855xbbGa/NbPzw6fdaGaLzexj4Ebg6mCirULLg3w3U8nOvY9v/AwsA1p3DSSsBi0zCy6a5JPFa7fAq2P3VJuNeSvWfRUtg0lnwJyH/FawV74Ezdsl6hN86+zxflbc9DHVjrtIIxdJFClQ6ynQBXfOuenA9H2O/brS7+OAccmOKyYtw42db9buPXC9folfXZwCzcUGKTMEFz8Cr/0KPvg7FC2FYffVbitW8Ht1f/SAn7LapLkfOzji3OR9jpzWcMatMPUG+OQ56HtZ8u4t6aFkh39s0jzYONDK7LrbkyjW7Z0ovl4CBx4ZTEyNRUYmDL4d2vWAGePgvhPgrN8zrO/lNa8lcQ4+ew1e/7VPMj0G+0HrFh2SE3tlR1/p60W9/r9+NXp2i+THIKmrNHUSRarPekpZb672/+hunDj92zn7JTt91dgOvQOOrpHIvwZ+9K4fD3r5J/D378IH98GW/+x9nnO+/tYH98G9x8NTl/qBwuFPweXPBJMkwI+7DP4TbP/axyZSWaRFEWoWbByoRVEnU+YX8sc3NvBRJnSwzXvm7OduNk7DqUWRTO26w3Wvw5KX4f2/+NpKM8ZC8w5+wZ5zsK0Qijf78zv29esZel8IWU2CjR3g4GOh5zk+9mOvS9pGV5IGIuOfTYJvaSpR1MGEmctYX5rNzoxsDjT/BVRcWs5777/PaaAWRbKZ+bpKvYf56q6fv+1nE+0ogowsP4PpoD6+WF8qLoQc+Cu470S/wO+s3wUdjaSKku3+sYlaFGnJz8031rnWdLSNe44fsmsZNGsBbboFF1xj16GX/0knBx4JR10GH02EE2+EFklYMCqpr3Snn0GZAhsXaYyiDiJz8z93HTnMvl1L8d3QZ37DmozM6t4qUrVT/tsXgfvw/qAjkVRRsgNCzX2LOWBKFHUQmbP/mcvjUFtLBhW0C+3mcPcfOPi7QYcn6ahddzjiPJjzIOzaFnQ0kgpKdqTEjCdQoqiTyP4PG5t2I9tKOf6AzdxzUhlGhV/ZK1IXJ/3ML+Kc+2jQkUgqKNmREuMToERRZ8P65XHLD30toKcHw4CMReGB0/yAI5O0ldcfDj3N17KK7EUgjVfpTrUoGoR23f0GRctnwqKX/N7JTQ8IOipJZyf9DLavg4+fCToSCVrJdj9GkQKUKOrDzFc0XToNtv7HLwATqY9up8KBffwMKNWAatxK1KJoOE7+BfQ8h+Xdr2fAlGbx2V1NGi8zOH4UfL0Ivnw/6GgkSBqjaECat2VKrwkMXXomhVt31X93NZE+l/iigR89EHQkEqTSHSmxKhuUKOIibrurSa3EdY/sVBLKgX5XwafTYOvqoKORoJTsTIk6T6BEERdx2V1NaiWue2SnomOvBxwUPBx0JBIUraNoWOq1u5rUSYNvxbXuAj3O9msqSncFHY0kW0U5lBUrUTQkdd5dTeqsUbTijh8FOzfCkilBRyLJVhqpHKtE0WBEVmrn5eZgQF5uDrdf2KfmTXSkzhpFK67bqX63xLmPBR2JJFukxHiKjFEEWj3WzAYD9wCZwEPOufH7vJ4NPA4cA2wELnPOfZHsOGMxrF+eEkMSjRnUk3GTF+7V/bRvK27K/MLYt0ZNRWbQfwS8cSsULYf2PYKOSJJlT4nxRj7rycwygXuBs4EjgcvNbN8df64DNjvnDgfuBv6U3CglVdXUimswg91HX+FLw8xTq6JR2dP1lCYtCjO7AXjCObc5zvc+DljhnFsVvs8zwFBgSaVzhgK3hn9/AfibmZlzWrIq0Vtx0Qa706pV0aID9DwbPn4avvdryMoOOiJJhkgF4eyWwcYRFkuL4kBgjpk9Z2aDzeJWHD0P+KrS89XhY1We45wrA7YCbau6mJmNMrMCMysoKiqKU4iSrhrUYHf/q/2g9tJXgo5EkmV3JFGkRu24GhOFc+5XQHdgEnA18JmZ/dHMUmpPSefcROdcvnMuv3177RDW2DWowe7DTodWB6v7qTGJtCiatgo2jrCYxijCXT3rwj9lQGvgBTP7v3rcuxA4uNLzzuFjVZ5jZllAK/ygtkhUDWrKckYm9LsSVv0LNn0edDSSDOnWojCzm8xsLvB/wHtAH+fcj/EzkS6qx73nAN3NrJuZNQGGA1P3OWcqMDL8+8XALI1PSCwa3JTlflf6/ZPn/yPoSCQZdm31jymybUEs02PbABc6576sfNA5V2Fm59b1xs65MjMbDczET4992Dm32Mx+CxQ456biu7v+YWYrgE34ZCISkwY1ZblVZzj8DJj/JJx2M2QGOrNdEm33NsjMTpnJCzX+1+ac+02U1z6tz82dc9OB6fsc+3Wl33cBl9TnHiINRv+R8OwV8Nlr0GtI0NFIIu3aljKtCdDKbJH00WMQNO+g7qfGYPe2lBmfACUKkfSRGYKjv++33t22NuhoJJHUohCROus/Alw5fPxU0JFIIqlFISJ11vYw6HISzPsHVFQEHY0kiloUIlIv/UfA5s/hy9lBRyKJsnsbZKfGYjtQohBJP0ee779E5j0edCSSKGpRiEi9hHLgqEthyVTYuSnoaCTeKsqh5BuNUYhIzabML2TA+Fl0G/sKA8bP2rtEev8RUL4bFj4fXICSGLu/8Y9qUYhINDXup9HxKOh4tN/9TlVtGpYUq/MEShQiKSnafhp79B8B6xfDmnlJjk4Sak+dJw1mi0gUMe2n0ediyMrRoHZDszNcILtZlVvvBEKJQiQFxbSfRtNW0PsCWPgC7N6epMgk4ZQoRCQWMe+n0X8ElGyHJVOSF5wkVmQmmxKFiEQT834ah3wX2nZX91NDsmODf8xpHWwclaiovUiKimk/DTPofxW8/msoWgbt03AHP9nbzo0+SaTQniNqUYiku76XQ0bWfq2KqOswJHXt3JhS3U6gFoVI+mvRAXqeze65T3LWvFP4z9YyWuWE2FFSRmm5X2MRWYcBNJxd/xqqFEwUalGINADvtzqH7JLN9P5mNg7YUly6J0lE7LcOQ1LTzk1KFABm1sbMXjezz8KPVY7amFm5mS0I/0xNdpwi6eKX89ux2rVjZNZrUc+rbn1GYJyDtR/Du3fCiz+A50bCtJ/Dgqe/XXjWiEyZX8j6r9fw7OKdKdVdGFTX01jgTefceDMbG37+yyrOK3bOHZ3UyETS0OqtJTyWeRa3hJ7iSPuCJa5rledVtz4jECtnwRu3wdoF/nmrg33Bw2/WQcEkCDWHY66G08amVN2jRPFlWz5hQcY2NtMypboLg+p6Ggo8Fv79MWBYQHGINAidcnN4tvw0drpsRmZW3aqoch1GEHZt9a2Hf1wAxZtgyB0wZhX8bBGMngO//AKunwVHDoUP/g73HgdfNPy9NybMXEZG6Q6yrZSNriWQOt2FQSWKA51zkU1/1wEHVnNeUzMrMLMPzGxYtAua2ajwuQVFRUXxjFUk5Y0Z1JPSUCsml5/EsMz3aMM2QhlG62ah6Oswkm3jSph4Gix6EU67GUYXwHE/gOaV+uQzMqHzMXDBffCDN6FJC3h8KBQ8EljYybBmSzFtzBcE3EKLvY4HLWFdT2b2BnBQFS/dUvmJc86ZWXXlL7s45wrN7FBglpktdM6trOpE59xEYCJAfn6+ymlKoxJJAM+/OpQrS97khy3e5cBzbgk+MVRWOBeevBRcBVz9CnQ5oeb35B3jk8WL18O0n0JpMZzwXwkPNQidcnM4cOtmANa5NnsdD1rCEoVz7ozqXjOzr82so3NurZl1BNZXc43C8OMqM/sX0A+oMlGINHZ+gd5IePwlflj0Fhx1T9AhfatwHjw2FJq1gSsnQ7vDY39v01Zw2ZPw4rUwcxxkt/ClSxqYMYN68u5LbwGwNpwoUqW7MKiup6nAyPDvI4GX9z3BzFqbWXb493bAAGBJ0iIUSVff/TF8swYWTwk6Eu/rxfDEhT5JXDujdkkiIqsJXPwIHDYQpv0MVr4V/zgDNqxfHtf0aQLA165N6nQXElyiGA+caWafAWeEn2Nm+Wb2UPicI4ACM/sYeAsY75xTohCpyeFnQvteMPsuqKgINpZv1sET4XLoI6fCAZ3qfq3MEFzyKLTrAc+N8OMdDcx3WuyAJi1ZOP4S3hs7MCWSBASUKJxzG51z33POdXfOneGc2xQ+XuCcuz78+/vOuT7Oub7hx0lBxCqSdjIy4ORfwPolsPzV4OIo3QXPXOFnOV3xPLTuWv9rNm0F33/WD3g/P9LfoyHZVli/ZJogWpkt0hD1vtB/Mb8zIZitUp2Df94EhQVwwf1w0Hfid+3cQ2DY/bBuIcy8OX7XTQXb1ihRiEiSZGbBST+HNfNh5ZvJv//7f4VPnvFTYI88P/7X7zkYTrzRL8xbNDn+1w/KtjVwQGp0N1WmRCHSUPW93H/p/OtPyW1VLH/Nlz0/chic+j8Ju83Lba9jsfVgy/OjGXr7CzWWu0h4Nd2yEj9uEtlPorbKy2D7upRsUah6rEhDldUEThnj1x8snQZHnJf4exYthxev811Nw/7u98uIsynzC7ntn4vZvLOUbvZDpjcZx093/oX/mtwcqLrchS+PsZDi0nIgftV0p8wv5P4ZBVy+4wkuyXqHZoTHTDofC2f9AQ45PvaLfbPGrzFRohCRpOp3FXxwH7z+G+gx2M8cSpTizfD0cMhsAsOfhibN436Lfb/wP3cdGV92ObeFHmNo6RtMmJld5Rf/hJnL9rxnT7jh8hixJoop8wuZMHMZa7YU0yk3h9N7tWfJ3Nk8mvEn2mVuZUr5Scy13lzZO5tuXzxL9sODuaX0Ot5tOYQxg3rWfJ8Nn/nHdt1jiieZlChEGrLMLDjzNv8FPvdRXy4jEcrL4PlrYMt/4OppkHtwQm5T1Rf+4+VnclZGAf+b9Q/O3tq7yvdVVwYj1vIYVbVIln/4Go83+T+20IKhJb9jsesGwKtLQ2SV9eUO/sz40IP8ZFsO4yaXADW0XjYs94/tesQUUzJpjEKkoesxGLqeDLN+D9urLIJQfzPGwqq34Ny7/D7eCVLVF7sjgzGlP6ScDP7S9EGoKN/vnOrKYMRaHmPfBNXZiri/yV2sc224cPdte5IE+L1ANpRmM6r05xRU9GBC6AE6lf2n5uJ+G5ZD01xo3j6mmJJJiUKkoTODc+6E0p0wY1z8r//RgzDnQThhdMJLa1T3xb6GdtzuruZotwT+fe9+r48Z1JOcUOZex2pTHqNygsphFw+G7iSTCq4v/QXrqXI7HUoI8eOSm9hFiAmhB1i3ZUf0mxQt93ueJ2Bcp76UKEQag/Y9/SK8RS/4WUnxsuJNePWXvtVy5m/jd91qVPWFD5CbE+L4YaOh17kw63d+jUUlw/rlcfuFfcjLzalTNd1vE5TjjtD99LCvuKH0Br50Hfc6LyeUSetm344DFdGaW0tH0j9jBTe1eCP6TTYsS8nxCQBzQSzGSbD8/HxXUFAQdBgiqaVsNzxwih90/tFsv9d2faxZAI+d5xfAXTsDslvGJcya7DuovNdA8Y6NcN8JkNMGRr3lN0KK0z3HTV7ItRUvMib0HH8o/T5PZAzlomPyeGtp0V6xAHuNZ4Djkew7OTnrU7JunAutqkhOOzbChEPhzN/BgBvjEnNtmdlc51x+Va9pMFukscjKhosfhgcHwuQf+CquGfv/dR6TomW+0F/TXPj+c0lLEhCpkltNS6B5Wz8t94mL/O55Z4+P2z0PWjuL7370HFPKBzC9xcXcPrhX1BbJt8msGaUnjSfrX+fDzHFM6X77/omu+SL/prz+cYk33tSiEGls5v0Dpo6G/GvhnLtq3ye+bqH/Iga45lVoe1j8Y6yv6f8DHz3gk+Hh36v/9dYvhYfO8J/12hl1a6m8PQHe+j0jSsbyTsVRew7nhDKZ0utNeq6YBOO+Ssi04lhEa1FojEKksel/FQz4KRQ8DG/eVrtV25+/C48MgYwsGDE1NZME+CnB7XvBSz+ErfVcgV28GZ653CeH4U/VuTtravOL+bziIG7NepQmlH57+dJy3Io3oXN+YEmiJkoUIo3RGbfCMdfA7Lvh5Z9Ayc7o51eUw7t3+S1JW3aE616DDr2SEmqdhHLgksf8jnjPXln3KrPlpeH1IV/BZf+oenwhRn9643N+XXY1h2asY1TmtD3HO7GBXm4ldD+zztdONCUKkcbIDM69G079JSx4Eh44GZbN2H//iooK+Ox1mHiqb30ccR5c/zq06hxM3LXRoRdc8ACsmQf/vDHmvTm+rQk1jWl/vNSvDznvz/VeH7JmSzHvVhzFtPLjGZ01hd72OQBXZIVnQ/W5tF7XTyQNZos0VmZw+s3Q5USYeiM8fRm0OsR/Iea0hh1F8NWH4T0SOvuB8N4XpuQ8/2odcS4M/JVfbJjTBgbfHjX+b1dgl/GLrOc5t3wWf6u4mM6czrB6htIpN4fCLcXcWno1R2ev5JEmE5hUdjbXZb7K6ryz6dy6Sz3vkDgazBYR38WyeAosmeL3ty7dAdmtoFNf6H2BX5+QlR10lHXjHMy8BT6413e3DbnDlzYJqzzdNsOMClfOuKynGZX1Ck+Xnc64suvJy23Ge2MH1iuMymVADrfVPBi6k24ZX7Ou6aEcNPo1aBHsimxNjxWR6DJDcNQl/qehMYNBf/CJbvZdsHEFDL0XWnfZr4ZTG7eZO0L3c2rmJzxadha3lY0ALOaaUNFEptJOmLmMlVs6MzLnXv73pGaceeLxeyWuVBRIi8LMLgFuxe+LfZxzrso//81sMHAPkAk85JyLaVK0WhQiUqX5T/qV5BVl0Hc4N3zSlUXbW9CBLZyROZcrMt8kkwpuKxvBU+UDAd9NlWlGhXP7L/BrQFKxRbEIuBB4oLoTzCwTuBc4E1gNzDGzqc65JckJUUQanH5XQLdT4O3xsOAp/lq+G8I9auXOmFZxAn8uu4jP9ynNUR7+gzpe+1ikm0AShXPuUwCLPih2HLDCObcqfO4zwFBAiUJE6i73YN/1NOiP/OzOSbgdRWyhBfMrDmcrLYBvWxAZZnuSRERt97FoCFJ5emwe8FWl56vDx6pkZqPMrMDMCoqKihIenIikuaatOHXIZczMPJV/VRy9J0nkhDK589K+fD7+HCqq6ZqPx5hFOklYojCzN8xsURU/QxNxP+fcROdcvnMuv3371KvnLiKpp6aqsvXdx6KhSFjXk3PujHpeohCovE1W5/AxEZG4iVZkcMygnvtUgq3dPhYNRSrPyZoDdDezbvgEMRz4frAhiUhjUnlKa5VlzRuJQBKFmV0A/BVoD7xiZgucc4PMrBN+GuwQ51yZmY0GZuKnxz7snFscRLwi0nhFLWveSAQ16+kl4KUqjq8BhlR6Ph2YnsTQRERkH6k860lERFKAEoWIiESlRCEiIlEpUYiISFRKFCIiEpUShYiIRKVEISIiUSlRiIhIVEoUIiISlRKFiIhEpUQhIiJRKVGIiEhUShQiIhKVEoWIiESlRCEiIlEpUYiISFRKFCIiEpUShYiIRBVIojCzS8xssZlVmFl+lPO+MLOFZrbAzAqSGaOIiHiB7JkNLAIuBB6I4dzTnXMbEhyPiIhUI5BE4Zz7FMDMgri9iIjUQqqPUTjgNTOba2ajgg5GRKQxSliLwszeAA6q4qVbnHMvx3iZk5xzhWbWAXjdzJY6596p5n6jgFEAhxxySJ1iFhGR/SUsUTjnzojDNQrDj+vN7CXgOKDKROGcmwhMBMjPz3f1vbeIiHgp2/VkZs3NrGXkd+As/CC4iIgkUVDTYy8ws9XACcArZjYzfLyTmU0Pn3YgMNvMPgY+Al5xzs0IIl4RkcYsqFlPLwEvVXF8DTAk/PsqoG+SQxMRkX2kbNeTiIikBiUKERGJSolCRESiUqIQEZGolChERCQqJQoREYlKiUJERKJSohARkaiUKEREJColChERiUqJQkREolKiEBGRqJQoREQkKiUKERGJSolCRESiUqIQEZGolChERCQqJQoREYlKiUJERKIKJFGY2QQzW2pmn5jZS2aWW815g81smZmtMLOxSQ5TREQIrkXxOvAd59xRwHJg3L4nmFkmcC9wNnAkcLmZHZnUKEVEJJhE4Zx7zTlXFn76AdC5itOOA1Y451Y550qAZ4ChyYpRRES8rKADAK4Fnq3ieB7wVaXnq4Hjq7uImY0CRoWfbjezZXWMpx2woY7vTQXpHj+k/2dI9/gh/T9DuscPyf8MXap7IWGJwszeAA6q4qVbnHMvh8+5BSgDnqzv/ZxzE4GJ9b2OmRU45/Lre52gpHv8kP6fId3jh/T/DOkeP6TWZ0hYonDOnRHtdTO7GjgX+J5zzlVxSiFwcKXnncPHREQkiYKa9TQY+B/gfOfczmpOmwN0N7NuZtYEGA5MTVaMIiLiBTXr6W9AS+B1M1tgZvcDmFknM5sOEB7sHg3MBD4FnnPOLU5CbPXuvgpYuscP6f8Z0j1+SP/PkO7xQwp9Bqu610dERMTTymwREYlKiUJERKJSoghL93IhZvawma03s0VBx1IXZnawmb1lZkvMbLGZ3RR0TLVlZk3N7CMz+zj8GW4LOqa6MLNMM5tvZtOCjqUuzOwLM1sYHv8sCDqe2jKzXDN7IVzm6FMzOyHwmDRGsadcyHLgTPzCvjnA5c65JYEGVgtmdgqwHXjcOfedoOOpLTPrCHR0zs0zs5bAXGBYmv07MKC5c267mYWA2cBNzrkPAg6tVszs50A+cIBz7tyg46ktM/sCyHfOpeWCOzN7DHjXOfdQeMZnM+fcliBjUovCS/tyIc65d4BNQcdRV865tc65eeHfv8HPdMsLNqracd728NNQ+Cet/hIzs87AOcBDQcfSGJlZK+AUYBKAc64k6CQBShQRVZULSasvqYbEzLoC/YAPAw6l1sLdNguA9cDrzrl0+wx/xq9xqgg4jvpwwGtmNjdc2ieddAOKgEfC3X8PmVnzoINSopCUYmYtgBeBnzrntgUdT20558qdc0fjKwkcZ2Zp0w1oZucC651zc4OOpZ5Ocs71x1ee/km4WzZdZAH9gfucc/2AHUDgY6ZKFJ7KhaSAcL/+i8CTzrnJQcdTH+HugreAwQGHUhsDgPPDffzPAAPN7IlgQ6o951xh+HE98BK+azldrAZWV2qJvoBPHIFSovBULiRg4YHgScCnzrm7go6nLsysfWQTLjPLwU+OWBpoULXgnBvnnOvsnOuK/39glnPuyoDDqhUzax6eDEG4y+YsIG1mAjrn1gFfmVnP8KHvAYFP6EiFMuOBc86VmVmkXEgm8HCSyoXEjZk9DZwGtDOz1cBvnHOTgo2qVgYAVwELw338ADc756YHF1KtdQQeC8+iy8CXnUnLKaZp7EDgJf93B1nAU865GcGGVGs3AE+G/2hdBVwTcDyaHisiItGp60lERKJSohARkaiUKEREJColChERiUqJQkREolKiEBGRqJQoREQkKiUKkQQzs2PN7JPwfhXNw3tVpE0NKBEtuBNJAjP7PdAUyMHX8rk94JBEYqZEIZIE4XIMc4BdwInOufKAQxKJmbqeRJKjLdACaIlvWYikDbUoRJLAzKbiS3d3w2/5OjrgkERipuqxIglmZiOAUufcU+HKsu+b2UDn3KygYxOJhVoUIiISlcYoREQkKiUKERGJSolCRESiUqIQEZGolChERCQqJQoREYlKiUJERKL6f7Z/eBIUUOA4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the data together with the fit.\n",
    "\n",
    "x_pts = np.arange(0,2*np.pi,0.01)\n",
    "y_fit = yhat(x_pts, beta_list)\n",
    "\n",
    "plt.plot(x_data, y_data, 'o')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.plot(x_pts, y_fit)\n",
    "plt.ylim(-2,2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we decide how good a fit like this is? To do this, we can use the mean squared error, which is the average squared distance between the data points and the curve at each x value in the dataset:\n",
    "$$\n",
    "\\mathrm{MSE} = \\frac{1}{n} \\sum_{i=1}^n (y_i - \\hat{y}(x_i))^2,\n",
    "$$\n",
    "where $(x_i, y_i)$ denotes one of the data points, and $\\hat{y}(x)$ is the curve that fits the data. If the MSE is close to zero, then the curve is doing a good job of fitting the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1** (3 points): Create a function `mse(x, y, p)` that fits a polynomial of order `p` to data `x` and `y` and returns the mean squared error from the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Solution ##\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2** (3 points): For the data that we generated above, make a plot of the MSE for different values of $p$. Make sure to label the axes appropriately. You should find that the MSE decreases as $p$ increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Solution ##\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In class, we found that, although higher-order polynomials may fit the data very well, they typically won't do a great job at describing new data that was not used for the fit. The idea of **cross validation** is to split the data into a **training set** and a **testing set**. The training set is used for fitting, and the testing set is used to quantify the goodness of the fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3** (5 pts): Create a new version of the function from the previous exercise. Call this function `mse_cv(x, y, p)`, which takes the same arguments as before. In this version, use the first half of the `x` and `y` data to perform the fit with a polynomial of order `p`. Compute the MSE using this first half of the data (the \"training set\"), then separately compute the MSE using the \"held out\" second half of the data (the \"testing set\"). Have the function return both MSE values. Finally, make a plot showing both the training MSE and the testing MSE for different values of $p$, using a legend to show which curve is which, and report the order of the polynomial that gives the optimal fit using cross-validation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Solution ##\n",
    "\n"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}