{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear basis function models with PyMC3" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.8\n" ] } ], "source": [ "import logging\n", "import pymc3 as pm\n", "import numpy as np\n", "\n", "print(pm.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a [PyMC3](https://docs.pymc.io/) implementation of the examples in [Bayesian regression with linear basis function models](https://nbviewer.jupyter.org/github/krasserm/bayesian-machine-learning/blob/dev/bayesian-linear-regression/bayesian_linear_regression.ipynb). To recap, a linear regression model is a linear function of the parameters but not necessarily of the input. Input $x$ can be expanded with a set of non-linear basis functions $\\phi_j(x)$, where $(\\phi_1(x), \\dots, \\phi_M(x))^T = \\boldsymbol\\phi(x)$, for modeling a non-linear relationship between input $x$ and a function value $y$.\n", "\n", "$$\n", "y(x, \\mathbf{w}) = w_0 + \\sum_{j=1}^{M}{w_j \\phi_j(x)} = w_0 + \\mathbf{w}_{1:}^T \\boldsymbol\\phi(x) \\tag{1}\n", "$$\n", "\n", "For simplicity I'm using a scalar input $x$ here. Target variable $t$ is given by the deterministic function $y(x, \\mathbf{w})$ and Gaussian noise $\\epsilon$.\n", "\n", "$$\n", "t = y(x, \\mathbf{w}) + \\epsilon \\tag{2}\n", "$$\n", "\n", "Here, we can choose between polynomial and Gaussian basis functions for expanding input $x$. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from functools import partial\n", "from scipy.stats import norm\n", "\n", "def polynomial_basis(x, power):\n", " return x ** power\n", "\n", "def gaussian_basis(x, mu, sigma):\n", " return norm(loc=mu, scale=sigma).pdf(x).astype(np.float32)\n", "\n", "def _expand(x, bf, bf_args):\n", " return np.stack([bf(x, bf_arg) for bf_arg in bf_args], axis=1)\n", "\n", "def expand_polynomial(x, degree=3):\n", " return _expand(x, bf=polynomial_basis, bf_args=range(1, degree + 1))\n", "\n", "def expand_gaussian(x, mus=np.linspace(0, 1, 9), sigma=0.3):\n", " return _expand(x, bf=partial(gaussian_basis, sigma=sigma), bf_args=mus)\n", "\n", "# Choose between polynomial and Gaussian expansion\n", "# (by switching the comment on the following two lines)\n", "expand = expand_polynomial\n", "#expand = expand_gaussian" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, to expand two input values [0.5, 1.5] into a polynomial design matrix of degree 3 we can use" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.5 , 0.25 , 0.125],\n", " [1.5 , 2.25 , 3.375]])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "expand_polynomial(np.array([0.5, 1.5]), degree=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The power of 0 is omitted here and covered by a $w_0$ in the model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example dataset\n", "\n", "The example dataset consists of N noisy samples from a sinusoidal function f." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8dfJQjJsBiSWTTOsFlQ2CSBERREFXCi41BQFbSSigkKpEhWQRTYpIqsKJsjSRvzRymL9It9QFhGoIJCIrAETTKFmZAkEJuuc3x8JfBFDCGRmzszcz/PxyMPMncu970PCxzPnnnuu0lojhBAi8AWZDiCEEMI7pOALIYRFSMEXQgiLkIIvhBAWIQVfCCEsIsR0gMupU6eOttvtpmMIIYRf+fbbb3/WWkeW9Z7PFny73c727dtNxxBCCL+ilMq83HsypCOEEBYhBV8IISxCCr4QQliEFHwhhLAIKfhCCGERUvCFEMIipOBbgMPhYNu2bTgcDtNRhBAGScEPcMnJyURFRdG9e3eioqJITk42HUkIYYjy1fXw27dvr+XGq8pxOBxERUXhdDovbLPZbGRmZhIZWeaNeEIIP6eU+lZr3b6s99zSw1dKJSmlspVSuy/zflelVI5Salfp12h3nFeULyMjgypVqvxiW2hoKBkZGWYCCSGMctfSCh8Ds4FF5ezzldb6ITedT1SA3W6noKDgF9sKCwuRNYqEsCa3FHyt9UallN0dxxLuExkZSWJiInFxcYSGhlJYWMjMmTPZvn07+/fv56677qJdu3YcPHiQ/v37U1hYSO3atfnNb35D3bp1iY2NpV27dqabIYRwE29etL1DKZWqlPofpdQtZe2glIpXSm1XSm2XGSXuERsbS1paGn379qVFixY8//zz9OrVi2HDhpGSkgKUDPPUrFmTG264gdOnT7Np0yZmzZp1YegnNTWVJ554gk8++YQzZ84YbI0QojLcdtG2tIf/udb61jLeqwm4tNa5SqlewAytdbPyjicXbSvnzJkz7N27lw4dOlBQUIDdbqdx48Z069aNrl27cuutt1KnTh2UUmX++eLiYlwuF6GhoaxYsYL4+Hiys7MJCwvj97//Pa+88or0/oXwQeVdtPVKwS9j3wygvdb658vtIwX/2uTk5DBr1iymT59OjRo1OHz4MEFBQRQUFPzqAu7VKC4uZsuWLSQnJ7Nw4UIAjh07Ro0aNdwVXQjhBh6fpVOBAHVVaVdSKdWh9LzHvXFuqygsLOTdd9/FbrczatQounTpwqeffkpQUMmPuDLFHiA4OJiYmBjmzJlDVlYWK1asoEaNGmitee211/juu+/c0QwhhAe5a1pmMrAFuFkplaWUilNKDVJKDSrd5TFgt1IqFZgJPKl99QYAP7VmzRqGDx/OHXfcwY4dO1i5ciUdOnTwyLkiIiLo1q0bAIcPH2b+/Pm0bt2al19+mbNnz3rknEKIypMbr/zYuXPn2LZtG3fffTdaazZt2kRMTMxlx+U95cSJE7z11lvMnj2bRo0akZiYyD333OPVDEKIEsaHdIT7paam0r59e3r16oXD4UApxZ133un1Yg9Qu3ZtZs2axcaNGwkODqZ///6cO3fO6zmEEOWTgu9hnli47MMPP6RDhw6cOnWK5cuX+8wyCXfeeSepqamkpKRQtWpVioqK+O9//2s6lhCilBR8D3L3wmVaa1588UUGDRpEt27dSEtLo3v37m5K6x5Vq1bl5ptvBuCtt96idevWbNy40XAqIQRIwfcYh8NBXFwcTqeTnJwcnE4ncXFxlerpK6WoXr06r732GqtWraJOnTpuTOx+Tz31FLVq1aJ79+787W9/Mx1HCMtz11o64hLnFy67eKXK8wuXXe0QzI8//sjJkydp1aoVU6ZMMTJOfy1atGjBli1b+N3vfke/fv04cuQII0aM8Jv8QgQa6eF7iLsWLtu/fz9dunTh8ccfp6ioyO+KZa1atVizZg2xsbGMHTuWAwcOmI4khGVJwfeQ8wuX2Ww2atasic1mIzEx8ap69zt37iQmJob8/HyWLl1KSIh/fiALCwtjyZIl7Nix48L4vq9OBxYikMk8fA9zOBxkZGRgt9uvqtinpqZy7733Ur16dVJSUmjWrNylh/xKYmIiX3/9NfPnzyc4ONh0HCECSnnz8P2zy+hHIiMjr2na5OTJk6latSrr16+nUaNGHkhmztGjR1mwYAEul4ukpKQLyz8IITxLCr6PSkpKIjs7m6ioKNNR3G7UqFG4XC7GjBlDREQE06dP97trE0L4I+la+ZCjR48SGxvLiRMnsNlsAVnszxs9ejRDhw5lxowZjB8/3nQcISxBevg+Iicnh549e3L48GESEhKoXbu26UgepZRi2rRpnDp1itDQUNNxhLAEKfg+ID8/nz59+rBnzx6++OILWrdubTqSVwQFBZGUlHRhOKeya/YLIconQzqGaa157rnnWLduHUlJST63VIKnnS/2W7dupWnTpuzcudNwIiEClxR8wxwOB5s2bWLcuHE8/fTTpuMYc9NNNwHw0EMPkZWVZTiNEIFJ5uH7gFOnTnHddddZfqZKWloaMTExNG/enK+++gqbzWY6khB+R9bD90G7d+9m6NChFBQUEBERYYlif6Wlolu1asVf//pXvv32W1544QW5G1cIN5OCb0BOTg59+vRh6dKlHD9ujUf7VnSp6IcffpgxY8aQl5dHYWGhl1MKEdhkSMfLtNb07duXVatWsX79emJiYkxH8jiHw0FUVNQvVg612WxkZmaWeReyy+VCKWWJTz1CuJsM6fiQadOmsXz5cqZOnWqJYg//t1T0xc4vFV2WoKAglFIcOHCAnj17kp2d7YWUQgQ+KfhedPLkSd5++20ee+wxhg4dajqO11zrUtF5eXmsX7+ep59+GpfL5cGEQliDFHwvqlWrFlu2bCExMdFSwxXXulR0q1atmDlzJmvWrGHy5MleSitE4JIxfC/QWrNmzRruv/9+SxX6S13LUtFaa/r168fSpUtZt24dd911l4dTCuHfZAzfsLlz59KjRw+++OIL01GMioyMJDo6+qqWi1ZK8eGHH9KkSROmTJniwXRCBD5ZS8fD0tLSGD58OL169aJXr16m4/ilGjVqsHr1aurXr286ihB+TXr4HpSXl0e/fv2IiIhgwYIFlh7OqazGjRsTHh7O6dOn2bRpk+k4QvglKfgeNHr0aHbv3s2CBQu44YYbTMcJCIMGDeLBBx+87JROIcTlScH3oC5duvD666/Ts2dP01ECxoQJE9Ba079/f5mqKcRVkoLvQb1792bixImmYwSURo0aMXPmTL766itmzZplOo4QfkUKvgcMGTKEqVOnmo4RsAYMGMBDDz3E66+/zoEDB0zHEcJvSMF3s5SUFGbPnn3ZFSFF5Z2fqtmzZ0/Cw8NNxxHCb8iNV26Um5vLbbfdRpUqVdi1a5es5y6E8Dq58cpL3njjDTIzM0lKSpJi7yVHjx7l0UcfJT093XQUIXyeFHw3yczM5P3332fIkCF06dLFdBzL0FqTkpJCfHy8PDBFiCuQgu8mUVFRbN68mQkTJpiOYikNGjRg6tSpFx4CL4S4PLcUfKVUklIqWym1+zLvK6XUTKVUulIqTSnVzh3n9RXnH7odHR1N9erVDaexnueee467776b4cOHc+zYMdNxhPBZ7urhfwz0KOf9nkCz0q944H03nde4PXv20LRpUxYsWGA6imUFBQUxb9488vLyGDVqlOk4QvgstyyeprXeqJSyl7NLb2CRLhlk3aqUilBK1dNa+3V3zOVy8fzzz1OtWjUeeugh03EsrXnz5ixfvpyOHTuWu9+1LNEsRKDw1hh+A+DHi15nlW77BaVUvFJqu1Jquz/MY//oo4/YtGkTf/nLX6R4+IAePXpQq1YtioqKyMvL+9X7FX2QuhCBylsFv6xlIn81pUJrPU9r3V5r3d7XC2h2djYjRoyga9euPPPMM6bjiFJOp5Po6GjGjBnzi+0Oh4O4uDicTic5OTk4nU7i4uLkBjlhKd4q+FnAjRe9bggc9dK5PWLHjh0opZg7d64se+xDbDYbbdu2Zdq0aeze/X9zCK72QepCBCJvFfyVQP/S2TqdgBx/H7/v0aMHP/74Iy1atDAdRVzinXfe4brrrmPQoEEX5uZf64PUhQgk7pqWmQxsAW5WSmUppeKUUoOUUoNKd/kCOAykA/OBF91xXhMKCwtZvnw5WmuqVatmOo4oQ506dZgyZQpff/01ixcvBq79QepCBBJZS+cqTZs2jT//+c9s2LBBHqjtw1wuF507d6ZatWqsXbv2wnaZpSMCXXlr6UjBvwpHjx7l5ptvpmvXrqxatcp0HHEFR48eJTIyktDQUNNRhPAaWTzNTUaMGEFhYSHvvfee6SiiAurXr09oaCinT58mMzPTdBwhjHPLjVdW8PXXX7NkyRLefPNNmjRpYjqOqCCtNTExMURERLBhwwaZUSUsTXr4FXTu3LkLz6gV/kMpxZAhQ/jqq69YunSp6ThCGCVj+CLgFRcXEx0djcPhYN++fTK7SgQ0GcOvhBMnTvCXv/ylzFv1hX8IDg5m5syZZGVlMWXKFNNxhDBGCv4VjBkzhhEjRnDw4EHTUUQlxMTEEBsbS1pamjwoRViWXLQtx549e5g7dy7x8fHcdtttpuOISkpKSiIsLEwu3ArLkh7+ZWitGTZsGDVq1GDcuHGm4wg3CA8PRylFZmYmaWlppuMI4XXSw7+Mf/7zn6xZs4bp06fLHZkBxOVycf/991OtWjW2bdtGcHCw6UhCeI308C+jXr16/OEPf+Cll14yHUW4UVBQEGPGjGHnzp0sXLjQdBwhvEqmZQrL0VrTuXNnfvjhBw4ePEiNGjVMRxLCbWRa5lU4ceIEQ4cOJTs723QU4SFKKd577z1++uknJk+ebDqOEEDJwn7btm3z6EN5pOBfYvz48cyaNYuffvrJdBThQR07dmTAgAHk5uaajiKE1x6/KUM6Fzlw4AC33HILzz77LPPmzfPquYX3aa1liqYwzuFwEBUVhdPpvLDNZrORmZl5TRNGZEingkaMGEF4eDjjx483HUV4wfliv2XLFlJTUw2nEVblzcdvyrTMUhs2bGD58uVMmDCB3/zmN6bjCC/Jy8ujT58+NGvWjI0bN0qPX3idNx+/KT38Una7nZdeeolhw4aZjiK8KDw8nLFjx7Jp0yY+++wz03GEBXnz8Zsyhi8sr6ioiNatW5Ofn8+ePXt+9fFaCG9w1+M3ZQy/HHl5eTzzzDPs3bvXdBRhSEhICFOnTuXQoUN88MEHpuMIi4qMjCQ6Otqjd/ZbvuDPnDmThQsXcuzYMdNRhEE9e/akT58+stSCCGiWHtL5+eefadKkCXfeeSeff/65R88lhBDeIEM6l/H222+Tm5srD8UQF2itWbp0KUeOHDEdRQi3s2zBP3ToEHPnzuWPf/wjt9xyi+k4wkccPXqUAQMGMHLkSNNRhHA7yxb8yMhIRowYwdixY01HET6kQYMGDB06lCVLlrBr1y7TcYRwq4Aew3fXNCdhLadOnaJJkya0b9+eL7/80nQcIa6KJcfwL7cYkdaa+Ph41qxZYzih8FURERGMHDmSNWvWyO+JCCgBWfAdDgdxcXE4nU5ycnJwOp3ExcXhcDj45z//yfz580lPTzcdU/iwF198kXvvvReXy2U6ihBuE5Br6ZxfjOji1edCQ0M5dOgQCQkJNG/enIEDBxpMKHxdWFgYa9euNR1DCLcKyB7+5RYj+ve//83333/PxIkTCQ0NNZRO+BOn08mMGTPIy8szHUWISgvIgl/WYkRz585l6tSpdOrUib59+5qOKPzE1q1bGTp0KHPnzjUdRYhKs8wsnVq1avHxxx/TsmVLOnfu7KaUwgp69OjBtm3bOHToEBEREabjCFGu8mbpBHTBF8Iddu3aRbt27XjttdfkGbjC51lyWubFJkyYIB/JxTVr06YN/fr1Y8aMGWRlZZmOI8Q1C/iCn5GRwbhx4/j2229NRxF+bPz48XTu3JkzZ86YjiLENXNLwVdK9VBK7VdKpSulEsp4/xmllEMptav06zl3nLciRo8eTVBQkCyhICrFbrezdu1aWrRoYTqKENes0gVfKRUMzAF6Ai2BWKVUyzJ2Xaq1blP69VFlz1sRqampLFmyhFdeeYWGDRt645QiwB07dow5c+aYjiHENXFHD78DkK61Pqy1LgA+AXq74biVlpCQQEREBAkJv/rQIcQ1WbRoEYMHD2bjxo2mowhx1dxR8BsAP170Oqt026UeVUqlKaWWKaVuLOtASql4pdR2pdR2h8NR6WDDhw9nzpw5MpVOuM2QIUOoX78+I0aMwFdnuAn/k5eXR0xMDKtWrfLoedxR8FUZ2y79l7AKsGutWwEpwMKyDqS1nqe1bq+1bu+O1S3vu+8+YmNjK30cIc6rWrUqY8eOZevWrXz22Wem44gAMXfuXL7++muqVq3q0fNUeh6+UuoOYIzW+oHS168DaK0nXWb/YOCE1vq68o4r8/CFryoqKqJVq1YUFxeze/duWaZDVMr55bijo6NZvXp1pY/n6Xn424BmSqlGSqkqwJPAyksC1Lvo5SPAXjecVwgjQkJCmDx5MtHR0Zw9e9Z0HOHn/vvf/9KoUSOvPGrVLXfaKqV6Ae8BwUCS1nqCUmocsF1rvVIpNYmSQl8EnABe0FrvK++Y0sMXQliF1hqlyhodv3qytIIQHpKamsq+ffv4/e9/bzqK8EOffvop9913H7Vr13bbMS2/tIIQnjJ27FgGDhyIO2aVCWv5/vvviY2N9cpQznlS8IWohIkTJ3L27Fnefvtt01GEn0lISKB69eq8+uqrXjunFHwhKuG3v/0tcXFxvP/++xw6dMh0HOEnNm7cyOeff05CQgJ16tTx2nllDF+ISjp27BhNmzbl4Ycf5pNPPjEdR/g4rTV33HEHWVlZHDhwwO1z78sbww/IZ9oK4U316tVj5MiROJ1Ot862EIHpzJkz1K1bl4EDB3r8RqtLSQ9fCCEM8FTnQGbpCOEFWmtWrVrFhg0bTEcRPmr16tUcPHgQwMgnQRnSEcJNioqK+NOf/kRYWBipqakEBwebjiR8yOnTp+nfvz+tW7fmf//3f41kkB6+EG4SGhrKpEmT+P7771m4sMz1AYWFTZ06FYfDwcSJE41lkDF8IdxIa03nzp05cuQIBw4coFq1aqYjCR/wn//8h2bNmtG7d2+Sk5M9ei4ZwxfCS5RSTJ06laNHjzJ9+nTTcYSPeOuttygqKmLChAlGc0jBF8LNYmJieOGFF7Db7aajCB+gtaZWrVoMGzaMxo0bG80iQzpCCBFAZEhHCAMKCgqYOXMm+/aVuxK4CGCbN29m3bp1pmNcID18ITzE4XDQpEkTunbtysqVK6/8B0RAKS4upl27dpw5c4b9+/d77clo0sMXwoDIyEjeeOMNVq1a5VO9POEdixYtIi0tjUmTJvnMYzClhy+EB+Xl5XHzzTdz/fXXs337doKCpI9lBWfPnqV58+Y0bNiQrVu3evWuWunhC2FIeHg4kyZNYufOnSxevNh0HOEl56fmvvvuuz61mJ4srSCEhz355JOkpKQYn5InvKd+/foMGjSILl26mI7yCzKkI4QQAUSGdITwASdPnuTPf/4zWVlZpqMID9m5cydJSUkUFxebjlImKfhCeMmpU6eYPXs2b7zxhukowgO01gwdOpSEhARyc3NNxymTFHwhvKRRo0YMGzaMxYsX880335iOI9zs73//Oxs3bmTcuHFcd911puOUScbwhfCi06dP07x5c5o0acKmTZt8agaHuHZ5eXm0aNGCGjVqsGPHDkJCzM2HkTF8IXxEzZo1mTBhAps3b2bixIk4HA7TkYQbvPvuu2RkZPDee+8ZLfZXIgVfCC8LDw8nODiYyZMnExUV5fH10YXntW3blmHDhnHvvfeajlIuGdIRwoscDgdRUVE4nc4L22w2G5mZmURGRhpMJgKFDOkI4SMyMjKoUqXKL7YFBweTkZFhJpColK1btzJy5EjOnTtnOkqFSMEXwovsdjsFBQW/2Hbu3Dl5WIofKi4uZvDgwSxYsACXy2U6ToVIwRfCiyIjI0lMTMRms1GzZk1CQkJwuVx89913pqOJq5SUlMS3337L1KlTqV69uuk4FSJj+EIY4HA4yMjIoG7dunTt2pWwsDBSU1N9Zhldqzn/87Db7RW6lnLy5EmaN29OixYt2LBhg09Nr5UxfCF8TGRkJNHR0dx4443MmDGDvXv3MmfOHNOxLCk5OZmoqCi6d+9e4VlTo0eP5sSJE8yaNcuniv2VSA9fCB/wwQcf8OSTTxIREWE6iqVc66yptLQ0Nm7cyODBg70R86qU18P33TsEhLCQQYMGAeByueQhKV50ftbUxQU/NDSUjIyMcgt+q1ataNWqlTciupX8ZgnhIw4fPky7du3417/+ZTqKZZQ1a6qwsPCys6Y++ugj/vCHP3D27FkvpHM/txR8pVQPpdR+pVS6UiqhjPfDlFJLS9//t1LK7o7zChFI6tWrx5kzZ3jhhRfIz883HccSLp01ZbPZSExMLLN3//PPPzNixAiysrKoWrWqgbSVV+mCr5QKBuYAPYGWQKxSquUlu8UBJ7XWTYHpwJTKnleIQGOz2ZgzZw4HDhzgnXfeMR3HMmJjY8nMzCQlJYXMzExiY2PL3C8hIYGcnBzmzp3rVxdqL+aOHn4HIF1rfVhrXQB8AvS+ZJ/ewMLS75cB3ZS//o0J4UE9evTgiSeeYMKECaSnp5uOYxnnZ01dbtx+06ZNJCYmMmzYMG699VYvp3MfdxT8BsCPF73OKt1W5j5a6yIgB7jeDecWIuBMnz6dsLAw3n33XdNRRKmEhASioqIYM2aM6SiV4o5ZOmX11C+d61mRfVBKxQPxADfddFPlkwnhh+rXr8+6deu47bbbTEcRpZYtW0ZWVhbVqlUzHaVS3NHDzwJuvOh1Q+Do5fZRSoUA1wEnLj2Q1nqe1rq91rq9rBworKxdu3aEhoZy6tQpTp48aTqOZZ08eRKXy0XdunVp377Mqe1+xR0FfxvQTCnVSClVBXgSWHnJPiuBAaXfPwb8S/vqHV9C+Ii8vDzatGnD0KFDTUexJK01ffv2pXfvSy9J+q9KF/zSMfnBwJfAXuBTrfX3SqlxSqlHSndLBK5XSqUDfwJ+NXVTCPFL4eHhPP300yxatIgvv/zSdBzLmT9/PuvXr+fhhx82HcVtZGkFIXxYfn4+bdq0wel0snv3br9ZldHf/fjjj9xyyy1ER0eTkpLiV9MwZfE0IfxUWFgY8+fPJzMzk9dff910HEvQWvP8889TXFzM/Pnz/arYX4kUfCF8XExMDC+//DKHDx+mqKjIdJyA99NPP7Fv3z4mTZpE48aNTcdxKxnSEcIPFBYWEhISElC9TV929uzZCw+b9zcypCOEnwsNDUUpRUZGhqyb7yFaaz788EOcTifVqlXzy2J/JVLwhfAj77//PoMHDyYlJcV0lIAzb948Bg0axKeffmo6isfIkI4QfsTpdNK2bVtyc3NJS0ujdu3apiMFhPT0dFq3bk2XLl1YvXq1Xz+TQIZ0hAgQNpuNv/71r/z00088//zz+GqHzZ8UFRUxYMAAQkNDSUpK8utifyWB2zIhAtTtt9/O+PHjWbZsGUuXLjUdx+9NnjyZzZs3M2fOHBo2bGg6jkfJIw6F8EOvvvoqQUFBPPjgg6aj+L3HH38cl8tFv379TEfxOBnDF8LPnTt3jpCQEKpUqWI6il/Jz8+nSpUqATfVVcbwhQhQZ8+epWPHjiQkyPJUV0NrzVNPPUX//v0tdR1ECr4QfqxatWrcc889TJ8+nRUrVpiO4zfmzJnDsmXLuO222wKuh18eGdIRws/l5+fTpUsXDh06xM6dO7Hb7RX6cw6Hg4yMDOx2+2Uf7ReItm7dyl133cUDDzzAihUrAm5WjgzpCBHAwsLCWLp0KS6Xi8ceewyn03nFP5OcnExUVBTdu3cnKiqK5ORkLyQ1z+Fw8Pjjj9OwYUMWLVoUcMX+SqzVWiECVJMmTVi8eDEAp06dKndfh8NBXFwcTqeTnJwcnE4ncXFxOBwOb0Q16tChQ7hcLpYtW0atWrVMx/E6mZYpRIB45JFHePDBB6+4BkxGRgZVqlT5xSeB0NBQMjIyAn5op1OnThw+fJiwsDDTUYyQHr4QASQ4OJjTp0/z+OOPs3HjxjL3sdvtFBQU/GJbYWFhhcf+/dGiRYuYNGkSWmvLFnuQgi9EwNFa891339G3b1/S09N/9X5kZCSJiYnYbDZq1qyJzWYjMTExYHv3W7ZsYeDAgaSkpFBcXGw6jlEyS0eIAJSenk6nTp2oXbs2W7Zs4frrr//VPlaYpZOZmUnHjh2pXr0633zzjSUWm5NZOkJYTNOmTVm+fDmZmZn07duX/Pz8X+0TGRlJdHR0wBb748eP06NHD/Lz81m5cqUliv2VSMEXIkDFxMTw8ccfs2/fPo4cOWI6jtdt2rSJI0eOsGLFClq2bGk6jk+QIR0hAtypU6eIiIgASsb3rXRnaXZ2NjfccIPpGF4lQzpCWFhERARaa15//XVGjhxpOo5HFRcXM3DgQFatWgVguWJ/JVLwhbCI48ePM3HiRKZMmWI6ike4XC7i4+P56KOP+O6770zH8Uly45UQFqCU4v333yc3N5eEhAQKCwsDqrevteall14iKSmJUaNG8cYbb5iO5JOk4AthEcHBwSxatIiQkBBGjRpFUVERY8aMMR2r0oqLi4mPjycpKYkRI0YwduxY05F8lhR8ISwkJCSEBQsWEB4eTosWLUzHcYugoCDCw8MZPXo0Y8aMsdRF6asls3SEsLgNGzbQsWNHwsPDTUe5KqdPn8bhcNCkSRPLzT4qj8zSEUKU6T//+Q8PPPAA3bt35/jx46bjVNj5O4m7d+9Ofn6+FPsKkoIvhIU1aNCARYsWsW3bNqKjo9m5c6fpSFe0du1aOnToQHZ2NomJiZZeDO1qScEXwuKeeOIJ1q9fT2FhIXfccQcfffSR6UhlKi4uZuLEiTzwwAM0aNCAb775hnvuucd0LL8iBV8IQadOnRY3iCsAAAhcSURBVNixYwd33XUXmZmZpuOUSSnF2rVreeyxx9i8eTONGzc2HcnvyEVbIcQF55cPDg4OJiUlhdzcXH73u98Zy+NyuVi4cCHdunXjpptu4ty5c9hsNhmzL4dctBVCVEhwcPCFJ2bNmDGDPn368Oijj5a5rr6n7dq1i5iYGP74xz8yb948AKpWrSrFvhKk4AshyvSPf/yDt99+m9WrV9OyZUteeeUVfv75Z4+f9+DBgzz77LPcfvvtpKens2DBAsaNG+fx81qBFHwhRJlCQ0N58803SU9P59lnn2X27NmsX78eKFnKwJ201heOOXv2bD755BNefvll9u/fzzPPPENQkJQqd6jUGL5SqjawFLADGcATWuuTZexXDJxfzeiI1vqRKx1bxvCF8C379u2jadOmhISEMH78+AsXUPv06UODBg2u6Zg//PADn332GQsWLGDatGncf//9OBwOiouLqVu3rptbYA2eHMNPANZqrZsBa0tfl8WptW5T+nXFYi+E8D2//e1vCQkpWY2lTp06ZGdnM2TIEBo2bEjbtm1/sWDZsWPHOHfuHABFRUXk5uaSm5sLQEFBAfHx8TRp0oTGjRszfPhwqlWrdqEXHxkZKcXeQyrbw98PdNVaH1NK1QPWa61vLmO/XK119as5tvTwhfB9e/fu5R//+Afr16+nbt26LF68GIBGjRqRkZFBUFAQLpcLgH79+rFkyRK01jRq1Ig2bdpw3333cf/999O8eXOTzQgo5fXwK1vwT2mtIy56fVJrXauM/YqAXUARMFlrvfwyx4sH4gFuuumm2311PrAQonxLliwhKyuL3NxcwsLCCA8Pp3nz5vTu3Ruw3pO3vKlSBV8plQKU9fnqTWBhBQt+fa31UaVUY+BfQDet9aHyzis9fCGEuHrlFfwrLo+stb6vnAP/pJSqd9GQTvZljnG09L+HlVLrgbZAuQVfCCGEe1X2ou1KYEDp9wOAFZfuoJSqpZQKK/2+DtAF2FPJ8wohhLhKlS34k4HuSqmDQPfS1yil2iulzq/A1ALYrpRKBdZRMoYvBV8IIbysUk+80lofB7qVsX078Fzp95uB2ypzHiGEEJUnt68JIYRFSMEXQgiLkIIvhBAWIQVfCCEsQgq+EEJYhBR8IYSwCCn4QghhEVLwhRDGOBwOtm3bhsPhMB3FEqTgCyGMSE5OJioqiu7duxMVFUVycrLpSAGvUssje5KslilE4HI4HERFReF0Oi9ss9lsZGZmEhkZaTCZ//PkE6+EEOKqZWRkUKVKlV9sCw0NJSMjw0wgi5CCL4TwOrvdTkFBwS+2FRYWYrfbzQSyCCn4Qgivi4yMJDExEZvNRs2aNbHZbCQmJspwjodVarVMIYS4VrGxsdx3331kZGRgt9ul2HuBFHwhhDGRkZFS6L1IhnSEEMIipOALIYRFSMEXQgiLkIIvhBAWIQVfCCEsQgq+EEJYhM+upaOUcgCZbjhUHeBnNxzHX0h7A5u0N7C5o71RWusy57r6bMF3F6XU9sstJBSIpL2BTdob2DzdXhnSEUIIi5CCL4QQFmGFgj/PdAAvk/YGNmlvYPNoewN+DF8IIUQJK/TwhRBCIAVfCCEsIyAKvlKqh1Jqv1IqXSmVUMb7YUqppaXv/1spZfd+SvepQHv/pJTao5RKU0qtVUpFmcjpTldq80X7PaaU0kopv57KV5H2KqWeKP05f6+U+pu3M7pTBX6nb1JKrVNK7Sz9ve5lIqc7KKWSlFLZSqndl3lfKaVmlv5dpCml2rnt5Fprv/4CgoFDQGOgCpAKtLxknxeBD0q/fxJYajq3h9t7D1C19PsX/Lm9FW1z6X41gI3AVqC96dwe/hk3A3YCtUpf32A6t4fbOw94ofT7lkCG6dyVaO9dQDtg92Xe7wX8D6CATsC/3XXuQOjhdwDStdaHtdYFwCdA70v26Q0sLP1+GdBNKaW8mNGdrtherfU6rfW50pdbgYZezuhuFfkZA4wH3gHyvBnOAyrS3oHAHK31SQCtdbaXM7pTRdqrgZql318HHPViPrfSWm8ETpSzS29gkS6xFYhQStVzx7kDoeA3AH686HVW6bYy99FaFwE5wPVeSed+FWnvxeIo6S34syu2WSnVFrhRa/25N4N5SEV+xs2B5kqpr5VSW5VSPbyWzv0q0t4xwFNKqSzgC2CId6IZcbX/xissEB5xWFZP/dK5phXZx19UuC1KqaeA9sDdHk3keeW2WSkVBEwHnvFWIA+ryM84hJJhna6UfIL7Sil1q9b6lIezeUJF2hsLfKy1nqaUugNYXNpel+fjeZ3H6lUg9PCzgBsvet2QX3/cu7CPUiqEko+E5X2k8mUVaS9KqfuAN4FHtNb5XsrmKVdqcw3gVmC9UiqDknHPlX584baiv9MrtNaFWusfgP2U/A/AH1WkvXHApwBa6y1AOCULjQWiCv0bvxaBUPC3Ac2UUo2UUlUouSi78pJ9VgIDSr9/DPiXLr064oeu2N7S4Y0PKSn2/jy2e165bdZa52it62it7VprOyXXLR7RWm83E7fSKvI7vZySi/MopepQMsRz2Ksp3aci7T0CdANQSrWgpOA7vJrSe1YC/Utn63QCcrTWx9xxYL8f0tFaFymlBgNfUnK1P0lr/b1SahywXWu9Ekik5CNgOiU9+yfNJa6cCrZ3KlAd+H+l16aPaK0fMRa6kirY5oBRwfZ+CdyvlNoDFAOvaq2Pm0t97SrY3uHAfKXUMEqGN57x106bUiqZkqG4OqXXJN4CQgG01h9Qco2iF5AOnAOeddu5/fTvTAghxFUKhCEdIYQQFSAFXwghLEIKvhBCWIQUfCGEsAgp+EIIYRFS8IUQwiKk4AshhEX8f+oUVWECG9ZVAAAAAElFTkSuQmCC\n", "text/plain": [ "