diff --git a/examples/Basic Examples/Bayesian Binomial Regression/Bayesian Binomial Regression.ipynb b/examples/Basic Examples/Bayesian Binomial Regression/Bayesian Binomial Regression.ipynb new file mode 100644 index 0000000..65f42b4 --- /dev/null +++ b/examples/Basic Examples/Bayesian Binomial Regression/Bayesian Binomial Regression.ipynb @@ -0,0 +1,471 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bayesian Binomial Regression " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook is an introductory tutorial to Bayesian binomial regression with `RxInfer`." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "using RxInfer, ReactiveMP, Random, Plots, StableRNGs, LinearAlgebra, StatsPlots, LaTeXStrings" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Likelihood Specification\n", + "\n", + "For observations $y_i$ with predictors $\\mathbf{x}_i$, Binomial regression models the number of successes $y_i$ as a function of the predictors $\\mathbf{x}_i$ and the regression coefficients $\\boldsymbol{\\beta}$\n", + "\n", + "$$\\begin{equation}\n", + "y_i \\sim \\text{Binomial}(n_i, p_i)\\,,\n", + "\\end{equation}$$\n", + "\n", + "where:\n", + "\n", + "$y_i$ is the number of successes, $n_i$ is the number of trials, $p_i$ is the probability of success. The probability $p_i$ is linked to the predictors through the logistic function:\n", + "\n", + "$$\\begin{equation}\n", + "p_i = \\frac{1}{1 + e^{-\\mathbf{x}_i^T\\boldsymbol{\\beta}}}\n", + "\\end{equation}$$\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Prior Distributions\n", + "We specify priors for the regression coefficients:\n", + "\n", + "$$\\begin{equation}\n", + "\\boldsymbol{\\beta} \\sim \\mathcal{N}_{\\xi}(\\boldsymbol{\\xi}, \\boldsymbol{\\Lambda})\n", + "\\end{equation}$$\n", + "\n", + "as a Normal distribution in precision-weighted mean form.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Model Specification\n", + "\n", + "The likelihood and the prior distributions form the probabilistic model\n", + "\n", + "$$p(y, x, \\beta, n) = p(\\beta) \\prod_{i=1}^N p(y_i \\mid x_i, \\beta, n_i),$$\n", + "\n", + "where the goal is to infer the posterior distributions $p(\\beta \\mid y, x, n)$. Due to logistic link function, the posterior distribution is not conjugate to the prior distribution. This means that we need to use a more complex inference algorithm to infer the posterior distribution. Before dwelling into the details of the inference algorithm, let's first generate some synthetic data to work with." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "function generate_synthetic_binomial_data(\n", + " n_samples::Int,\n", + " true_beta::Vector{Float64};\n", + " seed::Int=42\n", + ")\n", + " n_features = length(true_beta)\n", + " rng = StableRNG(seed)\n", + " \n", + " X = randn(rng, n_samples, n_features)\n", + " \n", + " n_trials = rand(rng, 5:20, n_samples)\n", + " \n", + " logits = X * true_beta\n", + " probs = 1 ./ (1 .+ exp.(-logits))\n", + " \n", + " y = [rand(rng, Binomial(n_trials[i], probs[i])) for i in 1:n_samples]\n", + " \n", + " return X, y, n_trials\n", + "end\n", + "\n", + "\n", + "n_samples = 10000\n", + "true_beta = [-3.0 , 2.6]\n", + "\n", + "X, y, n_trials = generate_synthetic_binomial_data(n_samples, true_beta);\n", + "X = [collect(row) for row in eachrow(X)];" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We generate `X` as the design matrix and `y` as the number of successes and `n_trials` as the number of trials. Next task is to define the graphical model. RxInfer provides a `BinomialPolya` factor node that is a combination of a Binomial distribution and a PolyaGamma distribution introduced in [1]. The `BinomialPolya` factor node is used to model the likelihood of the binomial distribution. \n", + "\n", + "Due to non-conjugacy of the likelihood and the prior distribution, we need to use a more complex inference algorithm. RxInfer provides an Expectation Propagation (EP) [2] algorithm to infer the posterior distribution. Due to EP's approximation, we need to specify an inbound message for the regression coefficients while using the `BinomialPolya` factor node. This feature is implemented in the `dependencies` keyword argument during the creation of the `BinomialPolya` factor node. `ReactiveMP.jl` provides a `RequireMessageFunctionalDependencies` type that is used to specify the inbound message for the regression coefficients `β`. Refer to the ReactiveMP.jl documentation for more information." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "@model function binomial_model(prior_xi, prior_precision, n_trials, X, y) \n", + " β ~ MvNormalWeightedMeanPrecision(prior_xi, prior_precision)\n", + " for i in eachindex(y)\n", + " y[i] ~ BinomialPolya(X[i], n_trials[i], β) where {\n", + " dependencies = RequireMessageFunctionalDependencies(β = MvNormalWeightedMeanPrecision(prior_xi, prior_precision))\n", + " }\n", + " end\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " This example uses the precision-weighted mean parametrization (`MvNormalWeightedMeanPrecision`) of the Gaussian distribution for efficiency reasons. While this is less conventional than the standard mean-covariance form, the example would work equally well with any parametrization. The choice of parametrization mainly affects computational efficiency and numerical stability, not the underlying model or results." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Havin specified the model, we can now utilize the `infer` function to infer the posterior distribution." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32mProgress: 100%|█████████████████████████████████████████| Time: 0:00:03\u001b[39m\u001b[K\n" + ] + }, + { + "data": { + "text/plain": [ + "Inference results:\n", + " Posteriors | available for (β)\n", + " Free Energy: | Real[21992.9, 16235.8, 13785.0, 12519.7, 11800.1, 11366.2, 11094.1, 10918.7, 10803.3, 10726.3 … 10558.2, 10558.2, 10558.2, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1, 10558.1]\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "n_features = length(true_beta)\n", + "results = infer(\n", + " model = binomial_model(prior_xi = zeros(n_features), prior_precision = diageye(n_features),),\n", + " data = (X=X, y=y,n_trials=n_trials),\n", + " iterations = 50,\n", + " free_energy = true,\n", + " showprogress = true,\n", + " options = (\n", + " limit_stack_depth = 100, # to prevent stack-overflow errors\n", + " )\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now plot the free energy to see if the inference algorithm is converging." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXwUVb738VPV3Vk7ndUkkLAlQcANhaCyCAwQ9QpBQRgQwXFDHF/3+uCjw+iAOj6My1zEG/W6juCGuIDKCIqKgGIQIqCIsgkGQSCQhJCt03vV80dhG0ISOqZJdXd93n/4qjpUqn9ddve365yq05KqqgIAAKOS9S4AAAA9EYQAAEMjCAEAhkYQAgAMjSAEABgaQQgAMDSCEABgaAQhAMDQCEIAgKERhAAAQwunIHz22Wf37NnjX/V4PDoWE444Ym3l9XqZg7BNfD6foih6VxFOFEXx+Xx6VxFmvF5vcHcYTkG4cuXK3bt3+1edTqeOxYQjjlhbuVwugrBNPB4PH+tt4vP5+IbaJoqiuFyu4O4znIIQAICgIwgBAIZGEAIADI0gBAAYGkEIADA0ghAAYGgEIQDA0AhCAIChRWwQfnJQPeLQuwgAQMiL2CBcsFtZV8ZUTwCA0zDrXcCZkhglatx6FwEAequtrb399tuDPj9nx7NarQsXLjwTe47YIEyKEtUEIQDDq6ioWLNmzdNPP613Ie01adKkl1566UzsOWKDMDFKqnEzXTIACKvVOnHiRL2raK/JkyefoT1H7BghXaMAgEAQhAAAQ4voIORHvgAApxOxQZgUJVW7GCMEAJxGxAYhXaMAgEBEdBDSNQoAOJ2IvX2CrlEACCOHDx9et27d9u3bc3Nzb7zxxo586IgNQluUqPMIVQhJ70oAAKe1aNGitWvX2u32b7/9toODMGK7Rk2SiDOLOnpHASCU7NixY+XKlf7VXbt2ffjhh0KIWbNmrVy5cuzYsR1fUsQGoWByGQAIPfHx8dddd11dXZ22ev/99+/cuVPfkiK2a1T8euFol3i96wCAUPKLXb1lna9jzhKiZLFkpDmuUdR069Zt0KBBb7311vTp048ePfrxxx8/88wzHVJLiyI/CAEAjXWKlf7a16R2SBJGm0TcKTnz5z//ee7cudOnT1+4cGFhYWF6enpHlNKySA5CfoACAE5llsXIznpeRzh69Og777xzy5YtCxcuPEM/KNEmjBECADqULMu33nrrzTffbDabhw4dqnc5ER6EdI0CQCi67bbbdu/efdttt0nSiXPTf//737m5uY899tjatWtzc3PvvffeDiuGrlEAQEdzOBxRUVE33HCDv+U//uM/Gp8dRkdHd1gxkRyEiVFSNV2jABBili1b9vTTT998882pqan+xqioqKioKF3qiewgFPvr9S4CAHCy7du3FxYWzpgxQ+9CTojwIGSMEABCzezZs/Uu4SSRfLFMEl2jAIDTieQg5IwQAHBaBCEAwNAIQgCAoUXyxTKMEQKAyWQ6ePBgfn6+3oW0lyRJkiSpZ2CO1EgOQqtFOH3CqwhzJJ/3AkBrunfvXlJS4vV69S6kveLi4gjCNpOESLCIWo9I6bgJCgAg5PTt21fvEkJahJ8rMe82AKB1kR6EFq6XAQC0JsKDMCmaebcBAK2J8CCkaxQA0LpID0K6RgEArYr0IOSeegBAqyI8CBkjBAC0LsKDkDFCAEDrIj0IGSMEALQq0oOQMUIAQKsiPAiTopl3GwDQmggPQrpGAQCti/QgpGsUANCqCA9CukYBAK2L8CCkaxQA0LoID8JYs1CFcPn0rgMAEKoiPAiFEDZOCgEALYv8IGRyGQBAK4wQhKLGo3cRAIBQFflBmBQlql16FwEACFWRH4R0jQIAWmGEIKRrFADQosgPQrpGAQCtiPwgpGsUANCKQIPQ6XRu2LBh5cqVBw8ebNyuKMrGjRs/+uijmpqaxu2VlZUrVqz45ptvmuxny5YtK1asOHbsWOPG6urqjz76qKSkRFGUtj+F06BrFADQioCCcPfu3Z06dZo5c+azzz573nnnzZs3T2v3+XyFhYXTp09//vnne/Xq9cMPP2jtxcXFffr0WbBgweTJk//0pz/59zN16tTrrrtuwYIFffr0+eqrr7TGbdu29erV6/nnn7/llluuvvpqny/I08Aw7zYAoDVqACorK3/66Sdtef369bIsV1VVqaq6bNmynj17NjQ0qKo6e/bs8ePHa9sMHjz4ySefVFX1+PHjGRkZGzZs0P4wIyPj+PHjqqoWFRUNHTpU2/iaa66ZM2eOqqoNDQ15eXkffPBBS2WMGTNm+fLl/tXa2tpAin9vn2/cKm8gW0a8AI8Y/Orr630+n95VhBOHw+F2u/WuIpy43W6Hw6F3FeHE5/PV19cHd58BnRGmpqbm5ORoy7169VIUxW63CyHef//98ePHx8bGCiGmTJnywQcfeL3eo0ePrl+/fsqUKUKIpKSkq6666v3339c2Hj16dFJSkrbxunXrKisrPR7PihUrrr/+eiFEbGzsuHHjtI2DKDFKqnYxRggAaJ65rX/wxBNPDB06NDs7Wwhx8ODBAQMGaO3dunXzer1HjhwpLy+Pi4tLS0vT2rt27bpnzx4hxC+//NK7d2+t8ayzzoqNjT148KDD4fB6vV27dvVvvHXr1pYeura2tri4uKGhQVuVZXn8+PGnLTjBrNZ4xJkYfQw7iqJwHNqEI9ZWiqJIksRBC5zyK70LCRttPWKSJEmS1Po2bQvCN99885VXXvnyyy+1VbfbbbFYtGVtweVyuVwus/m33VosFqfT2WRjIURUVJS2sRDCv71/42YdP3583bp1WqwKIWw221VXXXXammNUUe1qbbfG4XK5Gv8vwGk5nU5JkmQ58i+uDhan02kymYI+0h/BPB4Ph6tNFEXRXmYBbh8VFdU4kprVhiB877337r777k8//dTfTZqZmVlZWaktV1RUCCE6deoky3JdXZ3L5YqOjtbaO3XqpP2Tf2OXy1VbW9upU6eUlBQhRGVlZefOnbWNtYVmdevWbcaMGWPGjNFW6+rq4uLiTlt2J1nUuD2BbBnxfD4fx6FNVFWNjY0lCAMny7LJZOL7VuC0IIyJidG7kLCh9ToE96Ms0Hf4ypUr77jjjuXLl5933nn+xkGDBq1du1ZbXrt27UUXXRQXF9e1a9fs7OwvvvhCa//8888HDx7cZOPPP/9c28xqtfbt27fxTgYNGhSUJ+aXGCVquX0CANCCgM4Id+7cOW7cuOHDhy9ZsmTJkiVCiNtvv7179+5/+tOf/vnPf/7lL3/p06fP7Nmzi4qKhBAmk+mee+654447HnjggeLiYrvdro3kTZgw4cEHH5wxY8bgwYMfeuihe+65R/uiPWvWrLvvvtvtdm/fvn3nzp3vvvtucJ+hRRYWWdi9Ir7N46EAgMgXUDjEx8c/9NBDJ/2Z2SyESE5O3rhx43PPPff1118vWLDAP2J35513ZmZmfvbZZxkZGcXFxVofaXR0dHFx8TPPPFNcXPzYY49NnDhR23jKlCmJiYnLly9PSkrauHGjdllpcCVGiRq3Gm8+zXgpAMCAJFUNm1sLCgsLm4wRJiQkBPKHfZZ63xtl6pNk9CAM/IhBY7fbGSNsE+0qBsYIA8cYYVspiuJwOOLj44O4T0O8w5OYXAYA0AJDBGFilKgmCAEAzTFIEPIDFACA5hkkCOkaBQA0zxBBmETXKACgBYYIQrpGAQAtMUgQ0jUKAGgeQQgAMDRDBGFSlKimaxQA0BxDBGFilMQZIQCgWQYJQrpGAQDNM0QQcvsEAKAlhghCbp8AALTEEEFos4h6j1CIQgDAKQwRhLIk4s2ijt+pBwCcwhBBKOgdBQC0wDhByIWjAIBmGCUIk6K5cBQA0AyjBGGihTNCAEAzDBOEjBECAJpjlCCkaxQA0CyjBCFdowCAZhkmCOkaBQA0xzhByBkhAKAZRgpCZpYBAJzCKEGYFCVVu+gaBQA0ZZQgpGsUANAsghAAYGhGCUJ+mxcA0CyjBCG3TwAAmmWUIIy3CJciPIredQAAQoxRglASwmYRtdxBAQA4mVGCUNA7CgBojqGCkAtHAQBNGSgIuXAUAHAqAwUhXaMAgFMZKgjpGgUANGWgIEyKEtUuvYsAAIQYAwUhP0ABADiVoYKQMUIAQFOGCkLGCAEATRkoCLl9AgBwKgMFIV2jAIBTGSoI6RoFADRloCBMIggBAKcwUBAmRknVdI0CAE5mqCDkjBAA0JSBgjDGJIQQTp/edQAAQomBglBwUggAOIXRgpA7KAAAJzFWEHJPPQCgCWMFIV2jAIAmjBaEdI0CAE5irCDknnoAQBPGCsJExggBACczWhDSNQoAOInBgtBC1ygA4CQGC0LGCAEAJzNWECZFM0YIADiJsYKQMUIAQBMGC0LGCAEAJzNWENI1CgBowlhBSNcoAKCJjg5Cp9NZUVFxantlZaXD4TjTj55oEbUeQRICAPwCDcIHH3ywX79+KSkpCxYs8DfOnDkzpZE+ffpo7aNHj/Y3Dh482L99UVFRZmbmhRde2K9fv/3792uN5eXlgwcPPv/88zt16vTQQw8F6Xk1zyyLaFnYPWf0QQAA4STQIOzZs+e8efN69+7tdDr9jUVFRVW/uuqqq0aNGqW119XVvfzyy1r7+vXrtcY9e/Y88MADJSUlhw4duuyyy+655x6t/f777+/Ro0dZWdkPP/zwv//7vyUlJcF7ds2gdxQA0Jg5wO2mTp0qhHjkkUea/deampr333+/uLi4lT288cYbV155Za9evYQQM2fO7NWrV11dXXx8/OLFi9esWSOEyM7OnjBhwqJFiy655JK2PYm2SIwSNR6RdeYeAAAQVoIzRrh48eK8vLyLLrrI3zJlypTY2NgLLrjg448/1lpKS0u1FBRC9OjRQ5KkX375pby8vL6+3t/eq1evffv2tfQoXq/3yJEjpb86cuTI7yg1KUpUu37H3wEAIlOgZ4StW7BgwfTp0/2rTz75ZJ8+fUwm06uvvnrttdd+//33OTk5NTU1cXFx/m3i4+Orq6vNZrMQwt+uNbb0KHv37p0zZ05sbKy2mpiY2Po5aLOsJsuRGnd9vNLWP4wAdrtdkiS9qwgnDQ0NPp9Plo11cXV7OJ1Ok8lksVj0LiRseDwen8/n9Xr1LiRsKIridDpVNdARrpiYGC1oWhGEIPz+++9/+OGHyZMn+1v8p4a33nrrggUL1qxZk5OTk56e7g85RVFqa2szMjJSUlKEEDU1NampqUKI6urq9PT0lh6od+/eM2bMGDNmjLZaV1dntVrbWm1KrM9tslitRvxoU1X1dxwxI5MkKTY2liAMnNlsJgjbRAvCmJgYvQsJG4qimEym+Pj4IO4zCO/wl156afz48Wlpac3+a11dnfb/+Pzzz9+0aZPWuGXLloSEhOzs7OTk5OzsbH/7pk2bzj///PaX1IokfpIQANBIoGeEJSUlBw4cKC8v//bbb5csWTJw4MDs7GwhhNvtXrx48ZtvvunfsqKi4rXXXhs6dKjZbF64cGFFRcUVV1whhJg6deqDDz740ksvDRky5K9//estt9wSHR0thPjzn/88Z86cLl26bN++/ZNPPpk3b94ZeJq/4QcoAACNtSEIi4uL+/TpU1tbu2TJki5dumhBuHPnzjFjxowYMcK/pcVi2bZt21tvvaUoSt++fb/88suzzjpLCJGSkrJy5coHH3zw6aefvvzyy+fOnattP2vWLJfLdf3116empr733nvdunUL9nM8CbdPAAAakwIfctRdYWFhkzHChISEtu7kmR3Kjmr1mUGmYFcXBn7fETMyu93OGGGbcLFMWzFG2FaKojgcjpAbIwwv3D4BAGjMcEFI1ygAoDEDBqGoYa5RAMCvDBmEXDUKAPiV4YKQMUIAQGOGC0LGCAEAjRkuCBMsosEnfEQhAEAIYcAglCURbxZ1XC8DABBCGDAIBb2jAIBGjBiEzLsNAPAzYhByBwUAwM+wQUjXKABACKMGocQZIQBAY8QgZIwQAOBnxCBkjBAA4GfMIOT2CQDACcYMQs4IAQAnGDEIGSMEAPgZMQjpGgUA+BkzCOkaBQCcYMQgpGsUAOBnxCBkZhkAgJ8xg5CZZQAAJxgxCOPNwqMIt6J3HQCAEGDEIBRC2CyilpNCAIBhg5A7KAAAGoMGYVI0F44CAIQwbBAmWriVEAAghHGDkK5RAIAQwrBBSNcoAEBj0CBkljUAgMaoQcgYIQBACGHYIOwcLx20M0YIADBqEObZpL21BCEAwLhBKPbW6l0EACAEGDQIu8RLx1yqw6t3HQAAvRk0CGVJdLdKP9XROwoARmfQIBRC5NoEw4QAAOMGYZ5N+olhQgAwPOMGYa5N+okzQgAwPOMGIXdQAACEsYOQOygAAAYOwu5W6XCD6vLpXQcAQFfGDUKzLLrESz/X0zsKAIZm3CAU9I4CAAwfhNLeGs4IAcDQDB2EuTYmlwEAozN0EHIHBQDA4EHIGCEAGJ2hg7BHgnSgXvUoetcBANCPoYMw2iQ6xUkHuIMCAAzM0EEo6B0FAMMzehDmJnC9DAAYmuGDkN+gAABjM3oQ5tnET3V6FwEA0A9ByOQyAGBoBKG0r15ViEIAMCqjB2GsWaRESwftJCEAGJTRg1BwBwUAGBtByIyjAGBoBCF3UACAoRGEdI0CgKERhHSNAoChEYSip036qVYlCQHAmAINwiVLltx+++0FBQUffPCBv3Hjxo0FjXz99ddae319/W233Zabmztw4MC1a9f6t1+9evWll16am5s7Y8YMu92uNbrd7rvuuisvLy8/P7/xzjuM1SKsFlHWQBQCgBGZA9xuy5Yt3bt3Ly4u/uWXX/yNFRUVZWVlRUVF2mpOTo62cO+99x46dGjdunVffvnluHHj9u7dm5aWVlFRMX78+BdffHHIkCHTp0+/7777nnrqKSHEI4888vXXX69evXrHjh2TJk369ttvc3Nzg/ocTy/PJu2tFZ3jOvhhAQD6CzQIH3vsMSHEqlWrmrQnJSWNGjWqcYvD4Xj11Ve//PLLrKysyZMnL1y4cNGiRTNnznz99dcvueSSSZMmCSEefvjh4cOH//d//3d0dPQLL7zw+uuvd+vWrVu3bmPGjFm4cOHDDz8cjKfWBrk2aW+tOjRT6uDHBQDorr1jhLt37x44cODo0aNfe+01VVWFEPv373c6nX379tU26Nev344dO4QQ27dvz8/P1xr79u1rt9t/+eWXqqqqI0eO9O/fv8nGHYw7KADAsAI9I2xWbm7uM888k5ubu2vXrnvuuaeqqmrmzJnHjh2zWq2SdOLsKikpadeuXUKIysrKvLw8rVGWZZvNVlFRoWVnQkKC1p6YmFheXt7Sw+3YsWPy5MkWi0VbTU1N3bp1a3vq98uOkj86ZKqvdwRlbyHLbrf7/78gEA0NDT6fT5a5pixQTqfTZDL536Q4LY/H4/P5vF6v3oWEDUVRnE6nGvAFjjExMWbzaZKuXUF4zjnnnHPOOUKI/v37u93uoqKimTNnJiUl1dfXq6qqfebW1tampqYKIZKTk+vr67U/VFW1vr4+JSUlKSlJCFFfX68t1NXVaRs3q3fv3g8//PAVV1yhrbpcLqvV2p76/c5LV5/d47Nao4Oyt5ClqmqwjphBSJIUGxtLEAbObDYThG2iBWFMTIzehYQNRVFMJlN8fHwQ9xm0d3hycrLD4RBCdOnSRZKkn376SWvftWuXdhFNbm7uzp07tca9e/dKkpSdnZ2Wlmaz2bRTxsYbN1+rLFut1uRfBfFA9LRJe/gxJgAwpECD0G63Hz9+3OPxOByO48ePu91uIURxcfGxY8eEEAcOHHjkkUe0czWbzTZu3Lh58+YpirJ169ZVq1Zdf/31QoipU6euWrVq69atiqLMmzfv2muvtVqtsixPmzbt8ccf93q9e/fufffdd2+44YYz9mRblBwtzLKocHb8IwMAdBZoEM6ZMyc/P//QoUPPPfdcfn7+6tWrhRCffPJJ9+7drVZr3759L7rookcffVTb+H/+53927dqVkpJSUFDw1FNPde/eXQjRo0ePoqKiUaNGpaSk/Pjjj/Pnz9c2njt3rt1uT0tLu/jii+fMmdOvX7/gP8sAML8MABiTFPiQY0scDkdsbGw7251OZ3R0dOuXchQWFs6YMWPMmDHaal1dnf8qm/a7fq3vyi7StLxIHg0K7hEzArvdzhhhm3CxTFsxRthWiqI4HI7gjhG262IZTbNp19Z23V8HeTbBHRQAYEB81T0h1ybxGxQAYEAE4QmMEQKAMRGEJ+TZpL3cQQEAxkMQnpAeK7yqOO7Suw4AQMciCH+TS+8oABgPQfgbhgkBwIAIwt/kJgguHAUAoyEIf8OPMQGAARGEv8mzST/VEYQAYCwE4W/ybII7KADAaAjC33SOl+o8os6jdx0AgA5EEP5GEiKHYUIAMBiC8CTcQQEARkMQniTPxh0UAGAsBOFJuIMCAIyGIDwJXaMAYDQE4UnoGgUAoyEIT9IlXjrmUh1evesAAHQUgvAksiS6W5lfBgAMhCBsimFCADAUgrAphgkBwFAIwqa4gwIADIUgbIogBABDIQibomsUAAyFIGyqu1U63KC6fHrXAQDoEARhU2ZZdImXfq6ndxQADIEgbAa9owBgHARhM/JsEj9VDwAGQRA2I9fG5DIAYBQEYTPOTpR2VROEAGAIBGEzBqZLJeVcOAoAhkAQNiM5WvROkr4q56QQACIfQdi8y7OlVQcVvasAAJxxBGHzCrLkTw9xRggAkY8gbN6gdOmnWrXCqXcdAIAzjCBsnlkWl2XKaw7TOwoAEY4gbFFBlrSK3lEAiHQEYYsKsqRPDhKEABDhCMIW9U6SzLLgznoAiGwEYWtGdpa4dhQAIhtB2JqCLGnVIa6XAYBIRhC2ZlSW/EUZc60BQCQjCFuTGi16JUobmWsNACIXQXgal2fTOwoAkYwgPI2CLJm7CQEgghGEpzEoQ9pVrVa59K4DAHBmEISnESWLIZnSauZaA4AIRRCeHr2jABDBCMLTuzxb+pS51gAgQhGEp3dOkqSo4scashAAIhBBGJBR/BIFAEQogjAg/CQTAEQqgjAgBVny52WKh0tHASDiEIQBSYsRuTaphLnWACDiEISB4pcoACAiEYSBKsiS+W1CAIg8BGGgLsuUdhxnrjUAiDQEYaCiZDEoQ/q8jN5RAIgoBGEbMNcaAEQegrANLs+WPmauNQCILARhG5ybLLl94qdashAAIgdB2AaSEKOyJK4dBYBIQhC2DXOtAUCEMQe43Y4dOzZu3Lhnz56xY8cOHDhQa9yzZ88777yzc+fOhISESZMmDR8+XGt//vnnf/75Z205IyPjrrvu0pb3798/f/78ioqKUaNG3XzzzZIkae1vvvnm8uXLk5KSZs6cefbZZwfpqZ0RBVnyf37l8SgmC18hACAiBPpxPmvWrBUrVrzxxhvffPONv3H+/Pnl5eVXXXVVXl5eYWHhsmXLtPbFixeXlZUlJycnJyfbbDat0eFwXHbZZRaLZeLEiY8//vj8+fO19oULF957773jxo1LT08fMmTIsWPHgvfsgi8jVnS3SpsqOCkEgAgR6BnhihUrhBAjR45s3Pjcc8/5z+oqKirefPPNa665RlsdP3781Vdf3Xjjt99+Oz09Xcu/pKSkG2644a677jKZTPPnz583b97EiRMnTpy4cePGV1555e67727nszqjLs+WVh1SB2VIehcCAAiCdnXw+VNQCHH48OGMjAz/6htvvHHnnXf+61//crlOzMXy9ddfDxs2TFseMmTI4cOHDx48WFdXt2PHjqFDh2rtQ4cOLSkpaU9JHaAgS2bSUQCIGIGeEbZu9erVy5cv/+6777TVYcOGJSUlmc3mBQsWvPjii+vXr4+KiiorK8vPz9c2iIqKstlsZWVlHo9HCJGamqq1p6WlHTlypKVHKS0tnT179hNPPKGtJiYmvv7660Gpv00usoofa6K2ltnzEsKsg9Rutzf+7oLTamho8Pl8ssyAcKCcTqfJZLJYLHoXEjY8Ho/P5/N6vXoXEjYURXE6naoa6MdvTEyM2XyapAtCEG7atGnKlCnvvPNO165dtZa5c+dqC7fddlvPnj0//PDDcePGxcXFud1u/1+5XK64uLjY2FhtWXvnOJ3OuLi4lh4oIyNjxIgRl1xyib+llY3PnDgh/vMc9ek98r+GhNnno8/n0+WIhS9VVWNjYwnCwMmyTBC2iRaEMTExehcSNhRFkSQp8I+yQN6/7Q3CrVu3jh079qWXXiooKDj1X2NjY/Py8srKyoQQ2dnZ+/fv19qPHj3qcrmys7NtNpvFYtm/f/+5554rhNi/f392dnZLjxUfH3/hhRf6H6iurk6vT6j/c57Ie8fzQIPUzRpOJ1iyLPOZ3ibyr/QuJGxwxNpKlmVVVTlibRL011i79vX9999feeWVRUVFhYWF/kan01ldXe3fYNOmTVqP6IQJEz788MPKykohxKuvvvqHP/whJSXFbDZfffXVr776qhCirq7uvffemzBhQntK6hiJUeKWXvIT3zNSCADhTw3MzJkzc3JyYmNj09LScnJyPvroI1VVx4wZExsbm/OrCRMmqKq6f//++Pj4Sy+9dNCgQVarde7cuf6d3H777dnZ2SNHjszMzNy8ebPWuHPnzuzs7OHDh/fo0WPy5Mk+n6+lGsaMGbN8+XL/am1tbYDFnwlHGtTk19xlDTqW0Gb6HrFwVF9f38oLEqdyOBxut1vvKsKJ2+12OBx6VxFOfD5ffX19cPcpqYENOZaXl9fX1/tXMzIy4uPjjxw50tDQ4G+MiYnp3LmzEKKmpmbHjh1CiN69eycnJzfez969e8vKyvr16xcfH+9vdDqdmzdvTk1N7dOnTys1FBYWzpgxY8yYMdpqXV1dQkJCIMWfIf/1lS/BIh4ZYNKxhjbR/YiFHbvdzhhhm3CxTFsxRthWiqI4HI7GCdJ+gY4Rpqenp6enN2nMzMxsduPExET/7DNN5OXl5eXlNWmMiYkZMmRIgJWEjll95X7ve2f1NSVF6V0KAOD34qvu79clXhrTVX5mByOFABDGCMJ2mfG9ij4AAA97SURBVH2h/NR2X71H7zoAAL8XQdgueTZpWKb80m5OCgEgXBGE7TX7Inn+94qbKASA8EQQtlffFOmCFPHaHpIQAMISQRgE919kenSr4iUKASAMEYRBcGm61MUqluwjCQEg/BCEwXFfX9M/vlWUMPs5CgAAQRgkV2RL8Rax4gAnhQAQZgjCoLm3r/zwVoIQAMIMQRg047rLDV6x+jDdowAQTgjCoJGEmNVXfnSrT+9CAABtQBAG05RceX+9+OooJ4UAEDYIwmAySeL/ni//cxsjhQAQNgjCILvpbPmbSnUNI4UAECYIwiCLMYnXhpuuW+vdV0cWAkAYIAiD7w+dpFkXmMZ/5nN49S4FAHA6BOEZcff58gUp0oxiriAFgFBHEJ4pzw82ba9Wn9vJhTMAENIIwjMl1izeHWn6f9/4vjzCYCEAhC6C8AzqniC9Ntx83VrfITtZCAAhiiA8swqypNt7yxNX+/gJewAITQThGTf7IjkrXvq/G7lwBgBCEUF4xklCLBxqWntYXfgjZ4UAEHIIwo6QYBHvFZju2+TbVMFgIQCEFoKwg/RKlF4cYpq42lfh1LsUAEAjBGHHubqbPCVXmrzG66WLFABCBkHYoebmm2JM4o9rfLUevUsBAAghCMIOZpLEvwvMZyeKAcu8248zXggA+iMIO5pZFo8NMD1wkTziI+/SfXSSAoDOzHoXYFDX58nnJkvXfubbXKk+nG8ySXoXBABGxRmhbi5Mlb6+xvxtpTrqI2+5Q+9qAMCoCEI9pUaLj640D+skXfxvL7cYAoAuCEKdmSTx936mpwbKhZ96mXoGADoeY4QhYWw3uVeSNH6Vr6RcfXqQKYrvJwDQUfjEDRW9EqWvxpqPOkT+Mu+y/Qr9pADQMQjCEJIYJZYVmJ64xPTIVqXve94l+4hDADjjCMKQMypL+vpq8xOXmB77Thn4gXf5AQYOAeAMIghD1KgsafM15tkXyg9sUQYRhwBwxhCEoUsSorCrvOUa813ny38pUYYs9645TF8pAAQZV42GOlkSE3vI47rJb/ykTP/S19Uqrs+Tx3WXU6P1rgwAIgJnhOHBLIs/9ZR3TTT/17nyZ4fU3Lc9V37sXbBbqXLpXRkAhDnOCMOJRRbju8vjuwunz7TqkLKkVL2nxHNusjSxh3xdrpweq3d9ABCGCMKwFGMShV3lwq7C7jV9eEB5Z5/64DeewRnStT3kYZlSro05vAEgUARheIs3iz/myH/MEXavafl+Zdl+9cEtikdRL02XB6ZLgzKk/mlSHP+TAaBlfEZGiHizmJwrT84VQoiDdvWro+qGcnXW18q2KvXcZGlgunRpunRBvHRugt6FAkCIIQgjUHa89Mcc6Y85Qgjh9IktleqGcnXJPvXuo9H1Xs/ZidLZiVLvJKlXotCW43kVADAwPgIjXIxJDM6QBmdI4nxRV9fgi07YXa3urlF316hL94kfa5Qfa9SzYqSzE8XZiVKOTeoUK7Ljpcw40SWePlUAhsBHnbEkRYlL0qVL0n+7mkZRxYF69cda8WON+nOduvWY+KVeOeIQB+2qRRZZcVLneNE5TsqKExmxUnK0SImWUqJFyq8LZm7AARDmCEKjkyXRPUHqniAuz2p6rWmNWxxqUA/ZRVmDetAuSuvUqkpR5VKqXKLKJapcapVLWM0iJUZKjRbJ0cJmkeLNIt4iEqNEgrZsFknRwmqW4szCahFRsojX/muWok2CM04AoYCPIrQoMUokRknnJAkhWrwfo9otqlzqMaeocYsat2r3CrtX1LpFrUetcIoGr6h2iXqv0uAV9R7hVoTdI1yKaPCqLp9o8IoYk4g1ixiTiDVJQghblDBJQhIiKVoIISyysJolIUSUSfgHMi2ysFpOLJslkWD5rTaTLGyWk8rTdn4qq1mytHwu6w9ph0OKjlZl+TQz2yV31Cw/JknYLCF9b4zLLcmysDAzbsA8HqEoUnRI/18NLYoiZEXEB3WfBCHaJSlKJEVJOScuRm3zu9npEw6vcPhUp08IIWrcQlGFoooatxBCuBVh96pCCLdP2L0n/sSjiHrPiWWvKo67f0spnyL2eBrvXji8QttzE/VexdPyh7UW0kIIRTHLsipEc7to5HhHze/jU0WtJ6Tnm1VVSZKEEJ7Tbgk/VZUkiSPWBjflmuYPCuYOCULoKcYkYkwiubUE1fOrst1uj42NlWUGQgPldDpNJpPFYjn9phBCCOHxeHw+JSYmRu9CwoaiKA6HO7j75B0OADA0ghAAYGgEIQDA0MI4CIuKitzuIPcUR7YFCxZUVlbqXUU4Wbp06d69e/WuIpx89tlnmzZt0ruKcLJly5ZVq1bpXUU4KS0tXbJkSXD3GcZB+Mwzz1RVVeldRThZtGhRaWmp3lWEk2XLln333Xd6VxFOVq1atX79er2rCCfr16//9NNP9a4inHz33Xfvv/9+cPcZxkEIAED7EYQAAEMjCAEAhiapakhPVNFY7969a2trY2NjtdUDBw5kZ2dzs3PgDh8+nJaWFhUVpXchYaO8vNxqtcbFxeldSNioqqoymUyJiYl6FxI2amtrPR5Pamqq3oWEjYaGhvr6+vT09AC3nzJlyty5c1vfJpyC8NixY9XV1doMTkIIl8sVHd1RkzxGBI5YW7ndbovF4n/J4bS8Xq8kSSZTcxO8ojk+n09VVbOZSb4Cpaqqx+MJ/At9p06d/KdPLQmnIAQAIOjoVwQAGBpBCAAwNIIQAGBoBCEAwNDC8lIlVVXffPPNzZs35+Tk3HrrrfyUV7Oqqqq2bNmyb9++Sy+99IILLvC3l5eXL1y4sKKiYvTo0SNGjNCxwlCzf//+FStWlJaWZmZmTps2LTMzU2tXFGXRokVbt27t2bPnLbfcwv0nfqWlpStWrPj5559tNtvo0aMHDBjg/6d33313/fr1Xbt2nT59enx8cH9OPBLs3bt3zZo1o0ePzsrK0lq+/fbbt956Kzo6+sYbb8zJydG3vNBRU1Pz9ttv+1cHDhx4/vnna8slJSVLly61Wq033XRT165d2/MoYXlG+Le//e3RRx/t2bPnBx98MH78eL3LCVHXXHPN3/72t7///e+Np/RtaGgYOHDg7t27u3XrNmXKlMavMEyaNGnLli1dunTZvn37Oeecs2/fPq39rrvuKioq6tmz5zvvvHP99dfrW2RIWbduXWlpaY8ePdxu98iRI5cuXaq1P/bYY/fdd19eXt7atWuvuOIKLk1vwuv1Tps2bebMmTt37tRaSkpKhg0blpyc7Ha7BwwYcODAAX0rDB1Hjx6dOXNm6a9qamq0du2llZmZWV1dPWDAgPLy8nY9jBpuampqrFbrDz/8oKqq3W5PTEz85ptv9C4qFGn3J40ePfrxxx/3Ny5YsGDAgAGKoqiqunjx4gsuuEC3+kKPw+HwL48YMeLhhx9WVbWioiImJqa0tFRV1Zqamri4uF27dulWYgh74IEHxo8fr6qq0+lMS0tbv369qqput7tz585r167VubgQ88gjj8yePTsrK2vVqlVay4QJEx544AFteerUqffee69+1YWW3bt3p6Wlndp+xRVXzJs3T1u++uqrtXfr7xZ+Z4Rbtmyx2WznnnuuECIuLm7w4MFffPGF3kWFombn3Fm3bl1BQYF2h3hBQcG2bduOHz/e4aWFqMZ97C6XKyEhQQhRUlLSpUuXHj16CCFsNtvFF1+8bt063UoMVQ6HY/PmzVqf1fbt210u18CBA4UQFotl+PDhvEMb271791tvvTVnzpzGjevWrbv88su15YKCAo5YYy6Xq6io6Nlnn921a5e/8YsvvgjiEQu/IDxy5MhZZ53lX83IyDh8+LCO9YSXsrIy/9FLTU01m81lZWX6lhSC3nrrrX379k2bNk3wejudDRs29OjRIyUlRVEU7cP9yJEjaWlp/ul4OGKNKYoyffr0J554ovG3Lo/HU1lZ6X+Zpaen8670M5vNw4YNq6ys3LRpU35+/uLFi4UQx48fdzqdQTxi4ReEZrPZ5/P5Vz0eD9OGBc5sNnu9Xm3Z5/P5fD4u/WiiuLj4zjvvfPvtt5OSkkRzrzeOWGMDBgzYsmVLcXFxTU3NvffeK05+jQneoSd78skn+/TpM3LkyMaNJpNJlmX/QfN6vbzG/HJycpYvX/6Pf/zj5ZdffvbZZ++55x4hhMViEUIE8YiFXxB27ty5rKxMURRt9dChQ506ddK3pDCSlZXl/3quLXD0GisuLr722muXLFkyZMgQraVz586HDh3yb3Do0KHOnTvrVF0oMpvNKSkp/fv3v//++7Vrrzp37lxRUeF2u7UNeIc2tnjx4nXr1uXn5+fn55eXl99xxx0vvviiLMuZmZn+lxmvsZYMHjy4rKzM4XBYrVabzdb4iLXzNRZ+QXjxxRdbLJa1a9cKIQ4fPrxx48bRo0frXVTYKCwsXL58udPpFEIsXbp0xIgRXNrut3HjxnHjxr3yyivDhg3zNw4ZMsRut2/YsEEIsW/fvm3btl155ZX61RhaGhoa/MubN2/WLmE/55xzsrKyVqxYIYSoqqpas2ZNYWGhbiWGmNdee23RokUvvPDCCy+8kJycfNddd40ZM0YIUVhYuGTJEiGEqqpLly7liPlpH1aa5cuX5+bmajNojx07VjtiiqK89957Y8eObdfDtOdKG728/PLL6enpN910U05Ozt133613OSHqoYce6t+/f2JiYnZ2dv/+/T/66CNVVb1eb0FBQf/+/adNm5aamvrVV1/pXWYIyc3N1U5uNP/4xz+09meffTYjI+Pmm2/u1q2b/9I+qKp6+eWXjxgxYtq0aUOGDMnIyNiwYYPWvnTp0rS0tBtvvLFXr1633nqrvkWGrMZXjZaWlnbu3Pnaa68dOXLkeeedV11drW9toeP+++/Pz8+fOnXq8OHDU1JSPv30U619586d6enpkyZNGjp0aH5+vt1ub8+jhOuvT+zevfubb77Jzc29+OKL9a4lRO3fv7+ystK/ql3RIITw+Xyff/55ZWXlsGHD/PeMQwjx/fff+zv0hBBpaWndunXTlnfs2PHdd9+dffbZ/fv316m6UGS320tKSo4ePZqenj5w4MDGP9xYWlpaUlLStWvXwYMH61hhKNu2bVv37t1tNpu2WlNT89lnn0VHR48aNYpJQvwcDsemTZsOHTqUmpp68cUXayP3mqqqqtWrVyckJIwYMaKdY4ThGoQAAARF+I0RAgAQRAQhAMDQCEIAgKERhAAAQyMIAQCGRhACAAyNIAQAGBpBCAAwNIIQAGBoBCEAwNAIQgCAof1/ye1XNulvviUAAAAASUVORK5CYII=", + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/html": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plot(results.free_energy)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Free energy is converging to a stable value, indicating that the inference algorithm is converging. Let's visualize the posterior distribution and how it compares to the true parameters." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dZ2AU1d4G8LMlu8lms9lN7703QgokoUkNggRQOioElKIX9Cpipagor15Brwiv0hUQBaIg0gQjgVAvBAKEkEYS0nuyu9le3g9z32VNIEDaJpnn92l39szMf8fIszNz5hyGXq8nAAAAdMU0dQEAAACmhCAEAABaQxACAACtIQgBAIDWEIQAAEBrCEIAAKA1BCEAANAaghAAAGgNQQgAALSGIAQAAFpDEPZQarX67bffNnUVPYharTZ1CT0IjoYxHA1jOBrtgCDsoaRS6ZYtW0xdRQ+iUChMXUIPgqNhDEfDmFKpNHUJvQ+CEAAAaA1BCAAAtIYgBAAAWkMQAgAArSEIAQCA1hCEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0xjZ1AQDQ+b7++uv09HRTV9FNtFoti8UydRWdb9SoUQsWLDB1FbSAIATog44ePRoTE9OvXz9TFwLtdPny5RMnTiAIuweCEKBvGjJkSGJioqmrgHZisVh37941dRV0gXuEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0hl6jALSwaPGin3/6uXv2xWKztm3dNnHixO7ZHUAHIQgBaOFeyb3lny6PeyquG/b19cdfl5eXGy/R6XRpaWlXr16tra194403HBwcHn9rWq32p59+OnHihEwm8/X1feGFF8LCwp6ontLS0sbGxiddC+gDQQhAFxaWFgJrQTfsiMPhtFjS1NT0xhtvxMTEbN269cUXX3z8INRoNJMmTcrPz//nP//p6up669atsWPH3rhxw8bG5vHrOXTo0Llz53788cfHXwVoBUEIAF1OJBJdu3aNELJt27YWH92+ffvatWvNzc3U21GjRvn4+Bg+/e677zIyMnJycqysrAghzzzzzMKFCy0tLQkhUqn0xx9/LCsri42NfeaZZ6j2O3bsmDBhwo8//tjY2Dh58uTw8PDa2tpz587l5+dv3rxZKBROmzaNEJKamnr27Fk7O7sZM2bY2toSQtLS0kQiUUREBCGkvLz8woULzz33nFar3bZt26RJk3bt2sXn8xcuXJiWlvbXX3/p9frg4ODp06czGIyuPnTQDdBZBgBMZuPGjSNHjszOzj59+vTChQv37dtXXV1t3ODAgQPz5s2jUpAiEok4HI5cLh8wYMCZM2ecnZ1XrFjx+uuvU58uXbp09uzZcrlco9EMHjy4tLS09U4///zzhQsX2tvbZ2VlRUZG1tTUEEJ27dr1xx9/UA1yc3M//vhjQoharV64cOFzzz2nUql4PN7Ro0fnz5/v5OTk6el56tQpnU7XRYcFuhnOCAHANDQazXvvvXf48OGhQ4dSbwMDA+Pi/nYX8969ey+88ELrdb///nuRSLR7925CyIQJE/z8/JYtW+bm5kYIWbRo0eTJkwkh169fP378+EsvvTRo0CBCCDVup1wu//DDDy9cuECd/FVVVW3YsOGjjz5qo84VK1aMGTOGELJ69epRo0YtWrSIEDJv3rzOOQrQA+CMEABMo6SkRC6XDx48mHo7YsSIK1eutGhjbm4ukUhar5uVlWVY0dXV1cPDIzs7m3obGRlpWE6d7RkrLCxksVhUChJChgwZcuvWrbbrjI2NpV7MmDHj+PHjgYGBb7755o0bNx7nO0KvgCAEANPg8XharVatVlNvFQoFj8dr0SY2NvaB80nxeDy5XG54K5fLDesapmRiMBh6vb71ikql0nBVUyaTUSuyWCytVmtYaLwKl8ulXgQFBRUUFHz//fcMBmPQoEE5OTlP9oWhp0IQAoBpODo6BgYG7tmzhxCiUql++umnp556qkWbN9544/fff6cugRJCdDrdhg0bmpqahg0bdujQIbFYTAj566+/mpubDSd5rQmFwrq6Ouq1u7u7q6vrzz//TAhRKBT79u2jdurl5ZWRkUG1SUlJeeB2ampqWCxWXFzcF1984e/vn5eX1/4vDz0JghAAukNcXJyvr69er09MTPT19W1sbCSEbN++ffXq1WPHjg0PD3d2dl64cGGLtSIiIn755Zf333/f19d32LBhTk5Op06d4nA448aNS0pKCg8PHzdu3MyZM7du3WrcoaaFxMTE4uLi4ODg5557jsVi7dix46233ho7dmxYWFhgYODcuXMJIXPnzr148eKQIUNiYmIMJ6ktfPbZZ8HBwc8++2xCQoK1tfWIESM67eiAST3g0gH0BA0NDb6+vvX19aYupKeQSCRt/EtHN488GmPHjv3nP/9pPB/huGfGqRlqN0+3rq+OXDpzacniJYsXLzZe2NTUZNzNUigUUs8eKJXK27dvi0QiLy+vh21Qp9NlZ2dLJBI/Pz87OzvD8oqKirKysqCgID6fTy1pbGwUCARMJpMQIpPJmEymubk59ZFUKtXpdAKBgBAil8uzs7Pt7Ow8PDwMW5PJZFlZWa6urg4ODjKZTCAQ6PX6xsZGQ6mEkNLS0pKSEjs7O39//w4do0f55Zdf9uzZ87Bz0zZIpVLD0YDHhF6jALSw7I1l//nPf7pnX6Evho4fP77FQmtr6wc25nK5/fv3b3uDTCYzNDS09XJnZ2dnZ2fjJUKh0PC6xR1H43iwsLCIiopqsTUej2foF0PlJYPBEIlExm3c3NyovqnQlyAIAWhhxIgRuJQH8EAIQgBaeOmlBQcOHOiefTEYZMuWzVOmTOme3QF0EIIQgBbKyyvefntTXFzio5t22NdfL2v9AN/27dtPnTpVW1vr4eGxePHi6Ojox9/gf/7znx07duTn5zs4OIwYMeL5559vPZxpG/R6/Z9//jl8+HDDkxUAxhCEAHRhYcEXCESPbtdhHA639cLGxsZJkybZ2tpeuHBh2LBh165de8z+Jj/88MMrr7yyfPny5557rra2lhqG7Z133nn8ejQazejRo8ViMfpbwQMhCAGgO7zxxhvUi5EjR/76668XLlyggrCmpiY1NbWkpIR6nt3Ly2v69OmGterq6hYvXvz9998bLrROnz7d0Jv6t99+O3funJ2d3Ysvvujo6EgISU1NZbPZDQ0Np0+fDgkJmT9/PpPJ3LVrFyHkyy+/5HK5s2bNcnd3Lyws/Omnn5qamhITE4cPH04Iqaio+P33319++WVqy1u2bBk/fryLi8uRI0ccHBxyc3MvXbr0wQcfqFSq77//vqqqysnJadasWW30dIVeBM8RAkA3kclkdXV1R48eLSsrowZIKyoqCg0NTU1NVavVX3/99Y4dOyoqKoxXOXHihEAgeO6554wXUnMwrVq16r333gsODq6srOzfvz81WvfRo0eTk5NPnjzZr1+/jRs3UmNnUyeCQqFQJBKx2eycnJwBAwYolUpfX9/k5ORvv/2WEFJSUvLFF18YdrFu3bqSkhJCyL59+6ZPn37p0qWwsDC1Wp2QkKDX64cOHcrhcO7cudPVRwy6B84IAaCbrF69esuWLRKJ5F//+hc119L69evHjx//3XffEUISEhKmTZu2ZMkS41VKSkq8vLxaz3YklUo/++yz69evBwUFEULu3bv3zTffUGNnBwcHf/PNN4QQR0fHlStXrlq1atKkSYSQ5ORkKhFXrFgxe/bs1atXE0J8fX1nzJhBjaP9MAMGDPj6668JIeXl5Y2Nja+88soTzYYIPR/OCAGgm3z++ecNDQ3Z2dlffvkl1YU1Ozt7yJAh1KeDBg2qq6trMXESj8drampqvam7d+/yeDwqBQkh8fHxt2/fpl6Hh4dTL5ydnWtra1uve+fOHcMcF3FxcTU1NS3mfmrB8MShi4vLrFmzPDw8nnnmma1btz5sABrodRCEANCt/P39x40bd+bMGUIIn8+XSqXUcrlcrtVqW/RnGThwYEFBQetpBfl8vkwmMwxVYzycyiO7hhrvtLm5mcFg8Pl8Nput0WgMbQwTBRNCjHuofvvtt0VFRdOmTdu4ceMTddiBngxBCABdTqVSFRUVUa9ramr+/PNP6rxt5MiRe/bsUSqVhJBt27b179+/xVXHAQMGDBs2bP78+YZRsy9dupSSkuLp6enm5kYN2C2VSn/66aeRI0c+bO9mZmZWVlaG074RI0bs2rWLOp/bunVrfHw8j8fz8PAoLy8vLy8nhJw+fbqsrKz1dqRSaVNTE9U356WXXsrPz+/ocYGeAfcIAaDLKRSKAQMGCIVCCwuLwsLCmTNnUhPbLliwID09PSgoyMXFpaKi4pdffmm97r59+xYvXuzl5eXv719XV2dpabl582YWi/X999/PnDlz8+bNRUVFiYmJs2fPbqOAZcuWDR48mMfj7d27d8mSJRcuXAgODra3t6+rq6N2amdnt2TJkv79+wcEBHh6enp7e7feSElJybBhw4KCgiwsLHJzcw1zYkBvh0G3eygMut0CBt021p5Bt8dNEIs1Tk4ebazVWa5fT1+2bEmLHig6ne7u3bsKhcLDw4MaydOgpqaG+oNv46pmU1NTYWGhg4ODi4uLYaFGo6GesjecR1KTFFpYWBBCtFqtVCo1HuO0sbGRugpKCKmqqhKLxb6+vtQI3ZTS0lKZTBYQECAWi3k8HpvNbm5uZrPZhikJVSpVQUEBIcTPz8/MzKydB+gxYNDt7oQzQgBaWLHivZs3b3bPvsaMiW496DaTyfTz83tge3t7e3t7+7a3aW1tbZh63oDNZhv6y1CoCKSwWKwWI30bD8nt6OhIPXpozDCgtiGqLS0tjRtwOJzg4OC2S4VeB0EIQAvx8fHx8fGmrgKgJ0JnGQBoqba29p133vP0DDYz45qb86Ojh2zfvp0a+QWg70EQAsDfXL9+PTAw4t///s+9e2s0mutK5dmMjOmvvrpm6NAxxg8VPBGpVDpt2rT333/fsEShUMycOfPNN9/spKrb48qVKzk5OYa3AoHg3r17HdngsWPHcF+/N0IQAsB9Uql09OhnGhsXKBQnCXmOkGBC+hPyD4UiMyNDPX/+4kdv4kFUKtX+/fu/++47w2MJv/322/Hjx0+cONF5tT+x77777tChQ4a327dvt7Oz68gGlyxZQnWlgd4F9wgB4L6tW7c2N7vodKtafWKlUOw6cCAgP3/lw/q8PNKUKVN27dpFPYe+Y8eOGTNmnD171vDpgQMH0tPTBQLB/PnzPT09CSH5+fkHDhwoKipycHAwLLx582ZeXp5IJNq/f7+Tk9Prr7/eog8qIUStVu/cuTMzM9PNzW3RokVUH5mcnJwdO3Y0NDQ4OzsnJyc3NTVlZmaWlZXp9frg4OCkpKTi4mLqmfqUlBQfH5+zZ8/evn07MTFx8uTJ+/fvT01NjYyMXLBgAYPB0Ol0+/fvv3jxolarTUhImD59OoPBOHz4cGNj4w8//JCampqYmBgZGSmVSrdu3ZqXlxcYGLhw4UJD11PoaXBGCAD3/frrCbl8JiEtx/YkhBDiaW6ecPLkyXZvfN68eTt37tTr9eXl5VlZWaNHjzZ8tGTJkm+//Xbo0KHW1tYJCQnUg+3nz5/n8Xjjx483NzcfOHAgddXxypUrS5cu3bFjx1NPPZWRkfHCCy+02Iter3/66afPnDlDTb2UkJCgVColEsmwYcNcXV0nT55sZ2d37949KysrKysrW1tbHx8fqvvohx9+2NjYSAj54Ycfpk+frlQqBwwY8PLLL7/wwgsXL14cNmzY119/vWnTJkKISqVKS0uLj48fOnToV199tWbNGkKIg4MDm812dXX18fERCAQymWzAgAGlpaVjxozJzMycMGFCu48bdDWcEQLAfeXlFYQ89FlDlcqLiqj2CQgIEIlEly5dOn369OzZsw1PDRYWFu7Zs6ekpIR6VqG0tHTr1q0rV6588cUXCSEymWzw4MGXLl06evTo888/Twjh8Xg7d+5kMpmDBg3y8/PT6XTGzwIeO3astrb25MmTDAZj4sSJV65c+f333yMiInQ63cyZM40vfvr4+Pj7+0+dOrV1qZMmTXrrrbcIIRkZGXfu3KEmcpLJZAcPHnz11VfNzc03bdqk0+mampqcnJzmz5+/YsWKgQMH8vn8kSNHxsbGEkI2bNgQHBxMzWiRlJQUEBCQmZnZr1+/dh896DoIQgC4j8+3IkT8sE/Z7CYrqw49RTdnzpydO3eeOXMmJSUlNzeXWnjr1i2lUjls2DDqbXV19VNPPUUIOXHixJIlSywtLe3s7LKzswcNGkQ1CAoKopLP0dFRoVA0NzcbDy9w48aNkpISKo0IIffu3cvPz3/22WefffZZT0/PhISEyZMnv/zyy20/Dm94PNHOzs7w4KCdnR11VqpUKufPn5+enu7m5sZkMh84HtuNGzcuXLgQExNj+FL5+fkIwp4JQQgA9w0dGpuVlapWJz/oQ5VOlx4bu+RBHz2uGTNmvP322yEhIcHBwYYgFAgErq6uV65cadH45Zdf3rNnDzU9xZQpUwxDbBuf/7UmEAiGDBly8ODBFsu//fbbzz777OjRo59//nlZWdknn3zSxkaMd9F6d7t3766uri4sLGQwGDdu3HjgA5oCgWDWrFnGcxxCj4V7hABw38KFL+n1Bwi53PojJvNTV1eHoUOHdmT7QqHw119/3bx5s/HCmJgYuVy+f/9+6q1MJqMuwIrFYupiaX5+/uP3L01MTExLS7t69Sr1tr6+vqGhob6+vr6+3traeubMmVOmTKEm3bW1tW0xD/BjampqsrS0ZDAYer1+48aNhuW2traGS8cTJ07cs2eP4XmM0tJShULRjn1BN0AQAsB9ISEhn3zyMYeTSMhOQgzz7dWwWK9xuV/u3//DIyc5eqQRI0YYpgykWFpapqSkrF69OiIiIi4uLigoKCMjgxCyatWqMWPGDB06dObMmYYLp4/k6+u7Y8eOZ599NjY2NioqKiYmpry8vLi4OCQkZMCAAQkJCbt27Vq+fDkhZM6cOSdPnvT09Hz99def6CvMnj07Ozu7f//+4eHhPB7PsPzNN9987bXXvLy89u7dO3To0FWrVg0cODAhISE8PHzkyJHUOKjQA2HQ7R4Kg263gEG3jbVj0O0n8tNPP/3zn+/W1dVxOP6EqGWy23Fxw7Zu3RASEtK+DT6m6upqhULh4uJCjYtNCGloaGhoaPD29m49Sf0jlZaWMplMJycn6tqmVqulTtdcXFw6Hucajaa0tFQoFBqPX9qaTqcrKSkxNzdvPa5p2zDodnfCPUIAaGnGjBlTp069fv16Xl6eubl5ZGSkl5dXN+zXwcGhxRKRSCQSidq3NcMI2hQWi+Xu7t7Oylphs9mPc0yYTCb1+CP0ZAhCAHgAFosVHR0dHh5uPD87QJ+Ee4QA0FJ5eflrS193c3bncrksFisyvP8333xDTekO0PfgjBAA/ubSpUvjx44Pcej3evz7XjY+Kq3qRlnGuo+/3PPDjydOHW89ntnjkEqlc+bMab18zZo1HZ/eLzk5mRoN3MnJafr06YbHDbuTRCL566+/kpKSun/X0HEIQgC4r6mpacL4pGnhL86KmWdY6CHyGh00/t0jS16a9/K+Az+3Y7NcLpcaJoYQMmnSpA0bNlC36560C8kDHTp0aOXKlcHBwRkZGSNGjDhx4gT1PH53qqiomD9/fk1NTTfvFzoFghAA7tuyZYsz3804BSlcNvfdkWtm7Hw6JycnMDDwSTdrZmY2ceJEQgjVTX348OGhoaFpaWmnT59WKpVHjx5dsmRJfn6+r6/vwIEDCSF37tz566+/Fi9eTAhRqVTfffddRkaGq6vr0qVLW3eoIYTExMQMHjw4MTHx2rVrKSkpYWFhO3fuvHnzpqWl5ZQpU0aMGEEIqa2t3bRp08SJE7/99tvQ0FBq4NMrV66w2ezRo0dTA63p9fp33313zpw53377rVQqXbJkiY+Pz7p164qLi2fOnGnohZuXl7d9+/aqqqqhQ4fOmTOHwWBs2rSpubmZGk/8nXfeEQqFN2/e3LVrV11d3ahRo2bOnEkIKS4uTklJiYuL27Fjx9ChQ6dOnbpp06Zr165xudwRI0bMmjXrSY8qdBbcIwSA+04c/eMp79EP/Mie79DPI+rUqVOdta+LFy++8sorZ8+enTJliq2t7W+//Xbt2jXqo/z8/D179hBC9Hp9YmJiZmbmjBkzLC0t4+Li2p4TUSwWm5ubFxYWEkJmzZo1aNCguXPn/vnnn4SQhoaGtWvXLl++fMSIEREREbW1tdXV1VOnTh07duyaNWu2bNlC7e6zzz5bunTpsGHDvL29ExMT582b5+3tPWzYsKlTp1JTLN24cWP48OGenp5Tp07du3cvFX5BQUFsNjs6Ojo6OprD4Zw/f378+PFBQUHPPffcpk2b1q5dSwgpLS396KOP1q5dO3bs2ICAgA8++CA9PX3evHmTJk2qra3trKMK7YAzQgC4r7Ky8infh16udLB0at9QLA/j7u5OzefwMCdOnKivr09NTWUwGImJiZcuXTp48ODs2bNbNNuzZ096enpmZua5c+fWr18fEhISGxvb1NRUV1eXnJy8b9++kSNHEkIUCsWWLVs8PP47qvjq1atlMlllZeXSpUv37t378ssvU8vXrFkzcODAZ599dsuWLSNHjpw7dy4h5NChQ2lpab6+vmvWrFm2bNmiRYsIIbGxsZ6enmvXrh0xYgSXyzWM3/3hhx9+9NFH1IqBgYGDBg169913CSFKpXLXrl3Uo4effvppfHz8sGHD2h4xDroBghAA7hMIBFKV9GGfSlUSa2vrTtxdWFhY2w2ysrIKCgoMMyDW1dVFR0c/sCWDwRg1atSXX37p5OSUk5MzY8YMJpNpY2NTUVHh4+NDtREIBIYUrKmpmTp1am1trbOzc1NTk/GJZkBAAPXCxsbG8NrW1rahoYEq6ezZsxs2bKCWKxSK1j8OsrKyPvjgg48//pgQotfrq6qqqO27ubkZHsB/9913X3zxxQ0bNjz99NNLly6NiIh41NGCroIgBID7EgbHXz12cXzo5NYfqTTKzNKMjwau7MTdGU8BYWZmZnhCQyz+7wwYIpFoyJAhx44da3s7s2fPHjx4sOHtqlWrZs+evWzZMkLIl19+mZqa2np3X375ZUhICHU+eujQIeqMjWI8io3xa+oGp1AoXLVq1bRp04wLyMvLM34rFArXr18/ZsyYNr5vXFxcbm5udnb27t27hw4dWlZWRg2sCt0Pp+QAcN/LC14+k/9nZllG64+2X/pfD0/3rns4wd/fPy0tTa/XKxSKHTt2UAvHjBlz+fLltLQ06m1lZeXj9MyUy+XU9cbGxsZt27Y9rA0Vckql0njg7EeaMmXK+vXrqSl8CSG3bt0ihNjZ2YnFYsPCKVOmfP7551KplBCi1+upNi1kZWXp9frg4OBly5bJZDKZTPb4NUDnQhACwH0BAQHr1q17+/Crv2TuVWqU1MIaadW/Uj/8PfvA7r27OzhKJ4PBMDc3p1LK3NzceMTqRYsW5ebm+vn5RUZGRkVFUYOpurm5/fzzzwsWLKBmbho6dGhlZWWLbQqFQsPYpJTly5d/9tln0dHRCQkJo0aNosbeZLFYxuOCvvrqq0ePHo2IiIiIiIiKijI8HykSiQxngQKBwHAOx+PxzM3NCSFLly596qmngoKCoqKi3NzcPvzwQ2qtf/7zn5GRkb6+vuXl5e+++25ISIifnx/VZv369YQQNptt/BTmmjVr3NzcBg4c2K9fv48++sje3r4jBxY6RA89Un19vUgkMnUVPYhYLDZ1CT3II49GYmLi8ePH2739Q4cOBfgGcM24vs7+bvYeTAZz3Njx+fn57d7g46uoqFCpVK2X19XV1dfXP/521Gp1WVmZWq1uo41Wqy0rK1MqlU9c5f+vq1Ao2mij0Wja3r5MJistLX3g901JSXn22WfbUZhEImnHWjSHe4QA0FJSUlJSUhLVUYXL5UZGRnbKk++Pw8nJ6YHLbWxsnmg7bDbbxcWl7TZMJvORbTqyLovFaruNhYWFq6tr+wqAToQgBIAHCw0NDQ0NNXUVAF0O9wgBAIDWEIQAAEBrCEIAAKA1BCEAANAaghAAAGgNvUYB+iAul5ucnGxhYWHqQrqDXq83Hgitb2hubqZmj4JugCAE6IP27NlTXV1t6iq6SXNzc58cpdPZ2dnUJdAFghCgD+Lz+dS4YnQgkUio8dgA2gf3CAEAgNYQhAAAQGsIQgAAoDUEIQAA0BqCEAAAaA1BCAAAtIbHJ9rp8OHDx48fr6io8PX1/cc//uHp6dm6jVwu37RpU3p6upWV1axZs8aOHUsI2bhxY0lJCdXAzc3tH//4R7fWDQAAf4czwnZauXKln5/f3LlzxWJxXFxcbW1tiwZarXbcuHGnTp2aNWvWmDFjDA1++OEHuVzu4+Pj4+ODOTkBAEwOZ4TtdO3aNerFhAkTTp06dfr06SlTphg32LNnT2Vl5a1bt1gsVot1n376aersEAAATA5B2FFyuby2ttbFxaXF8rS0tIkTJ27duvXWrVuRkZFz5841JOK2bdt+//330NDQ5ORkc3Pzbi8ZAADuY+j1elPX0LvNnz+/srLyyJEjLZaPHDny5s2bc+bM6d+//1dffRUWFrZ9+3ZCyMqVKx0cHNhs9vfff6/T6dLT083MzFpvtqamxtXVdfDgwYQQLpe7d+9eNpvWv1qkUil9xgx7JBwNYzgaxtRqtUgkMnUVvQyCsENWrFhx+PDhv/76q/Vf3jPPPCOTyVJTUwkhubm5QUFBjY2NAoHA0EAmk/n4+Gzfvn3cuHGtt9zQ0ODp6fnLL78QQiwsLAYNGtSV36MXwHiSxnA0jOFoGMPPgnag9UlGB33yySe//vrrA1OQEOLh4dHc3Gx4rdfra2pqjIOQx+N5enpWVlY+bPtsNnvUqFGdXjYAABhDr9F2Wr9+/Q8//HDy5El7e3vDQplMtm3bNplMRgiZMWPGmTNnqCw8cuSIg4ODp6dnc3NzfX091fjSpUs3btyIjY01Sf0AAEDBGWF7SKXSN998k8fjhbIkNmgAACAASURBVIaGUks+/fTTRYsW1dfXv/TSS4mJiTweb+jQoU8//XRISEhAQMDNmzd37NjBZrOLiooiIyMDAgJYLFZubu7atWvDw8NN+10AAGgO9wjbQ6/XNzY2Gi/h8XhcLlen09XV1dna2jKZ/z3VLioqqqurCwoKMkwcKpVKc3NzCSH+/v5t3NhoaGjw9fU1nD4C7gMZw9EwhqNhDPcI2wFnhO3BYDAeeF+QyWQaXyklhHh5eXl5eRkv4fP5UVFRXVoeAAA8PtwjBAAAWkMQAgAArSEIAQCA1hCEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0hiAEAABaQxACAACtIQgBAIDWEIQAAEBrCEIAAKA1BCEAANAaghAAAGgNQQgAALSGIAQAAFpDEAIAAK0hCAEAgNYQhAAAQGsIQgAAoDUEIQAA0BqCEAAAaA1BCAAAtIYgBAAAWkMQAgAArSEIAQCA1hCEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0hiAEAABaQxACAACtIQgBAIDWEIQAAEBrCEIAAKA1BCEAANAaghAAAGgNQQgAALSGIAQAAFpDEAIAAK0hCAEAgNYQhAAAQGsIQgAAoDUEIQAA0BqCEAAAaA1BCAAAtIYgBAAAWkMQAgAArbFNXQBAL6NWq6VSqVQqVSgUzc1KuVyp1erVag2DQczMzFgshoUFh8fjmpub83g8Pp/P5XJNXTIAtAVBCPBoCoWiurq6qqq+rKxeKlWx2Xw224rFMmez+RyOLYPBYLHYhBClUq3X62tqVBqNUqut12juaTTNHI7ewUHo5CS0sRHZ2tqyWCxTfxsA+BsEIcBDqdXqkpLSgoLy6upmDseex7O1tfVzd+c/0UZUKqVU2pid3ahSFWg0Gc7O1h4eDi4uzjwer4vKBoAngiAEeACJRJKTU5CXV8VmO4pEAYGBdgwGo32b4nC4NjaONjaOhBCdTtvYWHf9etXFi+fs7Mz9/V3d3d04HE6n1g4ATwZBCPA3Mpns5s07BQV1fL63t3eImVlnphSTybKxcbCxcdDrw8Ti+oyM0kuX/vL2tgsI8LK1te3EHQHA40MQAvyXTqfLzc2/dq3I0tLHzy+Cuu3XRRgMhrW1rbW1rVYbWl1dduzYTTs7Vni4j4uLS7tPPQGgfRCEAIQQIpFIzp+/1tjI8/IayuGYd9t+WSy2s7Ons7NnfX1VWlqBtXVe//4Bzs7OiEOAboMgBCClpWVnz2YJhcF+fu6mqoG6j9jYWHv6dI6NTX5MTIidnZ2pigGgFQQh0N2dO7lXrpR6eibweE/WHbQrCIV2QqFdXV3lsWM3fH0FkZEh6FwK0NUwsgzQWmbmratXq/38BveEFDSwtXUKCnqqpkZ46FB6fn6BXq83dUUAfRmCEOjr5s2srKwmf/+4zu0a2imYTKabm5+Hx+DLl2v++CNdIpGYuiKAPgtBCDRVUHD3xo06P7+BXdo7tIPMzXmBgXEajedvv10oLCwydTkAfROCEOiopqbm4sW73t6xPTkFDZycPDw9B50/X3L+/BW1Wm3qcgD6GgQh0I5SqUxLu+7iEsXlWpi6lsdlYWEZGDi4osLi+HFcJgXoZAhCoJ0rV26wWO4CgY2pC3kyDAbDyyuUzQ44cuRCVVWVqcsB6Dt6wXUhgE5UVVVVVNQcGBhj6kLaycHBlcfj//XXWSaT7ePjZepyAPoCnBECjeh0usuXbzs6hvbqcVv4fGtv7/jz54tv3co2dS0AfQGCEGiktLS0uZknEtmbupCO4nIt/P0TMjPrr1+/iacMAToIQQh0odfrr1/Pd3T0N3UhnYPNNvP3H3j7tvTq1UxkIUBHIAiBLqqrq2Uy817XR6YNLBbb339Abq7i+vWbpq4FoBdDEAJd5Offs7b2MHUVnYzJZPn7x2ZnS2/ezDJ1LQC9FYIQaEGlUhUX19vZOZu6kM7HZLJ8fWMzM2sLCu6auhaAXglBCLRQXV1tZmbHZLJMXUiXYLPNfH0HXrhQWFlZaepaAHofBCHQQklJlaWlg6mr6EIcjrm7e8zp0zfFYrGpawHoZRCEQAtlZfVCYR+f55bPtxaJwk6fxnikAE8GQQh9n0wmUyoZvWhk0Xazs3PWaJwuX76OByoAHh+CEPq+pqYmDkdo6iq6iYdHcFGRCnM2ATw+BCH0fc3NzSyWpamr6CYMBsPLq//Fi3mYpALgMSEIoe8Ti2VcLl2CkBBibs4TCoMuXMAFUoDHgtkn2kMsFu/cufPMmTNyuTwqKuqNN94QiUStmxUUFPzrX/8qKChwdHR8++23w8PDCSFVVVWrV6/Oz8/v16/fypUrBQJBt5dPOxKJgsPhduceFQqFWCxWKBRqjZoQYsY2Mzc35/P5PB6vg1vW6bSlpQVlZUWBgf1sbBwf1szR0SM3tyI//66/v28H9wjQ5yEI2+PWrVvnz5+fNm2aSCRat27dhAkT0tPTW7QpLCwcNGjQ/Pnzp06dWlFR0dzcTC1PSkoKDw9fvXr1unXrkpOTU1JSur182lEq1WZmnG7YkU6nq6ysLCwpbJI1cSw5LA6LyWRSy7UqrUauMWOauTq6uji5tO8H0O2s//zrk1cdLZxszG23165y8fVb9u4GK6sH3/708Ai/ejXd1dW54+kL0LcxcPGkg8rKytzc3Kqrq+3t/zanwQsvvGBtbf3NN98YL7x48eK4ceOqqqrMzMzq6+udnZ1zc3M9PT1bb7ahocHX17e+vr5rq+89JBKJlZVV+9Y9fPi0UBjN47Vz9cckFoszszLlDLmNg42l1YOvxKqUqqb6JlmjTGghDPANsLF5goFPq6vLXl849pMx65wELtSS1LwTp2pSv9x45GGrlJXl29k1JiT01skXH1NH/jb6HqlUyufzTV1FL4N7hB2Vm5srEAhaXxpNT0+Pi4t75513XnrppV9++YVaePXq1djYWDMzM0KIjY1NQEDA9evXu7timuraCQirq6vPZ5zn2HLcfd0floKEEA6XY+9s7xHkQazJpVuXMjIzFArFY+7i0C9bZ4S/YEhBQsgI/0QiUd+7l/uwVZydfQoLxfg5BdA2XBrtELFY/Morr6xZs4bN/tuR1Ol0paWlq1at+uCDD/r16/faa69VV1cvWrSoqqrKODJtbGweNiaWRqORSCRRUVGEEIFAcPDgwRa7oJvm5uZ2z6Yrlyu4XFnX/exramq6fPOyg5cDx4LzmMHGteA6ejnWVdedOnMqLCDMweHRo94U5mUN8Xq+xcIA28CcnBt2dm4PW8vCwiM9PWPEiLjHqaqX6sjfRt+D4RTagdb/tnZQc3Pz+PHjR4wYsWTJkhYfMZlMc3PzefPmJScnE0IUCsW///3vRYsW8fl8438oZTLZwy7psNlsHo+3ZcsWQoi5ublQSJfH4B5Gr9e3+4IPj2dhYWHRRbfKtFptdkG2i58L3+qJy3PzdFM4KG4X3dZoNb4+vm3/ay4U2TXKG1tcRm9UNtrZObbx1Tw9/XJyKpqbmx0dH9qzprfryN9G3yOVSk1dQu+DS6PtJJfLk5KSAgICNmzY8MAGnp6eTk5O1GtnZ+eGhgZq4d27/50iQKfTFRUVPfAGIYXFYkVHR0dHR4eGhnZ2+fRibm6m0XTVz+R7Jff05vp2pCDF3MLcM8Azvyo/+0522zfsh42a/NudA8ZLmhSNNyuvh4UNbHsXdnb+mZkPvXwKAAjC9lCpVFOmTLGzs9u8eTPVLZBy586dnTt3Uq9nz5596NAhnU5HCPnll1/i4+MJIePGjSsqKjp//jy10NLSMi6uL1+z6iG4XDO1WtUVW9br9XdL7to62nZkI0wW08PX417Dvdy8tuJqYNxo2wCvFSffyizLKG2892fuiWVHXl361heP7BBrY+NYU6OrqanpSJEAfRgujbbHH3/8cfToUWtra0NP0bS0tPDw8CtXrnz88cdz584lhCxZsiQ1NTUoKIjL5XK53F9//ZUQYmVltXHjxgkTJvj4+BQXF+/atYvF6psTA/UofL55Y+Pj9kl5IhKJRMPQcM07+pAig8nw8PW4m3PX0tLSzfWhN/zefOfr69fSTx37uSav3C8w/N9LT9jaOj164wyGSORz505hi47NAEDB4xOdSafTaTQaDuf+L/SioiIGg9Hi+qdEIikuLvbx8Wnj1g4en2ihI13kCwoKrl9XenmFdG5JhJDS0tLc6lxn986Z71etUpfllQ2OHvzIbyqTyZ7olqdOp8vN/XPSpARLyz44wg4enzCGxyfaAZdGOxOTyTROQUKIl5dX67uAVlZWYWFheMy52/B4PJ1O1hVblsllbE6nXVYx45gJXYQ3sm9QV9Q7EZPJtLBwKyy817mbBegbEITQ9wkEArW6S6arVWvUbNZjBWHx3Tunft+XeuRAWXF+G82EImGztrmioqKTCrzP3t49N7cMV4AAWsM9Quj7LC0tmUyVWq3q9IHWWCyWTvOIszeNWvXvNR8W3NE2S8cTorO0+iY00m7x8neZzAffHrZ3tb9z946zs7NxP6yO4/H4SiW3rq7Ozq6PT1AM8KRwRgi04OBgLZU2dfpmzTnmj3x++cet32ZdHyIRH9LpFuh0iyRNJzIu+h388YeHbtPCXGum7YpOnpaWLqWlnX+uCdDbIQiBFlxdbcXi2k7frKWlpVahbbvN5bNn1aqlxktUynfPnDzZxioCG8G98s6/nycSOd69W9XpmwXo7RCEQAv29nZKZecHobW1tUqmauPGm16v0+u4hLS4CmquUbd1r04gFNQ21Go0mk4q8794PL5CwRKLu+R2KUDvhSAEWhCJREymTKVSdu5mORyOkC9sljY/rAGDwWQwlYS0uI+obHvgWAaDwTZnd0Vicbn2tbWd/4MAoFdDEAItMBgMb2+HuroHD3HeEZ4unk21bd19jB002Mzsb+Pwcbn/M3jUqLY3yzJnyeXyTqjv73g8UWUlHk4F+BsEIdCFh4dzc3PndxVxdHTUyXUq5UOHcJu9YHFwv1QrwbMMxvdM5jYrwfiI2OzJs+e2vVk2my1XdH4QWlvblJc3dPpmAXo1PD4BdGFvb89g3FCpFByOeSduls1mB3gF5JXnuXk/eGg0thnnrY//pzAv687N62w2Kyj8ZXfvwEdvl0k6/bF6QgiXa6FQ6JVKJZfb0WHhAPoMBCHQBYvFCgx0KSwscXf379wte7h7FJUWSZokVtYPHejL2z/U2/8JZhHR6XSP+aj+kzIz40ulUgQhgAEujQKNeHm5NzeXdProKkwms19Iv7rSuk7s56lVaS0sLDpra8ZYLCuJRNIVWwbopRCEQCPW1tZOTtyu6DIjEokCPQPLCjttDDONTNNFA0mz2RYyWZfMxQHQSyEIgV5CQnzq6gq6YsveXt6uQtfSwtKOZ6FCruCyuF00UwSXay6VIggB7kMQAr04OTlZWambmuo6fcsMBiM0KNSR51h6t1Sv61AW1lXW+Xj4MBiMzqrNmJkZVybr5OcpAXo1BCHQC4PB6N/fv6oqpys2zmQyI8Ii3IRuxbnFbTxQ0TZJk4StZru6uHZubQYsFkujecSwcAC0giAE2nF1dbW2Vjc0dP6o1oQQBoMREhQS4RtRWVBZV/PE551ymby+tD4qPIrFevDcFB3HZCIIAf4GQQi0w2AwoqICKytvd93kfK4urkNih3DknOLcYon4cbtoNjU01RTVxITFCASCLiqMgkkJAYwhCIGOnJycPD3NKyuLu24XPB5vQPSAKP8oVY2q+E5xbXVtGxM2NUuaSwpK1HXq+P7xXT1foE6nY7HwPz7AfXigHmgqMjLk8OGLtrbOHE4XPlru4ODg4ODQ0NBQXllellemZ+jZFmymGZPFZhFCtFqtVqlVy9QCC0GYR5iTk1PnTsb7QDqd1sysq667AvRGCEKgKSsrq4gIt9u3s3x9o7p6XyKRSCQShQSFNDc3S6VShVKhUqkIIRwzDo/Hs7Ky6qJn5x9IrVYJhZxu2x1Az4cgBPoKCQksKTlbW1thZ+fcDbtjMBh8Pp/P53fDvtqgVMqtrLovdwF6PtwqAPpiMplxcRE1Nbc6fZ7CnkyjkfN4nTnsOEBvhyAEWhOJRDExnoWFGV3Xg7Sn0WjEXTR4G0AvhSAEugsM9Hd21peXd8m4az2QWi1BEAIYQxAC3TEYjPj4KI2mqIsese9RZDKJlZUZ5mACMIYgBCDm5ubDh0dXVV2Tyfr4/ESNjXXu7l37nCJAr4MgBCCEEJFINHhwUHHxVbW6nWOE9goyWbWTE4IQ4G8QhAD/5enpER3tXFBwWafrm0NxajRqvb7BwcHB1IUA9CwIQoD7goMDg4Ks8vOv6nQ6U9fS+WprKzw97dhsPD0M8DcIQoC/iYqK8PZmFRT0wSxsairy9/c0dRUAPQ6CEOBvGAzGwIFRnp6Mu3cz+lIWNjbWCoX6rh7RG6A3QhACtMRgMOLioj09Gfn5/9FqNaYup3NUVeVGRPiaugqAnghBCPAA1HlhUJBFXt7FPtCPtL6+SihUu7p21az3AL0aghDgwRgMRv/+Ef372+flpffq5wt1Om1VVVZsbAiDwTB1LQA9EYIQoC0hIYEjRgSWlV2sr682dS3tVFKSExAgsre3N3UhAD0UghDgEVxdXRMTY6TSG/fu5fS6sbkbGmoYjPLIyFBTFwLQcyEIAR5NJBKNGzfE1rYhN/eiSqUwdTmPS6mUV1VdHzYsisPBTLwAD4UgBHgsXC53yJCB0dF2hYVnq6vLTF3Oo6nVqrt3L8fH+9nY2Ji6FoAeDUEI8LgYDEZgoP/48QMIyc/Pv9qTp/PVaNQFBZf793f08fE2dS0APR2CEODJWFtbJyYOCQ+3LCpKq6go6oF3DVUqZV7ehfBwUWhokKlrAegFEIQAT4zJZIaEBE2YkMDnV+TknO1RExk2N4sLCs7FxDhFRKCDDMBjwfC7AO3E5/Ofeiq+srIyIyMrL8/C2TmIz7c2bUlVVSWNjdkjRoS5uLiYthKAXgRBCNAhTk5OTz/tWFx87/r1KxUVfEdHf4HABJ1TVCpFcfFNkUg+YUK8lZVV9xcA0HshCAE6isFgeHl5eni4l5aWZmZmVlSwRSJvOzsXJrM7bj3odNqysgKZrCgqyiswMAbDxwA8KQQhQOdgMpkeHh7u7u41NTU5OUU5ObctLFxsbFwFAlEX7VGjUVdUFEmlRQEBdmFhg3k8XhftCKBvQxACdCYGg+Hg4ODg4BAdLSstLcvNzSwv11lYOFlbO1hb23bK6Zper29qqquszFWrqwIDnQIC4nAtFKAjEIQAXYLH4wUE+AcE+IvF4srKqsLCO3fuSDkcGy7XxspKZGkpYLPNnmiDcnmzRNIgldYolTWWliQmxt/FJZjL5XZR/QD0gSAE6FoCgUAgEAQE+KtUqoaGhurquqqqnLt3xTqdmZkZn8m0YLEsOBwum23GYrENp4wajVqjUatUSq1WptXK1GqxpSXb2VkUGmpnbx+k0WhwFgjQWRCEAN2Ew+E4Ojo6OjqGhxNCiEwmk0qlcrlcJpPLZI1KpVqpVOt0hBDCZBIrKzMu14zH41haing8V4FAYDxeqETSi6eFAuhpEIQApsHj8dC9BaAnwMgyAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0hiAEAABaQxACAACtIQgBAIDWEIQAAEBrCEIAAKA1BCEAANAaghAAAGgNQQgAALSGIAQAAFpDEAIAAK0hCAEAgNYQhAAAQGsIQgAAoDUEIQAA0BqCEAAAaA1BCAAAtIYgBAAAWkMQAgAArSEIAQCA1hCEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0xjZ1Ab3V5s2bT5w4UV5e7uPj89Zbb0VGRrZocPHixfXr1xverly5MiwsjBDy3nvv5efnUwt9fX3Xrl3bbTUDAEBrCMJ2Onjw4Jw5c7y9vY8ePfrUU09lZWW5uroaNygpKcnOzl65ciX11sHBgXrx559/jho1igpOkUjUzWUDAEALCMJ2Onr0KPViwIABe/fuPXfu3LRp01q0cXR0nDp1aut1hwwZMnbs2C4vEQAAHgPuEXZUQ0NDaWmpn59f64/u3LmTlJQ0f/78U6dOGS//4osvkpKSVq5cWV9f311lAgDAg+GMsEN0Ot28efMmT54cFRXV4iMvL6/33nvP29s7Kyvr2Wef3bx584wZMwghs2bNcnd3ZzKZ27ZtS0hIyMjI4PF4rbes0WjEYrG3tzchRCAQnD59ms2m9X8sqVRq6hJ6EBwNYzgaxjQajalL6H0Yer3e1DX0Vnq9fuHChXl5eUePHrWwsGij5fr161NSUs6dO2e8UK1W+/r6btiwYeLEia1XaWho8Pb2zsjIIIRwOBw3N7fOLb7XkUgkVlZWpq6ip8DRMIajYUwqlfL5fFNX0cvQ+iSjI/R6/dKlS2/fvn38+PG2U5AQ4u7u3tjY2GKhmZmZo6NjQ0PDw9ZiMpk+Pj6dUCsAADwc7hG20zvvvHP+/PnDhw8b//iSSCQff/yxRCIhhNy4cUOtVhNC6uvrN2zYMGzYMEKIWCw2PDtx+PDhW7duJSQkmKJ8AAD4LwRhe0gkks8//zwjI8PGxobBYDAYjK+++ooQIhaLP/roI7FYTAjZuHGjSCTy8fFxc3NzcnKinhesr68fOHCgi4uLi4vLggULtm3bFhAQYOIvAwBAb7hH2IWam5urq6tdXFy4XK5hoU6nq6qqYjAYTk5Obazb0NDg6+uLbqUGuA9kDEfDGI6GMdwjbAfcI+xClpaWVLdPY0wm09nZ2ST1AABAa7g0CgAAtIYgBAAAWkMQAgAArSEIAQCA1hCEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0hiAEAABaQxACAACtIQgBAIDWEIQAAEBrCEIAAKA1BCEAANAaghAAAGgNQQgAALSGIAQAAFpDEAIAAK0hCAEAgNYQhAAAQGsIQgAAoDUEIQAA0BqCEAAAaA1BCAAAtIYgBAAAWkMQAgAArSEIAQCA1hCEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0hiAEAABaQxACAACtIQgBAIDWEIQAAEBrCEIAAKA1BCEAANAaghAAAGgNQQgAALSGIAQAAFpDEAIAAK0hCAEAgNYQhAAAQGsIQgAAoDUEIQAA0BqCEAAAaA1BCAAAtIYgBAAAWkMQAgAArSEIAQCA1hCEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0hiAEAABaQxACAACtIQgBAIDWEIQAAEBrCEIAAKA1BCEAANAaghAAAGgNQQgAALSGIAQAAFpDEAIAAK0hCAEAgNYQhAAAQGsIQgAAoDUEIQAA0BqCEAAAaA1BCAAAtIYgBAAAWkMQAgAArSEIAQCA1hCEAABAawhCAACgNQQh9AIKheLrr782dRU9hVgs/vbbb01dRU9RW1u7bds2U1fRU5SVle3evdvUVfQ+CML2qK2tfeeddwYNGhQRETF//vzS0tLWbTZs2DDaiF6vp5bn5+dPnjw5JCTk+eefr6qq6t7Ce6uamhr8029QWlq6fft2U1fRUxQUFOzZs8fUVfQU2dnZ+/btM3UVvQ+CsD2KioqUSuX//M//7N27V6FQjBs3rnWbnJwcV1fXt/8ftVCv10+cODEkJOT333/n8XjPP/989xYOAAAtsU1dQK8UExMTExNDvV63bp2zs3NVVZWjo2OLZj4+PqNGjTJecvr06bq6uo8//pjJZK5bt87BwSE3NzcgIKCb6gYAgFZwRthRV65csbe3t7e3b/3Rvn37hg8fPm/evMzMTGrJrVu3+vfvz2QyCSFWVlYBAQG3bt3q1nIBAODvcEbYIZWVlYsXL/7yyy+pbDNG3Rq0tbU9duxYfHz8lStXQkJCqqurhUKhoY1IJKqurn7glvV6fWNjo5mZGSGEwWC4uLgwGIyu+yI9nFarrays9Pb2NnUhPYJGo6mursbRoKhUqvr6ehwNikKh4HK5pq6i90EQtl9tbe3o0aMXLFgwe/bs1p9OnDiRejF48OCcnJxdu3atXbtWKBQ2Nzcb2kgkEpFI9MCN29jYXLx4US6XE0J4PJ6trW0XfIPeRKlU4v9wAxwNYzgaxox/asNjQhC2U319/ZgxY5KSklasWPHIxvb29hKJhBDi7e29efNmaqFGo7l7924bv2QHDBjQWdUCAMDD4B5hezQ1NSUmJkZHRy9btqyhoaGhoUGr1RJCrl279tlnn1Ft0tLSdDodIeTKlSs///wz1Wtm3LhxdXV1R44cIYTs2LHDyckpNjbWdN8DAAAIw/B8Gzy+EydOzJw503hJWlpaeHj4/v37165dm5GRQQiJjY29deuWpaUlIWT58uXLly+nWh47dmzevHlMJpPD4fz888847QMAMC0EYRdSKBQqlUogELRYrtVqxWKxUCikc/8XAIAeAkEIAAC0hs4y0OOcPXv2zJkzdXV1Li4uzz//vJOTU+s29+7d27Fjh1QqnTp1quHy8oEDB+rr66nXTk5OSUlJ3Vd0lzl06FBGRoZEIgkICHj++ef5fH7rNpmZmXv37uVwOC+++KKfnx+1UKPR7NixIzs7OywsbM6cOSwWq3sL73x6vf7y5ct//vlnfX19eHj4zJkzORxO62bnzp07cuQIdTR8fHyohT/88INCoaBetx7pojdSq9WnT58+d+6cXC5PSEhISkp64BWmgwcPnj9/3s7Obv78+YbO53K5fMuWLUVFRQMHDpw2bRouTaGzDPQ4Bw4cUCqVPj4+mZmZ4eHhZWVlLRpUVVXFxsY2NDQ4OTmNGTPm9OnT1PIPP/wwLS3t7t27d+/eLS8v7+66u8auXbs4HI63t3dKSkp8fLxSqWzR4OrVq0OGDLG2ttZoNAMGDCgqKqKWJycn79y509/ff/PmzYsXL+7uurtAUVHRjBkzGhsb3d3dN23ahKvgAAAAB+BJREFUNGbMGKqTmrHvv/9+2rRprq6uer0+Pj6+sLCQWv7GG29cvXqV+tt42MO7vcvx48eXL1+u0+ns7e3feOONV199tXWbt956a+XKlf7+/kVFRfHx8TKZjFo+YcKEY8eO+fv7f/TRR6tXr+7WunsmPUAPFhUVtX379hYL16xZk5SURL3+4osvxo4dS70OCwtLS0vr1vq6kUql4vF4ly9fbrF8xowZ77//PvV67ty5y5Yt0+v1hYWFXC63pqZGr9eXl5dzudyysrJuLrjTqdVqjUZDvaaGm7hx40aLNn5+fnv37qVev/LKK6+//jr12tbWNjc3t9tK7QZyudzw+tKlSxwOR6lUGjeQSqUsFuvWrVvU2yFDhmzdulWv11+8eNHGxoZa/dq1a9bW1lKptBsL74lwRgg9V3FxcWlpaVhYWIvlZ86cGTNmDPV69OjRZ86cMXx08ODB9evXnzx5svuq7C5Xr15ls9mtHzxNS0sbPXo09Xr06NFpaWmEkPT09H79+tnZ2RFCnJ2dg4KCzp8/380Fdzo2m224wKtWq7VabesLxRUVFYbLob6+vn/99Zfho927d3/11Vfnzp3rnmq7mrm5ueE1NaAMm/23W13V1dVarbb10UhLSxsyZAi1emRkJJfLvX79ejcW3hMhCKEn+uSTT9zc3Pz9/d99993Wj1pWVFQYBnd1cHCQyWSNjY2EkMjISCaTWV5ePm/evFmzZnV30V1mzpw5zs7OI0eO3LNnD5VtBlqttrq62vhoVFRUkL8fIkKIo6Njn7lWTHnttdemTp3a+mdBeHh4amoqIUSv16emphq+dUJCgkKhKCwsnDRp0ptvvtnd5XYlpVL55ptvLl++vMVAj+7u7kKhkDoaSqUyPT2dOhqVlZXGfxsODg597G+jPUx9Sgp0lJyczGolPj7e0EAmk1VWVv7222/29vanTp1qsXr//v337NlDvS4pKSGEtLi2U1FRwefzz58/39VfpFPExsa2PhovvfSSoYFEIikrK/vuu+9sbGzu3r1rvK5Op+NwOJmZmdTbY8eO+fj46PX69evXG64Y6/X64cOH/+///m+3fJuOGjt2bOujMXnyZOM277//fkRERH19fevVL1y44Ojo+PTTT8fGxo4aNcrDw6NFgzt37jCZzMLCwq77Cp3Izc2t9dFYuXKloYFGo5k+ffrEiRPVanXr1ffu3Wtraztp0qTQ0NBhw4aNGzdOr9e/9dZbycnJhjaBgYEHDx7shu/Sk6HXKJjA9u3b255a1sLCwsLCYsKECVOnTk1JSRk5cqTxp66urobfsGVlZUKhkBq4wMDJycnX17ewsDA+Pr7Ti+90ly9fbrsBn8/n8/kLFizYs2fPsWPHXnnlFcNHDAbD2dm5rKwsIiKCEFJWVubi4kIIcXV1Ne5kZFje8x07dqztBh9//PGhQ4dSU1MfOE5vXFxcfn7+jRs3XF1dDx8+fPDgwRYNAgMDbW1tCwsLvby8OqvmrkP9znsYrVb7wgsvSCSSX3/9tcV1UcqMGTNGjx6dnZ0dGBi4YsUK6nKoq6urYT4cajj73vK30XVwaRR6Fp1OZ+jmrtVqr1275uHhQQiRy+WpqakajYYQMmHChJSUFGoEu/3790+YMIEQolKpqCWEkNzc3JycnNDQUNN8h86jUCgMX6qpqSkvL486GjU1NYZ7XUlJSQcOHCCE6PX6AwcOUEdj1KhRBQUF2dnZhJDMzMyKiorhw4eb5jt0qnXr1u3evfuPP/4wvrjX2Nho6DlMCOHz+QkJCUKhcOPGjdSA+AqFQv//D0xfvHixoaEhMDCwewvvfDqdLjk5uba2NiUlxfgxkrKysitXrhje2traDh48uLm5+aeffqLGwxo/fnx6enpVVRUh5OTJk1ZWVv379+/++nsWU5+SAvyNRCKxtbVNSkp64YUXfH194+LixGKxXq/Pzc0lhNTW1ur1+ubm5ujo6KFDh06fPt3R0TE7O1uv12dkZLi7u0+ZMmXKlCkCgeCdd94x8TfpDOfPn/f09JwyZQr1TWfOnKnVavV6fUpKipOTE9WmuLjY1dV18uTJo0ePDgkJMVww/PTTT93c3ObNm+fi4rJ+/XqTfYfOc/v2bUKIt7d39P87c+aMXq8/ffq0mZkZ1WbDhg2jRo2aMWOGq6vr/PnzqcN15MgRPz8/6hKilZVV3zga+/btI4SEhIQYjgbVMfibb76Jioqi2ixdunTChAlTpkyxtbX9/PPPDesuWbLEz88vOTnZwcFh9+7dpvkCPQlGloEep7y8/OrVqzKZzMfHJyYmhnraV6lUZmZmRkVFUZeAlEplamqqVCodNWoUdYlMr9dnZWVRd4AiIyMNneV6u7y8vFu3bul0uuDg4JCQEGphQ0NDcXFxZGQk9VYsFp86dYrD4YwcOdLCwsKw7vXr16kH6sPDw01QemeTy+VUFhr4+flZW1tLJJKcnJyYmBhCSHNz87lz5+rq6kJCQvr160c102q1N27cyMnJ4fF40dHRrq6uJqi+s9XX1xuekqSEh4dzOJzq6uqamhrqckhdXd358+elUmlcXFyLjkUXLlwoKiqKiYnx9/fv1rp7JAQhAADQGu4RAgAArSEIAQCA1hCEAABAawhCAACgNQQhAADQGoIQAABoDUEIAAC0hiAEAABaQxACAACtIQgBAIDWEIQAAEBr/wcVvMM5UvckTgAAAABJRU5ErkJggg==", + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/html": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "m = mean(results.posteriors[:β][end])\n", + "Σ = cov(results.posteriors[:β][end])\n", + "\n", + "p = plot(xlims=(-3.1,-2.9), ylims=(2.5,2.7))\n", + "covellipse!(m, Σ, n_std=1, aspect_ratio=1, label=\"1σ Contours\", color=:green, fillalpha=0.2)\n", + "covellipse!(m, Σ, n_std=3, aspect_ratio=1, label=\"3σ Contours\", color=:blue, fillalpha=0.2)\n", + "scatter!([m[1]], [m[2]], label=\"Mean estimate\", color=:blue)\n", + "scatter!([true_beta[1]], [true_beta[2]], label=\"True Parameters\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References\n", + "[1] Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using Polya-Gamma latent variables. *Journal of the American Statistical Association*, 108(1), 136-146.\n", + "\n", + "[2] Minka, T. (2001). Expectation Propagation for approximate Bayesian inference. *Uncertainty in Artificial Intelligence*, 2, 362-369." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.11.2", + "language": "julia", + "name": "julia-1.11" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/Basic Examples/Bayesian Binomial Regression/Project.toml b/examples/Basic Examples/Bayesian Binomial Regression/Project.toml new file mode 100644 index 0000000..0416ec2 --- /dev/null +++ b/examples/Basic Examples/Bayesian Binomial Regression/Project.toml @@ -0,0 +1,9 @@ +[deps] +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +ReactiveMP = "a194aa59-28ba-4574-a09c-4a745416d6e3" +RxInfer = "86711068-29c9-4ff7-b620-ae75d7495b3d" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" diff --git a/examples/Basic Examples/Bayesian Binomial Regression/meta.jl b/examples/Basic Examples/Bayesian Binomial Regression/meta.jl new file mode 100644 index 0000000..50eb895 --- /dev/null +++ b/examples/Basic Examples/Bayesian Binomial Regression/meta.jl @@ -0,0 +1,15 @@ +meta = ( + title = "Bayesian Binomial Regression", + description = """ + An introductory tutorial to Bayesian binomial regression with RxInfer. + Learn how to model binary outcomes using logistic regression with proper Bayesian inference. + The example demonstrates the use of Expectation Propagation (EP) algorithm and Polya-Gamma augmentation. + """, + tags = [ + "basic examples", + "regression", + "multivariate", + "expectation propagation", + "polya-gamma", + ] +) diff --git a/examples/Basic Examples/Bayesian Linear Regression Tutorial/Bayesian Linear Regression Tutorial.ipynb b/examples/Basic Examples/Bayesian Linear Regression/Bayesian Linear Regression.ipynb similarity index 100% rename from examples/Basic Examples/Bayesian Linear Regression Tutorial/Bayesian Linear Regression Tutorial.ipynb rename to examples/Basic Examples/Bayesian Linear Regression/Bayesian Linear Regression.ipynb diff --git a/examples/Basic Examples/Bayesian Linear Regression Tutorial/Project.toml b/examples/Basic Examples/Bayesian Linear Regression/Project.toml similarity index 100% rename from examples/Basic Examples/Bayesian Linear Regression Tutorial/Project.toml rename to examples/Basic Examples/Bayesian Linear Regression/Project.toml diff --git a/examples/Basic Examples/Bayesian Linear Regression Tutorial/hblr-matrix-completion.jpeg b/examples/Basic Examples/Bayesian Linear Regression/hblr-matrix-completion.jpeg similarity index 100% rename from examples/Basic Examples/Bayesian Linear Regression Tutorial/hblr-matrix-completion.jpeg rename to examples/Basic Examples/Bayesian Linear Regression/hblr-matrix-completion.jpeg diff --git a/examples/Basic Examples/Bayesian Linear Regression Tutorial/hbr/osic_pulmonary_fibrosis.csv b/examples/Basic Examples/Bayesian Linear Regression/hbr/osic_pulmonary_fibrosis.csv similarity index 100% rename from examples/Basic Examples/Bayesian Linear Regression Tutorial/hbr/osic_pulmonary_fibrosis.csv rename to examples/Basic Examples/Bayesian Linear Regression/hbr/osic_pulmonary_fibrosis.csv diff --git a/examples/Basic Examples/Bayesian Linear Regression Tutorial/meta.jl b/examples/Basic Examples/Bayesian Linear Regression/meta.jl similarity index 84% rename from examples/Basic Examples/Bayesian Linear Regression Tutorial/meta.jl rename to examples/Basic Examples/Bayesian Linear Regression/meta.jl index 2555faf..b98bac8 100644 --- a/examples/Basic Examples/Bayesian Linear Regression Tutorial/meta.jl +++ b/examples/Basic Examples/Bayesian Linear Regression/meta.jl @@ -1,5 +1,5 @@ return ( - title = "Bayesian Linear Regression Tutorial", + title = "Bayesian Linear Regression", description = """ An extensive tutorial on Bayesian linear regression with RxInfer with a lot of examples, including multivariate and hierarchical linear regression. """,