Skip to content

Instantly share code, notes, and snippets.

@mlaves
Last active March 8, 2024 05:37
Show Gist options
  • Save mlaves/607d5252325d44fcea02d42179811d2e to your computer and use it in GitHub Desktop.
Save mlaves/607d5252325d44fcea02d42179811d2e to your computer and use it in GitHub Desktop.
Simple Classification Example in Pyro with SVI and MCMC
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simple Classification Example in Pyro with SVI and MCMC"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example, we train a simple Bayesian neural network on the task of binary classification using Pyro.\n",
"A sample of our toy data consists of two features $ X_{i} = [x_{i}, y_{i}]^{T} $ and a class label $ c_{i} $.\n",
"\n",
"We construct a linear model $ z = wX + b $ producing logits and define a categorical distribution for our final output\n",
"\n",
"$$ p(c | z) = \\mathrm{Categorical}(c | z) . $$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import torch\n",
"from torch import nn\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"from torch.distributions import constraints\n",
"\n",
"import pyro\n",
"import pyro.distributions as dist\n",
"import pyro.optim as optim\n",
"\n",
"pyro.set_rng_seed(1)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# generate toy data\n",
"num_samples = 200\n",
"data_1 = np.random.multivariate_normal((1,1), np.array([[1,0.5],[0.5,1]]), size=num_samples)\n",
"data_2 = np.random.multivariate_normal((5,1), np.array([[2,0.5],[0.5,2]]), size=num_samples)\n",
"data = torch.FloatTensor(np.concatenate([data_1, data_2], axis=0))\n",
"\n",
"y_1 = np.repeat([0], num_samples, axis=0)\n",
"y_2 = np.repeat([1], num_samples, axis=0)\n",
"y = torch.FloatTensor(np.concatenate([y_1, y_2], axis=0))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f890ccf9d90>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2de5Qc1X3nv7/p6UE9wtZIoMSrGRHJgRVGGFA0ZHHw8S4oQTYPWYY9wmadhOBj1hv8wOYIpHWCZU6yFtYuEK+dzRJwcs6alwLyIENAECB/mBwwkkePCFBwjEEz4BMhMbKNWkzPzN0/amqmuvreqluP7qru+n7O4YiprseveqTf797fU5RSIIQQUjy6shaAEEJINtAAEEJIQaEBIISQgkIDQAghBYUGgBBCCkp31gJE4eSTT1ZLlizJWgxCCGkrdu3a9ZZSaqH/eFsZgCVLlmDnzp1Zi0EIIW2FiLymO04XECGEFBQaAEIIKSg0AIQQUlDaKgZACCFZUKvVMDIyguPHj2ctSiBz5szBwMAAyuWy1fk0AIQQEsLIyAje8573YMmSJRCRrMXRopTC4cOHMTIygqVLl1pdQxcQIYSEcPz4cZx00km5Vf4AICI46aSTIu1SaAAIIcSCPCt/l6gy0gAQYsPercDtZwKb+pw/927NWiJCEkMDQEgYe7cCP/gicPQgAOX8+YMv0giQlvP4449j2bJlOPXUU7F58+bE96MBICSMp24BatX6Y7Wqc5yQFjE5OYnrrrsOjz32GF588UXcd999ePHFFxPdk1lAhIRxdCTacVJ4hoZHsWXHAbwxVsWivgrWr16GtSv6E93zRz/6EU499VS8//3vBwB88pOfxMMPP4wzzjgj9j25AyAkjHkD0Y6TQjM0PIqN2/ZhdKwKBWB0rIqN2/ZhaHg00X1HR0exePHimZ8HBgYwOprsnpkaABHpE5EHReRlEXlJRD6UpTyEaFl1M1Cu1B8rV5zjhPjYsuMAqrXJumPV2iS27DiQ6L66+e1JM5OydgH9BYDHlVL/WUR6APRmLA8hjZy1zvnzqVsct8+8AUf5u8cJ8fDGWDXScVsGBgZw8ODBmZ9HRkawaNGiRPfMzACIyHsBfATA1QCglBoHMJ6VPIQEctY6KnxixaK+CkY1yn5RX0Vztj3nnnsuXnnlFbz66qvo7+/H/fffj3vvvTfRPbN0Ab0fwCEAfyMiwyJyl4jM9Z8kIteKyE4R2Xno0KHWS0kIIRFYv3oZKuVS3bFKuYT1q5clum93dze+/e1vY/Xq1fjABz6AdevWYfny5cnumejqZHQD+C0AX1BKPS8ifwFgA4A/9Z6klLoTwJ0AMDg42OgEI4SQHOFm+6SdBQQAF198MS6++OLE93HJ0gCMABhRSj0//fODcAwAIYS0NWtX9Kei8JtNZi4gpdTPARwUEXdftApAsqoGQggh1mSdBfQFAPdMZwD9FMAfZSwPIYQUhkwNgFJqN4DBLGUghJCiwkpgQggpKDQAhBBSUGgACCGkTbjmmmvwa7/2azjzzDNTuR8NACGEtAlXX301Hn/88dTuRwNACCFp06QJch/5yEewYMGCVO4FZJ8GSgghnYU7Qc4dIuROkANy10+KOwBCCEmTNpogRwNACCFp0kYT5GgACCEkTdpoghwNACGEpEkTJ8h96lOfwoc+9CEcOHAAAwMDuPvuuxPdj0FgQghJkyZOkLvvvvsS38MLDQAhhKRNm0yQowuIEEIKCg0AIYRYoFT+BxJGlZEGgBBCQpgzZw4OHz6cayOglMLhw4cxZ84c62sYAyCEkBAGBgYwMjKCQ4cOZS1KIHPmzMHAgH26KQ0AIYSEUC6XsXTp0qzFSB26gAghpKDQABBCSEGhASCEkIJCA0AIIQWFBoAQQgoKDQAhWdGkqVGE2MI0UEKyoI2mRpHOhQaAtD1Dw6PYsuMA3hirYlFfBetXL8PaFf1ZixVM0NQoGgDSImgASFszNDyKjdv2oVqbBACMjlWxcds+AMi3EWijqVGkc2EMgLQ1W3YcmFH+LtXaJLbsOJCRRJa00dQo0rnQAJC25o2xaqTjuaGJU6MIsYUGgLQ1i/oqkY7nhrPWAZd9C5i3GIA4f172Lfr/SUthDIC0NetXL6uLAQBApVzC+tXLMpTKkjaZGpU6e7c2ZVwiiQ4NAGlr3EBv22UBFRWmv+YKyfOAAz+Dg4Nq586dWYtBCInL7Wc6St/PvMXAl/+59fIUBBHZpZQa9B/PPAYgIiURGRaRR7KWhRDSZJj+misyNwAAvgTgpayFIIS0AKa/5opMDYCIDAC4BMBdWcpBCGkRTH/NFVnvAO4AcCOAKdMJInKtiOwUkZ15n8dJCAmB6a+5IrMsIBG5FMC/KaV2ich/Mp2nlLoTwJ2AEwRukXiEkGZR1PTXHJLlDuB8AGtE5GcA7gdwoYh8L0N5COls2H6a+MjMACilNiqlBpRSSwB8EsDTSqlPZyUPIR2Nm39/9CAANZt/TyNQaLKOARBC4hB1NR/UfpoUllwYAKXUPyqlLs1aDkLagjireWP+/cHoLqGsXEl0YaVOLgwAISQCcVbzgXn2EVxCWbmS6MJqCjQAhLQbcappdfn3fmxcQlm5kujCago0AIS0G1Grad3um7UqIKXge4e1ZMiqlQNbSDQFGgBC2o0o1bR1rhMAatI5t7JAf++wlgxZtXJgC4mmQANASLsRpZrW5DoB4rVkyKqVA1tINAXOAyCkHbGtpjW5SKpvA5ffGX0wi/t52gNdwobENOu5BYfzAEgDQ8OjHLDSKbRD/33/kBjAWd2zR1Bq5HYeAMkXQ8Oj2LhtH0bHqlAARseq2LhtH4aGR1v2/PM3P42lGx7F+ZufbtlzO5Z2cJ0wwyczaABIHVt2HKibrwsA1doktuw40PRnZ218OpJ26L4ZJcOHxWCpwhgAqeONsWqk416Suo6CjA9dUAnIe/fNeQMGN5Uvw4fzhFOHOwBSx6I+fbGQ6bhLGqt3k5EZHavmzy3ElWh62Lqp0nQV8fcHgAaA+Fi/ehkq5fpioUq5hPWrlwVel4bryGRkBMiXW4htCdLF1k2VVjEYf38z0ACQOtau6Mc3Lv8g+vsqEAD9fRV84/IPhrpgkriOXHTGRwD489RaFZMwwqBl+py1zslK2jTm/HnWusZVemW+/tqoxWB5/P1ltCNhDIA0sHZFf4PCD/PvL+qrYFSj7MNcR/7nAqh7ju6eQDTDkjp5aksQlj+f1b2SovP3d5WBUg8wOT57XpyMpjz9/oBMYxs0ACQU17/vunhcNwwwq7TXr15Wdw7Q6DqyCRL7jc/5m59ObFhSxzZo2WzSVBx5C7DqVulTNaeFRc/cZEYqL78/l6AdCQ0AyRqb7Bzd6t2r4G2MiBfXWJh2ABecvlB7fkuK11bd3Fi4BADj7ziKtFUK06Q4tn3W+ey0i4BXnrBTlqZ7ff9zzv+32ggEVTDf9Gqye+t+f1nWRmS4I6EBIKHY+vd1riOXKCmefmOh45mXDxnPDzMuiXGV4WM3AdUjs8erR1q7ag5SEEcPAjvvrv85SDbTvdRkNjuBZq7S89ZWIsMdCYPAJJS4qaFeogSJdcYi6LpMitfOWue4IvwkCSZGDQRGVRBBspkCrGHXNYtmVzDrgs5ZkWG1Ng0ACSVuaqiXKEbEJsDrvS6NDKRYJN26exX+rUuBh6+LlppoM+QlrmxpXReXPFUwNztDJ8N3pQuIhBLm37fBJkjsEpT9o7sujQykWIRt3YOyavxBV68rySUsEFjnytDIYZJZR/Xt4Osq86cby1m4TMKyiWyzjfJQwdyq4HhG78puoKRl2AZqdTEAtx6gX3Od7vxKuWRVv5CIoC6WQHCHS1OXzgbEcVPEkcVPUIfNIHlKPYBSThZO2L3COnu2W+fPduimaoGpGygNAEmVtLJxot4nsxbWptVskEKdtzjCit2jaKKurE1ZQLr7AHoD4k4O0+1QdEowSGGuutnJKlKa+E5eFeqmPjSWIgLWhjkn0AAQa+Iq07CVeKHmDBgVh4uuxtlHM1bOYbsWnYGJogSD3rtcCdih5FShdvgOgDGAAhKkiJOkVAZl4+x87Qjuee71GdXQ9FTNrDHFB2YIUf7uitnr50+jWCjoPqZsGOO7KEdBeuU0nSulYPdUXmf75q1mIGWYBVQwwrp2JkmpDOrm6VX+Ue+bKXEzQOJk6MwgjcrYJuPIRlaTUTp60HxN0Lv4s5VMKY06t4/38zwpVO/3+NQtwNlX5SMbqQnQABSMMAWfJKXSlHVTEjGud0fHqto2z7mYDJaka2Rdap8BKemP63LyTStkb8aRjaymZwLma8LexVsnYEppNF0rpXwpVN33uOdex0DloWYgZWgACkaYgg/K1w9TyqZunpMhcSb/TiQ3k8Hido10V5DbrnV+HvyMflW88monw8bPu79sVMJhxUK2sgatxP3XzKyE5znB2yCXlncnoiuyMsn/ib/Kl0LNY6fQJkIDUDDCCrJMRV8XnL4wVCl7W0kDVmHOOqq1SdywdQ+uf2B3ZmMp64hT6GVaQercCJfeBvSc2HiPqVqjwgkrFrKVNWhH4r2m7j0QbjjCfPh5KuwKIm+dQpsMg8AFI6wgy1T0tWn7fqtePm4/IFMXzzCCdgupV/aGpVVG6dEycy/N+bWqk46pyxoxFWDpFE5QsZCtrKZGdv5rdCthEyYf/t6t9f2SKguAj92aXOmb0ljT6O2Tt06hTYY7gIJhM/Bl7Yp+PLvhQry6+RI8u+FCAMBYtaa9X5yYQX9fBfN7y5FlT7Wy18Znbtujxb9a1mFaQYb59m0Jk9XrluquzOb3Q/TXWK14A1bye7cCQ3/c2Czv4euStVLQ/d6G/jh6Gw0T2oC3BAfJm0ULhsRwB9CmJMmpD+raqSPI9aJTykPDo+gS0a7m5/eW8eyGC606fnqx7T1k/b3YpFXado20WS0bFPq/9p2P9x+936+GnSKuKATJqm07IU5s4pTz9NeEtZcIy4N/6pb6ymGXyfFkfe5NcwL8xO2n39Bew+PIbOWMhBa1oKABaEPi5urHNRpBq3m/UnZlM7lyfnV8AkPDozPP/foP9uPtY/rdhYuu/YOOSN+Lra/XpkdL2GrZ4CIZGh7FuT/7B0iD9ofjMvIT5rLyGwE3jqA1UArY+V3HAOgUeZCryCZtM7BVdQJ/epRr3VhGVNeQ+zvXFYG1aFBLq4bEZOYCEpHFIvKMiLwkIvtF5EtZydJuxMnVD8uscTN8lmx4FL+58e+xxJPpE+R62bLjQF0gOKyVc21Kzci5dkU/ghKEyl2CO648B89uuNBK+d+wdU/o9+K+58jUSfobxfH1Bl0TEOzcsuMAFuEt/XV+xWPjsjKdY1zJK3N2iz/1000ftQ3eBn4nCfzpUa6tzE82/L3ZAeEgF0+LgtFZxgAmANyglPoAgPMAXCciZ2QoT9sQx+8eZDS8xgGYDcS6RuKC0xei3KVbpjYaEptArfccU2wBQIN72kTYrsN9nvc9n5o6R298TK6XoH+sJv/75X8dmDceKahtk55oOicw9z9Aocykcx4FvnbE+XPVzc5zwvzSq252Zvj6KfUkK/rSfdfurGAv7jlJUjpt4jNx/fRhBj2t2FAImRkApdSbSqkfT///LwG8BKADewIkx59/32cIoAat1IOqdL/+g8YMH5dqbRKP7HkzUBm76ZthuwUbOb3UJpVV6mfYrsN9nve8S7uei+Z6CfrHGjPFMVJQ22ZFGDTVy/QLjKJQohTGnbUOWPuXnmAznP//+HeC20iHKVPdd732L537+r//KBlWOmwC63F3GGEGvUVDYnIRAxCRJQBWAHhe89m1AK4FgFNOOaWlcuUBnV+73CUolwS1ydklbFiQNKjHfpgPPnCVPs2kUti4bR+uWNmPh3aNGhWyX875veXA50fdUQQ9zz1vTdcPsUB+pb9ApxxsA8YRfbPrVy/D5FAXujHV+KF/1W6TnhjUf6jcC9TeqT8WdTVu65f2+91tUj+jBD1N37UuOJ8kpTMsCSCJnz7MoLdobGXmaaAiciKAhwBcr5T6hf9zpdSdSqlBpdTgwoULG2/Q4ehWt7Uphbk93YGpnH50BV5pU61N4pmXD9WlmfZVypjfWzbK+bXLlqNcMm8vbCqQg1pQeJ/nnndj91b96h/QKwdLf6xV+wrPKnftP67G4ZN/W18st/Lq+p9tVoRBmUN+5Q8gMACjw7YfUZxVcTMqcNNYRQeNjkzip7dx8bRgbGWmOwARKcNR/vcopbZlKUteMa1uj1Zr2P01+1RBVwle/8DuSM8vdznpnFOWumJ0rBop02jtin7sfO0Ivvfc69pnuxXIQZk9puI2v7Fxz1skhsAroFcOFqtv/05t5S+exLlDn4V6+DDE1HP/6EH8evktYOl/BH72Q8dVIyVH+V96W/2zbFaEOvdVEG7Fsa1isdmFxF0VNyPo2exVdJKisZx0Gc3MAIiIALgbwEtKqdvCzi8qaY47XLuiH1t2HAit0BVxFod9lTLeGZ/AlMZDYbwWmLm/bXrqMy8f0h4/cU43nnn5kDZ4vWn7/jpDc8XKfjzz8qFAw+P+/G8PL8T7oHlmZYFeOVj8Y/Xu1NZ0/RCby3ehV8adD91VcLemH36tChz5qRNkDSPMzRRHWfqvCUqbtFFacRV5sypwmzlqMYkSb5GLJ4wsXUDnA/h9ABeKyO7p/y7OUJ5cksZA9rD7+Vk0r4I7rjwHvzw+URdncCmJ4NPnnaLNDIrT8tm0yxk7VjN/Vq3VpbQ+tGsUd5zxCl799Zvw7PHLsfYfV2vdDmtX9ON9l/8PvWvgY7fqBdQFHs++qi4bZvAXT86cfmP31lnl71Kr6qdqAXol7A+G2gRI46awepu+bbs2WbBb18nUPR4kf4uCnqmStL9RC1w8YXAiWBuQ9iQt935hg9dNgVwBcPuV52D9g3u0BkJ3/qubLzF+buob5DaVs+kptKbrh7i1525U8O7swaCJWXEKhLzX+lZ+VZyAm8Y/g+1TH8ZPT7gKhqxZPf6xj/5VpWkm79lX1Y98PO0iYPj/OdW2Nrj32HNvSG8gy+lXbvsHf2VuVwlAV/3xUo/TCK/6dvr9fEgDHAlJGojbsC2KYnbPX796mfUUMsAxQFes7Meje98MzVICgB/2fBEDXRrffjNG9xnGBI6qk3H+u98yy1JZAExUg8c6Wg+LBxr6rZZ6gMkJQJdV5DJvcb2C9TZrC3qOzbhGo+wRx1+S1DEZgMyzgEh2xMkMEgAXnL7QuojJ1Er6yw/snqk2BtDQoM5NJ/UrfzeryI8xsGvjF49azGO45yI5jKtP/BF65XijunNdTGEug0h+fN9TJscRqvy9LgfAQvnD7NbxY5TdYpHZwT3380wu6gBIcxgaHq3rtdNXKWPTmuUzK29/62dTAzcvCsBDu0bRF5C/XxLBlFIzK31dKqt/NvAVK+tdWo/ufVPrgpp7Qrc26+dNnIx+XVuFML94nKZbhoClVOZj08T/BcRnHP1tkINWuaGzhOPi6WjpbfiWJkHzgMPmCQDOtZv66P5pIaE7ABH5vIhYLgFIXhgaHsX6B/fUKemxag3r/25PwxAXt/XzlKU7sFqbDEwhn1JqppX02hX9obuFam0S33vu9bodgsm4vDFW1ba0fmPljfGCiHHyz0359pPv6n3pPXPtlZkuGFrq0bRViBJkABo6Wu7dar/bqB5p3B3pdk2mQO7KqyPMR07YzplEwsYF9D4AL4jIVhH56HT6Jsk5W3Yc0AZoa1MK1z+wO1JBlY6j1Rr6KnYtKeYZzouDe2//zIJz1/zXeBkZcdIWTfn245piq7B7+dFllnz8O067A++xwWviGwrXwEXKGvIo5ke+oi/2AvS/g0tvqz9eWaDvE6STkTQVqyDwtNK/CMAfARgEsBXA3Uqpf22uePV0chA47UyfpRseDfW8+oulovTodwO7/vMFwH857xT82doPztzTNlsoDDeUaNse2gpT4DIoeLypD5GGXTYjEA3YTcYyupMEuPzO8Olg2ksNLp0o7+mV3fhdWgafSSimILBVDEAppUTk5wB+DqeL53wAD4rIk0qpG9MVtXjE7e8fRFDvHxdvjr5rfOZVyphT7grMvHHrENwq3nuee33mn7AbIxj8jQUzhWdxlX9fpYy5J3TPvIc/bgDE/35cXvjNL+DMH/9pY/pokOsoTLHWKTSJPtzFNkXVpieO0cANNBYjVeYDE+/q20Z4Mfnzo+50wjKfOnQMY56wiQF8UUR2AfgmgGcBfFAp9d8ArARwRZPlKwRR+/vb9JxZv3pZYI8dl9GxKq5/YPeM/32sWsPx2hQ+fd4p2gyhvkq5btfwyJ43A4u/4s7xFQCb1izH+tXLtN7uNIbEDw2P4g9e+A3cNP4ZjEydjCklGFUn44UPfj3YdRQYV1Cod7soJ88+SZvgbdc6bpc4hBVYeYuRbnoV+OobTtvnTUfNA+RN7aXjKux2LALrEEJdQCJyCxx3z2uazz6glHqpWcL5aQcXUBxXjsldoyug0rlpyl2CE+d0Y+xYre6ZfzK0D/c8/3rknl9AeO6+K0tQbyEBAjOLLLLDAwkrMAsjqADNnYVs5Nal+hTKpO6RoFz6y++MlxkTtKMI+0zX6kBXPObP449aaOc930079RaJMSMoEbFdQEopoxlupfJvB+K6cqL0+zF1B3VdNu4zd752BA/tGo2l/IHZbJsguTdt3x94DwUEppUqOCmjYamngFPpe2P3ViySt/CGOhnfnFiHXe/9vdDrgog0WMev0JZ/Qq8ETf50W/dIUC69TeM2k+I1VUMHpcAG9asxzRK2ua8O93ktmoVLHFgHkCJBrpywVs26Slhdvx8bl0q1Non7nj9opVhNzKuUcc7Xn5iZBTC/t4yvXbYcAKwaytny3ko3fnV8ArWAdqP+5moD8hY2l+/C/jOWAAhZqQdgbXh1SmnPvY2tGNxpWUn82UHxhTAjElV5mlJgv/85x+3kvpN/5xK2uk/SJ79Fs3CJAyuBUyTOqEbA2R1csbIfpekM25IIrlipX33bpmomUf4A8IvjtbpBMG8fq+ErW3dj/YN7UlP+7n0hTmzBFLHQNVfrlXGc+6//O9GzrRvtmZTSK080NvMyjSwcf8eu0njVzYg9vStqTUPgBLGIc4fTmmfbolm4xIEGIEVMyjlMaQ8Nj+KhXaMzSntSKTy0a9QY3G32YJcTuru0/f+nFFJJ5/RTm1SYe0I3Xt18yUyfIS+J2jwEoCso0w7WiaKU/Hn8lQVOf+3qEVgVOZ21zsnx9xsBm6BokJy6wi2bXUmtCmz77Ow1NkYmyTzbFs3CJQ40ACkSt3VzlCwgr9IC9GvFpJV6705EGACQEm+MVTE0PIpj4xONn6mT9RdFnWer6ffjLyjTuuqiKiVvZk3P3MbunGFFTpfe5gR8oxa1meSpzNev2k+7yL5C173Gxj2VJKuHGUEthQYgJdzsn2ptcsaVYzOqEYjuOlq7on/G2PjX4/N7y4kya7Kir7eMjdv2aesPvjmxDuNyQv3BKEohyfBuIJlSiuvSiNIr3jVuRw9Cu3MA9Kv2XX/rxDFcQ2NK7/ReY5MCmqRPftIe+yQSDAKngD/7Z1KpumKpMMKCkbrUUt2uAQB6e7rR29Md6qcviaCnW1Ct2a/2uwQQEUxazIfsEliPkayUS1AKxgrk7VMfxoJSDzbNfSher3iT2+Kxm+wLrtz7RH1+syZduTSkarp1CMpRnqtudgK6OtSkMz/g49/RZ+CYrvFnO7nGMMmMBS/NnOJF6uA8gBSwzSc31QiY+uF/43KnnYKu3UJA8Txuv/Kc0JYOAgR29PQTRaFHuYebXfTlB3YH7lwS5fzbtm5oRk96Uy59Ws+xaWURNmOgssApAnPlNWUyufedyXbyKPrXnwN2fhd13zN7/OcGzgNoIjYuHFfJeztebty2D0PDo4HByKBWyjoW9VUa4gSm88YslT+QXPn77zG3p4Q7rjwHwzdfhLUr+kMD5XFmIM9gu9puRgOyZrs0bFxMOheWF29Bm+t6uvyvzW4vv3sKaFT+ABu6tQF0AaWATT55WI2AqegqSisFU8DZv2Nwz/POCmg174xPYudrR2beWVcL4ZJkBjIA/fBuE81IN2ymS8PGxeQ+e9tn7e8bxe311C0wLkvc7zMt9xBJFe4AUsAm+ydujYDtyte7a/DuNoDGf5pzyl3Y+doR/Op4Y8ZNK7nv+VnF5a+FcLENpGtxg6PbrgW6K46rw12FVxbor8ki3TDKRDL/ubpMHl2A+qx15nc2HbcNRAcZTXfofJIgPGka3AGkgH+ylq5vTpR2D4Dd4HaX+b3luliDKUDs8vaxGr733Ouh92023mI1fy0EgEiB9Ab8vvfqEUcxuv10TL75VqcbBlXvAvWr5tMuqm8/EVSRrFPWH7tVP7S9eqR+UlhUjNXLMhsvYHVvLqEBSImwvjlR2j1E6csPoKHfT9wOnFGx7eMTdL1L3DYaRsKUjo2LY+/W+qHp/tGOaRCUoeQdIH/0oNnP7lYkh1H3zm7KqG9SmPc8L0EuHK2LTYClHwkeOs/q3syhC6hFWFedInwF72esWqurGk4UMLXkZ5svwf9ad3aie0wqNdPOOq6LzIhNcDTIxbF3K/DwdfXKq3rEWUGn6bowyVk9oolZmPzsB+0H2rvvPG9x4/1MQdswF44u0D14DfDaPwUPnWd1b+ZwB9BCwnYJLnGUnrfraFBANU3WrujHpu3763oGRcXNhppXKWvvE9uYJc2/f+qWxgpewHGfpOm6SG0IvApfxXuJUqBm48LxB7pvP7PR1eSF1b25gDuAHGJSevN7y8Y+QN7WEe5uwzSzNymVctfMQBoRZx5BEqq1SYggVhsNI0lbCgS5J9JyXezdqp8jXOoBJOY/TdvUyyjtLeJUM4d9R6wPyAU0ADnElFX0tcuWzxSH6fDuHNau6MemNcubYgQmptRMPcPbx2qoTSlUys5fJW9HU++fYYwdq1m7yKxImn8ftFNIw3XhulX8LpKeuU5QR+kqtC0NrY2BimIg4zRoC/xsMZV/TqALKIeEZRWZsoO8O4c/GdpXN6s3DqaKY11H0OO1Kdxx5TkNCsp9aewAABNESURBVHvphketnuUWsBkVfpw88iT596tuBoY+B0xp3GhRZ/y6eN9BuvSTw2rH9celBKy8Gth5d/hzbAxUlDx/XZA3bDe16mZ9xlGph66fHMFWEG2Ct41EX29ZO0TFO7QlrLVCM9CNUjS1yTBdr037bHY7BROmsY+24x292PTZCUScYLVJJpdmfS9xDHArsqiIFaZWEDQAbcDQ8CjW/92eBoWvW6FXyiXMKXdlUuFrO8M4CLcHUp0RsOl34yWtqlNjD6FpZRyEX4bxd4IV98ytQ2YKaw2Jr/kbFSzxEXsmMMmeTdv3a0cm6lRTtTbZ9OwfE329jfEGvzsrbLmhzf2PEoRMc6Zs3EwinQw2BA1cd90mSTqTEuKDQeA2IEmaZSv51fEJ7RQzd+jK7VeeYxXGbEiDjRKEjDoWMYi4mUQ6GUxICXVB6ktvCw9eR5kV0E5EaYlBUoEGoAPpqzSmiwqA3nI6v+7+voo2u6g2pXDD1j1YuuHRmQIvL1t2HLCKSzSkwVoo4qHhUZy/+WlMjcUcqK4jbiaR7bPKFeATf+W0pwCcnkW3n+n8fycq+CDYLygTMjUAIvJRETkgIj8RkQ1ZypJn5mtcKyYq5RI2rVnekFJ5+5Xn4FiE4S9B97/g9IXGXcmkUg3trl1sCty0uf8hitjb/C6V8ZH+Z0dVxsbRjAsa3+H15xzF3yrFl8Uq2+aZae7ciDWZxQBEpATgOwB+D8AIgBdEZLtS6sWsZGoWpkEwtnztsuVY/+Ae40B2163iv7f/GTds3RO5d8/83jJ6e7pnZL/g9IV4aFejm0eH359vaohXEsGUUsHfTUBKp7d1xjcn1mFz+S70iqeKt9VVp6ddpE/XXP4Jx8UDNGbIePFX2aYV1E4zPpL2M+OOziSJyDII/NsAfqKU+ikAiMj9AD4OoG0MgI1i92fBuCtjoFFBm3DPu/6B3cZzbKZlxWncphTq3uv8zU9HCjJ7V/2mhniJCr58z9g+9WGgBtzYvRWL5DC6+nwK06tMK/OdY9W30w2mvvJE8HGblNCjI3ojEVVpP/IVZ/avLrMIaH5XTttxnJX5hpRb9gtqJlkagH4AXoftCID/4D9JRK4FcC0AnHLKKa2RzAJbxW7qcrlp+34r4+E9py9hv5x+wwo8iLFqre69ovYp8spm0zY7Dv6dxfapD2P7+IeduoQve+oSdC2iXdJaDe/das76cVezNkHiynyzkbBV2o98xa5wrJmr7KBmd+73f/Qg0FV2isQmM9y5FZAsYwC6hJCGJapS6k6l1KBSanDhwoUtEMuOoPbFXkwKc6xa046HdNGNkHxnfKKh706Ufjm6FhM2eN8rSnM2nWxuRtCrmy/BsxsuTKz8AbuBPADCFW9Sn7NrYEy4q9lQhSuz8pg4ejDch7/rb0Oe45MrjDjxA9t7T9WAnhObNzqTaMlyBzACYLHn5wEAb2QkS2Rs2xeb/N5+qrVJ3LB1DwAYZwHXJlWDTz5sBe3fRVyxsj/WMBj3vYI6jZZLgrk93TharaW2urfBemdhs9JNshoOMjDe1WxYB9DBa6Z7/4fg37X4YwUmt49JriDixg+ijOOsvj07nJ60hMwqgUWkG8C/AFgFYBTACwCuUkrtN12Tp0pgU4sDfzuEqD15yl2CE+d0Gyt5vdW2YTEIXRVupVzCCd1dkWsLvO/lnVbmDoUxtnEwkDQwHgtTRbGXOG0eXIyVwwAGPzM7tcvk7wacTKGbXrWT1cWtAI7UakJmp4zZTBOLWo3txbYqOsl3TwLJXSWwUmpCRD4PYAeAEoDvBin/vGEz4csdc+hVCQKnnbIpJbM2pQLbOLguGJsYhMlNNafchUq5FKk9g/e9bOcamEgjMB6LMCWZ1Ofc06tv71xZUF/dWz0CdJUaG82VK06vHBtZvRwdiVZ8NvgZJxspyqo+SZaOP4Mr63GcHFA/Q6Z1AEqpv1dK/Xul1G8qpf48S1miYjPhS6eAFYATyqVYvnivIraJQRjjD77Wy32VMsql+tiC+1NYW2a3AMtU/KXDNn6SOv56gsqC+kHxSXzOj3xFr/wBYPLdRuU8NamvC/AOWDn7qulKYTh/lufq7x/mUvLew1X+QLTc+zgtoU0kbdWdBBac1cFeQAkIWwkHKeDbrzwnUl6+38ViE4MIGkTvlz2OSybuSj718Y9RSNIiOoiggKvJMAT5vPdudXYNrh9fTQI1w31MhgEIdqtEWdXHaQkdRLN+D2FwQH0dNAApoVOgYQoYgFWnTF2b5aB7u0QZRB/HrRN3kLuN7G2HTcDVT9DqOYpL562XzZ8FKegoze46pQkdC87qYC+gFNClbG7ctg8XnL4wMD3R70bSuWJMCluX+ijTz3ZdMUFuqjiuGz9xV/LWaZvthAS49CoLojeVS0shhWXolHrqjwUNbOmEJnRpurI6AO4AUsC0En7m5UP4xuUfDHStxHXFuMe8Q9ldZ5LfFdOM6mQg/kq+WQVhmWKc1tU1G9iNsnpObVh8CH4XZBvNB4lF2q6sNocGIAWCVsJRXStRz393Qp9NFOSKieu68RPFxeQnaSZR7nADq7v+Znaeb3kucNkd9YFdWyKndWqoLAj+/KlbGkc2TtU62x/eKa6slKABSIGsfNo6Re4lqosmahC2I1fySbj0tllDkBS/oqrMB8Z/Vd8qwaXU42QVeeMQpZ7ZnYeJovrDswpA5xAagBRIshKOitdFFLZZNxmgNA1W267kbXLBs84X1+XPP3WL4xpyR0e6RWBA9PdhA7bCQwOQAq1aCUeZrxtkgFppsHKJTQFUFq2TwwhbuQZ9pnufUo/ThM3rBiqwP7yIcCh8G2FqP+Gnr1LGpjXLI/UIyoPrpmUy2bQ1SNL6wE+UnUSzdh2m96ksAHrm0h/e4eSuFQSJTpCPXtA4ECaIZrpuWllUFgsb33da/vEoO4lm7jqMbZnZgK3IsA6gjZinmcMLOCv+NNsrJ8FUExFWZ9DS9hA2ueBp5YtHabfQzLGIzH8nGmgA2gjRTVAIOJ4FcRV5S9tDWAyZtzrHBuNOQtPPP2zXkWSeb1rv4yWL+cIkVWgA2ogxQ5dQ0/EsiKvIgzKWUkfXjOzsq5yVtqvMgHQalgWtsP2NyIzD5OcDty4Ftn02fhOztBuwmZqqPfIVGoU2gjGANqIdeujElbHlmUnejBqT7/2yb8XvT+9N2YTAOCfA24hMV/xV6gHe/WVjwZb/WhvSzH83uat2fhcz75qHzCkSCHcAbUQ79NCJK6NNe+2mkbbvvW51DDgKMcBP57p4dKv0nhP1yt9/bVKiunOMz/UZurRiGKQpcAfQRrRD5W0SGTMrKku7IlbbyVPNFm/58bp+/Kv0TX3Bz0ojiBsn+yhKryL/95h1gR2ZoeMNQB7z3U3YyGpq7pand2y76uAobZFtMBkONekEXqM0IgtStGkVbcXpka/tVWRwdXm/xzwW2BWYjnYBxU1JzIK4sqb9jmm0ic4FUVwaaWfIGFMuF0cPxK662anW9VNZMHtt0mycODsgnbtq8Jrw77GZqa4mmK1kpKN3AGl1vbQh6So8rqxpvmNms3rTJuoqM+0OkUEth+MEYv15vm6jN1f5J11Rp7EDGn8H2P99Rw5/nyKvHEncbXFcR9xxBNLRO4BW5ZansQqPK2ua75jZrN60ibPKTHPYSZopl0/d0tgBdHJ89l3SWFHH2QH500CrR2Yby7muLp2CjluQFneWbxY7jjaiow1Aq3LL01CccWVN8x0zndWbJnloc5yWQQl7F5t3DXOBxDFYYSMrTUo2rrstriLPw9+FHNPRBqBVaZNpKM64sqb5ji0txmomndT2IOxdwj63XTlHNVg2ClR3TtzdUVxF3kl/F5pARxuAVuWWp6E4o8rqBmu//MBuzCl3oa9STvyO7VBnYEUz2h7EIY3g42kXBR8Pe9dmuUBsFKjXCHm/ByD67iiuIs/L34Wc0tFBYKA1KYlpVbEGyeoNMs+rlPHO+ARqk07K3dvHaqiUS7j9ynMSvWs71BlYkYexf2kFH195Ivh42Ls2ywUSNrLSVbJpfQ9xZ/nm4e9CjuE8gJRoZi6+7SCY/r4Knt1wYSrPJDHwZqlIl6HoK+I8gU190LeREGcFHUaacw38eN+3Mt85Vn27XslmNVeB1MF5AE2mmTuNsNm/Lm0XrO0k/CtdnfIHoq+8k6Zoxl0522CT0prmDoSzfFOno2MAnYKtYm+7YG0nEZYV42JS3KZ4QVIfdtpdQKPCIGyu4Q6gDTB12PTSlsHaTsJmRWtS3DZ+8iSujyxXzs3cgZDE0AC0Abogc7lLcOKcbowdq7VvsLaTMLlqpASoqWDFHdaLp51dHwzC5hoGgdsEU5A5b43gCot/FQ84K10bd0vSQC8hITAI3OaYuoB2RO+eTiDJSjftbqSEWEID0Ma0stkdsSCuq4Z+cpIRmWQBicgWEXlZRPaKyPdFJGTqBdHRMb17ik7WmTqksGS1A3gSwEal1ISI3ApgI4CbMpKlbWmHGcHEklYHellURZDRDkAp9YRSamL6x+cA0NkZg47p3UNaS9zWyqTjyEMh2DUAHjN9KCLXishOEdl56NChFoqVfzIdpE7aF/bIJ9M0LQ1URP4BwPs0H31VKfXw9DlfBTAI4HJlIUiR00AJSQ2mnRaOlqeBKqV+N0SgPwRwKYBVNsqfEJISTDsl02SVBfRROEHfNUqpY1nIQEhhYY98Mk1WMYBvA3gPgCdFZLeI/FVGchDSWtIYEpMUpp2SaTJJA1VKnZrFcwnJlLSGo6RBO/cXIqmRhywgQooBs29IzqABIKRVNGs8IyExoQEgpFVwOArJGTQAhLQKZt+QnEEDQEirYPYNyRlsB01IK2H2DckR3AEQQkhBoQEghJCCQgNACCEFhQaAEEIKCg0AIYQUFBoAQggpKDQAhBBSUGgACCGkoNAAEEJIQaEBIISQgkIDQLIlDxOyCCko7AVEsiNPE7IIKSDcAZDs4IQsQjKFBoBkBydkEZIpNAAkOzghi5BMoQEg2cEJWYRkCg0AyQ5OyCIkU5gFRLKFE7IIyQzuAAghpKDQABBCSEGhASCEkIJCA0AIIQWFBoAQQgoKDQAhhBQUGgBCCCkoopTKWgZrROQQgNcS3OJkAG+lJE6e4Hu1D534TgDfK+/8hlJqof9gWxmApIjITqXUYNZypA3fq33oxHcC+F7tCl1AhBBSUGgACCGkoBTNANyZtQBNgu/VPnTiOwF8r7akUDEAQgghsxRtB0AIIWQaGgBCCCkohTMAIrJFRF4Wkb0i8n0R6ctapiSIyEdF5ICI/ERENmQtT1JEZLGIPCMiL4nIfhH5UtYypYmIlERkWEQeyVqWtBCRPhF5cPrf1Usi8qGsZUqKiHx5+u/fP4vIfSIyJ2uZmkHhDACAJwGcqZQ6C8C/ANiYsTyxEZESgO8A+BiAMwB8SkTOyFaqxEwAuEEp9QEA5wG4rgPeycuXALyUtRAp8xcAHldKnQ7gbLT5+4lIP4AvAhhUSp0JoATgk9lK1RwKZwCUUk8opSamf3wOQDtPIP9tAD9RSv1UKTUO4H4AH89YpkQopd5USv14+v9/CUeZ9GcrVTqIyACASwDclbUsaSEi7wXwEQB3A4BSalwpNZatVKnQDaAiIt0AegG8kbE8TaFwBsDHNQAey1qIBPQDOOj5eQQdoiwBQESWAFgB4PlsJUmNOwDcCGAqa0FS5P0ADgH4m2nX1l0iMjdroZKglBoF8D8BvA7gTQBHlVJPZCtVc+hIAyAi/zDtu/P/93HPOV+F4264JztJEyOaYx2R1ysiJwJ4CMD1SqlfZC1PUkTkUgD/ppTalbUsKdMN4LcA/B+l1AoA7wBo61iUiMyHs5NeCmARgLki8ulspWoOHTkUXin1u0Gfi8gfArgUwCrV3oUQIwAWe34eQAdsVUWkDEf536OU2pa1PClxPoA1InIxgDkA3isi31NKtbtiGQEwopRyd2kPos0NAIDfBfCqUuoQAIjINgC/A+B7mUrVBDpyBxCEiHwUwE0A1iiljmUtT0JeAHCaiCwVkR44gartGcuUCBEROP7kl5RSt2UtT1oopTYqpQaUUkvg/J6e7gDlD6XUzwEcFJFl04dWAXgxQ5HS4HUA54lI7/Tfx1Vo88C2iY7cAYTwbQAnAHjS+d3iOaXU57IVKR5KqQkR+TyAHXAyFb6rlNqfsVhJOR/A7wPYJyK7p4/9d6XU32coEwnmCwDumV6E/BTAH2UsTyKUUs+LyIMAfgzHTTyMDm0JwVYQhBBSUArnAiKEEOJAA0AIIQWFBoAQQgoKDQAhhBQUGgBCCCkoNACEEFJQaAAIIaSg0AAQkgAROXd6tsQcEZk73UP+zKzlIsQGFoIRkhAR+TM4/X0qcPrifCNjkQixggaAkIRMt0B4AcBxAL+jlJrMWCRCrKALiJDkLABwIoD3wNkJENIWcAdASEJEZDucaWxLAfw7pdTnMxaJECuK2A2UkNQQkT8AMKGUund6RvM/iciFSqmns5aNkDC4AyCEkILCGAAhhBQUGgBCCCkoNACEEFJQaAAIIaSg0AAQQkhBoQEghJCCQgNACCEF5f8DMrVDWwbpZJAAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.scatter(data[:num_samples,0], data[:num_samples,1], label=\"0\")\n",
"plt.scatter(data[num_samples:,0], data[num_samples:,1], label=\"1\")\n",
"plt.xlabel('x')\n",
"plt.ylabel('y')\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"from pyro.nn import PyroSample\n",
"from pyro.nn import PyroModule\n",
"\n",
"\n",
"class BayesianNet(PyroModule):\n",
" def __init__(self, in_features, out_features):\n",
" super().__init__()\n",
" self.linear = PyroModule[nn.Linear](in_features, out_features)\n",
" self.linear.weight = PyroSample(dist.Normal(0., 1.).expand([out_features, in_features]).to_event(2))\n",
" self.linear.bias = PyroSample(dist.Normal(0., 1.).expand([out_features]).to_event(1))\n",
"\n",
" def forward(self, x, y=None):\n",
" logits = self.linear(x)\n",
" with pyro.plate(\"data\", x.shape[0]):\n",
" obs = pyro.sample(\"obs\", dist.Categorical(logits=logits), obs=y)\n",
" return logits"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from pyro.infer.autoguide import AutoDiagonalNormal, AutoMultivariateNormal\n",
"\n",
"\n",
"model = BayesianNet(2, 2)\n",
"guide = AutoDiagonalNormal(model) # restrict variational distribution to diagonal Gaussian\n",
"#guide = AutoMultivariateNormal(model)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"from pyro.infer import SVI, Trace_ELBO\n",
"\n",
"\n",
"adam = pyro.optim.Adam({\"lr\": 0.03})\n",
"svi = SVI(model, guide, adam, loss=Trace_ELBO())"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[iteration 0001] loss: 1.1963\n",
"[iteration 0101] loss: 0.1787\n",
"[iteration 0201] loss: 0.1413\n",
"[iteration 0301] loss: 0.1327\n",
"[iteration 0401] loss: 0.1425\n",
"[iteration 0501] loss: 0.1311\n",
"[iteration 0601] loss: 0.1340\n",
"[iteration 0701] loss: 0.1330\n",
"[iteration 0801] loss: 0.1259\n",
"[iteration 0901] loss: 0.1299\n"
]
}
],
"source": [
"pyro.clear_param_store()\n",
"\n",
"for j in range(1000):\n",
" # calculate the loss and take a gradient step\n",
" loss = svi.step(data, y)\n",
" if j % 100 == 0:\n",
" print(\"[iteration %04d] loss: %.4f\" % (j + 1, loss / len(data)))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"from pyro.infer import Predictive\n",
"\n",
"\n",
"guide.requires_grad_(False)\n",
"predictive = Predictive(model, guide=guide, num_samples=1000,\n",
" return_sites=(\"linear.weight\", \"obs\", \"_RETURN\"))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# slice through x\n",
"\n",
"x_ = [[x,1] for x in np.linspace(0, 7)]\n",
"x_svi = torch.FloatTensor(x_)\n",
"\n",
"samples_svi = predictive(x_svi)\n",
"\n",
"y_new_svi = []\n",
"\n",
"for i in range(50):\n",
" y_ = torch.bincount(samples_svi['obs'][:,i], minlength=2).float()/samples_svi['obs'].shape[0]\n",
" y_new_svi.append(y_.numpy())"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxV9Z3/8dfn3puVJBAgEEjCHpaAG0ZwQUVxQUWoVlsc1NqxtZ0ZrV1+nS4z2k6nWjtjO+1U22pdqrjg1iqjuGAL7iAg+x72ACEBQsi+3c/vj3OQEJIQyD05N7mf5+NxH3c5557zviz3c8/3+z3nK6qKMcaY2BXwO4Axxhh/WSEwxpgYZ4XAGGNinBUCY4yJcVYIjDEmxlkhMMaYGGeFwJh2EpFZIvKO3znaS0T+LCI/dx9fKCIbT3E7fxSReyKbzkQTKwTmlIjIJBH5WETKROSgiHwkIueIyHkiUikiqS28Z7mI3CkiQ0RERSTUiXl/KiLPdGQbqvqsql4RqUydSVU/UNVRJ1pPRG4TkQ+bvfebqvqf3qUzfrNCYE6aiKQBrwO/A3oDWcB/ALWq+glQCHyx2XvGAXnA852bNjI6UrTE0aH/a51ZNE3ssUJgTsVIAFV9XlUbVbVaVd9R1VXu8qeAW5u951bgDVU9cKKNi8hCEfmFiHzqHnG8JiK9myyfLiJrReSQu+6YJst+ICK7RaRcRDaKyBQRmQr8GPiyiFSIyEp33Z4i8riI7HXf83MRCbrLbnOPcv5HRA4CP23+a1lEzheRJW7GJSJyfrPPcJ+IfARUAcNa+JzbReRHIrJOREpF5EkRSXSXTRaRQvfzFAFPuq9PE5EV7mf/WEROb7K9s0TkM/ezvwAkNlk2WUQKmzzPEZG/iEiJiBwQkYfcP8c/Aue5f06H3HWbNjGtF5FpTbYTEpH9IjLefX6um+uQiKwUkckn+vs2UUBV7Wa3k7oBacABnC/8q4D0ZstzgHpgkPs8gHOU8AX3+RBAgVAr218I7AbGAT2AV4Bn3GUjgUrgciAO+FegAIgHRgG7gIFN9jPcffzTI9tosp9XgUfcffQDPgW+4S67DWgA7gJCQJL72ofu8t5AKXCLu/wm93mfJp9hJzDWXR7XwufcDqxx/7x6Ax8BP3eXTXb3/0sgwd3/eKAYmAgEga+420hwP/8O4Dvun8sN7t9B0+0Vuo+DwErgf9zPnghMavK5P2yW889NtnMv8GyTZdcAG9zHWTj/Lq52/84vd59n+P1v1m5t3+yIwJw0VT0MTML5Mv8TUCIic0Wkv7t8F/AecLP7lik4XzZvnMRuZqvqGlWtBO4BvuT+Wv8yzpHFfFWtBx7E+ZI8H2jE+VLME5E4Vd2uqlta2rib9Srg26paqarFOF+MM5ustkdVf6eqDapa3WwT1wCbVXW2u/x5YANwbZN1/qyqa93l9a18zodUdZeqHgTuwykoR4SBn6hqrbv/rwOPqOpidY7EngJqgXPdWxzwG1WtV9WXgSWt7HMCMBD4vvvZa1T1w1bWbe45YLqIJLvP/8F9DZy/73mqOk9Vw6o6H1iKUxhMFLNCYE6Jqq5X1dtUNRvnl/tA4DdNVmnaPHQL8FwbX4Yt2dXk8Q6cL7m+7n52NMkRdtfNUtUC4Ns4v/6LRWSOiAxsZfuD3W3udZsxDuEcHfRrJUNzx+RokjOrne9vaZ0d7naPKFHVmmaZv3ckr5s5x33PQGC3qja9imTzfEfkADtUtaEd+Y7h/hmvB651i8F0jhaCwcCNzfJNAgac7H5M57JCYDpMVTfgNB+Ma/LyX4AsEbkEuB54+iQ3m9Pk8SCcZo79wB6cLxzA6Yh1193tZnlOVSe56yhO0wru46Z24fya7quqvdxbmqqObfrR2sh3TI4mOXe38/1HNP+ce9p4/y7gviZ5e6lqsns0shfnz1uaba8lu4BBrXRAtyfz8zhHLjOAdW5xOLLd2c3y9VDVB9qxTeMjKwTmpInIaBH5nohku89zcL4YFh1Zx23SeRmnk3OHqi49yd3cLCJ57q/OnwEvq2oj8CJwjdsJHAd8D+cL/WMRGSUil4pIAlADVOM0FwHsA4YcGb2jqnuBd4BfiUiaiAREZLiIXNzOfPOAkSLyD26H6ZdxRkW9fpKf819EJNvtDP8x8EIb6/4J+KaITBRHDxG5Rpyhup/g9Cl8y81zPU4TUEs+xSkcD7jbSBSRC9xl+4BsEYlvI8cc4Argnzh6NADwDM6RwpUiEnS3O/nIvxMTvawQmFNRjtNhuVhEKnEKwBqcL+WmnsL51XyyRwMAs3GOMopw+he+BaCqG3Haon+Hc4RwLXCtqtbh9A884L5ehNPM82N3ey+59wdE5DP38a04nazrcDp6X6adzRjqjH6ahvOZD+B0Wk9T1f0n+TmfwylIW93bz9vY51KcfoKH3LwFOJ27uJ//evd5KU5fyl9a2U4jzp/bCJwO7UJ3fYC/A2uBIhFp8bO4RfQTnH6ZF5q8vgvnKOHHQAnOEcL3se+ZqCfHNika4z8RWYgzwucxv7N4SUS2A19T1Xf9zmJim1VqY4yJcVYIjDEmxlnTkDHGxDg7IjDGmBjX5S5k1bdvXx0yZIjfMYwxpktZtmzZflXNaGlZlysEQ4YMYenSkx2SbowxsU1EWjvT3JqGjDEm1lkhMMaYGGeFwBhjYpwVAmOMiXFWCIwxJsZ5VghE5AkRKRaRNa0sFxH5XxEpEJFVR6a6M8YY07m8PCL4MzC1jeVXAbnu7Q7gDx5mMcYY0wrPziNQ1fdFZEgbq8wAnnZnVFokIr1EZIB7iduI27a/ks37yltcFgwImT0TyU5PpmdSnBe7N6b7aqyHyhIoL3Lu6yqgoRYaao69DzeeeFumbaOmQtbZEd+snyeUZXHsNH2F7mvHFQIRuQPnqIFBg1qbdKlt76wt4hdvbjjhemmJIbLTk8npnURWr2SCAahtCFNbH6a2odF53BCmR0KI7PQkctx1c9KTGdgrifiQdbuYKFNRAkUroWgN1FdDKAFCic59XJJz31AHFUVQUex8oVcUO89rDrvrH3mPewsEoPIAVOyDqgO0b2IzOfEqpm2pmd2uELT0r6LFf02q+ijwKEB+fv4pXSXv+vHZXDCib4vLGsLK3kPVFJZWs6u0il0Hq9haUskHm515ORJCARJCQRLiAiSEAsSHAhwubuDN1XtpCB+NIwLpyfEkhgIkxAXd9znvTYoP0jclgX5pCWQcc59IdnoScUErIOYk1NdA+V7n13Z99bG/vOsqoGQD7F0FRauc9dorPgVS+ju3zNMgsadTJI75dV/j7D99MORMcNZNdd/Tox8kpB5fbIIJTvEwUcnPQlDIsfO1ZnPsfK0RlZGaQEZqQqvLz8zpddLbbGgMs6+8ll0Hq5wicrCKA5W17tFDmJr6I0cQjRysrGPTvnJKymuPKR4AoYAwuE8yI/qlMKJfCsMzUsjtl0rewDSCAfsVFZNUnV/aJRvhQAEc2gmHdjj3pTucX+ttkSD0HQlDL4LM02HA6c4Xe0JPaGyh2SYQcr7IE1I65/OZqOJnIZgL3Ckic3CmPSzzqn/AK6FggKxeSWT1Smr3e8Jhpay6nuLyWkrKayk6XMPWkgoKiivYXFzBu+uLaXQLxbisNH5145mMykz16iOYaLFnBWz/APZvgpJNsH8jVJceXS5B6JkFvQZD7mXOfVoWxPc4+su76S/w9CHOfUsCSa0vMzHJs0IgIs8Dk4G+IlII/ASIA1DVP+JM/n01zryrVcBXvcoSTQIBIb1HPOk94lv8gq9rCLPzYCXLdpTyX29tZNrvPuDuKbl88+LhhKz5qHsJN8LGN+GTh2Hnx85ryX2g7yjIm+HcZ4yEPiMgLRuCXe4akaaL6HIT0+Tn52usXH30QEUt9762ljdW7+X07J48eOMZjOxvRwddXl0lLH8WFv0eSrdBz0Fw7j/BaTdCSotXCTamw0Rkmarmt7jMCkH0e2PVXu55bQ0VNQ3cfVku37homB0ddBW15U6b/qGdzm3/JljzCtQcguxz4Lw7YfQ0+7VvPGeFoBtoenQQDMhxI5MS44KM7J/K/defRkqCfan46uA2eO1OKF57bDs/QFwPp43/vDudETfGdJK2CoF9Y3QRfVISeHjWeK5bt4/lu0o/H5l05NyGqrpG3li9lz2Hqnnyq+eQmmgnxvmibDc8Pd0Zfz/2OmeIZS/3lj7Y6QMQGwlmoosVgi7msrz+XJbXv8Vlb67ey13PL+fWJz7lqX+cQJoVg85VUewUgapS+MpcyLLLZ5muwRqau5GrThvAw7PGs2Z3Gbc8/ill1fV+R4odVQdh9nXOEcGsl6wImC7FCkE3c+XYTP4w62zW7SnjlscXU1ZlxcBzNYfh2RucjuCbnoPB5/mdyJiTYoWgG7osrz+P3HI2G/aWM+vxRRyqqvM7UvdVVwXPz3ROCLvxKRh+qd+JjDlpVgi6qUtH9+fRW89m074KvvzIIl5YspOishq/Y3UvDbXwws2w42O4/lEYfbXfiYw5JTZ8tJv7YHMJP3xlNbsPVQMwOjOVS0b3Y/LIDMYPTreL3Z2qQzvhpdtg9zKY/jsYf6vfiYxpk51HEONUlU37Kli4sZgFG4tZur2UhrCSlhjiT7fmM3FYH78jdi2b58Nfvu5cImLGw5A33e9ExpyQFQJzjPKaej4qOMB989aRGAry5t0X2pnK7RFuhAX3wwcPQv/T4EtPQZ/hfqcypl3aKgT2vz8GpSbGMXVcJvdck8fm4gqe+3Sn35GiX0UxzP6CUwTOugW+Nt+KgOk2rBDEsMvz+nP+8D78ev4mG2balt3L4I8Xwq4lMOP3MOMhu4yz6VasEMQwEeGeaXkcrq7nN3/b5Hec6FRzGF64FULx8LV34axZficyJuKsEMS4MQPSmDlhELM/2UFBcYXfcaLPO/8O5Xvghichc5zfaYzxhBUCw3cvH0lSXJD75633O0p0KfgbfPYUnH8XZLfYx2ZMt2CFwNA3JYG7pozg7xuKeW9Tid9xokNNGcy9y5klbPKP/U5jjKesEBgAbjt/KEP6JPPz19fR0Bj2O47/3v43KN8LX/gDxCX6ncYYT1khMADEhwL8+OoxNpwUYPO7sHw2nP8tyD7b7zTGeM4KgfmcDScFqg85TUIZo2Hyj/xOY0ynsEJgPtd0OOmv5m/0O44/3v43qNgHX/i9NQmZmGGFwBxjzIA0bj1vCLMX7WDZjtITv6E72fQOrHgGLrgbsqxJyMQOKwTmOP/vylEM7JnED15ZRW1Do99xOkdFMcy9EzLGwOQf+p3GmE5lhcAcJyUhxM+vG0dBcQV/WLjF7zjeCzc6VxOtKYMvPgahBL8TGdOprBCYFl0yqh8zzhzIwwsK2Lyv3O843vrg17B1IVz1Szt72MQkKwSmVfdOyyMlIcQPXllFY7hrXa683bZ/CAvvh9NuhPFf8TuNMb6wQmBa1SclgXum5fHZzkM8s2iH33Eir6IEXr4deg+Daf8DIn4nMsYXVghMm647K4uLRmbwX29t+Hy6y24hHIa/3gHVpXDjnyEh1e9ExvjGCoFpk4hw3xfGEVb497+upqvNaNeqD38NW/4OVz0Amaf5ncYYX1khMCeU0zuZ/3flKBZsLGHuyj1+x+m4HR/Dgvtg3Bfh7K/6ncYY31khMO1y2/lDOCO7J/fPW09NfRc+t6C+Bl75GqQPgWm/sX4BY/C4EIjIVBHZKCIFInLcWToiMkhEFojIchFZJSJXe5nHnLpgQPjhVWPYd7iW57vyRenWvQaHd8PVD0Jimt9pjIkKnhUCEQkCDwNXAXnATSKS12y1fwdeVNWzgJnA773KYzruvOF9OHdYb36/cAvVdV30qGDpE9B7OAy7xO8kxkQNL48IJgAFqrpVVeuAOcCMZusocORnWU+gGzRAd2/fuWwkJeW1PLu4Cw4n3bcWdi2C/K9CwFpFjTnCy/8NWcCuJs8L3dea+ilws4gUAvOAu1rakIjcISJLRWRpSYnNoOWnicP6MGlEX/6wcAtVdQ1+xzk5S5+AYAKcaRPQG9OUl4WgpV645mMPbwL+rKrZwNXAbBE5LpOqPqqq+aqan5GR4UFUczK+c3kuByrrmP1JFzoqqK2AlS/A2OsgubffaYyJKl4WgkIgp8nzbI5v+rkdeBFAVT8BEoG+HmYyEXD24N5cNDKDP763hYraLnJUsPolqCuHc273O4kxUcfLQrAEyBWRoSISj9MZPLfZOjuBKQAiMganEFjbTxfwnctyKa2q56mPt/sd5cRUYenj0H8cZJ/jdxpjoo5nhUBVG4A7gbeB9Tijg9aKyM9EZLq72veAr4vISuB54DbtNqeudm9nDUrn0tH9ePT9rZTXRPm0lruXQdFqyP9HO2/AmBZ4OnRCVeep6khVHa6q97mv3auqc93H61T1AlU9Q1XPVNV3vMxjIuvbl+VSVl3Pkx9t9ztK25Y8DvEpcPqX/E5iTFSyMXTmlJ2e3YvLxvTnsQ+2UlYdpUcFVQdh7V+cImAXljOmRVYITId8+7JcDtc08MSH2/yO0rKVc6ChxmkWMsa0yAqB6ZBxWT25cmx/Hv9wG8XlNX7HOZaqc+5A9gS7wqgxbbBCYDrsB1NHU9cQ5r431vsd5VjbP4ADm+1owJgTsEJgOmxYRgrfnDyc11bs4cPN+/2Oc9SSxyEpHcZ+we8kxkQ1KwQmIv558nAG90nmntfWRMdlqsv3wYbXnctJxCX5ncaYqGaFwEREYlyQ/5wxjm37K3nkva1+x4FPH4FwozULGdMOVghMxFw0MoNppw/g4YUFbNtf6V+Q2nJY8hjkTYc+w/3LYUwXYYXARNQ90/JICAa497U1/s1vvOwpqCmDC+72Z//GdDFWCExE9U9L5HtXjOSDzft5fdXezg/QUAeLfg9DLoSsszt//8Z0QVYITMTdct4QTsvqyc9eX8fhzr4O0ZpXnKko7WjAmHazQmAiLhgQ7rtuHPsravnV2xs7b8fhMHz0W+g3FkZc1nn7NaaLs0JgPHF6di9uPXcwTy/aQUFxeefstGA+lKx3jgbsKqPGtJsVAuOZu6bkEhDh5WW7O2eHH/0WeubAuOs7Z3/GdBNWCIxn+qYkcFFuX15bsZtw2OMRRLuWwI6P4Lx/gWCct/syppuxQmA8dd34bPaW1bBo2wFvd/TRbyCxF5x1i7f7MaYbskJgPHX5mP6kJIR4dbmHzUP7N8OGN2DC1yEhxbv9GNNNWSEwnkqKDzJ1XCZvri7y7hpEH/8OQgkw4RvebN+Ybs4KgfHcdWdlUV7bwLvr90V+4+VFsPJ55+JyKRmR374xMcAKgfHcucP6kJmWyF8/86B5aMlj0FjvdBIbY06JFQLjuWBAmHHmQN7bVMKBitrIbbixHj57GnKvsIvLGdMBVghMp7hufBYNYY3s9Yc2vAEV++Cc2yO3TWNikBUC0ylGZ6YxOjOVv0Zy9NDSJ6DnILuchDEdZIXAdJrrx2exYtchtpZUdHxj+wtg23tw9lcgEOz49oyJYVYITKeZfkYWIvDqij0d39iyJyEQshPIjIkAKwSm02T2TOT84X14dfnujk1aU18NK56FMddCav/IBTQmRlkhMJ3qurOy2Xmwis92lp76Rta+CtWlNh+xMRFihcB0qqnjMkmMC3Ss03jpE9An15mFzBjTYVYITKdKSQhxRV4mr6/aS11D+OQ3ULQaCj91jgZszgFjIsIKgel0152VxaGqehZuLD75Ny99AkKJcOZNkQ9mTIzytBCIyFQR2SgiBSLyw1bW+ZKIrBORtSLynJd5THS4MLcvvZLjeGtN0cm9sbYcVr0I474ISenehDMmBoW82rCIBIGHgcuBQmCJiMxV1XVN1skFfgRcoKqlItLPqzwmeoSCAaaM7s/8dUXUN4aJC7bz98iqF6GuAvLtTGJjIsnLI4IJQIGqblXVOmAOMKPZOl8HHlbVUgBVPYW2AtMVXTG2P4drGvh028H2vUHVaRbKPB2yxnsbzpgY42UhyAJ2NXle6L7W1EhgpIh8JCKLRGRqSxsSkTtEZKmILC0pKfEorulMF+VmkBgX4J217WweKlwC+9Y41xWyTmJjIsrLQtDS/9bmZxGFgFxgMnAT8JiI9DruTaqPqmq+quZnZNg157uDpPggF+Vm8M66fe07ueyzpyA+Fcbd4H04Y2KMl4WgEMhp8jwbaH5tgULgNVWtV9VtwEacwmBiwJVjM9lbVsPq3WVtrxgOw6a3YeSVNhWlMR7wshAsAXJFZKiIxAMzgbnN1nkVuARARPriNBVt9TCTiSJTxvQjGBDePlHzUNFKqCyB3Ms7J5gxMcazQqCqDcCdwNvAeuBFVV0rIj8Tkenuam8DB0RkHbAA+L6qHvAqk4kuvZLjmTi0N++sPcEUlpvfde6HT/E+lDExyLPhowCqOg+Y1+y1e5s8VuC77s3EoCvy+vPT/1vH1pIKhmW00uxTMB8GnmVzEhvjETuz2Pjq8rGZALyzrpWjgqqDzoihEdYsZIxXrBAYX2X1SuK0rJ6t9xNsXQAaduYlNsZ4wgqB8d0Vef1ZvvMQxYdrjl+4+V1I6m0nkRnjISsExndXjnOah+avb9Y8FA5Dwbsw/FKbjtIYD1khML7L7ZfCkD7JvN189FDRSqgstmGjxnjMCoHxnYhw5dhMPtmyn8M19UcX2LBRYzqFFQITFa4Y25/6RmXBhibXHbRho8Z0CisEJiqclZNO35SEo8NIbdioMZ3GCoGJCoGAcHlefxZuKKamvrHJsFErBMZ4zQqBiRpXju1PZV0jn2w54A4bTYess/2OZUy3Z4XARI3zhvchNSHE/LV73GGjU2zYqDGdwNNrDRlzMhJCQSbl9mX3+k+hwYaNGtNZ7IjARJVLRvVjXPUS54kNGzWmU7TriEBEEoF/BibhzDL2IfAHVW3hmgDGnLqLR2WwM7iC4tQ8+tmwUWM6RXuPCJ4GxgK/Ax4CxgCzvQplYlf/UBXjAwV8oGf6HcWYmNHePoJRqnpGk+cLRGSlF4FMjNu6gCBh5pSO4oqaelIT4/xOZEy3194jguUicu6RJyIyEfjIm0gmpm1+l/r4XixrHM5HBTZZnTGdob2FYCLwsYhsF5HtwCfAxSKyWkRWeZbOxJZwIxTMJ5g7hR4J8SzcWHzi9xhjOqy9TUNTPU1hDDiXlKgsITDmGibV9mXhxhJUFRHxO5kx3Vq7CoGq7vA6iDGs/z8IxsOIy5lcfYg31xSxcV85ozPT/E5mTLdm5xGY6KAKG96AoRdDYhoXj+wHwMKNJT4HM6b7s0JgokPxOijdBqOvASCzZyJjBqRZP4ExncAKgYkOG94ABEZd/flLk0dlsHR7KeVNJ6sxxkScFQITHdb/H+RMgNT+n780eWQGDWG1YaTGeMwKgfHfoZ1QtApGTzvm5fGD00lNCPHeJmseMsZLVgiM/za84dy7/QNHxAUDTMrty4INzjBSY4w3rBAY/214A/rlQZ/hxy2aPCqDosM1bNxX7kMwY2KDFQLjr8oDsOOj444GjrBhpMZ4zwqB8demt5y5iVspBJk9ExmdmWrDSI3xkBUC468Nr0NaNgxo/bLTk0f1s2GkxnjICoHxT10lbPm7czTQxvWELhllw0iN8ZKnhUBEporIRhEpEJEftrHeDSKiIpLvZR4TZQr+Bg01MGZam6vZMFJjvOVZIRCRIPAwcBWQB9wkInktrJcKfAtY7FUWE6U2vAFJ6TDo/DZXiwsGuHBkX+av20d9Y7iTwhkTO7w8IpgAFKjqVlWtA+YAM1pY7z+B/wJs/uNY0lgPm96EkVdB8MQXwb3x7Bz2V9Txztp9nRDOmNjiZSHIAnY1eV7ovvY5ETkLyFHV19vakIjcISJLRWRpSYkNI+wWdnwENWWtjhZq7qKRGWSnJ/HsYrsiujGR5mUhaKn37/PTQ0UkAPwP8L0TbUhVH1XVfFXNz8jIiGBE45v1r0MoCYZf2q7VgwHhpgmD+HjLAbaUVHgczpjY4mUhKARymjzPBvY0eZ4KjAMWutNfngvMtQ7jGKAKG+fBiCkQn9zut30pP4e4oPDc4p0ehjMm9nhZCJYAuSIyVETigZnA3CMLVbVMVfuq6hBVHQIsAqar6lIPM5losGsxHN593EXmTiQjNYErx2by8rJCauobPQpnTOzxrBCoagNwJ/A2sB54UVXXisjPRGS6V/s1XcCqF5xmoRMMG23JrImDKauu5/VVez0IZkxsau/k9adEVecB85q9dm8r6072MouJEg21sOYvThFISD3pt587rDfDM3rwzKId3HB2tgcBjYk9dmax6Vyb50PNITj9y6f0dhFh1sTBrNh1iDW7yyIczpjYZIXAdK5Vc6BHBgy75JQ38cXx2STGBXjWOo2NiQgrBKbzVJfCprdh3A3tOomsNT2T47j29IG8tmK3XYjOmAiwQmA6z9pXobEOzji1ZqGmbj53MFV1jby6fHcEghkT26wQmM6z6kXoO7LNS0631+nZPRmXlcazi3faNJbGdJAVAtM5SnfAzo+dTuI2LjndXiLCzRMHs6GonGU7SiMQ0JjYZYXAdI7VLzr3p90YsU1OP3MgqQkh6zQ2poOsEBjvqTrNQoPOh/TBEdtscnyI68dn8caqvZSU10Zsu8bEGisExnt7lsP+TRHpJG7uK+cPoT4c5ulPtkd828bECisExnurXoRgPOS1NB1FxwzLSOHKvEye/mQHlbUNEd++MbHACoHxVmMDrHkZRk51ZiPzwDcuHkZZdT1zluw68crGmONYITDe2roAKktO+ZIS7XHWoHQmDO3N4x9staksjTkFVgiMt1bOcY4Ecq/wdDf/dPFw9pTV8H8r95x4ZWPMMawQGO/UljsT1I+9DkLxnu5q8qgMRvVP5ZH3ttoJZsacJCsExjuL/wgN1XDGTZ7vSkS446JhbNxXzsJNNq+1MSfDCoHxxsFt8P6DMGY65EzolF1ee8ZABvRM5JH3tnTK/ozpLqwQmMhThXnfh0AIpj7QabuNDwW4fdJQFm09yIpdhzptv8Z0dVYITOStnwsF8+GSf4OeWZ2665kTBpGWGLKjAmNOghUCE1m15fDmDyDzNJhwR6fvPiUhxC3nDenrB1IAABFkSURBVOattUVs21/Z6fs3piuyQmAia8H9UF4E037ToclnOuK284cSFwzw6Ptbfdm/MV2NFQITOXtXOiOF8r8K2fm+xchITeCGs7N55bNCistrfMthTFdhhcBERrgRXv8OJPeBKff6nYavTRpKXUOYFz61y04YcyJWCExkLPsz7F4GV9zn2TWFTsawjBQmjejL85/upDFsJ5gZ0xYrBKbjKorh3f+AIRfC6V/yO83nbj53EHvKaliwodjvKMZENSsEpuMW3Af1VXDNryMyDWWkTBnTn36pCTyzeIffUYyJalYITMeUboflz8DZt0HGSL/THCMuGGDmhEG8t6mEXQer/I5jTNSyQmA65v0HQYJw4Xf9TtKimefkIMBzn9q8xsa0xgqBOXUHt8KK55zhomkD/U7TooG9kpgypj8vLtlFXYPNVWBMS6wQmFP33n9DMA4mfcfvJG2aNXEQByrreGttkd9RjIlKVgjMqTmwBVbNgXO+BqmZfqdp00W5GeT0TuLZRdZpbExLPC0EIjJVRDaKSIGI/LCF5d8VkXUiskpE/iYig73MYyLovV9CMAEuuNvvJCcUCAj/MGEwi7cdZPO+cr/jGBN1PCsEIhIEHgauAvKAm0Qkr9lqy4F8VT0deBn4L6/ymAgq2QirX4IJX4eUfn6naZcv5WcTFxSeXWydxsY05+URwQSgQFW3qmodMAeY0XQFVV2gqkfG9S0Csj3MYyLlvV9CKKlLHA0c0SclgavGDeCVzwqpqmvwO44xUcXLQpAFNL3QS6H7WmtuB95saYGI3CEiS0VkaUmJTUPoq+L1sOYvMPEb0KOv32lOyqyJgyivaeD1lXv9jmJMVPGyELR0immLF30RkZuBfOC/W1quqo+qar6q5mdkZEQwojlpCx+A+BQ4/y6/k5y0CUN7k9svxc40NqYZLwtBIZDT5Hk2sKf5SiJyGfBvwHRVrfUwj+moojWw7lU495uQ3NvvNCdNRJg1cRCrCstYVWhTWRpzhJeFYAmQKyJDRSQemAnMbbqCiJwFPIJTBOzKYNFMFf72M0hIg/P+xe80p+z6s7NJSQjx07lrqW+0E8yMAQ8Lgao2AHcCbwPrgRdVda2I/ExEprur/TeQArwkIitEZG4rmzN+W/tX2Pw2XPyvUXGZ6VOVlhjHfdeN47Odh3jwnY1+xzEmKng6l6CqzgPmNXvt3iaPL/Ny/yZCqg7Cm/8KA86Eif/kd5oOm3FmFou3HeSR97YycWhvLh3d3+9IxvjKziw2J/bOPU4xmPGQb/MQR9q90/IYnZnK915cyZ5D1X7HMcZXVghM27YsgBXPOOcMZJ7md5qISYwL8vCs8dQ1hPnW88utv8DENCsEpnV1VfD6t6H3cKdvoJsZnpHC/defxtIdpfx6/ia/4xjjGysEpnUL73cmnpn+vxCX5HcaT8w4M4ubJuTwh4VbWLjRBq6Z2GSFwLRsz3L45GFn5rEhk/xO46mfXDuW0ZmpfPfFlRSV1fgdx5hOZ4XAHK+xHubeBT36wWX/4XcazyXGBXnoH8ZTU9/Id15YQTjc4gnwxnRbVgjM8T55CIpWwzUPQlIvv9N0ihH9Urh3Wh6fbD1gl6AwMccKgTlWWSEs/CWMngZjrvU7Taf68jk5XDQyg1/M28DOAzbZvYkdVgjMsd79KaAw9QG/k3Q6EeGB608jFBC+//JKayIyMcMKgTlq16fOhDPnfwt65Zx4/W5oYK8k/n3aGBZvO8hsm9rSxAgrBMYRDsNbP4TUAV1qwhkvfCk/h4tHZvDAmxvYcaDS7zjGeM4KgXGsfgl2L4MpP4GEFL/T+EpEeOCLpxEKCt9/eZU1EZluzwqBgbpKp29g4Hg4/ct+p4kKA3omcc+0PD7ddpCnPtnudxxjPGWFwMBHv4XyPU4HccD+SRxx49nZXDIqg1++tYHt+62JyHRf9r8+1h3a5RSCcV+EQRP9ThNVRIRfXH86ccEAd89Zzv4Km0DPdE9WCGLd39wzh2PgDOJTkdkzkQdvPIMNReVc878fsGT7Qb8jGRNxVghimQ0XbZcrx2by13++gKS4IDMfXcSj729B1TqQTfdhhSAWqULJRmfWMRsu2i55A9OYe9ckrsjrz/3zNvCN2csoq673O5YxEdE9ppsyJ1ZXCdveh83vwOZ3oWwnIHDD4zE/XLS90hLj+P2s8Tz50Xbun7eea3/3Ib+fNZ5xWT39jmZMh0hXO8TNz8/XpUuX+h2j6yjeAG//CLZ/CI11EJ8CwybDiMucmzUJnZJlO0q587nPOFBZx5O3ncMFI/r6HcmYNonIMlXNb3GZFYJubN9aeGo6iMAZM2HE5TDoPAjF+52sWzhQUcusxxaz40AVs2+fQP6Q3n5HMqZVbRUC6yPoropWw5+nQTAe/vFtuOLnMOxiKwIR1Cclgdm3T2RAz0S++uQS1uwu8zuSMafECkF3tHclPHUtxCXDV9+APsP9TtRtZaQm8MzXJpKWFMctjy9m075yvyMZc9KsEHQ3e5Y7RSA+BW57HXoP8ztRtzewVxLPfX0iccEAsx5bbGchmy7HCkF3UrgMnpoBiT3htjeg91C/E8WMwX168OzXJtIYVmY9tpjdh6r9jmRMu1lncVdTstH51d9QAw21R+/rq2DJ45DcG77yuo0G8sma3WXc9KdF9OkRzy3nDSEnPYmc3slkpyeRmhjndzwTw2zUUHfQUAfv/zd88CvQxuOXB+Kg32i4aQ70zO78fOZzy3aU8o3Zy467NlGv5Dhy0pOZPCqDW84bTL/URJ8SmlhkhaCr27MCXv1nKF4LZ9wEk77rnAQWSoRQAgQTIGjnBkYTVaW0qp5dB6soLK1mV2kVhaVVbCmuZNG2A8QFAsw4cyC3XziU0Zlpfsc1MaCtQmDfHtGsoQ4+eNA5Ckju6/zaH3WV36lMO4gIvXvE07tHPGfk9Dpm2fb9lTzx0TZeWlrIS8sKuTC3L7dPGsrFIzMQEZ8Sm1hmRwTRqKEO9nwGb3wP9q2B02fC1F847f+m2zhUVcdzn+7kqY+3s+9wLenJcYzol8KIfikMz0j5/PHAnkkEAlYgTMdY01A0aqiDyhIo3wv7N8P+jVCyybk/uM3pB0jpD9N+A6Ov9jut8VBdQ5h5q/eyaOsBtpRUUFBcQWnV0QvaJYQCZKcnkZ2eTE5v9z49mYG9EkmOD5EQCpAQFyAhFHQehwKEgjYg0BzLt0IgIlOB3wJB4DFVfaDZ8gTgaeBs4ADwZVXd3tY2o64QqDpf6Id2Qul2qNzfZERP9dGRPXWVULEPKoqd+6oDx24nEILewyFjJPQdBX1HwsgrICndl49l/HWgopYtJZUUFFewbX9Fk36Gag5Vnfiqp4N6JzN2YJp768nYgWn0S2u5czocVkSwZqluzpc+AhEJAg8DlwOFwBIRmauq65qsdjtQqqojRGQm8EvAm0lzGxsg3NDysnA9lBW6X+Y74JB7Kyt0vuiPdMo2va855Kx/aKczdLMlgbij74lLhpR+zgleg86FlEzneWom9BkB6UMgaMMLjaNPSgJ9UhKYMPT45sDDNfUUHqym6HA1tfVhahoaqa0PU9sQprahkYraRgqKy1m75zBvrin6/H19UxJISwq56x59T11jmFBAyEhNcG4pCfRLc+7TkuKoawwft5/GcJj0HvH0S00kIzWBfkfem5pAj/gQQWvK6lK87CyeABSo6lYAEZkDzACaFoIZwE/dxy8DD4mIqBeHKYsehvn3tm/dUBL0GuQMwwyEjv7Cryk7+ks/PtX5Ah8+BdIHQ6/BzntSM49++QeCEf8YxqQlxpE3MI68gScebXS4pp71ew6zds9h1u09TG1D+PPmo4RQ0G1SClDbEKakvJaS8lr2lNWwsrCMA5W1NP2fGN/kfcEAlFbWU9cYbnG/cUH5vKkqMc65t+LQcd+aksu1ZwyM+Ha9LARZwK4mzwuB5pPifr6OqjaISBnQB9jfdCURuQO4A2DQoEGnlmbwBTDlJy0vCwQhLcv5Mk8fDD0ynCt2GtPFpSXGMXFYHyYO63PS721oDFNZ20hCXID4YOC4DmtVpay6npLyWorLaykur6GkvJaqukbn6OTIkUdDmJr6RsJdrD8yGvVM8qbVwMtC0NI3afN/Ce1ZB1V9FHgUnD6CU0qTne/cjDHtEgoG6JnceqeziNArOZ5eyfHk9k/txGQm0rwcWlAINL3OQTawp7V1RCQE9ARsdnBjjOlEXhaCJUCuiAwVkXhgJjC32Tpzga+4j28A/u5J/4AxxphWedY05Lb53wm8jTN89AlVXSsiPwOWqupc4HFgtogU4BwJzPQqjzHGmJZ5eokJVZ0HzGv22r1NHtcAN3qZwRhjTNvs9ENjjIlxVgiMMSbGWSEwxpgYZ4XAGGNiXJe7+qiIlAA7TvHtfWl21nKU60p5u1JW6Fp5u1JW6Fp5u1JW6Fjewaqa0dKCLlcIOkJElrZ29b1o1JXydqWs0LXydqWs0LXydqWs4F1eaxoyxpgYZ4XAGGNiXKwVgkf9DnCSulLerpQVulberpQVulberpQVPMobU30ExhhjjhdrRwTGGGOasUJgjDExLmYKgYhMFZGNIlIgIj/0O09bROQJESkWkTV+ZzkREckRkQUisl5E1orI3X5nao2IJIrIpyKy0s36H35nag8RCYrIchF53e8sbRGR7SKyWkRWiMhSv/OciIj0EpGXRWSD++/3PL8ztURERrl/pkduh0Xk2xHdRyz0EYhIENgEXI4zGc4S4CZVXdfmG30iIhcBFcDTqjrO7zxtEZEBwABV/UxEUoFlwBei8c9WRATooaoVIhIHfAjcraqLfI7WJhH5LpAPpKnqNL/ztEZEtgP5qtolTtASkaeAD1T1MXfOlGRVPeR3rra432W7gYmqeqon1h4nVo4IJgAFqrpVVeuAOcAMnzO1SlXfp4vM1Kaqe1X1M/dxObAeZy7qqKOOCvdpnHuL6l9CIpINXAM85neW7kRE0oCLcOZEQVXror0IuKYAWyJZBCB2CkEWsKvJ80Ki9MuqKxORIcBZwGJ/k7TObWZZARQD81U1arO6fgP8KxD2O0g7KPCOiCwTkTv8DnMCw4AS4Em32e0xEenhd6h2mAk8H+mNxkohkBZei+pfgl2NiKQArwDfVtXDfudpjao2quqZOHNoTxCRqG16E5FpQLGqLvM7SztdoKrjgauAf3GbOKNVCBgP/EFVzwIqgWjvO4wHpgMvRXrbsVIICoGcJs+zgT0+Zel23Pb2V4BnVfUvfudpD7cZYCEw1ecobbkAmO62vc8BLhWRZ/yN1DpV3ePeFwN/xWmSjVaFQGGTI8KXcQpDNLsK+ExV90V6w7FSCJYAuSIy1K2qM4G5PmfqFtwO2MeB9ar6a7/ztEVEMkSkl/s4CbgM2OBvqtap6o9UNVtVh+D8m/27qt7sc6wWiUgPd7AAbhPLFUDUjnpT1SJgl4iMcl+aAkTdAIdmbsKDZiHweM7iaKGqDSJyJ/A2EASeUNW1PsdqlYg8D0wG+opIIfATVX3c31StugC4BVjttr0D/NidrzraDACeckdeBIAXVTWqh2R2If2Bvzq/CwgBz6nqW/5GOqG7gGfdH4dbga/6nKdVIpKMM+rxG55sPxaGjxpjjGldrDQNGWOMaYUVAmOMiXFWCIwxJsZZITDGmBhnhcAYY2KcFQJjjIlxVgiMMSbGWSEwpoNE5BwRWeXOd9DDnesgaq9hZExzdkKZMREgIj8HEoEknGvY/MLnSMa0mxUCYyLAvUzBEqAGOF9VG32OZEy7WdOQMZHRG0gBUnGODIzpMuyIwJgIEJG5OJeKHoozdeedPkcypt1i4uqjxnhJRG4FGlT1OffKph+LyKWq+ne/sxnTHnZEYIwxMc76CIwxJsZZITDGmBhnhcAYY2KcFQJjjIlxVgiMMSbGWSEwxpgYZ4XAGGNi3P8HFsSwTISIC68AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"y_new_svi_np = np.array(y_new_svi)\n",
"x_svi_np = np.array(x_svi)\n",
"\n",
"plt.plot(x_svi_np[:,0], y_new_svi_np[:,0])\n",
"plt.plot(x_svi_np[:,0], y_new_svi_np[:,1])\n",
"plt.xlabel('x')\n",
"plt.ylabel('p')\n",
"plt.title('SVI posterior predictive');"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Sample: 100%|██████████| 1100/1100 [01:04, 17.06it/s, step size=1.22e-01, acc. prob=0.895]\n"
]
}
],
"source": [
"from pyro.infer import MCMC, NUTS, Predictive\n",
"\n",
"\n",
"nuts_kernel = NUTS(model)\n",
"\n",
"mcmc = MCMC(nuts_kernel, num_samples=1000, warmup_steps=100)\n",
"mcmc.run(data, y)\n",
"\n",
"hmc_samples = {k: v.detach().cpu().numpy() for k, v in mcmc.get_samples().items()}"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" mean std median 5.0% 95.0% n_eff r_hat\n",
"linear.weight[0,0] -1.27 0.77 -1.31 -2.60 -0.08 317.19 1.00\n",
"linear.weight[0,1] 0.34 0.74 0.36 -0.95 1.48 533.85 1.00\n",
"linear.weight[1,0] 1.35 0.76 1.31 0.03 2.45 302.17 1.00\n",
"linear.weight[1,1] -0.34 0.73 -0.35 -1.49 0.85 530.33 1.00\n",
" linear.bias[0] 3.01 0.81 2.98 1.50 4.24 680.84 1.00\n",
" linear.bias[1] -2.98 0.78 -2.95 -4.16 -1.62 794.84 1.01\n",
"\n",
"Number of divergences: 0\n"
]
}
],
"source": [
"mcmc.summary()\n",
"posterior_samples = mcmc.get_samples()\n",
"posterior_predictive = Predictive(model, posterior_samples, return_sites=(\"linear.weight\", \"obs\", \"_RETURN\"))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# slice through x\n",
"\n",
"x_ = [[x,1] for x in np.linspace(0, 7)]\n",
"x_mcmc = torch.FloatTensor(x_)\n",
"\n",
"samples_mcmc = posterior_predictive(x_mcmc)\n",
"\n",
"y_new_mcmc = []\n",
"\n",
"for i in range(50):\n",
" y_ = torch.bincount(samples_mcmc['obs'][:,i], minlength=2).float()/samples_mcmc['obs'].shape[0]\n",
" y_new_mcmc.append(y_.numpy())"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5xU9bnH8c8zM9tgWRbYpe3S+6I0ERSIooiCDUtUjIkx0ZhY0jS9meTm3iQ3uaZpjMaosSJ2FBSxB1BpAtJdmix9aUvbNvPcP84Bh2V3WWDOnCnP+/Wa10458zvPLMv5zvmdc34/UVWMMcakr4DfBRhjjPGXBYExxqQ5CwJjjElzFgTGGJPmLAiMMSbNWRAYY0yasyAwJoZE5DoRed3vOppKRB4Rkd+49z8nIitPsJ1/iMjPY1udiRcLAnMUEVknItUiUlDn+YUioiLSNeq5YSIyTUR2i8hOEZkjIl9xXxvtLv98nXYGus+/E/WciMi3RGSJiOwXkTIReUZETvX0wx5Z1y9F5PGTaUNVn1DV82NVUzyp6n9Utc+xlhORG0RkZp33fkNV/8u76oyXLAhMQ9YC1x564G6Qc6IXEJEzgbeAd4GeQBvgFmB81GLbgREi0ibquS8Dq+qs7y/At4FvAa2B3sCLwEUx+CxxISKhk3iviMhJ/X88mfWbNKeqdrPbETdgHfAzYG7Uc38Efgoo0NV9biZwbyPtjAbKgH8At7nPBd3nfgG84z7XCwgDw46jxneA3wJzgD3AS0DrqNcvBZYCu91l+0W99kNgI7AXWAmMAcYB1UANsA9Y5C7bEvgXsNl9z2+AoPvaDcAs4E/ATve1G4CZUesaAcx1a5wLjKjzGf7bbeMg0LOBf4sfA8uAXcDDQHad3+8PgS3AY+7zFwML3c8+GxgQ1d5gYIH72Z8GJgG/iW4vatlOwPM4Yb4DuAfoB1S6/177gN3uso9EtbMcuDiqnRBQDgxxH5/h1rUbWASM9vtvPt1vtkdgGvIBkCci/UQkCFwDHO42EZFmwJnAs01o61Hgevf+BTgb6E1Rr4/B2QDNOc4arwe+CnQEaoG/urX1Bp4CvgMUAtOAl0UkU0T6ALcDp6tqC7eedar6GvA/wNOqmquqA911/NttuyfORvR84KaoGoYDa4C2OBv1w0SkNTDVrasNcDcwtc7e0ZeAm4EWwPoGPud1bp09cPaUfhb1WnucPaguwM0iMgR4CPi6u877gSkikiUimTh7WY+573kGuLK+Fbr/5q+4NXUFioBJqroc+Abwvvt7yq/n7U8RtTfp1l6uqgtEpMj9nfzGreF7wHMiUtjAZzdxYEFgGvMYzsZ2LLAC5xvxIa1w/n42H6sRVZ0NtHY3wtfjBEO0Nk1pp776VHWJqu4Hfg5cHRVaU1V1hqrW4OzN5OB8Ow8DWUCJiGSo6jpVXV1f4yLSDqeb6zuqul9Vt+F8+58YtdgmVf2bqtaq6sE6TVwEfKKqj7mvP4Xze7wkaplHVHWp+3pNA5/zHlXdoKo7ccImeiMbAe5S1Sp3/V8D7lfVD1U1rKr/BqpwvoWfAWQAf1bVGlV9FmcvpT7DcAL2++5nr1TVmQ0sW9eTwKXulwWAL7jPAXwRmKaq01Q1oqozgHnAhU1s23jAgsA05jGc/8Q3cPTGexfORqjDcbR1O3AO8EKd13YcRzvRNkTdX4+zkSvA2YAd/natqhF32SJVLcXZU/glsE1EJolIxwba7+K2udk9GL4b5xt22wZqqOuIOqLqLGri++tbZr3b7iHbVbWyTs13HqrXrbmT+56OwEZVjR5psqG9kE7AelWtbUJ9R3B/x8uBS9wwuJTPgqALcFWd+kZxYv/+JkYsCEyDVHU9zkHjC3H6iqNfOwC8TwNdC/V4DLgV59vggTqvvQkUi8jQ4yyxU9T9zjj9++U43U5dDr0gIuIuu9Gt/UlVHeUuo8Dv3UXrDsW7AefbdIGq5ru3PFXtH7VMY8P3HlFHVJ3Re1ZNGf637ueM7larr+b/jqo3X1WbuXsjm4Ei9/cR3V59NgCdGzgA3ZSaD3UPTQCWueFwqN3H6tTXXFV/14Q2jUcsCMyx3Aic63a/1PUD4AYR+f6hfm/31NBJdRdU1bXA2TgHnOu+9gnwd+Ap95TTTBHJFpGJIvKjRmr7ooiUuN86fw08q6phYDJwkYiMEZEM4E6cDfpsEekjIueKSBbOQc+DON1FAFuBrofO3lHVzcDrwP+JSJ6IBESkh4icfaxfmmsa0FtEviAiIRG5BijB6Xs/HreJSLF7zOEnOAd5G/JP4BsiMtw9E6m5iFwkIi1wgrsW+JZbzxU4XUD1mYMTHL9z28gWkZHua1txgjuzkTom4RxPuYXP9gbAOc50iYhcICJBt93RIlJ8jN+B8ZAFgWmUqq5W1XkNvDYbONe9rRGRncADOBvA+pafqaqb6nsN57TRe4B7cc4mWQ1cDrzcSHmP4ZytsgXIdttAVVfi9EX/DWcP4RLgElWtxjk+8Dv3+S043Tw/cdt7xv25Q0QWuPevBzL57KydZ2liN4aq7sA5g+dOnO6vH+CcTVPelPdHeRInkNa4t980ss55OMcJ7nHrLcXp2sP9/Fe4j3fhHEt5voF2wji/t57ApzhnJ13jvvwWzgH/LSJS72dxQ/R9nOMyT0c9vwFnL+EnOGcjbQC+j22LfCVHdhcakxzci9EeV9UH/a7FSyKyDrhJVd/wuxaTuiyFjTEmzVkQGGNMmrOuIWOMSXO2R2CMMWku6QapKigo0K5du/pdhjHGJJX58+eXq2q9Q3kkXRB07dqVefPqPZvRGGNMA0SkoavIrWvIGGPSnQWBMcakOQsCY4xJcxYExhiT5iwIjDEmzXkWBCLykIhsE5ElDbwuIvJXESkVkcXuzErGGGPizMs9gkdw5oFtyHicuWp74UzVd5+HtRhjjGmAZ9cRqOp7ItK1kUUmAI+6syV9ICL5ItLBHb42LlSV/dVh9lXWog3MtZERDJAVCpAVCpIRFI6c08OYFBauhf3boL5haDQC4WqorXRuNe7PcDUEQhDKglD2kT9Vobbqs/fUVjqPIw1NgiZR749qKxiC6v1wcDdU7nFvu6GyouG26taUkeP8lADUVkPtwSNrC9dCKLPOZ8iGYCZEaj5btibqc2gkZr/6BvUZB0WnxbxZPy8oK+LIKfjK3OeOCgIRuRlnr4HOnRuaUKlxUxdv5sk566k4WEtFZQ0VB2uoqKwlHGn6WEsiHA6FNs0zKWqVQ6fWzShulUOnVs7PtnnZZIcCZGUEyQoFCAUsPEwSqD4AW5fClsXObfNi2LbM2cgllYb+r8VjTLU4/D9v0T7lgqC+31q9/1qq+gDOhCcMHTr0hP5Fa8IRDlaHaZObSffC5uRlZ5CXE6JlTgbNs0IE69lYK1AbjlBV695qwlTVRqisCVO+r5oNuw6w5OPN7DrQ0JzjEBDICgVp1SyDs/sUckH/9ozoUUBmyI7TmzioqYTNi6BsDpTNg/3ln32Djf5mfmDHZ99os1tC+wFw+k3QpgdI8Oh2RY7+xh/KgWCG084R35bdmwQ+Wy76W3Yg6LRXl0bcb+vR9VZBuAoycyEn36k1uyVk50NWnrO3UJ9IuIG9kfDRewmhbGcPIlxn3TWVzrqDmfV87kwIJO//aT+DoIwj52It5si5WGPqssFFXDa46NgLnoB9VbWU7TpA2c6DlO+rcoMjTFVN5PD9jbsPMmXhJp6as4EWWSHO7deWC/q35+zehTTPChGO6BHvqa6N0CE/m4xg8v5xmRhY/77zDb1kgvNt8FhqKqF0hvO+DR867w1XO6/ld4aWnZ0NZ92NX/NCaH+qEwD5nevfMCezQBAymzm3pgplObc04GcQTAFud+e3HQ7siefxgVjKzQrRt30efdvnNbpcZU2Y2avLeW3JFt5Yvo2XFm4iGBAEqK2niyovO8SYfu0OB0ZOZj3fzExq2rsVZvwcFruzPE7/CfS7xPmW3mXkkRtqVdi8ED56HD5+xukzD2VDxyFwxi1QPAyKT4cW7fz5LCbheTYfgYg8BYwGCnAmu74LyABQ1X+I03F+D86ZRQeArzQ0N260oUOHaioMOlcbjjB33S5mlZajKFmhoHv8wTm+EBCYs3YXb67Yyu4DNWRnBDirl9O1NKpXAW1bZNmxh1QUroV5/4K3fuN0R4z4FvS/HBY9BR895mzkC/s6gdBrLKx81QmArUucjX+/S2HwdU5YBDP8/jQmgYjIfFUdWu9ryTYxTaoEQVPVhCPMWbuT6Uu38PrSrWypcA7etcgK0aNtLj0Kc+nZ1rkNKG5Ju7xsnys2J2zDXJj6XdjyMXQ/By78IxT0/Oz16gOw9HmY809nD+CQjkNg8BfhlCudfnNj6mFBkCIiEWXxxj0s2rCb0m37KN22j9Xb97FtbxUAwYDw+SHF3H5uTzq1Po6+UOOv2iqY/lOY+09o0RHG/Q+UXNZ4P/3G+bBuJvQ8D9r1j1+tJmlZEKS4PQdrWL19H1MWbuLJDz9FUa4e2onbzulJx/wcv8szjanYDJO/BGVz4Yxb4ZyfQFYLv6syKciCII1s3nOQe98u5em5GxCEa4d14tZzelqXUSJa/z5Mvt65OOqyv0P/y/yuyKQwC4I0VLbrAPe+Xcoz88oIBYW/TBzMBf2bcPqh8Z4qzH0QXvsR5HeBiU9A235+V2VSXGNBYCepp6jiVs347RUDeOvO0fRtn8ctj8/nsQ8anKnOxEtNJbx0O0z7HvQYA197y0LA+M6CIMV1btOMp752Buf0acvPX1zCH6avINn2AlNG+Sfw0AWw8HE4+4dw7SQ7y8ckBAuCNJCTGeT+L53GtcM6ce/bq/neM4upCcdhgCzjONQV9I/Pwe71MPFJ56BwEg9JYFKLn1cWmzgKBQP8z+Wn0j4vhz+9sYrt+6r4+3VDyM2yPwFP7d0KU74Jn0yHHufChL9DXge/qzLmCPaVJI2ICN8+rxe/v/JUZpWWM/GB9ynfV+V3WalrxVS470xY+y6M/1+47jkLAZOQLAjS0DWnd+af159G6bZ9XPvAB2zbm2xDDSe4yj3OXsCkL0BeEdz8Lgz/unUFmYRlf5lp6ty+7XjkK8PYuPugEwYVFgYxsfwVuHe4M/7PqDvgpjehbV+/qzKmURYEaeyM7m145CvD2LKnkokPfMCWPRYGJ6xiM0y6Dp6+Dpq1gZvegPPucma5MibBWRCkuWHdWvPvrw5j294qJj7wPpv3HPS7pOQSicDcf8G9w6D0DRhzF9z8jiezSBnjFQsCw9CuThjs2FfNNfd/wMbdFgZNsqcMHrkQpt4BHQfBLbPhc3fY8M8m6VgQGABO69KKx24azq4D1Vxz//ts2HnA75IS2/5yePQy2LIEJtwL109xpnU0JglZEJjDBnXK54mbhlNxsIYJ985i9upyv0tKTJUV8PiVsGcDXDfZmQvAJgkyScyCwBxhQHE+L9w2klbNMvjSv+bw4H/W2JAU0WoqndNCty6Bqx+FLiP8rsiYk2ZBYI7SozCXF28byXn92vKbqcv5ztMLOVgd9rss/4Vr4bkbYd1/4LL7oPcFfldkTExYEJh6tcjO4L7rTuP7F/RhyqJNXHHf7PQ+bqAKL38bVrziXCU84Gq/KzImZiwITIMCAeG2c3ry0A2ns3HXAS65Z2Z6HjdQhdd/5o4a+iPnKmFjUogFgTmmc/q0Zcrto2jdPJM7nl5EOJJmxwzmPADv3wPDbobRP/K7GmNizoLANEnXgubcObYPWyoq02uvYNc6mHEX9B4H435vZweZlGRBYJpsTL+25GWHeG5+md+lxIcqvHIHBEJw0d02aJxJWfaXbZosOyPIxQM78trSLeytrPG7HO8teQ5Wvwljfg4ti/yuxhjPWBCY43LlkGIqayK8umSL36V468BOZ3L5otPg9Jv8rsYYT1kQmOMypHM+3Qqap3730Bt3OWFw8Z8hEPS7GmM8ZUFgjouIcMXgIj5cuzN1rytYPxsWPApn3godBvhdjTGesyAwx+3yIU5/+QsfbfS5Eg/UVjkXjrXsDKN/7Hc1xsSFBYE5bsWtmnFm9zY8v6As9cYhmvUXKF8FF98Nmc39rsaYuLAgMCfkiiFFrNtxgPnrd/ldSuyUl8J7f4T+V0CvsX5XY0zcWBCYEzL+1A7kZAR5bkGKdA/VHIQXb4GMbBj3O7+rMSauPA0CERknIitFpFREjro2X0Q6i8jbIvKRiCwWkQu9rMfETm5WiPGntOeVxZuorEnykUnDNfDMDVA2Fy75K7Ro53dFxsSVZ0EgIkHgXmA8UAJcKyIldRb7GTBZVQcDE4G/e1WPib0rhhSzt7KWGcu2+l3KiYtE4KXbYNVrcNH/Qf/L/K7ImLjzco9gGFCqqmtUtRqYBEyos4wCee79lsAmD+sxMXZmjzZ0aJnN8wuS9JoCVZj+E1j8NJz7Mzj9Rr8rMsYXXgZBEbAh6nGZ+1y0XwJfFJEyYBrwzfoaEpGbRWSeiMzbvn27F7WaExAMCJcPLuK9T8rZtrfS73KO33t/gA/vgzNuhc99z+9qjPGNl0FQ3zCNdc81vBZ4RFWLgQuBx0TkqJpU9QFVHaqqQwsLCz0o1ZyoK4YUE44oL32UZDtzc/4Jb/83DJgI5/+3jSpq0pqXQVAGdIp6XMzRXT83ApMBVPV9IBso8LAmE2M92+YysFM+z8zfQG044nc5TfPxszDt+9B7PEy4x0YVNWnPy/8Bc4FeItJNRDJxDgZPqbPMp8AYABHphxME1veTZL46siurtu7ju5MXJX4YVGxyThPtfCZc9TAEM/yuyBjfhbxqWFVrReR2YDoQBB5S1aUi8mtgnqpOAe4E/iki38XpNrpBU+5S1dQ3YVARm3ZX8vvXVhAKCH+8aiDBQIJ2tXxwH0Rq4fL7ICPH72qMSQieBQGAqk7DOQgc/dwvou4vA0Z6WYOJj1tG9yAcifDH11cRDAj/e+UAAokWBpV7YN7D0P9yaNXV72qMSRieBoFJL7ef24uasPKXNz8hFBD+5/JTEysM5j0E1Xth5Lf9rsSYhGJBYGLqO+f1ojYS4d63VxMKCv814RQkEc7Iqa1yuoW6nwMdBvpdjTEJxYLAxJSI8L3z+1AbUe5/dw2hQIC7LinxPwwWPw37tsLl9/tbhzEJyILAxJyI8KNxfampVR6atZbzS9oxoqePZwVHIjDrr9B+AHQf7V8dxiQoO4HaeEJE+P4FfcjJCPo/v/HKabDjExj1HbtwzJh6WBAYz+RkBhndp5DpS7cQifh0VrAqzPoz5HeBfnWHujLGgAWB8di4U9qzbW8VH23waQKbT993hpce8U0IWk+oMfWxIDCeOqdvWzKCwmt+dQ/N+gs0awODrvNn/cYkAQsC46m87AxG9izgtaVb4j+/8bblzjwDw74Omc3iu25jkogFgfHcuP7t2bDzIMs2V8R3xbP/BhnNYNjX4rteY5KMBYHx3NiSdgQEpseze6hiMyyeDEOuh2at47deY5KQBYHxXJvcLIZ1a81rS+MYBMunQKQGTr8pfus0JklZEJi4GNe/Pau27mP19n3xWeHKaVDQGwp6xWd9xiQxCwITF+f3bw/A9HjsFVTugXUzoc9479dlTAqwIDBx0TE/h4Gd8uNznKD0DWfOgT4Xer8uY1KABYGJm3H927OobA8bdx/0dkUrX3WuHSg+3dv1GJMiLAhM3FzQvx0Ar3vZPRSugU9eh97jIBD0bj3GpBALAhM33Qtz6dOuhbdXGX/6vnOMwI4PGNNkFgQmri44pT1z1+2kfF+VNytY+SoEs5wJaIwxTWJBYOJq/CntiSi8sWxr7BtXhRVTofvZkJUb+/aNSVEWBCau+rZvQZc2zby5uGzbcti93s4WMuY4WRCYuBIRxvVvz6zScioqa2Lb+Mppzs/e42LbrjEpzoLAxN0Fp7SnJqy8vGhTbBte+Sp0HAJ5HWLbrjEpzoLAxN2g4nyGdmnF715dEbtrCvZuhY3zrFvImBNgQWDiLhAQ/nTNICIR5c7JC2MzjeWq15yfdtqoMcfNgsD4olPrZtx1aX8+WLOTf81ce/INrnwVWnaGdv1Pvi1j0owFgfHNVacVc35JO/4wfSXLT2bSmuoDsOZtZ29AJHYFGpMmLAiMb0SE315xKnk5GXz36YVU1oRPrKE170BtpXULGXOCLAiMr9rkZvGHzw9gxZa9/N/rK0+skZXTIKsldB0V2+KMSRMWBMZ35/Rty3XDO/PgzLXMXl1+fG+OhJ0Dxb3Og2CGNwUak+I8DQIRGSciK0WkVER+1MAyV4vIMhFZKiJPelmPSVw/vagfXds053uTF7Hn4HFcaLZxPuzfbqeNGnMSPAsCEQkC9wLjgRLgWhEpqbNML+DHwEhV7Q98x6t6TGJrlhniT9cMYuveKv40Y1XT37hiKgRC0HOMd8UZk+K83CMYBpSq6hpVrQYmARPqLPM14F5V3QWgqts8rMckuEGd8jm3b1veWL4V1SZeW7ByGnQZCTmtvC3OmBTmZRAUARuiHpe5z0XrDfQWkVki8oGI1DtIjIjcLCLzRGTe9u3bPSrXJILP9SqgbNdB1u84cOyFyz+B8lXQ92LvCzMmhXkZBPWd0F33a14I6AWMBq4FHhSR/KPepPqAqg5V1aGFhYUxL9QkjlE9CwD4T2kTDhqvmOr87GvHB4w5GV4GQRnQKepxMVB3lLEy4CVVrVHVtcBKnGAwaapbQXOK8nOY9UkTg6DDQGhZ7H1hxqQwL4NgLtBLRLqJSCYwEZhSZ5kXgXMARKQAp6tojYc1mQQnIozqWcDs1eWEGxuDaO9WKJtr3ULGxIBnQaCqtcDtwHRgOTBZVZeKyK9F5FJ3senADhFZBrwNfF9Vd3hVk0kOI3sVUFFZy+Ky3Q0vtOpVQKHvRXGry5hUFfKycVWdBkyr89wvou4rcId7MwaAkT3aADDzk3IGd27gbKAVUyG/C7Qtqf91Y0yT2ZXFJuG0yc2if8c8ZjZ0wLhqrzO+UN+LbZA5Y2LAgsAkpFG9Cljw6S72V9Ue/WLpmxCutm4hY2LEgsAkpFE9C6gJK3PW7jz6xRVTIac1dBoe/8KMSUEWBCYhnd61NZmhAP+pexppuAY+me4MOR309BCXMWnDgsAkpOyMIMO6tmZW3eME62dB5R4bZM6YGLIgMAlrVK8CVm7dy7aKys+eXDEVQjnQ41z/CjMmxVgQmIR1aLiJw2cPqcKKaU4IZDbzsTJjUosFgUlYJR3yaN08k5mHjhNsXgQVZTa2kDExZkFgElYgIIzo0YaZpeXOsNQrpoIEoHe9g9QaY06QBYFJaJ/rVcC2vVV8sm2fM/dA5zOheYHfZRmTUiwITEIb6R4n+GjxIti6xM4WMsYDFgQmoRW3aka3guaEl73iPGHHB4yJOQsCk/BG9Syg287/ECnsC627+12OMSmnSUEgItkicoeIPC8iz4nId0Uk2+vijAE4q3suQ1jJ1oIRfpdiTEpq6h7Bo0B/4G/APUA/4DGvijIm2ojQJ2RJDbMZ4HcpxqSkpg7W0kdVB0Y9fltEFnlRkDF1NS97jxpCPL29M1f6XYwxKaipewQficgZhx6IyHBgljclGVPHmrfZ2nIgczdWsm1v5bGXN8Ycl6YGwXBgtoisE5F1wPvA2SLysYgs9qw6Y/Zthy0fk9n7XFThzeXb/K7ImJTT1K4hu5TT+GPtuwAUDhxHp6UVvL50C9cO6+xzUcakliYFgaqu97oQY+q15m3Ibol0HMz5JSt57IP17KuqJTfL5iIwJlbsOgKTuFRh9TvQ7SwIBBlb0o7q2gjvrdrud2XGpBQLApO4dpQ6o412PweAoV1a0apZBq8v3eJzYcakFgsCk7jWvOP87OEEQSgYYEy/dry1Yhs14Yh/dRmTYiwITOJa/TbkdzliWImxJe2oqKytf1J7Y8wJsSAwiSlcC+v+A91HH/H0Wb0Kyc4IWPeQMTFkQWAS08b5UFVxuFvokJzMIJ/rVciMZVudyWqMMSfNgsAkpjXvAALdzj7qpbEl7di0p5KlmyriXpYxqciCwCSmNW9Dx0HQrPVRL43p25aAYN1DxsSIBYFJPFV7oWzuUccHDmmTm8XQrq15fdnWuJZlTKqyIDCJZ91MiNQevn6gPueXtGPFlr1s2HkgjoUZk5osCEziWfMOhHKg8xkNLjK2pB2A7RUYEwOeBoGIjBORlSJSKiI/amS5z4uIishQL+sxSWL129BlBISyGlykS5vm9G3fwo4TGBMDngWBiASBe4HxQAlwrYiU1LNcC+BbwIde1WKSyJ6NUL6yweMD0caWtGPuup3s3F/teVnGpDIv9wiGAaWqukZVq4FJwIR6lvsv4H8Bm3HEHDWsRGPOL2lPROGtFTZHgTEnw8sgKAI2RD0uc587TEQGA51U9ZXGGhKRm0VknojM277dRp5MaWvegeaF0Lb/MRc9pSiPDi2zrXvImJPkZRBIPc8dvhRURALAn4A7j9WQqj6gqkNVdWhhYWEMSzQJZedaWDEVeoyBwLH/NEWEsSXteO+T7RysDsehQGNSk5dBUAZ0inpcDGyKetwCOAV4x53+8gxgih0wTlORMLzwDQgE4dyfNfltY0vaUVkTYVZpuYfFGZPavAyCuUAvEekmIpnARGDKoRdVdY+qFqhqV1XtCnwAXKqq8zysySSqmX+CDR/AhX+E/E7HXt41vFsbcrNCvLHcTiM15kR5FgSqWgvcDkwHlgOTVXWpiPxaRC71ar0mCW36CN75LfS/HAZcfVxvzQwFOLtPIW8s30YkYoPQGXMiPJ34VVWnAdPqPPeLBpYd7WUtJkHVHITnb4bmbeGiu0HqO7TUuLH92jF18WYWle1mcOdWHhRpTGqzK4uNv2bcBeWr4LK/1zvAXFOc06ctwYBY95AxJ8iCwPin9E2Ycz8Mv6VJ1w00pGWzDIZ1bc0by+x6AmNOhAWB8ceBnfDirVDYF86766SbO6+kHSu37uXTHTYInTHHy4LA+OOV78KBHXDFA5CRc9LNndevLQAzrHvImONmQWDib9tyWPYinPV96DAwJk12adOc3u1yecNGIzXmuFkQmPhb/jIgcNqXY9rsef3aMWfdTvYcqOAykOYAABQqSURBVIlpu8akOgsCE3/Lp0Cn4dCifUybPa+kHeGI8s4qO2hszPGwIDDxtXMtbPkYSmJ/TeGg4nwKcjOZYd1DxhwXCwITXyvcgWb7XhzzpgMBYUzfdry7cjvVtZGYt29MqrIgMPG1bIpzgLhVF0+aP6+kHXurapmzdqcn7RuTiiwITPxUbIayOdDvEs9WMapnAdkZAbvK2JjjYEFg4udQt1A/78YczMkMMqpnITOWbUXVBqEzpiksCEz8LH8ZCnpDYR9PVzO2pC0bdx9k+ea9nq7HmFRhQWDi48BOWDfT026hQ87t2w4RrHvImCayIDDxsfJV0HBcgqCwRRaDOuVbEBjTRBYEJj6WvwwtO0GHQXFZ3Xn92rG4bA9b9lTGZX3GJDMLAuO9qr2w+i1nb+AEJp45EeNPca5anjT307isz5hkZkFgvPfJ6xCuiku30CHdC3M5r19b/j17HQerw3FbrzHJyILAeG/5y9C80BlfKI6+fnYPdh2o4Zn5G+K6XmOSjQWB8VZNJax6HfpeBIFgXFc9tEsrhnTO54H31lAbtiEnjGmIBYHx1pq3oWZ/XLuFDhERvnF2D8p2HWTaki1xX78xycKCwHhr+cuQ1RK6nuXL6s/r147uhc25/93VdqWxMQ2wIDDeCdfAymnQZxyEMn0pIRAQvn5Wd5ZuqmBW6Q5fajAm0VkQGO/MewgO7vKlWyjaZYOLaNsii3+8u9rXOoxJVBYExhsLHoVXfwC9LoDe430tJSsU5KujujGztJwlG/f4WosxiciCwMTewidhyregxxi4+lEIhvyuiC8M70xuVoj731vjdynGJBwLAhNbiyfDi7dC97Nh4hOQke13RQDkZWdw3fDOTF28iU93HPC7HGMSigWBiZ0lz8ELX4euo2DiU5CR43dFR/jqqG4EA8KDM22vwJhoFgQmNpa9BM99DTqdAV94GjKb+V3RUdrlZXP54CImz9vAjn1VfpdjTMKwIDAnb8278OxXoXgoXDcZMpv7XVGDbj6rO5U1EX718jIiEbuuwBjwOAhEZJyIrBSRUhH5UT2v3yEiy0RksYi8KSLezGhuvDXzT5DbHq57FrJa+F1No3q2bcH3zu/NlEWb+PUry+wiM2PwMAhEJAjcC4wHSoBrRaSkzmIfAUNVdQDwLPC/XtVjPFKxCda8A4O+ANl5flfTJLed05OvjuzGI7PX8be3Sv0uxxjfeXle3zCgVFXXAIjIJGACsOzQAqr6dtTyHwBf9LAe44XFTwMKAyf6XUmTiQg/u6gfuw9Uc/eMVbRqlsGXzuzqd1nG+MbLICgCosf/LQMaG4f4RuDV+l4QkZuBmwE6d+4cq/rMyVKFhU85w0u36eF3NcclEBB+//kBVFTW8IspS8nLyWDCoCK/yzLGF14eI6hvKqp6O2RF5IvAUOAP9b2uqg+o6lBVHVpYWBjDEs1J2fQRlK+Egdf6XckJyQgGuOcLQzi9a2vunLyIt1du87skY3zhZRCUAZ2iHhcDm+ouJCLnAT8FLlVVO6cvmSx6CoJZ0P9yvys5YdkZQR788lD6tG/BLY/PZ/76nX6XZEzceRkEc4FeItJNRDKBicCU6AVEZDBwP04I2NexZFJbDR8/C30vhJx8v6s5KXnZGfz7q8Non5fNN5/8iP1VtX6XZExceRYEqloL3A5MB5YDk1V1qYj8WkQudRf7A5ALPCMiC0VkSgPNmUTzyetwcCcM/ILflcREQW4W/3f1QDbtqeTuGav8LseYuPJ0NDBVnQZMq/PcL6Lun+fl+o2HFj0FzdtCj3P9riRmTuvSmuuGd+bhWWu5bFARpxa39LskY+LCriw2x2//Dlg1HQZcnRAji8bSD8b1pU1uFj9+YbHNc2zShgWBOX5LnoNITdKeLdSYljkZ/PKS/izZWMEjs9f5XY4xcWFBYI7foieh/anQ/hS/K/HEhae259y+bbl7xio27j7odznGeM6CwByfbSuc6wdScG/gEBHh1xP6owq/eHGJjUdkUp4FgTk+i54CCcKpV/ldiaeKWzXjzvN78+aKbby2ZIvf5RjjKQsC03SRsDO2UK+xkNvW72o8d8OIrvTvmMddU5ZSUVnjdznGeMaCwDTdmndg7+akGmDuZISCAX57xamU76vi96+u8LscYzxjQWCa5uBueO8PkN0Seo/3u5q4GVCcz42juvHEh5/y8Ky1fpdjjCdS6yRw441ty2HSdbB7PVz6t4SZkD5efjiuL+t3HOBXLy8jv1kGlw8u9rskY2LK9ghM45a9BP8cA1V74cuvOBPQpJlQMMBfrx3MiB5t+N4zi3lrxVa/SzImpiwITP0iYXjjVzD5emhXAl9/F7qc6XdVvsnOCPLA9UPp3zGPWx5fwJy1NkqpSR0WBOZoB3bCE1fBzLvhtBvghqmQ19HvqnyXmxXi4RtOp6hVDjc+Mpelm/b4XZIxMWFBYI60pwweHAPr/gOX/MW5hbL8riphtMnN4vEbh9MiO8SXH5rD2vL9fpdkzEmzIDCf2f0pPHwh7C+HL7/s7A2Yo3TMz+HRG4cTUfjigx/aZDYm6VkQGMeudfDwRVC5G65/ETqf4XdFCa1n21z+/ZVhqCqf/8f7/PzFJey1i85MkrIgMLBzLTxyMVRVwPUvQdFpfleUFE4tbsmMO87mhhFdefzD9Yy9+z1eX2rDUZjkY0GQ7nashkcugup98OUp0HGw3xUlleZZIe66pD8v3DqS/GYZ3PzYfG55fD7bKir9Ls2YJrMgSGflpU4I1FY6xwQ6DPS7oqQ1qFM+L39zFD8Y14e3VmxjzN3v8tz8Mhu51CQFC4J0tWWJEwLhGudCsfan+l1R0ssIBrh1dE9e+85Z9Gufx53PLOK2Jxewa3+136UZ0ygLgnS09j14eDxIAG54xblgzMRMt4LmPHXzGfxwXF9mLNvKBX9+j/dWbfe7LGMaZEGQbj5+Fh6/0rlA7KYZ0Laf3xWlpGBAuGV0D164dSQtczK4/qE5/HLKUiprwn6XZsxRLAjSyex74LkboWgofPU1aGmDp3ntlKKWvPzNUdwwoiuPzF7HxX+bybJNFX6XZcwRLAjSQSQCr/0EXv8plEyAL70AOa38riptZGcE+eWl/XnsxmFUHKzhivtm8dLCjX6XZcxhFgSprnq/sxfwwb0w7Ovw+YfTbhjpRPG5XoVM/dbnGFCUz7cnLeQ3ryyjNhzxuyxjLAhSVvUBpyvozwNg6fMw9tcw/vcQCPpdWVorbJHFE18bzg0juvLgzLV86V9z2LGvyu+yTJqzIEg1NZXwwT/gr4OcrqD2p8CNM2Dkt0HE7+oMzmmmv7y0P3+8aiDzP93FJX+bycdlNpKp8Y8FQao4uBvmPgh/HQyv/RDa9IIbpjlDRnQa5nd1ph6fP62Y574xAhHhyn/M5okP1xOJ2AVoJv4k2a58HDp0qM6bN8/vMvxzYCdsXgTlq2D7Sudn+SrY586a1ekMOOcn0O0s2wNIEjv2VfGtSR8xq3QHQzrn81+XnUL/ji2P+Z6FG3ZzerfW5GVnxKlSk8xEZL6qDq33NQuCJLBnI6yYCsunwPrZoO656Fl5UNAbCvs4P4tPhy4jLACSkKry3IKN/HbacnYdqOb6M7tyx/m9j9jIqyoLPt3N4x+sZ+rizVSHI2RnBLjwlA5cNbQTw7u1JhCwf3tTPwuCRLV3i3OVr0acyV9C2Z/9lCCsnwUrXoGN853lC3pD34uh+9lQ2Bdy29lGP8XsOVDDH19fyeMfrqdN8yx+elFfzi9pz5RFm3js/fUs21xBblaIK4cUMbpPW2Ys38rLCzext6qWzq2bcdVpxVx5WjEd83P8/igmwVgQJJIdq52N+/JXoGzOsZfvOBj6XQJ9L4HC3t7XZxLCx2V7+NlLS1i0YTcZQaEmrPRt34IvndmFywYV0TwrdHjZg9Vhpi/dwuR5G5i9egci0LttC4Z0yWdw51YM6dyK7gXNbW8hzfkWBCIyDvgLEAQeVNXf1Xk9C3gUOA3YAVyjqusaazNhg6C2ypnmcd82ZzTP2qojf+4pc7p3ti11lu8w0Nm4974AsnKPXr62CtqWQH4nfz+X8U0kokyet4Elm/Zw2aAiTuvSCjnGHuCGnQeYsmgTc9ftZMH6XVRU1gLQMieDwZ3z6ZifQ1YoQFYoSFYoQHaG87NNbib9O+bRrSCXoAVGSvIlCEQkCKwCxgJlwFzgWlVdFrXMrcAAVf2GiEwELlfVaxpr94SD4JMZsPRFyG4JOfnOz+yWkJ0PWS3qP79eFcLV9WykK50N/u71zvSOu9bD3s1AY79Lcfrv+14MfS+CVl2O/zMYcxwiEWVN+T4WrN/Ngk938dGnu9mxv5qq2jBVtRGqa4++mC0nI0jfDi3o3zGP/h1b0r2geb3BIAKZwSBZGYHDwZKdESAzFKA2olTVRKiscdZzaH0Ch5fLcgMoKxQgI1j/yYsRVaprI1TWfNZGVU2E6nCE3KwQLXMyyMsJkZMRPGZA1iccUcIRJSMoTXq/qlIdjpARCCTl3lVjQRCq78kYGQaUquoat4hJwARgWdQyE4BfuvefBe4REVEv0mn3p7D6LajcAzUxmHBcApBXBPldoPtoZ8Oe39npt8/Iierzd/v9s/IgO+/k12tMEwUCQs+2LejZtgVXn370nmUk4mzYKmvCbNpdybLNFSzdtIelmyp46aNNPP7Bpz5UffwygkJedgYtskP1horibPSr3GA6FFC17qm6Ihyxl5SVESAjEHB/N5F6gzMzGDi8bFbICcRgHI7XfWtMLy4Z2DHm7XoZBEXAhqjHZcDwhpZR1VoR2QO0AcqjFxKRm4GbATp37nxi1Zx+o3MDZwz+ygpnft7K3c79hr7NB+scxA1lObecVhC00/ZM8goEhOxAkOyMIPnNMinpmMfnT3MGIoxElA27DrBh50G0nv8bEYWa2sjhb/zRG8xQQA5/4z/U9ZQZcjbQhzbGVTWfva8mUv8wG4I4ew9RG+isUJBQUDhQFWbPwRoqKmuoOFjj3q8l3EBbocCRG+5D7QYEZ6+jNvJZbbXOXkfdcDh0vyYcOfwZKmvDh3/G43hryxxvtjleBkF98Vj3N9WUZVDVB4AHwOkaOunKghnQvI1zM8YcJRAQurRpTpc2zf0uxcSBl1cWlwHR+6PFwKaGlhGRENAS2OlhTcYYY+rwMgjmAr1EpJuIZAITgSl1lpkCfNm9/3ngLU+ODxhjjGmQZ11Dbp//7cB0nNNHH1LVpSLya2Ceqk4B/gU8JiKlOHsCE72qxxhjTP28PEaAqk4DptV57hdR9yuBq7yswRhjTONs9FFjjElzFgTGGJPmLAiMMSbNWRAYY0yaS7rRR0VkO7D+BN9eQJ2rlhNcMtWbTLVCctWbTLVCctWbTLXCydXbRVUL63sh6YLgZIjIvIYGXUpEyVRvMtUKyVVvMtUKyVVvMtUK3tVrXUPGGJPmLAiMMSbNpVsQPOB3AccpmepNplohuepNplohuepNplrBo3rT6hiBMcaYo6XbHoExxpg6LAiMMSbNpU0QiMg4EVkpIqUi8iO/62mMiDwkIttEZInftRyLiHQSkbdFZLmILBWRb/tdU0NEJFtE5ojIIrfWX/ldU1OISFBEPhKRV/yupTEisk5EPhaRhSJyAhOLx5eI5IvIsyKywv37PdPvmuojIn3c3+mhW4WIfCem60iHYwQiEgRWAWNxJsOZC1yrqssafaNPROQsYB/wqKqe4nc9jRGRDkAHVV0gIi2A+cBlifi7FWeG8uaquk9EMoCZwLdV9QOfS2uUiNwBDAXyVPViv+tpiIisA4aqalJcoCUi/wb+o6oPunOmNFPV3X7X1Rh3W7YRGK6qJ3ph7VHSZY9gGFCqqmtUtRqYBEzwuaYGqep7JMlMbaq6WVUXuPf3Astx5qJOOOrY5z7McG8J/U1IRIqBi4AH/a4llYhIHnAWzpwoqGp1ooeAawywOpYhAOkTBEXAhqjHZSToxiqZiUhXYDDwob+VNMztZlkIbANmqGrC1ur6M/ADoP5Z2ROLAq+LyHwRudnvYo6hO7AdeNjtdntQRJJhguaJwFOxbjRdgkDqeS6hvwkmGxHJBZ4DvqOqFX7X0xBVDavqIJw5tIeJSMJ2vYnIxcA2VZ3vdy1NNFJVhwDjgdvcLs5EFQKGAPep6mBgP5Doxw4zgUuBZ2LddroEQRnQKepxMbDJp1pSjtvf/hzwhKo+73c9TeF2A7wDjPO5lMaMBC51+94nAeeKyOP+ltQwVd3k/twGvIDTJZuoyoCyqD3CZ3GCIZGNBxao6tZYN5wuQTAX6CUi3dxUnQhM8bmmlOAegP0XsFxV7/a7nsaISKGI5Lv3c4DzgBX+VtUwVf2xqharalecv9m3VPWLPpdVLxFp7p4sgNvFcj6QsGe9qeoWYIOI9HGfGgMk3AkOdVyLB91C4PGcxYlCVWtF5HZgOhAEHlLVpT6X1SAReQoYDRSISBlwl6r+y9+qGjQS+BLwsdv3DvATd77qRNMB+Ld75kUAmKyqCX1KZhJpB7zgfC8gBDypqq/5W9IxfRN4wv1yuAb4is/1NEhEmuGc9fh1T9pPh9NHjTHGNCxduoaMMcY0wILAGGPSnAWBMcakOQsCY4xJcxYExhiT5iwIjDEmzVkQGGNMmrMgMOYkicjpIrLYne+guTvXQcKOYWRMXXZBmTExICK/AbKBHJwxbH7rc0nGNJkFgTEx4A5TMBeoBEaoatjnkoxpMusaMiY2WgO5QAucPQNjkobtERgTAyIyBWeo6G44U3fe7nNJxjRZWow+aoyXROR6oFZVn3RHNp0tIueq6lt+12ZMU9gegTHGpDk7RmCMMWnOgsAYY9KcBYExxqQ5CwJjjElzFgTGGJPmLAiMMSbNWRAYY0ya+380b6S649Lz3wAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"y_new_mcmc_np = np.array(y_new_mcmc)\n",
"x_mcmc_np = np.array(x_mcmc)\n",
"\n",
"plt.plot(x_mcmc_np[:,0], y_new_mcmc_np[:,0])\n",
"plt.plot(x_mcmc_np[:,0], y_new_mcmc_np[:,1])\n",
"plt.xlabel('x')\n",
"plt.ylabel('p')\n",
"plt.title('MCMC posterior predictive');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Comparing Posterior Distributions"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deXxb1ZX4v0eSt+yOnd1JnJXsCUlIKFAIKVB22tJCMm2B0g5TCjN0usyvHaYBAl1n2kILLaXQ0lIgZStJKYQ9LA0EspOF7IudfU+cxLakd39/SE/R8mTLtvzelXS/n48+lt+7lq7Pe+/cc88591xRSmEwGAyG3MfndQcMBoPBkB2MQjcYDIY8wSh0g8FgyBOMQjcYDIY8wSh0g8FgyBMCXn1xZWWlqq6u9urrc4IlS5bsV0r1aM3fGvk2jZFt+9Ja+RrZNk9TsvVMoVdXV7N48WKvvj4nEJFtrf1bI9+mMbJtX1orXyPb5mlKtsblYjAYDHmCUegGg8GQJxiFbjAYDHmCZz50gGAwSG1tLfX19V52w3NKS0upqqqiqKjI664YDIYcxlOFXltbS+fOnamurkZEvOyKZyilOHDgALW1tQwaNMjr7hgMhhzGU5dLfX09FRUVBavMAUSEioqKgp+lGAyGtuO5D72QlbmNkYHBYMgGnit0g8FgMGSHglfoIsKXv/zl2O+hUIgePXpw+eWXx4699NJLTJ48mZEjRzJixAi+853vAHDnnXciImzcuDHW9pe//CUiElscUVdXx7/9278xZMgQRo8ezbnnnsuiRYtc+u8MBkMhkbFCFxG/iCwTkRcczpWIyF9FZKOILBKR6mx2sj3p2LEjq1at4uTJkwC8+uqr9OvXL3Z+1apV3HrrrfzlL39h7dq1rFq1isGDB8fOjx07ljlz5sR+f+aZZxg1alTs96997Wt0796dDRs2sHr1ah599FH2798PRGIIU6ZMYfz48VxxxRXccccdKf0TkRtEZJ+ILI++vpZ1IRgMhrygJRb6bcDaNOe+ChxSSg0Ffgn8tK0dc5NLLrmEf/zjHwA8+eSTzJw5M3buZz/7GbfffjsjRowAIBAI8I1vfCN2/jOf+Qxz584FYPPmzXTt2pUePSJlFjZt2sSiRYu455578Pkioh48eDCXXXYZACUlJbzxxhusWLGC5557jvnz5/P+++87dfGvSqkJ0dfD2f7/DQZDfpBR2qKIVAGXAT8EvuXQ5Crgzuj7Z4D7RURUC/a3u+vvq1mz82imzTNiVN8u3HHF6GbbzZgxg9mzZ3P55ZezcuVKbrzxRt555x0gYqF/+9vfTvu3Xbp0oX///qxatYq5c+dy7bXX8sc//hGA1atXM2HCBPx+v+PfigidOnUCIq6eYDBoAqQGg6HVZGqh3wv8F2ClOd8PqAFQSoWAI0BFciMRuUlEFovI4n379rWiu+3DuHHj2Lp1K08++SSXXnppi/9+xowZzJkzh+eff57PfvazLfrbcDjMhAkTOOecc7jwwguZOnWqU7OrRWSliDwjIv3TfZau8jUYDO7QrIUuIpcDe5VSS0RkWrpmDsdSrHOl1EPAQwCTJ09OOJ+JJd2eXHnllXznO99hwYIFHDhwIHZ89OjRLFmyhPHjx6f92yuuuILvfve7TJ48mS5duiT87YoVK7AsK+ZyScbv97N8+XIWLVrE9773PVatWsWYMWPim/wdeFIp1SAiXwf+BEx3+qym5GswGPKfTCz0s4ErRWQrMAeYLiJ/SWpTC/QHEJEA0BU4mMV+tjs33ngjs2bNYuzYsQnHv/vd7/KjH/2I9evXA2BZFr/4xS8S2pSVlfHTn/6U22+/PeH4kCFDmDx5MnfccQe292nDhg0xn3s8Xbp0Ydq0acyfPz/huFLqgFKqIfrr74FJbfg3C5JwOMzpp5+ekLlkk8sBfYMhmWYVulLq+0qpKqVUNTADeEMp9aWkZvOA66PvPx9to7eFeHx/wq9VVVXcdtttKc3GjRvHvffey8yZMxk5ciRjxoxh165dKe1mzJjBxIkTU44//PDD7N69m6FDhzJ27Fj+9V//lb59+wKwb98+Dh8+DEQyXl577bVY8NVGRPrE/Xol6QPThjTcd999jBw5Mt3pnA7oe0V8htbo0aNNhpYuKKUyfgHTgBei72cDV0bflwJPAxuBD4DBzX3WpEmT1Jo1a5Rn1O3z7rujrFixQk2YMEGNHTtWDR06VN11111KKaV+8IMfqLlz5ypgMfBjYDWwAngTGKEyuFaTJk3y5p/SjJqaGjV9+nT1+uuvq8suuyx2HFgc+cHLwCei7wPAfkCUkW2TWJaljh07ppRSqrGxUU2ZMkW99957sfPRe/cG4H7VAh1jZNs89r3r9GpRcS6l1AJgQfT9rLjj9cAXWvJZhoj1v2zZMgDWrl0bsyJnz54da6OU+j7wfS/6lw9885vf5Gc/+xnHjh1L1yQhoC8idkA/YQonIjcBNwEMGDCg/TqcI8RnaAWDQZOhpQmFuVJUc2+QITu88MIL9OzZk0mTmgw7ZBzQV0pNVkpNttcZFDp2hlbPnj3blKFlsrOyR2Eq9CO1ULfX614Y2pl//vOfzJs3j+rqambMmMEbb7zBl76UHP7J/YC+V9gZWrW1tXzwwQesWrUqucnfgWql1DjgNSIZWimYwTJ7FKZCP7Efwg3NtzPkND/+8Y+pra1l69atzJkzh+nTp/OXvyQnaOVgQF8zunXrZjK0NKEwFbqNeW4LklmzZkHEEgd4BKgQkY1EVkF/z6t+5RLxGVonT540GVqa4OmORZ6jLBDnZfmG/GLatGlMmzYNiASd77777iNgAvqtZdeuXVx//fWEw2Esy+Kaa67h8ssvZ9asWUyePNlu9h8iciUQIuLGusGr/hYKRqHj54c//CFPPPEEfr8fn89Hnz59mDBhAj/+8Y9jTZcvX87MmTNZu3Yt1dXVLF68mMrKSu/6bjB4SHyGVjwmQ8tbCl6hv/fee7zwwgssXbqUkpIS9u/fz+rVq/nKV76SoNDnzJnDv/zLv3jYWYPBYGiaAvehW+zatYvKykpKSkoAqKys5LzzzqNbt24JG1E89dRTzJgxw6ueGgwGQ7PoY6G/9D3Y/VF2P7P3WLjkJ+nPK8VFF13E7NmzGT58OBdccAHXXnst5513HjNnzmTOnDlMnTqV999/n4qKCoYNG5bd/hUi4RD49bntDIYWoRQs/DWMuAwqhnjdmxQK20JH0alTJ5YsWcJDDz1Ejx49uPbaa3n00UeZMWMGzzzzDJZlMWfOnIRNLwytZM08uLsCDm31uicGQ+s4vh9e/QE8/nmve+KIPqZSU5Z0exFNW/T7/bEsiLFjx/KnP/2JG264gerqat566y2effZZ3nvvPff7l28si+aA71kD5dWedsVgaBXHdkZ+Bk962480FLyFvm7dOjZs2BA7snz5cgYOHAjAzJkz+c///E+GDBlCVVWVV53MHyR6u6l0+6QYDJpzIrqIuKiDt/1IQ2ErdKWoq6vj+uuvZ9SoUYwbN441a9Zw5513AvCFL3yB1atXm2BotrAVuhX0th8GQ2uxQok/NUMfl4snKCZNmsTChQsdz/bo0YNgMFX5bN26tZ37lafY1fga6rzth8HQWsKNkZ9W2Nt+pKHgLXSDB5g6OoZcJRw18ERP1alnr1zDKHRX8UUnhGHjcjHkKLarRdPS780qdBEpFZEPRGSFiKwWkbsc2rR6qylPC9tpYqEXTHE/X7RuTshY6IYcxXa5aGqhZ+JDbwCmK6XqRKQIeFdEXlJKvZ/U7q9KqVtb8uWlpaUcOHCAiooKj3Y78V6RKqU4cOAApaWlXnel/bEfAmOhG3KV2L2rp4nerEKP1oa2o1hF0VdWNGFVVRW1tbW4vkvJ4ejmFvtCUOz9RhelpaUFkhYZfQiMD92Qq9gZWpput5dRlouI+IElwFDgAaXUIodmV4vIucB64D+VUjUOn5OwL2NRURGDBg1qdedbzZ1nRn5e/ksYf6P731+oqGhmgD1tNbSO+iNQuxiGfsrrnhQe+RAUVUqFlVITgCpgioiMSWqSm1tNhb3PJa2vr2fKlCmMHz+e0aNHc8cdd6S0EZESEfmriGwUkUUiUu16R7OBHVAyLpe28dL/g798Dg5u8bonhYfmLpcWDTNKqcPAAuDipOO5udWU8j6XtKSkhDfeeIMVK1awfPly5s+fz/vvJ4cn+CpwSCk1FPgl8FPXO5oN7Nxdl4KiGQ6WrQ7oe0bt4sjPhmPe9qMQ0dzlkkmWSw8R6RZ9XwZcAHyc1CY3t5rSYAm6iNCpUycAgsEgwWDQKUB8FadmPc8AnxJvoshtw3LX5ZLhYAmRgP6E6OthVzrXJqIhrEazQMt18sBC7wO8KSIrgQ+BV5VSL4jI7Oj2UhDZamq1iKwA/gOdt5qKTxHUZLVXOBxmwoQJ9OzZkwsvvJCpU6cmN+kH1AAopULAEaAiuZGI3CQii0VkseuB5kxw2Yee4WCZe/gjtftpPO5tPwqRmEL3PkPOiWYVulJqpVLqdKXUOKXUGKXU7OjxWUqpedH331dKjVZKjVdKna+U+rjpT/WQeCWugcsFItUely9fTm1tLR988AGrVq1KbuKkhVLuKO1iFMm4bKFDRoMlRAL6K0XkGRHp79RAq8HSZ6d/muCy69guF01ruegZqm1P4i+EJha6Tbdu3Zg2bRrz589PPlUL9AcQkQCRHesPuty9tuNBUDSDwTL3AvrivUIvqGB+PGGj0PUi3irXwIe+b98+Dh8+DMDJkyd57bXXGDFiRHKzecD10fefB95Qubi81Ja3BytF0w2WORnQtxV6yDuFXlDB/HhiCl0vY9Cm8BS6pZdC37VrF+effz7jxo3jjDPO4MILL+Tyyy9n1qxZELHEAR4BKkRkI/At4Hte9bdNuOxyyWSwzMmAvgYWekEF8+OJVVvU00IvvPK58Ra6BqPsuHHjWLZsWcrx2bNnc/fddx8BUErVA19wu29ZJ+ZycUcR7dq1i+uvv55wOIxlWVxzzTVOg+V/RIP7ISJurBtc6VxbiCl0b1fchsNhJk2axMaNG7nllluaDeaLiB3M3+9yV7NHrB6697rDicJT6FacVa5JULRgiGW5uONDz3Cw/D7wfVc6lC00qYljxycOHz7MZz/7WVatWsWYMQlrDjMK5ievINcaW+aa6o4CdLnoGxTNe2x5azpdzRliPnQ9auK0NZivVcC5OWIuF+/dtU4UnkLXLChaUMQUuhlI24QGFnpBBfPjsY0RTS30AnS5GIXuGcpY6FnB1oke+tCbik9MnjzZbvYI8Fg0mH8QyP3NeTXPcik8ha5ZULSgiFnopjhXm9CgamVT8QmbvAnmx2PLXFMLvfBcLhquFC0YNN8xPWewrUQP89ALFs2zXApcoRuXi6so40PPCpr7cfMak+WiGcbl4h2Wu2mLeYuZ6XhHvJtLw0yXwlPoxuXiHSZtMTsYhe4d8fEfDfVH4Sn0BAtdvxE2rzFZLtlBcz9uXhPWex1L4Sn0eGWi4Qib1xhFlB3CxkL3jHiXi4b6owAVevzSf2Ohu4pJW8wOZmD0jvh7V0P5F55CN0FR77AHUGNZtg3NN1nIa8I57kMXkVIR+UBEVkS3mbvLoU3uFLI3QVHviA/m5fgKcE8xQVHviFfoGsbgMrHQG4DpSqnxwATgYhE5M6lN7hSyN8W5vMOsAcgOYeNy8QwreKqWjoYGYSZ7iiqllL29eFH0lWxe5U4he1OcyztUmFhFVZOL3nqMhe4d4SD4iiLvNRxQM/Khi4hfRJYDe4FXlVKLkprkzq70JijqHVYIfIFT7w2twyh07wgHwR+9h3PRQgdQSoWVUhOAKmCKiIxJapI7u9KboKg32AOp3z2FnpcbGStlgqJeYcs+1y10G6XUYWABcHHSqdzZld6+COLTcoTNW2xZxx6G9ldGebmRcfysUkOFktfY96w9y9Rwhp9JlksPEekWfV8GXAB8nNQsdwrZ2xdF/OaBcJPkh8EFhZ6XGxknZFkYC91VbNn7o0ZJLip0oA/wpoisBD4k4kN/QURmRzfXhVzalT5mKfo9vyA1NTWcf/75jBw5ktGjR3PfffeltBGRaSJyRESWR1+zPOhq27EHT5d96OFwmAkTJtCzZ08uvPDCZjcyRuf4DyRlaRmF7iq2q0tjl0uzG1wopVYCpzscnxX3PncK2ccrFo8VeiAQ4Oc//zkTJ07k2LFjTJo0iQsvvJBRo0YlN31HKXW5F33MGrbycdGHDtnbyFgp9RDwEMDkyZO9m32a0hXekWKh6yf/AlwpGlXiPu9dLn369GHixIkAdO7cmZEjR7Jjxw5P+9RuxOQefRjC7lqXbd3IWBvMOgrviCl02yjRT/6Fp9DjfbkajbBbt25l2bJlTi4BgE9EV+q+JCKj032GNm4BJzxwueTlRsbG5eIdyS4XjfSHTeHtKRrLcvHeh25TV1fH1Vdfzb333kuXLl2STy8FBiql6kTkUuB5YJjT52jjFnDCA5dLUxsZE7HEIdc2MjZBUe9IdrloaKEXnkKPD4pqcEGCwSBXX301X/ziF/nc5z6Xcl4pdTTu/Ysi8hsRqVRK7Xe1o21FJVvo7b9StKmNjO++++4jkGPxH9DGQq+pqeG6665j9+7d+Hw+brrpJm677baENiIyDZgLbIkeek4pNZtcxVboGqctFp5CTwiKeqvQlVJ89atfZeTIkXzrW99ybCMivYE9SiklIlOIuMkOuNnPrGAl56F7P5jmJJYeBklBBfRt7FroGt/DhafQE4Ki3o6w//znP3nssccYO3YsEyZMAOBHP/oR27dvj2/2eeBmEQkBJ4EZWvt40+FBHnpeEu/H9VCGffr0oU+fPkBiQN9BoecPlv5ZLoWn0OMXFnk8ZTrnnHNoSjfffPPNKKXuB+53r1ftRGwgNQq9TcRiEd4q9HgyCegDO4HvKKVWJzcQkZuAmwAGDBjQrn1tE+Fko0Q/hV6AWS76uFwKClvutnVjqi22jng/rgYKPcOA/njg10QC+iloUeMpE2yXi8YWeuEpdI1WihYUKUFR/R6GnCDeIPFYhpkE9O3S20qpF4EiEal0u59Zw0oKimp4DxeeQo9PW9TwguQtxoeeHRJcLt7dv5kG9O26ODkd0LdJDuxraBAWoA/dWOiekLKwyLhcWoUmQdGCCujbpKyl0M8gLDyFblwu3mDL2u9e+dy8RJOgaEEF9G1Syufqp9AL2OXiMwrdTVIsdP0ehpwgOaifwwZvzhGfIQda3sOFp9BVOKLMjUJ3F+NDzw5mYPSOZNkbC10DrKhCR8zD4Cb2zW/7H03aYutIyRYyA6NrWEn3sIb6owAVeshY6F7g0QYXeYcJLntHnmxB119E3hSRtSKyWkRuc2iTO7vqKAtEIi8NL0jeYizL7BAf1ActrcS8JcVtqJ/sM8lyCQHfVkotFZHOwBIReVUptSapXW4U4Yl3uWjoA8tbrKQNLoxCbx0pflxjlLhGPmS5KKV2KaWWRt8fA9YS2YcxNzFBUW8wQdHsEF9cDrS0EvOWHAhIt8iHLiLVRPYXXeRwutlddbTYUSfmQxfPqy0WFMlBURcUel5uwm0rdDEDo+vkgIWe8cIiEekEPAt8M37ThSgZ7aqjxY46VhgwFrrreBAUzcua3VaSD11DpZK3xBS6LXv99EdGFrqIFBFR5o8rpZ5LPp9TRXjsoCgmKOoqsWCee5tE5+Um3Ca47B0pbkP99EcmWS5CZN/FtUqpX6RpkztFeOygqJigqKvEr9AVn+uKKBubcGtBsoWuoR83b8mBhUWZuFzOBr4MfCQiy6PH/hsYAKCUepBcKsJjgqLe4KFCz8Ym3NpswpBsoZt72D3yIW1RKfUuIM20yZ0iPAlBUf0uSN5iKyIRVxV6tjbh1iL+A3Hpn8bl4jpWCIjev6ClhV6AK0Xj8tBRpriRW3hgoedlzW6zsMg7rBD4osYgaCn7AiyfG7dSFCIKXZqcgBiygYpT6D6/Kwq9mZrd9l5nueMuBFNCwUtUOFJpUWMLvfAUuuV0UbyZqNTU1HDdddexe/dufD4fN910E7fdllhZIWo93gdcCpwAbrAXeuUUCRa6Owq9qZrdN9988z7IMXchpFroGiqVvMWKi7+BllkuBajQQ0kWuncXJcM86UuIBOmGAVOB30Z/5hbJLhcX0hbzkpSFRUahu4Ydf7NDihoOpoXnQ1fhxIvi4QORYZ70VcCfVYT3gW4i0sflrradhKCoSRltNSZt0TviEyo0Lb9deAo9edqkSdpXE3nS/YCauN9rSVNLR4vSCumI35xbfFo+DDmBSspyMQOje1ihU7sViU9L2ReeQk8Jinqv0JvJk3aK2Do6hpVSDymlJiulJvfo0cOpiXfEtu9yN20x70ix0L2RY4Z1ckREfiUiG0VkpYhM9KCr2cMKR7JcQFujpAB96NGgqCZ+sObypIlY5P3jfq8CdrrSuWyiknzoGlo3OYEKE8mF9tblUlDxH5uYDx1t91MoPAs9FhS1XS7eZahlkicNzAOui1o7ZwJHlFK73OtllrAzAmJ56EahtwornDjD9EiOBRX/sUlQ6Hrew4VnoceW/nsfFG0mT9rmRSIpixuJpC1+xfWOZoOYRS6m7EJbUMkxIO+VSiviPwkGiTZlFZojB3zohafQNQqKNpUnDXDzzTcTXeRyi3u9aifiMwSMD731pORCe6tUshH/0aasQnPEVpmjbemQwnO5JKctGkvRHWxXAWj7MOQEsaC+rdC9GxgLJv5jk+xy0dBCLzyFbiW5XDS8KHlJ/Ipc8Ru5txZlaeFyKaj4j028Qsf40PXACoO/SAuXS0FhWdpPV3MCy85y8dblUlDxH5uEtEU9s1wKT6FrtFK0oFD65/DmBMlBUY/kWFDxH5scyHIpTJcLegRFCworxCmXi57+x5wgJahv5OgasTUsaHsPF6ZCj69prHGl1LwiIShqslxajQprExQtOJIXFuWihS4i/UXkTRFZKyKrReQ2hza5s8TXfiA0WSlaMMRcXWg7Xc0JLEsLl0tBkgNZLpn40EPAt5VSS0WkM7BERF5VSq2Ja5M7S3yNy8UbkoOiGj4MOYEmPvSCxN6xCLQ1Spq10JVSu+wNFZRSx4C1pFb7y50lvhqtFC0oEqwbvyubA+RlAalYlouZYbpOvA+dPMhyEZFq4HRgUdKp3Fnim5KHrt9FyUsSXC7iiu83LwtIxQwSb6stFiRWKJLyDNqWr8g4KCoinYBngW/G75Run3b4E8clvp6Xd40t3zUuF1dJDoq6YFnmZQGpZIPEzDDdIx+CogAiUkREmT+ulHrOoUnuLPGNZQmYKaureBwUbesGItpsHhJb+m/2FHWdHAiKZpLlIsAjwFql1C/SNMudJb4pLheTtugKlncKPVsFpDyfXULqBi0aWol5S55scHE28GXgIxFZHj3238AAAKXUg+TSEl+zUtQbkivVuWTd5F0BKc2qLRYUCYvjcjQoqpR6F2crJr5N7izxjZVxNT50V4l3ubhU2KgFBaRuFZE5RIKh+s4u4VSRs9hGxSYo6ho5kLZYWLVclIqrVmeyXFwlPijqc8f/2EwBKdtvkjuzS/AkuGyIovRf+l9gCt3eBs2sFHWdFAu9/S3LpgpI3Xzzzfsgx2aX4HlwuaAxG1xohn0BjMvFfVKCokburSJhxa1R6K6SXA9dQ2OwsBR6ws7zJijqKvEK3SWXS16ijBw9w5TP1Qx7mp9goZu0RVfwwOWSl8T70I0c3cUKgc/2oeuZ5VJgCj3OQtdkT9Ebb7yRnj17MmbMGMfzIjJNRI6IyPLoa5bLXcwOZk/R7KDMzk+eYFlJsjcWuvfEB0U1WSl6ww03MH/+/OaavaOUmhB9zXajX1knPkPAZ/YUbTUqWpwLMHuzuki8u9b+qaHsC0uhaxgUPffcc+nevbunfXAFKxTnKohOV427q+XkQKZFXhLvrgVtZV9gCj3eh55TQdFPiMgKEXlJREZ73ZlWYVkkbEEHuSJ7vbCX/kNkpuORDJtzFQKd88JVaBPTHXrnoReWQk+YNulhoWfAUmCgUmo88Gvg+XQNtSkg5URy/rR9zNAy4i10D1eKFoyr0CbFQtcz9bawFHqCy0WPoGhzKKWOKqXqou9fBIpEpDJNWz0KSDmRvMLRPmZoGckDo0eDYsG4Cm2sZB+6nrtuFZZCj1nokjMKXUR6RyteIiJTiFyzA972qhVYocTpqn3M0DKSSyjoLcOMXIVazyxtbDnbtVxcqkfUUgpr6b89RRI/urhcZs6cyYIFC9i/fz9VVVXcddddBIPB+CafB24WkRBwEpih0q1n1xkVb6HrkWGUk6Tk82trkBwHJiml6kTkUiKuwmFODZVSDwEPAUyePFnPezvFQtfTh15gCt32g+lTT/rJJ59Me+7mm29GKXU/cL97PWonEpas29un6fdAaI+yiKUt6r1S1Ip3FYrIb0SkUim13+uOtQpHH7p+si9Mlwu540PPG5L3FAUtHwjtiR8Y9V4pGsgLV6FNSpaLnj70ArPQ9ctDLxiSi3OBlg+E9iTXcvFoUGzKVfj1r38doBxYlfOuQhsnl4uG7q5mFbqI/AG4HNirlEpJOhWRacBcYEv00HPapiglpC0aP66rJBc2AmOht4bkbCGPLPSmXIVR9imlJrvRF1dwcrloqDsysdAfJeLD/XMTbd5RSl2elR61JzmYtpg3mKBodoivJ2KvuDW0P8lZLrm6UlQp9TZw0IW+tD9WfNqisRJdJbmON7S77POy8Fl8LRcPV4oWHMk+dPT0oWcrKJob+aZO9dBz2K2XUziuFG1f6zIvVzNaydUWtQ2K5hfxxiBEXS76zY6yodAzXpru+UpGDcvnFgwJPnR3slzycjVjguvKVFt0DcdaLvoVmGuzQm/J0nTPcbTQzQPhCglZLnYeuhbWZW7MLm1Sqi1qIcP8x6naImhnELZZoefU0nTHHYv0uiB5i55B0dyZXcY6kpyHbu5fV4gFRZPLV3h+DyeQSdrik8A0oFJEaoE7gCIApdSD5NLS9NjSfyHmctHsguQlMbnrlbaolDoa9z43VjPGD4w+HwQbve1PoeCUhw46GCUJNKvQlVIzm+2yic4AACAASURBVDmfO0vTE1wuxkJ3jeTdXjSpoyMivYE9Siml/ezSJnmBlmYKJW9xykMHz42SZAp4paiePrC8JJ3/sZ0fhrwrfKYUoIjNLjWtJ5KXpPWh6yX/AlPo8RfFKHTX8Gi6mneFzxyXn+ulUPKWlDx0PS30wizOZVwu7qIccnhBu4dBe5RDLEIzCzFvSVkpqqf+KCyFHh8UNRX/3CNHAkra4zQwmrRFd3DasSj+uCYUlkJPDs5putor70in0I0yahlOSkUzhZK3pAuKamaUFJZCT34gNK3HkHekDKS2dWMG0xaRYqGblaKu4bRSFLQbUAtMoTtEqo2F3v7Ebnq9rRvtSZGjsdBdI0eyXApLoRuXizeYoGh2sDMq40soGBm6gy3n5KCoZvIvLIXu5HIx0/72J93DoJl1oz0pA6Op5eIaaX3oeumPwlLoKn7pP1q4XDKo2S0i8isR2SgiK0VkostdbDtOroKE44aMcCzhamToCk57ioJ293BhKXQNt5HKoGb3JcCw6Osm4Ldu9CurOLm6wPPBNOdwkqNmCiVvsaIrjFPKV+gl/wJT6E5ZLt4qlQxqdl8F/FlFeB/oJiJ93OldlkiRu57+R+1xTP9UnrgNm5tZAuT8zDKetKm3et3DhaXQndLn9LcS+wE1cb/XRo/lDslyt33pxv/bMlJcht5ZiRnMLLuS6zPLeEweuoY4pc9pNsI6IA7HHAtIabMJQzI58jBoT8z48H6BVgYzy27k+swyHnvHreSa/polVRSWQrcfCF9OWei1QP+436uAnU4NtdmEIRkrybI0tehbh1NQNP64XhSR4cxSW0MkHisUV5gLbY2SwlLoMUsmTrHor9DnAddFs13OBI4opXZ53akWkTYoqtfDoD3p5Jg7rivHmaW2hkg88XXoQdvBtMDK54YBQae0xQxqdr8IXApsBE4AX/Gin20ibUBJ+8FUL9IWOdNSjkEynFnmBFbo1MwetF0pmskWdH8ALgf2KqVSQtrR/UTvI6J0TgA3KKWWZrujWUGFky6K9ytFM6jZrYBb3OtRO5Ajy6a1J11QVDMrMcphIjPLOcBUcnFmGY/tQ4+hp9swE5fLo8DFTZzPnTzpZD8YphaGK6R1FbSv7PNu0ZZyCOqDJy6XmTNn8olPfIJ169ZRVVXFI488woMPPsiDDz5oNzkCbCYys/w98A3XO5lN0vrQ9ZodZbKn6NsiUt1Ek1ieNPC+iHQTkT5ajsZOfjDNLkheki6Y184W+g033MCtt97Kddddl65JvDEylYgxMrVdO9UWkoPLHsYimppZ2iilcntmGU+yha5pHCgbQdGM86Q9j2anXBRTPtcV0pYebV/LMu8WbaUtQ5wzQdHcJcUYzN+0xYzzpD2PZjsqdL0uSF4S274rWaF7LvvcMUbAISjqTzxuaD8KyELPOE/ac1ICG75TJUkN7UfYroORVNjI+4chd4wRcC4uB0ahu0FKloueAelsKPTcyZO2rFQLXbMLkpek22DXe9nnjjECcQNg3I5FCccN7UaOLCzKJG3xSWAaUCkitcAdRFaBoZR6kFzKk7ZCp6b9YIKibpFsoYMWlS6JGCO35kxqXY5sVJyXpPWh6yX7TLJcZjZzPnfypJ1ySb1XKvmPXXo0YTBt/9lR3i3aSgmK2j50ExRtd9L60PUyCAtspagJinpC8sIiIBK/aF+FnneLtlLSFrWJReQ/TvE30M5CL7xaLiYP3X3CSVkuEK10aWTfImJBUe1iEflPirtWz8G0sBS6MkFRT0jZ7YVIgNS4ClqG2WzbO3KkOFdhKXRHH7pJW2x3nIKiLrhc8o60Rc7MwNjuFFAeeu7g6HLR64LkJckLi8DMjlpDulou5h5uf5zib6DdPVx4Cj1ZqRgfevuTvPQfzGDaGmzloV8+f/4TDuZEca4CU+gOuaSaXZC8JBwkYfsuMEHR1pCy4tYodNcIB1OD+qCd7AtMoTv40DW7IHmJFUySO6YwWmtIVxPHyLH9CTeaLBftMOVzvcEKJwVEIUc26NaL5GwhnwmKuoYVBF/csh1joWuAKZ/rDeEgIYRFWw6eOmZ86C0nnByL0FOp5CVRH/qiLQcj97Gms6MCWykaTnW5mLTF9scKoogoIVupTzUWesuJlVAwFrrrhBsTLXRNB9MCU+gh8Bef+l3M4hZXCAdR4mPhwc48vqMnQUt4s4OfcivkWL/WkIaUfH49My3ykpSgqO1D10v2Behyic+0MEFRV7DCNKgAv9rSl26BENUd6tldX8T2/XVe9yy3SE7/9OlpJeYl0aCoUnAs5ENpaqEXoELXq3zu/PnzOe200xg6dCg/+clPUs6LyA0isk9ElkdfX/Ogm20iHGrkUKiYqtIG7jhtO98bWkuHItiy7yh7j9Z73b3cIW1deW9mmYVw78YIN3IoWMR311TztRXDueT1CpT4I4peIwpLoatw4q4jHpfPDYfD3HLLLbz00kusWbOGJ598kjVr1jg1/atSakL09bDb/Wwr2/Ye5IQq5ob+eyn2KXwCXQMhRIX55WsbvO5e7hAOAhIXB/Iuda5Q7l0gYoUri3k7OnI4FOCavvvYdRxOWEUcq9NrllkwCv2JRds5Xt/A3rrgqWwLjxcWffDBBwwdOpTBgwdTXFzMjBkzmDt3rmf9aQ/ClmL3gSP4/H7GdDkROx7wQbdSH08vruGRd7Z42MMcwkry4/q8q4deCPdujKgVvqu+mK8O2MPVfQ7wg+HbaSTAu+t2EArr40fPSKGLyMUisk5ENorI9xzO58TUSqwwm0904I51A7hzeScaLG9dLjt27KB//1M7oFVVVbFjxw6npleLyEoReUZE+js1AE02Mk7izY/3QvAkHYuSwp/io7zUR8hSLNl+yJvO5RrhZJehdzvPZ/ve1Zn6+ohbsEcH4cxuxwDoX9ZIwO/n6LHj/O7tzV52L4FmFbqI+IEHgEuAUcBMERnl0FTrqZVSimCwkQ0nOtBoCY9tKuOVXaXUnfDOh6scUiZFUvI+/g5UK6XGAa8Bf2ri87zfyDiJx97fRqdAmA5JCt0SP6V+iwHdO/DhloOOssgG+eDnfWLRdp5YtL2JjYrdt9Czee/qaIjE88wHmwAY2fFEQk5Fkd/HoPIA9762no179XC9ZGKhTwE2KqU2K6UagTnAVe3breyzdPthxApS3THIj0du43+GbacuXESwscGzKVNVVRU1NTWx32tra+nbt29CG6XUAaVUQ/TX3wOT3Oth26g5eIK3N+yjdweVlMMLlgTwW41MGljOvroGltcczvr354Ofd/eRep5ZUsNvF2zk7XU7OWkFWH4wQMjC08Ut2bx3dTREbI7WB/nzu5E4T++yRDkrn5+KUvD7hFlzV7WbUdISMlHo/YCauN9ro8eSaXZq5dVIXNcQYv6qXZRImMGdIrm8ozqf5LTOjfhVkN8u2ORaX+I544wz2LBhA1u2bKGxsZE5c+Zw5ZVXJrQRkT5xv14JrHW1k23g6cU1oKBINWJJokJXEsAfbmBsv64U+YWnl9Rm/ftz3c/7h3e3cOmv3uGjHUcoCvg4cOQ4R0IBPvNGd856oRuv7CqLNPQg0yLf712bh97azImoy0X5EstXWBKglCAXjerNwk0HmLdipxddTCAThe609iN5KMpoauXVSPzHd7dwvDFEMYklMHuXhSmWMA8s2MiOwydd649NIBDg/vvv59Of/jQjR47kmmuuYfTo0cyaNYt58+bZzf5DRFaLyArgP4AbXO9oKwhbimeW1DK0ZydHhW75IhZ6aZGfMX278vcVO6kPZtfSzKaf1wtj5LW1ezh0vJEbzx7E184ZzOl9SuhcLPz7oJ108FvctKgy0jDkvkLP53vXZu/Reh55dwufHtEdIJKmGIclfnxWI1MGdWdsv67c84+1HK0POn2Ua2SyUrQWiL/Jq4CEoUgpdSDu198DP21717LD4RONPPTOZsb07oQcVlhxo6wSf0TJAz9/ZR2/uGaC6/279NJLufTSSxOOzZ49O/ZeKfV94Psud6vNvLtxPzuP1DNzygD8Hzdg+YoSzlsSwGdFFNHEgeUsqznMy6t3c9UEp8lf62iBn/dJpVSDiHydiDEy3eGzHgIeApg8eXK7z60P1DXw4daDTBxQzsCKjgAUhU6Av5hzuh9lSrdj/GXPQBoOBNi99xAD27tDDuTrvWvzqzc20BAKM6S8GLZEFHg8SgL4rQZ8InxyWCW/XbCJX7++gdsvcwoxukMmFvqHwDARGSQixcAMYF58A52nVn/851aO1Ye46LRyIHGUtSSAD4svTali7vKdbD9wIt3HGFrIU4trKO9QxMjenfFbDc4ul6hCH1TZkaryMv76YY3TR7WaXI5RPPb+NoJhxSeHVcaOBcLHCfsipSuKfYov9tpGUIp49+NaDp/Qa4FLrrN1/3HmfFDDGdXdKS+NHFNJJaAtXwBftL5OVXkHJvTvxmPvb2PvMe8SLZpV6EqpEHAr8DIRRf2UUmq1iMwWEdtppuXUqq4hxKMLt3LRqF706Rz5V+MVi63ce3cUBPjNgo1edDPvOHi8kVdX7+Ezp/cj4PfhDzek+h99fvzhiB71iTC6b1cWbjrAqh1HstaPXPXzKhVxVw3r2YmeXUpjx4tCpxQ6QIlPEfD5UKEGvvqnxV50NW/51RsbKPL7mD6iJ/6o0k52uSjxx4wSgOkjehIMKx56y7s0xozy0JVSLyqlhiulhiilfhg9NkspNS/6/vtKqdFKqfFKqfOVUh+3Z6cz4YlF2/nu0ys4cjLIN84f6nhR7Pddi2FydTlPLa7hN28apd5Wnli0jcawRZfSIlAWfhXEkkSXi4pzuQBMqe5OScCX1Zzepvy8QNdoM+2MkaXbD1N76CTj+3dLOB4IncDylSQcE7+fAV38LNl2iPc2HcDQdn79+gaeX7aDSQPL6VxahKhIWmiqDz2Az2qI/V7RqYTRfbvw2PvbON7gTTmGvF0pGgxbvLthP+cMrWRC/26xqVGCyyVqNfqsIOcOiwRp/7lxv/udzSPqg2EeXbiN4b060atLKf7oDZ8SFJUAfhWMLewqK/YzZVB3Xlixk/vfyN6geumll7J+/Xo2bdrE7bffDsT8vEdAT2Pk7yt2UhLwMapPl4TjReFECx0iA2O/zn7KOxRx+/Mf0RDSq1hULrJg/b6YXxyI6Y5kH7rlKyIQbkg4dtaQShpCFs8tzX7WVibkrUJfuv0QxxpCnNa7M08s2h6bGlkOFrrfCtKtQzHjq7rxwdaDHDpu/JGtZe7yHeyva+CT0QHSH02ps5Ly0MPRIGkgfCq76Owhlfh8woJ1e13qrX6EwhYvrNzJ9BE9KS1KVCCB0HHC/kSFbkmAItXIZyb0Y/O+4/zmTW9ScPOFmoMnWLb9EGcM6k7n0sg9ms7lEvYVEwgfTzjWv7yMqvIyHl24FctyPy89LxV6KGzx9vp99C8vY3BlJEPglIWe6kO3z31yeA+CYcWf39vmco/zA8tS/P6dLYzu2yUmd78VCRAlW+hhfySHujh4NHasS1kRZ1R3Z+n2Q9QcLMwA9XubD7C/rpHyDomKG2VRFD5JOMnlYkX9uMN6deaqCX357YJNbNx7zMUe5xe/fWsTIhKbsQMx12CKy8VXTCCUeJ+KCJ8YXMGmfcd514PZfl4q9H98tItDJ4KcN7xnLE3N52Ch20rGpyIKvXeXUkb07syjC7dwotFsfNFSFqzfy8a9dYzu2zUmd9sCT05bDPkjwb54hQ5w3vAe+ESy6nbJJeYtj7hbTuvdOeG4rThSXC6+QMyt9YPLR9GpNMCtTyzLek5/IbDz8EmeXlzD5IHldC07db8GwhGjJJwUBwr7ivGrEL6khV1jq7rSo3MJf/in+0Xn8k6hK6X43Vub6dGphBF9Tj0UTj70Uxb6qQty3vAeHDoRzHoKXSHw0Nub6VpWxNh+XWPHAqHIlDTZsgynUehdo1b6s0trCy6NtD4Y5qVVuxndtytF/sRHsyg6tU8OikYWt0Tu7VdW7+HK8X35ePcx7vr7anc6nUf87q2Iu+q84YmLHu1ZpvKlKnQ4dW1sAj4fXz5zIAvW7WP9HndnS3mn0N/ZsJ81u47yyWGV+OIWkdg+dJWwsMgXPXdqddfAio5UV3Tk/15eZ3J7W8CKmsO8v/kgZw+pwO87JfeiqGVpJfl+Qz7b5ZKapnje8B74fMIDBZZx9NraPdQ1hJiQlN0C6S10S05Z6ADDe3XmvOE9ePKDGp5f5rgq1uDAnqP1PL5oO+OrutEtyd1lp9emWuiRwTXZ7QLwpTMHUlrk4+F33E1hzDuF/uBbm+jVpSTloShysBRjLhcrcbnuFeP7cDIY5s55q7UouJML3Pvaerp1KGJydfeE405yh1MulxIHhd6lrIhJA8t5eklNQVnpzy/bQa8uJQzu0THlXEyOyUFRXxGBUGLZigtG9mJgRQf++28faVMFUHcefGsTllJMO61nyrlTbsPk8hXF0fPHU/5m/qrdjK/qxrNLd7i6K1deKfR3Nuxj4aYD/OsnBxNImrLaQo+3cJxcLgB9upZx/oiePL98Jz+Z/7En0epcYsm2Q7y5bh9Tq7unZmY4yB0gFIgorZJG51ro5w2L+tLfLIwdjQ4eb2TBun1cNaFfwszSpigUUcxhX2nC8bCvNHbOxu8TZpwxAKXg3x5bzOPvmyB/U+w9Vs8Ti7Zzev9yuncsTjlvz4BUSmA/6nJxsNABzhlaiVKK37hY/C9vFHp9MMw9L6ylvENRiv8RnC3FmEJXqQV1pp/Wky9OHcDv3trMFfe/yx89CHDkAkop/u/ldXQs9nPmkIqU8+ksdMtXRNDfkdIG58UwXcqKmDKoO88u3cG2A6kWUL7x9OIaQpaiJOD8SJ5S6MmxiBKKQql+2q5lRUwf0ZNN+46zwVjpTfKtv66gMWQx7TTngoH+cD1hCSRuME+cDz3kfH9WdCph4oBynli03bXif3mh0JVS/PdzH7FuzzGuGNc3Y4VuxfLQU33lIsI9nxnD/1w2kjU7j/K7tzaz64j7FRl15/89u5L3Nh9g+oielAT8KedPyT3V8gn7SyhrTJ/ade7wHgR8+Z/xEgxbPLpwK4MrO9Kna5ljm5hC9ycpdF8JReGTiMMmF1MHd6e8QxGvrN5tXIdp2HusnkVbDjChfzcqOpU4tgmET6ZkaUGcDz2c3i04fURPwkpx25PLstPhZsgLhf7wO1t4btkOPjWiJyOSVtfZFAePYuFLuDDJeejJPPlBDR2KA1x/VjWHTjTyb48tMelgcew4fJK/r9hF//Iypg5Otc4BikLHUEjM3xhPMJDeQgfoUlrE5IHlPLu0Nqs1XnTj2SW17DpSzzlDK9O2STfTsRW8k5UY8Pn41Ihe7DxSz8urd2exx/nD/85fh2VFFG86ikInHO/fmA89jYUO0K1DMVOiays272v/mVLOK/Q75q7mRy+uZXTfLpzfxEUpbTwQ8dvGTZusZhS6zfBenfnCpP6srD3CjY9+mJ2O5zgNoTDfeHwpllJ8YVJ/R78vQGnDAYJJcrcJBjo1qdABpo/oRceSAN95OjItzjeOnAjys5fXMXlgeUrueTzFUbdKsoUeivrUndwuAOP7d6OyUwm/eHU9YRMLSuCDLQd5ekktZw2tSGudQ8QYDPlTZ0724JrO5WIz7bQeBPw+7vlH+9d9y2mFvnFvHXM+3E7vrqVNKhWA0oaDEcUSx6kodfOulFF9u3Dm4AoWbjrAwk2FXe9FKcXtf1vFiprDXD2xisrO6R+GsoYDBAOdHM8F/R0pbcLlApEaL5+Z0I+Pdx/jvtfXt6nfuqGU4r+ejRSQu+uq0U612mMEQscJSyB1+XnMQne2/vw+4VMje7J+Tx0vrPR+Rx1dOHi8kW/OWcbAig5Md8hsiac4eCSWlRVPOKY/ms7E6lxaxKdG9OSNj/fy+to9re90BuSsQj9yIshNf15MwCd86cyBFKcJJtmUNh5IUej24paiYGbJ/xeP7k1Fx2K++/RKjnm8M4mX/N8r63hmSS3TR/RkTNwiIifKGvamV+iBjpQEjzY7QxrZpwuTB5bzwJubmL9qV6v7rRNKKX46fx0vr97DxaN7s6KmaZdSUaguxTqHU1ZicRP38Nh+XenVpYS7X1jj2f65OnHkZJAvP7KI/ccb+dWM0ykpSo39xFPaeICQPzWVVIkfC1/aLJd4zhpSSY/OJXz76RXtuqtRTir0UNji3+cso+bQCb44dWBq3QsHShoPEky6KJYEUEjKSq90FAd8fGFyf3YdOcnsvztuNpzXWJbiyw8v4oE3N3FGdTmfasLFBYBSdDpRQ31xuePpxqiiL2nG7QJwxfi+TOjfjW89tYKPdx9ttr3OKKX45WsbePCtTUwZ1J2zHLKDkikO1aX4zyGuhEIovUx8Ilw0qjf76xr51euFkQaajv11DXz5kUWs33OMmWcMYPXOpu8lsUJ0PLHD+R4WwfIXp50dxeP3CVef3o+jJ4N89+kV7Taw5pxCV0pxzz/W8vb6fdx91RiqK1NHToc/OuXLjUckmvaVebBiQPcO3DxtCE8vqXV9FZiXNITCfOfpFbyzcT9nDu7OVRP6NekigMggWhyqo764u+P5UFShlzU2r9CL/D4uHt0bvwhf+eOHOVu86/CJRr791Ap+9foGJg0s58rxfZuVI0Ty9UP+DinHg4FIEkBZfdMVKkf26cLEAd349ZsbeapAy1ps3lfH536zkLW7jjLzjAFNxixsOtTvwodFQxqjJOTv4Lja2YkBFR25ZEwfXl69h5sfX8q+Yw3N/1ELyWRPUW2wLMXPX13Howu3cvaQCjKN8QTCJwhYDbHFLPGEfSXNBjWS6dO1jNF9u3DPP9ZSc/AEt10w3HFBQr5Qc/AE33h8KR/tOMIFI3tx/mk9MlJCw7Y/BcCJ0l6O5+0BtrnAqE2XsiJuOLuaPy3cyhcfXsRfvjqVARWpSk43GkMWb3y8h5dW7eaNtXs53hjimxcMo7JTSZNxn3hK08QigoGOhH3FdDzZvCvqyvH96FgS4L+eXcnfV+7kMxP6cd5pPahsIiCYL7y7YT///uRSfCJ87ZzB9O+e2X3T+Xhk8EtnlDQGOlLWkHlM7eyhlUyuLudHL67lrJ+8zoT+5cy+ajQj02TntZSMLHQRuVhE1onIRhH5nsP5EhH5a/T8IhGpzkrv4li76yhfefRDHnhzEzPO6M8lY/s0/0dRyhoiu7SnWOhEsgRKGg+3qC8+Ea49oz83nj2IP723jTN++Brn/98C7nttAws37m9Rpcb58+dz2mmnMXToUH7yk5+knHdDtukIW4rHF23j0l+9w4a9x/jS1IFMH9EzI2Ve2rCf0ZsfZkeP8zjWwXkLY9vlUtqCB6JP1zL+/NWpHDkZ5MoH3mXOB9ub3NTBS/nuPlLPjY9+yKS7X+Xrf1nKq2v2MLxXZ245fyg9O5dmrMwB5xlm5B/geFlfutY1P1ssDvi4aFRvLh7dm+U1h/n20yuYfM9rXHn/u/z8lXUs2XawRZkwOt+7NruOnOR/nv+ILz2yiIDPxw1nVWeszAE6n9gOpFfooQwytZL5ytmDePmb5zKuqhvLth/ikvveYdr/vskLK3cSbKMrplkLXUT8wAPAhUAt8KGIzFNKxTuRvwocUkoNFZEZwE+Ba1vambClaAxZHDrRyPaDJ9h24Dgrao+weOtB1u+po1NJgMvH9WFsv64ZKRWAjidquWDRVwDni1JfUkm/fW8zfOvjrK/+YsZ9Dfh8DO3Zids+NYyl2w6xaX8d9762HkXEXzaqTxcmDSxnbL+u9OlaSs8uJXQsCVDk91Hk91Hs9+HD4pZbbuHVV1+lqqqKM844gyuvvJJRoxJ2Dc+qbI81BDlyIsjhk0GOnowEZ3w+IeAT/D7BL8LJYJhVO47w/PKdbNxbx6DKjlw9sapFs5ChNc/gtxrZ33WMY8oiQDDQmbCvmIG7X2Zbn4sJBzJ70NbsPMrXzhnEm+v28r3nPuKef6xlVN8uDOvZiQHdO9CnWxn9upXSq3OxK/JVShGyFA0hi52HT/JR7RFeWrWbN9ftxbIUw3t15szBFQzt2SmhcFmmFAWPUBo8xN400/5dlWcxtOYZAsE6QkXOAWgbv084d3gPzhlWya7D9azbc5T1e+p44M2N/PqNjXQtK+KcYZVMHdSd03p1prxjMX6fYFmR/zFsKfw+oSwg3PyNb/Dyy69QPXBAu967SkW+N2QpgmGLxpBFY9iiIRj52RiysJRCEEQiQc/1e47xzob9vL1+H5ZSfGJIBZ8e1bvZ5Ilkuh9dQ9Bfljaw3xjoRPcTNZGdtySzz35iUWSQuHpiFZeM6c2SbYdYtOUgtz6xjB6dS/js6f24YGQvhvXsRNeyInwtuGcycblMATYqpTYDiMgc4CogXqFfBdwZff8McL+IiMpgedqHWw9y/R8+oCFkOVoHpUU++pd34LKxfZjQvxsdSzL3EvU6sIjzFt9KIFr+8kRJahBvf9cxVBxdzeS1P6FD/W6Wj/h2xp8P0KtLaWy2cLIxHBmIDh5n24ETzPlwO48uTD/iNuxYy2HVlYsfWYfIOg5VTOTOXz/KU7/9WXyzVssW4NyfvcmuIycJhlueg1xVXsaMM/q3aAAFGLr9KcZteIDDnYZQX5I+4Kd8fmp7nMeAPa/xuTemM3fafBqLUysNOlHRqYSrJ1Yxvqoba3YdZfeRej6qPcLJuIVfDTvWctjqykUPf4zft46yIWcxd+7cZKXTavn++MW1PLpwK41hi+TWXUoDnDW4gqmDK9rsjpu66i4AjpVVOZ7f0vcKTtv2BNM/vIlXznoio8/0idCvvIx+5WVMH9GLE40hNu6tY/2eOt5ev49/rGzahRO5d7tx0SPrKAls4FDFRB7/67P88K7syHbp9kN88feLCFlWq+5dgO4dizlzcAWfGFxBeSuuQZ997zCk9m8c6DIygNN0pwAABO9JREFUrVFSV1ZF74MfMmbj71g17OYWf0eH4gCfHNaDs4dW0q9bGY8v2s4f3t3CQ3F76xYHfPz3JSO44exBzX6eNCdbEfk8cLFS6mvR378MTFVK3RrXZlW0TW30903RNvuTPusm4Kbor6cB65r/lzOiEnArOTyb31UOdAHs6kndgU7A9ujvp0XPNSvb6LmWytctuXl1fZqS70ClVA+P71035ZJtmpJtJdAR2EPhyNbN7xuolHIsPJOJues0NCWPApm0QSn1EPBQBt/ZIkRksVJqcrY/t72/S0S+AHw6abCcopT6d/u7AKfiHo6jcEvl65bcvLo+zcnX/hOHj3Hl3nVTLtmmKdlG/69qEXHaZSMvZavLtczE6VML9I/7vQpIXnIWayMiAaArcDAbHcxzjGzbFyPf9sPIVkMyUegfAsNEZJCIFAMzgHlJbeYB10fffx54I1Mfb4FjZNu+GPm2H0a2GtKsy0UpFRKRW4GXAT/wB6XUahGZDSxWSs0DHgEeE5GNREbgGe3ZaQey7sZx47uak230u/5M+8nWLbl5cn1y4N51Uy5ZpRnZvh9tVkiy1eJaNhsUNRgMBkNukHNL/w0Gg8HgjFHoBoPBkCfkpEIXkf8VkY9FZKWI/E1EHFejiMhWEflIRJZHUwBb8h3tXu5ARPqLyJsislZEVovIbQ5tponIkej/sFxEZrX0e+I+q13l5obMop/jqtyyTXNyymXa8sy18XvzQie0GaVUzr2Ai4BA9P1PgZ+mabcVqGzF5/uBTcBgoBhYAYxKavMN4MHo+xnAX1vxPX2AidH3nYH1Dt8zDXhBd7m5JTMv5JbNVyZyyuVXa5+5LHxvXuiEtr5y0kJXSr2ilLIrYL1PJAc2m8TKHSilGgG73EE8VwF/ir5/BviUSAvWxwNKqV1KqaXR98eAtUC/NvW86e9rT7m5IjNwX25ZJhM5GVpIvuiEtpKTCj2JG4GX0pxTwCsiskQiy4szpR8QXzS6llSFEWsTvZGOAM3vVJCG6PTsdGCRw+lPiMgKEXlJREa39juSyLbcXJcZeCK3tpKJnHKZ1j5z2SQvdEJr0LYeuoi8BvR2OHW7UmputM3tQAh4PM3HnK2U2ikiPYFXReRjpdTbmXy9w7FWlTvIBBHpBDwLfFMplbyFylIitRvqRORS4HlgWBOf5ZXcXJUZZFduLpJVGWhIa5+5ZikkndBatFXoSqkLmjovItcDlwOfUlGnlcNn7Iz+3CsifyMybcrk4rVkWXOttGFZs4gUEVFKjyulnnP4H47GvX9RRH4jIpXKoThXtI1XcnNNZpB9ublIJnLKWdrwzGXy2QWhE9qE2077bLyAi4mU7+3RRJuOQOe49wuJVH7L5PMDwGZgEKcCIKOT2txCYgDkqVb8H0JkJei9TbTpzakFYFOIVLMT3eTmlsy8kFuW791m5ZSrr7Y8c1n47rzQCW2Wg9c3QSsv3kYivqrl0ZctxL7Ai9H3g6NCXwGsJjIta8l3XEoke2KT/bfAbODK6PtS4OloXz4ABrfi/ziHyJRsZdz/cinwdeDr0Ta3Rvu/gkiw5yxd5eaGzLyQWzvcvylyyodXW5+5Nn53XuiEtr7M0n+DwWDIE/Ihy8VgMBgMGIVuMBgMeYNR6AaDwZAnGIVuMBgMeYJR6AaDwZAnGIVuMBgMeYJR6AaDwZAn/H+h59KHlG6yDQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 4 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# visualize posterior on parameters\n",
"\n",
"fig, ax = plt.subplots(1,4)\n",
"sns.distplot(samples_mcmc['linear.weight'][:,0,0,0], kde_kws={\"label\": \"MCMC\"}, ax=ax[0])\n",
"sns.distplot(samples_svi['linear.weight'][:, 0,0,0], kde_kws={\"label\": \"SVI\"}, ax=ax[0])\n",
"\n",
"sns.distplot(samples_mcmc['linear.weight'][:,0,0,1], ax=ax[1])\n",
"sns.distplot(samples_svi['linear.weight'][:, 0,0,1], ax=ax[1])\n",
"\n",
"sns.distplot(samples_mcmc['linear.weight'][:,0,1,0], ax=ax[2])\n",
"sns.distplot(samples_svi['linear.weight'][:, 0,1,0], ax=ax[2])\n",
"\n",
"sns.distplot(samples_mcmc['linear.weight'][:,0,1,1], ax=ax[3])\n",
"sns.distplot(samples_svi['linear.weight'][:, 0,1,1], ax=ax[3]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that our variational family has diagonal covariance.\n",
"The SVI approximation of the posterior of the model parameters is over-confident compared to the true posterior from MCMC."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment