{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Gaussian process classification\n", "\n", "[Gaussian process regression](gaussian_process_regression.ipynb) introduced how we can use Gaussian processes for regression when we assumed normally distributed data. In that scenario derivation of the posterior GP was especially easy, because it has a closed form solution.\n", "\n", "In this section we will extend the GP framework with classification of data consisting of two groups, where we use most material 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", "We will again use GPy and scipy for demonstration.\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", "from scipy.special import expit\n", "from scipy.stats import bernoulli \n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.rcParams['figure.figsize'] = [15, 6]\n", "\n", "N = 5\n", "base = plt.cm.get_cmap(\"viridis\")\n", "color_list = base(scipy.linspace(0, 1, N))" ] }, { "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": [ "### GP prior\n", "\n", "For classification we will still use a Gaussian process prior, even though our data is not Gaussian:\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", "where $m$ is again the mean function and $k$ is a *kernel function* $k$ is a psd kernel with respective hyperparameters.\n", "\n", "For classification we will a data set $\\mathcal{D} = \\{(\\mathbf{x}_i, y_i)\\}_{i=1}^n$ with $y_i \\in \\{0, 1\\}$.\n", "\n", "Consequently, we need to *squash* a **latent** sample of the GP prior through a logistic function to receive a prior on the probabilities $\\pi(\\mathbf{x})$, i.e. the expeted values of $y_i$;\n", "\n", "\$$\n", "\\pi(x) = \\frac{1}{1+\\exp(-x)}\n", "\$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's visualize some prior samples $\\pi(\\cdot)$. First we generate some data:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "scipy.random.seed(23)\n", "\n", "n = 50\n", "x = scipy.linspace(0, 1, n).reshape((n, 1))\n", "beta = 2\n", "z = expit(scipy.sin(x) * beta)\n", "y = scipy.zeros(n)\n", "y[int(n/2):] = 1" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEcCAYAAAB53pugAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHLFJREFUeJzt3X20HVWZ5/HvkwBq0iGAQGuLJB0EoiPokshbWt4ypiM9Kg3SLyqK6Di8OLgYsrqd1m6II622LaZBGRE72hiUFl2C4yiyFHyDRbdhDWozBGgwAQ3KmyBkICI880fVlcPxntzzUuecfW++n7XO2qR21a5dlZv7Y1ftqhOZiSRJJZs17g5IkjQVw0qSVDzDSpJUPMNKklQ8w0qSVDzDSpJUPMNKklQ8w0qSVLztxt2BmWbFihV55ZVXjrsbkjSdxFQrOLJq2H333TfuLkjSjGNYSZKKZ1hJkopnWEmSimdYSZKKV2xYRcTrIuL8iPhuRPwyIjIi1vbZ1h4RsSYiNkXElojYEBGrI2LnrWzzooj4fETcExGPRcQtEbEqIp7V/1FJkvpR8tT19wAvAR4BfgIs7qeRiNgLuA7YHbgCWA8cCLwTWBERSzPz/rZtDgKuBrYHvgDcBRwF/A2wLCKWZeaWfvojSepdsSMr4AxgH2BH4JQB2rmAKqhOz8xjMvNdmXkU8BFgX+Cc1pUjYjbwKWAO8LrMfH1m/iVwEPBFYGndN0nSiBQbVpl5TWbelgN8lXFELAKWAxuAj7VVnwVsBk6IiLktyw8HXgh8JzO/3NKfJ4G/qP94ckRM+RCbtC245BJYuBBmzarKSy7pb3lTbY26jUHXbXLZqacOZ52p/tz6dzg0mVn8BzgCSGBtj9u9rd7uwg71X6/rl7Use1+97L932OaWun6vyeoPOOCAlLYVa9dmzpmTCU995szJPOWU3pavXdtMW9tvn7nDDqNrY9B1m17W/mlqnak+E3+HA5jy93lk9j1wGZmIOAK4BrgkM9/Yw3YfAlYCKzPzw5PUfxQ4DTg1M/9nvewy4HVUlwC/OMk2XwH+CDg6M7/WXr9kyZJct25dt12UprWFC2Hjxt9ePns2PPFE98sXLKjKJtqazDDbGHTdmWLBAtiwoe/Np7xSVfIEiybMr8uHOtRPLN9pwG2kbdKdd06+vNMv6k7LO7XTT1ujbmPQdWeKrf0dNqHYe1YjMpHmvQwv+9lGmpH23HPy5bNn97Z8zz2ba2vUbQy67kzR6e+vKTM9rCZGQfM71O/Ytl6/20jbpHPOgTlznr5szhx4+9t7W37OOc20tf32sMMOo2tj0HWbXtauqXWmMvF3OFTd3Nga9wcnWEjFWrs2c8GCzIiqnLjR3uvyptoadRuDrtvkslNOGc46U/15wMkVmU6wiL2Af6eaur5XVtPPJ+rmAXdTjS53y8zN9fKjgG9STV0/vK29RcDtwEbg93OSk+cEC0nq2bbxfVYRsX1ELK7D6Tcy83bgKmAh1ay/VquAucDFE0FV+zZwM3BYRLymZR+zgA/Wf/z4ZEElSRqOYkdWEXEMcEz9x+cAfwjcAXy3XnZfZq6s110I/BjYmJkL29ppf93SzVRvozgSuBU4NKd+3dKdwDJgCXAt1WXDSV+35MhKkno2raeuvxR4c9uyRfUHqktxK6dqJDNvj4glwHuBFcDRVJf/zgNWZeYDk2zzLxHxcqrR13JgXr2/9wIf6BRUkqThKHZkNV05spKknm0b96wkSTObYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqnmElSSqeYSVJKp5hJUkqXtFhFRF7RMSaiNgUEVsiYkNErI6Inbvc/oiIyC4+z2/bbmvrXj+co5UkdbLduDvQSUTsBVwH7A5cAawHDgTeCayIiKWZef8UzWwAVnWo2w84FrgpM++apH4j8OlJlv9kys5LkhpVbFgBF1AF1emZef7Ewog4FzgDOAc4eWsNZOYG4OzJ6iLic/V/fqLD5hsyc9JtJUmjVeRlwIhYBCynGhl9rK36LGAzcEJEzO2z/WcDfww8Cnym/55Kkkah1JHVUXV5VWY+2VqRmQ9HxLVUYXYw8M0+2j8ReAZwcWb+osM6O0XEScBzgIeAGzLT+1WSNAalhtW+dXlrh/rbqMJqH/oLq7fV5YVbWeclwD+2LoiIHwAnZOaP+tinJKlPRV4GBObX5UMd6ieW79RrwxFxOLCYamLFdR1WOxdYCuwGzANeDnyBKsCujojn9bpfSVL/Sg2rqURdZh/bvr0uO46qMvPMzLwuM+/LzEcyc11mHg98EdgVWNnHfiVJfSo1rCZGTvM71O/Ytl5XImIX4Dj6n1jx8bo8rI9tJUl9KjWsbqnLfTrU712Xne5pdfJmqokVn8/MB/vo17112dcsRElSf0oNq2vqcnlEPK2PETGP6n7So0Cvs/P+c112erZqKgfX5R19bi9J6kORYZWZtwNXAQuB09qqV1GNbC7OzM0TCyNicUQs7tRmRLwCeCHwb1uZWEFEvGyy57ciYn+qB5EB1nZ5KJKkBpQ6dR3gVKrXLZ0XEcuAm4GDgCOpLv+9u239m+symNzExIqpRlWnA8dGxNXAXcAWqtmDK4DZwEXA5zpvLklqWrFhlZm3R8QS4L1UQXE0cDdwHrAqMx/otq36xbevo7uJFZdTTeDYn+rh5GcC9wNfAy7KzC/3eCiSpAFFZj+zv9XJkiVLct26dePuhiRNJ52uiP1GkfesJElqZVhJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSimdYSZKKZ1hJkopnWEmSild0WEXEHhGxJiI2RcSWiNgQEasjYuce2vhWRORWPs/ssN2LIuLzEXFPRDwWEbdExKqIeFZzRyhJ6sZ23awUEd8D3pGZNw65P6373Au4DtgduAJYDxwIvBNYERFLM/P+Hppc1WH5ryfZ90HA1cD2wBeAu4CjgL8BlkXEsszc0sO+JUkD6CqsgEOB70fERcB7MvOBIfZpwgVUQXV6Zp4/sTAizgXOAM4BTu62scw8u5v1ImI28ClgDvDazPxyvXwW8HnguHr/H+h235KkwXR7GfCCuvwvwK0RcXJExJD6REQsApYDG4CPtVWfBWwGToiIuUPY/eHAC4HvTAQVQGY+CfxF/cehHr8k6em6CqvMfAewBLgW2IUqQG6IiKVD6tdRdXlVHRKtfXm47scc4OBuG4yIP42Id0XEf4uIV0XEM6bY95XtFZl5B3ArsABY1O2+JUmD6XqCRWb+IDMPA94E/Ax4KfCdiLg4Ip7bcL/2rctbO9TfVpf79NDmpcD7gQ8DXwXujIjXjWjfkqQB9DwbMDPXUv1CP5dqcsIbgPURsTIiur0HNpX5dflQh/qJ5Tt10dYVwKuBPYBnAYupQmsn4J8j4lVD3LckqQF9TV3PzEcycyXwEqpZc/OADwI/iojlDfavk4n7RTnVipn5kcz8Smb+NDMfy8xbMvOvgDOpjv9vh7VvSVIzBnrOKjPXZ+YrgeOBO6lGXF+LiC9FxO8P0PTE6GV+h/od29brxyepRoYvjYh5I963JKkHTT0U/GXgJKrZewG8BrgpIs7u9NDtFG6py073hfauy073laaUmY8BD9d/bJ1VOPR9S5J601dYRcSiiPjziPiHiLge+CXwDapZclAF1jOBvwb+LSKO6tBUJ9fU5fL6+abWfc8DlgKPAtf30/+6nX2BnakC676WqqvrcsUk2yyiCrGNwB397luS1Jtu32DxSuCgls+zW6vr8gngR1TTyq+jCpNVwH7AVRHxPzKz01skniYzb4+Iq6ietToNOL+lehXVSOjCzNzc0sfF9bbrW5YtArZk5k/bjmdXqgd/AS7NzNa3WHwbuBk4LCJe0/ZQ8AfrdT6emd6zkqQRiW5+50bEk1QTClofhH0I+BeeCqfrW8Oj3m4W1euRPkAVjG/NzE931bHfft3SzVRBeSTVJbhDW1+3FBEJkJnRsuxEqntT3wZuBx4A9gSOprontQ54ZWY+2Lbv9tct3Qks46lnzTq+bmnJkiW5bt26bg5RklSZ8iULvYTVHTwVTNcCN3U7uoiIPwcuAW7IzJd3s0293fOB91Jdkns2cDdwObCq/ZVPHcJqP6pZfwcAv0c1OeJh4CaqVyddmJm/6rDvF1GN4o6kmu24Efgc8IHMfLRTnw0rSepZY2G1e2beM1BPIh4EZmfmvClXnsYMK0nq2ZRh1e3rlgYKqtqDVK9IkiSpJ029caIbp1LN4pMkqScjC6vM/CrVO/kkSepJ0d8ULEkSGFaSpGnAsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVr+iwiog9ImJNRGyKiC0RsSEiVkfEzl1uPzci3hARn42I9RGxOSIejoh1EXFmROzQYbvcyuf6Zo9SkjSV7cbdgU4iYi/gOmB34ApgPXAg8E5gRUQszcz7p2jmFcBa4AHgGuByYBfg1cDfA8dGxLLMfGySbTcCn55k+U96PxpJ0iCKDSvgAqqgOj0zz59YGBHnAmcA5wAnT9HGz4A3Apdl5q9a2pgHfAs4FDgN+PAk227IzLMH6L8kqSFFXgaMiEXAcmAD8LG26rOAzcAJETF3a+1k5o2ZeUlrUNXLH+apgDqiiT5Lkoan1JHVUXV5VWY+2VqRmQ9HxLVUYXYw8M0+9/F4Xf66Q/1OEXES8BzgIeCGzPR+lSSNQalhtW9d3tqh/jaqsNqH/sPqpLq8skP9S4B/bF0QET8ATsjMH/W5T0lSH4q8DAjMr8uHOtRPLN+pn8Yj4h3ACuBGYM0kq5wLLAV2A+YBLwe+QBVgV0fE8/rZrySpP6WG1VSiLrPnDSOOBVZTTb44LjMfb18nM8/MzOsy877MfCQz12Xm8cAXgV2BlQP0XZLUo1LDamLkNL9D/Y5t63UlIo4BLgXuAY7IzDt67NfH6/KwHreTJA2g1LC6pS736VC/d112uqf1WyLieOAy4OfA4Zl5yxSbTObeutzqLERJUrNKDatr6nJ5RDytj/UzUkuBR4GuZudFxOuBzwGbqILqtj77dXBd9joikyQNoMiwyszbgauAhVQP7bZaRTWyuTgzN08sjIjFEbG4va2IeDPwGeBO4LCpLv1FxMsme34rIvanehAZqrdiSJJGpNSp6wCnUr1u6byIWAbcDBwEHEl1+e/dbevfXJcTky+IiCOpZvvNohqtvSUi2jbjwcxc3fLn06lew3Q1cBewBVhMNXtwNnAR1ShNkjQixYZVZt4eEUuA91IFxdHA3cB5wKrMfKCLZhbw1OjxpA7rbKSaHTjhcqoJHPtTPZz8TOB+4GvARZn55R4PRZI0oMjsefa3tmLJkiW5bt26cXdDkqaT37rk1a7Ie1aSJLUyrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFKzqsImKPiFgTEZsiYktEbIiI1RGxc4/t7FJvt6FuZ1Pd7h7D3rckaXDbjbsDnUTEXsB1wO7AFcB64EDgncCKiFiamfd30c6z63b2Aa4GLgUWA28B/igiDsnMO4axb0lSM0oeWV1AFRanZ+YxmfmuzDwK+AiwL3BOl+38LVVQfSQzl9XtHEMVPLvX+xnWviVJDSgyrCJiEbAc2AB8rK36LGAzcEJEzJ2inbnACfX6Z7VVf7Ru/w/r/TW6735dcgksXAizZlXlJZcMXjesdmdKf0ZZV/q+Bvk7lYYqM4v7AG8DEriwQ/3X6/plU7TzH+v1vt6h/sK6/q1N7fuAAw7Ifq1dmzlnTiY89Zkzp1reb92w2p0p/Rll3SmnlL2vfvsnNWDqXOhmpVF/gA/VgXBmh/qP1vWnTNHOafV653eoX1nXf7CpfQ8SVgsWPP2XwcRnwYL+64bV7kzpzyjrZs8ue1/99k9qAFN9Sp1gMb8uH+pQP7F8pyG009S+e3bnnb0t77ZuWO3OlP6Mqu6JJ8reV9P9k5pU5D2rLkRd5hjaaWrfv2XPPTsv77duWO3OlP6Msm727LL31W//pJHoZvg16g/b6GXAbeEeUWn98Z7V4P2TGjB1LnSz0qg/bKMTLDKrf/wLFmRGVGXrL4N+64bV7kzpzyjrSt/XIH+n0gCmzIXIbPxq1sDqh3L/nWr6+F6Z+WRL3TzgbqpLmLtl5uattPM7wD3Ak8BzM/PhlrpZwO3AwnofdzSx7yVLluS6dev6Om5J2kbFVCsUec8qM28HrqIKktPaqlcBc4GLW8MiIhZHxOK2dh4BPlOvf3ZbO++o2/96trzBop99S5KGq8iRFUz6yqObgYOAI4FbgUOz5ZVHEZEAmRlt7bS/bulfgRcCr6UadR1aB1Tf+27lyEqSejY9R1bwmxHOEuDTVEFxJrAXcB5wSKewmKSd+4FD6u1eULdzEPAp4ID2oGpy35KkZhQ7spquHFlJUs+m78hKkqQJhpUkqXiGlSSpeN6zalhE3AtsHHc/JGkauS8zV2xtBcNKklQ8LwNKkopnWEmSimdYSZKKZ1iNQETsERFrImJTRGyJiA0RsToidu6xnV3q7TbU7Wyq291jWH0ftkHPTUTMjYg3RMRnI2J9RGyOiIcjYl1EnBkROwz7GIalqZ+btjYPi4gnIiIj4n1N9neUmjw3EbFfRFwcEXfVbd0TEd+OiDcNo+/D1uDvmz+IiCvq7R+LiDsj4qsRsdWJEMPiBIshm+Q9g+uBA6neM3gLsLSb1zdN8o7D7wOLeeodh4e0vpB3Omji3NT/cL4GPABcQ/XG/F2AVwPPqdtflpmPDekwhqKpn5u2NucBPwR2BX4HOCcz39Nkv0ehyXMTEScCnwT+H/AVqm9b2Al4MbApM/+s4e4PVYO/b04BLgA2A18CfgLsARwLzAHek5nnDOMYOurme0T8DPTdXBPff/Vf25afWy//eJftTHz31rlty0+vl1857mMdx7kBXgq8Adihbfk84Aa28kWaJX+a+rlp23YNVaj/Vd3G+8Z9nOM8N8DBwK+BG4HnTFK//biPdRznBtgeeBB4FNi3re6FwGNU4f6MkR7buE/uTP4Ai+ofkB8Ds9rq5gGPUP2fy9wp2plb/3A8Asxrq5tVt5/AonEf86jPzRT7eH29j/817uMd97mhGoEn8EbgxOkaVk2eG+A7dVsvHvdxlXRugN+t2/lBh/of1vXPHuXxec9quI6qy6uy5UscAbL6IshrqYbUB0/RziHAs4Brs+ULJOt2nqT6/i2ohvrTRVPnZmser8tfD9DGODR6biJid+Ai4PLMXNtkR8egkXNT3+d9BbAOuCkijoyIlfV9zmX1l7NON0393NwD3AvsExF7t1ZExD7A3sCNOeJvn5iOfyHTyb51eWuH+tvqcp8RtVOSURzTSXV55QBtjEPT5+YTVP/WTx6kU4Vo6ty8vGX9q+vPh4C/B74B3BgRLxign+PQyLnJavh0GtXPzA0R8U8R8f6IuJjq0vpNwPEN9Lcn2416h9uY+XX5UIf6ieU7jaidkgz1mCLiHcAKqvsRa/ppY4waOzcRcRLVJcA/zcyfN9C3cWvq3Oxel38C3Ec1ceCbwG7AWcAJwP+OiP0y81f9d3ekGvu5yczLImIT8DmgdVbkz6m+C3Dkk7kcWY3XxHe4DDols6l2StL3MUXEscBq4GfAcZn5+BSbTDddnZuIWEh1Hi7LzM8PuU+l6PbnZnZL+bbM/FJm/jKrL159M9XlwX2A44bTzbHo+t9URLyRaoT5XapJFXPq8pvAR4FLh9THjgyr4Zr4P5n5Hep3bFtv2O2UZCjHFBHHUP1Dugc4IqfZdP5aU+dmDdWMrlOb6FQhmjo3v6jLLcBXWyvqy2BX1H88sNcOjlEj56a+L7WG6nLfCZm5PjMfzcz1VCPOG4DjI+KIwbvcPcNquG6py07XiCduXna6xtx0OyVp/Jgi4njgMqpLFYdn5i1TbFKqps7Ny6gud91bPwScEZFUl3EA3l0vu3yw7o5U0/+mHm6fjFCbCLNn9dC3cWvq3Cynmr7+7UkmajxJNYsS4IB+Otkv71kN1zV1uTwiZrX+xdcPaC6l+j/f66do5/p6vaURMa91RmA9a2l52/6mg6bOzcQ2rwcuBn4KHDlNR1QTmjo3F1Ndvmm3N3AY1f28G4D/M3CPR6epc/NDqntVu0bE705yP+/Fdblh8C6PTFPn5hl1uVuH+onlo72XN+5nA2b6hx4f0qN6K8XiSdqZeCj4w23Lt5mHgrdybt4MPEF103fBuI+rpHPToe0TmabPWTX8c/O+ev1/ouW5JGA/ql/qjwMvGPfxjvrcUF36TKpnO/dvq3tpfW6eBP7DKI/N1y0N2SSvP7kZOIjqmahbgUOz5XmF+jINmRlt7bS/bulfqW54Trxu6dCsbg5PG02cm4g4kupG8Cyq6+x3TbKrBzNz9ZAOYyia+rnp0PaJVJcCZ8rrlvr9NzWHasLAwVSjy29RjRqOo7r8d2Zmnjvkw2lUg+dmDfAWqtHTl6i+UHYhcAywA7A6M88Y8uE83bj/T2Bb+ADPp/rlcHf9l78R+Adgl0nWTep7vJPU7VJvt7Fu526qX9B7jPsYx3VueGqUsLXPhnEf5zh/biZZd+KcTcuRVZPnhuoy6dlU79DbQjX54BvAq8Z9jOM8N1QzB0+kCvBfUD1Y/wBVuP/ZOI7LkZUkqXjOBpQkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJIkFc+wkiQVz7CSJBXPsJJmmIhY2vLNwMd3WOegiHikXufvRt1HqVe+dV2agSLiCuA1VF998eLMfKKlbl/ge8CuVF88+Jb0F4EK58hKmpneRfXtyYuBN04sjIjfo/o22V2BrwBvM6g0HTiykmaoiPgk8Fbgx8C+wFzgO1Rf2/49YHlmPjq+HkrdM6ykGSoingfcRvUV7WcAfwwcBvwIOCwzHxxj96SeGFbSDBYRHwD+smXRBmBpZm4aT4+k/hhW0gwWEc8FfkJ1f/oB4ODMvG28vZJ65wQLaYaKiO2AT/DUv/M5gPeoNC0ZVtIMFBEBfBL4T8C9VJMsngmsGme/pH55GVCagSLiQ8BK4BHgKOAFwGepprPvn5n/d4zdk3rmyEqaYSJiJVVQPQ4cl5nfBy4FfgjMBt4/xu5JfTGspBkkIt4E/B2QwImZeRVA/eDvX9ervSYilo6pi1JfvAwozRARcTRwBbAdcEZmrp5kneuBg4BrM/MPRtxFqW+OrKQZICIOAS6jCqoPThZUtXfX5dKIeO1IOic1wJGVJKl4jqwkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxTOsJEnFM6wkScUzrCRJxfv/rPESI0aWNGAAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "_, ax = plt.subplots()\n", "ax.scatter(scipy.sin(x), y, color=\"blue\")\n", "plt.xlabel(\"$x$\", fontsize=25)\n", "plt.ylabel(\"$y$\", fontsize=25)\n", "ax.tick_params(axis='both', which='major', labelsize=20)\n", "ax.tick_params(length=0, color=\"black\")\n", "ax.spines[\"top\"].set_visible(False)\n", "ax.spines[\"right\"].set_visible(False)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we define the kernel and likelihood:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "m = scipy.zeros(n)\n", "kernel = GPy.kern.RBF(input_dim=1)\n", "\n", "k = kernel.K(x, x)\n", "lik = GPy.likelihoods.Bernoulli()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we sample some prior values. First we generate a sample from the latent GP. Then we transform these values using the sigmoid function. In the end we binarize the outcome of the last step." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "f_prior = mvrnorm(mean=m, cov=k)\n", "p_prior = lik.gp_link.transf(f_prior)\n", "y_prior = lik.samples(f_prior).reshape(-1,1) " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEACAYAAAByG0uxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X+QVOWd7/H3d9TIZGaCaLC4FaI4XBWIPyiYKP6IIHhdIug1ImWVa9RN4upKImSlsq5ubcS6JnETzfgriyaLqIkxRm7FbCoYIfgres1eYpEYRVQIKFGugEpgBKLyvX+c09I0fXpOP31O9/TM51XVdWbOec73ec7TZ/rbz/k15u6IiIiEaGl0A0REpHkpiYiISDAlERERCaYkIiIiwZREREQkmJKIiIgEUxIREZFgSiIiIhJs30Y3IG9Tp071hx9+uNHNEBFpJpa2YL8fiWzatKnRTRAR6bf6fRIREZH8KImIiEiw1EnEzM41s1vN7Ekz+4uZuZn9MKRSMxtuZgvM7HUz22lma82s28yGVFhnjJk9YGZvmtkOM1tlZvPMrDWkDSIiUrtqTqz/C3AssA1YD4wKqdDMRgJPAwcDDwEvAscBs4GpZnaSu28uWed4YBmwH/Ag8BowGfhXYIqZTXH3nSHtERGRcNUczvoqcATwMeAfaqjze0QJ5Ap3P9vdr3L3ycB3gSOB64sLm9k+wF3AR4Fz3f18d/8n4HhgEXBS3LZclD4pv7cn51dbPu84WcdqRPysVGpns2xDkmZvf4iBuM1JGtkXqZOIuz/q7i97Df+AxMw6gdOBtcDtJYu/DvQAnzeztqL5E4HRwBPu/vOi9uwCvhb/epmZpb4kLa2FC+H223e/Ie7R7wsXZlM+q3rrFasR8bNSqZ3Nsg1Jmr39IQbiNidpdF/U+8T65Hj6SJwEPuTuW4GniEYcE8qss9fNHu6+BngJOBTozLKh7rBtGyxatPsNuv326Pdt28pn/mrKZ1VvvWI1In5WKrVz69bo1de3IUmzvAdZGojbnKQv9EW9bzY8Mp6+lLD8ZaKRyhHAr6tY54j4tTqDNgJgBrNmRT8vWhS9AGbMiOaXjnuqLZ9VvfWK1Yj4WemtnYUyfXkbkjTLe5ClgbjNSfpCX9R7JDI4nm5JWF6Yf0CN62Si+A0qqPTGVFs+q3rrFasR8bNSqZ3Nsg1Jmr39IQbiNidpdF/0tftECptdzSAsZJ1UCkPDYsXHHmstn1W99YrViPhZqdTOZtmGJM3e/hADcZuTNLov6p1ECqOGwQnLP1ZSLnSdmhUfW5wxA5Yti6bFxx5rKZ9VvfWK1Yj4WanUzttui159fRuSNMt7kKWBuM1J+kJf1PucyKp4ekTC8sPjafH5j5B1amYG7e17HlssDBnb28ufE6mmfFb11itWI+Jnpbd2Qt/fhiTN8h5kaSBuc5K+0BcWcsWumU0CHgV+5O4XVLHeSOAVokt8RxZfoWVmHcAbRKOjoe7eE8+fTHSS/Ql3n1gSr5PoZPo64LBylx93dXX58uXLq9q+Yu57vhGlv9daPu84WcdqRPysVGpns2xDkmZvf4iBuM1JcuiLxj7F18z2M7NRcdL4kLuvBh4BRgAlp4KYB7QB9xQSSOxxYCVwipmdVVRHC3BD/Ov8Wu5fqaTciCPL8nnHyTpWI+JnpVI7m2UbkjR7+0MMxG1O0si+SD0SMbOzgbPjX4cBfwOsAZ6M521y97lx2RHAn4B17j6iJE7pY09WEt19firRIakTUzz25FVgCtBFdG9J4mNPah2JiIgMQKnTUDXnRMYCF5XM62T3TX7rgLm9BXH31WbWBVwHTAXOIDqMdQswz93fKrPOb83s00SjldOBjri+64Bv6blZIiKNEXROpJloJCIiUjX9Z0MREcmfkoiIiARTEhERkWBKIiIiEkxJREREgimJiIhIMCUREREJpiQiIiLBlERERCSYkoiIiARTEhERkWBKIiIiEkxJREREgimJiIhIMCUREREJpiQiIiLBlERERCSYkoiIiARTEhERkWBKIiIiEkxJREREgimJiIhIMCUREREJpiQiIiLBlERERCSYkoiIiARTEhERkWBKIiIiEkxJREREgimJiIhIMCUREREJpiQiIiLBlERERCSYkoiIiASrKomY2XAzW2Bmr5vZTjNba2bdZjYk5fqTzMxTvD5Zsl6lss9Usw0iIpKdfdMWNLORwNPAwcBDwIvAccBsYKqZneTum3sJsxaYl7DsaOAc4Hl3f63M8nXAwjLz1/faeBERyUXqJAJ8jyiBXOHutxZmmtlNwFeB64HLKgVw97XAteWWmdmP4x/vTFh9rbuXXVdERBoj1eEsM+sETicaSdxesvjrQA/weTNrC2mEmR0EfA7YDtwbEkNEROov7Uhkcjx9xN13FS9w961m9hRRkpkA/DqgHRcD+wP3uPvbCWUOMLMvAMOALcDv3F3nQ0REGihtEjkynr6UsPxloiRyBGFJ5Evx9I4KZY4F/qN4hpn9Hvi8uz8XUKeIiNQo7dVZg+PploTlhfkHVNsAM5sIjCI6of50QrGbgJOAoUAH8GngQaLEsszMPlFtvSIiUrus7hOxeOoB6/59PE0chbj7le7+tLtvcvdt7r7c3WcCi4CPA3MD6hURkRqlTSKFkcbghOUfKymXipkdCMwg/IT6/Hh6SsC6IiJSo7RJZFU8PSJh+eHxNOmcSZKLiE6oP+Du71S5LsDGeBp0VZiIiNQmbRJ5NJ6ebmZ7rGNmHUTnK7YD1V4tdUk8Tbo3pDcT4umawPVFRKQGqZKIu68GHgFGALNKFs8jGgnc4+49hZlmNsrMRiXFNLPPAKOBP1Y4oY6ZjSt3/4mZHUN0gyPAD9Nsh4iIZKuaO9YvJ3rsyS1mNgVYCRwPnEp0GOuakvIr46lRXuGEem+jkCuAc8xsGfAasJPoaq6pwD7A94EfJ68uIiJ5SZ1E3H21mXUB1xF9gJ8BvAHcAsxz97fSxoof2Hgu6U6o/4zoxP0xRDc9DgI2A4uB77v7z9PWKyIi2TL3kKtym0dXV5cvX7680c0QEWkmSUeQ9qL/JyIiIsGUREREJJiSiIiIBFMSERGRYEoiIiISTElERESCKYmIiEgwJREREQmmJCIiIsGUREREJJiSiIiIBFMSERGRYEoiIiISTElERESCKYmIiEgwJREREQmmJCIiIsGUREREJJiSiIiIBFMSERGRYEoiIiISTElERESCKYmIiEgwJREREQm2b6MbICJST++99x7r169nx44djW5Kww0aNIjhw4ez3377BcdQEhGRAWX9+vV0dHQwYsQIzKzRzWkYd2fz5s2sX7+eww47LDiODmeJyICyY8cODjrooAGdQADMjIMOOqjmEZmSiIgMOAM9gRRk0Q9KIiIidbbPPvswduxYjjrqKGbOnMm7774bHGv+/Pncc889GbauOkoiIiJ11trayooVK/jjH//IRz7yEebPnx8U5/333+eyyy7jwgsvrGqdLCmJiIhUsHTpUqZPn86YMWOYPn06S5cuzTT+Zz7zGV555ZW95re3t3PllVcybtw4pkyZwsaNGwGYNGkSV199NRMnTuTmm2/m2muv5Tvf+Q4AK1asYMKECRxzzDF87nOf4+233y67TpaUREREEixdupTZs2ezYcMGhg4dyoYNG5g9e3ZmieT9999n8eLFHH300Xst6+npYdy4cTz77LNMnDiRefPmfbjsnXfe4fHHH+fKK6/cY50LL7yQG264gT/84Q8cffTRqdaplZKIiEiC7u5uWltb6ejooKWlhY6ODlpbW+nu7q4p7vbt2xk7dixdXV0ccsghfPGLX9yrTEtLC+eddx4AF1xwAb/5zW8+XFaYX2zLli288847TJw4EYCLLrqIJ554ouI6WdB9IiIiCdasWcPQoUP3mNfW1saaNWtqils4J1KN4iup2traqq4zZJ00qhqJmNlwM1tgZq+b2U4zW2tm3WY2pIoYj5mZV3gNSlhvjJk9YGZvmtkOM1tlZvPMrLWabRARSauzs5Oenp495vX09NDZ2Zl73bt27eLBBx8E4L777uPkk0+uWH7w4MEMGTKEJ598EoB77733w1FJnlKPRMxsJPA0cDDwEPAicBwwG5hqZie5++Yq6p6XMH+vSwfM7HhgGbAf8CDwGjAZ+FdgiplNcfedVdQtItKrOXPmMHv2bCD6Jt/T08P27duZM2dO7nW3tbXx/PPPM378eAYPHsxPfvKTXte5++67ueyyy3j33Xfp7Ozkrrvuyr2duHuqF/ArwIGvlMy/KZ4/P2Wcx6JqU9e7D/BCXMdZRfNbiBKKA1clrT9+/HgXESl44YUXqiq/ZMkSnzZtmo8ePdqnTZvmS5Ysyalle2pra6tLPQn9kTo3mEcfyBWZWSewGlgLjHT3XUXLOoA3AAMOdveeskF2l38MmOjuqW6VNLPJwK+BJ9x9YsmyQrvWAYd5mY3p6ury5cuXp6lKRAaAlStXMnr06EY3o1ft7e1s27Yt93oS+iP1rexpz4lMjqePFCcQAHffCjwFfBSYkLZiMzvPzK4ys380s8+a2f691P1w6QJ3XwO8BBwK5H+QUkSkTuqRQLKQNokcGU9fSlj+cjw9ooq67we+CdwI/BJ41czOrVPdIiKSgbRJZHA83ZKwvDD/gBSxHgLOBIYDrcAoomRyAPATM/tsjnWLiEiGsrpPpHD8rNcTLO7+3ZJZq4Crzex14FbgG8DiPOoWEZFspR2JFL7tD05Y/rGSciF+QHR579j4ZH096xYRkQBpk8iqeJp03uHweJp03qJX7r4D2Br/WnxrZe51i4jUi7tz8skns3jx7gMuDzzwAFOnTg2K1+hHwac9nPVoPD3dzFrKXOJ7ErAdeCa0IWZ2JDCEKJFsKlq0DLgGmEp07qR4nU6i5LIOqO05BCIiZbhD8f9uKv29WmbG/PnzmTlzJqeeeioffPAB11xzDQ8/vNcFqL0qPAq+2nX23Te7J16liuTuq83sEeB0YBbRuYuCeUQjhzuK7xExs1Hxui8WzesEdrr7n4vjm9nHgcKtlfe7e/Fd648DK4FTzOwsd/95vE4LcENcZn65e0RERGqxcCFs2wazZkWJwx1uvx3a2+Hii8PjHnXUUZx55pnccMMN9PT0cOGFFzJy5Mg9yrS3t3PppZfy6KOPMmTIEO6//36GDh3KpEmTOPHEE3nqqac466yz2Lp1K+3t7cydO5cVK1Z8eMf6yJEjWbBgAUOGDNlrnSyf5FvNs7MuB94EbjGzn5nZN81sGfBVokNJ15SUXxm/ip0CrDOzX5vZnWb2LTO7D3gFOAFYDnyteAV3/wD4O+Bd4EEzu8/MvgX8FjiX6B6V0pP1IiI1cY8SyKJFUeIoJJBFi6L5tX5t/frXv859993H4sWL+drXvrbX8mZ5FHzqMU08GukCriM6tHQG0Z3qtwDz3P2tFGF+B/wQGA+MJTopvhV4DniAaDTz1zJ1/9bMPk006jkd6CA6hHUd8C3Xc7NEJGNm0QgEosSxaFH084wZu0cmtWhra+O8886jvb2d/fff+17r0kfBn3POOR8uS/so+JkzZ1ZcJwtVHRhz99eIRgVpyu7Vxe7+HHBxNXUWrfsCMLPXgiIiGSkkkkICgWwSSEFLSwstLekOCPWLR8GLiAwkhUNYxQqHtvLW7x4FLyIykBSfAykcwir8DtmOSMpplkfBp3qKbzPTU3xFpFg1T/HN6+qsNJrlKb4aiYiIJLj44j3vCymcI8lzBNJsdE5ERKSC0oRRrwTS3x4FLyIishclEREZcPr7ueC0sugHJRERGVAGDRrE5s2bB3wicXc2b97MoEGDaoqjE+siMqAMHz6c9evXs3HjxkY3peEGDRrE8OHDa4qhJCIiA8p+++3HYYcd1uhm9Bs6nCUiIsGUREREJJiSiIiIBFMSERGRYEoiIiISTElERESCKYmIiEgwJREREQmmJCIiIsGUREREJJiSiIiIBFMSERGRYEoiIiISTElERESCKYmIiEgwJRERkSa3dOlSpk+fzpgxY5g+fTpLly6tW91KIiIiTaJcsli6dCmzZ89mw4YNDB06lA0bNjB79uy6JRLr7/9nuKury5cvX97oZoiI1KSQLFpbW2lra6Onp4ft27fT0dHB+++/T0dHx4dlt27dyrBhw/jFL34RWp2lLaiRiIhIH1NuxNHd3U1raysdHR20tLTQ0dFBa2srzz//PG1tbXus39bWxpo1a+rSViWRXjTyWKOI9G/VHJ567rnnyiYLgJ6enj3m9/T00NnZWZdt0OGsCpKGjzfffDMA3d3drFmzhs7OTubMmcNpp5324TeG0vkiMnCV+1wAqjo89eqrr3LIIYfsNX+fffZh27ZtZT+navjsSX04S0mkgunTp7Nhw4bUb9pFF13E3XffXVXSEZH+I4tksWrVKsaNG0dLy+4DRbt27eJPf/rTh4ew6vD5oiRSUEsSGTNmDEOHDt3rzXz22Wc58sgjM/mmABrRiDSbPJNF0ufLsGHDmDNnTr0+F5RECvIYiSS9+c888wwTJkxIvVNoRCPStzUiWeR0eKpaSiIFeZwTqfaYZbU7UZYjmsJ2KPGIlJf091Ht33+WyQIa/mUxnyRiZsOB64CpwEHAG8DPgHnu/naK9duAs4FpwDjgk8AuYBXwY+BWd/9rmfUqNfK37j4haWGt94lU800kaQRR7U6X1YimsDPq4gAZSColhbR/yzfffDPd3d1VHYlo0mSRJPskYmYjgaeBg4GHgBeB44BTiZLASe6+uZcYU4HFwFvAo8ArwIHAmcCwOP4Ud99Rsp4D64CFZcKud/cfJNWZ182GWeyoeY9ohg0bBqCLA6RfyvML3rBhw1izZk1V50SbNFkkySWJ/Ao4HbjC3W8tmn8T8FXgDne/rJcYY4FPAT8tHnGYWQfwGNHoZK6731iyngOPu/ukVI0t0lfuWG/EiGbjxo0Afe7iAI10pJy+9MVs48aNdHZ2VvUFrEmTRZJsk4iZdQKrgbXASHffVbSsg+iwlgEHu3tP2SC913E+8CPgF+5+Zsmypk8iSfL8w6k0EmnUxQEhIx0lneaUxb7dqEPEhSuhqj0U3I9knkS+BHwfuNPdLy2zvDBKOc3df11FQ4tjzAQeAH7m7p8rWebA74FbiA57bQF+5+7P9Ba3ryeRalV7XDdpWaMuDqh2pJNl0qmUjJSodqu2//JMCo28WGWAf4HJPIl8G5hLmUNN8fLbgFnA5e7+71U0tDjGYqIT9pe5+x0ly5Ia+Xvg8+7+XFLc/pZEklT7AQmN+eZX7TfCrJJOb8koq4sP+sP8avaLvJOCLptvmMyTyJ3AJcAl5U5im9n1wNXA1e7+zSoaWlj/y8CtwArgOHd/r2T5jcAi4CVgBzAK+CfgXGATMNbd/1wu9kBJIiEacQw67w+RauvN8uKDZp9f6YqkvvZ+6gbe3NU9iXwD+Gfgn939W1U0FDM7h+gw1kaiK7xSP3rSzB4EZgDd7v7VcmWURLKTxYimUd9ck+ZnefFBs8+vdEVSXxtZagSRu9RJZN+U5bbE08EJyz9WUi4VMzsbuB94Ezi1mgQSm0+URE6pcj0JcNppp5X9Iy18gy39Yx43blzq+RAlI2CPD4u5c+dy99137zX/U5/6FD09PXt8GPX09DB48OCq5heedFr67bvwVNRyT03dsmVLv5xfeD/K9UW1/TpmzBi2bdv2Yeze3s9KSSFpPwKUNPqAhp1Yj0+k3wdsACa7+8vVNDyOcSzRIbBV7j6qXBmNRJpHniduQ86JZHWorlnmV7oiSVfbDTiZH84aSXRj4FqSL/FtAYamucQ3vpz3HuDPhI1ACnEuJRqNLHb3M8qVURLpv7I6mZwUC7JJVM0yv7crknTfz4CSOong7qlewK8AB75SMv+meP78kvmjgFFl4lwEfACsAQ5NUe84oK3M/GOITqo7cH7S+uPHj3eRUEuWLPFp06b56NGjfdq0ab5kyZJ+PV8kljo31PLYk5XA8USPPXkJONGLHntSuCzX3a1o3qnAUqJRywLgtTJVvePu3UXrLATOAZbF5XfGCWoqsA/RYbZLPWFDNBIREala5ifWcffVZtbF7gcwnkF0GOsWogcwvpUizKHs/pe8X0gosw7oLvr9Z0Qn7o8BJgODgM1Ez+D6vrv/PO02iIhItvQoeBERKZV6JNLSexEREZHylERERCSYkoiIiARTEhERkWBKIiIiEkxJREREgimJiIhIMCUREREJpiQiIiLBlERERCSYkoiIiARTEhERkWBKIiIiEkxJREREgimJiIhIMCUREREJpiQiIiLBlERERCSYkoiIiARTEhERkWBKIiIiEkxJREREgimJiIhIMCUREREJpiQiIiLBlERERCSYkoiIiARTEhERkWBKIiIiEkxJREREgimJiIhIMCUREREJpiQiIiLBlERERCRYVUnEzIab2QIze93MdprZWjPrNrMhVcY5MF5vbRzn9Tju8LzrFhGR7OybtqCZjQSeBg4GHgJeBI4DZgNTzewkd9+cIs5BcZwjgGXA/cAo4O+AaWZ2gruvyaNuERHJVjUjke8RfYhf4e5nu/tV7j4Z+C5wJHB9yjjfIEog33X3KXGcs4kSwsFxPXnVnRn3yr+nLd+oOCGx8m5r3vND9LVtqMf+0te2rVH7XZax+uK+nRXzFK0ws05gNbAWGOnuu4qWdQBvAAYc7O49FeK0ARuBXcB/c/etRcta4jpGxHWsyaLurq4uX758ea/bWI2FC2HbNpg1C8yiN/L226G9HS6+OH35FStg7Nj6xwmJlXdb855fqS+qfZ8btQ312F9gYG1zf+6LGlnagmlHIpPj6SPFH+IAcSJ4CvgoMKGXOCcArcBTxQkkjrMLeCT+9dQc6s6Ee7RjLVoUvYGFN3LRomh+uW8K5co/+CBs2VL/OCGx8m5r3vMr9UW173OjtqEe+8vWrdGrr2xbo/a7/tAX9ZR2JPJtYC4w191vLLP8NmAWcLm7/3uFOLOA24Db3P0rZZbPBb4N/Ju7/1MWdecxEil+AwtmzNj9DSFt+csvh+99r/5xQmLl3da851fqiyR9bRvqsb/AwNrm/twXNUodMW0SuRO4BLjE3X9QZvn1wNXA1e7+zQpxriY6f3G9u/9LmeWXAHcCd7r7pVnUnUcSgWhnnTx59+/LllV+I5PKNypOSKy825r3/BB9bRvqsb/0tW1r1H7XH/qiBpkfzkpbYa2DqpA4WdWdWuHbTrHCELOa8rt2NSZOSKy825r3/JDhfl/bhnrsLwNtm/tzX9SNu/f6IjrE5MCVCctvi5f/Qy9xZsXlbk1YPjdefkNWdY8fP96ztGuX+623uk+aFE3L/Z6m/MSJ7ueeW/84IbHybmve8yv1RbXvc6O2oR77yy23RK++sm2N2u/6Q19kIFVucPfU94msiqdHJCw/PJ6+lEOcrOrOhFl0NUTxscjC8dP29r2HlpXKr1gBEyfWN05orLzbmvf8pL5I0he3oR77Cwysbe6vfVFPac+JjAReofJlti3AUK98iW878CbVXeJbU915nhMpfuNKf09bvlFxQmLl3da854foa9tQj/2lr21bo/a7/tAXNUgdNdU5EXdfTXT57QiiQ1LF5gFtwD3FH+JmNsrMRpXE2QbcG5e/tiTOl+P4vyokkNC666Hct6CQ8o2KExIr77bmPT9EX9uGeuwvfW3bGrXfZRmrL+7bWUk1EgEo8+iRlcDxRPd0vASc6EWPHjEzB3B3K4lT+tiT/wJGA/+TaJRyYpw4gusultdIRESkH8t2JAIfjgi6gIVEH+BXAiOBW4ATkj7Ey8TZTHTT4S3Af4/jHA/cBYwvTSBZ1i0iItlKPRJpVhqJiIhULfuRiIiISKl+PxIxs43Auka3Q0SkiWxy96lpCvb7JCIiIvnR4SwREQmmJCIiIsGUREREJFi/TiJmNtzMFpjZ62a208zWmlm3mQ2pMs6B8Xpr4zivx3GH5113PdTaVjNrM7O/NbP7zOxFM+sxs61mttzMrjSzjySs5xVez2S7ldnI4n01s8d62fZBCeuNMbMHzOxNM9thZqvMbJ6ZtWa3hdnJYL+a1Es/FV6fLFmvafYrMzvXzG41syfN7C9xG38YGKvq/s5in+q3J9bL3OX+InAc0V3uq4CT0tykWOYO+/8LjGL3HfYnFD+mJcu66yGLtprZVGAx8BbwKNGzzg4EzgSGxfGnuPuOkvWc6Mq5hWXCrvcy/z+mkTLcpx4DJhI9tqec/+Xu75esczzR/rcf8CDwGtF//ewi+u+eU9x9Z/VblY+M9qsRwMUJi48GzgGed/ejStZrmv3KzFYAxwLbgPVEny0/cvcLqoxTdX9ntk9V88jfZnoBvyJ6RPxXSubfFM+fnzLOHXH5m0rmXxHPfzivupuln4CxwN8CHymZ3wH8joRH+cfzH2t0HzRgn3os+tNLXe8+wAtxHWcVzW+J//gduKrR/ZNHX1WI/+M4zhXNvF8RfcgfTnRz36S47T/Mu7+z3Kca3ok5vTGdcSf8CWgpWdZBlPV7gLZe4rQB78blO0qWtcTxHejMuu5m6qde6jg/ruM/yyxrpj/2zPoqIIlMjut+vEK71hIfWWj0K+/9CjgI2BH/bQ5p5v2qpN1BSSSkv7Pcp/rrOZHCP5B8xIseHQ/g0ePnnwI+CkzoJc4JQCvwlBc9tj6Os4vo6cIQfZvIuu56qEdb34un7ycsP8DMvmBmV5vZLDPrC/1STuZ9ZWbnmdlVZvaPZvZZM9u/l7ofLl3g0aHUl4BDif74+4K896uLgf2Bn7r72wllmmW/ykJIf2e2T/XXJHJkPE36R1Uvx9Okf3RVS5ys6q6HerT1C/F0r501dizwH8D1RP+l8v+Y2QozO7qGOvOQR1/dD3wTuBH4JfCqmZ1bp7rzlHd7vxRP76hQpln2qyw09HOqvyaRwfF0S8LywvwDcoiTVd31kGtbzezLwFRgBbCgTJGbgJOAoUTD7k8THY89FlhmZp8IqTcnWfbVQ0QXHQwnGumOIkomBwA/MbPP5lh3PeTWXjObSNRfz7v70wnFmmm/ykJDP6f6axLpTeEJlbVemhYSJ6u66yG4rWaEbeizAAADVklEQVR2DtANbABmuPt7pWXc/Up3f9rdN7n7Nndf7u4zgUXAx4G5NbS93lL3lbt/191/4e5/dvcd7r7K3a8m+hcHLcA38qq7j6ilvX8fTxNHIf1sv8pCrp9T/TWJFLLo4ITlHyspl2WcrOquh1zaamZnEx2qeROY5CWXQKcwP56eUuV6earH+/oDonNHYy3618/1rDtLee1XBwIzgO1E/yG1Wn1xv8pCQz+n+msSWRVPk47nHR5Pk44H1hInq7rrIfO2mtlM4KfA/wMmuvuqXlYpZ2M8bQtYNy+5v68e3UdTuICjeNubaZ+C/Np7EdEJ9Qfc/Z2AdvXF/SoLjf2cavRlbTldKjeS3i95e5feL/Ftp/pLfDOpu5n6qWid84m+Sa8r7pOAdl0at+uXje6jer6vRCc7HfgLsG/R/Ga7xDeXvmL3fQ0n9pf9qqR9kwi7xLfq/s5yn2p4x+X4hlR7880oYFSZOIWbDW8smT8gbzas0E8XAR8Aa4BDU9Q7rtyHCHAMsCmu+/xG90/WfRX/gX6iTOyPE91x7MCdJcsq3Rj2U/rBzYZJ+1XR8s/E6z3X3/arojZWTCJEd5aPAkZm0N+Z7VMD6bEnK4n+P/upREO0E73oMQDxoxJwdyuJU/rYk/8CRrP7sScnesn/ha+27kbKop/M7FRgKdEOuIDo8Qml3nH37qJ1FhI9tmJZXH4n0R/IVKId/PvApd6HdtCM+upionMfjwOriR4VcwhwBtHx6eXA//CSwzVlHlHxKjCF5nnsSdDfX9Hye4ELiO5Qv7VCvQtpov0qPn94dvzrMOBviL6IPRnP2+Tuc+OyI4hGG+vcfURJnKo/czLbpxqdfXPO7J8E7gLeAP5KdJjlZuDAMmWdhLuIiZ4DdXO8/l/jeAuA4VnU3ehXrf1EdPOX9/JaW7LO2cD/JnrO1l+K+vU/Kfpm1NdeGfTV0UTPdHoO2Ex0M+ZbRB8aX6Hk0TEl644h+pa4iejD8SWi52+1Nrpf8uiromVDiE6mvwsc0EudTbVfAdem/bsBRpT7Wwrp7yz3qX47EhERkfz116uzRESkDpREREQkmJKIiIgEUxIREZFgSiIiIhJMSURERIIpiYiISDAlERERCaYkIiIiwZREREQk2P8H0YbQSH4TEr4AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "_, ax = plt.subplots()\n", "ax.scatter(x, p_prior, color=\"k\", marker=\"o\", alpha=0.75, label=\"P prior\")\n", "ax.scatter(x, y_prior, color=\"blue\", marker=\"x\", alpha=0.75, label=\"Y prior\")\n", "ax.spines['top'].set_visible(False)\n", "ax.spines['right'].set_visible(False)\n", "ax.tick_params(axis='both', which='major', labelsize=20)\n", "ax.tick_params(length=0, color=\"black\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is of course fairly random, so our prior guess is pretty off charts." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### GP posterior\n", "\n", "We will use a binomial likelihood for constructing the posterior:\n", "\n", "\\begin{align*}\n", "p(\\mathbf{y} \\mid f, \\mathbf{x}) = \\prod_i^n \\mathcal{Bernoulli}\\left(\\text{logit}^{-1} \\left(f \\right) \\right),\n", "\\end{align*}\n", "\n", "where we used a *logit link function*. The inverse of the $\\text{logit}$ is the well-known $\\text{logistic}$ function which we already defined above. We will call the inverse of the link function the *mean function*. The necessity of a mean function is the same as in GLMs: relating the expected value of the response to the predictor.\n", "\n", "The posterior Gaussian process is then given by:\n", "\n", "\\begin{align*}\n", "\\text{posterior} & \\propto \\text{likelihood} \\times \\text{prior},\\\\\n", "p(f \\mid \\mathbf{y}, \\mathbf{x}) & \\propto p(\\mathbf{y} \\mid f, \\mathbf{x}) \\ p(f \\mid \\mathbf{x}).\n", "\\end{align*}\n", "\n", "Since the prior is not conjugate to the likelihood, the posterior does not have an analytical form and thus needs to bei either approximated deterministically, e.g using Laplace's approximation, expectation propagation or variationally, or stochastically using sampling. I'll go over two different approaches:\n", "\n", "- the Laplace approximation (which is one of the typical textbook examples),\n", "- a Markov Chain Monte Carlo sampler." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Laplace approximation to the posterior\n", "\n", "Since the Gaussian process prior distribution is multivariate Gaussian, the posterior is often also fairly close to a multivariate Gaussian, unimodal and roughly symmetric. Note that this implies that for more complex likelihoods the normal approximation should be a fairly bad choice.\n", "\n", "We shall use a Taylor expansion around the mode $\\bar{f}$ of the log-posterior $\\log p(f \\mid \\mathbf{y}, \\mathbf{x})$, i.e.\n", "\n", "\\begin{align*}\n", "\\log p(f \\mid \\mathbf{y}, \\mathbf{x}) = & \\ \\underbrace{\\log p(\\bar{f} \\mid \\mathbf{y}, \\mathbf{x})}_{\\text{constant}} \\\\ \n", "& + \\underbrace{(f - \\bar{f}) \\left[ \\frac{\\partial }{\\partial f} \\log p({f} \\mid \\mathbf{y}, \\mathbf{x}) \\right]_{\\mid \\;f = \\bar{f}}}_{0} \\\\\n", "& + \\frac{1}{2}(f - \\bar{f})^T \\left[ \\frac{\\partial^2 }{\\partial f^2} \\log p({f} \\mid \\mathbf{y}, \\mathbf{x}) \\right]_{\\mid \\;f = \\bar{f}}(f - \\bar{f}) \\\\\n", "& + \\dots\n", "\\end{align*}\n", "\n", "The first line is a constant in the mode $\\bar{f}$. The first partial derivative evaluated at $f = \\bar{f}$ is zero, because it is the mode of the function, so we can safely omit it. The other higher order derivatives can also be omitted due to the fact that they have less relative importance for large $n$. So let's take a closer look at the second partial derivative. We observe that\n", "\n", "\\begin{align*}\n", "\\frac{1}{2}(f - \\bar{f})^T \\left[ \\frac{\\partial^2 }{\\partial f^2} \\log p({f} \\mid \\mathbf{y}, \\mathbf{x}) \\right]_{\\mid \\;f = \\bar{f}}(f - \\bar{f})\n", "\\end{align*}\n", "\n", "is proportional to the exponent of a multivariate normal distribution with mean $\\boldsymbol \\mu = \\bar{f}$ and variance $\\boldsymbol \\Sigma = \\left( - \\frac{\\partial^2 }{\\partial f^2} \\log p({f} \\mid \\mathbf{y}, \\mathbf{x})\\right)^{-1}$. So we use this to approximate the posterior as:\n", "\n", "\\begin{align*}\n", "p(f \\mid \\mathbf{y}, \\mathbf{x}) & \\approx q(f \\mid \\mathbf{x}, \\mathbf{y}) \\\\\n", "& = \\mathcal{N}\\left(\\bar{f} , \\left( - \\frac{\\partial^2 }{\\partial f^2} \\log p({f} \\mid \\mathbf{y}, \\mathbf{x}) \\right)^{-1}\\right).\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to find the mode $\\bar{f}$ of the posterior $p(f \\mid \\mathbf{y}, \\mathbf{x})$, we need to maximize it (or its logarithm):\n", "\n", "\\begin{align*}\n", "\\log p(f \\mid \\mathbf{y}, \\mathbf{x}) = \\log p(\\mathbf{y} \\mid f, \\mathbf{x}) + \\log p(f \\mid \\mathbf{x}) - \\log(\\mathbf{y} \\mid \\mathbf{x}).\n", "\\end{align*}\n", "\n", "Maximizing the posterior w.r.t. $f$ is independent of the marginal likelihood, so it suffices to maximize\n", "\n", "\\begin{align*}\n", "\\Psi(f) = \\log p(\\mathbf{y} \\mid f, \\mathbf{x}) + \\log p(f \\mid \\mathbf{x}).\n", "\\end{align*}\n", "\n", "We can optimze this using a standard gradient based solver, but since we need the Hessian matrix anyways for the variance of the normal approximation, it makes sense to use Newton's method.\n", "\n", "The Jacobian and Hessian of $\\Psi(f)$ are given by:\n", "\n", "\\begin{align*}\n", "\\frac{\\partial }{\\partial f} \\Psi &= \\frac{\\partial}{\\partial f} \\log p(\\mathbf{y} \\mid f, \\mathbf{x}) - k(\\mathbf{x}, \\mathbf{x}')^{-1}f,\\\\\n", "\\frac{\\partial^2 }{\\partial f^2}\\Psi & = \\frac{\\partial^2}{\\partial f^2} \\log p(\\mathbf{y}\\mid f, \\mathbf{x}) -k(\\mathbf{x}, \\mathbf{x}')^{-1},\n", "\\end{align*}\n", "\n", "where the Hessian of the likelihood depends on the choice of the mean function. As above we will usea logistic mean function, and define a matrix $\\mathbf{W}$ as the negative Hessian:\n", "\n", "\\begin{align*}\n", "\\mathbf{W} & = - \\frac{\\partial^2}{\\partial f^2} \\log p(\\mathbf{y}\\mid f, \\mathbf{x}) \\\\ & = \\operatorname{diag}(\\sigma_1(1 - \\sigma_1), \\dots, \\sigma_n(1 - \\sigma_n))\n", "\\end{align*}\n", "\n", "Having laid the ground-work, we update $\\bar{f}$ in every iteration of Newton's method as:\n", "\n", "\\begin{align*}\n", "\\bar{f}_{t + 1} & = \\bar{f}_{t} - (\\nabla \\nabla \\Psi)^{-1} \\nabla \\Psi\\\\\n", "& = (k(\\mathbf{x}, \\mathbf{x}')^{-1} + \\mathbf{W})^{-1}\\left(\\mathbf{W}f + \\nabla \\log p(\\mathbf{y} \\mid f, \\mathbf{x})\\right)\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having solved this for $\\bar{f}$, the approximate posterior is given as:\n", "\n", "\\begin{align*}\n", "q(f \\mid \\mathbf{y}, \\mathbf{x}) = \\mathcal{N}(\\bar{f}, \\frac{\\partial^2}{\\partial f^2} \\log p(\\mathbf{y}\\mid f, \\mathbf{x}) -k(\\mathbf{x}, \\mathbf{x}')^{-1} ),\n", "\\end{align*}\n", "\n", "The implementation is fairly straight-forward, so we won't do that here. If you want to see a solution in R, you can for instance look at some [of my code](https://github.com/dirmeier/GPy)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's again demostrate this using GPy. First we set the model and the inference method of our choice: the Laplace approximation as explained above. For convenience we redefine the kernel and likelihood from above." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = GPy.core.GP(X = x.reshape((n, 1)),\n", " Y = y.reshape((n, 1)), \n", " kernel = GPy.kern.RBF(input_dim=1), \n", " inference_method = GPy.inference.latent_function_inference.Laplace(),\n", " likelihood = GPy.likelihoods.Bernoulli())\n", "m.optimize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we plot the data and the mean posterior process." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " /Users/simondi/anaconda3/envs/py36/lib/python3.6/site-packages/matplotlib/figure.py:2022: UserWarning:This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEYCAYAAAA9AaOpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4VMX6wPHvbEvdNEICBJDem4IiolRFkCtYsBfsXRAB9er9idgLWCgqFgT0XlHBhqKANCkigqBA6IYOKaT37O78/tjNkoQENskmm/J+nifPSU6ZMwvZfTMz78xRWmuEEEKI2sbg6woIIYQQFSEBTAghRK0kAUwIIUStJAFMCCFErSQBTAghRK0kAUwIIUStJAFMCCFErWTydQXKIzIyUrdo0cLX1RBCCFHC5s2bk7TWDavznrUqgLVo0YJNmzb5uhpCCCFKUEodrO57SheiEEKIWkkCmBBCiFpJApgQQohaSQKYEEKIWkkCmBBCiFpJApgQQohaSQKYEEKIWsnjAKaUGqWUmq6UWqOUSldKaaXUZxW5qVKqqVJqtlLqmFIqTyl1QCn1tlIqvCLlCSGEqH/KM5H5P0B3IBM4AnSoyA2VUq2B9UAU8B2wC7gAGAsMVUr11VqfrEjZQtQmWmuUUsV+Bqp9n8PhwGAwlPu6mrSvpr+Goj8L7ylPF+I4oB0QAjxYiXu+izN4jdFaX6W1fkprPQh4C2gPvFSJsoWoFZ577jnGjRvn/sDTWtOnTx/69OlTrfsmTZpEz549mTRpkk/rUZdfw7hx43juuecq9osizkxrXe4vYACggc/KeV0r13VxgKHEMSvO1l0WEFTa9T179tRC1HYOh0OPHTtWA3rs2LHa4XDoMWPGaNd7Q48ZM6Za9tntdt2jRw8N6B49emi73e6TetTl11Dy/7kuAzbpCsSTynwp7fpLoTyUUgOAlcB/tda3luO6e4APgQ+01veXcnwJMAS4VGu9vOTxXr166fKuhai1JiMjg/T0dLKzs7Hb7eW6XoiqkpycTEZGhvtnq9UKUO37LBYL+fn5VXpPm83GsmXLWLJkCZs2bWLMmDEATJs2zX1eZfb16NGDrVu3eqUsb+8bO3Ysb731Vp3vRlRKbdZa96rWe1ZzAHsDmABM0FpPLeX4DOBh4CGt9Xslj3frca5e+MMKj+uptaYgJw0TNsIjIggMDMJoNNb5XyRRe8TGxrq/79Spk2/2dexE7M6qK79jx47Y7XZ2795FSkoKC7/5gccmPo1Sil49e7rP27R5M0DF9m3aRK9evcp/XTXtqw+fOe1bRG932PK7Vuc9q3s1+lDXNq2M44X7w0o7qB3g7+/n8c3yc7JA2WnW7ByMRqPntRSiGpyIjy+WeJCQkABQ7fsOHDxQpeUnJiYCYDSaaNAgkuuuuZKvv5rPicSTBAQGuM+bOXMGQIX23XPv3V4ry9v73n13JuPHPw7U7SBmKPqfXk1q2uNUCv+Hy98sLEVBXhaR4RESvESNcyI+nuTkZCIiImgUHe3+Gai2fdHR0cTFxZGbm4u/vz8tW7YkvhrqcfDgQX7buJn58+dz4403MmHCeKZMmcr8+fMByrVv/PjHufXWW9m9ew/t27fjs88+Y+rUNytUVlXsmzr1TT7//HOAehHEqlt1B7DCFlZoGcdDSpxXKbb8XAKDGnujKCG8ymg0uj/QARpFR5OTk+P+vrr2FY5ZWa1WVDXdMyYmhl7n9aBLly5MmDAeUEyYMJ7t27cDlHtfv379AejXrz9KGSpVlrf3OYMWWK3BSPDyvlqVxNG127l68Yp1Htcz9UQc7dq1rxf9z0JUlKZ6P1q11uzZs5uwRi1K3Lnws6j8+7R2oJThjOf4dl/d/wxq2SQi1paf27k671ndLbCVru0QpZRBa+0oPKCUsgJ9gRxgg7duKMFLiDOr7nfIqfdkyTuXVhPP9hUPXpUrq3r2CW+okkE3pZRZKdXBteqGm9Z6P7AUaIEz27CoyUAQME9rnVUV9RJCCFF3eNwCU0pdBVzl+rGRa9tHKTXH9X2S1nqC6/sYYCdwEGewKuohnEtJTVNKDXad1xsYCOwBninfSxBCCFEflacLsQcwusS+Vq4vcAarCZyF1nq/UqoX8DwwFLgCOA5MAyZrrZPLUSchhBD1lMcBTGv9HPCch+ce4Awdv1rrw8Cdnt5bCCGEKEmeByaEEKJWkgAmPPLC85Pxs5jcX19+8cVZrxk58spi1xw4cKDqKyqEqDckgIkKmTtvzhmPHzt2jGVLl1ZPZYQQ9ZIEMFEukZGRBAUFsWL5cg4fPlzmef/97FPsdjvntGhRfZUTQtQrEsBEuQQFBXHNNdficDj47NN5ZZ43d+5cAG6/7fbqqpoQop6RACbK7fbRztkU8z6dR2lLka1bt5a9e/fQslUrLr7kkrOWt3r1Km679RbatG5JiDWIhpER9L3oQqZMeYOsrNLntOfk5LBo0SIefOB+zu/Vk5gmjbAGB9LinGaMuvYafv75pzLvN2/eXPwsJtq1dc6z//PPzdx8042c07wp1uBA2rdvy8SJ40lJSfHkn0MI4SMSwES5XXJJP1q1bs0/+/ezdu2a047PK9L6OtNSXjabjQcfuJ8hl13Kl19+weHDhzGbzWRlZbFp0yaeefrfXNj7Ag4ePHjatV999SWjrr2a2bM/5u+//yIrKwuTycTx48dZtOh7Ro64kiefnHjW1zL/88/pd8nFLFy4gJycHGw2Gwfi4pj2zjsMGjiAzMzMcvzLCCGqkwQwUW5KKXfX4Nw5c4ody8rKYsGCrzAYDNx2e8l578U9+eQTzJ79MdHR0UybPoPjJxI4mZxKWnomS5f9Qo8e57Jnz25uuP46HA5HsWtDQ8O4+557WbrsF44djyc1LYOU1HQOHDzMs89Owmw28/Zbb7Fo0aIy75+YmMh9993Dbbfdzr79cSQknuRkcipvvzMNs9lMbOwOpk55o2L/SEKIKlfTngdW7fo+PN/XVfCqdTNvrJb73Hrb7Tz//GS+/nohb78zjeDgYAAWLPiKzMxMBl96Kc2aNeOff/aXev2O7duZOWM6gYGBLF78M126nnqQq9lspn//AfyyfAU9undly5Y/WbRoESNHjnSfM3LkyGI/F2rcuDHP/Of/CAgM5N9PPcnMmdO58sorS61DdnY2t912O++9P8u9LzAwkAcffIi4uH945+23+eLLL5j03OQK/RsJIaqWtMBEhTRr1oxBgwe7W1yF5s2dA8Ado8+80Monc2ajtWbYsCuKBa+irFYrV45wBqlly8qXkj9s2BUA/L5hA3a7vczznvr306Xuv/LKEQDs37eP7Ozsct1bCFE96n0LrLpaLHXR6Nvv4Jdly5g75xPuuONO9u3bx9q1awkLC2NEKa2jotavWw/AkiU/07xZTJnnFY5BHTp0+jhYfHw8s95/j19+WcbevXtJS0s7LVhlZ2eTkpJCZGTkaddHRETQpk2bUu/buHET9/cpKSkEBgae8fUIIapfvQ9gouJGXnUV4eHhrF+/nj179rjT6m+48Sb8/f3PeO3x48cAZ4DyJFEip0QraMOG3xg54kpSU1Pd+4KDgwkMDEQphd1uJykpCXCOy5UWwAqfRlwak+nUW6OgoOCs9RNCVD/pQhQV5ufnx/U3OFuwcz6ZzX//+xkAo0efOXkDcLeUXnrpZfLybWf9WvbLCve1NpuN22+7ldTUVLp378F33y8i6WQKJ5NTOXzkGIcOH+XXNaee3F2Rp44LIWo+CWCiUgqD1fTp0zhy5AidO3ehZ89eZ70uOtr5SLnt27eX+54bNvzGwYMHMRqNfPPtdwwdOuy01lR8/IlylyuEqF0kgIlK6dmzF126dCU/Px+A0Xfc4dF1fS66CICfflpc7rlWRw4fAaBhw4bExJQ+frZi+fJylSmEqH0kgIlKe+nll3ls3DgeGzeOm2++xaNr7r7rbpRSpKam8tRTT57x3IKCgmJBLiQ0BHAmccTHx592/pEjR5g5c0Y5XoEQojaSACYqbejQYbz22hu89tobNGzY0KNruvfowaNjxgDw4QezuOnGG/hr61b3eJXdbufvv/7i5ZdepGOHdvz111b3tX37XkxQUBBaa265+Ub27Nnjvmbp0iVcdtngM64AUhs5NOQV2MnOs5GdZyOvwI5DxvZEPSdZiMJnXn31dbTWTJ82ja+/XsjXXy/E39+foKAg0tLSsNls7nOLBqTQ0FBefe11Hn3kYdasWUPXLp0IDg7GZrORm5tLZGQkH3z4Eddec7UvXpbXaCA7t4C0rHxy8gooLV6ZTQb8LSYC/c0E+Zsw1LHALcSZSAATPmM0Gpky5U1uveU2PvjwA9b8+itHjx4hLS2N8PBw2rZty+DBlzJixEi6de9e7Nr77rufZs2a8dabU9m8eTM2m40mMTEMHTqMiROfcI/J1VYaSEzNJj3L+TqUAovZiMmoQIPNobHZ7BTYHBTY8snIzkcpCPQ3ExxgIcjfjEFimajjVG1KMe7a7Vy9eMW6s5/oknoijvbtO1RhjYTwPg0cP5lFdm4BBqUID/EnJNCCsURE0kB+gZ2cfBtZOQXk5tvcrTSjUREa6EdIsB+mGhjJdu/eRVijlr6uhvCilk0iYm35uZ2r857SAhOihknLyiM7twCjUdG4QTD+ZmOp5ynAz2zEz2wkLMgPm0OTlVNAenYeefl2kjNyScnMxRpgITzEH7NRhrxF3SIBTIgaxO7QpKTnAhAVFlhm8CqNyaAIDbIQEmQhN99GWmYeWbkFpGfnk5GTjzXQjwirHyYJZKKOkAAmRA2SnJGL3aEJ9HMmZlSEAgIsJgIiTBTYHaSk55KRk096Vh4Z2XmEBfsTbvWThA9R68mfYkLUEA4N6Vl5KAUNQgPwRngxGw1EhQfSLCoEa6AFrSElI5dD8emkZ+dTe0bAhTidBDAhaojCJAyLa1zLmywmA9HhgTSNsuJvMWGzaxJSsjmamEleQdmPmxGiJpMAJkQNkZ3rXPU+0K9iXYee8DcbiWkYTHREECajgdx8G0cSMziZnotDmmOilpEAJkQNkZ3nnLgd6F+1Q9MKsAaYaR5tJTTYD3B2Kx5JyCAn33bmi4WoQcoVwJRSTZVSs5VSx5RSeUqpA0qpt5VS4eUs52Kl1Heu63OVUoeUUouVUkPLV30h6gabXZNfYMdgUPhbqie3yqAUDUMDiIm0YjEbybfZOZqYSWJqjixTJWoFjwOYUqo1sBm4E9gIvAX8A4wFflNKNfCwnAeBNcBg1/YtYDXQH/hJKfVMeV6AEHVBdp6z+zDAYvJK8kZ5+FuMNIuyEmH1RynnPLTDCRnk5svYmKjZytMCexeIAsZora/SWj+ltR6EMwC1B146WwFKKTPwCpAL9NRa36a1/rfW+jagF5AHPKOU8ivvCxGiNstxdR8GVHH3YVkUEBHiT9OGVvzMRgpsDo4mZZCckSuZiqLG8iiAKaVaAUOAA8DMEocnAVnAbUqpoLMUFQGEAnu01ruLHtBa7wT2AAFAsCf1EqKuKMwEDKim7sOy+JmNNI2yEm71ByA5PZejiZkU2B0+rZcQpfG0BTbItV2qtS72m6y1zgDWAYHAhWcpJwFIBNoppdoWPaCUage0BbZqrU96WC8haj0N2FwBoiaskqGABiH+NGkQ7M5UPJyQQXp27V4gWdQ9nr5b2ru2e8o4vte1bXemQrRz5eCHXffdrJSaq5R6RSk1D+f42g7gOg/rJESd4HBoHA6NwaAw1KCFdwP8TDSLshIcYMHhcM4bS0jJlnR7UWN42l8R6tqmlXG8cH/Y2QrSWn+llDoGfA7cXuRQPPAJzsQQIeqNoq2vmhO+nIwGRXREIIHZJpJSc0jPzie3wE6jiCAsJt+3FkX95q3fwML33Vn/NlNK3Qr8gjMDsSPOrseOwHJgBjDfS3USolYoHF+qqavFKyAk0EJMlBWLyUh+gZ0jiRlk5BT4umqinvP0HVPYwgot43hIifNK5Rrnmo2zq/A2rfUurXWO1noXcBvObsTrlFIDPKyXELWeuwVWw1s0fiYDTaOsWAOdXYrxyVkkpOZIl6LwGU/fMYUZg2WNcRUmZJQ1RlZoCGAGVpeSDOIAfnX92NPDeglR69lsNSeB42wMCqLCA2kYFohSzsWHjyZmSJai8AlP3zErXdshSqli1yilrEBfIAfYcJZyCud3NSzjeOF+SXeqYV54fjJ+FlOxL38/M5ENwmnV8hz697uYMWMe5euFC8nPr5r/vtTUVF54fjIvPD+Z1NTUKrmHL9T0LsSSFBAaZKFpQytmk4G8AjuHEzLIzJUuRVG9PHrHaK33A0uBFjizCIuaDAQB87TWWYU7lVIdlFIdSpy7xrUdpZTqVvSAUqoHMArnONoKT1+AqH7R0dFER0cTFRWFUopjx46xYcMGZr3/HjfddAMtzmnGrFnvo728HFFqaiovvvgCL774Qp0KYLWlC7GkwjljwQFmHA7NiZNZnEyXic+i+pRn1uRDwHpgmlJqMLAT6A0MxNl1WHIJqJ2urTuxSmu9USn1Cc7lqP5QSn0DHMQZGK8CLMDbWusd5X8porocOny02M92u52dsbH8svwX3nvvXQ7ExTHm0UdYv24dc+bOQ8mDE8+oNnUhlmRUiuiIIPwy80hOzyElI5e8AjvR4YEYa9CUAFE3efyOcbXCegFzcAau8UBrYBrQpxyTj+/GGcB+Ay53lXMZsBa4SWs9ztM6iZrBaDTSpWtXHntsHFu2/MX1198AwPz5n/PG66/5uHY1m0Nr7A6NQVFrP/AVEB7sR+MGwRgNiuzcAo4kZshzxkSVK9effFrrw1rrO7XWjbXWFq31OVrrsVrr5FLOVVrr096R2mmO1nqA1jpca23SWkdorQdrrSWFvpYLDAzk49mf0KPHuQC88cbrJCef+vVwOBysX7+OZ57+N5dcfBGtWp5DcFAAjRtFcenggXzwwSwKCk4fS7ns0kG0b9fG/XP7dm2KjcdddumgSt/DFwrszg43k9FY4+aAlVegn4mmUVb8LK61FBMzJdVeVCnfLrwm6iSLxcKTTz7FTTfdQHp6Ot9//x133HEnAIcOHWLggP7uc00mE4GBgSQnJ7NmzRrWrFnDF/M/54cffyIgIMB9Xnh4BJGRkSQlJQEQGRmJ0WgsdrxQRe/hC+7uw1o2/lUWs9FATKSVpLRs0rPyiU/OIi/YjwahAbU+QIuap268a0SNM+Tyy90BZs2vv7r3m0wmrrxyBP/97+fEHThERmY2iUnJnExO5cOPPqZJkyasXbuWZ5/9T7HyvvxqAevWn0pyXbd+A4cOH3V/ffnVgkrfwxdq0hqI3mJQ0DDsVKp9amYex5IyscmEMeFldeddU0OVzMTzdmZeTRUcHEzLlq0A+Oef/e79TZs2ZcHCrxl13XU0adIEg8HgPv/220ezYOHXAHz80Ufk5uZW6N7VcQ9vsdWyFHpPFabax0RaMRkN5OTZOCLPGBNeVrfeNTXMC89PZsKE8e6gpbVmwoTxvPD8ZB/XrHpERDgf1J2cnOLxNT179iIqKoqsrCz++mtrldSrOu7hKbvDGcCMxrrZweZvcabaB/iZsNmdzxhLk1XthZdIAKsiWmtS09KYMX2aO4hNmDCeGdOnkZqWVi9aYmW9xvz8fD74YBZXXDGUFuc0I8QaVCwhIyEhAYCjR46Wer0nquMe3uBwdasZ6vBUA5NB0SQymLBgP7SGxJRsMrLzyZcsRVFJksRRRZRSTJkyFYAZ06cxY/o0AB55dAxTpkytF3OjUlKck40bNDiVYJGQkMCwoZezffs29z5/f/9iSRmJiYk4HA6ysrKoiOq4h7e4A1gtTaH3lAIiQwPwsxhJTMkhN9/OMx/8ysSbexMVHujr6olaSlpgVahoECtUX4JXZmYmcXHOJ+O0atXavX/ihPFs376NBg0a8MGHH3Hw0BHS0jM5euyEOyGjSZMmAOgKrulQHffwFruuHwGskDXAQkxD53yx/UdTmThzJVv3Jvi6WqKWkgBWhQq7DYsqOiZWly1dsgS73dlF1K+/M6W9oKCAb7/9BoC3357G6NF30KhRo2LX2e12d6p8RVTHPbypcEnrutyFWJKf2Ui41Z/z2kWTkZ3Pi3PW8fXq3fXifSG8SwJYFSk65vXIo2PIzSvgkUfHFBsTq6vy8/N57bVXAQgNDWXEiJGAs9uuMOuve48epV67bt3aMjMDC7MJoezxtcreo7o5XBGsvrTACikFT9/eh+sHdcCh4bMlsbz+39/JlgWBRTlIAKsiSinCQkOLjXlNmTKVRx4dQ1hoaJ3tRszJyeGeu+9i69YtADzxxJOEhTkf1B0SEuJ+3dv+/vu0a202G5Oe/b8yyw4JCXF/n1bGYr6VvUd1qw9JHGUxGBQ3XtqRp2+/kEB/M7/HHufJd1dxOD7d11UTtYQEsCr0f89OKjbmVRjE/u/ZST6umXc5HA52bN/O22+/RY8e3fjiC+eKYLfccivjJ0x0nxccHMxFF10EwBNPTGDlyhU4XGnkO7ZvZ8SIf7F582aCgoJKvU9YWBgxMTEAzJ03F5vNdto5lb1HddKAQztbI/Uwfrn16tCYNx4ewDmNQjialMmT761i/TbfZoeK2kECWBUr2dKqCy2v5s1i3F/RUZEEBfpz3nk9ePKJiRyIiyMyMpIZM99l9idzTnu9U6a+SVBQEEePHmXo5UMIC7US2SCc887rwepVq3jvvVlERkaWee97770PgHdnziAiPJQ2rVvSrm1rbr3lZq/do7oUbX3V/t+KymncIJhXHujPJd2bkptvZ8rnG5n70zbs8qBMcQYSwES5xcfHEx8fT0JCAjabjUaNGtG7d2/uu/8B5s//krgDh9yBpqTzzuvJunW/MWrUdURGRuJwOLBarYwadR2rf13DLbfeesZ7P/nUv5n65lv07NkLs9nMkSNHOHjwIPHxJ7x2j+riqGcZiGfjbzHx2PW9uGt4V4wGxXdr9jH5k3WkZub5umqihlK1KZmga7dz9eIV6zw+P/VEHO3bl3ymphA1Q+GTjP3MRppFWX1dnWq1e/cuwhq1LPN4bFwSUz7fSGpmHg1CA5h48wW0axZR5vnC91o2iYi15ed2rs57SgtMCB+pL5OYK6JTy0imPDKQ9s0jOJmWw38++JVF6/bV6exdUX4SwITwkfo2ibm8IkICeP6eSxh+UWtsds0nP27jtc9+JzNH1lIUThLAhPCR+pxC7ymzycDd/+rGE7f0JtDfzMadxxk/fSV7Dp/2DF1RD0kAE8JHJInDcxd2bsLURwbSpmk4ianZ0qUoAAlgQviMtMDKJzoiiJfu6yddisJNApgQPiJJHOVXWpfi49NXEBtXM9a2FNVLApgQPiJdiBVXtEsxKTWHZz9aw/+WxrqfcC3qBwlgQviIdCFWTnREEC/f349rB7RDAwtW7ebpWb9yLCnT11UT1UQCmBA+4opfGKUFVmEmo4FbhnTm+XsuITI0gH1HUpgwYwW/bDogCR71gAQwIXxEWmDe07llJG+NGcTF3ZxrKb779Rbe+N9GMrJlGaq6TAKYED4iSRzeFRRgYdwNvRh7XU8C/Exs2HGMcdNWsHn3ibNfLGolCWBC+Ij7YZbSAvMapRT9z23Om48OosM5ESSn5/LS3N+YvmAzWZJuX+dIABPCR1yPKcMg70Kvi44I4oV7+zF6WBcsJgMr/zzE2HeWS2usjpG3jhA+4HyYpXY9zFJaYFXBaFCMvKQtUx8dRPvmxVtjMvm5bihXAFNKNVVKzVZKHVNK5SmlDiil3lZKhZf3xkqprkqpeUqpw66yEpRSq5VSt5e3LCFqG3mYZfWJaWjlxfuKt8bGvPULa/46LJmKtZzHAUwp1RrYDNwJbATeAv4BxgK/KaUalKOsO4AtwFXAGmAqsABQwBWeliNEbSWTmKtX0dZYx3MakJqZx1tfbOKFOes5cVLmjdVWpnKc+y4QBYzRWk8v3KmUehMYB7wEPHC2QpRSFwIfAduBoVrrEyWOm8tRJyFqJUmh942YhlZeuPcSVvx5kHk/bWfr3gQee2c5owZ1YOTFbTGbZFSlNvEogCmlWgFDgAPAzBKHJwH3AbcppcZrrbPOUtzrgBG4tWTwAtBaF3hSJ286mpBOboG9um/rdf5mIzFRIb6uRo21fds2Xn75JdatW0tiYiJ2u51u3brzx6bNrF69iiGXXQpAXr6t3GWX93ppgfmOwaC4tFcLenVozNzF21i99TD/WxrLr1sOc9e/utGjbZSvqyg85GkLbJBru1RrXWyxMa11hlJqHc4AdyGwvKxClFJNgUuATcAOpdRAoCfOMe2twMqS5VeH3AI7/v7+1X1br8vNza2W+9jtdr75+msWL/6R33//ncTEBLKzswkLC6Nt27b07XsxN910M527dKmW+ngiLi6OAQP6kZGRAUBERARms5kGkR73fHuVtMB8LyzYj7HX92LAec354LutHEnM4PlP1nF+x0bccUVXGjcI9nUVxVl4GsDau7Z7yji+F2cAa8cZAhhwfpHzVwADShzfppS6Rmu9z8N6iWr2++8buPuuu9i799Svgtlsxmq1cvLkSRITE1m/fj1vvPE6V111NZ9+9l8sFosPa+z00YcfkJGRQes2bVi69BeaNm1a7HhgYCDt2rUv42rvs8sk5hqje5so3h47mEXr9rNg5W7+2HmCLXsSuLJva0YNbE+An4xq1FSedviGurZpZRwv3B92lnIK2+bXAx2Ba1xltwE+BboCPyqlfP+JJ07zww+LuOzSwezdu4cGDRrw4osvsX3HTjKzcjh+IoHMrBzW/7aBiROfICQkhG+//Ybs7GxfVxuA7Tu2A3DllVeeFrwAzj//ArZt38G27TuqpT6FyW/SAqsZzCYj1/Rvx4zHL2Pgec2x2R188+teHp66jKUb42SV+xqqPEkcZ1L4LjxbTqqxyPYerfUPrp/TlVKjcQa1XsC1wOdeqpvwgr1793LnHaPJy8ujY8dO/PDj4tMCgdFopGfPXvTs2YvHx0/gvnvv8VFtT5fjCqTBQTWjW6hwDEziV80SEeLPo6N6MrR3Sz7+4W/2HE7h/W+38t2avdx0WScu6hIjreYaxNMWWGELK7SkK/99AAAgAElEQVSM4yElzitLimubBywuekA7J2R85/rxAg/rJarJc5OeJT09HX9/f778akGprZiiIiIiWLDwa0JDi//KnDhxgqeeeoIe3bsRER5KeFgIPbp3499PPUl8fHypZR04cAA/iwk/i4kDBw4QHx/P44+Po127NoRYg2jWtAm33nIzu3btOu3adm1b42cxsXr1agBefPEFd1nO/asAZxJG4b6y7Nq1i9G330bzZjGEWINo374tjz02tsx6l5Sbm8uMGdO5dPBAOrZpRvd2jenRuTWjrr2GJUt+LvO6onXNyMhg0rP/R9cunQkNCaZxoyiuumoEGzf+ftb7L1u2lFtvuZm2bVoRGhJMo+iG9DzvXB57bCwbNvxW6jV2u5158+YyfPgwmjVtQnBQADFNGjF8+DC+/OKLOjuPqm2zCF6+vz+P33g+jRsEcfxkFm/O/4OJ765ky574Ovu6axtPW2C7Xdt2ZRxv69qWNUZWspyMMpI1CgNcgIf1EtUgPj6er79eCMBNN91Mu3Zl/RqcrugqE7/+uprrRl1Lamoq4Bx3Ukqxc2csO3fG8skns1n49Tf07XtxmeXFxu7g/vvuJSEhgcDAQAASEhL46qsvWbLkZ5YvX0m37t3d50dGNiQ3N5fk5GQKCgoICgoiOPhUK8zT8bklS37mulHXkpfnXN08ODiYE8eP8967M/n2m6+Z/PwLZ7x+7969XDVyBPv27XX/uwQHW0lMSGDRou9ZtOh77rv/AaZPn1FmGSeOn6B37/PZv28f/v7+GAwGkpOT+WnxYn5Ztoyvv/mWIUMuP+267Oxs7rn7LhYuXODeZ7Vayc7OZvv2bWzfvo11a9fyx6bNxa6Lj49n1LVXs3HjRve+0NBQkpKS+GXZMn5ZtowvvpjP5/O/qBHjnN5mMCgu7taUCzs3YeWfB/li+S7ijqXxwpz1dGrRgFED29O9TZSspOJDnrbAVrq2Q5RSxa5RSlmBvkAOsOEs5fwNJAGRSqnoUo4Xpq0d8LBeohqsXrUKh2vhvhEjR1aojMOHD7uDV8eOnVi5ajUpqekkp6SxfMVK2rVrT0pKCqOuvYajR4+WWc5dd95BmzZtWP/bBvf1i3/6mcaNG5Oens64cWOLnb/+tw0cOnyUPn36ADBu3OMcOnzU/dWnz0VnrfuRI0e49ZabycvLo2vXbqxdt56TyamkpKbz/aIfMBqNPDFxQpnXp6amMnz4MPbt28uAgQNZvmIl+w4lsuHvf9h34BivvzGF4OBgPpj1PtOnTyuznLFjH8VitrBk6TL3a1+3/jfatWtPQUEBjzz8kPv/qah777mbhQsXYDAYmDBhIvv/OUDSyRTS0jP5J+4gc+d+Su8LLyx2TX5+PtdcPZKNGzdy7rnn8e1335OSmk5C4kmSU9L4+ONPiIqK4ocfFvH000+d9d+wNjMZDVx2fktmjh/C7UM7ExxgJvbASZ7/ZD1PvreajbHH3Fmlonp5FMC01vuBpUAL4OEShycDQcC8onPAlFIdlFIdSpRjA2a5fny9aDBUSnUF7gBsOFflEDVEbOypxIYePc6tUBmvvfYqqamphIeH8/OSpVx0UV/3sYsvvoSffl5CSEgIycnJvP76a2WWExUVzeKfltCzZy8ATCYTgwdfyoyZ7wKwdu1ajhw5UqE6nqnu6enpNGjQgMU//cz55zt7uA0GA5dfPpTvv/+BrKyypz+++urLHDxwgAEDB/Ljjz9x8cWXuFssoWFhjB37GLNnz3Ge+8rL2GylzyMzmUwsXfYLAwYMxGAwoJSiV6/z+fzz+QAcPHjwtK7AFSuWs2DBVwC8/c40Xnr5FXf3r1KKmJgYbrzpJmbMKD698+OPP2LTpk106tSZZb8sZ9iwK9wt3qCgIG697Ta++34RSilmvf8+CQkJ5fknrZX8zEau6teO9ydezq2XdyIkyMK+Iym8+tnvPD59BWv+Ooxdkj2qVXmmnT8EJADTlFLfKqVeUUqtwLkKxx7gmRLn73R9lfQyzpba7cAmpdSbSqlPgd8Bf+BJSaOvWU4mJ7u/j4iIKPf1WmsWuj5E7733Pho1anTaOU2bNuXee+8D4KsvvyizrMfGjSMg4PQe5qFDh7mDwvbt28pdx7JorVnw1ZeAs+5RUadPcu3cpQvXXHNtmdfPnTMHgMceG4fJ5Oy1P5XE4ex+GjFyJCEhISQlJfHnn5tLLevuu+8p9f5dunalRcuWAGzbVvy1F967U6fO3H//WRfKcftk9mwA7n/gAaxWa6nnnHdeTzp16kx+fj6rV63yuOzaLtDfzDX92zNr4uXcNbwrESH+HIpP560vNvHAlKV8vXq3PEizmnichai13q+U6gU8DwzFuWbhcWAaMFlrnXym64uUk62UGgw8AdyIs0WXC6wHpmqtfyrfSxBVrbID1nFxcSS7guCgwYPLPG/wpZcydeoUTp48SVxcHC1dH8pFXXB+6fk9JpOJhg0bcvToUVKSU0o9pyKK1n3AwIFlnjdg4EC++GL+aft3xsa6r7/3nrsxuJ6dYndotHau0Vc4hJKZ6VyT79DBQ1xwQe/Tyjr/grJzm5o0bsyBuDhSkou/DX9ztciGDx9e5rUlZWRksG3b3wBMfm4SL7/0YpnnFr62g4cOelx+XeFnMfGvvm24vHdLVm05xHdr9nEsKZPPlsTy5fJd9OvRjOEXteacRmXlvonKKlcavdb6MM7FfD05t8yRTa11NvCc60vUcJENTq1WkZycTJMmTcp1fWLiqe6lJk1iyjwvJuZUZmNiYkKpASy4jNYA4G7dFNi8txqZ53Uv/dix48eKlJXo0T2zc0qfO1dWSwiKvPaC4q89/oRztbbm55zj0b3BmSlaOJaWnOzR36XuaQr1kdlk5LLzWzK4Zwu27kvgx/X72bInnl82HeSXTQfp1KIBg3qeQ58uMQT4eWvmkgDvzQMTdVinTp3d32/duqXcAawoTzO2amJmV0XqZLefWmPz0OGjREc7c5cOJ2aQl2+naZQVf7OxrMsrrbDO5al70TqvWbuu1NagOJ3BoDivXTTntYvmaGIGi3/7h5V/HiT2wEliD5zko0V/07dbDIN6nkOH5hE18ne8tpGll8VZ9R8wwN319f13353l7NM1bHhq3ObIkcNlnnf06Knki8jIhuW+T1UoWvei9Svp2LFjpe5vFH1qvK/o2Fxhr2xVf4hFu8YbDx444Pk10acShLdv3+7tKtULMQ2t3DuiOx//exgPXXMu7ZtHkJtvY/mmgzwz61ceeXMZ/1sWy8ETaTKnrBIkgImzio6O5uqrrwFg/vzP2bPnbNP9TtFa07JlS3fyx8oVK8o8d8Vy5zKaDRo0KLX70BeK1v1MiQqrVq4sdX/nLl0ICXHO8//qyy/d+92r0VfxH+F9LnROH/jxxx89viY8PJyOHTsB8OUZEmrE2QX4mbm0VwteeaA/08ddytX92hJm9eP4ySwWrNzNuGkrGPPWL3wuwaxCJIAJjzw3+XmCg4PJycnhhuuvO+NcLYCUlBRuuP460tLSUEox6rrrAfjoow85ceK0p+hw7NgxPvroQwCuv+FG77+AClJKce2o6wD48MMPSEpKOu2cnbGx7oneJZlMJkbfcQcAn346j3Xr1gKgHcWzEAt5OubkqTvudA5Zx8buYNas9z2+7u57nMuArVyxgi+/OHMQ83ad66qYhlZuG9qFD58YyqS7+nLZ+S2wBlo4mpTJV65g9tDUZXz4/V9s3n2iQo/1qW8kgAmPtGvXjtmfzMVisRAbu4MLzu/JlDdeZ9++UzMe7HY7W7dsYfJzk+jQvi3ffvuN+9iTTz5FWFgYycnJDBt6Ob/9tt59bP36dQwbejmpqalEREQwceIT1frazuaJJ57EarWSlJTEFcOGsnnzJsDZuly2bCkjRvzLPUeqNE8//R9atW6NzWbjyn8N5+233+KkKxAalCItLY0lS37m7rvuZNDAAV6t+4ABA7n++hsAeGzsGP7zzNPueXJaa44dO8bs2R9z/333Frvuvvvu5wJX1uOdd45m0rP/x+HDp7p/s7OzWb16FWPHjqFjB89XZhFgNBro3iaKB68+l9n/HsazdzqDWUighfjkLH7a8A8vzf2N0S/+yItz1vPj+v3EHU+TydKlkCQO4bGRI0eyZOky7rnnbvbv28czzzzNM888jcViITg4mNTUVHf2mlKKG264kaCgIMA5z+urBQsZde01xMbuYED/fu5jhZOAw8LC+GrBwjIz+nylefPmzPv0M2684Xr++msrF/W5EKvVis1mIycnh8aNG/P6G1PKXLw4IiKCxYt/5vrrRvH333/x5BMTgYmEhIQCmvT0dPe5rdu08Xr9Z33wIfn5+Xz77Te88cbrvPHG64SEhJCXl+deGqtbt+7FrvHz8+Obb7/nlltuYtXKlbz66iu8+uorhISEYDAYSEs71d1VmAEpys9oNNCjbRQ92kZx38ge7DuSwpY98WzefYL9R1P5c088f+5xrrUZ5G+mY4sGdG4ZSeeWkbRoHIrJWL/bIPKbh/NJxtX1MMiqVJXZbIUuuqgv27btYOGCBSxe/CMb/9hIYkICGRkZRERE0L59ey65pB8333Ir7dsXf75Wv379+XvbDt56ayo///QzBw8eQClFhw4dGXbFMB577PFSJznXBFdcMZzff/+Dl19+iVWrVpKWlkajxo254orh/PvfT7NrV2lz9k9p2bIlv234nS/mz2fBgq/4Y9MmUlOSMRqNtGjZku7dujP8X/9i+PB/eb3ugYGBfPHlVyxe/CNzPvmEP/7YSFJSElarlXbt2tOvf39uuumm066LjIzk55+X8sMPP/C//37GH39sdK+4ERMTQ+cuXRg27ApGjKjY8mKiOKNB0b55BO2bR3DjpR1Jzchly954/t6fyI64JJJSc9i06wSbdjm74M0mAy0ah9ImJpw2TcNpExNGk4ZWjPVotXxVmwYNu3Y7Vy9esc7j81NPxNG+fYeznyhENbI5NAeOp2E0KFo2rp+TXHfv3kVYo5qRqFNbJKRksyMuidi4JGIPJHH85OnLl/mZjTSNstI8OoTm0SE0iw6heZSViJCAKn8MTMsmEbG2/NzOZz/Te6QFJkQ1K/yjUZ4rJcojKjyQqPDmDDyvOQCZOfnsP5rK/qMp7DuSyr4jKSSl5bj2pRa71mIyEBUeRKMGQUQXbiOCaBQRSIPQwFo7wbp21lqIWsxRTXPARN0WHGChe5sourc5NVcxMyefw/HpHIrP4FB8Oofi0zmckE56Vj5HEjM4kphRalkBfiYirP6Eh/gTERJAhNWfiBB/QoL8sAZaCA60YA20YA0wE+BnrjF/fEkAE6KaaXkas6giwQEWOraIpGOLyGL7c/IKOJGcTfzJTE4kZ7m/4pOzSUnPISfPxtG8TI4mZZ71HgblvE9woIXgADPBARafteAkgAlRzQrngBkkgolqEuBnpmXj0FLHXLXWZOUWkJyeS3J6LikZOe7v07PyyMwpICM7n8ycfDKy88nJs5GenU96dr4PXklxtSqAJaZl88nibQT7O6N+UICZ4AAzQQGn/hIIDrTUqywcUfucWoVDfk+F7ymlnJ+dARaaR4ec9Xyb3UFmTgGZ2XlkZOeTlWsjO7eA0XOqvq4l1aoAlpFdwKK1Z35UmEFBSLAf4cH+jLkihoTUbIwGA0ajwmQ0YDYaMBkNGAwK+fgQvnBqDMy39RCiIkxGA2HBfoQF+/m6KrUrgEWGBnD70M5k5hSQlVNAZk4+WbkFrp/zXc3cAlIz8kjNyCO/oBHpWaU3cw0G5QxmJoN7azEZsZgMGI0GCW6iymhpgQnhFbUqgIUEWriq35mXrbHZHaRl5pGamYvVkUxUeCB2u8bmcGCzO7DZHBTYHTgcmjyHnbwC+2llGAwKi8mI2WTAz2zEz2LEYjZilA8c4QUln8YshKiYWhXAPGEyGmgQGkCD0ABST6RiDTCf9kGhAYdDY7M7KHAFNJvNQb7NTn6BHbtDk5tvIzcfiiadngpoJvwtRvzMpipfTVzUPYVrB9SUVOTqVpsWTxA1W50LYEUpgxG73X7aWm0K57ItRoMRv1KWX7I7tCuYOcgvcLbS8gpszmBncw5ggnMMw89sIsDPGdD8/UzSShNndWoleh9XxEfsdjvKUPXLnom6r04HMJPFn+ysLEJCy7dcj9GgCLCYCLCc2qeBfJudvHznV26+jXyb3dVScz72oDCgBfqbCPQz42cxyliaOE1970LMzsrCZPH3dTVEHVCnA5jZL4jklGSCgoMxGiv3F58C/ExG/ExGcD05w641efl2cvJs7kBW+JVMLgaDItDPTFCAmSB/kwzaC6BIF2I9/H2w2+0kpyRj9q+fa0AK76rbAcw/EFtBLocPHyIiPILAoCCMRqPX/vI1KkWgn4lA1yx0h9bk5NnIznPOi3B2NzonACqFO5gF+psx1dPxD1G0BebjilQTrTV2u53srCySU5LRRj/M/mU/P00IT9XpAKaUIsAaQUFuNkmp6dgSEtCO07MOq4rdod1jaAU2R7FjFrMRf1d2Yz35HBMu6Vn55BXYyU6xlDoGWxcpgxGTxR+zfyhm/8B6230qvKtOBzBwBjFLQBCWgCCf1iM1I5c/dh3n9x3H+WtfAnbXQH6An4k+XZow4NzmdGoRWW8z0+qTaZ+sY+veBJ4Z3Yee7Wvm88+EqA3qfACrKcKs/lx2fksuO78l6Vl5rPv7KKu2HmLv4RRWbD7Eis2HiI4I4vLeLRncsznWQN/PchdVIy/f2Qvgb5a3nxCVIe8gHwgJ8mNYn1YM69OKo4kZrN5ymFVbDxGfnMW8n7bz+bJYLu7WlKG9W9K2WYSvqyu8LLfAmbXqZ6kf3YdCVBUJYD4W09DKzUM6ccOlHdmy5wQ/bYhjy554Vv55iJV/HqJ1TBhXXdKWC7vEyCLFdYS7BSYBTIhKkQBWQxgNil4dGtOrQ2NOnMxkycY4lm86yP6jqUyd/wfREbGMvLgNA3ueU28G/uuqwuXLLNKFKESlyDuoBmrUIJjRw7py46WdWPXnIb5ds5f45Cw++P4vvli+i+EXtWbYhS0JKjrTWtQaea6J79ICE6JyDOU5WSnVVCk1Wyl1TCmVp5Q6oJR6WykVXtEKKKX6KaXsSimtlHqxouXURX5mI5f3bsmMxy9j/I3n06pJGGlZefxvWSwPvLGUBat2k5Nn83U1RTnluroQ/Szy96MQleHxO0gp1RpYD0QB3wG7gAuAscBQpVRfrfXJ8txcKWUF5gLZQHB5rq1PjAZF325NuahrDH/vT2Thyt1sj0vif0tj+WHtPq7u346hF7aSrsVawG53PhVBKbCYyvX3oxCihPK8g97FGbzGaK2v0lo/pbUeBLwFtAdeqsD93wFCgVcqcG29o5Sie5sonr/3EibffTHtm0eQnp3P3J+289CUpfy04R9sdsfZCxI+Uzj+5Wc2yWReISrJowCmlGoFDAEOADNLHJ4EZAG3KaU8ni2slBoJ3AmMAY55ep1w6tq6IS/f34//jO5DqyZhpGTk8uH3fzHuneVs2nVcHllRQ7m7D6W1LESledoCG+TaLtVaF/sTX2udAazDucTthZ4UppSKAj4EvtVaf+ZhHUQJSinOa9+INx4ewBM3X0DjBkEcTcrk5XkbmDx7HXHH03xdRVFCXoEkcAjhLZ4GsPau7Z4yju91bc/8uORTPnDd+wEPzxdnoJTiwi4xvD32Uu4a3pXgADN/709kwowVzPz6T1Izcn1dReEiCRxCeI+nAazw2Qdl/UlfuD/sbAUppe4CRgIPaa3jPby/8IDZZOBffdswc/wQhl/UGoNSLN90kEfe+oWffvvHvf6i8J2cXOfDUAP8JIAJUVneSoMqHI0+4yekUqoF8Dbwldb6Sy/dW5RgDbRw97+68c5jgzmvXTTZuQV8uOgvnnpvFXsPJ/u6evVatmvaQ6C/2cc1EaL28zSAFbawynoKXUiJ88oyG8gBHvLwvqISmkRaeWZ0H564pTcNQgPYfzSVp95fzfvfbiEjO9/X1auXsl0tsEBpgQlRaZ4GsN2ubVljXG1d27LGyAqdhzMVP9E1cVkrpTTwiev4M65933pYL3EWSiku7NyEaY9dylX92mJQiqUbDzDmrV9Yv+2or6tX77gDmLTAhKg0T/8MXOnaDlFKGYpmIromI/fF2bLacJZy5uHMViypLdAP2ApsBrZ4WC/hoQA/E7cP7cKAc5vzwXdbiT1wkimfb+TCv5tw74juhFv9fV3FeiE7zxnAgiSACVFpHgUwrfV+pdRSnHPBHgamFzk8GQgCZmmtswp3KqU6uK7dVaScMaWVr5S6A2cA+1Fr/Z9yvgZRDs2jQ3j+nktYujGOeT/vYMOOY2z/J5E7h3dlwLnNZXJtFcvOlTEwIbylPB3xD+FcSmqaUmowsBPoDQzE2XX4TInzd7q28olYwxgMiqEXtqJnh0a8/+1WtuyJZ/qCP1n791EevPpcIkMDfF3FOivL3YUoY2BCVJbHWYha6/1AL2AOzsA1HmgNTAP6lHcdROF7DcMC+c/oPjw6qifBAWa27Iln3DvLWfv3EV9Xrc6SMTAhvKdcfwZqrQ/jXP7Jk3M9bnlprefgDIyimimlGHhec3q0jeLdr/9k8+543pz/B3/sPM59I7rLI1u8TAKYEN4jy2ELAMKt/jx9ex8euKoHfmYja/46wmPTVrBtf6Kvq1anFI6BBflJABOisiSACTelFEMuaMnURwfRtlk4J9NymPTxWj5ZvI181yrqonIKsxBlDEyIypMAJk7TJDKYl+/rx42DO2IwKBat3ccT767i4AlZHLiysqQLUQivkQAmSmU0Grh+cAdeub8fjRsEcSg+nSffXcWS3+PkUS2VIGn0QniPBDBxRm2bRTD10UEM7nUO+TYHs77bytT5f7hbEsJzDocmJ0+WkhLCWySAibPyt5h4+JrzeOz6XvhbTKzfdpQJM1ay70iKr6tWq+Tm29Da+Swwo1HeekJUlryLhMf69WjGlEcG0rJJKPHJWTw9azWL1u2TLkUPFa5EHyAZiEJ4hQQwUS5NIoN59YH+XNGnFTa75pMft/HqZxvIyM7zddVqvGxZhUMIr5IAJsrNbDJyz5XdefLW3gT5m/lj5wnGz1jJXulSPCOZxCyEd0kAExXWu1MT95yxpNQcnpn1K0s3SpZiWQoDmKxEL4R3SAATlRIVHsiL9/Zj2IWtsNkdvP/tVmYs/JO8fJuvq1bjuFPoZQxMCK+QACYqzWwycO+I7oy9vhcWs5GVfx7i37N+5cTJTF9XrUaRleiF8C4JYMJr+vdoxmsP9qdxgyAOHE9jwsxV/LHzuK+rVWOcWkZKWmBCeIMEMOFV5zQK5fWHB3JBp8Zk5xbwyqcb+O/SHdgdMi4mSRxCeJcEMOF1Qf5mnrylN7cN7YxBwcJVe3hhzjrSMut3qr17JXoJYEJ4hQQwUSWUUlzdrx2T7rqY0CA//t6XyMSZ9TvVXuaBCeFdEsBElerauiFTHhlI++YRJKU5U+2X/XHA19XyCVmJXgjvkgAmqlyD0ACev+cSd6r9e99s4d2v/6x3zxhzt8AkjV4Ir5AAJqpFYar9o6N6YjEZ+GXTQf7z4RqSUrN9XbVqU7gWorTAhPAOCWCiWg08rzkvP9CfqPBA9h1JYcKMlWzbn+jralWLUytxyBiYEN4gAUxUu1ZNwnjj4QH0aBtFenY+k2ev5ds1e+v8ElSSRi+Ed0kAEz5hDfTjmdEXMWpAexwa5v20namf/+F+4GNdo7WWx6kI4WUSwITPGA2Km4d04qlbexPoZ2L99qM8+e5qjiZm+LpqXpeamYfDoQkOMGM2ydtOCG+Qd5LwuQs6NeH1hwfQLMrKkcQMnnh3Fb/HHvN1tbwq0ZWsEhUe5OOaCFF3SAATNUKTSCuvPjiAi7rEkJNn47XPfq9TS1AlpjgDWMOwAB/XRIi6QwKYqDEC/EyMv+l8Rg/r4l6C6qW56+vE054TXC2whuGBPq6JEHWHBDBRoyilGHlJWybddTEhQRa27k1g4sxV7D+a6uuqVUphCywqTAKYEN4iAUzUSF1bN2TKwwNp2zSchJRsnpm1mpV/HvR1tSosUVpgQnhduQKYUqqpUmq2UuqYUipPKXVAKfW2Uircw+uDlFK3KKX+p5TapZTKUkplKKU2KaXGK6UsFXsZoi6KDAvkxfsu4bLzW5BvczB9wZ/M+m4rBTaHr6tWbgnuMTAJYEJ4i8cBTCnVGtgM3AlsBN4C/gHGAr8ppRp4UMwlwGfA5cB2YDrwORADTAFWKqX8y/MCRN1mNhl58OpzefDqczEZDSz5PY5nP1rDybQcX1fNY1prElOd9Y2SFpgQXlOeFti7QBQwRmt9ldb6Ka31IJyBrD3wkgdlnABuBRprrUe5yrgPaAf8CVwEPFyuVyDqhcvOb8FL9/cjMjSA3YeSmThzJTv+SfJ1tTySmVNAbr6NAD+TPAtMCC/yKIAppVoBQ4ADwMwShycBWcBtSqkzTnLRWm/VWv9Xa51fYn8GMNX14wBP6iTqn7ZNw3nj4YF0a92Q1Mw8Jn28hi+X76rxqfZFuw+VUj6ujRB1h6ctsEGu7VKtdbEBCFfwWQcEAhdWoi6FawjZKlGGqONCg/34vzsu4toB7dDA/OU7mTx7Lcnpub6uWpncCRwyB0wIr/I0gLV3bfeUcXyva9uuEnW5y7X9uRJliHrAaDRwy5DO/N8dfQkL9mP7P0k8Pn05f+6J93XVSiWrcAhRNTwNYKGubVoZxwv3h1WkEkqpR4ChwFZgdkXKEPVPj7ZRTH10EN3aNCQ9K58X56xn3s/ba1yWYoKswiFElfDWPLDCjv1yD0Yopa4B3saZ4HGt1rpuLkcuqkS41Z9n7+jLzUM6YTAovv11L0+9t4qDJ9J9XTW3U12IkoEohDd5GsAKW1ihZRwPKXGeR5RSV3N7I6cAABQXSURBVAHzgQRggNb6n/JcLwSAwaAYNaA9L9x7CdERQcQdT2PizJV8u2ZvjUjwOJ6UCUgKvRDe5mkA2+3aljXG1da1LWuM7DRKqeuAr4B4oL/WevdZLhHijDqe04A3Hx3IZee3wGZ3MO+n7Uz6aA3xyVk+q1Nyei6HEzKwmI20aFzW339CiIrwNICtdG2HKKWKXaOUsgJ9gRxggyeFKaVuxjmB+RjO4LX3LJcI4ZEAPzMPXn0uT9/eh7BgP2IPnGTctBUs+T0Ohw9aY3/vTwCgS8tIzCZjtd9fiLrMowCmtd4PLAVacPpE48lAEDBPa+3+U1cp1UEp1aFkWUqp0cCnwCGgn3QbiqrQq0Mj3h47mD5dmpCbb2PWd1v5vw/XcCShesfG/trnDGDd20ZV632FqA9M5Tj3IWA9ME0pNRjYCfQGBuLsOnymxPk7XVv3zE2l1ECcWYYGnK26O0uZ2JmqtX67HPUSolQhQX5MuOkCftt+jI8W/cXOgyd5fPpKrunfjmv6t8NirtoWkdaav/a6AlgbCWBCeJvHAUxrvV8p1Qt4HmfK+xXAcWAaMFlrnexBMedwqtV3VxnnHMSZlShEpSmluKhrDN3aNOTTn3ew7I8DfLliF6u2HOKOYV3o3blJla2OcfBEOqmZeUSE+NMsylol9xCiPitPCwyt9WGci/l6cu5pnwpa6znAnPLcUwhvCA6w8ODV59KvRzM+WvQXB0+k8/r/NtKlVSR3Du9GyypIsNi61zmxunubKFlCSogqIM8DE/VK55aRTHl4IPeN6E5wgJnt/yQxfvqK/2/vzKOsqM4E/vsaBhfCjuAgQgsKGARbIaIwKu4aNXI0TlxHXEYzLqDCGR3mBHRGjRltcOIKLiExCSTRJI4e4oqAisa0W2wVUEMjoh420Ub2ft/8ce+zX8q3VL2u915X8/3OeafOu7fu992vvrr3q7p16xa1c17l49WNsenZur2Jea+4x7sjBu8Zm1zDMJqJdAdmGG2Bdu2qOPHQAYwZ3pffzl/CU39ezktvr+Ll+lUcUbM3444YRL/enQsLysP/vfA+azdspvofuzBqaJ+Yam4YRiYWwIydlk67d+DiU4Zz2j/tyyMLlvJc3QoWvLGSBW+s5MB9e3Hy6IEcPKg3VVXRhv/WbNjE7xe6VyIvOnkY7SKWNwwjHBbAjJ2enl1354fjDmLcEYN47IX3WfD6R7z1wWre+mA1vbvtzuhhezF6WF8G9OlS8FnWis++5OZfLGbr9iYOHdqHAwbsUSYrDGPnQ1Qrv9ROWIYNP0jnzX+p0tUw2jgbN2/j2b80MO+Vv7F2Q/OXn3t370jNfr0Y0r8He/fqRK9uu9OhfTu2bNvBqjWNvPjXVTz/+gq2bGticL/u/Mf5h9K54y4VtMQwysc+fbq/u2PblqHl1GkBzDBy0JRSljSsY3H9Kl5+ZxUbGreGKnf48L5cfsbB7FLi98wMozVRiQBmQ4iGkYN2VcLQAT0ZOqAnF50ynGUr1/Pu8rW8v/JzPl23kbVfbGb7jhQd2lexZ49vMbhfd477TrWteWgYZcICmGGEoF2VsH//Huzfv0elq2IYhsfeAzMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FYADMMwzASiQUwwzAMI5FECmAi0ldEHhKRT0Rkq4g0iMgdItItopzuvlyDl/OJl9s3WvUNwzCMnZX2YXcUkYHAYqAX8BiwBDgEmAicKCJjVHVdCDk9vJxBwHxgLjAEuBA4WUQOU9W/RTWkbaGABP5TME01hUhV3n1aU1rS6hu/DZn/DcOISpQ7sHtwwWuCqo5T1etV9WhgBjAYuDmknFtwwWuGqh7j5YzDBcJeXs9Oy6xZM6mtnU5zp6eMH38h48dfmDdt5syZnHfeecycOTNSuUqlJa2+cdtQWzudWbPS5QzDKIZQAUxEBgDHAw3A3YHsacBXwPki0rGAnI7A+X7/aYHsu7z8E7y+nRClsXEjc+bM+TqI3X57LfX19dTX13P77bVZ01RTLFq0kKVLl7Fo0UJUU6HKVSotafWN24ba2unMmTOHxsaNNAc1wzCiIqqFG5CIXALcD8xS1cuy5D+FC3DHqupzeeQcCzwDPK2qJ2TJnwlcClyiqg8G84cNP0jnzX+pYH2TTXMHl+ass84CYO7cuXnTBg8exNKlyyKXq1Ra0uobpw1nn302kyZdiw0jGm2Fffp0f3fHti1Dy6kzbAC7DZgMTFbV2iz5dwFXAJer6r155FyBu9O6S1WvypI/GbgN+B9VvS6YP7zmIH30ifkF65t0VJWRI0Z8/b/utdcACqfV1TFy5Mjo5SqVlrT6xmyDiAUvo+0wuLp3fWrHtmHl1Bl2EkcXv/0iR346vWsp5XRo3479+kaa8Jg4VJVrrrmGxjUNX6fdffuNAAXTfnDasUWVq1Ra0uobpw331P4XM2bMsCBmtBm0afvW8itVLfgDZuEG6y/JkX+Lz7++gJwpfr+bcuRf6vPvy5Y/YsQIbcukUimdOHGiAjpx4kRNpVI6YcIE9cdEJ0yYkDWtqalJa2pqFNCamhptamoKVa5SaUmrb9w2BH1sGG0BoE5DxJM4f2ED2G2+EU7KkX+Xz/+3AnKu8PvdmSN/ss//Sbb8th7AVFWnTZv2dx1bKpXSUaNG6ahRo/KmTZ06VWtqanTq1KmRylUqLWn1jduGiRMn6rRp04o7SQyjFVKJAJaoSRwjR47Uurq6gvVNOqr6d0NLaR8VSkulUlRVVUUuV6m0pNU3bhts+NBoS4jIa6o6svCe8RH2Gdjzfnu8iFSpaiqdISKdgDHAZuCVAnJe8fuNEZFOqtqYIacKFwQz9e2UBDu2bB1dtrTMjjRKuUqlJa2+pbbBMIxohHoPTFU/BJ4GqnHDgJncCHQEfqGqX6UTRWSIiAwJyNkIPOz3vyEg50ov/ynd6VfiMAzDMAoReikp4HLcElA/FZFjgPeAUcBRwDLgPwP7v+e3wUvNKcBY4FoRqQFeBfYHTgNW880AaRiGYRjfIPRSUv4ubCQwGxe4JgEDgZ8Ch2mIdRC9nHXAYb7cvl7OKOBnwAivxzAMwzDyEuUODFVdiVt0N8y+OQf5VXU9bu3DiVH0G4ZhGEYa+x6YYRiGkUgsgBmGYRiJxAKYYRiGkUgsgBmGYRiJJNRKHK0FEVkDrKh0PQzDMIxv0F9V9yinwkQFMMMwDMNIY0OIhmEYRiKxAGYYhmEkEgtghmEYRiKxAGYYhmEkkooGMBEZLSLzRGS9iGwSkb+KyNUi0i6CjL1E5CoR+ZOINIjIVhFZJyLPiMjpOcqMFRHN87u1BTb1FZGHROQTX5cGEblDRLpFlNPdl0vb9ImX27fUuiPWs0U6RaSjiJwrIr8WkSUi8pWINIpInYhMEpEOOcrl81+hz/oUTRzHWEQWFKj/rjnKfVtEfisiq0Vki4gsFZEbRWS3+Cz8hs6W+rdQW0v/9g6UK7t/ReT7InKniLwgIl96Xb8sUlbk41Zu/8Zhr4j0EJFLROQPIvKBiGwWkS9E5EURuVjcZ7KCZaoL+HduaP2VmoUoIqcBjwJbgN8A64FTgcHAI6p6Zkg5twLXAcuBhcBnQH/gdGAXYIaqXhsoMxb3zbGFwIIsYl9U1WeLsGkgbsX+XsBjwBLgENyK/UuBMWEWPRaRHl7OIGA+8BdgCM0r9h8W/ORMXLqjEIdOETkR+BPO/88DHwDdcefCnl7+Maq6JVBOca9UzM4i9mNVfaBow3LXNS7/LgCOxH2KKBs3qeqOQJlRuHPhH4BHgJXA0bgFtl/CHaOt0a3KW884/FsNjM+RPQzXTt9R1QMC5Srh3zeBA4GNwMe4NvcrVT0vopzIx61C/m2xvSLyQ+Be4FNc+/0I6I3zaxdcH3+mZgQaf04sB94C/phFbL2qPhKqAuX+BLS3ozOuI94KjMxI3xXneAXOCinrdODILOn7A194WSMCeWN9+g0x2/WUl3tVIH26T78vpJyZfv/pgfQJPv3JUukut71ADXAu0CGQ3gl4zcuZlKWcAgvKfN7G5d8FrumF1tsOeNfr+F5GehWus1Pg+tZqbx75c7ycCa3Ev0cB++E+AZXuI35Z6uNWQf+22F5ckD0VqAqk74kLZgqcEcir9umzW2xDOU+QDAMu8gb8PMcBUWBhDHpmZesAKUEAAwZ4mcuzOLMT7irnK6BjATkdgU1+/06BvCovX4EBceuuhL0FdJzjdTyeJa+sHVyc9hI9gOVsExn1asCPqLQ2e3PI74EbfdkEdKu0f7PoL7ZDj3zcKuHfuOwtIHOKl3lnIL2amAJYpZ6BHe23T2bJW4Q7qUeLyC4t1LPdb3fkyN9XRK4UkSkicpGI7NcCXWmbnlbVVGaGqjbihgF2Bw4tIOcwYDfgJV8uU04K92VscFdPceuOQjl0FvJfV++3KSJyhYjEaV+Q2O0VkR+IyPUicq2InJTnfM/ZXtQNJS/DDZsPCKs7BKX273jcEP/vVPXzHPuU079xUcxxq4R/y0Gh9ttHRC7z/r1MRIZHVVCpADbYb5cFM9SN/S/HfausaIeJSGfgDFykfzrHbucCdwI3Aw8Cy0TkkbAPqAPktMnzvt8OKoGcuHRHoRw6L/LbbBc64MbvH8T57y7gZRF5U0SGtUBnLkph71zgx0AtMA/4SES+XybdhSi1zkv8dmaefcrp37hISvstKSLSHvgX/zdX+z0OuA/n3/uAt0TkeRHpF1ZPpQJYF7/9Ikd+Or1rMcJFRIAHcA8T71XV9wK7rAGuxz1E7gTsAZwEvIELeo9nmz1TgLhsKkZOSY9nDkrtwyuBE4E3gYey7DIdGIPzXSfgO7jnBQcC80Vkr2L05iFOex/DPTfoi7vbHoILZF2B34jISSXUHZaS6RSRI3E2v6Oqi3PsVm7/xkVS2m+puRU4AJinqk8F8jYB/w2MALr535G4SSBjgedEpGMYJUUHMD8tNMz02PQvyvTM9Nectcjq1QJnAi8A1wYzVfUdVf2Jqtar6kZVXauqT+IO3nJcwzm1SN25aKlNLZETl+4oFK1T3OsPd+BmlJ6hqtuD+6jqJFVd7H23UVXr1M1cfRToCUxuQd2LIbS9qjpDVZ9Q1VWqukVVl6rqFGASrk3eUirdMdISnZf6bc67r1bo37hISvstGhGZgDuXlwDnB/NVdbWqTlXV11V1g/8tAo4H/gzsS/Mdel5acgf2IW5KaNjfJxll01cUXchO58B+oRGR24BrcM/SvqsRpp6q6pfAr/3fIyKqjsumYuSU7HjmoSQ6RWQcbmhtNTBWA68LhOA+v43qv0KU4xg/gHteUCMincqsO0ip/NsdN8qxGXi4iHqVyr9xkZT2WxJE5Argf3GzKo9S1fVhy/rHR+nXI0L5t33kGjYrO6bYsriANhI3pvtaZoYfO90H15AjdV4iMgO4Gncreoqqbiqibmv8NtQtbAZL/TbXOHV6gkiuce6WyIlLdxRi1ykiZ+IuID4DjlbV9wsUyUax/itEyY+xqm4RkUbckEpHID2Jp03413MBbvLGz1V1QxH1KpV/4yIp7Td2RORqYAZQj3tvbXURYqL5t1TTMgtMr4x1Gj3uFvtumids7NaCuqXfTfn3iOUGUnj67CYKT6P/FtGn0ceiuxL2ZpQ5B3fRsiLTtiLqdZmv17yYz9mSH2Pcw3wFvgTaZ6RXYhp9Seyl+X2n0a3Jv1n0jKW4afSRj1sl/BuXvRnlr/Pl3wB6tqAeP/Zy7gm1fylPgjyV7IyLtKFfZMZNPR0C9AukC3B/+qQGdg2hf0zw5PLp5wEpX6/qIuyK+gLjEGBIFjnpF5lrA+mJfpE5j70XAE24O+7+IfQenK3jBIYDa73uc1qjvb5D2iuL7J4Z5/6sQF6+F11/Ryt5kTmXfzPyD/fl3m6N/g3oGkueDh23YsYQYGAMx60i/o3R3h/5snVA9xC6RhFYuMCnH417NzD0BU4ll5Iah5tVtAX3zGM98D38UlLAP2tG5TKXf1LVsRnp04AbcGPqdwDbsqh7U1X/mFGmAXdyLMYtobIrbpbTIbi7gH9V1dlF2BRcQuY9nLOOwt3+j9aMJWT8cjmoqgTkBJeSehW3skh6KanRqvphS3THQRz2ishRwLM4fzyEW0InyAZVvSOjzGzcCizz/f5bcY3rRFxncD9wmcZ8csdk73jcOP9C3HPk9UA/4Lu4ZyB1wHEaGF7LstTQR8AxlHcpqaLO54z8h3EXiRNU9c48emdTGf+OA8b5v3sCJ+Auql7waWtVdbLftxp3l7VCVasDciK3xQr5t8X2isgFuOW+mnCvJGV7TteQ2Z/6pdSG4l7o/9gnD6f5fbgfqepNoYwoZVQPEYnH4O6aPscFoLdxEzDa5blCWBBIn+3T8/1mB8pcBzyDaxybcUH0Q+BnwIEttGlvL+dTXDBdgXuo+Y0rk3T9csjp7sut8HI+xXXwfePQHaMPW2Qv7oXWQv5rCJQZB/wet27ilxnH53EyrmBbqb3D/Dn7NrAO97LnelyncRVZrkwzyn4bd0W+FtepL8Otp1j0kHkZz+duvq1tAroW0FkR/+IuhEOdhzSvJtGQQ1bktlhu/8ZhbwgZ2frsi4EncMOiG72tH+HWxD08ig0VuwMzDMMwjJZg3wMzDMMwEokFMMMwDCORWAAzDMMwEokFMMMwDCORWAAzDMMwEokFMMMwDCORWAAzDMMwEokFMMMwDCORWAAzDMMwEsn/A+xz1kOoqWIiAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m.plot()\n", "plt.tick_params(axis='both', which='major', labelsize=20)\n", "plt.tick_params(length=0, color=\"black\")\n", "plt.legend(fontsize=25)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Laplace approximation of the posterior predictive" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the approximation of the posterior, we can calculate the predictive distribution of the latent variable $f^*$ for new data ${x}^*$:\n", "\n", "\\begin{align}\n", "q(f^* \\mid x^*, \\mathbf{x}, \\mathbf{y}) & = \\int p(f^* \\mid {x}^*, \\mathbf{x}, {f}) q({f} \\mid \\mathbf{x}, \\mathbf{y})d{f} \\\\\n", "\\mathbb{E}[f^* \\mid x^*, \\mathbf{x}, \\mathbf{y}] & = \\kappa({x}^*, \\mathbf{x} )(\\mathbf{y} - \\sigma(\\bar{{f}})) \\\\\n", "\\mathbb{V}[f^* \\mid x^*, \\mathbf{x}, \\mathbf{y}] & = \\kappa({x}^*, {x}^*) - \\kappa({x}^*, \\mathbf{x} ) (\\mathbf{K} + \\mathbf{W}^{-1} )^{-1} \\kappa(\\mathbf{x}, {x}^*)\n", "\\end{align}\n", "\n", "The class of a novel data-point $y^*$ is then calculated as:\n", "\n", "\\begin{align}\n", "P(y^* = 1 \\mid x^*, \\mathbf{x}, \\mathbf{y}) & = \\int \\sigma(y^*) q(f^* \\mid x^*, \\mathbf{x}, \\mathbf{y}) d f^* \\\\\n", "& = \\int \\sigma(y^*) \\mathcal{N}(\\mathbb{E}[f^*], \\mathbb{V}[f^*] ) d f^*,\n", "\\end{align}\n", "\n", "which can be calculated using the Gaussian error function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Prediction is fairly simple using GPy. We just plugin some new data into the predict function of the GP." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "x_new = rnorm(size=(100, 1))\n", "y_hat = m.predict(x_new)[0].reshape((100, 1))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X18VdWd7/HPL88QQwhPggJiFAi0giJDaK1PEJW2XnW0rZ2LFb3a18vauUBHXrS+nFaoMl65jA04pePUEXW0ra2dam9n2pG0oFIfccSH8igxKiryFEOMEEKy7h/rbM7JSQ7ZJ5CcnHO+79frvFbO2vvsvfbJPvu319p7rW3OOUREJHvlpLoAIiKSWgoEIiJZToFARCTLKRCIiGQ5BQIRkSynQCAikuUUCEREspwCgYhIlstLdQGSMWvWLPeHP/wh1cUQEUkn1tUMaVUj2LNnT6qLICKScdIqEIiIyPGnQCAikuUUCEREspwCgYhIllMgEBHJcgoEIiJZToFARCTLKRCIiGS5tOpZnIzm5mb27dtHY2Mjra2tqS6OSDu5ubmUlJQwaNAgCgsLU12co1q6FJYtg/p6yMmBAQOgshLmz4eqqlSXTo4HS6dnFk+dOtWtX7++y/mam5t59913KSsrY8CAAeTn52PWZS9rkV7hnKOlpYX9+/dTX1/P6NGj+2wwWLoUfvAD//ehQ9H84cOhrAyWL1cwSAOZNcREWPv27aOsrIwhQ4ZQUFCgICB9iplRUFDAkCFDKCsrY9++fakuUkLLlvlaQGurT3NzwQz27IF+/aC6OtUllOMhIwNBY2MjAwYMSHUxRLo0YMAAGhsbU12MhBoaID8f2tra5x8+DMXFUFubmnLJ8ZWRgaC1tZX8/PxUF0OkS/n5+X36GlZpKbS0+NpArLw8aGqC8vLUlEuOr4wMBICagyQt9PX9dMECXxvIzfVpays4B0OGwIED/oKxpL+MDQQicuwWLoQf/tDfKZSb65uJBg+GKVNgzhx/jWDiRLj0UqipSXVppbsUCETkqBYuhF27fBNRczPs3u1rAg89BDt3wtChPp03T8EgXSkQSErU1dVhZlx33XXt8q+77jrMjLq6uh5Z79q1azEzFi1a1CPLzzQ1Nf5sP/6sv7ra3zVUUuKvH5SU6C6idKZAkMHMrN0rNzeXIUOGMGPGDB599NFUF69HJAowkryaGn+W39lZf22tv2solu4iSl8Z27NYom6//XYAWlpa2LJlC0888QRr1qzhlVde4Z577klx6dq76667+N73vsfJJ5/cI8ufNm0amzZtYsiQIT2y/EwSe9YP0bS62t8ttHNnNA90F1E6UyDIAvHNIH/84x+56KKLqK6uZu7cuYwZMyYl5erMiBEjGDFiRI8tv3///lRUVPTY8jNJba2vCcQKzvpXrPC1gyCvqUl3EaUzNQ1loZkzZ1JRUYFzjpdffhlo36SydetWrr76aoYNG0ZOTg5r16498tl9+/Zx6623MmHCBPr160dpaSkzZ87kqaee6nRdjY2N/N3f/R0jR46kqKiIiooK7rnnHtrieyhFHO0awUsvvcTVV1/NySefTGFhISNGjODiiy/ml7/8JeAD3qmnngrAQw891K5Z7MEHHwSOfo1g27ZtXHvttZx88skUFBRw0kknce2117Jt27YO8y5atAgzY+3atTz++ONMmzaN/v37M2jQIL7+9a/z/vvvJ/r600Z5uT/AxwrO+quq/PASw4dDXR28+y40Nvragi4Ypx/VCLJUMMZU/H3s27dvp7KyknHjxjF79mwOHDhwpJf2O++8wwUXXEBdXR3nnnsus2bNoqmpid/97nfMmjWL++67j29+85tHltXc3MzMmTN5+eWXmTx5MrNnz+bjjz/mjjvu4Omnn06qvD/96U/51re+RW5uLpdddhljx45l165drF+/npUrV/K1r32NCy64gI8//pjly5czefJkrrjiiiOfP/PMM4+6/JdffpmqqioaGxu57LLLmDhxIps3b+bRRx/lySef5I9//CNTp07t8LmVK1fy29/+lssuu4zzzz+fF198kccee4zXXnuNDRs29NkxhMKYP//oZ/3BGEPz5vmaQ3Fx9DqCxiBKM865tHmdffbZLoyNGzeGmi/TAc7/i9tbvXq1MzNnZq6urs4559zbb799ZP5bb7210+Wdf/75zszcz3/+83b59fX1bvLkya6oqMjt3LnzSP6SJUsc4K688krX2tp6JL+2ttaVlZU5wM2ZM6fdsubMmeMA9/bbbx/J+8tf/uLy8vJcWVmZe/PNNzuU67333jvyd7Ad8csNrFmzxgHu9ttvP5LX1tbmKioqHOAeeeSRdvP/4he/cIAbP358u224/fbbHeBKSkrc66+/3u4zf/M3f+MA99hjj3Vahnh9eX9dvdq5L3/ZuQkTfLp6dfvpX/6yc2ef7dykSc4NGuRc//7OlZY6N21aasornery2JpU05CZjTSzB8zsAzNrNrM6M6s2s7Ikl/MFM3sy8vmDZvaumf2nmc1KZjmplujWur5m0aJFLFq0iNtuu42vfOUrzJo1C+cc8+fP55RTTmk374knnnjk4nKs1157jaeffpqrrrqKr3/96+2mDRw4kMWLF3Pw4EF+/etfH8lftWoVOTk5LF26lJyYMQpOPfVU5s6dG7r8P/nJTzh8+DDf//73+cxnPtNh+siRI0MvqzPPPfccmzdv5nOf+xyzZ89uN+3qq6/mC1/4Alu2bGHdunUdPjt37lzOOOOMdnlBreill146pnL1BVVV8LvfwcaNPo0/y6+t9X0L3nrLj06an+97H2/Y0Hd/D9JR6KYhMzsNeA4YBjwJbAamAfOAWWZ2jnNub4jlfAtYCTQBvwF2ACOBK4EvmtnfO+eWJLshvS24ta5fv/a31vXFKvHixYsB3ww0cOBAzj33XG644QauueaaDvNOnjy50+aM559/HoCGhoZO29d3794NwKZNmwB/beCtt95i1KhRnHbaaR3mv+CCC46UqysvvPACAF/84hdDzZ+s//7v/wZgxowZnU6fMWMG69at49VXX+W8885rN62z5qJRo0YBUF9ff5xL2veUl8O6dX74iUOHfGoGhYX+ekFf+y1I55K5RrASHwTmOufuDTLN7B7gO8AS4KajLcDM8oG7gIPA2c65LTHT/gF4FbjNzJY555qTKFuvO9qtdX1t53dJPHNi+PDhnebv3etj/OrVq1m9enXCz3/yySeADxjgaxjJrKczH3/8MUCP3VIalDXR3UpBflCOWAMHDuyQl5fnf1Z9eTC542X+fPiv//K1ADP/cs6PTvrGG6kunYQVqmnIzMqBi4E64Mdxk2/Hn91/w8ziuph0MAgoBbbGBgEA59wmYCvQDzghTLlSKVM71CQaBK20tBSA5cuXH7WtcdWqVe3m/+ijjzpd3s6dO0OXKTjY9tSdOEFZE5Xpww8/bDefRFVV+f0+2G3M/AlSbq6/sCzpIew1gqDO/JRzrt19f865RuDPQH9gehfL2QXsBsaZ2djYCWY2DhgLbAjTxJRqR7u1LhNNn+7/tc8++2yo+UtKSjj99NN5//332b59e4fpsbekhl3373//+y7nzc3NBZI7Gz/rrLOOWqYgf8qUKaGXmU0GDICiIujfPxoUzHyepIewgWB8JN2aYHpwo/W4oy3E+TaKb0fW+4qZPWRmd5nZw8ArwF+Ar4YsU0rNn+/PeBobfbtoY2Nmd6iZOnUq5557Lv/+7//OAw880Ok8b7zxBrt27Try/vrrr6etrY3vfve77foNvP3226xYsSL0ur/1rW+Rl5fHHXfcwcaNGztM37Fjx5G/y8rKMDPefffd0Ms/55xzGD9+POvWrePxxx9vN+3xxx/nmWeeYdy4cXzhC18IvcxsMmkSnHQSFBT4gekKCvz7SZNSXTIJK+w1gqBO3JBgepDfscE0jnPuV2b2AfBz4NqYSR8Bq4C0aFwJOtRUV/vmoPLyzH+Y989+9jNmzJjBDTfcwIoVK6isrGTgwIHs2LGD119/nTfffJPnn3+eYcOGAXDLLbfwxBNP8Otf/5opU6ZwySWX0NDQwGOPPcZ5553Hb3/721DrnThxIitXruSmm27irLPO4vLLL2fs2LHs3buX9evXU1JSwpo1awA44YQTqKys5Nlnn2X27NmMGzfuSN+DSQmOTGbGQw89xEUXXcTVV1/N5ZdfTkVFxZHhOEpKSnj44Yfb3fkkUUF/g1NPVS/jtBXmHlPgX/D3mN+YYPo/RKZ/L8SyrgEOAI8CFfhrAhWR9w74ZaLPqh9BckjQj6AzXd1/H9i/f79bsmSJmzJliisuLnZFRUVuzJgx7ktf+pK777773CeffNJu/oaGBved73zHnXTSSa6wsNCNHz/eLVu2zG3fvj10P4LAc88956688ko3dOhQl5+f70aMGOEuueQS96tf/ardfNu2bXOXXnqpGzRokDMzB7hVq1Y55zrvRxDYvHmzu+aaa9zw4cNdXl6eGz58uJs9e7bbvHlzh3mDfgRr1qzpMC3sdxnIhP21q/4GklJdHuPNhbijxMz+L7AAWOCc+8dOpv8TvsnnZufcT46ynHHAm8DrwDQXc73BzHKAl4CzgQudc2vjPz916lS3fv36Lsu7adMmJkyY0OV8In2B9lfpYV0+Bi9sXTe4wyfRNYDgwm+iawiBi4F84GnX8aJzG/BM5O3ZIcslIiLHKGwgWBNJL46cuR9hZiXAOfjmnhe6WE7QU2logulB/qGQ5RIRkWMUKhA457YDTwFj8E1AsRYDxcDDzrkjN1SaWYWZxY/3G9x7+BUza3flzszOBL6Cb9f+U9gNEBGRY5NMz+Kb8UNMrDCzmcAmoBK4EN8kdFvc/Jsi6ZH2KefcS2a2CrgeeNnMfgO8gw8wVwAFQLVz7i/Jb4qIiHRH6EDgnNtuZlOBHwKzgC8BHwIrgMXOuX0hF3UD/lrAdcAlQAmwH1gH/NQ594vQpRcRkWOW1PMInHPv4c/mw8zb6ZVq529TejDyEhGRFFMPGRGRLJexgSBM/wiRVNN+Kn1BRgaC3NxcWlpaUl0MkS61tLQcGShPJFUyMhCUlJSwf//+VBdDpEv79++nJHiYhUiKZGQgGDRoEPX19ezZs4dDhw6p+i19inOOQ4cOsWfPHurr6xk0aFCqiyRZLqm7htJFYWEho0ePZt++fdTV1WXFk6IkveTm5lJSUsLo0aM7fTSoSG/KyEAAPhiMGDEi4eMHRUTEy8imIRERCU+BQEQkyykQiIhkOQUCEZEsp0AgIpLlFAhEpFfU1MCll8LEiT6tqUl1iSSgQCAiPa6mBubNg507YehQn86bp2DQVygQiEiPq66Gfv2gpARycnzar5/Pl9RTIBCRHldbC8XF7fOKi32+pJ4CgYj0uPJyaGpqn9fU5PMl9RQIRKTHzZ8PBw5AYyO0tfn0wAGfL6mnQCAiPa6qCpYvh+HDYfduny5f7vMl9TJ20DkR6VuqqnTg76tUIxARyXIKBCIiWU6BQEQkyykQiIhkOQUCEZEsp0AgIpLlFAhERLKcAoGIpIyGpu4bFAhEJCU0NHXfoUAgIilRXQ0tLfD227B+vU9bWjQ0dSpoiAkRSYk33vDjDuXmQn4+HDoEO3ZAc3OqS5Z9VCMQkZQ4cMCneXlg5tPYfOk9qhGISEoUFflnEhw8CIcPQ2urDwjxD7CRnqcagYikxKRJUFrqm4RaW/0jLPPzfWDQBePepUAgIikxf75/QE1hIQwY4J9hnJsLI0bognFvU9OQiKREVRUMGuSDQXOzbyoaNcrXEvQs496lQCAiKXPGGb7/QElJNK+xUc8y7m1JNQ2Z2Ugze8DMPjCzZjOrM7NqMytLdsVmdoaZPWxm70WWtcvMnjaza5NdloikJz3LuG8IHQjM7DTgFeB64CXgR0AtMA943swGJ7Gs64BXgSuAZ4F/BB4HDPhS2OWISHrTs4z7hmSahlYCw4C5zrl7g0wzuwf4DrAEuKmrhZjZdOB+4E1glnNuZ9z0/CTKJCJpTs8yTr1QNQIzKwcuBuqAH8dNvh1oAr5hZmHuAF4K5ALXxAcBAOdcS5gyiYjI8RG2RjAjkj7lnGuLneCcazSzP+MDxXTgj4kWYmYjgXOB9cBfzOxC4GzAARuANfHLFxGRnhU2EIyPpFsTTN+GDwTjOEogAP4qZv4/ARfETX/DzK50zr0VslwiInKMwl4sLo2kDQmmB/kDu1jOsEj6NWACcGVk2acD/wacAfyHmRWELJeIiByj49Wz2CKp62K+3Jj0Rufcb5xz+51z24E5+CajccBVx6lcIiLShbCBIDjjL00wfUDcfInUR9Jm4D9jJzjnHPBk5O20kOUSEZFjFDYQbImk4xJMHxtJE11DiF9OY4KLwkGg6BeyXCIicozCBoI1kfRiM2v3GTMrAc4BDgAvdLGc14E9wBAzO7GT6Z+NpHUhyyUiIscoVCCItOE/BYwBvh03eTFQDDzsnGsKMs2swswq4pZzGLgv8nZpbFAxszOA64DD+F7GIiLSC5LpWXwz8BywwsxmApuASuBCfJPQbXHzb4qkFpf/D8BM4FrgDDNbCwzFXyAuAm7R7aMiIr0n9F1DkVrBVOBBfAC4BTgNWAF8zjm3N+RyPsUHgsVAf3wN4zJ8kPmSc+6eJMovIiLHyPzNOulh6tSpbv369akuhohIOolvlelATygTEclyCgQiIllOgUBEJMspEIiIZDkFAhFJOzU1cOmlMHGiT2tqUl2i9KZAICJppaYG5s3zD70fOtSn8+YpGBwLBQIRSSvV1dCvH5SUQE6OT/v18/nSPQoEIpJWamuhOO6huMXFPl+6R4FARNJKeTk0NbXPa2ry+dI9CgQiklbmz4cDB6CxEdrafHrggM+X7lEgEJG0UlUFy5fD8OGwe7dPly/3+QHdVZQcjTUkIhkluKvo0CHYtw8+/RTy8uD734eFC1NdupTQWEMikl2qq30Q+OADn+bk+GsI3/2uv7to+nTVEOIl8zwCEZE+r7bW1wRycsA5f/0gaPhoaYHXXoMbb4Sbb4ZnnvHzl5fDeee1fz9/fvvmpkyW8TWCsG2FybQpdjZvTQ1UVvp7mktKju2sY+lSGDYMCguhrAwqKjqWK7YM06f7dY8e7T93yimJt6Gr7Yxd97Bh/n38tk2Y4N93Vqbp0/1ZV16eTysrk/vO4/OWLu34vrKy/To6+66XLoWiIjDzr9xc+OY3O363o0Z1/J91tr0VFdF1FhTAwIFH/54ldcrLfXNQbi40N0eDAPi/8/L8tYU77oh2Stu6FX7wA9iy5dg6qcXuO/36+f1k9GifN3164t9dsN+fcorfP4PPVFQcn2NKVzL6GkHQVtivn7/PuKnJnx10dmEpzHyJ5t2zx+94jY1+J3MOWlvhxBPh/vuTO6tYutTvkDk5/gB24IDPHz0ahgzx7+fMgYce8mVoafE7cWtr9IAHMHIk5Oe334autjN23fn5ftktLX7e5ma/bYcPw8GDfvr48f6AGpRp5UpfHT98OLo9eXkwYgT86792/Z3v3u23YcgQn/fhh/Duu/5gfdJJ/v3bb/s7RWJ327w8Pz34rpcuhVtv9fPFy8nxZTbz/zPw21JQ4L/DAQN83v797be3MwUF/of78cd+WZ9+6pd/+unwT/+UPWeTfU1NDfz1X/v/Z3NzdD8Ifh/Fxf63WlDgD64Ar7/u9+N+/WDSJJ/X2OgvRP/ud+HXe8MNsGuXf9/c7NPc3Oi6x43z+1vs7y74LbS0wI4d/jNtbf69c35/zc3t/jGFbL9GELYHYjI9FTubt77ev/Ly/Cs/36f19cn3dly2LHqwCto3zfwBNijXsmXRMuzY4dfV1uZ3lGCn2bu34zZ0tZ2x6w7StjZoaIhu2+HDflpbW8cy1df7/JwcX4ZgvoaGcN95Q4NfRpC3d69Pg2r+3r1+eZ2tI/a7Xras8yAAPj/4bgOHD0f/Z3v3+lfs9saymJ9Uayu8954PYMF97c75s8qLLvLLDGpV0nuqqvyF4eCkCKI1w8JC/38D6N8/Oj04uYkN+sl2Uquujv5Wgt9JTo7/u63N5+/Y0fF3F/wW9u71ZS4sjH7GrP3+2Z1jShgZHQjC9kBMpqdiZ/O2tPidK3bHy831+cn2dmxo8P90aH8wCw5IxcV+nqAMBw/6dcWeJefm+vz4behqO2PXHYhdZlAmM58GP5qgTMEZTPznO/seEn2PLS3R9/E/zkRn5vHraGjofL5Ysd9t7DYGgSZ2e+PXFft3bJli///B9E8+8bUsBYPetXAh/OY38Fd/5Q+gZr4GAP63VFgIgwZF5y8q8vtQUVE0L9lOarW1fhnBfhQr2KdifzPB/hr8FoLfcjB/7Geh+8eUMDI6EITtgZhMT8XO5s3Pj1bdAq2tPj/Z3o6lpdGDYU7MfycvL1qu0tJoGYqK/LqCmkOw7qKijtvQ1XbGrjsQu8ygTM75NPjRBGXKz29/xhx8vrPvIdH3GBuI4n+csT/So62jtLTz+WLFfrex2xicxcVub/y6jvZ3fIAJalfLlnVdJjm+qqrghRfg97+HqVOjAWHyZFi0yAeGoFPa4ME+HTSo+53Uysv9fhjsR7GCfSr2NxPsr8FvIfgtB/PHfha6f0wJI6MDQdgeiMn0VOxs3rIy/zp82L9aWnxaVpZ8b8cFC/xym5v9jhqc6Z90UrRcCxZEyzByZLQaGlwca231O3b8NnS1nbHrDtKcHH9gDbYtaIbKyelYprKyaFNNa2t0vtLScN95aalfRqIf5+DB0QN1/Dpiv+sFCzr+EAM5OdHvNpCXF/2fDR7sX7HbGyu2RhA0TUHHABisC/yPN0wtRXpGVRW8+KLfrxobfXBYuLB9p7Rx4+CHP/TXvRJ1UuvK/PnR30rwOwmahIImopEjO/7ugt/C4MHR6xrBZ4KL28dyTAnFOZc2r7PPPtsla/Vq5778ZecmTPDp6tXHNl+ieVevdm7aNOdOOMG/KiuPvoyjuftu54YOda6gwLmBA50bP75juWLLUFnp1z1qlP/c6NGJt6Gr7Yxd99Ch/n38tlVU+Pedlamy0rmiIudyc306bVpy33l83t13d3w/bVr7dXT2Xd99t3OFhc75n5JzOTnO3Xhjx+925MiO/7POtnf8+Og68/OdKy313/P48T7PzL+C9YFzxcV+Hf36+XWG+f4lvcXuO0VFfj8ZNcrnVVYm/r8H+8Xo0X5fCT4zfvxxOaZ0eWzN6LuGRHrD0qVw113Rs/6gSaiw0J/JtbX5s80pU8LfnSZyHGX3XUMivWHhwugdU21tcPfd/jbUgwfhhBN8EFi4UOPoS9+lGoFIL5k40XdWir1+0dbm26Q3bkxduSTjqUYg0ldoHH3pqxQIRHqJxtGXvkqBQKSXdDaO/pw5/hqBxs2XVNI1ApEUSWaMK5FjoGsEIn2V7iKSvkKBQCRFYsdb2rfPj4C5aRM8/bSaiKR3KRCIpEhwF9G+ffDWW35EVDNfO+jOWPgi3aVAIJIiwV1E77wTHZTs4EE/1syOHfD3f5/qEkq2UCAQSZHgLqLWVl8bOHTID4YXjEL56quqFUjvUCAQSaGqKjj/fH/w79cv+vQ0M104lt6jQCCSYvPnRx9Y4lz06VSjRvXMQ0hE4ikQiKRYVRWceWb0CVQFBf65x4WFGn5CekdSgcDMRprZA2b2gZk1m1mdmVWbWVl3C2Bm55lZq5k5M7uzu8sRSWdLlsDJJ8OECfDZz0YfcK7hJ6Q3hA4EZnYa8ApwPfAS8COgFpgHPG9mg5NduZmVAA8Bnyb7WZFM0tnwE+phLL0lr+tZjlgJDAPmOufuDTLN7B7gO8AS4KYk178cKAXuinxeJGtVVenAL6kRqkZgZuXAxUAd8OO4ybcDTcA3zKw47IrN7HJ87WIu8EHYz4lki5oaPxCdBqSTnha2aWhGJH3KOdcWO8E51wj8GegPTA+zMDMbBvwUeMI590jIMohkjWBAup07/cNsdu7s3d7GCkLZJWwgGB9JtyaYvi2Sjgu5vH+JrDvZpiSRrBAMSNfSAm++6ccgev99uO22nl93qoOQ9L6wgaA0kjYkmB7kD+xqQWb2v4DLgZudcx+FXL9IVqmt9UNNBGMQ5ef73sYbNvT8AVmjomaf49WPIBjv+qgPNzCzMUA18Cvn3C+P07pFMk55Obz3nj8Qt7TA/v3wySf+77/9255dd+yoqIHiYnVuy2RhA0Fwxl+aYPqAuPkSeQA4ANwccr0iWSkYkK652b9inx+1ZQssXdpz69azlbNP2ECwJZImugYwNpImuoYQmIK/BXV3pAOZMzMHrIpMvy2S90TIcolkpKoqOOssP9xELOf8OETLlvXcuvVs5ewTth/Bmkh6sZnlxN45FOkUdg7+TP+FLpbzMP7uonhjgfOADfhOa6+GLJdIxrrzTrjooo75zkF9fc+tN+jcVl3tm4PKy30QUB+HzBUqEDjntpvZU/i+BN8G7o2ZvBgoBu5zzh2pUJpZReSzm2OWM7ez5ZvZdfhA8B/OOY3CLoI/8Obm+ovEEB2V1Dl/7aCn160Df/ZIZne6GdgFrDCzJ8zsLjP7E75X8VYg/sa2TZGXiHRTWVn0qWXgm2qc8y/dzinHS+hA4JzbDkwFHgQqgVuA04AVwOecc3t7ooAi2ayy0o87FNQEAm1tcMMNCgZyfJhzR73js0+ZOnWqW79+faqLIdJrgs5ddXX+gm0QEAoL/fRJk+DFF1NaROn7rKsZ9DwCkT4suHDb3OzfB81Ehw75PgWvv57a8klmUCAQ6eOqqqKPsQwq8Ga+eai5Wc1DcuwUCETSwGc+42sB8TT0gxwPCgQiaeDOO32TUFATAP9Iy/JyDf0gxy6ZB9OISIpUVcGUKX54iZYWKCqCkSP9YHTDh6e6dJLuVCMQSRPBc41POsm/37rVB4bzzkttuST9KRCIpImqKpgzxz8f4NNP/fWBESPgoYd0wViOjZqGRNLIM8/A+PH+GQGBxkZ/wVhDQkh3qUYgkkZinxWwb5/vR7BpEzz9tGoF0n0KBCJpJHhWwL590aeXBZ3M9DhJ6S4FApE0Ejwr4J13ogPROQennqo+BdJ9CgQiaSQYcqK11b8KCmAOzRLTAAAKZklEQVTsWD9KqR4nKd2li8UiaaaqCs4/3989VF8PGzf6J5nl5sLpp6e6dJKOVCMQSUPz5/uH29fW+iBg5msIb73Vs88zlsykQCCShqqq4ODB6LATOTnQv7/vadyTzzOWzKSmIZE0dfAgnHBC+8dWtrVBQ0PqyiTpSTUCkTRVWurHHYrV0uLzRZKhQCCSphYsiD6TIDZdsCDVJZN0o6YhkTS1cKFPly3zzUGlpT4IBPkiYemZxSIimU3PLBbJBjU1cOmlMHGiTzXUhCRDgUAkzdXU+HGGdu6EoUN9qnGHJBkKBCJprrrajzPU0gJvvulHI33/fbjttlSXTNKFAoFImqut9XcMBaOR5uf7XsYbNqhWIOEoEIikufJyP9xETg7k5fmexuCfa6zRSCUMBQKRNBcMTe2cfx0+7PsTjBmj0UglHAUCkTRXVQVnneVHH21piQ5NnZ/vawsiXVEgEMkAd94JI0fChAnw2c/6JqIDB3xtQaQrCgQiGSB4YM3w4bB7t0+XL9cD7SUcDTEhkiGqqnTgl+5RjUBEJMspEIhkIA05IclQIBDJMBpyQpKlQCCSYYIhJ0pKfCezkhL/Xp3LJBEFApEMU1sLxcXt84qL1blMElMgEMkw5eXQ1NQ+r6lJncsksaQCgZmNNLMHzOwDM2s2szozqzazspCfLzaz2Wb2MzPbbGZNZtZoZuvN7BYzK+jeZohIIBhyorHRDzXR2KjOZXJ0oQOBmZ0GvAJcD7wE/AioBeYBz5vZ4BCLORd4BLgEeBO4F/g5cDKwDFhjZkXJbICItKfOZZKsZDqUrQSGAXOdc/cGmWZ2D/AdYAlwUxfL2AlcA/zKOXcoZhklwFrg88C3gX9MolwiEkedyyQZoWoEZlYOXAzUAT+Om3w70AR8w8ziLlG155zb4Jx7NDYIRPIbiR78LwhTJhEROT7CNg3NiKRPOefaYidEDuJ/BvoD04+hLC2R9PAxLENERJIUNhCMj6RbE0zfFknHHUNZ/lck/cMxLENERJIUNhCURtKGBNOD/IHdKYSZ/S0wC9gAPNCdZYhI5zTchHTlePUjiDwcD5f0B82uBKrxF5Kvcs61dPEREQlJw01IGGEDQXDGX5pg+oC4+UIxsyuAXwC7gAucc+r7KHIcabgJCSNsINgSSRNdAxgbSRNdQ+jAzL4K/Ar4CDjfObeli4+ISJI03ISEETYQrImkF5tZu89E+gCcAxwAXgizMDP7n/iOZB/gg8C2Lj4iIt0QO9xEfT28/jq89BLs2aPmIYkKFQicc9uBp4Ax+A5fsRYDxcDDzrkjI5yYWYWZVcQvy8zmAP8GvAucp+YgkZ4TDDfx/vuwcSN8/DE0N/vgcMMNCgbimXPhru9Ghph4Dt+7+ElgE1AJXIhvEvq8c25vzPwOwDlnMXkXAjX4APQA8F4nq/rYOddpC+bUqVPd+vXrQ5VXRLyaGrjqKj/mUE4OFBb69PBhmDQJXnwx1SWUHmZdzhA2EACY2Sjgh/hbPQcDHwJPAIudc/vi5u0sEFwHrOpiNe8458Z0NkGBQKR7Skp8mhczqMzhSNfNxsbeL4/0qi4DQVIPr3fOvYcfdC7MvB1W7px7EHgwmXWKyPERf87nHFiXhwjJBnoegUgW+MxnoLXV1wKc82lrq88XUSAQyQJ33gknnuivDRw65NMTT/T5IgoEIlmgqgruvx/OPRfGjPHp/fdrqGrxkrpGICLpS88okERUIxARyXIKBCICaJTSbKZAICIapTTLKRCIZLmlS2HWLD8ExSuv+LGIWlo0Smk2USAQyWJLl8Jtt/k+BYGDB31QaG7WKKXZQoFAJIstWwZtMU8hD3oaHz4M773nRy+VzKdAIJLFGhp8T+OcyJEgdhiKAwf86KWS+RQIRLJYaamvBZj5YBDUCMzgrLPU7yBbKBCIZLEFC3wAiG8eGjpUw09kEwUCkSy2cCEsWQIDBkSbiMaNg0cfVW0gmyT1PIJU0/MIRESS1uVg46oRiIhkOQUCEZEsp0AgIpLlFAhERLKcAoGISJZTIBARyXIKBCIiWU6BQEQky6VVhzIz2w28k+pyiIikkT3OuVlHmyGtAoGIiBx/ahoSEclyCgQiIllOgUBEJMspEPQCMxtrZt81sz+Z2XtmdsjMPjKzJ83swlSXL5XMLN/M5pnZKjPbEPlunJndmOqy9RYzG2lmD5jZB2bWbGZ1ZlZtZmWpLlsqmdlXzOxeM3vWzPZH9otHUl2uvsDMBpvZjWb2GzN7y8wOmFmDma0zsxvMLKljuy4W9wIz+wVwNbARWAfsA8YDlwG5wDzn3IrUlTB1zGwgUB95+xFwCBgFfNM5d3/KCtZLzOw04DlgGPAksBmYBlwIbAHOcc7tTV0JU8fMNgCTgU+AHUAF8Khz7pqUFqwPMLObgJ8AHwJrgHeBE4ErgVLg18BXXdgDvHNOrx5+AdcBZ3WSfz7+wNcMjEh1OVP03RQAXwy2H1gEOODGVJetl7b/vyLb+7/j8u+J5P9zqsuYwu/mQmAsfjz9CyLfxyOpLldfeAEzgP8B5MTlD48EBQdcFXZ5ahrqBc65B51zr3aS/zSwFn8w/Hxvl6svcM4dcs793jn3YarL0tvMrBy4GKgDfhw3+XagCfiGmRX3ctH6BOfcGufcNhc5wkmUc+5Pzrn/55xri8vfCfxz5O0FYZenQJB6LZH0cEpLIakwI5I+1ckPuhH4M9AfmN7bBZO0lvQxRYEghczsFGAm8CnwTIqLI71vfCTdmmD6tkg6rhfKIhnAzPKAayNv/xD2c3k9UxzpipkVAo8ChcBC51x9Fx+RzFMaSRsSTA/yB/ZCWSQz/B/gs8B/Ouf+K+yHVCMIKXJLn0vilfA2NzPLBf4NOAd4DFjWW9vRE47ndyPtBA8dVxu5dMnM5gK34O88+0Yyn1WNILztwMEk5v+gs8xIEHgE+CrwS+CaDLgYdly+mywUnPGXJpg+IG4+kU6Z2beB5fhb1Gc65/Yl83kFgpCcczOPdRmR9ruf4YPAz4BrnXOtx7rcVDse302W2hJJE10DGBtJE11DEMHM5gM/At7EB4FdyS5DTUO9xMwKgMfxQeBh4BuZEATkmKyJpBfH9wQ1sxJ80+EB4IXeLpikBzP7Lj4IbAAu7E4QAAWCXhG5MPwb4HLgX4Hr428XlOzjnNsOPAWMAb4dN3kxUAw87Jxr6uWiSRows+/jLw6/gq8J7On2stK/ebrvM7NV+N7Fe4CVdH7xb61zbm0vFqvPMLPv4YcPADgTP6zAc0Rvn1znMnS4iU6GmNgEVOJ71W4FPu+yd4iJK4ArIm+HA5cAtcCzkbw9zrkFqShbqpnZHOBBoBW4l86vI9U55x4MszxdI+gdp0bSIcAPjjLf2p4vSp80Cz/cRqzP0763dUYGAufcdjObCvwQ/z18CT9+zApgcbIX/TLMmcCcuLzyyAv80wqzMhAQPabkAvMTzPM0Plh0STUCEZEsp2sEIiJZToFARCTLKRCIiGQ5BQIRkSynQCAikuUUCEREspwCgYhIllMgEBHJcgoEIiJZToFARCTL/X94uODgSQ4ujAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "_, ax = plt.subplots()\n", "ax.scatter(x_new, y_hat, color=\"blue\", alpha=0.75, label=\"Prediction\")\n", "ax.spines['top'].set_visible(False)\n", "ax.spines['right'].set_visible(False)\n", "ax.tick_params(axis='both', which='major', labelsize=20)\n", "ax.tick_params(length=0, color=\"black\")\n", "plt.legend(fontsize=20)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Metropolis-hastings sampling" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }