{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction to Variational Bayesian Methods\n", "\n", "### Written by D. Schleuter and C. Fonnesbeck" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Bayesian analysis, the most common strategy for computing posterior quantities is through Markov Chain Monte Carlo (MCMC). Despite recent advances in efficient sampling, MCMC methods still remain computationally intensive for more than a few thousand observations. A more scalable alternative to sampling is Variational Inference (VI), which re-frames the problem of computing the posterior distribution as a minimization of the Kullback-Leibler divergence between the true posterior and a member of some approximating family. \n", "\n", "In this lecture, we provide a basic overview of the VI framework as well as practical examples of its implementation using the Automatic Differentiation Variational Inference (ADVI) engine in PyMC3." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are a couple of alternatives to evaluating complex posterior distributions without using an MCMC algorithm:\n", "\n", "1. Divide and Conquer approaches: See [this paper](https://arxiv.org/abs/1311.4780) or [this paper](https://arxiv.org/abs/1508.05880) for more details\n", "\n", "2. One alternative approach is to build an approximation to the posterior $p(\\theta|X)$ using some other distribution $q(\\theta)$\n", "\n", "Let's focus on the latter option, and try to use a simple distribution to approximate a complex one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Distributional Approximations\n", "\n", "We have already seen one example of a distributional approximation: the Laplace (normal) approximation. In that scenario, we constructed a Taylor series of the posterior and threw away the terms of higher than quadratic order. \n", "\n", "Variational inference is another distribitional approximation method. Here, rather than a Taylor series approximation, we choose some class of approximating distributions, and choose the member of that class which is **closest to the posterior**.\n", "\n", "Thus, we shift from a stochastic sampling approach to approximation to a deterministic approximation that places bounds on the density of interest, then uses opimization to choose from that bounded set." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: Approximating a Student's t with a gaussian \n", "Let's approximate a Students-t distribution with $\\nu =3$ with a Gaussian distribution of some mean and variance.\n", "\n", "One naive approach would be to build a set of test points and minimize the MSE between the $\\log p(z)$ and $\\log q(z)$. \n", "\n", "$$\n", "{\\hat \\phi} = \\underset{\\phi}{{\\rm arg\\,min}} \\frac{\\sum_{i} q(\\theta_{i};\\phi)\\left[\\log q(\\theta_{i};\\phi) - \\log p(\\theta_{i})\\right]^{2}}{\\sum_{i} q(\\theta_{i};\\phi)}\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline \n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import numpy as np\n", "import scipy as sp\n", "import copy\n", "import time\n", "import warnings\n", "warnings.filterwarnings(\"ignore\", module=\"mkl_fft\")\n", "warnings.filterwarnings(\"ignore\", module=\"matplotlib\")\n", "warnings.filterwarnings(\"ignore\", module=\"theano\")\n", "warnings.filterwarnings(\"ignore\", module=\"scipy\")\n", "\n", "sns.set(style='ticks')" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "ν = 3\n", "logp = sp.stats.t(ν).logpdf" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEBCAYAAACXArmGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd0XGeZ+PHvqPfemyVb0usuW7bjkhgncXojlJCYELPAZklggbAL57flEFg4ZLMbfmw2xHuSH9WQYLKhpEBIJcWprpL7o2ZblmRZXZasPjO/PzQSsiJbI3tGd8rzOecezbz33rnPtSU9est9X5vT6UQppZQCCLE6AKWUUr5Dk4JSSqlxmhSUUkqN06SglFJqnCYFpZRS4zQpKKWUGqdJQSml1DhNCkoppcZpUlBKKTVOk4JSSqlxmhSUUkqNC7M6gOkYYyKBVcBJwG5xOEop5S9CgWxgp4gMunuSzycFRhPCdquDUEopP7UeeNvdg/0hKZwEePLJJ8nKyrI6FqWU8gvNzc3ceeed4Pod6i5/SAp2gKysLPLy8qyORSml/M2Mmt21o1kppdQ4TQpKKaXGaVJQSik1TpOCUkqpcZoUlFJKjdOkoNRFcDqdVoeglEf5w5BUpXxK38Awv3u9hjf2NNDW1U9mcgxXrMznY5fPIypCf6SUf9OaglIz0NTay9d++AZPv1bFnKx4PrZhHtlpsfz6pSPc98M3aW4/Y3WISl0U/bNGKTc1t5/h/zz6Ng6nkwe/fBkLi1LH91VWt/Lg1p3885a3+cHXPkJqYrSFkSp14bSmoJQbBoZGeOAXOxi2Oz6UEADKStL5/r2X0ts/zINbdzI8onM3Kv+kSUEpNzzx5yMcbTrNN+5cQX5m/JTHzM1N5Gt3LOfI8U6efq16liNUyjM0KSg1jeoTnTy/vZbr1xayckHmeY+9rCyXy8vzePq1KuqbT89ShEp5jiYFpc7D6XTy0+cOkhgXyWdvXOjWOX/70cVERYTxs+cPejk6pTzPrY5mY0wpsBVIBdqBzSJSPemYzwFfBxyMLu7wYxF5xLXvO8CXgCbX4e+IyJc9cQNKedNeaeVgXTv3fHwpsdHhbp2TGBfJbRtL+PkfD7G/to0l89K8HKVSnuNuTeExYIuIlAJbgMenOOZ3QJmILAPWAf9ojFk6Yf8vRWSZa9OEoHye0+nkV38+REZKDNesnjOjc2+8bC4pCVE8+eIRL0WnlHdMW1MwxmQA5cDVrqJtwKPGmHQRaR07TkQmNqDGAOHAjB73NMYkAUmTinURBWWJiqpWahq6+eqnlhEeNrOW1sjwUD52+Tx++txBquo7KS1I9lKUSnmWO9/p+UCjiNgBXF+bXOVnMcbcYow5CBwHHhKR/RN232GM2WeMedkYs/Yc17oPODpp06U4lSWefauWpPhILl9xYX+XXLN6DrFRYTzzZq2HI1PKezza0Swiz4nIIqAUuMsYY1y7HgOKRGQp8BDwrDEmdYqPeBgomrSt92SMSrnjxKkedh9p4cZLiwgPC72gz4iJCue6tYW8U9nIqY4+D0eolHe4kxROALnGmFAA19ccV/mURKQe2AHc5HrfLCLDrtevuM5dPMV5XSJybOIGNMzslpS6eH965yjhYSFct6bwoj7nxkvnAvDKB8c9EJVS3jdtUhCRFqAC2OQq2gTsndifAGCMmT/hdRpwBbDf9T53wr5lQCEgFxm7Ul4xNGznjT0NrFuSQ1J85EV9VnpyNOXzM3llRz12u8NDESrlPe42H90DfMUYUwV8xfUeY8wLxpiVrmO+aIw5aIypAF4DHhWRl137HjDGHDDGVAI/Bu4SkWbP3YZSnvP+gZOc6R/mqks+1G12Qa5ZPYeO0wPsPtLikc9Typvcek5BRI4Aq6cov2HC66+f5/zPXlB0Slng1R31pCdHs7Q43SOft2phJsnxkbz8wXEuWZTlkc9Uylv0iWalJmjt7KeiupWNKwsICbF55DPDQkPYUJ7H7iOn6Okb8shnKuUtmhSUmuDNvQ04nbBxlWeajsZsWJ7HiN3Ju/uapj9YKQtpUlBqgrcrGyktSCIrNdajnzsvL5Hc9Fje2tvo0c9VytM0KSjlcrLtDLUN3VxWljv9wTNks9n4yPI89te20d7d7/HPV8pTNCko5fJ25ehf8ZcuzfHK539keS5OJ2yv0CYk5bs0KSjl8s6+JkxBMhkpMV75/LyMeIpyErRfQfk0TQpK8demo0vLvFNLGLN2cTZHjnfQ2TPg1esodaE0KSjFaC0BvNd0NGbNkmycTthx8JRXr6PUhdKkoBSw42Az8/ISvdZ0NKYwO4GMlBjeP3DSq9dR6kJpUlBBr7t3kCPHO7hkofefNrbZbKxZnEVldSt9A8Nev55SM6VJQQW93UdO4XQyK0kBYM3ibIZHHOw9e05JpXyCJgUV9HYcPEVKQiRzcxNn5XoLC1OIj4nQJiTlkzQpqKA2POJgj7SwamGWx+Y6mk5oaAgrF2Sw+0gLdseMVqxVyus0KaigdrCujf7BEVYtyJzV666Yn0lP3xA1Jzpn9bpKTUeTggpqOw+dIiIshLJSz0yT7a7lJoMQG7rGgvI5mhRUUNt56BRLS9KJinBraRGPSYiNoKQgmV2H9XkF5Vs0Kaig1dx+hpPtZ1gxP8OS669ckElNQxfdvYOWXF+pqWhSUEGrsnp0SGhZyew2HY1ZMT8DpxP2iDYhKd/hVp3ZGFMKbAVSgXZgs4hUTzrmc8DXAQcQCvxYRB5x7QsFHgGuA5zAgyLyE0/dhFIXorK6jZSESPIy4iy5/rzcJJLiItl9uIUrVnh2UR+lLpS7NYXHgC0iUgpsAR6f4pjfAWUisgxYB/yjMWapa9+dQDFQAqwFvmOMKbyYwJW6GA6Hk301rSwtScdmm52hqJOFhNhYbtLZIzo0VfmOaWsKxpgMoBy42lW0DXjUGJMu8tdHMkXk9ITTYoBwRmsFALczWnNwAK3GmGeA24CHJl0rCUiaFEKe+7ejlHuON5+mu3eIsmJrmo7GrJifyeu7G6ht6KK0INnSWJQC92oK+UCjiNgBXF+bXOVnMcbcYow5CBwHHhKR/a5dBa6yMfVTnQ/cBxydtG1371aUcl9ldRtgXX/CmLHrj/VvKGU1j3Y0i8hzIrIIKAXuMsaYGX7Ew0DRpG29J2NUCkZ/CeekxZKeHG1pHEnxkRTlJFBRpUlB+QZ3OppPALnGmFARsbs6jXNc5VMSkXpjzA7gJkAYrRnMAXa6Dplccxg7rwvomlg287yi1PmN2B0crGvj8nLf6NwtK0nnj28fZWBoZNafl1BqsmlrCiLSAlQAm1xFm4C9E/sTAIwx8ye8TgOuAMaaj54G7jbGhBhj0oFbGe2YVmrWVdd30T9ot7zpaMzy0gxG7A4OHe2wOhSl3BuSCtwDbDXG3A90ApsBjDEvAPeLyC7gi8aYa4BhwAY8KiIvu87/FbAaGBvG+l0RqfPQPSg1I5U1rdhssKQ4zepQAFhYlEJYaAiVVa2UG2sepFNqjFtJQUSOMPpLfXL5DRNef/0859uBey8kQKU8rbK6laKcRBJiI6wOBYCoyDAWFKZQoZ3NygfoE80qqAwMjXDkWKfPNB2NKStNo66xW6e8UJbTpKCCyqGjHYzYHZSV+EbT0ZhlriS1zzVUVimraFJQQWVfdSuhITYWFaVaHcpZivOTiY0K0yYkZTlNCiqoVNa0YeYkExXpW0M/Q0NsLC1Jp6KqBadTp7xQ1tGkoIJGb98QtQ1dPtefMKasJJ2Wzn6a2/usDkUFMU0KKmjsr23D6bR+aotzWeZa/a2iSqfSVtbRpKCCRmV1G5ERoT478VxOWixpSdHj8zIpZQVNCipoVFa3smhuKuFhvvltb7PZWFqcxr6aNhw6lbayiG/+dCjlYe3d/TS09Fo+VfZ0ykrS6ekb4tjJ09MfrJQXaFJQQeGvU2X71vMJk43Fp1NpK6toUlBBobK6lfiYCIpyEq0O5bxSE6PJy4jTpKAso0lBBTyn08m+6laWFqcREmLN0pszsbQ4jYN17QyPOKwORQUhTQoq4DW1naGte8Dnm47GlJWkMzBkp6q+0+pQVBDSpKAC3lhTjK8+nzDZkuI0bDbYV6NDU9Xs06SgAl5ldStpSdFkp8VaHYpb4mMimJebqP0KyhKaFFRAczic7K9po6wkDZvN9/sTxpSVpCPHOxgYHLE6FBVkNCmogFbX1E1P3zBLffz5hMmWFqczYnfqEp1q1mlSUAFt33h/gn90Mo8ZXaLTxr4abUJSs8ut+YONMaXAViAVaAc2i0j1pGO+BdwBjLi2fxGRl1z7fgFcBYz1nD0tIt/3xA0odT6VNW3kZcSRmhhtdSgzEhUZhpmTov0Kata5W1N4DNgiIqXAFuDxKY7ZAawSkTLg88BTxpiJP4kPisgy16YJQXnd8IiDg3XtfjPqaLKyknRqG7vp6RuyOhQVRKZNCsaYDKAc2OYq2gaUG2PO+kkTkZdEZGwi+H2AjdGahVKWqKrvZHDI7ndNR2OWFqfhdMJ+HZqqZpE7zUf5QKOI2AFExG6MaXKVn6tuuxmoFZGGCWX/YIz5IlAL/LOIHJ58kjEmCUiaVJznRoxKfUhldSshNlgyzz+TQmlBMlERoVRWt7JuaY7V4agg4fGOZmPMBuB7wKYJxf8KFIvIEuD3wIvGmNApTr8PODpp2+7pGFVwqKxuZW5eEnExEVaHckHCw0JYNDdVH2JTs8qdpHACyB37Je76muMqP4sxZi3wBHCriMhYuYg0iojD9fqXQBxT1wAeBoombetnckNKAfQPjiDHOykr9s9awpiyknQaWnpp7+63OhQVJKZNCiLSAlTw17/8NwF7ReSspiNjzCrgKeCTIrJn0r7cCa+vBexA4xTX6hKRYxM3oGHycUpN52BdO3aH0287mceMxa+rsanZ4taQVOAeYKsx5n6gk9E+A4wxLwD3i8gu4H+AaOBxY8zYeXeJyH7XuZmAAzgN3CIi+qim8prK6lbCQkNYUJRidSgXpTA7gfiYCCqrW7lyZb7V4agg4FZSEJEjwOopym+Y8HrVec6/6oKiU+oC7atuY0FhClER7v7d45tCQlxLdFa34nQ6/WqqDuWf9IlmFXC6ewepa+r226Gok5WVpNHWPcDJtjNWh6KCgCYFFXD2144tvenf/Qlj/tqvoE83K+/TpKACTmV1G9GRYZTkT37kxT9lp8WSlhStnc1qVmhSUAGnsrqVRXNTCQ0NjG9vm83Vr1DThsPhtDocFeAC46dGKZeWzj5Otp0JmKajMWUl6fT0DXHs5GmrQ1EBTpOCCij+OlX2dMbuR/sVlLdpUlABpbK6jcS4COZkJVgdikelJkaTlxGnSUF5nSYFFTCcTicV1a2UlaQTEhJ44/mXFqdxsK6d4RGH1aGoAKZJQQWM+uYeunoGWRZg/QljykrSGRiyU1XfaXUoKoBpUlABo2KsP6E0MJPCkuI0bDZ01lTlVZoUVMCoqGolJy2WjOQYq0PxiviYCOblJmq/gvIqTQoqIIzYHRyobQvYWsKYspJ05HgHA4M6n6TyDk0KKiDI8U4GhuwB258wZmlxOiN2J4eOdlgdigpQmhRUQBhbenOpny+qM52FRSmEhdrYV6NNSMo7NCmogFBR1Upxvv8uvemuqMgwzJwU7VdQXqNJQfm9voFhpL4z4Ka2OJeyknRqG7vp6RuyOhQVgDQpKL93oK4dh8PJsgDvZB6ztDgNpxP269BU5QWaFJTfq6hqJSI8lPlz/HvpTXeVFiQTFRGqTUjKK9xaq9AYUwpsBVKBdmCziFRPOuZbwB3AiGv7FxF5ybUvBvg5sMK17xsi8kdP3YQKbhVVrSwqSiEiPNTqUGZFeFgIi+am6kNsyivcrSk8BmwRkVJgC/D4FMfsAFaJSBnweeApY0y0a983gB4RKQZuBn5ijIm7uNCVgvbufk6c6gmapqMxZSXpNLT00t7db3UoKsBMmxSMMRlAObDNVbQNKDfGnPVTKCIviUif6+0+wMZozQLgdkYTC64axi7g+ouOXgW9sdXIgqWTecxfl+jU2oLyLHeaj/KBRhGxA4iI3RjT5Co/V6PmZqBWRBpc7wuA4xP217vOP4sxJgmYvIZinhsxqiC1V1pIjIugKCfR6lBmVWF2AvExEVRWt3Llyg/9KCl1wdzqU5gJY8wG4HvA1Rdw+n3Atz0bkQpUDoeTPdJC+fyMgJwq+3xCQlxLdFa34nQ6sdmC6/6V97jTp3ACyDXGhAK4vua4ys9ijFkLPAHcKiIyYVc9MGfC+4KpzgceBoombevdiFEFoZqGLk6fGWKFybA6FEuUlaTR1j3AybYzVoeiAsi0NQURaTHGVACbGP2FvwnYKyJnNR0ZY1YBTwGfFJE9kz7maeCLwC5jTAmwyvU5k6/VBXRN+lz370YFlT3Sgs0Gy4M2KYz2K+yVFnLSddyG8gx3Rx/dA3zFGFMFfMX1HmPMC8aYla5j/geIBh43xlS4tiWufQ8BScaYGuCPwN+JSI/H7kIFpd2HT1Gcl0RiXKTVoVgiOy2W7NRYdh1psToUFUDc6lMQkSPA6inKb5jwetV5zj8D3HYhASo1lZ6+IarqO/nUVcFbk7TZbKxYkMHLH9QzOGwnMkie01DepU80K79UIa04nLBiQXA2HY1ZuSCToWE7B2p1aKryDE0Kyi/tOnKK+JhwSvKTrQ7FUovnpRERHsquw6esDkUFCE0Kyu+MDUVdXppBaJANRZ0sMjyUpcVp7Dp8CqfTaXU4KgBoUlB+52hTN109g5TPD+6mozErF2TS3N5Hkw5NVR6gSUH5nT0yOtqmPEiHok62wpUctQlJeYImBeV3dh9pYW5uIskJUVaH4hOyUmPJz4zTpKA8QpOC8iunzwxx+FgHKxdkWh2KT1kxP5MDte30D45YHYryc5oUlF/ZdfgUDoeT1YuyrA7Fp6xckMmI3cE+XXhHXSRNCsqv7DjYTEpCJMV5kyfTDW4Li1KJjgxlpzYhqYukSUH5jeERO3vkFKsWZgXdrKjTCQ8LYVlpBjsPjdaklLpQmhSU39hf007/oJ01i7OtDsUnrVmcTcfpAWoauqY/WKlz0KSg/MYHB08SGTH6sJb6sFULMwkJsfH+gZNWh6L8mCYF5RecTic7DjZTbjKI0InfphQfE8HiuamaFNRF0aSg/EJdYzdt3QNcslBHHZ3PmsXZnDjVS0OLzkyvLowmBeUXPjjYjM022kSizm314tGk+cGBZosjUf5Kk4LyCx8cbGb+nJSgXVDHXRnJMczLS9QmJHXBNCkon9fcfoa6xm4ddeSmNYuzkfpOOk4PWB2K8kOaFJTPe3dfEwCXluVYHIl/WLM4G6dz9EE/pWbKreU4jTGlwFYgFWgHNotI9aRjrgEeAJYAPxKRb0zY9x3gS0CTq+gdEfnyRUevgsLblU0U5yeRmRJjdSh+YU5WPFmpMby3/yTXrS20OhzlZ9ytKTwGbBGRUmAL8PgUx9QBdwMPneMzfikiy1ybJgTllpaOPqpPdHHZUq0luMtms7FuSQ6V1a2cPjNkdTjKz0ybFIwxGUA5sM1VtA0oN8akTzxORGpEZC+g0zQqj3nH1XS0TpPCjKxflovd4eS9/drhrGbGneajfKBRROwAImI3xjS5ymcyJeMdriamZuDbIvLe5AOMMUnA5JnO8mZwDRVg3tnXxNzcRLLTYq0Oxa/My0skOzWWtysauXbNHKvDUX5ktjqaHwOKRGQpo81LzxpjUqc47j7g6KRt+yzFqHxMa2c/cryTy7SDecZsNhuXLcthX00r3b2DVoej/Ig7SeEEkGuMCQVwfc1xlbtFRJpFZNj1+hXXuYunOPRhoGjStt7d66jA8u5+16gjbTq6IOuX5eJw/nX0llLumDYpiEgLUAFschVtAvaKiNtNR8aY3AmvlwGFgExxrS4ROTZxAxrcvY4KLNv3NlKUk0BOepzVofilwuwEctPj2F6hSUG5z60hqcA9wFZjzP1AJ7AZwBjzAnC/iOwyxlwG/AZIAGzGmDuAL4jIS8ADxpgVgB0YAu4SER1Erc6pqbUXqe/kczctsjoUv2Wz2Vi/LJenXhU6Tw/omtbKLW4lBRE5AqyeovyGCa/f5hydwiLy2QsNUAWnN/c0YLPBhvLc6Q9W57R+WQ6/eUV4u7KJm9fPtToc5Qf0iWblc5xOJ6/vaWDJvDRSE6OtDsevFWQlUJSTwBt73O4CVEFOk4LyOVX1nZxsO8MVK3Q0sidcuTKfqvouTpzS6bTV9DQpKJ/zxu4GIsJCWLtERx15wobleYSE2Hh9t9YW1PQ0KSifMmJ3sL2ykUsWZREbHW51OAEhOSGKcpPB67tOYHc4rQ5H+ThNCsqn7JUWunuHuLxcm4486cqV+bR1D3Cgps3qUJSP06SgfMorO+pJjIugfL6usOZJqxdlERsVxmu76q0ORfk4TQrKZ3T2DLDjYDNXriwgPEy/NT0pIjyUy5bl8u7+k/QP6pyV6tz0J0/5jL/sHG3zvmZ1gdWhBKSNKwsYHLKzvaLR6lCUD9OkoHyC0+nk5Q+Os2huKnkZ8VaHE5DmFyZTkBXPi+8dszoU5cM0KSifcKCunaa2M1pL8CKbzcZ1awqpPtFFTUOX1eEoH6VJQfmElz84TmxUmC6m42VXrMwnIjxUawvqnDQpKMv19g3xbmUTG8rziIpwd45GdSHiosP5yLJc3tzTQN/AsNXhKB+kSUFZ7uUP6hkacegi87Pk+nWFDAzZeXOPzkqvPkyTgrKU3eHkT+8eZfG8VIpyEq0OJyiU5CcxNyeRF949htOpTzirs2lSUJbaeaiZlo4+brpMp3WeLTabjevXFXLs5GkO1LVbHY7yMZoUlKWe315HWlI0axZlWR1KULliZT7xMRE8+2at1aEoH6NJQVnm+MnT7Ktp44Z1hYSG6rfibIoMD+WGSwvZcaiZptZeq8NRPkR/EpVlnn2rloiwEK5ZPcfqUILSjeuKCA0J4bntdVaHonyIJgVlifbufl7f3cDGSwpIjIu0OpyglJwQxeXleby6s56eviGrw1E+wq1B4caYUmArkAq0A5tFpHrSMdcADwBLgB+JyDcm7AsFHgGuA5zAgyLyE4/cgfJLz71Vh8Ph4OOXF1sdSlD76IZ5vLqznhffO8ZtG0utDkf5AHdrCo8BW0SkFNgCPD7FMXXA3cBDU+y7EygGSoC1wHeMMYUzjlYFhN6+If783lEuK8slKzXW6nCCWmF2AstL03luex0DQzp7qnIjKRhjMoByYJuraBtQboxJn3iciNSIyF5gqu+s24Efi4hDRFqBZ4DbprhWkjGmcOIG6GorAeaFd4/RP2jnE1eWWB2KAj51VSldPYO8/P5xq0NRPsCdmkI+0CgidgDX1yZXubsKgInfcfXnOP8+4OikbfsMrqN83MDQCM9vr6N8fgZzc/VhNV+weF4ai+el8rvXqxkatlsdjrKYr3U0PwwUTdrWWxqR8qgX3jlKV+8gt1+l7de+5I6rDR2nB3llh67MFuzc6Wg+AeQaY0JFxO7qNM5xlburHpgD7HS9n1xzAEBEuoCz5vQ1xszgMsqX9Q0M89u/1FBuMlhYlGp1OGqCpcVpLChM4bevVXHN6gLCw0KtDklZZNqagoi0ABXAJlfRJmCvq2/AXU8DdxtjQlx9EbcCv5tpsMq/Pf92HT19Q9x53XyrQ1GT2Gw27rjG0NY9oLWFIOdu89E9wFeMMVXAV1zvMca8YIxZ6Xp9mTGmAfgH4IvGmAZjzLWu83/F6OikauB94Lsiok/MBJHe/mH+8EYtlyzMorQg2epw1BSWl6azoDCF37wsuo5zEHPrOQUROQKsnqL8hgmv3+YcI4VcndP3XmCMKgD8/vVqzvQP8+lrtTnQV9lsNj5/8yK++aPtPPNmLZuu0f+rYORrHc0qALV09vHsm7VcXp7HvLwkq8NR5zG/MIV1S7P5/evVdPYMWB2OsoAmBeV1v3rhMAB33bDA4kiUOzbfsJChEQe/eVmsDkVZQJOC8qqq+k7e2NPARzfMIyM5xupwlBty0+O4bs0cXnz/OCdO9VgdjpplmhSU1zgcTn7y7AGS4iL5pD697Fc2XTOf6MgwHvv9Pl2dLchoUlBe89rOeg4f62DzDQuIiQq3Ohw1A0nxkdx1/QL21bSxvaLR6nDULNKkoLyiu3eQn//xIAuLUti4qsDqcNQFuG5tIcV5ifz0uQP0DQxbHY6aJZoUlFf84o+H6BsY4UufKCMkxGZ1OOoChIbYuPcTZXT2DPLrl7TTOVhoUlAet7+mjVd31nPrhnnMyU6wOhx1EUoLkrluTSHPb69FjndYHY6aBZoUlEf1DQzz8FN7yU6N5Y6r9eGnQPDZGxeSkhjNf23by6DOohrwNCkoj/rZ8wdp7ezjvk3LiYp064F55eNio8P52u3LaGzt5Yk/H7Y6HOVlmhSUx+w+coqX3j/OxzYU6yyoAWZZaQbXryvk2bdqOVjXbnU4yos0KSiP6Dw9wMO/2Ut+ZrzOghqgPnfTIjJTYvjhr3fT0zdkdTjKSzQpqItmdzj5wZO76RsY4f/ctZKIcJ2LPxBFR4bxzc+spOP0AP/9m736UFuA0qSgLtpvXhb21bRx78eX6GijAFdakMznblrEBwebeebNWqvDUV6gSUFdlD3SwlOvCleuzOeqS+ZYHY6aBTevn8vaJdls/dMhDh3V/oVAo0lBXbATp3r4z1/upCAznns/vtTqcNQssdlsfPX25WSkxPDAL3ZwqqPP6pCUB2lSUBfk9JkhvvfTDwgLC+FbX1ijw0+DTFx0ON/6/GpGRhx876fv6zQYAUSTgpqx4REH/751B23d/fzr36wmM0WnxA5G+Znx/NNnV3GipZeHntiN3e6wOiTlAW79eWeMKQW2AqlAO7BZRKonHRMKPAJcBziBB0XkJ6593wG+BDS5Dn9HRL7siRtQs8tud/B/n9zNgdp2/vHT5SwoSrE6JGWhZaUZ3PPxpfzPbyv50dMVfPVTy3WuKz/nbp3/MWCLiDxhjPkM8Dhw5aRj7gSKgRJGk8deY8yrInLMtf+XIvIND8SsLOJwOHn06Ure2dfEF25ZzOUr8q1Pn8gxAAAOOUlEQVQOSfmA69cW0nV6gF+/LMREhXP3Rxdjs2li8FfTNh8ZYzKAcmCbq2gbUG6MSZ906O3Aj0XEISKtwDPAbZ4MVlnH6XTy0+cO8OrOeu642nDrhnlWh6R8yB3XGG75yFye317Hky8e0WcY/Jg7NYV8oFFE7AAiYjfGNLnKWyccVwAcn/C+3nXMmDuMMdcAzcC3ReS9yRcyxiQBk1d2z3MjRuVFDoeTx/+wjxfePcYt6+fy6Wt1ojt1NpvNxt/espj+gRGeerWKEbuDz964UGsMfmi2OpofA4pEZCnwEPCsMWaqyXHuA45O2rbPUoxqCna7g/9+ai8vvHuMj11ezN9q04A6B5vNxt/ftozr1hbyu9drePwP+3E4tMbgb9xJCieAXFdH8liHco6rfKJ6YOLTSwVjx4hIs4gMu16/4ipfPMW1HgaKJm3r3b0Z5VkDQyP8x6928ZddJ/j0tfP53E36l586v5AQG1/6xFJu3TCPP71zlP/atofhEZ1u259M23wkIi3GmApgE/CE6+teV7/BRE8Ddxtjfs9oR/OtwEcAjDG5ItLoer0MKAQ+tJSTiHQBXRPLjNGmCit0nB7gez/7gNqGLu7+6GJu+Yj2ISj32Gw2Pn/zIuJiwnniz0do6ezjX/7mEhLjIq0OTbnB3dFH9wBbjTH3A53AZgBjzAvA/SKyC/gVsBoYG6r6XRGpc71+wBizArADQ8BdItLsoXtQHlbX2M33fvo+vf3D/OvfXMLqxdlWh6T8jM1m4/arDNmpsTz8m71885HtfOsLq8nPjLc6NDUNm6+PEjDGFAJHX3vtNfLytM/Zm5xOJ6/sqOfx3+8jPjaC+7+whrm5iVaHpfzckeMdfP9nOxgcHuHLn1zGhnL9OZ4NDQ0NbNy4EUb7c4+5e54+0awA6B8c4Yfb9vCj/61gQVEK/3XfBk0IyiPmz0nhv76+gcLsRH7w5G4efbpCl/X0YTphjeJgXTuPPLWXk+1n+PS18/nUVaWE6lOpyoPSkqJ54EuX8uSLR/jtX6o5fKyD++5YTkl+stWhqUm0phDEBoZG+PGz+/nn/3mbEYeT799zKZuuMZoQlFeEhYbw2RsX8p2713Cmf5hvPLKdX75wSEcn+RitKQQhp9PJ+wea+elzBzjV0ccN6wr5m5sWEa0znapZsGJ+Jo9+80p+9twBnn6tmnf3neTuWxezYn6m1aEpNCkEnROnevh/z+ynoqqV/Mx4Hrj3UpYUp1kdlgoycdHhfPX25VxWlsvjf9jHd378PqsXZfGFWxaTnRZrdXhBTZNCkGjp6OOpV6t4dWc90ZFh3H3rYm5YV0RYqLYgKuuUz8/g0W9ewXNv1fHUq8KX/vMvXLtmDrdtLCE1Mdrq8IKSJoUA19bVz/++VsUrHxwHbNywrpA7rjb6IJHyGeFhoXziyhKuWJnPr186wovvHePlD45z/dpCPnllCckJUVaHGFQ0KQSo6hOdPPdWHdsrGrHZ4OrVc/jUxlLSkvSvL+WbUhKi+PvblvHJK0v431er+OM7R3nh3WNcXp7HRzfMozA7weoQg4ImhQAyPGLn/f3NPP92HYePdRAdGcaNlxXx0fXzyNDV0ZSfyEqN5au3L+eTG0t49s1aXtt1gld31rOsJJ0bLi1k5YIswsO02dNbNCn4OafTSW1DN6/urOfNPQ309g+TmRLD3350MVdfUkBMVLjVISp1QXLS4rj3E2V85voFvPjeMf70zlEe+MVOEuMiuLw8n6svKWCO1h48TpOCH3I6nRxtOs27+5p4d38TJ071Eh4Wwtol2WxcVUBZSbo+a6ACRnxMBLdtLOXjlxezt6qVV3Yc50/v1PHsW7UUZiewbkk268pyKMiM11l8PUCTgp8YGrZz+GgHu46c4r39JznV0UeIDRbPS+Pmy+ayfnkecdFaK1CBKzQ0hJULMlm5IJPu3kHe2tvIO/ua2PaK8OuXhbyMOFYvyqJ8fgYLClO1iekCaVLwUQ6Hk+PNp6moamWvtHDwaAdDw3bCQm2UlaTzqatKWb0oS0cRqaCUGBfJzevncvP6uXScHuC9/Sd5d18Tz7xZy+9eryEqIpQlxWmUmwwWz0ujIDOeEK09u0WTgo/o7R+m6ngnR453cPhYB1X1nfQNjACQnxnHtWvmsKw0ncVzU7WfQKkJUhKiuPHSIm68tIi+gWH21bSxR1rYKy3sPHQKgNioMOYXprCwKJUFRSkU5yXpE/znoP8qs8zpdNLWNcDRk90cazrN0aZujjadpqmtF6cTQmxQkJXAhuV5zC9MZmlxug4jVcpNMVHhrFmczRrXGiDN7Wc4dLSDQ0fbOXS0g91HDgNgs0FOWixzc5OYm5vI3NxE5uUmas0bTQpeMzhsp7ntDE1tvTS2nqGptZemtjMcP3ma3v7h8eMyU2IozE5gQ3keCwqTKS1I1pqAUh6SlRpLVmosV67MB6Cnb4jDxzqobeimrrELOd7B9orG8eMTYiPIy4gjNz2OvIx48jLjyMuIIzM5htAgefpfk8IF6hsYprWrn9bOftq6RrdW19eT7Wdo6+pn4vpFSfGR5KTFcmlZDkXZCRTmJFKYnUCsdg4rNWviYyK4ZGEWlyzMGi/r6RuirnG0xt7Q0kNDSy87D53ilR3148eEhNhITYwiIzmGzJQY0pOjyUyOISMlhvSkaJITogKmOSow7sIDnE4n/YMjdPUO0t0zRFfv4Ojr3kG6ewbPet/RPcAZV3v/mBDbaNtmenIMCwtTyc2IIyctltz0OHLSY/Wvf6V8VHxMBGUl6ZSVpJ9V3tM3RGNLLw0tPTR39NHS0UdLZz/7atro6O7HMWnRyujIUJLjo0hOiCI5PpKUhCiSXF8TYiOIj40gISaCuJgIYqPDfXbYuFtJwRhTCmwFUoF2YLOIVE86JhR4BLgOcAIPishPptvnTcMjDj44eJLu3iHO9A/T2z/s+jr6fmLZmf7hD/0nj4mNDicpLpKk+EjyM+MpK04nPTmatKTRLT0phpSEyKCpXioVDOJjIphfmML8wpQP7RuxO2jr6qels4/Wzn46ewbpPD1AZ88gHacHqGvsZveRFvoHR6b45NE+jdio8AmJYvR1bFQ4MVFhxESFExsVRmJcJKsXZc3q7xZ3awqPAVtE5AljzGeAx4ErJx1zJ1AMlDCaPPYaY151rQ16vn1es/NQM//xy13j7yPCQ4mLDiM2OoK46HCS4qPITY8nLiac2OhwYqPCSYqPIDEucjwJJMRG6nhnpdRZwkJDxvsrzqd/cITOngF6zgzR0zdMT9/Q2a9d77t7B2lo6aVvYJgzAyM4JvyF+m9/t5Zyk+HtWxo3bVIwxmQA5cDVrqJtwKPGmHQRaZ1w6O3Aj0XEAbQaY54BbgMemmaf16xbmsNP/vVqIsJDiIsOJzws1JuXU0qps0RHhhEdGQczWLLE6XQyOGSnb3CEEbuDjOTZnbfMnZpCPtAoInYAEbEbY5pc5ROTQgFwfML7etcx0+0bZ4xJApImFee5EeM5ZepEcEopP2Kz2YiKDCPKoo5rX+tovg/4ttVBKKVUsHKnsfwEkOvqLB7rNM5xlU9UD8yZ8L5gwjHn2zfRw0DRpG29GzEqpZTygGlrCiLSYoypADYBT7i+7p3UnwDwNHC3Meb3jHYm3wp8xI19E6/VBXRNLDPGzOiGlFJKXTh3h9XcA3zFGFMFfMX1HmPMC8aYla5jfgXUAdXA+8B3RaTOjX1KKaV8hFt9CiJyBFg9RfkNE17bgXvPcf459ymllPIdOgBfKaXUOF8bfTSVUIDm5mar41BKKb8x4XfmjB7Q8oekkA1w5513Wh2HUkr5o2yg1t2D/SEp7GR0WOpJwG5xLDORB2xnNPYGi2OZLcF2z8F2v6D37E/3HMpoQtg5k5N8PimIyCDwttVxzNSEobQN3p7jyVcE2z0H2/2C3rMf3rPbNYQx2tGslFJqnCYFpZRS4zQpKKWUGqdJwXu6gH9j0rQdAS7Y7jnY7hf0ngOezek8x3JjSimlgo7WFJRSSo3TpKCUUmqczz+nECiMMZcDrwFfE5FHLQ7Hq4wxW4CNwCDQy+g97zr/Wf7HGFMKbGV0Ovh2YLOIVFsblfcYY1IZnfF4HqP/tzXAF6eYRj/gGGO+DXwHWCIiBywOx6u0pjALjDHxwH8Af7Y6llnyZ0Z/eMqAfweesjgeb3kM2CIipcAW4HGL4/E2J/CfImJEZCmjD0Y9aHFMXmeMKQfWMLpYWMDTpDA7fgg8BLRZHchsEJE/isiw6+17QJ4xJqC+14wxGUA5sM1VtA0oN8akWxeVd4lIh4i8MaHofc5eUTHgGGMiGU34X2I0KQa8gPpB9UXGmOuBJBH5rdWxWOTvgT+JiMPqQDwsH2h0rRUytmZIk6s84LmS/L3Ac1bH4mXfBZ4QkaNWBzJbtE/hIhlj9jC65vSUuxmtXl89exF53zT3nDn2i9IYcwfwaaZYelX5vR8x2l8UsP1jxpi1wCrgn6yOZTbpcwpeZIy5DPg90OcqSmO0g+6/ReS7lgU2C4wxHwN+AGz0w0nEpuVqPqoCUkXEbowJZbSzuSTQO16NMT8AlgI3uyasDEjGmH8CvgoMuYrygFPA50TkZcsC8zJNCrPIGPMLYFcQjD66idG/JK8WkRqr4/EWY8wbwE9E5AljzGeAL4jIFRaH5VXGmO8D64AbRaRvuuMDiTHmGHBToI8+0uYj5Q0/Z/Svq99OmHZ4o4i0WxeSV9wDbDXG3A90ApstjserjDGLgH9htIb0ruv/9qiIfMzSwJRHaU1BKaXUOB19pJRSapwmBaWUUuM0KSillBqnSUEppdQ4TQpKKaXGaVJQSik1TpOCUkqpcZoUlFJKjfv/QnJ+KMAYHqgAAAAASUVORK5CYII=\n", "text/plain": [ "