{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Regression with Outliers\n", "\n", "In the standard __Gaussian process regression__ setting it is assumed that the observations are __Normally distributed__ about the latent function. In the package this can applied using either the `GP` or `GPE` functions with which *exact Gaussian process* models.\n", "\n", "One of the drawbacks of exact GP regression is that by assuming Normal noise the GP is __not robust to outliers__. In this setting, it is more appropriate to assume that the distribution of the noise is heavy tailed. For example, with a __Student-t distribution__,\n", "$$\n", "\\mathbf{y} \\ | \\ \\mathbf{f},\\nu,\\sigma \\sim \\prod_{i=1}^n \\frac{\\Gamma(\\nu+1)/2}{\\Gamma(\\nu/2)\\sqrt{\\nu\\pi}\\sigma}\\left(1+\\frac{(y_i-f_i)^2}{\\nu\\sigma^2}\\right)^{-(\\nu+1)/2}\n", "$$\n", "\n", "Moving away from the Gaussian likelihood function (i.e. Normally distributed noise) and using the Student-t likelihood means that we can no longer analytically calculate the GP marginal likelihood. We can take a Bayesian perspective and sample from the joint distribution of the latent function and model parameters." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAHU5JREFUeJzt3W9sXfV9+PHPiSLdEIJBTteONmRtVPLPInYwk6JVLKDFmRaQiO/VHoFMpqIgpCpEqxRXYiwhBKuTULSIPSiqpqxRJx5MPuKB24n8EaR7MCGlQOU6xCGIQKIRVS0uUYBYYZzfAw//6hI7Dv7ec33t10uK1Ovc3Hz6aeS9973H52ZFURQBAEAyCxo9AADAXHPNwBodHY3vfe97cfvtt8cdd9wRDz30UBlzAQA0rYXXesIPfvCDyLIsTp8+HVmWxYULF77wnI8//jhOnToVq1evjsWLF9dlUACAZpFNdQ3WRx99FLfeemucP38+WlpaJn2R1157LTo7O+OXv/xl3HnnnV9qkA8//DBuvvnmL/VnuT52XR67Lo9dl8euy2PX5Um96ynfInz77bejtbU1+vr64q677oq77747jh07NunzL126FBcvXhz/NTo6Ou1B/vd//3f6UzMjdl0euy6PXZfHrstj1+VJvesp3yL89NNP49133421a9fGD3/4w3j99dejq6srhoaG4mtf+9oXnr9x48YJj3ft2hW9vb3TGmRkZOQ6xmYm7Lo8dl0euy6PXZfHrssz0123trZOeDxlYC1fvjwWLFgQDz74YERErF+/Pr71rW/F4ODgVQPr+PHj0dHRMf64UqlEpVL50sNRP3ZdHrsuj12Xx67LY9flSbnrKd8i/MpXvhJ/9Vd/FS+99FJERLzzzjvxzjvvxJo1a676/CVLlkRLS8v4r+uJKwCAueKaP0X4ox/9KL773e9Gb29vLFiwIJ5//vn4xje+UcZsAABN6ZqBtWLFinj55ZfLmAUAYE5wJ3cAgMQEFgBNLc/zaF/fGTcsvjHa13dGnueNHgkEFgDNK8/zqNVqMXhlaVy+f08MXlkatVpNZNFwAguApvXU089E1tYVxY6BiK6dUewYiGztpti7r6/RozHPCSwAmtbp4VNRrO2KyLKxL2RZFG2bY/jUm40djHlPYAHQtFauWh3ZySMRn3+sblFENnQ4Vk1yv0YoyzVv0wAAs9XuJ5+IWq0W2YH7omjbHNnQ4ShOHo3drsGiwZxgAdC0qtVq9Pf3x7rKSCwa2BPrKiOR53l0d3c3erSm56czZ0ZgAdDUqtVqvPHaifjk44/ijddOiKsE/HTmzAksAGACP505cwILAJjAT2fOnMACACbw05kz56cIAYAJ/HTmzDnBAgAm8NOZM+cECwD4gmq1GtVqtdFjNC0nWAAAiQksAIDEBBYAQGICCwAgMYEFAJCYwAIASExgAQAkJrAAABITWAAAiQksAIDEBBYAQGICCwAgMYEFAJCYwAIASExgAQAkJrAAABITWADAvJXnebSv74yvL7st2td3Rp7nSV5XYAEA81Ke51Gr1WLwytIYvX9PDF5ZGrVaLUlkCSwAYF566ulnImvrimLHQETXzih2DES2dlPs3dc349cWWADAvHR6+FQUa7sismzsC1kWRdvmGD715oxfW2ABALPC59dD3bD4xqTXQ01m5arVkZ08ElEUY18oisiGDseqNWtm/NoCCwBouD+8Hupy4uuhJrP7ySeiGDoS2YH7Io78c2QH7ovi5NHY/Q9PzPi1BRYA0HD1vB5qMtVqNfr7+2NdZSQqA3tiXWUk8jyP7u7uGb/2wgTzAQDMyOnhU1Hcv+eL10MN7Knr31utVqNarcYHH3wQra2tyV7XCRYA0HD1vB6qEZxgAQANt/vJJ6JWq41dB9W2ObKhw2PXQ9X5Qvd6cYIFADTcH14PtSjx9VCN4AQLAJgVPr8eai5wggUAkJjAAgBITGABACQmsAAAEptWYB08eDCyLIsXX3yx3vMAADS9awbW2bNn48c//nFs2LChjHkAAJrelIH12WefxSOPPBLPPfdcVCqVsmYCAGhqU94Ha//+/fGd73wnOjs7p/Vily5diosXL44/rlQqwgwAmHcmDaxf//rX0d/fH7/4xS+m/WIbN26c8HjXrl3R29s7rT87MjIy7b+HmbHr8th1eey6PHZdHrsuz0x3/ccfFD1pYP3Xf/1XnD17Nm6//faIiLhw4UJs37493n///Xjssceu+meOHz8eHR0d44+v9wQr5adYMzW7Lo9dl8euy2PX5bHr8qTc9aSB9dhjj00IqXvuuSd27twZW7dunfTFlixZEi0tLcmGAwBoRu6DBQCQ2LQ/7PmVV16p4xgAAHOHEywAgMQEFgBAYgILACAxgQUAkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkJjAAoAmkOd5tK/vjBsW3xjt6zsjz/NGj8QUBBYAzHJ5nketVovBK0vj8v17YvDK0qjVaiJrFhNYADDLPfX0M5G1dUWxYyCia2cUOwYiW7sp9u7ra/RoTEJgAcAsd3r4VBRruyKybOwLWRZF2+YYPvVmYwdjUgILAGa5latWR3bySERRjH2hKCIbOhyr1qxp7GBMamGjBwAAprb7ySeiVqtFduC+KNo2RzZ0OIqTR2O3a7BmLSdYADDLVavV6O/vj3WVkVg0sCfWVUYiz/Po7u5u9GhMwgkWADSBarUa1Wq10WMwTU6wAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACSyfM82td3xg2Lb4z29Z2R53mjR4KGmDKwLl++HFu3bo2VK1dGe3t7dHV1xZkzZ8qaDYAmkud51Gq1GLyyNC7fvycGryyNWq0mspiXrnmCtX379hgeHo5f/epX8cADD8QjjzxSxlwANJmnnn4msrauKHYMRHTtjGLHQGRrN8XefX2NHg1KN2VgLVq0KLZs2RJZlkVExIYNG+Ls2bNlzAVAkzk9fCqKtV0R//d/MyLLomjbHMOn3mzsYNAA13UN1oEDB+KBBx6Y9PcvXboUFy9eHP81Ojo64wEBaA4rV62O7OSRiKIY+0JRRDZ0OFatWdPYwaABFk73iX19fXHmzJk4duzYpM/ZuHHjhMe7du2K3t7eab3+yMjIdEdhhuy6PHZdHrsuz2S7/vudO2Lbtm2RHdgSRdtfRzb0UhQnj8Xf/+Qn8cEHH5Q85dzg33V5Zrrr1tbWCY+nFVjPPvts5HkeR48ejcWLF0/6vOPHj0dHR8f440qlEpVK5UsPR/3YdXnsujx2XZ6r7frhhx+Om266Kfbu64vhgT2xavWa2J3n0d3d3YAJ5w7/rsuTctfXDKz9+/fHCy+8EEePHo1bbrllyucuWbIkWlpakg0HQHOpVqtRrVYbPQY03JSBdf78+fj+978fK1asiHvvvTcixk6lXn311VKGAwBoRlMG1rJly6L4/GJFAACmxZ3cAQASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYMEsl+d5tK/vjBsW3xjt6zsjz/NGjwTANQgsmMXyPI9arRaDV5bG5fv3xOCVpVGr1UQWwCwnsGAWe+rpZyJr64pix0BE184odgxEtnZT7N3X1+jRAJiCwIJZ7PTwqSjWdkVk2dgXsiyKts0xfOrNxg4GwJQEFsxiK1etjuzkkYiiGPtCUUQ2dDhWrVnT2MEAmNLCRg8ATG73k09ErVaL7MB9UbRtjmzocBQnj8Zu12ABzGpOsGAWq1ar0d/fH+sqI7FoYE+sq4xEnufR3d3d6NEAmMI1T7DeeuutePjhh+O3v/1t3HzzzfFv//Zv0dbWVsZsQIxFVrVabfQYAFyHa55gPfroo7F9+/Y4ffp09Pb2xrZt20oYCwCgeU0ZWL/5zW/ixIkT8dBDD0VERK1Wi3PnzsWZM2dKGQ6gXtzAFainKQPr3Llzceutt8bChWPvJGZZFsuXL4/33nvvqs+/dOlSXLx4cfzX6Oho+okBZsgNXIF6S/pThBs3bpzweNeuXdHb2zutPzsyMpJyFKZg1+Wx6/Jcz67/cc/eyNo2jd3ANcui2PR4ZAe2xD8+9XTcc8899RtyjvDvujx2XZ6Z7rq1tXXC4ykD67bbbov3338/Pv3001i4cGEURRHvvfdeLF++/KrPP378eHR0dIw/rlQqUalUvvRw1I9dl8euyzPdXb995q0o7n/wj27g+tfx9sAe/3tNkz2Vx67Lk3LXU75F+NWvfjXuvPPO+OlPfxoREf39/bFs2bL49re/fdXnL1myJFpaWsZ/XU9cAZTFDVyBervmW4TPP/98bNu2Lfr6+qKlpSUOHjxYxlwAdeMGrkC9XfM2DatWrYr//u//jtOnT8eJEyfijjvuKGMugLpxA1eg3nxUDjAvuYErUE8+KgcAIDGBBQCQmMACKIm7x8P8IbAASuDu8TC/CCzgC5y0pPfU089E1tY1dvf4rp1R7BiIbO2m2Luvr9GjAXUgsIAJnLTUx+nhU1Gs7fqju8dvjuFTbzZ2MKAuBBYwgZOW+nD3eJhfBBYwgZOW+tj95BNRDB2J7MB9EUf+eewu8iePxu5/eKLRowF1ILCACZy01Ie7x8P84k7uwAQ+p69+3D0e5g8nWMAETloAZk5gAV9QrVbjjddOxCcffxRvvHZCXDUxt9yAxhBYAHOUW25A4wgsgDnKLTegcQQWwBzllhvQOAILrpNrWmgWbrkBjSOw4Dq4poVm4uam0DgCC66Da1rqx8lgem65AY0jsOA6uKalPpwM1o9bbkBjCCy4Dq5pqQ8ng8BcI7DgOrimpT6cDNJsvKXNtQgsuA6uaakPJ4M0E29pMx0+7Bmukw/sTc8HTNNMJrylnWVRbHo8sgP3xd59fb43MM4JFtBwTgZpJt7SZjqcYAGzgpNBmsXKVatj8OSRKDY9PhZZ3tLmKgQWAFwHb2kzHd4iBIDr4C1tpsMJFgBcJ29pcy1OsAAAEhNYAACJCSwAgMQEFk3Nx1UAMBsJLJqWj6sAYLYSWDStCR9X0bUzih0Dka3dFHv39TV6NADmOYFF0/JxFQDMVgKLprVy1erITh6JKIqxL/i4CgBmCTcapWn5uAoAZisnWDQtH1cBwGzlBIum5uMqAJiNnGABACQmsErihpgAMH8IrBK4ISYAzC8CqwRuiAkA84vAKoEbYgLA/CKwSuCGmAAwv7hNQwncEBMA5pernmBdvnw5tm7dGitXroz29vbo6uqKM2fOlD3bnOGGmAAwv0x6grV9+/b4m7/5m8iyLP7lX/4lHnnkkXjllVdKHG1ucUNMAJg/rnqCtWjRotiyZUtk/3dR9oYNG+Ls2bNlzgUA0LSmdQ3WgQMH4oEHHrjm8y5duhQXL14cf1ypVKJSqXz56QAAmtA1A6uvry/OnDkTx44du+aLbdy4ccLjXbt2RW9v77QGGRkZmdbzmDm7Lo9dl8euy2PX5bHr8sx0162trRMejwfWoUOHYv/+/RER8fjjj8ff/d3fxbPPPht5nsfRo0dj8eLF13zx48ePR0dHx/jj6z3B+uPhqB+7Lo9dl8euy2PX5bHr8qTc9Xhg9fT0RE9Pz/hv7N+/P1544YU4evRo3HLLLdN6sSVLlkRLS0uy4QAAmtFV3yI8f/58fP/7348VK1bEvffeGxFjp1GvvvpqqcMBADSjqwbWsmXLovj8ruMAAFwXH5UzR+V5Hu3rO+OGxTdG+/rOyN01HgBKI7DmoDzPo1arxeCVpXH5/j0xeGVp1Go1kQUAJRFYc9BTTz8TWVtXFDsGIrp2RrFjILK1m2Lvvr5GjwYA84LAmoNOD5+KYm1XxP/diT+yLIq2zTF86s3GDgYA84TAmoNWrlod2ckjEZ//oEJRRDZ0OFatWdPYwQBgnpjWR+XQXHY/+UTUarXIDtwXRdvmyIYOR3HyaOx2DRYAlMIJ1hxUrVajv78/1lVGYtHAnlhXGYk8z6O7u7vRowHAvOAEa46qVqtRrVYbPQYAzEtOsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkNi8DKw8z6N9fWfcsPjGaF/f6UOQAYCk5l1g5XketVotBq8sjcv374nBK0ujVquJLAAgmXkXWE89/UxkbV1R7BiI6NoZxY6ByNZuir37+ho9GgAwR8y7wDo9fCqKtV0RWTb2hSyLom1zDJ96s7GDAQBzxrwLrJWrVkd28khEUYx9oSgiGzocq9asaexgAMCcMe8+i3D3k09ErVaL7MB9UbRtjmzocBQnj8Zu12ABAInMuxOsarUa/f39sa4yEosG9sS6ykjkeR7d3d2NHg0AmCPm3QlWxFhkVavVRo8BAMxR8+4ECwCg3gQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxAQWAEBiAgsAIDGBBQCQmMACAEhMYAEAJCawAAASE1gAAIkJLACAxK4ZWAcPHowsy+LFF18sYx4AgKY3ZWCdPXs2fvzjH8eGDRvqNkCe59G+vjO+vuy2aF/fGXme1+3vAgAow6SB9dlnn8UjjzwSzz33XFQqlbr85XmeR61Wi8ErS2P0/j0xeGVp1Go1kQUANLVJA2v//v3xne98Jzo7O6f9YpcuXYqLFy+O/xodHZ3y+U89/UxkbV1R7BiI6NoZxY6ByNZuir37+qb/3wAAYJZZeLUv/vrXv47+/v74xS9+cV0vtnHjxgmPd+3aFb29vZM+f3j4VBT374nIsrEvZFkUbZvj1MCe+OCDD67r72b6RkZGGj3CvGHX5bHr8th1eey6PDPddWtr64TH44F16NCh2L9/f0REPProo3H27Nm4/fbbIyLiwoULsX379nj//ffjsccem/TFjx8/Hh0dHeOPK5XKlG8vrlq1OgZPHoli0+NjkVUUkQ0djtVr1nxhUNKy3/LYdXnsujx2XR67Lk/KXY8HVk9PT/T09Iz/xh+G1D333BM7d+6MrVu3TvliS5YsiZaWlmn/5buffCJqtVpkB+6Lom1zZEOHozh5NHa7BgsAaGINvQ9WtVqN/v7+WFcZicrAnlhXGYk8z6O7u7uRYwEAzMhVr8H6Y6+88krdBqhWq1GtVuODDz5wDAoAzAnu5A4AkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYlMG1ujoaHzve9+L22+/Pe6444546KGHypoLAKBpLZzqN3/wgx9ElmVx+vTpyLIsLly4UNZcAABNa9LA+uijj+Jf//Vf4/z585FlWURE/Omf/mlpgwEANKtJ3yJ8++23o7W1Nfr6+uKuu+6Ku+++O44dOzbli126dCkuXrw4/mt0dDT5wAAAs92kJ1iffvppvPvuu7F27dr44Q9/GK+//np0dXXF0NBQfO1rX7vqn9m4ceOEx7t27Yre3t5pDTIyMnIdYzMTdl0euy6PXZfHrstj1+WZ6a5bW1snPJ4QWIcOHYr9+/dHRMSDDz4YCxYsiAcffDAiItavXx/f+ta3YnBwcNLAOn78eHR0dIw/rlQqUalUvvRw1I9dl8euy2PX5bHr8th1eVLuekJg9fT0RE9Pz/jjI0eOxEsvvRRbtmyJd955J955551Ys2bNpC+2ZMmSaGlpSTYcAEAzmvKnCH/0ox/Fd7/73ejt7Y0FCxbE888/H9/4xjfKmg0AoClNGVgrVqyIl19+uaxZAADmBHdyBwBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgQUAkJjAAgBITGABACQmsAAAEhNYAACJCSwAgMQEFgBAYgILACAxgUUyeZ5H+/rOuGHxjdG+vjPyPG/0SADQEAKLJPI8j1qtFoNXlsbl+/fE4JWlUavVRBYA85LAIomnnn4msrauKHYMRHTtjGLHQGRrN8XefX2NHg0ASiewSOL08Kko1nZFZNnYF7IsirbNMXzqzcYOBgANILBIYuWq1ZGdPBJRFGNfKIrIhg7HqjVrGjsYADTAwkYPwNyw+8knolarRXbgvijaNkc2dDiKk0djt2uwAJiHnGCRRLVajf7+/lhXGYlFA3tiXWUk8jyP7u7uRo8GAKVzgkUy1Wo1qtVqo8cAgIZzggUAkNisCKzR0dH4p3/6pxgdHW30KHOeXZfHrstj1+Wx6/LYdXnqseusKD7/sa8v77XXXovOzs745S9/GXfeeed1//mLFy/GzTffHB9++GG0tLTMdBymYNflsevy2HV57Lo8dl2eeux6VpxgAQDMJQILACCxJD9F+Mknn0RExJtvfrm7dl+6dCkiIt54441YsmRJipGYhF2Xx67LY9flsevy2HV5Uu169erVsXjx4ohIdA3Wv//7v8dDDz0005cBAGhaf3gtepLA+u1vfxsvvfRSfPOb34wbbrhhxgMCADSb5CdYAAD8fy5yBwBITGABACQ2KwJr8+bNsW7duujo6Ii77747Xn/99UaPNCddvnw5tm7dGitXroz29vbo6uqKM2fONHqsOWvHjh3xzW9+M7IsizfeeKPR48xZb731VvzFX/xFrFy5Mv78z/88hoaGGj3SnOXfdDl8ry5X3RqkmAVGRkbG/3Oe58W6desaOM3c9cknnxQ/+9nPis8++6woiqJ47rnnio0bNzZ2qDns+PHjxblz54o/+7M/K15//fVGjzNn3XvvvcXBgweLoiiK//iP/yjuuuuuxg40h/k3XQ7fq8tVrwaZFSdYt9xyy/h//vDDDyPLsgZOM3ctWrQotmzZMr7fDRs2xNmzZxs71Bz2l3/5l7Fs2bJGjzGn/eY3v4kTJ06M3yamVqvFuXPn/H/7deLfdDl8ry5XvRokyY1GU+jp6YmXX345IiJ+/vOfN3ia+eHAgQPxwAMPNHoM+NLOnTsXt956ayxcOPatLMuyWL58ebz33nvx7W9/u8HTQRq+V9dfPRpk1gTWoUOHIiLiJz/5SfT29oqsOuvr64szZ87EsWPHGj0KAJPwvboc9WiQhrxFeOjQoejo6IiOjo44ePDghN97+OGH4+WXX47f/e53jRhtzrnarp999tnI8zz+8z//c/yGaMzcVP+uqY/bbrst3n///fj0008jIqIoinjvvfdi+fLlDZ4MZs736vKlbJCGnGD19PRET09PRET8/ve/j//5n/+Jr3/96xER8eKLL8bSpUujtbW1EaPNOX+464iI/fv3xwsvvBBHjx6d8L4zM/fHu6b+vvrVr8add94ZP/3pT2Pbtm3R398fy5Yt8/YgTc/36nL8/ve/j48//rguDdLwO7m/++678bd/+7fxySefxIIFC+JP/uRP4tlnn42Ojo5GjjUnnT9/Pm677bZYsWJF3HTTTRERUalU4tVXX23wZHPTo48+Gj/72c/iwoULsXTp0rjppptcfF0Hw8PDsW3btvjd734XLS0tcfDgwbjjjjsaPdac5N90OXyvLk89G6ThgQUAMNfMits0AADMJQILACAxgQUAkNj/A5G7EVey9DyRAAAAAElFTkSuQmCC" }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Load functions from packages\n", "using GaussianProcesses, Plots\n", "using Distributions:Normal, TDist\n", "using Random\n", "using Statistics\n", "\n", "#Simulate the data\n", "Random.seed!(112233)\n", "n = 20\n", "X = range(-3,stop=3,length=n);\n", "sigma = 1.0\n", "Y = X + sigma*rand(TDist(3),n);\n", "\n", "# Plots observations\n", "pyplot()\n", "scatter(X,Y;fmt=:png, leg=false)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We fit a standard (exact) Gaussian process model to the Student-t data and compare this against the Monte Carlo GP which is applicable for non-Gaussian observations models." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "GP Approximate object:\n", " Dim = 1\n", " Number of observations = 20\n", " Mean function:\n", " Type: MeanZero, Params: Float64[]\n", " Kernel:\n", " Type: Mat32Iso{Float64}, Params: [0.0, 0.0]\n", " Likelihood:\n", " Type: StuTLik, Params: [0.1]\n", " Input observations = \n", "[-3.0 -2.68421 … 2.68421 3.0]\n", " Output observations = [-3.98707, -2.7855, -1.74722, -2.27384, -7.12662, -1.43191, -0.840899, 0.305618, -0.444213, -0.924319, -0.0303161, -0.368707, -0.706642, 2.07517, 5.68693, -1.02634, 3.25863, 5.46091, 4.45861, 4.60741]\n", " Log-posterior = -77.863" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Build the models\n", "\n", "gpe = GPE(X,vec(Y),MeanZero(),Matern(3/2,0.0,0.0),0.5) #Exact GP assuming Gaussian noise\n", "\n", "l = StuTLik(3,0.1)\n", "gpa = GPA(X, vec(Y), MeanZero(), Matern(3/2,0.0,0.0), l) #Approximate GP with student-t likelihood" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimate the parameters of the exact GP through maximum likelihood estimation" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " * Status: success\n", "\n", " * Candidate solution\n", " Minimizer: [6.57e-01, 1.71e+00, 1.41e+00]\n", " Minimum: 4.578633e+01\n", "\n", " * Found with\n", " Algorithm: L-BFGS\n", " Initial Point: [5.00e-01, 0.00e+00, 0.00e+00]\n", "\n", " * Convergence measures\n", " |x - x'| = 2.08e-08 ≰ 0.0e+00\n", " |x - x'|/|x'| = 1.21e-08 ≰ 0.0e+00\n", " |f(x) - f(x')| = 1.42e-14 ≰ 0.0e+00\n", " |f(x) - f(x')|/|f(x')| = 3.10e-16 ≰ 0.0e+00\n", " |g(x)| = 4.69e-12 ≤ 1.0e-08\n", "\n", " * Work counters\n", " Seconds run: 0 (vs limit Inf)\n", " Iterations: 13\n", " f(x) calls: 38\n", " ∇f(x) calls: 38\n" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "optimize!(gpe)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Taking a Bayesian perspective, we can add prior distributions to the model parameters and sample from the posterior distribution using Markov chain Monte Carlo through the `mcmc` function which uses a Hamiltonian Monte Carlo sampler." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of iterations = 10000, Thinning = 2, Burn-in = 1000 \n", "Step size = 0.100000, Average number of leapfrog steps = 10.015000 \n", "Number of function calls: 100151\n", "Acceptance rate: 0.461500 \n" ] } ], "source": [ "set_priors!(gpa.lik,[Normal(-2.0,4.0)])\n", "set_priors!(gpa.kernel,[Normal(-2.0,4.0),Normal(-2.0,4.0)])\n", "\n", "samples = mcmc(gpa;nIter=10000,burn=1000,thin=2);" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xd8U/X6B/DPSUeatrRQ9mwBKeICFHCAICiIXmWKA1Cm4GIoCk4Qrlf5ObgiXhGQJaLIlYIMF8r2qogMEaEDuvfMaHbO8/sjprY0bUaTnJP2eb9efUmSc5Kntfn2yff7nOcrEBGBMcYYY4z5jELqABhjjDHGGhtOsBhjjDGZmDp1KgRBQEZGhtShsAbiBKsRyMjIgCAI9X4lJCRIHaZTCQkJDYotPz8fr7zyCm688Ua0bNkSYWFhiIuLw0033YTnn38ef/75Z61zbrvttho/G4VCgRYtWuDWW2/Fpk2bwKvmrKmoPna0a9cOVqvV6XHnz5+XbCw5dOgQBEHAq6++GrDX/OGHHzBx4kQkJCRApVIhKioKvXr1wuzZs/HLL78ELA5/4rHT/0KlDoD5Tvfu3TF58mSnjzVv3jzA0fjftm3bMGPGDOj1elx33XWYMGECWrZsCY1Gg9OnT+Odd97Bm2++iR07dmDs2LG1zl+wYAGio6Nhs9lw6dIlJCUl4dixY/jtt9+watUqCb4jxqQRGhqKwsJCfPXVVxg1alStx9evXw+FovF/HjcYDJg+fTq2bduGyMhI3HHHHUhMTAQApKSkYOvWrVi7di0+/vhjPPzwwxJH6z0eOwOEWNBLT08nAHTnnXdKHYrH4uPjKT4+3uPzvvrqK1IoFNSqVSv65ptvnB6Tk5NDTz75JK1fv77G/UOGDCEAlJ+fX+P+33//nVQqFQmCQJcuXfI4JsaCjWPsGDx4MMXGxtLo0aNrHWOxWKht27Y0YsQIUiqVXr1fG+LgwYMEgJYsWeL313rooYcIAA0fPpwKCgpqPV5eXk4LFy6k9957z28xTJkyhQBQenq6X56fx87A4QSrEfA0wXrjjTcIAM2ePbvOxx577LGq+3Jzc2nx4sV04403UuvWrSk8PJzi4+Pp8ccfp8LCQqevYTKZaMWKFdSvXz+Kjo6mqKgo6tWrFz399NNUVlZWFbOzL1cDqcVioa5duxIAOnjwoMvv12Kx1Lhd1yBBRDRy5EgCQNu3b3f5vIwFu+pjx+zZsyk0NLTWe3rnzp0EgLZt21ZngqXT6Wjx4sXUs2dPUiqV1KJFC7r77rvp2LFjtY5dsmRJ1Xt369at1Lt3b4qIiKB27drR3LlzSa/X1zrW2Vf1BMRkMtE777xDffv2pcjISIqOjqZBgwbRl19+6fbP4sCBAwSAEhMTqbKyst5jjUZj1b9PnDhBTz75JF199dUUExNDERERdM0119Abb7xBZrO51rmOD5Xl5eX05JNPUqdOnSgkJIQ2btxIRPUnWBs2bKABAwZQVFQURUVF0YABA6rOcwePnYHFS4RN0MKFC7F//36sWbMGI0eOxJgxYwAAx48fx+LFi3HVVVdhxYoVVccfOXIE77zzDm6//XbceOONCAsLw6lTp7B69Wp8++23OHnyJGJjY6uONxgMGD58OH788Uf06NED06ZNg1KpRGpqKtasWYNHHnkECQkJWLJkCd59910AwPz586vOv+222+qN/+DBg0hPT8egQYNcHgvYlz88JQiCx+cwFsymT5+ONWvWYMuWLViwYEHV/Rs2bEBcXFzVOHE5o9GIYcOG4fjx47j++usxf/58FBYW4vPPP8e3336Lzz77DBMmTKh13vvvv49vvvkGo0ePxrBhw/DNN9/gvffeQ0lJCbZu3QrAPhZkZGRg8+bNGDJkSI33u6PswWQyYeTIkTh06BD69OmDGTNmwGKxYN++fRg9ejRWrVqFp556yuX3v379egDAs88+i8jIyHqPVSqVVf9et24d9uzZg8GDB+Puu++GXq/HoUOH8MILL+DXX3/Fjh07ap1vMpkwbNgw6HQ6jBo1CqGhoWjbtm29rzl37lysWrUKHTt2xIwZMwAAO3bswLRp03Dq1CmsXLnS5ffIY2eASZ3hsYZzfArt3r07LVmyxOnX119/XeOcnJwcatmyJcXFxVFOTg5pNBrq3r07KZVKOnPmTI1jCwsLSavV1nrdzZs3EwB67bXXaty/YMECAkAPP/wwWa3WGo9VVFTUeC5vlgiXLl1KAOiVV17x6DyHuj6F/fHHH1XT3P6anmdMTi6f/b7mmmvo6quvrno8Pz+fQkNDac6cOURETmewHO/HSZMmkSiKVfefPHmSwsPDqXnz5qTRaKrud8xKxcbG0oULF6ru1+v1lJiYSAqFgnJzc6vud7VE+OKLL1aNB9VfX6PRUL9+/Sg8PLzG89UlISGBAFBaWprLY6vLzMysNc6JokjTp08nALVm8eLj46t+5tVn6xyczWAdPnyYAFCvXr2ooqKi6v6ysjJKTEwkAHTkyBGXsfLYGVg8g9WIXLx4EUuXLnX62Lx58zBy5Miq2x07dsT69esxZswYTJ48GZ06dcLFixexcuVKXHfddTXObdOmjdPnfPjhhzFnzhx8//33eOmllwAAVqsVa9euRWxsLFauXImQkJAa51Sf6fJWQUEBAKBDhw61HsvIyMCmTZtq3JeQkICpU6fWOvbtt9+uKtRMT09HUlISDAYD5s6dK9urLhnzp+nTp+OZZ57BL7/8ghtvvBGbN2+G1WrF9OnT6zxn8+bNCAsLw/Lly2vMXvTt2xdTpkzBunXrsGvXrlpF4fPmzUPPnj2rbqtUKjz00ENYunQpfvvtN6fv78uJoojVq1eje/fuWLp0aY3Xb9asGRYvXoxRo0YhKSnJ5SyWY1zp1KmTy9etrkuXLrXuEwQBTz75JDZs2IDvv/8eAwcOrHXMm2++CZVK5dZrbN68GQDw6quv1hhDW7RogSVLlmDSpEnYtGkTbr311nqfh8fOwOIEqxG588478c0337h9/OjRo/HYY4/hww8/BADcfffdmDt3rtNjk5KSsGbNGpw8eRLl5eWw2WxVj+Xl5VX9+8KFC9BqtbjjjjvQokULL78T72VkZNRKMocMGeJ0kHjnnXcA2AfDmJgY9OvXDzNmzMAjjzwSiFAZk53Jkydj0aJF2LBhA2688UZs3LgRffv2RZ8+fZwer9FocOnSJfTq1ctpYjJ06FCsW7cOp0+frpVg3XDDDbWOdzxHRUWFW/EmJyejvLwcHTp0cPrhsri4GIB9XPIXs9mM999/H9u2bcOFCxeg0+lqtCuoPj46RERE4Nprr3X7NU6dOgXAefnE0KFDAQCnT5/2MPKaeOz0PU6wmrixY8dWJVh1fcJ755138Oyzz6J169YYMWIEOnXqVPXJ691334XJZKo6Vq1WA7DPkPmLo1bB2cB12223VQ1uBQUFaN++fZ3Pk5+fj3bt2vknSMaCUOvWrXHvvfdi27ZtmDBhApKTk+u97F6j0QBAnfVDjvef47jqYmJiat3nqPmp/gGuPmVlZQCAc+fO4dy5c3UeV1lZ6fK52rVrh4yMDOTm5qJbt25uvT4A3HfffdizZw8SExPxwAMPoE2bNggLC0NFRQVWrlxZY3x0aNOmjUe1ShqNBgqFAq1bt671WNu2bSEIgtOfsbNjAR47A4UTrCasoqICjz76KKKiomCz2TBnzhycOnUKzZo1qzrGarXin//8J9q3b4/Tp0/XWC4kIrz55ps1ntNReJqbm+u3uG+55RYA9oJNxphvzZgxA0lJSZg6dSoiIiIwadKkOo91JEmFhYVOH3csSTlLpnzB8bzjx4/HF1980aDnGjhwIDIyMvDDDz+4nWD9+uuv2LNnD+68807s27evRknEzz//XGfhuaeF4DExMRBFEcXFxbVKNoqKikBEbv2MeewMrMbfOY7VadasWcjKysLKlSvx1ltv4eLFi3jyySdrHFNSUgK1Wo2bb7651hv7xIkTMBgMNe7r2bMnYmJi8Ouvv6K8vNxlDCEhIW5/WnUYOnQounbtimPHjuHIkSMencsYq9+dd96Jjh07Ijc3F2PGjKl3qT8mJgbdunVDWlqa0w9Vhw4dAoA6lxjd4UhanI0TvXr1QkxMDE6cOAGLxeL1awCoujLvnXfeqTWuXc4xK3Xx4kUAwD/+8Y9a9aZHjx5tUDzV9e3bF8DfP8/qPPkZ89gZWJxgNVHr16/Hf//7X0yYMAEzZszAU089hXvuuQdbtmzBp59+WnVcmzZtoFKpcPLkSej1+qr7y8vLMWfOnFrPGxoaitmzZ0OtVmPevHm1BkW1Wg2dTld1Oy4uDiUlJTAajW7HHhoaivfffx8KhQL33Xcf9u/f7/Q4d+s4GGN/CwkJwa5du7Bz50688cYbLo+fMmUKLBYLXnjhhRq1R7///js2bdqE2NjYOls8uCMuLg4AkJ2dXeux0NBQPP7448jMzMSzzz7rNMn6448/UFRU5PJ1hg4dioceegjJyckYN26c03M0Gg1efPFFrF27FgAQHx8PADh27FiN486dO+fWz85dU6ZMAQAsXbq0xlKgWq2uqptyHFMfHjsDi5cIG5G0tLR69+t6/vnnERERgZSUFMybNw+dO3euGigAe7+b6667Do8//jhuvvlmdO3aFQqFAk888QTeeecd9O7dG/feey80Gg2+/vprxMfHO70aZdmyZfj555+xZcsW/Pzzz7jrrrugVCpx6dIlfPPNNzh27FjVp61hw4bhxIkTuOuuu3DrrbciPDwcgwcPxuDBg+v9Xu+++2588sknmDlzJkaMGIHevXvj5ptvRlxcHCoqKnDp0iX88MMPEATB6RU8jLG69evXD/369XPr2IULF2Lfvn3YsmULzp8/j9tvvx1FRUX4/PPPYbVasW7duhplB5668sor0aFDB2zbtg1KpRKdOnWCIAiYM2cOYmNjsXTpUpw8eRLvvfce9u3bh8GDB6NNmzbIzc3F2bNncebMGfz00091Xg1d3fr160FE2LZtG7p27YoRI0YgMTERRITU1FT88MMP0Gq12LJlCwBgwIABGDBgALZv3478/HzcdNNNyMrKwu7du/GPf/yjwcuWDoMHD8acOXOwatUqXHPNNRg/fjyICDt27EBOTg7mzp3rcsx04LEzgCRrEMF8pr6u6NW/ysvLyWQy0fXXX08KhYIOHz5c67m+++47EgSBbrrppqouvmazmf71r39Rjx49SKlUUpcuXWjBggWk1Wrr7GNlNBrp7bffpj59+pBKpaLo6Gi66qqraMGCBVReXl51nFarpUcffZTat29PISEhHm+JkZeXRy+99BL179+fmjdvTiEhIdS8eXPq378/Pffcc3Tu3Lla59TXjZixpsTTXSDq6+T+yiuvUGJiYlXvq7vuuouOHj1a69jqndwvt3HjRgJQqzv5zz//TEOGDKFmzZo57eRutVppzZo1NHDgQIqJiakap0aOHEmrV68mnU7n1vfnsH//fnrooYcoPj6eIiIiKCIignr06EEzZ86kX375pcaxRUVFNH36dOrQoQNFRETQtddeS//5z3/o0qVLBICmTJlS43hXvf9cdXLv378/RUZGUmRkJPXv3582bNjg0ffmwGOn/wlEvP01Y4wxxpgvcQ0WY4wxxpiPcYLFGGOMMeZjnGAxxhhjjPkYJ1iMMcYYYz7GCRZjjDHGmI9xgsUYY4wx5mM+SbD0en2tTt+MMSZnPG4xxvzJJwnWhQsXcMMNN+DChQu+eLqAUqvVUofgFY47sDjuxsedcSuYfn7BFCsQXPFyrP4TTPF6GmuTXyL0dKNhueC4A4vjbpqC6ecXTLECwRUvx+o/wRSvp7E2+QSLMcYYY8zXOMFijDHGGPMxTrAYY4wxxnyMEyzGGGOMMR/jBIsxxhhjzMc4wWKMMcYY8zFOsBhjjDHGfIwTLMZY0Jk7dy4SEhIgCAJOnz5ddX9qaipuueUWJCYmon///jh37pyEUTLGmjKXCZbJZMJTTz2FHj164Nprr8XkyZMDERdjrIkyWl0387vvvvtw7NgxxMfH17h/9uzZmDVrFlJSUrBo0SJMnTrVT1EyxoKdSIQLFeS35w91dcDzzz8PQRCQkpICQRBQUFDgt2AYY02b1mxFUaUJ3VtE1Xvc4MGDa91XVFSEEydO4LvvvgMAjB8/Hk899RTS0tJwxRVX+CVexlhwKjUSkjII+XrC/GsUiAkXfP4a9SZYlZWVWL9+PXJyciAI9hdv166dz4NgjDG1yYI8rRFhCu8qF7Kzs9G+fXuEhtqHNUEQ0KVLF2RlZdWbYOl0Omg0mqrbSqUSSqXSqxgYY/L3a7GI73IIFtF++7cSwtAOAU6wLl68iLi4OLz++uv4/vvvoVKp8Oqrr+L22293ejwPVIwxb1QYLcjXGeG/yfq6DRkypMbthQsXYtGiRQCA8vJyCSLyTjDFCgRXvByr/wQyXp0V+DovFJd0NT/EHc0iXKO0IMRFjuUq1ri4uBq3602wrFYrMjMzcdVVV2H58uU4deoUhg8fjnPnzqFt27a1jq9voJKrYPtldOC4A4vj9h+NRUSZWaxKriwCgLj6lwid6dy5M/Lz82G1WhEaGgoiQlZWFrp06VLveYcPH0afPn2qbl/+wfDyQVPOgilWILji5Vj9J1Dx7kkRkW8jqFQ17xcBFAgKXBvnehbLk1jrTbC6dOkChUKBSZMmAQD69u2Lrl274uzZs04TLFcDlVwF2y+jA8cdWBy375XozdCLJkRUG/DCvVwibNOmDa6//np88sknmDp1Knbs2IFOnTq5rL+Kjo5GTEyMV6/JGJPWiWIR/Vq7N2YUG+ueIz9eRG4lWJ6oN6pWrVrh9ttvx7fffgsASE9PR3p6Onr16uX0eMdA5fgKhuSKMSaNokoTivQmr86dPXs2OnXqhJycHNx5551VSdSaNWuwZs0aJCYmYvny5di4caMvQ2aMychvxYSDee4VFhitBJ2l7sezKwkFet8WKbi8ivDDDz/EjBkzsGjRIigUCqxZswYdO3b0aRCMsaaDiFCoN6HMUM9o58KaNWuc3t+zZ0/89NNPXj8vYyw4FOgJX2eLsBKgtxIiQ+uffSoxun7O48WEUfG+m8VymWB169YNBw8e9NkLMsaaLiJCfqUJFUbvkyvGWNNmshH+e8meXAH25KlLdP3nlLgxWX62jDC8I0HlIllzF3dyZyxAkpKS0Lt3b6hUKvTu3RtJSUlShxRQRIQ8nZGTK8ZYg+zOJJRWS5iKDa7PKamn/srBIgKnSn23TMgJFmMBkJSUhPHjx+Ps2bMwGo04e/Ysxo8f32SSLJEIOVoj1Car1KEwxoLY8SIR58prJkHuJE/uLBECwIliApFvkixOsBgLgKVLl0IQhKo3LhFBEAQsW7ZM4sj8TyRCjsYArZmTK8aY9/IqCd/m1E5+it1IntxJwgCgzASkaVwf5w5OsBgLgJSUlFqfiogIycnJEkUUGDaRkKUxQGdxvb8gY4zV52A+weYkT3KVPIlEKPfgguXjxTyDxVjQSExMrNpuykEQBPTs2VOiiPzPKhIyNXroOblijPlAbqXzxEdtBszOMq+/lJvgNDGryyWNb5YJOcFiLACWLFlStSwIoGq5cMmSJRJHVreGFOVbRBGZaj2MVtGPETLGmopyE0FfR5UBATWK3i/nbv2Vg40Agw8+F3KCxVgAjBs3Djt27MB1112HiIgIXHfddUhKSsLYsWOlDs2phhTlm20iMtUGmGycXDHGfCO3sv7H67uS0N36q+q0Zo9PqcVlHyzGmG+MGzcO48aNkzoMt9RXlF/f92Cy2pClMcAiSrFtM2OssapredDBnkQ571/l6QwWYN8YuvaGgJ7hBIsxVos3RfkGqw3ZagOsPrrEmTHGHPL09T9e35WE3iRYWh+06+MlQsZYLZ4W5estVmRxcsUY8wORCPku9gmsbxnQmyVCnYWL3BljfuBJUb7ObEWWxgAbJ1eMMT8oNgBmFyWdZSZ7Ina5Sgt5VbDOM1iMMb9wtyhfY7IgW2MAl1wxxvwl18XyIGC/8q/MyZWE3iwPAoDOBwkW12AxxpxyVZRfYbQgX2cE51aMMX9yVeDuUGwAWkXUvM/bBMsXM1icYDHGPFZqMKOo0sTJFWPM7/Jc1F85OLuS0Jv6K8A3NVicYDHGPFJUaUKJwQdNYhhjzAWrSCisp8dVdc6uJCzxYIuc6rgGizEWMESEAp2RkyvGWMDk6+F2jaez5UBvZ7AsImDyZH8dJzjBYoy5RETI05lQZvTBxzrGGPtLhYvPa+7WXwH2ZKp6/z6rSKjwcgYLaPgsFidYjLF6iUTI1RqhNnFyxRjzrd/KQpClqzuJcucKQgezCGiqDVOlRjSoTrShVxJygsUYq5NIhGyNARpzHbusMsZYA1zUKfBTYT0JlgczWEDNPQnr2wDaHTyDxRjzC5tIyFQbUGnxwbbyjDF2mXITodws4EIFodxUO5EyWMlpb6v6VK+58rb+yqGhVxJygsUYq8UqishU62GwcnLFGPOPNM1fm8kDOF5UO5lxtf+gM9WvJPS2B5YDLxEyxnzKbBORoTbAaHOxNwVjjDVAmvrvf58spVpX7Xm6PAjUTKoammDxEiFjzGdMVhsy1XqYOblijPmRTSSka/9OoEw24FTJ5QmW589b7NMlwgadzgkWY8zOYLUhU22AhTcWZIz5Waau9gbOvxTXbLOQ62YH9+r0VkBvJWjM5HKDaFe0XIPFGGuoSosVWWoDrE52o2eMMV9z1F9VV24CLlTY/602k9czSMWGhi8PAjyDxRhrIK3ZimyNATZOrhhjAeIswQKAn/4qds/zYnnQocTY8OVBADDY7M1KvcV7ETLWhKmNFuTpjLxpM2MsYDRmQlEd+wtm6Qh5leTV8qBDsZHc3l7HFZ0FaK707lyewWKsiSozmGWXXO3fuxv3DL4JKpVK6lAYY36Sqq7/8Z+LyKsCdwf7DJb351fXkCsJeQaLsSaoRG9CkV5emzbv37sb86dNhiAINQpdGWONS13Lgw7nygkhgvfPX2wk+GoIaUgdFidYjDUhRIQivRmlBnklVwDwwVvLOblirJETqWZ7BmdsZP/ylsaHw5v9SkLvsj1OsBhrIogI+ToTKmS6aXPGxVROrhhr5LJ1gNHPG0T4chRpyAwW12Ax1gSIRMjVGgOSXO3fuxtjh9yCvp1aY+yQW7B/7263zkvo3gOC0IB1AcaY7LlaHpSbhtRgcYLFWCMnEiFbY4DGbPX7aznqqFLPn4PZZELq+XOYP22yW0nWE889DyLiJIuxRixVHVwJlq4BwyYnWIw1YjYiZKoNqLQEZtPmy+uoHAnT6reXuzx3+D2j8O7GT9DzqmsQERHh71AZYwGmsxAK62jPIFdas/cJISdYjDVSFpuIAqMIgzUwyRXgvI6KiJCelurW+cPvGYU9R36CwRBkozBjzKU0jW/rowKBZ7AYYzWYbCIy1HqYA7yvoLM6KkEQ0PWKxIDGwRiTlwoT4ceC4NtEvtJiL7PwBidYjDUy9k2b9ZJs2nx5HZVjufCJ554PeCyMMXnI0RE+ShZR7KPmn4FEsCdZ3uAEi7FGpGrTZgmSK+DvOqrEq65GuFKJxKuuxspNW3HHP+6VJB7GmLTOlRE2p4oN3jhZSt4uE3IfLMYaCa3JilydwWd7cHlr+D2jMPyeUdIGwRiT3NF8EQfyKOjqri6nNQPtIz0/jxMsxhqBCqMF+TLbV5Ax1jRZRcK+LMKp0sYxIvEMFmNNVKnejEK9SeowGGMMRQbCjnQx6Nox1Mfb7XI4wWIsSMl5X0HGWNNzvEjEdzkEa+OYuKribf2YW0XuGzduhCAI2LVrl3evwhjzKce+go0puarUabHjk48xZRwXxDMWTPRWwmdpIr7KbnzJFeD9djkuZ7AyMjKwbt063HTTTd69AmPMpxz7CmoDsPVNIJw99Ru2b9qAr3ftgNFgwMDbhkkdEmPMTSVGwuYUsUF79smdX2awRFHEzJkzsWrVKiiVSu9egTHmMzaRkKUxBH1yZTaZsGvbVtx/xxA8OGIofjpyCDPnPoP9p85h4xdfSh0eY8xNR/OpUSdXgKMGy3P1zmCtWLECAwcOxA033ODWk+l0Omg0mqrbSqWSEzPGfMRiE5GlMcBkC75uyA7lpaX4bMNabNv4EUqLizFo2B34z9bPcevtIxASEiJ1eIwxD2jMhD/KG+Ga4GW8ncGqM8H6448/sGPHDhw5csTtJxsyZEiN2wsXLsSiRYu8iyxAysvLpQ7BKxx3YEkdt1kkFBptHtc3GA3yaJ2cn5ONT9auxq5tn0IQgHsnPIgHps1EQvcrAABm89+1ZBYBQFyURJEyxtz1SxHB1vjzK9jIXmfmqToTrKNHjyIjIwM9evQAABQUFGDWrFnIz8/H448/7vScw4cPo0+fPlW3g2UGKy4uTuoQvMJxB5ZUcestNpRoDAiLIIR5cb5KpfJ5TO7KzkjH2n+/jd3bP0N0TAxmzp2PiTNmoXlcyzrPCVfwBhOMyZ3JRvitpAlkV3/RWTxvu1Dn8Y8//niNROq2227D/PnzMWbMmDqfLDo6GjExMR6GwBiri9ZsRa5W+u7snsrLzsIHby3H7u2foXlcSzy9eBnuf2QaIqN4ZoqxxuBkCcFokzqKwNFagBYensN9sBiTqWDszl5SVIS1/34Ln2/egNjmzfHsq69hwiPToIr0Yp8JxpgsiUT4uSiYRqaG01mAFh72GnU7wTp06JCH4TDGvFWiN6MoiLqzG/R6bF79PtavehehoaF4auGLmDhzNqKio6UOjTHmY+fKAXXjacHnFq2FgHDPzuEZLMZkhIhQWGlCmTE4rnsmIuz94nP8+5+vorSkGJNmzsasp59F8xbBWWfHGHPtf4XBeyWzt3QWcILFWLAiIuTqjNCYgqPH1R+nT+KNFxfi9K/HMeLeMXhm8VJ0TugqdViMMT9K1xDy9VJHEXje9PriBIsxGbCJhBytAZUW+VeNqivKsfJfy7B98wb06HUVNuzcixsHDZY6LMZYAPyvidVeOXjTC4sTLMYkZhXtDUSNVnlPuzuWA99a8hKMBiOef205Hpz+KEJDeRhhrCnjCLChAAAgAElEQVQoMhDS1E0zwfKmmzuPjIxJyPxXd3azzLuz52RmYOmz8/G/Qwdw19jxWLTsDbRu107qsBhjAaIxE/ZkUlBd1exLPIPFWBAxWG3I1hhglXGTK1EU8cna1XjvjX+ieVxLfPjZF7j1jhFSh8UYC6A/ywl7MkUY5F/B4Ddm0f7lCU6wGJNApcWKHI0RNpJvcpV56SJenvsETv7yEybNnI15Ly1GVHQzqcNijAWI2Ub4OptwqlS+41QgeTqLxQkWYwGmMVmQpzPKtjs7EeHT9WuxYtlitG7bFpu//Br9bhkodViMsQDKqyTsSBdRGjzt+PxOb/Os0ygnWIwFkNy7sxcV5OOlOY/jf4cO4KHpj+KZxct4exvGmpgzpYTdmWKT2MjZE57+PDjBYixASg1mFFbK9+Pg9/v2YPHTTyE8XIk1nydh0LA7pA6JMRZgh/JEHMrnzMoXeNt6xgKgqNIkSXJ14Ot9GDvkFvTt1Bpjh9yC/Xt31zrGoNdj6bPzMW/qJPS/ZRB2HfmJkyvGmhibSNiVwcmVL/EMFmN+REQoqDShXIKtb/bv3Y2Fs6dDEAQQEVLPn8P8aZPx7sZPMPyeUQCAtAvn8czMKcjNysSSt9/FhEemQRA83NGUMRYQROSX96fRSvj8EiFdy8mVL/EMFmN+QkTI0xklSa4A4IO3llclV454BEHA6reXg4iQtHULHhhxGxQKBT7ffxj3T5nOyRVjMpas9v1zqs2EDSkiJ1d+wDNYjPmBSIQcrRE6s3T7CmZcTK1KrhyICOlpqXjxqdnYvX0b7nt4Kl741/8hQqWSKErGmLvOlRNaRwAtI3zzQajMSNicKkJt9snTscvwDBZjPmYTCVkag6TJFQAkdO9Ra0bKcXv/3j34vw8/wtIV73FyxZo0IkKFKThmb0qMhAsVvom1UG+fueLkyn84wWLMh6x/JVd6GWza/MRzz9eo2XAsF7aIa4lt3x3EPePvlzhCxqRXbkJQ9HoiIpQYgfMVDX+u3ErCphTRq+1fmPs4wWLMR6yiiEy1Hgar9MkVAAy/ZxTeXLMBPXpdDUVICIgIfW+8CXv+dwJX9LxS6vAYk4UiI6ANgkRDbQYsoj050pq9n8XK0BI+Tm3a294ECidYjPmA2SYiQ22ASWabNvfu1x+xzZtDIQh4/rXl2LLnW0RFR0sdFmOyUWggaBqQsARKidH+XwJw3stlwjQ1YWuaCBMnVwHBRe6MNZDJJiJLrYdFZnvfnD31G+ZNmQSr1YoNSXtxw823SB0SY7JTZAAig+AvYbHx7/HlQgUwoI1n56epCdsuirDKa5hq1HgGi7EGMFptyJRhcrVr21Y8cu9ItG7XDl8cOMrJFWN1KDIQNEFQ6O2YwQKADB3B4EGmdEnjWXJ1/sBOrJt4A94YGI11E2/A+QM7PYyWAZxgMeY1g8WGTLUBVhklV1arFctffh4vzXkc94y/H2u370Kbdu2lDos1MefLCfuy5LVc7oxNJJSagqMGq7hagiWS+z2x0jWEzzxMrrYvnIDuplwsGhCP7qZcbF84gZMsLwTBxChj8qO3WJGtMcJG8kmuKsrL8Nysafjl6BG8+MZbmDhjFoxGo+sTGfORAj3hmxxChpagCgHu7uyfzuO+UmK0Jysai3zex3UpMdaM8Xw5oU/L+n+2GVrCpxdFWDzIdY999BpujW+FraOvhyAIeLRvPCbuOokf17+GXsPGehN6k8UJFmMe0pmtyNEaIKOJK1xMScZTkx+ARl2BdV98iRsHDZY6JNaEVFoIB/IIJ0sIjreFwWZvf9AqQtLQ6lX01+ePSou9ObBCpslgpYWgv6yt3kUNwWwjhIc4jzlTS/g0zbPkCgBKMpMxY0B8jfYuQ+Nb4vjxZG9Cb9J4iZAxD2hlmFwd/f47TBx5O5QREdj27UFOrljA/Ted8Fu15MohWydJOG4rMvy1jRTkvUxY4mQi2kpAmsb58dk6+9WCZi9WaVvF98ShrNIaW2wdzCxF64Senj+ZHwVDnRgnWIy5SW2yIEcjn+SKiLBp9ft4YtL96HfLQGz9aj86J3SVOizWBJUanb8pcipl8mapQ6Hh739rZVzoXlzHz9dZu4Z8g+B1cgUAg2a+jKOZJZi46yTWnszAxF0ncSyrBANnvuLdE/pBsNSJcYLFmBu0FhF5WmOtT+hSMZvNWPL0HLy1+EVMe3Ie3tv8KaKim0kdFmuCLCLVOfuTLfMEyzGDBQCaIJvBAoAUNcFW7RNfgZ6wPSsMxgb0ueo1bCzuf/O/SFd1xP8dz0S6qiPuf+sL9Bo6xvsn9bHqdWKzrk/A1tHXY1CXVvhx/WtSh1YD12Ax5kKZwYxSs4gImWzZV1FWivnTHsapX3/Bv1atxpgHJ0kdEmvCyuq5jqLYAJhsBGUddUJSMtuoxj58WgsBkF+cQM0rCKsz2YBLWqBHrD1Z/DhVhNEGNHSo6jVsrKwL2oOlToxnsBirR4nejIJKk2xmrtLTUvHgncOQlnweG5L2cnLFJFdWzz5+BCC3MmCheKTIgBrvazn3wrr8CsLqLlQQSoyEzSlirUL4xipY6sQ4wWKsDkWVJhTp5bML7M9HDmHiyNsRrlTis28O4IabbpY6JMZQZqr/44dclwmLLpsVkmuR++UzbZc7X2FPriqbSHIFBEedGMAJFmNOFVaaUGKQz0fa/27ZhNkPjMM1fa/nYnYmK/XNYAFAjmxnsGomfnKtwaqr/spBb5VvcugvwVAnBnANFmM1EBEKKk0oN8pjxBJFESuWLcbG/7yHB6fNxAuvv4nQUH7bupKQkAClUgmVyl6N8sILL+CBBx6QOKrGyWWCpSMQya/haJGh5m2tTDd8rqv+qqmTe50YwAkWY1WICHk6I9Qmecy16ysr8cKTs3Dg63144V//h0mPPia7P1Jy9vnnn6NPnz5Sh9HouVoilGvD0cKgmcGSZ+LHXOMEizHYuzjnao3QmuWRXBUXFODJyQ/gUmoKVm35DLeNuEvqkBirxSq6t1Fytk5eCValhWrVLFlEwGglRITK60MMz2AFL67BYk2eSIRsjUE2yVXKn+fw4MhhKC4swJa933Jy5aVHHnkE1157LWbMmIHi4uI6j9PpdNBoNFVfJpN8LmyQu3IT3LrCVm4NRy9fHnSQ4ywWz2AFL57BYk2aTbQnV3prAzrz+dCPB3/A09MfQeeErvjP1s/RrkNHqUMKSkeOHEGXLl1gsVjw8ssvY8qUKfjqq6+cHjtkyJAatxcuXIhFixYBAMrLy/0eq69IEeslrQCDIczlcclFhIHNamYvUv5sU8sUMBhq//nLLrYgNLp2QiNVrDYCcivCPdo9wmAIrimvYIpXrVajLKru/xlxcXE1bnOCxZosqygiS2OA0erlnhI+tuOTj7H02XkYOPR2vL1uI3dmb4AuXboAAMLCwjB//nwkJibWeezhw4dr1GoplUoolcqq25cPmnIW6FhFiwiVyvVf/0oAUbGKWg1HpfrZmrTO41ZERiEuzvkSoRSxFhsIygjPxyfHxR3BIljijY2NRVxcC7eP5wSLNUkWmz25MtmkT66ICO+98U+s/ffbeGDqDLz4xlt8pWADVFZWwmKxoHnz5gCAzz77DH379q3z+OjoaMTExAQqvEbF1RWEDo6Go91k8mOua4lQbt3cuf4quPEozpock1VElkYPiwx2bbbvKfgUdm/fhmcWL8P0p+bxlYINVFhYiPHjx8Nms4GI0K1bN3z88cdSh9UouZtgAfaGo91i5PG7XVRHXZPcurlz/VVw4wSLNSkGqw3ZagOsJP3ApdWoMX/aw/jt5//hrbUbcPfY+6QOqVHo1q0bTp06JXUYTUKpBwmAXBqOVpgIpjpKLuXWsJNnsIIbX0XImoxKixVZMkmuigry8ciou/DnmdNYt30XJ1cs6NhE8uiqO0fDUanVtTwIABqL9PFVxzNY8nD+wE6sm3gDBiZ2wA19+yApKcmt8zjBYk2CxmRBtsYAmwwG+LTkC5h41x3QlJdjy95v0X/gIKlDYsxj5WZ4dnXbXw1HL2cTCXpr4N6XdS0PAvKawSIil9vkMP87f2Anti+cgO6mXDzbPx7R5XkYP368W0kWJ1is0Ss3WpCrNXr0x6A++/fuxtght6Bvp9YYO+QW7N+72+1zT/7yMx6+ZwSimzXD1q+/xxVX9vJNUIwFWJkXf/yzdTVvp2sIq8+LOFEcwASrnhmsSos94ZMDtdne/JT5lmM26o2B0Vg38QacP7Cz3uOPffQabo1vha2jr8es6xPwyei+uDW+FV5/7Z8uX6veBMtoNGLMmDFITExE7969MXz4cKSlpXn23TAmoRK9Gfk6o1vNEN2xf+9uzJ82Gannz8FsMiH1/DnMnzbZrSTr4LdfY+Z9o5DY62p8vOebRtXjKkQQEBkaghYRYWgbqUTHZhGIj1WhW/NIXNEiCj3i/vpqEYXuLaKQEBuJzjEqdIiOQCtVOGKVoYgIVcjo+i3miqstcpxxNByttBD25IZic6qIEiOQrPZ1dHW7fJPn6giATh79hrn+yg+qz0YtGhCP7qZcbF84od4kqyQzGbd1aVl18ZEgCLitS0ucv3DB5eu5nMGaNWsWkpOTcebMGYwePRozZ8704NthTBr2TZuNKNL7tiv3B28thyAIVbUkjk1sV7+9vN7zdn76CeZNmYhbbx+Btdt3Iia2uU/jCrQwhYDmyjB0iI5A9+ZR6NkyGgnNI9E+OgItI8MRqwxDVFgoIkJDEB6iQJjir68QBZQhCkSGhaBZeCiaR4ShTZQSHZup0K15FK5sGY3OMTLaU4XVyZMrCB2yKwknikW8f07En+q///zkVVJANlsW3Vh2k8uVhFx/5XuXz0ZtHX09BnVphR/Xv1bnOa3ie+JQVmmNMf9QVil6XXmly9erN8GKiIjA3XffXZW53XTTTcjIyPDg22Es8BybNpcZfF9QkXExtVahLhEhPS21znPWr3oXL897AuMmPYIV6zdDGRGcCYQyRIHWkeHo1jwSPeKi0aFZBJpHhEEZ6rtKA0EQoAwN8dnzMf/xJsEqMgB7swiGy67iIwApAZjFyq0EXJV7yaUOi2ewfM/ZbNTQ+JYozkiu85xBM1/G0cwSTNx1EmtPZmDSlydxNLMEL72y2OXreTQyrly5EqNHj67zcd7Ti0nNJhKyNAaoTf6Z50/o3qNWnypBEND1itqdwokIb7/6MlYsW4zHFizEkrffRUhIcCUPIYKAlip7UtW9RRRaRyoRwQkQg3dLhPVJVvtvxkZrJuzNFLEpxXVRkyYAM2n1ISL8WizijzKewfI1Z7NRBzNL0TqhZ61jRSJoLYTYG8dg4JLtOBfWEW8ez4Q+rhOSkpIwduxYl6/ndh+s119/HWlpafjhhx/qPKa+Pb3kKpj2GquO467NSoQiowiTH4pUjX/tlzVz3jNYOHt61TKh478z5z0Dg+Hv6lmr1Yp/Pb8Ae7Zvw7OvvoYHpz8KozHwH0mNXu7zFa4QEBsmIDJEgMJkht4E6H0cW32CaXuapkgkgtrHS2npWoJFJIQpfFeJZ7ASjhUQjheT2wXjUs5gqc2ELzMIl7ScXPnDoJkvY/vCCZi46ySGxrfEwcxSHMsqwT9eX40LFYRCA1BoIBQYgIrqG5l3HAPFojEY0taKA6Pc39bHrQTr7bffRlJSEr7//ntERkbWeZyrPb3kKlgHc477byabiGyNAQqlCH/taqVSqfCPcfchPDwcq99ejvS0VHS9ogeeeO4F3PGPe6uOM5tMeH7OYzj4zVdY/sFa3DvhQT9F5B5P9vmKCgtBS1U4osO5BzGrW4XJvhGxL1lE4JIG6NnA8kSNmVCgt9d7/VpMMHq4j7snvb186VQJ4Zscsc4mqKzheg0bi/vf/C+OfvQafjmejPAOPRH25AfY12o0cImgCgHaRQI9Y4E4pYCYMKBZOBATBkSFAmPaGgEP/sK4HEVXrFiBzz77DN9//33V3l514T29mBT0FhtyNIFrIDr8nlEYfs8op49V6nSYO2UiTh3/GSs3f4qhd94VkJgaKiosBK0jlYgM4+U/5po39VfuSFYTejZ3fwarwkTIrQTyDYR8PVCgJ1Q2sDog0DNYajNhXxYhxY9LpMy+pP1HGZDacTTyn7WXOrWOBK6IATpFCWirApqFwadbldWbYOXk5GDBggXo1q0bhg4dCsA+K/XLL7/4LADGGkJjsiBP57seVw1RUV6Gxx+6DxeTk7FmW1JQNBCNCFGgTZSSZ6yYR3xdf+WQoqaqpXdndBbCbyX2pCq3suHJlDOBqsGyiIQfCwg/Frq/fMk8YyNCcgXwWwnhkhYIVwDdY4BR8QJ6xADRYf5tDFPvqNqpUydZbG3AmDNlBjMKK00+63HVEMWFhZh1/xgUFeRj4669uLp3X6lDqleoQkCbSCVilaG8uTTzmL9msHQW+5V+naKdP/5VFuHPCv++4+ubwUpVE9pHNvwP87kywv5cERUyaQnR2KjNhBPFhFOlQKUV6BQFjI4XcHUL+LTGzxX+2MqCDhGhUG/ySxsGb+RlZ2HG+FEwGgzYvPsbXNHTdX8UqQgAWqjC0FqlREgABxrWuPgrwQLsy4Sdomv/bqaq/Z9cAfZaMIOVoAqtHcOPhQSbCExNhFfvn3w94dscQgYXsfuF1kw4Wkj4rQQIUwDXxQE3tBLQViXNWMcJFgsqIhFytUZozfJot5xxMRUzxo9GWFgotuz9Fp3iE6QOqU6q0BC0j+Y2C6zh/LVECNiXCW+/bJMDi0j4Kjtw62haC6C67K9juYmQqSUQgH3Z9mUmd+VWEo7kk19bUTRllRb7UuuvxUCoAhjaXsCA1kB4iLQfIjnBYkHDKorI1hhhsMrjMpvkc3/g0Qmj0TyuJT764ku0adde6pCcUggC2kSGo0VEGC8HsgYjIlT4cQar0GAvXm+u/Pt39XA+oTyAbRU1ZqDNZReLnSmlqnKEkyWE9pEi+reuv5Vkto5wOJ+QpuHEyh+sor0Nx09F9tsD2wI3tRUQIXFi5cAJFgsKRqsNORojzKI8qkF/P3kCsx8Yh05dErB2+060aNlS6pCcUoUI6NY8EuEhvK878w212XU39IZKVhNubGP/I1lkIPxUGNgE5fI6LCLgzGWNP7/JJrSJIMQ3q/3H/KLGXsDO/az8p1BPSMoglJqAG9sAA9sKiHSyrCslTrCY7OnMVuRqjbDJ5IKLX388hicm3Y/Eq67Gh9u+QLOYWKlDqkUhAG0ilSDBzMkV8yl/1l85JFfY/2gSEfZmkc97brmisRBQbevxbL1QawbNRsD2SyJm9VIgNlyASPY2AP8rFFFgAPMTkQg/FwEH8ggtlcCjPQW0jZRXYuXACRaTNTldKQgAR3/Yj3lTJ6HvgJuw6uPPEBkVJXVItahCQ9ChWQSUIQqUGeQ58LDg5c/6K4dMHcFkI5wrB7J0gX/3ay+7uu8PtfO6xUorsO2iiOviBPxc5Pvu9qymChNhVyYhUwfc3AYY1kFAqIwv1uEEi8kSEaGw0oQyozyuFASA7/ftwYJHp2LQsDuw4iP5bdosAGipCkfryHCutWJ+E4gZLBvZa54O5knz0ap6N3ezjXBBo0BoHZuS5OvtVwcy/zpfTvgyk6AMAab0EJDgZGlWbjjBYrJjEwm5WgN0FnkUswPA1zu/wKvPzMXwe0Zj+ep1CAsLkzqkGsIUAjpERyCKG4YyPwtEggUA3+YEfmnQQWv5+4XPV9hbN/A7SxpEhEP5hCMFQK/mwKguAiJkVmtVF/6dYbJi/mtPQZNNHsXsAPDfLZuwdME8jH5gIpa9+z5CQuTV5iA6zL4kGKrgWivmf4FYIgR8v9ehJzTVlvpOl/LslFRMNnshe4ravhw4qK1vt7LxN06wmGxUWqzI1RgDtqegOz7+8D/4v1dewIQp07H4zRVQyCiJEQC0jlSipYrbL7DAIApsuwSp6K32mXSdFdwUVCKlRsK2iwStBXiou4DE2OAb4zjBYrJQbrSgQGeUTTE7AKxZ8Rbee+OfmDHnaTz27CJZJVdhCgEdm0UgMozfwixwUtRoEvvmEeytGs6UkazGpKYiTU34IoPQLBSYeaWAVhHBl1wBnGAxicmxmJ2I8O5rS/HReysw5/mXMfuZ52A0GqUOq0pUWAg68pIgk8BPRU0n3dBY7IX2LLB+KyHsyyJcEQOM6yqfpqHe4ASLScYmEnK0BlTKqJhdFEW88eJCfLp+LRb+83VMeewpqUOqoaUqHG34KkEmgXx909pD71w5BaygvylLOfwljn/8FkoykxHRMRGaES+h34ixuKuzAEWQj3P8EZhJwmi1IV2tl1VyZbPZsHj+U/hswzoseWelrJIrhSCgU7MItI1ScnLFJPFzgLupS+1kSdP6fqVw/sBO7Hp5MrqbcrFoQDyus+YBHz6Arqm7gj65AngGi0lAY7IgT2eCKKNidrPZjBeemIX9e7/E8g/W4p77HpA6pCrhIQp0ahbBmzQzyWjNhD/K5fN+DYSmUGsmtWMfvYZb41th6+jrIQgCHu0bj4m7TuJ/G/6Fq24fJ3V4DcYzWCxgiAjlZhG5WqOskiujwYD5Uyfhh6/34t8btsgquYoOC0HX2EhOrphPGKwEm+j5e+94sXQ9qVjjVZKZjNu6tKyalRcEAUPjW6I4I1niyHyDEywWEPZ6KyMqLKKsrsqp1Gnx+MQJOP7jUXywdTtuv/seqUOq0lIVjs4xKoTIeCsIJn9GK+F0KWFrmoi3fxeRovbsfItIOFEsp3ctawyKDQS0TcShzFLQXx+4iQgHM0vROqGnxNH5Bi8RMr8zWW3I0Rpl1TwUACrKy/DYg+ORnpqKNZ/vxA033Sx1SADsGzW3j4pAbIS8usWz4FJiJOzPIaRpas4+nSkj9GrhftJ+upRgkE+pJGsEsnSEzy4SIka/hKOrHsDEXScxNL4lDmaW4lhWCe5/60OpQ/QJTrCYX2lMFuTrTLDJaEkQAIoLCvDo/WNQUlSIjTv34qrefaQOCYC9v1WnZiqownhJkDXM6VJCsrr2+y5VTdBbCZFubDdCRE2uuJ3514UKwo50Qsco4MFJ45Dc9hP8uuUtHD+ejNYJPXH/Wx+i19AxUofpE5xgMb8gIhTrzSg1mGW1JAgAOZkZmHnfaJhNJny851t065EodUgAAFVoCDo1i0BYCK/cs4ZLc5JcAfYtaP4oIwxo4zrBSlUDpdyqgPnIiWLCV9mEXs2BsQkCQhUCEoeMRu+RD0odml/wSM58zioSsjUGlMgwuUq7cB4P33MnBEHAlr3ySa5iwkMRH6vi5Ir5hNZMKDDU/fjpUveepyk1FmX+Q0T4PlfEvmxCv9bA+K725Kqx49Gc+ZTBakNGhR46GfW3cvj95Ak8MmokmrdsiY/3fIuOXeKlDgkA0EoVjo7NIhpF3xcmD2ma+h/P05O9yLgeBXpCehNqLMr8w2wjbL9E+LEQGN5RwF2dgr+BqLs4wWI+U260IFOth1mUVzE7APx0+CCmj70XXa9IxOYvv0Lrtm2lDgkCgPbREWjDzUOZj6XWsTxY3Zmyuo8hInydzckVaxiNmbAxhXBRCzzYTcAtbYUmNdZxgsUaTCRCns6IfJ0RXrTY8btvd+/CYw/dh34334J1/92FmNjmUoeEEEFA5xgVWvCVggyoukzdF0QiXHJj5un3UqrzdU+VApk6Gb6ZWdDIqySsSyborcD0RAE9mzedxMqBEyzWIGabiAy1HhUy2qy5uu2bN2DBzCm4c9QYrNqyDZFRUVKHhDCFgPhYFaLD+RoTZk+uvsvxXTKTrQOMbqzQayxAurb2/ToL4bsc+c1Cs+DxR5l95io2DHj0SgHtIptecgVwgsUaQGu2Ir1CD6NVfoMxEeGDt5dj6bPzMXHGLCz/YB3CwqSfLYoIUSCBO7Ozas6V24vJS4y+SbLSNO4/z5nS2sd+nU1uJWiMXc5sI+zOFLEjw36l4JREAdFhTTO5ArhNA/OCnFswAPZNm19/4Tls2/gR5r24GI/OXyCLdf+osBB0asad2dnfRCIczLN/QPm1mHBX54b/brhTf+VwvoLwDxshPMT+uikVhHNNbM9B5hsFent/K7UFGNVFQJ+WkMW4KyVOsJhHrKJ9L8FKGV4lCAAmoxGLHn8UP3y1B0tXrMJ9D0+ROiQAQKwyFB2iI5r8gMNqOlP6d5+p06WE2zv8nex4w1V7hsuZReDPCqBPS/vsw75s+c1GM3kjIhwvBvbnElpFALOuFNAqgsc5gBMs5gG9xYZcrQEWOVayA9CoKzD3kYn4/eQJrNz8KYaNvFvqkAAAcRFhaMtXCrLL2ETC4fy/ExqTzX5lX//W3v+euGrP4MyZUkKflgJ+yCOozV6/NGuCtBbC3kxCigYY0NrehqEp9LdyFydYzC2lejOK9CZZLgkCQGF+Hh57cDwK8nLx0Re7cf2NN0kdEgCgTaQSrSLDpQ6DydDJUkLFZQnNr8WE/q29f05P6q8cMrSEc2WE49xUlLmJiHC23F6vFyLYWzA0xasEXeEEi9XLJtpbMGjNVqlDqVPahfOY/cA4QBCwZe93uKLnlVKHBAFAu+gIbsPQhNlEgkJwXodiFQlH8msnNEUGe8KT0MzzP1YiARe9SLAIwI4MUbYfnpi8aC2EfVmEZDVwTQvgrs6CW/taNkWcYLE6Gaw25GqNMNvkW5fx64/HMOeRh9C+c2d8+NkXaNu+g9QhQSEAHaIjEKPk5Kop+73MXld1TxegtarmH6BfiwnaOjqbHC/yLsHKMwheX/0n01V/JiOXz1rd31VArxacWNWH2zQwp8oMZntXdhknV1/t/AKP3j8GV/Xug493fy2T5EpA52YqTq4YzpQRMnWED8+LOJgnwvpXFmO2EY4V1J3RXFATNGbPM55LOh7OmX8YrPYrBHdmEK6IAZ64ipMrd/AMFqvBJhIKKo1Qm+S7JEhE2PD+SqxYtkMIgPoAACAASURBVBij7n8QS//9PsLDpa9zCv2rO7sqjHtcNVYmG4EIiHCxJKI2EzL/6qZuI+Bwvr39wT1dFMjSESrreXuJBPxWQhjawbM/YJd0CvvaNGM+lKklJGUQTCIwPkHANXFN95fM0++cEyxWxfjXkqBJxrNWVqsVr7/4HD7fuB6PLViIpxa9JIur88IU9uSKG4g2PsUGQqqGkKoGsnSEfq0Fl/2qfi+lWjVNJUZgU4oId8pVfishDG5HbvdM01kIRUYBESq3DmdNxPkDO3Hso9dQkpmMVvE9MWjmy+g1bKxb59qIcCiPcKwQ6BINjEsQEBsu/VgrJaWHwzvPKTdSSUlJ6N27N1QqFXr37o2kpKR6jy83WpCh1ss6uarUafHU5AfwxcebsOzf72PO8y/LIrkKD1EgnruzNzomG2FNahj+86eI73II6VqCjYBTJQSDtf4lvN/r2UjZxakAAJ3F3p/KXWkacJE6q+H8gZ3YvnACuptysWhAPLqbcrF94QScP7DT5bnlJsLGZMKPhcDQ9gKm9ODkCgCiQj17l3GC1QglJSVh/PjxOHv2LIxGI86ePYvx48c7TbJEIuRqDbLdqNmhIC8XD98zEqeO/4IPt+3A+MmPSB0SAEAZokB8rArhIfxWamxsBFRYav9RMYvAieK63yx5lYRiY8Nf/8cCEZlagljPRtAVJsKxAhGH8uT7wYhJ49hHr+HW+FbYOvp6zLo+AVtHX49BXVrhx/Wv1XteTiXhI8cmzT0FDG4vQCGDD7JSUwhAlIefoXmJsBFaunQpBEEA/TUwExEEQcCyZcswbty4quOCYUkQAP48cxpPTn4AipAQbNn7LRKvulrqkAAAqtAQdI5RcWO9JuiXYsLNbcnp//v6Zq88UWAANqaIiAgBujYT0CNWwBUxQIgAnCsnnC2z/zGU8eciJqGSzGTMGBBfNcsvCAKGxrfE8ePJdZ7zZ7m9kL1DJPBAd26/UF1kKOBpnskJViOUkpJSlVw5EBGSk/9+Y1UYLSiolPesFQAc+HofFj42A917Xon3t3yO1m3bSh0SACAyLASdeV/BJktnsbdhuL5VzftFIpz1UYLlYLTZ9ww8X2F/XgG8HMhcaxXfE4eycvFo3/iqD9wHM0vROqFnrWOJCP8rAr7PJVzTAhgdzx3ZLxftRbbE6xqNUGJiYq3aJEEQ0LNnT4hEyNMakSfzJUEiwqYPVmHulIkYNGw4Nu36SjbJVXRYCLrEcHLV1P1UKNb6IHNRg3qvEPQFGb9tmYwMmvkyjmaWYOKuk1h7MgMTd53EsawSDJz5So3jRCLsyyZ8n0u4tZ29mJ2Tq9qaeVGDxglWI7RkyZKqZUEAVZ9eXnz5FaRX6FFhqqPDoUyYzWYseWYu3lryEmbMeRor1m+GKjJS6rAAAM3CQ9EpRsU1CQzFRiBVXfO+M6Wc/jB56DVsLO5/879IV3XE/x3PRLqqI+5/6wv0Gjqm6hiLSNh2kXCqBBjVRcCwDgpZXDgkR97MYPESYSM0btw47NixA8uWLUNycjJ69uyJZ194CdcNvVP29VYVZaV4esYUnPzlJ/xr1WqMeXCS1CFViVWGokN0BA9ArMr/igiJf+3BZrIRLlRwgsXko9ewsXW2ZTDZCJ9dJOTpgYlXCOgew+NafaK96B3tMsFKTU3FlClTUFJSgtjYWGzatAlXXy2PImNWt3HjxmHcuHEQiVBQaUKF0VLv1UhycDElGU9Ouh9ajRrrv9iNfrcMlDqkKs1CBU6uWC0ZWkJeJaFDlIA/y91rwcCY1AxWwtY0QokRmHyFgC7RPK654k2C5XKJcPbs2Zg1axZSUlKwaNEiTJ061YvQmBRMVhEZFXpUGOW9JLh/726MuOFajBrYH4X5eZjz/MuySq7iIsLQMpynzplzPxbas6ozPi5uZ8wfdBbC5lRCmQmYksjJlbuahfm4BquoqAgnTpzA5MmTAQDjx49HdnY20tLSvIuQBYzaZEG6Wg+jzJcEv9vzJeZPm4zcrEwAgMVsxj8XPoP9e3dLHJldK1U42vHMFavH+Qr7tjiOrXEYkyu1mbApxb5V09REAe0jeVxzl89nsLKzs9G+fXuEhtpXEgVBQJcuXZCVleX0eJ1OB41GU/VlMpk8j4g1iEiEfJ0RuVqj7JcETUYjljw9p8Z9juL81W8vlyiqv7WJVKJNlFLqMJjMiQRsvyTy1X1M1ipM9uTKSsC0RAFtVJxceULyIvchQ4bUuL1w4UIsWrTIly/hc+Xl5VKH4BVncVtEQpFJhFnG/ReMBnuL65LCQjw3exo06tr7gRAR0tNSYTAYAh0eAHufobhwBRRGC8r+6sjdmH5PgkFcXJzUIXjE360ZGGuICpN9WVAQgKm87Y1XosMAnYfn1Jtgde7cGfn5+bBarQgNDQURISsrC126dHF6/OHDh9GnT5+q20qlEkql/GcAgm0wd6get9pkQbHOhBAlQe77vaaeP4d5U/6/vTsPj7K8Gj/+vZ/ZJ5PJvkFIQoAEkB0XBBHQ0rpVIKlalSJU1Eqt2loBW1vUWqUVebX2VUEtSrH+pCVYX5e2Kpv7hrZIkUXZRJAtJJA9mef3x5AQQiaZfeZJzue6ckEgeeaQDM+cnPvc5/buDszrXcjuHdtPmieklKJ33yIcjuj/SxSQ47KTbD+1HtwVnidCiO6lol5n6Vbv/VXOFAyOVQOrKcw9WJmZmYwYMYJly5YBsGLFCnJzc+nbt2+7H+9yuXC73S1vRkiujK71kmBTnC8JArz01+e55tILycnNZflra/nZr+9pd2bXrNvnRj02BfRIbD+5EkIIo6ms91auPEhyFYrEIF8SOt1FuGjRIhYtWkRRURHz589nyZIlwT2SCLu6Jg87Kqopj/NdggANDQ3c/4vZ3HXbzVxSerl3Mnt2NhMvuZSHliyjaOBpWG02igaexsNPP8u3Lv5uVOPTFOS6HSTZJLkSQhjf0eO7BT26N7lKtklyFSxXEDsIwY8erOLiYt59992gLi4i51ijhwNHqg1RtTp88CA/m3kNn7z/LnPunc8Prr/xpF15Ey+5lImXXBqz+DSl6JVoJ8Eqc3eFEMZ3rEFn6RadRo93FEOKJFchCWYHIcgkd8Px6DrfVNVxoM6DwxH/ydVnn67nlulTaaiv409lLzFw2PC4GnlgUopebgdOiynWoQghRMiqGmH5dp06j7dylSrJVciCTbDkLEIDqWtsYscRYywJAqx87ll+cMl3yMjKYvnr6xh59uiwP8ZrL73IlHGjGZ6bwZRxowOan2XWFPlJklwJIbqG6kad5bssVDfCtH6KNLskV+EQsR4sER+O1DawvaIm7geHgvew5ntm/5Q7b76R737vCpa++A+ye/QM++O89tKL3DpjKls3baS+ro6tmzZy64ypfiVZFk2Rn+TEbpbkSghhfDWNOn/eqlPVqLimnyJdkquwCbYHSxKsOOfRdfYcreHrY9EZHBpKRQjgm71fM33ShZQ9+2fuWvgH7nnoj1gjtJv00Qfmt+w6BP+HlNpMGgVJTmwmefoLIYyvtkln2Tadinq4PK+BDBkiGlbBDBkFSbDiWk1jE9uPVFNRF50phqFUhAA+ePtNLjt/LN98/TVL/+8fXPaD6RGNd8cXW0+anwUnhpT6Yjdr5Cc5sEhyJYToAuqavAc3H66DH/RTZNrjvzfXaKQHqwvRdZ1DNfXsOFJNXRSXBIOtCOm6zp/++DAzSy+lb/+BLH99HUNGnB7xeAv69DulYb55SGl7nBYT+W4nZk2e9kII46s9nlwdqIGpfeVswUiRBKuLaPR42H20lm+q6qJ+tlkwFaGjlRXcOuMHPHj3r5jx41tYvHwlaRkZkQ4VgFm3z/V7SGmi1Uye24FJkxuQEML4qo6PYjhQC1P7KXomyL0tEjQFCbJEaHzH6hv58kg1x+pjc7BZoBWhLf/dyBUTx/PeujX84Zm/8NNf3dVyMHg0+DukNMlmITfRjhZH4yGEECJYlfU6T2/VqWzwjmLIleQqYhLMBD1aSOZgxQGPrrO/uo7ymoaoV61am3X7XG6dMbWlEtRRRejF5c9x989vJa+wD8tfX0t+YZ8YRNz5kNI0h5VMpzWuZm8JIUSwyuu8Zwt6dJhRJKMYIi3Y5UGQClbM1R5vZD8c4+QK/KsI1dXWctdtt3DHj2/ggkkl/OWV12OWXHUmy2kjK8EmyZUQoks4UKOzZIuOpuCHxZJcRUOwIxpAKlgxo+s6h2sa2F8d/V6rjnRUEdq9Yzs/u/Yatm3exD3/80dKp06LcnT+UUCOq+NDm8vKyrj77rvZsmULRUVFzJs3j5KSkugFKYQQAdhxVGf5lzqJFu9uwVBe+IX/gh0yClLBion6Jg87K2v4Js6Sq4688cpLXHb+uRytPMJfXnk9bpMrTUEvt6PT5Kq0tJQNGzZQW1vLhg0bKC0tpaysLIqRCiFE53Rd5/393mXBbCdML5LkKpqCnYEFkmBFla7rlNfW8+WRaqobmmIdjl8aGhpYcNed3HzNVZw19lz++sabDBgyNNZhtcusFHluJ65ODm2+++672x1Hcc8990QjTCGE8EujR+fvO3X+8ZXOWZneUQwOsyRX0SQ9WAZQ3+RhV2UNe4/V+T2RPdSp6qHa9/Uepk+6iD8vepTb77mPh5YsI9GdFNUY/GU1aeQnO/06V3DLli3tjqPYvHlzpMITQoiAVNR7+602lsOUAsV3cjXZCR0DoVQLJcGKMG+vlbdqVRVA1SrUqeqhemvV65ROGMO+PV/x9N9fZfqNN8Vts7jdrFGQ5PD76JuioqJ2x1EUFxdHIjwhhAjItkqdJz7XqWr0NrMPSY3Pe293IBWsOFXX2MTOihr2VflftWoW7FT1UDU2NvLwffdwwxUlDBo2gr+teovhZ57l1+fGouLmspgoSApsOvu8efPaHVA6b968SIUphBCdqm3SeXGnh2e36WQ74Lpimc4ea9LkHmc8us7+qjpvr1VjcL1WwUxVD9WBffuY+b1JPPnwQm795Twee+5vpKSl+fW5sai4Jdst9HI7Ai6bl5SUsGLFCoYMGYLdbmfIkCGUlZUxZcqUCEUqhAjEplUreeKqkdw/xsUTV41k06qVsQ4p4rZV6jz2X++S4CV5iqv7KhKkmT3mpIIVR5qnsR+sqQ9ph2CgU9VD9d66NZSeN4Yd27byp7KXuO7W29ACqApFu+KW6bTSw2UPetmypKSETz/9lJqaGj799FNJroSIE5tWrWT57MvoU7eHOWfm06duD8tnX9Zlk6zWVat0O8waqBiZruK2JaM7sZnAEsLxapJghUl9k4fdlTXsqqyhPgwHNAdyzl4ompqaePSB+5n5vUkUDRzEitVvc8aYcwK+TrQqbgro6bKT7rSF9brxrqysjKFDh+JwOBg6dKiMlBBd1ltP3svY/HSenTSC60cU8OykEZyTl87bT90b69DCStd1Pjmo88eNJ6pWU/sqkqySWMWLUEY0gCRYIWvy6HxTVceXR6o4GsYzBP09Zy8UB/fv5/rLp/DoA/OZNfsOFj1fFvRBzdGouJmUIi/JQVIHM666IpnbJbqTgzs3Mz4v7aQfLifkp3FgR9fZ5bvrmM4Tm3Ve3KXTO1GqVvEqMcRkVxKsIDXvDvyivIpDNfV4IjAxdOIll1K25h0++eoAZWveCWty9cHbb1I6YQxbN/2XJ1e8yKyfz8Vk6nzEgS+RrrhZTRoFSU4SLN3v8AGZ2xW4rVu3Mnr0aIqKijjjjDPYuHFjrEMSfkrPL2bNrkMnPd9X7zxERoHxd/lW1Ous2O5hyRbvv21GkaK0tyZVqzglFawo8w4LbWBbeRX7qupoDHB3YKw1NTXx+IO/59qS71JYVMyK1W8zauy4kK8byYqb8/hOQZu5ez5dZW5X4G644Qauv/56tmzZwpw5c5g+fXqsQxJ+Omfmnby58yBXvbCexet3cNUL63lr10HGzPxVrEML2tF6nVd3e3hko86OozApX3FdsSLPJYlVPAulwR3kLEK/eXSdiroGDtU0hKXHKhYOHTjA3FkzeXftGm78+Rx+dNuckKpWbXV0jmGwkmxmclz2bj1gr6ioiA0bNpyUZMncLt/279/PRx99xL/+9S8ASktLuemmm9i2bRt9+/aNcXSiMwPOm8Llv/8rbz91Lx98sJmMgmIuf+BxBkyYHOvQTrJp1UreevJeDu7cTHp+MefMvJMB5528WeZovc5b3+h8fBAsGozNVozKBJup+97PokEBShHyypIkWBHW6NE5UtvA4dp6GiOxDhglH7z9JrNvuBaPp4kn//Z3Rp07PtYhdUgB6U4rGd2smb098+bNo7S0tGWZUOZ2dWz37t3k5ORgNntvb0op8vLy2LVrV7sJ1rFjx6isrGx532azYbPJ8y6WBpw35ZRkJZ4073Qcm5/OtWfms2aXd6fj5b//KwPOm0JFvc47bRKrszLBLolVVCRZwazBwdrQrhPqmY9hTbCqGxpP6sMxspqGJsprG6isb4hIf1W0eDweFv/PAv739/cx8uwxPPD4U2RkZ8c6rA5pCnJcdpJs3auZ3ZfmuV333HMPmzdvpri4mHnz5sloiTAZN+7kJfLZs2czZ84cqhuhpsY41eqamhBfTaLMSPG2jfXNJ+5p2emolOK64flc9cJ6Vi/+DRvyLuHzSg2rBmenNTEytQmbCfR6qIlBrPEuEvFmaB4cJqipCa2tpKmqgcPqRAJQXl7e4cenpqae9H5YE6x9x+r44kg1yTYzSTYLFj+PLokXjR4PlXWNlNc2UGfQZcDWvEuC1/Hu2tX86LbZ3BhiI3s0WDRFbqIDhx9nCnYnJSUllJSUxDoMQ+jVqxd79+6lsbERs9mMruvs2rWLvLy8dj9+7dq1DBs2rOX95gqWvVHH4ajC4XBEK/SQGSlWMFa8rWM9tGsrM8/MP2Wn4zvvb6Gh1sR3chXD08Aao/utkb6uEP5481IVbivsagitOpKbnkBqm0n6bZOojoR9ibC+ycP+6noOVNfjtJhItJpx28wBHWUSTY0enfLaeirrGqluaAppOGg8+eidt7n9hh/S2NjA4uUrGT3+vFiH1CmH2URuot1wibmIL5mZmYwYMYJly5Yxffp0VqxYQW5urs/+K5fLhdvtjnKUwsi8Ox33cN3w/JYl+1U7D5GSV8ys01S37hmNByk2SLcrCPEVPW57sHSgqqGJqoYmvqmqw2Ex4bKYSLCasZu0mC0j6rpOdWMT1Q1NHKtvorymCTt1MYklEjweD0/94X/4w/2/YeSo0fx+0VNkZufEOqxOuW1menTzZnYRPosWLWL69Oncd999uN1ulixZEuuQRBfg0XW2VQKX/II3F17BlS+s57z8NFbvPMTbuw5y+QOPyz0sDqTaFFkhFsU0Bc4QM6SoNLnrQHWDN6mhuh6TUjgtJhxmEw6zht1swhTCOHqfj6vrNHh0ahubqGn0HP+16aSeqmhUrF576UUefWA+O77YSkGffsy6fW7Yd9sBHD54kLmzruOdNau47tbb+PHsX7Q0+sYraWYXkVBcXMy7774b6zBEF6DrOl9VK7Yd9LDpCBxtgOzhUxhx53K+XP5bPgxwp6M/uw9FaFJtkGwFqwb1QXb7uMyEXAiKyatvk65ztL7xpMnnFk1hNWneN03DrCnMmsKkecutmgKFovnfq+venyY8eKepN3l0GnWdhiYPDR4P9U06dU0ePDGeU9V8CHJzGbn5EOSHliwLa5L18Xvvcvv1M2hoqGfR82WMmXB+2K4dKZpS9HDZcEszuxAiDMKVvHh0nd1V8N9y/XhSZSXRAgOSYWiqIscJakAJTA6sL7Kz3YciPFJs3uQo06H4qiq4HCDUHYQQR2MaGjw6DR7vkmJX0tEhyOFIsDweD0898hCP3P8bhp85igcW/8kQS4JWk0Zuoh27WZrZhRChCzV5OdbgXf7bVqnzRSXUNkGiBQYmQx9nPX1TbSFXNFqfs9h69+HbT90rCVaYJJhPzBnLcsBXVcFdJ9T+K4ijBKuriuQhyIcPHuQXN93AW6teN8ySIECCxUTPRAfmCCwLCyG6p0CSF13XOVIPe6pgT7V3uvq+4zMUejjhzAzo61bkJngrITU14Rk/dHDnZq5tZ/fhBx/IqQzhkmI78X0KpQ9LEiwDKOjTj62bNp4yhTvUQ5A/ef897rz5Rurr63j8/63gnPO+FWqoUeG2aOS5HV1iVpoQIn50lLwcqNHZXwv7a3T2VsOeaqg+3qGSbIVeCXB2lqJPIiSEYWnIl/Z2H3aVcxbjRWqrdt4sR/A7CUM9hxAkwYq4WbfPPakHK9RDkJt3CT4y/16GnzmK3y96iqycHmGOOvw0BTkJdpqqGyS5EkKEXVp+EWt2fX3K6ISmrGIe3eR9kU0wQ44TTk+H3ARFD2dkE6q2zpl5J8tnX8ZVL6xnwvHdh28d330owqN1gpUZQgUr0Rp6LJJgRVjzIciPLZjP9m1b6d23H7NuvyOoQ5BbDw6dcdMt3PKLXxtiSdCiKXLdDhxmE4erYx2NEMLIdF3ncB18XQ1fV3srUvtroOa8X7L38ROjE1YdH50w/M7HGNJPkWGPbjLVnmDPWZSdh/5rvUToMCvcFqhsCPw6LnMXanLvysJxCPL7b61jzo9m4vE0sXj5SoafdbYhkivptxJChKq8TmdLBWyt0PmqCuqOb71Ptnp7pkZlKtJKS6jouZzPlgU+OiGaAj1nUXYeBia1zcSfLKeisiLwZcJu24MVrblS8aCpqYnHFszn8Qd/zxljxvL7x54kIzubmpponGoVmjSHlUynVZYEhRAB0XWdPdXw+RFvYnWg1ttmUOCCMdnepb0cJzjbVhkuLGH0hV3rSCnZeRiYlLYJlgO2VgR+nW6ZYEVrrlQ8+Gbv18z50Uw+fu8dZs2+gxt+envcnyUIMt9KCOEV6NJWbZPOfw7BRwd1DtR6J2n3c8P4HEUf94nt992J7Dz0n1U7dX5VsI3u4UiwDHfoW0dzpbqSNf96lZLxo9n55Rf8qewlZhngoGYAm0mjd5JDkishurnmpa0+dXuYc2Y+feq8S1ubVq085WP3Vev83y4PCzfo/OMrnQw7/KCv4rbBiskFGgNTVMSSq02rVvLEVSO5f4yLJ64a2W58seTdeXjopNc82XnYvrbVK4BMe+DXcVm8vcOhMlyCFcm5UvGgvq6O+XfO5cdXX8HQ08+kbM07nDHmnFiH5Re31UxBkhObDA8VottrvbR1/YgCnp00gnPy0nn7qXtbPmZPlc6yrR4Wfa6ztQLGZCl+OkhxWaFGoTvyhyYHkgTGyjkz7+TNnQe56oX1LF6/g6teWM9buw4yZuavYh1a3Em1nfp8SbdDoLn5aSnhed4ZLsEq6NPvlJ6ecMyVigfbt23lygvO5/8teZI7fvs7/nfZ86SkpcU6rE4pIMtpI9ftiMiZkkII4zm4czPj89JOWdo6sGMz+6p1Vuw28+RmncoG+F5vxS2DFONyFInW6N1D/EkCY6155+F2R09+98FOtjt6cvkDf4u75v14qAS2bXAHMGmK9ACrWENSI5hg1dbWMnnyZIqKihg6dCgTJ05k27ZtYXnAUM26fW7LsiAQ8lypeKDrOmXP/pnLzh9LXW0Nz/1jFVOvv9EQzeFmTZGX5CDNGYahIUKILsPX0palRzGLPtc5VKcoKVD8aIDitBSFKQb3u46SwHgy4LwpzHz2Y+54+xgzn/04LpOreKgEtrdECM19WP5Jt0PPhAhXsK6//no2b97Mv//9byZNmsTMmTPD8oChap4rVTTwNKw2G0UDT+Php58Naq5UPDhSfpifXXsNv7r1x1w4uZTlr69jwOAhsQ7LL06LicJkJwkWw+2VEEJEmK+lLf2SX3JpnmJmnwYGp0Z+GbAj0t8UHvFSCWxviRACGzgaruoV+NhFaLfbueiii1reHzVqFAsWLAjbg4YqHHOl4sEHb7/JHbOup7qqioVPPsN3Jhlny62MYBBCdGTAeVMo/d1y3lj8W955fzOm7GKG/vIxLpk0GbOmiIdJMzJZPTziZadjxxWszncSKqKQYLX18MMPM2nSpE4/rr6hnrq6uhMXN5kwGWAYZrTV19Xxh/n38vT//oHTzx7D/Y8uJqdnbqzD8otJKXq47CTa5PsqhGifrut8XgGrsidTMWcyI9O9oxZiPUm9rWAnq4uTxcMZiyYFST46Vfw99Dk/UZHsowoWjE5fJe+77z62bdvGG2+80enFnl7yNCmZWS3vjxkzhjFjxoQWYYTV1tRG9fG2bd7Er2/5MV9u3cJP7vgVV1/3I0wmU8CDQ6MdN4BNU6TZNBqq6jlcFdw1ysvLwxtUlEjc0ZWamhrrEESQDtR4Ry18eRT6uuHKQkVGAD0w0RboZHVxqnBXAl0WSLcrdhz1f35VshWfy81uq8Jhgpqmjq8RzuoVtEqwli5dysKFCwG45ZZbmDFjBgsWLKCsrIzXX38dp9PZ6cWmz5jOgMFDT1zcIBUshyOEEyH91NTUxDOP/ZE/3P8b8gv78Py/1tB/0OCQrhmNuJul2C1kJdjC0i9h1BdPiVsI3+qadNbt1XlvPyTZ4Mo+iqKk+E2sRPiEuxJY4FL0csGOo/5/TkonlacsZ8cJm1nBaSn+P54/WrKfadOmMW3atJa/WLhwIc899xyvv/46ycnJfl3MarFis/lYBO3Gdm3/kl/+5EY++eA9ps/6CT+Zeyc2exDTz2JAU4ocl40kGRwqRLfV0UR2Xdf5rBxe26NT0wjjchSjs5DzR7uZcFYC8xOhKEnx6m7/K1jtjWhoLdPeccLWPzn8w2zbLS999dVX3HbbbRQWFjJhwgQAbDYb77//flgfvKvzeDw899Ri/ufeu0jLyOSZv7/KyLNHxzosv9lMGrmJDmxmw41LE0KESUeHDaedPZlXduvsPAYDkuHbPcPbwyK6pwKXIsmq6OFUfF3t/7/y9gAAIABJREFUX5Llq8G9WWeN7uFeHgQfCVZubu4p09K7gvYOiT7n/IkReaydX2zj1z/9CR+9+zZXXXs9t955FwkuV0QeKxKSbRayXeFZEhRCGJevw4b/8fi9HE2eRKoNpvZV9HHLvUKEzmWhpWdvQDJ8Xe3f5/ka0dCso0b3BDP0TfI3Qv91m9JE8yHRWzdtpL6uruWQ6FWvvhzWx2lsbOSpRx5iyvjRfLN3D0teeJlfzl9gmORKU9DDZadHol2SKyGEz2Gcx77azIQe3kGhklyJcMl3nXgu9U/2/3nVWQUr0+Edw9CeQRGax9ZtEixfh0Q/8dCDYXuMzz5dz5XfOY+H7r2LK394HSvXvseZY8aG7fqRZjNpFCQ5SbZLv5UQwsvXMM7MgmLGZivptRJhVZB44vcZDv+OuVF03oNlNSmSfXzM0AgsD4Kfc7C6Al+HRO/8MvQjgI4dreSR++/lL08tpt+A0/jLP95g8PCRIV83mpJsZnJcUrUSQpxs1A9/Sdncy7nyhfWcJ8M4RYS1rmABDEhWvLmv45alRIt/myrG5WjsOqZTUQ+V9d5fk6zQI0xH47TVbRKsgj792Lpp40lJllKKgj59g76mruu8XPZXFsy7k2NHK7ntrnuZet2PMBtgNEUzTUFWgp0UqVoJERMd7dCLJV3X2XQEXs+ajHbj83z2r/v4UIZxighKMENmm5lp/f1IsDob0dBsWJpiWNrJH9voiVy/uXEygRDNun0ut86Y2rJM2Pzrdbf+PKjrbd74Gff/YjYfvvMWEy+ZxOx7fkuPXnlhjjqybCaNnol27GZTrEMRolvqaIdeLJOsQ7U6r36l80UlFCXBBVNLSLm2NGbxiO6hbfUKoIcT3BaobPD9eZ0tD3Ykkkvc3aYHy9ch0RMuuKjzT27l4P79zPvZzXzvvHM4uP8bFi9fyUNL/my45CrJZqZ3slOSKyFiKF4OyW1W16Tz2h4Pj27SOVQL3y9UXNlH87tCIEQoWvdfNVNKddrs3lmDe6x0mwoWtH9ItL9H1FQdO8afFz3Kn/74MCaziTm/uZ8rZszEYjHW0posCQoRP+LlkFxd19lwfFhobSOcm+0dFmqRBnYRRQWJ7T/fBiQrPjjgeymvsxENsdKtEqxg1NfV8bdlz/D4g7+nsuIIV/7wOm742e0kpxjv+BFZEhQivsTDIbm7jum8sUdnVxUMTIaJMixURECeSzEyXbFyh6fdv3eaIcPHjsH8RDo8SzCUJcJIkgTLh/q6Osr+8meeeOhBvtn7NZdefiU3zfmF4ZYCm8kuQSHiT7gPyfWXruvsOAbr9np/zbTDD/oqCmWelYiAYWmK7+YpTJris3LF1opTq1H5LtVSyW1LU4riZMWnh9qvYskSoUFUV1Xxtz8/zTOP/ZFv9n7NhVNKufHncynsVxTr0IIiS4JCxK9wH5LbGV33Nq6v26ezuwqyHXB5oaJ/Ej5f3IQIlgIm5ipGZ51o9/5OruLLSp2mNrlSe/1XrfX3kWA5TOAwx+dzVxKs4w4dOMDzTz/Js08u4mhFBReXXsa1N/+MvsX9Yx1a0GRJUIj4F85Dcn05XOs9kPmzcp0DtdDTCVf2UfRzS2Il/OeywLEOdvO1ZjNBaYFGUZsG9XS74owMxXv7T06WCtrZQdhaHzdYNGhos8IYzxswun2CtW3zJv629GleXP4cmqYx5cqpzPjxzfTMy491aCGRJUEhurfKep2Nx5Oqr6u9L079k+CCXEXvREmshP/MCs7vqTgzQ/HKbp2PD3Y8OyrNBlf00U6ZadVsfI7iP4d1qhu97zvN3qNsOmLRFONzFBX14DCfqFylxenyIMQwwWrv4OW2O/wipaGhgVWvvsxzTy3mw3feIiMrm1k/n8tl18wwZPN6a7IkKET3pOuwp0pnS4XOlgrYVwMmBf3cMLq3oihJdgWKwGXYobS3RrbT+9z5br4i2eph1dc67aVZQ1MVF+cprCbfzzW7WXFeD8VLu7xX6Kj/qrUx2caaLBWTBKv54OXmXTPNBy8/tGRZRJOs3Tu2s2LZUsr+8mcOHdjPyFGjuf9/F3PhlFLDjVtojywJCtF9bFq1kjef/A0Hd27B1qOI+gt/Sf2wKdhNx5OqLEW/JLB38EInREfOzFBMzFWnJOZjczSSbTov7PC09FJZNbgoTztlUrovI9IVHx7Q2VED+a5wRx4fYpJg+Tp4+bEF88OeYFUdO8a//u/v/P3/PcuH77xFojuJSy//Pt/7wXSKBp5GTU1Nl0iuZElQiK5P1709VG+/upL/3Hc55+SlM/PMfNbs/Jo3H72CCfcs55wLp8h9QITEboIpBRrFHQz4HJyqSLRoPP+FB7cVLivUSLf7/7zTlOKCXI3HD/uef2V0MUmwfB28vH3b1rBcv6GhgXdWv8FLK5az+h+vUFtTw1ljxzH/0cV86+JLcTidYXmceCBLgkJ0bY0ene1HYUuFztZKqKgH87J7GZuXzrOTR6CU4rrh+Vz1wnq2PPtbzr2oJNYhCwNLtsLVfTUyfPRPtVaQqLh+gOb3Yctt9XYrBid7yOqk/8qoYpJg+Tp4uXff4EchNDQ08OHbb/KPv6/k9ZdfpKK8nL79B3DDT2/n4tLLDDu/qiOyJChE11TVoLO5wptUfXnUu3Mq2QrFSdAvSfHXb7YwPg4mwIuupWeC4so+CpfF/2Qp1F1838lp7LIbLmKSYPk6eHnW7XMDuk5NdTXvrlvDGy//H6v/+QoV5eX0KujN5dN+yAWTSyg+bVCX/ca5zIreyU5ZChCii6io19l0BD4/orPrmPfPchO8x9YUJ0G6/cTOv3iYAC+6lgHJipLep/ZbRVpXbhGMSYLVfPDyYwvms33bVnr37ces2+/gWxd/t9PP3b9vL+te+xdr/vUq765dTW1NDYVFxVwx/VomXnwpA4YM7bJJFZxYEtSrGyS5EsLgKuq986k2luvsrfa+2BQmwiV53qQqwUclIVYT4EXXNDpLMbGnfzv5hP9iNqahvYOX21NfX8+nH77PO6tX8eYbr/H5Z/9B0zSGnXEWP579C8Z/50LDTlkPVOslwcPVsY5GCBGMqgad/x6Bzw57z/8zK+iXBKMzvbv+bH78SN/eBPjJv304YhPgRdd1Qa5iVJaxxh8YRdwNGvV4PGz570bee3Mt769bw0fvvkN11TFS0tIYPf48Ztx0M+dMOJ/k1LRYhxpVyTYL2S6bVK2EMKAGj87mI/CfwzrbKr1HiBS6YXK+on+yf0lVW20nwNfU1IQxYtHVKeCSPI2RGfKaEikxT7AaGhrYtOHffPLBe3z0ztt8/N47VJSXY7PbGX7mKK6/9TbGTDif/oOHoGndL8vWlCI7wUay7BIUwlB03Xve36eHdP5bDnUe6JUAF/VSDEwBZ5yenya6Pk3B5HyNIX7OrBLBiXqCtX/fXjas/5h/f/Qh//n4QzZ88jG1NTXY7HaGnn4mV8+8gdNHn8Ow08/EZrdHO7y4Yjdr9HQ5sJm7X2IphBFtWrWSdU94h39q2UXUX/RLkkZN4axMGJKqSAtgTpAQkWBS8L3eGgNS5LkYaRFLsDweD3t27WTzxs/4/LP/8N///Jv//vtTDnyzD4CsnB4MGXkGP5l7JyPOGkX/wUOxWq2RCsdwUuwWshJkSVAII2jy6Kx6aSXv3OMd/nndWceHfz5+Bd/uvZyBg2Q2lYg9iwZXFGr0TZLXlWgIa4L19+XP8dyfnmDr5/9l2+efU1NdBUBaRgbFpw1mylVTOW3ocAYNH0F2j57hfOguw6QUOS4bbpssCQoR7/ZV6/z7sM5/DkPdkvaHf77zp98y8HxJsERsWTS4qo9Gb7ckV9ES1gTrb0ufpk///vQtHsB3Lp1C3/4D6D9oCBlZWeF8mC7LaTHR02XHYpIlQSHiVVWDzoZy+PchnX014DTD4FRYv1+Gf4r4ZFZwpSRXURfWBGvpS//gtKHDw3nJbkEBaQ4rGU6rzCERIg5V1Hsnq28+orPjKKC8U9XH5yj6Jnkrz7tl+KeIQ2YF3++jUSjJVdSFNcFSyDcwUBZN0TPRjtMS8w2dQohWdh/TefuAiS+rPOytAQ0oSITv9FIMamcXoAz/FPHGpOCKPtJzFSvyqh5DbquZHJcdU5SPJhBCdO7u9R4+PGyiXxKcnaXo5wZ7B6MV2hv+efkDj8vwTxETJgWXF2r0k+QqZiTBigFNKbISbKTIbCsh4tbisRo9LLW4nA6/P6ft8E8hYkE7PoqhOFmSq1iSBCvKHGYTPRLt2KSRXYi4pinVpQ+iFV2TpqC0QOZcxQNJsKJEAelOK+kOaWQXQggRfpqCi3s0clqqvMbEA0mwosBq0ujhsuO0mGIdihBCiC5IAZPyNXopT6xDEcfJOlWEpdgtFCY7JbkSQggREQq4NF9jqJwtGFekghUhFk2R47LjssqXWAghRGQo4JI8jeHpklzFG3n1j4Akm4XsBJuMXxBCCBFRGQ4YmSGvNfFIEqwwsmiKbJedRKlaCSGEiIIeTkmu4pVkAmGSbLOQJVUrIYQQIdAUnJ2pqG6ETw7pnX58jjMKQYmgSIIVIqtJIzvBJr1WQgghQpLpgMn5Gj0SFJ8e0v1KsKSCFb8kKwiSAlKPH9CsyVwrIYQQQTIpOCdbcW62alkFyU3o/PMUkC0VrLglCVYQHGYTOS4bdrOMXhBCCBG8XgmKS/IUWW0qUWk2sJugtsn356bbvb2/Ij5JghUAk1JkOK2k2C0yjV0IIUTQkq3wrZ4ag3xMXVdK0TNB8UWl72VCWR6Mb50OGl2yZAlKKV544YVoxBOXFJBst9AnJYFUOepGCCFEkKwanNdDcdNpvpOrZp0tE0qDe3zrMMHasWMHTzzxBKNGjfLrYrddN4PXXnoxLIHFC4fZREGSkx4uO2YpxQohhAiCpmBEuuInp2mcm6P59XrSs5MKVY8EeU2KZz4TLI/Hw8yZM3nkkUew2Wx+XWzX9i+4dcbULpFkWTRFD5edgiQHDjnmRohupaysjDGnD+ehb2XwxFUj2bRqZaxDEgalKRiaqrhpoMal+RqJVv+TolyX779TQLYj9PhE5PhMsBYuXMiYMWMYOXJkQBdUSvHoA/dTV1dHU2NjyAFGm6Ygw2mlT0oCydJrJUS3U1ZWRmlpKe4je5hzZj596vawfPZlkmSJgChgcKrixwM1pvTWSLUH/lriNCtSfNQ30u1gNcnrUzxrt8n9s88+Y8WKFaxbty7gC+q6zrbNnzN//nzGjBnDmDFjQg4ykmpragHvfwaXWZFs0TDVNnCktiq2gXWivLw81iEEReKOLqPGnZqaGrPH/u1v7uHc/AyWTRqOUorrhudz1Qvrefupexlw3pSYxSWMwaxgaJri7CxFehBJVVu5CYryulMb3XOkwT3utSRYS5cuZeHChQDccMMN7Nixg379+gGwb98+rr/+evbu3cuNN97Y4QWVUvQt7s/cuXMxm0yYzPG/UTHD7SIzwXhjF2L5IhQKiTu6jBp3rHy+eTO3nZ7fUr1WSjEhP40PPtgc48hEPEu0wBkZitMzFE5z+JKfnk7YcPjUP+8hDe5xryX7mTZtGtOmTWv5i9aJ1Pjx47n11luZPHlypxfUdZ0fz/6F331bsWQ3ayTZTfRMkmeqEMKrf3Exa3d9zXXDvUmWruus3nmIjILiWIcmosiiQYOn848zKfhuvsbgFCJyVFpuggKkgmVEnY5pCER+YR8efvpZvnXxd8N52bCzmjR6uuz0TnLikDVsIUQrv/zVr1m38wBX/309i9fv4KoX1vPWroOMmfmrWIcmouicbIU/rw7D0xTD0lTEzqHNdnqTuNYUMqLBCPxKsNasWeNX9WrB4iVxnVyZlSI7wUafZCdJ0sAuhGhHSUkJK1as4FhKLr/7YCfbHT25/IG/MWBC5/dA0XUUuFSnc6pMCsbmRPZ1xKwpsttOeZcGd0OI/wapMNAUpNqtpDmsEfspQwjRdZSUlHDBpVO4+70qHA7ZC98dJVpgXI7is8O+J6kPS1MkBTB2IVi5CbCn1b4rWR40hi6dYDVPYM9wWjFrYV0NFUII0YW5LN4q0aBUxQd7Tv17TcHY7OgkOt6BoycSPWlwN4Yum3UkWs0UJieQ47JLciWEEMJvNtOJJbhzffRiDU1VJNuik2C1PTJHKljG0OUqWE6zicwEG06Zvi6EECIIiZYTv89wKPq7PexoOPFnmoJzI9x71VqqXeE0Q3WjNLgbSZcp7VhNGrmJdgqSnZJcCSGECFqi5eTkaXRG00lVrCGpipQoVa+a9Tx+7mCqDWzS4G4Ihk+wTK12Brptls4/QQghhOiAq81LSbpNZ2CKN6mJZu9Vaz2PV63kgGfjMGyCpYBUu4W+KQmkOqwyckEIIURYJLbzs/q4HG8v1qAURVoYjsAJVO7xxEqWB43DkD1YLouJrAQ7NrNh80MhhBBxqr0EK9OhOC1FRbX3qrWeCd7CQg9pcDcMQyVYVpNGVoKNRKuhwhZCCGEgLkv7Scyl+SpmAz4dZkWaXSpYRmKITEVTkO6wkeqwoMlSoBBCiAhqr4IFsZ+ePjhVSYO7gcR9guW2mslKsGExyXKgEEKIyPOVYMXayHRJrowkbhMsq0kjO8GGS5YDhRABKCgowGaztRxxc8cdd3DFFVfEOCphJPGaYPlauhTxKe6yFwWkOaykO62yHCiECMrzzz/PsGHDYh2GMCCrFvulQNE1xFWC5TCbyHHZsJtlUKgQQojoi9fqlTCeuEiwNKXIdFpJsVtknpUQImTTpk1D13XOPPNM5s+fT0ZGhs+PPXbsGJWVlS3v22w2bDZbNMIUcSjRKq9BIjxinmC5LCayXXas0sQuhAiDdevWkZeXR0NDA3feeSfXXHMNr7zyis+PHzdu3Envz549mzlz5lDdCDU1nkiHGzY1NbWxDiEgcRuv1cPhw40n/VF5eXmMggmckWIFY8XbWaypqaknvR+zBMukFJkJNlLsUo8VQgRv6dKlLFy4EIBbbrmFGTNmAGCxWLj11lspKirq8PPXrl17Ur9WcwXL3qjjcFS1NMsbgZFihfiMNydZkZp66g/8bV8845mRYgVjxRtIrDFJsFwWEzkuu4xeEEKEbNq0aUybNg2Aqqoqjhw5QnJyMgDPPfccw4cP7/DzXS4Xbrc74nGKwJgVNOrRf9xEa/QfU3RNUU2wNAWZThupDnkGCyHC75tvvqG0tJSmpiZ0XaewsJClS5fGOiwRhAt7aWyp0NlcEd0sK1FGIYgwiVqC5TSb6JEovVZCiMgpLCzkk08+iXUYIgwK3TAsTfH8l7AlikmW7CIU4RLxbEcBmU4r+UkOSa6EEEJ0ym2BFJvCpCkuL1T0dUevquSSBEuESUQrWDaTRg+XHYdF5loJIYTwT37iiYTKrCm+3wee+wK+qPRdydIUuMzeMQsuM9hNcKgOvqnRaQhgM6hUsES4RCzBSrZbyE6wyTR2IYQQAcl3nfx+S5K1DXYc00m3Q45DkeOEbKciww5OM+3OUdR1ncN1sK/a+7kfHvCdpFk15DBlETZhT7BMSpHjsuG2yY8BQgghApfvOjXJsWiKq/qCjsKi+Z8EKaVIs0OaHYqS4eODOh4fOZYsD4pwCmtTlN2sUZjslORKCCFEUBLMkOFoP4Eya4ElV21ZNEVqB0P6ZQehCKewJlgy20oIIUQo8tqpXoVTto/kDaT/SoRXWLMhOUdQCCFEKNr2X4VbttP330mCJcJJyk1CCCHiRusdhJHQUQVLerBEOEmCJYQQIi7YTZAd4eMJO65gySqMCB9JsIQQQsSFXi4V8VYTl0X5rFTJEqEIJ0mwhBBCxIVI918187VMKAc9i3CSBEsIIURcKIjwDsJmvpYJXVE7nVd0B5JgCSGEiDmrBjkd9EeFU3sVLIsGdrP0YInwkQRLCCFEzOUmeA93job2KljSfyXCTRIsIYQQMZefGL3HSrN5K1atyQ5CEW6SYAkhhIi59s4fjBSlFJltlgllBpYIN0mwhBCim5tSoHFutsIdoyTDpCA3IbqP2XbeliwRinCTPRNCCNGNWTUYlAImTWN8D52tFbD+oM7WSh2PHp0YeiYozFHqv2rWtg9LEiwRbpJgCSFEN1aQeKK5XFOK4mQoTlYcrtV5ZquHivrIxxCt+VeteXcSnsggXdKDJcJMlgiFEKIb6+Nu/89T7Ypr+mkRr+zYTVCUFP3kJssBrR9VKlgi3CTBEkKIbqyP23dyk2pXTOunkRCBtY5Mu8538zR+NlijVxQb3JtZTYpU24n3JcES4SZLhEII0U25LZBu7zi5yXAoftBP45ktHmqaQns8TcGgFMUZGYqE+gZSU2O7LJftVByq8y4TSoIlwq3DClZdXR033XQT/fr1Y/DgwUydOjVacQkhhPDBboK+bsWEHt5lvIk9g0tUOqpetZbtVFzdV8NmCuphWmTaoaR3bCpW7WneSShT3EUkdFjBmjt3LkoptmzZglKKffv2RSsuIYQQbeS5FBf3UmQ6vLOcmmU44KUvAr+evwkWQK5LcVUfjae3eAh2c2GWj0OWYyXreKO7zMASkeAzwaqqquKpp57iq6++avmPnJ2dHbXAhBBCnDAwWVHSu/1xBi6LopfTw36P/9dTQKGPBndf8hMVSVY4EuTOwkxH5x8TTc2jGmSKu4gEn0uEX3zxBampqdx3332cfvrpjB07ljfeeKPDix07dozKysqWt7q6urAHLIQQ3c3ZmYrLCjueFdXfHUB2hTe5cAaxLJYRQhUq3ipYbqsiwSz9VyIyfFawGhsb2blzJwMHDmT+/Pl88sknTJw4kY0bN5KVldXu54wbN+6k92fPns2cOXPCG3GYlZeXxzqEoEjc0SVxR1dqamqsQ4gLCrigl+KszM43fBe7PbxdAU1+rt8FsjzYWqYdtlYE9alkxVkFC7xJnywRikg4KcFaunQpCxcuBODqq69G0zSuvvpqAIYPH07v3r3ZsGGDzwRr7dq1DBs2rOV9m82GzWZr92PjiVFv5hJ3dEncIprMCkp7awxI8S8Rcpig0K3YWuFfhlWYGGSC1WZAp7+cZki0xlcFC5orebGOQnRFJz2tpk2bxrRp01ref+211/jnP//JRRddxPbt29m+fTsDBgzweTGXy4XbHeCivhBCiFMMTVN+J1fNBqX4l2BZNMgLcnp6hj24z2t7uHK8yHaooJv2hehIh3n7448/zrXXXsucOXPQNI1FixbRs2fPaMUmhBDd1sAAkyuA/sneyldjJxlDviv4s/8yjk9ADzQpicflQfBWsI41xDoK0RV1mGAVFhayevXqaMUihBAC75JV78TAP89mUvRLUmw60nH64+t4HH9YNEWKDQ4HuIcp3hrcm6XbvRU9IcJNnlZCCBFn+icrNBVcQjKok8rXplUrueHbI0hwOhg5fBhlZWUBP0ZGJ9Pf25MZ5NJipGlKkWKLz+RPGJskWEIIEWcGJgf/gl+UDFYfd/ZNq1ayfPZlpFR8zW2n5+Mq/5rS0tKAk6xA51mpID5HCKOTBEsIIeKIwxTc8mAzi6YoSmo/QXvryXs5Nz+dZZOGc/2IApZNGs7Y/HTuu/c3AT1GoA3rKTbv4cpCdCeSYAkhRIiKfSQ0weifrDAF2YDebJCPQ5QP7tzMuLy0ltM5lFKMz0tj0+efB3T9QHcSxmv/lRCRJAmWEEKEQAGX5Cl6JYQniQhm92Bb/dzeA6Hbyu5dzNpdh9B1bxO8ruus2XWIAf37B3T9dDsEkgPK8qDojmS8mhBChCA/UZFoVUzoAUu3hjZRyW6CwhCWB5uZNEX/ZMWnh7wHGQ9KUQxOVQz57a8pLS3l6r+vZ3xeGmt2HeLNnQcp+5/FAV3frClSbXCw1r+PlwqW6I6kgiWEECEYlOL9tdCtyHd1nEjYTJDXwceEY3mw2dmZih/00/jZYI0Lemn0TFCUlJSwYsUKqlNzefCjnVSn5lJWVsaUKVMCvn4gOwnjdQaWEJEkFSwhhAiSpmBAqx1/43MUz3RQxbq4l0a/JFj8uU55O3OkQtk92FaWU9HeoWYlJSWUlJSEfP1MB2w60vnHWTRIjf8T04QIO6lgCSFEkHonKhIsJ5Ki3m5FgY8z/oamKoakKRxmxfcLtVNGKdhNoQ0AjbZMPytYGXbV0lQvRHciCZYQQgTptHYa0ifknPpnaTa4OO/En2c5FZMLTr79FieFb3kwGjL8XPaT5UHRXUmCJYQQQTApGJB86p/nJyp6t6pimRR8r1A7ZQ7UwBTFudnqpPeNJM3m/bd1RhIs0V1JgiWEEEEodHuX+9ozoceJP5/YU5Hj9P1xRUkKm8GWB8G7UzHNj3lYgQ4lFaKrkCZ3IYQIQnvLg83yXIo+boVZwags3z/HKqUo7Q3v7/eOPjCaDLtif03HoymkgiW6K0mwhBAiQGYF/ZM6/piJPRVua+fXspkU57bTt2UEmQ7YWO77710WTtoEIER3IgmWEEIEqI9bYfexPNgs28eyYFfi3Unou4IlA0ZFdyY9WEIIESBfZ/11N53tJMwM8MxCIboSSbCEECIAFg2KO1ke7C5Sbd7lUl+yukEVTwhfJMESQogA9HOrU0YudFeaUqR3UKWSBnfRnUmCJYQQfjIpGJomyVVrGT76rDQFGbJEKLoxaXIXQogOmBX0TVIMTFYUJdFpc3t346vPKtVmzNETQoSLJFhCCNEOiwaTchs5vdepU9jFCd5BoqfuJOwp/Veim5MlQj+VlZUxdOhQHA4HQ4cOpaysLNYhCSEiyKIp+rs9klx1ou0yoEWD83soLs2Xr5vo3qSC5YeysjJKS0tRSqHrOhs2bKC0tJQVK1ZQUlIS6/CEECJmUmzepKrBA/2SFBf1UqTYJLkSQipYfrj77rtbkisAXddRSnHPPffEODIhhIgtpRSFiYrLCzWu7qtJciXEcVLB8sOWLVtakqtmuq6zefPmGEUkhBDx4/uN/THIAAAFQ0lEQVR9FEpJYiVEa1LB8kNRUdEpNw+lFMXFxTGKSAgh4ockV0KcShIsP8ybN69lWRBoWS6cN29ejCMTQgghRDySBMsPJSUlrFixgiFDhmC32xkyZAhlZWVMmTIl1qEJIYQQIg5JD5afSkpKZMegEEIIIfwSlgpWfX39Sb8aRV1dHb/73e+oq6uLdSgBkbijS+KOrrq6Ou66666Yx22kr5+RYgVjxSuxRo6R4g0mVqW33R4XhHXr1jFu3DjWrl3LueeeG+rloqayspKkpCQqKipwu92xDsdvEnd0SdzRFa24169fz8iRI/n4448ZMWJEzOIIByPFCsaKV2KNHCPFG0ys0oMlhBBCCBFmkmAJIYQQQoRZWJrca2trAe9ATpfLFY5LRsWxY8cA+PTTTyXuKJC4o8vocVdXV0d02aCmpgaATZs2dRiHEb5+RooVjBWvxBo5RorX31j79++P0+n0vqOHwdNPP63jPU5d3uRN3uQtLG/vvvtuOG5PPi1btizm/0Z5kzd561pvH3/8ccs9JixN7gcPHuSf//wnBQUFOByOUC8nhBAn/yQYAXLfEkKEW+v7VlgSLCGEEEIIcYI0uQshhBBChJkkWEIIIYQQYRaWBOvb3/42Q4YMYdiwYYwdO5ZPPvkkHJeNqNraWiZPnkxRURFDhw5l4sSJbNu2LdZh+eXmm2+moKAApRSffvpprMPxy9atWxk9ejRFRUWcccYZbNy4MdYh+cWIX2sjP7fj6V4ST7F0xIjfb6P8vzLSfcsoX1Mw5nM2qPtBOHbjlJeXt/y+rKxMHzJkSDguG1E1NTX6yy+/rHs8Hl3Xdf2RRx7Rx40bF9ug/LR27Vp99+7den5+vv7JJ5/EOhy/TJgwQV+yZImu67r+17/+VT/99NNjG5CfjPi1NvJzO57uJfEUS0eM+P02yv8rI923jPI11XVjPmeDuR+EpYKVnJzc8vuKigqUUuG4bETZ7XYuuuiillhHjRrFjh07YhuUn84991xyc3NjHYbf9u/fz0cffcTUqVMBKC0tZffu3XH/EwsY72sNxn5ux9O9JJ5i6YgRv99G+H9ltPuWEb6mzYz4nA3mfhCWQaMA06ZNY/Xq1QC88sor4bps1Dz88MNMmjQp1mF0Sbt37yYnJwez2ft0U0qRl5fHrl276Nu3b4yj6/qM9tyOp3tJPMXiL6N9v+OV3LeixyjP2UDvB2FLsJYuXQrAM888w5w5cwxzMwK477772LZtG2+88UasQxEirIz43I6ne0k8xeIPI36/RfdmpOdsoPeDoJYIly5dyrBhwxg2bBhLliw56e+uueYaVq9ezaFDh4K5dES1F/eCBQsoKyvj1VdfjehQw1B09PU2gl69erF3714aGxsB0HWdXbt2kZeXF+PIujYjPLc7Eu17iZHua0a7lxnxHib3rciL5+dsR/y+H4Sj8WvPnj0t769cuVLv2bNnS/NaPHvwwQf1ESNG6IcPH451KEExQjNjs3Hjxp3ULDpy5MjYBhQgI32tdd2Yz+14upfEUyz+MOL3W9fj//+VEe9b8f41bWak52yw94OQJ7nv3LmTyy67jJqaGjRNIyMjgwULFjBs2LBQLhtxX331Fb169aKwsJDExEQAbDYb77//fowj69wNN9zAyy+/zL59+0hLSyMxMTFuGy+bbd68menTp3Po0CHcbjdLlixh8ODBsQ6rU0b8Whv1uR1P95J4iqUzRvx+G+X/lZHuW0b5moLxnrPB3g/kqBwhhBBCiDCTSe5CCCGEEGEmCZYQQgghRJhJgiWEEEIIEWb/HxrsbUSjTDYOAAAAAElFTkSuQmCC" }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Plot posterior samples\n", "xtest = range(minimum(gpa.x),stop=maximum(gpa.x),length=50);\n", "\n", "#Set the parameters to the posterior values the sample random function\n", "fsamples = [];\n", "for i in 1:size(samples,2)\n", " set_params!(gpa,samples[:,i])\n", " update_target!(gpa)\n", " push!(fsamples, rand(gpa,xtest))\n", "end\n", "\n", "#Predict\n", "p1=plot(gpe,leg=false, title=\"Exact GP\") #Exact GP (assuming Gaussian noise)\n", "\n", "\n", "sd = [std(fsamples[i]) for i in 1:50]\n", "p2=plot(xtest,mean(fsamples),ribbon=2*sd,leg=false, title=\"Monte Carlo GP\") #GP Monte Carlo with student-t noise\n", "scatter!(X,Y)\n", "\n", "plot(p1,p2;fmt=:png)" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.0.0", "language": "julia", "name": "julia-1.0" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.0.0" } }, "nbformat": 4, "nbformat_minor": 2 }