Bayesian_Inference_pymc3.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Bayesian_Inference_pymc3.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "y_u0Yvlg-cig",
        "outputId": "bd6daabb-9641-45f0-fb78-eae716d25c17",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 462
        }
      },
      "source": [
        "! pip install pymc3 --upgrade"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already up-to-date: pymc3 in /usr/local/lib/python3.6/dist-packages (3.9.3)\n",
            "Requirement already satisfied, skipping upgrade: typing-extensions>=3.7.4 in /usr/local/lib/python3.6/dist-packages (from pymc3) (3.7.4.3)\n",
            "Requirement already satisfied, skipping upgrade: scipy>=0.18.1 in /usr/local/lib/python3.6/dist-packages (from pymc3) (1.4.1)\n",
            "Requirement already satisfied, skipping upgrade: dataclasses; python_version < \"3.7\" in /usr/local/lib/python3.6/dist-packages (from pymc3) (0.7)\n",
            "Requirement already satisfied, skipping upgrade: fastprogress>=0.2.0 in /usr/local/lib/python3.6/dist-packages (from pymc3) (1.0.0)\n",
            "Requirement already satisfied, skipping upgrade: h5py>=2.7.0 in /usr/local/lib/python3.6/dist-packages (from pymc3) (2.10.0)\n",
            "Requirement already satisfied, skipping upgrade: pandas>=0.18.0 in /usr/local/lib/python3.6/dist-packages (from pymc3) (1.1.2)\n",
            "Requirement already satisfied, skipping upgrade: contextvars; python_version < \"3.7\" in /usr/local/lib/python3.6/dist-packages (from pymc3) (2.4)\n",
            "Requirement already satisfied, skipping upgrade: arviz>=0.9.0 in /usr/local/lib/python3.6/dist-packages (from pymc3) (0.10.0)\n",
            "Requirement already satisfied, skipping upgrade: theano>=1.0.5 in /usr/local/lib/python3.6/dist-packages (from pymc3) (1.0.5)\n",
            "Requirement already satisfied, skipping upgrade: numpy>=1.13.0 in /usr/local/lib/python3.6/dist-packages (from pymc3) (1.18.5)\n",
            "Requirement already satisfied, skipping upgrade: patsy>=0.5.1 in /usr/local/lib/python3.6/dist-packages (from pymc3) (0.5.1)\n",
            "Requirement already satisfied, skipping upgrade: six in /usr/local/lib/python3.6/dist-packages (from h5py>=2.7.0->pymc3) (1.15.0)\n",
            "Requirement already satisfied, skipping upgrade: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.18.0->pymc3) (2018.9)\n",
            "Requirement already satisfied, skipping upgrade: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.18.0->pymc3) (2.8.1)\n",
            "Requirement already satisfied, skipping upgrade: immutables>=0.9 in /usr/local/lib/python3.6/dist-packages (from contextvars; python_version < \"3.7\"->pymc3) (0.14)\n",
            "Requirement already satisfied, skipping upgrade: netcdf4 in /usr/local/lib/python3.6/dist-packages (from arviz>=0.9.0->pymc3) (1.5.4)\n",
            "Requirement already satisfied, skipping upgrade: packaging in /usr/local/lib/python3.6/dist-packages (from arviz>=0.9.0->pymc3) (20.4)\n",
            "Requirement already satisfied, skipping upgrade: matplotlib>=3.0 in /usr/local/lib/python3.6/dist-packages (from arviz>=0.9.0->pymc3) (3.2.2)\n",
            "Requirement already satisfied, skipping upgrade: setuptools>=38.4 in /usr/local/lib/python3.6/dist-packages (from arviz>=0.9.0->pymc3) (50.3.0)\n",
            "Requirement already satisfied, skipping upgrade: xarray>=0.16.1 in /usr/local/lib/python3.6/dist-packages (from arviz>=0.9.0->pymc3) (0.16.1)\n",
            "Requirement already satisfied, skipping upgrade: cftime in /usr/local/lib/python3.6/dist-packages (from netcdf4->arviz>=0.9.0->pymc3) (1.2.1)\n",
            "Requirement already satisfied, skipping upgrade: pyparsing>=2.0.2 in /usr/local/lib/python3.6/dist-packages (from packaging->arviz>=0.9.0->pymc3) (2.4.7)\n",
            "Requirement already satisfied, skipping upgrade: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=3.0->arviz>=0.9.0->pymc3) (0.10.0)\n",
            "Requirement already satisfied, skipping upgrade: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=3.0->arviz>=0.9.0->pymc3) (1.2.0)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ryM293nQ91Or",
        "outputId": "404f2401-d8be-405e-e2af-89fef4f0ce44",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 289
        }
      },
      "source": [
        "! pip install arviz"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: arviz in /usr/local/lib/python3.6/dist-packages (0.10.0)\n",
            "Requirement already satisfied: netcdf4 in /usr/local/lib/python3.6/dist-packages (from arviz) (1.5.4)\n",
            "Requirement already satisfied: setuptools>=38.4 in /usr/local/lib/python3.6/dist-packages (from arviz) (50.3.0)\n",
            "Requirement already satisfied: numpy>=1.12 in /usr/local/lib/python3.6/dist-packages (from arviz) (1.18.5)\n",
            "Requirement already satisfied: pandas>=0.23 in /usr/local/lib/python3.6/dist-packages (from arviz) (1.1.2)\n",
            "Requirement already satisfied: matplotlib>=3.0 in /usr/local/lib/python3.6/dist-packages (from arviz) (3.2.2)\n",
            "Requirement already satisfied: scipy>=0.19 in /usr/local/lib/python3.6/dist-packages (from arviz) (1.4.1)\n",
            "Requirement already satisfied: packaging in /usr/local/lib/python3.6/dist-packages (from arviz) (20.4)\n",
            "Requirement already satisfied: xarray>=0.16.1 in /usr/local/lib/python3.6/dist-packages (from arviz) (0.16.1)\n",
            "Requirement already satisfied: cftime in /usr/local/lib/python3.6/dist-packages (from netcdf4->arviz) (1.2.1)\n",
            "Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.23->arviz) (2018.9)\n",
            "Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.23->arviz) (2.8.1)\n",
            "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=3.0->arviz) (1.2.0)\n",
            "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=3.0->arviz) (0.10.0)\n",
            "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=3.0->arviz) (2.4.7)\n",
            "Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from packaging->arviz) (1.15.0)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fe4PPz3VGN5k"
      },
      "source": [
        "\n",
        "Bayesian prior beliefs act as regularization from making inference from the data\n",
        "\n",
        "### Excercise 1 - Fitting a Normal distribution\n",
        "\n",
        "1. Generate data drawn from Normal distribution\n",
        "2. Fit a model to this data and estimate the parameters (mean & std) of this distribution\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z1z5824h8qRS"
      },
      "source": [
        "import pymc3 as pm\n",
        "import numpy as np\n",
        "import scipy\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uS5JTNtjDeOM"
      },
      "source": [
        "# Generate data from normal distribution with mean=0.5 & std=0.35\n",
        "data = scipy.stats.norm(loc=0.5, scale=0.35).rvs(1000)"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5CtP5uQvHJ06",
        "outputId": "0bf723cf-0584-4863-a76e-0310f2b2425f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "data.shape"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(1000,)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p_ZM54S_bG5c"
      },
      "source": [
        "#### Estimating parameters of the data generated above.\n",
        "\n",
        "Since standard deviation is a positive value,we use Half-normal distribution for std since its range is [0,inf]."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xGQsUJTZYVAQ",
        "outputId": "0cd8c443-1925-4857-96f1-eb9a4ec8ea1b",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 142
        }
      },
      "source": [
        "# Instantiate model with context. This model computes Log-likelihood and its gradients\n",
        "n_post_samples = 1000\n",
        "\n",
        "with pm.Model() as model:\n",
        "  # Define priors for model parameters\n",
        "  mu = pm.Normal('mu', mu=0.0, sd=0.5)\n",
        "  std = pm.HalfNormal('std', sd=0.5)\n",
        "  \n",
        "  # Likelihood function. Use observed data to compute likelihood\n",
        "  y = pm.Normal('y', mu=mu, sd=std, observed=data)\n",
        "  \n",
        "  # Get posterior samples from NUTS sampler\n",
        "  trace = pm.sample(n_post_samples)\n"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Sequential sampling (2 chains in 1 job)\n",
            "NUTS: [std, mu]\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "\n",
              "    <div>\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",
              "      <progress value='2000' class='' max='2000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
              "      100.00% [2000/2000 00:02<00:00 Sampling chain 0, 0 divergences]\n",
              "    </div>\n",
              "    "
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "\n",
              "    <div>\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",
              "      <progress value='2000' class='' max='2000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
              "      100.00% [2000/2000 00:01<00:00 Sampling chain 1, 0 divergences]\n",
              "    </div>\n",
              "    "
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 4 seconds.\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dJlSPxltHvF5",
        "outputId": "b9eb05da-69dd-4bc1-8a9b-f92332ce06de",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 266
        }
      },
      "source": [
        "plt.hist(data, histtype='stepfilled', bins=120)\n",
        "plt.show()"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD5CAYAAAA+0W6bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAOcUlEQVR4nO3db4il5XnH8e+vxprShKrZYbOom01TSWtfdA3DNsFSrCZFDaihEvSF2RTD2BJpQvNGUqgh9MW2NBEKrckmituSasT80TamrTUWKzTSWWv9S9CI0l1Wd62JGhoSVq++mGey43hmztnzb+ae+X7gMOc8z3POueaesz/ufc517pOqQpLUnp9b6wIkScMxwCWpUQa4JDXKAJekRhngktQoA1ySGvWmfgckeTNwH3BSd/ztVXVdkncCtwJvA/YDV1bVT1d7rC1bttSOHTtGLlqSNpP9+/e/UFUzy7f3DXDgJ8B5VfWjJCcC9yf5NvDHwPVVdWuSLwBXATes9kA7duxgfn5+iPIlafNK8myv7X1PodSCH3U3T+wuBZwH3N5t3wdcOoY6JUkDGugceJITkjwEHAbuBr4P/LCqjnaHHABOm0yJkqReBgrwqnq1qnYCpwO7gF8d9AmSzCWZTzJ/5MiRIcuUJC13XF0oVfVD4F7gfcDJSRbPoZ8OHFzhPnuraraqZmdm3nAOXpI0pL4BnmQmycnd9V8APgA8wUKQX9Ydthu4Y1JFSpLeaJAulG3AviQnsBD4t1XVPyZ5HLg1yZ8B/wXcOME6JUnL9A3wqnoYOLvH9qdZOB8uSVoDfhJTkhplgEtSowY5By5tejuu/RYAz+z54ND3Hfb+0kqcgUtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RG2UaoDcN2PW02zsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIa5WqEEq5kqDY5A5ekRhngktQoA1ySGtU3wJOckeTeJI8neSzJJ7rtn0lyMMlD3eWiyZcrSVo0yJuYR4FPVdWDSd4K7E9yd7fv+qr6y8mVJ0laSd8Ar6pDwKHu+itJngBOm3RhkqTVHVcbYZIdwNnAA8A5wDVJPgLMszBL/0GP+8wBcwDbt28fsVxtRLbwvdHimDgeWs3Ab2ImeQvwNeCTVfUycAPwLmAnCzP0z/W6X1XtrarZqpqdmZkZQ8mSJBgwwJOcyEJ4f6Wqvg5QVc9X1atV9RrwJWDX5MqUJC03SBdKgBuBJ6rq80u2b1ty2IeAR8dfniRpJYOcAz8HuBJ4JMlD3bZPA1ck2QkU8Axw9UQqlCT1NEgXyv1Aeuy6a/zlSJIG5ScxJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNOq5vpZdaM61vvPdb5LUWnIFLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRtlGqHVvaStgL6O07vV67Em2BI7jsafVGqn1zxm4JDXKAJekRhngktSovgGe5Iwk9yZ5PMljST7RbT81yd1Jnux+njL5ciVJiwaZgR8FPlVVZwHvBT6e5CzgWuCeqjoTuKe7LUmakr4BXlWHqurB7vorwBPAacAlwL7usH3ApZMqUpL0RsfVRphkB3A28ACwtaoOdbueA7aucJ85YA5g+/btw9YpjU2/tsRB72sLn9bawG9iJnkL8DXgk1X18tJ9VVVA9bpfVe2tqtmqmp2ZmRmpWEnSMQMFeJITWQjvr1TV17vNzyfZ1u3fBhyeTImSpF4G6UIJcCPwRFV9fsmuO4Hd3fXdwB3jL0+StJJBzoGfA1wJPJLkoW7bp4E9wG1JrgKeBT48mRIlSb30DfCquh/ICrvPH285kqRB+UlMSWqUqxGqeaO0BY7yuON+3kn9Htq4nIFLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ytUItSG5sp82A2fgktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1yj5wjWSx3/qZPR9cl4/X67HXg/VUi9rlDFySGmWAS1Kj+gZ4kpuSHE7y6JJtn0lyMMlD3eWiyZYpSVpukBn4zcAFPbZfX1U7u8td4y1LktRP3wCvqvuAF6dQiyTpOIxyDvyaJA93p1hOGVtFkqSBDBvgNwDvAnYCh4DPrXRgkrkk80nmjxw5MuTTaRx2XPutn13W6rkljc9QAV5Vz1fVq1X1GvAlYNcqx+6tqtmqmp2ZmRm2TknSMkMFeJJtS25+CHh0pWMlSZPR95OYSW4BzgW2JDkAXAecm2QnUMAzwNUTrFGS1EPfAK+qK3psvnECtUiSjoOfxJSkRhngktQoVyPUwCa5UuBGtJ7aJpfW4t9v43AGLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplG6HWzHpqs1uvHCOtxhm4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapRthHqDYVau63WfzdgCt5a/s6tFbj7OwCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjbCPU2K3WSjdom91mbEEcF8du83AGLkmNMsAlqVEGuCQ1qm+AJ7kpyeEkjy7ZdmqSu5M82f08ZbJlSpKWG2QGfjNwwbJt1wL3VNWZwD3dbUnSFPUN8Kq6D3hx2eZLgH3d9X3ApWOuS5LUx7BthFur6lB3/Tlg60oHJpkD5gC2b98+5NNprdiStr7599ncRn4Ts6oKqFX2762q2aqanZmZGfXpJEmdYQP8+STbALqfh8dXkiRpEMMG+J3A7u76buCO8ZQjSRrUIG2EtwD/Abw7yYEkVwF7gA8keRJ4f3dbkjRFfd/ErKorVth1/phrkSQdBz+JKUmNMsAlqVEuJ9swv4VckzTM66tXX7qvz8lxBi5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZRvhJre07ct2r83Fv337nIFLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRtlGuEFNskXMb0Jvm3+/jcMZuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUbYSbVK9WMtvL1I9fpL2+OAOXpEYZ4JLUKANckho10jnwJM8ArwCvAkeranYcRUmS+hvHm5i/U1UvjOFxJEnHwVMoktSoUWfgBfxLkgK+WFV7lx+QZA6YA9i+ffuIT7cx9GrFGmX1wH73tT1Q4+YXIq8Po87Af6uq3gNcCHw8yW8vP6Cq9lbVbFXNzszMjPh0kqRFIwV4VR3sfh4GvgHsGkdRkqT+hg7wJL+Y5K2L14HfBR4dV2GSpNWNcg58K/CNJIuP8/dV9U9jqUqS1NfQAV5VTwO/McZaJEnHwTZCSWqUqxGuQ73a/hZbtWwJ1LT5mlu/nIFLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRtlGOKJBV2Ub9UuEBz3Wli8NYxqvm2mvYLgZVkx0Bi5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqPsA2e4flH7raUF4/q3MOgyyr3+jS7u36j93itxBi5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIalaqa2pPNzs7W/Pz8UPcdpdWv3/GrtS/1O66Xpfe13VAar16thYPeZxSD5sRK9xmlhiT7q2p2+XZn4JLUKANckho1UoAnuSDJ95I8leTacRUlSepv6ABPcgLw18CFwFnAFUnOGldhkqTVjTID3wU8VVVPV9VPgVuBS8ZTliSpn1EC/DTgf5bcPtBtkyRNwdBthEkuAy6oqo91t68EfrOqrll23Bww1918N/C94csdyRbghTV67vXI8TjGsXg9x+OY9TIW76iqmeUbR1lO9iBwxpLbp3fbXqeq9gJ7R3iesUgy36uPcrNyPI5xLF7P8ThmvY/FKKdQ/hM4M8k7k/w8cDlw53jKkiT1M/QMvKqOJrkG+GfgBOCmqnpsbJVJklY10jfyVNVdwF1jqmXS1vw0zjrjeBzjWLye43HMuh6Lqa6FIkkaHz9KL0mN2rABnuTUJHcnebL7ecoKx72a5KHusqHehO231EGSk5J8tdv/QJId069yegYYj48mObLk9fCxtahzGpLclORwkkdX2J8kf9WN1cNJ3jPtGqdlgLE4N8lLS14XfzrtGleyYQMcuBa4p6rOBO7pbvfy46ra2V0unl55kzXgUgdXAT+oql8Brgf+fLpVTs9xLP3w1SWvhy9Ptcjpuhm4YJX9FwJndpc54IYp1LRWbmb1sQD49yWvi89OoaaBbOQAvwTY113fB1y6hrWshUGWOlg6RrcD5yfJFGucJpd+WKKq7gNeXOWQS4C/rQXfBU5Osm061U3XAGOxbm3kAN9aVYe6688BW1c47s1J5pN8N8lGCvlBljr42TFVdRR4CXjbVKqbvkGXfvi97pTB7UnO6LF/s3CpjNd7X5L/TvLtJL++1sUsGqmNcK0l+Vfg7T12/cnSG1VVSVZqt3lHVR1M8svAd5I8UlXfH3etasI/ALdU1U+SXM3C/07OW+OatPYeZCEnfpTkIuCbLJxaWnNNB3hVvX+lfUmeT7Ktqg51//U7vMJjHOx+Pp3k34CzgY0Q4IMsdbB4zIEkbwJ+Cfjf6ZQ3dX3Ho6qW/u5fBv5iCnWtVwMtlbEZVNXLS67fleRvkmypqjVfI2Ujn0K5E9jdXd8N3LH8gCSnJDmpu74FOAd4fGoVTtYgSx0sHaPLgO/Uxv1gQN/xWHaO92LgiSnWt97cCXyk60Z5L/DSklOSm0qSty++N5RkFwu5uS4mOk3PwPvYA9yW5CrgWeDDAElmgT/oVlH8NeCLSV5j4Y+yp6o2RICvtNRBks8C81V1J3Aj8HdJnmLhTZzL167iyRpwPP4oycXAURbG46NrVvCEJbkFOBfYkuQAcB1wIkBVfYGFT1hfBDwF/B/w+2tT6eQNMBaXAX+Y5CjwY+Dy9TLR8ZOYktSojXwKRZI2NANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RG/T8beRfxqelWXwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MkvGVHCDH35O",
        "outputId": "6a7e7148-8843-4127-e97a-44a02506dcce",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "trace['mu'].shape"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(2000,)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BC7__LM79sev",
        "outputId": "ef7dd4ce-f20e-4725-b373-f07b26e75618",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 410
        }
      },
      "source": [
        "pm.plot_posterior(trace)\n",
        "plt.show()"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/arviz/data/io_pymc3.py:91: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n",
            "  FutureWarning,\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+kAAAFTCAYAAABF801FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3hUVf4/8PedmfTeeyVAQk0ghNCRjqhI8UtQKWLFXtayooisK7o/Vt0VcVdRUVBEKUEERJBeQgkJLQTSe+89mZn7+wPISg8hmTPl/XqePJg7d2beEyFnPnPO/RxJlmUQERERERERkXgK0QGIiIiIiIiI6CIW6URERERERER6gkU6ERERERERkZ5gkU5ERERERESkJ1ikExEREREREekJFulEREREREREeoJFOhERERGRnpIkKVySpEWSJDm38fyRkiTJkiSN7ORoRNRJWKQTEREREemvcADvAGhTkU5Eho9FOhEREREREZGeYJFOZIAuLXuTJUkKlSRpuyRJdZIkZUuS9Mil22dJkpQsSVKtJEm7JUnq8qf7ypIkLbrq8QIvHZ+r21dCREREkiR1kyRpoyRJxZIkNV4a03+WJOkxAN9cOi3l0lgtS5IUeOl+bpIk/SBJUrUkSZWSJH0HwFHQyyCiDqISHYCI7sjPAL4EsBTA0wC+liSpK4CRAN4AYAbgXwB+ADBQUEYiIiK6uS0AKgDMB1AKwAfA3QA2A3gPwFsAHgCQe+n8gkt/bgDQF8CbAFIAzADwqc5SE1GnYJFOZNj+nyzL3wGAJEnHAdwL4EkAQbIsV1867gXgX5IkBciynCUuKhEREV1NkiRXACEAJsuy/Mufbvrh0u1pl75PlGU59U/3GwtgKICZsiz/eOnwdkmStgHw7fzkRNRZuNydyLBtu/wfsixXACgGEHe5QL8k+dKffroMRkRERG1SBiAdwAeSJD1+aUVcWwwCoAGw/qrjP17nXCIyICzSiQxbxVXfN9/gGABYdn4cIiIiuh2yLMsAxgI4DmAJgAuSJKVLkjT/Fnf1AlAhy3LLVceLOiEmEekQl7sTmZ4mAOZXHXMREYSIiIgAWZbTAcyWJEnCxWvMnwWwXJKkzJvcrQCAkyRJZlcV6h6dl5SIdIEz6USmJwtAr6uOTRIRhIiIiP5HvigRwMuXDvXCxQ/XAcDqqtMPA1ACmHbV8ZjOS0hEusCZdCLT8yOAtyRJWgAgDsAwADPFRiIiIjJNkiT1wcWdWNYCSMXFwnsuADWAXZf+BIBnJEn6FkALgFOyLO+QJOkAgP9eaj53ubv71R/EE5GB4Uw6kelZAmAZLi6liwUQBmCW0ERERESmqxBANi7Onv8CYA0AbwD3yLIcL8vySQCLcHEHlwMAjl26HQCmAtiKi2P7WlycgHtWl+GJqONJF3tVEBEREREREZFonEknIiIiIiIi0hMs0omIiIiIiIj0BIt0IiIiIiIiIj3BIp2IiIiIiIhIT7BIJyIiIiIiItITt9onna3fiYiI9IN0h/fnmE5ERKQfbjqmcyadiIiIiIiISE+wSCciIiIiIiLSEyzSiYiIiIiIiPQEi3QiIiIiIiIiPcEinYiIiIiIiEhPsEgnIiIiIiIi0hMs0omIiIiIiIj0BIt0IiIiIiIiIj3BIp2IiIiIiIhIT7BIJyIiIiIiItITKtEBiEg/yLKM+mYN6prVcLY2h0rJz/CIiIhuRpZlVDeqIcsyZPniMVtLFcw4hhLRHWCRTmTiDqeV4ef4HOxMKkJ1oxoAoJCAYDdbjA51xwORvghxtxOckoiISDyNVkZCdgV2ny/G/pRSpBXXoq5Zc8U5CgnwcrBCuJ8jBgQ6ISrIBaGedlAoJEGpicjQSPLlj/2u76Y3EpHhKqlpwl83nMbOc0Wwt1RhbA9PdPOwhbW5EiU1TUjIqURcehnUWhmT+3rjzbvD4G5vKTo2kSm703f4HNOJ2qm2SY3v47Lw3eEs5FU2QKmQ0N/fCT197OHtYAWV8n//PCvqmpFRVo/4zHLkVzUCAOwtVRjW1Q339vXCqFAPmKs4005k4m46prNIJzJBh9JK8dwPCahpUuPlsd0wd3AgLM2U15xXVtuEL/anY+XBTFiZK7FkSm9M7O0lIDERgUU6kc5ptTLWHs/BP38/j9LaZkQHO2NmlD9GdneHg5XZLe+fW1GPY5nlOJJejp3nilBa2wxXWws8GOWHh6ID4MEPv4lMFYt0IvqfP84VYf73JxDgbI3PHuqHbh63XsqeVlKLl9cm4mRuFV4d3x1Pj+wCSeKyPSIdY5FOpEOFVY14dd1J7E8pxYBAJ7x5dxgi/J3a/XhqjRb7U0ux6nAWdp8vhplCgVmDAvDMXSFwtjHvwOREZABYpBPRRZtP5uOltYno4W2Pbx+JgtNtvCloVmvx6rqT2JSYj7mDA/H2PT2g5PV1RLrEIp1IRxJzKvHYt8dQ16TBW/eE4cEo/w79cDqrrA6f7U7FuvhcWJurMH9kFzw6NOi6q9qIyCixSCci4FBqKR7+6ggiA5zx1dxI2Fneepne1bRaGUu2ncOX+zMwJcIHSx/oy0KdSHdYpBPpwO7zxZi/Oh5udhb4es4AdG3DirP2Si2uwYe/nceOpCL4O1vj7Xt6YGwPj057PiLSGyzSiUxdfmUD7vn0AJxtzBH7zBDYWtzZxg7LdqVg6e8XEDPAD0um9ubSdyLdYJFO1MmOZpRj1ldH0NXDFisfiYKrrYVOnnd/Sgne3ZyE1OJaTInwwbuTe8K+HR+mE5HBYJFOZMqa1Br833/jkFZci03PDkEXN9sOedyl289j2e5UzB0ciHfu7cFCnajzsUgn6kRJ+dWY8cVhuNlZYN1Tg3V+nXiLRovPdqfi012p8LS3xL9iwhEZ6KzTDESkMzcd07n/A5GR+/uWcziZU4mlD/TtsAIdAF4Z1w2PDg3CykOZ+Gx3aoc9LhERka7lVzZgzjdHYWuhwqpHBwpp5GamVODFMd3w81ODoFJKiPkiDt8fydJ5DiISj0U6kRHbfb4Y3x3OwqNDgzChl2eHPrYkSXhrUtjFa9N/v4D18bkd+vhERES60NiiwZOr4tHYrMF386Lg42glNE8/fyf88uxQDO3qigUbz+Ct2NNo0WiFZiIi3WKRTmSkyuua8dq6U+juYYdXx3fvlOeQJAkfTuuDwV1c8Pr6UziQUtopz0NERNRZ3oo9g9N5Vfh4RninNom7HQ5WZvhqzgA8OSIYq+Oy8ei3x1HXpBYdi4h0hEU6kZH6YNs5VNQ14+MZ4Z26pYu5SoH/zOqPLm62eGp1PM4VVHfacxEREXWkTYl5WBefi+dHd8UYPeuqrlRI+OvEMHw4rTcOppZixheHUVzTKDoWEekAi3QiIxSfVYGfjufi0WFB6OFt3+nPZ29phm8eGQBbCxUe+/Y4SmubOv05iYiI7kROeT3e2ngGkQFOeH5UiOg4NzRjgD++nN0facV1mP75YRRUNYiORESdjEU6kZHRaGW8HXsGnvaWeH5UV509r7ejFb6Y3R+ltU14evUJNKt5/RwREeknrVbGq+tOQgbw8YxwqJT6/ZZ4VKgHvn98ICrqmjHzizgUVXNGnciY6fdvJCK6bd8fyUJSQTXevqcHbO5wP/Tb1cfXEf+Y3gdHM8vx7uazOn1uIiKitvrhaDbi0svx1qQw+Dlbi47TJv38nbByXhRKapow88s4Ln0nMmIs0omMSGV9M5ZuP48hIS64u3fHdnNvq8nhPnhyRDC+P5KN1XHcOoaIiPRLXmUDlmw9hyEhLpgxwE90nNvSP8AJ3zwShYLKRjz05RGU8fIyIqPEIp3IiHy2OxU1TWq8fU8PSJIkLMdr40MxsrsbFv1yFkfSy4TlICIi+jNZlvHmhtPQysAHU/sIHSvbKyrIGV/PHYDs8no89t1xNLZoREciog7GIp3ISORW1OPbQ1mY1s8XoZ6d3yzuZpQKCf+KiYC/izWe/v4EcivqheYhIiICgF9PFWDvhRK8NqG7wSxzv55BXVzwr5hwJOZU4qW1idBqZdGRiKgDsUgnMhIf/X4BkICXx3YTHQXAxT1ev5wdiWa1Fk98F4/6Zu7vSkRE4tQ3q/H+1nPo4WWP2YMCRce5YxN6eWHB3WHYdqYQH/yWLDoOEXUgFulERiApvxobE/PwyOBAeDtaiY7TqoubLf49MwLnCqvx6rpTkGV+0k9ERGJ8vicNBVWNeHdyTygVhrfM/XoeHRqEOYMC8MW+dKw6nCk6DhF1EBbpREbgw9+SYW9phqdH6t8+r3eFuuP1CaHYcqoAy/ekiY5DREQmKLusHv/dl477w70xINBZdJwOI0kSFt7bE6ND3fHu5iTEZ5WLjkREHYBFOpGBO5RWir0XSvDMXV3gYG0mOs51PTk8GJPDvbH09/PYmVQkOg4REZmYv21Jgkoh4Y2JYaKjdDilQsJHM8Lh7WiFZ75PYMd3IiPAIp3IgMmyjI93XICnvaVeX18nSRI+nNYHvbwd8OLaRKQU1YiOREREJmLvhRLsSCrCc6O6wtPBUnScTuFgZYblD/VDeX0zXlybCA0byREZNBbpRAbscFoZjmVWYP7ILrA0U4qOc1OWZkp8Mbs/LM2UePy746iqbxEdiYiIjJxGK+P9LecQ4GKNeUMDRcfpVL18HPDufT2xP6UUy3alio5DRHeARTqRgZJlGZ/sTIGHvQVmDPATHadNvBys8N9Z/ZBX2YCXf+KWMURE1LnWn8jF+aIavD4hFBYq/f4wuyPEDPDD1AgffPLHBcSll4mOQ0TtxCKdyEAdTivD0cxyPD0yRO9n0f+sf4Az3r6nB/5ILsay3fykn4iIOkdjiwYf77iAvn6OmNjLU3QcnZAkCe9N6QV/Z2u8tu4Utz8lMlAs0okMkCzL+OQPw5pF/7NZ0QGYEuGDj3dewJ7zxaLjEBGREfrmYCYKqhrx14mhkCTj2HKtLazNVfhwWh9kl9fj/20/LzoOEbUDi3QiA3Q4vQxHM8oxf4T+X4t+PZIk4f0pvdHdww4v/JiInPJ60ZGIiMiIVNY3Y/meVIwKdUd0sIvoODoXHeyC2YMCsPJQJo5lcls2IkPDIp3IAP1rZwrc7SwQE+UvOkq7WZkr8d9Z/aGVZTy5Kh6NLRrRkYiIyEh8tjsVtU1qvD4hVHQUYV6fEAofRyu8tu4UGpo5xhIZEhbpRAbmcFoZjmSUG0RH91sJcLHBJzPCkVRQjb/9miQ6DhERGYHcinp8eygL0/r5orunneg4wthYqPCPaX2QUVqHj3Zw2TuRIWGRTmRgPt11cRZ9pgHPov/Z6DAPPD4sCN8fycYf54pExyEiIgP30e8XIEnAy2O7iY4i3OAQV8QM8MPXBzNxoahGdBwiaiMW6UQGJDGnEofSyvD4sGCDm0VvamrCK6+8And3d9jY2GDSpEnIzMwEAPxlfHeEednjtXWnUFLTdMX9Ro4cCUmSrvlqbGwU8CqIiEifJeVXY2NiHuYOCYS3o5XoOLflZuPkzTz11FMIDQ2Fra0tnJycMHz4cOzcubP19tcmhAL5Z9Hd0/6asfSNN97oxFdERO3FIp3IgPxnTxrsLVWYOdDwZtGff/55rFy5EkuXLsW6detQWlqKsWPHorGxERYqJf4VE47aJjVeW3cSsnzl/ul33XUXDh8+fMWXhYWFoFdCRET66oPfkmFvaYanR4SIjnLbbjZO3kxDQwOeffZZbNy4EatXr4arqysmTpyIuLg4AICzjXnrTjCvLFl2xVj6zDPPdPrrIqLbpxIdgIjaJq2kFtuTCvHMyBDYWhjWP93c3Fx89dVX+PrrrzF79mwAQJ8+fRAUFITVq1fjscceQzcPO/x1YigWbU7C6rgszBoU2Hp/Z2dnREdHC0pPRESG4GBqKfZdKMGbd4fCwdpMdJzb0pZx8ka+/fbbK76fMGECgoKC8MMPP7SOnaPDPLAEwI4CCyzuFwlrc8N6H0FkajiTTqRH5s6di8jISGzZsgU9evSAtbU1Jk2ahPLycnz44x4UrXkT70zth8jISJw6dar1flqtFh988AFCQkJgYWGBbt26XTNob9myBWPHjoW7uzvs7e0RHR2N33///YpzFi1aBFdXVyQkJCA6OhrW1taIiIjA/v377+h1XX6eqVOnth7z8fHB0KFDsW3bttZjcwYHYkQ3N7y/NZnbshER0TVuNE6Wlpbh7e92ouLnBXhpYh+jHSfbQqlUwtHREc3Nzf87pri4T3xJbRNW7M+4o6xE1PlYpBPpmezsbCxcuBDvvfcevvjiCxw6dAizH3kU373/Eu66ezLWrVsHtVqNmJiY1mXhzz33HN577z088cQT2LJlC6ZMmYJ58+bh119/bX3cjIwM3HvvvVi1ahXWr1+PwYMHY+LEiTh48OAVz19fX485c+bgySefxPr162FhYYGpU6eivv5/RbNWq4Varb7pl0bzv+1ekpOT4evrC1tb2yueKywsDMnJya3fS5KEJVN7QyEBb8WeaX19v//+O6ytrWFtbY3x48df8caLiIhMy/XGyckzZ+Pgl28jZkaMUY+TNyLLMtRqNcrKyvDxxx8jJSUF8+bNu+a8ip/fxgtjQ+EfEID33nvvigxEpEdkWb7ZFxHp0Jw5c2SlUimnpqa2Hnv11VdlALLbPS/L2WV1sizL8pYtW2QAclJSkpySkiJLkiSvXLnyiseaNWuWHBkZed3n0Wg0cktLizxu3Dj5kUceaT3+zjvvyADkP/74o/VYQkKCDEDetm3bFTkB3PRrxIgRrec/9thjct++fa/JsWDBAtnLy+ua418fSJcDXv9V3pSYJy9cuFD++uuv5X379smrVq2SQ0NDZXt7ezkjI+PmP0wi43OrMZtjOhm9642Tr7zyFxmA3GvmX2W1RivLsvGPk1dbs2ZN6+Pa2NjImzZtuuL2EydOyG+88Yb85fc/y54x78nR9z0sKxQK+fnnn7/lYxNRp7jpmM0LUoj0TGBgILp06dL6vY9/EADgnglj4edsDQAICbnYECcvLw9paWlQKBSYMmUK1Gp16/1Gjx6NNWvWQKPRQKlUIjc3FwsWLMDOnTtRUFDQOrswZMiQK57f3NwcI0eObP2+R48eAC5eL3fZokWL8Oyzz970ddjZtX9v2tmDAhGbkIfFm89i5+sL4GhtDgAYNmwYxowZg9DQUHzyySf45JNP2v0cRERkmK4eJ6vNXQAArz86vXVZt7GPk1cbP348jh07htLSUnz//feIiYnB1q1bW3NGREQgIiICAHDB/CRiE/rhyUA3LF/+Kd5++224urp2WBYiunMs0on0jKOj4xXfH8uuBgDMH9en9Zi5+cWitbGxEaWlpdBoNHBwcLju4xUUFMDb2xv33XcfampqsHjxYoSEhMDGxgYLFy5EcXHxFefb2dlBofjflTB/fq7L/P394evre9PXIUlS6387OTmhqqrqmnMqKirg5OR0zXGlQsKSqX1w77IDWLI1GR9O/99r9/T0xJAhQ3DixImbPj8RERmnP4+TTWoNdqeUAwAmR3VtPW7s4+TVnJycEBkZCeBi47j8/HwsXLgQ+/btu+bcF8d0Q2xiPqo8+kOtVuPUqVMYNWrULZ+DiHSHRTqRHmts0eBAaikAIMzL/rrnODs7Q6VS4eDBg1e8abjM3d0dqampSEhIwLZt2zBhwoTW2xoaGtqVa968edc03LnaiBEjsGfPHgBAaGgocnJyUFdXBxsbm9ZzkpOTERoaet379/C2x+PDgvGfvWmY2s8HA4NdWm+7vL8rERGZtrXHclBR3wIANxwXjHWcvJmIiAj8+OOP173N29EKs6MD8Pm68wBu/HMjInFYpBPpsZ+P56C2SX3Tc0aNGgWNRoOqqiqMHTv2uudcfpPx573Fs7KycPDgQfTp0+e697mZ213GN27cOADAxo0b8fDDDwMA8vPzsX//fixfvvyGj/HC6K74JTEPi39Nwi/PDoVSIaGwsBAHDhy4bkMcIiIyHQ3NGny6KxXBrjYou8l5xjxOXo8syzh8+DCCgoJueM4TI4Lx0fuHoFCq2vX6iKhzsUgn0lNqjRb/3ZcOf2frm7756N69O5566inExMTgtddeQ2RkJBobG3H27FlcuHABK1asQGhoKHx9ffHKK6/gb3/7G2pqavDOO+/Ax8enXdkCAwMRGBjY5vN9fX3x6KOP4sUXX4Qsy3Bzc8OiRYsQEBDQ+mYEABYvXozFixe3XjOYknwWzVvfxxH7PnjXPB9drBqxZMkSKBQKvPjii+3KTkRExmFVXCZKapowq6cnjt3kPGMeJ/fv34+PPvoIU6ZMgb+/P8rKyvDtt98iLi4Omzdvbr3f/Pnz4ebmhgEDBsDc3Bxbt25F5ZH1cIicjEaldbteIxF1HhbpRHpqy+kC5FY0YHqoOxJuce5nn32Gbt264csvv8TChQthb2+PHj164NFHHwVwcWZgw4YNeOaZZzB9+nT4+vpiwYIF2LNnD86cOdP5LwbAv//9b9jY2ODll19GfX09RowYgTVr1sDS0rL1HK1We8V2MC4uLnC3s8CZA6vw3vZlcHKwx8iRIxEbGwt/f3+d5CYiIv1T26TG53vSMKyrK7ooS255vrGOk35+flAqlXjzzTdRUlICNzc3hIeH48CBAxg0aFDreWFhYVixYgU+/vhjNDc3IyQkBO+89wG+q+6BL/am4d3JvXTyGomobaTLnStv4KY3ElHnkGUZk/59AM0aLX5/cTgUCtO+XuxEdgWmLj+E50aF4JVx3UXHIRLlTn8RcEwno/HvP1Lw0Y4L2PTMEPT1c7z1Hei6Xl93ChsT83Dg9bvgbmd56zsQUUe56Zh+bfcMIhLucFoZkgqq8fiwIJMv0AGgn78TJod744t96cirbF8THyIiMg6V9c34cl86xvbwYIF+h+aP7AK1RouvD2SKjkJEf8IinUgPfbk/Ha625pgc3r5r4YzRaxMudrf9cFuy4CRERCTSF/vSUdusxivjuomOYvACXW0wvqcn1hzNRn3zzRvVEpHusEgn0jOpxTXYfb4Es6IDYWmmFB1Hb/g4WuHJ4cH45WQ+4rMqRMchIiIBSmub8M3BTNzTxxuhntffmpRuz7yhQahqaMGGE3mioxDRJSzSifTMVwcyYKFS4OFoNka72pMjusDdzgJLtp7DLfppEBGREVq+Ow3NGi1eGtNVdBSjERnghN4+DvjmYAa0Wo6tRPqARTqRHimtbcL6E3mY1t8XLrYWt76DibGxUOG50V1xPKsCey7cupsvEREZj4KqBqw+koVp/XwQ7GYrOo7RkCQJjwwJRFpJHfanloqOQ0RgkU6kV1bHZaFZrcW8IUGio+itGZF+8HWywj9/P8/ZdCIiE7JsVypkWcZzoziL3tEm9fGCm50Fvj6QIToKEYFFOpHeaGzRYNXhLIwOdUeIO2cIbsRcpcCLY7rhTF41fjtTKDoOERHpQE55PdYey8HMKH/4OVuLjmN0LFRKPDwwAHsvlCC9pFZ0HCKTxyKdSE9sTMhDWV0zHhsWLDqK3psS4YMubjb4544L0PD6OSIio7dsVyoUCgnP3BUiOorRmhnlB5VCwpqj2aKjEJk8FulEekCrlfHVgQz09LZHdLCz6Dh6T6mQ8PLY7kgtrsWmRHajJSIyZjnl9Vh/IhcPRvnDw95SdByj5W5viTFhHlgXn4vGFo3oOEQmjUU6kR7Ye6EEqcW1eHxYMCRJEh3HIEzs5YkeXvb4ZGcKWjRa0XGIiKiTLN+TBoUk4akRXURHMXoPDvRHRX0Ltp/l5WREIrFIJ9IDKw6kw8vBEpP6eImOYjAUCgmvjOuG7PJ6bOTerkRERimvsgHr4nMQE+UHTwfOone2oSGu8He2xvdHuOSdSCQW6USCnS+swcHUMsweFAgzJf9J3o5Roe7o5WOPZbtToeZsOhGR0Vm+OxUAOIuuIwqFhJlR/jiaUY7U4hrRcYhMFisCIsFWHsqEhUqBmAF+oqMYHEmS8Pyorsgur8emxHzRcYiIqAPlVzbgp+M5+L9IP3g7WomOYzIeiPSFmVLCD0dyREchMlks0okEqqpvQWxCHu4P94GTjbnoOAZpbA8PhHldnE1np3ciIuPxn71pAICn2dFdp1xtLTCupyfWn2ADOSJRWKQTCfTT8Rw0tGgwZ3Cg6CgG6+JseggySuvw6ynOphMRGYPS2iasPZaDqRG+8OEsus49FOWPqoYWbD1dIDoKkUlikU4kiEYr49vDmYgKckYPb3vRcQza+J6e6OZhi093cTadiMgYrDyYiWaNFk+OCBYdxSQN6uKCQBdr/HScS96JRGCRTiTIruRi5FY0YC5n0e+YQiHhuVFdkVpci21n+Kk/EZEhq21S47vDmRjfwxPBbrai45gkSZIwtZ8v4tLLkVNeLzoOkclhkU4kyLeHMuHlYIlxPTxERzEKd/f2Qhc3G3z6Ryq0nE0nIjJYPx7NRnWjGk+NZEd3kab28wEAbEzgNqdEusYinUiAlKIaHEgtxcPRAVBx27UOobw0m36+qAa/JxWKjkNERO3QrNZixf4MRAc7I9zPUXQck+brZI1BwS5YfyIXsswPv4l0idUBkQDfHs6EuUqBmVH+oqMYlXv6eMHf2Rqf703nGwoiIgMUm5iHwupG7ouuJ6b390VWWT2OZ1WIjkJkUlikE+lYdWMLNpzIw319veHMbdc6lEqpwOPDg3EypxJHMspFxyEiotug1cr47940hHnZY0Q3N9FxCMCEXp6wNldifXyu6ChEJoVFOpGObTyRh/pmDeYMChQdxSg90N8XLjbm+O+l/XWJiMgw7D5fjLSSOjw1IhiSJImOQwBsLFSY2MsLv54qQEMz90wn0hUW6UQ6JMsy1hzNRh9fB/T2dRAdxyhZmikxd3Agdp8vQXJhteg4RETURl8fzICXgyXu7u0lOgr9yfT+vqhtUrPfC5EOsUgn0qET2ZVILqzhteidbNagAFibK/HF3nTRUYiIqA3OF9bgYGoZZg8KhBkbquqVgUHO8HG0wjoueSfSGf4WJNKhNUezYWOuxL19vUVHMWqO1uaIGeCPX07mI6+yQXQcIiK6hW8OZsDSTIGZUX6io9BVFAoJ0/r74kBqKQqqOKYS6QKLdCIdqSfTsaAAACAASURBVGpowa+n8nFfuA9sLVSi4xi9R4cFQQbw9YEM0VGIiOgmymqbsCEhD9P6+cLRmg1V9dG0fj6QZe6ZTqQrLNKJdGRTYh4aW7R4aCCXuuuCj6MV7uvrjTVHs1FZ3yw6DhER3cCao9loVmvxyJBA0VHoBgJcbDAg0Anr47lnOpEusEgn0gFZlvHDkWz09nFALx82jNOVJ0cEo75Zg9VxWaKjEBHRdag1WqyOy8awrq4IcbcTHYduYkqEL9JK6nA2n01ZiTobi3QiHUjIYcM4EUI97TGyuxtWHspEYwu3jiEi0jd/JBejsLoRs7ktqd67u7cnzJQSNiVyyTtRZ2ORTqQDa45kw9pcifvC2TBO154c3gWltc3sSktEpIdWx2XBy8ESd3V3Ex2FbsHR2hwjurnjl5P50Gi55J2oM7FIJ+pk1Y0t2HwqH5PDvdkwToDoYGf09XPEl/vT+aaCiEiPZJbWYX9KKWZG+UPFbdcMwuRwbxRVN+FIRpnoKERGjb8RiTrZpoSLDeMejAoQHcUkSZKEp4YHI6usHtvPFoqOQ0REl/xwNBtKhYSYAdx2zVCMCfOAjbkSvyTmi45CZNRYpBN1sp+O56KHlz16+7JhnCjjenrC39ma27EREemJxhYNfj6eg3E9POBubyk6DrWRlbkS43t6YuvpAjSp2euFqLOwSCfqRMmF1TidV4UHIn1FRzFpSoWEOYMDcTyrAqdzq0THISIyedvOFKCivgUPR3OVmaG5L9wb1Y1q7DlfIjoKkdFikU7UiX4+ngszpYTJ4T6io5i8ByJ9YWOuxDcHOZtORCTa6rhsBLnaYFCwi+godJuGhrjCxcacS96JOhGLdKJO0qLRIjYhD6NDPeBsYy46jsmztzTDA5F+2HwqH8U1jaLjEBGZrOTCasRnVeChgf5QKCTRceg2qZQK3NPHCzvPFaGmsUV0HCKjxCKdqJPsTi5GWV0zl7rrkTmDA6HWyvg+Llt0FCIik7Xu0iqzqf04Phqq+8J90KTWYvvZItFRiIwSi3SiTvJzfC7c7Cwwohv3ftUXQa42uKu7O74/ksWGN0REAqg1WsQm5uOu7u5cZWbA+vk7ws/ZCpsS80RHITJKLNKJOkFpbRN2JxdjaoQP937VM48MCURpbTN+PVkgOgoRkcnZl1KC0tomTO/PWXRDJkkSJvf1wcHUUl5CRtQJWD0QdYLYhDyotTLfhOihoSGuCHG3xTeHMiDLsug4REQmZX18HpxtzDGyu7voKHSH7o/whlYGtpzih95EHY1FOlEHk2UZ6+Jz0dfPEV097ETHoatIkoS5gwNxJq8ax7MqRMchIjIZlfXN2JFUhPv6esNcxbeghi7E3Q49vOyxiV3eiTocf0MSdbCkgmokF9ZwFl2PTe3nA3tLFbdjIyLSoc2nCtCs0XJ8NCKTw72RmFOJzNI60VGIjAqLdKIOFpuQBzOlhHt6e4mOQjdgba7CjAF+2H62CEXVvJaOiEgX1sfnItTTDj297UVHoQ5yX7g3JAn45SRn04k6Eot0og6k0cr45WQ+RnZ3hxO71uq1BwcGQKOVsfZYjugoRERGL62kFok5lZjWzxeSxL3RjYWXgxWiAp0Rm5jHPi9EHYhFOlEHiksvQ1F1E+4P9xEdhW4hyNUGw7q6Ys3RbKg1WtFxiIiM2vr4XCgVEiZHeIuOQh1scrgP0kvqcCavWnQUIqPBIp2oA21MyIOdhQqjw9i11hA8NDAABVWN2JVcLDoKEZHR0mhlbDiRh+FdXeFuZyk6DnWwu3t7wkwpcc90og7EIp2ogzQ0a/DbmUJM7O0JSzOl6DjUBmPC3OFpb4nVR7JFRyEiMlqH0kpRWN2IaWwYZ5QcrS9uqffLyXxotFzyTtQRWKQTdZCd54pQ26TG/RFc6m4oVEoFYqL8sO9CCbLK2JmWiKgzrI/Phb2lCmPCPERHoU5yf7gPimuaEJdeJjoKkVFgkU7UQWIT8uDlYInoIBfRUeg2xAzwh1Ih4QfOphMRdbiaxhb8drYQ9/b15iozIzY6zB22FioueSfqICzSiTpAWW0T9l4owX3h3lAo2LXWkHg6WGJsmAd+Op6DxhaN6DhEREZl2+lCNLZoudTdyFmaKTG+p+el/98cS4nuFIt0og6w5XQB1FoZU7jU3SA9FO2PivoWbD9bKDoKEZFRWRefi2BXG0T4OYqOQp3s/ghv1DSpsZvNWInuGIt0og6wMSEPoZ52CPW0Fx2F2mFIF1f4Ollxz3Qiog6UXVaPo5nlmNafe6ObgsFdXOFqa4FYLnknumMs0onuUGZpHRKyK9kwzoApFBJmRPrhUFoZG8gREXWQ9SdyIUnA1H4cH02BUiHh3r5e2J1cgqqGFtFxiAwai3SiO7QpMR+SBNzX11t0FLoD0yN9oZCAn45zNp2I6E5ptTI2JORiSBdXeDlYiY5DOnJ/uA+aNVr8dqZAdBQig8YinegOyLKM2MQ8RAe5wNuRb0IMmZeDFUZ2d8fPx3Oh1mhFxyEiMmjHMsuRU96Aaf05i25K+vg6IMjVBrEJ+aKjEBk0FulEd+BkbhUySutwfwRn0Y3BjAF+KK5pwp7zJaKjEBEZtHXxubAxv9jxm0yHJEmYHO6NuIwyFFY1io5DZLBYpBPdgdiEPJirFJjY20t0FOoAo0Ld4WprgR/ZQI6IqN3qm9XYeroAk/p4wdpcJToO6djkcB/IMrD5JGfTidqLRTpRO6k1Wvx6Kh+jQ91hb2kmOg51ADOlAtP7+2L3+WIUV3MGgIioPbafLURdswbT+nFvdFMU5GqDvr4O7PJOdAdYpBO104HUUpTWNmNyOK+3MyYzBvhBo5Wx7kSu6ChERAZpfXwe/JytMCDQWXQUEmRyuA/O5lcjtbhGdBQig8QinaidNiXmw95ShbtC3URHoQ4U5GqDgUHOWHssB7Isi45DRGRQ8isbcDCtFFMjfKFQcG90U3VPXy8oJLCBHFE7sUgnaof6ZjW2ny3EpD5esFApRcehDhYT5YessnrEpZeLjkJEZFA2JuRBlsGl7ibO3c4SQ0JcselkHj/wJmoHFulE7bAjqQj1zRoudTdSE3t5wc5ShbXHskVHISIyGLIsY318LqKCnOHvYi06Dgk2OdwHOeUNOJFdKToKkcFhkU7UDpsS8+HtYIkoXm9nlCzNlJgS4YOtZwpRVd8iOg4RkUFIyKlEemkdpnMWnQCM7+kBC5UCm9hAjui2sUgnuk1ltU3Yd6EE94Z783o7IzZjgB+a1Vp2pyUiaqP18bmwNFNgYm/ujU6AnaUZxoR54NdTBWjRaEXHITIoLNKJbtOW0wVQa2Xcz6XuRq2ntwN6+zhgzdFsXk9HRHQLjS0abD6Zjwk9PWHHbUnpksnh3iiva8aB1FLRUYgMCot0otsUm5CH7h52CPOyFx2FOtmMAX5ILqzB6bwq0VGIiPTaznNFqG5UY3p/P9FRSI+M7O4OByszxCZwVRrR7WCRTnQbssvqcSK7EvdHcBbdFNwX7g1LMwV+PJYjOgoRkV5bH58LLwdLDOriIjoK6RFzlQL39vXCb2cKUd3IHi9EbcUineg2XG5+cl+4t+AkpAv2lmaY1NsbvyTmo75ZLToOEZFeKq5uxN4LJZjazwdK9mqhq0zv74cmtRZbThWIjkJkMFikE7WRLMuITcxDVJAzfBytRMchHYmJ8kNtkxq/8s0FEdF1bUjIg1YGl7rTdfX1dUBXd1usi88VHYXIYLBIJ2qjs/nVSCupY8M4ExMZ4IQubjZYyyXvRETXkGUZ6+JzERnghCBXG9FxSA9JkoTp/X0Rn1WBtJJa0XGIDAKLdKI2ik3Ig5lSwt3cWsakSJKEmAH+iM+qwIWiGtFxiIj0ysncKqQW12J6f+6NTjc2JeLipRDrOZtO1CYs0onaQKOV8cvJfNzV3R2O1uai45COTe3nAzOlxNl0IqKrrIvPgaWZAnf38RIdhfSYu70lRnRzw4YTedBoua0p0a2wSCdqg8NpZSiuaWJXdxPlYmuBcT08seFELprUGtFxiIj0QmOLBr8kXtwb3Z57o9MtTO/vi8LqRu6ZTtQGLNKJ2iA2MQ92FiqMCnUXHYUEiYnyQ0V9C7afLRIdhYhIL+xI4t7o1Hajwy7umc4GckS3xiKd6BYaWzT47UwhJvTyhKWZUnQcEmRIF1f4Ollh7bFs0VGIiPTCuvhceDtYYjD3Rqc2sFApMTncG9vPFqKqgXumE90Mi3SiW9h+thC1TWpM4VJ3k6ZQSJgR6YeDqWXIKqsTHYeISKjCqkbsTynBtP6+UHBvdGqjB/r7oVmtxeaT+aKjEOk1FulEt7AuPhe+TlaIDuZMgambHukLhQT8yAZyRGTiNl7aG31aP3Z1p7br5WOP7h52XPJOdAss0oluIr+yAQdSSzGtH2cKCPBysMLoMA/8fDyHDeSIyGTJsoyfj+dgQKATArk3Ot2Gy3umJ+ZUIrWY25oS3QiLdKKb2HAiF7IM7v9KrWZFB6C0thm/nSkUHYWISIgjGeVIL61DzAB/0VHIAN1/ac/0nzmbTnRDLNKJbkCWZayLz0V0sDP8nK1FxyE9MTTEFYEu1lh1OEt0FCIiIdYczYadpQp39+be6HT73OwscFd3N2w8kQe1Ris6DpFeYpFOdAPHsyqQWVbPrWXoCgqFhIejA3A8qwJJ+dWi4xAR6VRFXTO2nSnE1AgfWJlzxxNqn+n9fVFc04T93DOd6LpYpBPdwLrjubA2V2JiL0/RUUjPTO/vCwuVAquPcDadiEzLhoQ8NKu1mDmQS92p/UaFesDJ2gzrjnPJO9H1sEgnuo76ZjV+PZWPSb29YGOhEh2H9IyjtTnu6+uN2IQ8VDdyr1ciMg2yLGPN0WyE+zki1NNedBwyYOYqBSaH+2BHUhEq65tFxyHSOyzSia5j2+lC1DVr2DCObmj2oEDUN2uw8USe6ChERDoRn1WB1OJaPBjFWXS6c9P7+6JZwz3Tia6HRTrRdaw9loMAF2tEBTmLjkJ6qrevA/r6OWJVXBZkWRYdh4io0/1wNBu2Firc05cN4+jO9fJxQJiXPbu8E10Hi3Siq6QU1eBoZjlmRvlDkrg3Ot3Y7OgApBbXYl8KG98QkXGrqm/BllMFmBzuDWtzXgZGHWN6f1+cyq3C+ULumU70ZyzSia6y5mgOzJQSl7rTLd3b1xvudhZYsT9ddBQiok617kQumtRaPMiGcdSB7g/3hkoh4afjOaKjEOkVFulEf9LYosH6E7kY39MTrrYWouOQnjNXKTBncCD2p5RyFoCIjJZWK2PV4Uz0D3BCT28H0XHIiLjYWmBsDw9sTMhDk1ojOg6R3mCRTvQnW08XoKqhhU1xqM0ejPKHpZkCXx3gbDoRGaf9qaXILKvH7EEBoqOQEZoxwA/ldc3YkVQkOgqR3mCRTvQna45mI8jVBoO6uIiOQgbCycYc0/v7IjYhHyU1TaLjEBF1uO8OZcLV1gITe7FhHHW8YV3d4ONohbXHuOSd6DIW6USXJBdW41hmBWZG+bFhHN2WeUOC0KLV4rvDmaKjEBF1qJzyeuw6X4wHo/xgruLbRup4SoWEByJ9sT+lFDnl9aLjEOkF/rYluuSbA5mwNFPg/yL9REchAxPsZovxPTyx8lAmqhpaRMchIuowq+OyoJAkPDiQS92p8zwQ6QdJAn5mAzkiACzSiQAA5XXNiE3Mw9R+vnC0NhcdhwzQc6NDUNOoxreHMkVHISLqEI0tGqw9noPxPT3g6WApOg4ZMR9HK4zo5oafjudCrdGKjkMkHIt0Ily8Fr1JrcUjgwNFRyED1dPbAWPCPPDVgQzUNHI2nYgM3y8n81FZ34JZ0YGio5AJiBngh8LqRuxLKREdhUg4Fulk8lo0F68lHtbVFV097ETHIQP2/OgQVDW04LvDWaKjEBHdEVmW8c3BTHTzsEV0sLPoOGQCRoV6wNXWHD8e5ZJ3IhbpZPK2ni5AUXUT5g0JEh2FDFwfX0fc1d0NK/ano65JLToOEVG7HUorw7mCajw6NIjNVEknzFUKTOvniz+Si1Fc0yg6DpFQLNLJpMmyjBX7MxDkaoMR3dxExyEj8Pzorqiob8HqOM6mE5HhWrE/Ha625pgc7iM6CpmQGQP8oNHKWBefKzoKkVAs0smkHUgtxem8KjwxPBgKBWcK6M5F+DtheDc3fLGPs+lEZJhSi2uw+3wJZg8KhKWZUnQcMiHBbraICnLG2mM5kGVZdBwiYVikk0lbvjsNHvYWmNqPMwXUcV4a0xVldc1YsT9DdBQiotu2Yn8GLFQKPDTQX3QUMkExA/yQVVaPuPRy0VGIhGGRTiYrIbsCh9PL8NjQYFioOFNAHSfC3wkTenrii31pKK1tEh2HiKjNSmubsCEhD9P6+8LF1kJ0HDJBd/f2gr2lCmuOZouOQiQMi3QyWcv3pMHBygwzOVNAneDVCd3RqNZi2a5U0VGIiNps1eEsNKu1bKZKwliaKTG1ny+2nSngB91kslikk0k6X1iDHUlFmDM4ELYWKtFxyAh1cbPFjAF++P5IFrLK6kTHISK6pcYWDVbHZWF0qDtC3G1FxyET9tBAf7Ro2ECOTBeLdDJJn+y8AFsLFR4ZHCg6ChmxF0d3hUqhwNLfL4iOQkR0Sz8dz0FZXTMeHx4sOgqZuK4edogKcsYPR7Kh1bKBHJkeFulkcs7kVWHbmULMGxoEJxtz0XHIiLnbW+LRoUHYfDIfp3IrRcchIrqhFo0W/92bjn7+jhgY5Cw6DhEejg5Adnk9DqSWio5CpHMs0snk/PP383CwMsNjw3i9HXW+J0cEw9nGHO9vPcftZIhIb/2SmI+8ygY8c1cIJIlbkpJ443t6wMXGHKvjskRHIdI5FulkUuKzyrH7fAmeHBEMe0sz0XHIBNhZmuGlMV0Rl16OHUlFouMQEV1Dq5Xx+d40hHraYVSou+g4RAAAC5USD0T64Y/kYhRUNYiOQ6RTLNLJZMiyjKXbL8DV1hxzeS066dDMKH90cbPBkm3JaFZrRcchIrrC70lFSC2uxfyRXTiLTnrlwSh/aGUZa4/liI5CpFMs0slk7EouxuH0Mjx7VwiszdnRnXRHpVRgwaQwZJTW4fsjXLZHRPpDlmV8vicVAS7WmNTbS3Qcoiv4u1hjeFc3/Hg0B2oNP+Qm08EinUxCi0aL97eeQ7CrDR6KDhAdh0zQXd3dMTTEFZ/sTEFlfbPoOEREAICDqWU4mVuFJ4d3gUrJt4Wkfx4a6I/C6kb8kVwsOgqRzvC3MZmENUezkVZShzfvDoMZ34SQAJIkYcGkMFQ3tuDTXami4xARAQCW70mFu50FpvX3ER2F6LpGhbrDy8ES3x/JFh2FSGdYrZDRq2powcc7LmBwFxeMDmNDHBInzMse/9ffD98dzkRmaZ3oOERk4hKyK3AorQxPDA+GhUopOg7RdamUCsQM8Me+CyXILqsXHYdIJ1ikk9H7eMcFVDa0YMGkMDbEIeFeGdcNZkoFPtiWLDoKEZm45XvS4GhthplR/qKjEN3UjAF+UCok/HCUs+lkGlikk1FLzKnEt4czMSs6AD29HUTHIYK7vSXmj+iC384W4kh6meg4RGSizhfWYEdSEeYODoSNBZupkn7zdLDEmDB3/HQ8B01qjeg4RJ2ORToZrRaNFn/dcBrudhZ4dXx30XGIWj02LBheDpZ4b8s5aLWy6DhEZIL+vSsFthYqbklKBuOhgQEor2vGb2cKRUch6nQs0slofX0gA+cKqvHufb1gZ2kmOg5RKytzJV6b0B2n86oQm5gnOg4RmZgLRTXYeroAcwYHwNHaXHQcojYZGuKKABdrfHeYW5mS8WORTkYpu6weH++8gHE9PDChl6foOETXmNzXB318HfCP386joZlL94hId/79RwqszZR4bGiw6ChEbaZQSJg9KBDxWRU4nVslOg5Rp2KRTkZHlmUsiD0NlUKBdyf3FB2H6LoUCgkL7g5DYXUjVuxPFx2HiExESlENtpwuwJzBgXCy4Sw6GZYHIn1hba7EykOZoqMQdSoW6WR0fjmZj/0ppXh1fHd4OViJjkN0QwODXTChpyc+35uG4upG0XGIyAT8e1cqrMyUeGwYZ9HJ8NhbmmF6f19sPpmP0tom0XGIOg2LdDIqFXXNWLw5CeF+jng4OkB0HKJbemNiKFo0Wvzz9wuioxCRkUstrsGvp/IxZ3AgnDmLTgZq9qBANGu0WHOE27GR8WKRTkbl/a3nUNXQgiVTe0Op4J7opP8CXW0we1AgforPQVJ+teg4RGTEPr00i/44Z9HJgIW422J4NzesistCi0YrOg5Rp2CRTkbjUFopfo7PxWPDghHmZS86DlGbPT+qKxyszPD3rUmQZW7JRkQdL7W4Fr+czMfsQZxFJ8P3yOBAFNc0YRu3YyMjxSKdjEJjiwYLNp6Bv7M1XhjdVXQcotviYG2GF0Z3xcHUMuw+Xyw6DhEZoWW7UmCpUuLxYUGioxDdsRHd3BDoYo2VBzNERyHqFCzSySh8tjsVGaV1+PuUXrAyV4qOQ3TbHo4OQLCrDd7fmgyNlrPpRNRx0kouzaIPDoCLrYXoOER3TKGQMGdwIE5kV+JkTqXoOEQdjkU6GbwLRTX4z940TInwwbCubqLjELWLmVKBl8d1Q2pxLbacLhAdh4iMyLJdqbBQKfEEr0UnIzK9vy9szJX4ltuxkRFikU4GTauV8dcNp2FjocJbk8JEx2m32NhY9OnTBxYWFggKCsJHH3100/NfeuklSJKEv/zlL1ccT05OxsCBA+Hg4ICYmBjU1tZecfu+ffvg4+NzzfHrWblyJSRJuu65ixYtgqura+v3mZmZkCSp9cvGxgZdunTBQw89hP37919z/7lz5yIyMvKWGUzN3b280NXdFv/+I4Wz6UTUIdJLarEpMQ+zB5n2LDrHWeNjZ2mGByL9sPlUPopruI0pGRcW6WTQvj+ajfisCiy4O8xg33wcPHgQU6dORVRUFDZv3ox58+bh9ddfxyeffHLd85OSkvDVV1/B3v7a5nhz585FSEgIfvrpJyQlJeH9999vvU2r1eKFF17AkiVLYGtr2ymvZenSpTh8+DC2bt2Kt99+G2VlZRg+fDjefffdTnk+Y6NQSHhhTFfOphNRh7k8i/74cNOdRec4a7xmDwpAi0bGmiM5oqMQdSiV6ABE7VVY1YgPtyVjSIgLpvf3FR2n3RYvXowhQ4ZgxYoVAIBx48ahsrISixcvxtNPPw1z8yu78D733HN44YUXsGrVqiuO19bW4siRI9i8eTPc3NxQWVmJpUuXtr6B+Prrr2FmZoZZs2Z12mvp3r07oqOjAQAjRozA3LlzsXDhQixatAgjRozAyJEjO+25jcXF2fQUfPpHCib19uJWgkTUbukltYhNzMOjQ4PgaqAfZHcEjrPGK9jNFiO7u2H1kSzMH9kF5irOP5Jx4N9kMlgLN51Bi0aLv9/fG5JkuIVMYmIixo4de8WxcePGoaKiAocPH77i+Lp165CcnIw33njjmsdpbm4GAFhZWQEArK2tW49VV1fjrbfewr/+9S+d/6zeeecdeHt74z//+Y9On9dQXZ5NTymuxVbOphPRHVi2OxXmKgWeGN5FdBShOM4at7mDA1FS04RtZzhmkvFgkU4G6bczBfg9qQgvjumGQFcb0XHuSGNj4zWf4l/+/ty5c63HGhoa8Morr+CDDz6Ajc21r9nZ2RmBgYH49NNPUV5eji+++KL1erS//e1vGDNmDAYNGnTb+TQaDdRq9RVfWq22zfdXKpUYNWoU4uLibvu5TRWvTSeiO5VRWofYhDw8PDAAbnamO4sOcJw1dsO7uiHY1QYr2UCOjAiXu5PBqWpowcJNZxHmZY/HjGC/15CQEBw7duyKY0ePHgUAlJeXtx5bsmQJvLy88PDDD9/wsZYvX44HHngAb775Jrp27YrPPvsMqampWLFiBU6fPt2ufI6Ojtc97uLi0ubH8PX1RVFRUbue3xRdnk1/9ocE/HamEJP6eImOREQGZtmuS7PoI0z3WvTLOM4aN4VCwuxBAVi0OQmJOZUI97v+z5PIkHAmnQzOh78lo7S2CR9O6w0zpeH/FX7qqacQGxuLL7/8EhUVFdi+fXtr11mF4uLry8jIwNKlS2+5jG7ixIkoLi7G+fPnce7cOfj7++Pll1/GSy+9BF9fX3z22Wfw9/eHv78/li9f3qZ8+/btw7Fjx674evzxx2/rNcoyZ4Nv18ReXgh2tcHne1P58yOi25JZWofYxDw8NDAA7naWouMIx3HW+E3r7ws7CxVW7E8XHYWoQ3AmnQzK0Yxy/HAkG48NDUIfX+P4pHTevHk4efIk5s+fjyeeeALW1tb48MMP8dxzz8HT0xMA8MYbb2DixIno3r07KisrAVzsItvU1ITKyko4ODi0vqmwtrZGt27dAAA7duzAyZMnsXbtWpw8eRJvv/02Dh06BAAYNGgQhg4dij59+tw0X0RExDVdan/99dfbeo15eXnw8PC4rfuYOqVCwpMjgvH6+tM4kFqKYV3dREciIgOxbHcqVJd+hxDHWVNgZ2mGBwf648v96cgpr4efs7XoSER3xPCnIclkNKk1+OuGU/B1ssLL47qJjtNhlEolli1bhpKSEpw6dQpFRUWtnVsv/3n+/Hls2LABTk5OrV85OTlYtmwZnJyckJeXd83jajQavPTSS/jHP/4BKysr7NmzB6NGjUJoaChCQ0MxevRo7N27t9Nfn1qtxq5du9p1nZ6puz/CBx72Fvh8T5roKERkINJLarHhRC4ejuYs+mUcZ03DI0OCoFRInE0no8CZdDIY/92bjrSSOqx8ZACszY3vr+7lNwXAxWveBg8ejNDQUADAihUrUFtbe8X5MTExGDFiBObPnw83t2tnWT///HM4OTlhxowZrcfq6+tb/7uurk4ny+MWL16M/Px8PPXUU53+XMbGiLXfswAAIABJREFUQqXEY0OD8fet53idHRG1ycc7U2BppsT8kabd0f16OM4aN08HS9wf7oO1x3PwwphucLYxv/WdiPSU8VU6ZJRyyuvx2e5UTOrthZHd3UXH6VBxcXE4cOAAwsPDUV1djTVr1mD79u04cOBA6zmXu8f+maWlJfz8/n979x0fRbX/f/x1dtN7SEIJJQkJXQQEAUFQFAUs2AtXr3ptV7lgx64E9Xft196uispVka9dsKACiiACgnQIBEIIIISQRnrZ+f2xASH0kM1sNu/n47GPJLszy+dzdtgzn5kzZ9oe8J6oubm5TJgwgenTp+95bvDgwdx9991MnDgRy7KYOXMmTzzxRL3mkpaWRmxsLBUVFWRkZPDRRx/x3Xff7bl/qxy9Uf3a8dLMdbz+03pe/3tvu8MRES+2+s9Cpi7dyr+GJDfp+6LXpn626bhxcHs+XrSZSfM2cttQ3xl1KU2PinRpFFK/WonTYXjwnC52h1Lv/P39mTJlCqmpqTgcDgYNGsTcuXPp3r17nd8zNTWVkSNHcsIJJ+x5rlevXjz11FM88MADADzzzDP06NHjmOPf21133QW4d2xatWrFSSedxOzZsxk0aFC9/jtNSVigH1cPSOTlWemkZxeR0jzs8CuJSJP03A9rCQ/y48ZBOou+N/WzTUeHFuEM7dKc937dyD8HJxMc4LQ7JJE6MYcZhtO0p4oUr/DDqu3cMOl37j+rMzcO1o6HND07i8oZ+ORMzj0+nqcvqd8dPmlUDj7l9JFRn+7Dlmblc94rc7njjI7ccnoHu8MRsc3Cjblc8vo8HjmvG1edlGh3OCIHc8g+XRPHiVcrragm9auVdGwRxj8GNv57oovURUxYIJf0bsuXS7aSvavM7nBExAs9+8NaokP8ufZk9ZXStPVJiOaEdlG8+csGqqpddocjUicq0sWrvfpTOlvyS3nkvON84p7oInX1j4GJVFS7eP+3TXaHIiJeZkFGLrPX7uDmU5MJC9SVjNK0GWP45ynJZOWW8u2KbXaHI1InqnrEa2XuLOaNnzdwfs94+rePsTscEVu1jwvj9M7N+eC3TMoqq+0OR0S8hGVZPPN9GnHhgfy9f6Ld4Yh4hTO6tKB9bCiv/bS+QWbYF6lvKtLFaz329Wr8nIb7zvK9yeJE6uLak5PYWVzBV0u22h2KiHiJOek5LMjIZexpKZokS6SGw2G46dRkVv1ZyKy0bLvDETlqKtLFK/2ybgc/rNrOmNNSaBERZHc4Il5hQHIMnVuGM3Fuhs4MiEjNWfS1tI4K5rIT29odjohXuaBXa9pEB/PijHT1mdLoqEgXr1NZ7WLC1FUkxIRwnSbAEdnDGMO1JyexZtsu5qbvtDscEbHZtyu2sTQrn1tP70Cgn86ii+zN3+ng5lOTWZKVz5z0HLvDETkqKtLF60yal0l6dhEPnd1VOx0itYzsEU9sWAAT52bYHYqI2Kiy2sXT09Po2CKMi3q3sTscEa90ce82tIwI4qUZ6XaHInJUVKSLV8kpKuf5H9cyuGMcp3dpbnc4Il4nyN/Jlf0TmLkmm/U7iuwOR0Rs8tGCTWTkFHPP8M44HYe83a5IkxXo5+Sfp7RnwcZcftugEWjSeKhIF6/yzPQ0SiuqeficrhijnQ6RA7myfwIBfg7e0dl0kSapqLyKF2aso29SM07rrAPaIocyqm87YsMCeWnmOrtDETliKtLFayzfXMCU37O4ZkAiKc3D7A5HxGvFhgVyfs94Pl20hfySCrvDEZEG9ubsDeQUVXDfiM46oC1yGEH+Tm4cnMTc9J0sysyzOxyRI6IiXbyCZVmkTl1JTGgAtwztYHc4Il7v2pOTKK2s5sMFm+wORUQaUPauMt78ZQNndW9Jr3bRdocj0ihc0S+B6BB/nU2XRkNFuniFL5dsZVFmHncP60xEkL/d4Yh4vc4tIxiYEsOkXzOprHbZHY6INJAXZ6yjosrFuGGd7Q5FpNEIDfTj+kHt+SltB0uy8u0OR+SwVKSL7YrLq3j829Uc3yaSizVDrcgRu+7kJLYVlvHN8j/tDkVEGsD6HUVMXpDFqL7tSIoNtTsckUbl6gGJNAsN4Onpa+wOReSwVKSL7V6amc72wnImjOyGQzPUihyxUzs2p31cKG/PycCyLLvDEREPe2zaKkL8ndxyui4LEzlaYYF+/GtICnPTdzJnne6bLt5NRbrYasOOIt6es4GLe7fRtXUiR8nhMFw7MIllmwtYuFGT4Yj4sllrspmVtoNbTu9AXHig3eGINEpX9GtH66hgnp6+Rge3xaupSBfbWJbFhKmrCPJzcs9wXVsnUhcXndCG6BB//jt7g92hiIiHVFS5ePTrVbSPDeXqAYl2hyPSaAX5O7l1aAeWbi5g+sptdocjclAq0sU2P67O5ue1O7jtjI46KyBSR8EBTv7eP4EZa7azYUeR3eGIiAdMmreRDTuKeeicrgT4addN5Fhc2Ks1yXGhPPP9Wqo08ap4KX3Tiy3KKqt5dNoqOjQP46qTEuwOR6RR+/tJifg7Hbw9J8PuUESknuUUlfPCj+s4tVMcQzo3tzsckUbPz+lg3LBOpGcX8dkfW+wOR+SAVKSLLd6cvYFNuSVMGNkNf6c2Q5FjERceyIW9WvPJos3kFlfYHY6I1KNnv0+jtLKaB8/uancoIj5jWLeW9GgTyQs/rqO8qtrucET2o+pIGtzmvBJe+Smds7q3ZEBKrN3hiPiE6wclUV7l4v3fMu0ORUTqyYotBXy0MItrBiSS0jzM7nBEfIYxhruHd2ZLfin/m6d+U7yPinRpcP/+ZjUAD+isgEi9SWkezpBOcUyat5GySp0VEGns3JOrrqRZSABjdcs1kXo3MCWWwR3jeHHGOo1CE6+jIl0a1Nz0HL5Zvo1/nZpC66hgu8MR8Sk3DGpPTlEFX+gaO5FG76ulW1m4MY9xwzoRGexvdzgiPunBs7tQXFHN8z+utTsUkX2oSJcGU1ntIvWrlbRrFsINg9vbHY6IzzkpOYaurSJ4a04GLpfu/yrSWBWUVPLotFUc3yaSS/q0tTscEZ/VsUU4V/RrxwfzN7F2+y67wxHZQ0W6NJj3ft3IuuwiHjqnK0H+TrvDEfE5xhhuHNye9Owifly93e5wRKSOnvhuNXkllTx+YXecDmN3OCI+7bahHQkNcPLotFVYlg5wi3dQkS4NYkt+Kf/5YS1DOsUxtItuISPiKecc34p2zUJ4eVa6djZEGqEFGblMXpDFdScn0S0+0u5wRHxes9AAbh3akV/W5fDDKh3gFu+gIl08zrIsHv5iBZYFj5x3HMborICIp/g5HYw+NZllmwuYvS7H7nBE5CiUV1Vz/+fLaR0VzG1DNVmcSEO56qQEOrUIZ8LUVZRWaPJVsZ+KdPG4b1dsY8aabO48syNtm4XYHY6Iz7vwhDbERwbx0ox1Opsu0oi8/tMG0rOLeOyC4wgJ8LM7HJEmw9/p4NHzj2NLfimvzEq3OxwRFeniWQWllYz/aiXHtY7gmgGJdocj0iQE+Dm46dRkfs/MY35Grt3hiMgRWL+jiFdmpXNuj3iGdNJlYSINrW9SMy7s1Zr/zt7Ahh1FdocjTZyKdPGoJ79bw86ich6/4Hj8nNrcRBrKpX3aEhceyMszdUZAxNu5XBb3fbacIH8HD5/T1e5wRJqs+87qQqC/gwc+X6G7pIitVDWJx/y6PocP52/iHwOT6N5Gk9+INKQgfyc3DmrPnPQcFm/KszscETmE9+ZtZEFGLg+e05W48EC7wxFpsuLCA3ngrC7M27CTDxZssjscacJUpItH7CqrZNzHy0iMCeHOMzvaHY5Ik/S3fu1oFhrAcz+stTsUETmIjJxinvxuDad1bs4lvdvYHY5Ik3fZiW0Z1CGWx79ZTVZuid3hSBOlIl084t/frObPglKevbSHJr8RsUlooB+jT03ml3U5/Jqumd5FvE21y+Kuj5cS4HTw+IXddfcTES9gjOGJi47HYQz3fLpME7CKLVSkS72blZbN5AVZ3Dg4md4JzewOR6RJu7J/AvGRQTw5PU07GiJeZuKcDBZl5jHhvG60iAiyOxwRqdE6Kpj7z+rCr+t38sF8DXuXhqciXepVfkkF93yyjI4twrj9DN3jVcRuQf5ObhvakaVZ+Uxfud3ucESkRnp2EU9/n8YZXVtwfs/WdocjIrWM6tuWk1Pcw9435hTbHY40MSrSpd5YlsUDX6wgt7iC/1zak0A/p90hiQhw4QmtSY4L5Znv06iqdtkdjkiTV1Xt4s6PlxIa4OTfF2iYu4g3Msbw5MXuuxP968PFlFVW2x2SNCEq0qXefLQwi6+X/cntZ3TkuNaazV3EW/g5HYwb1pn07CImL8yyOxyRJu+/v2xgaVY+j5x3nGZzF/FiraOCefaSHqzcWsj/+3q13eFIE6IiXerF2u27SP1qJYM6xHLzKcl2hyMitQzr1oKT2sfw7Pdp5BVX2B2OSJO1cmsBz/+wjrO6t+Sc41vZHY6IHMbQri24YVAS//stk2nLttodjjQRKtLlmJVWVDPmw8WEB/nzn0t74nBo2J6ItzHGkDqyG7vKqnj2hzS7wxFpksoqq7n1oyVEhfjz/87XMHeRxuLu4Z3p1S6Kez9druvTpUGoSJdj9si0VazdXsRzl/XQsD0RL9apZTh/75/Ah/M3sXJrgd3hiDQ5j3+zmvTsIp69tAfRoQF2hyMiR8jf6eClUb1wOgw3TPqdwrJKu0MSH6ciXY7Jx79nMXnBJm4+NZlBHeLsDkdEDuP2oR2JCgngoS9WUO3SLdlEGsqstGzem5fJtQOT1F+KNEJtokN47YoTyMgpZsyHf2giVvEoFelSZ0uz8nngixUMTInhzjM62h2OiByByBB/Hjy7C4s35fO/eRvtDkekScgpKmfcx8vo3DKcu4d3sjscEamjASmxPHb+ccxeu4NHpq2yOxzxYSrSpU527Crnn/9bRPPwQF4edQJ+Tm1KIo3FBb1ac0rHOJ6ankZWbond4Yj4NMuyuPfTZRSWVvL85T0J8tftSUUas8v7tuOGQUlMmpfJaz+ttzsc8VGqrOSoVVS5+NcHi8kvreC/f++j6+pEGhljDP++sDsGuP/z5ViWhr2LeMr7v2Xy4+ps7h7eic4tI+wOR0TqwX0jujCyRzxPfreGD+dvsjsc8UEq0uWoWJbFQ1+sYMHGXJ66uAdd47XDIdIYtY4K5t4RnfllXQ4fLtAOhognrNxawKNfr+bUTnFcOzDJ7nBEpJ44HIZnL+3BqZ3ieOCL5UxdqluzSf1SkS5H5YUZ65jyexa3nJbCyB7xdocjIsfgin4JDOoQyyNTV5G2bZfd4Yj4lOLyKsZ++AfRIf48e0kP3Z5UxMf4Ox28dkVvTkxoxm1TlvDlki12hyQ+REW6HLH/W5jF8z+u4+LebbhdE8WJNHoOh+E/l/YkPMifMR8uprSi2u6QRHzGQ1+sYOPOYp6/rBcxYbo9qYgvCg5w8s4/TuTExGhum7KE/1uYZXdI4iNUpMsRmZWWzX2fL2dQh1gev7A7xuiMgIgviAsP5D+X9mBddhETpq60OxwRn/DJos189scWxp7WgZOSY+wOR0Q8KDTQj3eu6cvJKbHc/ekyXvtpveZ6kWOmIl0Oa0FGLqPfX0znluG8dmVv/DWTu4hPGdwxjtGnJvPRwiwmzdtodzgijVp6dhEPfbGC/u2bccvpHewOR0QaQHCAkzev6sPZx7fiye/WcPcny6io0n3Upe787A5AvNuizDz+8c4C4qOCePcffQkL1CYj4ovuPLMTadt2MWHqKpJiQxnUIc7ukEQanaLyKm56fxHBAU5euLwXTl2HLtJkBPk7eenyXiTHhvLizHQ25Zbw+pW9dRckqROdEpWDWpqVzzUTFxAXHsiHN/QnLlzX1In4KqfD8MKoXnRoHsboDxaTnq2J5ESOhmVZjPt4KRt2FPHyqF60iAiyOyQRaWAOh+GOMzvx/GU9+WNTPiNfmcPyzQV2hyWNkIr0RmDVqlWcfvrphISEEB8fz8MPP0x19ZFP8ORyuejTpw/GGKZNm7bPaxUVFTzyyCOkpKQQHBxMSkoK48ePZ86aP7nirflEhfrz4Q39tbMh0sh9+eWXdO/enaCgILp27cqUKVP2WyYs0I+3ru5DoJ+Tv7+9gP+88gbGmP0er7/++p51/vzzT8aNG0ePHj0ICwujbdu2XH311WzdqtvRSONSl7525cqVDB8+nPj4eAICA3lz9AhaL59EUmjlPsuNHz+e7t27ExERQXh4OH369Dng/0ERaZxq97Hla+fw0T/7U11tcdFrv/K/3zL3u079k08+YcCAAcTExBAUFESnTp147LHHqKio2O/9ly9fzjnnnENkZCTh4eH07duXRYsWNVR6YgONXfZyeXl5DB06lK5du/Lll1+yfv167rzzTlwuF4899tgRvcdbb73F5s2bD/javffey+uvv85jjz1Gr169WLx4Mfc/8CAvfruEEy67nQ+u70/LSBXoIo3ZnDlzuOiiixg9ejQvvvgi33zzDaNGjSI6Opozzzxzn2XbRIcw6dq+XP7fefx39gYAZs6cSXBw8J5l2rdvv+f3RYsW8fnnn3P99dfTr18/tm/fTmpqKgMGDGDFihWEhYU1TJIix6CufW1BQQFJSUmceMb5vLkoj27h5ayf/i5nnXUWCxcuxM/PvZtVWFjINddcQ9euXXE6nXzyySdcfvnlOJ1OLr744oZKU0Q84GB97Hfffce0W07ljv9bwkNfrGBBRi7/74LjiAjyB2Dnzp2cdtppjBs3jqioKBYsWEBqairbtm3j5Zdf3vP+S5YsYdCgQZx33nl7Du4tXLiQ0tJSW/KVBmJZ1qEeYrN///vfVlRUlFVQULDnuSeffNIKDg7e57mDyc3NtWJjY6233nrLAqypU6fu83qLFi2sO+64Y8/fny7KsiL7nm8FhEVbO3aV1V8iImKbM8880xoyZMg+z40YMcIaOHDgQdf5fWOu1WrkHRZgZfyZc9Dl8vLyrMrKyn2eS0tLswDr3XffPbbApbbD9dnq0+voWPraDTuKrONTp1vDnvvZKi6vtL7//nsLsBYtWnTI9QYMGGCde+659RK/iNjncH1sdbXLennmOivp3mlW/3//aP2cln3Q97r//vutyMhIy+Vy7XmuX79+1qhRozwTvNjpkH22hrt7uW+//ZZhw4YRERGx57nLL7+c0tJSfv7558Ou/9BDDzFw4EBOP/30A75eWVlJZGQklmXxnx/Wcsf/LaVdy1gigv2I1X1dRRq98vJyZs2axaWXXrrP85dffjnz5s2joODA18r1TojmHwMSAbjq7QVsyT/wEfuoqKg9Zwt369ixIyEhIRryLo1GXfva/JIKrnt3IU6H4b9/70NIgB8xMe5brh1oyOreYmJiDruMiHi3I+ljHQ7Dv4ak8NnogYQEOLlq4gLu+2w5ReVV+71f7e+FVatWMX/+fMaOHevxXMS7qEj3cmvWrKFz5877PNeuXTtCQkJYs2bNIdddtmwZEydO5JlnnjnoMtdffz2vv/EGF6e+w/PfLqN/SDY7Fkxj7Jgx9RK/iNhr/fr1VFZW7vc90qVLF1wuF2vXrj3ouh1ahAPw6+N/o01MGO1TOvDGG28c9t9ctmwZJSUldOzY8diCF2kgdelrK6pc3Pz+YjbnlfLaFb1oGe5HWloa9957LyeeeCJ9+/bdb52qqiry8/P54IMP+P7777nppps8ko+INIyj6WN7to3i61sGcePg9ny0cBPDnpvNz2t3UF1dTUlJCXPmzOHFF1/k5ptvxhj3nSHmz58PuC/J6dGjB35+fiQnJ/P22283XJJiCxXpXi4vL4+oqKj9no+OjiYvL++Q644dO5YxY8aQkpJy0GX+cdsDONr347NHriPruUuYMv5aLrzwQh5++OFjjl1E7Lf7e6L290h0dPQ+rx9Iq1atePTRR3n1zYmkXPEI+WGJ3HTTTTz33HMHXcflcnHrrbfSoUMHRo4cWQ8ZiHje0fa1lmXxwOfLmbdhJ09c1J3xo68kMDCQzp07k5uby7Rp03A49t3F+u233/D39yc6OpprrrmGF154gfPPP99jOYmI5x1tHxvk7+T+s7rwyU0nEejn4OqJCwgMDiE0NJRBgwZxyimn8PTTT+9Zftu2bQBcddVVXHHFFfzwww8MHz6c66+/nm+++caTqYnNNHGcj/roo49IS0tj6tSpB3zdsiw+XbyFm+9+mLxFP3L7w49z/ukDWLp0KQ899BAxMTE88sgjDRy1iHiTYcOGMWzYMACGnzmM6947mXlvlfNw6iPceuut+xUhAPfddx/z5s3j559/xt/fv6FDFmkQT01P4+NFm7n19A5ceEIbur/0Erm5uaxbt47HHnuMESNGMHfuXIKC/pp4tXv37ixcuJD8/Hy+/vprxowZQ0REBKNGjbIxExGxQ++EZnx72yDenL2Bp3c8g6kup09oLl9+8Cpjxozh1VdfBdgzI/z111/P3XffDcCQIUNYvXo1jz/+OGeddZZtOYhnqUj3ctHR0Qe8ZjQvL2/PUbraKisrGTduHPfccw8ul4v8/HwKCwsBKC4uZv2WHTw5I5Nvf19L7s//46nnXuDOsaMBGDx4MAEBAYwZM4YxY8bQvHlzzyUnIh63+3ui9vfI7qP7B/seqa1dTAifjh7AuetGMPv1X7jptW95/oZhhAT81Y28+uqrPP3000yePJl+/frVUwYinnc0fe3bczJ47af1/K1fO24b2gGADh3cP/v168egQYNISkriww8/5Nprr92zXmhoKH369AFg6NChFBQUcM8996hIF2nEjqWPDfRzMua0DpzXszUTpq7ix9XbiR9h8dprT3DHHXeQkpKyZ/0hQ4bss+5pp512yFFt0vhpuLuX69y5837Xw2VlZVFSUrLf9S+7FRcXs3nzZu644w6io6OJjo6mR48egHsii67de/Dz2h1c0TUIV3UVg/ufuM/6vXr1oqqqiszMTM8kJSINJjk5GX9///2+R9asWYPD4Tiq68YjgvwZPSQZgG9WbGPky3NZs819APDTTz9l7NixPPXUU1x22WX1l4BIAzjSvvbD+Zt4dNoqhndryaPnHbfnutG9JSQk0KxZMzZs2HDIf/OEE04gKyuLqqr9J48SkcahPvrYts1CeOvqPky8pg/hrd3LX/fSN/yxKY8uXboAf51R382yrAOOZhPfoU/Xy40YMYLp06eza9euPc9NmTKF4OBgTjnllAOuExYWxqxZs/Z5PP7SmwBEDb6KQTc+wne3Deb2CwYCsHjx4n3WX7RoEQCJiYkeyEhEGlJgYCBDhgzh448/3uf5KVOmcNJJJxEZGXlU7/f5Z58RGxvL5NvPpaC0kvNensu9L03miiuuYOzYsdx11131Gb5IgziSvvbduRnc//lyhnSK4/nLe+J07F+gA6SlpbFz506SkpIO+W/OnTuXNm3a7Hd3BBFpPOqzjz2tcwuuTi4DINcZxQWv/spHmcFERkUzc+bMfZadMWPGnhNw4pucqamph3r9kC+K5x133HG88cYbzJo1i/j4eH788Ufuu+8+br/9dkaMGLFnuZSUFJYuXcp5552Hw+EgMTGRxMREnBHNeWtJEZN+386u379i7C238t4D19IsNICwsDCWLl3Km2++SXBwMKWlpXzxxRc8+OCDjBw5kmuuuca+xEWk3iQlJZGamkpubi5BQUG8+uqrvPfee7z55pskJ7vPjGdmZhIbG0tiYuKejv+iiy5i48aN7Nq1i7S0NB599FEmT57ME088wSVnncYFvVrz66KlvPPQ9YTEtuafY+6gOC+bzZs3s3nzZsrLy2nWrJmdqfuaCce4fmp9BOGLDtfX/nf2eh79ejV579xEl4A8Lr7wAgDuuusuZs6cya5du9i2bRvffvstN998M7GxsbzyyisEBASQmZnJRRddREVFBQUFBSxfvpwnnniC999/nyeffHLPEHgRaZzq2scOHz6cbdu2UVBQQEZGBpMmTeKRCRO48MIL+erl8fg7HXz2x5+UVRtmTn6NwvJqHK4qnn76aT777DPeffddEhIS7Exdjs2h+/TD3EhdvMDKlSutIUOGWEFBQVbLli2tBx980KqqqtpnmYSEBOvqq6/e8/emncXWvZ8us1Lu/9rq9OA31gPv/WgB1tSpU/dZr6CgwLrzzjut9u3bW0FBQVZycrI1btw4q7CwsCFSE5EG8vnnn1vdunWzAgICrE6dOlmTJ0/e5/WMjAwLsN555509z913331Wx44dreDgYCsoKMg64YQTrEmTJu2z3sSJEy3ggI+9v5OkXhyuz1affgwO1NdWVlZaz36fZiXcM80a/cGi/frayZMnWwMGDLCio6Ot4OBgq1OnTtYdd9xh7dixY88y+fn51pVXXmklJiZagYGBVosWLawhQ4ZYX3/9tQ1Ziogn1KWPffDBB61u3bpZoaGhVmRkpNWrVy/rxRdftCoqKvYss72g1Jrw1Uorbuj1ll9kC8vh9LM6dO5qffrppw2VmnjOIftsY9W6xqF2DV8fhwmk4WzYUcQbP2/g08WbcRjDZSe25V9DUmgZGXT4lUVE6ii/pIL//LCWD+dvwt/p4PpBSdw4uD3hQZrhvR4deHz1kVOffhTKKqu56+OlTFv2Jxf3bsMTF3bHz6mrBEWk4eUUlTNxTgaT5mVSVF7F6Z2bc8Pg9vRLanbAuTGkUTjkB6ci3QdYlsWv63cycU4GM9ZkE+Dn4G9923HTKckqzkWkQW3MKebZH9YydelWmoUGcP2gJEad2I7o0AC7Q/MFKtIbyPbCMm6c9DvLthRwz/DO/HNwe+0Ii4jtCkoqeW/eRt6Zm0FeSSXd4iO47uQkzjk+ngA/HURsZFSk+6qyymq+WrKViXMzWLNtFzGhAVzZP4Er+ycQFx5od3gi0oQt31zAU9PX8Mu6HAJe14cxAAAVPklEQVT9HFzQqzVXD0ikS6sIu0NrzFSkN4D5G3Zy60dLKCyr5PnLenJmt5Z2hyQiso/Simo+/2MLE+dmkJ5dRFx4IFf1T+CK/gk000HxxkJFuq/Zsauc93/L5IP5meQUVdC5ZTjXnpzEyB7xBPk77Q5PRGSPtG27ePfXjXz+x2bKKl10i4/ggl6tOef4eI30OXoq0j2oqtrFCzPW8cqsdBJiQnnlbyfQNV4HlUTEe1mWxex1Obw9J4PZa3cQ6Ofg7O6t+Fu/dvROiNYIIO+mIt1XrNpayMS5GXy1ZCsV1S5O79yc605O4qTkGP0nFBGvlldcwed/bOGLJVtYtrkAgC6tIhjSKY5TOsbRo22UDjIenop0D1m+uYAHvljOss0FXNy7DRNGdiM0ULdGE5HGY932Xbw3byNf/LGVovIqOrUI52/92nF+r9ZEBmt+GC+kIr0xc7ksZq7JZuLcDH5dv5NgfyeX9GnDNQMSaR8XZnd4IiJHLT27iBmrtzNzTTa/Z+ZR7bLwdxq6toqgV7toerWLokurCBJjQnWN3b5UpNezwrJKnp2exv9+yyQmLJDx53blnOPj7Q5LRKTOisurmLp0Kx8u2MSyzQUE+Ts49/h4/tavHT3bRunEnvdQkd4YFZdX8enizbwzdyMZOcW0igzi6gGJjDqxHZEhOhomIr6hoLSS+Rt28kdWPosz81i2uYDSymoAnA5DQkwIHZqH0aF5OEmxobSLCaFdsxDiwgJxOJrcjoaK9HpSUeViyu9ZvDhjHTlF5VzVP4E7h3UiQncjEBEfsmJLAR/M38SXS7ZQUlFNl1YR/K1vW87tEU9UiK5dt5mK9MZka34p783byOT5mygsq6Jn2yiuOzmJ4ce1xF+3fhERH1dV7SJt+y7WbS8iPbuIddm7WJddRObOEqpdf3VJgX4O2jVzF+xtm4WQUFO87/7bR4fOq0g/RtUuiy/+2MJzP65lc14pJyZG89A5XTm+TZTdoYmIeExReRVfLtnCh/M3sXJrIQFOB6d1bs5Fvdtwaqc41Rj2UJHeGCzJyuftORl8s/xPLMtixHGtuPbkJHonRNsdmoiI7SqqXGzJLyVzZzFZuSVsyi0hc6f756bcEkoqqvdZvkVE4F8FfLNQ2sUE1xTxocSGBTTW4X4q0usot7iCKQuzeP+3TLbkl3Jc6wjuOrMTp3SMa6zbgohInazcWsBni7fw5ZIt5BRV0Cw0gJE94rm4dxu6xUfoO7HhqEj3VpXVLn5YtZ2352SwKDOP8EA/Lu/blqsHJNImOsTu8EREGgXLssgtriAzt8RdwO8sIbOmeM/KLeHPgrJ9lg/2d7oL9r3Ovu/+u010MIF+XnsWXkX6UbAsiyVZ+bz/2yamLttKRZWLk9rHcM3ARM7s2kI7oiLSpFVWu5i9dgefLd7CD6u2U1HtolOLcEb2jOfMri1IaR6m70nPUpHubbJyS5iyMIv/+z2L7F3ltGsWwj8GJnJJn7aEaTZZEZF6VVZZzea80v3OwO/+e/c18ADGQMuIoD2Fe3xUMNEh/kSHBhAd4n5EBvsT5O8g0M9JoL+DQD9HQ+3IqEg/Aut3FPHlkq18tWQLG3eWEBLg5KIT2vD3kxLo2CLc7vBERLxOQUkl05Zv5bPFW1iUmQdA+9hQzuzWkjO6NqdHmyj8NCS+vqlI9wbF5VXMWJPNp4s2M3vdDgwwpFNzRvVtx5DOzXE2vQmQRERsZ1kWO4rKD1rAby8sP6L3CfBzF+sBTgctIoL45tZBnghXRfoBuFwWy7cUMCstmx9Xb2fFlkKMgQHJMZzXozUjurckXBPCiYgckW0FZfywejvfr9zGvPU7qXJZhAf60a99DCenxDAwJZbkuLCmOHlrfVORbpeyympmrclm2rI/mbFmO2WVLlpFBnFpn7ZcdmJb4qOC7Q6xTlJTU5kwYYLdYYhILePHjyc1NdXuMHxKVbWLgtJK8koqyCupJK+4goLSSsqrXDWPasor//q9ospFZLA/dw/v7IlwVKTjPrCyfkcxCzJyWZCxkznpOeQUVWAM9GobxdnHx3Pu8a1oHhFkd6jHRH2tiHdpin1sQUklc9JzmJOew9z0HDbllgAQEeRHj7ZR9GobRc92UXSLj6R5eKCGxx+dQzaWxlbXo2qXxZpthcxZ596YF2TkUl7lIjYsgEt6t+Wc41txYmIzHXkSEWkk/JwOYsICiQkLtDuUJsmyLHbsKmfl1kJWbi1gxZZCfs/MJaeoAoDYsEAGJMcypHMcp3RsTrNQ3VJIRKS+RIb4c/bxrTj7+FaA+5Ldeet38kdWHn9syuflWensvvFKZLA/HVuE0bFFOB1bhJMcF0ZCTAitIoM0VL4OVKTXUVF5FRk7itmQU8SqPwtZsimf5VsK9sww3LFFGFf0S+D0Ls3pl9RMG6eIiEgtxeVV5BSVk1NUUfOznK35pWzcWULmzmIyc0rYVV61Z/mEmBAGdYijb1Iz+iU1Iyk2VGduREQaSNuau6ZcemJbAEoqqli+uYA123bV3D51F1OXbqWw7K/vbT+HoU10MO1iQklo5p6gtUVEUM0jkJaRQYQEqCStrUkPd7csi4pqF6UV1ZRUVFNaWV3r9yp2lVWRU1TBjl3lZO8qI3tXOZk7i/e5TjHA6aBLfAQ920TSo20UA5JjaRnZuIfZiYiI1/G64e5Tl25l4cZcql0WLst9fXi1ZeGyrJrfobzS3a8WV1RRUl7zs6KaovIqKqpc+72nn8PQtmbivsSYEBJjQ+naKoIu8RFE6NpyERGvZlkW2wvL2ZBTxKaaeV4yd995ZWfxPgX8buGBfrSIdBft0SEBRAT7ExHkT0SwX81Pf8KD/AjycxJQM/9LgJ/7YXB3bpZl1ZzVt6istiirdNdz5ZUuyiqrKauqprTir9/LKve+ZO2vS9fKKqspr3JRUeXCZVkkxYby9CU9PNFU3jXc/d25Gbw/fxNOY3A4DE4Hf/2+10+nY/fv7HeUvHZGLsuiymVRWe2iqtr9e5XL/Xtltcv9d83v1a6/niuvcv99JEIDnDSPCCIuLJBBHeJIig0lOS6UpNgwEmNDvPmWPSIiIh6xNCufr5ZuxWkMxuzbpztq+vJAPwchAU7CAv1oHh5IaKAfoQF+hAQ4iQoJIDYsgNjwQOLCAokNCyQ2LECjz0REGiljDC0jg2gZGcSA5P1fLyqvYltBGdmFZWwrLGN7YTnbC8vYXvP3n/mFFJZVUlBaSWW1Z88XB/g5CPJzEOjvJNDPQVDNz8CaAwD+DvdPOzT4mfRpy7by7fJtVO8+2l7zs9r118O1+++ao/L7BrR/SAaDn9Pg73Dg53TvFPg7Hfjt/uk0+Dkc+Nd6LdDfQUiAH8H+ToIDnIQEOAnyd//c/XtYoB+xYe6dChERERt53Zl0ERERT7As9wnVwrJKCkurKCyrpLzSRUW1+yx3RZWLiupqLMt9+1RHzUldYwz+DuMuuP0dBPu7a7ogf2fN7+5iPMDpsHueMM3uLiIi4gNUpIuIiPiGQ/bpGk8mIiIiIiIi4iVUpIuIiIiIiIh4CRXpIiIiIiIiIl5CRbqIiIiIiIiIl1CRLiIiIiIiIuIlVKSLiIiIiIiIeAkV6SIiIiIiIiJeQkW6iIiIiIiIiJdQkS4iIiIiIiLiJVSki4iIiIiIiHgJFekiIiIiIiIiXkJFuoiIiIiIiIiXMJZlHfTFCRMmfAfENlw49S4e2Gp3EB7m6zn6en6gHH2FcvQN3pxjzvjx44fXdWUf6NPt4M3bg7dT29Wd2q7u1HZ1o3aru7q23aH7dMuyfPaRmppq2R2DclR+ylE5+spDOerR1B7aHtR2arvG9VDbqd18pe003F1ERERERETES/h6kT7B7gAagK/n6Ov5gXL0FcrRNzSFHOXIaXuoO7Vd3ant6k5tVzdqt7rzSNsd8pp0EREREREREWk4vn4mXURERERERKTRUJEuIiIiIiIi4iVUpIuIiIiIiIh4iUZTpBtjRhtjMowxZcaYRcaYQUe43snGmCpjzIoDvBZhjHnRGLPVGFNujEk3xlxa/9EfmfrO0RjzkzHGOsBjpWcyOKJYPfE53mqMWWOMKTXGbDbGvGKMCav/6I+MBz5Hf2PMw8aY9TXvudQYU+d7JdeHo8nRGHPqQbbDzrWWu8gYs6rm/+IqY8wFns/koDHXa37GmG7GmE+MMRtqXkttkEQOwQM53mCM+cUYk2eMyTfGzDLGnNww2Rw07vrO8RJjzO81+RUbY5YYY65umGykPhzlNnGKMeZXY8zOmv5ljTHmrlrLNJltor7brtbyo2r+v03zTPT28cA2d81BvquCPJ9Nw/LENme8bN/fUzyw3XldTeEpHtrujr5WsfveckfyAC4DKoEbgC7AS0AR0O4w60UDG4DpwIpar/kD84FvgZOBxJqfJ/pQjs2Alns9EoBCYLwP5fg3oBz4e81neBqQAbztQzk+CfwJnA20B24GSoFejSFH4FTAArrW2h6dey1zElAFPFDzng/U/N3PR/I7EXimZnvdAKTa8dl5OMcPgDFAL6AT8DpQDHTwoRxPA84HOgPJwK012+lZdn6eenhsm+gNXA50A5KAK2u26dFNbZvwRNvttWx7YDMwG5hmd67e3m7ANTXP7f091dLuXBtJ23nVvn8jazuvqikaWdvVqVaxvTGOsMHmA2/Wem4d8Phh1vsMGA+ksn/hcyPuneUAu/PzVI4HWPYK3DsPbX0lR+Bl4Odaz004XFs0shy3ArfWeu5T4P3GkCN/FT+xh3jPKcAPtZ77EZjsC/nVWn4F9hfpHs2xZh0DbAPG+mqONestPtz/bz2841HX7+day392uO8lX9wmPNV2/FU0XQ28i+8V6fXebriL9CK7c2ukbedV+/6Nqe0O8LqtNUVjajvqWKt4/XB3Y0wA7qMU39d66XtgwCHWGw20AB47yCLnA3OBl4wx24x7eG2qMca/HsI+Kh7MsbYbgO8sy8qqS5zHwoM5zgF6GmP61yzfDhgJfHOsMR8tD+YYCJTVeq4U99HfBlXXHGv8boz50xgzwxgzpNZrJx3gPacfwXvWKw/m5zUaMMcAIAjIq1ukddcQORq303GPGph9TAGLxx3jNrH7PXrVLPvzQV73yW3Cw233/4CNlmW9d6xxehsPt1uwMSazZtjstJrlfIYH285r9v09pSG+62rYVlN4igfbrk61itcX6UAs4AS213p+O+7hFvsxxnTHfVbySsuyqg/yvu2BS3AfxT0beAi4CXi8HmI+Wp7Kce/lOwKnAG8eW6h15pEcLcv6CLgfmG2MqQQygeXAPfUU99Hw1Oc4HbjNGNPJGOMwxpwBXAi0qp+wj8pR54h7qP7NwEW4404DZtS6xqflUb6np3gqP2/SUDk+hnuI2FfHFG3deCxHY0ykMaYIqAC+Bm6xLOvbeoxdPKMu2wQANYVQOfA78KplWa/Xet3XtwmPtJ0x5kzgUuCf9Ruu1/DUNpcGXAucB4zCfRB/rjGmQ30F7gU81XbetO/vKR77rttrObtrCk/xSNvVtVbxO+rwvZwxJhD30Nm7LMvKOMSiDiAbuKGmOFpkjIkBnjPGjLNqxiJ4o6PIcW834N4J/dpjgdWjI83RGHMK7i/Z0biHqKQAL+AeRvJwA4RaZ0fxOd6K+4twFe7huOuBd3B30l7Psqw03DsVu80zxiQC44Bf7IipPvl6fnD0ORpjbsW94z3UsqzChojxWB1FjruAnkAYcDrwH2PMRsuyZjRQqNLwBuH+vPsDTxpjMizL+t9er2ubOLgDtp0xJg738PZRlmXl2xmglzroNmdZ1jxg3u4FjTG/AkuAscAtNsTqbQ71/7XR7vs3kMN91+3WqGqKBnLQtqtrrdIYivQcoBr3cOC9tcB9vWNtrXBf6P+OMeadmuccuEei7Z7M5XvcG1dlrbOXq4EQ3EdSdtRfCoflqRyBPcM3rsZ9jUVVfQd/hDyV42O4r/t4q2aZ5caYUOAtY8wjDZyvR3K0LGsHcL5xz9wag/sa9SdwX1fV0I42x4OZj3uijd221cN71gdP5edNPJqjMeY24FFghGVZC+oa5DHyWI6WZbmA9Jo/lxhjuuA+Qq6CzLvVeZvY6yDqcmNMC9xzh/xvr9d9fZvwRNt1w90HzjDG7F7cAVDT/3WrOXDWmHlsm6u1bLUx5nfAl86ke6rtvGnf31M8ut15SU3hKZ5quzrVKl4/3N2yrApgEXBGrZfOAH49wCpbgO64j2rvfryOuwPtudc6c4EUY8zebdARKMH9ITUYD+a42/m4v3zerr+oj44HcwzB/R9qb9W4J61qUJ7+HC3LKrMsawvug2sXAV/WZ/xHog45HkxP3J3lbvPq4T2PmQfz8xqezNEYcwfuAv1sy7LmHEucx6KBP0cH7nkjxIvV4zZxJJ+3T20THmq7hezf/32Fe9RKT9wzHzdqDbXNGfdRjuPx0j6nLjzYdl6z7+8pDbDd2V5TeIoH265utUpdZ79ryAfu6fArgOtxn3l8Afe1jgk1r08CJh1i/VT2nzG7Le5bB7yEe5KXYbhvAfK0r+S412s/UmvmbF/Jsea5Qtxnu5Jw/0dKBz71oRz74b5Gtj3u4TQzcJ9Fj2oMOQK34f5S74D77MnjuIftX7jXMgNwzxJ6L+5bGd2H+xYYdt2Crb7zC+CvHdF03AdjegIpPvQZjqt5z0vZ9zYtkT6U4wPA0Jr/i12AO2u205vsyFEPj28TY4FzaraJDsB1uPubJ5raNuGJtjvAv/Euvje7uye2ufG491nb4+5HJtZsc33tzrcRtJ1X7fs3prbba1mvqCkaU9tRx1rF9sY4ikYbDWzEfZ+5RcDgvV77CfjpEOumcoACFvd1A7/inik7A3gEG2/L4KEc2wMu4FK7P0NP5Ij7rPJ43LdHKAWygFeBaB/K8RTc16OX4T7SOwmIbyyfI3D3Xp9PLu4zJfvdQxi4GFiD+8txNXsVR409P9z3xbQO8DjottAIc9x4kBzf9aEcH6+1zK+4r6m1JT89PL5N3AasxH3P2wLct1YbDTia4jZR3213gPd/Fx8r0j20zT2He+KpctzXV08HTrI7z8bQdjXLedW+fyNrO6+qKRpL21HHWsXUrCwiIiIiIiIiNvP6a9JFREREREREmgoV6SIiIiIiIiJeQkW6iIiIiIiIiJdQkS4iIiIiIiLiJVSki4iIiIiIiHgJFekiIiIiIiIiXkJFuoiIiIiIiIiXUJEuIiIiIiIi4iVUpIuIiIiIiIh4if8Pn73/Tauj7FgAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 993.6x331.2 with 2 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IPQ64pVs9v-C",
        "outputId": "84c98081-32b4-448c-8aea-079eed366d61",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "# Estimates of parameters mu & std from the posterior samples of the fitted model\n",
        "\n",
        "trace['mu'].mean(), trace['std'].mean()"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(0.49726704003022826, 0.3482576656808353)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 12
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Og5yO1SJ-usd"
      },
      "source": [
        "### References \n",
        "* https://www.youtube.com/watch?v=j5QyTDUww6U&t=1520s"
      ]
    }
  ]
}