{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Gaussian process regression\n", "\n", "[Bayesian linear regression](bayesian_regression.ipynb) derived linear regression in a Bayesian context. Here, we discuss Gaussian process regression using GPy and scipy. Most of the material is from [Rasmussen and Williams (2006)](http://www.gaussianprocess.org/gpml/). I also recommend Michael Betancourt's [Robust Gaussian Processes in Stan](https://betanalpha.github.io/assets/case_studies/gp_part1/part1.html) as a resource, for instance to learn more about hyperparameter inference which won't be covered here.\n", "\n", "As usual **I do not take warranty for the correctness or completeness of this document.**" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import GPy\n", "import scipy\n", "from sklearn.metrics.pairwise import rbf_kernel\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.rcParams['figure.figsize'] = [15, 6]" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "rnorm = scipy.stats.norm.rvs\n", "mvrnorm = scipy.stats.multivariate_normal.rvs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Priors on functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Bayesian linear regression we assumed a linear dependency\n", "\n", "\\begin{align*}\n", "f_{\\boldsymbol \\beta}& : \\ \\mathcal{X} \\rightarrow \\mathcal{Y},\\\\\n", "f_{\\boldsymbol \\beta}(\\mathbf{x}) & = \\ \\mathbf{x}^T \\boldsymbol \\beta + \\epsilon,\n", "\\end{align*}\n", "\n", "which was parametrized by a coefficient vector $\\boldsymbol \\beta$. \n", "In order to model uncertainty, regularize our coeffiencts, or what reason whatsoever, we put a prior distribution on $\\boldsymbol \\beta$ and by that introduced some prior belief to the model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When we use Gaussian Processes, we instead set a prior on the function $f$ itself:\n", "\n", "\\begin{align*}\n", "f(\\mathbf{x}) & \\sim \\mathcal{GP}(m(\\mathbf{x}), k(\\mathbf{x}, \\mathbf{x}')) ,\\\\\n", "p(f \\mid \\mathbf{x}) & = \\mathcal{N}(m(\\mathbf{x}), k(\\mathbf{x}, \\mathbf{x}')) .\n", "\\end{align*}\n", "\n", "So a Gaussian process is a distribution of functions. It is parametrized by a *mean function* $m$ that returns a vector of length $n$ and a *kernel function* $k$ that returns a matrix of dimension $n \\times n$, where $n$ is the number of samples. For instance, the mean function could be a constant (which we will assume throughout the rest of this notebook), and the kernel could be a radial basis function, i.e.:\n", "\n", "\\begin{align*}\n", "m(\\mathbf{x}) &= \\mathbf{0},\\\\\n", "k(\\mathbf{x}, \\mathbf{x}') &= \\exp\\left(- \\gamma ||\\mathbf{x} - \\mathbf{x}' ||^2 \\right),\n", "\\end{align*}\n", "\n", "where $\\gamma$ a hyperparameter we have to optimize. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The parameters $\\mathbf{m}$ and $\\mathbf{k}$ apparently do not have a fixed dimensionality as in Bayesian linear regression (where we had $\\boldsymbol \\beta \\in \\mathbb{R}^p$), but have possibly infinite dimension. That means with more data, the dimensions of $\\mathbf{m}$ and $\\mathbf{k}$ increase (in Bayesian regression $\\boldsymbol \\beta$ was independent of the sample size $n$). For that reason we call this approach *non-parametric* (the turn itself sounds confusing, because we apparently *have* parameters)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's look at some *prior functions* $f$. A prior function is just a sample from a $n$-dimensional multivariate normal distribution with mean $\\mathbf{m}$ and variance $\\mathbf{k}$. The kernel must be postive definite, which the *RBF* kernel is." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "n = 50\n", "x = scipy.linspace(0, 1, n).reshape((n, 1))\n", "beta = 2\n", "y = scipy.sin(x) * beta + rnorm(size=(n, 1), scale=.1)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAGnJJREFUeJzt3X+QXWd93/H3R2sbRuASS9qk1NLuuqlDY1zwjzsK6Q+wh2CEZ2qlE9JKiMS0JhpITNOkw4wzOwMdM56hpQ2QBAdE0dgBxSalJVFTg3EC1AlBxFeNcWx3DIqw7K2Y8cZy3CFybSR/+8c5N77evefuubvnOeeeez+vmTt7zznPvfscSz5fPb++jyICMzOztWxqugJmZtYODhhmZlaKA4aZmZXigGFmZqU4YJiZWSkOGGZmVooDhpmZleKAYWZmpThgmJlZKec0XYEqbdu2LRYWFpquhplZaxw9evQvI2K2TNmJChgLCwt0u92mq2Fm1hqSTpQt6y4pMzMrxQHDzMxKccAwM7NSHDDMzKwUBwwzMyvFAcPMzEpxwDAzs1KSBQxJByU9IenBguvvlXR//npQ0llJW/Jrj0r68/yaF1aY2dQ5dAgWFmDTpuznoUNN1yhtC+M2YFfRxYj4UERcFhGXAb8C/M+IONVX5Or8eidhHc3Mxs6hQ7B/P5w4ARHZz/37mw8ayQJGRNwLnFqzYGYvcEequpiZtcniIpw+/eJzp09n55vU+BiGpM1kLZH/2nc6gC9JOipp/xqf3y+pK6m7vLycsqpmZrV47LHi8012VTUeMIB/CnxtRXfUP4qIK4C3AL8g6fVFH46IAxHRiYjO7Gyp/FlmZmNtbm7w+S1bmu2qGoeAsYcV3VERcTL/+QTweWBnA/UyM2vELbfA5s0vPtc7brKrqtGAIekVwBuA3+s79zJJ5/feA9cAA2damZlNon374MABmJ8HKft54ACcKhgVLurCqlqy9OaS7gCuArZJWgLeD5wLEBEfz4v9M+BLEfHXfR/9IeDzknr1++2I+GKqepqZjaN9+7JXv8XFrBtqpaIurKolCxgRsbdEmdvIpt/2nzsOvDZNrczM2uuWW7Ixi/5uqc2bs/N1GIcxDDMzK6Goq2plSyQVBwwzs3Voanrrvn3w6KPw/PPZz7qCBUzYFq1mZnXorcTudQ31prdCvQ/wurmFYWY2onFdiZ2aA4aZ2YiGrcQuMo7JBEflgGFmNqKiaaxF58c1meCoHDDMzEZUtBK7aHrrsC6sNrU8HDDMzEY06vTWoq6qXkujLS0PRUTTdahMp9OJbtf7LZnZeFlYGLxCe2YGzp5dfX5+PpsyWwdJR8vuO+QWhplZYkVdWIOCBdSXG2pUDhhmZokVdWHNzw8uX1duqFF54Z6ZWQ0GJROEZnNDjcotDDOzhjSdG2pUDhhmNpHaMl21ydxQo3KXlJlNnGnN9ZSaWxhmNnGqzPXUlpZKHdzCMLOJs55cT4O4pfJibmGY2cQZNddTkWnNSlskWcCQdFDSE5IeLLh+laSnJd2fv97Xd22XpEckHZN0U6o6mtlkGjXXU5FpzUpbJGUL4zZg1xpl/igiLstfNwNImgE+BrwFuATYK+mShPU0swlT1XTVac1KWyRZwIiIe4FT6/joTuBYRByPiOeAO4HdlVbOzCZeFdNVq8xKOwmaHsP4cUnflPQFSa/Oz10IPN5XZik/Z2ZWq6qy0o5rbqhRNRkw/hcwHxGvBX4d+N38vAaULUypK2m/pK6k7vLycoJqmlnTmhwXGKWlUtVg+7hqLGBExP+NiO/l7+8CzpW0jaxFsaOv6Hbg5JDvORARnYjozM7OJq2zmdWvTeMCVQ22j6vGAoakvy1J+fudeV2eBO4DLpZ0kaTzgD3A4abqaWbNatO4QNtyQ40q2cI9SXcAVwHbJC0B7wfOBYiIjwNvBd4t6QzwDLAnst2czki6EbgbmAEORsRDqeppZuOtbeMCRVlpJ0GygBERe9e4/hvAbxRcuwu4K0W9zKxd5uYG71Y3KeMCbdL0LCkzs6EmfVygTRwwzGysTfq4QJs4+aCZjb1JHhdoE7cwzMysFAcMM2utSU70N44cMMyslape0OfgszYHDDNLJuVDuOpd9dqymrxJytbKTYZOpxPdbrfpapgZq3erg2w6bFUznDZtyh7uK0lZ3qdRLCwMXusxP5/lj5pkko5GRKdMWbcwzCyJ9bQARmmRVJnor22ryZvigGFmSQx7CA8KDKN2C1W5oG/Ss8xWxQHDzJIoethu2TI4MPziL47WIqlyQZ9Xk5fjgGFmSRQ9hGFwYHjyycHfM6xbqIpd9Xrf49Xka3PAMLMkih7Cp0bcuLmubqGqgs8kc8Aws2QGPYSLAsDWre4WGncOGGZWWhXrKoq6qj760Xq6hbxAb/2cfNDMSlm5rqI3WA2jPdR7ZRcXs/GJubksiPTOp+wKquoeppUX7pnZKocOrX6gLy62f3HbNC/QKzLKwj23MMzsRYr+Fb5yZlNPmxa3eYHexngMw8xepGiF9szM4PJtWtzmBXobkyxgSDoo6QlJDxZc3yfpgfz1J5Je23ftUUl/Lul+Se5jMqtR0b+2z55t/ywmL9DbmJQtjNuAXUOufwd4Q0S8BvgAcGDF9asj4rKyfWtmVo2if233Zi21eXGbF+htTLKAERH3AoVLdCLiTyLiqfzwCLA9VV3MrLxh/wovWtzWpqmqXqC3fuMyhnED8IW+4wC+JOmopP0N1clsKo36r3DvJTE9kk6rlbQA/H5EXDqkzNXArcA/jogn83N/JyJOSvpB4B7gPXmLZdDn9wP7Aebm5q48MWjOnJkl46mq7daa/TAkvQb4z8DuXrAAiIiT+c8ngM8DO4u+IyIOREQnIjqzs7Opq2xmK3iq6vRoLGBImgP+G/AzEfGtvvMvk3R+7z1wDTBwppWZNc9TVadHymm1dwBfB14laUnSDZLeJeldeZH3AVuBW1dMn/0h4I8lfRP4U+B/RMQXU9XTzDbGU1WnR7KV3hGxd43r7wTeOeD8ceC1qz9hZuNordxQNjmcGsTMNmzfPgeIaTAu02rNzGzMOWCYmVkpDhhmZlaKA4aZmZXigGE2ptqUn8mmgwOG2RiqKz+Tg5KNwgHDbAwVbWK0uFjd73DSQBuVA4bZGBqWn6mqVkEdQckmiwOG2RgqysO0ZUt1rQInDbRROWCYjaGi/ExQXavASQNtVA4YZmOoaBOjUwV7WK6nVeCkgTYqBwyzMTVoK9EqWwXe39pG5YBhlmvDFNNhrYL11N/7W9sonK3WjBemmPbGB3qDyTBeD9GiVOLQjvpbuyXd07tunU4nut3u2gXNVmj7vtRtr781pzV7epuNi7ZPMW17/a0dHDDMqHYwuYmxEE+RtTo4YJhR3RTTptJteIqs1SFpwJB0UNITkh4suC5JvybpmKQHJF3Rd+16Sd/OX9enrKdZVVNMm0q34SmyVoekg96SXg98D/itiLh0wPVrgfcA1wI/Bnw0In5M0hagC3SAAI4CV0bEU8N+nwe9rWmbNmUti5Uk+PSnV89u8gPdmjY2g94RcS9QsDYVgN1kwSQi4gjwA5JeCbwZuCciTuVB4h5gV8q6mlWhjhxQZk1pegzjQuDxvuOl/FzR+VUk7ZfUldRdXl5OVlGzMurIAWXWlKYDhgaciyHnV5+MOBARnYjozM7OVlo5s1HVkQPKrClNB4wlYEff8Xbg5JDzZhuWetpr6hxQZk1pOmAcBn42ny31OuDpiPgucDdwjaQLJF0AXJOfM9sQT3s1W7/U02rvAL4OvErSkqQbJL1L0rvyIncBx4FjwCeBnweIiFPAB4D78tfN+TmzDfG0V7P1cy4pmyrDpr0+/3z99TFr2thMqzUbN+M4ltCGtOpm4IBhU2Y9YwkpH+hNjamYrYcDhk2VUccShj3QqwgkTY2pmK2HxzDMhijaZ2LrVnjmmRc/7DdvHn0g22Mq1jSPYZhVpGhh3ZNPVtMyGMcxFbMihQFD0l2SFuqritn4GfXBPerK7brWZ3hg3aowrIVxG/AlSYuSzq2pPmZjpeiBvnXr4PKjBpg61md4YN2qMnQMQ9LLgPeRZYr9NPA3vaoR8avJazcij2FYCocOrU5LDtlDd6NjGHXwft82zChjGOescf37wF8DLwHOpy9gmE2LffuKg0Ab9rfwft9WlcKAIWkX8Ktk+Z6uiIjTRWXNptGwQDJO5uYGtzA8sG6jGjaGsQj8dETc5GBh1l5OfGhVKQwYEfFPIuKhOitjZtVz4kOrylpjGGY2AdrSfWbjzQv3zMysFAcMMzMrxQHDzMxKccAwM7NSHDCsEc5tZNY+qff03iXpEUnHJN004PqHJd2fv74l6a/6rp3tu3Y4ZT2tXs5tZNZOyfbDkDQDfAt4E7AE3AfsjYiHC8q/B7g8Iv5Vfvy9iHj5KL/TuaTawbmNzMbHuOyHsRM4FhHHI+I54E5g95Dye4E7EtbHxoRzG5m1U8qAcSHweN/xUn5uFUnzwEXAl/tOv1RSV9IRST+ZrppWN28aZNZOKQOGBpwr6v/aA3wuIs72nZvLm0lvAz4i6YcH/hJpfx5YusvLyxursdWiytxGHjw3q0/KgLEE7Og73g6cLCi7hxXdURFxMv95HPgqcPmgD0bEgYjoRERndnZ2o3W2GlSV28iD52b1SjnofQ7ZoPcbgf9DNuj9tpUJDSW9CrgbuCjyyki6ADgdEc9K2gZ8HdhdNGDe40Hv6eLBc7ONq3IDpXWLiDOSbiQLBjPAwYh4SNLNQDcielNl9wJ3xosj148Cn5D0PFkr6INrBQubPh48N6tXshZGE9zCmC5uYZht3LhMqzVLyhsDmdXLAcPGStGsp0HnvTGQWb28gZKNjd6sp9P5hsC9WU9f+xrcfvvq8+CNgczq5DEMGxtFYxIzM3D27OrzHqsw2ziPYVgrFc1uGhQshpU3szQcMGxsFKUGmZkZrbyZpeGAYWOjaNbT/v2eDWU2DhwwbGwUzXq69VbPhjIbBx70NjObYh70tg1zFlgzW8nrMGyVovUQ4G4gs2nmFoatsrj4QrDoOX06O29m08sBw1ZxFlgzG8QBw1bxFqpmNogDhq0yCVlgPWhvVj0HDFul7VlgvXWrWRpeh2ETxxsrmZXndRg21Txob5aGA4ZNHA/am6WRNGBI2iXpEUnHJN004Po7JC1Luj9/vbPv2vWSvp2/rk9ZT5sskzBobzaOkq30ljQDfAx4E7AE3CfpcEQ8vKLoZyPixhWf3QK8H+gAARzNP/tUqvra5OgNzi8uZt1Qc3NZsGjLoL3ZuEqZGmQncCwijgNIuhPYDawMGIO8GbgnIk7ln70H2AXckaiuNmG8datZ9VJ2SV0IPN53vJSfW+mnJD0g6XOSdoz4WTMzq0nKgKEB51bO4f3vwEJEvAb4A+D2ET6bFZT2S+pK6i4vL6+7spaGF9CZTY6UAWMJ2NF3vB042V8gIp6MiGfzw08CV5b9bN93HIiITkR0ZmdnK6m4VcML6MwmS8qAcR9wsaSLJJ0H7AEO9xeQ9Mq+w+uA/52/vxu4RtIFki4ArsnPWYs4663ZZEk26B0RZyTdSPagnwEORsRDkm4GuhFxGPjXkq4DzgCngHfknz0l6QNkQQfg5t4AuLWHF9CZTRanBrFknKLDbPw5NYiNBS+gM5ssDhiWTNuz3prZi3lPb0vKC+jMJodbGGZmVooDhlXCC/TMJp+7pGzDegv0emsuegv0wN1RZpPELQwbyaCWhBfomU0HtzCstKKWxMpg0eMFemaTxS0MK62oJTEzM7i8d7gzmywOGFZaUYvh7Fkv0DObBg4YVlpRi6G3IM8L9MwmmwOGlTYs1ce+fVl+qOefz346WJhNHgcMK82pPsymmwPGFKhyUZ1bEmbTy9NqJ5wX1ZlZVdzCaKFRWgxeVGdmVXELo2VGbTF41zszq4pbGC0zrMUwqOVRNBXWi+rMbFQOGBuQOkProO8vahn0WhonTkDEC8fXXutFdWZWjaQBQ9IuSY9IOibppgHXf1nSw5IekPSHkub7rp2VdH/+OpyynuvR6xpa+YCuKmgUff+WLYPLz8wMbnncdVfxVFinJDezUSgi0nyxNAN8C3gTsATcB+yNiIf7ylwNfCMiTkt6N3BVRPyL/Nr3IuLlo/zOTqcT3W63snsYZmEhe4ivND+fTTdN9f1bt8Izz7w4OGzeXJwAUMqmwK60ciyk9z1eV2E2XSQdjYhOmbIpWxg7gWMRcTwingPuBHb3F4iIr0RE75F1BNiesD6VSj2YXPQ9p04NbjHMzw8uXzRW4dlTZjaqlAHjQuDxvuOl/FyRG4Av9B2/VFJX0hFJP5mighuRejB52PcPWjw3LG3HIJ49ZWajShkwNODcwP4vSW8HOsCH+k7P5c2ktwEfkfTDBZ/dnweW7vLy8kbrXNqoD+jU3z9q2g7PnjKzUaUMGEvAjr7j7cDJlYUk/QSwCFwXEc/2zkfEyfznceCrwOWDfklEHIiITkR0Zmdnq6v9GqrMqzRo8Hk93z9K2o7UAc/MJlBEJHmRLQo8DlwEnAd8E3j1ijKXA38BXLzi/AXAS/L324BvA5es9TuvvPLKGAef+UzE/HyElP38zGeGl928OSKbC5W9Nm8e/pkm6mlmkwnoRsnnerJZUgCSrgU+AswAByPiFkk35xU8LOkPgH8AfDf/yGMRcZ2kfwh8AnierBX0kYj41Fq/r85ZUkVGnX2UeraVmdkwo8ySShow6jYOAWPUALBpU9auWKloOqyZWZXGZVrtVBp19pEHn82sLRwwKjZqAPDgs5m1hQNGxVJPhzUza4rTm1es96BfXMy6oebmXtjzethnHCDMbNw5YCTgAGBmk8hdUmZmVooDhpmZleKAUSPvP2FmbeYxjJqMuhe3mdm4cQujJt5/wszazgGjJt5/wszazgGjJk4BYmZt54BRE6cAMbO2c8CoiVOAmFnbOWDk6pjyOsqOeGZm48bTavGUVzOzMtzCwFNezczKcMDAU17NzMpwwMBTXs3MykgaMCTtkvSIpGOSbhpw/SWSPptf/4akhb5rv5Kff0TSm1PWc60pr84BZWaWMGBImgE+BrwFuATYK+mSFcVuAJ6KiL8HfBj49/lnLwH2AK8GdgG35t+XxLApr70B8RMnIOKFAXEHDTObNilbGDuBYxFxPCKeA+4Edq8osxu4PX//OeCNkpSfvzMino2I7wDH8u9LpmjKqwfEzcwyKQPGhcDjfcdL+bmBZSLiDPA0sLXkZwGQtF9SV1J3eXm5oqq/wAPiZmaZlAFDA85FyTJlPpudjDgQEZ2I6MzOzo5YxbV5QNzMLJMyYCwBO/qOtwMni8pIOgd4BXCq5Gdr4RxQZmaZlAHjPuBiSRdJOo9sEPvwijKHgevz928FvhwRkZ/fk8+iugi4GPjThHUt5BxQZmaZZKlBIuKMpBuBu4EZ4GBEPCTpZqAbEYeBTwGflnSMrGWxJ//sQ5J+B3gYOAP8QkScTVXXtezb5wBhZqbsH/STodPpRLfbbboaZmatIeloRHTKlPVKbzMzK8UBw8zMSnHAMDOzUhwwzMysFAcMMzMrxQHDzMxKmahptZKWgRMb/JptwF9WUJ228P1ONt/v5KrqXucjolRepYkKGFWQ1C07J3kS+H4nm+93cjVxr+6SMjOzUhwwzMysFAeM1Q40XYGa+X4nm+93ctV+rx7DMDOzUtzCMDOzUqY2YEjaJekRScck3TTg+kskfTa//g1JC/XXsjol7veXJT0s6QFJfyhpvol6VmWt++0r91ZJIanVM2vK3K+kf57/GT8k6bfrrmNVSvxdnpP0FUl/lv99vraJelZF0kFJT0h6sOC6JP1a/t/jAUlXJKtMREzdi2x/jr8A/i5wHvBN4JIVZX4e+Hj+fg/w2abrnfh+rwY25+/fPen3m5c7H7gXOAJ0mq534j/fi4E/Ay7Ij3+w6XonvNcDwLvz95cAjzZd7w3e8+uBK4AHC65fC3yBbGvr1wHfSFWXaW1h7ASORcTxiHgOuBPYvaLMbuD2/P3ngDdKGrTXeBuseb8R8ZWIOJ0fHiHbFretyvz5AnwA+A/A/6uzcgmUud+fAz4WEU8BRMQTNdexKmXuNYC/lb9/BQ1t71yViLiXbIO5IruB34rMEeAHJL0yRV2mNWBcCDzed7yUnxtYJiLOAE8DW2upXfXK3G+/G8j+xdJWa96vpMuBHRHx+3VWLJEyf74/AvyIpK9JOiJpV221q1aZe/13wNslLQF3Ae+pp2qNGfX/73VLtkXrmBvUUlg5XaxMmbYofS+S3g50gDckrVFaQ+9X0ibgw8A76qpQYmX+fM8h65a6iqz1+EeSLo2Iv0pct6qVude9wG0R8Z8k/TjZNtCXRsTz6avXiNqeVdPawlgCdvQdb2d1s/Vvykg6h6xpO6xZOM7K3C+SfgJYBK6LiGdrqlsKa93v+cClwFclPUrW73u4xQPfZf8+/15EfD8ivgM8QhZA2qbMvd4A/A5ARHwdeClZ3qVJVer/7ypMa8C4D7hY0kWSziMb1D68osxh4Pr8/VuBL0c+wtRCa95v3kXzCbJg0db+7Z6h9xsRT0fEtohYiIgFsjGb6yKirRvCl/n7/LtkExuQtI2si+p4rbWsRpl7fQx4I4CkHyULGMu11rJeh4GfzWdLvQ54OiK+m+IXTWWXVESckXQjcDfZrIuDEfGQpJuBbkQcBj5F1pQ9Rtay2NNcjTem5P1+CHg58F/ysf3HIuK6xiq9ASXvd2KUvN+7gWskPQycBd4bEU82V+v1KXmv/xb4pKRfIuuaeUeL/7GHpDvIuhK35eMy7wfOBYiIj5ON01wLHANOA/8yWV1a/N/RzMxqNK1dUmZmNiIHDDMzK8UBw8zMSnHAMDOzUhwwzMysFAcMs0Qk7ZD0HUlb8uML8uNWZwK26eWAYZZIRDwO/CbwwfzUB4EDEXGiuVqZrZ/XYZglJOlc4ChwkCxj7OV5llWz1pnKld5mdYmI70t6L/BF4BoHC2szd0mZpfcW4LtkCQ/NWssBwywhSZcBbyLLiPtLqTa2MauDA4ZZIvkOjb8J/JuIeIwsweN/bLZWZuvngGGWzs+RZf29Jz++Ffj7ktq8OZVNMc+SMjOzUtzCMDOzUhwwzMysFAcMMzMrxQHDzMxKccAwM7NSHDDMzKwUBwwzMyvFAcPMzEr5/zwPx+2TJnGHAAAAAElFTkSuQmCC\n", "text/plain": [ "