examples/diagnostics_and_criticism/Bayes_factor.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Bayes_factor)=\n",
    "# Bayes Factors and Marginal Likelihood\n",
    ":::{post} Jun 1, 2022\n",
    ":tags: Bayes Factors, model comparison \n",
    ":category: beginner, explanation\n",
    ":author: Osvaldo Martin\n",
    ":::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running on PyMC v4.0.0b6\n"
     ]
    }
   ],
   "source": [
    "import arviz as az\n",
    "import numpy as np\n",
    "import pymc as pm\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.ticker import FormatStrFormatter\n",
    "from scipy.special import betaln\n",
    "from scipy.stats import beta\n",
    "\n",
    "print(f\"Running on PyMC v{pm.__version__}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "az.style.use(\"arviz-darkgrid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"Bayesian way\" to compare models is to compute the _marginal likelihood_ of each model $p(y \\mid M_k)$, _i.e._ the probability of the observed data $y$ given the $M_k$ model. This quantity, the marginal likelihood, is just the normalizing constant of Bayes' theorem. We can see this if we write Bayes' theorem and make explicit the fact that all inferences are model-dependant. \n",
    "\n",
    "$$p (\\theta \\mid y, M_k ) = \\frac{p(y \\mid \\theta, M_k) p(\\theta \\mid M_k)}{p( y \\mid M_k)}$$\n",
    "\n",
    "where:\n",
    "\n",
    "* $y$ is the data\n",
    "* $\\theta$ the parameters\n",
    "* $M_k$ one model out of K competing models\n",
    "\n",
    "\n",
    "Usually when doing inference we do not need to compute this normalizing constant, so in practice we often compute the posterior up to a constant factor, that is:\n",
    "\n",
    "$$p (\\theta \\mid y, M_k ) \\propto p(y \\mid \\theta, M_k) p(\\theta \\mid M_k)$$\n",
    "\n",
    "However, for model comparison and model averaging the marginal likelihood is an important quantity. Although, it's not the only way to perform these tasks, you can read about model averaging and model selection using alternative methods [here](model_comparison.ipynb), [there](model_averaging.ipynb) and [elsewhere](GLM-model-selection.ipynb). Actually, these alternative methods are most often than not a better choice compared with using the marginal likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian model selection\n",
    "\n",
    "If our main objective is to choose only one model, the _best_ one, from a set of models we can just choose the one with the largest $p(y \\mid M_k)$. This is totally fine if **all models** are assumed to have the same _a priori_ probability. Otherwise, we have to take into account that not all models are equally likely _a priori_ and compute:\n",
    "\n",
    "$$p(M_k \\mid y) \\propto p(y \\mid M_k) p(M_k)$$\n",
    "\n",
    "Sometimes the main objective is not to just keep a single model but instead to compare models to determine which ones are more likely and by how much. This can be achieved using Bayes factors:\n",
    "\n",
    "$$BF =  \\frac{p(y \\mid M_0)}{p(y \\mid M_1)}$$\n",
    "\n",
    "that is, the ratio between the marginal likelihood of two models. The larger the BF the _better_ the model in the numerator ($M_0$ in this example). To ease the interpretation of BFs some authors have proposed tables with levels of *support* or *strength*, just a way to put numbers into words. \n",
    "\n",
    "* 1-3: anecdotal\n",
    "* 3-10: moderate\n",
    "* 10-30: strong\n",
    "* 30-100: very strong\n",
    "* $>$ 100: extreme\n",
    "\n",
    "Notice that if you get numbers below 1 then the support is for the model in the denominator, tables for those cases are also available. Of course, you can also just take the inverse of the values in the above table or take the inverse of the BF value and you will be OK.\n",
    "\n",
    "It is very important to remember that these rules are just conventions, simple guides at best. Results should always be put into context of our problems and should be accompanied with enough details so others could evaluate by themselves if they agree with our conclusions. The evidence necessary to make a claim is not the same in particle physics, or a court, or to evacuate a town to prevent hundreds of deaths."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian model averaging\n",
    "\n",
    "Instead of choosing one single model from a set of candidate models, model averaging is about getting one meta-model by averaging the candidate models. The Bayesian version of this weights each model by its marginal posterior probability.\n",
    "\n",
    "$$p(\\theta \\mid y) = \\sum_{k=1}^K p(\\theta \\mid y, M_k) \\; p(M_k \\mid y)$$\n",
    "\n",
    "This is the optimal way to average models if the prior is _correct_ and the _correct_ model is one of the $M_k$ models in our set. Otherwise, _bayesian model averaging_ will asymptotically select the one single model in the set of compared models that is closest in [Kullback-Leibler divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence).\n",
    "\n",
    "Check this [example](model_averaging.ipynb) as an alternative way to perform model averaging."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Some remarks\n",
    "\n",
    "Now we will briefly discuss some key facts about the _marginal likelihood_\n",
    "\n",
    "* The good\n",
    "    * **Occam Razor included**: Models with more parameters have a larger penalization than models with fewer parameters. The intuitive reason is that the larger the number of parameters the more _spread_ the _prior_ with respect to the likelihood.\n",
    "\n",
    "\n",
    "* The bad\n",
    "    * Computing the marginal likelihood is, generally, a hard task because it’s an integral of a highly variable function over a high dimensional parameter space. In general this integral needs to be solved numerically using more or less sophisticated methods.\n",
    "    \n",
    "$$p(y \\mid M_k) = \\int_{\\theta_k} p(y \\mid \\theta_k, M_k) \\; p(\\theta_k | M_k) \\; d\\theta_k$$\n",
    "\n",
    "* The ugly\n",
    "    * The marginal likelihood depends **sensitively** on the specified prior for the parameters in each model $p(\\theta_k \\mid M_k)$.\n",
    "\n",
    "Notice that *the good* and *the ugly* are related. Using the marginal likelihood to compare models is a good idea because a penalization for complex models is already included (thus preventing us from overfitting) and, at the same time, a change in the prior will affect the computations of the marginal likelihood. At first this sounds a little bit silly; we already know that priors affect computations (otherwise we could simply avoid them), but the point here is the word **sensitively**. We are talking about changes in the prior that will keep inference of $\\theta$ more or less the same, but could have a big impact in the value of the marginal likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing Bayes factors\n",
    "\n",
    "The marginal likelihood is generally not available in closed-form except for some restricted models. For this reason many methods have been devised to compute the marginal likelihood and the derived Bayes factors, some of these methods are so simple and [naive](https://radfordneal.wordpress.com/2008/08/17/the-harmonic-mean-of-the-likelihood-worst-monte-carlo-method-ever/) that works very bad in practice. Most of the useful methods have been originally proposed in the field of Statistical Mechanics. This connection is explained because the marginal likelihood is analogous to a central quantity in statistical physics known as the _partition function_ which in turn is closely related to another very important quantity the _free-energy_. Many of the connections between Statistical Mechanics and Bayesian inference are summarized [here](https://arxiv.org/abs/1706.01428)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using a hierarchical model\n",
    "\n",
    "Computation of Bayes factors can be framed as a hierarchical model, where the high-level parameter is an index assigned to each model and sampled from a categorical distribution. In other words, we perform inference for two (or more) competing models at the same time and we use a discrete _dummy_ variable that _jumps_ between models. How much time we spend sampling each model is proportional to $p(M_k \\mid y)$.\n",
    "\n",
    "Some common problems when computing Bayes factors this way is that if one model is better than the other, by definition, we will spend more time sampling from it than from the other model. And this could lead to inaccuracies because we will be undersampling the less likely model. Another problem is that the values of the parameters get updated even when the parameters are not used to fit that model. That is, when model 0 is chosen, parameters in model 1 are updated but since they are not used to explain the data, they only get restricted by the prior. If the prior is too vague, it is possible that when we choose model 1, the parameter values are too far away from the previous accepted values and hence the step is rejected. Therefore we end up having a problem with sampling.\n",
    "\n",
    "In case we find these problems, we can try to improve sampling by implementing two modifications to our model:\n",
    "\n",
    "* Ideally, we can get a better sampling of both models if they are visited equally, so we can adjust the prior for each model in such a way to favour the less favourable model and disfavour the most favourable one. This will not affect the computation of the Bayes factor because we have to include the priors in the computation.\n",
    "\n",
    "* Use pseudo priors, as suggested by Kruschke and others. The idea is simple: if the problem is that the parameters drift away unrestricted, when the model they belong to is not selected, then one solution is to try to restrict them artificially, but only when not used! You can find an example of using pseudo priors in a model used by Kruschke in his book and [ported](https://github.com/aloctavodia/Doing_bayesian_data_analysis) to Python/PyMC3.\n",
    "\n",
    "If you want to learn more about this approach to the computation of the marginal likelihood see [Chapter 12 of Doing Bayesian Data Analysis](http://www.sciencedirect.com/science/book/9780124058880). This chapter also discuss how to use Bayes Factors as a Bayesian alternative to classical hypothesis testing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytically\n",
    "\n",
    "For some models, like the beta-binomial model (AKA the _coin-flipping_ model) we can compute the marginal likelihood analytically. If we write this model as:\n",
    "\n",
    "$$\\theta \\sim Beta(\\alpha, \\beta)$$\n",
    "$$y \\sim Bin(n=1, p=\\theta)$$\n",
    "\n",
    "the _marginal likelihood_ will be:\n",
    "\n",
    "$$p(y) = \\binom {n}{h}  \\frac{B(\\alpha + h,\\ \\beta + n - h)} {B(\\alpha, \\beta)}$$\n",
    "\n",
    "where:\n",
    "\n",
    "* $B$ is the [beta function](https://en.wikipedia.org/wiki/Beta_function) not to get confused with the $Beta$ distribution\n",
    "* $n$ is the number of trials\n",
    "* $h$ is the number of success\n",
    "\n",
    "Since we only care about the relative value of the _marginal likelihood_ under two different models (for the same data), we can omit the binomial coefficient $\\binom {n}{h}$, thus we can write:\n",
    "\n",
    "$$p(y) \\propto \\frac{B(\\alpha + h,\\ \\beta + n - h)} {B(\\alpha, \\beta)}$$\n",
    "\n",
    "This expression has been coded in the following cell, but with a twist. We will be using the `betaln` function instead of the `beta` function, this is done to prevent underflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def beta_binom(prior, y):\n",
    "    \"\"\"\n",
    "    Compute the marginal likelihood, analytically, for a beta-binomial model.\n",
    "\n",
    "    prior : tuple\n",
    "        tuple of alpha and beta parameter for the prior (beta distribution)\n",
    "    y : array\n",
    "        array with \"1\" and \"0\" corresponding to the success and fails respectively\n",
    "    \"\"\"\n",
    "    alpha, beta = prior\n",
    "    h = np.sum(y)\n",
    "    n = len(y)\n",
    "    p_y = np.exp(betaln(alpha + h, beta + n - h) - betaln(alpha, beta))\n",
    "    return p_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our data for this example consist on 100 \"flips of a coin\" and the same number of observed \"heads\" and \"tails\". We will compare two models one with a uniform prior and one with a _more concentrated_ prior around $\\theta = 0.5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.repeat([1, 0], [50, 50])  # 50 \"heads\" and 50 \"tails\"\n",
    "priors = ((1, 1), (30, 30))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABOhUlEQVR4nO3deXhU9cH+//ucM5OVJQmbhB2RXVRQQVxwbyu1tkLRS2utbd1abbXt861t7bebT/vU3+9prbWLa126ah9b6r70USoCQlllC4shBAgQSAjZSGbOOd8/kpkQCWSbmc8s79d19RJnJjN3xtNwzyefxfJ93xcAAACAmLNNBwAAAADSFWUbAAAAiBPKNgAAABAnlG0AAAAgTijbAAAAQJxQtgEAAIA4oWwDAAAAcULZBgAAAOIk0NUHVldXxzNHp/r376+amhqjGWAe1wG4BiBxHYBrAMlxDRQWFnb6mJQZ2bbtlImKOOI6ANcAJK4DcA0gda6B1EgJAAAApCDKNgAAABAnlG0AAAAgTijbAAAAQJxQtgEAAIA4oWwDAAAAcULZBgAAAOKEsg0AAADECWUbAAAAiBPKNgAAABAnlG0AAAAgTijbAAAAQJxQtgEAAIA4oWwDAAAAcULZBgAAAOKEsg0AAADEScB0AABAFzRUKbjpH7KOHJIkeQPHKzzuMskJms0FADghyjYAJLPGQ8pe8gsF1v9NltvU7i6v71A1z7xV4VMXSJZlKCAA4EQo2wCQpKzafcp5/otyDm6TJLlDpsoderosLyRn25uyayuU8+b3FarcrKaL7pVsx3BiAKbV1dXpiSee0KJFi7Rv3z7l5ORo8uTJuv322zVp0iTT8TISZRsAkpBVXabcv94ku7ZCXv5gNX3sfrkjzm4bwZ5zj4Jrfq+sd36m4No/S42H1PSx+5lWAiShV155RWvXrtXmzZu1fft2hUIh3Xvvvfr4xz8e09epq6vTzTffrNLSUs2YMUMXXnih9u7dq7feekt33XWXFi5cqJycnJi+Znck6n1INpRtAEg2oUbl/OPLLUW7cLQa5z0mv9+w9o8J5ih01hfl9xum7Fe+qeCWV+X3H6Hm879mJjOA43r44Ye1d+9eFRQUaMCAAdq7d29cXuexxx5TaWmp7rzzTl1//fXR2++77z69+OKLKisr04QJE+Ly2l2RqPch2bAbCQAkmey3fyLn4HZ5+QPVuODpY4v2UcITPqYjV9wvScpa8aicsncTFRNAF33729/W3/72N7366qu6+uqr4/IaruvqpZde0vDhw9sVbUkKBFrGVvv27RuX1+6qRLwPyYiyDQBJxNnyqoLvPydflpo+dr/8/EGdfo07/qMKTbtGkpT9yj2y6g/EOyaQVBoaGvToo49qwYIFuuCCCzR//nz9/e9/lyStWLFCs2bN0rvvmvsgevbZZ2vo0KFxfY2tW7eqtrZW5513Xrvbw+Gw/vWvf+mkk05ScXFxXDN0JhHvQzKibANAsmiuV/b/3idJCp19s9yR53T5S5vm3CN3wCmyGw4oa/HP45UQSDr79+/XTTfdpCeffFITJ07U1Vdfrfr6ev30pz9VSUmJnnzySU2ePFnnnnuu6ahxtXnzZknS5MmTo7c1NTXpRz/6kaqqqnTttdeaipbxmLMNAEkia8VjshsOyisYqeZzvty9Lw7mqOnyHynvT9cqsOFvCk2/Qd6gifEJipTg+76OHDGd4vhyciSrl1tWep6ne+65R2VlZbr//vt1wQUXSJJmz56tr3zlK/r973+vlStX6mc/+1mnz/Xoo49267WvvfZa49MyjlZSUiJJmjhxot544w0tXbpUS5cuVXV1ta644gotWLCgS8+T6u9DMqJsA0ASsGr3KbjySUlS0/lfl5ysbj+HN/Q0hSZ8TMGSV5S16H4dmfc4+29nKN/39aU7fb2/3nSS4zt1qvTrX/aucL/zzjvauHGjLr300mjRlhTd4u6NN97QlClTNHv27E6f6/HHH+/Wa8+dOzepSmZJSYny8vI0YsQIffe7342W72AwqOLiYvm+36XnSfX3IRlRtgEgCWQteVBW+Ijc4ulyx13W4+dpPu9rCmx7U4GdS+XseEfumAs6/yIgRb3++uuSdMyobTDYtgXmzTff3KXnWrZsWeyCJZjrutq+fbsmTZoky7L0+OOPq6amRuvXr9fjjz+uxx57TKFQSLfffnunz5XK70OyomwDgGFWzW4FNv5dktR0wX/0ajTa7z9codOvV9bKJ5W17DdqpGxnJMuy9OtfKu2nkaxevVr9+vXT1KlT290eGcU99dRTNWvWrF69RirYsWOHmpqaNHFiy9SxQCCgAQMGaM6cOZo2bZrmzp2rxYsXd6lsI/Yo2wBgWHD107J8T+FRs+UVn97r5wud+XkF1/xBTsUa2XtWyys+o/chkXIsy1JurukU8VNbW6uqqipNnTpVtt1+v4clS5ZIUrcWRabyXOXIlJGO9tDu27evHMdRXl5el54rld+HZEXZBgCTjhxW8P2/SpJCM26KyVP6+YMUnnilghueV9bKp3SEso001NjYKEnHFO2mpib9+te/liQ5jtPl50vlucqRsj1+/Phj7lu+fLlCoZCmT5/epedK5fchWVG2AcCg4Pq/ygo1yB0wTu6o2G1NFpp+o4Ibnpez7Q1ZNbvk9x8es+cGkkFRUZGys7NVUlKiioqK6P7NDzzwgHbv3i1Jqq6u7vLzJdNc5dtvv12rV6/Wr371K82YMaPTx2/ZskVSyxz2o6eK1NTU6Be/+IWys7N11VVXdem1k+l9SBeUbQAwxQ0puOoZSVJoxudiunOIN2i8wqPOVaDsXQVXPa3mi74ds+cGkkEgENDll1+uF154QbfddpvmzJmj7du3a+XKlbrlllv09NNP64UXXlBWVpZuvPFG5RqcU7Nw4UKtXbtWkrR9+3ZJ0j/+8Q+tWrVKkjRnzhzNmTMn+vjInPOujMz7vh8t20899ZTWrVunKVOmqKamRosWLVJdXZ2+853vGD/QRur++5AuKNsAYIiz4x3ZdXvl5RYpPPHjMX/+0PQbW8r2xoVqPv/rUiA75q8BmPS1r31N2dnZeuutt/T888+rqKhId999t6655hoNHjxYDz30kJ599lndeuutRnOuXbtWL7/8crvb1q1bp3Xr1kmShg4dGi2Zvu+rtLRUQ4cOPWbhZ0fKy8tVX1+v888/X57nac2aNVq3bp0KCgp05pln6rrrruvS8yRCd96HdGL5Xdx4sTu/iomHwsJC4xlgHtcB0ukayPnHnQpse1PNMz6n5jnfjP0LeK7yHr9Mdm2Fjsz9b4UnXBH71zAkna4D9Ey6XgPbt2/X9ddfr2984xuaP39+p49/8803de+99+qee+7RJz/5yfgHTCLJcA0UFhZ2+hiOawcAExqq5HzwtiQpPOVT8XkN21F4css8zcD6v8XnNQDE1Jo1a1RUVKQrr7yyS48/0eJIJAfKNgAYENz8giwvLHfIVHkD4/eXZKi1yDtl78qq3Ru31wEQG/PmzdPLL7+s7OyuTfsqKSmR4zg6+eST45wMPUXZBoBE8/3oSHMoXqPakZcqGCl32Jmy5CuwcWFcXwtA4m3ZskWjRo3qcjlH4rFAEgASzK7cJOdAiXwnS+GJc+P+eqGpV8vZ/W8FN/xNobNviemuJwDMevXVV01HQCcY2QaABAuUtPzl6I69UMrpH/fXC59yufxAjuxDZbIrN8f99QAAbSjbAJBIvq/AltckSeHxH03Ma2blyx19viRFXxsAkBiUbQBIILtys+yanfKdbIXHXJCw1w2P/4gkKbDlValrO74CAGKAsg0ACRQZWXbHXCBl5SfsdcNjL5TvZLVMJTmwJWGvCwCZjrINAIni+y0jy2obaU4YppIAgBGUbQBIEPvAFtmHylp2IRl7YcJfPzqVZOtrTCUBgAShbANAgjjb3pQkuaPPS+gUkoiWqSRB2VUfyKouTfjrA0AmomwDQIIEPlgkSQqffLGZANl95Q4/u10WAEB8UbYBIAGs+ko5+96X1Lo40hB3zBxJkvPB28YyAEAmoWwDQAI4pe9IktwhU+TnDzKWIzy2tWzvWSU11RrLAQCZgrINAAkQKG2ZthEZWTbFLxgpr2isLC8sp+xdo1kAIBNQtgEg3tzmaLE1sQvJh4VbC3+AqSQAEHeUbQCIM2f3SlnN9fLyBsobMsV0HLmthd8pfUfyPbNhACDNBUwHAIB055T+S5LkjjlfssyPcbjFZ8jP7iu7sUr23vflDT3NdCQAhtTV1emJJ57QokWLtG/fPuXk5Gjy5Mm6/fbbNWnSJNPx0gJlGwDizNnRMoXEHX2e4SStnKDckecosPV1OWVLKNtAHNXW1uqRRx7Rpk2btGfPHtXW1qqgoEAjR47UvHnzdNFFF8myrGO+buPGjXr00Ue1fv16hUIhjR07Vtdcc40+8pHYnT5bV1enm2++WaWlpZoxY4YuvPBC7d27V2+99ZbuuusuLVy4UDk5OTF7ve7oyvt2PIl477qDsg0AcWTVV8o5uFW+LIVHnmM6TlS4tWwHdi5VaNbtpuMAaevQoUN68cUXNWXKFM2ZM0f9+vVTVVWVFi9erG9/+9u66qqr9K1vfavd16xcuVJ33XWXgsGgLr30UvXp00dvv/22vve976miokKf+9znYpLtscceU2lpqe68805df/310dvvu+8+vfjiiyorK9OECRNi8lrd1ZX37f777z/m6xL13nUHZRsA4sjZuUyS5A2eJOUWGk7Txm0t/vaeNVJzvZETLYFMUFxcrDfeeEOBQPvKVV9fry9+8YtauHChrrnmGo0dO1aSFA6H9ZOf/ESWZek3v/lNtOx+4Qtf0M0336xHH31UF198sUaOHNmrXK7r6qWXXtLw4cPbFW1J0ax9+/bt1Wv0Rlfet5tvvlkDBw6M3peo9667zE8eBIA05pQtkdRWbpOFXzBSXr9hsryQnN0rTccBeqWhoUGPPvqoFixYoAsuuEDz58/X3//+d0nSihUrNGvWLL37rpmtLh3HOaYwSlJ+fr5mzpwpSdq1a1f09pUrV2rXrl26/PLL240q5+fn66abboqW5N7aunWramtrdd557ae3hcNh/etf/9JJJ52k4uLiXr9OT3XlfSsrK2t3X6Leu+6ibANAvPi+nJ1LJUnuqNmGw3yIZUU/AEQ+EACpaP/+/brpppv05JNPauLEibr66qtVX1+vn/70pyopKdGTTz6pyZMn69xzzzUdtZ2mpiatXLlSlmVp9OjR0dtXrVolSdFCebTIbZHH9MbmzZslSZMnT26X6Uc/+pGqqqp07bXX9vo14uHo923cuHHt7kvUe9ddTCMBgDixqktl1+2T72TJLZ5uOs4x3FHnKLj+r9EPBEgzvi+FG02nOL5ArtTBwsDu8DxP99xzj8rKynT//ffrggsukCTNnj1bX/nKV/T73/9eK1eu1M9+9rNOn+vRRx/t1mtfe+213ZpmUVtbqz//+c/yfV/V1dVasmSJ9u3bpy984QvtpjWUl5dLkkaMGHHMc/Tr108FBQXtRsJ7qqSkRJI0ceJEvfHGG1q6dKmWLl2q6upqXXHFFVqwYEGXnyue792J3rfRo0eruro6+thEvXfdRdkGgDgJRKaQDJshBc2s6D+RyIJN58AWWfWVRo+RR4z5vnL/cr2cPatNJzkut3i6Gq/5fa8K9zvvvKONGzfq0ksvjRZtSdEt69544w1NmTJFs2d3/pulxx9/vFuvPXfu3G6X7aNfIxAI6M4779R1113X7nF1dXWSWqY+dCQ/P1/79+/vVtaOlJSUKC8vTyNGjNB3v/vdaPkOBoMqLi6W7/tdfq54vnddfd+kxL133UXZBoA4iU4hSbL52lG5hXIHT5Kzf5OcncsUnnSl6USIqd6NGqeC119/XZKOGYUNBoPRP998881deq5ly5bFLlgHiouLtWzZMrmuq3379unNN9/Ub3/7W61bt07/+Z//2eH85HhxXVfbt2/XpEmTZFmWHn/8cdXU1Gj9+vV6/PHH9dhjjykUCun227u2U1E837sTvW+//vWv4/a6sUTZBoB48Fw5u/4tSXJHzjIc5vjcEbNaynb5csp2OrGsllHjNJ9Gsnr1avXr109Tp05td3tkVPbUU0/VrFnJ9f8/x3FUXFysz372s7JtWw899JAWLlyoefPmSZL69OkjqWXXjY7U19dHH9NTO3bsUFNTkyZOnCipZbR4wIABmjNnjqZNm6a5c+dq8eLFXS7bidDR+/bss8/qYx/7WPQxiXjveoKyDQBxYB/YIqvpsPys/JZt/5KUO+IsaeXv5OxaYToKYs2ypGCe6RRxU1tbq6qqKk2dOlW23X6/hyVLWqZwdWdRZLznbHdk5syZeuihh7Rq1apo2Y7MNy4vL4+W4YjDhw/r0KFDOvXUU3v1upEpIx3tod23b185jqO8vK5fO4l+7yLv2/Lly9uV7US8dz1B2QaAOIiUV7d4umQn749at3iGfFmyD5XJqtsvv89g05GALmlsbBm1/3DRbmpqik4vcByny88X7znbHamsrJTUPucZZ5yhp556Su+9954uu+yydo9/7733JEnTp/duwXWkbI8fP/6Y+5YvX65QKNSt10j0e9fR+yYl5r3rieT9GwAAUli0bA8/y3CSTuT0kzd4YstUkl0rFJ4413QioEuKioqUnZ2tkpISVVRUaOjQoZKkBx54QLt375akdjtVdCZe8463bNmi4uLiY6Yv1NTU6Le//a0k6Zxz2tZ1nHnmmRo2bJhef/11LViwIFqI6+vr9bvf/U6O42ju3GP/f3r77bdr9erV+tWvfqUZM2Z0mklqmfN+9FSRmpoa/eIXv1B2drauuuqqLn+P8XjvuvK+Hb0oVur5exdvlG0AiDXfa5uvnexlWy0ZKdtINYFAQJdffrleeOEF3XbbbZozZ462b9+ulStX6pZbbtHTTz+tF154QVlZWbrxxhuVm5trJOdLL72kf/zjH5o+fbqGDh2qnJwc7d27V0uWLFFDQ4MuuugifeQjH2n3fX3729/WV7/6Vd1222267LLLlJ+fr7ffflt79uzRrbfe2uEJiJF56p2N5vu+Hy3bTz31lNatW6cpU6aopqZGixYtUl1dnb7zne8YPdBG6tr7duWVV6qmpib6NT197+KNsg0AMWYf3CbryCH5wTx5Q6aYjtMpd/hZ0qqnox8QgFTxta99TdnZ2Xrrrbf0/PPPq6ioSHfffbeuueYaDR48OLqI7tZbbzWW8aKLLlJdXZ3Wr1+vNWvW6MiRI+rXr5+mTZumK664QpdddpmsDy0UnTFjhh5++GE9+uij+uc//6lQKKSxY8fqlltu0Uc/+tFjXsP3fZWWlmro0KHHLBb9sPLyctXX1+v888+X53las2aN1q1bp4KCAp155pm67rrrOn2OROjK+/bhKURS99+7RLD8Lm6k2J1fxcRDYWGh8Qwwj+sAqXANBNf8Qdn/e5/Co2bryLzuzWU0orFafX7Tsg9x/W2L5ecNMByoc6lwHSC+uAbabN++Xddff72+8Y1vaP78+Sd87Jtvvql7771X99xzjz75yU8mJmCcJMM1UFhY2OljOK4dAGLMTqEpJJJa9tse2DK30WZ0G0g5a9asUVFRka68svPtO0+0OBLxQdkGgFjy/ZSarx0RycoWgEDqmTdvnl5++WVlZ2d3+tiSkhI5jqOTTz45AckgUbYBIKasQztlNxyQ72TJG5L4/Vx7yhvWsnuBs2eV4SQA4mnLli0aNWpUl4o5YoMFkgAQQ86e1ZIkb8hUKZBlOE3XucVnSJLsyhKpuV7KyjecCEA8vPrqq6YjZBxGtgEghiIjw27x6WaDdJPf9yR5fYfK8j05FWtNxwGAtEHZBoAYsqNlO/GnlPVWJLPdOjoPAOg9yjYAxErjITkHt0tqm5aRSrxhLZmZtw0AsUPZBoAYcSrWSJK8wtFSXpHRLD0RGdl2KtZKnms4DQCkB8o2AMSIs2eNpNScQiJJ3sDx8rPyZTXXyz641XQcAEgLlG0AiJG2xZGpN4VEkmQ7coee1vLH3UwlAYBYoGwDQCy4Idl732/5Y6qWbR01lYRFkgAQE5RtAIgBu7JEVviI/Oz+8ovGmI7TY97Q0yWJ7f8AIEYo2wAQA87elnLqDp0mWan7o9U9aaokya4plxqrDacBgNSXun8jAEASsSsiZfs0w0l6Kad/y24qkpzWaTEAgJ6jbANADDgV6yRJ3knTDCfpPbf1e7Ap2wDQa5RtAOitxmrZh8okSe5JpxoO03uRDwzO3nWGkwBA6qNsA0AvRUe1C8dIuQVmw8RA5AODU7FO8n3DaQAgtVG2AaCXIiPA7tDUn0IiSd6gifLtoKwjh2TV7DIdBwBSGmUbAHopbRZHRgSy5A2aKImpJADQW5RtAOgN34vu2pEOiyMjIqP0LJIEgN6hbANAL1jVO2Q1HZYfyJE3cLzpODHjReZtU7YBoFco2wDQC9FR7cGTJSdoOE3sRLf/27dBckOG0wBA6qJsA0Av2Ps2SJLcIVMMJ4ktv3CU/Oy+stwm2Qe3mo4DACmLsg0AveDsWy9J8oZMNZwkxixbbuv3xLxtAOg5yjYA9JQXlr1/syTJPSnNyrYkr3WRZGQfcQBA91G2AaCH7KoPZIUb5Qfz5BeONh0n5ji2HQB6j7INAD0Uma/tDZkiWen34zQyNcY+uE1qrjecBgBSU/r97QAACWLvbZmv7abbfO1Wfp/B8voOlSVfTusHCwBA91C2AaCH2hZHptdOJEeL7Ldtc5IkAPQIZRsAesINya5sXRyZxmXb5XAbAOgVyjYA9IBdtV2W2yw/u6/8gpGm48QNiyQBoHco2wDQA9H52oPTc3FkhDdkinxZsmsrZNVXmo4DACknff+GAIA4yoT52pKkrHz5RWMlSfa+jYbDAEDqoWwDQA9Ej2lPw8NsPswdPFmSZO+nbANAd1G2AaC73GbZB0okSd7gNB/ZluQNaSnbDmUbALqNsg0A3WQf2CrLDcnP7i+//3DTceLOY2QbAHqMsg0A3dQ2hWSKZFmG08SfO2iiJMk+vEdqrDacBgBSC2UbALqpbXFk+s/XliTl9JPXv2V7Q2f/JsNhACC1ULYBoJvsfZFj2tN/vnaEO4SpJADQE5RtAOiOcJPsA1slZdDItpi3DQA9RdkGgG6wK0tkeWH5uYXy+w41HSdhImXbYa9tAOgWyjYAdIO9v3Vx5JCpGbE4MsIdPEmSZB8qk5rqDKcBgNRB2QaAbnD2RhZHTjacJMHyiuS1juTblSySBICuomwDQDfYlZsltZ2qmEmiU0mYtw0AXUbZBoCuckOyD7YujszAsh2dSsL2fwDQZZRtAOgiu+qD1pMj+8rvN8x0nISLTJ2xWSQJAF1G2QaALoqM6HqDJmbU4sgIb3DLvuJ21XYp1Gg4DQCkBso2AHRRZI/pyPHlmcbPHyQvb6As35N9YIvpOACQEijbANBFTuviyEycry1Jsix50XnbTCUBgK6gbANAV/h+dCcSb3BmjmxLbfO2OdwGALqGsg0AXWAd3i2rqVa+E5RXdLLpOMa4HNsOAN1C2QaALoiUS2/AKZITNJzGnOgiyQNbJbfZcBoASH6UbQDoAufonUgymN+vWH52f1leSPbBbabjAEDSo2wDQBfYmb44MsKy5A5pXSS5b4PhMACQ/CjbANAFkT223QxeHBnhDWot25UlhpMAQPKjbANAZxqqZNftk8Q0EknyBk2Q1LYVIgDg+CjbANAJp7J1vnbBSCkr33Aa8yIfOOzKEsn3DacBgORG2QaATrRNIcnw+dqtvKIx8p2grOY6WYd3m44DAEmNsg0AnbD3RxZHTjKcJEk4WfKKxklqWzgKAOgYZRsAOhGdRsJ87ajIvG0WSQLAiVG2AeBEQg2yqkolMbJ9NBZJAkDXULYB4ATsyi2y5MvLGyg/f5DpOEmj3SJJAMBxUbYB4ATsyBQSRrXbcSPTSGrKpaY6w2kAIHlRtgHgBJzI4kjma7eXWyivz0mSJPsAo9sAcDyUbQA4gci2f4xsH4t52wDQOco2AByPF5Z9YIskjmnvCPO2AaBzlG0AOA7r0E5ZbpP8QK78/iNNx0k6bP8HAJ2jbAPAcTitJdIbeIpkO4bTJB83MrJ9YIvkuYbTAEByomwDwHFERmwjI7hozy8YKT+QKyt8RNahMtNxACApUbYB4Dgiu2x4AynbHbKdllF/sUgSAI6Hsg0AxxEZ2XYZ2T4uFkkCwIlRtgGgI0dqZNdWSJK8geMNh0le0cNtGNkGgA5RtgGgA5Et/7y+Q6WcfobTJC9GtgHgxCjbANCB6E4knBx5QtHt/+r2SY3VhtMAQPKhbANAB9p2ImEKyQll5ctr3YOcRZIAcCzKNgB0oG1xJCPbneFwGwA4Pso2AHyY58o+uLXlj2z71ykWSQLA8VG2AeBDrJqdssJH5Ady5BdwTHtnWCQJAMdH2QaAD4nO1x7AMe1d4Q1uLdsHt0tus+E0AJBcKNsA8CEOiyO7xe9bLD+7nywvJLuq1HQcAEgqlG0A+BCbbf+6x7KiB//YlZsMhwGA5ELZBoAPsQ+07kTCyZFd5jJvGwA6RNkGgKMdOSz78B5JbVvaoXMeO5IAQIco2wBwlPbHtPc3nCZ1eIMnSWqd7+77htMAQPKgbAPAUZwDkfnajGp3hzdgnHzLkdVYLau+0nQcAEgalG0AOEp0cSSH2XRPIFte0WhJkr2fRZIAEEHZBoCj2Gz712PewNZFkgdYJAkAEZRtAIjwXNkHWo5pd9n2r9tYJAkAx6JsA0Arq6ZcVrhRvpPNMe090G6RJABAEmUbAKLsytadSAaeItkBw2lST2Rk26reIYWOmA0DAEmCsg0ArZzW6Q8eh9n0iJ83UF5ukSzfk31wm+k4AJAUKNsA0Cq6xzbztXvGspi3DQAfQtkGgFaRgsge2z0XLdutH1wAINNRtgFAkppqZR/eLUlymUbSY5Gy7TCyDQCSKNsAIEnRLf+8PidJuQVmw6SwyGFAduUWjm0HAFG2AUDS0VNIGNXuDa/oZPl2QFZTjay6vabjAIBxlG0AkOSwODI2AlnyCsdIajuNEwAyGWUbAHTUyDbztXutbUcSyjYAULYBwPc4pj2G2nYkoWwDAGUbQMazanbJCjXId7LkF44yHSflRRZJcmw7AFC2AaBtCsmAcRzTHgMc2w4AbSjbADKeU8niyFjy8wfJzy3k2HYAEGUbANj2L9YsSy7ztgFAEmUbAKJHi3sDGdmOlbbDbSjbADIbZRtAZmuul11TLklyGdmOmciUHMo2gExH2QaQ0aKj2n2GSLmFhtOkj8iUHKeyhGPbAWQ0yjaAjBYZeeUwm9jyisbJt5zWY9v3mY4DAMZQtgFkNCe6OHKC4SRpJpAlryhybPtmw2EAwBzKNoCMFp1GwrZ/MRedt936HgNAJqJsA8hcvnfUTiRMI4m1yHvKIkkAmYyyDSBjWTW7ZTXXy3eC0SkPiJ3IyDbHtgPIZJRtABkrcuAKx7THR9ux7aVSuMlwGgAwg7INIGNFdyJhvnZccGw7AFC2AWQwh23/4suy5HKSJIAMR9kGkLGi00jY9i9uIofbsP0fgExF2QaQmZrrZR1qPaZ9IGU7Xtj+D0Cmo2wDyEj2ga2y5MvLHyTlFZmOk7YiU3Scys0c2w4gI1G2AWSktikkLI6MJ29A67HtR2pk1e03HQcAEo6yDSAj2SyOTIxAtryi0ZKYtw0gM1G2AWQkJ3pMO/O1480byLxtAJmLsg0g8/j+UXtsU7bjLfIeM7INIBNRtgFkHOvwHlnNdS3HtBdyTHu8Rcv2AfbaBpB5KNsAMk5khNUrGic5QcNp0p8XOdimagfHtgPIOJRtABmnbScSFkcmgt9nsPycAlm+y7HtADIOZRtAxoke0862f4lhWXIHcWw7gMxE2QaQcaLTSNj2L2GYtw0gU1G2AWSWo49pZ2Q7YaLzthnZBpBhKNsAMop9YEvrMe2DOaY9gSIj205lCce2A8golG0AGSU6hYT9tROq7dj2QxzbDiCjULYBZJTo4sjBkwwnyTBHH9vOvG0AGYSyDSCj2Ps3SWJk2wTmbQPIRJRtAJnDc2Uf2CqJxZEmeGz/ByADUbYBZAyrZqescKP8QI78glGm42Sc6Mj2gc2GkwBA4lC2AWQMZ/9R+2vbjuE0mSc6ss2x7QAyCGUbQMawoydHMl/bBL/PEPk5/Tm2HUBGoWwDyBiRbf+Yr22IZcmNTiXZYjgMACQGZRtAxmjbY5uybUrbIknmbQPIDJRtAJmhsVp23T5JTCMxqa1sM7INIDNQtgFkhOhhNv1HSln5htNkrshvFZzKzRzbDiAjULYBZASOaU8OXtHJ8i275dj2+krTcQAg7ijbADJC2+JIyrZRwRz5hWMkMW8bQGagbAPICNFt/wZPMpwE7qDxkjhJEkBmoGwDSH9us+yD2yUxjSQZROZtU7YBZALKNoC0Z1d9IMsLyc/uJ79vsek4Gc8b2DqyfYCyDSD9UbYBpD17/1GLIy3LcBpER7arSjm2HUDao2wDSHucHJlc/D5D5Ge3Httetd10HACIK8o2gLQXXRzJfO3kYFnRXWGYtw0g3VG2AaQ335dTuUkSx7QnE4+yDSBDULYBpDWrbp+sIzXyLUfegHGm46BVdJEke20DSHOUbQBpLXpyZNFYKZBtOA0iIvudO/s5th1AeqNsA0hrbce0M4UkmXgDTpFvB2Q11ciqrTAdBwDihrINIK2xODJJBbLkFZ0sSbL3bzIcBgDih7INIK05kZHtwYxsJ5u2qSSUbQDpi7INIH2FGmRVl0liGkkyinwAYpEkgHRG2QaQtuzKLbLky8sfKD9vgOk4+BB3UMvItr1/o+EkABA/lG0Aacs+EJmvPclwEnQkMo3Erq2QGqsNpwGA+KBsA0hbkbnALI5MUtl95fUfLklyONwGQJqibANIW207kTBfO1l50akkLJIEkJ4o2wDSk+/JPrBFkuRStpOWO5iyDSC9UbYBpCXr0E5ZoQb5Trb8wlGm4+A42nYkoWwDSE+UbQBpKTpfe+B4yQ4YToPj8QZNliTZVaVS6IjhNAAQe5RtAGkpsp2cN2Sy4SQ4Eb/PYHm5RbJ8V/bBrabjAEDMUbYBpKXIHGB3MGU7qVlW2xaAzNsGkIYo2wDSj+/LiYxsU7aTXmS3GIfDbQCkIco2gLRj1e2V1Vgt3w7IG3iK6TjoRHRkm2PbAaQhyjaAtGPvax3VLjpZCmQbToPORKb62JVbJM81nAYAYouyDSDtOCyOTCl+wUj5gVxZ4UZZ1TtMxwGAmKJsA0g7NvO1U4vtyBs0QZLkMJUEQJqhbANIO+xEknraTpJkkSSA9ELZBpBWrPoDsuv2yZcVHS1F8ovsSMIiSQDphrINIK1ERrX9ojFSVr7hNOiqyJQfZ/8myfcNpwGA2KFsA0grkWkI7qBJhpOgO7yBp8i3A7Iaq2XV7TUdBwBihrINIK2wE0mKCmTLGzBOkmTv22A4DADEDmUbQFphJ5LUFZ1KQtkGkEYo2wDSx5Ea2TW7JLXtboHU4Z40VRIj2wDSC2UbQNpw9rfsZOH1Hy7l9DecBt3lDZ4iqbVss0gSQJqgbANIG0whSW3eoAny7YDsxioWSQJIG5RtAGkjuhMJU0hSUyBb3oBTJEn23vWGwwBAbFC2AaQNh5HtlBfZRcbZz7xtAOmBsg0gPTTXy6oqlUTZTmXukMi8bY5tB5AeKNsA0oJdWSJLvrz8wfLzB5qOgx6KLJJ09q1nkSSAtEDZBpAWOMwmPUQWSVqN1bJqK0zHAYBeo2wDSAuRvZk9jmlPbUcvkmS/bQBpgLINIC3Y+96XJLknnWo4CXorukiSsg0gDVC2AaS+5nrZBz+QJHmtpxAidblDOEkSQPqgbANIefa+DS2LI/ucJD9/kOk46CWvdUcSZz8nSQJIfZRtACnP2ddyAAqj2unBGzieRZIA0gZlG0DKs/dG5mtPM5wEMdFukSQnSQJIbZRtACmPke30E51KwuE2AFIcZRtAamusll2zS5Lkth6IgtTXdpIkiyQBpDbKNoCU5uxtHdUuHC3l9DMbBjHTNrLNSZIAUhtlG0BKi+6vPYT9tdOJN3CCfDso68ghWbV7TMcBgB6jbANIadGRbeZrp5dAlryB4yRJ9l6mkgBIXZRtAKnL9zk5Mo15rb+tcPauM5wEAHqOsg0gZVl1+2TXH5BvOfIGTTQdBzHmDm3ZypGyDSCVUbYBpKzI/trewFOkYK7hNIg1r3XfdHvfBslzDacBgJ6hbANIWdH9tYcwXzsdeUVj5Wflywo1yD64zXQcAOgRyjaAlNV2ciTztdOS7UR3mbEr1hoOAwA9Q9kGkJp8X07rgSfsRJK+POZtA0hxlG0AKck6VCar6bB8J1vegFNMx0GcuJF52xWUbQCpibINICVF99cePFFygobTIF4iI9v2wW1Sc73hNADQfZRtACmJkyMzg58/SF7fobLkRz9gAUAqoWwDSEmcHJk53KGnSZJs5m0DSEGUbQCpx21u2XtZbXN6kb4i+2077EgCIAVRtgGkHLuyRJbbJD+nv/zC0abjIM4iJ0nae9dJvm84DQB0D2UbQMpx9qyRJLlDT5csy2gWxJ83eLJ8y5FdXymrbq/pOADQLZRtACnHrlgtSXKLTzcbBIkRzJU3aLwktgAEkHoo2wBSTmRk2xt6utEcSJy2eduUbQCphbINIKX4hytk11bIt2yOac8gkR1JnL0skgSQWijbAFKKv/PfkiRv4AQpK99wGiRK9CTJfRslL2w4DQB0HWUbQErxy1vLdvFphpMgkfyiMfKz+8oKN8o+sNV0HADoMso2gJTil6+UJLlDzzCcBAll2XKHtBxgxOE2AFIJZRtA6gg3y9/TUrTYiSTzeENZJAkg9VC2AaQMe/8GyW2Wl1skv/8I03GQYG7r7jPOntVmgwBAN1C2AaSMyHHdXvEZHGaTgdzilqlDdnWp/PoDhtMAQNdQtgGkjMiIZmQbOGSYnP5yB5wiSfLLVhgOAwBdQ9kGkBp8X3bkmHbma2csb9h0SZK/8z3DSQCgayjbAFKCVVshu36/ZAfkte5KgczjDpshSfLLlhtOAgBdQ9kGkBIiU0isoVOkYK7hNDDFLW4d2d6zTgo1Gk4DAJ2jbANICXbr4khrxJmGk8Akv1+xvD5DJC8sZ+/7puMAQKco2wBSgtM6X5uyneEsS27rvG1790rDYQCgc5RtAMkvdER25SZJkjWSsp3pvOKWedvO7lWGkwBA5yjbAJKevXedLC8sL3+Q1H+46TgwLDKy7VSsljzXcBoAODHKNoCk5+xq2VPZHX6mLA6zyXjewPFSdh9ZzfWyD2wxHQcAToiyDSDpObv+LUlyh59lOAmSgu3IGtFyLTCVBECyo2wDSG5us5yKNS1/pGyjlTXqbEmSvYdFkgCSG2UbQFKz966XFT4iP7dQftHJpuMgSUTKtrNrpeT7htMAwPFRtgEkteh87WFnSszXRitr2Bny7YDs+v2yDu8xHQcAjouyDSCpMV8bHbGy8uQNnixJcthvG0ASo2wDSF5eWM6elgVw7gjKNtpzh7Xut72HRZIAkhdlG0DSsvdtlBVqkJ/dv2W7N+AonCQJIBVQtgEkrbb9tWdIFj+u0J5b3Hq4zcFtUmO14TQA0DH+9gKQtNrKNlNI0IG8IrkDWnaoiVwrAJBsKNsAkpPnRhe+UbZxPO7wmZIkp5yyDSA5UbYBJCW7crOs5jr5WX3kDZpoOg6SlDuidb/t8vcMJwGAjlG2ASSltv21p0u2YzgNklXktx7Owa1SQ5XhNABwLMo2gKTE/trokrwiuQNOkcS8bQDJibINIPn4npzdlG10jTsiMm97ueEkAHAsyjaApGMf2CrrSI38YNspgcDxROdt72LeNoDkQ9kGkHSi87WLz5CcoOE0SHbu8DPly5JzcLus+gOm4wBAO5RtAEnH2blMUtuIJXBCuYXRHWvYlQRAsqFsA0guXljOrpa5t+7IcwyHQaqIXCtO2RLDSQCgPco2gKRi79sgq6lWfnY/5mujy9xRrWV751LJ9w2nAYA2lG0ASSUyMumOmMn+2ugyd9gM+U5Qdm2FrENlpuMAQBRlG0BSCexcKkkKj5ptOAlSSjBXbvF0Sa2j2wCQJCjbAJJHqEH2njWSmK+N7otcMwHmbQNIIpRtAEnD2bVSlheS169YfsFI03GQYqKLJMvfkzzXcBoAaEHZBpA0Ir/+d0eeI1mW4TRINd6QKfKz+8lqqpW9b73pOAAgibINIIk4OxZLklzma6MnbKft6Paydw2HAYAWlG0AScGq3Svn4Fb5lq0w87XRQ+HR50mSAq0f3ADANMo2gKTg7HhHkuSddKqUW2g4DVKVO+Z8SZJdsVZqPGQ2DACIsg0gSURGIsOjzzecBKnM7ztU7oBxsnwvuo0kAJhE2QZgnhtqO8yGso1eilxDDlNJACQByjYA4+yKtbKa6+TnFMgbMsV0HKQ4t3XetrPjHY5uB2AcZRuAcdEpJKPO5Yh29Jo77Ez5gVzZ9ZWyD5SYjgMgw1G2ARjn7PiXpLbFbUCvBLLkjmzdArD0HcNhAGQ6yjYAo6zavXL2b5IvK/rrf6C3wmMukCQFPnjbbBAAGY+yDcAop7UMeUOnyc8bYDYM0oY7Zo4kya5YIzVWmw0DIKNRtgEYFfjgLUlSeOxFhpMgnfj9iuUOmtiyBWDpItNxAGQwyjYAc0INcnYukyS5lG3EmDv2QklSYPvbRnMAyGyUbQDGOGVLZbnN8voVyxt4iuk4SDOR35Y4ZYslt9lwGgCZirINwJh2U0gsy3AapBvvpKny8gfKaq6Xs2uF6TgAMhRlG4AZviendS4tU0gQF5YdXSjpMJUEgCGUbQBG2BVrZdcfkJ+VL3f4WabjIE2FT75EkhTY/k9OkwRgBGUbgBGBrW9Kap1CEsgynAbpyh01W34wT3Zthex9603HAZCBKNsAEs/3Fdj2uiQpPO4yw2GQ1gLZbQfcbH3DcBgAmYiyDSDh7MrNsmt2yXey5Y7h1EjEl3tKywe6wNbXmUoCIOEo2wASLjLC6I4+TwrmGU6DdBceM0e+kyX7UJnsg1tNxwGQYSjbABLO2dZStsOnXG44CTJCVr7cUedKkhymkgBIMMo2gISyqj6Qc3CbfDugcOsJf0C8hY+eSgIACRQwHaArfN9XQ4Ovxkbm2mW67Gyug1SXt/5lSVLz8HPU6PeVuvnfk2sAUvevA6v4QmXbATkHtqh591a5RePimA6JwM8CFBSkxn//pC/bvu/rS3f6en99lekoSApcB6nN1/MXvaKxfaUf/v0KvfRQT35Qcg1A6v510F+/OPtczTlpkZ774Sv6TckdcUmFROJnQaY74/TDevDnvqwkP4GYaSQAEmZ8vxKN7fuBmtwsvb33YtNxkGFe232FJOkjw16RlBojYgBSn+X7XdsHqbq6Ot5Zjsv3feXkFOrQIXMZkBwKCrgOUln+sp8pb/VjahpzmQ5/9Bc9eg6uAUg9uw6s5noNePI8WW6Tquf/VeFBk+OUDonAzwIMHVqoQ4cOGc1QWFjY6WOSfhqJJFmWpbw8S01Nyf1rAsQf10EK833lbH9FkuRNuUK5uT3778g1AKmH10FuH7ljL1Rg62vK3/GKmkdOiU84JAQ/C5Ds00cimEYCICHsirWyD++WH8yTO2aO6TjIUKEJLVNJAiUvS75nOA2ATEDZBpAQwU0LJUnhky+RgrmG0yBTuWMukJ/VR3Zthexd/zYdB0AGoGwDiL9wswKbW6aQhKdcZTgMMlowR+HxH2n548aFhsMAyASUbQBx53zwtqymGnl9hsgdMct0HGS40ORPSpICW16VQg1mwwBIe5RtAHEX3Ph3SVJ40ick2zEbBhnPGzZDXv8RskINCmx903QcAGmOsg0grqyGg3J2vCNJCk1mCgmSgGVFr8VA6wdBAIgXyjaAuApsekGWF5Y75FT5A042HQeQJIVby7azc5ms2grDaQCkM8o2gPjxfQXf/6skKTT1U4bDAG38/sMVHn62LPkKrH/edBwAaYyyDSBu7N0rZVdtlx/IVXjilabjAO2Epy2QpJYPhF7YcBoA6YqyDSBugu8/K0kKT7xCyu5jOA3QXnjcZfJzC2XX7Y2uKwCAWKNsA4iPxmoFtrwmSQpNu8ZwGKADgSyFJrdMbwque9ZwGADpirINIC6CGxfKcpvlDp4kb8hU03GADoWmzZckOaX/YqEkgLigbAOIPd9TcO2fJEmhUxdIlmU4ENAxv3CMwiNmyvI9Bdf+xXQcAGmIsg0g5pzSRbIP7ZSf3U/hSSyMRHILnX69JCm47i9S6IjhNADSDWUbQMwFVz0jSQpNnS9l5RtOA5yYe/LF8voNk3XkkAKbXzAdB0CaoWwDiCn7wFYFdi6Vb9kKnXGd6ThA52ynbXR71TOS7xsOBCCdULYBxFRwdcuotnvyJfL7DTOcBuia0NR58oN5cg5ulVO+zHQcAGmEsg0gZqz6SgU2LpQkNU+/wXAaoBty+ik0+ZOSpOCKJ8xmAZBWKNsAYia46umW7f6GniZv2Jmm4wDdEprxOfmWrUDZYtn7NpiOAyBNULYBxMaRw9Ht/prPvoXt/pBy/IIRCk+4QpIUXPGY4TQA0gVlG0BMBNf+SVZzvdwBp8gde6HpOECPhM66WZIU2PKarOpSw2kApAPKNoDeCzUouOrplj+efbNk8aMFqckbNF7hsRfJkq+s5YxuA+g9/kYE0GvB1X+Q3Vglr/9IhSd8zHQcoFeaZ94qSQpsXMjoNoBeo2wD6J2mWmWteFyS1HzOlyU7YDgQ0Dve0NMUHnuhLN9V1pJfmY4DIMVRtgH0StbKJ2U11cgdcLLCE+eajgPERPPsr0iSAiUvy67cYjgNgFRG2QbQcw1VCq56SpLUfM5XJNsxHAiIDW/wJIXGf7Rl7vaSX5iOAyCFUbYB9Fj20l+27EAyeLLcUy41HQeIqebZd8q3HAW2/6+c8uWm4wBIUZRtAD1iH9iqwLpnJUlNc77JDiRIO37RWIWmLZAkZS36L8lzDScCkIr42xFAj2Qtul+W7yk87lJ5I842HQeIi+bZd8rP7itn/yYFNi40HQdACqJsA+g254O3FChbLN8Oqun8b5iOA8RPbqGaZ94uScp69+dSU63hQABSDWUbQPeEGpT9zx+1/HH6DfILRxkOBMRX6Izr5RWOll1/QFmLHzAdB0CKoWwD6Jaspb+SXVshr19xy77aQLpzstR0yfclScG1f5Jdsc5sHgAphbINoMvsys0KrmzZ6q/p4nulYJ7hREBiuCNnKjT5Klnylf3m9yQ3ZDoSgBRB2QbQNeFmZb9yjyzfVfiUy+WOvch0IiChmi74P/Jz+sup3Kzg8kdMxwGQIijbALoka9mv5BwokZ9bqKaLv2s6DpB4eUXRaz9r2W9k733fcCAAqYCyDaBT9p7VCq54TJJ05NIfyM8faDgRYEZ44lyFJnxMlu8q55VvSqFG05EAJDnKNoATa6xWzktfl+V7Ck2+Su4pl5lOBBjVdPH/lZc/SHZ1qbL/9z7TcQAkOco2gOPzPeW88s2W3UcKRqrpou+YTgSYl1ugpo/dL9+yFdzwvALrnzedCEASo2wDOK7gew8rsOMd+U62jlz5oJTd13QkICm4I2ep+Zw7JEnZ//yh7P2bDCcCkKwo2wA65Gx5TdlLHpQkNV3yXXmDJhhOBCSX0MxbFR59viy3STkLvyyrbr/pSACSEGUbwDHsinUti78kNZ9xg8JT5xlOBCQhy9aRK/4/eYVjZNdWKGfhl1kwCeAYlG0A7VjVpcpZ+CVZbpPCYy5Q85xvmo4EJK+c/mr81G/l5xTI2bdeOS99jQNvALRD2QYQZdXsVu5zn5fdcFDuoEk6Mvdnku2YjgUkNb9gpBqveki+k63AB28r+5X/I3mu6VgAkgRlG4AkyTq8W7l//bzsur3yisaqcd5jUla+6VhASvCGzdCRTzwo3w4quOVVZb/2bckLm44FIAlQtgHIqvpAuX/+jOyanfL6D1fjvCekvCLTsYCU4o65QEfm/v/yLUfBTf9Qzotfk8LNpmMBMIyyDWQ4e88a5f7lhrYR7QXPyO87xHQsICW5p1yuI1c+IN8JKrDtDeX8/VbpyGHTsQAYRNkGMlhg4z+U+9yNshur5A6erIYFz8jve5LpWEBKc8ddqiOfelh+ME+BncuU96drZVWXmo4FwBDKNpCJws3KeuvHynn1m7LcZoVPvkSNC55m6ggQI+7Ic9R4ze/l9R0qu7pUeX+8Rs7W103HAmAAZRvIMC3zs69V1upnJEnNZ9+iI594kMWQQIx5gyep8bq/yC0+Q1ZTrXJf+Kqy/vlDKdRgOhqABKJsA5nCDSm4/BHlPfMpOfs3yc8pUOMnf6Pm8+6WLH4UAPHg5w9S46efUvNZX5QkZa39k/Ke/qScncsMJwOQKAHTAQDEn7PjXWUt+i85B7dJksKjzlPT5fexEBJIBCeo5vO/LnfELGW/8V3ZNeXK/etNCk24Qs3nf11+v2LTCQHEEWUbSGN2xVplLf2VAjvekST5Of3VdOG3FJ70CcmyDKcDMos7+lw1fPYfylr8MwXX/lnBkpcV2PZPhU6/XqEzb5KfP9B0RABxQNkG0o3vydmxWMFVTylQtqTlJjug0OnXqXnWl6Sc/oYDAhksu4+aL/m/Ck+dp+xFP5Wza4WyVj6h4No/KjR1nkJnfEZ+4WjTKQHEEGUbSBNW9Q4FNr+k4KYXZB8qkyT5lqPw5KvUPPNW+QUjDScEEOENmaLGTz8lp/Rfylr2azl71ylrzR+UteYPCo+arfDEjys87lIpu6/pqAB6ibINpCrfbynYH7ytQMlLcvZtaLsrq0/bKFn/4QZDAjguy5I7do4ax1wgZ+dSBVc9rUDpIgXKlihQtkT+m9+XO+YChSd8TOGR50i5haYTA+gByjaQKnxP1qGdcnavklO+TM7O92TX72+723Lkjpqt8MQrFB53GVv5AanCsuSOmi131GxZh8oV2Pyigptfkl21XYFtbyqw7U35suQNmih35Ey5I2bKO+lU+XkDTCcH0AWUbSDZuCFZtRWya8pbyvWBLbL3b5Z9YIusD+3P6ztZcovPUPiUyxUe/1EOpQFSnF8wQqFZtys08zbZB0oU2PySnA/elnNwm5zKTXIqN0krn5QkefmD5A2aKG/wJHlFJ8srGCG//3D5eQNZAA0kEco2EC++L4WbZIXqpeaGlqLcXC8r1CDrSI2shoOyGqpa/tnY+s/6Slm1e2X5XsdP6WTLGzypdXRrltyhp0vBnMR+XwDiz2oZyW4eNFE6/+uy6ivllC+XU/6enF0rZFWXya6vlF1fKbXuNhThB3Ll9yuWlzdAfl6R/Nyiln+2/ruy+8oP5snPypeO+qecoKFvFkhvlu/7flceWF1dHe8sx2Xv36j89/+k5iP1x97ZtfjtWCf8muPcF/Ov6eZznfD5EvU1J9CT1+nB+xYMBBQKh7r3Mj14HetEX+OFW/7nhmR5YckLSW5I8sKyWv/Z8u+h45bmzvhOtryC4fL7DZc3YFxLwR40UX7hKMnO7M/IhYWFRn8eITlk/HXQXC/7wBbZlZvl7N8s61BZy2/DTvBhvTO+HWz5+eIE5duB6J9lB+Q7QckORv9dtiNfVusIuiVZaj0cK/LvrSPrltV2e+Sxan1s67/7PRyFz8rKUnNzc4++Nmlk5G8gYvc9ZxcN1+EZtxidNllY2PlaipT4Wzv4/nPy1/5VfOaGrxS5aD/ED+QeNYqUJ2X3bR11GtA24pRbJD9/oPx+w+TnD8rQH8IAuiQrX17xGfKKz1D46NvdZlmHK2TX7mn9zVmVrMaDbX9uOCg117X8hq25Xgo1tAwSSLK8lkEChRtjWIfix5foBRnOk+QMP1fuiLNNRzmhlOgtzefcoeyRp6ux9nDHDzhhKTnOfSf4Ev+4X9OT1+nB15zIcZ/vBM+VsK/p9h3Hv+84r5+fn6/6+oYO7+tROe3Bfx8/OtoTbBsBajf603a/n5UvBXIk2+l+NgDoLidLfuEouYWjuv41brMUapQVajzqt3Ztv7GL/hbv6H/3w22/HfR9SX7rP73WXyb67W+Pjra33W75Rz22B3LzctXY0Nijr00KPfrtd6qL7fecN3i03OIZMX3OeEiJsu3nDZBz9ucUyuRfGUKSZBcWKsx1AACx42S1lPSjDrxKhRrYp7CQXpDh7MJCKQWuAdt0AAAAACBdUbYBAACAOKFsAwAAAHFC2QYAAADihLINAAAAxAllGwAAAIgTyjYAAAAQJ5RtAAAAIE4o2wAAAECcULYBAACAOKFsAwAAAHFC2QYAAADihLINAAAAxAllGwAAAIgTyjYAAAAQJ5RtAAAAIE4o2wAAAECcWL7v+6ZDAAAAAOmIkW0AAAAgTijbAAAAQJxQtgEAAIA4oWwDAAAAcULZBgAAAOIkYOqF161bp1/+8pdas2aNQqGQxo0bpxtvvFFXXnlll5/D8zz98Y9/1F/+8heVlZUpLy9PM2fO1N13363Ro0fHLzxiorfXwL///W+9+eabWr58uXbv3q2GhgYNGzZMl1xyiW699Vb169cvzt8BYiEWPwuOFgqFNH/+fG3evFljxozRq6++GuPEiLVYXQN1dXV64okn9Prrr6u8vFzBYFAjRozQJZdcojvuuCNO6RErsbgODh8+rN/97nd68803tWvXLmVlZWn48OH61Kc+pU9/+tPKzs6O43eA3li4cKFWrlyp9evXa8uWLQqFQvrJT36iq6++ulvPk4zd0MjWf++9956+8IUvKBgMau7cuerbt69ef/117dq1S3fffbduu+22Lj3Pd7/7XT377LMaN26c5syZo4MHD+rll19Wdna2/vznP2vcuHFx/k7QU7G4Bs4991xVV1drxowZmjRpkizL0vLly7Vx40aNHDlSf/7znzVgwIAEfDfoqVj9LDjaL37xCz355JNqaGigbKeAWF0De/bs0Y033qjy8nLNnj1bkyZNUnNzs3bu3Kk9e/bohRdeiPN3gt6IxXVw+PBhXX311SovL9eMGTN02mmnqbm5Wf/617+0c+dOzZo1S7/73e9k2/xSPxldfPHF2r17twoLC5WXl6fdu3f3qGwnZTf0EywUCvmXXnqpP3XqVH/Dhg3R22tra/25c+f6kydP9ktLSzt9nqVLl/rjx4/3r7vuOr+pqSl6+5IlS/wJEyb4119/fTziIwZidQ08/PDD/r59+9rd5nme/73vfc8fP368//3vfz/W0RFDsboOjrZ+/Xp/8uTJ/tNPP+2PHz/e/8hHPhLj1IilWF0D4XDYnzdvnj9t2jR/6dKlHb4OklesroNHHnnEHz9+vP/jH/+43e1NTU3+vHnz/PHjx/vLly+PdXzEyLvvvuvv2rXL9/2Wv9/Hjx/v/8///E+3niNZu2HCP94tW7ZMO3fu1Mc//nFNnjw5enufPn30pS99SeFwWM8//3ynz/Pcc89Jku666y5lZWVFbz/nnHN03nnnacWKFSotLY39N4Bei9U1cMstt2jw4MHtbrMsS1/60pckSStWrIhtcMRUrK6DiObmZt1zzz067bTT9JnPfCYekRFjsboGXnvtNb3//vv6/Oc/r1mzZh1zfyBgbMYkuiBW10F5ebkkac6cOe1uz8rK0rnnnitJOnjwYAyTI5Zmz56tYcOG9eo5krUbJrxsL1++XJJ03nnnHXNf5P8MkcecyHvvvae8vDxNnz79mPsiz03ZSk6xugaOJ/IXq+M4PX4OxF+sr4OHHnpIZWVl+s///E9ZlhWbkIirWF0DL7/8siTpox/9qCoqKvSnP/1JjzzyiF555RXV19fHMDHiIVbXwSmnnCJJeuedd9rdHgqFtGTJEuXk5OiMM87obVwksWTthgn/uL9jxw5J0qhRo465r3///iosLFRZWdkJn6OhoUGVlZUaP358h4UqMgE+8lpILrG4Bk7kf/7nfyS1/ZBGcorldbBu3To99thjuvvuuzVmzJhYxkQcxeoaWL9+vSRp5cqV+slPfqLm5ubofUVFRXrggQc0c+bM2IRGzMXqOvj0pz+thQsX6oknntD69es1depUhUIhvfPOO6qpqdF///d/a8iQIbGOjySRzN0w4SPbdXV1kqS+fft2eH+fPn1UW1t7wueI3N+nT5/jPsfRr4XkEotr4Hg2bdqkX/3qVxowYIC++MUv9jgj4i9W10Fzc7O+9a1vadKkSfr85z8f04yIr1hdA5GpAffdd59uvPFGLVq0SEuXLtW9996r2tpaffnLX9b+/ftjFxwxFavrICcnR88884w+8YlPaPny5XriiSf0zDPPRKeodDTaifSRzN2QJblIG+Xl5br11lvluq5+9rOfqaioyHQkJMADDzygsrIy/fjHP2bqUIbyWzfVuvDCC/WNb3xDJ510koqKinTDDTfoc5/7nGpra/XXv/7VcErEW1VVlW666SatXbtWjzzyiP7973/r3Xff1Q9+8AM9//zzWrBggWpqakzHRAZKeNmOfLI43qfUurq64366jYjcf7xPJ5Hbj/fpBmbF4hr4sN27d+vGG29UVVWVHnzwwQ4XSSG5xOI62LBhg5588knddtttmjBhQswzIr5i9bMg8jwXX3zxMfdddNFFktqmmiD5xOo6+K//+i+tXr1aDz74oObMmaO+fftq4MCBWrBggf7jP/5D5eXleuqpp2KaHckjmbthwst2ZM5MR/OvampqVF1d3eG8raPl5eVp0KBB2rVrl1zXPeb+yHwcDrZJTrG4Bo62a9cu3XDDDdq/f78eeOCB6F+uSG6xuA5KSkrkuq5++ctfasKECe3+J0mlpaWaMGGCzjzzzJjnR+/F6mdBZJ5+RwdZRW5ramrqRVLEU6yug0WLFqmgoEATJ0485r7IAMyGDRt6FxZJK5m7YcLL9llnnSVJWrx48TH3vfvuu5Kks88+u9PnOfvss9XQ0KBVq1Ydc1/kuSOvheQSq2tAainan/3sZ7V//379/Oc/16WXXhq7oIirWFwHo0eP1vz58zv8n9Qy0jF//nx98pOfjG14xESsfhZEitS2bduOuS9yW2+3FEP8xOo6aG5uVl1dXbsFshFVVVWS1G47OKSfpO2Gid7YOxQK+Zdccok/depUf+PGjdHbj968/oMPPojefvDgQX/btm3+wYMH2z1Psm5cjs7F6hooLy/3L7roIn/y5Mn+a6+9lrD8iI1YXQfHw6E2yS9W18DOnTv9qVOn+uecc46/d+/eds9z1VVX+ePHj/eXLFkS/28IPRKr6+Dzn/+8P378eP/nP/95u9ubmpqi9z3zzDNx/V4QG50dapNq3dDIce3Lli3TF7/4RQWDQX384x9Xnz59osey3nXXXbr99tujj/3lL3+phx56SHfccYfuvPPOds9z77336rnnnkuuIznRJbG4BiJHu55++ukd7s8q6ZhrBsklVj8LOjJhwgSOa08BsboGnnnmGd13330qKCjQZZddpqysLL399tvavXu3rrnmGv3whz9M9LeGbojFdbBp0yZdf/31qq+v17Rp0zR9+nQ1NTVp8eLFKi8v15QpU/SnP/1J2dnZJr5FdOK5557TypUrJUlbtmzRhg0bNH369OgUoksvvTT62+tU64ZGjtWaNWuW/vjHP+rBBx/UK6+8olAopHHjxumrX/2qPvGJT3T5eX74wx9qwoQJ+stf/qJnnnlGeXl5uuiii9hrNwXE4hrYvXu3JGnNmjVas2ZNh4+hbCe3WP0sQOqK1TVwww03aNiwYXr88cf10ksvyXVdjRs3TrfddpsWLFgQx+8AsRCL62DSpEl6/vnn9fDDD2vZsmX6wx/+IMdxNHLkSN155536whe+QNFOYitXrtTf/va3dretWrUqOiVk2LBhXZoqmozd0MjINgAAAJAJ2GcbAAAAiBPKNgAAABAnlG0AAAAgTijbAAAAQJxQtgEAAIA4oWwDAAAAcULZBgAAAOKEsg0AAADECWUbAAAAiBPKNgCkkJdeeknXXXedpk+frrPOOktf/vKXtXPnTtOxAADHwXHtAJACwuGwvvnNb+rFF1/UyJEjNXXqVJWVlWnDhg0aNGiQXnrpJfXv3990TADAhzCyDQAp4Mc//rFefPFF3XnnnXrttdf085//XM8//7zmz5+vyspKPfPMM6YjAgA6QNkGgCS3fPly/eEPf9DVV1+tO+64Q7bd9qP7s5/9rCTp3XffNRUPAHAClG0ASHIPPviggsGg7rrrrmPuGzBggCRpz549CU4FAOgKyjYAJLHS0lKtWLFCF198sYYMGXLM/Y2NjQZSAQC6KmA6AADg+F5//XVJUkVFhe65555j7j906JAkqV+/fomMBQDoIso2ACSx5cuXS5LWrVundevWHfdxo0aNSlQkAEA3ULYBIIlt2LBBubm5WrNmTYf3f/3rX9eLL76oKVOmJDYYAKBLmLMNAEnq8OHDqq6u1qBBgzq83/f96Mj32WefnchoAIAuomwDQJKqqamRJOXn53d4/7p167R//34NHjxYZ5xxRiKjAQC6iLINAEnKsixJUigU6vD+Z599VpI0f/78dntvAwCSBz+dASBJDRkyRIFAQLt27VJzc3O7+7Zt26aFCxeqoKBAn/vc58wEBAB0irINAEkqGAxq5syZOnLkiJ566qno7Xv27NEdd9yhUCikH/zgB+rfv7/BlACAE7F83/dNhwAAdGzlypW64YYb5LquZs6cqT59+mjJkiVqamrSd77zHX3mM58xHREAcAKUbQBIcosXL9aDDz6okpIS5eTk6Mwzz9Stt96qadOmmY4GAOgEZRsAAACIE+ZsAwAAAHFC2QYAAADihLINAAAAxAllGwAAAIgTyjYAAAAQJ5RtAAAAIE4o2wAAAECcULYBAACAOKFsAwAAAHFC2QYAAADihLINAAAAxAllGwAAAIiT/weO5PmKFzEPIgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 720x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for a, b in priors:\n",
    "    distri = beta(a, b)\n",
    "    x = np.linspace(0, 1, 300)\n",
    "    x_pdf = distri.pdf(x)\n",
    "    plt.plot(x, x_pdf, label=rf\"$\\alpha$ = {a:d}, $\\beta$ = {b:d}\")\n",
    "    plt.yticks([])\n",
    "    plt.xlabel(\"$\\\\theta$\")\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell returns the Bayes factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5\n"
     ]
    }
   ],
   "source": [
    "BF = beta_binom(priors[1], y) / beta_binom(priors[0], y)\n",
    "print(round(BF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the model with the more concentrated prior $\\text{beta}(30, 30)$ has $\\approx 5$ times more support than the model with the more extended prior $\\text{beta}(1, 1)$. Besides the exact numerical value this should not be surprising since the prior for the most favoured model is concentrated around $\\theta = 0.5$ and the data $y$ has equal number of head and tails, consintent with a value of $\\theta$ around 0.5."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequential Monte Carlo\n",
    "\n",
    "The [Sequential Monte Carlo](SMC2_gaussians.ipynb) sampler is a method that basically progresses by a series of successive *annealed* sequences from the prior to the posterior. A nice by-product of this process is that we get an estimation of the marginal likelihood. Actually for numerical reasons the returned value is the log marginal likelihood (this helps to avoid underflow)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Initializing SMC sampler...\n",
      "Sampling 4 chains in 4 jobs\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='100' class='' max='100' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [100/100 00:00<00:00  Stage: 2 Beta: 1.000]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/osvaldo/anaconda3/lib/python3.9/site-packages/arviz/data/base.py:220: UserWarning: More chains (4) than draws (3). Passed array should have shape (chains, draws, *shape)\n",
      "  warnings.warn(\n",
      "Initializing SMC sampler...\n",
      "Sampling 4 chains in 4 jobs\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='100' class='' max='100' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [100/100 00:00<00:00  Stage: 0 Beta: 1.000]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/osvaldo/anaconda3/lib/python3.9/site-packages/arviz/data/base.py:220: UserWarning: More chains (4) than draws (1). Passed array should have shape (chains, draws, *shape)\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "models = []\n",
    "idatas = []\n",
    "for alpha, beta in priors:\n",
    "    with pm.Model() as model:\n",
    "        a = pm.Beta(\"a\", alpha, beta)\n",
    "        yl = pm.Bernoulli(\"yl\", a, observed=y)\n",
    "        idata = pm.sample_smc(random_seed=42)\n",
    "        models.append(model)\n",
    "        idatas.append(idata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block !important;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.DataArray &#x27;log_marginal_likelihood&#x27; ()&gt;\n",
       "array(5.)</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'>'log_marginal_likelihood'</div></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-48256fa5-6b5d-4e8f-bd2a-dae1b1a59ec4' class='xr-array-in' type='checkbox' checked><label for='section-48256fa5-6b5d-4e8f-bd2a-dae1b1a59ec4' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>5.0</span></div><div class='xr-array-data'><pre>array(5.)</pre></div></div></li><li class='xr-section-item'><input id='section-3a1b2b71-5877-4243-a7d1-28311a0154be' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-3a1b2b71-5877-4243-a7d1-28311a0154be' class='xr-section-summary'  title='Expand/collapse section'>Coordinates: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'></ul></div></li><li class='xr-section-item'><input id='section-cbd49c3b-648e-4f41-a218-c2488c73baca' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-cbd49c3b-648e-4f41-a218-c2488c73baca' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.DataArray 'log_marginal_likelihood' ()>\n",
       "array(5.)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BF_smc = np.exp(\n",
    "    idatas[1].sample_stats[\"log_marginal_likelihood\"].mean()\n",
    "    - idatas[0].sample_stats[\"log_marginal_likelihood\"].mean()\n",
    ")\n",
    "np.round(BF_smc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see from the previous cell, SMC gives essentially the same answer as the analytical calculation! \n",
    "\n",
    "Note: In the cell above we compute a difference (instead of a division) because we are on the log-scale, for the same reason we take the exponential before returning the result. Finally, the reason we compute the mean, is because we get one value log marginal likelihood value per chain. \n",
    "\n",
    "The advantage of using SMC to compute the (log) marginal likelihood is that we can use it for a wider range of models as a closed-form expression is no longer needed. The cost we pay for this flexibility is a more expensive computation. Notice that SMC (with an independent Metropolis kernel as implemented in PyMC) is not as efficient or robust as gradient-based samplers like NUTS. As the dimensionality of the problem increases a more accurate estimation of the posterior and the _marginal likelihood_ will require a larger number of `draws`, rank-plots can be of help to diagnose sampling problems with SMC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes factors and inference\n",
    "\n",
    "In this example we have used Bayes factors to judge which model seems to be better at explaining the data, and we get that one of the models is $\\approx 5$ _better_ than the other. \n",
    "\n",
    "But what about the posterior we get from these models? How different they are?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>a</th>\n",
       "      <td>0.5</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.4</td>\n",
       "      <td>0.59</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   mean    sd  hdi_3%  hdi_97%\n",
       "a   0.5  0.05     0.4     0.59"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(idatas[0], var_names=\"a\", kind=\"stats\").round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>a</th>\n",
       "      <td>0.5</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.42</td>\n",
       "      <td>0.57</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   mean    sd  hdi_3%  hdi_97%\n",
       "a   0.5  0.04    0.42     0.57"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(idatas[1], var_names=\"a\", kind=\"stats\").round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may argue that the results are pretty similar, we have the same mean value for $\\theta$, and a slightly wider posterior for `model_0`, as expected since this model has a wider prior. We can also check the posterior predictive distribution to see how similar they are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:00<00:00]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:00<00:00]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ppc_0 = pm.sample_posterior_predictive(idatas[0], model=models[0]).posterior_predictive\n",
    "ppc_1 = pm.sample_posterior_predictive(idatas[1], model=models[1]).posterior_predictive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA48AAAJjCAYAAACsmCRCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAApIElEQVR4nO3dfZTWdZ3/8Rc6Q7CDErreQdiNOhS4BmqihmxR631GyXqzCipkYoLm6WbtbO3+1tPaqdQKLewk2aaVyubJFGKozlIqcrMkIbSonAwQdTEhvAlkhrl+f+wyiTP4QblmLhgfj3M6J77fr995X/M518w85/u95upRqVQqAQAAgFexR60HAAAAYNcnHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoKhuRw9cv359Z85BB/r27ZsNGzbUegyqzLp2P9a0e7Ku3Y817Z6sa/djTWujX79+xWNcedyF7bGH5emOrGv3Y027J+va/VjT7sm6dj/WdNdlZQAAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoqqv1AAAAQMem3dJa6xG2a8JF3fM61KJFi3LZZZdlwoQJufjii1/3eUaPHp0k+clPflKdwXYB4hEAAGAX09ramh//+Mf5yU9+kieeeCK9e/fOUUcdlUsuuSQHH3xwTWbqnr8uAAAA2I19+ctfznXXXZfW1taMGTMmxx9/fO67776MHz8+jz/+eE1mcuURAABgF7Jo0aLcfffdGTp0aKZMmZKePXsmSU455ZRcfvnl+cpXvpKpU6d2+VziEQAAqImXv75w+PDhuemmm7J8+fL07Nkzo0aNyuWXX55evXrlwQcfzHe/+9089thj6d27d0477bRceuml2XPPPdvO1dLSkunTp2fGjBlZvXp16urqMmjQoJx33nl573vf2+5jb9q0KdOmTUtTU1P+9Kc/5S1veUvOOuusDBw4cLvzPvnkk7nllluyYMGCrFu3LnvvvXeGDx+eiy++OAcddFDVPi933313kuSSSy5pC8ckec973pPhw4dn3rx5WbVqVZffvuq2VQAAoKaWLVuWyZMnp0+fPjnnnHNywAEH5K677sqXvvSl/PKXv8xVV12VAw44IKNHj85ee+2V2267Ld///vfb/vtKpZLPf/7z+cY3vpHNmzfnzDPPzIknnpgVK1bkU5/6VO64445tPl5ra2s+85nP5NZbb81ee+2Vs846K0OGDMk3vvGN/PCHP+xwxqVLl2bcuHGZOXNm3vnOd+ass87K0KFD09TUlPHjx2fNmjVV+3z85je/Se/evXPEEUe023fssce2HdPVXHkEAABqat68efnKV76SkSNHpl+/fnnmmWdy4YUXZvbs2Zk/f36mTp2awYMHJ0kuvvjijBkzJrfffnvGjh2burq6zJo1K3PmzMmwYcMyZcqU1NfXJ0kuuuiiXHDBBbnhhhtywgknpH///kmSmTNnZuHChTn22GNz3XXXtV3BPPvss3PRRRe1m6+lpSVf+MIX0tramu9973s57LDD2vYtXrw4l112Wa6//vpcd911O/252LhxY/74xz/mkEMO2ebK6lZbr4yuXr16pz/WayUeAdhpPefeWOsRqm5L717puXFTrceoqs3HT6r1CAAdOvLIIzNy5Mi2f9fV1WXUqFFZsWJFRowY0RaOSdLQ0JARI0bknnvuydq1a9O/f//MmDEjSTJp0qS2cEyS/fffP+eee26+9a1vZdasWRk/fnyS5Gc/+1mSZOLEidsE2qGHHpqTTz4599xzzzbz3X///Xnqqafy8Y9/fJtwTJKhQ4fmhBNOyK9//eu8+OKLaWho2KnPxQsvvND2ODuydfuLL764Ux/n9RCPAABATTU2Nrbbtu+++yZJu1h7+b5nnnkm/fv3z6OPPpo3velNGTJkSLtjjzzyyCTJY4891rbtscceS69evfLOd76z3fFDhw5tF49Lly5NkqxcuTLf+c532v03zz77bFpbW7Nq1aq8613v2u7j3N2JRwAAoKY6usq29Yrgq+1raWlJ8r9X4fbff/8Oz73PPvsk+csVvR09/uWee+65JElTU9N2H0Pyv7ec7qw+ffok2f6Vxa3bd/YK5+shHgEAgN1aQ0ND1q1b1+G+rdtfHlsNDQ1Zv379qx7/yvMnybXXXpsRI0bs7Livqnfv3vnrv/7rPPnkk9myZUu71z1ufa3jq/1V2M7ir60CAAC7tcbGxrz00ktZtmxZu30PPfRQ2zFbHXbYYdm0aVOWL1/e7vjFixe327b1dtiHH364ShO/umHDhmXjxo1ZsmRJu33z5s1L8pfbcbuSeAQAAHZrp556apJk6tSpbbeyJsnatWvzox/9KHvuuWdOOumktu2nnHJKkuSmm27Kli1b2ravWLEis2bNanf+kSNH5sADD8ztt9/eFqMv19LS0mF0vl6jR49Oknz7299Oc3Nz2/aFCxdm/vz5GTZsWJe/x2PitlUAAGA3d8opp2TOnDn59a9/nfPOOy8jRozIxo0b88tf/jIbNmzI5ZdfngEDBrQdf+qpp6apqSnz5s3LuHHjctxxx+W5557Lz3/+8xxzzDF54IEHtjl/z549c8011+TKK6/MpZdemqOPPjqHHHJIkuTpp5/Ob3/72/Tt27fd+0m+XkcddVTOOOOM/PSnP824ceNy/PHHZ/369fnFL36RhoaGfPazn63Kx3mtxCMAAOyiJlzkRsEd0aNHj1xzzTW54447MnPmzEyfPj11dXUZNGhQzj333G3eBiRJ9thjj3z1q1/NzTffnKamptx5550ZMGBArrjiigwcOLBdPCbJ4MGDc+utt+YHP/hB5s6dmyVLlqS+vj777bdfRo4cmRNPPLGqj+mqq67KoYcemp/85CeZPn16evfunREjRmTixIk1ueqYJD0qlUplRw7c3gtK6Tz9+vXzee+GrGv3Y0275/s89urdK5u8z2O34rnaPVnX7sea1ka/fv2Kx/hVBgAAAEXiEQAAgCKveQQAAOhEjz76aH71q18VjzvooINy+umnd8FEr494BAAA6ESPPvpopk2bVjxu2LBh4hEAAOCN6vTTT9+lo3BHec0jAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAMD/WbRoUY499th85zvf2anzjB49OqNHj67OULuIuloPAAAAdKzn3BtrPcJ2bT5+Uq1H6LbWrVuXe+65J8uXL8/y5cvz1FNPJUnmzZtX07nEIwAAwC7k8ccfz9SpU9OjR48MHDgwvXr1yqZNm2o9lngEAADYlbztbW/L1KlT09jYmIaGhpx99tlZuXJlrccSjwAAQG0sWrQol112WSZMmJDhw4fnpptuyvLly9OzZ8+MGjUql19+eXr16pUHH3ww3/3ud/PYY4+ld+/eOe2003LppZdmzz33bDtXS0tLpk+fnhkzZmT16tWpq6vLoEGDct555+W9731vu4+9adOmTJs2LU1NTfnTn/6Ut7zlLTnrrLMycODA7c775JNP5pZbbsmCBQuybt267L333hk+fHguvvjiHHTQQVX7vOy7777Zd999q3a+ahGPAABATS1btiy33XZbhg8fnnPOOScPPPBA7rrrrrz44osZOXJkrr766pxwwgkZMmRI5s6dm9tuuy0NDQ256KKLkiSVSiWf//znM2fOnBx88ME588wzs3Hjxvzyl7/Mpz71qVx55ZU5++yz2z5ea2trPvOZz2ThwoU55JBDcuKJJ2bDhg35xje+kSOPPLLDGZcuXZpPfvKT2bhxY0aMGJG3vOUtefrpp9PU1JQHH3wwN998cwYMGNAln69aEY8AAEBNzZs3L1/5ylcycuTI9OvXL88880wuvPDCzJ49O/Pnz8/UqVMzePDgJMnFF1+cMWPG5Pbbb8/YsWNTV1eXWbNmZc6cORk2bFimTJmS+vr6JMlFF12UCy64IDfccENOOOGE9O/fP0kyc+bMLFy4MMcee2yuu+66tiuYZ599dluQvlxLS0u+8IUvpLW1Nd/73vdy2GGHte1bvHhxLrvsslx//fW57rrrOvtTVVPeqgMAAKipI488MiNHjmz7d11dXUaNGpVKpZIRI0a0hWOSNDQ0ZMSIEdmwYUPWrl2bJJkxY0aSZNKkSW3hmCT7779/zj333LS0tGTWrFlt23/2s58lSSZOnLjNra+HHnpoTj755Hbz3X///Xnqqady/vnnbxOOSTJ06NCccMIJefDBB/Piiy/uzKdhl+fKIwAAUFONjY3ttm19zd8rY+3l+5555pn0798/jz76aN70pjdlyJAh7Y7dehvqY4891rbtscceS69evfLOd76z3fFDhw7NPffcs822pUuXJklWrlzZ4fs/Pvvss2ltbc2qVavyrne9a7uPc3cnHgEAgJpqaGhot23rFcFX29fS0pIkefHFF7P//vt3eO599tknSfLCCy+0bduR41/uueeeS5I0NTVt9zEkycaNG191/+5OPAIAALu1hoaGrFu3rsN9W7e/PEIbGhqyfv36Vz3+ledPkmuvvTYjRozY2XF3W17zCAAA7NYaGxvz0ksvZdmyZe32PfTQQ23HbHXYYYdl06ZNWb58ebvjFy9e3G7b1tthH3744SpNvHsSjwAAwG7t1FNPTZJMnTq17VbWJFm7dm1+9KMfZc8998xJJ53Utv2UU05Jktx0003ZsmVL2/YVK1Zs84d1tho5cmQOPPDA3H777W0x+nItLS0dRmd347ZVAABgt3bKKadkzpw5+fWvf53zzjsvI0aMaHufxw0bNuTyyy/f5j0YTz311DQ1NWXevHkZN25cjjvuuDz33HP5+c9/nmOOOSYPPPDANufv2bNnrrnmmlx55ZW59NJLc/TRR+eQQw5Jkjz99NP57W9/m759++aOO+6o2mO6+uqr2/7/H//4x3bbLr/88rz5zW+u2sfbEeIRAADYrfXo0SPXXHNN7rjjjsycOTPTp09PXV1dBg0alHPPPXebtwFJkj322CNf/epXc/PNN6epqSl33nlnBgwYkCuuuCIDBw5sF49JMnjw4Nx66635wQ9+kLlz52bJkiWpr6/Pfvvtl5EjR+bEE0+s6mOaOXPmq2772Mc+1uXx2KNSqVR25MDtvaCUztOvXz+f927IunY/1jTpOffGWo9Qdb1698qmjZtqPUZVbT5+Uq1HqCnP1e7JunY/1rQ2+vXrVzzGlUcAdtpvFu/Q7yF3K/V1zWlu6V6P6/Djaz0BALszfzAHAACAIlceAQAAOtGjjz6aX/3qV8XjDjrooJx++uldMNHrIx4BAAA60aOPPppp06YVjxs2bJh4BAAAeKM6/fTTd+ko3FFe8wgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoKiu1gMAbyzTbmmt9QhV17v3n7NxY/d6XBMu8rtFAGBbfjoAAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi7/MIdKnhG79Z6xGqrr65Ls0tLbUeo8om13oAasx7su4evCcr0JV8xQEAAKDIlUcAoB13Cewu3CUAdB1XHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUFRX6wEAAOga025prfUIVde795+zcWP3elwTLnJ9h12TeGSX5pvc7sE3OQCA7s9PfAAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAiupqPQC8muEbv1nrEaquvrkuzS0ttR6jyibXegAAADqZK48AAAAUufIIAPAG4Y6e3YU7etg1ufIIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAUiUcAAACKxCMAAABF4hEAAIAi8QgAAECReAQAAKBIPAIAAFAkHgEAACgSjwAAABSJRwAAAIrEIwAAAEXiEQAAgCLxCAAAQJF4BAAAoEg8AgAAUCQeAQAAKKqr9QDVNO2W1lqPUFW9e/85Gzd2r8c04SK/rwAAqJbu9vNv4mfgXVm3isfhG79Z6xGqqr65Ls0tLbUeo8om13oAAIBuo7v9/Jv4GXhX1qNSqVRqPQQAAAC7tu5x/RQAAIBOJR4BAAAoEo8AAAAUiUcAAACKxCMAAABF3eqtOnZ1S5YsyQ033JDFixenubk5hx56aC644IJ86EMf2qH//r/+67/yi1/8IgsWLMiaNWvy5z//OQMGDMgHPvCBXHLJJdl77707+RHQkZ1d1/nz5+fOO+/M7373uzzzzDNpbm7OgQcemCOPPDIXX3xx3vGOd3TyI+CVdnZNX6m5uTljxozJ8uXL8/a3vz2zZs2q8sTsiGo8V8eNG7fd/XfccUeGDh1apWnZEda0e6rW1+AXXngh3/3udzN79uysXr069fX1GThwYD7wgQ9k0qRJnTQ9HdnZNR07dmwWLFjwqsd8+ctfzujRo6swLa9GPHaR+fPnZ8KECamvr89pp52WvfbaK7Nnz86nP/3prFmzJhMnTiye44orrsj69etz1FFH5cMf/nB69OiRBQsW5Oabb87s2bNz++23Z9999+2CR8NW1VjXuXPnZtGiRTniiCMyYsSI1NfX5/e//33uvvvu3HvvvfnOd76TY489tgseDUl11vSVvvWtb2XVqlWdMC07qprreswxx+SYY45pt/3AAw+s5sgUWNPuqVrr+uSTT+aCCy7I6tWrc/zxx+dv//Zvs3nz5qxatSpNTU3isQtVY00/8pGPdPgcbWlpybe//e3sscceOe644zpjfF6pQqdrbm6ufPCDH6wcfvjhlWXLlrVtf/755yunnXZaZfDgwZXHH3+8eJ5vf/vblf/5n//ZZltra2vlX/7lXyqNjY2V//f//l+1R+dVVGtdN23a1OH2uXPnVhobGysf/ehHqzUyBdVa05dbunRpZfDgwZXvf//7lcbGxspJJ51U5akpqda6zps3r9LY2FiZMmVKJ07LjrCm3VO11rWlpaVy5plnVo444ojKgw8+2OHHoWt0xvfVl5s1a1alsbGxcskll1RhWnaE1zx2gXnz5mXVqlU5/fTTM3jw4Lbtffr0ySc+8Ym0tLTkrrvuKp7n4x//ePbff/9ttvXo0SOf+MQnkiQLFy6s7uC8qmqt65ve9KYOtx933HHp27evK1ZdqFprutXmzZtz1VVX5d3vfnfOP//8zhiZHVDtdaX2rGn3VK11bWpqysMPP5zx48d3eOdOXZ0b77pKZz9Xp0+fniQZM2bMTs/KjvHs6QJb79EeMWJEu33vfe97tznm9dj6RXDPPfd83efgtevsdX3ooYeyYcOGHHXUUa/7HLw21V7TG2+8MStXrszdd9+dHj16VGdIXrNqr+sf/vCHfP/738+mTZvSv3//HH/88dlnn32qMyw7xJp2T9Va15kzZyZJTj755Dz11FOZM2dOnn/++QwcODAjR45MQ0NDFafm1XTmz0pPP/10Hnjggey333553/ve97pn5LURj13gD3/4Q5LkrW99a7t9ffv2Tb9+/bJy5crXff4f//jHSf7yJKRrVHtd58+fnwULFmTz5s1ZuXJl/vM//zP9+vXL5z73uWqNTEE113TJkiW5+eabc+WVV+btb397NcfkNar2c/Xee+/Nvffe2/bvXr16ZfLkyfnYxz6207OyY6xp91StdV26dGmSZNGiRfnSl76UzZs3t+3bZ5998vWvfz3Dhw+vztC8qs78GfjHP/5xWltb85GPfMTV5C7kM90FXnjhhSTJXnvt1eH+Pn365Omnn35d5/7v//7vfPOb38y+++7rm1wXq/a6LliwIDfeeGPbv9/61rfm+uuvz+GHH75zg7LDqrWmmzdvzuc+97m8613vyvjx46s6I69dtdZ1n332yWc/+9m8733vS//+/fPcc89l/vz5ufbaa/PVr341ffr0yTnnnFPV2emYNe2eqrWuzz77bJLki1/8YsaPH5/zzz8/PXv2zIwZM/LlL385l112WWbOnNnupUBUX2f9DFypVNpud3XLatfymsfd2OrVq3PJJZdky5Ytuf76691is5ubPHlyHnnkkTz00EOZPn163vGOd+Tcc8/NPffcU+vReI2+/vWvZ+XKlbnmmmvcTt6NHHbYYZkwYUIOOeSQ9O7dOwcccEDOOOOM3Hzzzamvr88NN9yQ1tbWWo/Ja2BNu6dKpZIked/73pdPf/rTOfDAA7PPPvtk7NixufDCC/P888/nP/7jP2o8JTtj3rx5eeKJJ3LMMcd0eFWTziMeu0CfPn2SJM8//3yH+1944YXt/kZme9asWZMLLrgg69aty5QpU7yVQw10xromyV/91V/liCOOyI033ph3vOMd+ed//uesW7dup2Zlx1RjTZctW5bvfe97mThxYgYNGlT1GXntOuu5ulVjY2Pe/e53549//ONOvQSBHWdNu6dqrevW84waNardvve///1J/nJrK52rs56r/lBO7YjHLvC2t70tSTr8BrRhw4asX7/+Nf3W5IknnsjYsWOzdu3afP3rX2/7QkjXqva6vlJdXV2GDx+eP//5z3n44Ydf93nYcdVY00ceeSRbtmzJDTfckEGDBm3zvyR5/PHHM2jQoBx99NFVn5+OdfZzNUn69euXJNm0adNOnYcdY027p2qt69bXme+9997t9m3d9tJLL+3EpOyozniubtiwIT//+c+z995756STTqrGmLwG4rELvOc970mS3H///e32PfDAA0nS4RufduSJJ57IuHHjsnbt2nzta1/LBz/4weoNymtSzXXdnrVr1ybxZ8W7SjXW9G1ve1vGjBnT4f+S/33dx5gxYzJ69OjqDs92dfZztaWlJb/73e/So0ePHHTQQa/7POw4a9o9VWtdt96NtWLFinb7tm4bMGDA656THdcZz9Wf/vSn2bx5cz70oQ+lV69eOz8kr02t32jyjaC5ubnygQ98oHL44YdXfve737Vtf/kbpP7+979v2/7ss89WVqxYUXn22We3Oc/q1asr73//+yuDBw+uNDU1ddn8dKxa67pgwYJKa2tru/Pfd999lSFDhlSOOuqoyosvvth5D4Q21VrT7WlsbKycdNJJVZ+bV1etdf3Nb37T7rna3Nxc+bd/+7dKY2NjZcKECZ37QGhjTbunaq3rqlWrKocffnjluOOOqzz99NPbnOfDH/5wpbGxsTJ37tzOf0B0yvfVM844o9LY2FhZtmxZp85Ox3pUKv/3qmI61bx58/Kxj30s9fX1Of3009OnT5/Mnj07TzzxRD75yU/m0ksvbTv2hhtuyI033phJkyZl8uTJbdtHjRqVNWvWZOjQoR2+X06SbY6n81VjXY8++uj069cvf/M3f5MDDzwwL730Uh555JEsXLgw9fX1ufbaa3PyySfX4uG9IVVjTbdn0KBBefvb355Zs2Z15kOgA9X6Gpwkw4YNywEHHJDnn38+CxcuzOOPP57+/fvntttuczWjC1nT7qlaX4NvvfXWfPGLX8yb3/zm/N3f/V169uyZOXPmZM2aNTn77LNz9dVXd/VDe8Oq5vfVpUuX5swzz8yQIUPa/toqXcu9cF3k2GOPzQ9/+MNMmTIlP/vZz9Lc3JxDDz00V1xxRc4444wdOseaNWuSJIsXL87ixYs7PEY8dq1qrOvkyZNz3333ZdGiRVm3bl3bbVJ///d/nwsuuCCHHXZYJz8KXq4aa8qupxrres455+S+++7LggULsn79+tTV1eXggw/OxIkTM378+PTt27eTHwUvZ027p2p9DR47dmwGDBiQadOmZcaMGdmyZUsOPfTQTJw4MWeddVYnPgJeqZrfV7f+lVx/KKd2XHkEAACgyB/MAQAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCwCvMmDEj//AP/5Ajjzwy73nPe3LZZZdl1apVtR4LAGqqR6VSqdR6CADYFbS0tOQf//Efc++99+bggw/O4YcfnpUrV2bZsmXZb7/9MmPGjPTt27fWYwJATbjyCAD/55prrsm9996byZMnp6mpKV/72tdy1113ZcyYMXnmmWdy66231npEAKgZ8QgASRYsWJAf/OAH+ehHP5pJkyZljz3+8i1y3LhxSZIHHnigVuMBQM2JRwBIMmXKlNTX1+eTn/xku3377rtvkuTJJ5/s4qkAYNchHgF4w3v88cezcOHCjBo1KgcccEC7/Rs3bqzBVACwa6mr9QAAUGuzZ89Okjz11FO56qqr2u3/05/+lCTZe++9u3IsANiliEcA3vAWLFiQJFmyZEmWLFmy3ePe+ta3dtVIALDLEY8AvOEtW7YsvXv3zuLFizvc/6lPfSr33ntvhgwZ0rWDAcAuxGseAXhDe+6557J+/frst99+He6vVCptVyaPOeaYrhwNAHYp4hGAN7QNGzYkSRoaGjrcv2TJkqxduzb7779/hg0b1pWjAcAuRTwC8IbWo0ePJElzc3OH+++8884kyZgxY7Z570cAeKPxXRCAN7QDDjggdXV1eeKJJ7J58+Zt9q1YsSJ333133vzmN+fCCy+szYAAsIsQjwC8odXX12f48OHZtGlT/v3f/71t+5NPPplJkyalubk5//qv/5q+ffvWcEoAqL0elUqlUushAKCWFi1alLFjx2bLli0ZPnx4+vTpk7lz5+all17KP/3TP+X888+v9YgAUHPiEQCS3H///ZkyZUoeeeSR9OrVK0cffXQuueSSHHHEEbUeDQB2CeIRAACAIq95BAAAoEg8AgAAUCQeAQAAKBKPAAAAFIlHAAAAisQjAAAAReIRAACAIvEIAABAkXgEAACgSDwCAABQJB4BAAAoEo8AAAAU/X9t2mb+5F5RwQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_, ax = plt.subplots(figsize=(9, 6))\n",
    "\n",
    "bins = np.linspace(0.2, 0.8, 8)\n",
    "ax = az.plot_dist(\n",
    "    ppc_0[\"yl\"].mean(\"yl_dim_0\"),\n",
    "    label=\"model_0\",\n",
    "    kind=\"hist\",\n",
    "    hist_kwargs={\"alpha\": 0.5, \"bins\": bins},\n",
    ")\n",
    "ax = az.plot_dist(\n",
    "    ppc_1[\"yl\"].mean(\"yl_dim_0\"),\n",
    "    label=\"model_1\",\n",
    "    color=\"C1\",\n",
    "    kind=\"hist\",\n",
    "    hist_kwargs={\"alpha\": 0.5, \"bins\": bins},\n",
    "    ax=ax,\n",
    ")\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"$\\\\theta$\")\n",
    "ax.xaxis.set_major_formatter(FormatStrFormatter(\"%0.1f\"))\n",
    "ax.set_yticks([]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example the observed data $y$ is more consistent with `model_1` (because the prior is concentrated around the correct value of $\\theta$) than `model_0` (which assigns equal probability to every possible value of $\\theta$), and this difference is captured by the Bayes factor. We could say Bayes factors are measuring which model, as a whole, is better, including details of the prior that may be irrelevant for parameter inference. In fact in this example we can also see that it is possible to have two different models, with different Bayes factors, but nevertheless get very similar predictions. The reason is that the data is informative enough to reduce the effect of the prior up to the point of inducing a very similar posterior. As predictions are computed from the posterior we also get very similar predictions. In most scenarios when comparing models what we really care is the predictive accuracy of the models, if two models have similar predictive accuracy we consider both models as similar. To estimate the predictive accuracy we can use tools like PSIS-LOO-CV (`az.loo`), WAIC (`az.waic`), or cross-validation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Authored by Osvaldo Martin in September, 2017 ([pymc#2563](https://github.com/pymc-devs/pymc/pull/2563))\n",
    "* Updated by Osvaldo Martin in August, 2018 ([pymc#3124](https://github.com/pymc-devs/pymc/pull/3124))\n",
    "* Updated by Osvaldo Martin in May, 2022 ([pymc-examples#342](https://github.com/pymc-devs/pymc-examples/pull/342))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: Tue May 31 2022\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.9.7\n",
      "IPython version      : 8.3.0\n",
      "\n",
      "sys       : 3.9.7 (default, Sep 16 2021, 13:09:58) \n",
      "[GCC 7.5.0]\n",
      "matplotlib: 3.5.1\n",
      "pymc      : 4.0.0b6\n",
      "numpy     : 1.21.5\n",
      "arviz     : 0.12.0\n",
      "\n",
      "Watermark: 2.3.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv -w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{include} ../page_footer.md\n",
    ":::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "d4ca51fc2fdee62b1a00ff5126f64ae66836e25d3ba6f45d8551026256283997"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}