{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Homework 10 - Spline Regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this homework, you will implement a general spline regression.\n", "Recall that a spline of degree $d$ is a piecewise polynomial of degree $d$ with continuity of derivatives of orders $0, 1, \\ldots, d-1$.\n", "\n", "In the case of a cubic spline ($d = 3$) with $K$ knots, this regression function can be expressed by\n", "\n", "$$ f(x_i) = \\beta_0 + \\beta_1 \\, b_1(x_i) + \\ldots + \\beta_{K+3} \\, b_1(x_i) $$\n", "\n", "using appropriate basis functions.\n", "\n", "From Slide 352, you know that we can start off with monomials $x, x^2, x^3$ and then add for each knot $\\xi$ one **truncated monomial**\n", "\n", "$$ h(x,\\xi) = (x-\\xi)_+^3 = \\begin{cases} (x - \\xi)^3 & \\text{if } x > \\xi, \\\\ 0 & \\text{otherwise.} \\end{cases} $$\n", "\n", "\n", "A one-dimensional example is provided below." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def y_func(x):\n", " return np.sin(8*x) / np.exp((5*x))\n", "n = 100\n", "eps = 0.1\n", "np.random.seed(1)\n", "\n", "x = np.random.rand(n)\n", "y = y_func(x) + eps * np.random.randn(n)\n", "plt.scatter(x, y, label='Samples')\n", "xlin = np.linspace(0,1,100)\n", "plt.plot(xlin, y_func(xlin),'r--', label='Population line')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Task**: Write a function that implements spline regression of degree $d$ and try to estimate the function `y_func` using **cubic** spline regression with $K$ equidistant knots (use $K=4$ in your tests)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def generateSplineRegressionMatrix(x, xi, d = 3):\n", "\n", " # Get number of samples\n", " n = len(x)\n", " \n", " # Get number of knots\n", " K = len(xi)\n", " \n", " # Initialize matrix with predictor variables 1\n", " X = np.ones((n,1))\n", "\n", " # Append columns for the monomials x^1, ..., x^d\n", " for i in range(d):\n", " X = np.hstack((X,np.power(x,i+1).reshape(n,1)))\n", " \n", " # Append one column for each knot using the \n", " # truncated monomial (x - xi_k)_+^d\n", " for k in range(K):\n", " h = np.power(np.maximum(x-xi[k],0), d)\n", " X = np.hstack((X,h.reshape(n,1)))\n", " return X" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.2 0.4 0.6 0.8]\n" ] } ], "source": [ "K = 4\n", "xi = np.linspace(0,1,K+2)\n", "xi = xi[1:-1]\n", "print(xi)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd8FNX6+PHPyaaSQoCEFkroBEhMIDRDU2mCFBE7ihcU9Ip8bfzECnYUFSwoolfBhl4BEUQvKgERpBOIQECQmlADJCSkJ+f3xyYxCbupu9n2vF8vXmRnZ2fOJLPPnHlOGaW1RgghhGtxs3UBhBBC1D4J/kII4YIk+AshhAuS4C+EEC5Igr8QQrggCf5CCOGCJPgLIYQLkuAvhBAuSIK/EEK4IHdbF8CcoKAgHRoaautiCCGEQ9mxY0ey1jq4ovXsNviHhoayfft2WxdDCCEcilLqWGXWk7SPEEK4IAn+QgjhgiT4CyGEC7LbnL8Qwvpyc3NJTEwkKyvL1kURVeTt7U2zZs3w8PCo1ucl+AvhwhITE/H39yc0NBSllK2LIypJa8358+dJTEykVatW1dqGpH2EcGFZWVk0aNBAAr+DUUrRoEGDGt2xSfAXwsVJ4HdMNf27SfAXQggXJMHfRpbHJREzK5ZW01cRMyuW5XFJti6SEDZhMBiIjIykS5cu3HzzzWRkZFh0+wsXLmTKlCnlrrNu3Tr++OOP4tfz58/ns88+q/G+jx49SpcuXQDYvn07U6dOrfE2LUWCvw0sj0viyWV/kpSSiQaSUjJ5ctmfcgEQLsnHx4ddu3axZ88ePD09mT9/fq2XoWzwv//++7n77rstuo/o6Gjeeecdi26zJiT428Ds1QfIzM0vtSwzN5/Zqw/YqERC2Ie+ffty6NAhAN566y26dOlCly5dmDt3LmCsSXfs2JHx48cTERHB2LFji+8UQkNDSU5OBoy17AEDBlyx/ZUrV9KzZ0+ioqIYOHAgZ86c4ejRo8yfP585c+YQGRnJ77//zsyZM3njjTcA2LVrF7169SIiIoIbb7yRixcvAjBgwACeeOIJevToQfv27fn999/LPbZ169Zxww03ADBz5kwmTJjAgAEDaN26damLwhdffEGPHj2IjIxk8uTJ5Ofnm9tkjUjwt4GTKZlVWi5ErRkw4Mp/779vfC8jw/T7Cxca309OvvK9KsjLy+Onn34iPDycHTt28Omnn7JlyxY2b97MRx99RFxcHAAHDhxg0qRJxMfHExAQwPtF5auEPn36sHnzZuLi4rjtttt4/fXXCQ0N5f777+eRRx5h165d9O3bt9Rn7r77bl577TXi4+MJDw/n+eefL1XmrVu3Mnfu3FLLK2P//v2sXr2arVu38vzzz5Obm0tCQgLffPMNGzduZNeuXRgMBr788ssqbbeyJPjbQNNAnyotF8KZZWZmEhkZSXR0NC1atGDixIls2LCBG2+8EV9fX/z8/BgzZkxxzbp58+bExMQAMG7cODZs2FDpfSUmJjJkyBDCw8OZPXs2e/fuLXf91NRUUlJS6N+/PwDjx49n/fr1xe+PGTMGgG7dunH06NGqHDbDhw/Hy8uLoKAgGjZsyJkzZ1izZg07duyge/fuREZGsmbNGg4fPlyl7VaWDPKygWlDOvDksj9LpX58PAxMG9LBhqUSAli3zvx7deqU/35QUPnvm1GU8y9Ja212/bJdHIteu7u7U1BQAGC2//tDDz3Eo48+ysiRI1m3bh0zZ86scnlL8vLyAoyN1nl5edX6bMnPa60ZP348r776ao3KVRlS87eB0VEhvDomnJBAHxQQEujDq2PCGR0VYuuiCWEX+vXrx/Lly8nIyODy5ct89913xemY48ePs2nTJgAWL15Mnz59AGPOf8eOHQAsXbrU5HZTU1MJCTF+zxYtWlS83N/fn7S0tCvWr1u3LvXq1Su+6/j888+L7wKs4brrrmPJkiWcPXsWgAsXLnDsWKVmaK4yqfnbyOioEAn2QpjRtWtX7rnnHnr06AHAvffeS1RUFEePHiUsLIxFixYxefJk2rVrxwMPPADAjBkzmDhxIq+88go9e/Y0ud2ZM2dy8803ExISQq9evThy5AgAI0aMYOzYsXz//fe8++67pT6zaNEi7r//fjIyMmjdujWffvqp1Y67U6dOvPTSSwwePJiCggI8PDyYN28eLVu2tPi+VHm3V7YUHR2t5WEuQlhXQkICYWFhti5GpR09epQbbriBPXv22LoodsHU308ptUNrHV3RZyXtI4QQLsgiwV8pNVQpdUApdUgpNd3MOrcopfYppfYqpb6yxH6FEK4lNDRUav0WUuOcv1LKAMwDBgGJwDal1Aqt9b4S67QDngRitNYXlVINa7pfIYQQ1WeJmn8P4JDW+rDWOgf4GhhVZp37gHla64sAWuuzFtivEEKIarJE8A8BTpR4nVi4rKT2QHul1Eal1Gal1FAL7FcIIUQ1WaKrp6lJpct2IXIH2gEDgGbA70qpLlrrlFIbUmoSMAmgRYsWFiiaEEIIUyxR808Empd43Qw4aWKd77XWuVrrI8ABjBeDUrTWC7TW0Vrr6ODgYAsUTQhh70pO6TxixAhSUlIq/lAtGjZsmEXKVHKyuOeee45ff/21xtusCUsE/21AO6VUK6WUJ3AbsKLMOsuBawCUUkEY00DWmbBCCOFQSk7pXL9+febNm2eR7VZ1ugVzfvzxRwIDAy2yrSIvvPACAwcOtOg2q6rGwV9rnQdMAVYDCcB/tdZ7lVIvKKVGFq62GjivlNoHrAWmaa3P13TfQgjn0rt3b5KS/nmuxezZs+nevTsRERHMmDGjePmLL75Ix44dGTRoELfffntxjXrAgAE89dRT9O/fn7fffptz585x00030b17d7p3787GjRsB+O2334iMjCQyMpKoqCjS0tI4deoU/fr1K74LKZrSoeRU0eammQ4LC+O+++6jc+fODB48mMzM8mfoveeee1iyZEnx9mfMmEHXrl0JDw9n//79AFy+fJkJEybQvXt3oqKi+P777y3xKy5mkekdtNY/Aj+WWfZciZ818GjhPyGEHXr4fw+z6/SuilesgsjGkcwdOrdS6+bn57NmzRomTpwIwM8//8zBgwfZunUrWmtGjhzJ+vXrqVOnDkuXLiUuLo68vDy6du1Kt27direTkpLCb7/9BsAdd9zBI488Qp8+fTh+/DhDhgwhISGBN954g3nz5hETE0N6ejre3t4sWLCAIUOG8PTTT5Ofn3/FE8VKTjOttaZnz57079+fevXqcfDgQRYvXsxHH33ELbfcwtKlSxk3blylf09BQUHs3LmT999/nzfeeIOPP/6Yl19+mWuvvZZPPvmElJQUevTowcCBA/H19a30dssjc/sIIWyqaErno0eP0q1bNwYNGgQYg//PP/9MVFQUAOnp6Rw8eJC0tDRGjRqFj49xCvQRI0aU2t6tt95a/POvv/7Kvn3FQ464dOkSaWlpxMTE8Oijj3LnnXcyZswYmjVrRvfu3ZkwYQK5ubmMHj2ayMjIUtstOc00UDzN9MiRI2nVqlXx+tWZ3rnk1NDLli0rPv4VK1YU39VkZWVx/Phxi03HIcFfCAFQ6Rq6pRXl/FNTU7nhhhuYN28eU6dORWvNk08+yeTJk0utP2fOnHK3V7JmXFBQwKZNm4ovFEWmT5/O8OHD+fHHH+nVqxe//vor/fr1Y/369axatYq77rqLadOmlXqUY3nzoJWdnrmitI+5z5ecGlprzdKlS+nQwTpTvcvcPkIIu1C3bl3eeecd3njjDXJzcxkyZAiffPIJ6enpACQlJXH27Fn69OnDypUrycrKIj09nVWrVpnd5uDBg3nvvfeKXxc9N+Dvv/8mPDycJ554gujoaPbv38+xY8do2LAh9913HxMnTmTnzp2ltlXeNNPWMGTIEN59993ii07Rk8wsRWr+VrA8LonZqw9wMiWTpoE+TBvSQaZvFqISoqKiuOqqq/j666+56667SEhIoHfv3gD4+fnxxRdf0L17d0aOHMlVV11Fy5YtiY6Opm7duia398477/Dggw8SERFBXl4e/fr1Y/78+cydO5e1a9diMBjo1KkT119/PV9//TWzZ8/Gw8MDPz8/Pvvss1LbKm+aaWt49tlnefjhh4mIiEBrTWhoKD/88IPFti9TOlvY8rgkk0/pkoe1CHvkaFM6F0lPT8fPz4+MjAz69evHggUL6Nq1q62LVetkSmc7Mnv1gVKBHyAzN5/Zqw/YqERCOJ9JkyYRGRlJ165duemmm1wy8NeUpH0s7GSK6YYec8uFEFX31VcyK3xNSc3fwpoG+lRpuRBC2IIEfwubNqQDPh6GUst8PAxMG2Kd7lpCCFEdkvaxsKJGXentI4SwZ1LzF0IIFyTB38KKunompWSigaSUTJ5c9ifL45Iq/KwQrujll1+mc+fOREREEBkZyZYtW6y2rwEDBuCIXcitQdI+FlZeV8/aTP3IQDPhCDZt2sQPP/zAzp078fLyIjk5mZycHFsXyyVIzd/C7KGrp9x9CGtZHpdEzKxYWk1fRcys2BqfU6dOnSIoKKh4bpugoCCaNm3KCy+8QPfu3enSpQuTJk0qnuJgwIABPPLII/Tr14+wsDC2bdvGmDFjaNeuHc888wxgnGK5Y8eOjB8/noiICMaOHXvFDJ1gnDitd+/edO3alZtvvrl4Gonp06fTqVMnIiIiePzxx2t0fPZMgr+FFH0pzI2Xrs2unjLQTFiDNSoVgwcP5sSJE7Rv355///vfxVMxT5kyhW3btrFnzx4yMzNLTWvg6enJ+vXruf/++xk1ahTz5s1jz549LFy4kPPnjY8JOXDgAJMmTSI+Pp6AgADef//9UvtNTk7mpZde4tdff2Xnzp1ER0fz1ltvceHCBb777jv27t1LfHx88QXFGUnwt4CSXwpTPNwUGTl5FqstVcQe7j6E87FGpcLPz48dO3awYMECgoODufXWW1m4cCFr166lZ8+ehIeHExsby969e4s/M3Kk8RlR4eHhdO7cmSZNmuDl5UXr1q05ceIEAM2bNycmJgaAcePGsWHDhlL73bx5M/v27SMmJobIyEgWLVrEsWPHCAgIwNvbm3vvvZdly5ZRp06dah+bvZOcvwWY+lIUCfTx4HJOHhczcoF/aktApXLw1cndNw30MXkhkoFmoiasVakwGAwMGDCAAQMGEB4ezocffkh8fDzbt2+nefPmzJw5k6ysrOL1i1JEbm5upaZSdnNzK54OWSlVah9lX2utGTRoEIsXL76iPFu3bmXNmjV8/fXXvPfee8TGxtbo+OyV1PwtwNzJrwBfL3dy80sngypbW6rubbYMNBPWYI3R6wcOHODgwYPFr3ft2lU8f31QUBDp6enFjzusiuPHj7Np0yYAFi9eTJ8+fUq936tXLzZu3MihQ4cAyMjI4K+//iI9PZ3U1FSGDRvG3Llzi6eAdkZS87eA8mraNaktVbfnkAw0E9YwbUgHkzPW1qRSkZ6ezkMPPURKSgru7u60bduWBQsWEBgYSHh4OKGhoXTv3r3K2w0LC2PRokVMnjyZdu3a8cADD5R6Pzg4mIULF3L77beTnZ0NwEsvvYS/vz+jRo0iKysLrXWFD45xZDKlswWUN43z7NUHTF4YQgJ92Dj92nK322r6KpMNyAo4Mmt4DUstRNWndHaELsRHjx7lhhtuYM+ePbYuitXVZEpnqflbQHk17Us5F3ly1WKy87Mw6EAMOhB/9xaVqi1J7l7Ym9FRIXYX7EX1SPC3kJJfitz8XN7e8ja9Pl7C1qStaIOGEin4NI+6bDh3L10u3E/b+m3NbtMat9lCOLvQ0FCXqPXXlAR/C9t3bh93fXcXO0/tpGdIT2b0n8GQtkNo4NOAs5fPcjLtJEsSlvD2lrd5c9Ob9AkZRfa58ZxL9SSwjgdaQ2pmbvHdQ1HqyJ5vs4Vj01pf0RtG2L+apuwl529B7255l2m/TMPfy58FNyzgxrAbza57Ku0UU1fOYslf72MggAY5j+BTEFVqnZKPf3SEXKtwPEeOHMHf358GDRrIBcCBaK05f/48aWlptGrVqtR7kvOvZe9ueZep/5vKDe1v4OMRH9PIr1G56zfxb8LJ46Nokt2GZM83OOv1LAG5YwnMG4/C+CUs2SW0ZPqnqmMFhDCnWbNmJCYmcu7cOVsXRVSRt7c3zZo1q/bnJfhbwHcJ3/F///s/RnUYxdJblmJwM1T8IYxB3JM2NM6ey0WPj7jksYQCdYn6uQ+iChsJTqZk2s1kccL5eHh4XFFzFK5Bgn8NbTqxiTuW3UHPZj356qavKh34AQxKka81bnhRP/dBDDqQVI+vKSCToNzHULjXeKyAEEKYIiN8a+Dc5XOM/HokzQKaseK2FdTxKGcekJwc+N//4NIl4+vYWD79+hle+d+7XL9/A3WzLhOYN47A3H+R4f475zxn4e1h7PEjzwUWQlia1PxrYNov00jNSuW3e34j2DfY9EoXL8Izz8CXX0JqKnzxBdx5J+TkEJSXQWTC79yxezX5yo31raL4f9c/zN91PbjouYA2rZYwKnIYgHT5FEJYlAT/avrt6G8s2r2Ip/o8RafgTqZXWrUKJk2CM2eMAf/mm2HgQON7Q4fy14o13PptHO2PJzDg8A76Ht1JdkBdPhn7HBvO+vLm5jnM3RzBI70fAWS6BiGE5UhXz2rIyc8hcn4kmXmZ7P33XtPpnoICiImB9HRYuBC6dTO5rVJdOOt6M21oR0aH1qGgbx9uGV+HZRk7WHLLEsaEjbHuQQkhnIJ09bSiN/94k4TkBFbdserKwJ+fD3l54OUFy5ZB/frGn80wOVz+zBncgoL5/KnfSJoWzF3f3UXHoI7m7zCEEKKKpMG3ik6lneLF9S8yJmwMw9oNK/2m1vDAAzB8OOTmQpMm5QZ+sxo1gl9/xeeJp1n63jl803O4efEYLudctsxBCCFcngT/Knp7y9tk52fz2sDXrnzzscfgo4+gZ0/w8KjZjtzd4aWXaLr4B776zo2ECwd4YNUDNR7SLYQQIGmfKknNSuWD7R9wU9hNV07ItnQpzJkDDz0EL71Upe2WO3XD8OEMbPAbMy5+x8ytr9O3RV/u63afhY5ICOGqLBL8lVJDgbcxzl35sdZ6lpn1xgLfAt211vbZmlvIVED+K+NLLmVf4omYJ0qvfOoUTJ4M0dHw5ptQhTlSyj4LwOTUDb168UxBdzae3cHUFQ/Q16sdHbsMsMRhCiFcVI3TPkopAzAPuB7oBNyulLqiZVIp5Q9MBbbUdJ/WZurxiU8s28Gs39/iulbX0a1pmZ47yckQEgKff17ldE9lH4ptcDOw6KoZ+GYXMO6j68k9L3OxCCGqzxI5/x7AIa31Ya11DvA1MMrEei8CrwNZJt6zK6YCcnLBGi5mn7my1g8QHg67dkHHjlXeV1WmbmgS2ZcPo55lR/0sXnqsOxQ+fk4IIarKEsE/BDhR4nVi4bJiSqkooLnW+gcL7M/qygZeTQGX3JfhWdCGga0H/vPGqVPGRt6MjCqlekqq6tQNN936PHf79+XllsfYPHm4sYeREEJUkSWCv6moVxyRlFJuwBzgsQo3pNQkpdR2pdR2W04xWzbwZrnFk+eWSEvPW0rPef7ii/DOO8aLQDVNG9IBH4/Sk8FVNHXDO/9eSYghkLvrrSXjUEK19y2EcF2WCP6JQPMSr5sBJ0u89ge6AOuUUkeBXsAKpdQVI9C01gu01tFa6+jgYDNz5dSCazoGl7qiXTbE4qZ9eXHIhH8WHjpk7NY5aRK0aVPtfY2OCuHVMeGEFF5wDEoV5/yXxyWZ/Exd77p8Om4JBwMLmHHk02rvWwjhuiwR/LcB7ZRSrZRSnsBtwIqiN7XWqVrrIK11qNY6FNgMjLTX3j7L45JYuiOp+NalgEwyDH8Q0WAot0aX6N753HPg6WmctK2GRkeFFN8B5BemcYp6/Zi7AFzb+jomd5vMW5vfYuurU+DIkRqXQwjhOmoc/LXWecAUYDWQAPxXa71XKfWCUmpkTbdfW5bHJREzK5aHv9lVqrE3w7AJrbLIS+v7z8q7dsHixfDww8ZRvBZQ2V4/Jb0+6HWa1mnEhJMfkH3zjZBl923pQgg7YZERvlrrH7XW7bXWbbTWLxcue05rvcLEugPsrdZfsmtnWZcNsbgXNCLtUonUjq8v3HEHTJtmsTJU54EtsfvS8M9+iL1BBbwcsJu/J0yxWHmEEM5NpnfAdK0bII9kstx245t/LSGBJSZwa9fOOD9/YGDxHUOr6auImRVrNk1Tkar2+im6YGWkReCbdy2v9FVcjv0Pm979rFr7F0K4Fgn+mK9dX3ZfB0rTQA38p/fNwoWwZw9gejBYeXn68lS110/JC1a93HsBf8aN8aT1M49I+kcIUSEJ/piuXWs0lw2x+KsuvDlmiHGqhbNnjdM4zJ8PVC9Pb07JXj8KCAn04dUx4WYf2FLygmUggHp5k/mzUQ6D7+oL3t5V3r8QwrXIxG4Ya91lH5Po5nGcXLfjvDts/j8BeMEC47N4pxhz65Z+sLrJuf3NaBroU6qNok5+P7zzY9kTvJrjqcdpkekJjRtXqxxCCOcnNX9M17p7dz6CQnFj2I3GlXJz4YMPYNCg4mkcbPlg9bJpIoUiRD+Ep8GNBz4Yjm7bxjgWQQghTJDgX2h0VAgbp1/LkVnD2Tj9Wg6l/UbPZj1p6NvQuMKyZXDyJEydWvyZ6ozOtWR5y16w3hgzkFcHvsyP2Xv4bycN//qX8cliQghRhqR9TDiVdoptJ7fx8rUv/7Pw9GmIiIDrry9eVJSisdWD1U2lifILHuLLP79k6oi/GPzKBuq9/TY8+mitlEcI4TjkAe4mfLzzY+5beR+7799NRKOIf94oKAA3+79ZijsVR/ePujPhbAgLPj4D8fHQvr2tiyWEqAWVfYC7/UcyG/jhrx9oUbcF4Q3DjQvOnDHOnukAgR8gqkkUj/Z+lI+Cj7M+PAD27rV1kYQQdsYxolktysrL4pfDvzCi/QjjDJ4FBcZn8k6eXKnPW2rQV03N6D+DVoGtmHR3PbJHDKv4A0IIlyLBv4zYI7Fk5GYwov0I44I//oBjx6Bfvwo/a8lBXzXl6+nLB8M/4MCFv3j191eMI5KTbHMhEkLYHwn+Zaw8sBJfD1/6h/Y3LvjqK/DxgdGjK/ysJQd9WcKQtkO4vcvtvLphFglPTIAHH5SHvwghAAn+pWit+eHgDwxuMxhvd29j3/7//hdGjQI/vwo/b+lBX5YwZ8gcfD19mTw5hIIV38N339msLEII+yHBv4T4M/EkXkr8J+Xz889w/jzceWelPm/LQV/m2hoa+TVi9qDZ/F5whE9GNDOOU0hLs8q+hBCOQ4J/CbFHYgEY1GaQccE118DXX8PgwZX6vK0GfVXU1jAhagL9WvZjWo9UzqQmwbPPWm1fQgjHIMG/hHXH1tG2fluaBTQzLqhTB2691fjErkqo6uRsllJRW4NSig9v+JAMnc0jj3eB/v2tti8hhGOQEb6F8gvyWX9sPWPDxhoX/PabsafP1KnGh7dUUlUmZ7MUc20KSSmZLI9LYnRUCB2DOvJUn6eY+dtM7g73YaiF92XLdg0hRNVJzb9Q/Jl4UrJSGBA6wLhg4UJ4/XXw8rJlsSqlvDaFkimZ6X2m0zGoIw/88ACXX5phnKXUQvuqjXYNIYTlSPAvtPboWgBj8C8ogFWrjPP4uNv/zZGptoYiJVMyXu5efHjDhxxNPcrzxxYZH0N5+nSN91Vbk9kJISxHgn+hdUfX0a5+O0ICQmDrVjh3DkaMsHWxKqWorcGckimZfi37cW/UvbzVPJFdARnw+OPV2ldtt2sIISzL/qu1NbA8LqlSM24W5ftv6XyLccHKlWAwwNDqZsZr3+ioEGavPmDyIfRlUzKvD3qdlX+t5L6Jbmx+4UsMEycaezZVYV8S7IVwbE5b869Kl8TdZ3aTmp36T74/JQWuuw7q1avVMtdUZVMy9Xzq8fbQt9muTvHe9fXhoYeMqS4hhMtw2uBflS6Ja4+UyPcDzJsHP/1k7SJaXFVSMrd0voVh7YbxdO9Mjn8wy2FmLBVCWIbTpn2q0iVx3bF1tG/Qnqb+TY1PvjIYHDYYVjYlo5Ti/WHv0+n9Tvz7xHxW6uGoggLjsVP5lJkQwjE5ZoSrhMp2SSzK9w9oOcC4YNQoGDfOyqWzDy0DW/LSNS+x6uAqvp02DG67DZBRvEK4AqcN/pXNf+86vYtL2ZeMKZ/0dPjlF2jYsBZLalsP9XyIbk26MdV/Axd/WAKrV8soXiFcgNMG/8rmvzclbgKgT4s+sG4d5OTA8OG1X2AbcXdz5+ORH5OsMpk2NgAeeojk5Esm15VRvEI4D6fN+UPl8t9bkrbQxK+JcT6f2DnGEb0xMbVUQvsQ2TiSx3o/xuv6dcZtuMSjf67k1agxV6wno3iFcB5OW/OvrC2JW+gR0sP4yMbYWGPg9/a2dbFq3YwBM2hdrzWTbvfj1p3LCSSv1PsyilcI5+LSwf9C5gUOXjhIz5CexidcTZwIDzxg62LZRB2POnx4w4cc9Epn9kd3MPPWaBnFK4QTc+q0T0W2JW0DoGeznqCUcbCTCxvYeiD3RN7D6/Hz2d5zIhvv6QSNG1fqs9I1VAjH4tI1/y1JW1AooptGw7ZtcOqUrYtkc28OfpMGPg2Y+N4g8vpcDVlZFX5GuoYK4XhcPviHBYcR4BUAd90F991n6yLZXH2f+rx7/bvs8ExmbqMj8NprFX5GuoYK4XhcLu1TlJ5ISskgyWcD/ZpfD0lJcOAATJpk6+LZhbGdxjKqwyieK/iBke+/wjPpbdhmqGc2nSMPeBHC8bhUzb9keiJXnSaPS+w92pDtny41rnDttbYtoJ1QSjFv2DyUhzf3Dcvj/mVz0FqbTefIA16EcDwuFfxLpidy3P4CwC23HaeX/wT160NEhC2LZ1dCAkJooiazPrSAHY0PEHz5ImA6nSMPeHFty+OSiJkVS6vpq4iZFSttPQ7CIsFfKTVUKXVAKXVIKTXdxPuPKqX2KaXilVJrlFItLbHfqiqZhsh2O4DSXnjolkQe3Mnapp1ZvlsafEvKS7sGn/wIpg8q4JRWKDsbAAAgAElEQVTfP1M+l03n2PoBLxJ8bEca+x1XjXP+SikDMA8YBCQC25RSK7TW+0qsFgdEa60zlFIPAK8Dt9Z031XVNNCn+GEn2W4H8Cxog8LAbbe/gldeLieX/QkgXRQLhQTWITd1Kqe8HuSix7uMPDiaTaFRJtM5tnrAS1HwKbqjKwo+RWUS1lVeY7/8/u2bJWr+PYBDWuvDWusc4GtgVMkVtNZrtdYZhS83A80ssN8qK0pPaHLJUYfxKjCmJRIDG/N3UHPpoVLGtCEdCHAPITD3bjLcdzA04TkiLh63q3SO9DSyLWnsd1yWCP4hwIkSrxMLl5kzEbDJk1KK0hN1A06CysWzoAM3x//CmD1riteRk/YfRb+vjn4345vbnoev1zy0bg6Pfr3TbtIrEnxsSxr7HZclgr8ysUybXFGpcUA0MNvM+5OUUtuVUtvPnTtngaJdaXRUCA8ONjZOeuq2TNy2nJH71he/LydtaaOjQvjjyUG8NuQj0r0MvNf1IGPjf7ab3K4EH9uSxn7HZYngnwg0L/G6GXCy7EpKqYHA08BIrXW2qQ1prRdoraO11tHBwcEWKJppu07voo57AE3zAmiffJwdIR0BOWnL89XGPPzzxrOiI7Q7+xH1M1LtIr0iwce2bN3YL6rPEoO8tgHtlFKtgCTgNuCOkisopaKAD4GhWuuzFthnjew+s5vokEheDczBDc3OkDBCZD6acp1MycSfUeRnr2Xa4GN0upgEderaPL1S9PeSeYVsx1aN/aJmahz8tdZ5SqkpwGrAAHyitd6rlHoB2K61XoExzeMHfKuUAjiutR5Z031XR4EuIP5MPBOiJnD1loPg5saXH00Ff39bFMdhFPWU8tP/j5NeU9nVcCnBOWGEBNaxddEk+AhRDRaZ3kFr/SPwY5llz5X4eaAl9mMJf1/4m8u5l7mq0VWQtNk4sEsCf4WmDelg7EKZ25x6eXdz0eM/XHfoeW58dLGtiyaEqAaXGuELxpQPwFWNr4KPPoItW2xcIsdQMrcbkDeSJtmtWNZuO5Er59i6aEKIanC5id12n96NQRnoHNzZuMDT07YFciAl0yt/X+hMxNsdmXjkfVbvug8VGWnj0gkhqsLlav67zuyiQ1AHfL5ZCsOGQVqarYvkkNrUb8Psa17hlzYw/8VRkJdX8YeEEHbD5YL/7tO7jfn+X36BnTvBz8/WRXJY9/d/jEG+ETze8Th/vfW0rYsjhKgClwr+FzIvcOLSCSIbR8Iff0Dv3sbHN4pqcVNufHrfKrw8fbjL92fyCqT2L4SjcKngv/t0YWOvd0s4dAiuvtrGJXJ8IXWbMX/sQrYm7+KV9S9DQUHFHxJC2JxrBf+inj7HCgcY9+5tw9I4j1s638K4jrfywrrn2fr2/7N1cYQQleBywb+RbyMa+zaCAQOga1dbF8lpvDvyA0JyvLjz2FukJ+y2dXGEEBVwqeC/6/QuY//+IUNg7VqoY/vRqc4i0KceX9z0JYfraqbOHQz5+RV/SAhhMy7Tzz83P5d95/YxqNVAyM4GLy9bF8np9O02hqc3jOJFt+/p9NQtfFfvQZlvRwg75TI1//3J+8nJzyHSq4Wxe+eXX9q6SE7puSnf0i21Li8YviPxwt/yaD8h7JTLBP+ixt6IU9o4IKlVKxuXyDm5Gzzw832Vy57enPGeg8bY/dMepn8WQvzDZYL/vnP7cHdzp/2+M+DmBjIdgdUcy21B/dwpZBsSCEh/o3i5rad/FkL8w2WCf0JyAm3rt8Vzxy4IC5PGXitqGuiDb35/eiW2ZU/wBlomLy1eLoSwDy7T4Lvv3D66BHeBHRuNvX1c3PK4JKs9AKVo+ueL/jNpn3wPexsspFlBL6YNGW6R7QvrsNY5Yc1zTVSfSwT/7Lxs/r7wNzd3vAke6wXh4bYukk0tj0viyWV/kplr7I5Z1CALWORLWfLpWt4HHyUx4HW8M59kePg9Nd52RSTQVI+1zglrn2v2xNHOPZdI+xy8cJB8nU+nRl1g2jQYOtTWRbKp2asPFH8Zi1iyQbbklyCr3RCeyx7BlroXeGrODRbZfnn7fXLZnySlZEovoyqy1jlh7XPNXjjiuecSwT/hXAIAYZfrwMkrni3vcsw1vFa1QXZ5XBIxs2JpNX0VMbNiWR6XZPJLsMj7fu682JE3Mn5hWcIyCxyBaa4SaCxteVwSSRY6Jyr7eWdr/HfEc88lgv++c/tQKDq88iEMHmzr4ticuYbXqjTImqvpzFyx94ovQUae5u9Gc+gR0oN7lt/DX6f31qT4ZrlKoLGkor+jOTVtpLfEueYIzJ1jSSmZpSpH9sQlgn9CcgKhgaHU2bYLunWzdXFsbtqQDvh4GEot8/EwMG1Ih0pvw1xNJyUz1+T6p1Pz+Xbsf/HMzOHGd64mLetS1QteAVcJNJZk6u9YpKrnhCmWONdqk6m72coo7xyz1zSQywT/MP9WcPq0BH9KP49XASGBPrw6JrxKjVNVrU03DfShRWBLvvG4gwMel7jrzRgKtGWnf3a0QGMPyvs7VvWcMMUS51ptqUne3tS5V5a9pYGcvrdPfkE+B5IPMLhBG+MCCf5A6efxVkfTQB+TeeJ6dTzIyi0oVZssGYCve/pj3py8jYeb7eGFzyYyc/yn1S5DWSV7GTlKjwtbM/d3DAn0sdjvrabnmrWU7Z1zOTvPbN6+ovKXPfe0mfXsKQXp9MH/SMoRsvOzCTtX+Oe46irbFshJFPXlLxvkZ4zoDJQTgN3cmDr7d3ZNbcnzLCRiw9WM6XOfxcplr4HGXpn7O9rb3ZKlu1Ga6oJqTlJKJjGzYivcZ8lzL2ZWrMlt2lMK0umD/75z+wDodP3dEH63PLPXQiqqZZf3JVGBgXwwbR0JH/TkrrVTadm6K92ayh2ZLTjC3ZI1xgqU19ZhSlX36QgXVaW1uRsU24qOjtbbt2+v8XZe2/Aa09dMJ+WJFOp617VAyYSlnE4+Ss8v+pObn8uWe7fQvG7zWtu3ow3IcWXmatEhgT5snH5ttbbZavoqs6mZ8lRln7Y6x5RSO7TW0RWt5/Q1/4TkBJr6NaHul0uM3Tyb116AEeVrHBTKo5H/4YnYYfR68SqaN/iM6UOjrP4FcaVRp87AGl14zbV1ACiwSM7e3lOQTt/bZ9+5fYR5N4d774WtW21dHFHCM8v/ZO5PWTz2RxfO+F7k9LmpTF8WZ/XucI44IMdVLY9Lwk0pk+/VJH8+bUgHPAymt6sBgxX2aW+cOvhrrUlITqBTlr9xQUSEbQskii2PS+LLzcfRSvFt5HNM29iEY3WPkH35OV7/X4JV9y2DwRxD0R1avonUtEXy5+XkffK1rrVuw9UdW1BTTp32SbyUSHpOOmEpBcYpnFu3tnWRRKHZqw8Uf/dy3D34sfNbTN18P+/0ikOdmwUMtNq+zd3y27pWZ6/tELYql7lGWYNSV4wVqGoZZ68+QG6B+egfUrgNax+3LVOQTh38E5IL5/Q5mAKdO4Oh/EEYovaUrWWn+vizpdUcxu96kEWRa3jjjzd4/OrHrbJve+yJYa/tEM8s/9N4h1b4ujbLZe5OrEDrKwJ/VX935d3lFZ0LtZGzLy8Fae19O3Xap2hCt047jkvKx86YqmWfDmjIllYLiWk6gmm/TOP9rfOssm97HHVqj+0Qxam5Mstrq1yVna6jOr87c9s2dVdhTbZMQTp18P9X1L/4Y8IfBO/YD88/b+viiBJMDYdXwE39woidsIQR9Xvz4E9T+HjdW1bZ/+ioEDZOv5Yjs4azcfq1Nk+v2OPEYCVTc2XVRnCq7HQd1Qmg5rb95i1X1eq5YMv5qJw67RPgFUDv5r1tXQxhQkWDi76NeJHRnwxmkn4MT4Mnd/edUurz9pofr67yuh6WnGcGai8NVF7wrI3gVNkBaNVpw6ntwW3mzldbpiCdfpAXK1dCXBw884zxwe3CYWT+8hMjvrqB2JYFzO8/m0nXGNsAlsclMW3JbnLz/zl3PQyK2WNrt9ZmSWXz1ubUZGBTVZkbXKWAObdG2s3v2tTvzsfDYPNUXpGKymfpiowM8iry7bcQGwvPPWfrkogq8hl0PSsLVjD2y5FMVtNIy7vMY4Nm8PzKvaUCP0Buvub5lcbnBDjiHYE9TgxmqlaqgDt7tbCr36k9TFGRm59Lek46l3Mvk5GbweUc4/8ZuRk8+eNmkgvS0IZsCshGq1xSdDZTVml+P9uYrLwsOoZlEZqfRXZeNh8nZLHuTFvmDp1r1TJbJPgrpYYCbwMG4GOt9awy73sBnwHdgPPArVrro5bYd4Xi46Wx14H5DBnOd76xjFs1kcf/mMkl9wIuZHRHceUgnIsZuXbZY6ay7G1iMHsIqpVV3Z45WmvSc9K5kHmBi1kXuZh5kZSslOJ/W46eYP3fx0nLTsXTM5sWDdzw9somLTuNtJw00rLTuJx7mZz8nPJ35HnlotQ8NxbsrIO3uzdeBi+83L3wdvfG292bhr4Nq3wsVVXj4K+UMgDzgEFAIrBNKbVCa72vxGoTgYta67ZKqduA14Bba7rvCuXmwr59Lv/MXkfn2ac/i68+gP/KSbyw/gVa5EajPZ5G4XHFurbqNmdp9tId1d6nKCiWmwvZ2WhfXy5mXeTMb6s4k3yMM2mnOJuZzNms85wLMJDc0J/kjGSS//6TZDK4YMghV5X/XAnfHDd8c9zQ2X7sz29Mx6AGdNh3Cv98d/zyvfEr8MNPu+Mb1RO/mGvxzYU68/+Dr5sXddy82H4kndxsAxtDe7C9WVfqZuYwYcdqfH19mTy4E3h5Gf/16QNhYXDpknFckpVZoubfAziktT4MoJT6GhgFlAz+o4CZhT8vAd5TSilt7QaHAweMJ0V4uFV3I6zP4Gbg45Ef02LV78wM2U7IpYdx83gNNyqepbVsDdoRGoutXeu2+9+B1sYgmJhofO726dNkexpIHNSTxEuJnHjzORLPHiKJNJI8Mjnpk8epYG9O1ykwWQt3K4AGZ7wIdmtDUJ0gOpzM4eoMAw3yAqif70X9Ak/qd42h3t2Tqeddj8B/PcCRI2mobDcUxl5Bq9v35suoYQRedGPpwbeMZczPh4IC4/8x10G3e+HcOdjxDuSmQG4uURlZXE7PJKVOOHHNAmmYkcijG78xFuznEoX88ENj8P/rL2PMcrduVt4SWw8BTpR4nQj0NLeO1jpPKZUKNACSLbB/8xITwcdH0j5OQinFjFf+IHhCDx6OPEJg1mQ83F7FgxZ4uCl8vdzNPkZyeVxSceOao6SGrFXrruh3UCsXBq3h/Hn4+284fBiOHCE3N4tjD93N4YuHOTLzYY6eSuBIPTgaCMfrwil/YH/h54OM/wLyPQghgBBDIO3rNqNJ51409mtM47OXaeQTTKN6zWnYoAX164dg8PP/J6DeU0H5fvqDAWZm/jx6uQB++cXsR5cn5jD7zrdL/f4A1q4+gErJJLtNO5bvOMHozsGQkwPZ2ZCVBXULZx3u2BE8TeSJLMwSwd/UDEhlf2eVWQel1CRgEkCLFi1qXrKhQyEtDcxM0iQcUFAQ//7vPoLvHsiDoRu56D6Fth6P88LIqQA8/M0ukx+buWIvo6NCbDqi0l5UNCjKohfHrCzjHfj+/XDsGLmPPcLhi4c5+Oy/+Wt3LIfqw8H6cKi+McAXvPui8XNtwaOtgZYewbT0C+H6eq1o0agdLYLb0bxuc5oFNKNZQDP8PK33fI7qdCE1d2F9dUy46V5aXl7g7196WS09c8QSwT8RKDlPcjPgpJl1EpVS7kBd4ELZDWmtFwALwNjV0wJlkykdnJG3Nzd/8ztXvzWTm0+8xaZ6r7P2dBavDXoNvjH9kaI7ApnUrfzfgbkLw2P/3Q2UcwHIz4dDh4zzZ3l4cH7B2yQseoP92Unsb6DZHwR/NYBDLz6JVgXQGGgM9Qx+tA0IpXejMMYFd6BN/Ta0rtea1vVa08SvCQY3231/q9Pu4kiVC0sE/21AO6VUKyAJuA24o8w6K4DxwCZgLBBr9Xy/cG5KEfLY86zLmc60NdN5Z+s7rNm5hBym4Uk7sx+zxaRu9pZfL+93YO7CkK916TuApCT48UfO7trI3qNb2Zt6iH2Buey7Ppp9Gcc4l3EOBhs/aygw4JPXCFRL/PNa4OvWjP/r3497e8fQoE4Dqx1nTZXX7mLub+pIlYsaB//CHP4UYDXGrp6faK33KqVeALZrrVcA/wE+V0odwljjv62m+xUCwNPTh7evf5vr4zOYePZjzvo9QuOM0RgM40v1BqpXx/hzbfeiscc2hvJ+B7NXHyh1YVC6gLbJJ+h49k/qZu5idkYov55qwJ7969l7+k+SGwKFvRLrKh86uStGdhhJWFAYYcFhPLvkAsmpAcWNpkVWbPPhievsN/AXMdXuUt7f1F5njDXFaUf42lttS1iZ1lz8cC4Pxk5jced86mfWxUM9jLfufsXo39o8N6rzCMLaKJ+5fXy77W/+34qfURnx9D62lHM+50kILiCxxBNQ/T396RLUic51WtK5dS+6NAqnU3Anmvg1QZVpXzP3uEQFHJk13KLHVFvK+5uau7DW5mhjlx7ha4+1LWFlSlHv/kf4auAIRv/faJ5tvZe/gp6nntvVzOz3Yqm/e232Xa9qGqC2zt1RkU3p0VoTv30V8bt/ZslXu3huUSIJgXnkuRdAAJzoBI3T/fHLC6F1akfSvcJpFhDG9um3XRHkzXGkmnBllfc3daSBcU4Z/B2p0UVYWNu23PLDn4xe/AVz6u7j5fh5PLxuIFuPjWTG8Nm0a2C+PcAaygt+pmrf1jh3s/KySDiXwO7E7ez++w/ic44Tfyae5Ix/elo3125EuAUzvF4PIq69nbMXGvLBr5lkeSoyC3sd1nc38OzQ8EoHfrCfwWqWVNEFzVEGxjll2qe8W805t0Y6xFVZWMb5jPPMfrwX79Q/RLa7YnToUB6+5kn6tOhTpSBWXeYm9bqpWwhLdyRdsdzcxG6VSZNorTmZdpL4M/HsPrOb3ce3En9iOwcyE8lXxm+ET56iS4toIhpFcNXBS0T4tiai71jqdep2RZdoS6WfnC0Fa+8Tybl02sfclbmuj4ekg1xMgzoNmPX0Oh558QnmHv6KBZk/sezYT0QFhjGh97+5rcttBNUJstr+zaUBzNXwDUqZfGZt2TRJRm4Ge8/uJf5MPH+e/ZP4M/HEn4zjfE5K8Tot8/yI+DudMefciPBvR3i7PrSLGYFhxChQiuVxSTy1+gAnPz9D08C1VgvKjlITrixHSu2Uxylr/uauzN4eblzMuHIEaG1Okyts6NgxMua8zhdbP+aDIfXZxWnc3dwZGjqQ0Z3HMqzdMJr4N6mVopi7O4XSdwCafAweZ7m5VwF1fE8VB/pDFw6hC7dQJ99A57MQeTKfiDOwvN8T/Ovmu7nTOw9SU6F7d/D2LrWPykwzbM+1W0uo6I7EUe9YKlvzd8rgD6b/cI98s8vpeh6IarhgHF/4Z24Sn3/7DIuPrCQxwHhmdGsURf9W19CnRR9iWsTUeHZFcwGkbI+RAjLJUyfx9ztL1zaZ/HowjtS8I+S5JaExVlgUirY5fkS06E54h36EH8sk4tHXaFC/PWv8WrGlSUd2hIRxuH4IPp7u5QbqinohVaeXkiNx5oufywd/U5z9hBbVEB+Pfu9d4jcsZVXQRf7XTrG1GWS7Gb8XTfya0KVhF7o07ELreq1pUbcFLeq2ILhOMPV96uPjUf5Q/+nLdnI5N40ClUa+SsHgnsawq7xIzjzJmoMJZHOaPHWGAvVPukahaBnYks6+rej0cxydDqYQfgbCkqFOcFN47z248Ubj1Ak5OcS8v73K53VFXTCdsYtmSc588XPpnL85ztjzQNRQRARqwUdclT+fq/74g6e+/57s7X+y44Pn2JS4mT+XzGPP8d+ZHxhLptuVjbFeBi98PHzwNHji4eaBRpNXkEdufi4pWelo99wrvmUf7zF+rp5fY1qeuEzoeQOdkwPoejKdsLMFtB3/f/jMmAMZGbB6PIyJgG7doGtXaNz4nw15e4O3t9nHP1b0GMbyeqw4YxfNkirqgutII3Wry6WCv7M01AgrMBigb1/o2xcv4Grg6hYxsNML/vc/CuJ2ci71FMcC4XjXNpx/bhoXsy5yYeEHZKddJN2QTQpZ6AJI92uAR7/+bPgrnXFx6wi5lE7IpRyaphcQnAG7G/fm7p0bjb2NWreG+vWhTRvo3o7tPo2452Jdtk5fZTw/n5pb7vm5PC4JhYlZEik/UFdUEXL2ipKrX/zAxYI/OF/PA2FlU6bAlCm4AY0uX6bR4cP0uHwZonsZ3199meNn9rPv79N4ZmfhpjXxjdvyQZNbaeXhRuS5POrkZpFax48tDQJI8QkgJbQd44u6VR4+XLyr4jyze+V7o81efcBsw/E1HYPNHlZFFSFnryi5+sUPXCznL4Q1mMsP16vjQVZuQaUbDauTZy6v15Aj5KdtydV7+7hczV84r7Jf1ms6BrN2/zmrf3nN5YFTMnKrNKiwOnlmc+mJij4nKs4COHuWQIK/cAqm5sT5YvPx4vetOaCvvPxwVQJIdfLM5XVhtpf8tKPWoC3FXo/fzdYFEMISTI2YLavk06osadqQDvh4lJ6yuDr54epsZ3RUCHf2anHFo/LsJT9ddFFOSslE889FeHlckq2LVivs+fgl+AunUNkUhzVSIaOjQnh1TDghgT4ojLn26gwGqu52XhodzpxbI2u8f2uo6JGRzs6ej1/SPsIplJf7LrueNVgqP1zd7dhrftoR+8tbMk1jz8cvNX/hFEylTMqyl1SIKzF3sbWX9oiyLJ2msefjl+AvnIKplMm4Xi3sMhXiSizVHlJbLJ2msefjl7SPcBr2mvpwZY42WMzSaRp7Pn4J/kIIq3Kki7I1pnWw1+OXtI8QQhSy5zSNpUnNXwghCtlzmsbSJPgLp2GvIymFY7HXNI2lSfAXTsHU9A7yfGbbkwuy/ZKcv3AK9jyS0lXZ89QGQoK/cBL2PJLSVckF2b5J8BdOwZ5HUroquSDbNwn+winYYxe95XFJxMyKpdX0VcTMinW5dIdckO2bNPgKp2BvXfQcuQHaUo201XkUojQQ1x55jKMQVlCdRzLag7IXLSj/0ZOV2V5lg7ml9+2q5DGOQtiQo+a7y2uktfZU05betyifBH8hrMBSc8TUdhrElhctR71gOioJ/sJmHC2/W5XyViffbWp/td1uYI2JzRxh365IevsIm3C0AUBVLa8lHu1oi37ytuw1ZY89tpyZ1PyFTThafrc65a3pHDG2SIPYsteUvfXYcnYS/IVNOFp+1xbltVUaxJYTm7nKpGr2oEZpH6VUfaXUL0qpg4X/1zOxTqRSapNSaq9SKl4pdWtN9imcg6MNALJFeSUNIqyppjn/6cAarXU7YE3h67IygLu11p2BocBcpVRgDfcrHJyjBTZblNcS7QZCmFPTtM8oYEDhz4uAdcATJVfQWv9V4ueTSqmzQDCQUsN9CwfmaPldW5VX0iDCWmo0wlcplaK1Dizx+qLW+orUT4n3e2C8SHTWWheYeH8SMAmgRYsW3Y4dO1btsgkhhCuy2AhfpdSvQGMTbz1dxQI1AT4HxpsK/ABa6wXAAjBO71CV7Qshap+jjdUQ/6gw+GutB5p7Tyl1RinVRGt9qjC4nzWzXgCwCnhGa7252qUVQtgNR568TtS8wXcFML7w5/HA92VXUEp5At8Bn2mtv63h/oQQdkIe1uLYahr8ZwGDlFIHgUGFr1FKRSulPi5c5xagH3CPUmpX4b/IGu5XCGFjjjZWQ5RWo94+WuvzwHUmlm8H7i38+Qvgi5rsRwhhf2QuHscmc/sIIarF0cZqiNJkegchRLU42lgNUZoEfyFEtckgNMclaR8hhHBBEvyFEMIFSfAXQggXJMFfCCFckAR/IYRwQRL8hRDCBUnwF0IIFyTBXwghXFCNHuZiTUqpc4AlnuYSBCRbYDuOxNWOWY7X+bnaMdfkeFtqrYMrWslug7+lKKW2V+apNs7E1Y5Zjtf5udox18bxStpHCCFckAR/IYRwQa4Q/BfYugA24GrHLMfr/FztmK1+vE6f8xdCCHElV6j5CyGEKMNpgr9SaqhS6oBS6pBSarqJ972UUt8Uvr9FKRVa+6W0nEoc76NKqX1KqXil1BqlVEtblNOSKjrmEuuNVUpppZRD9w6pzPEqpW4p/DvvVUp9VdtltKRKnNMtlFJrlVJxhef1MFuU01KUUp8opc4qpfaYeV8ppd4p/H3EK6W6WrQAWmuH/wcYgL+B1oAnsBvoVGadfwPzC3++DfjG1uW28vFeA9Qp/PkBRz7eyh5z4Xr+wHpgMxBt63Jb+W/cDogD6hW+bmjrclv5eBcADxT+3Ak4auty1/CY+wFdgT1m3h8G/AQooBewxZL7d5aafw/gkNb6sNY6B/gaGFVmnVHAosKflwDXKaVULZbRkio8Xq31Wq11RuHLzUCzWi6jpVXmbwzwIvA6kFWbhbOCyhzvfcA8rfVFAK312VouoyVV5ng1EFD4c13gZC2Wz+K01uuBC+WsMgr4TBttBgKVUk0stX9nCf4hwIkSrxMLl5lcR2udB6QCDWqldJZXmeMtaSLGGoQjq/CYlVJRQHOt9Q+1WTArqczfuD3QXim1USm1WSk1tNZKZ3mVOd6ZwDilVCLwI/BQ7RTNZqr6Pa8SZ3mGr6kafNluTJVZx1FU+liUUuOAaKC/VUtkfeUes1LKDZgD3FNbBbKyyvyN3TGmfgZgvLP7XSnVRWudYuWyWUNljvd2YKHW+k2lVG/g88LjLbB+8WzCqjHLWWr+iUDzEq+bceUtYfE6Sil3jLeN5d1y2bPKHC9KqYHA08BIrXV2LZXNWio6Zn+gC7BOKXUUY450hQM3+lb2nP5ea52rtT4CHMB4MXBElTneicB/AbTWmwBvjHPgOKtKfc+ry1mC/zagnVKqlVLKE9LekwsAAAE1SURBVGOD7ooy66wAxhf+PBaI1YWtKg6owuMtTIF8iDHwO3IuuEi5x6y1TtVaB2mtQ7XWoRjbOUZqrbfbprg1VplzejnGhn2UUkEY00CHa7WUllOZ4z0OXAeglArDGPzP1Wopa9cK4O7CXj+9gFSt9SlLbdwp0j5a6zyl1BRgNcZeA59orfcqpV4AtmutVwD/wXibeAhjjf8225W4Zip5vLMBP+Dbwnbt41rrkTYrdA1V8pidRiWPdzUwWCm1D8gHpmmtz9uu1NVXyeN9DPhIKfUIxvTHPQ5cgUMptRhjyi6osB1jBuABoLWej7FdYxhwCMgA/mXR/Tvw704IIUQ1OUvaRwghRBVI8BdCCBckwV8IIVyQBH8hhHBBEvyFEMIFSfAXQggXJMFfCCFckAR/IYRwQf8fFffvnMmIE4gAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "d = 3\n", "X = generateSplineRegressionMatrix(x, xi, d)\n", "beta = np.linalg.solve(X.T @ X, X.T @ y)\n", "plt.scatter(x, y, label = 'Samples')\n", "Xspline = generateSplineRegressionMatrix(xlin, xi, d)\n", "plt.plot(xlin, y_func(xlin), 'r--', label='Population line')\n", "plt.plot(xlin, Xspline @ beta, 'g-', label='Regression line')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Hint**: If you have no idea how to get started, the following code cells might help you.\n", "They perform a **global cubic regression** and can be adapted to perform a **cubic spline regression**." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def generatePolynomialRegressionMatrix(x, d = 3):\n", " \n", " # Get number of samples\n", " n = len(x)\n", " \n", " # Initialize matrix with predictor variables 1 (for the intercept)\n", " X = np.ones((n,1))\n", "\n", " # Append columns for the monomials x^1, ..., x^d\n", " for i in range(d):\n", " X = np.hstack((X,np.power(x,i+1).reshape(n,1)))\n", " \n", " return X" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Set degree of spline regression\n", "d = 3\n", "\n", "# Generate Spline matrix\n", "X = generatePolynomialRegressionMatrix(x, d)\n", "\n", "# Solve the normal equation, remember the @ sign\n", "# performs the ordinary matrix multiplication for numpy arrays.\n", "# You should also avoid to invert the matrix X^T * X!\n", "\n", "beta = np.linalg.solve(X.T @ X, X.T @ y)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XlcVNX7wPHPYccNVDQVNdx3BHENF8q1zCXLNLPsq2lWZmlZ2qJlm2WmLfYrK9M2ra+ae19zSXPJXMkdNVfQ3EGRRZbz++MCAs7AADPMMPO8Xy9fMjN37j13GJ577lmeo7TWCCGEcC1u9i6AEEKI4ifBXwghXJAEfyGEcEES/IUQwgVJ8BdCCBckwV8IIVyQBH8hhHBBEvyFEMIFSfAXQggX5GHvApgTEBCgg4KC7F0MIYQoUXbu3HlRa10pv+0cNvgHBQWxY8cOexdDCCFKFKXUSUu2k2YfIYRwQRL8hRDCBUnwF0IIF+Swbf5CCNtLSUkhOjqapKQkexdFFJCPjw/Vq1fH09OzUO+X4C+EC4uOjqZs2bIEBQWhlLJ3cYSFtNZcunSJ6OhoatWqVah9SLOPEC4sKSmJihUrSuAvYZRSVKxYsUh3bBL8hXBxEvhLpqL+3iT4CyGEC5LgbyeLd8cQPmUdtcavIHzKOhbvjrF3kYSwC3d3d0JCQmjatCm9evUiNjbW3kXK4Z577rFKmV5//XU++OADACZOnMiaNWuKvM+ikOBvB4t3xzBh0V5iYhPRQExsIhMW7ZULgHBJvr6+REZGsm/fPipUqMDMmTOtst/U1FSr7GflypX4+/tbZV+ZJk+eTJcuXay6z4KS4G8HU1dFkZiSluO5xJQ0pq6KslOJhHAM7dq1IybmZiVo6tSptGrViuDgYCZNmpT1/JtvvknDhg3p2rUrDz30UFaNOiIigpdffplOnTrx0UcfceHCBe6//35atWpFq1at2Lx5MwAbNmwgJCSEkJAQQkNDuXbtGmfPnqVjx45ZdyEbN24EjFQzFy9eBODDDz+kadOmNG3alBkzZgBw4sQJGjVqxPDhw2nSpAndunUjMTExz/N87LHHWLBgQdb+J02aRIsWLWjWrBmHDh0C4Pr16wwdOpRWrVoRGhrKkiVLrPERZ5GhnnZwJtb0F8Pc80IUh+f+9xyR/0ZadZ8hVUKY0WOGRdumpaWxdu1ahg0bBsBvv/3GkSNH2LZtG1prevfuzR9//EGpUqVYuHAhu3fvJjU1lRYtWhAWFpa1n9jYWDZs2ADAoEGDGDNmDO3bt+fUqVN0796dgwcP8sEHHzBz5kzCw8OJj4/Hx8eHWbNm0b17d1555RXS0tJISEjIUb6dO3fyzTff8Ndff6G1pk2bNnTq1Iny5ctz5MgR5s2bx5dffsmDDz7IwoULGTx4sMWfU0BAALt27eKzzz7jgw8+4KuvvuLtt9/mrrvuYvbs2cTGxtK6dWu6dOlC6dKlLd5vXiT420E1f19iTAT6av6+diiNEPaVmJhISEgIJ06cICwsjK5duwJG8P/tt98IDQ0FID4+niNHjnDt2jX69OmDr6/x99KrV68c+xswYEDWz2vWrOHAgQNZj69evcq1a9cIDw9n7NixPPzww/Tr14/q1avTqlUrhg4dSkpKCn379iUkJCTHfjdt2sR9992XFXz79evHxo0b6d27N7Vq1craPiwsjBMnThToM+jXr1/WexctWpR1/kuXLs26q0lKSuLUqVM0atSoQPs2R4K/HYzr3oAJi/bmaPrx9XRnXPcGdiyVcHWW1tCtLbPNPy4ujnvvvZeZM2cyevRotNZMmDCBJ554Isf206dPz3N/2WvG6enp/Pnnn1kXikzjx4+nZ8+erFy5krZt27JmzRo6duzIH3/8wYoVK3jkkUcYN24cjz76aNZ7tNZmj+nt7Z31s7u7e77NPube7+7untVXobVm4cKFNGhgm7ggbf520Dc0kHf7NSPQ3xcFBPr78m6/ZvQNDbR30YSwGz8/Pz7++GM++OADUlJS6N69O7NnzyY+Ph6AmJgYzp8/T/v27Vm2bBlJSUnEx8ezYsUKs/vs1q0bn376adbjyEijWeuff/6hWbNmvPTSS7Rs2ZJDhw5x8uRJKleuzPDhwxk2bBi7du3Ksa+OHTuyePFiEhISuH79Or/88gsdOnSwwSdh6N69O5988knWRWf37t1W3b/U/O2kb2igBHshcgkNDaV58+bMnz+fRx55hIMHD9KuXTsAypQpw/fff0+rVq3o3bs3zZs35/bbb6dly5b4+fmZ3N/HH3/M008/TXBwMKmpqXTs2JHPP/+cGTNm8Pvvv+Pu7k7jxo25++67mT9/PlOnTsXT05MyZcrw7bff5thXixYteOyxx2jdujUAjz/+OKGhoQVu4rHUa6+9xnPPPUdwcDBaa4KCgli+fLnV9q/yupWxp5YtW2pZzEUI2zp48KDV2pCLU3x8PGXKlCEhIYGOHTsya9YsWrRoYe9iFTtTvz+l1E6tdcv83is1fyFEiTNixAgOHDhAUlISQ4YMccnAX1RWCf5KqR7AR4A78JXWeoqJbR4EXgc08LfWepA1ji2EcD0//vijvYtQ4hU5+Cul3IGZQFcgGtiulFqqtT6QbZt6wAQgXGt9RSlVuajHFUIIUXjWGO3TGjiqtT6mtb4BzAf65NpmODBTa30FQGt93grHFUIIUUjWCP6BwOlsj6MznsuuPlBfKbVZKbU1o5lICCGEnVijzd9UUuncQ4g8gHpABFAd2KiUaqq1zpEqTyk1AhgBULNmTSsUTQghhCnWqPlHAzWyPa4OnDGxzRKtdYrW+jgQhXExyEFrPUtr3VJr3bJSpUpWKJoQwtFlT+ncv3//W3LqFNWcOXMYNWpUntusX7+eLVu2ZD3+/PPPbxnnXxgnTpygadOmAOzYsYPRo0cXeZ/WYo3gvx2op5SqpZTyAgYCS3Ntsxi4E0ApFYDRDHTMCscWQpRw2VM6e3l58fnnnxd7GXIH/5EjR+ZI7WANLVu25OOPP7bqPouiyMFfa50KjAJWAQeBn7XW+5VSk5VSvTM2WwVcUkodAH4HxmmtLxX12EII59KhQweOHj0KmE+f3LBhQ4YMGUJwcDAPPPBA1p1C9tTLO3bsICIi4pb9L1u2jDZt2hAaGkqXLl04d+4cJ06c4PPPP2f69OmEhISwcePGHAuvREZG0rZtW4KDg7nvvvu4cuUKYKSPfumll2jdujX169fPSgFtzvr167n33nsBY2GXoUOHEhERQe3atXNcFL7//ntat25NSEgITzzxBGlpaeZ2WSRWye2jtV6pta6vta6jtX4747mJWuulGT9rrfVYrXVjrXUzrfV8axxXCGFlERG3/vvsM+O1hATTr8+ZY7x+8eKtrxVAamoqv/76K82aNcuRPnnr1q18+eWXWbltoqKiGDFiBHv27KFcuXJ8llk+C7Rv356tW7eye/duBg4cyPvvv09QUBAjR45kzJgxREZG3pKv59FHH+W9995jz549NGvWjDfeeCNHmbdt28aMGTNyPG+JQ4cOsWrVKrZt28Ybb7xBSkoKBw8e5KeffmLz5s1ERkbi7u7ODz/8UKD9WkoSuwkh7CozpXPLli2pWbMmw4YNy5E+uUyZMlnpkwFq1KhBeHg4AIMHD2bTpk0WHys6Opru3bvTrFkzpk6dyv79+/PcPi4ujtjYWDp16gTAkCFD+OOPP7Jez56KuaA5fnr27Im3tzcBAQFUrlyZc+fOsXbtWnbu3EmrVq0ICQlh7dq1HDtmmxZySe8ghLhp/Xrzr5UqlffrAQF5v25GZpt/dnnlHFNKmXzs4eFBeno6YOS+N+WZZ55h7Nix9O7dm/Xr1/P6668XuLzZmUrFXND3Zn+/1pohQ4bw7rvvFqlclpCavxDC4eSVPvnUqVP8+eefAMybN4/27dsDRpv/zp07AVi4cKHJ/cbFxREYaExDmjt3btbzZcuW5dq1a7ds7+fnR/ny5bPuOr777rusuwBb6Ny5MwsWLOD8eWMe7OXLlzl58qRNjiXB3wYW744hfMo6ao1fQfiUdbIwuxAFlD19cps2bbLSJwM0atSIuXPnEhwczOXLl3nyyScBmDRpEs8++ywdOnTA3d3d5H5ff/11+vfvT4cOHQgICMh6vlevXvzyyy9ZHb7ZzZ07l3HjxhEcHExkZCQTJ0600VlD48aNeeutt+jWrRvBwcF07dqVs2fP2uRYktLZyhbvjjG5Spcs1iIcUUlL6XzixAnuvfde9u3bZ++iOISipHSWmr+VTV0VlSPwAySmpDF1VZSdSiSEELeS4G9lZ0wszJ7X80IIywUFBUmt30ok+FtZNX/fAj0vhBD2IMHfysZ1b4CvZ87OJl9Pd8Z1b2CnEgkhxK1knL+VZXbqTl0VxZnYRKr5+zKuewPp7BVCOBSp+QshhAuS4G9lmUM9Y2IT0UBMbCITFu2Vsf5CmPH222/TpEkTgoODCQkJ4a+//rLZsSIiIiiJQ8htQZp9rCyvoZ7F2fSzeHeMND0Jh/fnn3+yfPlydu3ahbe3NxcvXuTGjRv2LpZLkJq/lTnCUE+5+xC2Yu3Z62fPniUgICArz01AQADVqlVj8uTJtGrViqZNmzJixIisXD8RERGMGTOGjh070qhRI7Zv306/fv2oV68er776KpB32ufsfvvtN9q1a0eLFi3o378/8fHxAIwfP57GjRsTHBzMCy+8UKTzc2QS/K0k84/C3Hzp4hzqKRPNhC3YolLRrVs3Tp8+Tf369XnqqafYsGEDAKNGjWL79u3s27ePxMREli9fnvUeLy8v/vjjD0aOHEmfPn2YOXMm+/btY86cOVy6ZCwTkl/a54sXL/LWW2+xZs0adu3aRcuWLfnwww+5fPkyv/zyC/v372fPnj1ZFxRnJMHfCrL/UZji6aZIuJFabLl+HOHuQzgfW1QqypQpw86dO5k1axaVKlViwIABzJkzh99//502bdrQrFkz1q1blyP1cu/exhpRzZo1o0mTJlStWhVvb29q167N6dOngfzTPm/dupUDBw4QHh5OSEgIc+fO5eTJk5QrVw4fHx8ef/xxFi1aRKlSpQp9bo5O2vytwNQfRSZ/X0+u30jlSkIKcLO2BFjUBl+Ytvtq/r4mL0Qy0UwUha0qFe7u7kRERBAREUGzZs344osv2LNnDzt27KBGjRq8/vrrOVI0ZzYRubm55UiL7ObmlpVW2Vza50xaa7p27cq8efNuKc+2bdtYu3Yt8+fP59NPP2XdunVFOj9HJTV/KzD35VdAaW8PUtJyNgZZWlsq7G22TDQTtmCL2etRUVEcOXIk63FkZCQNGhjf04CAAOLj41mwYEGB92su7XOmtm3bsnnz5qwlIxMSEjh8+DDx8fHExcVxzz33MGPGjFvWGXAmUvO3grxq2kWpLRV25JBMNBO2MK57A5MZa4tSqYiPj+eZZ54hNjYWDw8P6taty6xZs/D396dZs2YEBQXRqlWrAu83M+3zE088Qb169bLSPmeqVKkSc+bM4aGHHiI5ORmAt956i7Jly9KnTx+SkpLQWjN9+vRCn5ujk5TOVpBXGuepq6JMXhgC/X3ZPP6uPPdba/wKkx3ICjg+pWcRSy1EwVM6l4QhxK6U9rkoKZ2l5m8F+dW0C1tbkrZ74Wj6hgY6XLAXhSPB30rM/VEUpQnGFrfZQjg7SftsGQn+xSCv2lL222j/Up5oDXGJKVkXicymI0e+zRYlm9b6ltEwwvEVtcleRvvYUe7RPFcSUohNTMkxsgdg8/i7mD4gBIAxP0XKusDCanx8fLh06VKRA4koXlprLl26hI+PT6H3ITV/O8prfgDkHBKavfmnoHMFhDCnevXqREdHc+HCBXsXRRSQj48P1atXL/T7JfjbkbkZwdmdiU10mGRxwvl4enpSq1YtexdD2IE0+9iRuwXtrEWdKyCEEKZI8C8uN27A//4HV68aj9et45v5r/LO/z7h7kObKJcUf8tbMkf2yLrAQghrk+Bva1euwNNPQ+XKcPfdsGyZ8fyNGwSkJnDvwY3835Ip7P54EN/8dxKV4q8AxiSwd/s1o29ooKRrEEJYnbT529KKFTBiBJw7Bw8/DP37Q5cuxms9enB46VoG/Hc39U8dJOLYTjqc2EVyOT9m9A+hb3AVcDcCvqRrEEJYm6R3sJX0dAgPh/h4mDMHwsJMbpZjuryfD+N6NKRvUCno0AHGjYMhQ4q33EKIEk3SO9hLWhqkpoK3NyxaBBUqGD+bYXIC2LlzEBAAjz0Gf/wBn3wCTpxXXAhR/KTN35q0hiefhJ49ISUFqlbNM/CbddttsGYNvPIKzJ4N7dvD+fPWL68QwmVJ8Lem55+HL7+ENm3A07No+/LwgLfeguXL4dAhGDvWOmUUQgik2cd6Fi6E6dPhmWeMoF0AeabJ7dkT1q2DevVsUGghhKuySs1fKdVDKRWllDqqlBqfx3YPKKW0Uirfzgh7y1yQ3aJ1d8+ehSeegJYtYdo0KECSLItW62rbFipWhORkGDoUoqMLf2JCCIEVgr9Syh2YCdwNNAYeUko1NrFdWWA08FdRj2lrBV4+8eJFCAyE774rcHNPgRbFPnECFiww5gvExhboOEIIkZ01av6tgaNa62Na6xvAfKCPie3eBN4Hkky85lAKFJABmjWDyEho2LDAxypQ6oYGDeCXXyAqCvr2Ne4EhBCiEKwR/AOB09keR2c8l0UpFQrU0Fovt8LxbM7igHz2rNHJm5BQoKae7AqcuqFzZ2PewIYNxhwAB52nIYRwbNYI/qaiXlZEUkq5AdOB5/PdkVIjlFI7lFI77Jli1uKA/Oab8PHHxkWgkAqVumHQIHjvPVi71mgKEkKIArJG8I8GamR7XB04k+1xWaApsF4pdQJoCyw11emrtZ6ltW6ptW5ZqVIlKxStcO5sWOmWK9otAfnoUWNY54gRUKdOoY/VNzSQd/s1IzDjwuKuVFYTU56dzOPGwf79IOl4hRCFYI3gvx2op5SqpZTyAgYCSzNf1FrHaa0DtNZBWusgYCvQW2vtkLkbFu+OYeHOGLI3pijg/rBcM3EnTgQvL3j11SIfM3vytrSMZpx8O5mVMpLFpacbdwHHjxe5HEII11Hk4K+1TgVGAauAg8DPWuv9SqnJSqneRd1/cckc2vncT5G3dPZq4PdD2ZqhIiNh3jx47jljFq8VFLiTOdOZMzBlCtx/PyQ5fF+6EMJBWGWcv9Z6pda6vta6jtb67YznJmqtl5rYNsLRav3Zh3aak6Ozt3Rpo9193DirlaEwC7Ys3h1D+PeHGdZ5NOzezT9DR1mtPEII5ybpHch/LV3I1dlbrx788AP4+xdsMpil+7fg+ewXrLV1W/NNWC/qzPuaPz/5tlDHF0K4FqcO/gkpCVxLvkZ+aavzWw4xR2fvnDmwbx9QiMlgeSjoqJ/cF6wpEf/hYKUg6r0yVpp/hBD5curcPt/9/R0jV4zEy92LgFIBVCpVierlqlOjXA1q+tWkboW61KtYj9v84N840/sIzJ5r5/x5I43D8OHw6adWXVi9oAu25L5gJXt4Mar3S/glx7PIx6dAxxZCuB6nDv5tq7dlatepXLh+gYsJFzl3/RzRV6PZGr2VS4mXcmzr6VMFj/Tb8UwPwju9DuXcG/FBvztzBt9Zs4y1eEcZbevWXljdZG5/M6r5+97SR/FPQI2sIaP8+y9UqVKocgghnJ9TB//mVZrTvEpzk69dv3Gdo5ePcvjSYaIuRbEqajs7z0ZyNX07qHQuAE+tq8r3R+6gQ80OdKjWltD/+z9U165ZaRxMBeDM521tXPcGTFi0N8edR1Yz0axZRgroyEioW9fmZRFClDyyjGMuiSmJRP4byfYz29kWs41NpzZxMu4kD+6DnxbAuy/eQfn+j9Cjbg8ij3uaDMCZC6/bmtlU0NHR0LSpkXNo/fqstYCFEM7P0mUcJfhb4HTcac698zK3/byCDqNKczLeSKncpFIT6vt15OjJRly9GkSgf2nHWVh97lxjGchp02QhGCFciAR/W0hPRytF1KUoVh5ZyYojK/jj5B+kpqdStUxV7mt4H/2b9KdDzQ64u9m5tq019OkDv/0Ge/ZA/fr2LY8QolhI8Lemc+eMVAomMnfGJcWx4sgKFh1cxK9HfyUhJYGqZaoyoMkAHg5+mLCqYahCZvwssrNnoXVrI/ncfffZpwxCiGIlwd9a0tOhdm3o1s3oSM1DQkoCk1Z9y+zd33M57S9QqdQoW59RbYbxSPAjVC1rnVQQBZKcXLhF5IUQJZKlwd+pJ3lZxZYtcPIkdOyY76a/7bvC8r+CKHt9AjWSvqfCjVFcvOrOS2teosb0GvSd35cVh1eQlp73bGKr8vY2moB++AFiCjf7WAjhfJx6qKdV/Pgj+PoaK2flI/ukLzfKUDatB2XTelCh3EW6tjrAnL/nsCRqCTX9ajIybCTDWgyjcunKtj4DI+g//jh0726sBGavZighhMOQmn9eUlLg55+NjtMyZfLd3NzkritXA3iv63ucHnOaBf0XUK9CPV5e9zI1ptfg0V8eZdfZXdYueU7Vq8PkybBkiRH8hRAuT4J/Xn77DS5dgocftmjz/JKzebl7cX/j+1nz6BoOPn2QES1G8MuhXwibFUanOZ1YfGgx6Tq9UEXNN8Hcc89B8+YwejRcu1aoY1h8LCGEw5Pgn5c774T5843OXgsUJDlbw4CGfHLPJ0SPiWZat2mcjD3JfT/dR+OZjflq11ckpVqenM2iBHOenvDFF0b+/9des3jfhTqWEMLhSfDPS6lSMGCAsWKXBbIvyagwksLlN9vXz8ePse3GcnT0UebfP5/SXqUZvmw4tT6qxbQt04i/EZ/vcS1eCKZNG5g0CTp1suh8inQsIYRDk6Ge5mzYYIz0GT3aWLylmGitWXt8Le9uepd1x9dRwbcCz7V5jtFtRuPn42fyPbXGr8Dcb3HGgBCrzjg2dywFHJ/S02rHEUIUjgz1LKo5c+D994t9jLxSii61u7D20bX8OexP7qhxBxPXTyTooyAmb5hMbFLsLe/JK5GcySaZtDR4++185y2YUtBFZ4QQjkmCvynp6bBiBdx9N3jYbzRs2+ptWfbQMnaO2Emn2zsxaf0kan1Ui3c2vpOjOchUX0Mmk00ybm7Gnc24cUbq5wIo6KIzQgjHJMHflG3b4MIF6NXL3iUBoEXVFiweuJidI3bSvmZ7Xln3CrU/qs30P6eTlJqU1ddgzi1DUJWCTz81Vvx64YUClaUw/RpCCMfj1G3+ZlMe5+eVV+C994wLQPnyRSqDLWyN3sprv7/GmmNrqFGuBm9EvMGjzR+l4/sbTK4vEOjvy+bxd926o4kT4c03Yd06Y2STEKLEc/k2/yINSYyNhc6dHTLwg9EctPqR1ax9dC1VylRh6NKhBH8eTETz0/h45vyV5tkkM2EC1KoFzzxjNHUJIVyG09b8w6esK1gtOLf0dKNt3MFprVl0cBET1k7gyOUjNK3YDrerj3D1ak3L7nY2bjTSV7TMt6IghCgBLK35O21un0Kvr5uWZqx8VQICPxijg+5vfD+9G/Tmy11f8vr617mQ8hSD2wzm3c7vUr1cPs1cHTrc/Dnz3ClCk5kQokQoGRGuEAo9JLFPHxg82AYlsi1Pd0+eavUUR0cf5eX2L/Pf/f+l/if1mfj7RK7fuJ7/DkaPhoEDAZnFK4QrcNrgX6ghifHxsHq1sXBLCVXOuxxvd36bqFFR9GnYhzf/eJMGnzbghz0/kGcT3223wYIFsGqVzOIVwgU4bfAv1JDE9evhxg3oWfJnqt7ufzvz7p/Hpv9sokqZKgz+ZTB3zL6DHWfM9KO88ALUrQvPPMPFi1dNbpJvk5kQosRw2uAPxgVg8/i7OD6lJ5vH35V/m/W6dcaM3vDw4ilgMQivGc624duY3Xs2x68cp/WXrRmxbAQXEy7m3NDb2xj7f+QIY/cuM7kvmcUrhPNw6uBfYOvWGYHfx8feJbEqN+XGf0L/Q9SoKJ5r+xyzd8+m3if1+L/t/5dzVbHu3eH++3ls+xL8Sc2xD5nFK4RzkeCfSWsYNgyefNLeJbEZPx8/Puz+IXue3ENolVCeWvkUbb5qw7aYbTc3+vhjvPdE8vqAljKLVwgn5rTj/EXetNb8tP8nxq4ay7/x/zK8xXDe7fIuFXwrZG4A585BlSoW7U+GhgrhGFx+hm+Bbd8OZ8/auxTFRinFwKYDiRoVxZi2Y/h699c0/LQhcyPnGqOCnnzSaAJLyn9RGRkaKkTJI8E/0yOPwPDh9i5FsSvrXZZp3aexc8RO6laoy2NLHuPOuXdyqlsbOHbMyHGUDxkaKkTJ47QzfM0x2TxRGYiKghEj7F08u2lepTmbhm7iq11f8dKal6h7+gl23tmExu+8w0PxddjuXt5sc06hZ1MLIezGpWr+5pondnyz0NjgLgty/jgxN+XGiLARHHr6EA82eZDuofu5TgpPLHoTrbXZ5hxZ4EWIkselgr+55ol/F/8KFSpAcLCdSuZYbitzG9/3+57yFd/jjYgyNLhwGPfkqaRxzWRzjizw4toW744hfMo6ao1fQfiUddLXU0JYJfgrpXoopaKUUkeVUuNNvD5WKXVAKbVHKbVWKXW7NY5bUOaaIUKO7OL3ak1Y/LfrdPhaIuFqExY1/4o2T/ThmP8mzvg8yXX3jcTEJuTYzt4LvEjwsR/p7C+5itzmr5RyB2YCXYFoYLtSaqnW+kC2zXYDLbXWCUqpJ4H3gQFFPXZBVfP3NZnmeeBD7+CdmsKZRXsBZIhiBuPzAk+Gc/v1jjT8dxqr6r5HmtsfRF8NoXq56lnb9g0NtMvnlhl8Mu/oMoNPZpmEbeXV2S+fv2OzRs2/NXBUa31Ma30DmA/0yb6B1vp3rXVmdXErUB07MLfWbbR/Ff4JqCEjVHLJ/nk9v3kby3/4l45n+nFd7abxzMZ8vuNz0rV9F4GRkUb2JZ39JZc1gn8gcDrb4+iM58wZBvxqheMWWPbmiUz996ym3761WY/lS3tT9s9rdqs+XPMtw0e/HqJy/Me4pdTlyRVPcufcOzly6YjdyijBx76ks7/kskbwVyaeMzltWCk1GGgJTDXz+gil1A6FQumfAAAgAElEQVSl1I4LFy5YoWi3ykz2lnkBGLZ9Mb0P/JH1unxpc8r8vCY91pH3Oz9OyOkDDPp7L34Jk6mS9hw7z/xN8OfBTN08ldT01Px3aGUSfOxLOvtLLmsE/2igRrbH1YEzuTdSSnUBXgF6a62TTe1Iaz1La91Sa92yUqVKViiaeeO6N6ByWiL1L55iZ2BDQL60eZm6Kop5je7krxpNmbD+GyomXMX7RhcaqS/pUbcHL655kXZft2Pvub3FWi4JPvZl785+UXjWmOS1HainlKoFxAADgUHZN1BKhQJfAD201uetcMwi6xsaSOUtybih2RXYiEDJR5OnM7GJoBSvdHuKqSs/wi8pnsul/LgYV4ptDy7ivwf+y6iVowibFcarHV9lQvsJeLp72rxcmb8vyStkP/bq7BdFY5XEbkqpe4AZgDswW2v9tlJqMrBDa71UKbUGaAZkjqU8pbXundc+iyWx2+uvw5tvQmwslC1r22OVcOFT1t0cKaU1KKO1L9Dfl83jjclxF65f4Nn/Pcu8ffNofltzvunzDaFVQ+1VZCFcUrEmdtNar9Ra19da19Fav53x3ESt9dKMn7torW/TWodk/Msz8BebmBhjYpcE/nzlaF5RCr/Ea7z6xxzGd6yZtU2l0pX48f4fWTxgMeeun6PVl614bd1rJKeabOUTQtiRS83wvcWXX8Jff9m7FCVC7rbd9olnefzPBfRa9vUt2/Zp2If9T+3n4eCHeWvjW4TNCjO/fKQQwi4kn78ovGHDYO5c2LkTmjc3ucmKwyt4YvkT/Bv/Ly+Gv8ikTpPw9vAu5oIK4Tokn39+vv8e7rkHrl2zd0lKrqlToWJF4yKQanqYZ8/6Pdn31D6GNB/Cu5vepcWsFmyP2V7MBRVC5Oa6wX/1ati1C8qUsXdJSq4KFYxF33fuhBkzzG7m7+PP132+5teHf+Vq8lXaft2WCWsmSF+AEHbkusF/yxZo1y5r1IoopAcegClTYED+qZp61O3Bvif38Z+Q/zBl8xRazGqRc/1gIQR/Rf/F7rO7bX4c1wz+Fy7A0aNwxx32LknJpxS89BLUqGEMAU3PO9ePn48fX/X+KusuoN3X7eQuQIgMFxMu8sB/H+CRXx6xed4s1wz+W7ca/7drZ99yOJOrV6FrV/jsM4s2z7wLeKz5Y1l3AdIXIFxZuk5nyOIhnL9+nm/v+xY3Zdvw7JrB38sLIiKgRQt7l8R5lC1rfK4vvWTcVVnAz8cvR19Au6/b8fLal+UuQLikD7Z8wMojK5nefTotqto+NslQT2E9MTHQtKnxb/16cL81fbY5cUlxjF01ltmRs2lSqQlz+s6hZbV8R6sJ4RQ2n9pMpzmd6NeoHz898BOqCH2RMtTTHK0hWWqWNhEYCB9/DJs2sW/cGwVaXSvzLmDloJXEJsXS9qu2vLL2FbkLEE7vXPw5BiwYQJB/EF/2+rJIgb8gXC/4x8QYwzt/+MHeJXFOgwdzNqI7pefO5tylawVe2u/uenez76l9PNr8Ud7Z9A5hs8KkL0A4rZS0FAYsGMDlxMssfHAhfj5+xXZs1wv+O3caE5Jq1bJ3SZyTUgxr/wS9Hp1OqvvNpLEFWV3L38ef2X1ms2LQCmKTYqUvQDitl9a8xIaTG5jVaxbNq5ieJW8rrhn83dwgJMTeJXFaB1O8iPcuhVdqChH/3Oy3KejqWvfUu+eW2cEyL0A4i/n75jN963Seaf0Mg4MHF/vxXTP4N2oEpUrZuyROK3MVrSf+WsDsBW/Q5tTeHM8XRObs4JWDVhKXFEe7r9sxfs14klKTrFpmIYrTrrO7GLpkKOE1wvmg2wd2KYNrBX+tjeAfFmbvktjd4t0xBeqQLYjM9M9ft+rLyfJV+HD5h1ROSyzS6lp317ub/U/tZ2jIUN7b/B6hX4SyNXqr1cosbPedsOV3rSQ6e+0sfeb3IaBUAAsfXIiXu5ddyuFaQz1TU2H6dGjWDHr0sO6+S5DFu2OYsGgviSlpWc/5erpbdfm9xbtjmLoqikoHIlnwwzjOdb6HwFVLrZJO47d/fmP4suFEX41mTNsxTL5zMqU8S+U4rqzqVTC2+k4Ux3fNUVjy3UtKTSJiTgR7z+9l89DNhFSxfvOzDPU0xcMDxo1z6cAPxpKH2f8YoWAdsvnJ/kdwoXEIUU+/SODq5TBrllX2361ON/Y9uY8nwp5g2p/TaP55czac2JAVaGJiEws8ysjV2eo7YevvmqOw5LuntWb4suH8FfMX3933nU0Cf0FYYw3fkuPwYWOYZ7Vq9i6JXZnreC1oh6ypmg6Qo6YXE5tI/3IdWDbwHHU6dChawbMp612Wz3p+xoNNHmTY0mFEzI3gNrfeeKUMxo2b/TmZgcbZapnWtHh3zM0lOnMp6HfC0vcXdb+OJq+LXOZ37/X1r/P9nu+ZHDGZfo362aOYObhW8B87Fk6cgH377F0Su6rm72vyj70gHbK5b+czazreHm63/BEkpGoebT6YzY0bG/0uqangaZ3F3SOCItgzcg+v/f4a0/+cgbv3n1RMGYVv+s1+HWcLNNaU+Xs0pzCd9LnfX9TvWklg7jsWE5tIrfEr8Ci7gaOpU3ks5DFe7fhqMZfONNdq9pHOXiDXerwZfD3dC9Qha66mE5uYYnL7M7GJRuAfNgyGDDF+tpLSXqX5sPuHNPP6CDd8OO89iYueH5LGVcD5Ao01mfo9Zirod8IUa3zXilNhO6fz+o4luEVyNOVDSqWH0rP6pGKbwZsf1wn+Z87Av/9K8OfW9XgD/X0L3AFX0Np0NX9fo7O3Th2YN89IA2Flk+/uR630T/FLGcB19w2c8XmKFK9NvNCtvtWP5Szy+j1ao1PWGt+14lKUPiNTFzmAG+ofLni9jaeuTsXk8UxffcwGJS8c12n22bnT+F+CP2D8URblD9Dc7Xz5Up4kpaTfMrojq6Y3YQJs3w7PPw+hodCxY6HLkFvm+UxdVYYTce2J8/2UM2oK30Ttp1Xdz6herrrVjuUszP0eA/19rRagi/pds5XcfVbXk1Pzbbc35+Z3z9ifBlJUDOe8J+Kmy1A5+Q3cKO1QTZCuU/P/+2/jfzMLjYuCMXc7P6lXk7xrem5uxqLvderAgw8auZasqG9oIJvH30XMu6OIfXUf07pNY82xNTSe2ZiZ22bafIGMkqakNMtYe66AqVq+uSbLmNhEi46Z+d07PqUnlfyuc85rIqC57cabeBAAOFYTpOvU/AcPhiZNZM1eK8ld08k9rjnPmpKfH/zyC3TrBsePG9lAbcDDzYOx7cbSt2FfRi4fyahfR/HD3h+Y1WsWTSs3tckxS5r8fo+OwNzgAsjne5aHvPo6TCnIMS8mXOS890R08lUqJ7+DpzbuOB3toupak7yEY0lKAh+fYjmU1pof9v7AmFVjiE2KpU+dJ4k+eTfn4tIdMuCJm8KnrDPbNLV5/F2F2met8SsoTOTL75hXEq/Q+dvOHLx4kAmt5rJqV0CxX1QtneTlGjX/xET48Uejplmjhr1LIzL5+LB4VzSnX3uLlKvxLOw51GZ/IEopBgcPpkfdHjw4/ykWHvkEj/QFVHB7ipjY0CLXJIXt2GKugLm+DgAFZi8MeR3zavJVevzQg/0X9rNk4BJ61O3BxG6FLqLNuUab/4ED8PjjsE0yQjqSVxfvZcxPkVQ59Q9jN/1A2Jb/2XxGbkCpAJLPj6Ry8tuAG+e9X+OC51TiUy463axTZ7B4dwxuZoZGFqX9fFz3Bni6m96vBtwLeMxryde454d72HV2F//t/1961HX8LAKuUfPfs8f4PzjYvuUQWRbvjuGHrafQSvFK96epGXuWqSun81jp8kxd5WXTGviZ2ER8aU615E+I81hAnMfPJLrvIOHaI6Sld8LdzfLlJ4XtZLb1p5lomrZK+3ke7T5pWuPr6W5+1Fo2sUmx3P3D3WyP2c78B+bTu0HvAhXDXvmoXKPmv2ePkcK5dm17l0RkmLoqKutv74aHJyP6vcpJ/2rMWvQmflH7bXrszNqbwgv/1EFUS56Jd3o9Lnt9Tpuv2tht5TBHzX5pr3KZ65R1V+qWuQIFLePUVVGkpJuP/pmj1PKbn3A58TJdvu3CzjM7WfDgAh5o/ECBztGe+ahcJ/g3aVKgBcWFbeVuO43zLcuQB98gzqcsHa6etOmxcw9v9NSB3K7f4fmwzzhz7QxtvmrDE8ue4FLCJZuWIztHTUqX2TRnj3KZa19P1/qWwF/Qzy6vtvvMGn72oZubx991S+A/F3+Ou+bexb7z+/hlwC/0bdi3YCeIfRPfuUbw37tXmnwcjKm207PlKtFt2Gc0enWM8US6bcbkm5p1OqVfMB/c+ySHRh1iTNsxfL37axp82oAvdnxBWrrlQwILyxGzX2Y1zeV6vrjKZa59PffzhfnszO3b1F2FKceuHCN8djhHLh9h6UNL6Vm/Z57bm2PPxHeuEfwPHoQ33rB3KUQ2piYXKaBfx/rGH97atdCqFZw7Z5Pjm6vVlfMux7Tu04gcGUmTyk0YuWIkrb9qzZbTW2xSjkz5JQazRzNQ9qa53IojOFk6Aa0wAdTcvqc92DzfwL/n3B7CZ4dzJekKax9dS7c6hR/SY+kFzhZcI/hXrGiziUSicEzVvqcPCOGtvs2MDby84NAh6NwZLly45f22boduWrkp64esZ9798zgXf47w2eE88ssjxFy1TQDO64/dXs1AeQXP4ghOluYFKkwALWzOoXXH19Hxm464K3c2/mcjbau3tehczH1f7TnD2vkneS1bBrt3w6uvGqkFRMmxbh307An16sHq1XDbbYDxhzRuwd+kpN387nq6K6Y+kH+trTDib8TzzsZ3+PDPD3F3c2dC+wk83+55fD2tFwBNrXhlSlEmNhWUuclVCpg+IMRh5kQU12ph3+z+hhHLR1C/Yn1WDlrJ7f63W6V81h7tY+kkL+cP/o8+agSR6Oii70sUvzVroE8f485t82aoVInQyb9xJeHWPCzlS3kyqVcTmw2bO3blGONWj2PRwUXU9KvJlM5TGNh0oNVS9GYPAub+KhVwfErh2pcLU57cQUsBD7etefMOzUHYcrhkuk5n4u8TeXvj23Sp3YUF/Rfg5+Nn8fttMUM5L8Ua/JVSPYCPAHfgK631lFyvewPfAmHAJWCA1vpEXvu0WvAPCTFW7lq5suj7EvaxZYsxQ/vjj8HNjaDxK8xuampstrVrgOtPrGfsqrHs/nc3bQLbMK3bNMJrhltt/1D8AcMcV18Ted72KEb9bxiX0zdT2e1uZt7zGQ+EBRVoH+ZSSdjqQl5s6R2UUu7ATKArEA1sV0ot1VofyLbZMOCK1rquUmog8B4woKjHzldKijG718XX7C3x7rjD+Adw7BgtYg6yK7CRyU0Lm5K3ICKCItgxYgff/f0dL697mfbftOe+hvfxbud3aRBgnbbacd0bmGwqKO7EYI6ajvkWKSmQnHwzcePGjXDpEly7BtevQ3w81K0LfTOGY44eDVevwo0bxvtSUqBrV3jmGeP1Tp04GxvL7VcOsUzfwFNXZmWDQF67cQivlDR6j33EaEbO/OfuDo88Ao89BrGxMHSosVqdpyefHrpIXKri1wZ3sLFWC/wSrzF0x1J8yvjCtEPg7W38a98eGjUyylWqlLHmuA1ZY++tgaNa62MASqn5QB8ge/DvA7ye8fMC4FOllNK2bnOKijJ+qc0c6xZVFMHo0cz732+M7TmWFY0sWxM4dw3aGrVZN+XGkJAhPND4AWZsncGUzVNYGrWUx1s8zsROE6lWtmjrRNs626bD1+i1NoJgdPTNhZh8fKB/f+P1oUONNO2XLsHly0aQv/vum3f4gwbd2tT7wAM3g/+qVUbOr8zA6+lpBO0MF5KvcOjqAdI8wEPXIZVy3PDwIjEljemrD9O7dGmjjGlpxpDk5GRjeVIwYs6RI8b/KSncmZDE9fhEDlSuxcZaUDEhjme3zDO2/S1b+b74wgj+hw8bMasEBP9A4HS2x9FAG3PbaK1TlVJxQEXgohWOb150NPj6yhh/ZzJnDvFd72Hm0ve4PfYsn7XtD0rh6aYo7e1hNif74t0xWZ1r1kwPXNqrNK90fIXhYcOZvGEyX+z8gm///pZn2zzLi+EvUt63fKFP1Va17vw+g2K5MGhtBO5//oFjx4zU3ikpMGmS8Xr37kYnf3bBwTeDv1JQpYoxebNiRahQwQicmRYtMoJnmTLGv9Kljdp0pijTcwBS0lKYsHYC0+7ei1d6HQJuTMBTV8mxzYnr6beWLZvF0TeY+vBHOT4/gN9XRaFiE0muU4/FO0/Tt0mlm3ceSUlGqnOAhg2N0W42VuQ2f6VUf6C71vrxjMePAK211s9k22Z/xjbRGY//ydjmUq59jQBGANSsWTPs5EkrzPRMSzO+KDLSx3kkJXH6voHU+N8S/le/HdMGvczTvUMBeO6nSJNv8ff1JHJSN5u3pR+7coyJv0/kx70/Us67HC/c8QLPtnmWst5li7xva8nrMzDX3FTofpOkJCPQHjoEJ0/Ciy8azw8ZAt9+m3PbunWNGjMYC/5cuADVqxt9dlWrGsG+rO0+xxOxJxi0cBB/Rv/J062eZuff93I27tbRV3l9V4pr5FFeiq3DVynVDnhda9094/EEAK31u9m2WZWxzZ9KKQ/gX6BSXs0+ks9f5ElrmDEDfv4Zfv89a12AvDqDT0zpWWydb3vO7WHi7xNZErWEir4VeTH8RZ5q9RRlvOy/mFBen4G5VMfuSuU9ASotDY4eNfJneXrC11/De+8ZNfuMmdrpKHpMWsxTfcLoG73LeK1OHeM9QUFG7dwOtNZ8+/e3PPPrMyil+LLXlzzY5MFCBXJH6Kgvznz+24F6SqlaQAwwEBiUa5ulwBDgT+ABYJ3N2/uFc1MKxowxOu7c3eHKFfj5Z5Suhlbm7/LMBTdrT1oKvi2YxQMXsy1mGyOXvMhLa15iwuq3qe7xIFO6j+OhVvZbVD6vz8DcxK40rXM2j8XEGO3ru3ZBZKSRPyshwfi/WTMoXx6aN+dQRE9mXfThgH8gx8tXIznJ3dhPvxb07dXLpudpiQvXL/DkiidZeHAhHW/vyLd9v80av59Xv4u5pjF7pmsoqCIH/4w2/FHAKoyhnrO11vuVUpOBHVrrpcDXwHdKqaPAZYwLhBBFl5msb/ZseOEF5tUO5bnuz/JvuYAcm5Uv5QkU/yiaM+cDSTw7nipp9xLrMZ9TaV8xeMVPLP1nODP7vEIF3wo2OW5e8voMpq6KynFhUDqduhdPE3L2MMH/HmHZ5Z70nfmU0dk6YoTRTh0aCsOHG+tjV8vo6O7XD/r1Y9iUdcRUyBn4bDECq6C01szbN4/Rv47mavJV3uvyHs+3e/6WdN6m+l3y6jMprsqFNVilO1lrvRJYmeu5idl+TgL6W+NYlnL40QzCusaOhbJlafncGH6b/TTvd3yUH0N6kO7mjqe7YlKvJkDxr1mbmXTMm0bcduMNklUUcZ4/Mz/qQ5bPmMXIsJE81/Y5AsvlzFJpy/Ll+RmkpzNh8X7KXLnAjOUfEHz2CGVvGMHsqlcpdmYOse3U6WYzTx6T3ByxJnwy9iSjfh3F8sPLaRPYhq97f02Tyk0sfn9eieQcZYiuJZxyMRdbLPgsHJxSMGIEHnfdRdLD/+Gt1f9H9avn+a7vU7cEz+Icu547yHnrBlS+8Rop6jgdm27mw60f8tFfHzGo2SBeuOMFjsaUL5bvbt/QQPo2r2okPdyyBT7+DrZsoW/PnvDI84yftwPflGR+aXIXkdXq83fV+hyrEEi18hnt8qVLG+31+XCkmnBSahIfbPmAdza+g1KK6d2n80zrZwq8eE9eF7TirlwUhVMG/7yuzI74SxBWVLculbf+AfPnM7J9e0bWqGHUUI8fh1q1ir045oJfkF9jBtS+j8OH7+VIwk989/dPzP17Ln6qBV5pvfAlDJWRd9Fq392kJDh92siVBMYwyUOHjJ8rVYJ27SAsLOs4D3vPKHIN1hFqwlprlkYt5YXVL3D08lHub3Q/H3b/kJp+NQu1v/wuaCVlYpxTBv+8rszSHOQClIKHHrr5+IUXjM7JkSNh3DioUaPYimIu+N3ZsFLG8/5U4An8Uh4i2Ws1l92Wkub9Bh7pVSmbdg+lU7viTpnCNZNcvGjkQ9q0yfh/506oWfPmcMpnnjHGvoeHG8MsszXfWKsGa++a8NborYxbPY5NpzbRMKAhvw3+ja51uhZpn45wQbMGpwz+5q7Mfr6e0hzkimbONNZz+Owz+L//My4ML71k1HxtzFzwy3136k45St24nzKqL9fcNnHNfQVXPL8m1uN7SqV1oLZvb7TW5pPIpacb4+m3bjXG0Lu5wfjxxpBLLy9o2RKefRY6dDCGySrF4nZ9jHJ9fZhq/qdtFpTtURPefXY3b2x4gyVRS7it9G183vNzhrUYhodb0UOevS9o1uKUWT3Njc/18XQzmQ2yuJNlCTs5eRKmT4cvvzQC42uv3ZoTppiYG2sPN5PT3VDHuOaxguvuG9AqiaaVmzIsdBiDmg2icunKRrD/6SfYupUbm7fgdTUOgEFjZvPgI93o634J4uKMRXEy5kFksiTNsL0nKxXGtphtvLPxHZZELcHfx58xbccwpu0Yk5Ps8msFKKmtBC6f0tnUL27MT5HFml1POKjLl43/K1SABQvgP/8x1g3o18/ID2PFWaTmAkh+s2yzv+e5Nn74HJ3PkTU/UfngaWbe4U6lTvfw0qVG3DF6Klfr1Gd1mSD+qtqQnYGNOFYhEF8vjyJNRnKEyUqWStfpLD+8nA+2fMDGUxuzgv7oNqPx9/E3+R5nvfhB8U7yckimbjVzj2HO5IhjcIUNVcg2tr5+faMZaMkSoxbt5WVkEF261LgIZDSRFEZeo85MtRuXc0tncl3oXCGVvuPvMmbA9ugBE45mbZNSpTLpt93BhLPb6Hx5GRVeK02qRw1IbI1veigKIydMfp3E+Q3BdMQhmrmdiz/HN5HfMGvnLI7HHqemX02md5/OsNBh+abTyG9QiCsMGnHa4G+Ks3TUCCsKDoZZs4y+gC1bjIvAgQM3m4EGDzYmNLVoAQ0aGB2jjRpZlCwwRwDRmrI3Eih/5Srfzb/MwvcGoVJSSH32WQLORVP76jkCY//FLS3N6JR+/30jn01IiLEgUVgYtGiBZ5UqPAo8nJ7GhpMbmLd3HrN3zSfdew1K++Kb1pJS6e3wTWvBmVjzZctvxIojDdHMLjk1mZVHVvL93u9ZFrWMlPQUIoIieLfzu9zf+H6L2/Sd4eJXVC4V/J2lo0bYgLu70RnaIVea6HbtjHbztWvhu++M51q1gm3bjJ87d4Zjx7jq5sWZZEhO10TVaorXp59wJjaRX759nhpx/1Iu6Tpe6UbK3yWNOsF7g+jTOgjO7TXuRFqFQ7167PC9jalX/Ng2foXx/Xx5hsnvp7ubO3fVuoursQ1YvaUniW57SHDfQoL7XyR4bATtRjm3pry3aR9danchpEpIjvHs+VWEHKmilJSaxOp/VrP40GIWHVpEbFIst5W+jVGtRzEibAQNAxoWeJ8l9eJnTS4V/KHkjMEVDmLUKOMfGIuCHDtm/J/pjjs4VaoCB/75Fy+ScNOaf7UX/7doL36+nuyo3ogDt9UizqcMl33LEetbjtigevTJfP+xY1m7ymom8rB8NNrUVVGAJ77pYfimh1Eh5SmS3Q6T6LYD37J/M37teFgLFXwrcGfQnXSo2YH2Ndtzb/PmWe83VRGyZ0VJa80/V/5h9T+r+e3Yb6z+ZzXXU65Tzrscver3YnDwYLrU7lKkkTsl6eJnK07b4StEcTHXOVq+lCdJKekWdxoWppM1r1FDgf6+LHi6EeuOr2PN8TWsO76OU3GnACjtWZqwamG0qtaKltVaEnxbMPUq1MPT3dOCM7auhJQE9p3fx1/Rf7ElegubT23m9FVjiZAg/yB61OnBfY3uIyIoAi936+W5d/XRPi5X8xfOK/cf650NK/H7oQs2/+M11w4cm5DC9AEhFgeQwrQzm2ueyHxf1bJVeTj4YR4OfhiA03Gn2Xx6M5tPbWb7me18uu1TktOSAfBy96JhQEPqV6xPvQr1qFuhLjX9alKjXA2ql6tOaa/Cp1xOTk3mzLUzRF+N5tiVYxy5fITDlw6z7/w+oi5Fka6NtM/Vy1UnvEY4HW/vSLc63ahTvo75uQ1FlF8rgLO3EkjwF07B1Mia77eeynrdlhP68mofLkgAKUw7c15DmE29r4ZfDQb6DWRgUyOxbkpaCvsv7Gfvub3sPb+X/Rf28/e/f7P40GJSM/ooMvl6+BJQKoCKpSpS1qsspTxLUdqrNB5uHigUbsqNlPQUklKTSEpNIi4pjitJV/j32kXiU3L2Prsrd2qVr0WjgEb0b9yfkCohhFULK3TKBUfmqHcQEvyFUzA1NC83Ww3Vs1b7cGH20zc0kB0nL/PD1lM5LgCWHt/T3ZOQKiGEVAnJ8Xxqeiqn4k5xOu400Vejib4azcWEi1xMvMjFhItcv3GdK0lXiL4aTZpOI12no7XG090THw8fvN29qViqIj6qGhcvpuGX5oeHroi7rkgZ92pM6duZ/mHFn2upuDlykkkJ/sIpWDoEzxZD9eydB+etvs1oeXsFq9YuPdw8qF2+NrXL1y70PsDox/BLzvmZp6XDjNXHXSL4O/J8AQn+wink1fadeztbsFb7cGH346jt0yVxvLw1m2kc+fxlVXPhFMZ1b4CvZ9552Z1tqF5JYO5i66jj5TObaWJiE9HcbKZZvDumUPtz5POX4C+cQt/QQN7t14xAf18UxjDHwW1r5nhcEvKyOBtTF2VHvgjn1UxTGI58/tLsI5yGozZ9uLKSNqve2s00jnz+EvyFEDZVki7Ktkjr4KjnL80+QgiRwZGbaaxNav5CCJHBkZtprE2Cv3AajjqTUpQsjtpMY20S/IVTcOSZlK5MLgQtiUcAAAcLSURBVMiOS9r8hVOw9hA9UXTWHjMvrEuCv3AKjjyT0lXJBdmxSfAXTsGRZ1K6KrkgOzYJ/sIpOOIQvcW7Ywifso5a41cQPmWdyzV3yAXZsUmHr3AKjjZEryR3QFurk7YwKaqlg7j4yDKOQthAYZZkdAS5L1qQ99KTluzP0mBu7WO7KlnGUQg7Kqnt3dbOP1+QMfOOnPveGUnwF8IGrJUjpribQex50SqpF8ySSoK/sJuS1r5bkPJaY2lHe/Qb2CKxWUk4tiuS0T7CLkraBKCCltfU+gIFbbu2xzh5e46acsQRW85Mav7CLkpa+25hylvUHDH2aAax56gpRxux5ewk+Au7KGntu/Yor72aQeyZ2MxVkqo5giI1+yilKiilViuljmT8X97ENiFKqT+VUvuVUnuUUgOKckzhHEraBCB7lFeaQYQtFbXNfzywVmtdD1ib8Ti3BOBRrXUToAcwQynlX8TjihKupAU2e5TXGv0GQphT1GafPkBExs9zgfXAS9k30FofzvbzGaXUeaASEFvEY4sSrKS179qrvNIMImylSDN8lVKxWmv/bI+vaK1vafrJ9nprjItEE611uonXRwAjAGrWrBl28uTJQpdNCCFckdVm+Cql1gBVTLz0SgELVBX4DhhiKvADaK1nAbPASO9QkP0LIYpfSZurIW7KN/hrrbuYe00pdU4pVVVrfTYjuJ83s105YAXwqtZ6a6FLK4RwGCU5eZ0oeofvUmBIxs9DgCW5N1BKeQG/AN9qrf9bxOMJIRyELNZSshU1+E8BuiqljgBdMx6jlGqplPoqY5sHgY7AY0qpyIx/IUU8rhDCzkraXA2RU5FG+2itLwGdTTy/A3g84+fvge+LchwhhOORXDwlm+T2EUIUSkmbqyFykvQOQohCKWlzNUROEvyFEIUmk9BKLmn2EUIIFyTBXwghXJAEfyGEcEES/IUQwgVJ8BdCCBckwV8IIVyQBH8hhHBBEvyFEMIFFWkxF1tSSl0ArLGaSwBw0Qr7KUlc7ZzlfJ2fq51zUc73dq11pfw2ctjgby1KqR2WrGrjTFztnOV8nZ+rnXNxnK80+wghhAuS4C+EEC7IFYL/LHsXwA5c7ZzlfJ2fq52zzc/X6dv8hRBC3MoVav5CCCFycZrgr5TqoZSKUkodVUqNN/G6t1Lqp4zX/1JKBRV/Ka3HgvMdq5Q6oJTao5Raq5S63R7ltKb8zjnbdg8opbRSqkSPDrHkfJVSD2b8nvcrpX4s7jJakwXf6ZpKqd+VUrszvtf32KOc1qKUmq2UOq+U2mfmdaWU+jjj89ijlGph1QJorUv8P8Ad+AeoDXgBfwONc23zFPB5xs8DgZ/sXW4bn++dQKmMn58syedr6TlnbFcW+APYCrS0d7lt/DuuB+wGymc8rmzvctv4fGcBT2b83Bg4Ye9yF/GcOwItgH1mXr8H+BVQQFvgL2se31lq/q2Bo1rrY1rrG8B8oE+ubfoAczN+XgB0VkqpYiyjNeV7vlrr37XWCRkPtwLVi7mM1mbJ7xjgTeB9IKk4C2cDlpzvcGCm1voKgNb6fDGX0ZosOV8NlMv42Q84U4zlszqt9R/A5Tw26QN8qw1bAX+lVFVrHd9Zgn8gcDrb4+iM50xuo7VOBeKAisVSOuuz5HyzG4ZRgyjJ8j1npVQoUENrvbw4C2YjlvyO6wP1lVKblVJblVI9iq101mfJ+b4ODFZKRQMrgWeKp2h2U9C/8wJxljV8TdXgcw9jsmSbksLic1FKDQZaAp1sWiLby/OclVJuwHTgseIqkI1Z8jv2wGj6icC4s9uolGqqtY61cdlswZLzfQiYo7WeppRqB3yXcb7pti+eXdg0ZjlLzT8aqJHtcXVuvSXM2kYp5YFx25jXLZcjs+R8UUp1AV4Bemutk4upbLaS3zmXBZoC65VSJzDaSJeW4E5fS7/TS7TWKVrr40AUxsWgJLLkfIcBPwNorf8EfDBy4Dgri/7OC8tZgv92oJ5SqpZSygujQ3dprm2WAkMyfn4AWKczelVKoHzPN6MJ5AuMwF+S24Iz5XnOWus4rXWA1jpIax2E0c/RW2u9wz7FLTJLvtOLMTr2UUoFYDQDHSvWUlqPJed7CugMoJRqhBH8LxRrKYvXUuDRjFE/bYE4rfVZa+3cKZp9tNapSqlRwCqMUQOztdb7lVKTgR1a66XA1xi3iUcxavwD7VfiorHwfKcCZYD/ZvRrn9Ja97ZboYvIwnN2Ghae7yqgm1LqAJAGjNNaX7JfqQvPwvN9HvhSKTUGo/njsRJcgUMpNQ+jyS4gox9jEuAJoLX+HKNf4x7gKJAA/Meqxy/Bn50QQohCcpZmHyGEEAUgwV8IIVyQBH8hhHBBEvyFEMIFSfAXQggXJMFfCCFckAR/IYRwQRL8hRDCBf0//p2Tvqqt8iYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.scatter(x, y, label='Samples')\n", "\n", "Xreg = generatePolynomialRegressionMatrix(xlin,d)\n", "\n", "plt.plot(xlin, Xreg @ beta, 'g-', label = 'Regression line')\n", "plt.plot(xlin, y_func(xlin), 'r--', label='Population line')\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 3.82794677e-02, 8.94074990e+00, -5.46596826e+01, 8.98472825e+01,\n", " -7.81583544e+01, -2.18941451e+01, 9.11085904e+00, 1.75846110e+00])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beta" ] }, { "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.6" } }, "nbformat": 4, "nbformat_minor": 2 }