{ "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": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHXRJREFUeJzt3W9sXfV9P/DPSSPdEFwXxaVSGaQhUuPErhK7qSY0iQW0OEAIEr73wR4QGdZJ2QOqFK1azMQgfwgRk1A02FRlVSVG1IkHk09VyftDEtTmEUWlEVWIsZMS0sZap6rFwjIlVrqd3wMX/+qSOA7+3nN97ddLspRzfXL84YN1efM53/s9WVEURQAAkMyyRhcAALDYXDNgTU5Oxte+9rX44he/GJ2dnbFz584y6gIAaFrLr3XC448/HsuWLYuzZ89GlmXxi1/84mPn/OY3v4nh4eFYv359rFy5si6FAgA0i2y2NVgffPBB/NEf/VGMjo5GS0vLVS9y6tSp2Lx5c/z4xz+OL3/5y5+okPfffz8+85nPfKK/y/XR6/LodXn0ujx6XR69Lk/qXs96i/Cdd96Jtra2OHjwYHzlK1+JO++8M1599dWrnj8xMRHj4+PTX5OTk3Mu5H//93/nXjXzotfl0evy6HV59Lo8el2e1L2e9Rbh5cuX4/z589HR0RHPPvts/OQnP4mtW7fG0NBQ3HzzzR87f8uWLTOO9+zZE/39/XMqZGxs7DrKZj70ujx6XR69Lo9el0evyzPfXq9atWrG8awB6wtf+EIsW7YsHnrooYiI2LRpU9x+++1x5syZuOuuuz52/smTJ6Orq2v6uFKpRKVS+cTFUT96XR69Lo9el0evy6PX5UnZ61lvEX72s5+NP/uzP4tXXnklIiJ+9rOfxbvvvhvt7e1XPL+lpSVaW1unv64nXAEALBbX/BThkSNH4qtf/Wr09/fHpz71qfjWt74Vn//858uoDQCgKV0zYK1duzZ+8IMflFAKAMDiYCd3AIDEBCwAmlqe57Gpe3PcsPLG2NS9OfI8b3RJIGAB0LzyPI9arRanL7fFpR374vTltqjVakIWDSdgAdC09j/9TGSdPVHsHozoeSyK3YORdWyNAwcPNbo0ljgBC4CmdXZkOIqOnogsm3ohy6Lo3BYjw283tjCWPAELgKa1rn19ZEPHIz56rG5RRHbmWLRv2NDYwljyrrlNAwAsVHuffCJqtVpkz98fRee2yM4ci2LoROy1BosGM8ECoGlVq9UYGBiIjZWxWDG4LzZWxiLP8+jt7W10aU3PpzPnR8ACoKlVq9V489Qb8eFvPog3T70hXCXg05nzJ2ABADP4dOb8CVgAwAw+nTl/AhYAMINPZ86fTxECADP4dOb8mWABADP4dOb8mWABAB9TrVajWq02uoymZYIFAJCYgAUAkJiABQCQmIAFAJCYgAUAkJiABQCQmIAFAJCYgAUAkJiABQCQmIAFAJCYgAUAkJiABQCQmIAFAJCYgAUAkJiABQCQmIAFAJCYgAUALFl5nsem7s1xy623xabuzZHneZLrClgAwJKU53nUarU4fbktJnfsi9OX26JWqyUJWQIWALAk7X/6mcg6e6LYPRjR81gUuwcj69gaBw4emve1BSwAYEk6OzIcRUdPRJZNvZBlUXRui5Hht+d9bQELAFgQPloPdcPKG5Ouh7qade3rIxs6HlEUUy8URWRnjkX7hg3zvraABQA03O+vh7qUeD3U1ex98okozhyP7Pn7I47/Q2TP3x/F0InY+3dPzPvaAhYA0HD1XA91NdVqNQYGBmJjZSwqg/tiY2Us8jyP3t7eeV97eYL6AADm5ezIcBQ79n18PdTgvrr+3Gq1GtVqNd57771YtWpVsuuaYAEADVfP9VCNYIIFADTc3iefiFqtNrUOqnNbZGeOTa2HqvNC93oxwQIAGu7310OtSLweqhFMsACABeGj9VCLgQkWAEBiAhYAQGICFgBAYgIWAEBicw5Y+/fvjyzL4q233qpnPQAATW9OAevUqVPxwx/+MFavXl3vegAAmt41A9bk5GQ8+uij8c1vfjOyj7avBwDgqq65D9ZTTz0VO3fujNtvv/2aF5uYmIjx8fHp40qlEpVKZX4VAgA0mVkD1muvvRY/+tGP4tlnn53TxbZs2TLjeM+ePdHf3z+nvzs2Njan85g/vS6PXpdHr8uj1+XR6/LMt9d/+KDoWQPWyZMnY3h4eHp6NTo6Gvfcc098+9vfjvvuu++K53d1dU0fX+8EK+VTrJmdXpdHr8uj1+XR6/LodXlS9nrWgPX444/H448/Pn28Zs2aGBwcjC996UtXPL+lpSVaW1uTFQcA0IzsgwUAkNh1Pez5woULdSoDAGDxMMECAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACgCaQ53ls6t4cN6y8MTZ1b448zxtdErMQsABggcvzPGq1Wpy+3BaXduyL05fbolarCVkLmIAFAAvc/qefiayzJ4rdgxE9j0WxezCyjq1x4OChRpfGVQhYALDAnR0ZjqKjJyLLpl7Isig6t8XI8NuNLYyrErAAYIFb174+sqHjEUUx9UJRRHbmWLRv2NDYwriq5Y0uAACY3d4nn4harRbZ8/dH0bktsjPHohg6EXutwVqwTLAAYIGrVqsxMDAQGytjsWJwX2ysjEWe59Hb29vo0rgKEywAaALVajWq1Wqjy2COTLAAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAEgmz/PY1L05blh5Y2zq3hx5nje6JGiIWQPWpUuX4sEHH4x169ZFV1dX3HvvvXHhwoWSSgOgmeR5HrVaLU5fbotLO/bF6cttUavVhCyWpGtOsHbt2hUjIyPx5ptvxo4dO2LXrl1l1AVAk9n/9DORdfZEsXswouexKHYPRtaxNQ4cPNTo0qB0swasFStWxPbt2yPLsoiIuOOOO+L8+fOlFAZAczk7MhxFR0/E7/6bEVkWRee2GBl+u7GFQQNc1xqsF154IR544IGrfn9iYiLGx8envyYnJ+ddIADNYV37+siGjkcUxdQLRRHZmWPRvmFDYwuDBlg+1xMPHToU586diyNHjlz1nC1btsw43rNnT/T398/p+mNjY3MthXnS6/LodXn0ujxX6/VfP7Y7Hnnkkcie3x5F5z2RnXkliqFX469feinee++9kqtcHPxel2e+vV61atWM4zkFrOeeey7yPI8TJ07EypUrr3reyZMno6ura/q4UqlEpVL5xMVRP3pdHr0uj16X50q9fvjhh+PTn/50HDh4KEYG90X7+g2xN8+jt7e3ARUuHn6vy5Oy19cMWIcPH46XX345Tpw4ETfddNOs57a0tERra2uy4gBoLtVqNarVaqPLgIabNWCNjo7GN77xjVi7dm3cfffdETE1lXr99ddLKQ4AoBnNGrBuvfXWKD5arAgAwJzYyR0AIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsAIDEBCwAgMQELACAxAQsWuDzPY1P35rhh5Y2xqXtz5Hne6JIAuAYBCxawPM+jVqvF6cttcWnHvjh9uS1qtZqQBbDACViwgO1/+pnIOnui2D0Y0fNYFLsHI+vYGgcOHmp0aQDMQsCCBezsyHAUHT0RWTb1QpZF0bktRobfbmxhAMxKwIIFbF37+siGjkcUxdQLRRHZmWPRvmFDYwsDYFbLG10AcHV7n3wiarVaZM/fH0XntsjOHIti6ETstQYLYEEzwYIFrFqtxsDAQGysjMWKwX2xsTIWeZ5Hb29vo0sDYBbXnGCdO3cuHn744fjVr34VN910U/zLv/xLdHR0lFEbEFMhq1qtNroMAK7DNSdYf/VXfxW7du2Ks2fPxp49e+Iv//Ivy6gLAKBpzRqwfvnLX8apU6di586dERFRq9Xi3XffjQsXLpRRG0Dd2MAVqKdZA9bFixfjlltuieXLp+4kZlkWq1evjp///OdXPH9iYiLGx8envyYnJ9NXDDBPNnAF6u2aa7Cyj/bf+Z3io4+LX8GWLVtmHO/Zsyf6+/vnVMjY2NiczmP+9Lo8el2e6+n1U/sORNa5dWoD1yyLYuvXI3t+ezy1/+m466676lfkIuH3ujx6XZ759nrVqlUzjmcNWLfddluMjo7Gb3/721i+fHkURREXL16M1atXX/H8kydPRldX1/RxpVKJSqXyiYujfvS6PHpdnrn2+p2fnotix0N/sIHrPfHO4D7/vuZIn8qj1+VJ2etZbxF+7nOfi+7u7vjOd74TEREDAwOxZs2aWLNmzRXPb2lpidbW1umv6wlXAGWxgStQb9e8RfjP//zP8cgjj8ShQ4eitbU1XnrppTLqAqgbG7gC9XbNbRra29vjtddei7Nnz8Ybb7wRnZ2dZdQFUDc2cAXqzaNygCXJBq5APXlUDgBAYgIWAEBiAhZASeweD0uHgAVQArvHw9IiYAEfY9KS3v6nn4mss2dq9/iex6LYPRhZx9Y4cPBQo0sD6kDAAmYwaamPsyPDUXT0/MHu8dtiZPjtxhYG1IWABcxg0lIfdo+HpUXAAmYwaamPvU8+EcWZ45E9f3/E8X+Y2kV+6ETs/bsnGl0aUAcCFjCDSUt92D0elhY7uQMzeE5f/dg9HpYOEyxgBpMWgPkTsICPqVar8eapN+LD33wQb556Q7hqYrbcgMYQsAAWKVtuQOMIWACLlC03oHEELIBFypYb0DgCFlwna1poFrbcgMYRsOA6WNNCM7G5KTSOgAXXwZqW+jEZTM+WG9A4AhZcB2ta6sNksH5suQGNIWDBdbCmpT5MBoHFRsCC62BNS32YDNJs3NLmWgQsuA7WtNSHySDNxC1t5sLDnuE6eWBveh4wTTOZcUs7y6LY+vXInr8/Dhw85L2BaSZYQMOZDNJM3NJmLkywgAXBZJBmsa59fZweOh7F1q9PhSy3tLkCAQsAroNb2syFW4QAcB3c0mYuTLAA4Dq5pc21mGABACQmYAEAJCZgAQAkJmDR1DyuAoCFSMCiaXlcBQALlYBF05rxuIqex6LYPRhZx9Y4cPBQo0sDYIkTsGhaHlcBwEIlYNG01rWvj2zoeERRTL3gcRUALBA2GqVpeVwFAAuVCRZNy+MqAFioTLBoah5XAcBCZIIFAJCYgFUSG2ICwNIhYJXAhpgAsLQIWCWwISYALC0CVglsiAkAS4uAVQIbYgLA0mKbhhLYEBMAlpYrTrAuXboUDz74YKxbty66urri3nvvjQsXLpRc2uJhQ0wAWFquOsHatWtX3HfffZFlWfzTP/1T7Nq1K44dO1ZmbYuKDTEBYOm44gRrxYoVsX379sh+tyj7jjvuiPPnz5daGABAs5rTGqwXXnghHnjggWueNzExEePj49PHlUolKpXKJ68OAKAJXTNgHTp0KM6dOxdHjhy55sW2bNky43jPnj3R398/p0LGxsbmdB7zp9fl0evy6HV59Lo8el2e+fZ61apVM46nA9bRo0fj8OHDERHx9a9/Pf7iL/4innvuucjzPE6cOBErV6685sVPnjwZXV1d08fXO8H6w+KoH70uj16XR6/Lo9fl0evypOz1dMDq6+uLvr6+6W8cPnw4Xn755Thx4kTcdNNNc7pYS0tLtLa2JisOAKAZXfEW4ejoaHzjG9+ItWvXxt133x0RU9Oo119/vdTiAACa0RUD1q233hrFR7uOAwBwXTwqZ5HK8zw2dW+OG1beGJu6N0du13gAKI2AtQjleR61Wi1OX26LSzv2xenLbVGr1YQsACiJgLUI7X/6mcg6e6LYPRjR81gUuwcj69gaBw4eanRpALAkCFiL0NmR4Sg6eiJ+txN/ZFkUndtiZPjtxhYGAEuEgLUIrWtfH9nQ8YiPPqhQFJGdORbtGzY0tjAAWCLm9KgcmsveJ5+IWq0W2fP3R9G5LbIzx6IYOhF7rcECgFKYYC1C1Wo1BgYGYmNlLFYM7ouNlbHI8zx6e3sbXRoALAkmWItUtVqNarXa6DIAYEkywQIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhsSQasPM9jU/fmuGHljbGpe7OHIAMASS25gJXnedRqtTh9uS0u7dgXpy+3Ra1WE7IAgGSWXMDa//QzkXX2RLF7MKLnsSh2D0bWsTUOHDzU6NIAgEViyQWssyPDUXT0RGTZ1AtZFkXnthgZfruxhQEAi8aSC1jr2tdHNnQ8oiimXiiKyM4ci/YNGxpbGACwaCy5ZxHuffKJqNVqkT1/fxSd2yI7cyyKoROx1xosACCRJTfBqlarMTAwEBsrY7FicF9srIxFnufR29vb6NIAgEViyU2wIqZCVrVabXQZAMAiteQmWAAA9SZgAQAkJmABACQmYAEAJCZgAQAkJmABACQmYAEAJCZgAQAkJmABACQmYAEAJCZgAQAkJmABACQmYAEAJCZgAQAkJmABACQmYAEAJCZgAQAkJmABACQmYAEAJCZgAQAkJmABACQmYAEAJCZgAQAkJmABACQmYAEAJDangLV///7IsizeeuutetcDAND0rhmwTp06FT/84Q9j9erVdSkgz/PY1L05brn1ttjUvTnyPK/LzwEAKMusAWtycjIeffTR+OY3vxlZliX/4XmeR61Wi9OX22Jyx744fbktarWakAUANLVZA9ZTTz0VO3fujNtvv31OF5uYmIjx8fHpr8nJyVnP3//0M5F19kSxezCi57Eodg9G1rE1Dhw8NPd/AgCABWb51b7x2muvxY9+9KN49tln53yxLVu2zDjes2dP9Pf3X/X8kZHhKHbsi/hoOpZlUXRui+HBffHee+/N+edyfcbGxhpdwpKh1+XR6/LodXn0ujzz7fWqVatmHM8IWEePHo3Dhw9HRMSf//mfx/Dw8PT0anR0NO6555749re/Hffdd98VL37y5Mno6uqaPq5UKlGpVK5aTHv7+jg9dDyKrV+fCllFEdmZY7F+w4aPFUpa+lsevS6PXpdHr8uj1+VJ2esZAauvry/6+vqmj//2b/92+s9r1qyJwcHB+NKXvnTVi7W0tERra+ucf/jeJ5+IWq0W2fP3R9G5LbIzx6IYOhF7rcECAJpYQ/fBqlarMTAwEBsrY1EZ3BcbK2OR53n09vY2siwAgHm56hqsP3ThwoW6FFCtVqNarcZ7771nDAoALAp2cgcASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBITMACAEhMwAIASEzAAgBIbNaANTk5GV/72tfii1/8YnR2dsbOnTvLqgsAoGktn+2bjz/+eCxbtizOnj0bWZbFL37xi7LqAgBoWlcNWB988EG8+OKLMTo6GlmWRUTE5z//+dIKAwBoVle9RfjOO+9EW1tbHDx4ML7yla/EnXfeGa+++uqsF5uYmIjx8fHpr8nJyeQFAwAsdFedYF2+fDnOnz8fHR0d8eyzz8ZPfvKT2Lp1awwNDcXNN998xb+zZcuWGcd79uyJ/v7+ORUyNjZ2HWUzH3pdHr0uj16XR6/Lo9flmW+vV61aNeN4RsA6evRoHD58OCIiHnrooVi2bFk89NBDERGxadOmuP322+PMmTNx1113XfHiJ0+ejK6urunjSqUSlUrlExdH/eh1efS6PHpdHr0uj16XJ2WvZwSsvr6+6Ovrmz4+fvx4vPLKK7F9+/b42c9+Fu+++260t7df9WItLS3R2tqarDgAgGY066cIjxw5El/96lejv78/PvWpT8W3vvUtC90BAK5h1oC1du3a+MEPflBSKQAAi4Od3AEAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7AAABITsAAAEhOwAAASE7BIJs/z2NS9OW5YeWNs6t4ceZ43uiQAaAgBiyTyPI9arRanL7fFpR374vTltqjVakIWAEuSgEUS+59+JrLOnih2D0b0PBbF7sHIOrbGgYOHGl0aAJROwCKJsyPDUXT0RGTZ1AtZFkXnthgZfruxhQFAAwhYJLGufX1kQ8cjimLqhaKI7MyxaN+wobGFAUADLG90ASwOe598Imq1WmTP3x9F57bIzhyLYuhE7LUGC4AlyASLJKrVagwMDMTGylisGNwXGytjked59Pb2Nro0ACidCRbJVKvVqFarjS4DABrOBAsAILEFEbAmJyfj7//+72NycrLRpSx6el0evS6PXpdHr8uj1+WpR6+zovjoY1+f3KlTp2Lz5s3x4x//OL785S9f998fHx+Pz3zmM/H+++9Ha2vrfMthFnpdHr0uj16XR6/Lo9flqUevF8QECwBgMRGwAAASS/Ipwg8//DAiIt5++5Pt2j0xMREREW+++Wa0tLSkKImr0Ovy6HV59Lo8el0evS5Pql6vX78+Vq5cGRGJ1mD967/+a+zcuXO+lwEAaFq/vxY9ScD61a9+Fa+88kqsWbMmbrjhhnkXCADQbJJPsAAA+P8scgcASEzAAgBIbEEErG3btsXGjRujq6sr7rzzznjzzTcbXdKidOnSpXjwwQdj3bp10dXVFffee29cuHCh0WUtWrt37441a9ZElmXx1ltvNbqcRevcuXPxJ3/yJ7Fu3br44z/+4xgaGmp0SYuW3+lyeK8uV90ySLEAjI2NTf/5u9/9btHd3d3AahavDz/8sPj3f//34v/+7/+KoiiKf/zHfyx6enoaXNXidfLkyeLixYvFF77wheL06dONLmfRuvvuu4sXX3yxKIqi+Ld/+7fijjvuaGxBi5jf6XJ4ry5XvTLIgphg3XTTTdN/fv/992PZsgVR1qKzYsWK2L59e2RZFhERd9xxR5w/f77BVS1ef/qnfxq33npro8tY1H75y1/GqVOnpreJqdVq8e677/q//TrxO10O79XlqlcGSbLRaAp9fX3x/e9/PyIi/uu//qvB1SwNL7zwQjzwwAONLgM+sYsXL8Ytt9wSy5dPvZVlWRarV6+On//857FmzZrGFgeJeK+uv3pkkAUTsI4ePRoRES+99FL8zd/8TfzHf/xHgyta3A4dOhTnzp2LI0eONLoUmJeP/i//I4WdZ1hEvFeXox4ZpCH34o4ePRpdXV3R1dUVL7744ozvPfzww/H9738/fv3rXzeitEXnSr1+7rnnIs/z+M///M/pDdGYv9l+r6mP2267LUZHR+O3v/1tREyFq4sXL8bq1asbXBnMn/fq8qXMIA2ZYPX19UVfX19ERIyPj8d///d/xy233BIREd/97nejra0tVq1a1YjSFp3f73VExOHDh+Pll1+OEydOzLjvzPz9Ya+pv8997nPR3d0d3/nOd+KRRx6JgYGBWLNmjduDND3v1eUYHx+PiYmJumSQhu/kfvHixajVavHhhx/GsmXL4uabb47nnnsuurq6GlnWojQ6Ohq33XZbrF27Nj796U9HRESlUonXX3+9wZUtTo8++mh873vfi//5n/+Jz372s9HS0hI//elPG13WojMyMhKPPPJI/PrXv47W1tZ46aWXorOzs9FlLUp+p8vhvbo89cwgDQ9YAACLjf0QAAASE7AAABITsAAAEvt/FZETtbpr7jMAAAAASUVORK5CYII=" }, "execution_count": 1, "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": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "GP Monte Carlo object:\n", " Dim = 1\n", " Number of observations = 20\n", " Mean function:\n", " Type: MeanZero, Params: Float64[]\n", " Kernel:\n", " Type: Mat32Iso, 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": 2, "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", "gpmc = GPMC(X, vec(Y), MeanZero(), Matern(3/2,0.0,0.0), l) #Monte Carlo 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": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Results of Optimization Algorithm\n", " * Algorithm: L-BFGS\n", " * Starting Point: [0.5,0.0,0.0]\n", " * Minimizer: [0.6566151884322414,1.7118340256199516, ...]\n", " * Minimum: 4.578633e+01\n", " * Iterations: 13\n", " * Convergence: true\n", " * |x - x'| ≤ 0.0e+00: false \n", " |x - x'| = 2.09e-08 \n", " * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: true\n", " |f(x) - f(x')| = 0.00e+00 |f(x)|\n", " * |g(x)| ≤ 1.0e-08: true \n", " |g(x)| = 4.84e-12 \n", " * Stopped by an increasing objective: false\n", " * Reached Maximum Number of Iterations: false\n", " * Objective Calls: 38\n", " * Gradient Calls: 38" ] }, "execution_count": 3, "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": 4, "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.016300 \n", "Number of function calls: 100164\n", "Acceptance rate: 0.377800 \n" ] } ], "source": [ "set_priors!(gpmc.lik,[Normal(-2.0,4.0)])\n", "set_priors!(gpmc.kernel,[Normal(-2.0,4.0),Normal(-2.0,4.0)])\n", "\n", "samples = mcmc(gpmc;nIter=10000,burn=1000,thin=2);" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XlcVGX7P/DPgYFhEAZF3BdQU9RS0TZ/paGlmaamoGamGbmXmUtB9ZSWWVkuj/X0uOVuZS7gkqWVuaSlqVE+roiyo7LDMMw+5/r9MTFfkR2HOWfger9evGJmzpy5wObmM/d9n/sWiIjAGGOMMcYcxk3qAhhjjDHG6hoOWIwxxpjM9OvXD4IgSF0GuwscsFxcUlISBEGo8CskJETqMssUFBSEoKCgGj//+vXrmDt3LkJCQtCwYUN4eHigSZMmCA0NxYcffoi0tLQyX/P23427uzsCAgLw5JNPYu/evXfx0zDmum5vR1q1agWr1VrmcefPn7cf17lzZ6fWePToUQiCgPfee88pr0dEiImJQVhYGFq3bg2lUglfX1/06NEDc+bMwaVLl5xSR237+++/MX36dHTt2hVqtRqenp5o0aIFnnzySaxYsQI5OTmlnnPn3xiFQoEWLVpgxIgR+PXXXyX4KeRJIXUBzDE6dOiA8ePHl/lY8+bNnVxN7Vu+fDmioqJgtVrRu3dvTJgwAWq1Grm5uTh9+jTeffddLFy4EGfOnEH37t1LPNfd3R3vvPMOAMBkMuHKlSvYt28ffv75ZyxduhTz5s2T4kdiTHIKhQI3btzAjz/+iCFDhpR6fP369VAoFLBYLBJU5zy5ubkYPXo0Dh8+jIYNG2LgwIFo3749TCYTLl68iJUrV+Lzzz/HL7/8gn79+kldbo2IoojIyEgsW7YMCoUCjz32GJ588kl4e3sjMzMTv//+O+bMmYP58+cjISEBAQEBJZ7fuHFjzJw5EwCg1+tx7tw57N27F/v27cOOHTswatQoKX4seSHm0hITEwkADRo0SOpSqi0wMJACAwOr/bzVq1cTAOrQoQP9+eefZR5z6dIlCg8Pp+PHj5d6TaVSWer4H3/8kQRBIG9vbyoqKqp2TYy5suJ25LHHHiM/Pz8KDw8vdYzRaKSAgAAaPnw4AaDg4GCn1njkyBECQAsWLKjV1zGbzfTYY48RABo/fjwVFBSUOubGjRsUERFBe/bsqbU6QkNDqTb/RL/55psEgB544AG6du1amcecPn2a+vXrR6mpqSXuL+/f/8svvyQAFBQUVCs1uxoOWC6uugHrgw8+IAA0c+bMUo8tWLCAANDs2bPt96Wnp9P8+fPp4YcfpiZNmpCnpycFBgbSjBkzKCMjo8zXMBqNtGLFCnrwwQfJx8eHGjRoQF26dKE5c+ZQbm6uveayviprPHNzc0mtVpOXlxddvXq10p/XbDaXuF1ewCIi6ty5MwGgM2fOVHpexuqS29uRqVOnkqenJ2VlZZU4ZufOnQSAdu/eXe4f2KKiIlqwYAEFBweTUqmkRo0a0ZAhQ+i3334rdWxxe3PkyBHavn079ezZk7y8vKh58+b06quvkk6nK3VsWV+JiYn244xGIy1btox69uxJ3t7e5OPjQ3369KG9e/dW+XexYcMGe9i0Wq0VHmswGOzfHz58mCIiIqhTp07UoEEDatCgAd1///20Zs2aMp8LgEJDQyktLY0mTpxIzZo1I0EQ6MiRI0RUfsAym820fPly6t69O3l5eZFaraZ+/frR/v37q/wzXr16ldzd3alp06al/p3vJIoiWSyWUrWX9e9vtVqpQYMGBKDS89YHHLBcXHUDltVqtX8627dvn/3+EydOkLu7O3Xv3r1Eo7Ft2zZq0KABDR8+nGbNmkXz5s2jxx9/nABQ+/btKT8/v8T59Xq9/fwdO3akV199lV5//XV65plnSKVS0V9//UV5eXm0YMEC8vPzIz8/P1qwYIH9q7hxKc/atWsJAE2YMKHqv6TbcMBirLTb25FTp04RAFqxYkWJYwYPHkxNmzYls9lc5h9Yg8FAvXv3JgDUq1cvioqKooiICPL29iaFQkHR0dElji8OTaNGjaIGDRrQuHHjaM6cOdSlSxcCQOPGjbMfe+TIEZo4caI9lNzeZuTl5dlfv1+/fgSAevbsSa+++ipNnz6d2rRpQwDoP//5T5V+F48++igBoJ9++qlav8NBgwZRhw4d6Pnnn6eoqCiaNm0aBQYGEgCaO3duqeMB0H333Udt2rShHj160KxZs2j69On2XvmyApYoihQWFkYAqFOnTjRv3jyaPn06+fv7EwD67LPPqlTr22+/TQDo3XffrdbPeHvt5QUsb29vDlj/4IDl4oobxg4dOpRodG7/OnDgQInnpKSkUKNGjSggIIBu3LhB+fn5FBQURCqVii5evFji2IyMDCosLCz1ups3byYAtGjRohL3v/HGG/YAdOennvz8/BLnqskQYUREBAGgDRs2VOt5t78mDxEyVtKdH9Tuvfde6t69u/3xtLQ0cnd3p3nz5hFR2X9gFy5cSADo+eefJ1EU7fefO3fO3pul0Wjs9xcHLD8/P7py5Yr9fp1OR506dSJBECg9Pd1+f2VDhMWh4b333ivx+hqNhh544AHy9PQscb6ymM1m8vDwIIVCQXq9vsJj75SQkFDm+QYOHEju7u6UnJxc4rHiHriIiIhSbSVR2QFry5Yt9pBpNBrt96emplLTpk3Jw8OjzDru1L9/fwJAhw8fruqPV6p2HiKsHAcsF1fRcFvx12uvvVbqebt27SIANGDAABo7diwBoJUrV1b5dUVRtHdNF7NYLKRWq8nPz49yc3MrPUdNAtbgwYMJAB08eLDUY3/++WepcPndd9+Vek13d3f742+//TaNHDmS3N3dCQAtX768WvUwVhfcGbCWLl1KAOjs2bNERLRo0SICYP8AVtYf2Pbt25OHh0ep+TpERNOmTSMAtHXrVvt9xQFr/vz5pY4vfuz2XvaKApbVaqVGjRrRPffcUyJcFdu3b1+VerFu3bpFAKh58+YVHlcd0dHRBIA2bdpU4n4AZQ7FFisrYBWPHvzxxx+ljv/4448JAH3wwQeV1lTcS3h7sC32yy+/lGpH75zLCoAaN25sfzwqKooGDRpEAMjNzY127dpVaQ31AV9FWEcMGjQIBw8erPLx4eHhmDx5MtatWwcAeOaZZzBjxowyj42JicGaNWsQGxuLvLy8Epdw37hxw/79lStXoNFoMGDAADRq1KiGP0nFqIKNB2JjY/H++++XuG/atGkYOnRoifusVqv9ODc3NzRq1AhPPPEEXnnlFQwfPtzxRTPmYiZMmIC33noLGzZswP33349Nmzbh4YcfRteuXcs8XqPRICEhAV26dEHr1q1LPd6vXz+sWbMGf//9d6mrnXv16lXq+OJz5OfnV6neuLg45OXloWXLlqXaAADIysoCYGujakthYSGWLl2KPXv24Pr16ygqKirx+O1tZbF27dqVujqvIn/99RdUKhUeeuihUo8VX834999/V3qeitrRw4cP48MPPyxxn5eXF/r06VPivpycHPvvuni5mxEjRmDu3Lno27dvpTXUBxyw6rGwsDB7wHrllVfKPGbZsmV4/fXX0aRJEzz55JNo3bo1VCoVAGDFihUwGo32Y4sbw1atWtVazc2aNQMApKenl3ps8uTJmDx5MgDbmjn9+/cv8xxKpRIGg6HWamTM1TVt2hRDhgzBtm3bMHz4cFy7dg2vv/56ucdrNBoA//f+vFPxUjEFBQWlHvPz8yt1n0Jh+9NU3npcd8rNzQUAXLx4ERcvXiz3uDtDz50aN24MDw8P5OTkwGg0QqlUVun1TSYT+vXrh9jYWPTs2RMTJkxA48aNoVAokJSUhM2bN5doK4uV9/sqj0ajQZs2bcp8rKLfcVmve+XKFaSnpyM4OLjEY4sWLcKiRYsAAJs2bUJERESZ5wgODq7VwFoX8EKj9VRubi6mTp0KHx8fKJVKzJw5s1TjY7FY8MEHH6Bly5a4ePEivv76a3zyySd47733sGDBAphMphLHN2zYEEDZ4cdRHnnkEQDAkSNHau01GGPASy+9hLy8PEyaNAkqlQrPPfdcuceq1WoAQEZGRpmPF99ffJyjFZ83PDwcZJv6UubXxo0bKzyPQqHAQw89BLPZXK0FM/fu3YvY2FhMnjwZsbGxWLVqFRYtWoT33nsPTz31VLnPq+5K7Wq12iG/Y25HnYMDVj01ZcoUpKWl4YsvvsDixYtx9epVvPbaayWOyc7ORkFBAXr37o0mTZqUeOzs2bPQ6/Ul7gsODoZarcaZM2eQl5dXaQ3u7u5V/oRabNSoUfD19cXOnTsRHx9frecyxqpuyJAhaN68OdLT0xEeHl7hH261Wo327dvj2rVrZX7AOnbsGADc1a4S7u7uAMru1erSpQvUajXOnj0Ls9lc49cAgEmTJgEAPvroowqH0gDYe6WuX78OAGVOMTh+/Phd1XO7nj17Qq/X4/Tp06Ueq87veOLEiXBzc8PatWuRnZ3tsPpYSRyw6qEvv/wSMTExePbZZzFx4kS89tprGDRoENavX49du3bZj2vatClUKhViY2Oh0+ns9+fl5eHVV18tdV6FQoFp06ahoKAAr732WqmGsKCgAFqt1n7b398f2dnZ1Rqu8/f3x8cffwyj0YjBgwfjr7/+KvO4qs7dYIyVTaFQYN++fdi9e3epOTllmThxIsxmM956660SweTChQvYuHEj/Pz8MGLEiBrX4+/vDwBlboGlUCgwY8YMJCcn4/XXXy8zZF24cAGZmZmVvs6ECRPQt29fHD16FBERESgsLCx1TEZGBqZMmWKf9xoYGAgAOHHiRInjjh07hi+//LLyH66KJk6cCAB46623SvyM6enpWL58ORQKBZ5//vlKzxMcHIy5c+ciMzMTgwcPtgfEO3E7end4DlYdce3atQr36Cp+LC4uDrNnz0bbtm2xevVqALZu6k2bNqF79+6YOnUqHn74YbRp0wZubm54+eWXsWzZMvTo0QPDhg2DRqPBgQMHEBgYiJYtW5Z6nYULF+LUqVPYunUrTp06hcGDB0OpVCIhIQEHDx7EiRMn7J+wHn/8cZw9exbDhg1D37594enpiT59+pSaTHmnV155BUVFRXj77bdx//33o3fv3njggQfg6+uLnJwcXL58GcePH4dSqcSDDz5Ys18oYwwPPvhgld9DkZGR+P7777F161ZcvnwZTzzxBLKysrB9+3aYzWZs2bIFvr6+Na6lc+fOaNmyJb799lt4e3ujdevWEAQBM2bMgJ+fH95//33Exsbi888/x/fff4/Q0FA0adIE6enpOH/+PM6dO4eTJ0+iadOmFb6OQqHAnj17MHr0aGzevBn79u3Dk08+iXbt2sFkMuHSpUs4evQozGazfcL+sGHDEBQUhE8//RQXLlzAfffdh7i4OOzfvx8jRoxAdHR0jX/u202YMAExMTHYu3cvunfvjqFDh6KoqAg7duxATk4Oli1bhvbt21fpXIsXL4bZbMZnn32G4OBghIaGonv37vatcv7++2+cPXsWarW61HZjrIokunqROUhVlmko/mc2Go3Uq1cvcnNzo2PHjpU613fffVdqBWOTyUQffvghdezYkZRKJbVt25bmzp1LhYWF5S6zYDAYaOnSpRQSEkIqlYp8fHyoa9euNG/ePPuigEREhYWFNGXKFGrRogW5ublVexuMuLg4mjVrFnXr1o3UajUpFAoKCAigvn370sKFCyklJaXUcypaaJSx+qq6CxajnHWQtFotvfvuu9SpUyfy9PSkhg0b0uDBg0td5k9UciX3O23cuJEA0MaNG0vcf+rUKQoNDSVfX98yV3K3WCy0Zs0aevTRR0mtVtvbrKeeeopWrVpFWq22Sj8fkW0pml27dtGIESOoZcuW5OnpSd7e3nTffffRrFmz6NKlSyWOT0hIoPDwcGrSpAl5e3vTgw8+SN9++225y0vgn/WsylPRSu5Lly6lbt26kVKpJF9fXwoNDa3WavW3O3v2LE2ePNm+Ar2Hhwc1a9aMBgwYQMuXLy9zGYny/v1ZSQJRJYPMjDHGGGOsWngOFmOMMcaYg3HAYowxxhhzMA5YjDHGGGMOxgGLMcYYY8zBOGAxxhhjjDkYByzGGGOMMQdzSMDS6XSlVvtmjDFXwu0YY8yRHBKwrly5gvvvv9/ldtauyq7jcuSKdXPNzuOqdUutKu2Yq/xuXaVOgGutLVyr41W3zno9RFjdjYblwhXr5pqdx1XrdgWu8rt1lToBrrW2cK2OV90663XAYowxxhirDRywGGOMMcYcjAMWY4wxxpiDccBijDHGGHMwDliMMcYYYw7GAYsxxhhjzME4YDHGGGOMORgHLMYYY4wxB6s0YBmNRsycORMdO3bEvffei/HjxzujLsZYPWG0ilKXwBhjZSIiZOkJf2YRbuqFaj1XUdkBb775Jtzc3HD16lUIgoCbN2/WuFDGGLudwWJFikaPTv4+UpfCGGMAgAwd4ZqGkKwFUrUE/T8LuD/h78CAVVRUhI0bNyItLQ2CYDtxixYtalYxY4zdxmCxIqVADwuR1KUwxuoxIkJqEXAln3Aln5BrdMx5KwxY169fR+PGjbFo0SIcOnQIKpUK7733Hp544okyj9dqtdBoNPbbSqUSSqXSMZUyxuoMvdnWc2XlcMUYk8hNHeHvHMLFPILW7PjzVxiwzGYzEhIS0LVrVyxevBjnzp3DgAEDcOnSJTRp0qTU8aGhoSVuR0ZGIioqyrEVO1BeXp7UJdSIK9bNNTuP3Os2WAkZRivEEtnKV6pyGGP1iM5COJ9L+CubcEtfu69VYcAKDAyEm5sbnn/+eQBAjx490K5dO1y8eBH9+vUrdfyxY8cQEhJiv+0KPVj+/v5Sl1Ajrlg31+w8cq27yGxBpsYApRf3XDHGnOeWjvBbBuFSHsHqpOanwoAVEBCAJ554Aj/++COGDBmC5ORkJCYmIjg4uMzjfXx8oFara6VQxphr05osSCvU39FzxRhjtSdRQziRQbiucX7DU+lVhKtXr8ZLL72EqKgouLu7Y+3atTzRnTFWLYUmC9I5XDHGnORyni1YpRdJ1+hUGrDat2+Po0ePOqEUxlhdpDGakV5oAGcrxlhtS9QQDt2QNlgV45XcGatFMTEx6NGjB1QqFXr06IGYmBipS3KqfAOHK8ZY7bulI3wVL2JzvCiLcAVUoQeLMVYzMTExCA8PhyAIICKcP38e4eHhiI6ORlhYmNTl1bo8gxm3tByuGGO1R2chHEy1XRkot7aGe7AYqyXvv/++PVwBtsXsBEHAwoULJa6s9uXqTRyuGGO1KqmQsPqSiP/JMFwB3IPFWK25evWqPVwVIyLExcVJVJFzZOtMyNQ5aClkxhi7g0iEX28Sjt2UZ7Aqxj1YjNWSTp062beYKiYIQrnLnNQFWTojhyvGWK3RmAhb4glHZR6uAA5YjNWaBQsW2IcFAdiHCxcsWCBxZWW72wn5GUVGZOlMtVQdY6y+S9QQVl8WkVQo92hlwwGLsVoSFhaG6OhodO/eHV5eXujevTtiYmIwcuRIqUsrpXhC/vnz52EwGOwT8qsSsogIN7UG5Og5XDHGase5HMJX10ToLFJXUnU8B4uxWhQWFuYSVwxWNCG/ovpt4cqIfGMt7JTKGGMAjt0UceSGa/Ra3Y4DFmOsRhPyiQjpWgM0Rhf6SMkYcxlWkfBdCuHvHNcLVwAPETLGUP0J+SIR0go5XDHGaofRSvj6muuGK4ADFmMM1ZuQLxIhVaNHoYnDFWPM8bRmwsY4EQkuMpm9PBywGGNVnpBvFQkpBXoUma0SVcoYq8vyjYQNcSJu6aWu5O7xHCzGGIDKJ+RbREKKRgeDRXRiVYyx+iJTT9gaL6KwjlwzwwGLMVYpsygipUAPo5XDFWPM8dK0hK+vidDXoc5xDliMsQqZrCJSNHqYZBSu8vPz0a9fP/ttnU6HhIQEZGZmwt/f337/0aNHMWTIEHTq1Ml+38mTJ6FSqZxZLmOsAklaAQdzRJjk08Q4BAcsxli5jFYRKQU6mEV5TTZt2LAh/v77b/vtpUuX4tixYyXCVbGuXbvi7NmzziyPMVZFcfmEXake8PSSuhLH40nujLEyGSxWJMswXJVl48aNmDRpktRlMMaq4VIeYUeCCKv8m5ga4YDFGCtFb7YiuUAPiwuEq5MnTyInJwdDhw4t8/G4uDj06tULDz74IFauXFnp+bRaLTQajf3LaOTNqxlztP/lEHYl1t1wBfAQIWPsDkVmC9I0BljJNVq+DRs24IUXXoBCUbo569WrF9LS0uDn54e0tDQMGTIEAQEBGDNmTLnnCw0NLXE7MjISUVFRAIC8vDzHFl9LXKVOgGutLXKu9VyeG368qUBxC6PXGyStp6o0Gg1yc8ufKHbnFAUOWIwxO63JgrRCPVyg4woAUFRUhO3bt+P06dNlPq5Wq+3ft27dGs899xyOHz9eYcA6duwYQkJC7LeVSiWUSqX9dlnzvOTIVeoEuNbaIsdaz2SJ+DWf4HXHdSaucOGJWq2Gv3/DKh/PQ4SMMQCAxmhGqsZ1whUA7Ny5E927d0fnzp3LfPzmzZsQRdsnzsLCQuzfvx89e/as8Jw+Pj5Qq9X2r9vDFWOs5n69KeL7FIILNTF3hQMWYwz5BjPSCw0u1/CtX7++1OT2yZMnY9++fQCA6OhodOvWDT169EDv3r0xcOBARERESFEqY/UWEeFgqojDN1ythbk7PETIWD2Xqzcho8jocuEKAI4fP17qvnXr1tm/nzlzJmbOnOnMkhhjtxGJsCeJ8L9cV2xh7g4HLMbqsWydCZk6vkqOMeZ4ZpGwI4EQX1D/whXAAYuxeiuzyIhsvUnqMhhjdZDOQth2jZBaVD/DFcABi7F6h4iQUWRErqGO7KjKGJOV9CLbAqIF9fzzGwcsxuoRIsJNrRH5Rg5XjDHHO5sl4mAqwVJ/O67s+CpCxuoJIkJ6oUFW4ern/fswMvQRl1gDhzFWPrNI2JMkYn8Kh6ti3IPFWD0gEiGt0ACtySJ1KXY/79+H2RHjIQgCyEVWjWeMlZZrsA0J3tJLXYm8cMBirI6zioTUQj10ZqvUpZSwcsliDleMubizWSJ+SiOYyt9Bpt7igMVYHWYRCakaPfQWeYUrAEi6Hs/hijEXpTUT9ibX3yUYqoLnYDFWR5lFEckFOqeEq+K5VD1bN8HI0Efw8/59lT4nqENHCIJQ67UxxhzrYi5h5SWRw1UlOGAxVgeZrCKSC/QwWmu/3754LlX85YswGY2Iv3wRsyPGVxqyXn7jTRARhyzGXITRSohJFLEzUYROPtM5ZYsDFmN1jEkkJBfoYHJCuAJKz6UqDk2rli6u8HkDhw7Hio1foVPXe+Hl5eWMUhljNXSjiLDmslgvt7ypKQ5YjNUheosVtwxWmEXnNYJlzaUiIiRei6/0uQOHDkfM0d+h1/PlR4zJ1akMEevjROTyrlrVwgGLsTpCZ7YgpUAPq5M/YJY1l0oQBLS7p5NzC2GMOZRtuxsRB9PI6e1KXcABi7E6oNBkQYpGD6sEV+XdOZeqeLjw5TfedHotjDHHyNDZhgTjeCJ7jXHAYszFFRjNSNPo4cRRwRJun0vlqVSiU9d78dmmrzHg6WHSFMQYuysJGsKGq7yX4N3idbAYc2F5BhNuaY2Q+jPmwKHDMXDocImrYIzdrXM5hH3JIg8JOgAHLMZcVLbOhEwdzzpljDnGrzdFHL7BycpROGAx5oIyiozI0XP/PWPs7olE+D6F8Gc2hytH4oDFmAshItwsMiLfYJa6FMZYHSASITqRcDGPw5WjVXmS+/vvvw9BEHDhwoXarIcxVg6RCOmFhjoTrsxmMw4f+B6zJo6TuhTG6iWrSNiRwOGqtlSpBys2NhanTp1C27Zta7sexlgZRLJt2lxklt+mzdV1Iy0VOzdvRMw3W5GdmYF7e/SUuiTG6h3LP+HqKi/DUGsqDVhGoxGvvPIKvvnmG/Tv398ZNTHGbmMRbeHKGZs21xYiwpnfT2DrmpU4+uMBeDfwwbDRz2LUhBfR+b5uUpfHWL1iFgnfXidc13C4qk2VBqz58+dj/PjxaNeuXaUn02q10Gg09ttKpRJKpfLuKmSsHjNbRaRonLNpc20wm804uCcam1f9F5fPn8M9nbvg3U+X4+nwMWjg4yN1eYzVOyYrYdt1QmIhh6vaVmHAOnnyJM6cOYPFiyvetLVYaGhoiduRkZGIioqqeXW1LC8vT+oSasQV6+aaq88kEjIMVliq2Q4a9IbaKag6NRj02Ld9G7auWYmbaan4f6H98cVXO/Bw38fsK76X3H/QV5pCGatHzCLhm+uEJA5XTlFhwDp27BiuXLli771KS0vDoEGDsG7dOgwePLjM40NCQuy3XaEHy9/fX+oSasQV6+aaq05vtiJbo4eHF8GjBs9XqVQOr6kqDHo9dmzegPX/WYHc7Cw8NSIM//1qO4LvvU+SehhjNpZ/hgU5XDlPhQHrzTffxJtv/t9+YkFBQdi/fz/uu6/sxtLHxwdqtdqxFTJWz2hNFqQVGiBKsK9gTZlMJuzauglr/70UudlZGD7mOUyZPQ+B7TtIXRpj9V7x1YI858q5eB0sxmSkwGjGjUKD5FvfVJUoivhh9y58sXgR0pKTMWz0WEyfF8nBijGZEImwK5GvFpRCtQJWUlJSLZXBGMvVm5BRJP2+glV19vff8On8t3Hx3F94fPDT+M+Wb9GxS1epy2KM/YOIEJNIuJzvKq1K3cI9WIzJQJbOiCyda2x9k56SjCUL3sHP+/fi3pCe2LLvIO7/f49IXRZj7DZEwN5kwgVeRFQyHLAYkxAR4VaREXkusDq7Qa/Hhv9+hnWfLYdfw0ZYvHItng4fAze3Km8IwRhzAiLCgZsKXDPWPFxdPrwbJ9YtQnZyHAICg9Fn8jvo8vhIB1ZZ93HLyJhEiAjpWoNLhKsThw9hxGMPY83yJRg/ZQb2n/wTw0aP5XDFmMwQEfanEM7n1/y9efnwbuyIHI0OxnREPRSIDsZ07IgcjcuHdzuw0rqPe7AYk4BIhDSNHlqZb32TlZGBT96JwoE9MXioz2NYtS3hQTMRAAAgAElEQVQa7e7pKHVZjLFy/JBK+DP77oYFT6xbhL6BAfj6mV4QBAFTegZi3J5Y/LZ+EfdiVQN//GTMySwiIblA3uGKiLDn268xvM+DOHX8GD7+7xpsiPmOwxVjMnYgVcSZrLufc5WdHId+bRvbFwUWBAH9AxsjKynurs9dn3DAYsyJzFYRyQU6We8reDM9DdOeDcO/Xp2B0IGD8N1vZzF8zHP2xpYxJj8/por4I9MxE9oDAoNxNCUH9M9afESEI8k5aBIU7JDz1xc8RMiYkxitIlIKdDCL8ryqh4iwd/s3+PjtKDTw8cGqb3bisYGDpC6LMVYBIsLBNHJYuAKAPpPfwY7I0Ri3Jxb9AxvjSHIOTqRkY8yS1Q57jbvlCpPwuQeLMScwWKxIlnG4ys3OxqyJ4/CvV2fg8cFDsOf4KQ5XjMlc8YR2R4YrAOjy+EiM+XQnElWt8MnpZCSqWmHMkl3o0n+EQ1+nplxlEj73YDFWy3RmK1I1elhluvXN8UM/4V+zXoYoWvHZpq8x4OlhUpfEGKsEEWFvMuHvnNppV7o8PlJ2PULFXGUSPvdgMVaLtCYLUjQ6WYYrk9GIj/8VhenPjUKXbt2x+9gpDleMuQCRCNGJtReu5M5VJuFzwGKslmiMZqRq9HD2qODhA99jZOgj6Nm6CUaGPoKf9+8rdUzS9XiMGzwA2zetx9sffYrV30ajSbNmzi2UMVZtFpGwM6F+r9DuKpPwOWAxVgvyDWakS7Bp88/79yFy2kuIv3wRJqMR8ZcvYnbE+BIh67ud32LU449BryvCtoOH8fyU6XyFIGMuoNBE2HSV9xbsM/kdHE/Oxrg9sVgbm4Rxe2JxIiUbj05+V+rSSuCAxZiD5epNuKl1frgCgJVLFkMQhBKf7ARBwKqli2HQ67Fg7iy8+fJUDBw6DDsO/You3bpLUCVjrLrStIS1V0SkFdXvcAXIfxJ+MZ7kzpgDZetMyNQZJXv9pOvx9nBVjIiQEH8V4wYPQNL1eCz89xcIe34C91ox5iL+yiZ8nyLCwtnKTs6T8ItxDxZjDpJZZJQ0XAFAUIeOpYKTIAiwWiww6HXYdvAwwse/UCfCVVBQEDp37oyQkBCEhIRg+/btZR63aNEidOjQAR06dMC778prCIGxiohEOJAqYm8yhytXxD1YjDnArSIDcvXSb9r88htvYnbE+FLDhN163Y+1O3bDV+0ncYWOtWvXLtx3333lPv7rr79i27Zt+N///geFQoFHH30Uffr0waBBvMYXk7dsA2FPEvGQoAvjHizG7gIR4aZWHuEKAAYOHY5P12xAh+DOENxsb+9ho8fimwO/1LlwVRXbt2/Hiy++iAYNGkCpVOKll17Ctm3bpC6LsXIREU5liFhzmedbuToOWIzVEBHhhtaIPIM8wlWxtu3awWjQQ+3nhy937sHilWvrxJBgWZ5//nl069YNkydPRlZWVqnHU1JSEBgYaL8dFBSElJSUCs+p1Wqh0WjsX0ajtMO+rP7INxK2xNu2vjGLUlfD7hYPETJWA0SEdK0BGqNF6lJK+Om7vXh75jS0bdcB63btQ+vAIKlLqjW//vor2rZtC7PZjHfeeQcTJ07EDz/8UOq428PlnRcAlCU0NLTE7cjISERFRQEA8vLy7rJq53CVOgGuFQBEAs7lu+FohgImBwUrvd7gmBM5gavUqtFokJtb/j+Qv79/idscsBirJpEI6YUGFJrkE65EUcTKJR9j1dJPMHDocHz0xRp4N2ggdVm1qm3btgAADw8PzJ49G506dSrzmKSkJPvt5ORk+/PKc+zYMYSEhNhvK5VKKJVK++07G1G5cpU6gfpda1w+4VC6iCwD4K4EVA48t0rlyLPVLleoVa1Ww9+/YZWP5yFCxqpBJEKaRi+rcFWk1WLOSxOwetmneO3t+fjov2vrfLgqKipCfn6+/fa2bdvQs2fPUseNHj0amzdvRlFREYxGIzZs2ICxY8dWeG4fHx+o1Wr71+3hijFHSdMSNsSJ2HbdFq5Y3cM9WIxVkUiEVI0eRWar1KXY3UhNwSsTxiItKQmfb9mGx58aAr1eL3VZtS4jIwPh4eGwWq0gIrRv3x5btmwBAAwZMgQLFy7EAw88gH79+mHMmDHo1q0bAGDs2LF46qmnpCyd1XOZesKRG7wae33AAYuxKrCKtnCls8gnXP11+g+89uI4eKm88fUPP6NT13ulLslp2rdvj7/++qvMx+6chzV//nzMnz/fGWUxVq5MPeHoP8GKo1X9wEOEjFXCKhJSZBau9u/ajoiRTyOoQ0d8++ORehWuGHMlmXrCjusiVl0ScYnDlUu6fHg3vhx3P57s2gL39wxBTExMlZ7HAYuxClhEQrJGB71MwpUoivj84w8QNWMKng4bjXW79sI/IEDqshhjd0jVEr7lYOXyLh/ejR2Ro9HBmI43HgyET94NhIeHVylk8RAhY+WwECG5QAejVR4L0hj0evzr1Rk4uDcGc959H5NenV1n17dizBUREa4WAL9lEFK0HKnqghPrFqFvYAC+fqYXBEHAlJ6BeH5vLD5a9AHCwsIqfC4HLMbKYLaKuGUQoVDKI1xlZ2bi1RfG4uqli1ix8SsMHDpc6pIYY/8wi4TzucDJDL4iUO4uH96NE+sWITs5DgGBwegz+Z0KN43OTo7DpIcC7R9mBUFAv7aNsezslUpfi4cIGbuDySoiWaOHWZTHJ9D4y5cwdlB/3ExLw+Z9BzhcMSYTeUbCT2kilv9PxL5kDldyd/twX9RDgehgTMeOyNG4fHh3uc8JCAzG0ZScEnu7Hk3JQZfOnSt9Pe7BYuw2JquI5AKdbMLV70cPY85LL6Bl27ZY+fUOtGjVWuqSGKvXiAgJhcAfmYT4Ap5b5UrKGu4btycWv61fVG4vVp/J72BH5GiM2xOL/oGNcTQlB8eTsxHz77WVvh73YDH2D6PFKqtwtXPrJkwfG46eDz2Mr/b/yOGKMQnpLITTOW74z0URW+NFXOVw5XKyk+PQr23jEsN9/QMbIysprtzndHl8JMZ8uhOJqlZYciYZOv/WiImJwciR5Q8rFuMeLMYAGCxWpGj0sMggXImiiBWL3sf6//wbYyMm462PPoVCwW9VxqSQqiWczSJczCMU6hRwgR1dWDlsw33pmNLTNqeKiHAkOQdNgoIrfF6Xx0eiy+Mj8YR/Ifq2q/pWOdxqs3pPb7EipUAPaxU2Aq5tRoMBb8+cjh/37UbkBx/hhWmv8JWCjDlZlp5wPpdwIY+Qa5S6GuYodw73HUnOwYmUbIxZsrpWXo8DFqvXdGYrUjXyCFd5OTmYOWEsrlz4H/69YStPZmfMiXIMhCv5tmB1q+7vNlUvFQ/3/bZ+EU6fjkOToGCMWbIaXfqPqJXX44DF6q0ikwWphQaIMghXyQnXMeO5USjUFGDj7v3ofv+DUpfEWJ1mFQlJWiC+gHC1gHuq6ovi4b6qKDQRUoqAFC0hRQtsdlMhIYiqPKrAAYvVS1qTBWmFeshgyhXOnT2NV8Y/i4aN/PHNgV/QJqid1CUxViflGAjXNbarABM0BJM8lrljMmERCUmFwJUCwnUNkG+y3e+vBNo2AEa3NkEQlFU+HwcsVu9ojGakFxpkcQXQ4QPf441pL6FLtx74Yus2NPRvLHVJjNUZFpFwTWPrpUooJORxLxW7g9FKiC+whar4AsAkAo08gU5+QKCPgLY+gI+Hrceqr3/1tkzjgMXqlXyDGTe18ghX2zZ8iY/eegMDnh6OxSvXQunlJXVJjLk8q2jrfbiYR7hSQDDKYxtRJiNEtuHhv7IJl/MBCwEtVMAjzQR0aQg08YJDLi7igMXqjVy9CRlFRsnDFRHhsw8X4svPlmHCtJcRufAjuLnxknSM1ZRIhMRCW6i6nEfQc6hiZcg3Ev7OBc7lEPJNtqG/0BYC7msENFQ6/mptDlisXsjWGZGpM0ldBsxmM+bPnol9O7bhjYUf4cUZM6UuiTGXRERI1gIXcgmX8wlFFqkrYnKVVkT4PcPWW+XpBnRtBPRsLKBNA8f0VJWHAxar8zKKjMjRSx+uirRazHlpAv448Ss+XbMeT4eNlrokxlxOvpHwZzbh7xxCoVnqaphcERGuFgC/Z9quAPRXAk+3EdDdH/B0d87aghywWJ1FRLhZZES+QfpWODc7GzPGjUJifDzWfBuN3o/1k7okxlyG+M8fy7NZtqsApR7mZ/JFRLiUDxy9Qcg2Am0aAM+2F9DJD3Bz8qLNHLBYnURESC80QGNyzLjBz/v3YeWSxUi6Ho+gDh3x8htvVnkh0BupKZgyegQKNQXYtOd7dO0R4pCaGKvrzCLhTBbhVAZBI/3nJCZzaUWEn9IIqUXAPWpgeKCANj7S7YRRYcAyGAwYO3YsLl26BG9vbzRv3hyrV69GUFCQk8pjrPpEIqRp9NCaHTPT9ef9+zA7Yrx976r4yxcxO2I8Vmz8qtKQde3KZUwdMxKeSk9s3f8TAtt3cEhNUhMAeLq7QenuBk93N3i4u8HDTYC7mwA3CHATAEEAiAACYCWCVSRYiWC2ijBZCSarCINVlMVCr0xerCIhNofw600eBmSVyzMSfrlBuJgHNFMBE+4R0F4t/RZjlfZgTZ06FYMHD4YgCPjiiy8wdepU/PTTT86ojbFqs4qEFI0eeovjLiNauWSxPVwBtt4xQRCwauniCgPWubOnMf25UWjeqjXWfhuDJs2bO6wmZxMAqDzc0cDDHd4Kd6g83B3S3U5kC1qMAbYPR+dygGM3Rfsij4yVRyTCiVvAr7cIKgUwvK2AHo2dPxRYngoDlpeXF4YMGWK/3bt3b6xYsaLWi2KsJsxWESkaPYwO/oOddD3eHq6KERESr8WX+5zfjx7GrInj0Pm+7lj5zQ6o/aq+A7tcuAmAj4cCvkoFfDwUcHdzfKMlCAKUCneHn5e5ngQN4WCaiEzeB5BVgcZEiEmyTWB/pBnwWHPBaZPXq6pac7A+//xzDBs2rNzHtVotNBqN/bZSqYRSWfVl5RmrKaNFRIpGB3Mt7H0T1KEj4i9fLBGyBEFAu3s6lXn8T9/txRvTXsL/e6wf/r1hK1Te3g6vqTapFO5o6OUBtWfthCrGbldgAn6+LuJyPg8Vs6q5kk/Yl0zwcANe6CggyFee7VSVA9ZHH32E+Ph4rF69utxjQkNDS9yOjIxEVFRUzaurZXl5eVKXUCOuWHdt1mywEjKNVlgd3D4b9AYAwOTX5iJy2kv2YcLi/05+bS70+pIft/ft2IZFkXPxxNPDsPDfXwCCUOqY2lZcd3UIABooBKgVblCKAkhnQIHO8bWVx9/f33kvxmTBLBJ+u0X4OckTHkoOV6xyZtE2if1sNtDZzzaJXaWQZ7gCqhiwli5dipiYGBw6dAjeFXwaP3bsGEJC/u8KKVfowXLVht0V666NmrUmCzIL9fCspV1mVCoVng4bBU9PT6xauhiJ1+LR7p6OePmNtzDg6ZK9uVvXrMTid97E6Bci8O6ny+HuLt3Ql0qlqtJxbgLQUOmBxipPeLjzavLMOeLyCQdSbfOsLCLgIXVBTPY0JsK264Rsg209q/sDaneRUEeoNGAtX74c27Ztw6FDh9CwYcXzSHx8fKBWqx1WHGMVcea+ggOHDi93QjsRYc3yJfjP4kWImPka5s1fKPs3vgCgkZcHGnt7woO36WFOkm8kHEglxBVwjxWruls6wjfXCW4CMDlYQDNvebevxSoMWGlpaZg3bx7at2+P/v37A7D1Sv3xxx9OKY6x8shl6xsiwvIPFmDDf1Zg1lvvYuqc12UfrvyUCjTxVsKTe6yYk1hFwslMwrGbBDNfNMqqIb6AsDOREOAFPNdBgK+HvNvX21UYsFq3bl3q6inGpEREuFVkRJ4MVmcXRREfvfUGtm34ElEffIwXpr8idUkVUinc0ayBEt4efNUec540LWFvsois6k8NZPXcmSxbj2cnPyAsSH5XCVaGV3JnLkP8Z3X2Qgetzn43rFYr5s+eib3bv8F7yz/H6AkvSl1SudwFAU0beKKh0kP2vWus7rCIhKM3bZvs1sLFvawOuXpsL05vWYLs5DgEBAbj0Un/QmrwCJzKBB5uCjzZSpDN2lbVwQGLuQSLKCJVY3DoAqI1ZTab8dbLU/DTd3uxeOVaDB31rNQllctPqUCzBkooeJ4Vc6IbRYQ9ybymFavc5cO7seed8egbGIBJDwXiaEo6dkaNAaZvx1PDw/BwU9cLVsU4YDHZM1pFpGr0sljx22Q0Yu7kiTj+y89Yvn5LqSsJ5cLDTUBzHy/4evJbnDmP9Z9eq9+414pV0Yl1i9A3MABfP9MLgiBgSs9APLcnFlcOfYSHJ4dLXd5d4daXyZrObEWaRg+LDOYC6nU6vPbi8zh78jd8sWUb+g54UuqSyuSjENC+YQNeJJQ5VYaOsDtJxC3utWLVkJ0ch0kPBdqnLwiCgMcDG+PM6TiJK7t7PG7AZKvAaEaKRieLcFWk1WLGuNGI/eMUVm3bJctwpXAT0MZXhSZKdw5XzGls+8GJWHuFwxWrvoDAYBxNySmx1+uR5Bw0CQqWuLK7xz1YTJaydSZk6YxOWeOqMoWaAkwfG474y5exdsdu9Hq4t9QlleLjqUBLH9tcq9wiqath9UWugbA7iZBaJId3KnNFvV78F/a/PQbP7YnF44GNcSQ5BydSsjFmSfm7xrgK7sFiskJEuFFoQKZMwlV+Xi4mhQ9HwtWrWB+zT3bhyk0AmjdQoo2vF09kZ071VzZh9WWRwxWrsXwj4beWI6B6ZTuuebXCJ6eTkahqhTFLdqFL/xFSl3fXuAeLyYZVJKQV6lFklv5KQQDIy8nGzPFjkXnrBjbs3o8u3bpLXVIJSnc3tPL1gpeC17VizmMRCT+kEmKzOVixmss1EDbHE9zdgBfCh6B5hGtPaC8LBywmC0ariDSNHkYZXCkIAFm3bmHamDAU5Odh054fcE/nLlKXVIKfUoHmDbx4rhVzqnwjYUcC4YaOwxWruSw9Ycs1gpcbMKGjAA95fKZ2OA5YTHJFJgvSCg2wymAyOwDcTE/DpLBh0Ot02PLdAQR16Ch1SXYCgGYNlPBXeUpdCqtnrmsI0YkidNKv88tc2C0dYes1go/CFq58PATo6+jFERywmKTyDGbcctKGzVWRlpyEl8KGgYiwduceWYUrDzcBrX1VUPFWN8yJiAjHbxGO3CDZvE+Za7pRZAtXjZTA+HsEeCvqdg88BywmCSJChs6IXL30ewoWS75+DS+FDYOHpyc27N6PRo0DpC7JzsfDHS19VVDwkCBzIq2ZEJNISCjkaMXuTmIh4dvrhKYq4PkOArzqeLgCOGAxCVhFQnqhHlqZTGYHgGtxVzApbBjUfn5YH/MdmjZvAb1M+q0DVJ5o4u3J+wgyp0rUEKKTRGjl8xmIuajLeYToJEKQLzCmnett2lxTHLCYU8ltMjsAXLlwHpNHDUeT5i2wbudeNG7SROqSAABugoCWPkqolR5Sl8LqESLbdje/3uQhQXb3YrMJ+1MIXRsBIwOFenVhDi+cw5xGa7IgKV8nq3B1/q8/ETHyabRs3RYbd++XTbjydHdDkJ+Kw1U5DAYDRowYgU6dOiEkJARPPfUUkpKSSh139OhReHt7IyQkxP4ll55JOcr559L5YxyumAP8dovwXQrh/gAgLKh+hSuAe7CYkxSYRRg0elk12n+e/B0zxo1Gxy5dsfrbXfBV+0ldEgCggYc7WvF8q0pNnToVgwcPhiAI+OKLLzB16lT89NNPpY7r2rUrzp49K0GFrsMs2nqsfs8gWOX0JmUuSSTCz+mEU5nAY82Bfi2EejnFgXuwWK0SiZBeaECuSZRVuDp57AimjQ3DvSE9sXbHbtmEK38vD7RVc7iqjJeXF4YMGWJvtHv37o2EhASJq3JNl/IIX1wUcfwWhyt297RmwtZ4wh+ZwODWAvq3dKuX4QrgHixWi8xWEWmFBugt8pnMDgBHfjyAOS9NwMN9HsNnm76Gl0oldUkQADT3UaKRF69vVROff/45hg0bVuZjcXFx6NWrF9zd3REREYGXX365wnNptVpoNBr7baVSCaVS6dB65SC9yLb0wjUNpyrmGKlaws5EAhEwsaOAQN/6GayKccBitUJntiKtUA+LKK/G+8DuaLz58hT0GzQES9ash6cM/nC6CwJa+3qhgSe/HWvio48+Qnx8PFavLr05bK9evZCWlgY/Pz+kpaVhyJAhCAgIwJgxY8o9X2hoaInbkZGRiIqKAgDk5eU5tvhaUl6dVgKuaNwQm+uOG3p5/PHT6w1Sl1BlXGvZiIDYPDccyVCgpYowvJUZPgpUeQFRV/m9ajQa5OaWP4fY39+/xG1u0ZnDyW3x0GIxX2/F/DkzMWz0s/jgs5VQKKT/39/T3Q1tfFVQKni0viaWLl2KmJgYHDp0CN7e3qUeV6vV9u9bt26N5557DsePH68wYB07dgwhISH223f2YN3ZiMrV7XUWmgh/ZhPOZpN92QUZdNzaqeRUTCW41pIMFtvelBfygN5NgQGt3OAuVL9tdYXfq1qthr9/wyofL/1fGFZnEBFuFRmRZ5DfwjlbVv8Xn7z7Fp59cRLe+WQZ3NykDzTeHu5o46uqd1fWOMry5cuxbds2HDp0CA0blt3o3bx5E82aNYObmxsKCwuxf/9+TJo0qcLz+vj4lAhmroqIcE0D/JlNuFpAkFlnMqsDEjSEvckEoxUIDxJwnz+3ZbfjgMUcwiyKSNPIb74VEWHVsk/w308+wkuvzsbcd9+XxYRLP6UCLX28ZFGLK0pLS8O8efPQvn179O/fH4Ctp+mPP/7A5MmTMXz4cAwfPhzR0dFYtWoVFAoFLBYLRo8ejYiICImrr10aE+G3LHdcSxdRYJK6GlYXmUXCoXTC6SygnS/wTKAAP09uy+7EAYvdNZ3Ztlmz3OZbERE+nf82tqz+L2a99S6mznldFoEmQOWJpg2kn/vlylq3bg0qZ3PwdevW2b+fOXMmZs6c6ayyJGMVCXEFwF85hOsaQpHOXVZDgKzuSC8i7E4iFJiAp1oLeKgJZNGuyhEHLHZXcvUmZBQZZTffymKx4L25s7B721f41+KlGDdpqtQl/XOloBcaefHiocwxMnSEv3II/8sl6CxSV8PqMpEIJ24BR28SWngD07oICPDiYFURDlisRkQi3NQaUGCUX6tuMhrxxrRJOHLwe3yy6ksMHfWs1CXB7Z8rBX34SkF2l0QiXM4D/sgipGjl9tGGyc3lw7txYt0iZCfHISAwGH0mv4Muj4+s1jk0JluvVZIW6NscCG0hwJ17rSrFrT2rNtM/+wkaZLTlTbEibSFefWEc/j7zBz7b/A36DxosdUlQuAlo46uCysNd6lKYCysy264CPJtFKJTfdSRMhi4f3o0dkaPRNzAAkx4KxNGUdOyIHI0xn+6scsiKy7dNZFe42da2Cqrna1tVBwcsVi2FJgtuFBpgLWf+i5TycnIwbWwYkq9fx9rtu/HAI49KXZJtGQa1Ckp36a9aZK7plo5wKpNwPpdXWmfVc2LdIvQNDMDXz/SCIAiY0jMQ4/bE4rf1iyoNWBbRtt3N6Sygk59tIru3gsNVdXDAYlVCRMjSmZCjN8luvhUA3EhNwdQxI6EpyMemvT+gS7fuUpcEL4XbP9vecLhi1UNEuFoAnMokJBbK8R3HXEF2chwmPRRon4QuCAL6BzbG6dNxFT5PbyF8e52QruOJ7HeDAxarlEUUkV5oQJFZXkswFLt25TKmjhkJD09PbN3/EwLbd5C6JF7jitWIyUo4l0s4lUHIMUpdDXN1AYHBOJqSjik9bSGLiHAkOQdNgoLLfU6+kfD1dUKR2TYk2MaH27Ca4oDFKqQzW5FeqIdZZkswFIv94xReGT8GLVq1wZpvo9GkeXOpS4KvpwKtfL3gxp/4WBVpTITTWYQ/swh6eX6OYS6oz+R3sCNyNMbtiUX/wMY4kpyDEynZGLOk9LZSgG04+uvrBIUATAoW0JivErwrPHbBypWrNyG5QCfbcHX44A+YPGo4grveh837fpBFuPJTeqA1hytWRTeKCNGJIj67IOLELQ5XzLG6PD4SYz7diURVK3xyOhmJqlYYs2QXuvQfUerY6xrCxqsEXw8OV47CPVisFDkvwVBs19bNeP/11/DEkGH4ZNWXUHp5SV0S/L080KyBkucqsEpd1xBO3OL5Vaz2dXl8ZKUT2v+XS9ibRGivBka3E+Dpzm2YI3DAYiUYLVakFRpglOESDIBt8u/KJR9j5ZLFGBsxGW9/vATu7tIvf8Crs7PKiES4lAf8liHipk7qahiz+SubsC+FENIYGNZW4N53B+KAxewKjGbc1BohynAJBgAwm8344I05iP56C2b/awEmvzZXFr1FzbyVaOztKXUZTKYsIuHvHMLvGYRcnrjOZOTPbML+FML9AcDTbQRZtKd1CQcsBiLCrSIj8gzyXb2wSKvF61NexO9HD+Pj/67B8DHPSV3SP1vfKNHIi8MVK81otS0KejKToJXvW4vVU2eyCD+kEh5sAgxuzeGqNvAk9zomJiYGPXr0gEqlQo8ePRATE1Ph8SariKQCvazDVdatW3jxmSH489RJrNq2SzbhqqWPF4crVorOQvglXcS/z4v4OZ3DFZOfPzJt4erhphyuahP3YNUhMTExCA8Pt693cv78eYSHhyM6OhphYWGljteaLEiX6arsxa7FXcH0seEQrVZs3f8jgu+9T+qSIABo5esFtZI3bWb/p9BkGwY8m00wy3MKI2M4mUH4KZ3w/5oCA1txuKpN3INVh7z//vv2cAXYhv4EQcDChQtLHEdEyCwyIlWjl3W4OnX8GMYPGQhftRrfHPxFFuHKTQDaqFUcrphdgYnwfYptqYWTmRyumHydyrSFq0ebcbhyBu7BqkOuXr1qD1fFiAhxcf+3LYLcV2Uvtufbr7Fgzqt4uG8olq/fDEAtCHYAACAASURBVB9ftdQlwU0Q0MbXCw08+W3DgDQt4UwW4UIe7xHI5O+PTMKPabZw9URLDlfOwH8p6pBOnTrh/PnzJUKWIAgIDrZti1Bktm3ULNeFQwFbIPzikw+xetmnGDXhRbzzyTJ4eEjfW+QuCGijVsHbQ/olIZh0LKLtsvYzWYQbOvm+jxi73ZkswsE027Aghyvn4YBVhyxYsKDEHKzi/86fPx/ZOhOydEZZbtRczGgw4F+zZuDA7mjMefd9THp1tiwaAndBQFs/FVQKDlf1WVIhYWO8JwRPHgNkruPP7H8mtDfhYUFn4zlYdUhYWBiio6PRvXt3eHl5oXv37ti1KxoPDhiMTJmHq5ysLLwUNgyHD3yPf2/Yismz5siiIVC4CQjkcMUA5JsAg7xH1hkrIfafda4ebAIM4qsFna7SHqz4+HhMnDgR2dnZaNiwITZt2oSuXbs6ozZWA2FhYfYrBvUWK9ILDSg0yXfLGwCIv3wJr4x/Fga9Hpv2/oDuvR6QuiQAgEIAAtUqKDlcMcZczIV8N/xwk/BAAC/FIJVKe7CmTZuGqVOn4urVq4iMjMSkSZOcURe7S8UbNZtkuuUNAPy8fx+evL8bRjzWG1kZt/BK1NuyCVcebgKae7lzuGKMuZwLuYQDNxXo2RgYwiu0S6bCgJWZmYnY2FiMHz8eABAeHo7ExEQkJSU5ozZWAyIR0gv1uFVkhIznsuOn7/ZidsR4pKckAwDMJhMWvj4bP+/fJ3FlgKebGwL9vOHhxo0SY8y1XM4nxCQRuvqJGNqWw5WUKgxYqampaNmyJRQK20iiIAho27YtUlJSyjxeq9VCo9HYv4xG3njLmYwWKxLzdSgwyntI0GQy4b25s0rcVzwpf9XSxRJVZaN0d0Ognwqe7jw9kTHmWuILCLsSCV0bAoNbWHjjZolVOgfrzvR75zpLtwsNDS1xOzIyElFRUTUsrfbl5eVJXUKNlFW31iIixyhCrgOCBr0BAJCfm4PI6ZNRkF/6ZyAiJF6Lh16vd3Z5AABPNwG+Xm4oLLB9MKhL/3/Inb+/v9QlMObSEjSE7QmEjmpgZDsBJoPUFbEKA1abNm2QlpYGi8UChUIBIkJqairatm1b5vHHjh1DSEiI/bZSqYRSqXRsxQ7mqg17cd0iETKKjNCKZihVEhdVidTEBMycMBZ6XRHaBLVDWnJSqTW72t3TCSqV838QL4Ub2qpVULiV7Lly9f8/GGN1X2IhYdt1QjtfYFQ7Ae7ccyULFY6DNG3aFD179sRXX30FAIiOjkZQUBCCgoLKPN7Hxwdqtdr+Jfdw5epsGzXrZL1Rc7GjP/6A54cMhI+vGtt/Oop5Cz6wDwsCsK/Z9fIbbzq9NpXCHYFq71LhijHG5O5aAeGba4S2PsCY9gIUPHdUNir9i7JmzRqsWbMGnTp1wuLFi7F+/Xpn1MUqoTGakZivg8Ei10FBG1EUsXLpYrw+JQKP9n8CW/f/iJZt2mLg0OFYsfErdOp6LzyVSnTqei8+2/Q1Bjw9zKn1eXu4o61aBXdulBhjLuZKPmFbAqG9Gniug8AX5shMpXOw/n97dx4fVXn2f/xzn9mTyWRPWEISAiTshEWgILIoVamiJLKINAIibRGr1kfQR4WKiFopdenPUgVRHi11IShSrGIVUIuiIoqCAQUSggsSIIGQDFnO748x0UD2zHaS6/165SUjkzNXhpNzvnOf+1x3Wloa27Zt80ctohF0XafgTCVllcF/gb341CnuvPF3bNrwCr+9dT43zLsD7WejRGMvG8/Yy8YHrL5Qi4lOLodMBBVCGM6uYzrrDur0iIAMuSwYlGSpHAMpq6gk/2QpRWWVOIL8Xy7vwH5+f+1UDufl8egz/2DY6AtrhKtAC7Oa6Rhml3AlhDCcT47qrM/T6RcF45OUHMeCVPCc8US9Tp4pZ/+J05SUB/9aHf/d/BaTfzkKd2kpa/79Hy4cd5lXt79pw3omjBxG/4RYJowc1uTeWS6bmQQJV0IIg9F1nfePeMLVwBi4QsJVUJOAFeT0H+8SzC8qoaKeFhnBQNd1Vj72ML+ZnEHfAYP45xtv07V7D6++xqYN67l5xjT27fmCM243+/Z8wc0zpjU6ZEXYLHR02qX5nhDCUMorPcHq9XydYXHwK+nQHvQkYAWxsopKDhaWUFByxi8LNbdkZOh0cTH/c/0Mli1awHU33sLj/3iR8IhIr9f4+EMPVN9xCE1rUBplt9DeaZODkhDCUE6W6TyzT2fXMbgySTE2QZPjmAEE+Uyetuuku5xvTpX6bdSqamSoKrxUjQw9vOrZBiei5+7/mpumX0N+bi7LVq7m4vFX+qzOg1/vO6fZbVWD0vrEOKzEhUrbECGEsRwu9jQQBZiRqugYKsHKKGQEK8hU6jrfnSrl0En/XhJs7sjQljf+zeSxozjjdvPP19/yabgCSO7S7ZxPblUNSusSF2KTcCWEMBRd1/mkQGfVXp1wK1zfXcKV0UjACiLuHxuHHgtA49CmjgxVVlby+EP3M+eaSQwaNpznN232+nyr2sy57fZGNyhVQHunjZgQq8/rEkIIbyku03lhv876XJ0+UXBtN0WYRcKV0UjAChInSgPbOLQpI0Mnjh/jhmsm8fhDD3Dj7Xfx6DP/IMwV7pc6G9ugVAEdwuxE2iVcCSGMY/dxncf36OQVw8TOiiuSNOnOblAyByvAKip1visupdBdHtA65tx2e405WHWNDO357FNunvlrigpP8Lc1LzHiwrF+r7WhBqWagoQwB06r7N5CCGM4Xa6z8ZDOF8ehR4TnLsFQGbUyNBnBCqDTZRUcOHE64OEKGjcytO4fz3LNr8YS5grnxTe3BiRcNcSkFImuEAlXQrQCe95ax5NTB3L/cCdPTh3InrfWBbokr6vUdXYc1Xl8t87XRZCRrJjYWcJVayABKwB0XeeH025yC09zpjJ41hIce9l4sjf/l0/yfyB783+rw5W7tJQ/3noTd900h8syJ/Hcxk0kJCUHtthamDVFUriDEIupzudkZ2fTr18/HA4H/fr1Izs7248VCiEaa89b63hh3kS6uA8zf3ASXdyHeWHexFYVsr4u0vn7Hp1X83Q6h8Gcnoo+UdLfqrWQj/l+5i6v5JtTpYboyA5wOC+XW2Zmse/L3dyz7DGu+vW1gS6pVlaTRqLLgdVU92eG7OxsMjMzqy9/7tq1i8zMTNauXUtGRoYfqxVCNOTdFYsZkRTDc1cMQCnF9f2TmPryDt5buZgeYyYEurwWOVKi88Zhz4hVYijMSpM7BFsjGcHyE13XOVZyhgOFxYYJV1s3vc5VF46g8MRxnvvXpqANV3azRnJ4/eEK4J577qm1FcWiRYv8UaYQogmO5uYwKjG6xh3Do5Oi+eFgToAra74fSnSyD1SyfI/OcTdMSlFMl95WrZYELD84U1FJXlEJ3xW7qWxka6uWrrfXEhUVFTyyZBG/mzqR/ucN4YU3t9CzX7rfXr8pQi0mklwhmBuxkPTevXtrbUWRk2PcA7YQrVVMUhqb8wpqfCB6O7eA2OS0AFfWdN+X6Ly0v5LH9+jknoJLOinm9FD0iJDLga2ZBCwfqhq12n/iNMVljR+1aul6ey1R8MMPzJ40gRWPLOPmOxfy12ef98mSN94QbjOT6HJgauQtzKmpqbW2okhLM94BW8C+ffsYNmwYqampDB48mN27d9f6vMWLF9OlSxe6dOnC3Xff7ecqRXOdP+su3sk9ytSXd/DEjoNMfXkH7+YdZfgs4/wbflOs8/zXnhGrw6fhskTFjb0Ug2NVo49bwrgkYPmIu6KS3OpRq6Z1ZG/JenstseOD95l44Qj27fmCJ196hetvvhWtESNDgRhti3JY6NDERZsXLlxYa5PShQsX+qpM4UO/+c1vmD17Nnv37mXevHlcd9115zxn69atrFmzhs8++4zdu3fz2muv8frrrwegWtFUPcZMYNKfXuSAoyMPbs/lgKMjkx56iR6jfbtaREvpus5XRTqr91XyZI7OkVIYn6SY20sxMEZJT6s2RAKWl+m6ztHTbg6cKOZ0E0atfq656+01l67rPLP8r8y4chwdE5N46a13GTpiZKO+19+jbQqID7XRLrRp4QogIyODtWvX0rdvX+x2O3379iU7O5sJE4w9YbYtOnLkCDt27GDatGkAZGZmcuDAAQ4ePFjjec8//zzTp08nNDQUm83GzJkzWbNmTQAqFs3RY8wEZj33MXe8d4pZz30clOHq560kHps8gEf+L5vnvtJxV3gahd7QU9E/WmGSS4FtjgQsLzpdVs7+E6c5cvpMo+da1aY56+0118miQm6e8Wv+dPf/Mm3273hq3Qbi2rVv9Pf7c7Stqjt7tKP53dkzMjLYuXMnJSUl7Ny5U8KVQR06dIgOHTpgNntuhFZKkZiYSF5eXo3n5eXlkZSUVP04OTn5nOec7dSpUxQVFVV/ud1u7/8AolWobiVR6mkl0bPsGwofncyob19mVpqiZ6RCk2DVZkmbBi8oq6zkSLHbaw1DG9tVvaX27PqMP1yXxfGCAh595h9cOO6yJm/DX6NtJqVIcNkJtbSdXTY7O5t77rmHvXv3kpqaysKFC6WdxM+c/SHk7P2wtufV9ZyfGzmy5ujtvHnzmD9/PoWFGiUlgW8K3BglJaWBLqHRjFprgVux8fF7OT8xhueurNlKIue5xQwec2kAKzXu+xrMioqKOHas7t6VUVFRNR63nbOVD1TqOsdKyigoOUNFE+dZ1aeqq/rflj7Aga/20blrN+bcdsc56+01l67rZD/3fyy+/Va6pHbn789nk9g5pVnbSu7SjX17vqhx4vL2aJtFUyS6HNjMdTcQbW2kZ1f9OnXqRH5+PuXl5ZjNZnRd59ChQyQmJtZ4XmJiYo3Lhrm5uec852xbtmwhPf2nu2ZtNhs2m41wXcfhKMbhcHj1Z/EVo9QJxqm1UodDZ+x89IPOviLQvtnL6CFJ57SS2L59b1D8TMFQQ2MZoVaXy0VUVESjny+XCJtB13UK3WV8fbyYI6fdXg1XVerqqt5SJadPc+eNv2PBLXO5YvJUntu4qdnhCjyjbbVNHPfWaJvdrJEcEdKmwhVIz66GxMXF0b9/f5599lkA1q5dS3JyMsnJyTWeN3HiRJ555hmKi4txu9089dRTTJkypd5tO51OXC5X9ZfNZvPVjyEM4uQZna3f6jzxlZU1X+ucLIMrkhTtklNbTSsJ4X0ygtVEJ8+U80Oxm9KK4FniprEOfLWPW2b+mvzcgzzw+BNcPrH+E01j+HK0LcxqpmOYvU3OYZCeXQ37+9//zvTp01myZAkul4tnnnkGgHHjxrFo0SIGDRrEqFGjmDRpEn369AFgypQpXHLJJYEsWxhEeaXOV0Xw2TGdL0+AWYPuYZUMaWemQ4jnw6Rt1t28MG8iU1/eweikaN7OLeDdvKNMemh5oMsXQUACViPous6pMxUcLTljmC7sZ3tt3VoW3HIj7Tp04J+vv03X7j28tu2xl41n7GXjvbY9gEi7hXahtjbbhC81NZVdu3adc+lVenb9JC0tjW3btp3z/zdu3Fjj8YIFC1iwYIG/yhIGVqnrHDgJnx/T2VMI7gqId3gag/aNAv1MOQ6Hpfr5Va0k3lu5mO3bc4hNTmPSQ8uD8m5H4X8SsOpRqesUucspKDmD24AjVgBn3G4eWngn/1j5BOMyruKPf36UUKcz0GXVqaoNQ1QL7hRsDRYuXFhjDpb07BLCN85U6Ow/CXsLdXIK4XQ5RNlgSCz0jlTEOn76kFdSy/f3GDOhyWsj7nlrHe+uWMzR3BxiktI4f9Zdhl9fUZzLqwHrdFl5jfk4RlVWUckJdxnHS8sob0m/hQA7nJfLH2ZdS84Xn3PXg39myoxZQf1voylFxzA7YVbJ/VU9uxYtWkROTg5paWksXLhQ2koI0QgNBZgTbp29RbCv0DNiVaFDtA36RUHvKEV7x7l3qXqzthfmTWREUgzXDU5ic95hXpg3kUl/elFCVivj1TPZd6fc7D9xmnCbhXC7GUsjuoAHi0pd59SZck6UllFcVoFxY5XHljf+ze03zCbMFc6z/3qD3ukDAl1SvSyaopPLgb2NTWavT0ZGhtwxKEQT1RVgzv/jC7j7XcnXRXDM7bnDKykMLuqo6OaCaLt/Pny+u2IxI5JieO6Kmq0d3lu5WAJWkAkxQ6JTkeiEJKfCVtq0K1leHypwV1Ry5LSbH067CbWYcNkshFnNQbnuUqWuc/JMOUXuck6dKffJ3YD+Vl5ezmP3L2bFo8sYdfGl3PfY34iIjGr4GwPIYTaR4LIbKpALIYJTbQHm6pd3sO2p+3AuupIuLriog6KzC+wm/5+XjubmcN3g2lo7yA0sgaYp6BSq6B4B3cIVMWeF7mNN7Dnss2sxOnCqrIJTZRUoIMRiwmk147SYAnrL/ZmKSorLKjh1ppyjpyuw6bVdVTemH777jtt+M5MdH2zj1oX3Mn3OjY1aSzCQXDYzHZxt805BIYR3FZfpHKklwIz5McDc1EsFfJpETFIam/MOc33/pOq5ldLaIXDMCrq4FN0jFGkREGL23v7hl8kuOlBcVkFxWQXf47kc5LCYCDGbsJtN2M2aT06wFZU67opKSssrKCmv4HRZBWU/m1Plj2nrmzas5/GHHuDg1/tI7tKNObfd7vU77gA+eHcrt82eiaZprFr3Lwb+YpjXX8PbYkOsxDisAT/gCSGMy12hs+eE586//SdBi0tlc+435wSYuOS0oDjWnD/rLmntEAQ6hir6R3tuZLB7MVT9XEBmE5dV6pS5PZfmwHPnmNWkYTVpWEwKq6Zh1hRmTcOkPJOftepPI55tVOo6lbrnvxWVOuW6TnmlzpmKSsoqddzlFZRX6gGdS1W1EHLVL3nVQsgPr3rWayGrsrKSFY8s47EHFnPe8BH8aflKYuLivLJtX9EUtHfaCbdZGn6yEEKcpaLS00l91zGdvYVQrkOSE8Z1Upjn3MUrd0wK2gAjrR0CJ9QM/aI9i2///O5QXwmK27V0PHO3jNoKoS71LYTsjYB14lgBt8+ZzbtvvcnsW/6HG+b9LyZTcE8St2iKBJcDh0xmF0I00TG3zsdHdXYWeNopxDtgVAdF70gIt/54whybgc0U3AGmOa0dQNo7NFe8A34Rr9EnEr/OBw+KgNVa+XIh5E8/2s4fZk2ntOQ0f1vzEiMuHNvibfqaTGYXQjRVha6TcwI+Puq5BGg3Qd8oGBijiKtjFKK5ASaYSXuHplFA13DFL+IUKa7AXBqWgOVDvlgIWdd1/rHi7zy6ZBG9+w9g6ZNP075jgjfK9SmnWZEU7pDJ7EKIRimv9IxUvfOdTlEZdAqFK5MUPSM9I+FtjbR3aLyuLsUvE+oO4P4iAcuH5tx2e405WC1dCPlkUSF3/f4G3vzXeqbPuZGb7/ojFktwz2NSQFyoDVVSJuFKiDassZe3zg5WfSJheLwiPqRtHz+kvUPDomxwcYJGWkRw7CtyrcaHqhZCTu3ZC6vNRmrPXjzy9HPNWgh596c7mXjhCD54ZwtLn1zFbffcF/ThyqQ8zUOj2/iyN0K0dVWXt7q4DzN/cBJd3J7LW3veWlf9nEpd56MfdB77Qudfh3SSnDCnhyKjs+aXcLXnrXU8OXUg9w938uTUgTVqCwae9g4FNeb0SnsHD6vmaRh7Q8/gCVcgI1g+19KFkHVd55+rVvDg3XeQ2rMXK15aT3RcvBcr9A27SSPB5cBqkgwvRFtX3+Wt5F9cQu5Jndfydb4vgd6RcEE7/9zlVcUI85ukvUPturoUVyQpwqzBE6yqGPLst2nDeiaMHEb/hFgmjBzGpg3rA12ST5wsKuTWWdNZPP9WJv56Os9ueIOEpORAl9Ugl81MckSIhCshBOC5vDUqMfqcy1s/HMzh1cNmnt6nY1YwK02R2Vnza7iCmgFw9oBknrtiAOcnxvDeysV+raM+Ve0dDjg68uD2XA44OjLpoZeC6u5If44CmhVckqC4pmtwhisw4AiWP3pLBYPdn+7kD7Ou5XhBActWrubi8cHzS1QXBcSG2IgJkUuCQoif1NW9vDI+jdxijSuSFP2ifLfAckOMMr8pmO+O9OcoYJwDMpP9c+m4JQw3xFBfb6nWQNd1nlvxd6aOu4gwVzgv/merIcJV1XwrCVdCiLOdP+su3sk9ytSXd/DEjoNc/fIO3s07SsrUO7m+yxnSowO7hI3Mb2o5f40CDo5VzO4e/OEKDBiwfNlbKtAKTxzn5hm/ZskdtzH52pk8t3ETiZ1TAl1Wg+xmjc4RITithhsQFUL4QY8xE7jy/hfYaerIfR/k8rHWkYsXv8jVV2VgC4Kew2cHwKk/BsDhs+4OdGmGUd9lYG8wK8jsrDEu0bPSixEY7ozoi95SwWDnhx9w22+u41RRUbPvNAyEcJuF9k6btGAQQtRK13U+Pw6b4q+k/M4ruaSjYmAMQXXMkOVrWs6Xi1iHmGFKF41EZ/DsM41Ra8AqLS1lypQp7N69m5CQENq1a8fy5ctJTk72c3nn8nZvqUCrrKxk5aN/4bEHFtNnwCCeeWUjHTolBrqsBikgPtRGlLRgEELUoeiMzoY8z7qBPSM8k5KDdUJyMM9vMgJf3eUYY4drumpE2oJzv6lPnZcIZ8+eTU5ODjt37uSyyy5j9uzZ/qyrTt7sLRVoR777lllXXcEjSxYxc+7NPG2QcGXRFEnhIRKuhBC10nWdT47qPL5b59sSmJKimJiiBW24Ei3ni7scO4cpZqUZM1xBHSNYdrudcePGVT8eOnQoDz/8sN+KakhLe0sFg7dff427fv87LBYrK9auZ+iIkYEuqVFCLSY6htkxy3qCQohaFJ7ReTVP5+siSI+CXyYoHGZjniBF03hzFDA9WnF5ovLr4sze1qg5WI8++iiXX97wCNGZsjO43e6fNm4yYTIbbpqXT5WcPs1Df7yT51etZNTFl7L4kceJjI4OdFmNEu2wEhdiDejdPkKI4KTrOp8UwOv5OjYTTO2i6BYuxwrRdEPjFBcnBPbOUm9oMP0sWbKEffv2sXx5w9dRn171NJE/6zI+fPhwhg8f3rIKfai0pNSvr/flrs+466Y5fJufz+33/YnMaVkopSgpKWnSdvxdt6YgxqphcZdx3F3crG0cP37cy1X5nhFrBmPWHRUVFegSRCPVtqZgpxFX8mquzt4iSI+GixMUdpOxT44iMEa1V4zq0DqukFQHrNWrV7Ns2TIAbrrpJmbMmMHSpUvJzs7mzTffJCQkpMGNTZ8xnR59+v20cQOMYDkcDp+/RkVFBSsf/Qv/709L6Nq9Jy/+ZytdUlt2Z4U/6gZPC4aEMO8seWPEk6gRawbj1i2CW13NJK03PI9l4ASmpKigWgtOGIfCE8yHxreOcAU/C1hZWVlkZWVV/8WyZctYs2YNb775JhEREY3amNVixWazeb9KA8vd/zX/O/e3fPbxh1x34y3MmXcHVqsxJodH2Cy0kxYMQogf1bam4NUv72DnxiX8bloGoRY5Voim0xSMT9LoH9O69p9ah5fy8/O59dZbSUlJYfTo0QDYbDY++OADvxZnZJWVlfxz1QqWLVpAdGwcT7/yGgOH/iLQZTWKpqBdqJ0IuyXQpQghgkhtS8qMSYrmw+05Eq5Es5gUXN6xvNWFK6gjYCUkJJzTLb012LRhPY8/9AAHv95HcpduzLrpD/wq4yqvv87hvFwW3HIj72/dzOQZ13HrgnsJdTq9/jq+YDVpJITZsZuDoL2yECKoxCSlsjn3G580kxRtj1nBpBSNmMrKQJfiE8E9QcqLalsket5vZmK1Wr3W8qGyspI1Tz3JX+79I+GRkTz54ssMGzXGK9v2B5fVTHun3dC3xQohfKOsUscy/k7eWTrZ680kRdtj0Tzd2bu4FMeOBboa32g9s8ka4OtFor/em8O14y9lyR23MX7yFF55533DhKuqruwJLoeEKyECbM9b63hy6kDuH+7kyakD2fPWukCXREGpzoovdb7pPoEhC17wajNJ0fbYTDCtqydctWZtZgTLV4tEn3G7efKRP/PEw3+mY6dEVr38LwYPH9GibfqTRVN0DHMQYpFLgkIEWl136U3604sBW8bli+M6r+bqOC0wK00Rl54B4zMCUoswPvuP4SrBYOsKNkebCVi+WCT6vbf/w323/w+HD+Vx3Y0385tbbsNmt3ujXL+QruxCBJfa7tKb+vIO3lu52O8Bq7xS543DOh/+AL0i4PIkhU16W4kWCDXDtG4a7UPaxn7UZs6sc267vfqyINCiRaK/yT/EH667ltmTJhDXvgNr336P399xt2HClQJiQ6wkuhwSroQIIkdzcxiVGF3jODU6KZofDub4tY5jpTorc3R2HIVfdVJkdpZwJVom3Aoz0tpOuII2FLBqWyT6oSdWNWmR6NPFxfz1wfu4fNggPn7/vzy4fAWr1m2ga1p3H1buXWZN0cnlIDbEZvhlCIRobWKS0ticV1Bjrqi/79L74rjO37/UOVPpuSQ4KNb4S5aIwIqxw8w0jRh729qP2swlQjh3kejGLlFTUVHB+ufX8NgDizlWcJRrfzuX2bfcSqgzzFel+kSIxURHpx2LF7qyCyG87/xZd/HCvIkBuUvvVJnOpsM6nx2DXpFweaKMWomW6xCiuKarapN90tpUwGoqXdfZsunfPLz4Hvbt2c3F4ydwy91/pFNy50CX1iQKiJKFmoUIej3GTGDSn17kvZWL2b49h9jkNCY9tNynd+lV6jrbf4DN3+g/dtRWpEchxwrRYslhiqu7tN2gLgGrFrqu897b/+GvD97Hrh0fc97wEfzzjbfp039goEtrMpNSdAizE2aVf2ohjKDHmAl+m9Cee1Jn4yGdI6UwKAbGdFA4zG3zZCi8q2+UYnySwtyGW//IWfdnqkasnvjLUj796EPSzxvMyrXrGTJipCE/zTnMnrsEvbFQsxCiOLnFNQAAH01JREFUddB1nbxi+OCIzp4T0DEEru+u6NCGJh8L3zEpz6LNg+PkvCMBCygrK+P19et46rFHyPliF/0HD2X5P9dy/piLDBmsAKLsFuJDZSK7EMKjvFLn8+OeYPVdCUTb5HKg8C6XxbP0TVvocdUYbTpgnSws5PlVK/i/J/7Gd4fzGTZqDE+/vJFBw4Yb9oCjKUV7p41wmyzULISA426dTwo8LReKy6GrC67poOjikmAlvKdzmOKqzm1zMntdAhawzl54ec5tt3ttTcCG7Nuzm3+sfIL1L6yhvLycX2VM5NrfzSWtV2+/vL6v2EwaCWEObGYZmhWiLdnz1jreXbGYo7k5xCSlMfDXt2EZMplPCnQOnASrBv2iYXCsanO3ygvf0hScH68Y1UGhSWCvISABq7aFl2+eMY2HVz3rs5BVWlLCG6++zAvPrOKT7e8TExfPtb+7katnXk9sfLxPXtOfImwW2jltsoML0cacs7xO7mFevXsa/NZC4gUTuCJJ0TMCrG30Ti7hO9E2mJAslwTrEpCAVd/Cy94MWLqu8/nOHaxb8ywb177EyaJChl4wimUrnmH0pb+ioqICh8PhtdcLBE1BfKidSLtcEhSiLXpnxb2MSDx3eZ19by1hxqzMQJcnWiEFnBerGJugsLThuwQbEpCA5auFl6t8cyiPDS+9wKsvPc/+vTnEt+/A1TOv58opU0nq0rX6eY1tNBqsrCaNhDA7drMs1CxEW6LrOvnF8PFRne8P7mXWkKRzltfZvt2/y+uItsFlgSuTNVJcEqwaEpCA5YuFl3/47jteX7+O117OZueHH+AICeHCcZcx/977+cXI0ZhMrSuEuKxm2jvtmOTTg2hjSktLmTJlCrt37yYkJIR27dqxfPlykpOTz3nu5s2bGTduHKmpPx1btm3bZtiR67JKnV3H4MMfPHcCRlohLCGVzXnfcH3/pOorA/5eXke0DT0jFJcnSa+0xgpIwJpz2+015mA1d+Hlw3m5vLnxVTa9up6dH36AyWzm/DEX8eDfnmT0Jb8i1On00U8QOAqItmokuIx5ghDCG2bPns2ll16KUoq//vWvzJ49mzfeeKPW5/bs2ZOPPvrIzxV617FSnQ+P6uwsgNIKSHXBhT/eCfjlb+8O2PI6om2waHBpJ40BMRKsmiIgAatq4eW/LX2AA1/to3PXbsy57Y4GF16urKzki08/YfPrr/HWvzey94vPsVitDBs1hsWPPs7oS8YRHhHpp5/C/6yaRscwOyUnywJdihABY7fbGTduXPXjoUOH8vDDDwewIt/QdZ2Dpzx9q3IKwWGCATEwKEYRafvpRFfb8jpX3veIT5fXEW1H+xDI7Nz2Fmr2hoC1aTh74eW6nDhWwLYtm3nv7f+w9c03KPjhCK7wCEZcNJbf3Pw/nH/hRTjDXH6oOLDCrGY6/HhJ0Ngzx4TwrkcffZTLL6/7w1lOTg4DBgzAZDIxY8YM5syZU+/2Tp06RVFRUfVjm82GzWbzWr0NqWoI+v4Rne9LIM7uWXi5TxR1Tig+e3kdo88vFcHhF3GKizoqmYrSTEHXaNRdWsrOj7bz/pbNbNv6Np9/sgNd1+navQfjJ1/NqLGXkD54CGZz0JXuEwqIC7UR7bAGuhQhgs6SJUvYt28fy5fXfjlswIAB5OfnEx4eTn5+PuPGjSMmJoZJkybVuc2RI0fWeDxv3jzmz59PYaFGSUm5V+v/uZIK2HncxI5jJoorFCnOCiYlVpAUoqMUlLuhsa9eUlLqszq9TWr1jZbUqikY266c9NBKCk94sag6HD9+3Pcv4gUN1RkVFVXjccBTyuniYj7b8REfb3uPD//7Hp9+tJ0zbjdRMTEMGTGSSVkzGTZ6DO06dAx0qX5n0RQdwxyEWFrXBH0hmmr16tUsW7YMgJtuuokZM2awdOlSsrOzefPNNwkJCan1+1yun0a3ExISuPrqq3nnnXfqDVhbtmwhPT29+nHVCFa4ruNwFHt9gvxxt877R3Q+KQBd9zQEHRqniLFbgOa3XzHSRH6p1TeaU6tFg6s6a6RF+HfU6uxwEqyaUqdfA5au6xw6eIBdn3zMpx9uZ+dH2/ly12dUVFQQHhnJoF8M5w9338N5w0eQ2rMXmtZ2O5KH/XiXYFteiVyIKllZWWRlZVU/XrZsGWvWrOHNN98kIiKizu/79ttviY+PR9M0Tp48yYYNG7juuuvqfS2n01kjmDXX2d3Vz591V43LeN+c1vnvdzq7T4DdBMPiPL2FZKkRESghZpjaRRqHeovPAlZFRQWHDu7ny1272LPrM3Z/tpMvPv2Ewh+H2JJSutBv0HlcNe1aBgwdRkq31DYdqKpUXRKMsltknTAhapGfn8+tt95KSkoKo0ePBjyjTB988AEAs2bNYvz48YwfP561a9fyt7/9DbPZTHl5ORMnTmTGjBk+r/Gc7up5h3lh3kQmPvgCtvMm8N73niVsIm1waSdFenTd86uE8IdIG0zrqhEtk9m9xqsB6+Xn/8FzK/7Ovi9383XOl5T+ONEyvn0Huvfpy7TZv6NP+gB69x9IZHS0N1+6VbBqGh3C7HJJUIh6JCQknNOo+OdWrFhR/ee5c+cyd+5cf5RVw7srFjMi6dzu6q8+fh+ld15J+xC4qrOiRwSyvJUIuDgHZHXTcMroqVd5NWCtfXY1Xbt3J7VHL36VOZFu3XvSvXdfomJivPkyrdLP7xIUQhjb0dwcrht8bnf1bR/k8Ouuis5hyAi1CArxP4YruTTtfV4NWKtffY1e/fp7c5OtnqYgLsRGlNwlKESrEZ2Uyubcc7urt+ucJkuMiKAR74BrUzVCpDO7T3g1YCnkH6kpbCZP41BZS1CI1qGkXGf7D1Bw0Z18+/hk6a4ugpaEK98LeJuGtirCZqGd0ybzL4RoBYrLdLYd0fnwB6jUof/YCcR1fIGdq++r7q4+6aHl0l1dBIV2DsiScOVzErD8zKQU7Zw2wm3N728jhAgORWd0/vu9zsdHPZf7B8XAL+KVZ7JwpwwGXZwR6BKFqEHClf9IwPKjELOJDmF2rCZpRyGEUem6Tl4xfHhEZ88JsJpgeDwMiVM45KQlgli7Hy8Lyn7qHxKw/EABMSFWYhxWuXNICIMqKdf57ITGJyc8awRG2+CXCZ4eVjaT/F6L4NY+xHO3oIQr/5GA5WNWk0YHp/S2EsLIKnWdSf+pIL/YTKoLxnZUpEirBWEQEq4CQwKWD0XYLbQLlYnsQhidphTz+ml8/n0J7cONsxadEBKuAkcClg+YNUX7UDthNnl7hWgtRrTTyD8e6CqEaLwOIYpfd5O5gYEiCcDLXFYz7Zw2zLKuohBCiACJtekSrgJMApaXSPsFIYQQwSDGDuNjyyRcBZgELC8Is5ppL6NWQgghAizC6plzVX4q0JUICVgtYNYU8aEyaiWEECLwXBZPnyuXVXEs0MUICVjNFW4zEx9qx6zJEKwQQojACjV7OrRH2uScFCwkYDWR1aTRLtSG0ypvnRBCiMALNXsuC8bYJVwFE0kJjaSAaIeVmBCr9LUSQggRFFwWz8iVhKvg06hZ2ffccw9KKT7//HNf1xOUQi0mUiJCiZOmoUIIIYJEpA1mpEm4ClYNBqwdO3bw/vvvk5iY2ODGbr1+Bps2rPdKYcHAoik6htlJCg/BZpY7BIUQQgSHGDvMkDlXQa3e1OB2u7nhhht4/PHHG7XmVt6Br7l5xjTDhyxNQYzDSpfIULlDUAghRFBp5/CEK5dVwlUwqzdgLViwgGnTptG5c+dGb1ApxeMP3Y/b7aaivLzFBfqby2aWy4FCCCGCUucwxfRUjVCLnJ+CXZ2T3Ldt28aHH37IAw880KQN6rrOVzlf8sADDzB8+HCGDx/e4iJ9pbSktPrPNk0RZdWwl5VxqrAkgFU17Phx4y2IJjX7jxHrjoqKCnQJQgS9oXGKXyYo+fBvEDUC1urVq1m2bBkAkydP5ssvv6wevcrPz+fiiy9mxYoVXHrppXVuUClF17Tu3H777ZhNJkzm4L5R0RUaQlyoDZfV3KjLoMHCiCckqdl/jFp3sMrOzuauhYv4el8OMUlpnD/rLnqMmRDoskQbYVZwWZJGerRxzlHirICVlZVFVlZW9eM77rij+s/Jycls2LCB3r1717tBXde5Yd7/YrPZvFyqd5mUZ8QqKTJUPg0IIeqUnZ1NZmYmFyTFcNXgJDbnHeaFeROZ9KcXJWQJnwuzwOQUjQSnnKeMxqu3xiWldOGRp5/jol9d7s3NepUCouwWzwR2iybhSghRr/vuXcQFSbE8e8UAZg9I5rkrBnB+YgzvrVwc6NJEK5fkVMzuLuHKqBp9/e7gwYMNPmfpE6vo2S+9JfX4lMtmJi7EhtUkLReEEI3zZU4Otw5Kqp5CoJRidFI027fnBLgy0VpZNLiwg2JInDLU1BVRU5tIGiEWE53DQ0gIc0i4EkI0Sfe0NLbkHUPXdcAzDeLt3AJik9MCXJlojRKdit/10Bgar0m4MrjgnoHeQjaTRlyojTBZN1AI0Ux33r2AzMxMrnllB6MSo3k7t4B3844y6aHlgS5NtCIyatX6tMrhHLOmaO+0kRIRIuFKCNEiGRkZrF27lgJXAg9uz+WAoyOTHnqJHqOvDHRpopXoGaG4oaeMWrU2rSp9aAqi7LIgsxDCuzIyMkgZOYE1e4pxOByBLke0Eu1D4JIEjaQwOV+1Rq0iYCnAZbMQF2LFInOshBBCBLEwC4zpoJEejYxYtWKGD1ghFhPxoTYcZlOgSxFCCCHqFGaBYfGKgTEKq0mCVWtn2IBl1TTiQq24ZDFmIYQQQSzKBsPjPSNWJk2CVVthuIClKYh2WIl2yDwrIYQQwSve4QlWvaOQ81UbZKiA5bKaiQuVRqFCCCGCV+cwxfB4RddwCVVtmSECls2kER9qwyktF4QQQgQhTUGPCE+w6hAqwUoEecCSy4FCCCGCmVWD/jGKoXGKSJucp8RPgjZghVpMtHfa5XKgEEKIoOO0wJBYxaBYhcMswUqcK+gClllTxIfYCLfL3YFCCCGCS0KoJ1T1jvScr4SoS1AFrHCbmfhQu+y0QgghgobdBP2iFSmmMtI6yFUV0ThBEbAsmqK90y6T2IUQQgQFk4IUl6J3pKJnpOc8deyYHuiyhIEEPNFE2i3Ehdik+ZoQQoiA0hQkOxW9fgxVMrdKtETAApZFU3Rw2gmVUSshRBMlJydjt9ux2+0A3HHHHUyePLnW5y5evJhVq1YBMHXqVO69916/1SmCn0WDLi5FargiNRycFglVwjsCkm4ibBbiQ2XUSgjRfC+99BK9e/eu9zlbt25lzZo1fPbZZ5jNZoYPH87555/PxRdf7KcqRTAKt/JjoFJ0DpPJ6sI3/BqwzJqifaidMJuMWgkhfO/5559n+vTphIaGAjBz5kzWrFkjAauNMStIClN0dUFXlyLWIYFK+J7fkk6Y1Ux7pw2zJndgCCFa7pprrqGyspIhQ4Zw//33Exsbe85z8vLyGDlyZPXj5ORkXnrppXq3e+rUKYqKiqof22w2bDab9woXfuEwQVqEokeEZ5TKapJQJfzL5wFLUxAXYiPKYfX1Swkh2oitW7eSmJhIWVkZd911F9deey0bN26s9bnqZ6tA6HrDd4H9PJABzJs3j/nz51NYqFFSUt6ywv2kpKQ00CU0mjdrdZigW1gFaa5KkkJ1TAqohFOF3tn+8ePHvbMhP5Bava+hOqOiomo89mnAsps1Ojrt2MwmX76MEKKVW716NcuWLQPgpptuYsaMGQBYLBZuvvlmUlNTa/2+xMREDh48WP04NzeXxMTEel9ry5YtpKenVz+uGsEK13UcjmIcDkcLfxr/MEqd0LJazcozUtUvStE1HJ8vq3b2STSYSa3e15Q6fRawohwW4kNsNT49CiFEc2RlZZGVlQVAcXExJ06cICIiAoA1a9bQv3//Wr9v4sSJzJ07lzlz5mA2m3nqqadYvHhxva/ldDpxuVze/QGE13UKVaRHK3pFgl3aKYgg5PWAZVaK9mF2wqT9ghDCB77//nsyMzOpqKhA13VSUlJYvXp19d+PGzeORYsWMWjQIEaNGsWkSZPo06cPAFOmTOGSSy4JVOmihRwmSI9WDIxVxNglVIng5tUU5DBrdI4MwSIT2YUQPpKSksInn3xS59+fPRdrwYIFLFiwwNdlCR9KcioGxniaf0pLBWEUXg1Y7Zx2CVdCCCFaTAE9IhQj2ivah0ioEsbj1YAl862EEEK0hKagd6RiRDvpVyWMTSZKCSGECDgF9I1SjGyviJL5VaIVkIAlhBAioNqHwLhOGp2cEqxE6yEBSwghREA4TDCifTljUjSZYiJaHQlYQggh/EoBA2IUF3ZUlBZVSrgSrZIELCGEEH4TY4fLEzWSwjyhyjiL+gjRNBKwhBBC+JymYHi8ZxK79LISbYEELCGEED7VIUQxPknRTvpZiTZEApYQQgifMCsY3UHxi3jl80WYhQg2ErCEEEJ4XZLTM2oVLT2tRBslAUsIIYTXWDW4qKPivFgldweKNk0ClhBCCK/oHKa4IkkRYZNgJYQELCGEEC1iVnBhR8XQOBm1EqKKBCwhhBDN1s4BGZ014mRhZiFqkIAlhBCiyRQwLF4xpoPCJH2thDiHBCwhhBBNEmqGzM4aKS4JVkLURavvL91uN3PnzqVbt2706tWLadOm+asuIYQQQSjJqfhtDwlXQjSk3hGs22+/HU3T2Lt3L0opvv32W3/VJYQQIohUXRK8sKM0DRWiMeoMWMXFxaxatYr8/Pzqu0Lat2/vt8KEEEIEB7sJJiRrpEVIsBKiseq8RPj1118THR3N4sWLGTRoECNGjOA///lPvRs7deoURUVF1V9ut9vrBQshhPCfDiGK3/SQcCVEU9U5glVWVsb+/fvp2bMnDzzwAJ9++ikXXXQRu3fvJjY2ttbvGTlyZI3H8+bNY/78+d6t2IuOHz8e6BKaxYh1S83+Y8S6o6KiAl2CqMWQOMUvO8pdgkI0R42AtXr1apYtWwbANddcg6ZpXHPNNQD069ePzp0788UXXzBq1KhaN7ZlyxbS09OrH9tsNmw2m49K9w6jHtiNWLfU7D9GrVsEB5sJrkjS6BkpwUqI5qoRsLKyssjKyqp+vGnTJl5//XXGjRtHbm4uBw4cIC0trc6NOZ1OXC6X76oVQgjhU+1DYGJnjShZpFmIFqn3LsLly5czc+ZM5s+fj8lk4oknnpCJ7kII0QpZNLignWJYvFwSFMIb6g1YKSkpbN682U+lCCGECISeEYqLOynCrRKshPAW6eQuhBBtVIwdLu2k0UWahgrhdRKwhBCijYlzwNA4jX5RyOVAIXxEApYQQrQBCugWrhgap2SZGyH8QAKWEEK0YiFm6BOlOC9WESN3BgrhNxKwhBCildEUpIUr0qMVqeFyGVCIQJCAJYQQrUSsHQbEKBLUGTrFOQNdjhBtmgQsIYQwMIsGvSMVA2IUnZyekapjxwJclBBCApYQQhhROwcMilX0iVLYTHIJUIhgIwFLCCEMwqygZ6RnwnrVaJUQIjhpgS5ACCFE/SKsMLaj4g99NTI6a14LV9nZ2Qzsn05oiIOB/dPJzs72ynaFEBKwhBAiaCU5FZNTNH7fW2N4O40Qs/dGrbKzs8nMzMR5/BtuHZSE8/g3ZGZmSsgSwkvkEqEQQgQRk/JMWh8ar2gf4rvLgPfdu4gLkmJ59or+KKW4vn8S17yygyWL7yUjI8NnrytEWyEBSwghgoDNBINiPJ3Ww/yw6PKXOTncOigJpTyvpZRiVGI0f/7oS5+/thBtgQQsIYQIIKcFhsR6Jq7bvXgJsCHd09LYkvcN1/f3hCxd19mcV0CP7t39VoMQrZkELCGECIAoGwyL93RbNweg0/qddy8gMzOTa17ZwajEaDbnFfBO7lGy//KE32sRojWSSe5CCOFHHUIUk1I0buylMShWC0i4AsjIyGDt2rWcjkrgzx/lcjoqgezsbCZMmBCQeoRobWQESwhhKCdOnGDUqFHVj0+fPs3+/fs5cuQIUVFRNZ67efNmxo0bR2pqavX/27ZtGw6Hw1/lVuvqUgyPV3R2BU//qoyMDJnQLoSPSMASQhhKREQEO3furH68dOlStmzZck64qtKzZ08++ugjf5VXgwK6RyhGtFN0CA2eYCWE8D0JWEIIQ1u1ahX33XdfoMuoQVPQK9ITrOIcEqyEaIskYAkhDGvbtm0UFBRw2WWX1fmcnJwcBgwYgMlkYsaMGcyZM6febZ46dYqioqLqxzabDZvN1qh6TAr6RSvOj1dE2SVYCdGWScASQhjWU089RVZWFmZz7YeyAQMGkJ+fT3h4OPn5+YwbN46YmBgmTZpU5zZHjhxZ4/G8efOYP38+hYUaJSXltX6PWYP0iArOi67AZQFOw7HTzf6xWuz48eOBe/Emklp9Q2r1vobqPHuaggSsRsjOzuaee+5h7969pKamsnDhQpkYKoQfrV69mmXLlgFw0003MWPGDIqLi3n++efZvn17nd/ncrmq/5yQkMDVV1/NO++8U2/A2rJlC+np6dWPq0awwnUdh6O4xgR5uwkGx3q6rntzGRtvqGtOWjCSWn1DavW+ptQpAasBVet1VTXi27VrF5mZmaxdu1ZClhB+kpWVRVZWVo3/9+KLL9K3b1+619MY89tvvyU+Ph5N0zh58iQbNmzguuuuq/e1nE5njWBWm1g7DIlT9I1SWE3BFayEEMFB+mA14J577qkOVwC6rqOUYtGiRQGuTIi2beXKlbWGpVmzZrF+/XoA1q5dS58+fejXrx9Dhw5l7NixzJgxo1mvp4CuYZVkddO4oZeJQbGahCshRJ1kBKsBe/furQ5XVXRdJycnJ0AVCSEA3nnnnVr//4oVK6r/PHfuXObOneuV1+sbBZ1UOVFB1MdKCBG8ZASrAampqdWLoVZRSpGWlhagioQQgXD2cUAIIeojAasBCxcurL4sCFRfLly4cGGAKxNCCCFEsJKA1YCq9br69u2L3W6nb9++sl6XEEIIIeolc7AaQdbrEkIIIURTeGUE68yZMzX+awRut5sHH3wQt9sd6FKaxIh1S83+Y8S63W43f/zjH4O+ZqO8t0apE6RWX5Fava85dSr97FvkmmHr1q2MHDmSLVu2cMEFF7R0c35RVFREeHg4hYWFDfa8CSZGrFtq9h8j1h0sNe/YsYOBAwfy8ccfM2DAgHP+PljqbIhR6gSp1VekVu9rTp0yB0sIIYQQwsskYAkhhBBCeJlXJrmXlpYCnqacTqfTG5v0uVOnTgGwc+dOw9QMxqxbavYfI9ZdVfPp06cDeomgpKQEgD179tT690Z5b41SJ0itviK1el9j6+zevTshISGeB7oXPP300zogX/IlX/LV7K9t27Z543DUbM8++2zA3wP5ki/5MvbXxx9/XH1M8cok96NHj/L666+TnJxcY6V5IYRorBqf/AJAjmNCiJb6+XHMKwFLCCGEEEL8RCa5CyGEEEJ4mQQsIYQQQggv80rA+uUvf0nfvn1JT09nxIgR7Ny50xub9ZnS0lKuvPJKUlNTSU9P55JLLuHgwYOBLqtRfv/735OcnIxSis8//zzQ5TRo3759DBs2jNTUVAYPHszu3bsDXVKDjPYeg3H3aSMdO4xQq9H2A6P8rhnlOGaU9xOMt6826/ffG3ffHD9+vPrP69at0/v37++NzfpMSUmJ/q9//UuvrKzUdV3XH3vsMX3s2LEBrqpxtmzZoh86dEhPSkrSd+3aFehyGjR69Gh91apVuq7r+osvvqgPHTo0sAU1gtHeY1037j5tpGOHEWo12n5glN81oxzHjPJ+6rrx9tXm/P57ZQQrIiKi+s+FhYVoWnBfebTb7YwbNw6lFABDhw5l//79Aa6qcS644AISEhICXUajHDlyhB07djBt2jQAMjMzOXDgQFB/SgFjvcdVjLpPG+nYYYRajbYfGOF3zUjHMSO8n1WMtq825/ffK41GAbKysnj77bcB+Pe//+2tzfrFo48+yuWXXx7oMlqdQ4cO0aFDB8xmz26mlCIxMZG8vDySk5MDW1wrZ6R92kjHDiPVCsbaD4KVHMf8wwj7alN//70WsFavXg3AM888w2233cbGjRu9tWmfWrJkCfv27WP58uWBLqVVqvp0UkWXriA+Z7R92kjHDiPVarT9IJjJccy3jLKvNvX3v1lj3KtXryY9PZ309HRWrVpV4++uvfZa3n77bQoKCpqzaZ+prealS5eSnZ3Na6+9FtAGh/Wp770Odp06dSI/P5/y8nLAc1A6dOgQiYmJAa6s9TLCPl2XYDt2GOU4Z6RjmxGPZ3Ic861g3Vfr0+jf/5ZO/CosLNQPHz5c/Tg7O1vv2LFj9cS1YPXnP/9ZHzBggH7s2LFAl9IsRpjEqOu6PnLkyBqTQ4cMGRLYgprAKO9xFaPt00Y6dhipVqPtB7oe/L9rRjuOBfv7WcUo+2pzf/9b3Mn90KFDZGZmUlJSgqZpxMbGsnTpUtLT01uyWZ/Kz8+nU6dOpKSkEBYWBoDNZuODDz4IcGUNu+GGG3jllVf47rvviImJwel08tVXXwW6rDrl5OQwffp0CgoKcLlcPPPMM/Tq1SvQZdXLaO8xGHOfNtKxwyi1Gm0/MMrvmlGOY0Z5P8FY+2pzf/9lqRwhhBBCCC8LvvuMhRBCCCEMTgKWEEIIIYSXScASQgghhPCy/w+KEaPfeJuVYQAAAABJRU5ErkJggg==" }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Plot posterior samples\n", "xtest = range(minimum(gpmc.x),stop=maximum(gpmc.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!(gpmc,samples[:,i])\n", " update_target!(gpmc)\n", " push!(fsamples, rand(gpmc,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 0.7.0", "language": "julia", "name": "julia-0.7" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.7.0" } }, "nbformat": 4, "nbformat_minor": 2 }