{ "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": { "collapsed": false }, "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": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "GP Monte Carlo object:\n", " Dim = 1\n", " Number of observations = 20\n", " Mean function:\n", " Type: MeanZero, Params: Float64[] Kernel:\n", " Type: Mat32Iso, Params: [0.0, 0.0] Likelihood:\n", " Type: StuTLik, Params: [0.1] 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": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Results of Optimization Algorithm\n", " * Algorithm: L-BFGS\n", " * Starting Point: [0.5,0.0,0.0]\n", " * Minimizer: [0.6566151884322408,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)|: false\n", " |f(x) - f(x')| = 3.10e-16 |f(x)|\n", " * |g(x)| ≤ 1.0e-08: true \n", " |g(x)| = 4.86e-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": { "collapsed": false }, "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.371500 \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": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XlcVGX7P/DPGQaGYRkU0dxBSVFTXNr89WhoZSalKaSZaeauZeZSUD2lZT5lpT7WU265W5kbLllqi0pampllroiKbCI7DMPsc67fH8R8RXacmXMGrvfrxStm5syZC4x7PnPf97lvgYgIjDHGGGPMYRRSF8AYY4wxVt9wwGKMMcZkpl+/fhAEQeoy2G3ggOXmrl27BkEQqvzq0aOH1GVWKCQkBCEhIXV+/pUrVzB79mz06NEDjRo1gqenJ5o2bYqIiAj85z//QVpaWoWvefPvxsPDA0FBQXj00Uexe/fu2/hpGHNfN7cjrVq1gs1mq/C4M2fO2I/r1KmTS2s8fPgwBEHA22+/7ZLXIyLExcUhKioKrVu3hkqlgr+/P7p3745Zs2bh/PnzLqnD2f766y9MnToVXbp0gUajgZeXF1q0aIFHH30US5cuRW5ubrnn3Poeo1Qq0aJFCwwdOhQ///yzBD+FPCmlLoA5RmhoKEaPHl3hY82bN3dxNc63ZMkSxMbGwmazoXfv3hgzZgw0Gg3y8vJw4sQJvPXWW5g/fz5+//13hIeHl3muh4cH3nzzTQCA2WzGxYsXsWfPHvzwww9YtGgR5syZI8WPxJjklEolrl+/jgMHDiAyMrLc42vWrIFSqYTVapWgOtfJy8vD8OHDcfDgQTRq1AgDBgxA+/btYTabce7cOSxbtgyffPIJfvrpJ/Tr10/qcutEFEXExMRg8eLFUCqVePDBB/Hoo4/Cx8cHWVlZ+PXXXzFr1izMnTsXV69eRVBQUJnnN2nSBNOnTwcAGAwGnD59Grt378aePXuwdetWPPXUU1L8WPJCzK0lJSURABo4cKDUpdRacHAwBQcH1/p5K1asIAAUGhpKf/zxR4XHnD9/nqKjo+nIkSPlXlOlUpU7/sCBAyQIAvn4+FBxcXGta2LMnZW2Iw8++CAFBARQdHR0uWNMJhMFBQXRkCFDCACFhYW5tMZDhw4RAJo3b55TX8disdCDDz5IAGj06NFUWFhY7pjr16/TuHHjaNeuXU6rIyIigpz5Fv3aa68RALrnnnvo8uXLFR5z4sQJ6tevH6Wmppa5v7J//88//5wAUEhIiFNqdjccsNxcbQPWu+++SwBo+vTp5R6bN28eAaCZM2fa70tPT6e5c+fS/fffT02bNiUvLy8KDg6madOmUWZmZoWvYTKZaOnSpXTvvfeSn58f+fr6UufOnWnWrFmUl5dnr7mir+oaz7y8PNJoNOTt7U2XLl2q9ue1WCxlblcWsIiIOnXqRADo999/r/a8jNUnN7cjkydPJi8vL8rOzi5zzLZt2wgA7dy5s9I32OLiYpo3bx6FhYWRSqWixo0bU2RkJP3yyy/lji1tbw4dOkRbtmyhnj17kre3NzVv3pxeeukl0uv15Y6t6CspKcl+nMlkosWLF1PPnj3Jx8eH/Pz8qE+fPrR79+4a/y7Wrl1rD5s2m63KY41Go/37gwcP0rhx46hjx47k6+tLvr6+dPfdd9PKlSsrfC4AioiIoLS0NBo7dizdcccdJAgCHTp0iIgqD1gWi4WWLFlC4eHh5O3tTRqNhvr160d79+6t8c946dIl8vDwoGbNmpX7d76VKIpktVrL1V7Rv7/NZiNfX18CUO15GwIOWG6utgHLZrPZP53t2bPHfv/Ro0fJw8ODwsPDyzQamzdvJl9fXxoyZAjNmDGD5syZQw899BABoPbt21NBQUGZ8xsMBvv5O3ToQC+99BK98sor9OSTT5JaraY///yT8vPzad68eRQQEEABAQE0b948+1dp41KZVatWEQAaM2ZMzX9JN+GAxVh5N7cjx48fJwC0dOnSMscMGjSImjVrRhaLpcI3WKPRSL179yYA1KtXL4qNjaVx48aRj48PKZVK2rFjR5njS0PTU089Rb6+vjRq1CiaNWsWde7cmQDQqFGj7MceOnSIxo4daw8lN7cZ+fn59tfv168fAaCePXvSSy+9RFOnTqU2bdoQAPrf//5Xo9/Fv/71LwJA33//fa1+hwMHDqTQ0FB69tlnKTY2lqZMmULBwcEEgGbPnl3ueADUtWtXatOmDXXv3p1mzJhBU6dOtffKVxSwRFGkqKgoAkAdO3akOXPm0NSpUykwMJAA0Mcff1yjWt944w0CQG+99Vatfsaba68sYPn4+HDA+gcHLDdX2jCGhoaWaXRu/tq3b1+Z56SkpFDjxo0pKCiIrl+/TgUFBRQSEkJqtZrOnTtX5tjMzEwqKioq97obNmwgALRgwYIy97/66qv2AHTrp56CgoIy56rLEOG4ceMIAK1du7ZWz7v5NXmIkLGybv2gdtddd1F4eLj98bS0NPLw8KA5c+YQUcVvsPPnzycA9Oyzz5Ioivb7T58+be/N0mq19vtLA1ZAQABdvHjRfr9er6eOHTuSIAiUnp5uv7+6IcLS0PD222+XeX2tVkv33HMPeXl5lTlfRSwWC3l6epJSqSSDwVDlsbe6evVqhecbMGAAeXh4UHJycpnHSnvgxo0bV66tJKo4YG3cuNEeMk0mk/3+1NRUatasGXl6elZYx6369+9PAOjgwYM1/fHK1c5DhNXjgOXmqhpuK/16+eWXyz1v+/btBIAeeeQRGjlyJAGgZcuW1fh1RVG0d02XslqtpNFoKCAggPLy8qo9R10C1qBBgwgA7d+/v9xjf/zxR7lw+c0335R7TQ8PD/vjb7zxBg0bNow8PDwIAC1ZsqRW9TBWH9wasBYtWkQA6OTJk0REtGDBAgJg/wBW0Rts+/btydPTs9x8HSKiKVOmEADatGmT/b7SgDV37txyx5c+dnMve1UBy2azUePGjenOO+8sE65K7dmzp0a9WDdu3CAA1Lx58yqPq40dO3YQAFq/fn2Z+wFUOBRbqqKAVTp68Ntvv5U7/v333ycA9O6771ZbU2kv4c3BttRPP/1Urh29dS4rAGrSpIn98djYWBo4cCABIIVCQdu3b6+2hoaAryKsJwYOHIj9+/fX+Pjo6GhMnDgRq1evBgA8+eSTmDZtWoXHxsXFYeXKlTh16hTy8/PLXMJ9/fp1+/cXL16EVqvFI488gsaNG9fxJ6kaVbHxwKlTp/DOO++UuW/KlCl44oknytxns9nsxykUCjRu3BgPP/wwXnzxRQwZMsTxRTPmZsaMGYPXX38da9euxd13343169fj/vvvR5cuXSo8XqvV4urVq+jcuTNat25d7vF+/fph5cqV+Ouvv8pd7dyrV69yx5eeo6CgoEb1JiQkID8/Hy1btizXBgBAdnY2gJI2ylmKioqwaNEi7Nq1C1euXEFxcXGZx29uK0u1a9eu3NV5Vfnzzz+hVqtx3333lXus9GrGv/76q9rzVNWOHjx4EP/5z3/K3Oft7Y0+ffqUuS83N9f+uy5d7mbo0KGYPXs2+vbtW20NDQEHrAYsKirKHrBefPHFCo9ZvHgxXnnlFTRt2hSPPvooWrduDbVaDQBYunQpTCaT/djSxrBVq1ZOq/mOO+4AAKSnp5d7bOLEiZg4cSKAkjVz+vfvX+E5VCoVjEaj02pkzN01a9YMkZGR2Lx5M4YMGYLLly/jlVdeqfR4rVYL4P/+Pm9VulRMYWFhuccCAgLK3adUlrw1VbYe163y8vIAAOfOncO5c+cqPe7W0HOrJk2awNPTE7m5uTCZTFCpVDV6fbPZjH79+uHUqVPo2bMnxowZgyZNmkCpVOLatWvYsGFDmbayVGW/r8potVq0adOmwseq+h1X9LoXL15Eeno6wsLCyjy2YMECLFiwAACwfv16jBs3rsJzhIWFOTWw1ge80GgDlZeXh8mTJ8PPzw8qlQrTp08v1/hYrVa8++67aNmyJc6dO4cvv/wSH3zwAd5++23MmzcPZrO5zPGNGjUCUHH4cZQHHngAAHDo0CGnvQZjDBg/fjzy8/MxYcIEqNVqPPPMM5Ueq9FoAACZmZkVPl56f+lxjlZ63ujoaFDJ1JcKv9atW1fleZRKJe677z5YLJZaLZi5e/dunDp1ChMnTsSpU6ewfPlyLFiwAG+//TYee+yxSp9X25XaNRqNQ37H3I66BgesBmrSpElIS0vDp59+ioULF+LSpUt4+eWXyxyTk5ODwsJC9O7dG02bNi3z2MmTJ2EwGMrcFxYWBo1Gg99//x35+fnV1uDh4VHjT6ilnnrqKfj7+2Pbtm1ITEys1XMZYzUXGRmJ5s2bIz09HdHR0VW+cWs0GrRv3x6XL1+u8ANWfHw8ANzWrhIeHh4AKu7V6ty5MzQaDU6ePAmLxVLn1wCACRMmAADee++9KofSANh7pa5cuQIAFU4xOHLkyG3Vc7OePXvCYDDgxIkT5R6rze947NixUCgUWLVqFXJychxWHyuLA1YD9PnnnyMuLg5PP/00xo4di5dffhkDBw7EmjVrsH37dvtxzZo1g1qtxqlTp6DX6+335+fn46WXXip3XqVSiSlTpqCwsBAvv/xyuYawsLAQOp3OfjswMBA5OTm1Gq4LDAzE+++/D5PJhEGDBuHPP/+s8Liazt1gjFVMqVRiz5492LlzZ7k5ORUZO3YsLBYLXn/99TLB5OzZs1i3bh0CAgIwdOjQOtcTGBgIABVugaVUKjFt2jQkJyfjlVdeqTBknT17FllZWdW+zpgxY9C3b18cPnwY48aNQ1FRUbljMjMzMWnSJPu81+DgYADA0aNHyxwXHx+Pzz//vPofrobGjh0LAHj99dfL/Izp6elYsmQJlEolnn322WrPExYWhtmzZyMrKwuDBg2yB8RbcTt6e3gOVj1x+fLlKvfoKn0sISEBM2fORNu2bbFixQoAJd3U69evR3h4OCZPnoz7778fbdq0gUKhwAsvvIDFixeje/fuGDx4MLRaLfbt24fg4GC0bNmy3OvMnz8fx48fx6ZNm3D8+HEMGjQIKpUKV69exf79+3H06FH7J6yHHnoIJ0+exODBg9G3b194eXmhT58+5SZT3urFF19EcXEx3njjDdx9993o3bs37rnnHvj7+yM3NxcXLlzAkSNHoFKpcO+999btF8oYw7333lvjv6GYmBh8++232LRpEy5cuICHH34Y2dnZ2LJlCywWCzZu3Ah/f/8619KpUye0bNkSX3/9NXx8fNC6dWsIgoBp06YhICAA77zzDk6dOoVPPvkE3377LSIiItC0aVOkp6fjzJkzOH36NI4dO4ZmzZpV+TpKpRK7du3C8OHDsWHDBuzZswePPvoo2rVrB7PZjPPnz+Pw4cOwWCz2CfuDBw9GSEgIPvzwQ5w9exZdu3ZFQkIC9u7di6FDh2LHjh11/rlvNmbMGMTFxWH37t0IDw/HE088geLiYmzduhW5ublYvHgx2rdvX6NzLVy4EBaLBR9//DHCwsIQERGB8PBw+1Y5f/31F06ePAmNRlNuuzFWQxJdvcgcpCbLNJT+M5tMJurVqxcpFAqKj48vd65vvvmm3ArGZrOZ/vOf/1CHDh1IpVJR27Ztafbs2VRUVFTpMgtGo5EWLVpEPXr0ILVaTX5+ftSlSxeaM2eOfVFAIqKioiKaNGkStWjRghQKRa23wUhISKAZM2ZQt27dSKPRkFKppKCgIOrbty/Nnz+fUlJSyj2nqoVGGWuoartgMSpZB0mn09Fbb71FHTt2JC8vL2rUqBENGjSo3GX+RGVXcr/VunXrCACtW7euzP3Hjx+niIgI8vf3r3Ald6vVSitXrqR//etfpNFo7G3WY489RsuXLyedTlejn4+oZCma7du309ChQ6lly5bk5eVFPj4+1LVrV5oxYwadP3++zPFXr16l6Ohoatq0Kfn4+NC9995LX3/9daXLS+Cf9awqU9VK7osWLaJu3bqRSqUif39/ioiIqNVq9Tc7efIkTZw40b4CvaenJ91xxx30yCOP0JIlSypcRqKyf39WlkBUzSAzY4wxxhirFZ6DxRhjjDHmYBywGGOMMcYcjAMWY4wxxpiDccBijDHGGHMwDliMMcYYYw7GAYsxxhhjzMEcErD0en251b4ZY8ydcDvGGHMkhwSsixcv4u6773a7nbVrsuu4HLlj3Vyz67hr3VKrSTvmLr9bd6kT4FqdhWt1vNrW2aCHCGu70bBcuGPdXLPruGvd7sBdfrfuUifAtToL1+p4ta2zQQcsxhhjjDFn4IDFGGOMMeZgHLAYY4wxxhyMAxZjjDHGmINxwGKMMcYYczAOWIwxxhhjDsYBizHGGGPMwThgMcYYY4w5WLUBy2QyYfr06ejQoQPuuusujB492hV1McYaCJNNlLoExhirkk0k5Jlq9xxldQe89tprUCgUuHTpEgRBQEZGRl3rY4yxMoxWG1K0BnQM9JO6FMYYq5BNJGy9SmjtocCdtXhelQGruLgY69atQ1paGgRBAAC0aNHidupkjDEA/4SrQgOsRFKXwhhjFRKJsD2JkFBIaB1Yu+dWOUR45coVNGnSBAsWLMA999yDvn374qeffqr0eJ1OB61Wa/8ymWrZn8YYaxAMFhuSOVwxxmRMJEJcEuFCQd3aqSp7sCwWC65evYouXbpg4cKFOH36NB555BGcP38eTZs2LXd8REREmdsxMTGIjY2tU2GukJ+fL3UJdeKOdXPNriP3uo02QqbJBrFMm+UvVTmMMVYOEWH3NcLZ/Lp/CKwyYAUHB0OhUODZZ58FAHTv3h3t2rXDuXPn0K9fv3LHx8fHo0ePHvbbKpUKKpWqzsW5QmBgLfv8ZMId6+aaXUeudRdbrMjSGqHy5p4rxpg8GayE79MIp/Nur52qMmAFBQXh4YcfxoEDBxAZGYnk5GQkJSUhLCyswuP9/Pyg0WhuqyDGWP2kM1uRVmS4peeKMcakJRIhVQdc0RKuFAHXiwmOaKaqvYpwxYoVGD9+PGJjY+Hh4YFVq1bxRHfGWK0Uma1I53DFGJMJIsIVLXA6j3CpkGCyOf41qg1Y7du3x+HDhx3/yoyxBkFrsiC9yOiQT4SMMXY7sg0lQ39/5xK0Fue+Fq/kzpgTxcXFoXv37lCr1ejevTvi4uKkLsmlCowcrhhj0rtUQFibIOKz8yKO3nB+uAJq0IPFGKubuLg4REdHQxAEEBHOnDmD6Oho7NixA1FRUVKX53T5Rgtu6DhcMcakQUQ4nw8cuSHihsH1r889WIw5yTvvvGMPV0DJH7sgCJg/f77ElTlfnsHM4YoxJokiM+GPbMKn50RsS5ImXAHcg8WY01y6dMkerkoRERISEiSqyDVy9GZk6XmRYcaYaxRbCNeKgKQiQlIRIVcmzQ8HLMacpGPHjjhz5kyZkCUIQqXLnNQH2XoTsvVmqctgjDUQf2QTvk0VZXmFMg8RMuYk8+bNsw8LArAPF86bN0/iyip2uxPyM4s5XDHGXCc+Q8Q3KfIMVwAHLMacJioqCjt27EB4eDi8vb0RHh6OuLg4DBs2TOrSyimdkH/mzBkYjUb7hPyahCwiQobOiFwDhyvGmPMREfYmizh0XabJ6h88RMiYE0VFRbnFFYNVTcivqv6ScGVCgckF1zwzxho8q0jYcRsbMLsSByzGWJ0m5BMR0nVGaE1WZ5fHGGMwWgmbrxCSdfIPVwAPETLGUDIhv3SuWKmqJuSLREgr4nDFGHONQjNhTYLoNuEK4IDFGEPtJuSLREjVGlBk5nDFGHO+G3rC6osiso1SV1I7HLAYYzWekG8TCSmFBhRbnLAzKmOM3eKqlrDukogiN5zmyXOwGGMAqp+QbxUJKVo9jFbRhVUxxhqq07mEPckibO4zKlgGByzGWLUsooiUQgNMNg5XjDHnIiLEZxAOZ7hpsvoHByzGWJXMNhEpWgPMMgpXBQUF6Nevn/22Xq/H1atXkZWVhcDAQPv9hw8fRmRkJDp27Gi/79ixY1Cr1a4slzFWQ0YrIe4a4VKhe4crgAMWY6wKJpuIlEI9LDJbKrlRo0b466+/7LcXLVqE+Pj4MuGqVJcuXXDy5ElXlscYq4NMPWHLVRF5MtlL8HbxJHfGWIWMVhuSZRiuKrJu3TpMmDBB6jIYY3X0dy5hdUL9CVcAByzGWAUMFhuSCw2wukG4OnbsGHJzc/HEE09U+HhCQgJ69eqFe++9F8uWLav2fDqdDlqt1v5lMtWjFp8xGTqU6YG4ayIs8pmF4BA8RMgYK6PYYkWa1ggbyT9cAcDatWvx3HPPQaks35z16tULaWlpCAgIQFpaGiIjIxEUFIQRI0ZUer6IiIgyt2NiYhAbGwsAyM/Pd2zxTuIudQJcq7O4S60/3fDAkXQL1GqD1KVUS6vVIi+v8hR46xQFDliMMTud2Yq0IoNsd6e/VXFxMbZs2YITJ05U+LhGo7F/37p1azzzzDM4cuRIlQErPj4ePXr0sN9WqVRQqVT22xXN85Ijd6kT4FqdRe61/pAm4pyBoFYb3OLCE41Gg8DARjU+nocIGWMAAK3JglSt+4QrANi2bRvCw8PRqVOnCh/PyMiAKJZ84iwqKsLevXvRs2fPKs/p5+cHjUZj/7o5XDHGHONguohfMt2osakDDliMMRQYLUgvMsLdmrs1a9aUm9w+ceJE7NmzBwCwY8cOdOvWDd27d0fv3r0xYMAAjBs3TopSGWP/OHxdxM833K21qT0eImSsgcszmJFZbHK7cAUAR44cKXff6tWr7d9Pnz4d06dPd2VJjLEq/Jwhuv0CojXFAYuxBixHb0aWnq+SY4w5l0iE71IIJ3MaRrgCOGAx1mBlFZuQYzBLXQZjrJ4zWAlbrxKSihpOuAI4YDHW4BARMotNyDO64fb0jDG3kmMkfHW5fi0gWlMcsBhrQIgIGToTCkwcrhhjznW5kLA9SYTRJnUl0uCrCBlrIIgI6UVGWYWrH/buwbCIB9xiDRzGWM0QEX7OEPHVlYYbrgDuwWKsQRCJkFZkhM5slboUux/27sHMcaMhCALITVaNZ4xVrchMiLvW8OZbVYQDFmP1nE0kpBYZoLfI66Pkso8WcrhirB65VEDYlSxCL5/PcZLigMVYPWYVCalaAwxWeYUrALh2JZHDFWP1gE0k/JBO+C2L3HI9PWfhOViM1VMWUURyod4l4ap0LlXP1k0xLOIB/LB3T7XPCQntAEEQnF4bY8x50nSElRdFHOdwVQ4HLMbqIbNNRHKhASZb5Tu/O0rpXKrEC+dgNpmQeOEcZo4bXW3IeuHV10BEHLIYc0NmG2Ffqog1CSKyDFJXI08csBirZ8wiIblQD7MLwhVQfi5VaWhavmhhlc8b8MQQLF33BTp2uQve3t6uKJUx5gCXCwnLzos8JFgNnoPFWD1isNpww2iDl7frmr2K5lIREZIuJ1b73AFPDMGAJ4agS5C/s8pjjDlIqo5wLJNwvoBjVU1wwGKsntBbrEjVGmFzcdsXEtoBiRfOlQlZgiCg3Z0dXVsIY8zhzDbC33mE37MJmTwUWCs8RMhYPVBktiJFa4BNgqvybp1LVTpc+MKrr7m8FsaYY4hE2J8qYvEZEXtTOFzVBQcsxtxcocmCNK0BokS99jfPpfJSqdCxy134eP2XeOTxwdIUxBi7bftSCcezCCb5rfDiNniIkDE3lm8044bOJPlE09K5VIwx9/dblojfs6VuVdwfByzG3FSO3owsfQPcop4x5jSXCwkH0jhcOQIHLMbcUGaxCbkGs9RlMMbqkSwDYVuSKNl0g/qGAxZjboSIkFFsQoHRInUpjLF6pNhC2HxF5DlXDlTjSe7vvPMOBEHA2bNnnVkPY6wSIhHSi4z1JlxZLBYc3PctZowdJXUpjDVoFpGw5Sohn2ccOFSNerBOnTqF48ePo23bts6uhzFWAZFKNm0utrj/x8vraanYtmEd4r7ahJysTNzVvafUJTHWYOkshK8uE67reVzQ0aoNWCaTCS+++CK++uor9O/f3xU1McZuYhVLwpUrNm12FiLC778exaaVy3D4wD74+Pph8PCn8dSY59Gpazepy2OsQcoyEL68LKKQp3M6RbUBa+7cuRg9ejTatWtX7cl0Oh20Wq39tkqlgkqlur0KGWvALDYRKVrXbNrsDBaLBft37cCG5Z/hwpnTuLNTZ7z14RI8Hj0Cvn5+UpfHWIN1ubBkQjvPuXKeKgPWsWPH8Pvvv2Phwqo3bS0VERFR5nZMTAxiY2PrXp2T5efnS11Cnbhj3Vxz7ZlFQqbRBmste+6NBqNzCqpNDUYD9mzZjE0rlyEjLRX/L6I/Pv1iK+7v+6B9xXeD4ealoXkvQsZc5WS2iO9Sia8WdLIqA1Z8fDwuXrxo771KS0vDwIEDsXr1agwaNKjC43v06GG/7Q49WIGBgVKXUCfuWDfXXHMGiw05WgM8vQmedXi+Wq12eE01YTQYsHXDWqz531Lk5WTjsaFR+OyLLQi7q6sk9TDGyjp0XUR8BicrV6gyYL322mt47bX/208sJCQEe/fuRdeuFTeWfn5+0Gg0jq2QsQZGZ7YircgIUYJ9BevKbDZj+6b1WPXfRcjLycaQEc9g0sw5CG4fKnVpjLF/HOZw5VK8DhZjMlJosuB6kVHyrW9qShRFfLdzOz5duABpyckYPHwkps6J4WDFmMzEZ4g4zOHKpWoVsK5du+akMhhjeQYzMoul31ewpk7++gs+nPsGzp3+Ew8Nehz/2/g1OnTuInVZjLFbxGeIOHTdXVqW+oN7sBiTgWy9Cdl697hWOj0lGR/NexM/7N2Nu3r0xMY9+3H3/3tA6rIYYxX4mcOVZDhgMSYhIsKNYhPy3WB1dqPBgLWffYzVHy9BQKPGWLhsFR6PHgGFosYbQjDGXOhIhoiDdQxXFw7uxNHVC5CTnICg4DD0mfgmOj80zMEV1m/cMjImESJCus7oFuHq6MEfMfTB+7FyyUcYPWka9h77A4OHj+RwxZhM/ZrtgZ9uI1xtjRmOUFM6Yu8LRqgpHVtjhuPCwZ0OrrJ+4x4sxiQgEiFNa4BO5lvfZGdm4oM3Y7FvVxzu6/Mglm/egXZ3dpDFKtekAAAgAElEQVS6LMZYFX7OEHEk2wN1Xa3l6OoF6BschC+f7AVBEDCpZzBG7TqFX9Ys4F6sWuCPn4y5mFUkJBfKO1wREXZ9/SWG9LkXx4/E4/3PVmJt3DccrhiTuZ9vY1iwVE5yAvq1bWJfFFgQBPQPboLsawmOKLHB4IDFmAtZbCKSC/Wy3lcwIz0NU56Owr9fmoaIAQPxzS8nMWTEM/bGljEmT/EOCFcAEBQchsMpuaB/1uIjIhxKzkXTkLDbPndDwkOEjLmIySYipVAPi0z3pyAi7N7yFd5/Ixa+fn5Y/tU2PDhgoNRlMcZqwJFLMfSZ+Ca2xgzHqF2n0D+4CQ4l5+JoSg5GfLTCIed3BHeYhM89WIy5gNFqQ7KMw1VeTg5mjB2Ff780DQ8NisSuI8c5XDHmBogI+1IduxRD54eGYcSH25CkboUPTiQjSd0KIz7ajs79hzrsNW6Hu0zC5x4sxpxMb7EhVWuATaZb3xz58Xv8e8YLEEUbPl7/JR55fLDUJTHGasAmEnZeI5zNd3zb0vmhYbLrESrlLpPwuQeLMSfSma1I0eplGa7MJhPe/3cspj7zFDp3C8fO+OMcrhhzEyYb4cvLzglXcucuk/A5YDHmJFqTBalaA1w9Knhw37cYFvEAerZuimERD+CHvXvKHXPtSiJGDXoEW9avwRvvfYgVX+9A0zvucG2hjLE60VkI6y+JuFrU8MIV4D6T8DlgMeYEBUYL0iXYtPmHvXsQM2U8Ei+cg9lkQuKFc5g5bnSZkPXNtq/x1EMPwqAvxub9B/HspKl8hSBjbuKGnrA2QUSGXupKpNNn4ps4kpyDUbtOYdWpaxi16xSOpuTgXxPfkrq0MjhgMeZgeQYzMnSuD1cAsOyjhRAEocwnO0EQsHzRQhgNBsybPQOvvTAZA54YjK0//ozO3cIlqJIxVlsiEY5kiPj8oog8k9TVSEvuk/BL8SR3xhwoR29Gll661u/alUR7uCpFRLiaeAmjBj2Ca1cSMf+/nyLq2THca8WYm8gzlkxmTy1umEOCFZHzJPxS3IPFmINkFZskDVcAEBLaoVxwEgQBNqsVRoMem/cfRPTo5+pFuAoJCUGnTp3Qo0cP9OjRA1u2bKnwuAULFiA0NBShoaF46y15DSEwVp0/sgkrLogcrtwQ92Ax5gA3io3IM0i/afMLr76GmeNGlxsm7NbrbqzauhP+mgCJK3Ss7du3o2vXrpU+/vPPP2Pz5s34+++/oVQq8a9//Qt9+vTBwIG8xheTN62Z8E0KIbGQg5W74h4sxm4DESFDJ49wBQADnhiCD1euRWhYJwiKkj/vwcNH4qt9P9W7cFUTW7ZswfPPPw9fX1+oVCqMHz8emzdvlrosxqr0Rzbhs/Mihys3xwGLsToiIlzXmZBvlEe4KtW2XTuYjAZoAgLw+bZdWLhsVb0YEqzIs88+i27dumHixInIzs4u93hKSgqCg4Ptt0NCQpCSklLlOXU6HbRarf3LZGrgM4qZyxSYCBsvifgmRYRJvtuVshriIULG6oCIkK4zQmuySl1KGd9/sxtvTJ+Ctu1CsXr7HrQODpG6JKf5+eef0bZtW1gsFrz55psYO3Ysvvvuu3LH3Rwub70AoCIRERFlbsfExCA2NhYAkJ+ff5tVu4a71AlwrQAgEnAqX4Gfs5SwiI45p8FgdMyJXMBdatVqtcjLq/wfKDAwsMxtDliM1ZJIhPQiI4rM8glXoihi2UfvY/miDzDgiSF479OV8PH1lbosp2rbti0AwNPTEzNnzkTHjh0rPObatWv228nJyfbnVSY+Ph49evSw31apVFCpVPbbtzaicuUudQINu9bLhYQDaSKyjYBS5dg3ZbVa7cCzOZc71KrRaBAY2KjGx/MQIWO1IBIhTWuQVbgq1ukwa/wYrFj8IV5+Yy7e+2xVvQ9XxcXFKCgosN/evHkzevbsWe644cOHY8OGDSguLobJZMLatWsxcuTIKs/t5+cHjUZj/7o5XDHmKDlGwpeXRXxxuSRcsfqHe7AYqyGRCKlaA4ot8pkccT01BS+OGYm0a9fwycbNeOixSBgMBqnLcrrMzExER0fDZrOBiNC+fXts3LgRABAZGYn58+fjnnvuQb9+/TBixAh069YNADBy5Eg89thjUpbOGrhCM+GXG4STOeTybbSYa3HAYqwGbGJJuNJb5ROu/jzxG15+fhS81T748rsf0LHLXVKX5DLt27fHn3/+WeFjt87Dmjt3LubOneuKshir1PViwq+ZhPMFHKwaCh4iZKwaNpGQIrNwtXf7Fowb9jhCQjvg6wOHGlS4YsxdEBESCgjrEkSsuijibD6HK3d04eBOfD7qbjzapQXu7tkDcXFxNXoeByzGqmAVCclaPQwyCVeiKOKT999F7LRJeDxqOFZv343AoCCpy2KM3eJaEeHzi4TNV0Qk6zhVuasLB3dia8xwhJrS8eq9wfDLv47o6OgahSweImSsElYiJBfqYbI56Lrp22Q0GPDvl6Zh/+44zHrrHUx4aWa9Xd+KMXeVZSD8kM4rsNcXR1cvQN/gIHz5ZC8IgoBJPYPx7O5TeG/Bu4iKiqryuRywGKuAxSbihlGEUiWPcJWTlYWXnhuJS+fPYem6LzDgiSFSl8QYu0mBiRCfQTidx8OAcnbh4E4cXb0AOckJCAoOQ5+Jb1a5aXROcgIm3Bds/zArCAL6tW2CxScvVvtaPETI2C3MNhHJWgMsMmklEy+cx8iB/ZGRloYNe/ZxuGJMJmwi4Xw+YVOiiI/Pivgzl8OVnN083Bd7XzBCTenYGjMcFw7urPQ5QcFhOJySW2Zv18MpuejcqVO1r8c9WIzdxGwTkVyol024+vXwQcwa/xxatm2LZV9uRYtWraUuibEGL8dI+DOH8FcuoVg+S+KxalQ03Ddq1yn8smZBpb1YfSa+ia0xwzFq1yn0D26Cwym5OJKcg7j/rqr29bgHi7F/mKw2WYWrbZvWY+rIaPS87358sfcAhyvGJJRtIBzN9sCy8zZ8ek7EL5kcrtxNTnIC+rVtUma4r39wE2RfS6j0OZ0fGoYRH25DkroVPvo9GfrA1oiLi8OwYZUPK5biHizGABitNqRoDbDKIFyJooilC97Bmv/9FyPHTcTr730IpZL/VBlztWwD4Wx+yTBgthEwGDzgBju6sEqUDPelY1LPkjlVRIRDybloGhJW5fM6PzQMnR8ahocDi9C3Xc23yuFWmzV4BqsNKYUG2GqwEbCzmYxGvDF9Kg7s2YmYd9/Dc1Ne5CsFGXOhfBPhbF5JsMqs/5siNCi3DvcdSs7F0ZQcjPhohVNejwMWa9D0FhtStfIIV/m5uZg+ZiQunv0b/127iSezM+YiOcaSBUEvFABpxdK3Bcw5Sof7flmzACdOJKBpSBhGfLQCnfsPdcrrccBiDVax2YrUIiNEGYSr5KtXMO2Zp1CkLcS6nXsRfve9UpfEWL1FREgtBhIKCBcLCLkmqStirlI63FdbV7SEX9JV+D6EajyqwAGLNUg6sxVpRQZZXFJ9+uQJvDj6aTRqHIiv9v2ENiHtpC6JsXopTVcy9Hcun1Bkkboa5i6SighfXyE8EIRaTdnggMUaHK3JgvQiI2SQrXBw37d4dcp4dO7WHZ9u2oxGgU2kLomxeuWGviRUnc0jFJilroa5m2QdYfMVQrA/8FZXE4CaX+XAAYs1KAVGCzJ08ghXm9d+jvdefxWPPD4EC5etgsrbW+qSGKsXbuhLeqnO5/PwH6u7VB3hq8uE1r7A0+0FqDxq93wOWKzByDOYkVlskjxcERE+/s98fP7xYoyZ8gJi5r8HhYKXpGPsduQYCX/nlgQrDlXsdl0vJnx5mdDcBxjZXoCnovZXc3PAYg1Cjt6ELL304wMWiwVzZ07Hnq2b8er89/D8tOlSl8SY29JbS4b+TucB6Xz1H3MAk41wRQt8k0JoqgZGhQrw8qjbUjkcsFi9l1lsQq5B+nBVrNNh1vgx+O3oz/hw5Ro8HjVc6pIYczt6KyGxELhQQEgsJNg4V7HbIBLhuh64qgWuFBHSdIAIoK0f8EyoAFUdwxXAAYvVY0SEjGITCozSXy6Ul5ODaaOeQlJiIlZ+vQO9H+wndUmMuY1sAyGhkHCpsGReDGcqdrsy9YS/8ghn8oBiK6BSAO38gcfaCAj1BxqranfFYEU4YLF6iYiQXmSE1uyYzcJ+2LsHyz5aiGtXEhES2gEvvPpajRcCvZ6agknDh6JIW4j1u75Fl+49HFITY/VdkpZw8DohlYf/mAMUWwhn8oHTuYQbBsBHCXQLBLo0EtDaF1A4eNeMKgOW0WjEyJEjcf78efj4+KB58+ZYsWIFQkJCHFoEY44kEiFNa4DOYnPI+X7Yuwczx422712VeOEcZo4bjaXrvqg2ZF2+eAGTRwyDl8oLm/Z+j+D2oQ6pSWoCAC8PBVQeCnh5KODpoYCnQoCHQoACAhQCIAgAEUAAbESwiQQbESw2EWYbwWwTYbSJsljolclLqq4kWCUV8f8b7PblGgm/ZBL+zitpj8ICgH4tBNwZAHg4cSuyanuwJk+ejEGDBkEQBHz66aeYPHkyvv/+e6cVxNjtsImEFK0BBqtjwhUALPtooT1cASW9Y4IgYPmihVUGrNMnT2DqM0+heavWWPV1HJo2b+6wmlxNAKD29ICvpwd8lB5Qe3o45NMeUUnQYgwoGbb58XrJ3CrGbld6cUmwulAA+CmB/i0E9AwCfJSu2d+1yoDl7e2NyMhI++3evXtj6dKlTi+Ksbqw2ESkaA0wOfgN+9qVRHu4KkVESLqcWOlzfj18EDPGjkKnruFY9tVWaAJqvgO7XCgEwM9TCX+VEn6eSnjU4TLl6giCAJWylovLsHrHIhIOXyccyyJZ7K7A3FeBiXClCDibR7imAwJVwBNtBXQPBJROaMOqUqs5WJ988gkGDx5c6eM6nQ5ardZ+W6VSQaVS1b06xmrIZBWRotXD4oTWOSS0AxIvnCsTsgRBQLs7O1Z4/Pff7MarU8bj/z3YD/9duwlqHx+H1+RMaqUHGnl7QuPlnFDF2M2SdAKOXheRz2tXsTqwiYSrRcBlbcnyCrmmkh73tn7A8HYCOjVy/NyqmqpxwHrvvfeQmJiIFStWVHpMREREmdsxMTGIjY2te3VOlp+fL3UJdeKOdTuzZqONkGWyOfxybaPBCACY+PJsxEwZbx8mLP3vxJdnw2AwlHnOnq2bsSBmNh5+fDDm//dTQBDKHeNspXXXhgDAVylAo1RAJQogvRGFesfXVpnAwEDXvRiTBb2VsD+V8Fu6J9Q1332EMTu9tWQbm7RiIMALCPUHHm4poJ0/4O2iYcCq1ChgLVq0CHFxcfjxxx/hU8Wn8fj4ePTo8X9XSLlDD5a7NuzuWLczataZrcgqMsDLSbvMqNVqPB71FLy8vLB80UIkXU5Euzs74IVXX8cjj5ftzd20chkWvvkahj83Dm99uAQeHtINfalr+I6lEIBGKk80UXvB04NXk2fOR0T4Mxf4MV2E3jEX+bIGKN9UstK6wQaM6yigje/tL6vgaNUGrCVLlmDz5s348ccf0ahR1fNI/Pz8oNFoHFYcY1Vx5b6CA54YUumEdiLCyiUf4X8LF2Dc9JcxZ+582f2h30oA0NjbE018vODJ2/QwF7mhJ+xNIaTxsgvsNmToS8KVlwcwIUxAoEqe7W2VASstLQ1z5sxB+/bt0b9/fwAlvVK//fabS4pjrDJy2fqGiLDk3XlY+7+lmPH6W5g86xXZh6sAlRJNfVTw4h4r5iImW8myC79n8yR2dnuuaAlbrxKCvEu2sfH1lG97W2XAat26dbmrpxiTEhHhRrEJ+TJYnV0URbz3+qvYvPZzxL77Pp6b+qLUJVVJrfTAHb4q+HjyVXvMdf7OJXyfLkIn/Z8sc3Oncwl7kgmhGuCpdnXfI9BVeCV35jbEf1ZnL3LQ6uy3w2azYe7M6di95Su8veQTDB/zvNQlVcpDENDM1wuNVJ6y711j9UeWgfBtCiFZxx/SWdUuxe/GiY0fISc5AUHBYegz8U10fmiY/XEiwpEbwKEMQo8mwOC2gmRXBtYGjxEwt2AVRSQXGmQRriwWC2KnTsA3277GwmWrZB2uAlRKhDb2QWNvLw5XzCVMNsL+VBErLogcrli1LhzciV1vjkaoKR2x9wUj1JSOrTHDceHgTgAlu0B8k0I4lEHo10LAEDcJVwD3YDE3YLKJSNUaZLHit9lkwuyJY3Hkpx+wZM3GclcSyoWnQkBzP2/4e/GfOHOdy4WE3ckiing4kNXQ0dUL0Dc4CF8+2QuCIGBSz2CM2nUKv6xZgPYRQ7EtiZCkBZ4MFtCjiXsEq1Lc+jJZ01tsSNMaYJXBXECDXo+Xn38WJ4/9gk83bkbfRx6VuqQK+SkFtG/ky4uEMpchIsRnlHxJ/5fK3ElOcgIm3Bds72EXBAH9g5vgxIkErL9EyDcBz94poL3G/dozHiJkslVosiBFq5dFuCrW6TBt1HCc+u04lm/eLstwpVQIaOOvRlOVB4cr5jIGK+GrK4TDHK5YHQQFh+FwSm6ZvV4PJecCzcOgtwLjwtwzXAHcg8VkKkdvRrbeJIsGu0hbiKkjo5F44QJWbd2JXvf3lrqkcvy8lGjpp4JSoUBesdTVsIYiQ0/YckVEgfQrpjA31Wfim9gaMxyjdp1C/+AmOJSci6MpOQh8eRnGhgnQeLlnuAI4YDGZISJk6EwoMMljEkdBfh4mjxiG1KQkrInbg24975a6pDIUAtDMR4XG3nyFIHMdm0g4kU34KZ1glcOnIOa2Oj80DEMXfIETmz7C8d8SQM3D0C5mGUY9NczlmzM7GgcsJhs2kZBWZECxxSZ1KQCA/NwcTB89Elk3rmPtzr3o3C1c6pLKUHko0MrfG95KXteKuU5CAeH7NBG5vDkzc5AWDzwJRbsRID0wqI2Au4PcO1iV4oDFZMFkE5GmNcAkgysFASD7xg1MGRGFwoJ8rN/1He7s1FnqksoIUCnR3Neb51oxl8nUEw6kEa4WcZcVc5wbesIXSV5QKIDnOwpo7Vt/2jQOWExyxWYr0oqMsMlgMjsAZKSnYULUYBj0emz8Zh9CQjtIXZKdAOAOXxUC1V5Sl8IaCItI+CGtZJsbefyFsvoi00DYmEjQeBJGd1DAT8bb3tQFBywmqXyjBTdctGFzTaQlX8P4qMEgIqzatktW4cpTIaC1vxpq3uqGuUiWgbA9SUSWQepKWH2T/U+4CvAChrexwM+z/sWR+vcTMbdARMjUm5BnkMdkdgBIvnIZ46MGw9PLC2t37kXjJkFSl2Tn5+mBlv5qt5/0ydzHqRzCvlQRFnmM2rN6JMdI2JBI8PcExnQQIMjnbcChOGAxl7OJhPQiA3QymcwOAJcTLmJC1GBoAgKwJu4bNGveAgaDPD62B6m90NSHt7phrmGyEb5JJpzNl0u/MqtPcv8JVz5KYMydAnyUAmT0OduhOGAxl5LbZHYAuHj2DCY+NQRNm7fA6m270aRpU6lLAgAoBAEt/VTQqDylLoU1ACIRzuQBhzNE5PMVgswJ8owlw4LeCuC5DgJ869mcq1txwGIuozNbkS6jyewAcObPPzB5xDC0CW6HVdt2olHjQKlLAgB4eSjQmpdgqJTRaMTIkSNx/vx5+Pj4oHnz5lixYgVCQkLKHHf48GFERkaiY8eO9vuOHTsGtVrt4orlSyTC33nAkQxeeoE5T3IRYcvVkp6r5zoK9W5Ce0U4YDGXKLSIMGoNspnMDgB/HPsV00YNR4fOXbDi6+3w1wRIXRIAwNfTA614vlW1Jk+ejEGDBkEQBHz66aeYPHkyvv/++3LHdenSBSdPnpSgQnmzif8Eqxsi8jhYMSf6M4ewN5XQ1hcY0V6AWtkw2jYOWMypxH9WZs8zi5BTp8Gx+EN46bln0K3XPfh009fw9fOTuiQAQKC3J+7wVfF8q2p4e3sjMjLSfrt3795YunSphBW5j3wT4Y8cwp85hGKr1NWw+kwkwo/phGNZQK8gILKNAI8G1LZxwGJOY7GJSCsywmCVz2R2ADh0YB9mjR+D+/s8iI/XfwlvGSQ/AUBzPxUae/P6VnXxySefYPDgwRU+lpCQgF69esHDwwPjxo3DCy+8UOW5dDodtFqt/bZKpYJKpXJova4mEiGhADiZQ7iq5fWsmPOZbIS4JEKiFnistYD7mqLBfXDkgMWcQm+xIa3IAKsor6Z8384deO2FSeg3MBIfrVwDLxm8cXoIAlr7e8PXi/8c6+K9995DYmIiVqxYUe6xXr16IS0tDQEBAUhLS0NkZCSCgoIwYsSISs8XERFR5nZMTAxiY2MBAPn5+Y4t3klK69RZgdP5Hjidr0CRVZ5vbgaDUeoSaoxrrZlso4A96UoUWQVEt7GgvR/BWEU57vJ71Wq1yMur/AKtwMCyc3i5RWcOJ7fFQ0vFfbkJc2dNx+DhT+Pdj5dBqZT+f38vDwXa+KuhUiqkLsUtLVq0CHFxcfjxxx/h4+NT7nGNRmP/vnXr1njmmWdw5MiRKgNWfHw8evToYb99aw/WrY2oHKXrBfxVGIDz+QQbAfAE1DK+GNWdLjrgWitHRDiZAxxIIwR5A5PuFBDk7V2j57rD71Wj0SAwsFGNj5f+HYbVG0SEG8Um5Bvlt6jJxhWf4YO3XsfTz0/Amx8shkIhfaDx8fRAG3817ydYR0uWLMHmzZvx448/olGjihu9jIwM3HHHHVAoFCgqKsLevXsxYcKEKs/r5+dXJpi5A5tISC0GrmgJlwoJ1/I8oVbL7SMOq88MVsKeZMLFQuDepsCjrYQGf6EOByzmEBZRRJpWfvOtiAjLF3+Azz54D+NfmonZb70ji3kAASolWvp5y6IWd5SWloY5c+agffv26N+/P4CSnqbffvsNEydOxJAhQzBkyBDs2LEDy5cvh1KphNVqxfDhwzFu3DiJq3cMo7VkMdDEQuCajmCS158ea0BSdIQdSQSzWHKVYOdG3K4BHLCYA+gtJZs1y22+FRHhw7lvYOOKzzDj9bcwedYrsgg0QWovNPOVfu6XO2vdujWokvXUVq9ebf9++vTpmD59uqvKcolMPeFENuHvPOJtbJikzDbCwQzCb1lAW18gqp2AAC/p21i54IDFbkuewYzMYpPs5ltZrVa8PXsGdm7+Av9euAijJkyWuqR/rhT0RmNvGU+GYbIkEuFCPnAim5Csk9tfG2uILheWrG1VbAEGtBLQu1nJ7hPs/3DAYnVSsr6VEYUm+S2kYzaZ8OqUCTi0/1t8sPxzPPHU01KXBMU/Vwr68ZWCrBbon+1r4nmVdVZHFw7uxNHVC5CTnICg4DD0mfgmOj80rM7n01sJ+1MJZ/KB9v7A2A4CGqs4WFWEW3tWa+Z/9hM0ymg/wVLFuiK89Nwo/PX7b/h4w1foP3CQ1CVBqRDQxl8NtSdve8NqhohwLr9kX8Ac97iCncnQhYM7sTVmOPoGB2HCfcE4nJKOrTHDMeLDbbUOWQUmwqlcwslsgAA8GSyge2DDW9uqNjhgsVopMltxXWb7CZbKz83FlJFRSL5yBau27MQ9D/xL6pJKlmHQqKHykP6qRSZ/BivhYgFwLEtElkHqapi7O7p6AfoGB+HLJ3tBEARM6hmMUbtO4Zc1C2oUsEQiXCoE/sghXNYCKgUQ3gR4sHnD2EvwdnHAYjVCRMjWm5FrMMtuvhUAXE9NweQRw6AtLMD63d+hc7dwqUuCt1KBtho1lDJYEoLJV7GlJFSdLyAkFRFkdq0Ic2M5yQmYcF+wvZdJEAT0D26CEycSqn3u5ULCnhRCkQVo6QMMbiuga2PAy4ODVU1xwGLVsooi0ouMKLbI8zrwyxcvYPKIYfD08sKmvd8juH2o1CXxGlesSlaRcKEA+CuXQxVznqDgMBxOScekniUhi4hwKDkXTUPCqnxepp6wLYnQ2hd4JlRACx9ux+qCAxarkt5iQ3qRARaZvgOc+u04Xhw9Ai1atcHKr3egafPmUpcEfy8lWvl78xU1rJxMfck8lr9zCQZ5fl5h9UifiW9ia8xwjNp1Cv2Dm+BQci6OpuRgxEflt5UqpbMQNl8lBKqAp9sL3GN1G3jsglUqz2BGcqFetuHq4P7vMPGpIQjr0hUb9nwni3AVoPJEaw5X7BaXCgirL4pYfkHEb1kcrphrdH5oGEZ8uA1J6lb44EQyktStMOKj7ejcf2iFx1tFwparBJtY0nPF4er2cA8WK0fOSzCU2r5pA9555WU8HDkYHyz/HKoa7nflTIHenrjDV8VX1TC7JC3h4HVCarE8P6Sw+q/zQ8NqNKGdqGSrmxt64PmOAjS8YOht44DFyjBZbUgrMsIkwyUYgJJGYNlH72PZRwsxctxEvPH+R/DwkH75A16dnd0sTVcSrK4WcbBi7uHIDeBMPhAdIqCVL4crR+CAxewKTRZk6EwQZbgEAwBYLBa8++os7PhyI2b+ex4mvjxbFr1Fd/io0MTHS+oymMQKTIRz+SX7A2bopa6GsZo7l084lEGIaCGga6D0bWp9wQGLgYhwo9iEfKNF6lIqVazT4ZVJz+PXwwfx/mcrMWTEM1KX9M/WNyo09uZw1VCZbIQ/cwhn84H0YpLlEiaMVeVEdsnK7OGBQIT001jrFZ7kXs/ExcWhe/fuUKvV6N69O+Li4qo83mwTca3QIOtwlX3jBp5/MhJ/HD+G5Zu3yyZctfTz5nDVgBmthI2JhP1phDQOV8zNiEQ4kCZiXyrh/mYlK7PLYUSgPuEerHokLi4O0dHR9vVOzpw5g+joaOzYsQNRUVHljteZrUiX6arspS4nXMTUkdEQbTZs2nsAYXd1lbokCABa+XtDo+JNmxsqk43wxWVCOk9eZ27IIhLikggJhcCgNgLua8rByhm4B6seeeedd+zhCigZ+hMEAfPnzy9zHBEhq9iEVK1B1uHq+JF4jI4cAH+NBl/t/0kW4eufblYAACAASURBVEohAG00ag5XDZjJRtiUWNJrxZi70VkIGy4RrhQBI0M5XDkTB6x65NKlS/ZwVYqIkJDwf9siWEURKVoDcmS65U2pXV9/iSkjhiH87nuwae8BNG/ZSuqSoBBKNm328+KO34bKZAO+4HDF3FSBibA2gVBoLlmKoWMAhytn4oBVj3Ts2LHcGLogCAgLK9kWodhiRVKBXrZb3gAlgfB/Cxfg3y9Nw9BnRuOzL7fCz18jdVnwEAS01ajhy+GqwTLbCNtTlbymFXNLBSbChsSSuYITwgS05O1vnI4DVj0yb948+7AgAPtw4dy5c5GjNyOlUL5b3gCAyWjEq1PGY8XiDzHrrXfw9uKP4ekp/VCchyCgbYAaPp7Sr7fFpHO+AEjTc5PJ3E++ibA+kSCgpOeqkYrDlStwa1GPREVFYceOHQgPD4e3tzfCw8OxffsO3PvIIGTpTbIeEszNzsb4qME4uO9b/HftJkycMUsWV7QoFQKCA9RQKzlcMcbcT56JsP4SwUMoCVcBvEK7y1Q73pGYmIixY8ciJycHjRo1wvr169GlSxdX1MbqICoqyn7FoMFqQ3qREUVm+W55AwCJF87jxdFPw2gwYP3u7xDe6x6pSwIAKAUgWKOGisMVY8wN5RlLhgU9FcBzHXj7G1ertgdrypQpmDx5Mi5duoSYmBhMmDDBFXWx21S6UbNZplveAMAPe/fg0bu7YeiDvZGdeQMvxr4hm3DlqRDQ3NuDwxVjzC3lm4H1/4SrsRyuJFFlwMrKysKpU6cwevRoAEB0dDSSkpJw7do1V9TG6kAkQnqRATeKTZDxdCt8/81uzBw3GukpyQAAi9mM+a/MxA9790hcGeClUCA4wAeeCm6QGGPuJ89E2JzsBS8FMLajAH8OV5KoMmClpqaiZcuWUCpLRhIFQUDbtm2RkpJS4fE6nQ5ardb+ZTKZHF8xq5TJakNSgR6FJnkPCZrNZrw9e0aZ+0on5y9ftFCiqkqoPBQIDlDDy4OnJzLG3E++qWSdK08FlYQrTw5XUql2DtatE41vXWfpZhEREWVux8TEIDY2to6lOV9+fr7UJdRJRXXrrCJyTSLkOiBoNBgBAAV5uYiZOhGFBeV/BiJC0uVEGAwGV5cHAPBSCPD3VqCosOSDQX36/0PuAgMDpS6BMbdXuhSDUgE83cYCf09eVkZKVf7227Rpg7S0NFitViiVShARUlNT0bZt2wqPj4+PR48ePey3VSoVVCqVYyt2MHdt2EvrFomQWWyCTrRApZa4qGqkJl3F9DEjYdAXo01IO6QlXysT2AVBQLs7O0Ktdv0P4q1UoK1GDaWibM+Vu///wRhrGAr+WYpBIZTMufKU73KHDUaV4yDNmjVDz5498cUXXwAAduzYgZCQEISEhFR4vJ+fHzQajf1L7uHK3ZVs1KyX9UbNpQ4f+A7PRg6An78GW74/jDnz3q1wza4XXn3N5bWplR4I1viUC1eMMeYOSnuuFOAJ7XJS7TvKypUrsXLlSnTs2BELFy7EmjVrXFEX+//t3Xl4FFXa9/Hvqd6Tzr5ASEgCgQRZwyIii4DKoIgoRBaRiWziDOKo4yPoMwqKiDgi4zKvMgqiPDqMCgEZxl1Z1EFREUGBgAoJQRHZEgjZu94/WqKBkLW3Su7PdeWSlk71naa66pdTp+5Ti4KSMvadOE1xeaBeFHRzuVw8vXAB/3PTJPoNvoz/W/c2rVonMmT4CB5f9hKpHTthtdlI7diJJ154mcuvutqn9QVZTCSGOjDJhHYhhAH9dFrn+T3uKwGZ7aXPVSCp9QJtWloamzdv9kUtog50XedoqYsyV7G/S6lV4alT/OXWP/Luutf5w52zuGXmPWi/GSUaMnwEQ4aP8Ft9wRYTrUMdaAHQ0FQIIepr/0mdf32nE2GDG9opnDKhPaDIDDgDKatwkXeymIIyF44A/5fL3fc9f7pxPAdzc3nyxX/Sd/BlVcKVv4VYzcSH2CVcCSEM6ZvjOqv36yQ5YUxbhc0kx7JAEzhnPFGjk6XlfH/iNEXlgT9z8b8bPmDs7wZRUlzMirfe57Jhwz26/XfXrWXkwL50T4hh5MC+9e6dFWozkyDhSghhUJ8c1lm5T6djOIxPkXAVqCRgBTj9l7sE8wqKqKihRUYg0HWdpU89zs1jR9G1Ry/+9c562nW4wKOv8e66tdw+aQJ7d31DaUkJe3d9w+2TJtQ5ZIXbLMQ77QGxzqEQQtRHmUvnjQMu3s7T6RsLI5OVzB8NYBKwAlhZhYv9+UUcLSr1yULNjRkZOl1YyP/cNIlFc2cz5dY7ePqfrxEWHuHxGp9+dEHlHYdQvwalkXYLcU6bhCshhOH8eFrn2d06W4/AsNaKIQmaHMsCXIDP5Gm+TpaU88OpYp+NWp0ZGToTXs6MDD2+7KVaJ6LnfP8dt028gbycHBYtXc7QEdd6rc793+09p9ntmQalNYl2WIkNlrYhQghjcek6Hx2CjT/qxDpgWgdFrEOClRHICFaAcek6h04Vc+Ckby8JNnRkaOM7bzF2yCBKS0r419sfeDVcASSntD/nt7YzDUrPJzbIJuFKCGE4x4p1lu3R2fCjTt8WMDVNwpWRSMAKICW/NA495ofGofUdGXK5XDz96MNMv2EMvfr245V3N3h8vlV1pt91d50blCogzmkjOsjq9bqEEMJTXLrOJ4d1Fu/WKSyDSamKy+I1mW9lMBKwAsSJYv82Dq3PyNCJ48e45YYxPP3oAm69+16efPGfhISG+aTOujYoVUCrEDsRdglXQgjj+LnIPWr1dp5O9yj4wwWK1k4JVkYkc7D8rMKlc6iwmPyScr/WMf2uu6vMwTrfyNCu7V9x++TfU5B/gmdWrGTAZUN8XmttDUo1BQkhDpxW2b2FEMZQoet8fAg2HdIJt7pHrRIlWBmajGD50emyCvadOO33cAV1Gxla/c+XuOGqIYSEhvHae5v8Eq5qY1KKxNAgCVdCNAG7PljNc+N78nA/J8+N78muD1b7uySvyD2l89xu91yrPrHuUSsJV8YnAcsPdF3n59Ml5OSfptQVOGsJDhk+gqwN/+XLvJ/J2vDfynBVUlzM/Xfexr23TWd4xhhefuNdEpKS/VtsNcyaIinMQZDFdN7nZGVl0a1bNxwOB926dSMrK8uHFQoh6mrXB6t5deZoUkoOMqt3EiklB3l15ugmFbIKSnWy9rlYtkdHU3BTB8Xl8RpmmWvVJMiv+T5WUu7ih1PFhujIDnAwN4c7Jmeyd/dOHlj0FNf9/kZ/l1Qtq0kjMdSB1XT+3xmysrLIyMiovPy5Y8cOMjIyWLVqFaNGjfJhtUKI2ny0ZB4DkqJ5+ZoeKKW4qXsS49ds5eOl87jg0pH+Lq9Ryl06mw/Dh4d0rBqMSFSkRyF9rZoYGcHyEV3XOVZUyr78QsOEq03vvs11lw0g/8RxXv7PuwEbruxmjeSwmsMVwAMPPFBtK4q5c+f6okwhRD0cyclmUGJUlTuGBydF8fP+bD9X1nDlLp0vj+g8vVNnww86PaNhRidF92gl4aoJkoDlA6UVLnILijhUWIKrjq2tGrveXmNUVFTwxPy5/HH8aLpfeBGvvreRjt3Sffb69RFsMZEUGoS5DgtJ79mzp9pWFNnZxj1gC9FURSelsSH3aJVfiNbnHCUmOc3PldXfyTKd9T+4+NvXOmtz3Q1D/9hRMTRBwy7rCDZZconQi3Rd53hxGYdPl+KqR9PQxnRVb6yjP//MzD9MYctHm7j9L3OY8qc70OoQXvwhzGamVT3WFUxNTWXHjh1VQpZSirQ04x2whWjq+k+9l1dnjmb8mq0MTopifc5RPso9wphHF/u7tDr7qUjnvz/pfH0cTAq6R0HvGEWUXUJVcxCYZ84moKTCRU7lqFX9OrI3Zr29xtj66SeMvmwAe3d9w3MrX+em2++sU7jyx2hbpMNSr3AFMGfOnGqblM6ZM8dbZQov2rt3L3379iU1NZXevXuzc+fOap83b948UlJSSElJ4b777vNxlaKhLrh0JGP++hr7HPE8siWHfY54xjy6kgsGe3e1CE/IL9V5fb+Lxbt0ck7B5a0Uf+6suLK1JuGqGZERLA/TdZ2jRaUcKSqt8+XAszV0vb2G0nWd5f/4fyx6YDZde17IY0teILZlXJ2+19ejbQqIDbYR5ah/A9FRo0axatUq5s6dS3Z2NmlpacyZM4eRI409Yba5uvnmm5k2bRoTJ05k5cqVTJkyhc2bN1d5zqZNm1ixYgXbt2/HbDbTr18/+vfvz9ChQ/1UtaiPCy4dGfAT2nd9sJqPlszjSE42UUmphI36C9+mjMRqgitbK3pGu9vHiOZHRrA86HRZOd+fOP3LJcGGb6ch6+011MmCfG6f9Hv+et//MmHaH3l+9bo6hyvw7Wjbme7sDQlXZ4waNYpt27ZRVFTEtm3bJFwZ1OHDh9m6dSsTJkwAICMjg3379rF///4qz3vllVeYOHEiwcHB2Gw2Jk+ezIoVK/xQsWiKzm4l0a74B3YvGEvqd6v5UydF7xgl4aoZk4DlAWUuFwdPFrE/v4iSisb3tarPenuNsWvHdsZcPpBPP9zIky/+k7seeAiLxVKvbfhqtM2kFIlhDsJs9avPyKRn1/kdOHCAVq1aYTa7B+GVUiQmJpKbm1vlebm5uSQlJVU+Tk5OPuc5Zzt16hQFBQWVXyUlJZ7/AYTh6brO+mcfZECiu5XEtB7JvHxtD/onRpO/ej42mbze7MklwkZw6TrHiso4WlRKRT3nWdXkTFf1ZxYuYN+3e2nTrj3T77rnnPX2GkrXdbJe/j/m3X0nKakd+McrWSS2adugbSWntGfvrm/OmTjuydE2i6ZIDHVgM5+/gWhTIz27anf2KO/ZQb+6553vOb81cODAKo9nzpzJrFmzyM/XKCry/6oLdVFUVOzvEurMaLWeLoc9JzW2HTdxNHcPgy5KOqeVxJYt2RQVFfm9VqMwSq0FBQUcO3b+QZTIyMgqjyVgNYCu6xSUlnO4sISyxlwLrEFt6+01VNHp0zw488+8/so/GZ05iXseegSb3d7g7dV1DcOGsps1Woc6sATonYzeUlPPLglY0Lp1a/Ly8igvL8dsNqPrOgcOHCAxMbHK8xITE6tcNszJyTnnOWfbuHEj6em/tiWx2WzYbDbCdB2HoxCHw+HRn8VbjFInBH6tBaU6u/Phm6MWDpx2H4vah4ErMZUNuT9wU/ekys/rmVYSgfAzBUINdWWEWkNDQ4mMDK/z85vXWcsDTpaWs+/EaQ6eLPZauPKWfd/u5forLuWdf69hwdPPcv9jTzQqXEHd1jBsqBCrmeSwoGYXrkB6dtUmNjaW7t2789JLLwGwatUqkpOTSU5OrvK80aNH8+KLL1JYWEhJSQnPP/8848aNq3HbTqeT0NDQyi+bzeatH0MEsIJSnc0/6SzNdvevevuAjknBVYmKO7sork/RuHTafXyYc4Txa7by7Nb9jF+zlY9yj9BvqtytKmQEq050XedUaQVHikoN04X9bG+uXsXsO26lZatW/Ovt9bTrcIHHtu2N0bYIu4WWwbZm291YenbV7h//+AcTJ05k/vz5hIaG8uKLLwIwbNgw5s6dS69evRg0aBBjxoyhS5cuAIwbN44rrrjCn2WLAHaqTGfncfjmuE5uobt3VbtQuDZJkRYGelkZDsevp80zrSQ+XjqPLVuyiUlOY8yjiw3RSkJ4nwSsGrh0nYKSco4WlXpk8ro/lJaU8Oicv/DPpc8ybNR13P/YkwQ7nf4u67wU0CLYRmQj7hRsCubMmVNlDpb07DpXWlraOW0ZAN54440qj2fPns3s2bN9VZYwmOMlOrtPQHa+Tu4p9zEo5UyoCqdKp/WisnO/vyGtJH7b2iE6KY3+U+8N+HYUov48GrBOl5VXufvNqMoqXJwoKeN4cRnlBrsM+FsHc3P489Qbyf7ma+595DHGTZoa0P82mlLEh9gJsUrul55dQniHruv8VAS7TriD1eFi90hV2xAYnqi4IBwcZu8dJ8+0dhiQFM2U3klsyD3IqzNHM+avr0nIamI8eiY7dKqE70+cJsxmIcxuNtTcGZeuc6q0nBPFZRSWVWDcWOW28Z23uPuWaYSEhvHSf96hc3oPf5dUI4umaB3qwN6M7hSszahRo2RCuxANUN0IUWzfa/n6OOw4pnO0BOwmSA2DgXGKdqFg9VFbhY+WzGNAkru1g1KKm7onMX7NVj5eOk8CVgCzaBBlrV8y8PhQQUmFi8OnS/j5dAnBFhOhNgshVjMmLfBGTly6zsnScgpKyjlVWu7RVgv+Ul5ezlMPz2PJk4sYNPRKHnrqGcIjImv/Rj9ymE0khNoNFciFEIHpnBGiHPcIEX94BWuvkXQIh6EJirah/umwfiQnmym9q2/tIAKTWcG4FI2Icj8HrDN04FRZBafKKlBAkMWE02rGaTH5tZ9RaYWLwrIKTpWWc+R0BTbdv71KPOnnQ4e46+bJbP10M3fOeZCJ028N2IWazwj9ZcFmLYAvXQohjGPjc782//ztCNGu9+Zz8+RRWPz8y350Uhobcg9W29pBBB6TgrEpGimhimPH6ve9PpnsogOFZRUUllXwE+7LQQ6LiSCzCbvZhN2seeUEW+HSKalwUVxeQVF5BafLKqq0VvDFtPV3163l6UcXsP+7vSSntGf6XXd7pb/Vpx9t4q5pk9E0jWWr/0PPi/t6/DU8LSbISrTDGtDzwoQQga+wTGf7Mdh+TOfnnD1MO0/zT3+HK4D+U+/l1ZmjGb9mK4OTolifc5SPco8w5tHF/i5NnEVTMLqtRvuwhu03fplNXObSKStxX5oD910bVpOG1aRhMSmsmoZZU5g1DZNyT37WKj8s7m24dB2X7v5vhUunXNcpd+mUVrgoc+mUlFdQ7tL9OpfKFwshu1wuljyxiKcWzOPCfgP46+KlRMfGemTb3qIpiHPam9WyN0IIz6rQdb7Nh21Hdfbku88NaWFQ1vr8zT8DgbR2MAZNQUayRofwhofygLhdS8c9d8uorRDOp6aFkD0RsE4cO8rd06fx0QfvMe2O/+GWmf+LyRTYk8QtmiIh1IFDJrMLIeqp3KWz/yTsKXD3qyosh5YO95yqzpEQZFbsuvm+gB8hakhrB5D2Dr4SYoHfJWh0imzciGdABKymypsLIX/1+Rb+PHUixUWneWbFSgZcNqTR2/Q2mcwuhKivwjKd7HzYk6/z/Ukoc0GYFTpFQHqUIi6o6kmwqY4QSXsH77FokByiaBsCKaGKWIdnLiVLwPIibyyErOs6/1zyD56cP5fO3Xuw8LkXiItP8ES5XuU0K5LCHDKZXQhRJyfLdD4+pPPFEajQobUTLmmpSA2DGPu5i33/VkNHiAKZtHfwLLOCtHBFt0hFSihe6XQgAcuLPL0Q8smCfO790y2895+1TJx+K7ffez8WS2DPY1JAbLANVVQm4UqIZqyul7dOlup89JM7WFk06N9ScWGM+/JfcybtHTyjdbAiPUrRKQLsXt6n5FqNF3lyIeSdX21j9GUD+PTDjSx8bhl3PfBQwIcrk3I3D41q5sveCNHcnbm8lVJykFm9k0gpcV/e2vXB6srnFFfovJXn4olv3HcEDmipuK2zYmCc8km42vXBap4b35OH+zl5bnzPKrUFAnd7h6NV5vQG0uT9QGZW0CNacUtHjSkdNHrGKK+HK5ARLK9r7ELIuq7zr2VLeOS+e0jt2IklK9cSFdvCgxV6h92kkRDqwGqSDC9Ec1fT5a3ki69g1wmdNw/olFS4g9VFsVXXAPQ2I8xvkvYO9ecwQa8YxUWxCqfF9yOghjz7vbtuLSMH9qV7QgwjB/bl3XVr/V2SV5wsyOfOqROZN+tORv9+Ii+te4eEpGR/l1WrUJuZ5PAgCVdCCMB9eWtQYtQ5l7d+3p/Nmjwzr36vE+eA6R3dI1a+DFdQNQBO65HMy9f0oH9iNB8vnefTOmpyZvL+Pkc8j2zJYZ8jnjGPrjT85H1P05R7wvqw1oo7umhcFq/5JVyBAUewfNFbKhDs/Gobf556I8ePHmXR0uUMHRH4HyIFxATZiA6SS4JCiF+dr3u5q0Uaeac1rmuj6Bhe88R1bzLK/KZAn7zvrzYSdhO0D1OkhrnXlfTmYt31Ybghhpp6SzUFuq7z8pJ/MH7Y5YSEhvHa+5sMEa7OzLeScCWEOFv/qffyYc4Rxq/ZyrNb93P9mq18lHuEhHF/YUpKKZ0ilF9XdJD5TY1Xl3l23tAtUnFXV42MNhpdIlXAhCswYMDyZm8pf8s/cZzbJ/2e+ffcxdgbJ/PyG++S2Katv8uqld2s0SY8CKfVcAOiQggfuODSkYx65FW2m+N56NMcPlfxDJr7GjeOGYUjAHoOnx0Ax/8SAPtNvc/fpRmGPy6zXtJSMbKN5pUWC55guDOiN3pLBYJtn33KXTdP4VRBQYPvNPSHMJuFOKdNWjAIIc7r23ydD1pcS+E913JJS0W/FmAOoJNiU21O6ku+vMyqKbiqtftuwEBWbcAqLi5m3Lhx7Ny5k6CgIFq2bMnixYtJTk72cXnn8nRvKX9zuVwsffJvPLVgHl169OLF19+gVetEf5dVKwW0CLYRKS0YhBDnkV+q806ezs4T0CYEJrRTRNkD86QY6PObAt355tl5+jKrVWvcAsy+dN5LhNOmTSM7O5tt27YxfPhwpk2b5su6zsuTvaX87fChH5l63TU8MX8uk2fczgsGCVcWTZEUFiThSghRrXKXzqYfdf7+jU5uIYxMVvw+gMOVaDxfXGZ1WmBiqjHCFZxnBMtutzNs2LDKx3369OHxxx/3WVG1aWxvqUCw/u03ufdPf8RisbJk1Vr6DBjo75LqJNhiIj7EjlnWExRCVGNPvs5bB3TyS6FPLFwSp7D5uO2C8D1vX2aNscMN7TTCbcbZl+o0B+vJJ5/k6qtrHyEqLSulpKTk142bTJjMhpvm5VVFp0/z6P1/4ZVlSxk09ErmPfE0EVFR/i6rTqIcVmKDrH6920cIEZiKK3TW5ujsOgFtQ2B8O0W0jFg1K966zNomRDG2rW+6r3tSreln/vz57N27l8WLa+8W+8KyF4j4TZfxfv360a9fv8ZV6EXFRcU+fb3dO7Zz723T+TEvj7sf+isZEzJRSlFUVFSv7fi6bk1BtFXDUlLG8ZLCBm3j+PHjHq7K+4xYMxiz7sjISH+XIOqoul5HURdfyyvf6xSW4/eeVqJpSY9SXJ2oAvZOwZpUBqzly5ezaNEiAG677TYmTZrEwoULycrK4r333iMoKKjWjU2cNJELunT7deMGGMFyOBxef42KigqWPvk3/t9f59OuQ0dee38TKamNm/jni7rB3YIhIcQzS94Y8SRqxJrBuHWLwHa+JWVMf3yF6L4jmZamiJRRK+Ehg1spBsYZdzpKZfrJzMwkMzOz8i8WLVrEihUreO+99wgPD6/TxqwWKzabzfNVGljO99/xvzP+wPYvPmPKrXcwfeY9WK3GmBwebrPQUlowCCF+Ud2agtev2cq2t+YzZdIoLAYcZRCBx6LBiCR341Ajq3Z4KS8vjzvvvJO2bdsyePBgAGw2G59++qlPizMyl8vFv5YtYdHc2UTFxPLC62/Ss8/F/i6rTjQFLYPthNst/i5FCBFAqut1dGlSFJ9tyZZwJTwiwgZj22q0DDL+/lRtwEpISDinW3pT8O66tTz96AL2f7eX5JT2TL3tz1w16jqPv87B3Bxm33Ern2zawNhJU7hz9oMEO50efx1vsJo0EkLs2M0B0F5ZCBFQIhNT2ZDzg9d7HYnmKSVUcV2bwFrupjECe4KUB1W3SPTMmydjtVo91vLB5XKx4vnn+NuD9xMWEcFzr62h76BLPbJtXwi1molz2g05mVAI4T0lFTofHdI5fPlfOPTMWMav2crgpCjW5xzlo9wjjHm09pughKhJvxaKy+JVk5qSYtzZY/Xk7UWiv9uTzY0jrmT+PXcxYuw4Xv/wE8OEqzNd2RNCHRKuhPCzXR+s5rnxPXm4n5Pnxvf0+mK5tdmbr/PUNzqfHIa+V4xk5IJX2eeI55EtOexzxDPm0ZWypIxoME3B8PhyhiRoTSpcQTMawfLWItGlJSU898RjPPv4Y8S3TmTZmv/Qu9+ARm3TlyyaIj7EQZBFLgkK4W/nu0tvzF9f8/kyLhW6zvsHdTYfhvahcFWiIsyqIH4UXS8f5dNaRNOkKRiZrBGPy9+leEWzCVjeWCT64/Xv89Dd/8PBA7lMufV2br7jLmx2uyfK9Qnpyi5EYKnuLr3xa7by8dJ5Pg1YJ0p0Vu7T+fE0/C5e0SdW+loJz1LANb/cKXjsmL+r8Y5mc2adftfdlZcFgUYtEv1D3gH+POVGpo0ZSWxcK1at/5g/3XOfYcKVAmKCrCSGOiRcCRFAjuRkMygxqspxanBSFD/vz/ZZDbuO6/xjt7tp6OQ0xcUtlIQr4VEKdxuGblFNe79qNmfX6haJfvTZZfVaJPp0YSF/f+Qhru7biy8++S+PLF7CstXraJfWwYuVe5ZZU7QOdRATZJODphABJjopjQ25R6vMFfXVXXqHTuu8+r2LV/fptAmBmzso4oPlGCE8SwHDEzW6Rzf9favZXCKEcxeJrusSNRUVFax9ZQVPLZjHsaNHuPEPM5h2x50EO0O8VapXBFlMxDvtWDzQlV0I4Xn9p97LqzNH+/QuvR9P62z6UWd3PoRb4ZokRbdIuSQoPE/hnsvXM6Z57FvNKmDVl67rbHz3LR6f9wB7d+1k6IiR3HHf/bRObuPv0upFAZGyULMQAe+CS0cy5q+v8fHSeWzZkk1MchpjHl3s8bv0dF0nrxA++klnT767ueM1SYoukWCSY4TwAqsGGW000sKbz/4lAasauq7z8fr3+fsjD7Fj6xdc2G8A/3pnPV269/R3afVmUopWIXZCrPJPLYQRXHDpSK9NaC+t0Nl+DD4/ovNTEUTZ75DrdgAAH6xJREFU4NpfglVTu0VeBI4wK1yf0jS6s9eHnHV/48yI1bN/W8hXn39G+oW9WbpqLRcNGGjIkR+H2X2XoCcWahZCGNfhIp3Pj+h8dRTKXNA+DC5rpWgXKpcChXfFByuuT1E4Lc1vP5OABZSVlfH22tU8/9QTZH+zg+69+7D4X6vof+nlhj34RNottAiWiexCNFflLp1dJ9yjVbmnINgMF8VCz+hf+lkJ4WWdIxTXJivMzbSBdbMOWCfz83ll2RL+79lnOHQwj76DLuWFNW/Qq28/wwYTTSninDbCbLJQsxDN0YkSnU8Pm9iR7261kOyE69ooOoQhKzUInzApuCxecXFs827x4beAdfbCy9PvuttjawLWZu+unfxz6bOsfXUF5eXlXDVqNDf+cQZpnTr75PW9xWbSSAhxYDPLJUEhmosKXefjN1fzxbJ5FB7cAy1SUcP/Qo/fjaRXtCLG0XxPcML3Im3uyezS4sNPAau6hZdvnzSBx5e95LWQVVxUxDv/XsOrLy7jyy2fEB3bghv/eCvXT76JmBYtvPKavhRus9DSaZOJqkI0Ay7dffff9mM6ezespvzpsQxIjGbQRUlsyP2BD58ZS3Kb14hp7dvldUTz1jlCcXWSwmaS8xD4KWDVtPCyJwOWrut8vW0rq1e8xBurVnKyIJ8+lwxi0ZIXGXzlVVRUVOBwODz2ev6gKWgRbCfCLpcEhWjqTpXpbD0CXxzRKSiDuCAIfushugbA8jqi+bJoMKx182geWh9+CVjeWnj5jB8O5LJu5av8e+UrfL8nmxZxrbh+8k1cO248SSntKp9X10ajgcpq0kgIsWM3y0LNQjRleYU6nxzW2XXc/UtVl0i4MEYRF6R4+OAeBvVOOmd5nS1bfLe8jmi+2oa4R60ibBKuzuaXgOWNhZd/PnSIt9eu5s01WWz77FMcQUFcNmw4sx58mIsHDsZkalohJNRqJs5pl0mrotkpLi5m3Lhx7Ny5k6CgIFq2bMnixYtJTk4+57kbNmxg2LBhpKb+emzZvHmzYUauDxbqbPhR59sCd8+qIQmK9Eiwm3/93LuX1znITd2TKq8M+Gp5HdF8OUwwtLVGehNfT7Ax/BKwpt91d5U5WA1dePlgbg7vvfFv3v33WrZ99ikms5n+l17OI888x+ArriLY6fTST+A/CoiyaiSEGuMEIYQ3TJs2jSuvvBKlFH//+9+ZNm0a77zzTrXP7dixI59//rmPK2ycH34JVnsLINoOGcmKThHV96zyx/I6onnrHKG4srUiuBn2tqoPvwSsMwsvP7NwAfu+3Uubdu2Zftc9tS687HK5+OarL9nw9pt88NYb7PnmayxWK30HXcq8J59m8BXDCAuP8NFP4XtWTSM+xE7RyTJ/lyKE39jtdoYNG1b5uE+fPjz++ON+rMgzSip0dp+AHcd1vvtlxGrUL8GqpptXqlte59qHnvD48jqiedMUtA9V9I5VpIRKsKoLv7VpOHvh5fM5cewomzdu4OP177PpvXc4+vNhQsPCGXD5EG6+/X/of9nlOENCfVCxf4VYzbT65ZKgsWeOCeFZTz75JFdfff5fzrKzs+nRowcmk4lJkyYxffr0Grd36tQpCgoKKh/bbDZsNpvH6v2t0gr33YDfHHePVlXo0DoYRiYpOtdj+Zqzl9cx+vxSETjCrdAjWtE9ShEiDWrrJeAajZYUF7Pt8y18snEDmzet5+svt6LrOu06XMCIsdczaMgVpPe+CLM54Er3CgXEBtuIclj9XYoQAWf+/Pns3buXxYurvxzWo0cP8vLyCAsLIy8vj2HDhhEdHc2YMWPOu82BAwdWeTxz5kxmzZpFfr5GUVF5o2vWdThYpNh+wkR2gUaZrmhpdzEgxkWH0ApCf7khuKS44a9RVNSIb/YxqdU7GlurRYMR8WWkOHWUgrJTcMxDtZ3t+PHjXtqyZ9VWZ2RkZJXHfk8ppwsL2b71c77Y/DGf/fdjvvp8C6UlJURGR3PRgIGMyZxM38GX0rJVvL9L9TmLpogPcRBkaVoT9IWor+XLl7No0SIAbrvtNiZNmsTChQvJysrivffeIygoqNrvCw39dXQ7ISGB66+/ng8//LDGgLVx40bS09MrH58ZwQrTdRyOwgZPkD9V5l4L8MujOkdL3CMD/Vu6F1qOsHn+UGyUifwgtXpLQ2tVwNgUjQ7hvhuxOjucBKr61OnTgKXrOgf272PHl1/w1Wdb2Pb5Fnbv2E5FRQVhERH0urgff77vAS7sN4DUjp3QtObbkTzkl7sEm+saTkL8VmZmJpmZmZWPFy1axIoVK3jvvfcIDw8/7/f9+OOPtGjRAk3TOHnyJOvWrWPKlCk1vpbT6awSzBorv1Rn4486246657F0DIerEhXJTlloWQSmy+KVT8NVU+W1gFVRUcGB/d+ze8cOdu3Yzs7t2/jmqy/J/2WILaltCt16Xch1E26kR5++tG2f2qwD1RlnLglG2i1y8BWiGnl5edx55520bduWwYMHA+5Rpk8//RSAqVOnMmLECEaMGMGqVat45plnMJvNlJeXM3r0aCZNmuSTOr98N4v1z7qXrzG1TKXrhHsZOmIkDrN8rkXg6hap6N9SzsWe4NGAteaVf/Lykn+wd/dOvsveTfEvEy1bxLWiQ5euTJj2R7qk96Bz955EREV58qWbBKum0SrELpcEhahBQkLCOY2Kf2vJkiWVf54xYwYzZszwRVmVisp11q3JYueCs5avmT+GtPDXpLu6CFitg91NQ4VneDRgrXppOe06dCD1gk5clTGa9h060qFzVyKjoz35Mk3Sb+8SFEIYS3GFTvYJ992A350E7eWHGJAYzcvXyvI1whjCrTAuRcm0FA/yaMBa/u836dStuyc32eRpCmKDbETKXYJCGM63BTqf/+zutF6hQ2IwDI1XvH9Ylq8RxuEwwfUpmjQO9TCPBiyF/OPUh83kbhwqawkKYTzbj+qsztFpFQSXxys6hkPoL32CvpLla4RBhFvhhnYaMQ45f3ua39s0NFfhNgstnbY6NxIUQgSOHcd01uTodI+CqxPVOTekyPI1wgjigmB8iiYNRL1EbhXwMZNSxIfYaRVil3AlhAF9fUxn9X6dbpHVhyv4dfmafY54HtmSwz5HPGMeXSnL14iA0T5MMSlVwpU3yQiWDwWZTbQKsWM1Sa4VwoiyCzTWHtTpEglXJ1Ufrs44e/kaIQJFj2jF8EQlv+R7mQQsH1BAdJCVaIdVelsJYVAfHHTx74NmOkfANUlychLGdElLxaXx8ku+L0jA8jKrSaOVU3pbCWFkLl3n/q0uUkNdXJtslnAlDOnSVopL4iRc+YoELC8Kt1toGSwT2YUwOk0p/nWpifU5JWjK4u9yhKi33yUo+raQcOVLErC8wKwp4oLthHhhAVchhH+0ClZID0ZhNAq4orXiolgJV74mCcDDQq1mWjptmGVdRSGEEH6kgOGJGj1j5DcDf5CA5SEmpWjptBFmk8sHQggh/EsBI5I0ukdLuPIXCVgeEGI1EyejVkIIIQKAAn4XVy7hys8kYDWCWVO0CJZRKyGEEIHjytaKdmaXv8to9iRgNVCYzUyLYLusPC6EECJgDE1Q9I7VOHbM35UICVj1ZDVptAy24bTKWyeEECJwXNZKcbG0YggYkhLqSAFRDivRQVbpayWEECKgDG6lGCBNRANKnf41HnjgAZRSfP31196uJyAFW0y0DQ8mVpqGCiGECCAKuCpRMVDCVcCp9V9k69atfPLJJyQmJta6sTtvmsS769Z6pLBAYNEU8SF2ksKCsJll5xVCCBE4LBqMTdG4MEbOT4Goxn+VkpISbrnlFp5++uk6LVKcu+87bp80wfAhS1MQ7bCSEhEsdwgKIYQIOEFmyGyv0SFcrqoEqhoD1uzZs5kwYQJt2rSp8waVUjz96MOUlJRQUV7e6AJ9LdRmlsuBQgghAlaEDaakabR2yjkqkJ13kvvmzZv57LPPWLBgQb02qOs632bvZsGCBfTr149+/fo1ukhvKS4qrvyzTVNEWjXsZWWcyi/yY1W1O378uL9LqDep2XeMWHdkZKS/SxDCENqGKEa1UTgtEq4CXZWAtXz5chYtWgTA2LFj2b17d+XoVV5eHkOHDmXJkiVceeWV592gUop2aR24++67MZtMmMyBfaNiaHAQscE2Qq3mOl0GDRRGPCFJzb5j1LoDVVZWFvfOmct3e7OJTkqj/9R7ueDSkf4uSzQjZgWXxysuilWGOlc1Z1XST2ZmJpmZmZWP77nnnso/Jycns27dOjp37lzjBnVd55aZ/4vNZvNwqZ5lUu4Rq6SIYLkUKIQ4r6ysLDIyMrgkKZrreiexIfcgr84czZi/viYhS/hEqyDFyGRFjEPOVUbi0VsPktqm8MQLL3P5VVd7crMepYBIu8U9gd2iSbgSQtTooQfncklSDC9d04NpPZJ5+Zoe9E+M5uOl8/xdmmjiNAWXtFRMSZNwZUR1vn63f//+Wp+z8NlldOyW3ph6vCrUZiY2yIbVJLe0CiHqZnd2Nnf2Sqq8LKOUYnBSFFu2ZPu5MtGUtXTAiCSNVsESrIyqWSSNIIuJNmFBJIQ4JFwJIeqlQ1oaG3OPoes64J4GsT7nKDHJaX6uTDRFZuVe8mbaBRKujC6wZ6A3ks2kERtsI0TWDRRCNNBf7ptNRkYGN7y+lUGJUazPOcpHuUcY8+hif5cmmpgkp+LqJEW0XYJVU9Akh3PMmiLOaaNteJCEKyFEo4waNYpVq1ZxNDSBR7bksM8Rz5hHV3LB4Gv9XZpoIhwmGJ6omJgq4aopaVLpQ1MQaZcFmYUQnjVq1CjaDhzJil2FOBwOf5cjmggF9IhWXBavCDLLOaupaRIBSwGhNguxQVYsMsdKCCFEgIsPVgxrrYiXeVZNluEDVpDFRItgGw6zyd+lCCGEEDWKskG/lhrdo5CGoU2cYQOWVdOIDbYSKosxCyGECGAKaB+m6B2jSAmVYNVcGC5gaQqiHFaiHDLPSgghROCyatAzWnFhjCJSJq83O4YKWKFWM7HB0ihUCCFE4NIU9IhSDGolizI3Z4YIWDaTRotgG05puSCEECKAdQhXXB4v7RZEgAcsuRwohBDCCNqGuEesEp1yrhJuARuwgi0m4px2uRwohBAiIFk1SI9S9I6VEStxroALWGZN0SLIRphd7g4UQggReGId7snr6VEKm0mClaheQAWsMJuZFsF2zJrssEIIIQKDAloFK+JDKuidqMlolaiTgAhYFk0R57TLJHYhhBABQQFJIYqO4e6J66FWxbFjFdJuQdSZ3xNNhN1CbJANk4xaCSGE8CMFJDoVnSKgY4S0WBCN47eAZdEUrZx2gmXUSghRT8nJydjtdux2OwD33HMPY8eOrfa58+bNY9myZQCMHz+eBx980Gd1isBnUpAcougQDh3CFCFWCVXCM/ySbsJtFloEy6iVEKLhVq5cSefOnWt8zqZNm1ixYgXbt2/HbDbTr18/+vfvz9ChQ31UpQhEFg3ahyo6hCtSw8BulnOR8DyfBiyzpogLthNik1ErIYT3vfLKK0ycOJHg4GAAJk+ezIoVKyRgNUOagpRQRZcI92iVVe7+E17ms6QTYjUT57Rh1qSvlRCi8W644QZcLhcXXXQRDz/8MDExMec8Jzc3l4EDB1Y+Tk5OZuXKlTVu99SpUxQUFFQ+ttls2Gw2zxUufEYBrZ2KzhHQKUIRLHOqhA95PWBpCmKDbEQ6rN5+KSFEM7Fp0yYSExMpKyvj3nvv5cYbb+SNN96o9rnqN6tA6Lpe67Z/G8gAZs6cyaxZs8jP1ygqKm9c4T5SVFTs7xLqzNO1mhS0DnKRGuqifYgL5y9nuZKTUNLIbR8/frzR9fmK1Op5tdUZGRlZ5bFXA5bdrBHvtGMzm7z5MkKIJm758uUsWrQIgNtuu41JkyYBYLFYuP3220lNTa32+xITE9m/f3/l45ycHBITE2t8rY0bN5Kenl75+MwIVpiu43AU4nA4GvnT+IZR6gTP1Joc4m782cHLc6rOPokGMqnV8+pTp9cCVqTDQosgW5XfHoUQoiEyMzPJzMwEoLCwkBMnThAeHg7AihUr6N69e7XfN3r0aGbMmMH06dMxm808//zzzJs3r8bXcjqdhIaGevYHEF4RbIZuUYoe0bJUjQg8Hg9YZqWIC7ETIu0XhBBe8NNPP5GRkUFFRQW6rtO2bVuWL19e+ffDhg1j7ty59OrVi0GDBjFmzBi6dOkCwLhx47jiiiv8VbrwALvJPVm9Y7h7srrcjS4ClUdTkMOs0SYiCItMZBdCeEnbtm358ssvz/v3Z8/Fmj17NrNnz/Z2WcKLYuyQGqZoH6ZIdIImV0aEAXg0YLV02iVcCSGEaDQFpIUrBrRUxAdLoBLG49GAJfOthBBCNIam3C0VBrRUxDrknCKMSyZKCSGECAhdIhWD45QsqCyaBAlYQggh/CrMCsMTNdqHSbASTYcELCGEEH6hgJ6RFVybpsnSNaLJkYAlhBDC52IdMCJRI6i0QsKVaJIkYAkhhPAZs4JL4hT9WihMmuLYMX9XJIR3SMASQgjhE0lOxdVJ0nVdNA8SsIQQQniV3QSXxyt6Ritp5yOaDQlYQgghvKZrpGJIvCLEKsFKNC8SsIQQQnhcqyDFla0VrZ0SrETzJAFLCCGExwSb4bJ4je5RsrqHaN4kYAkhhPCIzhGK4YkKu1mClRASsIQQQjSKpmBIvOLiFpq/SxEiYEjAEkII0WBOC1zXRiM5REathPgtCVhCCCEaJCFYMaatIlTuEBTiHBKwhBBC1ItZQZ8WisFx7m7sQohz1XjBvKSkhBkzZtC+fXs6derEhAkTfFWXEEKIAKOAbpGKWztrXB6vSbgSogY1jmDdfffdaJrGnj17UErx448/+qouIYQQAaRtiGJIgiIuSEKVEHVx3oBVWFjIsmXLyMvLq+xlEhcX57PChBBC+F+QGUYkaXQIl2AlRH2c9xLhd999R1RUFPPmzaNXr14MGDCA999/v8aNnTp1ioKCgsqvkpISjxcshBDCNxKdij9cIOFKiIY47whWWVkZ33//PR07dmTBggV89dVXXH755ezcuZOYmJhqv2fgwIFVHs+cOZNZs2Z5tmIPOn78uL9LaBAj1i01+44R646MjPR3CeI3FNC/pWJwK4Um3diFaJAqAWv58uUsWrQIgBtuuAFN07jhhhsA6NatG23atOGbb75h0KBB1W5s48aNpKenVz622WzYbDYvle4ZRj2wG7Fuqdl3jFq38L9gM4xM1mgXJsFKiMaoErAyMzPJzMysfPzuu+/y9ttvM2zYMHJycti3bx9paWnn3ZjT6SQ0NNR71QohhPAKs4Ie0YpL4hROi4QrIRqrxrsIFy9ezOTJk5k1axYmk4lnn31WJroLIUQTYtGgZ7SiXwtFiDQMFcJjagxYbdu2ZcOGDT4qRQghhK9YzwSrljJiJYQ3SCd3IYRoRmLscGGMoluUwmaSYCWEt0jAEkKIJk5T0CFM0TtWyaLMQviIBCwhhGiiImzQI0rRPVouAwrhaxKwhBCiCdEUpIa4GNRGo20IlStxCCF8SwKWEEIYnMME7cIU7UMV7cKguKCcyFAJVkL4kwQsIYQwoFALdItSpIYp4oOp0nG92I91CSHcJGAJIYRBKCAlVNErRpEahixjI0QAk4AlhBABzmGCnjGKntGKCJuEKiGMQAKWEEIEKIcJ+rRQ9ImVnlVCGI0ELCGECDB2E1zcQnFRjMJulmAlhBFp/i5ACCGEW1wQDE1Q3N5ZY2Cc5vVwlZWVRc/u6QQHOejZPZ2srCyvvp4QzYkELCGE8KNIGwyMU8zopHHzBSYubuH9YAXucJWRkYHz+A/c2SsJ5/EfyMjIkJAlhIfIJUIhhPAxmwk6Ryi6RykSnP65BPjQg3O5JCmGl67pjlKKm7onccPrW5k/70FGjRrll5qEaEokYAkhhI8kh7hDVccIsGj+nVu1OzubO3slVXZ6V0oxKDGKxz7f7de6hGgqJGAJIYQXKaBzpGJQnCLKHjgT1jukpbEx9wdu6u4OWbqusyH3KBd06ODv0oRoEiRgCSGEl7QPU1zWStEyKHCC1Rl/uW82GRkZ3PD6VgYlRrEh9ygf5hwh62/P+rs0IZoEmeQuhDCUEydOkJ6eXvmVmpqK2Wzm2LFj5zx3w4YNBAUFVXl+UVGR12tsHayYlKpxQzstIMMVwKhRo1i1ahWnIxN47PMcTkcmkJWVxciRI/1dmhBNgoxgCSEMJTw8nG3btlU+XrhwIRs3biQyMrLa53fs2JHPP//c63VF2aBThKJjRGCOWFVn1KhRMqFdCC+RgCWEMLRly5bx0EMP+eW1g8zQM1rRyUChSgjhGxKwhBCGtXnzZo4ePcrw4cPP+5zs7Gx69OiByWRi0qRJTJ8+vcZtnjp1ioKCgsrHNpsNm812zvOCzHBje40WEqyEENWQgCWEMKznn3+ezMxMzObqD2U9evQgLy+PsLAw8vLyGDZsGNHR0YwZM+a82xw4cGCVxzNnzmTWrFnk52sUFZUD7qVshrcsw1Ksc6zYcz+Ppxw/ftzfJdSZ1OodUqvn1Vbn2dMUJGDVQVZWFg888AB79uwhNTWVOXPmyLwFIXxo+fLlLFq0CIDbbruNSZMmUVhYyCuvvMKWLVvO+32hoaGVf05ISOD666/nww8/rDFgbdy4kfT09MrHZ0awwnQdh6OQCKeDzPYarYIDe+TqfHPSApHU6h1Sq+fVp04JWLU4s5zEmT4xO3bsICMjg1WrVknIEsJHMjMzyczMrPL/XnvtNbp27UqHGvo2/fjjj7Ro0QJN0zh58iTr1q1jypQpNb6W0+msEsx+y6rBhHaBH66EEP4nbRpq8cADD1SGKwBd11FKMXfuXD9XJkTztnTp0mrD0tSpU1m7di0Aq1atokuXLnTr1o0+ffowZMgQJk2a1KDXs2owJrHMb0vbCCGMRUawarFnz57KcHWGrutkZ2f7qSIhBMCHH35Y7f9fsmRJ5Z9nzJjBjBkzPPJ6HSMUx846FgghxPnICFYtUlNTK9fqOkMpRVpamp8qEkIIIUSgk4BVizlz5lReFgQqLxfOmTPHz5UJIYQQIlBJwKrFmeUkunbtit1up2vXrrKchBBCCCFqJHOw6kCWkxBCCCFEfXhkBKu0tLTKf42gpKSERx55hJKSEn+XUi9GrFtq9h0j1l1SUsL9998f8DUb5b01Sp0gtXqL1Op5DalT6WffItcAmzZtYuDAgWzcuJFLLrmksZvziYKCAsLCwsjPzz9vz5tAZMS6pWbfMWLdgVLz1q1b6dmzJ1988QU9evQ45+8Dpc7aGKVOkFq9RWr1vIbUKXOwhBBCCCE8TAKWEEIIIYSHeWSSe3Gxe7XTPXv24HQ6PbFJrzt16hQA27ZtM0zNYMy6pWbfMWLdZ2o+ffq0Xy8RFBUVAbBr165q/94o761R6gSp1VukVs+ra50dOnQgKCjI/UD3gBdeeEEH5Eu+5Eu+Gvy1efNmTxyOGuyll17y+3sgX/IlX8b++uKLLyqPKR6Z5H7kyBHefvttkpOTcTgcjd2cEKIZqvKbnx/IcUwI0Vi/PY55JGAJIYQQQohfySR3IYQQQggPk4AlhBBCCOFhHglYv/vd7+jatSvp6ekMGDCAbdu2eWKzXlNcXMy1115Lamoq6enpXHHFFezfv9/fZdXJn/70J5KTk1FK8fXXX/u7nFrt3buXvn37kpqaSu/evdm5c6e/S6qV0d5jMO4+baRjhxFqNdp+YJTPmlGOY0Z5P8F4+2qDPv+euPvm+PHjlX9evXq13r17d09s1muKior0//znP7rL5dJ1XdefeuopfciQIX6uqm42btyoHzhwQE9KStJ37Njh73JqNXjwYH3ZsmW6ruv6a6+9pvfp08e/BdWB0d5jXTfuPm2kY4cRajXafmCUz5pRjmNGeT913Xj7akM+/x4ZwQoPD6/8c35+PpoW2Fce7XY7w4YNQykFQJ8+ffj+++/9XFXdXHLJJSQkJPi7jDo5fPgwW7duZcKECQBkZGSwb9++gP4tBYz1Hp9h1H3aSMcOI9RqtP3ACJ81Ix3HjPB+nmG0fbUhn3+PNBoFyMzMZP369QC89dZbntqsTzz55JNcffXV/i6jyTlw4ACtWrXCbHbvZkopEhMTyc3NJTk52b/FNXFG2qeNdOwwUq1grP0gUMlxzDeMsK/W9/PvsYC1fPlyAF588UXuuusu3njjDU9t2qvmz5/P3r17Wbx4sb9LaZLO/HZyhi5dQbzOaPu0kY4dRqrVaPtBIJPjmHcZZV+t7+e/QWPcy5cvJz09nfT0dJYtW1bl72688UbWr1/P0aNHG7Jpr6mu5oULF5KVlcWbb77p1waHNanpvQ50rVu3Ji8vj/LycsB9UDpw4ACJiYl+rqzpMsI+fT6BduwwynHOSMc2Ix7P5DjmXYG6r9akzp//xk78ys/P1w8ePFj5OCsrS4+Pj6+cuBaoHnvsMb1Hjx76sWPH/F1KgxhhEqOu6/rAgQOrTA696KKL/FtQPRjlPT7DaPu0kY4dRqrVaPuBrgf+Z81ox7FAfz/PMMq+2tDPf6M7uR84cICMjAyKiorQNI2YmBgWLlxIenp6YzbrVXl5ebRu3Zq2bdsSEhICgM1m49NPP/VzZbW75ZZbeP311zl06BDR0dE4nU6+/fZbf5d1XtnZ2UycOJGjR48SGhrKiy++SKdOnfxdVo2M9h6DMfdpIx07jFKr0fYDo3zWjHIcM8r7CcbaVxv6+ZelcoQQQgghPCzw7jMWQgghhDA4CVhCCCGEEB4mAUsIIYQQwsP+P6xijLBtQ6orAAAAAElFTkSuQmCC" }, "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 }