{ "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" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XucTfX++PHXx5ib3BnFoJlyjdHIkJKakstJIXFSKRVHKpwoJ/3qnFTOSTknuuiUvqfoRh1p6KoipyIxYyRkSiIzJGRIM8NcPr8/1szYM/Z9r73X2nu/n4+Hh31Zs9dnzV7zXp/1/tyU1hohhBDRpY7VBRBCCBF6EvyFECIKSfAXQogoJMFfCCGikAR/IYSIQhL8hRAiCknwF0KIKCTBXwghopAEfyGEiEJ1rS6AK82bN9cpKSlWF0MIIcJKTk7OQa11kqftbBv8U1JSyM7OtroYQggRVpRSu73ZTtI+QggRhST4CyFEFJLgL4QQUci2OX8hROiUlpaSn59PSUmJ1UURXkpISKB169bExsb69fMS/IUQ5Ofn06BBA1JSUlBKWV0c4YHWmkOHDpGfn09qaqpfnyFpHyEEJSUlNGvWTAJ/mFBK0axZs4Du1CT4CyEAJPCHmUC/Lwn+QggRhST4WyQrt4A+s1aROv09+sxaRVZugdVFEsJSMTExpKen07VrV0aOHElRUZGpn79gwQImTpzodpvVq1ezdu3a6ufPPfccL7/8csD73rVrF127dgUgOzubyZMnB/yZgZLgb4Gs3ALuW/oNBYXFaKCgsJj7ln4jFwAR1RITE9m0aRNbtmwhLi6O5557LuRlqB38J0yYwE033WTqPjIyMnjqqadM/Ux/SPC3wOwVeRSXltd4rbi0nNkr8iwqkRD20rdvX3bs2AHAE088QdeuXenatStz584FjJp0p06dGDNmDN26dWPEiBHVdwopKSkcPHgQMGrZmZmZp3z+O++8w/nnn0/37t25/PLL2b9/P7t27eK5555jzpw5pKen8/nnnzNjxgz++c9/ArBp0yZ69+5Nt27duPrqqzl8+DAAmZmZ3HvvvfTq1YsOHTrw+eefuz221atXc+WVVwIwY8YMbr31VjIzMznrrLNqXBReffVVevXqRXp6Orfddhvl5eWuPtIvEvwtsLew2KfXhQi5zMxT/z37rPFeUZHz9xcsMN4/ePDU93xQVlbGBx98QFpaGjk5Obz00kt89dVXrFu3jhdeeIHc3FwA8vLyGD9+PJs3b6Zhw4Y8W1U+L1x00UWsW7eO3NxcRo0axeOPP05KSgoTJkxgypQpbNq0ib59+9b4mZtuuonHHnuMzZs3k5aWxkMPPVSjzOvXr2fu3Lk1XvfG9u3bWbFiBevXr+ehhx6itLSUb7/9ljfeeIM1a9awadMmYmJieO2113z6XE8k+FugVeNEn14XIhoUFxeTnp5ORkYGbdu2ZezYsXzxxRdcffXVnHbaadSvX5/hw4dX16zbtGlDnz59ABg9ejRffPGF1/vKz89n4MCBpKWlMXv2bLZu3ep2+yNHjlBYWMgll1wCwJgxY/jss8+q3x8+fDgAPXr0YNeuXb4cNoMHDyY+Pp7mzZvTokUL9u/fz8qVK8nJyaFnz56kp6ezcuVKdu7c6dPneiKDvCwwbWBH7lv6TY3UT2JsDNMGdrSwVEI4WL3a9Xv16rl/v3lz9++7UJXzd6S1drl97a6OVc/r1q1LRUUFgMt+8JMmTWLq1KkMGTKE1atXM2PGDJ/L6yg+Ph4wGq3Lysr8+lnHn9daM2bMGB599NGAyuWO1PwtMKx7Mo8OTyO5cSIKSG6cyKPD0xjWPdnqoglhKxdffDFZWVkUFRXx+++/8/bbb1enY3766Se+/PJLABYtWsRFF10EGDn/nJwcAN566y2nn3vkyBGSk42/t4ULF1a/3qBBA3777bdTtm/UqBFNmjSpvut45ZVXqu8CgqFfv34sWbKEX375BYBff/2V3bu9mqnZa1Lzt8iw7skS7IXw4LzzzuPmm2+mV69eAIwbN47u3buza9cuOnfuzMKFC7ntttto3749t99+OwAPPvggY8eO5R//+Afnn3++08+dMWMGI0eOJDk5md69e/Pjjz8CcNVVVzFixAiWLVvG008/XeNnFi5cyIQJEygqKuKss87ipZdeCtpxn3POOcycOZMBAwZQUVFBbGws8+bN48wzzzRtH8rdbZWVMjIytCzmIkRofPvtt3Tu3NnqYnht165dXHnllWzZssXqoljK2femlMrRWmd4+llJ+wghRBQyJfgrpQYppfKUUjuUUtNdbPNHpdQ2pdRWpdTrZuxXCBGdUlJSor7WH6iAc/5KqRhgHtAfyAc2KKWWa623OWzTHrgP6KO1PqyUahHofoUQQvjPjJp/L2CH1nqn1voEsBgYWmubPwHztNaHAbTWv5iwXyGEEH4yI/gnA3scnudXvuaoA9BBKbVGKbVOKTXIhP0KIYTwkxldPZ1NKl27C1FdoD2QCbQGPldKddVaF9b4IKXGA+MB2rZta0LRhBBCOGNGzT8faOPwvDWw18k2y7TWpVrrH4E8jItBDVrr+VrrDK11RlJSkglFE0KEk7///e906dKFbt26kZ6ezldffRW0fWVmZhLN3cnNqPlvANorpVKBAmAUcH2tbbKA64AFSqnmGGkgcyeqEEKEtS+//JJ3332XjRs3Eh8fz8GDBzlx4oTVxYpYAdf8tdZlwERgBfAt8KbWeqtS6mGl1JDKzVYAh5RS24BPgWla60OB7lsIYY1gLEa0b98+mjdvXj3XTfPmzWnVqhUPP/wwPXv2pGvXrowfP756vp/MzEymTJnCxRdfTOfOndmwYQPDhw+nffv2PPDAA4D7qZ8dffTRR1xwwQWcd955jBw5kmPHjgEwffp0zjnnHLp168Y999wT8DHaiSn9/LXW72utO2itz9Za/73ytb9prZdXPtZa66la63O01mla68Vm7FcIEXrBWoxowIAB7Nmzhw4dOnDHHXfwv//9D4CJEyeyYcMGtmzZQnFxMe+++271z8TFxfHZZ58xYcIEhg4dyrx589iyZQsLFizg0CGjfulp6ueDBw8yc+ZMPvnkEzZu3EhGRgZPPPEEv/76K2+//TZbt25l8+bN1ReUSCEjfIUQPgnWYkT169cnJyeH+fPnk5SUxLXXXsuCBQv49NNPOf/880lLS2PVqlU1pl8eMsRILqSlpdGlSxdatmxJfHw8Z511Fnv2GJ0QPU39vG7dOrZt20afPn1IT09n4cKF7N69m4YNG5KQkMC4ceNYunQp9erVC+j47EYmdhNC+CSYixHFxMSQmZlJZmYmaWlpPP/882zevJns7GzatGnDjBkzakzTXJUiqlOnTo2pkevUqVM9tbKrqZ+raK3p378/ixYtOqU869evZ+XKlSxevJhnnnmGVatWBXyMdiE1fyGET4K1GFFeXh7ff/999fNNmzbRsaOxxkXz5s05duwYS5Ys8flzXU39XKV3796sWbOmetnIoqIivvvuO44dO8aRI0e44oormDt37ilrDYQ7qfkHQVZuAbNX5LG3sJhWjROZNrCjTN8sIkawFiM6duwYkyZNorCwkLp169KuXTvmz59P48aNSUtLIyUlhZ49e/r8ua6mfq6SlJTEggULuO666zh+/DgAM2fOpEGDBgwdOpSSkhK01syZMyeg47MbmdLZZFWNYbX/MGSxFmFnvk7pHC4VnEif+jmQKZ2l5m8yd41hdvzjEMIfshhR+JOcv8mC2RgmhPCNTP3smgR/kwWrMUyIYLNrClg4F+j3JcHfZNMGdiQxNqbGa2Y0hgkRTAkJCRw6dEguAGFCa82hQ4dISEjw+zMk52+yqjxoODSGCVGldevW5Ofnc+DAAauLIryUkJBA69at/f55Cf5CCGJjY0lNTbW6GCKEJPibrHZXz6p5TwCp/QshbEOCv8ns0tUzXPphCyGsIcHfZHbo6il3H0IIT6S3j0mq5jd31VcilF09gzXrohAickjN3wTOpnRwFFtHUXSijNTp74UkBWOHuw8hhL1J8DeBs5p2lcaJsfx+oozDRaWA7ykYf3L3rRonUuAk0MtAMyFEFUn7mMBVjVoBp8XXpbS8ZjLI2xSMvysmyUAzIYQnEvxN4G5Kh0BSMP7m7od1T+bR4WkkN05EAcmNE2VWUSFEDZL2MYG7+c1nr8jzOwUTyIVDZl0UQrgjNX8TuKtpB5KCkUnihBDBIjV/k7iqaQcy10+wVkwSQghZyctijr15GteLRWs4UlxafZEAmSROCOE9b1fykrSPhWr35jlcVEphcWmNnj0Aa6Zfxpxr0wGY8sYm+sxa5bHHjxBCuCNpHwu5Gx8ANXv2yHQNQggzSc3fQs56AdW2t7BYpmsQQphOgr+FYpTyuE2gYwWEEMIZCf6hcuIEfPghHD1qPF+1ipcWP8A/PnyaP2z/goYlx075kaqePdLlUwhhNgn+wXb4MNx5J7RoAX/4A7zzjvH6iRM0Lyviym8/59/LZpH71PW89N8HSTp2GDBvrIAQQjgjDb7B9N57MH487N8PN9wAI0fC5Zcb7w0axHfLV3Ltf3Pp8NO3ZO7Moe+ujRxv2Ii5I9MZ1u0MiDECvqwLLIQwm/TzD5aKCujTB44dgwULoEcPp5vVmLWzUQLTBnViWEo96NsXpk2DMWNCW24hRFjztp+/1PzNVl4OZWUQHw9Ll0LTpsZjF5yODN6/H5o3h5tvhs8+g6efhnr1gltuIURUkZy/mbSG22+HwYOhtBRatnQb+F06/XT45BO4/3548UW46CL45RfzyyuEiFoS/M10993wwgtw/vkQGxvYZ9WtCzNnwrvvwvbtMHWqOWUUQggk7WOet96COXNg0iQjaPvA7WpdgwfDqlXQvn0QCi2EiFam1PyVUoOUUnlKqR1KqeluthuhlNJKKY+NEVarWpA9dfp7nufS2bcPbrsNMjLgX/8CLwZvOe7H42pdvXtDs2Zw/Djceivk5/t/YEIIgQnBXykVA8wD/gCcA1ynlDrHyXYNgMnAV4HuM9h8Xj7x4EFIToZXXvE53ePT1A27dsGSJcZ4gcJCn/YjhBCOzKj59wJ2aK13aq1PAIuBoU62ewR4HCgxYZ9B5fNcOmlpsGkTdOrk8758mrqhY0d4+23Iy4Nhw4w7ASGE8IMZwT8Z2OPwPL/ytWpKqe5AG631uybsL+i8Dsj79hmNvEVFPqV6HPk8dUO/fsa4gf/9zxgDYNNxGkIIezMj+DuLetURSSlVB5gD3O3xg5Qar5TKVkplHzhwwISi+cfrgPzII/DUU8ZFwE9+Td1w/fXw2GOwcqWRChJCCB+ZEfzzgTYOz1sDex2eNwC6AquVUruA3sByZ42+Wuv5WusMrXVGUlKSCUXzz6Wdkk65op0SkHfsMLp1jh8PZ5/t974c1/8FY6bPqhST20bmadNg61ZITfV730KI6GVG8N8AtFdKpSql4oBRwPKqN7XWR7TWzbXWKVrrFGAdMERrbcu5G7JyC3grpwDHZIoCrulRayTu3/4GcXHwwAMB79Nx8rbyyjSOx0ZmpYzJ4ioqjLuAH38MuBxCiOgRcPDXWpcBE4EVwLfAm1rrrUqph5VSQwL9/FCp6tp51xubTmns1cCn2x3SUJs2waJFcNddxiheE/i9YMvevTBrFlxzDZTYvi1dCGETpvTz11q/r7XuoLU+W2v998rX/qa1Xu5k20y71fodu3a6UqOx97TTjLz7tGmmlcGfBVuycgvo8+p3jO03GXJz+eHWiaaVRwgR2WR6BzyvpQu1Gnvbt4fXXoPGjX0bDObt53vxuuMFa2W7XrzU4yrOXvQfvnz6Zb/2L4SILhL88bwcYo3G3gULYMsWwI/BYG742uun9gVrVuYtfJuUQvv7p0r6RwjhkQR/3C+H6LiiFr/8Ykzj8NxzQAB5eicce/2o2vt1ovYF63jdOCYOuZfxQ6dDQoLP+xdCRBeZ2A2j1n3f0m9qBPLE2JhTg+/8+cZavBON3LrZC6s7ndvfhVaNE09po/iheZvqLqP8/DOccYZf5RBCRD6p+eNlrbu0FP79b+jfv3oaBysXVnebJpo/H9q1M8YiCCGEE1Lzr+Sx1r10qdGt8vnnq19ydccQioXV3a7rm3QF/OUvcMstsHp19VrAQghRRYK/t37+Gbp1M2bUrGT1wuouL1itW8OTTxrLQD75pCwEI4Q4hSzg7ouKCqgTJpkyrWHoUPjoI9i8GTp0sLpEQogQ8HYB9zCJZBbbv98IpuES+MGY/uH55yEpyZgDSAghHEjax5OKCmNN3gEDjIZUD9wuyRhqLVsajb7+LCIvhIhoYVSVtcjatbB7N1x8scdNzRz0ZZr4eOOu5bXXoMDCcgghbEWCvyevvw6JicbKWR6YOejLVAUFMG4c3HmnLP4ihAAk+LtXWgpvvmk0nNav73Fzswd9maZ1a3j4YVi2zFgGUggR9STn785HH8GhQ3DDDV5t7mzUbdXrweaxreGuu4zUz+TJxkC1Bg2Cty8hhO1Jzd+dSy+FxYuNxl4v+LUkowm8amuIjTV6/+zdC3/9a3D3JYSwPQn+7tSrB9dea6zY5QVfJ2czi9dtDeefDw8+CJdcEvx9CSFsTdI+rvzvf0ZPn8mTjcVbvOTL5GxmcdWmUFBYTFZuQc3yPPhgUPZlebuGEMInUvN3ZcECePzxsOgj765NwWlKprwc/v53r8YteLuvULRrCCHMI8HfmYoKeO89Yx6fuva/OXLW1lDFaUqmTh3jzmbaNGPOogD3FarJ7IQQ5pHg78z69XDgAFx1ldUl8UpVW4Mrp6RklIJnnjFW/LrnHr/2Fep2DSGEuexfrQ2A310S33nHmAZ50KDgF9Ikw7onM3tFnvddTTt0gHvvhUcegbFjjZ5NPuxLgr0Q4S1ia/4BdUksLIR+/aBJk6CX00w+p2Tuuw9SU2HSJCPVJYSIGhFb83fXJdFjrXXevLAMhj6vL5CYCAsXGv+H04ylQoiARWzw97tLYnm5kfIJ02Doc0qmb9+Tj6uOHRnFK0SkC88I5wW/uyQOHQqjRwehRDY3eTKMGgXIKF4hokHEBn+/uiQeOwYffwwtWgS5dDZ0+umwZAmsWCGjeIWIAhEb/P3qkrh6NZw4AYMHh6qY9nHPPdCuHUyaxMGDR51uIqN4hYgcEZvzBz/y36tWGSN6+/QJXqHsKj7e6Ps/aBBTv3mHR7sPP2UTGcUrROSI2Jq/X1atMgJ/QoLVJbHGwIFwzTXcvGEZjSmr8ZaM4hUiskR0zd8nWhuDnVq2tLok1nrqKeIrKphxQElvHyEimNI2XdYvIyNDZ2dnW12M6KU17N8PZ5zh1ebSNVQIe1BK5WitMzxtJ2mfKhs2wL59VpfCPm6/3UiBlZR43FS6hgoRfiT4V7nxRvjTn6wuhX2MHAk7d8Jjj3ncVLqGChF+oi7n7zQ90QLIy4Px460unn306wejRlH+j0e57tjZbIhp4jKdIwu8CBF+oqrm7yo9kf3SW8YGl11mafns5sNb7qGYOty+dC5aa5fpHFngRYjwE1XB31V64uesD6BpU+jWzaKS2dMjG4/wxEWj6fTLLpJ+Pww4T+fIAi/RLSu3gD6zVpE6/T36zFolbT1hwpS0j1JqEPAkEAP8n9Z6Vq33pwLjgDLgAHCr1nq3Gfv2has0RPr3G/m0bReOfL1Peqg42FtYzMIeV/JGt/78Hl+vxuuOfJ5N1GTS08g6VXfTVZWqqrtDQL4Dmws4+CulYoB5QH8gH9iglFqutd7msFkukKG1LlJK3Q48Dlwb6L591apxotPFTkZd9w/iy0rZKydtDVW/r9/j6xFXVkpG/lbWpqQ7TedYtcCLBB9rBTR1urCUGWmfXsAOrfVOrfUJYDEw1HEDrfWnWuuiyqfrgNYm7Ndnrta6zW98Bj80byM9VGpx/H1NWruYl9/8G90O/2SrdI70NLKWNPaHLzOCfzKwx+F5fuVrrowFPjBhvz5znOytysjNHzN8y8rq53LSnuT4+3qx51B+S6zPwx88w9TFG22T25XgYy1p7A9fZgR/5eQ1p8OGlVKjgQxgtov3xyulspVS2QcOHDChaKca1j2ZNdMvq74AjN2QxZBtn1W/LydtTVW/rwdvvpjH+40jfc82Rm7+2DYDuST4WEsa+8OXGcE/H2jj8Lw1sLf2Rkqpy4H7gSFa6+POPkhrPV9rnaG1zkhKSjKhaK5NG9iRFuXFdDj4EznJnQA5ad2ZvSKPRZ0v5as2Xblv9Us0LTpii/SKBB9r+TV1urAFM3r7bADaK6VSgQJgFHC94wZKqe7A88AgrfUvJuwzYMO6J9Ni7XHqoNmY3Jlk6SXi1t7CYlCK+wfcwez3n6RRyTF+rdfI8vSK1T2NhHWN/SIwAQd/rXWZUmoisAKjq+eLWuutSqmHgWyt9XKMNE994L9KKYCftNZDAt13oC488D3UqcNrL0yGBg2sLo6tVfX82dG8LVff+E8wvkdbpFck+AjhO1P6+Wut3wfer/Xa3xweX27GfkxXUGAM7JLA79G0gR1PdqlUikbFvzFpw1uc/s9HrS6aEMIPUTe3Tw0vvGAs2yg8qp1euah4H+O+XALvtIML5QIgRLiJ7uAPEBdndQnCRs30ymAo3wKzZ8OoUXDuuZaWTQjhm6ia26eGV1+FK66A336zuiTha/ZsaNbMWAGtrMzz9kII24je4P/xx7BxI9Svb3VJwlfTpsai7zk5MHeu1aURQvggeoP/2rVwwQXVvVaEn0aMgFmz4NqQT9UkhAhAdAb/Awdgxw648EKrSxL+lIJ774U2bYx1fysqrC6REMIL0Rn8160z/r/gAmvLEUmOHoX+/eHZZ60uiRDCC9EZ/OPiIDMTzjvP6pJEjgYNjN/rvfcad1VCCFuLzuA/cCB8+inUq+d5W+EdpYxxE3FxcMstUF7u+WeEEJaJvn7+WhsDu+LjrS5J5ElOhqeegptuYsu0h7itRabMtyOETUVfzb+gwOje+dprVpckMo0ezb7MgZy28EX2H/oNDbaZ/lkIcVL01fxzcowBSampVpckMinF2Itu46f0MspiTp5esrSfEPYSfTX/nByoUwfS060uScT6tjSOY5Xr/mb+kF39utXTPwshTorO4N+5szT2BlHVNM+3fbWEF5c8xPk/fVPjdSGE9aIr+GttBP8ePawuieWycgvoM2sVqdPfM3093qrVtf7Tcxi7m5zBE+8+QYvyYlldy+aCdU4E81wT/ouunH95Odx9N6SlWV0SS2XlFpycm5+TDbKAKTl5x+mfp1x5D0tem8bbXy8k+fFrAv5sT7JyC2RVLz8E65wI9rlmJ+F27kVXzb9uXZg2DQYNsroklpq9Iq/6j7GKmevxOv4RHDgnnbw7/0Lyx+/C/PmmfL67/d639BsKCoull5GPgnVOBPtcs4twPPeiq+b/3XdGN89WrawuiaVcNbz62iDrrKYDnFLTG9mwL++M2s/ZffsGVnAP3AUaO9fArJaVW0CBSeeEtz8faY3/4XjuRVfwnzoVdu2CLVusLomlqtbjdfa6t1zdzsfXrXPKH0FRmeamc0ez5pxzjHaXsjKIjQ3sIJyIlkBjpqrv0ZVAG+nNONfCgatzrKCwmNTp79kyDRRdaR9p7AVONsg6SoyN8alB1lVNp7C41On2ewuLjcA/diyMGWM8NpmrgBJpgcZMzr7HKr6eE86Yca6Fkr+N0+7OMbumgaIn+O/dCz//LMEfo6Ht0eFpJDdORAHJjRN5dHiaT7USX2vTrRonGvP/nH02LFpkTANhsnALNHbg7nv09ZxwxoxzLVQCyds7O/dqs1tbR/SkfXJyjP8l+AO11+P1navb+Sb1YikprahRm6wRgO+7DzZsMHpdde8OF1/sdxlqq73IvB1vte3G1feY3DjRtN9boOdasNRus/r9eJnfefva556r+1o7pSCVDsLttxkyMjJ0dna25w29NXMm/PWvxpq9snRjwGrn/MEI8o8ON7rRug3AR45Ar17G/zk5xoRwwhLuvkc7BWyzu1E6O253kn3cZ59Zq1xeVNdMv8ynsvpKKZWjtc7wtF301PxHj4YuXSTwm8RTLdvtH0mjRvD22zBgAPz4owR/C4XD3VIwxgq4a+twxtd9ThvY0elF1U4pyOip+Qv7KSmBhARLdh1uA3KiWTBq0anT33OZmnHHl31adY5Jzd9RcTG8/rpR02zTxurSiCoJCWRtzGfPX2dSevQYbw2+NSR/INE06jQSBKMLr6u2DgAFpuTs7drWUSU6evts2wbjxsH69VaXRDh4IOsbpryxiTN++oGpX7xGj7UfhqQ7XLSMOo0EWbkF1FHK6XuBdOGdNrAjsTHOP1cDMUHYp91ER81/82bj/27drC2HqJaVW8Br635CK8X9A++kbeE+Zr8/h5tPa8LsFXFBrTHJYLDwUHWHVu4kNW1K/txN3qdcaxJjY0KSs7cqPRQdNf/Nm40pnM86y+qSiEqzV+RV/+2dqBvL+OEPsLtxK+YvfYRGeVuDum+7Dgaz6+yXVpXLVaNsjFKn9EbytYyzV+RRWuE6+leNRwj2+AQr5wSKnpp/ly4Q434Qhgid2rXsI4kNGPPHh/jva/fS9+juoO7bjj0x7NoO8UDWN8YdWuXzUJbL1Z1YhdanBH5ff3fu7vKqzoVQ5OytnBMoOmr+33wjKR+bcVbL3tcwiQFjn6XzA1OMFyoqgrJvO446tWM7RHVqrtbroSqXt3do/vzuXH22s7uKYLIyBRkdwf/bb+Ghh6wuhXDgbDi8AoZf3MH4w1u5Enr2hP37g7L/Yd2TWTP9Mn6cNZg10y+zvFeGp4nBrEgDOabmagtFcPJ2ug5/Aqirz/7XH88N6blgZQoyOtI+zZpZXQJRi8fBRXFxsH079OsHn34KSUk1fj7S+um763romAuG0KWB3AXPUAQnbweg+TNzaKgHt7k6X61MQUb+IK933oHcXHjgAWPhdhE+Vq2CwYOhfXv4+GM4/XTA+EOatuRrSstPnruxMYrZI0JbazOTt9MNhGJ6gCquBlcpYM616bb5Xdt9igpP5TO7IuPtIK/ID/433WQEkfz8wD9LhN4nn8Dgi2Z1AAASvElEQVTQocYUEGvWQFIS3R/+iMNFp04d3aReLA9e1SVs7wgcg4Crv0oF/DhrcMjKUztoKeCG3m2ZOcxeS6Ha+U4w1PP8hHSEr1JqEPAkEAP8n9Z6Vq3344GXgR7AIeBarfUuM/bt0ebN0tgbzi6/3Kj1v/56dfrOWeCvet2OPWa85di7xFXACGV31HCY96dKsHrmmHFRseu4koCDv1IqBpgH9AfygQ1KqeVa620Om40FDmut2ymlRgGPAdcGum+PSkuN0b1RvmZv2LvwQuMfwM6dnFfwLRuTOzvdNNyW0nPFLt1R7T5FQbXSUjh+/OTEjZ9/DocOGbP4/v47HDsG7drBsGHG+5Mnw9GjcOKE8XOlpdC/P0yaZLx/ySUcPFpC24PHeKqiAoXm3U4Xc9/vw6lzvIQhU2800shV/2Ji4MYb4eabobAQbr3VWK0uNpZnth/kSJnig44X8nnqeTQq/o1bs5eTUD8R/rUd4uONfxddBJ07G+WqV89YczyIzPj0XsAOrfVOAKXUYmAo4Bj8hwIzKh8vAZ5RSikd7JxTXp7xpabZ6xZVBGDyZBZ9+BFTB0/lvc7erQlcuwZt5xRBlWDXum3/O9DaCIL5+ScXYkpIgJEjjfdvvRW+/toI8L/+agT5P/wB3n/feP/6609N9Y4YcTL4r1hhzPlVFXhjY42gXSUmht1HT1AcE0tFXaOt8HjdWIpLy5nz8XcMOe00o4zl5UaX5OPHjeVJwYg5339v/F9ayqVFJfx+rJhtLVL5PBWaFR3hz2sXGdt+5FC+5583gv933xkxKwyCfzKwx+F5PnC+q2201mVKqSNAM+CgCft3LT8fEhMl7RNJFizgWP8rmLf8Mc4s3MezvUeCUsTWUZwWX9flMpJZuQXVjWvhkhoKZirD3e8gJBcGrY3A/cMPsHOnMbV3aSk8+KDx/sCBRrrPUbduJ4O/UnDGGcbgzWbNoGlTI3BWWbrUCJ716xv/TjvNqE1XyfMwTmHVKka4mPlz1+8Vp5bNQVb+CWbf8GSN3x/ApyvyUIXFHD+7PVk5exjWJenknUdJiTHVOUCnTkZvtyALuMFXKTUSGKi1Hlf5/Eagl9Z6ksM2Wyu3ya98/kPlNodqfdZ4YDxA27Zte+zebcJIz/Jy40SRnj6Ro6SEPVePos2Hy/iwwwX86/r/x51DugNw1xubnP5I48RYNj04wNJFNuzC3e/AVbrJ754zJSVGoN2+HXbvhr/8xXh9zBh4+eWa27ZrZ9SYARYuhAMHoHVraNUKWrY0gn2DBr6XwU/+nCt26HkUygbffMBxnuTWwF4X2+QrpeoCjYBfa3+Q1no+MB+M3j4mlE2mdIhECQm0ef9tmDuXQW++yaB7L69eF8BV8K+6I7Br41soufsduBote/ebXwNu7o7Ky2HHDmP+rNhY+M9/4LHHjJp95UjtChSDijpxx9AeDBsxwljG8+yzjZ9JSTFq51XGjAn4OAPlT7uLldM1+MqM4L8BaK+USgUKgFHA9bW2WQ6MAb4ERgCrgp7vF5FNKZgyxWi4i4mBw4fhzTdRuhVaub7L82dAUKDsll939ztwdWEo17pmeqygwMivb9wImzYZveqKioz/09KgSRM491y2Zw5m/sEEtjVO5scmrTheEmN8zvDzGHbVVUE9zkC5a3dx9Z2GU+Ui4OBfmcOfCKzA6Or5otZ6q1LqYSBba70c+A/wilJqB0aNf1Sg+xUCOHln9+KLcM89LDqrO3cN/DM/N2xeY7Mm9WKB0PeisWMbg7vfwewVeTUuDEpX0O7gHtL3fUe3n7/nnV8HM2zeHUZj6/jxRp66e3f405/g3HONFA3A8OEwfDhjZ62ioGnNwGfXmrAzztpd3H2nVlQu/GVKc7LW+n3g/Vqv/c3hcQkw0ox9ectutS0RZFOnQoMGZNw1hY9evJPHL76J19MHUVEnhtgYxYNXdQFC33fdnzRAsM9dt7+Digruy9pK/cMHmPvuP+m273sanDCC2dG4euRUdbG95JKTaR4XC59AZKbZ3H2ndumi642InNvHjrUtEWRKwfjx1L3sMkpuuIWZH/+b1kd/4ZVhd5wSPEPZd93X4Beqc3dY92SGndvSmPRw7Vp46hVYu5ZhgwfDjXczfVE2iaXHebvLZWxq1YGvW3ZgZ9NkWjWpzMufdpqRr/cgnGrC3nL3nYbTwLiIDP7h1OgiTNauHS3WfQaLFzPhoouY0KaNUUP98UdITQ15cdwFP2c1/KCeuyUlsGePMVcSGN0kt283HiclwQUXQI8e1fu5IX5uwDXYcKoJe8vTBS1cBsZFZPB3d2WWdFAUUAquu+7k83vuMRonJ0yAadOgTRvXP2syV8Hv0k5JTmv4riZ28ytNcvCgMR/SF18Y/+fkQNu2J7tTTppk9H3v08foZumQvjGrBhtONWFvRcoFLSKDv6src6PEWEkHRaN584z1HJ59Fv79b+PCcO+9Rs03yFwFP1c1/BilnK5Z6zFNUlFh9Kdft87oJlmnDkyfbnS5jIuDjAz485+hb19jgJVSZF0w1CjXf76jVeM9QQvK4VIT9lakXNAiclZPVwMtEmLrOJ0ULJoG+ES13bthzhx44QUjMP71r6fOCRMiqS5GjwJOFw53OkgoLw/eeAPWrePEmrXEHT0CwPVTXuSPNw5gWMwhOHLEWBSnchxEFW+mGbZ6sFKwecoChGuWwNtBXhE57NXVMn2FLmaDDOeeB8IHZ54Jc+caee877zReW7bMGD06ahS8+aYxR4yJXC0s7qom72zh8NkDUhj22w9G2UePhg0bjI23b4cZMzjy/U6Wn9WbaX/4M5eNe461sUZKKau8mVHTrxX4wfPSh3ZcVtJMnhZOt3Jh9VCJyLQPOL/VrN2HuUo49zwQfmja9OTjDh2MNNCyZUYtOi7OmEF0+XJjKoHKFIk/3PXccZY3blingofbQb+mZQybfpkxOnbQILhvx8kPbdUKrrnGqM0PHAiFhVzxbPYp57WnRmJPvZAisYumI08N69HQaSRig78zkdJQI0zUrRvMn2+0Baxda1wEtm07mQYaPdoY0HTeedCxo9Ew2rmzV5MF1gggWtPgRBFNDh/llcW/8tZj16NKSyn7859pvj+fs47uJ7nwZ+qUlxuN0o8/btyRpKcbCxL16GGU4YwzTu4gIQESElwu/+hpGUZ3FaFI7KLpKNovfhBlwT9SGmpEEMTEGCmSvrWmib7gAiNvvnIlvPKK8VrPnrB+vfG4Xz/YuZOjdeLYexyOV2jyUrsS98zT7C0s5u2X76bNkZ9pWPI7cRXGlL/LOl8Cj13P0F4psP8b406kZx9o357sxNOZfbgR66e/Z5yf/2+u2/MzK7cABU7bD9wFak8VoUivKEX7xQ+iLPhD5PU8EEE2caLxD4xFQXbuNP6vcuGF/FSvKdt++Jk4SqijNT/rOP699BsaJcaS3boz205P5UhCfX5NbEhhYkMKU9oztOrnd+6s/qjqNFFd73ujzV6R57Lh+NJOSS7e8VwRivSKUrRf/CBCe/sIEUqupv5tUi+WktIKr3vM+DOFsLteQ9KLzb1o7+0TdTV/Eblq/7Fe2imJT7cfCPofr6s8cGFRKXOuTfc6gPiTZ3aVnvD0c8JzFiDSswQS/EVEcNaz5tV1P1W/H8wBfe7yw74EEH/yzNMGdmTKG5t8zvmHUrjWoM1i1+OPyH7+Ivo465pXW7D6qU8b2JHE2JqLBvk7D46vnzOsezI39G5L7c6odslPR0N/eXfsfPwS/EVE8DbFEYxUiKtBhf7Mg+PP58wclsaca9MD3n8wRPpgMU/sfPyS9hERwV3uu/Z2wWBWftjfz7Frfjoc+8ubmaax8/FLzV9EBGcpk9rskgqJJq4utnZpj6jN7DSNnY9fgr+ICM5SJqN7t7VlKiSamNUeEipmp2nsfPyS9hERw66pj2gWboPFzE7T2Pn4JfgLIYIqnC7KwZjWwa7HL2kfIYSoZOc0jdmk5i+EEJXsnKYxmwR/ETHsOpJShBe7pmnMJsFfRAR3C6dEwx+yXckF2b4k5y8igp1HUkYrO09tICT4iwhh55GU0UouyPYmwV9EBDuPpIxWckG2Nwn+IiLYsYteVm4BfWatInX6e/SZtSrq0h1yQbY3afAVEcFuXfTCuQHarEZaf5ZClAbi0JFlHIUIAn+WZLSD2hctcL/0pDef520wN3vf0UqWcRTCQuGa73bXSBvsqabN3rdwT4K/EEFg1hwxoU6DWHnRCtcLZriS4C8sE275XV/K60++29n+Qt1uEIyJzcJh39FIevsIS4TbACBfy2vG0o5W9JO3steUHXtsRTKp+QtLhFt+15/yBjpHjBVpECt7Tdmtx1akk+AvLBFu+V0rymtVGsTKic2iZVI1Owgo7aOUaqqU+lgp9X3l/02cbJOulPpSKbVVKbVZKXVtIPsUkSHcBgBZUV5Jg4hgCjTnPx1YqbVuD6ysfF5bEXCT1roLMAiYq5RqHOB+RZgLt8BmRXnNaDcQwpVA0z5DgczKxwuB1cC9jhtorb9zeLxXKfULkAQUBrhvEcbCLb9rVXklDSKCJaARvkqpQq11Y4fnh7XWp6R+HN7vhXGR6KK1rnDy/nhgPEDbtm177N692++yCSFENDJthK9S6hPgDCdv3e9jgVoCrwBjnAV+AK31fGA+GNM7+PL5QojQC7exGuIkj8Ffa325q/eUUvuVUi211vsqg/svLrZrCLwHPKC1Xud3aYUQthHOk9eJwBt8lwNjKh+PAZbV3kApFQe8Dbystf5vgPsTQtiELNYS3gIN/rOA/kqp74H+lc9RSmUopf6vcps/AhcDNyulNlX+Sw9wv0IIi4XbWA1RU0C9fbTWh4B+Tl7PBsZVPn4VeDWQ/Qgh7Efm4glvMrePEMIv4TZWQ9Qk0zsIIfwSbmM1RE0S/IUQfpNBaOFL0j5CCBGFJPgLIUQUkuAvhBBRSIK/EEJEIQn+QggRhST4CyFEFJLgL4QQUUiCvxBCRKGAFnMJJqXUAcCM1VyaAwdN+JxwEm3HLMcb+aLtmAM53jO11kmeNrJt8DeLUirbm1VtIkm0HbMcb+SLtmMOxfFK2kcIIaKQBH8hhIhC0RD851tdAAtE2zHL8Ua+aDvmoB9vxOf8hRBCnCoaav5CCCFqiZjgr5QapJTKU0rtUEpNd/J+vFLqjcr3v1JKpYS+lObx4ninKqW2KaU2K6VWKqXOtKKcZvJ0zA7bjVBKaaVUWPcO8eZ4lVJ/rPyetyqlXg91Gc3kxTndVin1qVIqt/K8vsKKcppFKfWiUuoXpdQWF+8rpdRTlb+PzUqp80wtgNY67P8BMcAPwFlAHPA1cE6tbe4Anqt8PAp4w+pyB/l4LwXqVT6+PZyP19tjrtyuAfAZsA7IsLrcQf6O2wO5QJPK5y2sLneQj3c+cHvl43OAXVaXO8Bjvhg4D9ji4v0rgA8ABfQGvjJz/5FS8+8F7NBa79RanwAWA0NrbTMUWFj5eAnQTymlQlhGM3k8Xq31p1rrosqn64DWIS6j2bz5jgEeAR4HSkJZuCDw5nj/BMzTWh8G0Fr/EuIymsmb49VAw8rHjYC9ISyf6bTWnwG/utlkKPCyNqwDGiulWpq1/0gJ/snAHofn+ZWvOd1Ga10GHAGahaR05vPmeB2NxahBhDOPx6yU6g600Vq/G8qCBYk333EHoINSao1Sap1SalDISmc+b453BjBaKZUPvA9MCk3RLOPr37lPImUNX2c1+NrdmLzZJlx4fSxKqdFABnBJUEsUfG6PWSlVB5gD3ByqAgWZN99xXYzUTybGnd3nSqmuWuvCIJctGLw53uuABVrrfymlLgBeqTzeiuAXzxJBjVmRUvPPB9o4PG/NqbeE1dsopepi3Da6u+WyM2+OF6XU5cD9wBCt9fEQlS1YPB1zA6ArsFoptQsjR7o8jBt9vT2nl2mtS7XWPwJ5GBeDcOTN8Y4F3gTQWn8JJGDMgROpvPo791ekBP8NQHulVKpSKg6jQXd5rW2WA2MqH48AVunKVpUw5PF4K1Mgz2ME/nDOBVdxe8xa6yNa6+Za6xStdQpGO8cQrXW2NcUNmDfndBZGwz5KqeYYaaCdIS2lebw53p+AfgBKqc4Ywf9ASEsZWsuBmyp7/fQGjmit95n14RGR9tFalymlJgIrMHoNvKi13qqUehjI1lovB/6DcZu4A6PGP8q6EgfGy+OdDdQH/lvZrv2T1nqIZYUOkJfHHDG8PN4VwACl1DagHJimtT5kXan95+Xx3g28oJSagpH+uDmMK3AopRZhpOyaV7ZjPAjEAmitn8No17gC2AEUAbeYuv8w/t0JIYTwU6SkfYQQQvhAgr8QQkQhCf5CCBGFJPgLIUQUkuAvhBBRSIK/EEJEIQn+QggRhST4CyFEFPr/raXDkyFE340AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "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": [ "**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": 3, "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": 4, "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": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XlcVNX7wPHPYccNVDQVNdx3BHENF8q1zCXLNLPsq2lWZmlZ2qJlm2WmLfYrK9M2ra+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()" ] } ], "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 }