Skip to content

Instantly share code, notes, and snippets.

@agramfort
Last active July 23, 2020 21:36
Show Gist options
  • Save agramfort/c1d0307f4545a54642f011f79d9966b8 to your computer and use it in GitHub Desktop.
Save agramfort/c1d0307f4545a54642f011f79d9966b8 to your computer and use it in GitHub Desktop.
generalized additive models with backfitting and splines
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Generalized Additive Models (GAM)\n",
"\n",
"Author : Alexandre Gramfort\n",
"\n",
"GAM using backfitting algorithm were proposed by Hastie et al. in\n",
"\n",
"`Hastie, T. J. & Tibshirani, R. J. (1990). \"Generalized Additive Models\". Monographs on Statistics and Applied Probability 43.`\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Backfitting is a kind of coordinate descent method where a loop is run over each feature.\n",
"\n",
"The model reads:\n",
"\n",
"$y_i = \\sum_{j=1}^p f_j(x_{i,j}) + \\epsilon_i$\n",
"\n",
"where $\\epsilon_i \\sim \\mathcal{N}(0, \\sigma^2)$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy.interpolate import UnivariateSpline"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"class UnivariateSplineSmoother(object):\n",
" def __init__(self, k=3):\n",
" self.k = k\n",
"\n",
" def fit(self, x, y):\n",
" order = np.argsort(x)\n",
" x, y = x[order], y[order]\n",
" self._spline = UnivariateSpline(x, y, k=self.k)\n",
" return self\n",
"\n",
" def predict(self, x):\n",
" return self._spline(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example of spline smoothing"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3hUxf7H8fdsSQ9JINTQkSJFBIKgWBBEio1iQRCwgQXRa8HCVdGrP/AKeBVRUcQCCqg06UUpCggYmvReJPSSRtqW+f1xQiBhU8Cc3U3yfT3PPknYs2e+WeXD7Jw5M0prjRBCiJLH4usChBBCmEMCXgghSigJeCGEKKEk4IUQooSSgBdCiBLK5usCLhYdHa1r1qzp6zKEEKLYWLdu3SmtdXlPz/lVwNesWZO4uDhflyGEEMWGUupgXs/JEI0QQpRQEvBCCFFCScALIUQJJQEvhBAllAS8EEKUUBLwQghRQpka8Eqp55RSW5RSW5VS/zKzLSGEEDmZFvBKqcZAf+A6oClwp1LqKrPaE0KIYmnOHPjgA3A6i/zUZvbgrwbWaK1TtdZOYDnQ3cT2hBCi+Bk3DsaMAVvR33dqZsBvAW5SSpVTSoUAXYBquQ9SSg1QSsUppeJOnjxpYjlCCOFnHA5YuhRuv92U05sW8Frr7cB/gUXAAmAj4PJw3Bda61itdWz58h6XUxBCiJJpzRpIToYOHUw5vakXWbXW47XWLbTWNwNngV1mtieEEMXK4sVgsUC7dqac3tTFxpRSFbTWJ5RS1THG31ub2Z4QQhQrixdDbCxERZlyerPnwU9TSm0DZgMDtdYJJrcnhBDFQ2IirF1r2vg7mNyD11rfZOb5hRCi2PrlF3C5TA14uZNVCCF8Yf58iIiA6683rQkJeCGE8DatYcECuO02U+a/nycBL4QQ3rZlC8THQ+fOpjYjAS+EEN62YIHxtVMnU5uRgBdCCG+bPx+aNIGYGFObkYAXQghvSk6GFStMH54BCXghhPCuJUuMNWhMHp4BCXghhPCu+fMhLAzatDG9KQl4IYTwlvPTI9u3h4AA05uTgBdCCG/ZsQMOHvTK8AxIwAshhPfMnm18veMOrzQnAS+EEN4yezZcey1Uu2TvI1NIwAshhDecOgWrVsFdd3mtSQl4IYTwhnnzwO2WgBdCiBJn9myoXBlatPBakxLwQghhtsxMWLgQ7rzT2KLPSyTghRDCbMuXG0sUeHF4BkwOeKXU80qprUqpLUqpyUqpIDPbE0IIvzR7NgQFGTc4eZFpAa+UigGeBWK11o0BK9DTrPaEEMIvaW0E/G23QUiIV5s2e4jGBgQrpWxACHDE5PaEEMK/bNkCBw54fXgGTAx4rXU8MBI4BBwFErXWi3Ifp5QaoJSKU0rFnTx50qxyhBDCN2bMAKVKVsArpaKAe4BaQBUgVCn1UO7jtNZfaK1jtdax5cuXN6scIYTwjenT4YYbjCmSXmbmEM1twH6t9UmttQOYDtxgYntCCOFf9u6FTZugRw+fNG9mwB8CWiulQpRSCmgPbDexPSGE8C/Tpxtfu3XzSfNmjsGvAaYC64HNWW19YVZ7Qgjhd6ZPh+bNoWZNnzRv6iwarfVQrXUDrXVjrXUfrXWGme0JIYTfiI+H1auhe3eflSB3sgohhBlmzDC++mj8HSTghRDCHNOnw9VXQ4MGPitBAl4IIYraqVPG+jM+HJ4BCXghhCh6M2YYa79LwAshRAnzww9Qty40a+bTMiTghRCiKB07BkuXQs+exhIFPiQBL4QQRWnqVGN4pqfvF8+VgBdCiKI0ZQo0aQING/q6Egl4IYQoMocOwcqVftF7Bwl4IYQoOj/+aHx94AHf1pFFAl4IIYrKlCnQsiXUqePrSgAJeCGEKBp79sC6dX7TewcJeCGEKBqTJhnTIu+/39eVZJOAF0KIf0prmDABbr0VqlXzdTXZJOCFEOKf+uMPY/emvn19XUkOEvBCCPFPTZgAISE+X3smNzM33a6vlNp40SNJKfUvs9oTQgifSE831p7p3h3Cw31dTQ42s06std4JXAuglLIC8cAMs9oTQgifmDMHEhL8bngGvDdE0x7Yq7U+6KX2hBDCOyZMgCpVoF07X1dyCW8FfE9gspfaEkII7zhxAubPh4ceAqvV19VcwvSAV0oFAHcDP+Xx/AClVJxSKu7kyZNmlyOEEEVn0iRwOqFPH19X4pE3evCdgfVa6+OentRaf6G1jtVax5YvX94L5QghRBHQGr78Eq67Dho39nU1Hnkj4B9EhmeEECXN6tWwdSv07+/rSvJkasArpUKBDsB0M9sRQgivGzcOwsL8ZmlgT0ybJgmgtT4HlDOzDSGE8LqkJGPue69eRsj7KbmTVQghLtekSZCa6tfDMyABL4QQl2/cOLjmGmPtdz8mAS+EEJdj/Xrj0b+/sTywH5OAF0KIy/H55xAUBL17+7qSAknACyFEYZ09C999Z4R7VJSvqymQBLwQQhTW118bF1efecbXlRSKBLwQQhSG2w2ffAI33gjXXuvragpFAl4IIQpjwQLYt6/Y9N5BAl4IIQrn44+hcmXo1s3XlRSaBLwQQhRk926jB//EExAQ4OtqCk0CXgghCjJmDNhsMGCAryu5LBLwQgiRnzNnjGWBH3zQGKIpRiTghRAiP2PHGlMjX3rJ15VcNgl4IYTIS3o6jB4NHTsaa88UMxLwQgiRl++/h+PHYfBgX1dyRSTghRDCE7cbRo40bmpq187X1VwRUzf8EEKIYmvePNixw+jF+/mqkXkxe8u+SKXUVKXUDqXUdqXU9Wa2J4QQRUJrGDYMqleH++7zdTVXzOwe/EfAAq31vUqpACDE5PaEEOKfW7IE/vgDPv0U7HZfV3PFTAt4pVQEcDPwMIDWOhPINKs9IYQoMv/5D8TEwKOP+rqSf8TMHnwt4CTwtVKqKbAOeC5rI+5sSqkBwACA6tWrm1iOEKK0m7khnhELd3IkIY0qkcEM7lifrs1ich60fDn89psxPTIw0DeFFhGltTbnxErFAquBNlrrNUqpj4AkrfUbeb0mNjZWx8XFmVKPEKJ0m7khntembybN4fL4fFSInaF3NaLr4H6wZQvs3w/BwV6u8vIppdZprWM9PWdmD/4wcFhrvSbr56nAqya2J4QQ2XL31lMznXmGO8DZVAeT/zeZrr/+CqNGFYtwL4hpAa+1PqaU+lspVV9rvRNoD2wzqz0hhDgvd289PiGtUK975rfvORsaSdQTT5hZnteYPYtmEPB91gyafcAjJrcnhBCMWLgzz966iyTSrGvJsOzAqY6gcWIhhNv22rnp4EbebvcwQ0NDvVyxOUwNeK31RsDj2JAQQpjliIceu0P9TaLtJ85Zl4NyYdFh2HQMFh2Ik9O8sewAByPg/et/Iu2Xigy+YTDlQsr5oPqiI3eyCiFKnCqRwdnDMm4ySLRNJsk2A4uyUdl6D/rczQToOiiMO1Q77VxJyyPDef7O7jStnsGIVSP4ZO3nVNSP4UxuS0xkiOcZN37OtFk0V0Jm0QghisL5Mfgk59+cDHgPh2UfEe4OfNh5BA+3bpp9zFuztpJ8Lp1F4weC1cLmeb/TNbY6Hy3/ldeW/os0tYUQ1w2Uy3yOUHsZerSIYemOk/lPs/QyX82iEUIIn+jaLIZtp9fy5soX0Boa2P6P4V365Qjjrs1ijJ/Hj4czh2HGDOrEGvfi/PiHonz6MJJsM0iwTcAR+DIVMt/i+9UuzneJ4xPSeG365uxz+SPpwQshSpx5u+fR48ceVI+ozoLeC6gVVcvzgcnJUL8+1KgBq1ZlLypW69W52UGeZtnIyYBhWAimYsYw7DpnmMdEBrPyVd+tNplfD16WCxZClCgL9yyk2w/daFS+ESsfXZl3uAMMHw5Hj8L//pdjxcgqkRfmwAe7r6VSxn/RODge+BoOFZ/jFJ4u6PoLCXghRImx5vAauv3QjYblG7KozyKiQ6LzPnjfPvjgA+jTB1q3zvHU4I71CbZbs38O0LWolDEMjZMTAW/hIjH7uYv/MfA3EvBCiBJh75m93Dn5TiqHV2bhQwspG1w2/xcMHgw2m9GLz6VrsxiGd29CTGQwCmMY5pFWN1PN9RYudZqTAe+iySTYbmVwx/rm/EJFQC6yCiGKvVRHKt1+6IZbu5nfez4VQivk/4IlS2D6dHj3XWPVSA+yL8JeJLZGWV6Zl8Au5384Fzqa/90xwW8vsIL04IUQxdyM9Yep+V43Nh/fQkXnYLYdKuAu1MxMePZZqFkTXnjhstrq2iyGnf9+m/fav8dp9zI2J3115YV7gfTghRDF1swN8Tw9cwQnrYuIcD5IanqTgqcujhoFW7fCrFlXvKDYy21eZtupbby9/G1urH4j7Wu3v9JfwVSF6sErpSoopboppQYqpR5VSl2nlJLevxDCp4bOn80xy2cEuZoR4ewJQJrDxYiFOz2/YO9eYzOPHj3grruuuF2lFJ92+ZSry19Nr+m9OJp89IrPZaZ8Q1opdatSaiEwF+gMVAYaAq8Dm5VSbyulyphfphBC5JTmSGN7xjCslCE68yUUF2a9eJy6qDU89ZSxBd/o0f+4/dCAUH667ydSMlPoNb0XLnfeSxH7SkG98C5Af611S631AK3161rrl7TWdwNNgQ1AB9OrFEKIXN5Y+gYOy9+Uy3wOKxE5nvM4dXHSJFi82Jg1U6VKkdTQsHxDPu3yKcsOLOOd394pknMWJbmTVQhR7Px+8Hdu+eYWOtbsw4E9D+ZYGjjYbmV49yY5x+CPH4fGjaFOHVi5EqxWD2e9cn1n9GXS5kmsemwV18VcV6TnLsg/vpNVKTUxaxPt8z/XVEr9WlQFCiFEYaVkpvDwzw9TK6oWPz34ySXz1S8Jd62hf39jWYKvvy7ycAcY3Xk0VcKr0GdGH1IdqUV+/itV2Fk0K4A1SqkXgBhgMPCiaVUJIUQehvw6hP1n97Ps4WWEBYTRtVlY/nPRv/kGZs827lq9+mpTaooMiuTre77mtom38eovrzK68z8f4y8KhQp4rfXnSqmtwFLgFNBMa33M1MqEECKXP+P/ZMzaMQxsOZCba9xc8AsOHoTnnoNbbjG+mqh97fY81+o5PlrzEXfXv5vbat9manuFUdghmj7AV0Bf4BtgnlKqaSFed0AptVkptVEpJYPrQogr5nQ7GTBnAJXDK/N/7f+v4Be4XNCvnzFE8/XXYDF/Zvfw9sNpEN2AR35+hIT0BNPbK0hhf+MewI1a68la69eAJzGCvjBu1Vpfm9dFACGEKIzRa0az8dhGRncaTZnAQszOfucdWL7cmBJZK58VJYtQsD2Yid0mcjT5KIPmD/JKm/kpVMBrrbtqrU9c9PNaoJVpVQkhxEUOJR7ijaVvcGe9O+l+dfeCX/Drr8YNTX37wsMPm17fxWKrxPLGzW/w3V/fMXXbVK+2nVtBNzq9rpTyuCSb1jpTKdVOKXVnPqfQwCKl1Dql1IA82higlIpTSsWdPHmy8JULIUqN5xYY4+djOo9BXbRuu0fHjkHv3tCgAXz6aY513r1lyE1DaFmlJU/MecKnd7kW1IPfDMxWSv2qlBqhlHpZKfVm1rTJzcBdwJp8Xn+j1ro5xl2wA5VSl1wV0Vp/obWO1VrHli9f/op/ESFEyfP6zM1U/ve7zNwxk8DU+xi3LCn/Fzgc0KsXJCXBTz9BaAELj5nEbrUzsdtE0hxpPDbrMXx1v1FBAX+v1roNsBDYCliBJOA74Dqt9fNa6zy73Vrr+KyvJ4AZgHfvABBCFFuvz9zMxNX7OGX7HJu7MuHOrny3+hCvz9yc94teeAGWLoWxY6FRI+8V60H96PqM6DCC+Xvm8/m6z31SQ0EB30IpVQXoDcwCPgcmAH8C+S7DppQKVUqFn/8euB3Y8o8rFkKUCpPX/E2ydQ5Oy2GiHI+jsGf/uUeffw5jxsCLLxpj737g6ZZPc3ud23lx0YvsPr3b6+0XFPBjgV+BBkDcRY91WV/zUxFYoZTaBKwF5mqtF/yzcoUQpUWmPkuCfRJBrhYEuy98+Hd5Gu5YtgyeeQY6d4b//td7RRZAKcVXd39FoDWQPjP64HQ7vdp+vgGvtR6ttb4a+EprXfuiRy2tde0CXrtPa90069FIa12IiatCCGFItE9Ak0FZR38UFy6UWnNfNN2xw1j+96qrYPJkU5Yi+CdiysTw2R2fsSZ+DcN+H+bVtgs7TfIpswsRQojz/oz/k2TbYso478Guq+Z47sFW1S78cPgw3H67sQTwnDkQEYE/eqDxA/Ru0pu3l7/N0v1LvdaubNohhPALMzfE0+a9JdR8dTbtxvcjMjCaAc1eyu6xW5XiodbVebdrE+MFZ85Ax46QkADz5xsrRfqxz+74jHrl6tFzWk+OJB/xSpuyZZ8QwudmbojntembSXO4SLEuI0Vvp1LqC7SpXZ0RPa6/9AXJyXDnnbBnDyxYAM2aeb/oyxQeGM7U+6Zy3ZfX0XNqT37t+yt2q93UNqUHL4TwuRELd5LmcOEmlQT71wS46xGQ2dbz1ntJScbF1LVrjTH3W2/1fsFXqFGFRoy7axy/H/qd1359zfT2JOCFED41c0M88Vlb7CXafsClzlI28wkUlku33ktMhE6dYM0a+OEH6F6IZQv8TK8mvRjYciCj/hjFl+u/NLUtGaIRQvjM+aEZAIeKJ8n2M6HO2wjU9YFcW++dPg1dusD69fDjj9Ctmy9KLhIfdvqQPWf28OScJ6kSXoUudbuY0o704IUQPnN+aAbgrH0cCjtRjn6AsfXe4I5G0LN/P7RpA5s2wbRpxTrcAWwWGz/e9yPXVLyGHj/2YMn+Jea0Y8pZhfChmRviGbFwJ0cS0qgSGczgjvVz7PhT0PPCe84PwaRa/iTNGkeU41GsRAFc2Hpv3Tq44w7IzDQ2zb7pJl+WXGTKBJZhUZ9FtP2mLff/dD/7n9tPeGB4kbYhAS+KpbxC+uLZGADxCWnZQwCFeV54z8wN8ViUwqkzOGv/Apu7KuHOuwBjb9WuzWKM3nq/fhAdbawxY9KWe74SHRLNL31/YdfpXUUe7iABL4qh12du5rvVh7J/jk9IY/DUTUDOj/znpTlc2bMxXvxx0yW3uqc5XPzrh42MWLjzinvz8qng8pz/h9alNUm2GTgtR6mQ8Q4KO8F2Ky+3rwOvvmosO9CqFcyYAZUr+7psU1QKq0SlsEqmnFsCXhQrMzfE5wj38xwuzduzt5KQ6vD4ur8TzvLC9LmcU8dwWZLQZKBVGho3SgegCGBXUij/mr6JU2k38kjr5lgthbvlXT4VXL7z/xA71QkSbT8S4rqBYHczrErxwS2V6DzkcfjlF3jySfjwQwgM9HXJxZIEvChWPM6LznI21UFMZDCHEk6SYdlChmU3DsteMi37cKkzxkEF/B9/Cui/GAYtDeKaitfQvFJzWsa0pEPtDlSLqObxNfl9apCA9+z82PtZ21cARDkeB6D9zlV0/vJz40am8ePh0Ud9VmNJIAEvipVL5kUDGjeZahep1tVYg3ZyOGgLKA3agl1XI1Q3w+KMwarLY9MVsOoIlA7GQjCg0GSiceBWyTjVKdzqFH1bB7Dh2AYmb5nM2HVjAWgQ3YBOdTpROaAdM9eEcDQxnSqRwdlzuHPL68+FMf1xT9IfpNpWEOF4iKj0UF7/9SMe2LzYuCv1u++gYUNfl1nsScALv5TXmPbFgZqh9nDO9gup1lVGD11bqWqNpWeD59myrxpJSTWoGhnB4I71GbFwZz6BmzXXWkcToGsRExnMBx3bGX+kNVtPbmXx3sUs2reIj9d+gkt/iM1dmVBbWxwJHbER7fGsl6x6KLI936E2vWc/jM1Vkb5/RfPGkicom5rEzkcHUf+zkRAQ4OsSSwQJeOF38hvTfurWijw/5xOSLIvItOwFbSfYHWuM4bpaEnmuPJMHtfN43sE/bcLhvnCB1aKMEL74z3LMvcZYz7txhcY0rtCYWkH3s3XzClKtf3DOuoxE25Ss8eM2lHHeTaBukKM9j+uWCwAOZU6n/vG/+XxRba7f9yHbqtZny+eTuLV3Z1+XVqJIwAu/42lMO9G5nyfmfESSZRnptnTs7tqUzXySEFdbrIRlH+dpCCdbrg611aJ4oGU1lu44WajZLyMW7sRCGGGuDoS5OuBQx0i2ziHFtphU228EuZoS4exNkNsYWoiJzHfTs1Lr6KaVVHn6FTZuAkvEGfjsMxr2709DP1vHvSQwPeCVUlaM3Z/itdZ3mt2eKP7Oh7RGk2HZSpJtGmnWP1GuQJ5o9gj9W/Rn0LcJHodcIkM8r843YuFOHE43Qc4MQjPTCXJmEujM5PCSw6zs3RQIBYsF3Mdgwwnje4sFQkIgPBzCwjhyNhUuGnax60qUdT5OpLM3KbYFJNqmcTzwZYJc11JJP8bgjsX7bssit20beuRIoid8w90WTdJzTxL55jCIivJ1ZSWWN3rwzwHbgTJeaEuUAJUjgtibvJIE+yQyLTux6DJEOHpRxnkHHWNuoXnlGAZ3jGfw1E04XMYwSEhmGrXOHqFG0gn+emUN1+hkYzOI+Hg4c4ZpB44QlZZMoMvDNMqPC1fXHmXhnD2I1IAgEoLCORMSkfUoQ0y9GqxL7cVfATvZXXYNB8sMYs6B9bS66v+oHJ5z/napmjPvchl3n378McybhysogLEtNLYhr/PUPe/4uroSz9SAV0pVBe4A/g94wcy2RMmw4tAKEsP/zYnMNVjdFSib+RShrvZYCALgozmb6eo6Std160hfPJeqJw5R+0w8VZJP5TxRQABUrQoxMXDVVfwZVI14SzCJweGkBASTbgsgwxZAaEQ4wx+MNXrmbveFh9bgdEJamjFlLyWFPbuPsPqvgwSmpxKVlkTZ1CQaHd9Lpcxkgtcn0z7X75Jh/Zr4iG84XL02VRq2wlK9OhstEczd7yQoJJrAiPLEJ1Ay58zv2AHffw/ffGP8Q1u+POlvDiHWOh57xcr8eddQX1dYKpjdg/8QeBko+ntwRYmy4egGhiwZwoI9C6gUVomymU8S5rydaolnaPX3Sloe3kqTY3uod+ogvGmMz3cJCGFfuar8Ub0J+8pWZV/ZGA5FVuJ4eDRxo3vlGE5xbYjno4su3IJxQXV49yZQyGCtD2z30Puu1SzGWCfl1Ck4csQItL//JnX3Zg7FzcESvxfLwr+plOTmWoeTcRed81RIBIcjKnBmZmXodB3UqAE1axpfa9SAMsXkg6/TCXFxMHeusbzA9u3G+9+xI2uffZ3BGTXZkPExyZaTvF93PDaLXP7zBtPeZaXUncAJrfU6pVTbfI4bAAwAqF69ulnlCD9zfpjiUMJhMkInccq9iKjgKD5t+m8ePV6F5ROn03jf1OyeeUJQGJsq12Ndw1b0HdgdWrSg85R9xCemX3LumMjgHOEOF3rH/3RopGuzGM+vCQiAKlWMR2wsAFHALVozZcsUeix4lqTUs9Q6fReNzrQmJuksVRNPZD2OU+PIPhizFtJz/T5RURcC39PXyMhLflevSEoylu3980/4/XdYvtz4M4sFbrkFnn4aunVj5gnjE8pZ13qSA+cS7rybr5ZaqRsVX7I+sfgppU2ayqWUGg70AZxAEMYY/HSt9UN5vSY2NlbHxcWZUo/wHzM3xPPK9D85rqeSYpnK9X876bH7Kvodh6g9xp2q6dEVWFKhAatiGrG2WiN2R1cnKMB+YYVBLp1OCRd65f4WHqdST/HioheZsGkCNndVojNfIlBflf28VSn2DusMJ07AgQNw8KDnr+fO5TxxePilwV+jhrE4V1TUhUd4+OX9Q5CaanwiOXUKjh0ztsbbvfvCY/9+YxgL4KqroH17aNfOeERfuC+gzXtL+DvhNEcCn0Fho3LGaCwEERMZzMpXPU9nFZdHKbVOax3r8TmzAj5XAW2BlwqaRSMBX/JprWn47pvE7P6YHtsT6bHdToVzDhwWK5trNKb5k72M7dgaN2bmxiMF9riL2wXLikPe4XTAR7hIINLZmzLOHiiM6YEK8v8dtDY2ms7vH4DERM8NWyxGbz842Pi0ERhofLVYwOEwhpjOPxITjYDPLTwc6tY1Ho0aQcuWxqeVaM83egHUenUup+yjSbH+QqXM/xLovjr7d93/3h2X8c6JvOQX8DIQJkyVHcBnU2mTHMddW8exbN0RKp6DVLudJbVbMb9+G36r3ZzkwFAOvHzhL32ewyEXKcwx/mLmhnjCdAsC0sdwxv4JCfYJpFniiHa8iE1XRFPAQmVKQblyxqNFC8+NJCTAoUPGPwRnz176SE/PGeYulxH05x92u/EPQXT0hUf58lCnDlSocNnDQUHhm0hxLqKM497scIdcOzUJ03gl4LXWy4Bl3mhL+I+ZG+IZNWE5nTcupPvWn2l4MpE0GyysW5uF9e5lWe3rSAsIyj6+JN/af/HyuFbCiXa8wjn3dZyxf8aRwGco63iaMJexefQ/WqgsMtJ4+IHjKcc5bv2IQEctIp29s/88993CwjzSgxdFT2tYtozIQUNZum0FNq1ZVRWe6VKXX+u9SFpgVY8vK8m39ue+O1ehCHO1I8jdiFP2UZwOGEW6cz1lHU9hIST/O3KLAZfbRe/pvUlzJfFe2+/4abW12AyjlSQS8KLoJCfDxInwySewbRuNg22MukEzsWlFTkY8S7C7ab4vL8m39ucV2DZdkYqZw0m0/UCibQoZlu1EZw6mdsS1Xq6waL3727v8uv9Xxt89nkebdeBfbX1dUekkm26Lf+7oUXjlFePGooEDOaVTefq+EKq9AMPbPkRy+Ngc4R4ZbOzac7GS/rE9rzHnqBA7IfYAIp29qJj5HhoXxwJfpnbtxbjcLo+v8Xe/7PuFt5e/Tb+m/Xjk2kd8XU6pJgEvrtyuXdC/vzE9b+RIzrW7keeGXkf5+w+w9KZGVNBjiHT2RHFhfRi7VfHW3Y0Y3r0JMZHBKIyeuz9ObSxKgzvW9/iP2tC7LrwXwe6GtAj6kjYxXfhu+3t0mNiB+KR4H1V8ZQ4mHKTXtF5cXf5qPunyCaoEX1cpDmSIRly+bdvg7bfhp58gIAD3I4/w3W0VeHrXByilGNN+DN//Wp80V8YlLw0NsGUHeUkO9NwKutHq4vdC67v4ZuM3DDSFlSgAABWpSURBVJo/iGvGXsNXd3/FPQ3uKVQ7vpw2mpCeQJdJXXC4HUy7fxqhAaFeaVfkTQJeFN6ePUawf/89hIbCa6+xq1cnHl79Cn9s/YNOV3Vi7B1jqRFZg5HT53o8RWKa5z1TS4PCTulUSvFIs0doU70ND057kK4/dOWp2KcYdfsogu15X6fw5d6wP63bT/+5PUh076JRwHvs+DucBnlPjxdeIkM0omAHDsBjj0GDBsY6I4MHk7lnJ+92DKbJjNvYeXonE7pOYF6vedSIrAHkPeYs858Lr165eqx6dBUvXv8in8V9RstxLdl8fHOex+e3N6yZZqw/TP9ZT5CoN1DOMYiUpAY8/8NGXp+Zd63COyTgRd5OnoRnnoF69Yxe+zPPwL59xD13Hy1ndOaNpW/QrUE3tg/cTp+mfXKMt+Y15lySL6SaIdAWyMjbR7LwoYWcSj1Fy3EtGbN2DJ7uQM9rpo6ZUy7d2s1TcweSaFlMhONBwlzGmpoa+H71IWZuKF7XEEoaCXhxqYwMGDnSuCV97Fij975nD6kjhjH4r1G0+rIVp1JP8XPPn5ly7xQqhFa45BRdm8WUugupZrq9zu389dRftK/dnkHzB3HPlHs4lXqKmRviafPeEmq9OhdLHhc0zfrU5NZuBswewHH3LMo4ehDh7JXjeQ2mf3oQ+fPKWjSFJWvR+JjWxhDMyy8bi0ndcQeMGAFXX83S/UvpP7s/e8/uZUDzAbzf4X0igiJ8XXGpo7Xm47UfM3jxYEJtUYSk9sea2QqVez/CLGYtvuZ0O3l81uN8u+lbqlofwpLygMcaZM0Z88laNKJgcXHw/POwYgU0aQKLFkGHDiSmJ/Ly7Cf4Yv0X1Imqw5K+S7i11q2+rrbUUkrxbKtnuaXGLdzwRQ/OWt8lOKAVZR1PYNPGJymrUri1LtJZNBfPzqkQ4SIz4gM2nlzO223f5poyj/H8Dxvx1FWUay6+JQFf2p05A0OGoL/4grOhEYzo+Awrbr6bF6Mbonb8zNPznuZYyjEG3zCYt9q+RYg9xNcVC6BppaaUTx1Fom0WibbvORL4NGWcPSjj7IbSQUXaa754dk6G2sP69P/izDjBwGtH8OYtLwEQd/AM368+lCPk5ZqL70nAl1Zaw7ffwuDBuM+eZcJ1XRl5/YOkBIbgSjpO35n/IdmygmsqXsPMB2bSMqalrysWucREhkNCd0JdbThrH0+i/XtSbPOoYX8Yh+t27FbPG5BfrhELd5LqSCfRNo1E2xSsOpKKmcPYsKN59jHvdm1CbI2yxWrp5tJAxuBLo82bjR13VqyA66+nX8tHWB5cBY3mnPUXztrH4yaD6rY+7HltbJEFhShauee9p1u2kxTwFWlqOzHhMTzX6jmq2O/g0yXHrjh0tdZU+vebnLF/g9NymBDnTZR1PIWVMjK+7idkDF4YUlNh6FD43/+MJWXHj4eHH+a3IfNxqGOcsX9MunUTga5GlHMMwpJeVcLdj+W+O7ZOmea8dHtPAsI2MeqPUbz8y8tY9FCCXa0JsdzI4YTmhb7xKTE9kZ+2/cToNaM5EbgZm7sq5TPeJMR9XfYxMr7u/6QHX1osXQqPPw779jGrZReGXv8QIZUr8kKHOrw4fziHnN8AFqIcjxDm6oTCYsrFOuE9TYeNZW/aNNKsf+BW51A6hEB3A8oF1Od/3bpRI7IGEYERKKVITE/kcNJhNp/YzG8Hf2PFoRVkuDJoVL4Rt8Y8yuK4eqQ7LsyS8detEUsjn2/ZV1gS8CZITITBg2HcOFKq1eTptk/xW5VGAGSqfZwN/Jh0tZswd2siMp7Ehuf7y+UvdPFT69W5aEDjIN2yiVTrH2RYduNQB0HlvVJl4wqN6VC7Az0b96RllZYopYrd1oiliU+GaJRSQcBvQGBWO1O11kPNak94MHs2PPmksWny4MHcHdaWfakaTSYJtikk2aZhIZx6tjd5r3N/Ri7axZGENCxKXbL5xj/aZUj4RJXIYOIT0lDYCXbHEuw2MqByhJVPHq7AkeQjJKYbe7iWCSxD5fDKNIhuQFhA2CXnKk5bI4oLzByDzwDaaa1TlFJ2YIVSar7WerWJbQqA06dh0CCYPNmY0/7zzxAby/5X55JmWccZ+1iclqOEOm8jyvEYmenhdGtelW7NjZ2War3qeaGw4r7LUGkzuGP9HBdhwfgk9kqnJjSvHEPzys3zebUoCUwLeG2M/aRk/WjPevjPeFBJNW+esbTA6dPGyo+vvgoBAcQnxZMcOoLT7uXY3DFUyHiXYLexa1Dui2Xne365yUW14qWgJYpFyWfqLBqllBVYB1wFfKK1XuPhmAHAAIDq1aubWU7JlpwML74I48YZvfYFC6BpU5xuJ2NWf8gbS98gEwfRrr6EZHbL3oTD080oefX85KaV4keGVko3r1xkVUpFAjOAQVrrLXkdJxdZr9Dvv0O/fnDgALv6PcWAOndxMMVFeJl9JAeP5UDSNjpf1ZkxXcbw14HAQvXo5KKaEMWDX8yiUUq9CaRqrUfmdYwEfMEuDt6aoVa+2D2Tut99AbVq8du/R/LE/mBSHGdIsE8gxbYQmy7HCy2H8V6X/rJ9mhAlUH4Bb9pywUqp8lk9d5RSwUAHYIdZ7ZUG5+9cjE9I4+rj+xg7+gnqTvyc/T0egk2beOV4KMfd04kPGkCKdTHhjq5UTv+MFZuvknAXohQycwy+MvBt1ji8BfhRaz3HxPZKvBELd5Ke6eCxP3/mleXfcjakDA/f+xa7W9xEuxUziEsfgjPgMEGuZkQ5+hOgjWsaMvtFiNLJzFk0fwHNzDp/aTNzQzwZ8Uf4Zu7/uGX/ehbVbc3LnZ/lZEgSZ1NfYdWaOGxUoXzGmwS7W+ZYm1tmvwhROslaNMXAzA3xzB3xNfN/HkV4Ziqv3/40E669iYSAKSRb56AIJMrxKOHOu7Jnx5wns1+EKL0k4P1EnrNWMjJIGfgs4/6Yzo7oGvTq+RbrKm0h0T4ANylE6I6EZ/TGSpTH88ryAkKUXhLwfiD3sq/xCWm8Nn0zYQf2EjtkIA/t2MK3zbvw7/Z1OBY8DJflRNY4+8N8en8PRizc6fHGpJjIYAl3IUoxCXgvyW9e+YiFO3PcVITW3BU3nzb//YIMewB97+3JlIarcVjmEeCuQ7mMZwl2X5sjwOXGJCFEbhLwXpBXDx2MOw0vnuVSJj2FYQvGcOfOFfxeow79utrYHzUFm7sS0ZmDCXHdhMqa3Xo+wOWWdCGEJ7JcsBe0eW+JxyGUyGA7oYG27OdaHN7GR7NHUDHlNP9pW4FhNx5DqSjKOO4n3NUpxwXUyGA7G4fe7rXfQQjhn2RHJx/Lax56QpqDhDQHVreLZ1b9wLOrJnMw0k6bR92si0knwvE4Ya7OWAjM8bpgu5W37m7kjdKFEMWYBLwX5LU6I0DlpOP8b/Z/aH34IBOvgUFdgrFY+xCT3hkLQR5fIzNjhBCFIQHvBZ5WZ9S4aL/ra0bNn4XN7aZf1xBmNepJGWcXLE7PwQ4yM0YIUXgS8F5wPpDfmrWVs2nncLsWMnTJBB5fn8afVew8dXcvTkZ0J8IZkO95ZGaMEOJySMCbwNOUyFZXWTht/Z7ohBlMmpZM/VMwunVrRt/4Ek5rEAqwWxQocLguXPhWGLukxMjMGCHEZZKAL2IXT4nUaPYlbqDfzHdItazg6dVO3l+sOBsczkMPvMKqmtdmv+58gINMdxRCFA0J+CI2YuFOzjlSSbWuINk2h0zLbiqkBDNtVjS37TrGL3ViebnLvzgTEpH9mpjIYFa+2i77Zwl0IURRkIAvQuuPrmfzuQ84F7QMtzqH3V2N7jvv4rNZvxORfppNL7/DIHsL0pzu7NfIuLoQwiymbfhRWpxNO8snaz+h2efNaPFFC87ZFhPsaknV1P8wekFLpk2eTWJQOP0HfkLT/77O8B7XEBMZjMLoucuURyGEWaQHfwXSnenM2z2PyVsmM3vnbDJcGTSr1IwxnccQodsyccIa3p/+Hk2P7eb7azsxsuMTDH2gJSCbIAshvEcCvpCcbidL9i9h0uZJzNgxg6SMJCqEVqB/8/480uwRDh2ryIgFO2i9YjzTfxmL02rlya5D2NyqPUPlQqkQwgdMC3ilVDVgAlARY6bfF1rrj8xqzwzpznSW7F/Czzt+ZubOmZw4d4IygWXofnV3Hmz8IO1qtcNmsTFzQzzvT/ydN+Z+TOddq1hTrTGvdn2Z5x5px1gJdiGEj5jZg3cCL2qt1yulwoF1SqnFWuttJrb5j51JO8PcXXP5eefPLNy7kJTMFMICwuhStws9G/Wkc93OBNly3mm68sNv+HnaSGMlyLaP8GXLrrgtVkYs3Ck9dyGEz5i5J+tR4GjW98lKqe1ADOBXAe90O/kz/k8W71vMlL/msv1MHODGTjlurdGV52/sza01byXQFnjpi5OT4YUXGDHhS7aXr0mf+99hR4Va2U/LZtdCCF/yyhi8Uqomxgbcazw8NwAYAFC9enXTa9Fas+v0LpbsX8LifYtZsn8JiRmJKBSB+irKOO8lxNWKAF2XQ3vtpDdt4jncV6yAvn3hwAEm3tKTd2IfINOWcz9U2exaCOFLpge8UioMmAb8S2udlPt5rfUXwBdgrAdf1O2nO9NZd2QdK/9eycq/V7Lq71WcSj0FQPWI6tzX8D461OnAyFl2TiTmXAsmzeG6dJglIwOGDoX334eaNeG33wgPrYV1+maQHZWEEH7E1IBXStkxwv17rfV0M9sCI8y3nNjChqMb2HDMeKw/up5MVyYAdcvW5c56d9KmWhturnEzdcvWRSkFwCsT53o8Z45hlj/+gMceg+3boX9/GDUKwsPpmvW0LDEghPAnZs6iUcB4YLvW+gOz2sl0ZdJ/dn82HN3A9lPbcbqdAJQJLMO1la7l2euepU31NtxQ7QYqhFbI8zx5rdleJTIYUlLg3/+Gjz+GqlVh3jzo3DnHcTK/XQjhb8zswbcB+gCblVIbs/5siNZ6XlE2EmANYMuJLVQtU5W76t1Fs8rNaFapGbWiamFRhb9R19Oa7cF2KyPCj0LjxnDoEAwcCMOGQXh4Uf4KQghhCjNn0azAWO3WdOsGrPvH58i9cXUDeyafb5xE9TlToUED+P13aNPmH7cjhBDeIneyXqRrsxi6XlMJvvwShgyBpCR4/XVjeCYo712WhBDCH0nAX2ztWmMYJi4ObrkFxowxhmeEEKIYktUkAY4dM2bFtG4N8fEwaRIsXSrhLoQo1kp3Dz4lBUaONB4ZGfDii/Dmm3IRVQhRIpTOgHc4YPx4eOstOH4c7rvPmB1z1VW+rkwIIYpM6Qr4zEz47jsYPhz27IGbboKff4ZWrXxdmRBCFLnSMQafng6ffQZ16xp3okZEwKxZsHy5hLsQosQq2T3448eNKY+ffgpHjsD118PYsdCpEyivTNEXQgifKXkBrzWsXAmffALTphnj7R06wMSJcOutEuxCiFKj5AT8/v0webIxxXHrVmMYZuBAeOopqFfP19UJIYTXFf+AT0kxeuirVxs/33gjjBsHDz4IoaG+rU0IIXyo+Ad8WBjUqQPdusEDD0CNGr6uSAgh/ELxD3gwpj4KIYTIoXRMkxRCiFJIAl4IIUooCXghhCihTAt4pdRXSqkTSqktZrUhhBAib2b24L8BOpl4fiGEEPkwLeC11r8BZ8w6vxBCiPzJGLwQQpRQPg94pdQApVScUiru5MmTvi5HCCFKDKW1Nu/kStUE5mitC7X3nVLqJHDQtIIKJxo45eMaLkdxqxekZm8pbjUXt3rBP2quobUu7+kJv7qTNa8ivUkpFae1jvV1HYVV3OoFqdlbilvNxa1e8P+azZwmORn4A6ivlDqslHrMrLaEEEJcyrQevNb6QbPOLYQQomA+v8jqh77wdQGXqbjVC1KztxS3motbveDnNZt6kVUIIYTvSA9eCCFKKAl4IYQooUplwCulOimldiql9iilXvXwfKBS6oes59dkzef3qULU/LBS6qRSamPW43Ff1HlRPfkuNqcMo7N+n7+UUs29XaOHmgqqua1SKvGi9/hNb9fooaZqSqmlSqltSqmtSqnnPBzjN+91Iev1q/dZKRWklFqrlNqUVfPbHo7xu8wAQGtdqh6AFdgL1AYCgE1Aw1zHPA2Mzfq+J/BDMaj5YWCMr9/fi+q5GWgObMnj+S7AfEABrYE1xaDmthg37vn8/b2opspA86zvw4FdHv7f8Jv3upD1+tX7nPW+hWV9bwfWAK1zHeNXmXH+URp78NcBe7TW+7TWmcAU4J5cx9wDfJv1/VSgvVJKebHG3ApTs1/RBS82dw8wQRtWA5FKqcreqc6zQtTsd7TWR7XW67O+Twa2AzG5DvOb97qQ9fqVrPctJetHe9Yj9+wUf8sMoHQO0cQAf1/082Eu/R8s+xittRNIBMp5pTrPClMzQI+sj+BTlVLVvFPaFSvs7+Rvrs/6qD5fKdXI18VcLGtYoBlGD/Nifvle51Mv+Nn7rJSyKqU2AieAxVrrPN9jP8kMoHQGfEk1G6iptb4GWMyF3oQoOusx1v1oCnwMzPRxPdmUUmHANOBfWuskX9dTkALq9bv3WWvt0lpfC1QFrlNKFWp9LV8rjQEfD1zcu62a9Wcej1FK2YAI4LRXqvOswJq11qe11hlZP34JtPBSbVeqMP8d/IrWOun8R3Wt9TzArpSK9nFZKKXsGGH5vdZ6uodD/Oq9Lqhef32fAbTWCcBSLt3MyN8yAyidAf8nUFcpVUspFYBxQWRWrmNmAf2yvr8XWKKzrp74SIE15xpTvRtjbNOfzQL6Zs3waA0kaq2P+rqo/CilKp0fV1VKXYfx98enf4mz6hkPbNdaf5DHYX7zXhemXn97n5VS5ZVSkVnfBwMdgB25DvO3zAD8bDVJb9BaO5VSzwALMWanfKW13qqU+g8Qp7WehfE/4ESl1B6Mi249fVdxoWt+Vil1N+DEqPlhnxVM9mJzbYFopdRhYCjGxSm01mOBeRizO/YAqcAjvqn0gkLUfC/wlFLKCaQBPf3gL3EboA+wOWuMGGAIUB388r0uTL3+9j5XBr5VSlkx/rH5UWs9x58z4zxZqkAIIUqo0jhEI4QQpYIEvBBClFAS8EIIUUJJwAshRAklAS+EECWUBLwQQpRQEvBCCFFCScALkQelVMusxduClFKhWWuBF4s1SIQAudFJiHwppd4FgoBg4LDWeriPSxKi0CTghchH1to/fwLpwA1aa5ePSxKi0GSIRoj8lQPCMHYfCvJxLUJcFunBC5EPpdQsjB20agGVtdbP+LgkIQqt1K0mKURhKaX6Ag6t9aSslQRXKaXaaa2X+Lo2IQpDevBCCFFCyRi8EEKUUBLwQghRQknACyFECSUBL4QQJZQEvBBClFAS8EIIUUJJwAshRAn1/2BkbxMmCF3mAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"n_samples = 60\n",
"sigma = 0.1 # noise level\n",
"x = 3. * np.random.rand(n_samples)\n",
"xx = np.linspace(-.2, 3.2, 200)\n",
"\n",
"def f(x):\n",
" return np.sin(x ** 2) + 2. + x\n",
"\n",
"y = f(x)\n",
"y += sigma * np.random.randn(n_samples)\n",
"\n",
"spline_smoother = UnivariateSplineSmoother(4).fit(x, y)\n",
"y_pred = spline_smoother.predict(xx)\n",
"\n",
"plt.plot(x, y, 'o', label='data')\n",
"plt.plot(xx, f(xx), 'g', label='true')\n",
"plt.plot(xx, y_pred, 'r', label='estimated')\n",
"plt.xlabel('x')\n",
"plt.ylabel('f(x)');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Now the GAM implementation"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"class GeneralizedAdditiveRegressor(object):\n",
" \"\"\"Fit Generalized Additive Model with backfitting\"\"\"\n",
" def __init__(self, smoothers, max_iter=10):\n",
" self.smoothers = smoothers\n",
" self.max_iter = max_iter\n",
"\n",
" def fit(self, X, y):\n",
" n_samples, n_features = X.shape\n",
" self.y_mean_ = np.mean(y)\n",
"\n",
" residuals = y.copy()\n",
" residuals -= self.y_mean_\n",
" for i in range(self.max_iter):\n",
" for j in range(n_features):\n",
" if i > 0:\n",
" residuals += self.smoothers[j].predict(X[:, j])\n",
"\n",
" self.smoothers[j].fit(X[:, j], residuals)\n",
" residuals -= self.smoothers[j].predict(X[:, j])\n",
" return self\n",
"\n",
" def predict(self, X):\n",
" n_samples, n_features = X.shape\n",
" y = np.ones(n_samples) * self.y_mean_\n",
" for j in range(n_features):\n",
" y += self.smoothers[j].predict(X[:, j])\n",
" return y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def f1(x):\n",
" return np.cos(3 * x)\n",
"\n",
"def f2(x):\n",
" return x ** 3\n",
"\n",
"def f(X):\n",
" return f1(X[:, 0]) + f2(X[:, 1])\n",
"\n",
"n_samples = 100\n",
"x1 = 2. * np.random.rand(n_samples) - 1.\n",
"x2 = 2. * np.random.rand(n_samples) - 1.\n",
"X = np.c_[x1, x2]\n",
"y = f(X)\n",
"y += 0.2 * np.random.randn(n_samples)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAes0lEQVR4nO3df5DcdZ3n8ecrQwMTzmKiyQkZGBNrKVg4lOgc4sayBFlBWAFBBXfvVu+0snu31N5aXq7G8kpYq67IHrVr3Z7euSnWWvbWUlA0Zi/sRTFY3rEXl4khhgjRiAoZWIlKsqeZJZPJ+/6Y7tDp+X67vz397W9/u/v1qEqlf3zT34/t8P5+5v15f98fRQRmZjb4lvV6AGZmVgwHfDOzIeGAb2Y2JBzwzcyGhAO+mdmQOK3XA2hm5cqVsWbNml4Pw8ysb+zateunEbEq6b1SB/w1a9YwPT3d62GYmfUNST9Oe88pHTOzIeGAb2Y2JBzwzcyGhAO+mdmQcMA3MxsSpa7SMTMbJlt2z3D39v08e3iW1WOjbLzmQm5aN57b5zvgm5mVwJbdM3zkS3uZnZsHYObwLB/50l6A3IK+UzpmZiVw9/b9J4N9zezcPHdv35/bORzwzcxK4NnDs229vhQO+GZmJbB6bLSt15fCAd/MLMGW3TOs37SDtVPbWL9pB1t2z3T1fBuvuZDRysgpr41WRth4zYW5ncOLtmZmDYpYQG1U+1xX6ZiZFajZAmq3Aj4sBP1ufr5TOmZmDYpYQO2FXAK+pM9Iel7S4ynvv0XSEUmPVf98LI/zmpl1QxELqL2Q1wz/L4BrWxzzvyPisuqfj+d0XjOzzLIuxBaxgNoLueTwI+Kbktbk8VlmZt2QthA7/eOf8/CThxIXSru5gNoLRS7avlHSHuBZ4N9HxL6kgyRtADYATExMFDg8MxtkaQuxn935NFF93liN0+8BvlFRi7bfBl4VEa8F/iuwJe3AiNgcEZMRMblqVeK2jGZmbUtbcI2G53m3MyiTQgJ+RPxDRPyi+vhBoCJpZRHnNjOD9hZc+70aJ00hAV/SOZJUfXx59bw/K+LcZmaQvBCrlGP7vRonTS45fEmfA94CrJR0ELgDqABExKeBdwH/RtJxYBa4LSIaf5MyM+uapIXYKy9axQO7Zk7J7Q9CNU4alTnuTk5OxvT0dK+HYWYDrNubjhRN0q6ImEx6z60VzGyoDWI1Thq3VjAzGxIO+GZmQ8IB38xsSDjgm5kNCQd8M7Mh4SodM+s7g1ZKWRQHfDPrK73YfnBQOKVjZn2l2faD1pwDvpn1lUHdfrAIDvhm1lcGdfvBIjjgm1lfGdTtB4vgRVsza0uvK2QGdfvBIjjgm1lmZamQGaaGZ3lywDfrY0XPtptVyDgAl58Dvlmf6sVsu50KmV6nfmwxL9qa9ale1KNnrZCpXYxmDs8SvHQx2rJ7pmtjs9ZyCfiSPiPpeUmPp7wvSX8q6YCk70h6XR7nNRtmvahHz1ohk3Yx+vD9e1g7tY31m3Y4+PdAXjP8vwCubfL+24ELqn82AP89p/OaDa1e1KPftG6cu26+lPGxUQSMj41y182XLkrVpF105iM84++hXHL4EfFNSWuaHHIj8JfVjct3ShqTdG5EPJfH+c2G0cZrLjwlhw/F1KNnqZBZPTbKTIvfNLzYW7yicvjjwDN1zw9WXzOzJco62+6FpNRPklYXBctX6ap0JG1gIe3DxMREj0djVm5lrUdvvDlqmcR8xKLjxMICbxn/Nwyiomb4M8D5dc/Pq762SERsjojJiJhctWpVIYMzs/zdtG6cR6au4oebrueP3/NalHBMgLtcFqiogL8V+O1qtc4VwBHn782Gx03rxlk8v1/gLpfFySWlI+lzwFuAlZIOAncAFYCI+DTwIHAdcAA4CvyrPM5rZv1jPGUh110ui5NXlc57W7wfwO/lcS4zK0bed8r2qqrIXlK6RVsz671utG1wl8vec8A3s0W61SStrFVFw8K9dMxsEW8jOJg8wzfrY93qSJl2p+zY8grrN+1wSqZPKRJuhiiLycnJmJ6e7vUwzEqpMc9es2J5hTvecUnbgbj+4nH2aIVfHjvO3PxL8aEyIgiYO/HSa6OVkZN397odcjlI2hURk0nveYZv1qeS8uwALxyda7rA2hiYr7xoFdu+8xwvHJ07eczh2Tkqy8SK5RUOH51j9dgov3zxOIdn5075rPp2zGXYCcuacw7frE81y6en9cVP6lP/VzufPiXY18ydCF6oBvuN11zIkdnFx9TG0Yve/NY+z/DN+lSrjpRJF4S03wqaqc3Wzx6tLJrhAyyTUseR1yKv00X58AzfrE+16kh59mhl0WtLDcCzc/NIJJ4vqSlaTR530Xr3rPx4hm/Wp2oz3Du37kucef/y2HH+45a9PPzkoZMz47RZehaHj87xiVsva9kBsyavu2i9cXp+HPDNeiSPNEXtRqZ1H//qojz83Hzw2Z1Pn2xaNnN4lsqIqCzTKZU2Wa0eGz3lxqm1U9tSjx3PMe3iewLy44Bv1gPtti5odXE4nLDoCizqUDk3H5x1+ggn5k40nZ03Spqtp60hjI+N8sjUVZk/u5W087jpWvucwzfrgXaqWrLksNsJfr88Np8p2I9Ii3bS2rJ7hvWbdrB2ahtHjx2nsuzULvfdaIaWdeN0a80zfLMeaCdNkSWHndSJUiye4bfjj9/z2lN+i2j8reSFo3NURsTYaIUjs3Ndq55x07X8OOCb9UA7aYosF4ekoLjmFaP87Q9+vqSgf9bpI4sCatKFZ24+OOuM03jsjrct4SzZuelaPpzSMeuBdtIUaema2uu1NMuH7nsMgE/cehkbr7mQbz995JRgL2C00vo/+cqI+E/vvHTR61487X8O+GY9cNO6ce66+VLG6mrlz0wJxs0uDmn5/T/8632LZuMBnFkZWfRZtRYKtXz93e96beJsutWFx8rPKR2zHnrx+ImTj9N64DTLYa/ftCMxv592N21jLX07+XDvWNX/8trT9lrgvwAjwD0Rsanh/fcDdwO1soJPRsQ9eZzbrF+1c0NRWg673XRKYy19O7x42v86DviSRoBPAb8OHAQelbQ1Ir7bcOh9EXF7p+czGxR55MRT+9aPVnjx+IncZ+NePO1veeTwLwcORMRTEXEM+DxwYw6fa9YT9bXm6zft6FrPljxy4mn5/TtvuIS7br6U8bHRRbX0NrzySOmMA8/UPT8IvCHhuFskvRn4HvChiHgm4RgkbQA2AExMTOQwPLPsurF5d+1z63verFhe4frXnMsDu2YW1c5fedGqzJ/bKs3iAG/1Ot7xStK7gGsj4oPV5/8SeEN9+kbSK4BfRMSLkn4HuDUiWt577R2vrGjrN+3IvV3Alt0zbPzCnkX9ayoj4vI1KxbVytfvImXWrmY7XuWR0pkBzq97fh4vLc4CEBE/i4gXq0/vAV6fw3nNcteNWvO7t+9PbFY2Nx/sfOqFRTdGeeMQ65Y8UjqPAhdIWstCoL8N+M36AySdGxHPVZ/eADyRw3nNcpe2CLpMYu3UtiVVpjS7WKT1tGn8N94AxPLQ8Qw/Io4DtwPbWQjk90fEPkkfl3RD9bDfl7RP0h7g94H3d3pes25I21RkPmLJm280W4QdkRJfr/833gDE8tJxDr+bnMO3XqifTadt8jE2WuGsM07LNONulsO/9Z+fv2jhtjGH3411BRtczXL4vtPWrEGWTT4Oz86drLhpVcmTtDPViuUV7njHJdy0bpzJV728abrGPWwsLw74Zk202ii8ptWWe81uWGp1M5M3ALG8uHmaDaWsN1e12ii8Xrdm3N4AxPLiGb4NnSw3V9Xn8ceWVzjjtGUnN/k4euz4ov1joXszbvewsbw44FtfWUp5YuO/OXrseNOmZUk7O41WRvjErZclvg/dn3G7h43lwSkd6xtLKU9M+jdJs3N4KSXTar/ZWi9796mxfuMZvvWNdtoJN/s3aWopmaxbCjrAW79xwLe+sZTyxKwLqfUpmaVUxfhOWOsHTulY3xhbXkl8vVkgTntvbLSSmpJptyrGd8Jav/AM3/rClt0z/OIfjy96vTKipouladvy3XnDJU3r4iF7VcxSUk1mveCAb30hrePkWaefllhK2Rik2023tJOj952w1i8c8K00mgXstOB5pNqqIKm2fuMX93Dn1n0n6+drZZV5852w1i+cw7dSaJUHb7UdYFJaZW4+ODw71/W8uu+EtX7hgG+l0Kr2PSmoioVAntZNslG3NhZxXb71C6d0rBRa5cHrc/Ezh2cRnNwpqvH5Us7TqN0yS9flWz/wDN9KoVXKBhaC6iNTVzE+NroouAcLM/6lnqeeyyxtUDngWym0kwdPm6UHnEyrrFheobLs1EtA1rx6s/RS1i6bZmXklI6VQjvlk2lVMSPSKf826+c1SrugzBye5Q/ue+yU5802PjErm1xm+JKulbRf0gFJUwnvnyHpvur735K0Jo/z2mC5ad04G6+5kNVjozx7ePbkjLpR1n1nAR6ZuoofbrqeR6auyhyU2ymn7NZCsFk3dBzwJY0AnwLeDlwMvFfSxQ2HfQB4ISJ+BfgE8EedntcGT9bceWNVTNJG4J0E4nY2PYHWC8FOA1lZ5DHDvxw4EBFPRcQx4PPAjQ3H3AjcW338ReCtUsJ/pTbUWpVm1qst4P5w0/WcSNhkHJZ+p2vjBaWVVk3VvABsZZFHDn8ceKbu+UHgDWnHRMRxSUeAVwA/bfwwSRuADQATExM5DM/KrL78Ma2sslXg7sadrvVlls3q/FstBLvPjpVJ6ap0ImJzRExGxOSqVat6PRzrosbZb5pWgbvbd7qmpXjGRistb7Bynx0rkzxm+DPA+XXPz6u+lnTMQUmnAWcDP8vh3NbHsmxOkiVwd3vP104+3312rEzyCPiPAhdIWstCYL8N+M2GY7YC7wP+L/AuYEdESuLVhkarWa6AW16f7Q7Wbt/putTPT2vP7D471gsdB/xqTv52YDswAnwmIvZJ+jgwHRFbgT8H/oekA8DPWbgo2BBJalWQNvutCeDhJw8VN8gu6PZvH2btUJkn2pOTkzE9Pd3rYViHGlsXw8Is95bXj/PArpmmaR0BP9x0fQGjNBsMknZFxGTSe6VbtLXBk1ap8vCTh06WP6ZZfnr2engza86tFazrmrUq+NB9j7F6bJTllWUcnTux6JhfHptnzdQ2xjtsl2BmDvhWgGa5+trNSK3MHJ5l4xf2gBY2Nqm95l42Ztk5pWNtWUqbgHZbFaSZOxEng32Ne9mYZecZvmWWtG9slhl2Y6VK3mUCvonJLBvP8C2zO7fuy9zrplF975u0RdrRyrJMvWsa+SYms2wc8C2TLbtnODw7l/heuzPstFYId938Gj5x62VNq3Ya+SYms+yc0rFMms3i251ht7oZqfb3ZX/41cSLzIjEiQhX6Zi1yQHfMm3Y3WwWv5QZdpZWBXfecEniDVutGpaZWTKndIZc1n7tabP4FcsrXQu+tb70Y6OVk6+dWfGPrNlS+b+eIZd105G0vPsd77ikq+Ob/vHPOVKX1nnh6Jw3EDFbIgf8IZe1X3vjLlDjY6NdT61s2T3DZ3c+vaiM07X3ZkvjHP6Qa6dfe7dbEDe6e/v+Je+CZWaLeYY/5LLuFtWLjbibBXXX3pu1zzP8IZelX3vWO2yzVPu0I+23D7G0yiCzYeeAby1TNVk24l5q24VmknaLEvBbV0y4LNNsCRzwraUsC7tpF4UP378HWFrQ925RZvnqKOBLejlwH7AG+BHwnoh4IeG4eWBv9enTEXFDJ+e1YmVZ2E27KMxHdDTTL3qh2GyQdbpoOwV8PSIuAL5efZ5kNiIuq/5xsO8zWRZ2my2iuozSrBw6Dfg3AvdWH98L3NTh51kJZanBb9Xz3mWUZr3XaQ7/lRHxXPXx3wOvTDnuTEnTwHFgU0RsSftASRuADQATExMdDs+yalVh0yq1Unvvw/fvYT4WV8+7jNKs91oGfEkPAeckvPXR+icREZLS7pN5VUTMSHo1sEPS3oj4QdKBEbEZ2AwwOTmZ914ZliCPCpvaBWM+AsEpN0y5hbFZObQM+BFxddp7kn4i6dyIeE7SucDzKZ8xU/37KUnfANYBiQHfipel7LKZxgtGwMmgP+7KGrPS6DSHvxV4X/Xx+4CvNB4gaYWkM6qPVwLrge92eF7LUdZ+OmmSLhi1YP/I1FUO9mYl0WnA3wT8uqTvA1dXnyNpUtI91WN+FZiWtAd4mIUcvgN+iaTl17Pm3Tu9YJhZMTpatI2InwFvTXh9Gvhg9fHfApd2ch7rro3XXMjGL+xh7sRLmffKMqXm3RsXeM8erSTuTOWFWrNy8Z22tqBx9/CU3cSTFngrI6KyTKdcMLxQa1Y+DvjG3dv3Mzd/akHU3HycvFmqfjZ/9NjxRfn6uflgxfIKy08/zS0QzErMAd9Sc+218sz62Xyaw0fn2P2xt3VlfGaWDwf8AdaYa7/yolU8/OShRbPwtF45I9Ki2Xwa5+vNys8boAyopM3J/2rn04mblaf1ykm6YzaJ8/Vm/cEBf0Al1cY3qr+5KqlXznjKrH1stFLo3rZmlg+ndAZU1hr42nFpvXIaNyAZrYxw5w2XOMCb9SHP8AdU1px6s+OydMk0s/7hGf6AStoesFGW3Ls3IDEbHJ7hD6ik2fm/uGLCs3WzIeYZ/gBr3BP24ScP+YYosyHmgD/A8uhzb2aDwymdAdasz72ZDR8H/AHmtsVmVs8pnQ7Vty8YW14hAo7MzpWigVhaywS3QTAbTg74HWjMkb9w9KWe8M3y5a02DM9LUmlmq1LMosZmZsVzSqcDrdoXJOXLk3rc1Hra5O2mdePc8vpxRrTQ3H5E4pbXp9fVFzk2MyueA34HsuTCG48pciF1y+4ZHtg1c7IJ2nwED+yaSQ3gXuQ1G2wdBXxJ75a0T9IJSZNNjrtW0n5JByRNdXLOMsmSC288psiF1HYDuBd5zQZbpzP8x4GbgW+mHSBpBPgU8HbgYuC9ki7u8LylkNRWuF5SvrzTDcNb2bJ7hvWbdrB2alvqhiVpATzL2Oo/f/2mHU73mPWRjgJ+RDwREa1+378cOBART0XEMeDzwI2dnLcsGtsXrFheYWy00rR1QVrv+Tz6yTfm4NOkBfYrL1q1aCvb+rE5x2/W34qo0hkHnql7fhB4Q9rBkjYAGwAmJia6O7IctNtcrLHdQZ6VMFl64KddXGr5/voLheCURd5mKSJX8piVX8uAL+kh4JyEtz4aEV/Je0ARsRnYDDA5OZlty6U+k0cHyqTyyWa5dkHTi0tSMA/g4ScPnXzuHL9Zf2sZ8CPi6g7PMQOcX/f8vOprtkRpPXLGlldOuRegZnxslEemrmr6mVmCuW/kMutvRZRlPgpcIGmtpNOB24CtBZx3YKWlViJY8vpAlgXbbq4/mFn3dVqW+U5JB4E3Atskba++vlrSgwARcRy4HdgOPAHcHxH7Oht2ORVVwZI2Gz8yO7fkHaqyBHPvgGXW3xRR3jT55ORkTE9P93oYmTSmWWAhYHYjIK7ftCMxtZIlddOM2yqY9T9JuyIi8b4o99LJSZEVLEvpkVPTLKh7O0OzweaAn5Nmi555z5yzlnY2nvfKi1bxwK4Zb4hiNqSc0slJWpplxfIK/zh3opBUT72kFJMg8YasTlNBZlYezVI6bp6Wk7RFzwh60pAsra4+ievozYaDA35O0ipYjswurouH7gfZdj7fdfRmw8E5/BwlLXrevX1/T25WSrtJqjGt4zp6s+HhGX6X9epmpbTz/tYVE66jNxtSnuF3WV7N0tqt9OlmkzYz60+u0mmhDDcjJVXcVEbEWaefVpoN082sHHzj1RKlNSmDYuvWkypu5uaDw9UFYdfTm1kWzuE3UZY9XrNU3HjvWTNrxQG/ibL0f89a0eN6ejNrxgG/iW7vP5tVq71za1xPb2bNOOA3UZb+70l751aWnbr7rOvpzawVL9o2UabSxsabuspQPWRm/cVlmWZmA8TN08zMrOMtDt8taZ+kE5ISryjV434kaa+kxyR5ym5m1gOd5vAfB24G/izDsVdGxE87PJ+ZmS1RRwE/Ip4AkNTqUOsxL/KaWVE5/AC+KmmXpA3NDpS0QdK0pOlDhw4VNLzBVmsRMXN4luClVgxbds/0emhmVqCWAV/SQ5IeT/hzYxvneVNEvA54O/B7kt6cdmBEbI6IyYiYXLVqVRunsDRlaRFhZr3VMqUTEVd3epKImKn+/bykLwOXA9/s9HMtm7K0iDCz3up6SkfSWZJeVnsMvI2FxV4rSFlaRJhZb3ValvlOSQeBNwLbJG2vvr5a0oPVw14J/B9Je4C/A7ZFxP/q5LzWnrK0iDCz3uq0SufLwJcTXn8WuK76+CngtZ2cxzpTphYRZtY77qUzJJI2WDez4eLWCmZmQ8IB38xsSDjgm5kNiYHL4buFgJlZsoEK+LUWArW7SmstBAAHfTMbegOV0nELATOzdAMV8N1CwMws3UAFfLcQMDNLN1AB3y0EzMzSDdSirVsImJmlG6iAD24hYGaWZqBSOmZmls4B38xsSDjgm5kNCQd8M7Mh4YBvZjYkHPDNzIZEp3va3i3pSUnfkfRlSWMpx10rab+kA5KmOjmnmZktTacz/K8B/ywiXgN8D/hI4wGSRoBPAW8HLgbeK+niDs9rZmZt6ijgR8RXI+J49elO4LyEwy4HDkTEUxFxDPg8cGMn5zUzs/bleaftvwbuS3h9HHim7vlB4A1pHyJpA7ABYGJiIsfhdcYbq5hZv2sZ8CU9BJyT8NZHI+Ir1WM+ChwHPtvpgCJiM7AZYHJyMjr9vDx4YxUzGwQtA35EXN3sfUnvB34DeGtEJAXoGeD8uufnVV/rG802VnHAN7N+0WmVzrXAfwBuiIijKYc9Clwgaa2k04HbgK2dnLdo3ljFzAZBp1U6nwReBnxN0mOSPg0gabWkBwGqi7q3A9uBJ4D7I2Jfh+ctlDdWMbNB0NGibUT8SsrrzwLX1T1/EHiwk3P10sZrLjwlhw/eWMXM+s/A9cPvBm+sYmaDwAE/I2+sYmb9zr10zMyGhAO+mdmQcMA3MxsSDvhmZkPCAd/MbEgouRtCOUg6BPy4x8NYCfy0x2NIU+axgcfXiTKPDTy+TnVzfK+KiFVJb5Q64JeBpOmImOz1OJKUeWzg8XWizGMDj69TvRqfUzpmZkPCAd/MbEg44Le2udcDaKLMYwOPrxNlHht4fJ3qyficwzczGxKe4ZuZDQkHfDOzIeGA30DSuyXtk3RCUmrZlKQfSdpb3fhlumRju1bSfkkHJE0VMbbqeV8u6WuSvl/9e0XKcfPV7+0xSV3d/azVdyHpDEn3Vd//lqQ13RzPEsb3fkmH6r6vDxY4ts9Iel7S4ynvS9KfVsf+HUmvK2psGcf3FklH6r67jxU4tvMlPSzpu9X/Zv9dwjHFf38R4T91f4BfBS4EvgFMNjnuR8DKso0NGAF+ALwaOB3YA1xc0Pj+MzBVfTwF/FHKcb8oaDwtvwvg3wKfrj6+DbivwP8/s4zv/cAni/w5qzv3m4HXAY+nvH8d8DeAgCuAb5VsfG8B/mePvrtzgddVH78M+F7C/7eFf3+e4TeIiCciYn+vx5Ek49guBw5ExFMRcQz4PHBj90cH1fPcW318L3BTQedNk+W7qB/zF4G3SlKJxtczEfFN4OdNDrkR+MtYsBMYk3RuMaPLNL6eiYjnIuLb1cf/j4XtXRs31Cj8+3PAX7oAvippl6QNvR5MnXHgmbrnB1n8g9Ytr4yI56qP/x54ZcpxZ0qalrRTUjcvClm+i5PHxML+y0eAV3RxTInnrkr7/+qW6q/8X5R0fjFDy6SXP2tZvVHSHkl/I+mSXgygmiZcB3yr4a3Cv7+h3PFK0kPAOQlvfTQivpLxY94UETOS/ikLm7g/WZ1xlGFsXdNsfPVPIiIkpdX8vqr63b0a2CFpb0T8IO+xDoi/Bj4XES9K+h0Wfhu5qsdj6hffZuFn7ReSrgO2ABcUOQBJ/wR4APiDiPiHIs+dZCgDfkRcncNnzFT/fl7Sl1n49bzjgJ/D2GaA+lngedXXctFsfJJ+IunciHiu+qvp8ymfUfvunpL0DRZmP90I+Fm+i9oxByWdBpwN/KwLY0nScnwRUT+We1hYJymLrv6sdao+wEbEg5L+m6SVEVFIUzVJFRaC/Wcj4ksJhxT+/TmlswSSzpL0stpj4G1AYqVADzwKXCBpraTTWViI7GolTJ2twPuqj98HLPqNRNIKSWdUH68E1gPf7dJ4snwX9WN+F7AjqitqBWg5voac7g0s5ILLYivw29VqkyuAI3UpvZ6TdE5tPUbS5SzEu0Iu5tXz/jnwRET8ScphxX9/vVjBLvMf4J0s5NJeBH4CbK++vhp4sPr41SxUVOwB9rGQbinF2OKl1f/vsTBrLmRs1fO+Avg68H3gIeDl1dcngXuqj38N2Fv97vYCH+jymBZ9F8DHgRuqj88EvgAcAP4OeHXBP2+txndX9WdsD/AwcFGBY/sc8BwwV/25+wDwu8DvVt8X8Knq2PfSpKqtR+O7ve672wn8WoFjexML63zfAR6r/rmu19+fWyuYmQ0Jp3TMzIaEA76Z2ZBwwDczGxIO+GZmQ8IB38xsSDjgm5kNCQd8M7Mh8f8BGEQtliX2Q8YAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"smoothers = [UnivariateSplineSmoother(), UnivariateSplineSmoother()]\n",
"gam = GeneralizedAdditiveRegressor(smoothers)\n",
"gam.fit(X, y)\n",
"y_pred = gam.predict(X)\n",
"\n",
"plt.scatter(y, y_pred); # plot goodness of fit"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAssAAAFlCAYAAAAd9qXYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3hURdvH8e+k94QSIIUeei+CShMLRaTbELHQRMXu86ivXR97wa5UQewCKlJEEKWD9BJqCC0hkJCQ3jbZef+YqJQEAmT3bJL7c117mew5Z88PhHDv7Mw9SmuNEEIIIYQQ4mxuVgcQQgghhBDCVUmxLIQQQgghRAmkWBZCCCGEEKIEUiwLIYQQQghRAimWhRBCCCGEKIEUy0IIIYQQQpTAw+oA51K9enVdr149q2MIIcQF27hx4wmtdajVOZxJfmYLIcqrc/3MduliuV69emzYsMHqGEIIccGUUoeszuBs8jNbCFFenetntkzDEEIIIYQQogRSLAshhBBCCFECKZaFEEIIIYQogRTLQgghhBBClECKZSGEEEIIIUpQJsWyUmqaUipRKbWjhONKKfWBUipGKbVNKdW+LO4rhBBCCCGEI5XVyPJ0oM85jvcFGhU9xgKfltF9hRBCCCGEcJgyKZa11suBlHOcMhD4QhtrgRClVFhZ3FsIIYQQQghHcdac5QjgyCnfxxU9dxal1Fil1Aal1IakpCSnhBNCCCGEEKI4LrfAT2s9SWvdUWvdMTS0Uu0UK4QQQgghXIyziuV4oPYp30cWPSeEEEIIIYTL8nDSfeYC45VS3wKdgTStdYKT7i2EkZ8PMTFw+DBobZ4LCIAmTSA0FJSyNp8QQgghLlpCWg77jmdyeYNqeHmU3XhwmRTLSqlvgKuA6kqpOOB5wBNAa/0ZsAC4HogBsoG7y+K+QpxTYSGsXg0//QQLF8Levea54lStCp07w6BBMHAg1Kzp3KxCCCGEuCS/bD3Kqwt2s/rJqwkP8S2z1y2TYllrPew8xzVwf1ncS4jzSk+HiRPhvffg6FHw8oKrr4YhQ6BpU6hfHzyK/uifPAl79sDOnbBkCdxzD4wbBzfcAE8+CVdeae2vRQghhBClsnR3Ik1rBZZpoQzOm4YhhOPl5sJbb8G770JqKlxzjfm6b18ICir5uj5FLcK1hh074Pvv4dNPoUsX6N4d3nkHOnZ0zq9BCCGEEBcsPdfGhoMnGdO9QZm/tst1wxDiovzxB7RpA889B1ddBX/9ZUaKb7nl3IXyqZSCVq3g5Zfh0CEzMr13r5me8eCDZsRaCCGEEC5n5b4TFNg1PZvUKPPXlmJZlG+5uXDvvWaaRWEhLFoEP/4Il112aa/r7w8PPQS7d8N998FHH0Hz5rByZdnkFkIIIUSZ+WN3IkE+HrSvE1Lmry3Fsii/Dh6Ebt3gs8/gscdg+3bo1ats7xEcDB9+CGvXgq+vGbWeMOHfbhpClENKqWlKqUSl1I4Sjiul1AdKqRil1DalVHtnZxRCiNKy2zV/7k2ie+NQPNzLvrSVYlmUT8uWQYcOZprEjz/C22+bYtZROnWCDRugf3949FEYNgzy8hx3PyEcazrQ5xzH+wKNih5jgU+dkEkIIS5K9NF0kjLyHDIFA6RYFuXRzz9D795Qo4YpYAcNcs59g4Nhzhx47TX47jvo1w8yMpxzbyHKkNZ6OZByjlMGAl9oYy0QopQKc046IYS4MH/sSUQpuCZlH9xxh+mEVYakWBbly/TpMHSoWcy3ciU0auTc+ytlWsrNmAF//mnmSp844dwMQjheBHDklO/jip47i1JqrFJqg1JqQ1JSklPCCSHEqf7Yk0jryBBCfvwBZs82g1tlSIplUX58/jncfbcpUH//HapVsy7LHXeY6R87dpg8KecapBOi4tJaT9Jad9RadwwNDbU6jhCikknJymfLkVSuiaoKs2aZfRL8/cv0HlIsi/Jh9mwYPdos4PvlF7NNtdX69zdZ9uwxUzIyM61OJERZiQdqn/J9ZNFzQgjhUv7YnYjWcEPKHkhKMi1jy5gUy8L1LVpkFtRdfrmZM+ztbXWif117LXz7renrPHiwLPoTFcVc4I6irhiXA2la6wSrQwkhxJmW7DpOzSBv6i+dbwbS+vYt83tIsSxc29atZo5ys2Ywb16Zf7RSJgYPhmnTzCYoo0ZJWznh8pRS3wBrgCZKqTil1Cil1Dil1LiiUxYAsUAMMBm4z6KoQghRoryCQpbvTeK6RlVRc+bAgAEO6Ywl210L15WYaP7gh4TAwoVQpYrViUp2550QHw9PPw0tW5pFgEK4KK31sPMc18D9ToojhBAXZW1sCln5hdyUutesHXLAFAyQYlm4qrw8GDLEzD9asQLCw61OdH5PPWUW/P3f/5mR8IEDrU4khBBCVFhLdh7H19OdFit/haAg01bWAWQahnBNDzwAq1aZVnEdOlidpnSUgqlTTd7hwyE62upEQgghRIWktWbJruP0bBCMx88/mT0XHLSmSYpl4Xq+/BImTzYjtTffbHWaC+PrazZN8fc32bOzrU4khBBCVDjRR9NJSMvltpO7IC3NYVMwQIpl4Wr27oVx46BbN3jpJavTXJzwcFPw79oFDz5odRohhBCiwvl9l9m177LVv5p9F667zmH3kmJZuI7cXDMa6+MDX38NHuV4Sv1115mR8alT4auvrE4jhBBCVChLdh3nyhpeeC+YZ2oHT0+H3UuKZeE6nnjCtIqbMQMiI61Oc+lefBG6dIF77oHYWKvTCCGEEBXC0dQctsencfeJbZCTY9YJOZAUy8I1/PEHfPCBWdjXr5/VacqGh4cZIXd3h5EjwW63OpEQQghR7v0WfQyAK9cugrp14YorHHo/KZaF9dLT4e67oVEjeP11q9OUrTp1YMIEWLYMPvzQ6jRCCCFEubco+jiX+drwW7bU7PDr5thyVoplYb3HH4cjR0ybOD8/q9OUvbvvNqPlTz1lFjAKIYQQ4qKkZOWz7kAy45I2QWEh3Habw+8pxbKw1uLFpk3c44/DlVdancYxlDK/Rh8fUzjLdAwhhBDioizZdRy7hivWLoJWrczDwcpxuwFR7uXkmDZxjRubxXBOlJVXwN7jGRxKziYxI5fj6Xlk5RVgK9QU2u34eLoT6ONBsK8nkVX8qFPNjwbV/Qnx87q4G4aFwXvvmW2xJ082i/6EEEIIcUF+iz5Gp8KT+G38C157zSn3lGJZWOd//zNdIpYuNaOuDnQkJZs1+5NZG5vMhkMnOZxy+mYhPp5uBPl44unuhrubIsdWSEaujVzb6aPAdav50a52CB3rVeXqpjUID/EtfYgRI0ynjyeeMFth16pVFr80IYQQolLIzCtg+b4TTI5bYz61dXAXjL9JsSysER0Nb75pRlp79nTILQ6cyGL+tqPM336MXQnpAFTz9+KyelW5sUMkTWoF0jDUn5pBPgR4e6CUOus1cm2FxJ3M4VByFvsSM9lyOJXV+5P5actRAJqFBdG3ZS2Gdogk4nyFs1Lw6afmI6NHHoFvvinzX7MQQghRUS3bk0S+rZDOq+bD1VdD7dpOua8Uy8L57HYz/SI4GN5+u0xfOr/AzqLoY3y17hBrY1MA6Fi3Cs/0a0b3xqE0qhFQbFFcEh9Pd6JqBBBVI4BrmtUEzH70+5OyWLr7OEt2JvLu4r1MWLKXLg2rc8cVdbm2WU3c3Eq4R+PG8PTT8PzzcNdd0Lv3pf6ShRBCiEphUfQxrkneh8+RQ/DKy067rxTLwvm++gpWrjS721WvXiYvmZ1fwNfrDjNpeSyJGXnUrurLE32aMrhdBLWCy3aKh1LqnwJ6bPeGHEnJZvamOH7YEMfYmRtpEOrPPd0bMKR9JJ7uxayhfeIJ83vwwAOwYwd4XeQ8aCGEEKKSyLUVsnR3IlMPrQJ/fxg82Gn3lmJZOFdGhikWO3UyI6uXKNdWyIzVB/ls2X5OZtu4okE13ryxNd0bhZY8ulvGalf14+FrGzO+ZxQLdxxj4vL9PDF7O5/+uZ/HezehX6uw00ezvb3h/fehb1+zEcvjjzslpxBCCFFeLd+bhC0ziw7rFsPQoRAQ4LR7S7EsnOu11yAhAebMuaQm4lpr5m49ypu/7iE+NYcejUN58JpGdKhbpQzDXhgPdzf6twnnhtZh/LEnkTcW7mH815uZXPsArwxqScuI4H9P7tPH9F5+6SW4/XZZ7CeEEEKcw4LtCQw8shGPjHS44w6n3lv6LAvniY2Fd94xXSEuv/yiXyYmMZNbJq3loW+3EOzryVejOzNjZCdLC+VTKaW4umlNFjzUjTdvbE38yRwGfLSSl+ftJCuv4N8T330XcnPNHGYhhBBCFCvXVsiSXYncHbPCLOpzUGOAkkixLJznscfA0/Oit7TOL7Dz3pK9XP/+CnYnpPPakFbMe6ArXaLKZt5zWXN3U9zcsTa/P9aDYZ3qMHXlAXpNWM5fB8zCQxo3hocegs8/hw0brA0rhBBCuKjle5PwP3GMpltXmQE3B29vfSYploVzrFgBP/1ktnwOD7/gy2MSMxny6SreW7KPPi1r8ftjVzGsUx2nzUu+FMG+nrwyuBWzxl2Bh7vi1klreGvRbmyFdnj2WQgNNfOWtbY6qhBCCOFyFmxPYPieZSi73eyE62RSLAvH0xr++1+IiDD9hS/oUs3MtYe44cMVxJ/M4bPbO/DBsHaEBno7KKzjdKxXlfkPdmNo+0g+/mM/N09cw3G84LnnYNkyWLDA6ohCCCGES8m1FbJk53Fui14CPXpAVJTTM0ixLBxvzhxYu9YsZvPzK/VlWXkFPPTtFp79aQed6ldj0cPd6dOyfC+EC/D24K2b2vDRbe3YcyyDfh+sZEOvG81f/iefhMJCqyMKIYQQLmPFvhM037+V6seOwMiRlmSQYlk4ls1mpl60aGF26yul/UmZDP5kFfO2HeU/vZsw/a7LqBHk2C2xnemG1uH8dH8XAn08uPXzjawc+ajpufzFF1ZHE0IIIVzG/G1HGR79Ozow0LSMs4AUy8KxJk+GffvMoj5391JdsnLfCQZ9vIoTmfl8MbIz9/eMKhdzky9U45qB/Dy+C10bVef21NocbdIa/eyzkJ1tdTQhhBDCcjn5hazaHEvf3StQw4aZzUgsIMWycJzsbHj5Zeje3fQULoWv1x3mzs//IjzYl7lFhWRFFuTjyZQ7OjKscx0eaT8MFR9PwcefWB1LCCGEsNzvu49zzdY/8crLtWwKBkixLBzpk0/g2DF45RVQ5x4Z1lrzxq+7+b8ft9M1qjqz7r2CyCqln99cnnm4u/Hq4FZ0H3Mjy+u1I+elV8hNSbU6lhBCCGGpuVuOMmLHYnSLFmbnX4tIsSwcIyPDTL3o3Ru6dj3nqQWFdp4s2h76ts51mHpnRwJ9PJ0U1DUopbi/ZxQ5zzxPYGYqP496kuz8gvNfKIQQQlRAaTk2ji9bS4v4PaixY8876OZIUiwLx/jgA0hONh0wziGvoJDxX2/muw1HeODqKF4Z1BIP98r7x7L3qIEc7XYtfRZ+yb0f/X76jn9CCCFEJbFoxzFu3jQfu7eP2YjEQpW3KhGOc/IkvPUW9O9/zo9Ncm2FjP1iI79GH+PZG5rzWK8mKAvfObqK8A/eIjgviw5zZjB6xgZybdJOTgghROXy27p9DN65DHXrLVCliqVZpFgWZe+99yAt7Zyjyrm2QsbO3MiyvUm8PqQVo7rWd2JAF9e2Ldx4I/du+YVdOw9y75cbyS+wW51KCCGEcIrEjFxC5/+IX34O6p57rI4jxbIoY2lp8P77MHiwKfqKkWsrZMwXG1ixL4k3h7bm1k51nByyHHj+eTyzMpmZuZY/9iTx8HebKbTLdthCCCEqvvnbEhi25VfymjaHyy+3Oo4Uy6KMffihKZifeabYw7ZCO+O/3syKfSd4Y2hrbr6stpMDlhMtW8LgwbSaNZ0Xe0SyYPsxXp63E62lYBZCCFGxbft5Ka2PxeB9/72WLuz7mxTLouxkZMCECaancvv2Zx222zX/+WErS3Yd5+WBLbi5oxTK5/TMM5CWxp2b5jGqa32mrz7IlBUHrE4lhBBCOExMYgaX//odNh9fuP12q+MAUiyLsvTJJ5CSAs8+e9YhrTXPz43mpy1m++oRV9Rzfr7ypn17uOEGePddnu4WSb9WYbyyYBe/bD1qdTIhhBDCIRYsi2bgrmUUDLsNQkKsjgNIsSzKSlYWvPMO9OoFnTufdfiTP/czc+0h7unRgPt7RlkQsJx69llIScFt4me8c3MbOtatwuM/bGV7XJrVyYQQQogyZbdr+Hw6PgX5+D78oNVx/iHFsigbU6dCUlKxc5Vnb4zjrUV7GNwugif7NLUgXDnWqZN5A/LOO/gU5PPZiA5UD/BmzBcbSEzPtTqdEEIIUWbWxiQxcM3PJLfvDK1bWx3nH1Isi0tns5lR5a5doVu30w6t3HeCJ2Zvo0tUNd4Y2lr6KF+Mp56CxESYMYPqAd5MuqMDaTk2xs7cKD2YhRBCVBg7P/+euqnHCHz0IaujnEaKZXHpvv0WDh+GJ5887emYxEzu/WojUTUC+PT2Dnh5yB+3i9Kjh5na8tZbUFBAi/Bg3r25DVuOpPLiL9FWpxNCCCEuWU5+IY1nzSC9SiheNw21Os5pyqR6UUr1UUrtUUrFKKWeLOb4XUqpJKXUlqLH6LK4r3ABdju88YZpdXb99f88nZqdz+gZ6/H2cGPKnR0J8vG0MGQ5p5R5IxIbC7NmAdC3VRj3XtWQb/46wqyNcRYHFEIIIS7N8oWr6RqzkfQRd4GXl9VxTnPJxbJSyh34GOgLNAeGKaWaF3Pqd1rrtkWPKZd6X+Ei5s+H6Gh44ol/eiHaCu3c//UmjqbmMnFEByKr+FkcsgIYMACaNjVvTIp6LT92XWMub1CVZ37azu5j6RYHFEIIIS6e/f0PKHR3J+KJh62OcpayGFnuBMRorWO11vnAt8DAMnhdUR688QbUrQu33PLPU6/M38WqmGReG9KKDnWrWhiuAnFzg//+F7Zsgd9+A8DD3Y0PhrUjyMeTe7/cREauzeKQQgghxIU7HBtP95Xz2H91P1R4uNVxzlIWxXIEcOSU7+OKnjvTUKXUNqXULKWU7EZREaxeDatWwWOPgaeZZvHj5jimrz7IqK71Gdoh0uKAFczw4RARAW+++c9TNQJ9+HBYOw4lZ/H8XJm/LIQQovyJff0D/G25VHv6v1ZHKZazVlz9AtTTWrcGFgMzSjpRKTVWKbVBKbUhKSnJSfHERXnnHahSBUaOBGDn0XSemrOdzvWr8lRfaRFX5ry84MEHYelSM8JcpHODajxwdSPmbIrn5y3xFgYUQgghLkxhvo0mP8xgd+N2hPa40uo4xSqLYjkeOHWkOLLouX9orZO11nlF304BOpT0YlrrSVrrjlrrjqGhoWUQTzjE/v3w449w773g709ato1xX24kxNeLj25rj4e7dL5wiLFjISDAvFE5xQNXR9GhbhWe+XEHR1KyLQonhBBCXJjdn35BWOpxsu4bb3WUEpVFRbMeaKSUqq+U8gJuBeaeeoJSKuyUbwcAu8rgvsJKEyaYqRfjx6O15j+ztpKQlsPHw9sTGuhtdbqKKyQERo0y7fri/u2C4eHuxnu3tAXg4e+2UGjXViUUQgghSs3nkw85UiWMVuNGWB2lRJdcLGutC4DxwCJMEfy91jpaKfWSUmpA0WkPKqWilVJbgQeBuy71vsJCKSnw+edw220QFsb01Qf5bedxnujTlA51q1idruJ76CHTsu/DD097unZVP14e1JKNh04yaXmsReGEEEKI0kn7fRkN924leuideHm7bovZMvmsXGu9QGvdWGvdUGv9StFzz2mt5xZ9/ZTWuoXWuo3WuqfWendZ3FdY5LPPIDsbHn2UbXGpvLpgF9c2q8GorvWtTlY51K8PQ4fCxImQkXHaoYFtw+nbshYTFu9lz7GMEl5ACCGEsF7Ki6+S6hNAo/970Ooo5yQTS8WFyc+Hjz6CXr3IaNSU8V9vJjTAm7dvaiNbWTvTY49BWhpMm3ba00op/jeoJYE+Hjz2wxZshXaLAgohhBAls+/aTd2Vi/m95400rB92/gssJMWyuDDffw8JCfDIIzw/N5r41Bw+GNaOED/X2m2nwuvcGa680kzFKCw87VC1AG9eGdySHfHpfPxHjEUBhRBCiJIlvvAqNjcPfB9zvU1IziTFsig9reH996FpU36p2ZI5m+IZ3zOKjvVk4xFLPPSQ6UqyYMFZh/q0DGNAm3A+/iOGvcdlOoYQQggXcuwY1eZ8y9x2vbi6e0ur05yXFMui9NasgQ0bSB09jqd/jqZdnRAeuDrK6lSV1+DBEBlp3sAU4/n+zQnw9uDJ2duwS3cMIYQQLiLr7XdxLyggccz9+Hi6Wx3nvKRYFqX3/vvokBAe8mpFoV3z/i3tpJ+ylYpa9/H777B9+1mHqwV48+wNzdl0OJUv1x2yIKAQQghxhtRUPD77lIVNrqTP4G5WpykVqXRE6Rw5ArNns6PPTSyLz+H5AS2oU83P6lRizBjw9YUPPij28OB2EXRrVJ03f91DQlqOk8MJIYQQp7N/9BHeWZmsuGksDUMDrI5TKlIsi9L5+GO01jxY/UquaVqDmzpEWp1IAFStCiNGwJdfwokTZx1WSvHq4FYU2O28MDfagoBCCCFEkcxMCt55lyUNL+OqW66zOk2pSbEszi8nBz1lCutadSWlejivDWklbeJcyYMPQm4uTJ1a7OHaVf148JpGLIo+zh97Ep0cTgghhCgycSJeqSf57roRXNusptVpSk2KZXF+332HSk7m/WZ9eGlgC2oE+VidSJyqRQvo2RM++eSsNnJ/G921AQ1C/XlhbjS5tuLPEUIIIRwmN5eCN99iZd02tL25b7la81R+kgpraE3ehPfZG1qX4L7XMaBNuNWJRHHGj4fDh2HevGIPe3m48dKAlhxKzpatsIUQQjjf1Kl4JB7ns27DGNapjtVpLogUy+Kc7KvX4L1tC993GsBLg1vK9AtXNWAA1K5tdlcsQddG1enXOoyP/4jhSEq2E8MJIYSo1HJzsb/6Ghtrt6BW/95U9S9fG5lJsSzO6dBLb5Lu5Ufz/95HjUCZfuGyPDzg3nthyRLYtavE057t1xx3N8Ur80s+RwghhChTkyfjdjSet7sM564u9a1Oc8GkWBYlOrbnAJFL5rO6e38Gd2tidRxxPqNHg5eXmbtcglrBPtx3VUN+jT7Gmv3JTgwnhBCiUsrORr/6KpsatMHWvQctI4KtTnTBpFgWxdJas/apN/C0F9Dmladk+kV5EBoKt94K06dDRslbXI/u1oCIEF9emreTQtnZr9JSSvVRSu1RSsUopZ4s5vhdSqkkpdSWosdoK3IKIcq5zz5DHTvGa5cPY3S3BlanuShSLIti/bY1ns6LZ3GkYxfCOrWxOo4orfvug8xM03e5BD6e7jzZtym7EtL5YcMRJ4YTrkIp5Q58DPQFmgPDlFLNizn1O61126LHFKeGFEKUf1lZ8PrrbGl6GYltO3Fd8/LTLu5UUiyLs2TmFbD0nWmEZSYT9uQjVscRF6JTJ2jXDj79FHTJo8Y3tA6jY90qvP3bHjJybU4MKFxEJyBGax2rtc4HvgUGWpxJCFHRfPABJCXxUsebGdmlPu5u5fNTaimWxVkmLN5Lv1U/kV8rDI+B8u9nuaKUGV3evh1Wrz7HaYrn+jfnRGY+E5dJK7lKKAI49WOFuKLnzjRUKbVNKTVLKVXbOdGEEBVCSgq88QZb23Vjf1RrbupYfnf+lWJZnCb6aBp//LKS7gc243XvONNlQZQvw4ZBcPA5F/oBtI4MoX+bcKasjOV4eq6Twoly5Begnta6NbAYmFHcSUqpsUqpDUqpDUlJSU4NKIRwYa+9hk5P57/tb2F45zr4eZXfekKKZfEPrTXP/RzN3TsWoz08THcFUf74+8Odd8KsWZB47u2t/9OrCYV2zYTFe50UTriIeODUkeLIouf+obVO1lrnFX07BehQ3AtprSdprTtqrTuGhoY6JKwQopyJi4MPP2RLjxs4ULMBd15Zz+pEl0SKZfGPOZvi2RFzjJt3LEYNGgThsltfuTVuHOTnw7Rp5zytTjU/br+8Lt9vOMK+4yV30BAVznqgkVKqvlLKC7gVmHvqCUqpsFO+HQBIc24hROm88AJaax5pOYShHSKoGVS+92mQYlkAkJ5r47WFu7k3aRPe6WlmgwtRfjVrBj17wsSJYLef89QHrm6Ev5cHb/y620nhhNW01gXAeGARpgj+XmsdrZR6SSk1oOi0B5VS0UqprcCDwF3WpBVClCs7d8Lnn7P++ls5HBDKPd0bWp3okkmxLAB4b/E+krPyGLNzMTRubAotUb7dcw8cPAi//XbO06r6ezHuqoYs2ZXIxkMnnZNNWE5rvUBr3Vhr3VBr/UrRc89precWff2U1rqF1rqN1rqn1lreTQkhzu/xx9GBgTwa1Y9+rcOpV93f6kSXTIplwb7jGcxYc5CHa+biv/EvU2TJJiTl3+DBZqOSiRPPe+rdXepRPcCLtxftcUIwIYQQFdKiRbBwIatuHUecuz/39ij/o8ogxXKlp7XmpXk78fdy555dS8Db2ywOE+WflxeMHAm//AJHj57zVD8vD+67Koo1scmsjjnhpIBCCCEqjIICeOwx7PUb8FjNbvRsEkrz8CCrU5UJKZYruT/2JLJi3wke6xKJz7dfw403QrVqVscSZWXMGCgshKlTz3vqbZ3rEBbsw1u/7UGfY0MTIYQQ4izTpkF0NL/f/RjH8+D+nlFWJyozUixXYvkFdv43bxcNQv0ZfnA1pKebKRii4mjYEK67DiZPNkXzOfh4uvPgNY3YfDiVpbvP3XJOCCGE+EdqKjz7LIVdu/IUUXRrVJ2O9apanarMSLFciX2x5iCxJ7J4tl9zPCZPhubNoWtXq2OJsjZuHBw5AgsXnvfUGztEUreaH+/8tldGl4UQQpTOCy9AUhI/3/kfTmTZeOiaRlYnKlNSLFdSJ7Py+eD3ffRoHErP3KOwfj2MHSsL+yqi/v2hVq1SLfTzdHfjgasbsTMhnSW7ZHRZCCHEeWzfDh99hG30GF495lvhRpVBihhjtxUAACAASURBVOVK68OlMWTmFfB0v2YwaRL4+MCIEVbHEo7g6WkW+i1YYHZVOo9BbcOpU9WP93+X0WUhhBDnoDU88AAEB/NN/7GcyMyvcKPKIMVypXTwRBYz1x7klstq0zjADb76yizsq1qx3gmKU4waZTYn+fzz857q4e7G+J5R7IhP5489MroshBCiBN99B8uWkffiS7y/JaVCjiqDFMuV0puLduPp7sYj1zaG7783C/vGjrU6lnCkBg3MQr8pU8670A9gcPsIalf15f0l+2R0WQghxNnS0uDRR6FdOyY1vobkrHweva6x1akcQorlSmbjoZMs2H6Msd0bUCPIx3RJaNpUFvZVBmPHwuHD593RD8zc5fuvimJrXBp/7k1yQjghhBDlyv/9Hxw/Ttp7HzFx1SF6t6hJuzpVrE7lEFIsVyJaa15bsIvQQG/GdGsAO3bAmjWmF68s7Kv4BgwwO/pNnlyq04e0jyQixJePlsbI6LIQQoh/rV0Ln34K48fzYXoI2fkF/Kd3E6tTOYwUy5XI77sS2XDoJA9f2wh/bw9TNHl5wR13WB1NOIOXF9x9N8ydCwkJ5z/dw42x3Ruw8dBJ/jqQ4oSAQgghXJ7NZj6pDA/n6ONP88XaQ9zYIZKoGoFWJ3MYKZYriUK75s1Fu6lf3Z+bO9aG3FyYOROGDIHq1a2OJ5xl9GgzZ3n69FKdfnPH2lTz9+KTP/c7NpcQQojy4Z13/mkXN2GtGXh56NqKOVf5b1IsVxI/bY5n7/FMHuvVGE93N5gzB06eNFMwROXRqBFcdZXZ/tpuP+/pvl7ujOxan2V7k9gRn+b4fEIIIVzXzp3w/PMwZAjRnXoya1Mcd15Rl4gQX6uTOZQUy5VAXkEh7y7eS6uIYK5vGWaenDLFdEi46ipLswkLjB4N+/fDsmWlOn3EFXUJ9PbgUxldFkKIyqugwEzlCwxEf/wx/5u3ixBfT8ZfXfH6Kp9JiuVK4Ku1h4lPzeGJPk1xc1MQEwN//GF677rJH4FKZ8gQCAkp9UK/IB9PRlxRlwU7EohNynRwOCGEEC5pwgT46y/46COWpCjWxCbzyHWNCfb1tDqZw0mlVMFl5xfwyZ8xXNGgGl0bFc1NnjoV3N3hrrsszSYs4usLt98Os2dDcnKpLhnZtT6e7m5MXnHAweGEEEK4nN274dlnYfBg8ofexKsLdtEw1J9hnepYncwppFiu4GasPsSJzHwe7100+d5mM4u7+vWD8HBLswkLjRkD+flm98ZSqB7gzdD2kczeFMeJzDwHhxNCCOEy8vNh+HAICIBPPmHmusMcOJHFM/2amzVQlUDl+FVWUum5Nj5btp+eTULpULdo+8n58+HYMTNvVVRerVvDZZeZqRil7KE8ult98gvsfLHmkIPDCSGEcBkvvgibNsHkyST6h/De4r10bxzKVU1CrU7mNFIsV2DTVh4gLcfGo9ed0ih8yhQIC4O+fa0LJlzD6NFmY5q//irV6Q1DA7i2WU1mrjlITv75t8wWQghRzq1cCa+/DiNHwuDBvL5gN3kFdl4c0AJViTYzk2K5gkrNzmfqigP0aVGLVpHB5sn4eFi40Kxm9fCwNqCw3q23gp+fmcNeSvf0aMDJbBuzNsU5MJgQQgjLpaXBiBFQrx689x5/HUhhzuZ4xnSvT/3q/lancyopliuoyStiycwv4JHrTmkUPn266a07cqRluYQLCQqCm2+Gb76BzNJ1uehYtwpta4cwZUUshXbZAlsIISokrc2nj0eOwJdfUuDnz3M/7yAixJf7e0ZZnc7ppFiugE5m5TN91UH6tQqjSa2i7Sftdpg2DXr2hIYNrQ0oXMeoUaZQ/uGHUp2ulGJMtwYcSs5mya7jDg4nhBDCEp9+CrNmwauvwhVXMH31QXYfy+DZG5rh51X5PpmWYrkCmrIylmxbIQ9ec0qj8D//hNhYUxwJ8bcuXaBJkwuaitG7RU0iQnz5fJW0kRNCiApn82Z45BGztunxxzmcnM3bv+3h2mY16N2iltXpLCHFcgVz6qhy45qB/x6YOtVsRDFkiHXhhOtRyryBWrXK9NEsBQ93N+64oi5rY1PYeTTdwQGFEEI4TWqqmZ4XGgpffIFWiqd/2o6HmxsvD2pZqRb1nUqK5Qqm2FHlkyfNBhTDh5sNKYQ41R13mAWfFzC6fOtldfD1dJfRZSGEqCjsdrNh1cGD8N13UL06czbFs2LfCZ7o04Sw4MpbP0ixXIGkZpcwqvzVV5CXJ72VRfFq1oT+/eGLL8ymNaUQ7OfJ0A4R/Lz1qGxSIoQQFcGLL5q9GN5/H7p0ISkjj5fn76RD3SoM71zX6nSWkmK5Apm26iBZ+WeMKoMZMWzXDtq2tSaYcH2jRkFiIsybV+pL7rrSbFLy9brDDgwmhBDC4ebOhZdegrvugnvvRWvNU3O2kZ1fyBtDW+HmVjmnX/xNiuUKIj3XxvRVpq/yaaPKmzfDli2ysE+cW+/eZvvzC5iKEVUjgB6NQ5m59hD5BXYHhhNCCOEw27aZaZodOpguGErxw8Y4luxK5L+9mxBVI/D8r1HBlUmxrJTqo5Tao5SKUUo9Wcxxb6XUd0XH1yml6pXFfcW/Zq45RHpuAeOvPqP/4dSp4O0Nt91mTTBRPnh4wJ13mk1r4uNLfdldV9YjKSOPRdHHHBhOCCGEQxw7BjfcYPru//wz+PgQdzKbl37ZSef6VRnZpb7VCV3CJRfLSil34GOgL9AcGKaUan7GaaOAk1rrKGAC8Mal3lf8Kzu/gKkrD9CzSSgtI4L/PZCTY+YrDx0KVapYF1CUDyNHmgUeM2aU+pIejUOpU9WPmWsOOTCYEEKIMpeTAwMHQnIy/PILRERQaNf854dtaK15+6Y2lX76xd/KYmS5ExCjtY7VWucD3wIDzzhnIPD3v8CzgGtUZe0/4gBfrztMSlY+468+Y67yjz+aNjAyBUOURlQU9OhhNq+xl25ahZubYsTldfnrYAq7EqSNnBBClAuFhabzxfr1ZlCtfXsAPlu2nzWxyTzfvwW1q/pZHNJ1lEWxHAEcOeX7uKLnij1Ha10ApAHVyuDelV5eQSGTlsdyRYNqdKh7xujx1KlQvz5cdZUl2UQ5NGoU7N8PK1aU+pKbOkbi7eHGFzK6LIQQrk9rGD8e5syBCRNg0CAANhxM4d3Fe+nfJpybOkZaHNK1uNwCP6XUWKXUBqXUhqSkJKvjuLw5m+JJzMg7e6/2Awdg6VK4+25wc7n/zcJVDR1q5q5dwEK/ED8vBrYN56fN8aTllK71nBBCCIu8/DJ89hk88QQ89BAAadk2Hvp2CxEhvrwyuPJuPlKSsqii4oHap3wfWfRcsecopTyAYCC5uBfTWk/SWnfUWncMDQ0tg3gVV6FdM3HZflpHBtMl6oyB+s8/N7uz3XWXJdlEOeXnZxaDzpoFaWmlvuyOK+qRYytk9sY4B4YTQghxST7+GJ5/3izofu01AOx2zeOztnI8PZcPh7UjyMfT4pCupyyK5fVAI6VUfaWUF3ArMPeMc+YCdxZ9fSOwVGuty+DeldrCHQkcTM7m3h4NT38XWFhoiuXevaF27ZJfQIjijBplFn58802pL2kZEUy7OiF8ue4Q8ldbCCFc0JQpZvrFgAEwebIZUAM+XbafxTuP89T1zWhTO8TikK7pkovlojnI44FFwC7ge611tFLqJaXUgKLTpgLVlFIxwKPAWe3lxIXRWvPJH/tpEOpP7xa1Tj+4eDHExcnCPnFxOnSAVq3MQr8LMLxzXWKTslgbm+KgYEIIIS7Kl1/C2LHQpw98/z14mtHj5XuTePu3PQxoE87ILvWszejCymQyq9Z6gda6sda6odb6laLnntNazy36OldrfZPWOkpr3UlrHVsW963Mlu87wc6EdMZ1b3h2a5dp06BaNbOFsRAXSinzRmv9eti+vdSX3dA6jCAfD77+S3b0E0IIlzFzppl2cdVVZlGftzcAR1KyefDbzTSpGcjrQ1vJPOVzkJVf5dRnf+6nVpAPg9qd0XjkxAn46ScYMeKfvxBCXLDbbwcvrwsaXfbxdGdI+0h+3ZFAcmaeA8MJIYQolcmT/y2Uf/kFfH0ByMi1MXrGBux2zWe3d8DPy8PanC5OiuVyaFtcKmtikxnVtT5eHmf8L/zyS7DZZAqGuDTVqplm9TNnQl7pC9/hnetgK9TMkoV+QghhrQ8++Hfqxbx54O8PQEGhnfFfb2Z/Uiaf3t6BetX9LQ7q+qRYLocmLo8l0NuDWzudsXhPa9Pyq1MnaNnSmnCi4hg1yuzsNPfM9bola1QzkMvqVeGbvw5jt8tCPyGEcDqt4emnTVu4QYPMBmVFI8paa178ZSfL9ibx8qCWdImqbnHY8kGK5XLmcHI2C7cncNvldQg8s73L+vWwY4fZtliIS3XttaabygX0XAa4rXMdDiZnsya22O6QQgghHOXvT5ZffRVGj4YffjhtSuak5bHMXHuIsd0bMKxTHQuDli9SLJczU1fG4u6mGNml/tkHp00z7x5vvdX5wUTF4+5u+nT/9hscLv2ivb4twwjx85SFfkII4UxpaaYt3Oefm17KkyaBx79zkb/fcITXFu6mX+swnujT1MKg5Y8Uy+XIyax8vt8Qx8C2EdQM8jn9YHa26Yt7000QHGxNQFHx3H23+Uhv+vRSX+Lj6c6gthEsjj7Oyax8x2UTQghhxMTA5ZfDkiVmUd8LL/zTRxngt+hjPDl7G90aVWfCzW1xP7OLljgnKZbLkS/XHiLHVsjY7g3OPjhrFqSny8I+Ubbq14drrjEjFXZ7qS+75bLa5Bfa+XHzmZt5CiGEKFNLlkDnzpCYaPZZGD36tMPL9yYx/pvNtIoM4bPbO5zdGECcl/yOlRN5BYXMWHOIHo1DaVwz8OwTpk6FRo2gWzfnhxMV26hRcPAgLF1a6kuahQXRJjKY79YfkR39hBDCEex2eOUV6NULatWCv/4yLeJOsWJfEmO+2EDD0ACm33UZ/t7SIu5iSLFcTszdcpQTmXmM7lbMXOV9+2D5crOwT5qKi7I2eDBUqXLBC/1uvqw2e45nsDUuzUHBhBCikkpONvOTn3kGhg2DdeugYcPTTlm57wSjZ2ygfnV/vhrdmSr+XhaFLf+kWC4HtNZMXXmAJjUD6Vpcm5dp08DNDe64w/nhRMXn4wPDh5v2Qyml38q6f5twfDzd+G79EQeGE0KISmbxYmjVyiy+/vhjs79CQMBppyyKPsbIGeupX92fr8dcTlUplC+JFMvlwKqYZHYfy2BUt/pnb0dZUAAzZsD110N4uDUBRcU3apTZnOSrr0p9SZCPJ/1ahfPL1qNk5xc4MJwQQlQC2dnwyCNm2kVIiJl2cd99Z32i/P36I9z75UaahwXxjRTKZUKK5XJgyspYqgd4M7BtMcXwwoWQkHDWhH4hylTbttChA0yZYrpjlNItl9UmM6+ABduPOTCcEEJUcMuWQZs28N57MH48bNxofi6fQmvNJ3/G8N/Z2+gSVV2mXpQhKZZdXExiBn/uSeKOK+ri7eF+9glTpkDNmmZkWQhHGjUKtm0zP6RL6bJ6VahbzY9ZG2UqhhBCXLCUFLjnHrNwT2uz0PrDD//Zke9veQWFPP7DNt78dQ8D2oQz9U5ZzFeWpFh2cZ+vOoiXhxvDOxez005CAsyfbzaO8PQ8+7gQZWnYMDN/+QIW+imluLF9JGtjUziSku3AcEIIUYEUFppNRRo3Nj9zH33UDFb07HnWqcmZedw+ZR2zN8Xx0DWNeP/WttIerozJ76YLS8u2MWdTPAPbhFMtwPvsE2bMMH+hZHtr4QwhIWbTm6+/NnPnSmlIh0iUgjmbpOeyEEKc1+LF0LGjGVFu0QI2b4Z33gE/v7NO3XgohX4frGRrXBofDGvHI9c1Pnttk7hkUiy7sG/XHybHVsjdxW1trbV5t9m9u3nnKYQzjBplNr+ZNavUl0SE+HJlw2rM3hQnPZeFEKIk69ebxXu9ekFqqhmY+PNP0/niDH93ybpl4lq8PNyYc++VDGgji/wdRYplF1VQaOeLNYfoXL8qzcODzj5h+XKzvaUs7BPO1L07REVdcM/lGztEcjglm/UHTzoomBBClFNr1kDfvtCpE2zaBBMmwO7dZupbMaPEiem53D19PS/P20nPpjX45YGutIwItiB45SHFsotasus48ak5xY8qg1nYFxQEQ4c6N5io3JQyo8vLl8PevaW+rHeLWvh7uctCPyGEADOF8scfza67V14JGzbA66/DgQPw8MPgffbUS60187Ydpdd7y1kbm8wL/ZszaUQHgn1lzZKjSbHsoqatOkhkFV+ua17z7IMnT5qPwYcPL3YOkxAOdeed4O5+QaPLfl4e9GsdxvxtCdJzWQhReR07ZoriRo1gyBCIi4N33zVF8hNPQGBgsZcdSclm9IwNjP96M/Wq+bPgwW7c1aWYvReEQ0ix7IJ2Hk3nrwMp3HFFXdzdivmL8NVXkJsLY8Y4P5wQYWHQvz9Mnw75+aW+bGj7SLLyC1m887jjsgkhhKvJzYXZs2HwYKhdG556CurUMYNe+/aZjUbO2IHvn0tthXzyZwy9JixnTWwyz/RrxqxxV9AgtPjzhWNIEz4XNHPtQXw83bilYzHt4rSGyZOhfXto18754YQAM1f+p59g3jwzOlIKl9WrSkSIr+nw0jbCwQGFEMJC2dlmO+off4Sff4a0NLMnwkMPmYGuJk3Oebndrvl5azxvL9pLfGoOvZrX5PkBLYgI8T3ndcIxpFh2MWnZNn7cHM/gdhEE+xUzD2nDBtNr8dNPnR9OiL/16QMREeaNWymLZTc3xaB24Xz6536SMvIIDSymHaIQwlIFhXaOpeeSmJFHYnoeKVn5ZOTaSM+1kWuzYys0DzB91D3cFN4ebvh6uuPr5UGQrwfBvp5U8fOiWoAX1fy9qervVfynpBWJ1mYdx6JFpkheuhRycqBKFRg0yEyb7NkTPM5ddhXaNQt3JPDR0hh2H8ugZUQQb93UmisbVnfSL0QUR4plF/P9hiPk2uyMuLxe8SdMmWJ27hk2zKm5hDiNu7vp7/2//8Hhw+YjxVIY1DaCj//Yz9ytRxnVtYTFq0IIh7PbNbEnsog+msbOhHT2HMvgUHI2R1KyKbCf3eLRTYGvpzteHm54uLuhALvWFNg1eTY7uQWFlNQZ0sNNUSPQm1rBPkRU8SOyii+RVXypW9WfutX8CA/xLX/FdG4ubN0K69bBypWwYoWZjwymY9CoUaZI7t69VJuG5eQX8vOWeCaviGV/UhYNQv1575a2DGgTjlt5+72pgKRYdiGFds3MtYfoVK+EdnGZmabv4s03Q7C0iREW+7tYnjYNXnihVJc0qhlIy4ggftocL8WyxZRSfYD3AXdgitb69TOOewNfAB2AZOAWrfVBZ+cUZUNrzd7jmSzbm8i62BQ2Hj5JarYNAE93RVSNQJqHBXF9q1rUruJHzSAfagR5U83fm0AfD/y83M+5mExrTY6tkIzcAtJybJzMyudEZj4nMvNIzMglIS2XhNRctsWlsnB7wmkFuZe7G3Wr+VG/uj8NawQQFRpAVI0AGtYIIMDqLZsLCuDgQTNqvHMnbN9uPt2Njgab+f2jdm245hrT2eK666BBg1K/fExiJj9sOMK364+QlmOjaa1APhzWjutbhZW/NxAVmBTLLmTZ3kQOp2Tz3z4lzGX6/ntTMMvCPuEK6tUz/zBMmwbPPmtGm0thcLtIXp63k5jEDKJqFL/yWziWUsod+Bi4DogD1iul5mqtd55y2ijgpNY6Sil1K/AGcIvz04qLVWjXrD+YwvxtCSzZdZyEtFwAGoT606t5TTrWrUqryGAahgZc8vbISin8vDzw8/KgZpDPeXMdS8/lUHIWh5OzOXAii9gTWexPymTp7sTTCunwYB+iagbSqEYAjWsG0KhmIFE1AgjyucR2abm5ZuOPlBQ4ccI8EhLMIz7efGJ28CAcOfJvUQwQHm42CenTBy67zDxq176gWyek5bBoxzF+3HKUrUdScXdT9G5RkzuvqEen+lWlw4ULkmLZhcxYfYgagd70blGr+BMmTYJmzUxPRiFcwZgxZgvsRYvg+utLdUn/NmG8Mn8nP26O5z+9mzo4oChBJyBGax0LoJT6FhgInFosDwReKPp6FvCRUkpp2YbR5e1PyuS79Uf4cXM8SRl5+Hi60aNxKA9d04geTUIJC3bgIrHcXFN4pqSYRW2pqZCRAVlZZrAnNxf33Fwi8vKIyM/nSpvNjN4WFoLdjr2wkMycfDJybOaRayM9p4CMPBuFds1B4CBmSkiQtztBPh4EersT6OVOgIfCS2mUzWYKXJsN8vJMppycfzNkZJTcycfd3SzEq1sXOneGW24xu+Q2bmwW5VWrdsG/JfkFdrbFpbIqJpnfdx9nW1waAM3CgnimXzMGtA2nRuC532AIa0mx7CIOJWexbG8SD1/bCE/3Yt7hb9tm5kZNmFDsjj5CWGLAAKhRwyz0K2WxXCPQh66NQvlp81Eeu66JzMezRgRw6g4xcUDnks7RWhcopdKAasCJU09SSo0FxgLUKeXcdVH2Cu2axTuPM23VAf46kIKHm+LqpjUY0Dacq5vWwM+rDP65t9nMaOu+faYv8KFD5nH0qHkcP24K0tLw8jJzef9+uLuDuztuShHk5kaQUuD277+FGiiwa7PA0K6xFZiv8+0aO5Cq3ElRyryGlxce3l54+njh4e+HV3A1fCL98A4KQgX4m17GISFmOmPVqhAaaorgWrXM16X8lKw4drsm7mQOu46ls/VIKlvjUtl8OJXs/EKUgjaRIfy3TxN6Na8pn6yVI1Isu4iv1x3G3U0xrFMJ/9hMnmx29BkxwrnBhDgXLy+4+254+23zj2V4eKkuG9Q2nEe/38qmwyfpWK+qg0MKR9JaTwImAXTs2FFGnZ0sr6CQHzbEMWVFLAeTs4ms4st/+zThxg6RlzZaefw4rF8PW7aYwZrt2yEmxowC/83b20xBiIgwo7B/F5vVq5suEH8XpEFB4O9vHr6+5ueG24VN+1CAZ9HjVHa7Jj41h/1JmexPMlM5Dp7I4uCJLI4WTTv5m5eHG2HBPtQK8qFGkA81Ar2pHuBNFT9PQjy8CMr0wN+Wgb+3O17u7nh6KNzdFFqbxYy2Ak22rYCsvELSc2ykZOWTkpXP0bQc4k/mEHcyh9gTmeTaTLcQDzdF07BAbuwQyZUNq3F5g2qE+Hld+P8LYTkpll1Arq2Q7zccoVfzmsXP9crOhpkzzdbWF/ERkBAONXo0vPEGfP45PP10qS7p1aIW3h7bmbv1qBTL1ogHTp1oGVn0XHHnxCmlPIBgzEI/4QIKCu3M2RTP+7/vIz41hza1Q/ikT1N6t6h14QvDtDaL15YvN4/Vq82c3b81aGDm6Q4ZYqYjNGpknqtR44KL3rLm5qaoXdWP2lX9uOqM5T65tkLiTuZwJCWbwynZHE3L4WhqLsfSctgWl0pieh45tsJLzuDv5U5EFV8iQny5smE1GhXNrW4eFoSP58WPUgvXIcWyC1iwPYGT2TZuv7xu8Sf88IOZ+zV2rHODCVEaUVFw9dWmreFTT5XqH88Abw+ubVaT+dsSeO6G5ngUN/VIONJ6oJFSqj6mKL4VuO2Mc+YCdwJrgBuBpTJf2TX8sSeRl+ftJDYpizaRwbw+tBVdo6pf2MKwtDT49VdYuND0BU5IMM+Hh0PXrmbzjMsug7ZtS9yC2dX5eLoTVcN01ijO3x08Tmab7h2ZeQVk5RWQlV+IrcD0ky6wa9yUwk2Bh7sbfl7u+Hq5E+zrSTV/L6r4exHo7SGL8io4KZZdwMy1h2gQ6s+VDUsYNZ482byb797ducGEKK2xY+HWW2HJEujVq1SX9G8TzvztCazan0yPxqEODihOVTQHeTywCNM6bprWOlop9RKwQWs9F5gKzFRKxQApmIJaWOhQchYv/bKT33cnUr+6PxNHdKBX85qlL9RSUsy2y7Nnm00zbDYzZ/faa83f2549oX79SrMu5tQOHrIznjgXKZYttiM+jc2HU3n2hubF/8CLjoZVq8yc0EryA0yUQ4MGmXmKEyeWuli+qkkogT4ezN1yVIplC2itFwALznjuuVO+zgVucnYucbZCu2bqylje+W0vHm6Kp/o25e4u9UvX7s1mM9vSz5gBCxaY7xs2NCPHgwebucaXsKBNiMpAimWLfbXuMD6ebtzYPrL4EyZONIsh7rjDucGEuBDe3nDXXaZbS0IChIWd9xIfT3f6tKjFwh3HeMXWUub2CVGMmMQMHvthG1uPpHJts5r8b1BLagWXYuHe4cPm349p08zOcrVqwfjxZtvl9u1l8EWICyATBS2UkWvj5y3x3NA6nGC/YhqsZ2fDF1/AjTeaFcZCuLKxY02v1GnTSn3JgLbhZOYV8MfuRAcGE6L80VrzzV+HueHDlRxJyeaDYe2YfEeH8xfK69fDsGFmAd7rr0PHjjB3rtlc4913oUMHKZSFuEBSLFto7tajZOcXclvnEtrFffedWYQxbpxzgwlxMRo1MnMfJ00yRXMpXNGgGtUDvJm79aiDwwlRfqTn2rjvq008NWc7HetW5deHujGgTfi55yYvX27+/nXqZKZbPPIIxMbCL79A//7gIR8kC3GxpFi2iNaar9cdpmmtQNrVDin+pIkTzY59Xbs6N5wQF+uee8zHv4sWlep0D3c3+rWqxdLdiWTlFZz/AiEquJjEDAZ9tIrFO4/zZN+mfDGyEzXOtX30unVmYV6PHrBjB7z1FsTFmf/WLaHDkhDigkixbJFtcWlEH01neOc6xY8WbN5sfgiOGycfmYnyY+BAMzfys89KfUm/1uHkFdhZsuu4A4MJ4foWRR9j4EerSM+18fWYyxnXo2HJO1zu22em6F1+OezaZdYLxMbC44+X21ZvQrgqKZYt8vW6w/h6ujOwXUTxJ0ycCD4+smOfmOAc/gAAIABJREFUKF88PWHkSJg/38yRLIWOdatQM8ib+dsSHBxOCNektWbS8v3cM3MjUTUD+eWBrnSqX8JmPZmZ8OST0KKF+QTnxRfNznoPPwx+fs4NLkQlIcWyBdJzbczdepQBbcIJ8ilmYV9GBnz1lelbW6WK8wMKcSnGjDE7gk2eXKrT3dwU17cK48+9SWTk2hwcTgjXUmjXvDA3mlcX7KZf6zC+G3s5YcEl9PydPRuaNjU7Zg4fbork556DgOI33RBClA0pli3w85aj5NjOsbDvyy/N6IEs7BPlUb16cP31pli2la74vaF1OPkyFUNUMnkFhdz31UZmrDnE2O4N+PDWdsW3UExIgKFD/+2MtHq12V6+Zk3nhxaiEpJi2QLf/nWY5mFBtI4MPvug1vDJJ6a9T6dOzg8nRFm47z7T2/Wnn0p1ervaIYQH+zBvq0zFEJVDTn4ho2dsYFH0cZ7v35z/u75Z8fOTv/oKmjc3U5tef920hrviCucHFqISk2LZybYXLey7tVPt4hf2rVxpVjTfe68s7BPlV+/eZoT5k09Kdbqbm6Jf6zCW70siLUemYoiKLTOvgDs//4uVMSd4c2hr7u5S/+yTUlP/v707j4+quv8//jpJyMIWCAlhlU12kMWAuKIsghugtYBL61K1xapt/em3Wuyi1drWXbRWRRGrVioWATeQzY1FAdkhECM72di3rHN+f9yhRpghCZmZO8v7+Xjkkcmdm8ybO8nhM3fO/RxnqsX11zvF8qpV8NvfqgWciAtULIfY219vJSkhjpG9/VzY98ILkJrqNJUXiVTx8c40ogULnCv1q+GyM1pQVmGZvTYvuNlEXHS4pJwbX/2KZVv28vSY3ozu1/rEnZYsgd69nV77Dz0En34KnTqFPqyIACqWQ+pIaTnTV+zksp7NSU3xcWFffj5MnQo33aSrmiXy3Xyzs1T7Cy9Ua/derVJp2SiFj9aoWJboVFxWwa2vL2X51r08O7bPiSdNrIUJE+D88513Fr/8En7/e51NFnGZiuUQ+mDVLg6VlDO2v58L+155xbkgShf2STTIyIDRo2HyZOeC1SoYY7i0ZzM+31TIAXXFkChTWu5h3BvLWJS7mydG9+KyM5r/cIcjR5xpF3fd5UxjWr4czjrLnbAi8gMqlkPo7a+30T6jHv3a+mgHV17uLOQweDB07hz6cCLBMG4cHDjgXKRUDcN7NKeswjJvfUGQg4mEjsdjuXfqSuZnF/LIqJ5c2afVD3fYvt05m/z22/DIIzB9utqGioQRFcshklNwkGVb9jK2n58L+2bOdBZxuOOO0IcTCZazz4Y+feC555y3mKvQp3UjmjVM5sPV6ooh0eNvszYwfcVO/m945xNbhn71FfTr56zIN3Mm/O53EKf/mkXCif4iQ2TK19tIiDNc1beV7x0mTIDTToPLLw9tMJFgMsZ5AbhmjXORUhXi4gzDezTj042FHC4pD0FAkeCavHAzL36ay08GtGHcwA4/vPP99+HCCyElBRYtgssucyWjiJyciuUQKKvw8N/lOxjctSnp9ZNO3GHNGpg/3+lNqws5JNpccw00aeK8IKyGS3o0o6Tcw/xsTcWQyDZ3fT5/mrmWod0y+dOI7j98V3HiRBg50lm2evFi57OIhCUVyyEwb0MBuw+XMsZXiyCA55+H5GS45ZbQBhMJhZQU53f7vfdg69Yqd89qm0Z6/SQ+Wq2uGBK5NuUf5Fdvr6BHi1SeHduH+MoLjvzlL86y8EOHOidKmjZ1L6iIVEnFcgj85+ttNG2QxAUdM068c98+eP3178++iUSjceOcz//8Z5W7xscZhvfIZH52AUdLK4IcTCTw9h4u5ZbXl5JcJ56XfnomKYneJaythQcegPHjnc4XM2dC/fruhhWRKtWqWDbGpBljPjHGbPJ+9nn5rjGmwhizwvsxozaPGWnyDxQzP7uAH53ZioR4H4d70iSnZZAu7JNo1qYNjBgBL78MxcVV7n5Jj+YcKa3g042FIQgnEjjlFR7u+Pdydu0r5sWfnEnz1BTnDmvhnnucbhe33OKcJKnjo9++iISd2p5Zvg+Ya63tCMz1fu3LUWttb+/HiFo+ZkR5d/l2PBZGZ/mYglFR4XQJOOcc6Ns39OFEQunOO6GoCP797yp3PatdGo3q1tFqfhJxnpqzkS9zdvPwqB6c2cZ7/shauPtuePJJ5+/gxRfV8UIkgtT2r3UkMNl7ezIwqpY/L6pYa3ln6Xb6t02jXXq9E3d4/33IzYVf/zr04URC7aKLoGdPePrpKtvIJcTHMbhLJnPW51NW4QlRQJHamZ9dwPPzv2VMVuvvl7G21mkH9/TT8KtfwTPPqFAWiTC1/YvNtNYea4iaB2T62S/ZGLPUGLPYGHPSgtoYc5t336WFhZH9FuyyLXv5rugwV2f5aRf39NPQujVceWVog4m4wRjnheGqVbBgQZW7D+ueyYHicpbk7gl+NpFa2rHvKL+ZsoKuzRvy4MhKnS3+/Gf461+dlVmfesr5OxCRiFJlsWyMmWOMWePjY2Tl/ay1FvB3uqiNtTYLuBZ42hjTwc9+WGtfstZmWWuzMjJ8XBAXQaYu207dxHgu69n8xDtXrnQKhjvvVLs4iR3XXgvp6c4LxSpc0CmDlDrxfLxWC5RIeCuv8HDnW8spr7D847q+JNfxXtD37LPwxz/CjTc6XY9UKItEpCqLZWvtEGttDx8f04F8Y0xzAO9nn41RrbU7vJ9zgQVAn4D9C8LU0dIK3l+1i0t7Nqdeko9i+JlnoG5dtYuT2JKc7JxhmzkTcnJOvmudeAZ2ymD22nw8nqpX/xNxy4R5OSzfuo+/XNXz+yl3b7/tTLu48kqnp7KmXohErNr+9c4AbvDevgGYfvwOxpjGxpgk7+104FxgXS0fN+zNWpvHoZJyrj7TxxSMggJ480244QZo7LOBiEj0GjfOeTelGouUDOuRScHBElZs3xeCYCI1t2zLHibM28RVfVsyolcLZ+OcOfDTn8L558Nbb0F8vLshRaRWalss/xUYaozZBAzxfo0xJssYM9G7T1dgqTFmJTAf+Ku1NuqL5XeWbaN1Wgr926adeOcLL0BpKdx1V+iDibitRQsYMwZefdXpM34SgzpnkhBnmKWuGBKGDhaX8au3V9CycQoPjvDOU169Gq66Crp0gRkznHdTRCSi1apYttbuttYOttZ29E7X2OPdvtRae4v39kJrbU9rbS/v51cCETycbd97hIXf7uZHfVsRF3fcHLWjR512cZdf7gymIrHo7rvh0CF46aWT7pZatw5nd2jC7LX52Co6aIiE2oMz17FrfzFPj+lDg+Q6kJ/vjO0NGsCHH0KjRm5HFJEA0CSqIJi2fAfWwo/6+piC8frrTq/Ze+4JfTCRcNGnDwwa5FwAVVp60l2HdW/Gd0WHySk4FKJwIlWbv6GAqcu284uB7Z1+ykePwqhRzvg+Ywa08tMFSUQijorlALPW8u7y7Qxon0brtLo/vNPjcZrSZ2XBBRe4E1AkXNxzD+zYAVOmnHS3IV2djpSz1+WHIpVIlQ4Ul3H/f1fTKbM+dw3u6PRSvuUWWLwY3ngDzjzT7YgiEkAqlgNs+da9bN59xPdZ5fffh40bnSJBLYQk1g0fDt26weOPn3SRkmapyZzRKpU561UsS3h45P31FBws5rGre5GUEO/0T37rLWcpa/XNF4k6KpYD7N3lO0ipE88lvnorP/44tGkDP/pR6IOJhBtjnBeOq1Y53QNOYmjXTFZs20fBweIQhRPx7fNNhUxZuo2fD+xAr9aNYN48uPdeZ1y//36344lIEKhYDqDisgreX7mT4T2aUf/43spLlsDnnzt9N7UIiYjj2muhWTN47LGT7ja0eybWwtz1Plu5i4REcVkF46etoX16PX41uCNs3ep0dunSBSZN0juGIlFKxXIAzdtQwIHicq7q2/LEOx991OmprEVIRL6XlOS8gPzkE1i2zO9unTMb0KpxCnM0b1lc9Ny8HLbuOcLDV/Yg2VbAj3/sXKA6bZrTAUNEopKK5QB6d9l2MhsmcU6H9B/esW4dTJ/uLG2tAVXkh8aNg4YN4W9/87uLMYah3TL5IqeII6XlIQwn4tiUf5AXP/uWq/q0dMb4+++Hr75y+oV36uR2PBEJIhXLAVJ0qIQFGwsZ1acl8cf3Vv7b35ylre+8051wIuEsNRVuvx2mTnUugPVjaLdMSso9fLaxKIThRMDjsYyftoZ6SQmMv6yrs1z7k0/CHXfoGhSRGKBiOUBmrNhJhcee2AVjyxbnKulbb4X0dN/fLBLrfv1rSEw86dzlfm3TaJicoK4YEnLvrdjBV5v3cN/wLjTZkw833AB9+zoXbYtI1FOxHCDTvtlBj5YN6ZR53DSLJ55wPt99d+hDiUSKzEy4+WaYPNnpvexDnfg4BnVpyrwNBVR4tJqfhMahknL++tEGerVuxOg+LeCnP4WyMqc/eFKS2/FEJARULAdATsFBVu/Yz5V9jjurXFAAEyfC9dfDaae5E04kUtx7r7Nwz7EXmD4M7prJnsOlrNi2N4TBJJY9Pz+HgoMl/OmKbsQ99SR8+ik89xycfrrb0UQkRFQsB8C0b3YQZ+CKXsf1Vn7iCSgpgfvucyeYSCRp1w6uuw7++U/nhaYPF3TKICHOqIWchMTmosO88vl3XNW3JX12b4YHHoCrr3bOLotIzFCxXEsej+W9b3ZyXscMmjZI/v6OoiJ4/nkYOxY6d3YvoEgk+d3vnBeYfs4up6bUoV/bNOZtULEswffwB+upE2+4b2Ab54VcRobzYk79lEViiorlWlq6ZS879h3lqj7H9VZ+6ik4cgTGj3cnmEgk6tzZeYH5/PPOC04fBndtyoa8g2zfeyTE4SSWLPy2iDnr87ljUEeaPv4XWL8eXnsNmjRxO5qIhJiK5Vqa9s126ibGc3H3zO837tkDEyY4Deu7dXMvnEgkGj/eeaH51FM+7x7UpSmAzi5L0Hg8lkc/3EDLRincHL/LaRP385/D0KFuRxMRF6hYroXisgreX7WLYd2bUTex0hLWzzwDBw8689tEpGa6dXNeaE6Y4LzwPE77jPq0T6+necsSNDNX7WT1jv3ce0Frkm69BVq1gr//3e1YIuISFcu1sCC7gIPF5VxZeQrG7t3w9NNw1VXQs6d74UQi2QMPwKFDfvvYDurSlEXf7uZwiVbzk8AqKa/gsVnZdGvekJHvvQzZ2fDyy84qkyISk1Qs18J73+wkvX4S53SoNIftscecs8oPPuheMJFI17MnjBnjvEvjozPGoK5NKa3w8EWOVvOTwPrXoi1s33uUR9qUYp54An72M7j4YrdjiYiLVCyfov1Hy5i3oYArejUnId57GPPy4Nln4ZproEcPdwOKRLoHH4TiYnj00RPu6tc2jQbJCczVan4SQAeLy3hufg4XdGhMn4d/63S/OMmqkiISG1Qsn6JZa/IorfAwsnelKRiPPgqlpTqrLBIInTo5ywq/8AJs3/6Du+rEx3FBpwwWZBdirVbzk8B49YvN7DtSxqM7P4Nly5x3Nho3djuWiLhMxfIpem/FDto2qUuvVqnOhq1bnf6bN92klZ1EAuUPf3BW9fvzn0+466LOTSk4WMLanQdcCCbRZt+RUiZ+nsvYZpaWjz8Cw4fD6NFuxxKRMKBi+RTkHyhmUe5uRvRuiTnWnP6hh5zPv/+9e8FEok3btnDbbfDqq7Bp0w/uGtgpA3AutBWprZc+y+VQaTnjZ70IFRXwj39o8RERAVQsn5KZK3diLYzs3cLZsGYNTJoEt98Op53mbjiRaPPAA5CUBPff/4PNGQ2SOKNVKvOzC10KJtGi6FAJk77czP/Z72jw8fvOOxrt2rkdS0TChIrlUzB9xU56tkylQ0Z9Z8NvfwsNGqivskgwNGsG994L774Lixb94K6LOjflm6172Xu41KVwEg3+ueBbbHExP3vnaejYEX7zG7cjiUgYUbFcQ98WHmL1jv3fn1WeNw8+/NBZdUzLoIoEx//7f07RfM89UOmCvou6NMVj4bNNOrssp6boUAlvLNnC4zvmkZj7rdPRKCnJ7VgiEkZULNfQjBU7MQau6NXCufDonnucqRd33ul2NJHoVb++c13AwoUwbdr/Np/RMpUm9RJZoKkYcoomfv4djfcUcOn0V2DkSOfCPhGRSlQs14C1lpmrdnJWuzQyGybDm2/CN9/AI49AcrLb8USi2003QdeuzrSnUmfaRVycYWCnDBZkF1DhUQs5qZm9h0v516LNPLtiCnEVFfDkk25HEpEwpGK5BtbuPEBu4WFG9GrprNL3299Cv35w7bVuRxOJfgkJTjGTk+P0v/W6sEtT9h4pY+X2fS6Gk0g06cvv6Lh5Hf0WfuS8S9i+vduRRCQMqViugZkrd5IQZ7ikRzN4+GHYtQsmTIA4HUaRkBg+HK64wpmSsWsXABd0TCfOwIINaiEn1XeguIxJX37HE4snQ2amc/JDRMQHVXnV5PFYZq7cyfkd02m8YzM89RTceCOcdZbb0URiy5NPOtMw7rsPgEZ1E+lzWmM+3ah5y1J9/1q0hQtWLKDDplXOyY8GDdyOJCJhSsVyNS3buped+4sZ0buF01YoOdlZ3lpEQuv0053uGK+//r9WcgM7ZbBqx352HypxOZxEguKyCt5csJE/LXwdzjjDmQ8vIuKHiuVqmrlyJ0kJcQzPWeK0ivvjH51WViISer/7HbRo4SwEVF7OwE4ZWAtf5BS5nUwiwLRvdnDZZ++SUbQLnngC4uPdjiQiYUzFcjWUV3j4cPUuLm9Tj5Tf/Mo5E3HXXW7HEold9es7/XBXrIBnnqFny1TS6iXyqVrIVckYk2aM+cQYs8n7ubGf/SqMMSu8HzNCnTNYKjyWt2at4q6v3sEOGwZDhrgdSUTCnIrlalicu4eiQ6XcNf812LkTXn4Z6tRxO5ZIbLvqKudivz/8gbgtmzm/YzqfbSrEoxZyVbkPmGut7QjM9X7ty1FrbW/vx4jQxQuuT9blM+yjf9HgyEGMptKJSDWoWK6G91ft5OyCjZz29mtwxx3Qv7/bkUTEGHj+eacbze23M7BjOkWHSlm364DbycLdSGCy9/ZkYJSLWULKWsuUGUu4Zel0PGPGQp8+bkcSkQigYrkKpeUePlmxjSfm/APTsqWzAImIhIfWrZ2/yY8/ZsiKeQDqilG1TGvtLu/tPCDTz37JxpilxpjFxpioKKiXbdnLoKkvkeipIO7hP7sdR0QihIrlKnyZU8RNc1+nxbYc+Mc/1F5IJNz88pcwYAAN7/0NF9Qv1bxlwBgzxxizxsfHyMr7WWst4G/eShtrbRZwLfC0MaaDn8e6zVtULy0sDO9jP+O/nzN25Sw8t9zqdFUREakGFctVWDX1I8YtnkrFjTc58yNFJLzEx8PkyVBczJ+mP8WyLXs4UFzmdipXWWuHWGt7+PiYDuQbY5oDeD/7XM3FWrvD+zkXWAD4nLNgrX3JWptlrc3KyMgIyr8nELbvPcIZkyZAQgIJf/y923FEJIKoWD6J4n0HGPXU79ifnkn8M0+7HUdE/OnUCR57jPZLP2fMNx+xUC3kTmYGcIP39g3A9ON3MMY0NsYkeW+nA+cC60KWMAjen/opV66dT8mtP4fmzd2OIyIRRMXySRT+4i5a79nF1idfgIYN3Y4jIiczbhyeIUN4YN5E1n26zO004eyvwFBjzCZgiPdrjDFZxpiJ3n26AkuNMSuB+cBfrbURWywfLimn1YTHKE9MpP4fxrsdR0QiTILbAcLWO+/QespkXj/nR1xzjaZfiIS9uDjiJk3C06UHl//lN3DzYEhJcTtV2LHW7gYG+9i+FLjFe3sh0DPE0YJmzn8XcMXqBRSMu4tmTZu6HUdEIozOLPuyaRP2Zz9jRcsuZP/qfurE6zCJRIRWrVjypyfptDOHg7/4pdtpJAx4PJYGf/sLJUnJNHvoAbfjiEgEUhV4vOJiGD2asrgEbh/xf1zSp43biUSkBtrfOIYXzrqaBq9PgrfecjuOuGzZrIVcuHIBW6+7BdLT3Y4jIhFIxXJl1sK4cbBiBa/c9iDFzVsxoH2a26lEpAbaNqnL2yNuY2OnPnDbbbBqlduRxEXlD/+FkjqJtH1Ec5VF5NSoWK7sscfgtdcof+D3TEjuyLDuzUjQFAyRiGKM4ZwuzbjtsnuwjRo5LR/z8tyOJS7I+2Yt/RZ9zJoR15LUzN/aKyIiJ6dK8Jhp0+C++2DMGOaM/gVHSiu4/Ay1FxKJRAM7pbM5MZV1L74BRUUwahQcPep2LAmx/PEPUhEXT8s/66yyiJw6FcsAX38N118P/frBpEl8sCaftHqJnNVOUzBEItHZHdKJMzCrTnN44w1YsgRuuAEqKtyOJiFS8t1mus2expcDR9Kiq8/FB0VEqkXF8qpVMGwYNG0K06dTnJDI3PX5DO+hKRgikSo1pQ69Wzfis01FcOWVzhSrd96BX/wCPB6340kI7PjdQ2AtdR+43+0oIhLhYrsa3LABhg6FunVh3jxo1owF2QUcKa3gsp6agiESyc7vmMGq7fvYd6QU7rkHHngAJk6EX//auZhXoldRES3ffYs5fYfQf6DPVbpFRKotdovlDRtgyBDn9ty50K4dAB+sztMUDJEocF7HdDwWFufudjY89BDcfTdMmOAUzzrDHLV2//0pkspKOHDHr4mLM27HEZEIF5vF8qJFcO65UFYGn3wCnTsDUFxWwbz1+eqCIRIFerduRL3EeL7IKXI2GAOPPw533glPPgk33uiMARJdjhwh5cV/MPf0/gy66iK304hIFKhVRWiM+bExZq0xxmOMyTrJfsONMdnGmBxjzH21ecxamzkTBg+GtDRYuBDOOON/d322sZDDpRVc2rOZiwFFJBDqxMcxoH0TvthU9P1GY+CZZ+Dhh+Ff/3Layh086F5ICbiyia9Q98A+lo+5jYwGSW7HEZEoUNvTp2uAq4DP/O1gjIkHngcuAboB1xhjutXycWuuogL++EcYORK6d4cvv4QOP7xC+qM1eTSqW4cB7ZuEPJ6IBN65p6ezefcRtu058v1GY2D8eHjlFZgzB/r3h9Wr3QspgVNeTtnfH2dZiy70u/5yt9OISJSoVbFsrV1vrc2uYrf+QI61NtdaWwq8DYyszePW2K5dzoV8Dz0EP/0pLFjgdL+opKS8gjnr8rm4WyZ1NAVDJCqc19FZ3njht0Un3nnzzc40rH37nIJ54kRd+Bfp/vtf6u7YyjuDruH8Tk2r3l9EpBpCURW2BLZV+nq7d5tPxpjbjDFLjTFLCwsLa/fI5eXw3HPQrZvTZ/W115yPevVO2PXLnCIOlpRzibpgiESNjk3rk9EgiS9ydvve4aKLYMUKOO88uPVWp41kdlWv/yUsWUvJ3x4jt3ELmv9kDPG6sE9EAqTKYtkYM8cYs8bHR1DODltrX7LWZllrszIyMk71h8DHH0NWlnMxT1YWLF/uLErgx4er82iQnMC5HdJPMbmIhBtjDOedns7CnCI8Hj9njTMznfHi2Wfhq6+gZ0+4914tkR1pFi0iaflSXus3kh/3P83tNCISRaoslq21Q6y1PXx8TK/mY+wAWlf6upV3W+AdPQovvww9esAll8CePc5CBLNn/6/jhS+l5R5mr81jaLdMEhM0BUMkmpx7ejq7D5eyPu+A/53i450X1tnZcN118MQT0KYN3HSTs3CRhD371NMcSK5P3sjRtGiU4nYcEYkioagMvwY6GmPaGWMSgbHAjKA80ptvwm23QVKSc6V7Tg5cfbVzQc9JLMrdzYHici7toSkYItHmvNOdd4u+zPExb/l4mZkwaZJTNN96K/znP3DppVomO9xt2QL/fZe3eg3jinM7uZ1GRKJMbVvHXWmM2Q6cDXxgjJnl3d7CGPMhgLW2HLgDmAWsB/5jrV1bu9h+XHstzJ8Py5bB9ddDYmK1vu2j1buolxj/v4uBRCR6NEtN5vSm9f3PW/alY0fneodt22DqVOfMs4Sv557Dg2Hq2aMY2i3T7TQiEmUSavPN1tppwDQf23cCl1b6+kPgw9o8VrXUrQsXXlijb6nwWGavy2dQ10yS6+g/RJFodG6HJkxZuo2S8gqSEmrwd56WBgMGBC+Y1N6hQ9iXX+bjLudy1gW9NI6LSMDF/ATdr77bw57DpVzSQwuRiESrc05Pp7jMw4qt+9yOIoE2eTJm/34m9h3B1We2cjuNiEShmC+WZ63NIykhjoGdTrHzhoiEvQHtmhBnYOG3NZiKIeHPWnj+eXLadOFA7zPp3bqR24lEJArFdLHs8Vg+XpPHwE4Z1Euq1YwUEQljqXXr0KNlqu/FSSRyzZ8P69fzQvdLuPrM1pgqLuYWETkVMV0sr9i+j7wDxQzXFAyRqHd2hyZ8s3UfR0rL3Y4igfL88xxt2Ij3u53PlX38rnUlIlIrMV0sz1qTR514w+CuunpaJNqd2yGdco/l68173Y4igbBtG/a99/hv3+FkdW5Os9RktxOJSJSK2WLZWstHa/I4p0M6qSl13I4jIkGW1bYxdeINC6vTb1nC34svAvBCl6GM7K2zyiISPDFbLK/fdZCte45oCoZIjKibmECf0xrrIr9oUFICL71EdtYFFDZprnFcRIIqZovlj9fmEWdQA3uRGHJOhyas2bmf/UfK3I4itTFtGhQW8lznixnSNZOGyXp3UESCJ2aL5dlr88hqk0Z6/SS3o4hIiJx7ejrWOkvcSwR78UWOtm7DB817MEoX9olIkMVksbxl92E25B3k4u46qywSS3q1akRKnXi1kItkGzfCggXMOXcEqfWS1CNfRIIuJovlWWvzABjWXfPcRGJJYkIc/dqlsVhnliPXSy9hExL4e7MBXNazOYkJMfnfmIiEUEyOMh+vyaN7i4a0TqvrdhQRCbGz2zdhY/4hig6VuB1Faqq4GF57jZ0XDmNbUqq6YIhISMSD1dtUAAARw0lEQVRcsVxwoJjlW/fprLJIjBrQPg1AZ5cj0bRpsHs3U/peQrOGyWS1aex2IhGJATFXLM9elw9oCoZIrOrZMpV6ifEqliPRiy/iadeeF+PbcmnP5sTFaXlrEQm+mCuWZ63No22TunTKrO92FBFxQUK8M295kfotR5ZNm+DTT1l72WhKPHB5r+ZuJxKRGBFTxfL+o2Us+nY3w7o3wxidkRCJVWe3b8K3hYcpOFjsdhSprkmTID6eie3Pp2WjFPq0buR2IhGJETFVLC/ILqDcY7lYUzBEYtrZHZoAsDh3j8tJpFrKy+G11yi9eBgfFMLlZzTXCQ8RCZmYKpZnr80nvX6SzkiIxLhuzRvSIClBUzEixaxZsGsXSy4cSbnHckWvFm4nEpEYkuB2gFApLqtgQXYBI3q31EUhIjEuIT6O/u3SWKKL/CLDq69C06a83LAbbU0Z3Vs0dDuRiMSQmDmzvOjb3RwurdCqfSICOFMxcosOk39A85bDWkEBzJjB0THX8sWW/VymKRgiEmIxUyzPXpdHvcR4zvHOVRSR2DagvTMWaCpGmHvjDSgvZ/45l+GxcEkPdcEQkdCKiWLZ47F8sq6AC7s0JSkh3u04IhIGujZvSMPkBJZ8p2I5bFnrdMHo358pRxrSOi1FUzBEJORiolj+Zts+ig6VcHE3TcEQEUd8nKFf2zSWqCNG+FqxAtas4ei1P+HLnCIu7aEpGCISejFRLM9el0edeMNFXZq6HUVEwshZ7dPILTpMgeYth6fXX4fERD7pOZByj2V4D7X9FJHQi4li+ZO1+Qxo34SGyXXcjiIiYeSsds685SXf6exy2Ckrg7fegiuuYMa2EpqnJtOrldp+ikjoRX2xnFNwiNyiwwzVFAwROU73Fg2pn6R5y2Fp9mwoKODoNdfx2aZChvdoprafIuKKqC+W56zPB2BIVxXLIvJDCfFxnNmmseYth6PJkyE9nblt+lBa7lEXDBFxTdQXy5+sy6dHy4a0aJTidhQRCUNntU9jU8Ehdh8qcTuKHLN3L8yYAddcw0fZe0ivn8SZbRq7nUpEYlRUF8tFh0pYvnWvziqLiF/H5i1/pXnL4eOdd6CkhJLrrmdBdgFDu2USrykYIuKSqC6W560vwFo0X1lE/DqjVSopdeJ1kV84efNN6NKFhQ3baOVVEXFdVBfLs9fl07JRCt2aq4m9iPhWxztveXGuLvILC9u2wWefwXXXMXt9vlZeFRHXRW2xfLS0gi9yChnStama2IvISZ3VLo3s/IPsO1LqdhT5978B8IwZq5VXRSQsRG2x/EVOEcVlHoZ2UxN7ETm5s9o3wVrNWw4Lb70FAwbwTWITrbwqImEhaovlOevyaZCUQP92aW5HEZEwd0arVBIT4vh6s4plV61dCytXwrXX8sm6fBLiDBd21sqrIuKuqCyWPR7L3A0FDOycQWJCVP4TRSSAkuvE07tVI51Zdttbb0F8PIwezex1eZzdoQmpKVp5VUTcFZWV5Mrt+yg6VKKWcSJSbf3bpbFm5wEOl5S7HSU2WesUy0OGkGPqkVuolVdFJDxEZbE8d30B8XGGCztnuB1FRCJEv3ZpVHgsy7fudTtKbFq8GDZv/t8UDNDKqyISHqKyWJ6zPp+sNo1pVDfR7SgiEiHObNOYOANfayqGO95+G5KSYNQo5m3QyqsiEj6irljetucIG/IO6u07EamR+kkJdG+RqsVJ3FBR4azad+ml7I1PZtmWvQzqojFcRMJD1BXLc9c7b98N1tt3IlJD/dulsWLbPkrKK9yOElu+/BJ27YLRo1mwsQCPhcFd1AVDRMJD9BXLGwrokFGPdun13I4iIhGmX9s0Sso9rN6+3+0osWXKFEhJgcsvZ+76AjIaJNGzZarbqUREgCgrlg8Wl7E4d7cuChGRU9KvbWMAvlK/5dCpqICpU+GyyyhLqcunGwsZ1LkpcXFaeVVEwkNUFcufbSyirMIyRPOVReQUNKmfxOlN60d1v2VjzI+NMWuNMR5jTNZJ9htujMk2xuQYY+4LWqBPP4WCAhgzhqWb93KwuJxBXTUFQ0TCR1QVy71Pa8QDl3Wl72mN3Y4iIhGqX9s0lm3eS4XHuh0lWNYAVwGf+dvBGBMPPA9cAnQDrjHGdAtKmilToF49uPRS5m3IJzE+jvNOTw/KQ4mInIqoKpZbNkrhlvPbE6+370TkFPVv15iDJeVk5x10O0pQWGvXW2uzq9itP5Bjrc211pYCbwMjAx6mvBzefReuuALq1mXuhgIGdGhCvaSEgD+UiMipiqpiWUSktvq1TQNg6ZbonYpRDS2BbZW+3u7ddgJjzG3GmKXGmKWFhYU1e5TNm6FBAxgzhu+KDpNbeFhdMEQk7KhYFhGppGWjFJo1TObrzZG7kp8xZo4xZo2Pj4CfHbbWvmStzbLWZmVk1HDV1NNPh9xcuOIK5m8oAGCQimURCTN6r0tEpBJjDFltG/P1d3uw1mJM5E3rstYOqeWP2AG0rvR1K++2wDMG4uNZsLGQDhn1aJ1WNygPIyJyqnRmWUTkOP3appF3oJgd+466HcUtXwMdjTHtjDGJwFhgRrAe7GhpBYtzd3NhZ51VFpHwo2JZROQ4Wd5+y0sjeCqGP8aYK40x24GzgQ+MMbO821sYYz4EsNaWA3cAs4D1wH+stWuDlWlx7m5Kyz1c2LmG0zhEREJA0zBERI7TpVlD6icl8PXmPYzq4/O6tohlrZ0GTPOxfSdwaaWvPwQ+DEWmBdkFpNSJp3+7tFA8nIhIjdTqzHINmttvNsasNsasMMYsrc1jiogEW3ycoW+bxlF5ZjkcLdhYyDkdmpCUEO92FBGRE9R2GkaVze0rucha29ta67eoFhEJF/3aNCY7/yD7j5S5HSWqfVd0mC27j2gKhoiErVoVy9Vsbi8iEnGyvP2Wl22N6X7LQbcg22kZp4v7RCRcheoCPwvMNsYsM8bcdrIda9XgXkQkQHq3bkRCnInofsuRYH52Ie3VMk5EwliVF/gZY+YAzXzcNd5aO72aj3OetXaHMaYp8IkxZoO11ufUDWvtS8BLAFlZWbaaP19EJKBSEuPp0TKVpZt1ZjlYjrWMu/6sNm5HERHxq8piOQDN7bHW7vB+LjDGTAP6U715ziIirrm4eya5hYcjdnGScLfnSClnt2/CkK6agiEi4SvoreOMMfWAOGvtQe/ti4GHgv24IiK1dfuFp7sdIaq1bJTC5Jv7ux1DROSkats6rsrm9kAm8IUxZiXwFfCBtfbj2jyuiIiIiEgo1OrMcnWa21trc4FetXkcERERERE3aLlrERERERE/VCyLiIiIiPihYllERERExA8VyyIiIiIifqhYFhERERHxQ8WyiIiIiIgfKpZFRERERPxQsSwiIiIi4oeKZRERERERP1Qsi4iIiIj4oWJZRERERMQPY611O4NfxphCYEsNvy0dKApCnFOhLCcKlxwQPlnCJQcoiy+nmqONtTYj0GHC2SmO2RD5z3UwhEuWcMkByuJLuOSA8MkS8DE7rIvlU2GMWWqtzXI7ByhLOOeA8MkSLjlAWcI5RzQLl2McLjkgfLKESw5QlnDOAeGTJRg5NA1DRERERMQPFcsiIiIiIn5EY7H8ktsBKlGWE4VLDgifLOGSA5TFl3DJEc3C5RiHSw4InyzhkgOUxZdwyQHhkyXgOaJuzrKIiIiISKBE45llEREREZGAiMhi2RjzY2PMWmOMxxjj94pHY8xwY0y2MSbHGHNfpe3tjDFLvNunGGMSa5ElzRjziTFmk/dzYx/7XGSMWVHpo9gYM8p732vGmO8q3dc7WDm8+1VUeqwZlbaH+pj0NsYs8j6Pq4wxYyrdV6tj4u95r3R/kvffmOP9N7etdN/93u3ZxphhNfuXn1KWu40x67zHYK4xpk2l+3w+V0HKcaMxprDS491S6b4bvM/lJmPMDbXJUc0sT1XKsdEYs6/SfYE8Jq8aYwqMMWv83G+MMc96c64yxvStdF9Aj0ksCJdxO1zG7Opm8e4X1HHb7THb+zPCYtwOlzG7mllCMm5rzAastRH3AXQFOgMLgCw/+8QD3wLtgURgJdDNe99/gLHe2/8ExtUiy9+B+7y37wP+VsX+acAeoK7369eAqwNwTKqVAzjkZ3tIjwnQCejovd0C2AU0qu0xOdnzXmmf24F/em+PBaZ4b3fz7p8EtPP+nPhaHIfqZLmo0u/CuGNZTvZcBSnHjcBzfn5fc72fG3tvNw5mluP2vxN4NdDHxPuzLgD6Amv83H8p8BFggAHAkmAck1j5IEzG7eqMT8ftH5QxuyZZ/P3eh/KYEKQxu6rnvdI+QR+3q5kj6GN2DbLcSJDH7erkOG7/qByzI/LMsrV2vbU2u4rd+gM51tpca20p8DYw0hhjgEHAVO9+k4FRtYgz0vszqvuzrgY+stYeqcVjBiLH/7hxTKy1G621m7y3dwIFQCAWcPD5vJ8k31RgsPcYjATettaWWGu/A3K8Py9oWay18yv9LiwGWtXi8U45x0kMAz6x1u6x1u4FPgGGhzDLNcC/a/F4fllrP8MpgvwZCbxuHYuBRsaY5gT+mMSEMBq3w2XMPpUs/xPqYxLEMRvCZ9wOlzG7WllOIpBjlMZsInQaRjW1BLZV+nq7d1sTYJ+1tvy47acq01q7y3s7D8isYv+xnPiL9Ij3LYOnjDFJQc6RbIxZaoxZfOxtRVw+JsaY/jivWL+ttPlUj4m/593nPt5/836cY1Cd762Jmv68n+G8Kj7G13MVzBw/8h7zqcaY1jX83kBnwfv2ZjtgXqXNgTom1eEva6CPiXwvFON2uIzZNckS7HHbzTEbwmfcDpcxuyZZgj1ua8wGEmodLUiMMXOAZj7uGm+tnR4uWSp/Ya21xhi/7UW8r3B6ArMqbb4fZ3BKxGl38lvgoSDmaGOt3WGMaQ/MM8asxhl0aiTAx+RfwA3WWo93c7WPSbQwxlwPZAEDK20+4bmy1n7r+yfU2kzg39baEmPMz3HO4AwK0mNV11hgqrW2otK2UB4TqaFwGbfDZcwOYJZaj9saswMrDMZsCL9xO2rH7LAtlq21Q2r5I3YArSt93cq7bTfOqfkE76vTY9tPKYsxJt8Y09xau8s7iBSc5EeNBqZZa8sq/exjr+ZLjDGTgHuCmcNau8P7OdcYswDoA7yLC8fEGNMQ+ADnP9LFlX52tY+JD/6ed1/7bDfGJACpOL8X1fnemqjWzzPGDMH5D2ugtbbk2HY/z9WpDDJV5rDW7q705UScOYzHvvfC4753wSlkqHaWSsYCvzwuZ6COSXX4yxroYxI1wmXcDpcxO1BZAjFuh/GYDeEzbofLmF2tLCEatzVmE93TML4GOhrnauFEnCdxhrXWAvNx5qEB3ADU5ozHDO/PqM7POmEuj3dgOjb/bBTg8yrPQOQwxjQ+9vaYMSYdOBdY58Yx8T4n03DmF0097r7aHBOfz/tJ8l0NzPMegxnAWONcdd0O6Ah8VYPHrnEWY0wf4EVghLW2oNJ2n89VEHM0r/TlCGC99/Ys4GJvnsbAxfzwLFvAs3jzdMG5EGNRpW2BPCbVMQP4qXEMAPZ7i4JAHxP5XijG7XAZs6uVJUTjtptjNoTPuB0uY3Z1s4Ri3NaYDRHbDeNKnDknJUA+MMu7vQXwYaX9LgU24ryKGV9pe3ucP6Yc4B0gqRZZmgBzgU3AHCDNuz0LmFhpv7Y4r27ijvv+ecBqnMHlDaB+sHIA53gfa6X388/cOibA9UAZsKLSR+9AHBNfzzvOW4IjvLeTvf/GHO+/uX2l7x3v/b5s4JIA/K5WlWWO93f42DGYUdVzFaQcjwJrvY83H+hS6Xtv9h6rHOCmYB8T79d/Av563PcF+pj8G+eK/jKc8eRnwC+AX3jvN8Dz3pyrqdTBIdDHJBY+CJNxmzAZs6ub5WS/96E8JgRxzPb3vOPCuF2NHCEZs6uZJSTjdlU5vF//iSges7WCn4iIiIiIH9E8DUNEREREpFZULIuIiIiI+KFiWURERETEDxXLIiIiIiJ+qFgWEREREfFDxbKIiIiIiB8qlkVERERE/FCxLCIiIiLix/8HyPtYcPFbNUMAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 864x432 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"xx = np.linspace(-1, 1, 100)\n",
"fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 6))\n",
"for j, (smoother, ax, f) in enumerate(zip(smoothers, axes.flat, [f1, f2])):\n",
" f_j = smoother.predict(xx)\n",
" ax.plot(xx, f_j - np.mean(f_j))\n",
" ax.plot(xx, f(xx) - np.mean(f(xx)), 'r')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.5"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment