{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "*This example is a Jupyter notebook. You can download it or run it interactively on mybinder.org.*\n", "\n", "# Linear regression\n", "\n", "## Data\n", "\n", "The true parameters of the linear regression:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "k = 2 # slope\n", "c = 5 # bias\n", "s = 2 # noise standard deviation" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true, "nbsphinx": "hidden" }, "outputs": [], "source": [ "# This cell content is hidden from Sphinx-generated documentation\n", "%matplotlib inline\n", "np.random.seed(42)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generate data:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x = np.arange(10)\n", "y = k*x + c + s*np.random.randn(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model\n", "\n", "The regressors, that is, the input data:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "X = np.vstack([x, np.ones(len(x))]).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we added a column of ones to the regressor matrix for the bias term. We model the slope and the bias term in the same node so we do not factorize between them:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from bayespy.nodes import GaussianARD\n", "B = GaussianARD(0, 1e-6, shape=(2,))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first element is the slope which multiplies x and the second element is the bias term which multiplies the constant ones. Now we compute the dot product of X and B:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from bayespy.nodes import SumMultiply\n", "F = SumMultiply('i,i', B, X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The noise parameter:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from bayespy.nodes import Gamma\n", "tau = Gamma(1e-3, 1e-3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The noisy observations:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Y = GaussianARD(F, tau)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inference\n", "\n", "Observe the data:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Y.observe(y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Construct the variational Bayesian (VB) inference engine by giving all stochastic nodes:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from bayespy.inference import VB\n", "Q = VB(Y, B, tau)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Iterate until convergence:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 1: loglike=-4.080970e+01 (0.003 seconds)\n", "Iteration 2: loglike=-4.057132e+01 (0.003 seconds)\n", "Iteration 3: loglike=-4.056572e+01 (0.003 seconds)\n", "Iteration 4: loglike=-4.056551e+01 (0.003 seconds)\n", "Converged at iteration 4.\n" ] } ], "source": [ "Q.update(repeat=1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Results\n", "\n", "Create a simple predictive model for new inputs:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "xh = np.linspace(-5, 15, 100)\n", "Xh = np.vstack([xh, np.ones(len(xh))]).T\n", "Fh = SumMultiply('i,i', B, Xh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we use the learned node B but create a new regressor array for predictions. Plot the predictive distribution of noiseless function values:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdUXPmd5/33BYQQykJZQkhCrUCsAokshDoppyaM7e62\ne7zt8azDzK7nzLTdveMOe/o8Ds/xrH3mmWc9O+PdmTP2cVNIQqDsbomcoUBkCSRAAokgQIgc6u4f\noBJFDkX+vs7hWFTd+7u3sPrL1ffe3++jqKqKEEKIuc9ipk9ACCGEeUhBF0KIeUIKuhBCzBNS0IUQ\nYp6Qgi6EEPOEFHQhhJgnrKb6AIqiyHORQggxAaqqKuPZflqu0FVVlS8zfX388cczfg7z5Ut+lvLz\nnG1ftbW1RERE8Nvf/nZCtXbKr9CFEEKMzGAwkJWVRXZ2NgaDYcLjSEEXQogZVFtbS2xsLPX19ZMe\nSwr6HBMUFDTTpzBvyM/SvOTnOT49PT1kZmaSk5ODqprnVqNiroGGPYCiqFN9DCGEmEuqq6uJjY2l\nsbFx2G2++93voo7zpqhcoQshxDTp7u4mPT2dvLw8s12V9ycFXQghpsHjx4+JjY2lqalpyo4hBV0I\nIaZQV1cXaWlp5OfnT/mxpKALIcQUqaysJC4ujufPn0/L8aSgCyGEmXV2dpKSkkJRUdG0HlcKuhBC\nmNHDhw+Ji4ujpaVl2o8tBV0IIcygo6OD5ORk7t69O2PnIAVdCCEmqby8nPj4eFpbW2f0PKSgCyHE\nBLW3t5OYmEhpaelMnwogBV0IISbk/v37JCYm0tbWNtOnYiQFXQghxqGtrY2EhAQePHgwJeMbDAbS\n09MntK8UdCGEGKN79+6RlJRER0fHlIxfWFiITqfD1tZ2QvtLQRdCiFG0tLQQHx9PRUXFlIz/5MkT\nIiIiePz4McHBwWi1Wv7yL/9y3ONIQRdCiBEUFRWRkpJCZ2en2cduaWnh8uXLpKamcuTIEb773e+y\naNGiCY8ny+cKIcQQmpubiYuL49GjR2Yfu6enh9jYWK5evYqHhwenTp1i+fLlAGy8exevyEg2lZbK\n8rlCCDFZBQUFpKam0tXVZdZxVVUlNzeXiIgI7Ozs+NGPfsTmzZsBWF5Xh/f586wrKyPtrbdgAo9C\nSkEXQog+TU1NxMXFUVVVZfaxKysr0el0NDY2EhYWhouLCwCL2trQXL/Ovvh4cl9/ndt//uf0WFvD\nv/zLuI8hBV0IseCpqkpeXh7p6el0d3ebdeympiaioqLIzs7mxIkTBAYGYmlpiWIwsDspif1RUTxy\nciLipz+lddWqSR1LCroQYkFrbGwkNjaW6upqs47b1dXFV199xc2bN/Hx8eHTTz9l6dKlAGwqLsZX\np6Pb2pob3/seddu3m+WYUtCFEAuSqqrk5OSQmZlJT0+PWcfNysriwoULbNmyhQ8++IANGzYAsLy2\nFp/z51lbUUHqW29x39MTlHHd9xyRFHQhxIJTX19PbGwstbW1Zh23rKwMnU5He3s73/zmN9mzZw/Q\n2yfXXr3K3sREcl9/nVvf/nZvn9zMpKALIRYMg8GAXq9Hr9djMBjMNm5DQwORkZEUFhZy+vRp/Pz8\nsLCwQDEY2JOYyP6oKB46O6P7+GPaVq4023EHkoIuhFgQ6urqiImJob6+3mxjdnR0cPPmTW7fvk1g\nYCCfffYZNjY2QF+fPDycLhsbrv/gB9Q5OJjtuMORgi6EmNd6enrIzMzkzp07ZrsqNxgMpKWlERkZ\niaOjIx9++CFr164F+vrkERHYPXxIanAwDzw8zNonH4kUdCHEvFVdXU1sbCyNjY1mG7OkpASdTgfA\nd77zHRwdHYHePrnH1avsSUzkzhtvcOv99+mZxDT+iZCCLoSYd7q7u8nIyCA3NxdzLT1SV1fHhQsX\nuH//PufOnePAgQMv++QJCeyPjqbC1XXK++QjkYIuhJhXHj9+TGxsLE1NTWYZr62tjWvXrpGQkMBr\nr73Ge++9h3XfEyqbi4rw1enotLHh2g9/yNNt28xyzImSgi6EmBe6urpITU2loKDALOMZDAYSEhKI\njo7G2dmZn/70p6zqm8m5oroan/PnWVNZ2dsn12qnrU8+EinoQog579GjR8TFxdHc3GyW8foHTfzw\nhz9kW9+V96K2NjyuXGFPUhI5b77JV9/5zrT3yUciBV0IMWd1dnaSnJxMcXGxWcYbKmhCURQUg4G9\nCQl4RkdT7uY2pX3y5cuX4+7uPqF9ZT10IcScVF5eTnx8PK2trZMe60XQRFpaGkeOHOHw4cPGoInN\nhYX46nR0LF1KcmjolPXJV69ejUajwdHRsfdmq6KMez10KehCiDmlvb2dpKQkSkpKJj3WSEET/fvk\nKcHBlE1Rn3z9+vVoNBq2D1igSwq6EGJeKy0tJTExkfb29kmNMzBoIiQkhC1btgBg3dqK9upVY588\n77XXpqRPvmXLFrRarTHgYiAp6EKIeam1tZWEhATKysomPdaLoImGhgZCQkJwcXHp7ZP39PT2yS9f\nptzNjYwzZ2hbsWLyJz/A9u3b0Wq1rFu3bsTtpqSgK4qyGIgDrOm9iRqhquqniqJsB/4IrAaygHdV\nVR20MrwUdCHEZNy9e5fk5GQ6OjomNc5wQRMAWwoK8NXpaF+2jOSwMJ7a25vj1I0sLCzYtWsXGo3G\n+OjjaKbsCl1RFFtVVVsVRbEEEoG/Bn5Eb3HXKYry/wPZqqr+doh9paALIcatubmZ+Ph4Hj58OKlx\nurq6uHXrFjdu3MDHx4cTJ04YgyZWVlfjExHBqsePSQkJodzd3ax9cktLS/bu3Yu7uzvLli0b175T\n3nJRFMWW3qv17wGXgY2qqhoURfEBPlFV9egQ+0hBF0KMizlCmgcGTQQHBxuDJqxbW/G4fJlXUlN7\n++SvvorBjH1ya2trnJyccHV1ZcmSJRMaYyIFfUzPoSuKYgFkAo7A/weUAo2qqr5YuuwRMHRnXwgh\nxshcIc39gybeffdd9u7dC2DSJy/TaNB9/DHtZuyT29jY4OrqirOzs3F5gOk0poLeV7i1iqKsAC4C\n+4bazJwnJoRYOMwV0jxc0AS87JO3LV/O1b/6K+rN2CdftmwZbm5u7N27FyurmZuvOa4jq6rapChK\nLOADrFIUxaKv2G8Fhv2V+sknnxj/HBQURFBQ0IROVggx/zQ2NhITE0NNTc2Ex+gfNHHw4EGToImV\nT5709smfPDF7n3zVqlVoNBp27dpl/MUxUTExMcTExExqjLE85bIW6FJV9ZmiKEuAG8DPgG8BF1RV\n/aLvpmiOqqr/c4j9pYcuhBjEYDCQk5NDVlbWhEOaDQYDqampREZGsmvXLs6dO2cMmrBuacHzyhV2\npaaSc+QIeYcPm61PvnbtWjQaDTt27ECZokW5puqxRVfg3wCLvq8vVFX9XFGUHbx8bFEPvKOq6qA7\nGFLQhRADPX36lNjYWOrq6iY8RklJCeHh4SiKQlhYmDFoQunpYV9cHB5XrlCm1ZJx6pTZ+uSbNm1C\nq9WydetWs4w3EplYJISY1QwGA5mZmeTk5Ew4Dm64oAmArfn5+Oh0tK5cSXJYGA19sz8na9u2bWi1\nWuNTMtNBCroQYtaqqakhJiZmwnFwA4Mm3njjDeOTJCufPMFXp2NFTQ0pISFUuLlNuk+uKAo7d+5E\nq9WyZs2aSY010eNLQRdCzCrd3d2kp6eTl5c3oTi4gUETZ8+eNc62XNzSgsfly+xKSyP76FHyDx/G\nMMmnTCwsLNi9ezcajYYVUzD1f6ykoAshZpWqqiri4uImHAfXP2giLCzMGDSh9PTgFBeH9soVHnh4\nkHnqFO19qyROlJWVFfv27cPNzc04k3QmSUEXQswKnZ2dpKSkUFRUNKH9hwuaANial4dvRAQtq1aR\nHBo66T65tbU1Li4uuLi4GB91nA2koAshZlxFRQXx8fG0tLSMe98XQROpqakcOXKEV1991Rg0serx\nY3wiInr75KGhVLi6TqpPbmtri6urK05OTsZjzCZS0IUQM2YywRP9gya0Wi2nTp0y9q8XNzfjefky\njunp6I8doyAoaFJ98hcRb3v27DGutjgbSUEXQsyI+/fvk5iYSFtb27j26x80sWbNGkJDQ41BE0pP\nD06xsXhcvcp9Dw8yTp+mY5wrFvY3MOJttpOCLoSYVpMJnhguaAJVxT4vD1+djud2dqSEhtIwTKrP\nWKxbtw6tVjso4m22m7LVFoUQYqDi4mKSk5Pp7Owc135NTU1cunSJnJycQUETq6qq8NXpWP70Kcmh\noTx0cZlwn3w6Z3XOFnKFLoQYl+fPnxMfH8+jR4/GtV9XVxdfffUVN2/eHBQ0sbi5Gc/oaBwzM9Ef\nO0Z+UBDqBPvbMzGrcypIy0UIMWVUVSU/P5+0tLRxLXE7UtCE0tODc0wM2mvXuO/pScapUxPqk7+Y\n1anRaLCzsxv3/rORFHQhxJRobGwkNjaW6urqce1XXl5OeHg47e3thIaGGoMmBvbJk0NDaZxAn/zF\nrE53d3dWrlw57v1nM+mhCyHMymAwkJ2djV6vH9cStyMFTayuqsLnRZ88LKy3Tz5OVlZWxqzO2TCr\nc7aQK3QhxJDq6uqIjY3l6dOnY96nf9BEYGAgR48eNc6+XNzczP6oKHZmZZF1/DgFhw6Nu09ubW2N\ns7Mzrq6us2pW51SQlosQYtJ6enrIyMjgzp07Y15My2AwkJaWRmRkJI6OjiZBExbd3TjHxKC5do3S\nAwfIPHWKjnFeVdvY2ODm5oaTk9OMZHXOBCnoQohJefz4MXFxcTx79mzM+wwXNIGqsu3OHXwiImha\nv56UkBAaN20a1/ksXboUd3f3Gc/qnAnSQxdCTEhnZydpaWkUFBSMeZ/+QRNnz57Fy8vrZZ+8shJf\nnY6ljY0k/dmf8WicffIVK1ag0WjYvXv3nJjVOVvIFboQC9x4F9MaKWjC5vlzPKOje/vkJ05QEBg4\nrj75mjVr0Gq17Ny5c8qyOucKabkIIcZsvItpjRQ0YdHdjfPt22iuX6fEy4uskyfH1Sdfv349Wq0W\nBweHCX2W+UgKuhBiTEpKSkhKSqK9vX1M2xcWFhIREYGNjQ1hYWEvC6+q4tDXJ2/csIGUkBCebdw4\n5vPYvHkzHh4ebJ7EWi3zlfTQhRAjam5uJiEhgYqKijFtX11dTUREBFVVVYOCJtY8eoSPTofts2ck\nfu1rPHJ2HvN5ODg4oNVqWb9+/YQ+hxiaXKELsUAUFBSQmppKV1cXAPa5uVQ7OtJpa2vcxrq1lQ2l\npRTt3Mnly5dJS0vjyJEjHD582BgCYdPUxP7oaLbr9WSdPEnhwYNj6pPPdOjyXCMtFyHEII2NjcTF\nxfHkyROT161bWzkQGUn62bN02tpi3drK/osX+ZWdHRFffolWq+X06dMs78vqtOjqwqWvT37Px4fM\nEyfoHEOf3MLCgldeeQWNRjPvpudPJSnoQgijsUzbf1HUs994A/s//pH3q6tZtG4dISEhxqAJVBWH\nnJzePvnGjWPuk1taWhqn5y+bRDDFQiUFXQgBQG1tLbGxsdTX14+6bUt+Pv/1N7/hwNq1+Hztay+D\nJujtk/uGh7Pk+XOSQ0OpdHIadbxFixbh5OSEm5sbS5YsmfRnWajkpqgQC1x3dzfp6enk5eWNOm2/\nqamJ2xcvcjY9nR+ePMn/bmoiw9GRTkXBpqmJA1FROGRnk3nqFEUBAaP2yRcvXoyLiwsuLi4sXrzY\nnB9LjJFcoQsxTzx69Ij4+HieP38+4nYvgiZSb9zgn1avpvJ738Nq7VqsW1vxunCBllWrcLl9m3ve\n3mSdPGly03QoS5YswdXVFWdnZ+ONUzF5coUuxALU0dFBcnIyd+/eHXG7gUET/++5czTs34+VrS2o\nKpuLi9laUEDbihVc+ru/o2mUxJ+FvM7KbCVX6ELMYaWlpSQlJdHW1jbidmVlZeh0usFBE4Ddw4f4\nhodj09xMclgYlfv2jTiWrLMyPeQKXYgFYqwThF4ETRQUFHDmzBmToIklTU3sv3QJhzt3evvk/v4j\n9slXrVqFVqtl165dC36dldlKrtCFmENe5Hqmp6cbJwgNpbOzkxs3bnD79m0OHjzIsWPHjIEQll1d\nuHz1Fe43b1Ls54f++PER++Rr165Fq9WyY8cOs38eMTx5bFGIeay+vp64uDhqamqG3WakoAlUlR16\nPd7nz1O/ZQspwcEj9snXr1+Ph4cH27ZtM/dHEWMgBV2Ieainp4esrCxycnIwGAzDbldSUoJOpwMw\nDZoA7Coq8NXpWNzaSnJoKFX9eugDbd68Ga1W+3JikZgR0kMXYp4ZS4JQ/6CJc+fOceDAgZd98mfP\nOHDpEttyc8k4dYrigADUYW5k2tvbo9Vq2TiO1RLF7CJX6ELMQp2dnaSkpFBUVDTsNm1tbVy/fp34\n+PhBQRMmfXJ/f7KOH6drmFmbDg4OeHh4sG7duin5LGJipOUixDxQWlpKcnIyra2tQ75vMBhITEwk\nKioKZ2dnzpw5w+rVq3vfVFV2ZGXhff48T+3tSQ0OpmmIJWoVRWHHjh14eHjIyoezlBR0IeawsTyK\nOGzQBL19cr8vvmBRezvJYWE83rNn0P6KorBr1y60Wq0xbUjMTtJDF2IOUlWVvLw8MjIyhn0UcaSg\niSXPnuEVGYl9Xh4Zp09T7O8/qE9uYWHB7t270Wg0rFixYso/k5gZcoUuxEy4cgX8/anr7iYuLo66\nujpjuMRDV1fjZi0tLVy+fJnU1FSOHDnCq6++alwvxbKzE9cvv8Ttyy8pCghAf+zYoD65paUle/bs\nQaPRyBK2c8yUtFwURdkK/DuwEegB/peqqr9RFGU18AXgAJQBYaqqDroVLwVdiMG66+qo++53uREY\nSMeSJYPCJnp6eoiNjeXq1atotVpOnTr18spaVdmZmYn3hQvUbdtGSnAwzwfc0LSysmLv3r1oNBps\nR1lcS8xOU1XQNwIbVVXNVhRlGZAJnAH+HHiqquovFEX5AFitquqPh9hfCrqYX/qurunfg25shMRE\nOHFi1N3Ly8tJTEyks6aGA5GR5Lz5Ju43b5J+9iwdS5aQl5eHTqfDzs7ONGgCWFtejm94OEsbG0n6\nsz+jws3N+J51ayubHzxgxde/jru7u6xFPsdNy01RRVEigX/s+zqkqmp1X9GPUVV10GwFKehi3mls\nhI8+gs8/7y3qA78fRmtrK4mJiTx48MD42rK6Or7x0Uf84fPPKe7oQKfT0dDQQEhIiEnQhG1jIwci\nI7HPzyf9zBnKNBr2R0UZr+iXdnXx+u3brPzHf8RGniOfF6a8oCuKsh2IAVyAh6qqru733lNVVe2G\n2EcKuph/XhTxv/1b+OUvRyzmL9ZfycjIoLOz0/j6izZLop8f6/7P/+Gvnj8n4ORJAgMDsexbJMuy\nsxO3L7/E9csvKTx4kOyjR419cuvWVryjo+n667/G+epVLH/2sxF/oYi5ZUqfculrt0QAf62qarOi\nKGOu0p988onxz0FBQQQFBY3jFIWYhVat6i3mO3bAgwfDFtK6ujri4+Opra01ed26tRXPCxf4f1as\nIPI3v+E1T0+ie3rI9vam09Kyt0+ekYH3xYvUOjhw8Sc/MemTW1tb4+rpiePx41jv2TPiOYi5ISYm\nhpiYmEmNMaYrdEVRrIDLwDVVVX/d91ohENSv5XJbVdVBCynLFbqYl0a5Qu/q6iI9PZ38/PxBUXCq\nqtIZGcn/SEtjub09wcHBbNiwwfiUS/vy5fiGh2PV2UlSWBhPdu827rt48WJcXV1xcXHBurV1zP9K\nEHPPlLVcFEX5d6BOVdUf9Xvt50C9qqo/l5uiYkEZpYf+4MEDkpKSaGlpGbRr/6CJkJAQ9vULk7Bt\naMArMpIthYVknDnDXV9f4/PkNjY2xkK+aNGiCffxxdwxVU+5+ANxQC6g9n19CKQB4YA9UAGEqqra\nOMT+UtDF/DLMUy6tf/oTccuXDznT80XQRGFhIadPnzYJmrDs7MTtT3/C9auvevvkx47R1bd2uY2N\nDW5uboPzOif5pI2Y/WTqvxAzwGAwcOfOHbKysuju7jZ5b6SgCVQVx/R0vC5epGbHDtLeeovnfWuX\nD1vIxYIhBV2IaVZVVUVCQgKNjab/OB0xaAJY9+ABvjodll1dJIeF8eSVVwAp5OIlKehCTJO2tjaS\nk5MpKSkZ9N5IQRNLGxo4EBnJlqIi0k+f5q6vL1hYSCEXg0hBF2KKDfdMOYwcNGHZ2Yn7zZu43LpF\nQWAg2UeP0m1jw+LFi3F3d5dCLgaRgi7EFKqpqSEhIYG6ujqT19va2rh27RoJCQmDgiZe9Mm9L1yg\n2tGR1HPnaF671vj4oaurqxRyMSQp6EJMgfb2dtLS0galBxkMBhISEoiOjsbZ2ZmzZ8+arDG+/v59\nfMPDUQwGksPCqN61C2tra9zc3HqfI39R9IUYghR0IcxIVVWKiopIS0ujo6PD5L2RgiaW1tfjdfEi\nm+7eJf3sWe55e7Oo74rczc1NCrkYEynoQpjJcO2VkYImrDo6cL95E+fbtyk4dIjsI0dg2TJcXFxw\nc3N7+biiEGMgBV2ISWprayMtLY3i4mKT11taWrhy5QopKSmDgiYwGNjV9zz5k127SDt3jrb163Fy\nckKj0cgytmJCpKALMUEGg4H8/HwyMzNNnl4ZMWiCfn1yVSU5LIzaV15h7969eHh4SLCEmBQp6EJM\nQFVVFYmJiTQ0NBhfU1WV3NxcIiIihgyaMOmTnztHibc3u/sK+fLly2fiY4h5Rgq6EOPw/PlzUlJS\nTAInACorK9HpdNTX1xMaGmoSNGHV0YH7jRs4x8SQHxREzpEjODg54enpafKEixCTJQVdiDHo7u4m\nOzubnJwcenp6jK83NTURHR2NXq/n+PHjHDp0yBg0gcHAK2lpHLh4kce7d5N27hx2Gg0HDhzAzm5Q\nrosQkyYFXYhRlJSUkJqaarK0bVdXF7du3eLGjRv4+Phw4sQJli5danx/Q2kpvuHhACSFhWEZEICX\nlxcbNmyY9vMXC4cUdCGGUVtbS1JSEtXV1cbXVFUlKyuLCxcusGXLFmPQxAtL6+vxvnCBjSUlpJ07\nR+OxY3j5+LB169aZ+AhigZGCLsQAra2tpKWlcffuXZPXRwqasGpvR3PjBk6xseQfPsyDkBC0AQEm\ni2wJMdWkoAvRp7u7mzt37pCdnW2yRvmLoImCggLOnDljEjSBwcArqakciIzk8e7d5L79Nntff529\ne/e+3EaIaSIFXQjg3r17pKWlmfTJBwZNHD161GTCz4aSEnzDw1EtLMh4+202v/UWrq6uWFmNOUdd\nCLOSgi4WtCdPnpCcnExtba3xtdGCJpY9fYr3hQtsKC0lPSSEJe+9h8bDQ6bpixknBV0sSM+ePSM1\nNZWysjKT14cKmrDPzaXa0RGDhQWa69dxioujMCAASx8fdvzN38ikIDFrSEEXC0p7eztZWVkUFBRg\nMBiMr48UNGHd3Mwb//zPrHryhMp9+6h69128MzJY8qtfmQYuCzHDpKCLBaGnp4fc3Fyys7NN1l1p\na2vj+vXrxMfH8+qrr/Lmm2+aLFW78d693j65otCxcSOWv/oVm3//e/j8cynmYtaRgi5mtytXwN/f\ntHg2NkJiIpw4Meruqqpy7949MjIyaG5uNr5uMBhITEwkKioKZ2dnzpw5w+rVq43vL6ur6+2T379P\n9te/ztof/IA9NjYoO3fCgwewfbs5P6UQZjGRgi638MX08feHjz56eUXc2Pjy+1E8fPiQ1NRU6uvr\nTV7vHzTx/e9/n+39ivOi9nY016+zLy6Ogtdf59mvf42PtzdWzc29x33wAH75S7lCF/OGXKGL6fWi\niP/t346pmNbW1pKamkpVVZXJ6yMFTWAwsDs5mQOXLlG5bx9Pf/Qj3E+c6F3Otv8vkYG/VKSoi1lE\nWi5ibigrgx07Rmx3PHv2jPT0dO7fv2/yektLC5cvXyY1NXVw0ASw8e5d/MLD6V60iHv/+T+z7733\nTB5TnGzbR4jpIgVdzH6jXKG3traSmZlJcXGxyZMrowVNLK+rw/v8edaVlZH79tts+i//hR07d07r\nRxPCnKSgi9lthHZHx5Il5OTkkJeXZzJVX1VV8vLy0Ol0rFmzhtDQUJOgiUVtbb198vh4Ct58E6u/\n+ztcvLxkqr6Y86Sgi9ltiHZHd10d5X/4A/ErVpg8gggvgyYaGhoICQkxCZpQDAZ2JyWxPyqKSicn\n6v/mb3A/flzyO8W8IQVdzBkGg4HCwkKysrJoa2szea+pqYmoqCiys7M5ceIEgYGBWFpaGmd52j16\nhG94ON3W1pR+61u4ODqy8hvfmKFPIsTUkIIuZj1VVbl79y6ZmZkmz5LD6EETayoqePO3v8Wipwf9\n22+z5b332Pm738kTKmJekoIuZrXS0lIyMzNpbGw0eX20oIlFbW1or11jb0ICRUFBbLK2xu4Xv8Dq\nH/5BirmYt2RikZiVysvLSU9PHzQpCEyDJt555x2ToIkXffIDly7x0NmZlH/+ZzxOnmRFff3Lxx6l\nmAthJAVdTJlHjx6RkZFBTU3NoPdGDJoANhUX4xseTpeNDQkffMC+d98lyN6+98mYX/5SZnkKMQRp\nuQizq6qqIiMjgydPngx6b2DQxLFjx0zWHl9eW4tPRAR2Dx+SERrKmr/4C9zc3XuLvczyFAuI9NDF\njHry5AkZGRmDpunD6EETi9ra8Lh6lT2Jidx54w2a338f76Agk5uiMstTLCRS0MWMqKmpISMjg0eP\nHg35fv+gidDQUHbt2mV8TzEY2JOQwP7oaCpcXbn77rt4njxpMnlIiIVICrqYVrW1tWRmZlJRUTHk\n+/2DJs6ePYvXgBmcm4uK8NXp6LSxIe0b38Dh3Dnc3NxklqcQSEEX06Suro7MzEzKy8uHfH+0oIkV\nNTV4nz+P3cOHpAYHQ3Awvn5+LFu2bLo+ghCznhR0MaVGK+SjBU0samvD48oV9iQlkfPmm5SfPYvv\n4cPY29tP10cQYs6YkoKuKMq/AieBalVV3fpeWw18ATgAZUCYqqrPhtlfCvoc9/TpU8r/6Z+4s3w5\nnba2xtfEnRm1AAAV1UlEQVStW1vZUFrKQ1dXk6CJ0NBQk6AJpaeHvYmJePb1ybPeeos9hw6h0Wiw\ntLScgU8kxOw3VQU9AGgG/r1fQf858FRV1V8oivIBsFpV1R8Ps78U9Dnq6dOnZGZmUlZWhnVrKwci\nI0k/e5ZOW1vj91f8/PjD1atUVlYODpoANhcW4qvT0WFrS3JYGEv8/AgICDBZ+lYIMdiUtVwURXEA\novsV9CLgkKqq1YqibARiVFXdO8y+UtDnmLq6OrKysigrKzN5/UURz3nzTZyuXuVjS0u+yswcMmhi\nRXU1PufPs6aykpTgYGr8/fH188PR0XGaP40Qc9N0FvR6VVXX9Hv/qaqqdsPsKwV9jhitRw5gW1PD\nO3//97gsXYqdp+egoAnr1lY8rlxhd3IyOUeOkP/aa+zVaNi/f7/JjVEhxMhm7Voun3zyifHPQUFB\nBAUFTcdhxRiN9vgh9C6gVZqZyb7f/57XHR353erVFJw7Z+ypKz097E1IwDM6mnJ3d3SffMLSnTs5\nHRhoGgEnhBhSTEwMMTExkxpjolfohUBQv5bLbVVV9w2zr1yhz1I1NTVkZmby8OHDEberrKzk+h//\nyF8+fEjhO+/g6OnJ4rY2Y099XXk5vuHhtC9bRnJYGM8dHdm/fz/Ozs4m/XQhxNhN5RW60vf1QhTw\nHvBz4FvApfEcVMys6upqMjMzh53Z+UJTUxPR0dHo9Xo+0mgw/MVfsGv5cgA6bW0p9vPj2K9/zZLm\nZlJCQijTaNixcydH/fxMp+wLIabFWJ5y+QMQBNgB1cDHQCSgA+yBCiBUVdXGYfaXK/RZ4vHjx2Rm\nZg651kp/L4Imbt68ibe396CgCeuWlt4+eUoK2UeOkPfqq9iuXk1AQADbtm2b6o8hxIIgE4vEkCor\nK8nKyuLx48cjbjda0ITS08O++Hg8Ll+m3N2djDNn6Fi1ChcXF/bv34+VlazGLIS5SEEXJh4+fEhW\nVhbV1dWjbts/aCI0NJS9e02fQt1SUIBveDhtK1aQHBpKvb0969atIzAwEDu7IR9wEkJMwqx9ykVM\nr/LycrKysqitrR1129GCJlY+eYJPRASrnjwhJSSEcnd3Fllb43fggNz0FGKWkSv0eeTBgwdkZWXx\n9OnTUbcdLWhicUsLHpcv80pqKvpjx8gPCsKwaBHbt2/H399fbnoKMcWk5bIAqapKaWkper2ehoaG\nUbfvHzSxc+dO3nrrLZPnxJWeHpzi4vC4coUHWi0Zp07RvmIFtra2BAQEmKzRIoSYOlLQFxCDwUBJ\nSQl6vZ5nz4ZcF22QkYImALbm5+Or09GyciXJYWE0bNmCoijs27cPLy8vmekpxDSSHvoCYDAYKC4u\nJjs7m+fPn49pn9GCJlY+eYKvTsfKmhpSgoMpd3cHRWH16tUEBgaaPOkihJi95Ap9juju7qaoqIic\nnBxaWlrGtM9oQRMv+uS70tLIPnqU/MOHMVhZYWlpiVarRaPRSHqQEDNEWi7zUFdXFwUFBdy5c4e2\ntrYx7WMwGEhKSuLSpUtDBk0oPT04xcbicfVqb5/89Gna+2aAbtq0icDAQFauXDkln0cIMTZS0OeR\nzs5O8vLyyM3NpaOjY8z79Q+aCAsLw8HBAfvcXKodHem0tcU+Lw8fnY62FSso9fSksG+hNGtra3x8\nfAY9fy6EmBnSQ58H2tvbyc3NJT8/n87OzjHvV11dTUREBFVVVYOCJqodHQn4/e+xaW5mWX09GadO\nsenePUq9vADYuXMnfn5+2PZLIxJCzD1yhT5LtLa2kpOTQ2FhId3d3WPer6WlhcuXL5Oamjpk0MTi\n5mY8o6PZlZ5Ow8aNxH7zm7jeukX62bMsWreOgIAAHBwcpuIjCSEmQVouc1BzczPZ2dkUFxfT09Mz\n5Db9WyYvWD5/Tt2lS3yenY1Wqx0UNKH09OAcE4P26lXue3qScfo0i9rb+cZHH/GHzz9nW2CgPIoo\nxCwmBX0OefbsGXq9npKSEgwGw4jb9s/z7FiypDdo4j/+g/9pb8/Rr32NLVu2vNxYVbHPy8NXp6N5\nzRqSQ0Np2LLFOMaDkBAOZ2ay9B/+AVatmuJPKYSYKCnoc0B9fT1ZWVk8ePCA8fxcrFtbcfr97/mo\nsZFvVFYagyb6r6WyuqoKH52O5U+fkhISQoWrKygK1q2teF26RMdPf4pbYCCWz5/DRx/B559LURdi\nlpKCPovV1NSg1+tHzOsczougiaeZmeS1tPAf//2/07p+vfH9xc3N7I+OZmdGBvrjxyk4dAhDv6Vs\nXSsq2Puf/hOrd+x4OWhjIyQmwokTk/pcQoipIQV9FqqqqkKv11NZWTnufV8ETdy4cYPXPD35tKeH\nguPHcb95k/SzZ+m2tsYpJgbttWuU7t9P5qlTdCxbZtzfysqKAwcO4OLiIqsiCjHHSEGfRSoqKtDr\n9WNai3yggUET3zh+nBNJSaSfPUunrS3WLS0c/t3vWFVdTdO6dSSHhtK4ebPJGFu3buXgwYMs75sw\nJISYW6SgzzBVVXnw4AHZ2dnU1dVNaIzy8nLCw8Npb28nJCSEffv2mTzlsrqqCt/wcJY9fUqJlxdZ\nJ09Cv6vvxYsX4+vry+7du831sYQQM0AK+gx5sfJhdnY2jY1DRquOqn/QxOnTp/H39zdZR2VxczP7\no6LYmZVFVl+fXLW0NBnD0dERPz8/lixZMqnPI4SYeVLQp1lPTw/FxcXk5OSMeeXDgQYGTRw9etSk\nIFt0d+McE4Pm2jXjFXnHgHCJpUuXygQhIeYZKejTpLu727hgVmtr64TGGC1oAlVlW24uPhERNK1d\nS0poKI2bNg0aZ9++fXh7e8sEISHmGSnoU6yjo4P8/PxxL5g1UGlpKeHh4cDQQROrKyvx1elY2thI\nSkgID11cBo2xcuVKAgMD2TREkRdCzH2yONcUaWtr486dOxQUFNDV1TXhcUYLmrB5/hzP6OjePvmJ\nExQEBg7qkyuKgpubG56enlhZyf99QoiX5Ap9BM3NzeTk5FBcXDyuBbMGGi1owqK7G+fbt9Fcvz5s\nnxzAzs6OQ4cOmbZmhBDzkrRczOTZs2dkZ2dz7969UddZGYnBYCAxMZGoqKghgyZQVRzu3MEnIoJn\n69eTHBrKs40bB41jaWmJh4cH7u7ukiAkxAIhBX2Snj59il6vH/c6K0MZKmiiv9WVlfiGh2P77Bkp\noaE8cnYecpyNGzcSGBjIKllzRYgFRXroE1RdXY1er6eiosIsY50/f57KyspBQRMANk1N7I+OZrte\nT9aJExQO0ScHWLRoEV5eXjgPU+iFEGKgBX2F/ujRI/R6PY8fP570WC0tLVy5coWUlJQhgyZM+uTe\n3mSeOEHnEH1yAHt7ew4ePMiyfuuyCCEWFmm5jFFZWRl6vZ7a2tpJj9XT00NsbCxXr14dMmgCVcUh\nJwef8+dp3LCBlJCQIfvkINP2hRAvSUEfgcFgoLS0lOzsbBoaGiY9nqqq5OXlodPpWLNmDaGhoaZB\nE8CaR4/wDQ9nyfPnJIeGUunkNOx4O3bsICAgQKbtCyEA6aEPyRzT8weqrKxEp9NRX19PaGjooOVp\nbZqaOBAVhUN2NlknT1J48OCQfXIAW1tb/P392dF/rXIhhJiAeXuF3tXVRUFBAbm5uROenj/Qi6AJ\nvV7P8ePHOXToEJb9CrVFVxcut26huXGDu76+ZJ04YZIDOtDu3bvx9fVl8eLFZjk/IcT8IS0Xeqfn\n5+bmkp+fP6np+f31D5rw8fHhxIkTLO1/Q1NV2Z6djU9EBA2bN/f2yTdsGHa8ZcuWERgYyNatW81y\nfkKI+WdBF/TW1lbu3LlDYWHhpKbn9zcwaCI4OJgNAwq13cOH+IaHY9PcTHJYGJX79o04prOzM15e\nXiZPwAghxEALsqA3NTUZZ3X29PSYbdyhgib6W9LUxP5Ll3C4c4fMkycpCggYtk8OvYtpHTp0iI3D\nPOEihBD9LaiCXl9fT3Z2NqWlpZOe1dlfY2MjkZGR5OfnDxk0YdHVheutW7jfuEGxnx/648dH7JMr\nioK7uzuenp4m/XYhhBjJgijoNTU16PV6ysvLzTYm9AZN3Lx5k1u3bg0ZNIGqsl2vx+f8eeq3bCEl\nOJimEfrkAGvWrCEoKEgW0xJCjNu8LuiVlZXo9XqqqqrMcFYvjRo0AdhVVOCr07G4pYXk0FCqRumT\nW1hY4OHhgUajkcW0hBATMu0FXVGUo8D/ACyAf1VV9edDbDOpgl5WVkZ2djY1NTUTHmM4JSUl6HQ6\nYOigiSXPnnHg0iW25eaSeepUb598lAK9fv16Dh06ZLqqohBCjNO0FnRFUSyAu8BrQBWQDnxNVdWi\nAduNu6Cbe1bnQKMFTVh2deHy1Ve437xJsb8/WceP0zXKDE4rKyv279+Pq6urySQjIYSYiOmeKeoF\n3FNVtbzv4H8EzgBFI+41gqmY1dnfwKCJ9957zzSLU1XZkZWF9/nzPLW3J/LHP6Zp/fpRx920aROH\nDh0yXcNFCCGm2WQK+hbgYb/vH9Fb5MdtKmZ19tc/aMLJyYm///u/H9QSsauowO+LL1jU3k7st77F\n4z17Rh130aJFeHt74zTCGi1CCDFdJlPQh/qnwLh6K+3t7eTl5ZGXl0dnZ+ckTmV4RUVF6HQ6bGxs\n+P73v8/27dtN3l/y7BlekZHY5+WRcfo0xf7+o/bJQZa4FULMPpMp6I+Abf2+30pvL32QTz75xPjn\noKAgDhw4QE5ODkVFRZPK6hxJdXU1ERERVFVVDRk0YdnVheuXX+L2pz9R5O/PF599NmqfHGSJWyHE\n1IiJiSEmJmZSY0zmpqglUEzvTdHHQBrwdVVVCwdsZ7wpaq6szpG0tLRw+fJlUlNTefPNN3nttddM\np9n375Nv20ZKcDDP160b09jbt28nICAA2xEmEgkhhDnM1GOLv+blY4s/G2Ibta6uzmxZncMZNWgC\nWFtejm94OIva20kOCxtTnxxgyZIl+Pv7s3Pnzqk4dSGEGGTWTiz67W9/O2XjjyVowraxkQORkdjn\n55N+5gx3/fzG1CcH2LVrF35+ftjY2EzF6QshxJAWXMDFaEETlp2duH35Ja5ffklRQMCY++TQGzxx\n8OBBHBwcpur0hRDCrOZkQR8taAJVZWdGBt4XL1Lr4MDFn/xkzH1ygL179+Lj42P6jLoQQsxyc6qg\nDwya+PTTT02DJoB1ZWX4hodj1dnJ7ffe48k4nkZZvnw5Bw8elOAJIcScNCcKuqqq6PV6zp8/z5Yt\nW/jggw8GBU3YNjbidfEiWwoLyThzhru+vmPuk0Nv8IS3tzdWVnPiRyKEEIPM+upVXl6OTqejra2N\nd955Z1DQhGVnJ25/+hOuX31F4cGDhH/2GV3juIEpwRNCiPli1hb0hoYGLl26NGzQBKqKY0YGXhcu\nULNjBxc//JDn41h3XFEU3Nzc2L9/vwRPCCHmhVlX0AcGTXz22WemQRPAugcP8A0Px7K7m9vf/jZP\nXnllXMdYvXo1QUFBrBvHjVIhhJjtZk1BHxg08eGHHw4KmrBtaMArMpItRUW9z5P7+MA4+uQWFhZo\nNBo8PDwkeEIIMe/MioLeP2ji/fffHxQ0YdnZifvNm7jcukVhYCDhn346rj45wNq1azl06BB2dnZm\nO28hhJhNZrSgjxY0garimJ6O94ULVDs6cuHDD2keZz6npaUlHh4euLu7y1W5EGJem5GC3t7ezrVr\n14YPmqC3T+73xRdY9PTw1fvvUz3gqn0sJA5OCLGQTGtBNxgMJCUlERUVxb59+4YMmlja0IDXxYts\nLi4m7exZ7nl7j6tPDhIHJ4RYmKatoPcPmvje9743KGjCsrMT9xs3cLl9m4LAQL749FO6J7AglsTB\nCSEWqmlZbdHd3Z3KykreeustPDw8TK+aDQZ2pafjdfEiT3btIu3cOZoncONy0aJFeHl54ezsbMaz\nF0KImTFrV1t0dHTkO9/5jmnQBLD+/n18w8NRDIYJ98kBtm7dSmBgoMTBCSEWtGkp6EeOHDH5fml9\nPV4XL7Lp7l3SJ9gnB7C2tsbHx4e9e/ea61SFEGLOmtabolYdHbjfuIFzTAz5QUHEv/32hPrkANu2\nbePgwYODVlsUQoiFanoKusHAK2lpHIiM5PErr3D+v/03WtasmdBQixcvxs/Pj1fGOd1fCCHmu2m5\nKVrd90RLclgY1Y6OEx5rx44dBAQEDFrbRQgh5ptZe1M07/BhSry8JtQnB7CxsSEgIEBCmoUQYgTT\nUtBLfHwmvK+joyP+/v4S0iyEEKOYFYtzDcXW1paAgIBBE5CEEEIMbVYW9N27d+Pn5ychzUIIMQ6z\nqqAvW7aMgwcPYm9vP9OnIoQQc86sKej79u3D29tbrsqFEGKCZrygL1++nEOHDrF58+aZPhUhhJjT\nZrSgOzs74+3tjZXVjP9eEUKIOW9GInxWrlzJ6dOn8ff3l2I+TjExMTN9CvOG/CzNS36eM29aC7qi\nKLi5uREcHMzGjRun89DzhvxHYz7yszQv+XnOvGm7PF61ahVBQUGsX79+ug4phBALyrQUdI1Gg6en\nJ5aWltNxOCGEWJCmZXGuKT2AEELMU+NdnGvKC7oQQojpMSNPuQghhDA/KehCCDFPTEtBVxTlY0VR\nHimKktX3dXQ6jjufKIpyVFGUIkVR7iqK8sFMn89cpyhKmaIoOYqi6BVFSZvp85lrFEX5V0VRqhVF\nudPvtdWKotxUFKVYUZQbiqKsnMlznCuG+VlOqGZO5xX6r1RV9ej7uj6Nx53zFEWxAP4ROAI4A19X\nFEWSsSfHAASpqqpVVdVrpk9mDvrf9P597O/HwJeqqu4BbgE/mfazmpuG+lnCBGrmdBb0cd2tFSa8\ngHuqqparqtoF/BE4M8PnNNcpSMtxwlRVTQAaBrx8Bvi3vj//G3B2Wk9qjhrmZwkTqJnT+Rf6+4qi\nZCuK8i/yT7Fx2wI87Pf9o77XxMSpwA1FUdIVRfnOTJ/MPLFeVdVqAFVVnwDrZvh85rpx10yzFXRF\nUf6kKMqdfl+5ff97CvgnwFFVVQ3wBPiVuY67QAz1m1qeN50cP1VV9wPH6f0PJ2CmT0iIfiZUM802\nU1RV1TfGuOn/AqLNddwF4hGwrd/3W4GqGTqXeaHvChJVVWsVRblIb1srYWbPas6rVhRlg6qq1Yqi\nbARqZvqE5ipVVWv7fTvmmjldT7n0X4nrLSBvOo47j6QDuxRFcVAUxRr4GhA1w+c0ZymKYqsoyrK+\nPy8F3kT+Tk6Egum/HqOA9/r+/C3g0nSf0Bxm8rOcaM2crsW5fqEoiobeJwvKgO9O03HnBVVVexRF\n+QFwk95fwv+qqmrhDJ/WXLYBuNi3LIUV8HtVVW/O8DnNKYqi/AEIAuwURakAPgZ+BugURfk2UAGE\nztwZzh3D/CwPT6RmytR/IYSYJ+SxLSGEmCekoAshxDwhBV0IIeYJKehCCDFPSEEXQoh5Qgq6EELM\nE1LQhRBinpCCLoQQ88T/BcpRV5LmFP/wAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import bayespy.plot as bpplt\n", "bpplt.pyplot.figure()\n", "bpplt.plot(Fh, x=xh, scale=2)\n", "bpplt.plot(y, x=x, color='r', marker='x', linestyle='None')\n", "bpplt.plot(k*xh+c, x=xh, color='r');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the above plot shows two standard deviation of the posterior of the noiseless function, thus the data points may lie well outside this range. The red line shows the true linear function. Next, plot the distribution of the noise parameter and the true value, 2−2=0.25:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8TXf+x/HX174z6tco0diXoJSSxBrL1BJbLEWprbYw\njdGiCx1m2popZgTVohOk1F5C1NIisReNWiuxdFRrL9UiJCLf3x+JVgm5Se6933vP/Twfjzwe97on\n57x7mnzyvZ/7Pd+jtNYIIYSwphymAwghhHAcKfJCCGFhUuSFEMLCpMgLIYSFSZEXQggLkyIvhBAW\nlmGRV0p5K6W2KKW+VUodVkqFprNNU6XUNaXU/rSvcY6JK4QQIjNy2bBNMvCq1vqAUqoQEKuU+kJr\nHffAdtu01h3sH1EIIURWZTiS11pf0FofSHt8AzgGlE5nU2XnbEIIIbIpUz15pVRZoDawJ52X/ZVS\n3yilPldK+dohmxBCiGyypV0DQFqrZgUwIm1Ef79YwEdrnaCUagNEApXtF1MIIURWKFvWrlFK5QLW\nAuu11tNs2P5/QF2t9dUH/l0WyhFCiCzQWmepJW5ru2Yu8O2jCrxSyuu+x/VJ/eNxNb1ttdbypTXj\nx4/P/n4scj7tci4s8iXnQs5Fel/ZkWG7RinVEOgFHFZKfQNo4C3AJ7XG6DlAV6VUCHAHuAV0z1Yq\nIYQQdpFhkdda7wRyZrDNTGCmvUIJIYSwD7ni1ZDAwEDTEVyGnIvfybn4nZwL+7Dpg1e7HUwp7czj\nWZ5SIOdTCMtTSqEd/MGrEEIINyRFXgghLEyKvBBCWJgUeSGEsDAp8kIIYWFS5IUQwsKkyAshhIVJ\nkRdCCAuTIi+EEBYmRV4IISxMirwQQliYzXeGEp4rOTmZ2NhYzp49y/Xr17l+/Tq5cuWiSpUqVKtW\nDS8vL5SSW/wK4YpkgTJ35sAFym7cuMGiRYvYuHEjW7ZsoUyZMlSoUIHChQtTuHBhkpKSiIuL49ix\nY2it6dixI7169SIwMJCcOR+7MrUQIpOys0CZFHl35oAif/PmTT788EOmTJlCo0aN6NSpEy1btuSp\np5565Pf8+OOPLFu2jE8//ZTz588TEhLCyJEjKVSokF2zCeGpZBVKkW1aa+bPn0/FihXZu3cvmzdv\n5rPPPuOll156bIEH8Pb25tVXXyU2NpZNmzZx7NgxKlWqxMyZM0lKSnLSf4EQIj0ykndndhrJ37hx\ng2HDhhEbG8vChQt59tlns73P/fv38+abb3L69GkWLFhA/fr1s71PITyVjORFlh0+fJh69eqRK1cu\n9u7da5cCD1CnTh02btzIu+++S/v27Xn33XdJTk62y76FELaTkbw7y+ZIfseOHQQHBzNlyhT69u1r\nx2B/9OOPP9KvXz9u3brF0qVL8fb2dtixhLAiGcmLTNu+fTvBwcEsWrTIoQUeUnv2X3zxBe3atcPf\n35/9+/c79HhCiN/JSN6dZXEkv23bNrp06cLixYtp2bKlA4I92meffcbQoUOZO3cu7du3d+qxhXBX\n2RnJy8VQHmbfvn107drVSIEH6NKlC97e3gQHB/PDDz8wbNgwp2cQwpNIkfcgFy9epEuXLsyePdtI\ngb/Hz8+PnTt30qJFC+7cucOIESOMZRHC6qTIe4ikpCS6detGv379CA4ONh2HcuXKER0dTbNmzUhJ\nSWHkyJGmIwlhSVLkPcSrr75K0aJFmTBhgukov/Hx8SEmJoZmzZqhtebVV181HUkIy5Ei7wEiIiL4\n4osv2Lt3LzlyuNaEqqeffpqYmBiaNGlC0aJFefnll01HEsJSZHaNO7Nhds2ZM2eoW7cuW7ZsoWbN\nmk4Klnnx8fE0bdqU8PBwgoKCTMcRwqXIAmWeKoMir7WmVatWBAYG8tZbbzkxWNbs2bOHdu3asXbt\nWvz8/EzHEcJlyMVQIl3h4eFcvXqVMWPGmI5iEz8/P+bNm0enTp04ceKE6ThCWIKM5N3ZY0by99o0\n0dHR1KhRw8nBsmfOnDn85z//Yc+ePRQtWtR0HCGMk3aNp3pEkdda06ZNGxo3bszYsWMNBMu+4cOH\nc/r0adasWSM3IREeT9o14g8iIyM5e/as27Rp0hMWFkZCQgLjxo0zHUUItyZF3mISExMZNWoUU6dO\nJXfu3KbjZFnu3LlZvnw5S5cuZfHixabjCOG2pMhbzIwZM/D19TW6bIG9lChRgsjISEJDQzl69Kjp\nOEK4JenJu7MHevKXL1/G19eXHTt2UKVKFYPB7Gv+/PlMmjSJvXv3yn1jhUdy6AevSilv4BOgJHAX\n+FhrPT2d7aYDbYCbQD+t9YF0tpEib08PFPmQkBDy5s1LWFiYwVCOMWDAAJKSkliwYAFKZelnXQi3\n5egiXxIoqbU+oJQqBMQCHbXWcfdt0wb4i9Y6SCnlB0zTWvunsy8p8vZ0X5E/evQozZo1Iy4ujuLF\nixsOZn8JCQn4+fkRGhrKoEGDTMcRwqkcOrtGa33h3qhca30DOAaUfmCzjqSO9tFa7wGKKqW8shJI\nZM348eN5/fXXLVngAQoUKMCKFSt46623OHTokOk4QriNTH3wqpQqC9QG9jzwUmngh/uen+XhPwTC\nQQ4dOsTOnTsJCQkxHcWhqlSpwpQpU+jZsye3bt0yHUcIt2DzKpRprZoVwIi0Ef0fXk7nW9Lty9y/\n1G1gYCCBgYG2RhCP8M477zBq1CgKFChgOorD9enTh/Xr1zN69Gg++OAD03GEcIiYmBhiYmLssi+b\nZtcopXIBa4H1Wutp6bw+C4jWWi9Nex4HNNVaX3xgO+nJ25NSHD50iD//+c+cOnWKggULmk7kFNeu\nXaN27dp88MEHtGvXznQcIRzOGVe8zgW+Ta/Ap1kD9EkL4w9ce7DAC8e4N4r3lAIPUKxYMRYsWMCg\nQYO4cOGC6ThCuDRbZtc0BLYBh0ltwWjgLcAH0FrrOWnbfQC0JnUKZX+t9f509iUjeXtSipJeXh41\nir/f22+/zf79+1m7dq1MqxSWJguUeSqlmDxpEqNHjzadxIikpCT8/f0ZPny43FFKWJoUeQ90+vRp\nypYrx/Vff6Vw4cKm4xhz5MgRmjVrxr59+yhbtqzpOEI4hKxC6YGmT0+96NiTCzxAjRo1GD16NP37\n9yclJcV0HCFcjozk3dAvv/xCuXLluPrzzxne49UT3L17lyZNmtC9e3dCQ0NNxxHC7mQk72HCw8Np\n1aqV6RguI2fOnERERPCPf/xDbhsoxANkJO9mkpOTqVixIsuWLaO+n5+M5O8zdepUIiMjiY6OJkcO\nGb8I65CRvAdZtWoV3t7e1K9f33QUlxMaGsqdO3eYNWuW6ShCuAwZybuZgIAARo8eTefOnR97I29P\nFRcXR+PGjWW2jbAUGcl7iL1793LhwgU6duxoOorLqlq1Kq+99hqDBw9GBhRCSJF3K7Nnz2bIkCHk\nzJnTdBSXNmrUKK5evcq8efNMRxHCOGnXuIlr165RtmxZ4uPj8fJKW6pf2jWPdODAAZ5//nkOHz78\n+/kSwk1Ju8YDLFy4kFatWknBslHt2rUZMGAAI0aMMB1FCKOkyLsBrfVvrRphu/Hjx/P111+zdu1a\n01GEMEaKvBvYtWsXSUlJNGvWzHQUt5I/f37mzJnDsGHDuH79uuk4QhghRd4NzJ49m8GDB8tyulnQ\nvHlzWrZsydixY01HEcII+eDVxV29epXy5ctz8uRJSpQo8ccX5YNXm1y9ehVfX1+ioqKoV6+e6ThC\nZJp88Gphn3zyCUFBQQ8XeGGz4sWLM2XKFAYPHkxycrLpOEI4lRR5Fzdv3jwGDhxoOobb69WrF088\n8QTTpj3qDpZCWJO0a1zYgQMH6NSpE9999136C25JuyZTTpw4QUBAALGxsfj4+JiOI4TNpF1jURER\nEbz00kuyoqKdVKpUiZEjRzJ8+HBZ8kB4DBnJu6g7d+5QunRpdu7cSaVKldLfSEbymZaUlETt2rV5\n9913Uxd5E8INyEjegtavX0/lypUfXeBFluTJk4dZs2YxYsQImTsvPIIUeRcVERFB3759TcewpCZN\nmtCyZUvGjx9vOooQDiftGhd05coVKlSowPfff0/RokUfvaG0a7Ls8uXL1KhRgw0bNvDss8+ajiPE\nY0m7xmIWL15M27ZtH1/gRbb83//9HxMnTiQkJISUlBTTcYRwGCnyLkhaNc7Rv39/cuXKxZw5c0xH\nEcJhpF3jYo4fP06TJk04e/ZsxjcHkXZNth0+fJgWLVpw5MgRnnzySdNxhEiXtGssZMmSJbzwwgty\n9ycnqVmzJn369GH06NGmowjhEDKSdyFaa3x9fZk7dy4BAQEZf4OM5O3ixo0bVKtWjYULF9K0aVPT\ncYR4iIzkLeLgwYPcvn0bf39/01E8SqFChQgLC2PYsGEkJSWZjiOEXUmRdyGLFy+mR48esm68AZ07\nd8bHx4epU6eajiKEXUm7xkWkpKRQrlw51q5dS82aNW37JmnX2NWpU6fw8/OTBcyEy5F2jQXs3r2b\nwoUL217ghd1VqFCB0NBQ/vrXv5qOIoTdSJF3EYsXL6Znz56mY3i8MWPGcOTIEbn5t7AMade4gOTk\nZEqXLs2uXbuoUKGC7d8o7RqH+OKLLxg6dChHjhyhQIECpuMIIe0adxcdHY2Pj0/mCrxwmOeff556\n9erxz3/+03QUIbItwyKvlApXSl1USh16xOtNlVLXlFL7077G2T+mtS1fvpzu3bubjiHu85///IeP\nPvqI+Ph401GEyJYM2zVKqUbADeATrfUz6bzeFHhNa90hw4NJu+YhycnJPPXUU+zbt4+yZctm7pul\nXeNQU6dO5fPPP+fLL7+Uaa3CKIe2a7TWO4CfM8qQlYML2Lp1K2XLls18gRcO98orr3D58mWWLl1q\nOooQWWavnry/UuobpdTnSilfO+3TI6xYsYJu3bqZjiHSkStXLj766CNee+01fvnlF9NxhMgSm2bX\nKKV8gKhHtGsKASla6wSlVBtgmta68iP2I+2a+9y9e5dSpUplflbNPdKucYqBAwdSsGBBpk2bZjqK\n8FDZadfkyu7BtdY37nu8Xin1oVKquNb6anrbT5gw4bfHgYGBBAYGZjeC29q+fTulS5eWWTUu7v33\n38fX15e+fftSp04d03GEB4iJiSEmJsYu+7J1JF+W1JH8Q5djKqW8tNYX0x7XB5Zprcs+Yj8ykr/P\n8OHD8fb25s0338zaDmQk7zTh4eHMmTOHXbt2yTLQwukc+sGrUmoRsAuorJQ6o5Tqr5QaopQanLZJ\nV6XUEaXUN0AYIHMBbXD37l1WrlxJly5dTEcRNrh3F6mPP/7YdBQhMkWueDVk+/bt/OUvf+HgwYNZ\n34mM5J3q3l2kDh8+jJeXl+k4woPIFa9uaPny5TKrxs3UrFmTvn37yl2khFuRkbwBWmuefvppNm7c\niK9vNmacykje6W7cuIGvry8RERE0a9bMdBzhIWQk72a+/vprChQoQLVq1UxHEZlUqFAhpk+fTkhI\nCImJiabjCJEhKfIGrFq1iuDgYLlU3k117NiRypUrM3nyZNNRhMiQtGsMqFatGvPnz8fPzy97O5J2\njTHff/89devW5auvvqJixYqm4wiLk3aNG4mLi+PXX3+lXr16pqOIbPDx8eGNN95g2LBhyMBFuDIp\n8k4WGRlJp06dyJFDTr27GzFiBBcvXmTJkiWmowjxSFJpnGzVqlV06tTJdAxhB7lz52b27Nm89tpr\n/PxzRgu1CmGG9OSd6OzZs9SsWZOLFy+SO3fu7O9QevIuYdiwYaSkpDBr1izTUYRFSU/eTURGRhIU\nFGSfAi9cxsSJE4mKimLnzp2mowjxECnyTnRv6qSwlmLFihEWFsbgwYNJSkoyHUeIP5B2jZP8/PPP\n+Pj4cP78eQoWLGifnUq7xmVorWnfvj0BAQGMHTvWdBxhMdKucQPr1q0jMDDQfgVeuBSlFDNnzmTq\n1KkcP37cdBwhfiNF3klWr15Nx44dTccQDuTj48PYsWMZOnSozJ0XLkPaNU6QmJiIl5cX8fHx9l2i\nVto1Lic5OZmAgABCQkIYMGCA6TjCIqRd4+JiYmLw9fWVNcg9QK5cufjvf//LG2+8wYULF0zHEUKK\nvDNIq8az1KpVi4EDBxIaGmo6ihDSrnE0rTVlypRh06ZNVK1a1b47l3aNy7p16xa1atViypQpdOjQ\nwXQc4eakXePC9u/fT8GCBe1f4IVLy58/Px9//DHDhw/nl19+MR1HeDAp8g62evVqGcl5qKZNm9K2\nbVvGjBljOorwYFLkHUz68Z5t0qRJrF+/nujoaNNRhIeSIu9Ap0+f5vz58wQEBJiOIgwpWrQoH330\nEQMHDuTmzZum4wgPJEXegdasWUNQUBA5c+Y0HUUYFBQUREBAAOPGjTMdRXggKfIOFBUVJf14AUBY\nWBhLlixh9+7dpqMIDyNTKB3kl19+oUyZMpw7d45ChQo55iAyhdKtrFixgnHjxvHNN9+QP39+03GE\nG5EplC5o48aNNGrUyHEFXridrl278swzz/C3v/3NdBThQaTIO8iaNWukVSMeMnPmTBYuXChtG+E0\n0q5xgOTkZLy8vDh48CDe3t6OO5C0a9yStG1EZkm7xsXs3LmTsmXLOrbAC7fVtWtXatWqxdtvv206\nivAAUuQdICoqivbt25uOIVzYzJkzWbRoEdu3bzcdRVicFHkHkH68yEiJEiWYNWsW/fr14/r166bj\nCAuTnrydxcfH06JFC3744QeUylILzXbSk3d7AwYMIHfu3MyePdt0FOHCpCfvQtasWUP79u0dX+CF\nJYSFhbFx40bWrVtnOoqwKCnydib9eJEZRYoUYf78+QwaNIgrV66YjiMsSNo1dnTlyhXKly/PxYsX\nyZcvn+MPKO0ay3jttdf4/vvvWb58ubwLFA+Rdo2LWLduHc2bN3dOgReW8t5773H8+HEiIiJMRxEW\nk2GRV0qFK6UuKqUOPWab6UqpE0qpA0qp2vaN6D6kVSOyKl++fHz66aeMHj2aU6dOmY4jLMSWkfw8\noNWjXlRKtQEqaK0rAUOAWXbK5laSkpL48ssvCQoKMh1FuKmaNWsyduxYXnrpJZKTk03HERaRYZHX\nWu8Afn7MJh2BT9K23QMUVUp52See+9i2bRtVqlTBy8vj/tOFHYWGhlKwYEHee+8901GERdijJ18a\n+OG+52fT/s2j3Js6KUR25MiRg4iICD766CN27NhhOo6wAHsU+fQ+8fWoKR9aa7lBiLCbUqVKER4e\nTq9evfj558e9iRYiY7nssI8fgTL3PfcGzj1q4wkTJvz2ODAwkMDAQDtEMOvo0aNoralRo4bpKMIi\ngoKCCA4OZuDAgaxYsUKmVXqYmJgYYmJi7LIvm+bJK6XKAlFa65rpvNYWGK61DlJK+QNhWmv/R+zH\nkvPk//nPf3Lu3DlmzJjh3APLPHlLS0xMxN/fn6FDhzJkyBDTcYRB2Zknn+FIXim1CAgEnlBKnQHG\nA3kArbWeo7Vep5Rqq5Q6CdwE+mcliDtbs2YNf//7303HEBaTN29elixZQqNGjWjQoAE1az40xhIi\nQ3LFazZdunSJypUrc/HiRfLmzevcg8tI3iN88sknTJw4ka+//lpuJ+mh5IpXg9auXcuf//xn5xd4\n4TH69OlDw4YNGTp0KFYbJAnHkyKfTWvWrKFjx46mYwiLmzFjBgcPHiQ8PNx0FOFmpF2TDQkJCZQs\nWZLTp09TvHhx5weQdo1HiYuLo3HjxmzatIlatWqZjiOcSNo1hmzevJk6deqYKfDC41StWpVp06bR\npUsXrl27ZjqOcBNS5LNBWjXC2V588UXatGlD3759SUlJMR1HuAFp12RRSkoKpUqVYufOnVSoUMFM\nCGnXeKSkpCQCAwNp3749b775puk4wgkcOk9epG/v3r088cQT5gq88Fh58uRh2bJl1KtXj3r16tGy\nZUvTkYQLk3ZNFkmrRpjk7e3NokWL6N27N6dPnzYdR7gwKfJZtHr1almQTBjVrFkz3njjDTp16sTN\nmzdNxxEuSnryWXDy5EkaNWrEuXPnyJHD4N9J6cl7PK01/fr1IzExkcWLF8tCZhYlUyid7N5t/owW\neCFI/eWfNWsWp06dYtKkSabjCBckVSoLIiMjpVUjXEb+/PlZtWoV06ZN4/PPPzcdR7gYaddk0uXL\nl6lYsSIXL14kX758ZsNIu0bcZ/fu3XTo0IEtW7bIipUWI+0aJ4qKiuL55583X+CFeEBAQABhYWF0\n6NCBS5cumY4jXIQU+UxatWoVnTp1Mh1DiHT16tWL3r17ExwczO3bt03HES5A2jWZcOPGDUqVKsWZ\nM2coVqyY6TjSrhHpSklJoXv37uTJk4cFCxbIBAELkHaNk2zYsAF/f3/XKPBCPEKOHDmIiIjgu+++\n4+233zYdRxgmRT4TIiMjCQ4ONh1DiAwVKFCANWvWsGzZMubMmWM6jjBI2jU2SkpKomTJkhw+fJjS\npUubjpNK2jUiAydOnKBx48bMmzePNm3amI4jskjaNU6wdetWKleu7DoFXggbVKpUiVWrVtGnTx/2\n7t1rOo4wQIq8jWRWjXBXAQEBhIeH06FDB+Lj403HEU4mSw3bICUlhdWrV7N582bTUYTIkg4dOvDT\nTz/RqlUrdu7cKe9IPYgUeRt89dVX/OlPf6Jq1aqmowiRZQMGDODSpUu0atWKbdu2yW0rPYS0a2yw\nYsUKunbtajqGENn2+uuv06pVK9q0acP169dNxxFOILNrMpCSkkLZsmVZt24dNWrUMB3nj2R2jcgC\nrTVDhw4lPj6edevWUaBAAdORRAZkdo0D7du3j4IFC1K9enXTUYSwC6UUH374Id7e3nTp0oXExETT\nkYQDSZHPwL1WjdyMQVhJzpw5mT9/Pvnz56dnz57cuXPHdCThINKueQytNeXLlycyMpJatWqZjvMw\nadeIbEpMTKRr167kyZOHJUuWkDt3btORRDqkXeMgsbGx5M6dm2eeecZ0FCEcIm/evKxYsYLExEQZ\n0VuUFPnHkFaN8AR58+bls88+49atW/Tq1UsKvcVIkX8ErbVMnRQe416hT0hIoFu3bvJhrIVIkX+E\ngwcPorXm2WefNR1FCKfIly8fK1euJHfu3HTo0IGEhATTkYQdSJF/hCVLltCtWzdp1QiPkidPHhYv\nXoyXlxdt2rTh119/NR1JZJMU+XSkpKSwePFievXqZTqKEE6XK1cu5s+fj6+vL82bN+fy5cumI4ls\nkCKfjl27dlGkSBG5473wWDly5ODDDz+kTZs2NGrUiNOnT5uOJLJIFihLx6JFi3jxxRdNxxDCKKUU\n77zzDk8++SSNGzdm/fr1rre0h8iQTSN5pVRrpVScUuq4Uur1dF7vq5S6pJTan/Y1wP5RnePOnTss\nX76cHj16mI4ihEt45ZVXmDRpEi1atCA6Otp0HJFJGRZ5pVQO4AOgFVAd6KmUSm/N3SVa6zppX3Pt\nnNNpNm3aRKVKlShXrpzpKEK4jJ49e7J48WJ69OjBggULTMcRmWBLu6Y+cEJr/T2AUmoJ0BGIe2A7\nS0xDkVaNEOlr3rw50dHRBAUF8b///Y+3335bZp+5AVvaNaWBH+57/mPavz2os1LqgFJqmVLK2y7p\nnCwhIYGoqCi6detmOooQLsnX15fdu3cTFRVF7969uXXrlulIIgO2jOTT+1P94KpYa4BFWus7Sqkh\nQATQIr2dTZgw4bfHgYGBBAYG2hTUGaKiovD398fLy8t0FCFcVsmSJdm6dSsDBgygadOmREZGUqpU\nKdOxLCUmJoaYmBi77CvDVSiVUv7ABK1167TnbwBaa/3+I7bPAVzVWhdL5zWXXoWyQ4cOdOnShb59\n+5qOYhtZhVIYpLXmvffeY/bs2axcuZJ69eqZjmRZjl6Fch9QUSnlo5TKA/QgdeR+f4CS9z3tCHyb\nlTAmXbhwge3bt9OlSxfTUYRwC0opxo0bx/Tp02nbti3h4eGmI4l0ZNiu0VrfVUr9BfiC1D8K4Vrr\nY0qpvwP7tNZrgVClVAfgDnAV6OfAzA6xYMECOnfuTKFChUxHEcKtBAcHU7VqVTp37syePXuYMWMG\nefPmNR1LpJGbhpD6trN69ep8/PHHNGzY0HQc20m7RriQ69ev079/f86cOcPSpUtlGrIdyU1DsmnP\nnj3cvXuXBg0amI4ihNsqXLgwy5cvp2fPnvj5+bFy5UrTkQQykgdg8ODBlC9fnjfeeMN0lMyRkbxw\nUXv37qVHjx4EBQUxefJk8uXLZzqSW8vOSN7ji/zNmzcpU6YMR44ccb9pYFLkhQu7du0agwYNIj4+\nnkWLFsm6N9kg7Zps+Oyzz2jQoIH7FXghXFyxYsVYtmwZI0eOpFmzZoSFhZGSkmI6lsfx+JF8YGAg\noaGhdO7c2XSUzJORvHATp06d4qWXXqJAgQKEh4fj4+NjOpJbkZF8Fp04cYJvv/2Wdu3amY4ihKVV\nqFCBbdu20bJlS5577jlmz56Nqw34rMqjR/IjRoygYMGCTJw40XSUrJGRvHBDR48epX///hQpUoQ5\nc+ZQvnx505Fcnozks+D69essXLiQkJAQ01GE8CjVq1dn165dPP/889SvX5/333+fO3fumI5lWR5b\n5CMiImjevDllypQxHUUIj5MrVy7GjBnD3r17iY6O5rnnnuOrr74yHcuSPLJdk5KSgq+vL3PmzKFJ\nkyam42SdtGuEBWitWbJkCaNGjaJVq1b861//4sknnzQdy6VIuyaTvvzyS/LmzUvjxo1NRxHC4yml\n6NmzJ8eOHeOJJ56gevXqTJs2TVo4duKRRX7GjBmEhobKXW2EcCFFihRh8uTJbNu2jXXr1lGjRg3W\nrFkjs3CyyePaNSdPniQgIIAzZ86QP39+o1myTdo1wsI2bNjAqFGjePLJJ5k0aRLPPfec6UjGSLsm\nE8LCwnj55Zfdv8ALYXGtW7fmwIED9OjRg06dOtG1a1eOHTtmOpbb8aiR/Pnz56levTrHjh2zxi3+\nZCQvPERCQgIffPABU6ZMoW3btowbN46KFSuajuU0MpK30eTJk+nTp481CrwQHqRAgQKMGTOG48eP\nU65cOfz3s0WBAAAJSklEQVT9/enXrx8nT540Hc3lecxI/tKlS1StWtU9V5t8FBnJCw917do1pk+f\nzowZM2jZsiWvv/46tWvXNh3LYWQkb4N///vfvPjii9Yp8EJ4sGLFivG3v/2N7777jueee46goCBa\nt27Npk2bZDbOAzxiJP/TTz9RuXJlDh48aK0rXGUkLwQAiYmJLFy4kKlTp5IjRw7++te/8uKLL1rm\nZiVy05AMjB07lp9++onZs2c7/dgOJUVeiD/QWrNp0yamTp1KbGws/fv3Z8iQIW5/v1kp8o9x7tw5\nnnnmGfbt2+f2/6MfIkVeiEc6ceIEs2bNIiIiAj8/PwYNGkRQUBC5c+c2HS3TpMg/Rp8+fShVqhT/\n+te/nHpcp5AiL0SGbt26xdKlSwkPD+fkyZP06dOH/v37U7VqVdPRbCZF/hF2795N165diYuLo3Dh\nwk47rtNIkRciU+Li4pg7dy4LFizA29ub3r1707NnT5dfEE2KfDpSUlLw8/NjxIgR9O7d2ynHdDop\n8kJkyd27d9m8eTMLFiwgKiqK+vXr0717d4KDgylevLjpeA+RIp+OuXPn8t///pedO3dadyEyKfJC\nZNvNmzf5/PPPWbp0KZs2baJBgwYEBwfTsWNHl7lwUor8A65du0a1atWIioqy9qJGUuSFsKvr16+z\nYcMGVq5cyYYNG6hevTrt27enXbt2+Pr6GhswSpG/j9aaXr168ac//YmZM2c69FjGSZEXwmESExOJ\njo5m7dq1REVFkSNHDlq3bk2rVq1o3rw5RYoUcVoWKfL3iYiIYPLkyezbt8/6K01KkRfCKbTWHD16\nlI0bN7Jx40Z2795N7dq1ad68OS1atMDf3588efI47PhS5NMcP36chg0bEh0dTY0aNRx2HJchRV4I\nIxISEti5cyebN29m8+bNxMXFUb9+fRo3bkyTJk3w8/OjYMGCdjueFHlS31oFBAQwaNAgQkJCHHIM\nlyNFXgiXcO3aNXbt2sXWrVvZtm0bhw4dolq1ajRo0ICAgADq169P+fLls9zT9/gir7Vm2LBhXLhw\ngZUrV1p3Ns2DpMgL4ZJu375NbGwsu3btYs+ePezZs4eEhATq1atH3bp1qVOnDnXr1sXHx8emeuXx\nRX7s2LGsX7+eLVu2UKxYMbvv32VJkRfCbZw7d459+/axf/9+9u/fT2xsLLdu3eKZZ56hVq1a1KpV\ni+rVq1O9evWHLt706CI/ceJEPv30U7Zu3UqJEiXsum+XJ0VeCLd26dIlDh48yMGDBzl06BBHjx4l\nLi6OEiVK4OvrS9WqValWrRpDhgzxzCJ/76YB27Zt46mnnrLbft2GFHkhLOfu3bucPn2aY8eO/fY1\nb948zyryt2/fZvTo0axbt44tW7bg4+Njh3RuSIq8EB7B4XeGUkq1VkrFKaWOK6VeT+f1PEqpJUqp\nE0qp3Uqpp7MSxhbHjx8nICCACxcuEBsb67kFXgghbJBhkVdK5QA+AFoB1YGeSqkH1+h8Gbiqta4E\nhAGT7B30xo0b/Pvf/6Zhw4YMHjyYZcuWufWHrDExMaYjuAw5F7+Tc/E7ORf2YctIvj5wQmv9vdb6\nDrAE6PjANh2BiLTHK4AW9gp45coVJkyYQLly5dizZw8xMTGEhIS4/TRJ+QH+nZyL38m5+J2cC/vI\nZcM2pYEf7nv+I6mFP91ttNZ3lVLXlFLFtdZXMxPm7t27XLx4kfj4eLZs2cLmzZs5fPgwL7zwAjt2\n7KBKlSqZ2Z0QQng8W4p8ekPmBz/te3Ablc42AAQFBaXuQGsSExO5ffs2iYmJ/PTTT5w/f55ixYpR\noUIFAgMDeeedd2jQoIH116ARQggHyXB2jVLKH5igtW6d9vwNQGut379vm/Vp2+xRSuUEzmutH7rV\nilJKpoIIIUQWZHV2jS0j+X1ARaWUD3Ae6AH0fGCbKKAvsAfoBmyxZ0ghhBBZk2GRT+ux/wX4gtQP\nasO11seUUn8H9mmt1wLhwAKl1AngCql/CIQQQhjm1IuhhBBCOJdNF0NllitdPGWaDedipFLqqFLq\ngFLqS6VUGRM5nSGjc3Hfdl2VUilKqTrOzOdMtpwLpdQLaT8bh5VSC52d0Vls+B0po5TaopTan/Z7\n0sZETkdTSoUrpS4qpQ49ZpvpaXXzgFKqtk071lrb9YvUPxwnAR8gN3AAqPrANiHAh2mPuwNL7J3D\nFb5sPBdNgXxpj4d68rlI264QsBXYBdQxndvgz0VFIBYokva8hOncBs/FbGBI2uNqwP9M53bQuWgE\n1AYOPeL1NsDnaY/9gK9s2a8jRvJGL55yMRmeC631Vq317bSnX5F6zYEV2fJzAfAO8D6Q6MxwTmbL\nuRgEzNRa/wqgtf7JyRmdxZZzkQLcu6FqMeCsE/M5jdZ6B/DzYzbpCHyStu0eoKhSyiuj/TqiyKd3\n8dSDhesPF08B15RSxR2QxTRbzsX9XgbWOzSRORmei7S3n95a63XODGaALT8XlYEqSqkdSqldSqlW\nTkvnXLaci78DLymlfgDWAq84KZurefBcncWGQaEtUygzy64XT7k5W85F6oZK9Qbqktq+saLHnguV\nuk7FVFKn4j7ue6zAlp+LXKS2bJoATwPblVLV743sLcSWc9ETmKe1npp23c5CUtfR8jQ215P7OWIk\n/yOpP5T3eAPnHtjmB6AMQNrFU0W01o97m+KubDkXKKVaAm8C7dPeslpRRueiMKm/uDFKqf8B/sBq\ni374asvPxY/Aaq11itb6NBAPVHJOPKey5Vy8DCwD0Fp/BeRTSnnYHYKA1HN1/8SMdOvJgxxR5H+7\neEoplYfUOfNrHtjm3sVT8JiLpywgw3OhlHoWmAV00FpfMZDRWR57LrTWv2qtn9Ral9dalyP184n2\nWuv9hvI6ki2/I5FAc4C0glYJ+M6pKZ3DlnPxPdASQClVDchr4c8oFI9+B7sG6AO/rURwTWt9MaMd\n2r1do+Xiqd/YeC4mAQWB5Wkti++11p3MpXYMG8/FH74Fi7ZrbDkXWuuNSqnnlVJHgWRglBXf7dr4\nczEK+FgpNZLUD2H7PnqP7ksptQgIBJ5QSp0BxgN5SF1GZo7Wep1Sqq1S6iRwE+hv037TpuMIIYSw\nIIdcDCWEEMI1SJEXQggLkyIvhBAWJkVeCCEsTIq8EEJYmBR5IYSwMCnyQghhYVLkhRDCwv4fn+sT\n+2lTV2MAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "bpplt.pyplot.figure()\n", "bpplt.pdf(tau, np.linspace(1e-6,1,100), color='k')\n", "bpplt.pyplot.axvline(s**(-2), color='r');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The noise level is captured quite well, although the posterior has more mass on larger noise levels (smaller precision parameter values). Finally, plot the distribution of the regression parameters and mark the true value:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAESCAYAAAD67L7dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdUVEf7B/Bno8nvzfvGBrvLwtLL0pGigGIBFVQUjA3F\nblQUxYIES0QQsSS2KJYodiyoIFZCrKBYUESxUiwISi8qIkjZ/f7+IN7jBjQYxaDO55w98d47e2eY\nc/Ls7tyZZ3gAiGEYhvlyfPVvN4BhGIb5uFjgZxiG+cKwwM8wDPOFYYGfYRjmC8MCP8MwzBeGBX6G\nYZgvTIMHfh6PN5XH49388zWloetjGIZh3q5BAz+PxzMmojFE1IaIzInIhcfj6TRknQzDMMzbNfQ3\nfkMiigdQAUBKRGeIqG8D18kwDMO8RUMH/ltE1InH47Xi8Xj/JSJnIlJr4DoZhmGYt2jakDcHkMLj\n8X4hopNE9JyIkoiouiHrZBiGYd6O9zFz9fB4vIVE9AjA+r+cZwmDGIZh3hEA3j9538eY1SP487/q\nVDO+H1ZXOQDs9QFeAQEB/3obPqcX60/Wn4319T4adKjnT/t5PJ4CEVUR0UQAzz5CnQzDMMwbNHjg\nB9CpoetgGIZh6o+t3P3M2Nvb/9tN+Kyw/vywWH82Dh/14e4bG8HjoTG0g2EY5lPB4/EIjfXhLsMw\nDNO4sMDPMAzzhWGBn2EY5gvDAj/DMMwXhgV+hmGYLwwL/AzDMF8YFvgZhmG+MCzwMwzDfGFY4GcY\nhvnCsMDPMAzzhWGBn2EY5gvDAj/DMMwXhgV+hmGYLwwL/AzDMF8YFvgZhmG+MCzwMwzDfGFY4GcY\nhvnCNHjg5/F43jwe7xaPx7vB4/F28Xi8bxq6ToZhGObNGjTw83g8FSKaTESWAMyoZnP3wQ1ZJ8Mw\nDPN2TT9CHU2I6H88Hk9GRP8louyPUCfDMAzzBg36jR9ANhEtJ6JMIsoioqcATjZknQzDMMzbNfRQ\nT0si6kNEGkSkQkTf8Xi8IQ1ZJ8MwDPN2DT3U042IHgAoJiLi8XiRRNSeiHb/teC8efO4f9vb25O9\nvX0DN41hGObTERsbS7GxsR/kXjwAH+RGdd6cx7Mmos1E1JaIKohoKxElAFj7l3JoyHYwH86jR49o\n//79dPbsWUpJSaG8vDx6+fIlNW3alFq0aEEqKiqkra1NRkZGZGlpSTY2NqSoqPhvN5thPjs8Ho8A\n8P7Rexs64PJ4vACqmclTRUTXiGgsgKq/lGGBv5ErKCggX19fOnLkCPXt25e6du1KJiYmJBKJ6Ntv\nv6Xq6mp68uQJZWdn0/379+nWrVuUmJhICQkJpKmpSY6OjtSzZ0/q3Lkzff311//2n8Mwn7xGHfjr\n1QgW+Bu1O3fuUI8ePcjNzY38/f2pefPm9X5vdXU1JSYm0vHjxykqKorS0tLIxcWF3N3dqVu3btS0\n6ceYWMYwnx8W+JkG8+TJE7KwsKD58+fTiBEj3vt+WVlZFBERQbt376aMjAwaPnw4jR07lvT19T9A\naxnmy8ECP9NgvL29qby8nNavX//B752amkpbtmyhbdu2kYmJCXl5eZGrqys1adLkg9fFMJ8bFviZ\nBlFWVkZisZhu375NKioqDVZPZWUl7d+/n4KDgyk3N5emTZtGY8aMoe+++67B6mSYT937BH6WpI15\no7Nnz5KZmVmDBn0iom+++Ybc3d3p4sWLFBYWRnFxcaSlpUWBgYH05MmTBq2bYb5ELPAzb3Tjxg2y\nsrL6qHXa2tpSREQEnT9/njIyMkhXV5f8/PyouLj4o7aDYT5nbEoF80b5+fmkrKz8j95bXV1NKSkp\nlJGRQc+fP6cmTZpQq1atSCwWk7a2Nv3f//3fW98vkUhoy5YtlJ6eTgsXLiSJREJTpkwhb29vatas\n2T9qE8MwNVjgZ96ourr6nadb3rx5k5YsWUJHjhwhJSUl0tbWpmbNmpFUKqWnT5/So0eP6NGjR6St\nrU22trbUqVMncnR0fONwkpaWFm3atIlmzZpFAQEBpKenR35+fuTh4UHffMMyfDPMP8GGepg3at68\nOT179qxeZQHQkiVLqFu3btS6dWtKTU2l1NRUio6Opn379tH+/fvp1KlTlJaWRs+ePaMdO3aQpaUl\nHT16lExMTMjS0pIWLlxId+/erfP+urq6tGvXLvrjjz/oyJEjZGxsTAcOHCA2KYBh3h2b1cO8UUhI\nCMXHx9OWLVv+tuzatWtp3bp1dPz4cRKLxe9UT3V1NZ07d472799P4eHhpKmpSaNHjyZ3d/c3LhY7\nfvw4+fj4EJ/Pp5UrV1Lr1q3fqU6G+dS9z6weAvCvv2qawTQ2Z8+ehbW19d+Wy83NhYKCAu7evfve\ndVZVVSEqKgr9+vVDy5YtMWHCBNy+ffuNZdeuXQuhUIiJEyeiqKjovetnmE/Fn3HzH8VcNtTDvJGF\nhQXdunWLKioq3louNDSU+vbtS7q6uu9dZ9OmTcnZ2Zn2799Pt2/fJiUlJerSpQv17NmTYmNj5YZ2\nmjZtShMnTqTk5GQCQEZGRrRlyxaSyWTv3Q6G+ZyxoR7mrWxsbOjnn38mBweHN5ZxdnYmDw8P+v77\n7//2fg8ePKDExETKysqiyspK+u9//0tKSkqkq6tLhoaG9J///KfWe16+fEk7d+6kJUuWkKKiIvn7\n+1OPHj2Ix5P/lZuYmEgTJ06kb775htavX0/Gxsbv/gczzCeCLeBiGkyPHj0oKirqrWUePHhAEonk\nrWWioqLIysqK7OzsaPfu3fTgwQMqLCyk5ORkCgsLo+HDh5OCggK1bduWfH196dSpU1RVVZPE9T//\n+Q+NHTuWkpOTadq0aeTr60u2trZ04sQJuV8AVlZWdOHCBRoyZAjZ29uTn58fvXz58v07gWE+N/90\njOhDvoiN8Tda165dg6amJmQy2RvLqKmp4eHDh2+8Pm/ePGhpaeHw4cOQSqVvLFdWVoazZ88iICAA\n1tbWUFBQwKhRo3Dy5Em590mlUoSFhUEikcDBwQGXL1+uda/s7Gz069cP+vr6OHv2bD3/Wob5dNB7\njPH/60EfLPA3ajKZDIaGhjh37twbyxgYGODWrVt1Xjty5Ai0tbWRn5//znU/evQIy5cvh7m5OTQ0\nNBAUFIScnBzuelVVFTZs2ABlZWW4u7vX+eGzf/9+qKiowMvLC8+fP3/nNjBMY/U+gZ8N9TBvxePx\naNSoUbRp06Y3llFTU6PMzMw6r82fP59+/fVXEggE71y3qqoqTZ8+na5du0aRkZGUmZlJhoaGNHLk\nSLpx4wY1bdqUPDw8KC0tjSQSCVlaWtLcuXPpxYsX3D369etHN2/epJKSEmrdujWdOXPmndvBMJ8b\n9nCX+VsFBQUkkUjo3r17dW6j6O3tTSKRiGbOnCl3Pj8/n/T19amgoOCNK4DT0tLo0KFDdOvWLSoq\nKqImTZqQgoICaWlpkZmZGdnY2MiljSguLqYNGzbQ6tWrydLSkvz8/MjW1paIaraFnDFjBp0/f55W\nrFhB/fv3l3sAfOTIEZowYQINGjSIFi5cSN9+++2H6B6G+Vewh7tMgxIIBNSnTx/asGFDnddtbW3p\n3Llztc6npaWRgYFBnUG/tLSURowYQR07dqSHDx9S586dafz48TR69Giys7Oj8vJyCgkJIWNjYzIy\nMqIff/yRzp8/Ty1btqTZs2fTgwcPqFevXjRo0CDq2bMnXblyhdTU1CgsLIx27txJgYGB1KNHD7p3\n7x5Xp4uLC924cYOysrLIysqKrl69+uE6iWE+Jf90jOhDvoiN8Td6N2/ehEgkQllZWa1r+fn5aN68\nOV6+fCl3/tixY+jatWut8tXV1ejSpQtGjhyJ0tLSt9YrlUqRkJAAf39/GBsbQ11dHX5+fkhPTwcA\nvHz5EmvXroWKigoGDBiAtLQ0AEBlZSWWLFkCRUVFLFy4EBUVFXL33bVrFwQCARYvXozq6up36QqG\naRSosT7cJSIJ1WywfvXP/z4joil1lGuovmE+IFdXV6xatarOax06dMCRI0fkzl28eBFt27atVXbn\nzp1o167dPwq4SUlJmDp1KhQVFdG7d2+cPHkSMpkML168wKJFi6CoqIipU6dyq3jT09PRs2dPmJqa\n1pr9k5GRgc6dO6Nz587IzMx857YwzL+p0QZ+uYpqhpWyiUitjmsN0jHMh5WYmAhlZWW8ePGi1rW1\na9di0KBBcudycnKgoKBQayqos7Mz9u7dW2cdVVVViI2NxfLlyzF79mz4+fkhODgY0dHRyMvL48q9\nePECISEhMDQ0hKWlJSIiIiCVSpGfn48JEyZAKBRi3bp1qK6uhkwmw65du6CkpIQZM2agvLycu091\ndTUWLlwIoVCIyMjI9+kehvmoPpXA70REcW+41hD9wjSA/v374+eff651vqioCC1atEBBQQF3TiaT\nQSQS4cGDB3Jl1dTUap0DgBMnTkBLSwuWlpaYPHkygoKCMG/ePHh6eqJLly5o0aIF9PX1MXXqVMTE\nxKC6uhpSqRQHDx5EmzZtYGJigv3790MmkyEpKQmdOnWChYUF4uPjAQB5eXkYOHAgDAwMcOnSJbm6\nL168CC0tLUyaNEnug4FhGqtPJfBvJqKJb7jWEP3CNICUlBTw+Xy5AP/KqFGjsGjRIrlzgwcPxqZN\nm+TOKSgo1Hr/2bNnIRQKcfz48TfWXV1djcTERAQFBcHc3BwqKiqYNWsW7t27B5lMhqNHj8Lc3Bxt\n27ZFTEwMZDIZduzYAZFIBE9PTzx9+hQAsGfPHgiFQvj5+aGyspK7/5MnTzBgwACYm5tzzwoYprFq\n9IGfiL4mogIiErzhOgICArhXTExMQ/QT84FMnjwZEydOrHX++vXrUFZWlnvIu2PHDvTu3VuunJaW\nVq1Mnm3btkV4ePg7tePWrVuYPn06+Hw+evXqhdOnT6O6uhq7d++GpqYmXF1dkZaWhuLiYowbNw5i\nsRgHDhwAULOyt2fPnrCyskJKSgp3T5lMhjVr1oDP579zeximIcXExMjFyU8h8LsS0R9vud4Q/cQ0\nkMLCQgiFQiQlJdW65uzsjHXr1nHHz549qzUE1K1bN0RFRXHHOTk5aNWqVZ0Pe/Py8rBu3TqMHz8e\n7u7uGDt2LIKCgnDo0CFuNXBZWRlCQkKgr68Pa2trHD16FOXl5fjll1+gqKiImTNn4vnz5zhz5gz0\n9PQwaNAg5OfnQyaTYd26dVBUVMSGDRvknkUkJCRAS0sLU6dOlftVwDCNxfsE/o81j9+diMI+Ul1M\nA1NUVKT58+fTpEmTaqVA9vf3p0WLFnHJ0Zo3b04uLi4UGhrKlbG2tqYLFy5wxw8fPiRdXV1q0qSJ\n3L1CQ0PJ0NCQLly4QKampuTs7ExWVlZUWlpK69atI11dXWrTpg2tWLGCunbtSnfu3KEff/yRZs2a\nRfb29tSmTRu6efMmPX78mIyMjKioqIiuX79OqqqqZGZmRgcPHiRPT0+Ki4ujdevW0cCBA+nJkydE\nRNSmTRtKTEyk+/fvk729PWVnZzdUdzLMx/dPPzHq+yKib6lmmKfZW8o0yCci03CkUimsra2xcePG\nWtdcXV2xbNky7vjChQvQ1tbmvtHHxMTA0tKSu37p0iW0adNG7h6XLl2CSCTCnTt33tiGyspKxMTE\nYOLEiRAIBLC3t8e+fftQUVGB3bt3Q1tbGz179sStW7cQExMDfX199OnTB48ePcL58+ehq6uLESNG\n4OnTpygvL8eUKVOgoaGBCxcuyP2dQUFBUFFRYcnemEaFGvtQz982ggX+T1JSUhIEAoFc4jQAuH37\nNgQCAQoLC7lz7du3x549ewDUTNkUCoXcOH9WVhYEAoHcPYYNG4bg4OB6t6WiogJ79uxBhw4doKGh\ngeDgYDx9+hS//vor+Hw+Jk+ejJycHPj7+4PP5yMkJATPnz/H+PHjoampyQX1gwcPQigUYtmyZXJD\nP9HR0RAKhVi9evVbM5UyzMfCAj/zr5k9ezb69etXKxhOnDgRkyZN4o6joqJgYmLCpVeeNm0aZs+e\nDaDmgaqSkhK3GhcAjIyMcP36dbl7SqVSbNu2DT169ICuri60tLRgY2OD0aNHIyQkBBkZGQCA+Ph4\n9OnTB8rKyggODsbjx48xfvx4iEQihIaG4vr167CysoKjoyMyMzNx5MgRiEQizJ07F1VVVUhPT4e1\ntTX69OmDJ0+ecPXfu3cPpqamGD16dK1VygzzsbHAz/xrysvLYWhoiLCwMLnzrx4AX716FUBNcLe1\ntcXOnTsB1EwLFQqFXAqIIUOGYP369dz7dXR05KZUymQyuLm5wdbWFhEREUhJScG9e/dw7tw5/Pbb\nbxgyZAgUFRXRtm1brFq1CkVFRbh69SqcnZ2hqamJsLAwXLp0CZaWlujSpQuSk5MRFBQEgUCAHTt2\nICsrC46OjujQoQMyMzNRUVEBLy8v6OjoyD3Efv78Ofr37w9bW9tav3QY5mNigZ/5V126dAlCoRDZ\n2dly5zdt2gRra2tubP/MmTPQ0NDgFki5uLhg7dq1AICIiAg4ODhw723fvr3ctN6oqCiYmpq+9Zt2\nVVUVjh07Bnd3d7Rs2RIeHh5ITU1FbGwsLC0t0a5dO8THx2PZsmVQVFTE8uXLkZCQACMjIwwaNAiF\nhYVYvHgxlJSUcPToUQDA7t27wefzsWPHDq4emUyGwMBAqKmpITEx8f06j2H+IRb4mX/d3Llz0aNH\nD7khH6lUio4dO8rl9+nbty8WLFgAoGZIRk1NDeXl5Xj58iUEAgE37j9p0iQsXbqUe9+UKVOwZMmS\nOuuWyWQoLS1FVVUVdy4vL48bzx88eDDu3LmDTZs2QUlJCRMnTsS1a9fQqVMntG/fHjdv3oSXlxc0\nNDRw7tw5xMXFQVVVFbNnz0ZVVRVu3LgBHR2dWlM7IyIiwOfzERER8WE6kWHeAQv8zL+usrISbdu2\nrfVANiUlBYqKilyKhgcPHkBRUZHbLcvFxQXLly8HAMyaNQteXl4Aah6yvv4LYMSIEdi6davcvYuK\nijB27Fi0bNkS3377LZo2bQoNDQ24ublh06ZNKCgoQElJCRYtWgQ+n49x48YhOTkZHh4eUFFRQURE\nBFauXAk+n48NGzbg0KFDUFJSwsKFC5Gbm4tu3brBwcEBeXl5KC4uRo8ePdClSxe5h9aJiYlQVVXF\nokWL2ENf5qNigZ9pFNLS0sDn82s9lF26dCk6d+7MPdgNCgpC7969IZPJ5GYAZWdno1WrVsjJyUFZ\nWRlatWrFZc308fHB4sWLuXvKZDJ06NABEyZM4IaYqqurkZaWhq1bt2LAgAFo0aIF+vbti5MnT6Ko\nqAg+Pj5QVFTEL7/8glOnTkEikcDNzQ1xcXEwNzdHnz59cOPGDXTo0AHdu3dHTk4O5syZAzU1NVy6\ndAnV1dXw9fWFtra23FaTjx8/hoWFBUaPHl0r/TPDNBQW+JlGIzQ0FAYGBnL721ZXV8POzo77Zl9R\nUQFjY2NueqeXlxfGjx8PAJg6dSomT54MoGZmkJ+fH4Ca/Pnff/89d8/ExETo6em9dfP2Z8+eYf36\n9TAyMoK5uTkiIiKQmpoKZ2dnGBoa4vjx4/Dx8YGysjL2798PHx8fqKqq4vTp05g5cybU1NQQHx+P\nAwcOgM/nc784QkNDIRAIuOcAAFBaWgoXFxd06dJFbiYQwzQUFviZRmXkyJEYPny43NDH/fv35X4N\nxMfHQ0lJiRtGEYlEuHjxIvLy8qCoqIjU1FSkpqaCz+ejpKQE+fn5aNGiBZcSOiIiQu6D4G1kMhkO\nHz4MKysrWFhY4Pjx49i/fz/EYjE8PDwQHR0NTU1NjB8/HpGRkVBSUsKCBQsQGRkJgUCA9evX4/bt\n29DT08PUqVNRVVWFCxcuQFlZGb/++iv3d1ZXV2Py5MkwNjbmppYyTENhgZ9pVEpLS2FsbIyQkBC5\n89u3b4ehoSEXvGfNmoU+ffpAJpNh9+7dMDExQUVFBZYsWQJnZ2fIZDK4u7sjKCgIQE0eoG3btgEA\nTp48iU6dOtWqe/v27WjTpg0EAgG0tbXh6uqKX3/9FRkZGZDJZAgPD4eOjg569+6Na9euYezYsVBX\nV8ehQ4cwdOhQ7peAnZ0dnJ2dcfnyZRgbG2Ps2LHIycmBk5MTnJyc8OTJEzx8+BAmJiaYOHEi92BZ\nJpNh+fLlUFVVrTXkxTAfEgv8TKOTnJwMPp+PK1eucOdkMhmGDRuGMWPGAKjZNrF169bYuHEjZDIZ\nevfujblz56KiogIGBgbYv38/7t69C0VFReTl5SEqKgqtW7eGTCZDYWEhmjdvLjemfvDgQejo6ODU\nqVPIzc1Famoq9u7dix9++AGKiopwdHTEkSNHUF5ejp9//hmKiopYsGABjhw5AlVVVUydOhUbN24E\nn8/H+vXrMW3aNGhpaeHcuXPo27cv2rdvj8ePH2Py5MkwMDDAvXv38OzZMzg5OcHZ2VlueGvPnj0Q\nCAQ4ffr0x+t05ovCAj/TKIWHh0NDQ0MuM2dJSQkkEglCQ0MB1KR34PP5SE5ORnZ2NoRCIS5duoSz\nZ89CLBajuLgY06dPxw8//ACZTAZzc3MutbKdnR0OHz7M3Xvo0KHYvHlznW0pLy/Hjh07YG5uDjMz\nMxw6dAjp6elwdnaGqakpTp8+jf79+8PMzAxHjhyBkZERRo4cie3bt4PP52PXrl0ICAiAuro6rl27\nhnXr1kEkEuH8+fOorKzEmDFjYGlpKbeo6/Tp0xAIBNi3b19DdC/zhWOBn2m0ZsyYga5du8rNsb9x\n4wb4fD5u3LgBANiwYQPMzMxQVlaGffv2QVdXF8+fP8ekSZMwfPhwPHv2DGKxGGfPnsXRo0dhYGCA\nyspKbN68GT179uTuO2DAgForiP9KJpPh0KFDMDExQadOnXD16lXuYW1gYCB+++038Pl8bNy4EUOG\nDIG5uTmioqKgqamJ2bNnIywsDHw+H4cPH0Z0dDQX2GUyGebPnw8tLS2kpqZy9SUlJUFFRUUuVTXD\nfAgs8DONVnV1NZycnDBt2jS586GhodDV1cWTJ08gk8kwePBgbgjohx9+wIgRI1BaWgpdXV2Eh4cj\nMjISEokEL168gJOTE5YvX47y8nIoKytzq2f9/f0xc+bMerdrw4YNEAqFmDp1KlJTU9GtWze0b98e\nf/zxByQSCTw8PLB06VIoKSkhMjISHTt2xPfff4+YmBgoKytj1apVSEpKgqqqKjdjadOmTRCJRHJb\nO96/fx/a2toICgpic/2ZD4YFfqZRKy4uhp6eXq1hGC8vL/Tq1QtSqRQlJSUwMDDA5s2bUVpaCkND\nQ2zZsgXx8fEQCoXIzMzEoEGDuCCtqKiIjIwMrFmzBo6OjpDJZLh8+TJ0dHRqTfEMDQ2Fqqoqvvnm\nG0gkEkyYMAFxcXGQyWQoKCjAyJEjoampiRMnTmDJkiUQCoUICwtD3759YWNjw23VuHr1aowaNQqW\nlpaIj4+HoaEhvL29kZ6eDiMjI0yfPh1SqRSHDx8Gn8/HsWPHuDZkZ2fD1NQU06dPZ8Gf+SBY4Gca\nveTkZAgEApw5c4Y7V1lZic6dO3NZOu/cuQM+n4+EhATcunULfD4fSUlJWLx4Mdq3b4/c3Fyoqqri\n2LFjWLhwIRwdHVFRUQFDQ0NERkZCJpPBwsICR44c4ep4tQPY1atXUVZWhuvXr2Px4sXQ19eHhYUF\ntzn777//DrFYDG9vb8TGxkJdXR0zZ87E/PnzIRaLER4eDn19fUyZMgXz58+Huro6zp8/j06dOsHN\nzQ05OTmws7PD0KFDUVFRgbi4OAiFQrnx/eLiYtjY2GDs2LF17jbGMO+CBX7mk3D8+HEoKSnJ7beb\nn58PTU1N7Nq1CwCwf/9+qKmpIScnh9tMpaCgAD179oSPjw9Onz4NZWVlPHr0CNbW1ggODkZMTAxU\nVVXx5MkTREREwMLCgvvW//DhQ6iqqtZqy6tv5ubm5rC1tcWVK1dQWFiIvn37wsLCAvHx8ejatSuc\nnJywa9cuCAQCbNy4EV26dEGfPn2wefNmbnP4AQMGwMHBAbm5uXBxcUHPnj3x4sULJCUlQVlZWW6z\n+efPn8PBwQHu7u5sS0fmvbDAz3wy1q9fD4lEgqKiIu7cjRs3IBAIcPHiRQBAQEAA2rVrh/Lycnh7\ne8PJyQm5ubnQ1NTE3r17ERAQgM6dO3O/EJKSkjBhwgSMHDkSMpkM7dq149YQVFdXQ0FBgUv98FdS\nqRSbN2+GkpISvL29UVpaitWrV0MgEGD//v3w9vaGnp4eDh06BHV1dQQGBmL48OGwsbFBREQEBAIB\nwsPD4enpCUtLS2RlZWH48OHo0KEDnj59irS0NGhoaODXX3/l6iwrK0PPnj3Rt29fluKB+cdY4Gc+\nKT/++CM6duzIpWcGgKNHj0JZWRkPHjyAVCrFgAEDMGzYMFRWVsLR0RFTp07F1atXwefzkZiYiO7d\nu8PHxwe7du2Cnp4eHj9+DD09PezZs4fbGSwrKwsA4OnpiTlz5ry1TQUFBXB3d4dEIkFiYiLi4+Oh\nqqqKefPmISQkBEKhEOHh4bCwsMC4ceMwc+ZMSCQSrt0hISGYO3cu9PX18fDhQ3h5ecHS0hIFBQXI\nyMiArq4uFi5cyNX38uVLfP/99+jdu7dcPzBMfbHAz3xSpFIpBg4cCDc3N7kHsatXr4aBgQGKiorw\n4sULtG3bFoGBgSguLoa+vj7Wr1+PPXv2QENDA8nJydDW1saOHTvg6emJPn36ICEhAXw+H6mpqQgI\nCED37t0hlUpx//59KCoq1mvjlLCwMAgEAqxevRrZ2dmwtbWFm5sbfv/9dwgEAmzevBmOjo7o06cP\nli1bBlVVVfz+++/Q1NTEsmXLsGzZMmhqauLevXuYOXMmTExMkJubi+zsbBgaGsLf3597uFtZWYkB\nAwagR48eLPgz76xRB34iakFE4USUTES3icimjjIN1DVMY1VeXo6OHTvC29tb7vz06dO5XwM5OTnQ\n1NTE9u3bcffuXSgpKSE6Ohr+/v6wsbFBQkIC98DYzs4Oc+fO5ZKyFRUVoV27dvj5558BADNnzoSb\nm1u92nZsGMnxAAAgAElEQVTv3j2YmZlhxIgRePLkCQYPHox27drh7NmzUFVVxbJlyzB48GDY29tj\nw4YNEIlEiIqKgkQi4dYCqKqqIjk5GYGBgdDX10dWVhby8vJgYmKCOXPmcMG/qqoKgwYNQvfu3Vnw\nZ95JYw/824ho9J//bkpEzeso0zA9wzRqxcXFMDY2lttgRSqVws3NDf3790d1dTVu377NPUQ9d+4c\nBAIBEhMTMWTIEPTv3x+HDh3i5s1raGhg586dGDduHFxdXZGeng6RSISTJ0+irKwMBgYGcjtpvU1p\naSkGDhyI9u3bIy8vDzNmzICBgQEuXLgAiUQCPz8/eHh4wNraGqGhoRAKhTh69CiMjY0xZ84cbN68\nGWKxGMnJyVi4cCEkEgmysrKQn58PU1NT/PTTT7WCf8+ePdlevky9NdrAT0TNiOh+Pco1RL8wn4BH\njx5BQ0NDbpOVly9fwsHBAZ6enpDJZIiLi4NAIEBCQgIiIiKgoqKCO3fuoHPnzpg8eTLWrFkDiUSC\nM2fOQCAQ4MSJE+jcuTN8fHwQExMDoVCItLQ0JCUlgc/n4+bNm/Vqm1QqxaxZs6Cnp4f09HQsX74c\nGhoauHjxIszMzODj44Np06bBwsKCGyI6evQozMzMMHv2bGzbtg0qKipITk7GokWLIJFIkJ2djfz8\nfJiYmMDf35+rq6qqCgMGDEDv3r3ZA1+mXhpz4G9NRJeIaCsRXSWiECL6to5yDdU3zCcgOTkZIpGI\ny8ED1My/t7Cw4ILjwYMHIRKJkJKSgrVr10JXVxepqakwNTXFwoULMWvWLNjY2ODIkSMQCoWIi4uD\ngYEBVq5ciZCQEOjq6iI/Px+hoaHQ1tZGXl5evdsXHBwMNTU1pKSkICQkBGKxGPHx8bCyssK0adPg\n6+sLc3Nz7Nu3DwKBANHR0TAzM4Ofnx+2bdsGsViMtLQ0LFiwAIaGhsjLy0NeXh4MDQ25bSiBmjH/\nPn36oF+/fnIpLhimLo058FsRURURtfnzeCURBdZRDgEBAdzr9U22mS/DlStXIBAIcPLkSe5cXl4e\nJBIJVqxYAQDYsmUL1NXVkZGRgXnz5qF169a4c+cOtLS08Ntvv2H06NFwcnJCaGgoxGIxzpw5A1VV\nVezatQtz5sxBmzZtUFJSAj8/P1hbW8tl0/w7W7ZsgVgsRmpqKrZu3QqxWIyEhARYWlrC19cXPj4+\nsLKywr59+yAUCnHy5EkYGRlhwYIF2LhxI9TV1ZGeno65c+fCzMwMRUVFyM7Ohp6eHpfuAaj5tdOj\nRw8MGTKELfJi5MTExMjFycYc+JWI6MFrxx2I6Egd5Rqmp5hPyquhmgsXLnDnMjIyoK6ujo0bNwIA\nVqxYAT09PWRnZ2Pq1Klo164drl+/DrFYjO3bt6Nfv37o168fgoODoa2tjZMnT0JJSQkHDhyAh4cH\nOnfujNLSUowdOxZdunRBWVlZvdu3efNmqKmpISMjAyEhIdDQ0MD169dhamqKwMBAeHl5wc7ODmFh\nYVBSUkJsbCx0dXWxatUqrF69Gjo6OsjKyoKPjw+sra1RUlKCzMxMaGpqyu1dUFZWBnt7e3h4eLD0\nDswbNdrAX9M2OkNEkj//HUBEv9RRpmF6hvnkREdHQygUconXgJq9fFVUVLBz504AQGBgIIyNjZGX\nl4cxY8bA3t4eCQkJEIlECAsLQ48ePeDu7o7FixdDIpFwWTQPHz6MYcOGoVu3bigpKcHQoUPRrVs3\nbmOY+lixYgUMDQ1RXFyMZcuWwdDQEHfu3IGOjg7WrFmDESNGoGfPnti6dStUVVVx7tw5qKmpYfv2\n7Vi4cCFMTExQWFgIDw8PODg4oLy8HHfv3oWKigr27t3L1VNSUgJra2v4+vp+uM5lPiuNPfC3JqIE\nIkoiokgialFHmQbqGuZTdODAASgpKcntYHXr1i2IRCIuBfLs2bNhbm6O/Px8DBs2DI6Ojrh48SKU\nlJSwZ88edOvWDcOGDcP8+fNhYGCAo0ePQiAQ4ODBgxg2bBgcHBzw9OlTjBgxAh06dHinfXKnTZuG\nbt26oaqqCj4+PujYsSPu3LkDZWVlhIeHo3fv3hgxYgR+/fVX6Ovr4/z581BSUsLRo0cxffp0tG/f\nHiUlJRg4cCD69euH6upqXL9+HUKhUC6xW2FhIYyMjLgpqQzzukYd+OvVCBb4mb/Ys2cPRCIRbt26\nxZ1LSkqCkpISl1jNx8cHlpaWyMvLg7u7O5ycnHD+/HkIhULs3r0bXbt2xZAhQxAUFAQ9PT0cPnwY\nQqEQe/bswZgxY2Bra4uCggJMnToVJiYm9d4nt7q6Go6Ojpg5cyakUin69u2L0aNH4/Lly+Dz+Th7\n9iysra3h7++PGTNmoH379oiJiYFAIEB8fDyGDRsGV1dXlJaWomvXrpgwYQI3e4nP5+Py5ctcXY8f\nP4ampuYbN5hhvlws8DOfpZ07d0JZWRm3b9/mzl29elUu+E+fPh0WFhbIzc2Fu7s7unXrhnPnzkFJ\nSQmhoaHo3r07BgwYgMWLF0NLSwuHDx/mEqf9+OOPMDIywsOHD7Fs2TKoqKhw+YL+Tn5+PsRiMU6d\nOoXnz5/D1NQUa9euRWRkJLffrqamJkJDQzF48GC4ubnhwIEDUFZWRkpKChwdHeHp6YmnT5/C3Nyc\nm93zal3C64nsUlNTIRKJ5HYbYxgW+JnP1o4dO6CiolJn8H817OPr6wszMzPk5ORgxIgR6Ny5My5c\nuMDl0OnTpw969uyJ4OBgiMViHDp0CJqamli0aBGXdiExMRGHDx+GQCDAhg0b6vVQNSoqClpaWigr\nK8Pdu3e5xWWBgYGws7PjcgvFxcWhXbt28Pf3R3BwMIyMjJCZmQlTU1MsW7YM2dnZ0NDQ4BaXbdiw\ngZt++sqlS5fA5/Pr/cHEfP5Y4Gc+a6+++b++8CopKQkikQi7du2CTCbDnDlzYGhoiEePHsHDwwM2\nNja4ePEi1NTUsGLFCgwfPhx2dnbYtGkTBAIBIiIiYGZmhgkTJmDv3r3g8/mIiIhASkoKTExM4O7u\nXq9x//79+2Px4sUAgN27d0NfXx+lpaXo1asXfHx8EBkZCTU1Ndy6dQvq6uqIiIjApEmT0KNHD6Sn\np0MsFuPgwYO4deuW3H4FP/30E2xtbeVmHUVFRdX6NcB8uVjgZz57YWFhEIlEuHbtGnfu1q1bUFFR\n4fLdL1q0CDo6Onjw4AGmT58OU1NTXLp0CXp6evDz8+PG8l+tst28eTO6d++O7t27c5uvzJkzB8+f\nP8fEiROhrq4u97C1Lnfu3IFQKORSLbi5ucHX1xeFhYVQU1NDdHQ0Zs2aBUdHR+4ZwPXr19GlSxfM\nmDGD+yZ/48YNbr+Ce/fuQSaTwd3dHQMHDpRLZPfq18DrG9gzXyYW+JkvQkREBIRCIeLj47lzqamp\nUFdXx6pVqwDUZPh89Q17wYIF0NbWxqVLl2BpaYmxY8di8eLFUFNTQ0REBNTU1DB//nxMnDiRy8Pj\n4OCAbt26ITc3F8eOHYOmpibc3d3x+PHjN7arY8eOiIqKAgDk5uaCz+fj9u3bOH36NMRiMfLy8tCx\nY0csWrQIW7ZsgYGBAdLT06GpqYl9+/Zh586d0NLSQlFREdauXQsjIyM8e/YM5eXlaN++Pfz8/OTq\nmzVrFuzs7Fheny8cC/zMF+PVtMzY2Fju3MOHD6Grq4vAwEDIZDKEhoZCSUkJ8fHxWL9+PZSVlXH2\n7Fk4Ojqid+/e2LJlCwQCAXbv3g0rKysMHz4ca9as4aZ7zpkzByoqKjh27BhKS0vx008/QUFBATNm\nzEBubm6tNnl7e2Pp0qXc8fLly+Hq6gqgZl/h0aNHIzMzEwKBAFeuXMGoUaMwatQoXLlyBXw+Hykp\nKfDx8YGTkxOqq6sxfvx4uLq6QiqVIi8vDxoaGggLC+Pu//p+BWyB15eLBX7mi3Ly5EkIBAL8/vvv\n3LmcnByYmZlhypQpkEqlOHLkCPh8PqKjoxEZGQk+n4+jR49i5MiRaNOmDbdWYOXKlXBzc4O1tTUO\nHToEVVVV/PTTTzh+/DhUVVXh5eWF0tJSZGZmwtPTEy1btsTgwYOxd+9e3L9/H6mpqWjTpo3c3rpl\nZWUQiUS4efMmSkpKIBaLcf78eezcuRPGxsYoKiqCRCJBWFgY1q9fDzMzMzx//hz29vbw9/dHRUUF\n7OzsEBgYCABccrnXF7W9ePECbdq0kdvchfmysMDPfHEuXrzIzcl/5cmTJ+jQoQOGDBmCiooKbuHU\ntm3buCmeGzduRGBgIDQ0NBAVFQVDQ0N4enoiKCgIysrKOHjwILp164aOHTvixo0bGD58OLS1tXH8\n+HEA4IZjevXqBXV1dWhoaMDX17dWXp25c+di2rRpAGpmJtnY2EAqlcLFxQWBgYHcXgKPHz/GwIED\n4eXlhZycHKioqOD48ePIzs6GiooKoqOjAQD79u2DhoaG3Nh+VlYWVFVV5ZLbMV8OFviZL9KNGzeg\noqKCdevWcefKysrg6uoKJycnlJSUIDk5GZqampg/fz5SUlKgo6OD2bNnY+fOnRAIBAgLC0OvXr3Q\nqVMnLsfOzz//jIULF0IoFGLv3r2IioqCpqYmBg4ciIcPH9arbUlJSZBIJABqhmZMTU1x9OhRZGRk\nQFFREQ8ePMDcuXPh6uqKoqIiqKmp4Y8//sCpU6egrKyM3NxcnDlzBkpKStx+wTNnzkTXrl3lPmRe\nPTB+faEb82VggZ/5Yt2/fx86OjqYN2+e3MYmHh4esLCwQHZ2NrKzs2FlZYXRo0fj8ePHsLOzQ//+\n/RETEwOxWIygoCDMmTMHampqOHDgAGxtbdGrVy8cO3YMEokEgwcPRkZGBgICAqCgoIAff/xRbo59\nXaqqqtCkSRMuSO/atQsODg4AgPnz52PQoEF4+fIlDA0NERERgZMnT0JVVRVPnjzB7Nmz4ezsDJlM\nhsWLF8POzg5VVVWoqqpCly5dau0fvH37dujq6qK4uLgBephprFjgZ75oOTk5MDc3h6enJxdoZTIZ\n5s+fD01NTdy5cwelpaVwcXGBg4MDt9DLwsICly9fhrW1NQYMGIA9e/ZAIBBgxYoV8PX1hVgsxuHD\nh+Ht7c2tGXj8+DE8PT3RqlUreHl5ITk5uc42FRYW4n//+x/3YVRRUQE+n4+HDx+itLQUIpEISUlJ\n3DTSFy9eYMKECRg7diwqKipgZWWF9evXQyqVwtHREXPnzgVQk6paLBZzQ0CvTJ48Gb169ZKb+sl8\n3ljgZ754z549Q9euXfH999/LLXravn07hEIhTp8+jerqanh7e0MikSAlJQVLliyBsrIyTp48ibFj\nx8LIyAjHjh2DpaUl+vbty6VfmDx5MmJjY2FmZoauXbvi9u3byMrKwk8//QSRSIS2bdti/vz5OH78\nOO7cuYO4uDj07NkTHh4ecm0cOXIkNyy1dOlSuLu7AwAGDhyIoKAgPH36FGKxGHFxcbh9+zb4fD7S\n09ORnZ0NJSUlnD9/HgAQGxsLkUiErKws7t6VlZWws7NDUFBQQ3c100iwwM8wqPlWPWTIELRr107u\nIejp06chFAqxZcsWAEBISAi3j++rNNBr165FSEgI+Hw+QkNDMWXKFKirq+Po0aMYMmQIdHV1cerU\nKaxcuRJ8Ph+TJk1CXl4eqqqqcOLECfj4+MDe3h76+vrcB8Ff59mvXbuW+zB4+vQpWrZsiby8PNy7\ndw+KioooKirC3r17YWZmhurqaixatAjdu3eHTCbDgQMHoKOjw20eExgYiC5dush9w8/KyuI+yJjP\nHwv8DPMnqVSKmTNnQk9PTy61QXJyMnR0dDBjxgxIpVLuW/OKFSuQlpYGExMTjBw5EufPn4eOjg7G\njx+P/fv3QyQSYdasWdi3bx9UVFTg4eGBu3fvYsqUKVBQUMDMmTPrvY3jzp07MXjwYO546NChWLNm\nDQDghx9+4J5TdOrUCSEhIaisrISZmRk3h3/EiBHw8vICUJMhtEOHDnIb1QM1U12VlZWRnZ39Xv3I\nNH4s8DPMX6xfvx4ikQjnzp3jzhUWFqJz585wdXVFSUkJ0tPTYW5ujmHDhiE/Px9Dhw6FmZkZEhMT\nMWjQIJiYmCA2Nhaurq4wMzPD2bNn4enpCZFIhC1btiA9PZ2b2z9mzBhcvnz5rQuqvL29ubF6oCb1\ntIuLCwD51A+XL1+GWCxGWVkZzp8/D7FYjJKSEhQXF0NFRQVxcXEAgPT0dC4FxOv8/f1rzf5hPj8s\n8DNMHV7tvLV7927uXEVFBcaNGwdjY2Pcu3cPL168wNChQ9G6dWvcu3cP69atA5/Px+7du7F582bw\n+XysWrUKW7duhUAggJ+fH+Li4mBjY4M2bdogNjYWeXl5WLhwIbS1taGvr4+ZM2ciOjoajx49wsuX\nL5GdnY2VK1dCKBTK5fx/8OAB1NTUuGN7e3uEh4cDAFxdXREcHAwAGDZsGDeTJzw8HIaGhqioqABQ\nsxdw69atuWOgZkZRx44dueRxzOeJBX6GeYPr169DQ0MDAQEB3LdxmUyGNWvWcOP8MpkMq1atglAo\nxJEjR3D16lXo6upi7NixSEpKgrW1NRwdHZGQkIC+fftCIpHgxIkT2LVrFzQ1NdGjRw/u2/6lS5cw\nd+5cODg4QElJCV9//TX4fD769esnl1oaAJ4/f45vv/2WO964cSMGDRoEAIiPj4eGhgaqqqqQmZkJ\nBQUFZGVlQSaTwdnZmduV69Xxq1W+r7xKEZGQkNCQ3cv8i1jgZ5i3yMnJgY2NDdzc3OT21301zv/L\nL79AJpPh/vTpMFdRwYwZM1BUVIThw4fDwMAAly9fRlBQEAwUFHBh6FAcOHAA6urqcHd3x/3797F2\n7VqoqqrCyckJx44dq3f+nLt378p943/06BEUFRW599va2uLQoUMAgOnTp3Pj+3fv3oWioiI3jp+Z\nmQk+n487d+7I3X/Pnj2QSCTvtKcw8+lggZ9h/kZ5eTmGDh0KKysrPHr0iDufmZmJtm3bYqO5OUCE\nKn19DHJwQPv27ZGRkcGt8A3280OZtjZAhGADA9y8eRNz5syBgoICAgICUFhYiC1btsDU1BQSiQS/\n/PLL327lOG3aNEyYMEHunJKSEpcJdOPGjejXrx+AmqyfrVq14pLE+fr6YsyYMdz71qxZg44dO9b6\n0BkyZAimTJnyzzuOabQadeAnoodEdJ2IrhHR5TeUaZieYZjXyGQy/Pzzz1BRUeHmxAPAy5cv4TN8\nONK++QYggszICMF+fhAKhYiMjETmlSt48N//AkR4qaOD1XPnQlFREYsXL0Zqairc3d2hrKyM4OBg\nlJWV4dy5cxg7diwUFRVhYWEBHx8f7Nq1C+fPn0dCQgIiIyMxcOBA6Ovr11oB3Lp1a1y9ehVATV6g\nZs2acesSJkyYAH9/fwA1eYkEAgH3Lb+6uhpWVlbYvn273P2KioogFovlspkyn4fGHvgfEFGrvynT\nEP3CMHWKioqCQCBASEiI3Pk9wcFIbtKk5n8LIyMkREWhrYYGslq1AohQqKQE/VatsHLlSqSlpcHZ\n2RkGBgY4fvw4rl69il69ekFVVRUrV67E8+fPUVVVhbi4OCxYsAADBw6EjY0NLCws4OzsjBUrVnBz\n8l+no6ODlJQU7tjW1hanT58GULPxjLKyMqqqqgAAixcvlpseGh8fDxUVFZSUlMjd8/Dhw9DR0WFD\nPp+Zxh7404lI8W/KNES/MMwbpaamwtDQEOPGjZNbaHUnNhZ3X33z5/Mh5fMBIqR98w0So6ORlpaG\njh07ol27drh58yYOHjwIbW1tuLi4ICUlBQkJCejfvz8UFRXx448/Ii0trd5tevDgARQUFLjADtTk\n81+xYgV3bG1tjT/++AMAUFJSAj6fL1fH8OHDa+XyAYDBgwfD19f3nfqIadzeJ/B/RfXA4/HUXvu3\nsD7veQ2I6BiPx0vg8Xjj3vG9DNMgJBIJXbp0iYqKiqhTp0706NEjIiIy7NyZlFNSqOQ//yFeYSF9\nVVhIJBBQ8tq11HPkSAoNDaXjx4/T8OHDycHBgS5fvkwJCQnUoUMHsrOzo02bNlFwcDBdvnyZiIjs\n7OyoQ4cOtHr1anr48OEb2/PkyRMaNWoUeXt7U9OmTbnzampqlJWVxR27ublRZGQkERE1a9aMxo8f\nT6tWreKuL1iwgH777TfKycmRu//KlStp27ZtdP369ffuO+bTx6v54PibQjzebiIaBaDyz8DfGUB4\nvSrg8UQAcnk8noCIThCRF4BzfymDgIAA7tje3p7s7e3f4c9gmH8GAC1dupR+/fVXCg0NJUdHR6L8\nfCITE6KCAiIiKvvuO/r23j3KlcnIw8ODHj16RFu3biUlJSXy9vamy5cv08qVK8nOzo5++eUX2rx5\nM40YMYJ8fX1JKBTSsWPHKDw8nKKjo6lZs2bUtm1b0tPTIz6fT1VVVZSSkkKHDh2ioUOH0vLly6lJ\nkyZc+5YsWUIFBQW0dOlSIiK6ffs2ubi40IMHD4iIKCsri0xNTSkzM5O+++47IiLy8fGhyspKWr16\ntdzfGhISQtu3b6e4uDj66qt6fedjGpHY2FiKjY3ljgMDAwkA7x/drD4/C4joh78c9/4nPy+IKICI\nptdx/sP+BmKYd3T69GkoKytjqa8vZEZGNaOgAgGqFBQAIjz87jvk3bwJmUyGbdu2QSAQ4KeffkJ5\neTlOnDgBAwMDODk54datW8jOzoa3tzdatWqFH374ATdv3gRQk07i9u3bCA0NRUBAALy8vODt7Y01\na9bg/v37dbbr9cRur+7RvHlzFBYWcudcXFywbds27jgvLw+tWrWSS+L26r1t27at9QCY+TRRQ4/x\nE5EzEe0hIhciMiOiWfV833+J6Ls///0/IjpPRE51lGvI/mGYesm5fp2bvVOlrw/k5QF5eZAZGgJE\nSGnSBEf/TPSWnZ2Nfv36QSKRICYmBpWVlVi5ciUEAgE8PDyQnZ2NgoICzJ8/H8rKyujUqRN27NiB\n0tLS+rcnJwcKCgrcRiyvvEon/UpYWBicnZ3lykyZMqXOMf2LFy9CLBbX+WCZ+bQ0eOCvqYP0iGgh\nEf1CRNr1fI8WESVRzVTOm2/6wGCBn/nX5eUBf37Tz+XzYSYScbNpXr+W9vXXmDhgALfpycGDB6Gm\npoYRI0YgLy8PRUVF8PHx4RK4FRQUoLKyEuHh4ejRowdatGgBNzc37Ny5Ezk5OW9sTnJyMszNzWut\nyAVqUjucOnWKO3769Cm+++47uXTU6enpUFBQqDPAu7u7IyAg4B92FNNYfJTA35AvFviZf92aNXg1\njRN5eTh27BiUlZUxd+7cmlk2rwX/vZ06QVVVFUePHgVQM7vGx8eHy+vzKs3C+PHjoaCggBkzZnBB\nPi8vDyEhIfj+++/RsmVLaGpqolevXvD09ISPjw88PDzQrl078Pl8BAcH17kK2NzcvFYqBhsbm1pz\n9V1dXbFx48Za7381e6i+WUWZxokFfob5ENasqQnwf8rJyYGjoyPat2+P9PT0mmt/plE+ffo0tLS0\nMGzYMC73/+3bt9GtWzcYGRlxO2Q9fPgQkyZN4jJ43rhxg7u/VCpFSkoKDhw4gNWrV2PJkiVYt24d\nTp06hfLy8jqb+OzZM3z33Xe1howmTpyIlStXyp07fPgw2rdvX+d9Xj1fYD5dLPAzTAORSqVYunQp\nBAIBdu7cKXettLQU06ZNg0gkws6dOyGTySCTyXDw4EHo6enByckJ165dAwDk5+cjKCgIKioqsLOz\nw9atW2sttKqPefPmcYncXrd06dJagbyyshKKiop1po7Izs6WSwHBfHpY4GeYBnb16lUYGhpi8ODB\ntTY1v3TpEszMzODo6MgtpqqsrMTq1ashEokwePBgbjVuZWUlIiMj4eLigubNm6Nfv37Yvn37W8f7\ngZqUDKtWrYKysnKth71Azf4D48aNq3V+5MiR3GYvfzVp0iTMmjWrXn8/0/i8T+Bnk3kZph4sLCwo\nMTGRBAIBtW7dmk6cOMFds7a2pitXrpCTkxO1a9eO/P39qaqqiry8vOju3btkampKHTp0oGHDhtHd\nu3epb9++dPjwYUpPT6fevXvToUOHyNDQkCQSCbm7u9O8efMoJCSEduzYQevWraMpU6aQnp4ehYeH\n09mzZ0lNTa1W+54/f07/+9//ap13cnKikydP1vk3+fj40MaNG6m0tPTDdRTzSajXAq4GbwSPh8bQ\nDoapjxMnTtCYMWOoV69etGTJEmrWrBl37fHjx+Tj40OXLl2iZcuWUf/+/YnH41FJSQmtWbOGVq1a\nRba2tuTj40MdO3YkHq9m/Y1UKqU7d+7QtWvX6O7du5Sbm0vl5eX03XffkY6ODnXt2pUsLCy48n81\nfPhw6tixI3l4eMidz8jIIFtb21oreV/p378/devWjTw9PT9Q7zAfC4/HIzTkAq6GfhEb6mE+MU+f\nPsWoUaOgqalZ5+bmMTExMDMzQ6dOneRm4Lx48QLr1q2DRCKBubk5QkJC/tFY/+uKi4uhoKAgl276\nFZlMhmbNmtUannrl5MmTMDMzq/ceAkzjQWyMn2H+HVFRUVBTU8O4cePw9OlTuWvV1dUICQmBsrIy\nt2nLK1KpFH/88Qe+//57tGjRAiNGjEB0dLTcFor1UV1djYEDB2LixIlvLGNkZCQ3m+h1UqkUmpqa\nXCpo5tPxPoGfjfEzzHtwdnammzdv0ldffUXGxsZ04MAB7lqTJk1o3LhxlJaWRvr6+tS2bVuaNGkS\nZWVl0VdffUXdu3enAwcOUGpqKllaWlJgYCCJRCIaNGgQbdq0iVJSUkgmk9VZLwCKj4+nrl270rNn\nz2jZsmVvbOP//vc/Ki8vr/PaV199RUOGDKGwsLD36wjmk8LG+BnmAzl79ix5eHiQgYEBrV69utZD\n2IKCAvrll19oy5YtNGzYMJoxYwapqqrKlcnJyaHo6Gg6deoUnT9/noqLi0lfX5/U1dWpRYsWBIAK\nCstOu+kAAAw+SURBVAro+vXr9H//93/k7e1NEyZMkEvs9leGhoYUERFBxsbGdV6/evUqDRo0iO7e\nvfv+ncB8NO8zxs8CP8N8QBUVFfTzzz/T6tWradasWTR16lT6+uuv5crk5ubSsmXLaMuWLdS3b1/6\n8ccfydDQsM77FRcXU9r/t3e3MVJVdxzHf3+k7i6xmAAJooAIEQlkXSgIEoqOgIA2VG2o7KbNNvjO\nWCRQjaKVB5P2DbEGWErS1G3QpOJiENdKiyYyERbFXR5kExC38gIWgkTALRSXJ/99sSOF2YeZgb2z\nM3u+n2TCHe7hnrMnh9+eOffeuV9+qcbGRjU1NcnM1LdvX40aNUrDhg1r92Tvle3p06ePjh8/3uZV\nP1LLp4dbb71VNTU1Gjp06LX94Mg6gh/IMQ0NDZo3b54OHz6sVatWacqUKa3KnDhxQqtXr9bq1as1\nduxYzZs3TzNmzOjUr0yurq7W8uXLtXXr1g7LlZWVafr06Zo7d26n1Y1oEfxADnJ3bdy4UQsXLtTY\nsWO1fPly3XHHHa3KNTc3680331RFRYVOnTqlJ554QuXl5Ro8ePB11X/hwgVNmDBBzz33nObMmdNh\n2VdffVVfffWVKioqrqtOZM/1BD8nd4GImJkee+wx7du3TyUlJRo3bpwWLVqkpqamq8oVFhZq7ty5\n2rlzp6qqqnTkyBGNGTNG9913nyoqKi4/HSwTzc3NKi8v16BBg/T444+nLD9y5Ejt378/43qQnwh+\nIGJFRUV66aWXtHfvXh07dkx33XWXVq1apfPnz7cqO27cOK1Zs0ZHjx7VM888o9raWo0ZM0bFxcWa\nP3++1q9fr4MHD7Z7tc/JkydVWVmpu+++W+6udevWpTwPILV+xCO6N5Z6gCz7/PPP9fzzz+vAgQN6\n+eWXVVZW1uFVOZcuXVJdXZ22bNmi7du3a/fu3Tp58qQGDx6sfv36qaioSM3NzTpy5IiOHz+uqVOn\n6qmnnmp5jGSaDh8+rIkTJ6qxsbEzfkRkAWv8QB6Kx+N68cUXderUKS1ZskSzZ8/u8BfAlU6fPq1D\nhw7pm2++UXNzswoLCzVgwAANHTr0qoe1p+vgwYOaMmVKhw+ER24h+IE85e7avHmzli1bpqamJi1a\ntEilpaWtLgGN2rZt2/Tss8/qk08+yWq9uHac3AXylJlp5syZ2r59u1asWKHKykrdeeedWrFihU6f\nPp21dtTW1qqkpCRr9aFrEfxADjAzPfjgg9qyZYuqqqpUU1OjIUOGaMGCBZHfUfvDSeBZs2ZFWg9y\nR1aC38x6mNkuM6vORn1APhs/fryqqqq0a9cuFRQUaNKkSZo2bZrWrVvX7nfuXI/q6mo1NTVp5syZ\nnX5s5KasrPGb2QJJYyX1dveft7GfNX6gHefOndM777yjyspK1dXV6dFHH1VpaakeeOCB6z4XsHfv\nXk2fPl3r16/X5MmTO6nFyIacXuM3s4GSHpb016jrArqjgoIClZaW6oMPPlB9fb1GjRqlxYsXq3//\n/iorK9PatWszvsnr/PnzWrNmjaZOnaqVK1cS+oGJfMZvZusl/UHSzZJ+x4wf6BxHjx7Vpk2btHnz\nZsXjcfXq1Uv33HOPSkpKNHz4cA0aNEj9+vVTYWGhLl68qBMnTqihoUE1NTXasGGDiouL9corr6i4\nuLirfxRcg5y9nNPMfibpIXf/rZnF1BL8rc4gmZkvWbLk8vtYLKZYLBZZu4Duxt3V0NCguro61dfX\nq6GhQY2NjZev8+/Zs6f69OmjYcOGafz48Zo1a5ZGjBjR1c1GBuLxuOLx+OX3y5Yty9ng/6OkX0u6\nKKlI0o8lbXD38qRyzPgBIAM5O+O/qiKz+8VSDwB0ipw+uQsAyC18ZQMA5CFm/ACAtBH8ABAYgh8A\nAkPwA0BgCH4ACAzBDwCBIfgBIDAEPwAEhuAHgMAQ/AAQGIIfAAJD8ANAYAh+AAgMwQ8AgSH4ASAw\nBD8ABIbgB4DAEPwAEBiCHwAC0zPKg5tZgaSPJd2YqOttd18WZZ0AgI5F/rB1M+vl7mfN7AZJNZKe\ndvfPksrwsHUAyEBOP2zd3c8mNgvUMusn4QGgC0Ue/GbWw8x2Szom6UN3r426TgBA+yJd45ckd/9e\n0hgz6y1po5mNdPd9yeWWLl16eTsWiykWi0XdNADIG/F4XPF4vFOOFfka/1WVmS2WdMbd/5T096zx\nA0AGcnaN38z6mdnNie0iSdMkfRFlnQCAjkW91DNA0loz66GWXzJvufumiOsEAHQgq0s97TaCpR4A\nyEjOLvUAAHIPwQ8AgSH4ASAwBD8ABIbgB4DAEPwAEBiCHwACQ/ADQGAIfgAIDMEPAIEh+AEgMAQ/\nAASG4AeAwBD8ABAYgh8AAkPwA0BgCH4ACAzBDwCBifph6wPN7CMz22dm9Wb2dJT1AQBSi/SZu2Z2\ni6Rb3H2Pmd0kaaekR9z9i6RyPHMXADKQs8/cdfdj7r4nsX1G0n5Jt0VZJwCgY1lb4zezIZJGS9qR\nrToBAK1lJfgTyzxvS5qfmPkDALpIz6grMLOeagn9N9z93fbKLV269PJ2LBZTLBaLumkAkDfi8bji\n8XinHCvSk7uSZGavS/rG3Rd2UIaTuwCQges5uRv1VT2TJH0sqV6SJ14vuPu/ksoR/ACQgZwN/rQb\nQfADQEZy9nJOAEDuIfgBIDAEPwAEhuAHgMAQ/AAQGIIfAAJD8ANAYAh+AAgMwQ8AgSH4ASAwBD8A\nBIbgB4DAEPwAEBiCHwACQ/ADQGAIfgAIDMEPAIEh+AEgMAQ/AAQm0uA3s9fM7Gsz2xtlPQCA9EU9\n4/+bpBkR14ErxOPxrm5Ct0J/di76MzdEGvzuvk3SqSjrwNX4j9W56M/ORX/mBtb4ASAwBD8ABMbc\nPdoKzG6X9J67391BmWgbAQDdkLvbtfy7np3dkDZY4tWua208ACBzUV/O+XdJ2yUNN7NDZjY3yvoA\nAKlFvtQDAMgtWTu5m87NXGa20swazGyPmY3OVtvyTaq+NLP7zexbM9uVeP0+223MJ2Y20Mw+MrN9\nZlZvZk+3U47xmUI6fcn4TJ+ZFZjZDjPbnejPJW2UudHM1iXG5idmNjjlgd09Ky9JP5U0WtLedvY/\nJOn9xPYESZ9mq2359kqjL++XVN3V7cyXl6RbJI1ObN8k6YCkEUllGJ+d15eMz8z6tFfizxskfSpp\nfNL+JyX9ObE9R9K6VMfM2ozfU9/M9Yik1xNld0i62cz6Z6Nt+SaNvpRSnFDH/7n7MXffk9g+I2m/\npNuSijE+05BmX0qMz7S5+9nEZoFaLshJXp9/RNLaxPbbkqamOmYuXcd/m6TDV7w/orYHDNJzb+Lj\n4ftmNrKrG5MvzGyIWj5N7UjaxfjMUAd9KTE+02ZmPcxst6Rjkj5099qkIpfHprtfkvStmfXp6JjZ\nuJwzXW3NADjzfG12Srrd3c+a2UOSNkoa3sVtynlmdpNaZkzzE7PVq3a38U8Yn+1I0ZeMzwy4+/eS\nxphZb0kbzWyku++7okjy2DSlGJu5NONvlDToivcDJR3torbkNXc/88PHQ3f/p6QfpZoBhM7Meqol\nqN5w93fbKML4TFOqvmR8Xht3/4+kuKSZSbsOKzE2zewGSb3dvcOl4GwHf0c3c1VLKpckM7tX0rfu\n/nW2GpaH2u3LK9eezWy8Wi7bPZmthuWpSkn73H1FO/sZn+nrsC8Zn+kzs35mdnNiu0jSNElfJBV7\nT9JvEtu/lPRRquNmbakncTNXTFJfMzskaYmkGyW5u//F3TeZ2cNm9m9J/5XEzV7tSNWXkmab2ZOS\nLkj6Ti1n+tEOM5sk6VeS6hNrqS7pBUm3i/GZkXT6UozPTAyQtNbMeqhlov5WYiwuk1Tr7v+Q9Jqk\nN8ysQdIJSaWpDsoNXAAQmFxa4wcAZAHBDwCBIfgBIDAEPwAEhuAHgMAQ/AAQGIIfAAJD8ANAYAh+\nAAgMwQ+0wcxmmFny198C3QLBD7Rtj1q+Phjodgh+oG0TJW3v6kYAUSD4gbZNlPSZmf3CzOoS3zEP\ndAsEP9C20ZJ+4u4bJE1294td3SCgsxD8QJLEU4zOSRpoZuXu/l1XtwnoTAQ/0FqJpDpJVZLuSTwX\nFug2CH6gtVGStko6rpaZPzN+dCs8gQsAAsOMHwACQ/ADQGAIfgAIDMEPAIEh+AEgMAQ/AASG4AeA\nwBD8ABCY/wH0dew5Kx3OKQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "bpplt.pyplot.figure();\n", "bpplt.contour(B, np.linspace(1,3,1000), np.linspace(1,9,1000),\n", " n=10, colors='k');\n", "bpplt.plot(c, x=k, color='r', marker='x', linestyle='None',\n", " markersize=10, markeredgewidth=2)\n", "bpplt.pyplot.xlabel(r'$k$');\n", "bpplt.pyplot.ylabel(r'$c$');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, the true parameters are captured well by the posterior distribution.\n", "\n", "## Improving accuracy\n", "\n", "The model can be improved by not factorizing between B and tau but learning their joint posterior distribution. This requires a slight modification to the model by using GaussianGammaISO node:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from bayespy.nodes import GaussianGamma\n", "B_tau = GaussianGamma(np.zeros(2), 1e-6*np.identity(2), 1e-3, 1e-3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This node contains both the regression parameter vector and the noise parameter. We compute the dot product similarly as before:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "F_tau = SumMultiply('i,i', B_tau, X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, Y is constructed as follows:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Y = GaussianARD(F_tau, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because the noise parameter is already in F_tau we can give a constant one as the second argument. The total noise parameter for Y is the product of the noise parameter in F_tau and one. Now, inference is run similarly as before:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 1: loglike=-4.130423e+01 (0.002 seconds)\n", "Iteration 2: loglike=-4.130423e+01 (0.002 seconds)\n", "Converged at iteration 2.\n" ] } ], "source": [ "Y.observe(y)\n", "Q = VB(Y, B_tau)\n", "Q.update(repeat=1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the method converges immediately. This happens because there is only one unobserved stochastic node so there is no need for iteration and the result is actually the exact true posterior distribution, not an approximation. Currently, the main drawback of using this approach is that BayesPy does not yet contain any plotting utilities for nodes that contain both Gaussian and gamma variables jointly.\n", "\n", "## Further extensions\n", "\n", "The approach discussed in this example can easily be extended to non-linear regression and multivariate regression. For non-linear regression, the inputs are first transformed by some known non-linear functions and then linear regression is applied to this transformed data. For multivariate regression, X and B are concatenated appropriately: If there are more regressors, add more columns to both X and B. If there are more output dimensions, add plates to B." ] } ], "metadata": { "celltoolbar": "Edit Metadata", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }