{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6a0dc113",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc as pm\n",
    "import numpy as np\n",
    "import arviz as az\n",
    "import pandas as pd\n",
    "import pytensor.tensor.subtensor as st\n",
    "from itertools import combinations\n",
    "\n",
    "%load_ext lab_black"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8a18872",
   "metadata": {},
   "source": [
    "# 9. Simvastatin*\n",
    "\n",
    "This one is about factorial designs (2-way ANOVA) with sum-to-zero and corner constraints.\n",
    "\n",
    "Adapted from [Unit 7: simvastatin.odc](https://raw.githubusercontent.com/areding/6420-pymc/main/original_examples/Codes4Unit7/simvastatin.odc).\n",
    "\n",
    "Data can be found [here](https://raw.githubusercontent.com/areding/6420-pymc/main/data/simvastatin_data.tsv).\n",
    "\n",
    "Thanks to [Anthony Miyaguchi](https://github.com/acmiyaguchi) for updating this example!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc5f944e",
   "metadata": {},
   "source": [
    "In a quantitative physiology lab II at Georgia Tech, students were asked to find a therapeutic model to test on MC3T3-E1 cell line to enhance osteoblastic growth. The students found a drug called Simvastatin, a cholesterol lowering drug to test on these cells. Using a control and three different concentrations, $10^{-9}$, $10^{-8}$ and $10^{-7}$ M, cells were treated with the drug. These cells were plated on four, 24 well plates with each well plate having a different treatment. To test for osteoblastic differentiation an assay, pNPP, was used to test for alkaline phosphatase activity. The higher the alkaline phosphatase activity the better the cells are differentiating, and become more bone like. This assay was performed 6 times total within 11 days. Each time the assay was performed, four wells from each plate were used.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fbbed4a",
   "metadata": {},
   "source": [
    "## Notes: \n",
    "\n",
    "A [good explanation](https://stats.stackexchange.com/questions/257778/sum-to-zero-constraint-in-one-way-anova) of STZ constraints.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "26142563",
   "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>apa</th>\n",
       "      <th>conc</th>\n",
       "      <th>time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.062</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.517</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.261</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     apa  conc  time\n",
       "0  0.062     1     1\n",
       "1  0.517     1     1\n",
       "2  0.261     1     1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv(\"../data/simvastatin_data.tsv\", sep=\"\\t\")\n",
    "data.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a091a38f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 0, 0, 0, 0, 1, 1,\n",
       "        1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2,\n",
       "        3, 3, 3, 3, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 0, 0,\n",
       "        0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 0, 0, 0, 0, 1, 1, 1, 1,\n",
       "        2, 2, 2, 2, 3, 3, 3, 3]),\n",
       " array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1,\n",
       "        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,\n",
       "        2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4,\n",
       "        4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5,\n",
       "        5, 5, 5, 5, 5, 5, 5, 5]),\n",
       " {'conc': Int64Index([1, 2, 3, 4], dtype='int64'),\n",
       "  'time': Int64Index([1, 2, 3, 4, 5, 6], dtype='int64'),\n",
       "  'id': RangeIndex(start=0, stop=96, step=1)})"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set up alternate coordinates, the ID3 or clusters column\n",
    "conc_idx, conc = pd.factorize(data[\"conc\"])\n",
    "time_idx, time = pd.factorize(data[\"time\"])\n",
    "coords = {\"conc\": conc, \"time\": time, \"id\": data.index}\n",
    "\n",
    "conc_idx, time_idx, coords"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de2cd6c1",
   "metadata": {},
   "source": [
    "## Model 1 with sum-to-zero constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "891c822f",
   "metadata": {
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [mu0, _alpha, _beta, _alphabeta, tau]\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:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\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='12000' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [12000/12000 00:07&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 8 seconds.\n"
     ]
    }
   ],
   "source": [
    "def differences(var, index):\n",
    "    \"\"\"Calculate differences between levels with names like \"alpha[low] - alpha[high]\".\n",
    "\n",
    "    var: aesara.tensor.var.TensorVariable\n",
    "    index: pandas.Index\n",
    "    \"\"\"\n",
    "    name = var.name\n",
    "    for i, j in combinations(range(index.size), 2):\n",
    "        a, b = index[i], index[j]\n",
    "        pm.Deterministic(f\"{name}[{a}] - {name}[{b}]\", var[i] - var[j])\n",
    "\n",
    "\n",
    "with pm.Model(coords=coords) as m:\n",
    "    apa_data = pm.Data(\"apa_data\", data.apa, mutable=False)\n",
    "    time_idx_data = pm.Data(\"time_idx_data\", time_idx, dims=\"id\", mutable=False)\n",
    "    conc_idx_data = pm.Data(\"conc_idx_data\", conc_idx, dims=\"id\", mutable=False)\n",
    "\n",
    "    mu0 = pm.Normal(\"mu0\", 0, tau=0.0001)\n",
    "    _alpha = pm.Normal(\"_alpha\", 0, tau=0.0001, dims=\"conc\")\n",
    "    _beta = pm.Normal(\"_beta\", 0, tau=0.0001, dims=\"time\")\n",
    "    _alphabeta = pm.Normal(\"_alphabeta\", 0, tau=0.0001, dims=(\"conc\", \"time\"))\n",
    "    tau = pm.Gamma(\"tau\", 0.001, 0.001)\n",
    "    sigma = pm.Deterministic(\"sigma\", 1 / tau**0.5)\n",
    "\n",
    "    # sum-to-zero constraints\n",
    "    # sets the first element of a dimension to the negative sum of the rest\n",
    "    sst_1d_0 = lambda var: st.set_subtensor(var[0], -var[1:].sum(axis=0))\n",
    "    sst_2d_0 = lambda var: st.set_subtensor(var[0, :], -var[1:, :].sum(axis=0))\n",
    "    sst_2d_1 = lambda var: st.set_subtensor(var[:, 0], -var[:, 1:].sum(axis=1))\n",
    "\n",
    "    alpha = pm.Deterministic(\"alpha\", sst_1d_0(_alpha), dims=\"conc\")\n",
    "    beta = pm.Deterministic(\"beta\", sst_1d_0(_beta), dims=\"time\")\n",
    "    _alphabeta = sst_2d_1(_alphabeta)\n",
    "    alphabeta = pm.Deterministic(\n",
    "        \"alphabeta\", sst_2d_0(_alphabeta), dims=(\"conc\", \"time\")\n",
    "    )\n",
    "\n",
    "    mu = (\n",
    "        mu0\n",
    "        + alpha[conc_idx_data]\n",
    "        + beta[time_idx_data]\n",
    "        + alphabeta[conc_idx_data, time_idx_data]\n",
    "    )\n",
    "    pm.Normal(\"apa\", mu, tau=tau, observed=apa_data, dims=\"id\")\n",
    "\n",
    "    # calculate differences between levels with appropriate names\n",
    "    differences(alpha, coords[\"conc\"])\n",
    "    differences(beta, coords[\"time\"])\n",
    "\n",
    "    trace = pm.sample(2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ef27e6a9",
   "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>mu0</th>\n",
       "      <td>0.239</td>\n",
       "      <td>0.025</td>\n",
       "      <td>0.192</td>\n",
       "      <td>0.286</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tau</th>\n",
       "      <td>17.233</td>\n",
       "      <td>2.874</td>\n",
       "      <td>12.020</td>\n",
       "      <td>22.775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>0.243</td>\n",
       "      <td>0.021</td>\n",
       "      <td>0.206</td>\n",
       "      <td>0.283</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[1]</th>\n",
       "      <td>0.049</td>\n",
       "      <td>0.043</td>\n",
       "      <td>-0.033</td>\n",
       "      <td>0.130</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[2]</th>\n",
       "      <td>0.068</td>\n",
       "      <td>0.044</td>\n",
       "      <td>-0.014</td>\n",
       "      <td>0.149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[3]</th>\n",
       "      <td>-0.075</td>\n",
       "      <td>0.043</td>\n",
       "      <td>-0.153</td>\n",
       "      <td>0.008</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[4]</th>\n",
       "      <td>-0.042</td>\n",
       "      <td>0.042</td>\n",
       "      <td>-0.121</td>\n",
       "      <td>0.037</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1]</th>\n",
       "      <td>0.046</td>\n",
       "      <td>0.055</td>\n",
       "      <td>-0.058</td>\n",
       "      <td>0.149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2]</th>\n",
       "      <td>-0.150</td>\n",
       "      <td>0.055</td>\n",
       "      <td>-0.258</td>\n",
       "      <td>-0.047</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3]</th>\n",
       "      <td>-0.019</td>\n",
       "      <td>0.055</td>\n",
       "      <td>-0.118</td>\n",
       "      <td>0.090</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[4]</th>\n",
       "      <td>0.230</td>\n",
       "      <td>0.055</td>\n",
       "      <td>0.123</td>\n",
       "      <td>0.329</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[5]</th>\n",
       "      <td>-0.022</td>\n",
       "      <td>0.056</td>\n",
       "      <td>-0.129</td>\n",
       "      <td>0.081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[6]</th>\n",
       "      <td>-0.086</td>\n",
       "      <td>0.057</td>\n",
       "      <td>-0.190</td>\n",
       "      <td>0.023</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 1]</th>\n",
       "      <td>-0.085</td>\n",
       "      <td>0.097</td>\n",
       "      <td>-0.267</td>\n",
       "      <td>0.099</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 2]</th>\n",
       "      <td>-0.083</td>\n",
       "      <td>0.096</td>\n",
       "      <td>-0.265</td>\n",
       "      <td>0.095</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 3]</th>\n",
       "      <td>-0.210</td>\n",
       "      <td>0.098</td>\n",
       "      <td>-0.393</td>\n",
       "      <td>-0.028</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 4]</th>\n",
       "      <td>0.501</td>\n",
       "      <td>0.095</td>\n",
       "      <td>0.316</td>\n",
       "      <td>0.673</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 5]</th>\n",
       "      <td>-0.019</td>\n",
       "      <td>0.097</td>\n",
       "      <td>-0.207</td>\n",
       "      <td>0.160</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 6]</th>\n",
       "      <td>-0.104</td>\n",
       "      <td>0.096</td>\n",
       "      <td>-0.276</td>\n",
       "      <td>0.081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 1]</th>\n",
       "      <td>0.150</td>\n",
       "      <td>0.097</td>\n",
       "      <td>-0.027</td>\n",
       "      <td>0.338</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 2]</th>\n",
       "      <td>-0.100</td>\n",
       "      <td>0.097</td>\n",
       "      <td>-0.283</td>\n",
       "      <td>0.081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 3]</th>\n",
       "      <td>0.354</td>\n",
       "      <td>0.097</td>\n",
       "      <td>0.173</td>\n",
       "      <td>0.535</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 4]</th>\n",
       "      <td>-0.176</td>\n",
       "      <td>0.095</td>\n",
       "      <td>-0.362</td>\n",
       "      <td>-0.003</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 5]</th>\n",
       "      <td>-0.187</td>\n",
       "      <td>0.094</td>\n",
       "      <td>-0.358</td>\n",
       "      <td>-0.006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 6]</th>\n",
       "      <td>-0.042</td>\n",
       "      <td>0.098</td>\n",
       "      <td>-0.225</td>\n",
       "      <td>0.141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 1]</th>\n",
       "      <td>-0.043</td>\n",
       "      <td>0.098</td>\n",
       "      <td>-0.226</td>\n",
       "      <td>0.140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 2]</th>\n",
       "      <td>0.056</td>\n",
       "      <td>0.097</td>\n",
       "      <td>-0.119</td>\n",
       "      <td>0.241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 3]</th>\n",
       "      <td>-0.046</td>\n",
       "      <td>0.097</td>\n",
       "      <td>-0.228</td>\n",
       "      <td>0.133</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 4]</th>\n",
       "      <td>-0.138</td>\n",
       "      <td>0.096</td>\n",
       "      <td>-0.310</td>\n",
       "      <td>0.051</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 5]</th>\n",
       "      <td>0.158</td>\n",
       "      <td>0.094</td>\n",
       "      <td>-0.020</td>\n",
       "      <td>0.335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 6]</th>\n",
       "      <td>0.013</td>\n",
       "      <td>0.098</td>\n",
       "      <td>-0.171</td>\n",
       "      <td>0.198</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 1]</th>\n",
       "      <td>-0.022</td>\n",
       "      <td>0.095</td>\n",
       "      <td>-0.196</td>\n",
       "      <td>0.158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 2]</th>\n",
       "      <td>0.126</td>\n",
       "      <td>0.094</td>\n",
       "      <td>-0.056</td>\n",
       "      <td>0.301</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 3]</th>\n",
       "      <td>-0.098</td>\n",
       "      <td>0.095</td>\n",
       "      <td>-0.283</td>\n",
       "      <td>0.076</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 4]</th>\n",
       "      <td>-0.187</td>\n",
       "      <td>0.097</td>\n",
       "      <td>-0.353</td>\n",
       "      <td>0.012</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 5]</th>\n",
       "      <td>0.048</td>\n",
       "      <td>0.095</td>\n",
       "      <td>-0.132</td>\n",
       "      <td>0.223</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 6]</th>\n",
       "      <td>0.133</td>\n",
       "      <td>0.095</td>\n",
       "      <td>-0.047</td>\n",
       "      <td>0.311</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[1] - alpha[2]</th>\n",
       "      <td>-0.019</td>\n",
       "      <td>0.071</td>\n",
       "      <td>-0.156</td>\n",
       "      <td>0.112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[1] - alpha[3]</th>\n",
       "      <td>0.125</td>\n",
       "      <td>0.069</td>\n",
       "      <td>-0.004</td>\n",
       "      <td>0.258</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[1] - alpha[4]</th>\n",
       "      <td>0.092</td>\n",
       "      <td>0.069</td>\n",
       "      <td>-0.034</td>\n",
       "      <td>0.226</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[2] - alpha[3]</th>\n",
       "      <td>0.144</td>\n",
       "      <td>0.071</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[2] - alpha[4]</th>\n",
       "      <td>0.111</td>\n",
       "      <td>0.070</td>\n",
       "      <td>-0.021</td>\n",
       "      <td>0.241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[3] - alpha[4]</th>\n",
       "      <td>-0.033</td>\n",
       "      <td>0.069</td>\n",
       "      <td>-0.164</td>\n",
       "      <td>0.098</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[2]</th>\n",
       "      <td>0.195</td>\n",
       "      <td>0.086</td>\n",
       "      <td>0.041</td>\n",
       "      <td>0.365</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[3]</th>\n",
       "      <td>0.064</td>\n",
       "      <td>0.086</td>\n",
       "      <td>-0.100</td>\n",
       "      <td>0.224</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[4]</th>\n",
       "      <td>-0.185</td>\n",
       "      <td>0.086</td>\n",
       "      <td>-0.343</td>\n",
       "      <td>-0.019</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[5]</th>\n",
       "      <td>0.068</td>\n",
       "      <td>0.086</td>\n",
       "      <td>-0.093</td>\n",
       "      <td>0.229</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[6]</th>\n",
       "      <td>0.131</td>\n",
       "      <td>0.087</td>\n",
       "      <td>-0.028</td>\n",
       "      <td>0.297</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2] - beta[3]</th>\n",
       "      <td>-0.131</td>\n",
       "      <td>0.085</td>\n",
       "      <td>-0.292</td>\n",
       "      <td>0.027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2] - beta[4]</th>\n",
       "      <td>-0.380</td>\n",
       "      <td>0.085</td>\n",
       "      <td>-0.542</td>\n",
       "      <td>-0.225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2] - beta[5]</th>\n",
       "      <td>-0.128</td>\n",
       "      <td>0.088</td>\n",
       "      <td>-0.293</td>\n",
       "      <td>0.039</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2] - beta[6]</th>\n",
       "      <td>-0.064</td>\n",
       "      <td>0.087</td>\n",
       "      <td>-0.227</td>\n",
       "      <td>0.099</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3] - beta[4]</th>\n",
       "      <td>-0.249</td>\n",
       "      <td>0.086</td>\n",
       "      <td>-0.404</td>\n",
       "      <td>-0.088</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3] - beta[5]</th>\n",
       "      <td>0.003</td>\n",
       "      <td>0.086</td>\n",
       "      <td>-0.160</td>\n",
       "      <td>0.158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3] - beta[6]</th>\n",
       "      <td>0.067</td>\n",
       "      <td>0.088</td>\n",
       "      <td>-0.093</td>\n",
       "      <td>0.235</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[4] - beta[5]</th>\n",
       "      <td>0.252</td>\n",
       "      <td>0.086</td>\n",
       "      <td>0.081</td>\n",
       "      <td>0.406</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[4] - beta[6]</th>\n",
       "      <td>0.316</td>\n",
       "      <td>0.087</td>\n",
       "      <td>0.149</td>\n",
       "      <td>0.475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[5] - beta[6]</th>\n",
       "      <td>0.064</td>\n",
       "      <td>0.087</td>\n",
       "      <td>-0.099</td>\n",
       "      <td>0.228</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       mean     sd  hdi_3%  hdi_97%\n",
       "mu0                   0.239  0.025   0.192    0.286\n",
       "tau                  17.233  2.874  12.020   22.775\n",
       "sigma                 0.243  0.021   0.206    0.283\n",
       "alpha[1]              0.049  0.043  -0.033    0.130\n",
       "alpha[2]              0.068  0.044  -0.014    0.149\n",
       "alpha[3]             -0.075  0.043  -0.153    0.008\n",
       "alpha[4]             -0.042  0.042  -0.121    0.037\n",
       "beta[1]               0.046  0.055  -0.058    0.149\n",
       "beta[2]              -0.150  0.055  -0.258   -0.047\n",
       "beta[3]              -0.019  0.055  -0.118    0.090\n",
       "beta[4]               0.230  0.055   0.123    0.329\n",
       "beta[5]              -0.022  0.056  -0.129    0.081\n",
       "beta[6]              -0.086  0.057  -0.190    0.023\n",
       "alphabeta[1, 1]      -0.085  0.097  -0.267    0.099\n",
       "alphabeta[1, 2]      -0.083  0.096  -0.265    0.095\n",
       "alphabeta[1, 3]      -0.210  0.098  -0.393   -0.028\n",
       "alphabeta[1, 4]       0.501  0.095   0.316    0.673\n",
       "alphabeta[1, 5]      -0.019  0.097  -0.207    0.160\n",
       "alphabeta[1, 6]      -0.104  0.096  -0.276    0.081\n",
       "alphabeta[2, 1]       0.150  0.097  -0.027    0.338\n",
       "alphabeta[2, 2]      -0.100  0.097  -0.283    0.081\n",
       "alphabeta[2, 3]       0.354  0.097   0.173    0.535\n",
       "alphabeta[2, 4]      -0.176  0.095  -0.362   -0.003\n",
       "alphabeta[2, 5]      -0.187  0.094  -0.358   -0.006\n",
       "alphabeta[2, 6]      -0.042  0.098  -0.225    0.141\n",
       "alphabeta[3, 1]      -0.043  0.098  -0.226    0.140\n",
       "alphabeta[3, 2]       0.056  0.097  -0.119    0.241\n",
       "alphabeta[3, 3]      -0.046  0.097  -0.228    0.133\n",
       "alphabeta[3, 4]      -0.138  0.096  -0.310    0.051\n",
       "alphabeta[3, 5]       0.158  0.094  -0.020    0.335\n",
       "alphabeta[3, 6]       0.013  0.098  -0.171    0.198\n",
       "alphabeta[4, 1]      -0.022  0.095  -0.196    0.158\n",
       "alphabeta[4, 2]       0.126  0.094  -0.056    0.301\n",
       "alphabeta[4, 3]      -0.098  0.095  -0.283    0.076\n",
       "alphabeta[4, 4]      -0.187  0.097  -0.353    0.012\n",
       "alphabeta[4, 5]       0.048  0.095  -0.132    0.223\n",
       "alphabeta[4, 6]       0.133  0.095  -0.047    0.311\n",
       "alpha[1] - alpha[2]  -0.019  0.071  -0.156    0.112\n",
       "alpha[1] - alpha[3]   0.125  0.069  -0.004    0.258\n",
       "alpha[1] - alpha[4]   0.092  0.069  -0.034    0.226\n",
       "alpha[2] - alpha[3]   0.144  0.071   0.010    0.275\n",
       "alpha[2] - alpha[4]   0.111  0.070  -0.021    0.241\n",
       "alpha[3] - alpha[4]  -0.033  0.069  -0.164    0.098\n",
       "beta[1] - beta[2]     0.195  0.086   0.041    0.365\n",
       "beta[1] - beta[3]     0.064  0.086  -0.100    0.224\n",
       "beta[1] - beta[4]    -0.185  0.086  -0.343   -0.019\n",
       "beta[1] - beta[5]     0.068  0.086  -0.093    0.229\n",
       "beta[1] - beta[6]     0.131  0.087  -0.028    0.297\n",
       "beta[2] - beta[3]    -0.131  0.085  -0.292    0.027\n",
       "beta[2] - beta[4]    -0.380  0.085  -0.542   -0.225\n",
       "beta[2] - beta[5]    -0.128  0.088  -0.293    0.039\n",
       "beta[2] - beta[6]    -0.064  0.087  -0.227    0.099\n",
       "beta[3] - beta[4]    -0.249  0.086  -0.404   -0.088\n",
       "beta[3] - beta[5]     0.003  0.086  -0.160    0.158\n",
       "beta[3] - beta[6]     0.067  0.088  -0.093    0.235\n",
       "beta[4] - beta[5]     0.252  0.086   0.081    0.406\n",
       "beta[4] - beta[6]     0.316  0.087   0.149    0.475\n",
       "beta[5] - beta[6]     0.064  0.087  -0.099    0.228"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace, var_names=\"~_\", filter_vars=\"like\", kind=\"stats\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f9ad9ec3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([<AxesSubplot: title={'center': '94.0% HDI'}>], dtype=object)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAGdCAYAAAAVEKdkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsaklEQVR4nO3de1xVZb7H8e8yEBQREUrJUmBES9EIJSw1mfJWdtSaLqRN4mUmp0Nmt6mxEm20PDUdj5XN6XSRXpaZNamVkyimlk1qijqOaVlKevKSF1RMRcTn/GHs45Yt1w1r6/N5v168Xu1nrf2s33r2A35bt+0YY4wAAAAsU8/tAgAAANxACAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAfC7Xbt26YEHHlBCQoJCQ0MVHR2tvn37Kicnp0r9vPbaa3IcR47jaMSIEdWuZ+PGjRo8eLBiYmIUGhqqX/3qV3r44Yd14MABn+uXlJRo7NixuvTSSxUSEqKOHTvqgw8+OGv/69atU1BQkB555JFq1bdkyRLPfpYnPz/fs15+fr7XsoyMDM+y0p+GDRsqJiZGXbp0UWZmphYtWqTyvilp3LhxchxHaWlp1doP4FxDCALgV+vXr1dSUpL+67/+S9u3b1diYqIiIyOVk5Ojvn37atKkSZXqZ8+ePXr00UdrXM/ixYvVqVMnzZgxQyUlJWrfvr127dql559/Xp06ddLu3bvLvGfMmDH685//rIMHD6pt27bauHGjbr31Vn344Yc+t5GZmamLLrpIY8eOrXG9NXXRRRepa9eu6tq1qzp27KiIiAjl5eVp6tSp6tmzp6688kqtX7/e7TKBgEAIAuA3J06c0K233qrdu3crLS1N27dv16pVq7R582YtWrRI4eHhGjNmjD777LMK+3rggQd04MAB9evXr9r1FBYW6o477tDRo0c1atQo/fjjj1q9erW2bdumrl27asuWLRo+fLjXe/bu3asXXnhBrVq10ubNm/XPf/5TixYtkuM4PkPO9OnTtWzZMj377LMKDw+vdq3+csMNN2jZsmVatmyZli9frk2bNungwYP629/+pg4dOmjdunXq0qWL1qxZ43apgOsIQQD8Zt68efr2228VEhKi7OxsXXjhhZ5l1113nR5//HEZYzR+/Phy+8nNzdXbb7+te+65R507d652Pf/93/+tPXv26PLLL9d//ud/Kjg4WJIUFRWlGTNmKCgoSPPmzVNeXp7nPevXr9exY8c0dOhQNWvWTJJ07bXXqlu3blq3bp0KCws96xYWFurRRx9Vt27ddNddd1W7ztrWoEED3XLLLVqxYoWuv/56HTlyRLfffrtKSkrcLg1wFSEIgN988cUXkqSUlBS1atWqzPLf/OY3kk5dA/PTTz/57OPYsWP6wx/+oIsuukhPP/10jeopvY4nIyNDF1xwgdeyli1bqmfPnpKk999/39NeWldpACoVExMjSTp06JCnbdy4cfrpp5/04osv1qjOutKgQQO99dZbCgkJ0Xfffaf33nvP7ZIAVxGCAPhNQUGBJKlFixY+l5e2nzx5Ul999ZXPdSZMmKDvvvtOzz33nJo0aVLtWk6cOKHVq1dLkrp27epzndL2FStWeNpatmwpSfr222+91v3mm28UFBSkqKgoSacutn7xxRd1zz33KCkpqdp11rXmzZtr4MCBkk4duQNsRggC4DcRERGSpB9//NHn8tPbv/nmmzLLN27cqOeee07du3fX3XffXaNa8vPzVVxcLEmKj4/3uU5p++bNmz1tV1xxhS666CK9/vrrys3NVWFhoaZMmaK1a9fq2muvVWhoqCTpvvvuU0REhCZMmFCjOt3QrVs3STprEAVsEeR2AQDOHykpKZKkVatWafv27br00ku9lp9+m3npUaNSxhjdc889OnnypF5++eUa13J6/5GRkT7XKW0/fd2GDRvqmWee0fDhw9WrVy9Pe6NGjfT8889Lkt577z0tWrRI//M//+Ppo7i4WHv37lVUVJTq169f7boruk3eH0o/l7OdkgRsQQgC4DcDBgzQxRdfrB07dmjQoEGaNWuW51qaefPmaeLEiZ51jx496vXe119/XZ9//rkefvhhJSYm1riWY8eOef77bKEkJCTEZy3Dhg3TxRdfrGnTpmnPnj1q06aNHnjgAbVt21ZHjhzRww8/rM6dO2v48OEyxuiJJ57QlClT9PPPPyssLEyjRo3SxIkTqxVoznbqTpKKioq0atWqKvd5prCwMEnyusgbsBEhCIDfhIaG6t1339WNN96oZcuWqWXLlmrbtq0KCgq0Y8cOtWzZUklJSfrss8/UqFEjz/tKnwl0ySWXKCsry2+1lDp+/LjX61JFRUWSTl0wfKa+ffuqb9++ZdonTpyo7du3a9asWapXr54mTJigp59+WjfddJNuvfVWffDBB3rmmWcUFhamxx9/vMp1L1u27KzL8vPzFRcXV+U+z3T48GFJUuPGjWvcF3Au45ogAH7VrVs35eXladiwYWrevLnnAuORI0dq1apVntuymzdv7nnPH//4R+3fv1+TJ0/2Ckc1cfopsDNPvZ3ZfrbTZWf6/vvv9fzzzysjI0OpqakqLi7W888/r9atW2vu3LkaMmSIZs+erdatW+v555/XiRMnar4jtWDbtm2STj1YEbAZR4IA+F3r1q31+uuvl2k/ceKE1q1bJ0nq1KmTp730wX2ZmZnKzMz0ek/pUYsZM2bo448/lnTqazkqEhsbq+DgYBUXF2vLli2e03Kn27JliyQpISGhMrul+++/X6GhoZ6nXm/atEkHDhzQoEGDVK/eqf+nrFevnnr37q2XX35Z33zzjdq3b1+pvutS6dGmq666yuVKAHcRggDUmZycHB0+fFgXX3yxkpOTyyz39RUWpY4ePVrm2p3yBAUFKTk5WStWrNAXX3zh81qb0ucapaamVtjfxx9/rHnz5mnKlCmeIyilAe3MJ0WXvj7bd5O5aefOnZ6v/6jJ07iB8wGnwwDUiePHj3u+duIPf/iD18ML165dK2OMz5/Sa4RKL0Iu7wtAz3TLLbdIkrKzs8s8HXnbtm3Kzc2V9P8PcTyboqIijR49WomJibr33ns97aV3WX3//fde65e+jo6OrnStdeHo0aP67W9/q6KiIrVp06bC/QbOd4QgAH7197//3evhg5K0fft2DRw4UHl5eWrXrl21v23dl/fff1+xsbGeZ9+cbuTIkYqOjtbGjRv14IMPep4btG/fPg0aNEgnTpzQDTfc4HVqzpdnn31W33//vV566SUFBf3/AfQWLVro0ksv1UcffaR//vOfkk597cZHH32k5s2bV/o0W207evSoZs+erdTUVC1atEhhYWGaNWtWmadoA7bhdBgAv1qwYIGmTJmiyMhIxcbG6tixY9q0aZOMMWrXrp0WLFjguTXdHw4fPqwffvjB57LGjRtr5syZuummm/TCCy/onXfeUcuWLbVx40YdOXJEsbGxeuONN8rtf9u2bZo0aZLS09PVo0cPr2WO42jcuHEaPny4UlJS1LZtW3377bcqKipSVlaW5zqhuvTJJ594AmFJSYkKCgq0ZcsWTwBMSkrS9OnT/fIYAuBcRwgC4FcDBw7Uzp07tXLlSm3cuFEhISFKSUnRHXfcoX//93/3awCqjOuvv16rVq3ShAkT9Omnn2r9+vVq0aKFbr75Zj3xxBMV3hn24IMPynEc/eUvf/G5fNiwYTp27JgmT56sTZs2qVWrVnrooYc0cuTI2tidCv3000+ehyCGhoYqIiJCycnJ6ty5s26++WZdf/31rtQFBCLHVOUEOwAAwHmCa4IAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJW+QDjDFGhYWFbpcBAMA5Kzw8XI7jVLgeISjAFBYWKiIiwu0yAAA4Zx08eFCNGzeucD2eExRgOBIEAEDNVPZIECEIAABYiQujAQCAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAVqq1ELRkyRI5jqNx48bVqJ9x48bJcRwtWbLEL3Wdrf/Sn8cee6xG/b322mte/WVkZPinUAAA4FdBbhcQKIYMGaLY2Fh169bNq33t2rWaNWuWVq9erby8PO3du1c9evQ4ayhLTk5WVlaWDhw4oClTptRB5QAAoDoIQb/IyMhQWlpamfY5c+bomWeeUf369dWmTRvt3bu33H6Sk5OVnJys/Px8QhAAq2zeXaj8fUcUG9VQCc3C3S4HqBAhqAK33Xab+vfvrw4dOmjfvn2KiYlxuyQA0L7DRW6X4FFw5Lgen/0vrdi639OWGtdUE29OVGTD+i5WVlZUoxC3S0AAqVIIOn78uF555RV9/PHH+vrrr/XTTz8pIiJC3bp105NPPqkrr7yywj5iY2MlSWvWrNEf//hHffjhhzp06JA6dOigJ554Qv379z/re2fNmqVJkyZp48aNatKkiW677Tb9x3/8hxo0aODXGk/Xvn37Kq0PoHqOHD/hdgnnlE4Tct0uwaOeI4WHBmvqoGSlxEXqq60FGjN7vXpP/kwnjdvVefv6qT5ulxDwGta35/hIlfZ0//79Gj16tLp3764bb7xRkZGR2rJliz788EN98skn+uyzz5SSklJhP8ePH1fPnj119OhRDRkyRAcOHNDMmTM1cOBATZ8+XYMHDy7znqlTp+qTTz7RgAEDlJaWpvnz5+vFF1/Uvn379Pbbb/u9RgB1q93YHLdLQDWdNNLTN3dQv46njpT36xgjI6PMGWtcrqws5lnF8if1c7uEOlOlEBQZGalt27apRYsWXu0bNmxQly5dNGbMGC1cuLDCfnbu3KnLL79cy5cvV3BwsCTpwQcfVHJysjIzM9W/f3+Fh3ufT164cKFWr16ttm3bSpImTpyopKQkvfPOO3ruued08cUX+7VGAEDlpcRFer2+Kq6pS5UAlVelEBQSElImXEinThn9+te/Vk5OjoqLiz3Bpjx//vOfvda77LLLNGzYME2dOlVz587VXXfd5bX+/fff7wlAktSgQQPdeeedGj9+vFavXu0JQf6sEUDd4TRF1QTaEY2vthZ4jgRJ0srTrg8KJMwznK7KJ/7Wrl2rZ599VsuWLdOuXbtUXFzstXzv3r0VXjwcHBysLl26lGnv3r27pk6dqrVr15YJQcnJyWXWv+SSSyRJBw4c8HuNAOqWTdch+MPqJ3q6XYLHvW/n6cm5/5KR0VVxTbVy636NnbtBqXFN9fLgsn+73cQ8w+mqNBv+8Y9/6LrrrpMk9e7dWwkJCWrUqJEcx9GcOXO0bt06FRVVfMdCVFSU6tUr+5zGZs2aSZIOHjxYZllERETZ4oNOlV9SUuL3GgEgkAXSXU5/vauT7p+5xusaoO4J0ZqSfqWahgXW3WHA6aoUgiZOnKiioiItW7ZMXbt29Vq2fPlyrVu3rlL97Nu3TydPniwThHbv3i3Jd+Cp6xoBAJXTNKy+pg9P5TlBOOdU6Wszvv/+ezVt2rRMuDhy5Ijy8vIq3U9xcbGWL19epv3zzz+XJCUlJVWlrFqpEQBQNQnNwtWrXTMCEM4ZVQpBrVq1UkFBgTZs2OBpKykp0cMPP6w9e/ZUacNPPvmk17U6mzZt0htvvKGIiAgNGDCgSn3VVo0AAOD8VaXTYffdd58WLFigbt266fbbb1doaKiWLFmiH3/8UWlpaZX+ktOYmBgdOHBASUlJ6tevnw4ePKh33nlHx44d06uvvlrm9ng3aiy1adMmTZo0SZJ09OhRT1vpF6NGR0frL3/5S7XrBQAA7qjSkaCbbrpJ77//vuLj4/XWW29pxowZuuyyy7Ry5Uq1atWq0v3Ur19fCxcu1DXXXKM333xTb775ptq2bas5c+b4fFCiGzWW2rVrl6fGWbNmSTp17VJp2/vvv1+jegEAgDscY0ydPtS89Gsz8vPz63KzZzVu3DiNHz9eixcv9vkFqtWVn5+vuLg4DRkyRNnZ2X7rFwAA+EeVjgSdz37961/LcRw99thjNerntddek+M4iouL81NlAACgNlj/1Kgzj/5069atRv0lJycrKyvL87omd7oBAIDaY/3pMAAAYKc6D0EAAACBgGuCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRD0i+zsbDmO4/lJT0+vdl+5ublefaWlpfmvUAAA4BdBbhcQaAYMGKCkpCQlJiZKkowxmj9/vj788EN98cUX+uGHH1RcXKyEhATdcccdevDBBxUaGurVR3x8vLKysiRJ48ePr/N9AAAAFXOMMcbtIgJBdna2hg4dqmnTpikjI8PTfuzYMTVo0EAhISFKS0tThw4ddOzYMeXk5Gjz5s1KSUnR0qVL1aBBA5/9Oo6jHj16aMmSJXWzIwBggc27C5W/74hioxoqoVm42+XgHMWRoApccMEFmjhxou699141adLE015cXKzf/OY3+uijj/TSSy/pkUceca9IAKgF+w4XuV1CGQVHjuvx2f/Siq37PW2pcU018eZERTas72JlvkU1CnG7BJSDEFSB4OBgjRkzxmf7n/70J3300UdaunQpIQg4Rxw5fsLtEs4ZnSbkul1CGfUcKTw0WFMHJSslLlJfbS3QmNnr1XvyZzoZgOc1vn6qj9slBLSG9d2NIYSgGggODpYkBQUxjMC5ot3YHLdLQA2cNNLTN3dQv44xkqR+HWNkZJQ5Y43LlfnGfCtf/qR+rm6fu8Nq4I033pAk9e7d2+VKAMAeKXGRXq+vimvqUiU413EIo5rmz5+vV155RZdffrmGDx/udjkAKonTE5UXqEcxvtpa4DkSJEkrT7s+KNAw3wIbIagaVq1apTvuuEMRERF67733FBLChW/AucLtaxDOJauf6Ol2CWXc+3aenpz7LxkZXRXXVCu37tfYuRuUGtdULw9Odru8MphvgY1Pp4rWrFmj3r17y3Ec5eTkqH379m6XBAC1IhDvbPrrXZ10/8w1XtcAdU+I1pT0K9U0LPDuDkNgIwRVQV5ennr16qWSkhItWLBAKSkpbpcEAFZpGlZf04en8pwg+AUhqJLy8vLUs2dPnThxQjk5OUpNTXW7JACwVkKzcMIPaowQVAmlAai4uFjz58/X1Vdf7XZJAACghghBFdi/f7969uypgoIC9e3bVwsXLtTChQu91mnSpIlGjx7tToEAAKBaCEEVOHTokAoKCiSdui1+/vz5ZdZp1aoVIQgAgHMMIagCsbGx4jtmAQA4//DE6DMMHTpUjuMoPT292n3k5ubKcRw5juPHygAAgD9xJOgXSUlJysrK8rxOTEysdl/x8fFefcXGxtakNAAAUAscw7keAABgIU6HAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRD0i+zsbDmO4/lJT0+vdl+5ublefaWlpfmvUAAA4BdBbhcQaAYMGKCkpCQlJiZ62j744APNmDFD69at0+7du3Xy5Em1atVKvXr10iOPPKIWLVp49REfH6+srCxJ0vjx4+u0fgAAUDmOMca4XUQgyM7O1tChQzVt2jRlZGR4Lfvtb3+r5cuXKyUlRTExMZKktWvXavHixYqIiNCyZcvUvn17n/06jqMePXpoyZIltbwHOJ9s3l2o/H1HFBvVUAnNwt0uBwDOSxwJqoRXX31VoaGhZdpff/11jRgxQuPGjdN7773nQmWoqX2Hi9wuwUvBkeN6fPa/tGLrfk9balxTTbw5UZEN67tYmbeoRiFulwAANUYIqgRfAUiSbrvtNo0YMULfffddHVdUdUeOn3C7hIDUaUKu2yV4qedI4aHBmjooWSlxkfpqa4HGzF6v3pM/08kAOmb79VN93C7hnNOwPn9ugUDDb2UNzJs3T5K8rh8KVO3G5rhdAirhpJGevrmD+nU8ddq1X8cYGRllzljjcmXemE9Vlz+pn9slADgDIagK5syZo7Vr1+rIkSPasGGDcnJyFBcXp6eeesrt0nAeSYmL9Hp9VVxTlyoBgPMbIagK5syZozfffNPzunPnzpo5c6bi4uJcrKpyOH3hWyAe0fhqa4HnSJAkrTzt+qBAwXwCcD4gBFVBdna2srOzdfDgQa1Zs0aPP/64OnXqpA8++EDXXXed2+WVi+sRfFv9RE+3S/By79t5enLuv2RkdFVcU63cul9j525QalxTvTw42e3yPJhPAM4H/CWrhoiICKWlpemTTz5R27Ztdffdd2vr1q0KDg52uzRUUaDd5fTXuzrp/plrvK4B6p4QrSnpV6ppWODcHQYA5wNCUA00btxYXbp00Zw5c/Tdd9/p8ssvd7sknOOahtXX9OGpPCcIAOoAIaiGduzYIUkKCmIo4T8JzcIJPwBQy/jusAoUFRVp+fLlPpdNmzZNK1euVOvWrZWQkFDHlQEAgJrg8EUFjh49qquvvlqJiYlKSkpSixYtdPDgQa1cuVJ5eXlq1KiRpk2b5naZAACgighBFQgLC9P48eO1ePFiLVq0SHv37lVwcLBiY2M1evRoPfDAA2rZsqXbZQIAgCoiBFUgODhYY8eO1dixY90uBQAA+BHXBJ1h6NChchxH6enp1e4jNzdXjuPIcRw/VgYAAPyJI0G/SEpKUlZWlud1Tb4PLD4+3quv2NjYmpQGAABqgWOMCaDvpgYAAKgbnA4DAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIegX2dnZchzH85Oenl7tvnJzc736SktL81+hAADAL4LcLiDQDBgwQElJSUpMTDzrOgcOHFD79u21Y8cO9enTR/Pnz/daHh8fr6ysLEnS+PHja7VeAABQPYSgMwwcOFAZGRnlrjNq1CgdPHjwrMvj4+M1btw4SYQg2GPz7kLl7zui2KiGSmgW7nY5AFAhQlAVffTRR5o+fbpeeOEFjRo1yu1yYKl9h4vcLsGj4MhxPT77X1qxdb+nLTWuqSbenKjIhvVdrKysqEYhbpcAIIAQgqpg//79+v3vf69Bgwbp3/7t3whBtezI8RNulxCwOk3IdbsEj3qOFB4arKmDkpUSF6mvthZozOz16j35M500blfn7eun+rhdQkBrWJ9/EmAXZnwVZGZmqqSkRC+88IIKCwvdLue8125sjtsloBJOGunpmzuoX8cYSVK/jjEyMsqcscblyspiTpUvf1I/t0sA6hQhqJJmz56td955R++++66ioqIIQcBpUuIivV5fFdfUpUoAoPIIQZWwd+9ejRw5UgMHDtTtt9/udjnW4NTF2QXaEY2vthZ4jgRJ0srTrg8KJMwpAKcjBFXCvffeq+LiYv31r391uxSrcH3C2a1+oqfbJXjc+3aenpz7LxkZXRXXVCu37tfYuRuUGtdULw9Odrs8L8wpAKfjL0IF5s6dq/fee0/Z2dlq3ry52+UAkgLrLqe/3tVJ989c43UNUPeEaE1Jv1JNwwLr7jAAOB0hqAJr1pz6w56RkeHz+UE5OTlyHEdXXHGF1q5dW7fFAQGgaVh9TR+eynOCAJxzCEEVSE5O1vDhw8u0Hz58WO+++64uueQS9enTRy1btnShOiBwJDQLJ/wAOKcQgirQv39/9e/fv0x7fn6+3n33XbVv316vvfaaC5UBAICa4AtUAQCAlQhBAADASpwOq6bY2FgZE2DfCQAAACqNI0FnGDp0qBzHUXp6erX7yM3NleM4chzHj5UBAAB/4kjQL5KSkpSVleV5nZiYWO2+4uPjvfqKjY2tSWkAAKAWOIZzOgAAwEKcDgMAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsFuV0AvBljVFhY6HYZAACcs8LDw+U4ToXrEYICTGFhoSIiItwuAwCAc9bBgwfVuHHjCtdzjDGmDupBJdXFkaBDhw7p0ksv1fbt2ys1SWzBuPjGuPjGuPjGuPjGuPhWW+PCkaBzlOM4dfYL0rhxY34ZfWBcfGNcfGNcfGNcfGNcfHNrXLgwGgAAWIkQBAAArEQIslBISIiysrIUEhLidikBhXHxjXHxjXHxjXHxjXHxze1x4cJoAABgJY4EAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiHoPLJr1y6NGDFCMTExCg0NVZs2bfTUU0/p+PHjVernxRdf1NChQ9WxY0cFBQXJcRwtWbLkrOtnZGTIcRyfP5dddlkN96rm3BoXf267NviztpycHKWlpalx48YKDw9XWlqacnJyfK4bCPPlq6++0o033qjIyEiFhYXpqquu0owZM6rUx8mTJ/XSSy+pY8eOatCggS688ELdfvvt2rx5c61utza5MS6xsbFnnQ8jR470x27VWE3H5aefftIzzzyjW2+9VXFxcZ79q+3t1jY3xsXf84WvzThP7Nq1S6mpqdq+fbsGDhyoNm3aaNmyZcrKytKXX36pefPmqV69ymXeUaNGSZJiYmJ04YUXateuXZV63/33368mTZp4tUVHR1dpP/zNzXHx57b9zZ+1vf3227rrrrsUHR2tIUOGyHEczZo1S3379tVbb72lwYMH+3yfW/NlyZIl6tOnj+rXr6/09HRFRETogw8+0ODBg5Wfn68xY8ZUqp+RI0fq1VdfVbt27XTfffdp9+7devfdd7VgwQL94x//ULt27Wplu7XFrXGRpIiICI0ePbpMe+fOnWu6WzXmj3H5+uuvNWbMGDmOo4SEBDVs2FBHjhyp9e3WJrfGRfLzfDE4L9x9991Gknn55Zc9bSdPnjRDhgwxkswbb7xR6b4+/vhjs3PnTmOMMffcc4+RZBYvXnzW9Uu3sXXr1uqWX2vcHBd/btvf/FXb/v37TZMmTUx0dLTZtm2bp33Hjh2mefPmpkmTJmb//v1e73FzvhQXF5tf/epXJiQkxOTl5XnaDx06ZNq3b2+CgoLMt99+W2E/n376qZFkunfvbo4dO+Zpz83NNY7jmGuvvbZWtltb3BoXY4xp1aqVadWqlV/2w9/8NS67du0yS5cuNYcOHTLGGNO2bVtT3j+/tsyXqo6LMf6fL4Sg88ChQ4dMSEiIiY+PNydPnvRatmPHDlOvXj1z9dVXV6vvczkEuTkutbntmvJnba+88oqRZMaPH19m2aRJk4wk88orr3i1uzlfcnJyjCQzdOjQMstmzpxpJJk//elPFfZz5513Gklm6dKlZZb17dvXSDLffPON37dbW9waF2MCOwTV1udW0T/2tsyXM7kRgjgddh748ssvVVRUpF69epU5nxoTE6MOHTpoxYoVOnbsmEJDQ2utjnnz5qmwsFAhISHq2LGj0tLSdMEFF9Ta9iri5rgEymdS27WVXhPVu3fvMsv69Omjxx57TEuXLtXvf//7MsvdmC/l1VvatnTp0kr1ExYWpq5du5ZZ1qdPH82fP19Lly5VmzZt/Lrd2uLWuJQqKirSm2++qR9//FGRkZG65pprdMUVV1RjT/zLrc/NlvlSXf6cL4Sg80DpBYcJCQk+lyckJGjdunXasmWLz/Px/pKZmen1uk2bNnrnnXeUnJxca9ssj5vjEiifiS/+rK28vkrbznZBrBvzpbx6IyMjFR0dXe6FzZL0888/a+fOnUpMTPQZ2nzttz+2W5vcGpdSu3btUkZGhldb3759NX36dFevK3Trc7NhvtSEP+cLd4edBw4ePCjp1MVivjRu3NhrPX/r0aOH/va3v2n79u06evSoNm7cqNGjR+v7779X7969tWPHjlrZbkXcHBe3P5Py+LO28voKCwvTBRdcUKYfN+dLZfa9ov2uzvj5Y7u1ya1xkaRhw4ZpyZIl2rNnjw4dOqTly5frhhtu0Pz589W/f38ZF7/e0q3PzYb5Ul3+ni+EoAASHR191lv/fP1UdHt2XRk6dKhuueUWXXLJJQoNDdVll12myZMn69FHH9W+ffs0efLkGvV/ro5LbTtXx6W25wvOLWPHjlWPHj0UHR2t8PBwpaam6uOPP1a3bt305Zdf6u9//7vbJSKA+Hu+cDosgNx5550qLCys9PrNmzeX9P9p/GzJ+9ChQ17r1ZXhw4fr6aef1hdffFGjfs7FcamLbQfCuJzeV1RUlNeyn3/+WSUlJZXeR3/Nl/JUZt8rqrc64+eP7dYmt8blbOrVq6ehQ4dq2bJl+uKLL9SvX78K31Mb3PrcbJgv/lST+UIICiAvvvhitd5X0bUXmzdvVr169RQfH1/t2qqj9NxsZZ77UJ5zcVzqYtuBMC4JCQlatWqVNm/eXCYEVXTt0Zn8NV/Kc/q+d+rUyWtZQUGB9u7dq2uuuabcPsLCwhQTE6OtW7eqpKSkzPUvvvbbH9utTW6NS3nqYj5UxK3PzYb54m/VnS+cDjsPdOnSRSEhIVq4cGGZ86E7d+7U+vXrlZqaWud3Ia1YsULSqSd8usHNcQnUz8TftfXo0UOStGDBgjLLSp8YXbpORepivpRXb2lbZert0aOHfv75Z59HrXztt7+2W1vcGpfyuP33Q3Lvc7NlvvhTteeL3262h6uq+vC7n3/+2WzcuNH88MMP5fZb0fNwdu7cab777rsy7f/7v/9rLrvsMiPJzJw5s+o75CdujUt1tl2X/DUu+/fvNxEREZV+WKLb86W4uNjEx8ebkJAQs2bNGk/76Q95O/05Nnv27DEbN240e/bs8ern9IcCFhUVedrLe1hiVbZb19walw0bNpiCgoIy9Xz++ecmNDTUhISEVPi7WJv8NS5nqszDEm2YL2eqaFxqY74Qgs4TO3bsMJdeeqlxHMfccsst5rHHHjNdu3Y1kkyfPn1MSUmJ1/qLFy82kkyPHj3K9PXMM8+YIUOGmCFDhpg2bdp4+iht+/zzz736cRzHdO/e3fzud78zjz76qLnjjjtMWFiYkWSGDBlS5oF8dcmtcanOtuuSP8dl+vTpRpKJjo42mZmZZtSoUaZZs2ZGkpk+fXqZftyeL59++qkJDg42jRo1Mr/73e/MQw89ZOLi4owkM2HCBK91s7KyjCSTlZVVpp8RI0YYSaZdu3bmkUceMXfffbcJCQkxERERZsOGDTXarhvcGJesrCzToEEDc9NNN5nMzEzz0EMPmT59+hjHccwFF1xgXn311drc5Urx17iU/p0YMmSIady4sWe+l/74CpQ2zJeqjEttzBdC0Hlkx44dZtiwYaZZs2amfv36pnXr1mb8+PFej68vVd4/aj169DCSzvozbdo0z7rbtm0zI0aMMB07djSRkZEmKCjIREVFmV69erl6BOh0boxLdbZd1/w1LsYY88knn5hrr73WNGrUyDRq1Mhce+21Zv78+WXWC5T5smLFCtO3b18TERFhGjRoYDp37mzeeuutMuuV98e7pKTEvPDCC6Z9+/YmJCTEREVFmVtvvbXc/0Ov7HbdUtfjsmTJEnP77beb1q1bm/DwcBMcHGwuueQSk56eblasWFEbu1gt/hiX8v526CxPULdhvlRlXGpjvji/FAEAAGAVLowGAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEr/B/29Z4Z6G8ngAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x440 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_forest(trace, var_names=[\"alpha\"], combined=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "73781c6c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([<AxesSubplot: title={'center': '94.0% HDI'}>], dtype=object)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGtCAYAAADuyIUrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzf0lEQVR4nO3de1xVdb7/8fcyboKIoKVkKpBoKRqhSI03bBJtrNRqinEaL2lT43G6WNaMlVs73k4z1djFZs6p5DyaSc3KnHKKwryMTV5RK1OzgnTyrqiYgAjf3x/+2McdoAh7s+G7X8/Hg0ft71r7uz7ruwHfrPVdaznGGCMAAIBGrom/CwAAAPAGQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwAqEGgAAYAVCDQAAsAKhBkCN7du3Tw8++KASExMVFhamVq1aafDgwcrOzr6gfl5++WU5jiPHcTRu3Lha17Nt2zb98pe/VGxsrMLCwnT55Zfr4Ycf1tGjR6tcv6ysTFOmTFG7du0UGhqq7t276+233662/y1btigoKEiTJk2qVX0rVqxw7+e55Ofnu9fLz8/3WDZ69Gj3soqv8PBwxcbG6pprrtGECRO0bNkyneuJN1OnTpXjOEpPT6/VfgCNhgGAGvjss89M69atjSQTGhpqevToYTp27GgkGUlm1qxZNernwIEDJiYmxv2+sWPH1qqejz/+2DRt2tRIMhdffLFJSUkx4eHhRpJJSEgw+/btq/SeRx55xEgykZGRplu3biYoKMg4jmOWLFlS5Tb69OljYmNjzfHjx2tV4/Lly937eS55eXnu9fLy8jyWjRo1ykgyl1xyiendu7fp3bu3SUtLM507dzbBwcHu91111VXms88+q7J/l8tlJJn+/fvXaj+AxoIjNQDO6/Tp07rtttu0f/9+paena/fu3dqwYYN27typZcuWKTIyUpMnT9aqVavO29eDDz6oo0ePasiQIbWup7CwUHfccYeKiop033336fvvv9fGjRu1a9cu9e7dW99++63Gjh3r8Z5Dhw7pueeeU4cOHbRz50599tlnWrZsmRzH0ZQpUypt47XXXtPq1av11FNPKTIysta1essNN9yg1atXa/Xq1VqzZo22b9+uY8eO6a233lK3bt20ZcsWXXPNNdq0aZO/SwX8hlAD4LyWLl2qr776SqGhocrKytLFF1/sXnbdddfpsccekzFG06ZNO2c/OTk5+tvf/qZ77rlHPXv2rHU9f/7zn3Xw4EFdeeWVeuaZZxQcHCxJatmypV5//XUFBQVp6dKlys3Ndb/n888/V3FxscaMGaPWrVtLkvr166c+ffpoy5YtKiwsdK9bWFioRx99VH369NGdd95Z6zp9rWnTprrlllu0du1a/fSnP9XJkyd1++23q6yszN+lAX5BqAFwXp988okkKTU1VR06dKi0/NZbb5V0Zg7JgQMHquyjuLhYv/nNb3TJJZdo5syZdaqnYh7M6NGjddFFF3ksa9++va6//npJ0ptvvulur6irItBUiI2NlSQdP37c3TZ16lQdOHBAzz//fJ3qrC9NmzbVX//6V4WGhurrr7/WokWL/F0S4BeEGgDnVVBQIElq27Ztlcsr2svLy7V+/foq15k+fbq+/vpr/eEPf1CLFi1qXcvp06e1ceNGSVLv3r2rXKeife3ate629u3bS5K++uorj3V37NihoKAgtWzZUtKZycfPP/+87rnnHiUnJ9e6zvrWpk0bDRs2TNKZI2tAICLUADivqKgoSdL3339f5fKz23fs2FFp+bZt2/SHP/xBffv21ciRI+tUS35+vkpLSyVJCQkJVa5T0b5z505321VXXaVLLrlEr7zyinJyclRYWKg5c+Zo8+bN6tevn8LCwiRJv/3tbxUVFaXp06fXqU5/6NOnjyRVGywB2wX5uwAADV9qaqokacOGDdq9e7fatWvnsfzsy6IrjupUMMbonnvuUXl5uebOnVvnWs7uPzo6usp1KtrPXjc8PFyzZs3S2LFjNXDgQHd7s2bN9PTTT0uSFi1apGXLlum///u/3X2Ulpbq0KFDatmypUJCQmpd9/ku6/aGis+lulOAgO0INQDOa+jQobr00ku1Z88ejRgxQm+88YZ7LsrSpUs1Y8YM97pFRUUe733llVf0z3/+Uw8//LCSkpLqXEtxcbH7/6sLGaGhoVXWctddd+nSSy/VvHnzdPDgQXXq1EkPPvigOnfurJMnT+rhhx9Wz549NXbsWBlj9Pjjj2vOnDn64YcfFBERofvuu08zZsyoVUCp7lSZJJWUlGjDhg0X3OePRURESJLHpGcgkBBqAJxXWFiYFi5cqJ/97GdavXq12rdvr86dO6ugoEB79uxR+/btlZycrFWrVqlZs2bu9x08eFCPPvqoLrvsMrlcLq/VUuHUqVMeryuUlJRIOjOB9scGDx6swYMHV2qfMWOGdu/erTfeeENNmjTR9OnTNXPmTN1444267bbb9Pbbb2vWrFmKiIjQY489dsF1r169utpl+fn5io+Pv+A+f+zEiROSpObNm9e5L6AxYk4NgBrp06ePcnNzddddd6lNmzbuCbf33nuvNmzY4L6MuE2bNu73PPLIIzpy5IieffZZj7BTF2efcvrxqa4ft1d3eurHvvnmGz399NMaPXq00tLSVFpaqqefflodO3bUkiVLNGrUKC1evFgdO3bU008/rdOnT9d9R3xg165dkqRLLrnEz5UA/sGRGgA11rFjR73yyiuV2k+fPq0tW7ZIknr06OFur7gR3IQJEzRhwgSP91QcVXj99df13nvvSTrzGIbziYuLU3BwsEpLS/Xtt9+6T4Od7dtvv5UkJSYm1mS3dP/99yssLEyzZ8+WJG3fvl1Hjx7ViBEj1KTJmb/9mjRpooyMDM2dO1c7duxQ165da9R3fao4GtSrVy8/VwL4B6EGQJ1lZ2frxIkTuvTSS5WSklJp+f79+6t9b1FRUaW5L+cSFBSklJQUrV27Vp988kmVc1Uq7quTlpZ23v7ee+89LV26VHPmzHEf4agIXD++k3DF6+qeLeVPe/fu1d///ndJqtPdmoHGjNNPAOrk1KlT7scM/OY3v/G4Gd7mzZtljKnyq2KOTcWkXHOOBzL+2C233CJJysrKqnT33F27diknJ0fS/90UsDolJSV64IEHlJSUpPHjx7vbK64i+uabbzzWr3jdqlWrGtdaH4qKivSrX/1KJSUl6tSp03n3G7AVoQZAjfzjH//wuJmdJO3evVvDhg1Tbm6uunTpUuunWVflzTffVFxcnPveK2e799571apVK23btk0TJ05037fm8OHDGjFihE6fPq0bbrjB41RYVZ566il98803euGFFxQU9H8Hrtu2bat27drp3Xff1WeffSbpzGMW3n33XbVp06bGp7V8raioSIsXL1ZaWpqWLVumiIgIvfHGG5XusgwECk4/AaiRDz/8UHPmzFF0dLTi4uJUXFys7du3yxijLl266MMPP3RfSu0NJ06c0HfffVflsubNm2vBggW68cYb9dxzz2n+/Plq3769tm3bppMnTyouLk6vvvrqOfvftWuXZs+erczMTPXv399jmeM4mjp1qsaOHavU1FR17txZX331lUpKSuRyudzzbOrT+++/7w54ZWVlKigo0LfffusOdMnJyXrttde8ctk80FgRagDUyLBhw7R3716tW7dO27ZtU2hoqFJTU3XHHXfoP/7jP7waaGripz/9qTZs2KDp06fr448/1ueff662bdtq+PDhevzxx8975dPEiRPlOI7++Mc/Vrn8rrvuUnFxsZ599llt375dHTp00EMPPaR7773XF7tzXgcOHHDfVC8sLExRUVFKSUlRz549NXz4cP30pz/1S11AQ+KYCzmRDQAA0EAxpwYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwArcp6aeGGNUWFjo7zIAAGi0IiMj5ThOtcsJNfWksLBQUVFR/i4DAIBG69ixY2revHm1y7n5Xj3hSA0AAHVzviM1hBoAAGAFJgoDAAArEGoAAIAVCDUAAMAKhBoAAGAFQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwAqEGgAAYAVCDQAAsAKhBgAAWIFQAwAArECoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAKhBoAAGAFQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwAqEGgAAYAVCDQAAsAKhBgAAWIFQAwAArECoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArXHCoWbFihRzH0dSpU31QjndNnTpVjuO4v373u9/Vqb+XX37Zo7/Ro0d7p1AAAFBnQf4u4Mfi4uIkSfn5+V7rc9SoUYqLi1OfPn082jdv3qw33nhDGzduVG5urg4dOqT+/ftrxYoVVfaTkpIil8ulo0ePas6cOV6rDwAA1F2DCzW+MHr0aKWnp1dqf+eddzRr1iyFhISoU6dOOnTo0Dn7SUlJUUpKivLz8wk1QADZub9Q+YdPKq5luBJbR/q7HADVCIhQU52f//znuvnmm9WtWzcdPnxYsbGx/i4JwAU4fKLEp/0XnDylxxZ/obV5R9xtafExmjE8SdHhIT7ZZstmoT7pFwgEdQo1q1at0hNPPKGNGzcqJCREgwcP1lNPPaXLLrvMY70DBw5o1qxZevfdd7V7925FRkaqf//+mjZtmpKSkiSdOd0UHx/vfo/jOO7/d7lcmjp1qk6dOqW//OUveu+99/Tll1/qwIEDioqKUp8+ffTEE0/o6quvvqD6u3btWoe9Bxq/k6dO+7uEOukxPcen/TdxpMiwYL04IkWp8dFan1egyYs/V8azq1RufLPNL58c5JuO60l4SED/rQw/q/V335o1azRr1iwNGTJE9913n3JzczV//nytXr1a69evV+vWrSVJ33zzjdLT0/X9998rIyNDw4YN04EDB/TWW28pOztby5YtU1pamlq0aCGXy6U//elPkqQHHnjAva2KU0dHjhzRAw88oL59++pnP/uZoqOj9e233+rvf/+73n//fa1atUqpqam1Hgwg0HSZku3vEhq0ciPNHN5NQ7qfOYo7pHusjIwmvL7JZ9ts7J9J/uwh/i4BAazWoSY7O1svv/yyxo4d62578skn5XK5NHnyZL3yyiuSpJEjR2rfvn3Kzs7WwIED3es+/vjj6tmzp+6++2599tlnatGihaZOnaqsrCxJqvLqqujoaO3atUtt27b1aN+6dauuueYaTZ48WR999FFtdwkAKkmNj/Z43Ss+xk+VADifWoeazp0766677vJomzRpkl544QXNnz9fL730krZu3ap//etfGjt2rEegkaROnTrp7rvv1jPPPKMvvvjCfRrqXEJDQysFGunMaaQBAwYoOztbpaWlCg4Oru1uAQGlsZ/qqI+jGuvzCtxHaiRp3Vnza3yhsX8mgD/VOtT07t3bY96LJDVt2lQ9evTQBx98oK+++kpr1qyRJO3bt6/KIy/bt293/7cmoUY6cxn2U089pdWrV2vfvn0qLS31WH7o0CEm/AI11NjnP2x8/Hqf9j/+b7l6YskXMjLqFR+jdXlHNGXJVqXFx2juL1N8ss3G/pkA/lTrn55LLrmkyvaKuTTHjh3TkSNn/qJZunSpli5dWm1fP/zwQ422+a9//UvXXXedJCkjI0OJiYlq1qyZHMfRO++8oy1btqikxLdXQwBoOHx9pdBLd/bQ/Qs2ecyh6ZvYSnMyr1ZMhG+ufgJQe7UONQcOHKiyff/+/ZKkqKgoNW/eXJL0/PPPa8KECbXdlNuMGTNUUlKi1atXq3fv3h7L1qxZoy1bttR5GwBQISYiRK+NTeM+NUAjUetnP33yyScyxvOaxqKiIm3cuFFNmzZVp06dlJaWJkn69NNPa9zvRRddpLKysiqXffPNN4qJiakUaE6ePKnc3NwL3AMAqJnE1pEa2KU1gQZo4Godanbs2KFXX33Vo+0Pf/iDDh48qF/84hcKCQlRr169lJaWpvnz52vhwoWV+igvL9fKlSs92mJiYnTo0CEVFxdXWr9Dhw4qKCjQ1q1b3W1lZWV6+OGHdfDgwdruCgAAsECtTz9lZGRo/PjxWrp0qa644grl5uYqOztb7dq108yZM93rzZ8/XwMGDFBmZqb+9Kc/qUePHgoLC9OuXbv06aef6uDBgx4B5rrrrtOGDRt00003qW/fvgoJCVGfPn3Up08f/fa3v9WHH36oPn366Pbbb1dYWJhWrFih77//Xunp6dU+s6k627dv1+zZsyWdOcpU0VbxoMpWrVrpj3/8Y22HCAAA1CdzgZYvX24kGZfLZVauXGn69u1rwsPDTYsWLUxmZqbZtWtXpfccOXLEPP744yYpKck0bdrUNGvWzCQmJpoRI0aYt99+22PdwsJCc/fdd5vY2FjTpEkT97YqvPnmmyYlJcWEh4ebVq1amdtvv9188803ZtSoUUaSycvLc6/rcrmMJLN8+fJz7kt1Xx06dKjyfXl5eUaSGTVq1AWOHgAA8BXHGOOjm33739SpUzVt2jQtX768ygda1lbFIx1GjRrlvlkgAADwr1rPqWlMBgwYIMdx9Lvf/a5O/bz88styHMfjGVUAAKBhsPouTz8+OtOnT5869ZeSkiKXy+V+nZycXKf+AACA91h9+gkAAASOgDj9BAAA7EeoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAKhBoAAGAFQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwAqEGgAAYAVCDQAAsAKhBgAAWIFQAwAArECoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWsD7UZGVlyXEc91dmZmat+8rJyfHoKz093XuFAgCAOgnydwH1ZejQoUpOTlZSUpIkyRijDz74QH//+9/1ySef6LvvvlNpaakSExN1xx13aOLEiQoLC/PoIyEhQS6XS5I0bdq0et8HAABQPccYY/xdhC9lZWVpzJgxmjdvnkaPHu1uLy4uVtOmTRUaGqr09HR169ZNxcXFys7O1s6dO5WamqqVK1eqadOmVfbrOI769++vFStW1M+OAACAcwqYIzU/dtFFF2nGjBkaP368WrRo4W4vLS3VrbfeqnfffVcvvPCCJk2a5L8i4RM79xcq//BJxbUMV2LrSH+XAwDwkoANNcHBwZo8eXKV7b///e/17rvvauXKlYSaenT4RIlP+y84eUqPLf5Ca/OOuNvS4mM0Y3iSosNDfLLNls1CfdIvAKCygA015xIcHCxJCgpq+MNz8tRpf5fgNT2m5/i0/yaOFBkWrBdHpCg1Plrr8wo0efHnynh2lcp9dBL2yycH+abjehIe0vB/BgCgAr+xqvDqq69KkjIyMvxcyfl1mZLt7xIajXIjzRzeTUO6x0qShnSPlZHRhNc3+Wybjf3zyZ89xN8lAECNWX9J94X64IMP9Je//EVXXnmlxo4d6+9y4GWp8dEer3vFx/ipEgCAt3Gk5iwbNmzQHXfcoaioKC1atEihoQ1/PkRjP71xtvo4qrE+r8B9pEaS1p01v8YXbPp8AKChI9T8f5s2bVJGRoYcx1F2dra6du3q75JqxKY5Dxsfv96n/Y//W66eWPKFjIx6xcdoXd4RTVmyVWnxMZr7yxSfbNOmzwcAGjp+40rKzc3VwIEDVVZWpg8//FCpqan+Likg+fpKoZfu7KH7F2zymEPTN7GV5mRerZgI31z9BACoPwEfanJzc3X99dfr9OnTys7OVlpamr9Lgo/ERITotbFp3KcGACwV0KGmItCUlpbqgw8+0LXXXuvvklAPEltHEmYAwEIBG2qOHDmi66+/XgUFBRo8eLA++ugjffTRRx7rtGjRQg888IB/CgQAABckYEPN8ePHVVBQIOnMZdwffPBBpXU6dOhAqAEAoJEI2FATFxcny5/lCQBAQAmYm++NGTNGjuMoMzOz1n3k5OTIcRw5juPFygAAgDdYf6QmOTlZLpfL/TopKanWfSUkJHj0FRcXV5fSAACAFzmGczAAAMACAXP6CQAA2I1QAwAArECoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAKhBoAAGAFQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwAqEGgAAYAVCDQAAsAKhBgAAWIFQAwAArECoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsYH2oycrKkuM47q/MzMxa95WTk+PRV3p6uvcKBQAAdRLk7wLqy9ChQ5WcnKykpCR329tvv63XX39dW7Zs0f79+1VeXq4OHTpo4MCBmjRpktq2bevRR0JCglwulyRp2rRp9Vo/AAA4N8cYY/xdhC9lZWVpzJgxmjdvnkaPHu2x7Fe/+pXWrFmj1NRUxcbGSpI2b96s5cuXKyoqSqtXr1bXrl2r7NdxHPXv318rVqzw8R4AAICaCJgjNVX5n//5H4WFhVVqf+WVVzRu3DhNnTpVixYt8kNlgH127i9U/uGTimsZrsTWkf4uB4CFAjrUVBVoJOnnP/+5xo0bp6+//rqeKwLq3+ETJT7tv+DkKT22+AutzTvibkuLj9GM4UmKDg/x2XZbNgv1Wd8AGqaADjXVWbp0qSR5zL8BaurkqdP+LuGC9Jie49P+mzhSZFiwXhyRotT4aK3PK9DkxZ8r49lVKvfhye8vnxzku879IDyEX9fA+fBTIumdd97R5s2bdfLkSW3dulXZ2dmKj4/Xk08+6e/S0Ah1mZLt7xIalHIjzRzeTUO6n5m3NqR7rIyMJry+yafbte1zyJ89xN8lAA0eoUZnQs3//u//ul/37NlTCxYsUHx8vB+rAuyRGh/t8bpXfIyfKgFgM0KNzlwhlZWVpWPHjmnTpk167LHH1KNHD7399tu67rrr/F0eGpnGdtqjPo5orM8rcB+pkaR1Z82v8ZXG9jkAqDtCzVmioqKUnp6u999/X507d9bIkSOVl5en4OBgf5eGRqSxzX3Y+Pj1Pu1//N9y9cSSL2Rk1Cs+RuvyjmjKkq1Ki4/R3F+m+Gy7je1zAFB3/NRXoXnz5rrmmmv0zjvv6Ouvv9aVV17p75IAn/H1VUIv3dlD9y/Y5DGHpm9iK83JvFoxEb67+glA4CHUVGPPnj2SpKAghgioi5iIEL02No371ADwOeuf/VSdkpISrVmzpspl8+bN07p169SxY0clJibWc2WAnRJbR2pgl9YEGgA+E7CHIYqKinTttdcqKSlJycnJatu2rY4dO6Z169YpNzdXzZo107x58/xdJgAAqKGADTURERGaNm2ali9frmXLlunQoUMKDg5WXFycHnjgAT344INq3769v8sEAAA1FLChJjg4WFOmTNGUKVP8XQoAAPCCgJlTM2bMGDmOo8zMzFr3kZOTI8dx5DiOFysDAADeYP2RmuTkZLlcLvfrujzPKSEhwaOvuLi4upQGAAC8yDHG+PCRcgAAAPUjYE4/AQAAuxFqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAKhBoAAGAFQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwAqEGgAAYAVCDQAAsAKhBgAAWIFQAwAArECoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAK1oearKwsOY7j/srMzKx1Xzk5OR59paene69QAABQJ0H+LqC+DB06VMnJyUpKSqp2naNHj6pr167as2ePBg0apA8++MBjeUJCglwulyRp2rRpPq0XAABcmIAJNcOGDdPo0aPPuc59992nY8eOVbs8ISFBU6dOlUSoAYBAtXN/ofIPn1Rcy3Alto70dzk4S8CEmvN599139dprr+m5557Tfffd5+9yAAC1cPhEic/6Ljh5So8t/kJr846429LiYzRjeJKiw0N8tt2WzUJ91rdtCDWSjhw5ol//+tcaMWKEbrrpJkINgIBy8tRpf5fgNT2m5/is7yaOFBkWrBdHpCg1Plrr8wo0efHnynh2lcqNzzarL58c5LvOfSQ8xD/xglAjacKECSorK9Nzzz2nwsJCf5cDAPWqy5Rsf5fQKJQbaebwbhrSPVaSNKR7rIyMJry+yafbbYyfT/7sIX7ZbsCHmsWLF2v+/PlauHChWrZsSagBAFQrNT7a43Wv+Bg/VYKqBHSoOXTokO69914NGzZMt99+u7/LAQC/aIynN6rj66Ma6/MK3EdqJGndWfNrfMWmz8fXAjrUjB8/XqWlpXrppZf8XQoA+I2/5j/4wsbHr/dZ3+P/lqsnlnwhI6Ne8TFal3dEU5ZsVVp8jOb+MsVn27Xp8/G1gB2pJUuWaNGiRcrKylKbNm38XQ4AwAt8eaXQS3f20P0LNnnMoemb2EpzMq9WTITvrn5CzQVsqNm06cw35ejRo6u8f012drYcx9FVV12lzZs3129xAIAGJyYiRK+NTeM+NQ1YwIaalJQUjR07tlL7iRMntHDhQl122WUaNGiQ2rdv74fqAAANVWLrSMJMAxWwoebmm2/WzTffXKk9Pz9fCxcuVNeuXfXyyy/7oTIAAFAb1j/QEgAABAZCDQAAsELAnn6qTlxcnIzx4f2uAQCATwTMkZoxY8bIcRxlZmbWuo+cnBw5jiPHcbxYGQAA8Abrj9QkJyfL5XK5XyclJdW6r4SEBI++4uLi6lIaAADwIsdwrgUAAFggYE4/AQAAuxFqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAKhBoAAGAFQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwAqEGgAAYAVCDQAAsAKhBgAAWIFQAwAArECoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVrA81WVlZchzH/ZWZmVnrvnJycjz6Sk9P916hAACgToL8XUB9GTp0qJKTk5WUlORuy8rK0pgxY6p9z969e9WmTRv364SEBLlcLknStGnTfFcsAAC4YAETaoYNG6bRo0dXuawi8PxYs2bNPF4nJCRo6tSpkgg1AAA0NAETas7lXIEHQN3s3F+o/MMnFdcyXImtI/1dDgCLEWqAAHX4RIlP+y84eUqPLf5Ca/OOuNvS4mM0Y3iSosNDfLbdls1CfdY3gIaNUCNp8+bNeuaZZ3T69GldfvnlysjIUGQkf1GiaidPnfZ3CV7RY3qOT/tv4kiRYcF6cUSKUuOjtT6vQJMXf66MZ1ep3Phuu18+Och3nftYeAi/koG64CdI0pw5czxeR0VF6YUXXtCdd97pp4rQkHWZku3vEhqFciPNHN5NQ7rHSpKGdI+VkdGE1zf5dLuN+fPJnz3E3yUAjZr1l3SfS0JCgubOnauvv/5aJ0+eVH5+vl588UU1adJEI0eO1Pvvv+/vEoFGLTU+2uN1r/gYP1UCIBAE9JGafv36qV+/fu7XHTp00Pjx49WxY0cNGjRIU6ZM0Q033ODHCtEQNebTG2erjyMa6/MK3EdqJGndWfNrfMWWzwfAhQvoUFOdjIwMtWvXThs3blRJSYlCQ5l4iP9jy7yHjY9f79P+x/8tV08s+UJGRr3iY7Qu74imLNmqtPgYzf1lis+2a8vnA+DC8dNfjVatWmn37t0qKioi1MBKvr5K6KU7e+j+BZs85tD0TWylOZlXKybCd1c/AQhchJoqHD9+XNu3b1eLFi0UFRXl73KARikmIkSvjU3jPjUA6k1ATxT+5JNPKrUVFRXp7rvvVlFRkTIzM+U4jh8qA+yR2DpSA7u0JtAA8LmAPlLTp08fdenSRT169NCll16qAwcOKCcnR7t379ZVV12lmTNn+rtEAABQQwEdaiZOnKg1a9YoOztbBQUFCg0N1ZVXXqkJEybot7/9rZo2bervEgEAQA0FdKh5+umn/V0CAADwkoCZUzNmzBg5jqPMzMxa95GTkyPHcZhnAwBAA2T9kZrk5GS5XC7366SkpFr3lZCQ4NFXXFxcXUoDAABe5BhjfPhoOQAAgPoRMKefAACA3Qg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAKhBoAAGAFQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsQagAAgBUINQAAwAqEGgAAYAVCDQAAsAKhBgAAWIFQAwAArECoAQAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAK1oearKwsOY7j/srMzKx1Xzk5OR59paene69QAABQJ9aHmgpDhw6Vy+XSbbfdVmnZgQMHNHHiRCUmJiosLEwtW7bUtddeq5deesljvYSEBLlcLrlcrvoqGwAA1JBjjDH+LsKXsrKyNGbMGM2bN0+jR4+utHzz5s3KyMhQQUGBhgwZoiuvvFInTpzQtm3bFBISon/84x9V9us4jvr3768VK1b4dgcAAECNBPm7AH8qLCzU0KFDJUkbN25U9+7dPZafPn3aH2UBVdq5v1D5h08qrmW4EltH+rscAGhwAjrUzJ07V7t27dIrr7xSKdBIUlBQQA8PaujwiRKf9l9w8pQeW/yF1uYdcbelxcdoxvAkRYeH+HTbLZuF+rR/APCmgP5Xe+HChXIcR7feeqt27NihDz/8UEVFRbriiis0ePBghYT49h8MnHHyVOM+ItZjeo5P+2/iSJFhwXpxRIpS46O1Pq9Akxd/roxnV6ncxyePv3xykG83UE/CQwL6Vx0QMAL2J/3UqVP67LPPdPHFF+uFF17QlClTVF5e7l6ekJCgd955R926dfNjlYGhy5Rsf5fQoJUbaebwbhrSPVaSNKR7rIyMJry+yefbtuWzyZ89xN8lAKgHAXP1048dOXJEZWVlOnz4sKZNm6annnpK+/fv17///W898cQTysvL00033aTi4mJ/lwooNT7a43Wv+Bg/VQIADVfAHqmpOCpTVlamCRMm6KGHHnIve/LJJ/XVV19p4cKFevPNN3XnnXf6q8yA0NhPcdTH0Yz1eQXuIzWStO6s+TW+1Ng/GwCBJWBDTVRUlPv/b7755krLb7rpJi1cuFAbNmwg1PhYY5/vsPHx633a//i/5eqJJV/IyKhXfIzW5R3RlCVblRYfo7m/TPHpthv7ZwMgsATsb6yIiAi1bdtW33//vVq0aFFpeUVbUVFR/RaGRsfXVwi9dGcP3b9gk8ccmr6JrTQn82rFRDCZHQAqBGyokaTrrrtOr732mr788kulpHj+xfvll19KkuLi4vxQGfB/YiJC9NrYNO5TAwDnEbAThSXp3nvvlSTNnj1bR48edbfv27dPc+bMUZMmTXTrrbf6qTrAU2LrSA3s0ppAAwDVCOgjNT/5yU80ceJEPfPMM+revbtuuukmlZaWasmSJTpw4IBmzpypTp06+btMAABQAwEdaiTp6aefVrdu3fTiiy+6n+h99dVX689//rOGDx/u7/IAAEANBXyokaTRo0dX+bBLAADQeATMnJoxY8bIcRxlZmbWuo+cnBw5jiPHcbxYGQAA8Abrj9QkJyfL5XK5XyclJdW6r4SEBI++uDIKAICGwzHG+PiReAAAAL4XMKefAACA3Qg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAKhBoAAGAFQg0AALACoQYAAFiBUAMAAKxAqAEAAFYg1AAAACsE+buAQGGMUWFhob/LAACg0YqMjJTjONUuJ9TUk8LCQkVFRfm7DAAAGq1jx46pefPm1S53jDGmHusJWIF8pOb48eNq166ddu/efc5vxkDF+Jwb43NujM+5MT7n1tjGhyM1DYTjOI3iG8aXmjdvHvBjcC6Mz7kxPufG+Jwb43NutowPE4UBAIAVCDUAAMAKhBr4XGhoqFwul0JDQ/1dSoPE+Jwb43NujM+5MT7nZtv4MFEYAABYgSM1AADACoQaAABgBUINAACwAqEGAABYgVADr9q3b5/GjRun2NhYhYWFqVOnTnryySd16tSpGvexc+dOzZw5U/369dOll16qkJAQtWvXTiNHjtT27dt9WL3veWN8JOn555/XmDFj1L17dwUFBclxHK1YscI3RXvZ+vXr9bOf/UzR0dGKiIhQr1699Prrr19QH+Xl5XrhhRfUvXt3NW3aVBdffLFuv/127dy500dV15+6js+BAwc0a9Ys3XbbbYqPj5fjOOe8A2tjU9fxWb16tR566CH16NFDLVu2VFhYmK644go9+uijOnr0qO8Kryd1HZ8VK1ZoxIgRuvLKK9WiRQuFh4erc+fOuuuuu7Rjxw4fVu4lBvCSvXv3mvbt2xvHcczw4cPNo48+anr37m0kmcGDB5uysrIa9XPHHXcYSSYpKcnce++95pFHHjE33HCDkWSaNm1qVq1a5eM98Q1vjY8xxkgykkxsbKxp06aNkWSWL1/uu+K9ZPny5SYkJMQ0a9bMjBs3zjz00EMmPj7eSDIzZsyocT933323kWS6dOliJk2aZEaOHGlCQ0NNVFSU2bp1qw/3wLe8MT7Lly83kozjOKZTp04mPDzc2PKr3hvj07p1a3PRRReZ/v37mwceeMA8+OCD5uqrrzaSzOWXX27279/v473wHW+Mz2OPPWbatWtnbrvtNnP//febSZMmmRtvvNFcdNFFJjQ01Hz88cc+3ou6seM7HQ3CyJEjjSQzd+5cd1t5ebkZNWqUkWReffXVGvUzb948s3nz5krt8+fPd/9D1hh5a3yMMea9994ze/fuNcYYc8899zSKUFNaWmouv/xyExoaanJzc93tx48fN127djVBQUHmq6++Om8/H3/8sZFk+vbta4qLi93tOTk5xnEc069fP5/U72veGp99+/aZlStXmuPHjxtjjOncubMVocZb4zN79myzZ88ej7by8nLzm9/8xkgy48eP93rt9cFb41NUVFRle05OjpFkevbs6bWafaHxf6ejQTh+/LgJDQ01CQkJpry83GPZnj17TJMmTcy1115b5+106tTJSDIHDx6sc1/1yZfj01hCTXZ2tpFkxowZU2nZggULjCTz+9///rz9/OIXvzCSzMqVKystGzx4sJFkduzY4ZWa65O3xufHbAk1vhqfCnv27DGSTNeuXetSpt/4enyMMSY6Otq0aNGiTn34GnNq4BWffvqpSkpKNHDgwErn72NjY9WtWzetXbtWxcXFddpOcHCwJCkoqHE9i7W+xqchq5jzk5GRUWlZRdvKlStr1E9ERIR69+5dadmgQYNq3E9D463xsZWvx6ex/m6p4Ovx+fTTT1VQUKCkpKRa91EfGuenhwanYoJmYmJilcsTExO1ZcsWffvtt+rSpUuttrFu3Tpt3bpVqampatGiRW1L9Yv6GJ+G7lxjEB0drVatWp13ou8PP/ygvXv3KikpSRdddFGl5RV9N8YJw94YH5v5enxeffVVSVWHgsbA2+OzYsUKrVixQiUlJdq5c6fee+89tWrVSs8++6zXavYFQg284tixY5KkqKioKpdXPNK+Yr3a9D9q1Cg1adJETz31VO2K9CNfj09jUJMx+Pe//13nPs5erzHxxvjYzJfjs3nzZk2bNk2XXHKJHnnkkVrX6E/eHp8VK1Zo2rRp7tcdO3bUggUL1KNHj7oV6mOcfoKHVq1auS8BrclXfVxGXFxcrFtuuUXbt2/Xf/7nfyo9Pd3n26xOQxwfALWXl5enG2+8UWVlZVqwYIFatWrl75IahKlTp8oYoxMnTmjdunW64oor1Lt37wu+/UJ940gNPPziF79QYWFhjddv06aNpP/766C6v5CPHz/usV5NlZSUaPjw4fr444/1+9//XpMnT76g93tbQxufxqQmY3C+/bd5HL0xPjbzxfh89913GjBggA4ePKi33npLAwYMqHOd/uKr75+IiAilpqZq8eLF6tmzp379619r4MCBuvjii+tUr68QauDh+eefr9X7zjeXYefOnWrSpIkSEhJq3GdxcbGGDRum7OxsPfLII5o5c2atavOmhjQ+jc3ZY/DjQ9gFBQU6dOiQfvKTn5yzj4iICMXGxiovL09lZWWV5tWcb+5SQ+aN8bGZt8cnPz9fAwYM0J49e7Ro0SLdeOONXq23vvn6+ycoKEgDBgzQli1btGHDBt1www11qtdXOP0Er7jmmmsUGhqqjz76SMYYj2V79+7V559/rrS0NIWFhdWov7MDzcMPP6z/+q//8kXZ9cbb49MY9e/fX5L04YcfVlpW0Vaxzvn6+eGHH/TJJ59UWpadnV3jfhoab42Prbw5Pvn5+UpPT9f333+vhQsXaujQod4r1E/q4/tnz549khr4FWJ+vqQcFrnQm8v98MMPZtu2bea7777zaC8qKjIZGRlGkpk4cWK91F4fvDU+P9ZY7lNTWlpqEhISTGhoqNm0aZO7/eybg519f5mDBw+abdu2Vbon0dk33yspKXG323DzPW+Mz4/Zcp8ab41PXl6e6dChgwkKCjJvvfVWfZXvc94an5UrV1a6l5YxZ+6DExwcbKKiosyJEyd8th911fi/09Fg7Nmzx7Rr1844jmNuueUW87vf/c79GIBBgwZVegxAxe3c+/fv79Fe8Y98mzZtjMvlqvIrLy+v/nbMS7w1PsYYM2vWLDNq1CgzatQo9w0JBw0a5G775z//WU97dWE+/vhjExwcbJo1a2buvvtuj9u4T58+3WNdl8tlJBmXy1Wpn3Hjxln5mARvjU/F98GoUaNM8+bNjSSPtsZ288oK3hifDh06GEnmmmuuqfb3S2PljfGJiooyl19+ucnMzDSTJk0yEyZMMP369TOSTHBwsFm0aFE97tGFI9TAq/bs2WPuuusu07p1axMSEmI6duxopk2b5nE7+wrV/aPdv39/97ONqvtq6EclquON8THm/GM0b9483+9MLa1du9YMHjzYREVFmaZNm5qePXuav/71r5XWO9c/2mVlZea5554zXbt2NaGhoaZly5bmtttua5R3Ev4xb4zP+X5+GuMfBRXqOj7nG5vGflSrruPzpz/9yQwePNhcdtllJjQ01ISFhZnExEQzbtw488UXX9TTXtSeY8yPTvADAAA0QkwUBgAAViDUAAAAKxBqAACAFQg1AADACoQaAABgBUINAACwAqEGAABYgVADAACsQKgBAABWINQAAAArEGoAAIAVCDUAAMAKhBoAAGCF/wdM02E+7HMrhgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 600x460 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_forest(trace, var_names=[\"beta\"], combined=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b9ba2d64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([<AxesSubplot: title={'center': '94.0% HDI'}>], dtype=object)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoQAAAI3CAYAAAAC6UqUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAACLC0lEQVR4nOzdfVxUZf4//tdRuZFbGUxFSGdUcg0URBH9oYVGqLDJOEuKxhaElZIZ9bE7VG5cMds+bh+tdLvTMb+SRnETGpCaN1AmKiBqummC22LqIuAd3oBzfn+4zDqBicw5DDPzej4e86hznTPX+80ZHd6ec53rEkRRFEFEREREVquLqRMgIiIiItNiQUhERERk5VgQEhEREVk5FoREREREVo4FIREREZGVY0FIREREZOVYEBIRERFZORaERERERFaOBSERERGRlWNBSERERGTlWBASEcnszJkzeOmll+Dt7Q17e3v07NkTkyZNQmFh4T318/HHH0MQBAiCgFmzZrU7n6NHj+KJJ56Ah4cH7O3tMXDgQMyfPx/19fWtHn/z5k0kJyfj/vvvh52dHYYNG4asrKw79n/w4EF069YNr7zySrvy27lzp/7n/D1VVVX646qqqgz2xcbG6vc1vxwcHODh4YHRo0dj7ty52L59O35v9dbU1FQIgoCQkJB2/RxE5oQFIRGRjA4dOgR/f3/83//9H3755Rf4+vrCzc0NhYWFmDRpEpYtW9amfv7973/jtddeMzqfHTt2YMSIEcjIyMDNmzfh4+ODM2fOYPny5RgxYgTOnj3b4j1JSUn4y1/+ggsXLmDw4ME4evQooqKi8NVXX7UaY+7cuejVqxeSk5ONztdYvXr1QnBwMIKDgzFs2DC4urqitLQU77//PkJDQzF8+HAcOnTI1GkSmRwLQiIimTQ1NSEqKgpnz55FSEgIfvnlF+zfvx/Hjx/H9u3b4ezsjKSkJOzevfuufb300kuor69HREREu/O5dOkSpk+fjqtXr2LevHmorq7GgQMH8M9//hPBwcE4efIk4uPjDd5TU1ODlStXon///jh+/DgqKiqwfft2CILQasG3fv16FBcX469//SucnZ3bnatUJk+ejOLiYhQXF+OHH37AsWPHcOHCBXz55ZcYOnQoDh48iNGjR6OsrMzUqRKZFAtCIiKZbNmyBT/99BPs7Oyg1Wpx33336fdNmDABCxYsgCiKSEtL+91+tm3bhg0bNuC5557DyJEj253P3//+d/z73//GkCFD8Le//Q02NjYAAHd3d2RkZKBbt27YsmULSktL9e85dOgQrl27hri4OPTu3RsA8NBDD2Hs2LE4ePAgLl26pD/20qVLeO211zB27FjExMS0O0+5de/eHRqNBnv37sUjjzyChoYGTJs2DTdv3jR1akQmw4KQiEgm3333HQAgMDAQ/fv3b7H/T3/6E4BbY+bOnTvXah/Xrl3DnDlz0KtXLyxdutSofJrH/cXGxqJr164G+/r164fQ0FAAwBdffKFvb86ruRhs5uHhAQC4ePGivi01NRXnzp3Du+++a1SeHaV79+74f//v/8HOzg4nTpxAZmamqVMiMhkWhEREMqmrqwMAeHp6trq/uV2n02Hfvn2tHrNkyRKcOHECb7/9Nnr06NHuXJqamnDgwAEAQHBwcKvHNLfv3btX39avXz8AwE8//WRw7D/+8Q9069YN7u7uAG49qPLuu+/iueeeg7+/f7vz7Gh9+vSBWq0GcOuKLpG1YkFIRCQTV1dXAEB1dXWr+29v/8c//tFi/9GjR/H2229j3LhxePLJJ43KpaqqCo2NjQCAAQMGtHpMc/vx48f1bX5+fujVqxc++eQTbNu2DZcuXcKKFStQXl6Ohx56CPb29gCAF154Aa6urliyZIlReZrC2LFjAeCORTmRNehm6gSIiCxVYGAgAGD//v345ZdfcP/99xvsv33qluaric1EUcRzzz0HnU6HVatWGZ3L7f27ubm1ekxz++3HOjg44M0330R8fDweffRRfbuTkxOWL18OAMjMzMT27dvx4Ycf6vtobGxETU0N3N3dYWtr2+687zb1jBSaP5c73bYnsgYsCImIZBIZGYm+ffvi9OnTmDlzJj7//HP92LstW7YgPT1df+zVq1cN3vvJJ5+gqKgI8+fPh6+vr9G5XLt2Tf//dyrQ7OzsWs3l6aefRt++fbF27Vr8+9//xgMPPICXXnoJgwcPRkNDA+bPn4+RI0ciPj4eoihi4cKFWLFiBa5cuQJHR0fMmzcP6enp7Sru7nR7GwCuX7+O/fv333Ofv+Xo6AgABg/IEFkbFoRERDKxt7fHpk2bEB4ejuLiYvTr1w+DBw9GXV0dTp8+jX79+sHf3x+7d++Gk5OT/n3Ncw56eXkhJSVFslya3bhxw2C72fXr1wHcetjityZNmoRJkya1aE9PT8cvv/yCzz//HF26dMGSJUuwdOlS/PGPf0RUVBSysrLw5ptvwtHREQsWLLjnvIuLi++4r6qqCiqV6p77/K3Lly8DAFxcXIzui8hccQwhEZGMxo4di9LSUjz99NPo06eP/uGM2bNnY//+/fqpTvr06aN/z6uvvora2lq88847BoWiMW6/Tfzb29O/bb/TLeXf+vnnn7F8+XLExsYiKCgIjY2NWL58OQYNGoTc3Fw89dRTyM7OxqBBg7B8+XI0NTUZ/4PI4J///CeAW5NYE1krXiEkIpLZoEGD8Mknn7Rob2pqwsGDBwEAI0aM0Lc3T5I8d+5czJ071+A9zVezMjIysHnzZgC3lsa7G6VSCRsbGzQ2NuLkyZP6W9e3O3nyJADA29u7LT8WXnzxRdjb2+tXWzl27Bjq6+sxc+ZMdOly63pDly5dEBYWhlWrVuEf//gHfHx82tR3R2q+Cjlq1CgTZ0JkOiwIiYhMpLCwEJcvX0bfvn0REBDQYn9ry8g1u3r1aouxfr+nW7duCAgIwN69e/Hdd9+1Ojaved7EoKCgu/a3efNmbNmyBStWrNBfWWsuVn+7Qknz9p3WSjalX3/9Vb8EnzGrwBCZO94yJiIygRs3buiXfpszZ47BRNHl5eUQRbHVV/OYwuYHOERRbHNMjUYDANBqtS1W5fjnP/+Jbdu2AfjvhNl3cv36dSQmJsLX1xcJCQn69uandX/++WeD45u3e/bs2eZcO8LVq1fx5z//GdevX8cDDzxw15+byJKxICQiktHXX39tMNEzAPzyyy9Qq9UoLS3Fgw8+iFdeeUWyeF988QWUSqV+br3bzZ49Gz179sTRo0fx8ssv6+clPH/+PGbOnImmpiZMnjzZ4PZ1a/7617/i559/xnvvvYdu3f57o8nT0xP3338/8vLyUFFRAeDW0nd5eXno06dPm29Fy+3q1avIzs5GUFAQtm/fDkdHR3z++ectVm8hsia8ZUxEJKNvvvkGK1asgJubG5RKJa5du4Zjx45BFEU8+OCD+Oabb/TTvUjh8uXLOHXqVKv7XFxcsHHjRvzxj3/EypUr8dlnn6Ffv344evQoGhoaoFQqsWbNmt/t/5///CeWLVuG6OhoPPzwwwb7BEFAamoq4uPjERgYiMGDB+Onn37C9evXkZKSoh9X2JHy8/P1xfHNmzdRV1eHkydP6othf39/rF+/XpKpfYjMGQtCIiIZqdVq/PrrrygpKcHRo0dhZ2eHwMBATJ8+Hc8//7ykxWBbPPLII9i/fz+WLFmCb7/9FocOHYKnpyemTp2KhQsX3vUJ45dffhmCIOB///d/W93/9NNP49q1a3jnnXdw7Ngx9O/fH//zP/+D2bNny/Hj3NW5c+f0E07b29vD1dUVAQEBGDlyJKZOnYpHHnnEJHkRdTaCeC8DUIiIiIjI4nAMIREREZGVY0FIREREZOVYEBIRERFZORaERERERFaOBSERERGRleO0M9SpiaKIS5cumToNIiIis+Xs7AxBEH73GBaE1KldunQJrq6upk6DiIjIbF24cAEuLi6/ewznIaROjVcIiYiIjNOWK4QsCImIiIisHB8qISIiIrJyLAiJiIiIrBwLQiIiIiIrx4KQiIiIyMqxICQiIiKyciwIiYiIiKwcC0IiIiIiK8eCkIiIiMjKsSAkIiIisnIsCImIiIisHAtCIiIiIivHgpCIiIjIyrEgJCIiIrJyshaEO3fuhCAISE1NNaqf1NRUCIKAnTt3SpLXb2m1WgiCAK1WK0v/UouNjYUgCPrX3//+d1OndFdeXl4GOVdVVZk6JSIiIvoPXiHsIIIgICQkRNI+X3zxRaSkpGDkyJEG7bt378b8+fMxfvx4uLq6QhAExMbGShr73XffRVxcHIYNG4Zu3brdtWB/+eWXkZKSAj8/P0nzICIiIuN1M3UC1H6JiYlQKpUt2tesWYN169bBwcEB/fr1w8WLFyWPPW/ePACAh4cH7rvvPpw5c+Z3j3/55ZcBAFVVVTh48KDk+RC15vjZS6g63wCluwO8ezubOh0iok6LBaEFmjt3Ll555RX84Q9/wL59+zBmzBjJY2zevBkjRoxAnz59MHv2bHzwwQeSxyDLd/7ydVn6rWu4gQXZh7G3slbfFqRSIH2qL9wcbCWP5+5kJ3mfREQd6Z4Kwhs3buCDDz7A5s2b8eOPP+LcuXNwdXXF2LFjsWjRIgwfPrxN/TRf1SorK8Orr76Kr776ChcvXsTQoUOxcOFCTJky5Y7v/fzzz7Fs2TIcPXoUPXr0wOOPP4633noL3bt3lyTP7OxsLF26FEeOHIGLiwumTp2KpUuXws3NrcWxlZWVSE9PxzfffIOzZ89CoVBg4sSJSEtLQ//+/QHcGkc5fvx4AMCuXbsgCIL+/WvXrkVsbCwuXLiAv//97/j6669x/Phx1NTUoGfPnggNDUVKSgoGDhzYpvPa7Le3kOUQEREhewxqm4YbTaZOod1GLNkmS79dBMDZ3gbvzwxAoMoN+yrrkJR9CGHv7IZOlD7ej4snSt9pB3Kw5bUBImt3T98CtbW1SExMxLhx4xAeHg43NzecPHkSX331FfLz87F7924EBga2qa8bN24gNDQUV69exVNPPYX6+nps3LgRarUa69evxxNPPNHiPe+//z7y8/MRGRmJkJAQFBQU4N1338X58+exYcMGo/P84osvsHXrVjz++OMIDQ3Frl278Pe//x179uzBnj17DIrOvXv3YuLEibhy5Qoee+wxDBo0CFVVVdiwYQPy8/OxZ88eDBgwAEqlEikpKfoi8faxfP7+/gCAo0ePIjk5GePHj8fUqVPh6OiIY8eOISMjA1u2bEFpaam+wCT6rQeTC02dQqejE4GlU4ciYpgHACBimAdEiJibUSZLPHP/DKqW8R94RNbungpCNzc3/POf/4Snp6dB+5EjRzB69GgkJSVh69atberr119/xZAhQ/DDDz/AxsYGwK1xZgEBAZg7dy6mTJkCZ2fDMT9bt27FgQMHMHjwYABAeno6/P398dlnn+Htt99G3759jcpzy5Yt2LZtGx555BF929NPP421a9fif//3f7Fo0SIAQGNjI6Kjo6HT6bB//36DByWKi4sREhKCF198EXl5eVAqlUhNTUVaWpr+/39ryJAh+PXXX6FQKAzad+zYgdDQUCxZsgQfffRRm84rEd0SqDK8qj9KpbjDkUREdE8FoZ2dXYsiCwB8fHwwfvx4FBYWorGxUV/g3c1f/vIXg2P/8Ic/4Omnn8b777+P3NxcxMTEGBz/4osv6otBAOjevTtmzJiBtLQ0HDhwQF8QtjfPRx991KAYBIAlS5bg//2//4d169bpC8LNmzejqqoKf/nLX1o8NTt27FhERkYiJycHFy9ehIuLy13Pg6ura6vt48ePh4+PD7Ztk+e2GlkGc75dKeeVtX2VdforhABQctt4QqmZ82dARAS046GS8vJy/PWvf0VxcTHOnDmDxsZGg/01NTXw8PC4w7v/y8bGBqNHj27RPm7cOLz//vsoLy9vURAGBAS0ON7LywsAUF9fb3Se48aNa9F/3759MXDgQBw7dgyXLl2Cs7MzfvjhBwDAsWPHWr3id+bMGeh0Ovz0009tHs+3c+dO/N///R/27t2LmpoaNDX9d1yYra30g+DJcpjz+K8DC0Nl6TdhQykW5R6GCBGjVAqUVNYiOfcIglQKrHqi5feIscz5MyAiAu6xIPz+++8xYcIEAEBYWBi8vb3h5OQEQRCQk5ODgwcP4vr1tj016O7uji5dWk6D2Lt3bwDAhQsXWuxr7Upat263foSbN28anWevXr1azbV37944duwYLl68CGdnZ9TW3rrScPu4xdZcuXLld/c3y8zMxPTp0+Hk5ISJEydCqVTCwcFBP1n2qVOn2tQPkbmR6+nc1TEj8OLGMoMxg+O8e2JF9HAoHPkPLCKi37qngjA9PR3Xr19HcXExgoODDfb98MMP9zS/3Pnz56HT6VoUhWfPngVw59uocuZ57ty5Vtubc2q+/dv837y8PPzxj39sd57NUlNTYW9vjwMHDsDb29tg38aNG43un8jaKBxtsT4+iPMQEhG10T2tVPLzzz9DoVC0KLIaGhpQWlp6T4EbGxv1t15vV1RUBOC/T+C2R3vzbI59u9OnT+Pnn3/GwIED9Q+5BAUFAQD27NnT5py6dOlicBXzt/kOGTKkRTHYHJuI2se7tzMefbA3i0Eioru4p4Kwf//+qKurw5EjR/RtN2/exPz58/Hvf//7noMvWrTIYGzfsWPHsGbNGri6uiIyMvKe+zM2z61bt2L79u0GbQsXLkRjYyOeeuopfVtkZCT69euHv/3tb9i9e3eLfhobG1FcXGzQplAo8K9//euO+Z44cUJ/JRIArl27hjlz5hiMJZRLVVUVBEFoddUTIiIisnz3dMv4hRdewDfffIOxY8di2rRpsLe3x86dO1FdXY2QkJDfXcv2tzw8PFBfXw9/f39ERETgwoUL+Oyzz3Dt2jV89NFHLaac6Yg8IyIiEB4ejscffxz3338/du3ahT179sDPzw/z58/XH2dnZ4cvvvgCkydPxsMPP4xHHnkEvr6+AIB//vOfKCoqgru7O44dO6Z/z4QJE/D5558jKioKw4cPR9euXREREYGhQ4fihRdewAsvvIDhw4cjKioKTU1N2Lp1K0RRhJ+f3z0v9VZcXIyPP/4YAPQFcHFxsX4OxD/84Q94/fXX9cfrdDoA/x2P2RbLli3T/3zNV0qXLVsGrVYLAJg1axbGjh17T3kTERGRiYj36IsvvhADAgJEBwcHsWfPnuK0adPEn3/+WXzqqadEAGJlZaX+2B07dogAxJSUFIM++vfvL/bv3188f/68OGvWLLFXr16inZ2dOHLkSDE3N7dFzJSUFBGAuGPHjhb71q5dKwIQ165d2+48b+8jKytLHDFihGhvby/26tVLfO6558Tz58+3ei7+9a9/iS+++KLo7e0t2tnZiS4uLuKQIUPEWbNmidu3bzc49tdffxWnTZsm9uzZU+zSpYtBzjqdTvz73/8u+vj4iPb29mKfPn3E+Ph48ezZs+LDDz8s/vZjau1naO2c3On18MMPGxyfm5srAhAXLFjQan+tac7rTq/ffh5tzZ2IiIg6niCKogwLOf2+5luTVVVVHR3aIsTGxmLdunWorKyU5Dbv/PnzsXr1apw6dQo9e/Y0PsHfIXXuREREZLx7GkNInYtKpYIgCPj73/9uVD9FRUV45plnZC0Gvby8IAgC1q1bJ1sMIiIiah/OpmqG1Gq1wdW1tk5+fSd79+41MqO7e/nll3Hx4kX9do8ePWSPSURERG3DW8ZEREREVs4kBSERERERdR4cQ2iGtFotBEHQv6Kjo02d0l3FxMQY5Nw8PQ0RERGZHscQmrHIyEj4+/vr50Bs9vPPP2P9+vUoLS3FgQMHcPr0afTv31+yW/Tnz5/Hl19+ic2bN+Pw4cOorq6Gs7MzAgMDkZiYiIkTJ7Z4j0ajwaBBg1BeXo7c3FxJ8iAiIiJpsCA0Y2q1Wj/Z9O2KioqQlpaGrl27YsiQIThz5oykcTMzMzFnzhx4enpiwoQJ8PT0xL/+9S98+eWXKCgowNtvv20wkTdwqyDUaDTQarUsCKnDcC1jIqK2YUFogR566CH9Civdu3eHvb29pP0/8MAD2Lx5MyZPnowuXf476mDhwoUICgpCUlISZs6cib59+0oalyzP+cvXZem3ruEGFmQfxt7KWn1bkEqB9Km+cHOwlSUmALg72cnWNxGRnFgQWqABAwZgwIABsvU/YcKEVtsHDx6M6dOn48MPP8T333+PqKgo2XKg/2q4If9613IZsWSbLP12EQBnexu8PzMAgSo37KusQ1L2IYS9sxs6GR+j+3Fxy+ESnZmDLX8FENEt/DYgSdnY2AC4t3WRyTgPJheaOoVORycCS6cORcQwDwBAxDAPiBAxN6NM1rjm9llULYswdQpE1EnwKWOSzKVLl/DFF1/A3t4e48aNM3U6ZOUCVW4G26NUChNlQkTU+fEyDklm9uzZOHv2LBYvXgx3d3dTp2M1zO025e3kvKK2r7JOf4UQAEpuG08oF3P+LIjIurEgJEkkJSUhIyMDkyZNQlJSkqnTsSrmPA7swMJQWfpN2FCKRbmHIULEKJUCJZW1SM49giCVAqueCJAlJmDenwURWTd+e5HR0tLS8Oabb2LChAnIyspC165dTZ0SmQm5nspdHTMCL24sMxgzOM67J1ZED4fCUb6njImIzBULQjJKWloaUlNTERISgry8PHTv3t3UKRFB4WiL9fFBnIeQiKiNWBBSu6WmpiItLQ0PP/wwtmzZAgcHB1OnRGTAu7czC0EiojbgU8YEAFAqlRAEoc3L26WkpCAtLQ3jxo1jMUhERGTmeIXQAtXU1BgsHdfY2IiamhqDZe60Wq3Be3Q6HYC2zR+o1WqxePFidOvWDaNGjcLbb7/d4piQkBCEhIS0K38iIiLqWCwILdDly5exbt06g7YrV64YtN1eENbX16O6uhrBwcHw8vK6a//NVxGbmpqwfPnyOx7HgpCIiMg8sCC0QEqlEqLY9vW5iouLodPp2jxdTGpqKlJTU9uZHREREXU2HENoxuLi4iAIAqKjo43qp6ioCH5+fggPD5cos5ZiYmIgCALi4uJki0FERETtI4j3cimJOoXy8nLk5OTot319fREVFWW6hNogKysLFRUV+m21Wg1/f3/TJURERER6LAiJiIiIrBxvGRMRERFZORaEZkir1UIQBP3L2DGEHaF5DGHz67fT3hAREZHp8CljMxYZGQl/f3/4+voatP/8889Yv349SktLceDAAZw+fRr9+/dv86TTdyOKIl599VXs27cPP/30E2pra+Hq6oqBAwciPj4eTz75JGxsbAzeo9FoMGjQIJSXlyM3N1eSPIiIiEgaLAjNmFqtNphsullRURHS0tLQtWtXDBkyBGfOnJE07s2bN/Huu+9i5MiRiIiIwH333Ye6ujoUFBRg1qxZyMzMxNdff40uXf57AVqj0UCj0UCr1bIgNANcA5iIyLqwILRADz30EPbs2QM/Pz90794d9vb2kvbfrVs31NfXt+i3qakJYWFhKCwsRH5+PiIiIiSNS7ecv3xdtr7rGm5gQfZh7K2s1bcFqRRIn+oLNwdb2eK6O9nJ1jcREd0dC0ILNGDAAAwYMEDWGK0Vmd26dYNarcaOHTtw4sQJWeNLoeFGk6lTaJcRS7bJ1ncXAXC2t8H7MwMQqHLDvso6JGUfQtg7u6GTcT6CHxdPlK9zmTjY8uuTiCwHv9FIMjqdDgUFBQDQYlxjZ/RgcqGpU+h0dCKwdOpQRAzzAABEDPOACBFzM8pkjWuOn0XVMl4BJyLLwYKQjNK8hF1NTQ22b9+OY8eOITY2Fo888ohpE6N2C1S5GWyPUilMlAkREXUUFoRklLS0NP3/C4KA+fPn48033zRhRm1njrcpAfmvpu2rrNNfIQSAktvGE8rFXD8LIiJLwYKQjCKKInQ6HU6fPo3NmzfjjTfewJ49e/D111/DxcXF1On9LnMdA3ZgYahsfSdsKMWi3MMQIWKUSoGSylok5x5BkEqBVU8EyBbXXD8LIiJLwW9hMlqXLl3g5eWF2bNnw93dHdOmTUN6ejreeustU6dmkeR8Ind1zAi8uLHMYMzgOO+eWBE9HApH+Z4yJiIi02JBSJIKCwsDAOzcudO0iVC7KBxtsT4+iPMQEhFZGRaEJKnTp08DuDUFDZkv797OLASJiKwI1zImAIBSqYQgCG1a3u7YsWM4d+5ci/aGhga8/PLLAIDJkydLnSIRERHJhJdxLFBNTQ3mz5+v325sbERNTY3BMndardbgPTqdDkDbruwVFBTgtddeQ0hICAYMGABXV1dUV1cjPz8f58+fR3BwsL4wJCIios6PBaEFunz5MtatW2fQduXKFYO22wvC+vp6VFdXIzg4GF5eXnftPzQ0FPHx8SguLsa+fftw6dIluLq6wtfXF9HR0Zg1axZvGRMREZkR/ta2QEqlEqLY9nXGiouLodPpkJSU1KbjfX19sWrVqvamR0RERJ0MxxCasbi4OAiCgOjoaKP6KSoqgp+fH8LDwyXKrKWYmBgIgoC4uDjZYhAREVH7COK9XEqiTqG8vBw5OTn6bV9fX0RFRZkuoTbIyspCRUWFflutVsPf3990CREREZEeC0IiIiIiK8dbxkRERERWjgWhGdJqtRAEQf8ydgxhR2geQ9j8+u20N0RERGQ6LAjNWGRkJFJSUlqMH/z555+RmpqKKVOmwNPTE4IgQKlUyppLZmamvtjbuHFji/0ajQYpKSmIjIyUNQ8iIiK6d5x2xoyp1WqDyaabFRUVIS0tDV27dsWQIUNw5swZWfM4d+4cEhIS4OjoiCtXrrR6jEajgUajgVarRW5urqz5EBFR58c10zsXFoQW6KGHHsKePXvg5+eH7t27w97eXtZ4c+bMgaOjI5566iksX75c1lhERNQxzl++Lku/dQ03sCD7MPZW1urbglQKpE/1hZuDreTx3J3sJO/TErEgtEADBgzAgAEDOiRWRkYGsrKyUFhYiO+//75DYhIRmYOGG02mTsEoI5Zsk6XfLgLgbG+D92cGIFDlhn2VdUjKPoSwd3ZDJ8O8Jz8unih9pzJxsDVdWcaCkNrtzJkzeOGFF/D0008jLCyMBSER0W0eTC40dQqdkk4Elk4diohhHgCAiGEeECFibkaZLPHM6XOoWhZhsth8qITa7bnnnoO9vT1vExMR0T0JVLkZbI9SKUyUCTXjFUJql08//RRfffUVcnNz0aNHD1OnQ0TU6ZjTrcrWyHllbV9lnf4KIQCU3DaeUGrm/jl0FBaEdM9Onz6NxMREREdHY8qUKaZOh4ioUzLleDApHFgYKku/CRtKsSj3MESIGKVSoKSyFsm5RxCkUmDVEwGSxzP3z6Gj8CzRPUtISEDXrl3x7rvvmjoVIiKSiVxP566OGYEXN5YZjBkc590TK6KHQ+Eo/VPG1DYsCOmelZeXo6amBvfdd1+r+2fMmIEZM2bgnXfeQWJiYscmR0REnZrC0Rbr44M4D2Enw4KQ7ll0dDRqampatJeWlqKsrAzjx4/HgAED4Ovra4LsiIjIHHj3dmYh2ImwICQAgFKpxKlTp1BZWXnXZe6WLVvWantqairKysrw7LPPmsX6ykRERHQLC0ILVFNTg/nz5+u3GxsbUVNTY7DMnVarNXiPTqcDAHTrxj8SRERE1oa//S3Q5cuXsW7dOoO2K1euGLTdXhDW19ejuroawcHB8PLy6qg0iYiIqJPgxNQWSKlUQhTF333drri4GDqdDklJSUbFTU1NhSiKvF1MRERkZlgQmrG4uDgIgmB0AVZUVAQ/Pz+Eh4dLlFlLMTExEAQBcXFxssUgIiKi9hHE314uok6vvLwcOTk5+m1fX19ERUWZLqE2yMrKQkVFhX5brVbD39/fdAkRERGRHgtCIiIiIivHW8ZEREREVo4FoRnSarUQBEH/MoeHOJrHEDa/fjvtDREREZkOp50xY5GRkfD392+xIsjPP/+M9evXo7S0FAcOHMDp06fRv39/VFVVSRY7JCQEu3btanXfxIkTUVBQYNCm0WgwaNAglJeXIzc3V7I8iIiIyHgsCM2YWq02mGy6WVFREdLS0tC1a1cMGTIEZ86ckS2HlJSUFm2DBg1q0abRaKDRaKDValkQksXhmqxEZO5YEFqghx56CHv27IGfnx+6d+8Oe3t72WKlpqbK1jeRFM5fvi5b33UNN7Ag+zD2Vtbq24JUCqRP9YWbg60sMd2d7GTpl4isGwtCCzRgwAAMGDDA1GmQBWi40WTqFIw2Ysk22fruIgDO9jZ4f2YAAlVu2FdZh6TsQwh7Zzd0Ms3f8OPiifJ0LCMHW/6qIers+LeUjLJx40ZUVlbC0dERgYGBGDNmjKlTIgk9mFxo6hQ6NZ0ILJ06FBHDPAAAEcM8IELE3Iwy2WKa42dStSzC1CkQ0V2wICSjzJgxw2A7MDAQmzZtgkqlMlFGRB0rUOVmsD1KpTBRJkRE7ceCkNpFrVbj9ddfh7+/P1xcXHD8+HH87W9/w6efforQ0FAcOnQIDg4Opk6TjGSOtyd/S+4ravsq6/RXCAGg5LbxhHKwhM+EiDofFoTULomJiQbbfn5+WLduHZqampCRkYG1a9fi+eefN01yJBlLGPt1YGGobH0nbCjFotzDECFilEqBkspaJOceQZBKgVVPBMgS0xI+EyLqfPjNQpKKj49HRkYGvvvuOxaE1CnI+VTu6pgReHFjmcGYwXHePbEiejgUjvI8ZUxEJAcWhCSpnj17AgAaGhpMnAmR/BSOtlgfH8R5CInI7LEgJEnt3bsXAKBUKk2bCFEH8u7tzEKQiMwa1zImALcKOEEQ2rS83cmTJ/Hrr7+2aD969CgWLFgAAGaxvjIRERHdwiuEFqimpgbz58/Xbzc2NqKmpsZgmTutVmvwHp1OBwDo1u3ufyR2796NZ555BuPHj8fAgQPh7OyM48ePY8uWLWhsbERycjJGjx4tyc9CRERE8mNBaIEuX76MdevWGbRduXLFoO32grC+vh7V1dUIDg6Gl5fXXfsPCAjAjBkzsH//fpSUlODKlStwd3fH5MmT8fzzzyMsLEyyn4WIiIjkx4LQAimVSohi29fNKi4uhk6nQ1JSUpuOHzZsGD799NP2pkdERESdDMcQmrG4uDgIgmD0eL2ioiL4+fkhPDxcosxaiomJgSAIiIuLky0GERERtY8g3sulJOoUysvLkZOTo9/29fVFVFSU6RJqg6ysLFRUVOi31Wo1/P39TZcQERER6bEgJCIiIrJyvGVMREREZOVYEJohrVYLQRD0L3OY8695DGHz67fT3hAREZHpsCA0Y5GRkUhJSWkxfvDnn39GamoqpkyZAk9PTwiCINvKITt37kRkZCR69eoFOzs73H///Zg6dSoOHjxocJxGo0FKSgoiIyNlyYOIiIjaj9POmDG1Wm0w2XSzoqIipKWloWvXrhgyZAjOnDkjS/z09HQsXLgQffv2hVqtRs+ePXH27Fl89913OHToEPz8/PTHajQaaDQaaLVa5ObmypIPWR+uIUxEJA0WhBbooYcewp49e+Dn54fu3bvD3t5e8hhfffUVFi5cCLVajYyMDHTv3t1gf1NTk+Qxybycv3xdtr7rGm5gQfZh7K2s1bcFqRRIn+oLNwdbyeO5O9lJ3icRUWfCgtACDRgwAAMGDJA1xuuvvw5nZ2dotdoWxSDQtiXwqH0abphHsT1iyTbZ+u4iAM72Nnh/ZgACVW7YV1mHpOxDCHtnN3QyzJvw4+KJ0nfawRxs+XeSiO6M3xB0zyoqKnD06FFoNBo4OTkhPz8fFRUVcHBwwEMPPWRwq5ik92ByoalTMDmdCCydOhQRwzwAABHDPCBCxNyMMlniWcI5r1oWYeoUiKgTY0FI92z//v0AAHd3d4wdOxY//PCDwf4nnngCa9asga2t9LfuiJoFqtwMtkepFCbKhIjI/LEgpHt27tw5AMCaNWugUqnw7bffIjAwEMePH8fzzz+PDRs2wNPTE2+99ZaJM7VM5nL7Uu6ravsq6/RXCAGg5LbxhFIzl3NORNReLAjpnul0Ov1/P//8cwwfPhwAMHz4cOTk5GDQoEF47733sHjxYtjZcTC+1MxlLNiBhaGy9Z2woRSLcg9DhIhRKgVKKmuRnHsEQSoFVj0RIHk8cznnRETtxW85umeurq4AAC8vL30x2KxXr14ICgrCtm3bcPToUa5XbMXkfDJ3dcwIvLixzGDM4DjvnlgRPRwKRw5VICK6VywI6Z4NHjwYANCjR49W9ze3X716tYMyImujcLTF+vggzkNIRCQRFoR0z0aPHo3u3bvj5MmTuHbtWot5Do8ePQoAsq2OQtTMu7czC0EiIglw6ToCcKt4EwQBVVVVdz3WyckJf/7zn3HlyhUsWbLEYN/69etx5MgRjB07Fh4eHnfogYiIiDoTXiG0QDU1NZg/f75+u7GxETU1NQbL3Gm1WoP3ND8o0tYJpZcuXYqdO3ciPT0dxcXFGDlyJI4fP468vDy4ubnhgw8+MPrnICIioo7BgtACXb58GevWrTNou3LlikHb7QVhfX09qqurERwcDC8vrzbFcHd3x549e5CWlobs7Gx8//33UCgUiImJQWpqquwrpRAREZF0WBBaIKVSCVFs+/pdxcXF0Ol0SEpKuqc4CoUCK1aswIoVK+41RSIiIupEOIbQjMXFxUEQBERHRxvVT1FREfz8/BAeHi5RZi3FxMRAEATExcXJFoOIiIjaRxDv5VISdQrl5eXIycnRb/v6+iIqKsp0CbVBVlYWKioq9NtqtZpzFBIREXUSLAiJiIiIrBxvGRMRERFZORaEZkir1UIQBP3L2DGEHaF5DGHz67fT3hAREZHpsCA0Y5GRkUhJSTEYPyiKIvLz8zFnzhwMGzYMrq6ucHBwgJ+fH5YuXYpr165JEjsvLw8vvPACgoOD4ejoCEEQkJqaesfjNRoNUlJSEBkZKUl8IiIikg6nnTFjarXaYLJpALh+/TrCw8NhZ2eHkJAQTJw4EdeuXUNhYSEWLFiAnJwc7Nq1C927dzcq9vLly7Fr1y64uLigb9++OHHixO8er9FooNFooNVqkZuba1RsIuqcuLY0kfliQWhhunbtivT0dCQkJKBHjx769sbGRvzpT39CXl4e3nvvPbzyyitGxfnLX/6CPn36YNCgQdi0aRNmzJhhZOZEJLfzl6/L0m9dww0syD6MvZW1+rYglQLpU33h5mArS0wAcHeyk61vImvDgtDC2NjYtDrBtI2NDd544w3k5eVh165dRheE48aNM+r9RJag4UaTqVO4JyOWbJOl3y4C4Gxvg/dnBiBQ5YZ9lXVIyj6EsHd2QyfjPBY/Lp4oX+cdxMGWv4apc+CfRCtiY2MDoO3rFRPR73swudDUKXQKOhFYOnUoIoZ5AAAihnlAhIi5GWWyxrWE81+1LMLUKRAB4EMlVmXNmjUAgLCwMBNnQkSWJlDlZrA9SqUwUSZE1B68VGQlCgoK8MEHH2DIkCGIj483dTpEFsHcblnKeUVtX2Wd/gohAJTcNp5QLuZ2/ok6MxaEVmD//v2YPn06XF1dkZmZCTs7DsQmkoK5jf86sDBUln4TNpRiUe5hiBAxSqVASWUtknOPIEilwKonAmSJCZjf+SfqzPi3ycKVlZUhLCwMgiCgsLAQPj4+pk6JiExErqdyV8eMwIsbywzGDI7z7okV0cOhcJTvKWMikg4LQgtWWlqKRx99FDdv3sQ333yDwMBAU6dERBZI4WiL9fFBnIeQyIyxILRQpaWlCA0NRVNTEwoLCxEUFGTqlIjIwnn3dmYhSGSm+JSxBWouBhsbG5Gfn48xY8bc9T1KpRKCIKCqqkr+BImIiKhT4RVCC1NbW4vQ0FDU1dVh0qRJ2Lp1K7Zu3WpwTI8ePZCYmGjQptPpALR9jsKcnBzk5OQAACorK/VtzQXl2LFjMWvWrPb/IERERNRhWBBamIsXL6Kurg7AralmCgoKWhzTv39/g4Kwvr4e1dXVCA4OhpeXV5vilJeXY926dQZtBw8exMGDB/XbLAiJiIjMAwtCC6NUKiGK97ZWVHFxMXQ6XatL3t1JamoqUlNT7zE7IiIi6ow4htCMxcXFQRAEREdHG9VPUVER/Pz8EB4eLlFmLcXExEAQBMTFxckWg4iIiNpHEO/1chKZXHl5uX78HgD4+voiKirKdAm1QVZWFioqKvTbarUa/v7+pkuIiIiI9FgQEhEREVk53jImIiIisnIsCM2QVquFIAj6l7FjCDtC8xjC5pdWqzV1SkRERPQffMrYjEVGRsLf3x++vr76NlEUUVBQgK+++grfffcdTp06hcbGRnh7e2P69Ol4+eWXYW9vb1Tc8+fP48svv8TmzZtx+PBhVFdXw9nZGYGBgUhMTMTEiRNbvEej0WDQoEEoLy9Hbm6uUfGJiIhIWiwIzZharUZsbKxB2/Xr1xEeHg47OzuEhIRg4sSJuHbtGgoLC7FgwQLk5ORg165d6N69e7vjZmZmYs6cOfD09MSECRPg6emJf/3rX/jyyy9RUFCAt99+G/Pnzzd4j0ajgUajgVarZUFIsuA6ukRE7ceC0MJ07doV6enpSEhIQI8ePfTtjY2N+NOf/oS8vDy89957eOWVV9od44EHHsDmzZsxefJkdOny31EHCxcuRFBQEJKSkjBz5kz07dvXmB+FLMz5y9dl6beu4QYWZB/G3spafVuQSoH0qb5wc7CVPJ67k53kfRIRmRoLQgtjY2PT6gTTNjY2eOONN5CXl4ddu3YZVRBOmDCh1fbBgwdj+vTp+PDDD/H99993+qlwzEnDjSZTp2C0EUu2ydJvFwFwtrfB+zMDEKhyw77KOiRlH0LYO7uhk2EOhR8XtxwSYQ4cbPl1T0R3xm8IK2JjYwOg7esVd9YY1ujB5EJTp9Bp6URg6dShiBjmAQCIGOYBESLmZpTJEs9cP4uqZRGmToGIOjE+ZWxF1qxZAwAICwuTpf9Lly7hiy++gL29PcaNGydLDKLWBKrcDLZHqRQmyoSIyDzxMo6VKCgowAcffIAhQ4YgPj5elhizZ8/G2bNnsXjxYri7u8sSw1qZ623K28l5ZW1fZZ3+CiEAlNw2nlBqlvBZEBH9FgtCK7B//35Mnz4drq6uyMzMhJ2d9IPik5KSkJGRgUmTJrU6hpGMYwnjvw4sDJWl34QNpViUexgiRIxSKVBSWYvk3CMIUimw6okAyeNZwmdBRPRb/GazcGVlZQgLC4MgCCgsLISPj4/kMdLS0vDmm29iwoQJyMrKQteuXSWPQeZPrqdzV8eMwIsbywzGDI7z7okV0cOhcJT+KWMiIkvEgtCClZaW4tFHH8XNmzfxzTffIDAwUPIYaWlpSE1NRUhICPLy8oya35CoPRSOtlgfH8R5CImIjMCC0EKVlpYiNDQUTU1NKCwsRFBQkOQxUlNTkZaWhocffhhbtmyBg4OD5DGI2sq7tzMLQSKiduJTxhaouRhsbGxEfn4+xowZc9f3KJVKCIKAqqqqNsVISUlBWloaxo0bx2KQiIjIzPEKoYWpra1FaGgo6urqMGnSJGzduhVbt241OKZHjx5ITEw0aNPpdADaNn+gVqvF4sWL0a1bN4waNQpvv/12i2NCQkIQEhLS7p+DiIiIOg4LQgtz8eJF1NXVAbg11UxBQUGLY/r3729QENbX16O6uhrBwcHw8vK6a4zmq4hNTU1Yvnz5HY9jQUhERGQeWBBaGKVSCVG8t/W6iouLodPp2jxdTGpqKlJTU9uRHREREXVGHENoxuLi4iAIAqKjo43qp6ioCH5+fggPD5cos5ZiYmIgCALi4uJki0FERETtI4j3ejmJTK68vBw5OTn6bV9fX0RFRZkuoTbIyspCRUWFflutVsPf3990CREREZEeC0IiIiIiK8dbxkRERERWjgWhGdJqtRAEQf8ydgxhR2geQ9j80mq1pk6JiIiI/oNPGZuxyMhI+Pv7w9fXV98miiIKCgrw1Vdf4bvvvsOpU6fQ2NgIb29vTJ8+HS+//DLs7e2NiiuKIl599VXs27cPP/30E2pra+Hq6oqBAwciPj4eTz75JGxsbAzeo9FoMGjQIJSXlyM3N9eo+ERERCQtjiE0Q1qtFnFxcVi7di1iY2MN9l27dg3du3eHnZ0dQkJCMHToUFy7dg2FhYU4fvw4AgMDsWvXLqPWHG5qaoKTkxNGjhyJIUOG4L777kNdXR0KCgpQVVWFiRMn4uuvv0aXLi0vQP9e7kREBK7LTSbBK4QWpmvXrkhPT0dCQgJ69Oihb29sbMSf/vQn5OXl4b333sMrr7zS7hjdunVDfX19iyuNTU1NCAsLQ2FhIfLz8xEREdHuGEREndX5y9dl6beu4QYWZB/G3spafVuQSoH0qb5wc7CVJaa7k50s/ZL5YUFoYWxsbFqdYNrGxgZvvPEG8vLysGvXLqMKQgCt3nbu1q0b1Go1duzYgRMnThjVPxFZtoYbTaZOod1GLNkmS79dBMDZ3gbvzwxAoMoN+yrrkJR9CGHv7IZOpnt5Py6eKE/HMnOwZfkiNZ5RK9I8rq8t6xW3h06n0y+Vd/u4RiKi33owudDUKXQ6OhFYOnUoIoZ5AAAihnlAhIi5GWWyxTTXz6FqGe9ASY0FoRVZs2YNACAsLEyyPpuXsKupqcH27dtx7NgxxMbG4pFHHpEsBhGRtQhUuRlsj1IpTJQJWRsWhFaioKAAH3zwAYYMGYL4+HjJ+k1LS9P/vyAImD9/Pt58803J+iciy2SutyoBea+q7aus018hBICS28YTysGcPweSFgtCK7B//35Mnz4drq6uyMzMhJ2ddIOIRVGETqfD6dOnsXnzZrzxxhvYs2cPvv76a7i4uEgWh4gsizmPATuwMFSWfhM2lGJR7mGIEDFKpUBJZS2Sc48gSKXAqicCZIlpzp8DSYt/EixcWVkZwsLCIAgCCgsL4ePjI3mMLl26wMvLC7Nnz4a7uzumTZuG9PR0vPXWW5LHIiIyNbmezF0dMwIvbiwzGDM4zrsnVkQPh8JRnqeMiZqxILRgpaWlePTRR3Hz5k188803CAwMlD1m8/jEnTt3yh6LiMiSKBxtsT4+iPMQkkmwILRQpaWlCA0NRVNTEwoLCxEUFNQhcU+fPg1AvieZiYgsnXdvZxaC1OG4lrEFai4GGxsbkZ+fjzFjxtz1PUqlEoIgoKqq6q7HHjt2DOfOnWvR3tDQgJdffhkAMHny5HvOm4iIiEyDl3EsTG1tLUJDQ1FXV4dJkyZh69at2Lp1q8ExPXr0QGJiokGbTqcD0LYrewUFBXjttdcQEhKCAQMGwNXVFdXV1cjPz8f58+cRHBysLwyJiIio82NBaGEuXryIuro6ALcKt+aJom/Xv39/g4Kwvr4e1dXVCA4OhpeX111jhIaGIj4+HsXFxdi3bx8uXboEV1dX+Pr6Ijo6GrNmzeItYyIiIjPC39oWRqlUQhTvbY2j4uJi6HS6Vpe8a42vry9WrVrVnvSIiIioE+IYQjMWFxcHQRAQHR1tVD9FRUXw8/NDeHi4RJm1FBMTA0EQEBcXJ1sMIiIiah9BvNfLSWRy5eXlyMnJ0W/7+voiKirKdAm1QVZWFioqKvTbarUa/v7+pkuIiIiI9FgQEhEREVk53jImIiIisnIsCM2QVquFIAj6l7FjCDtC8xjC5pdWqzV1SkRERPQfLAjNWGRkJFJSUgzGD4qiiPz8fMyZMwfDhg2Dq6srHBwc4Ofnh6VLl+LatWuy5JKZmakv9jZu3Nhiv0ajQUpKCiIjI2WJT0RERO3HaWfMmFqtRmxsrEHb9evXER4eDjs7O4SEhGDixIm4du0aCgsLsWDBAuTk5GDXrl3o3r27ZHmcO3cOCQkJcHR0xJUrV1o9RqPRQKPRQKvVIjc3V7LY1PlwHVYiIvPDgtDCdO3aFenp6UhISECPHj307Y2NjfjTn/6EvLw8vPfee3jllVckizlnzhw4OjriqaeewvLlyyXrl6R1/vJ1Wfuva7iBBdmHsbeyVt8WpFIgfaov3BxsJY3l7mQnaX9ERNaOBaGFsbGxaXWCaRsbG7zxxhvIy8vDrl27JCsIMzIykJWVhcLCQnz//feS9NlZNdxoMnUKRhmxZJus/XcRAGd7G7w/MwCBKjfsq6xDUvYhhL2zGzqJ5zL4cfFEaTuUkYMtv2aJqPPjN5UVsbGxAdC29Yrb4syZM3jhhRfw9NNPIywszOILwgeTC02dQqemE4GlU4ciYpgHACBimAdEiJibUSZ5LHP6LKqWRZg6BSKiu+JDJVZkzZo1AICwsDBJ+nvuuedgb2/P28SkF6hyM9gepVKYKBMiIroXvEJoJQoKCvDBBx9gyJAhiI+PN7q/Tz/9FF999RVyc3MNxipaMnO6Tdmajriqtq+yTn+FEABKbhtPKCVz/yyIiDobFoRWYP/+/Zg+fTpcXV2RmZkJOzvjBuSfPn0aiYmJiI6OxpQpUyTKsvMz97FgBxaGytp/woZSLMo9DBEiRqkUKKmsRXLuEQSpFFj1RICkscz9syAi6mz4rWrhysrKEBYWBkEQUFhYCB8fH6P7TEhIQNeuXfHuu+9KkCF1FLmfzF0dMwIvbiwzGDM4zrsnVkQPh8JR2qeMiYhIWiwILVhpaSkeffRR3Lx5E9988w0CAwMl6be8vBw1NTW47777Wt0/Y8YMzJgxA++88w4SExMliUmdn8LRFuvjgzgPIRGRGWJBaKFKS0sRGhqKpqYmFBYWIigoSLK+o6OjUVNT02rMsrIyjB8/HgMGDICvr69kMcl8ePd2ZiFIRGRmWBBaoOZisLGxEQUFBRgzZsxd36NUKnHq1ClUVlZCqVT+7rHLli1rtT01NRVlZWV49tlnzWJ9ZSIiIrqFBaGFqa2tRWhoKOrq6jBp0iRs3boVW7duNTimR48eLW7l6nQ6ANLNUUhERETmg7/9LczFixdRV1cH4NZUMwUFBS2O6d+/v0FBWF9fj+rqagQHB8PLy6ujUiUiIqJOghNTWxilUglRFH/3VVVVZfCe4uJi6HS6Vpe8uxepqakQRZG3i4mIiMwMC0IzFhcXB0EQjC7AioqK4Ofnh/DwcIkyaykmJgaCICAuLk62GERERNQ+giiKEi87T3IrLy9HTk6OftvX1xdRUVGmS6gNsrKyUFFRod9Wq9Xw9/c3XUJERESkx4KQiIiIyMrxljERERGRlWNBaIa0Wi0EQdC/zOEhjuYxhM0vrVZr6pSIiIjoP1gQmrHIyEikpKQYjB8URRH5+fmYM2cOhg0bBldXVzg4OMDPzw9Lly7FtWvXJIkdEhJiUODd/po0aVKL4zUaDVJSUhAZGSlJfCIiIpIO5yE0Y2q1GrGxsQZt169fR3h4OOzs7BASEoKJEyfi2rVrKCwsxIIFC5CTk4Ndu3ahe/fukuSQkpLSom3QoEEt2jQaDTQaDbRaLXJzcyWJTabFNYuJiCwHC0IL07VrV6SnpyMhIQE9evTQtzc2NuJPf/oT8vLy8N577+GVV16RJF5qaqok/ZD0zl++Lku/dQ03sCD7MPZW1urbglQKpE/1hZuDrSwx3Z3sZOmXiIhuYUFoYWxsbFqdYNrGxgZvvPEG8vLysGvXLskKQmvUcKPJ1Cm0yYgl22Tpt4sAONvb4P2ZAQhUuWFfZR2Ssg8h7J3d0Mk0Z8GPiyfK07HMHGz5FUtE5oHfVlbExsYGgLTrFW/cuBGVlZVwdHREYGAgxowZI1nfndWDyYWmTsGkdCKwdOpQRAzzAABEDPOACBFzM8pki2mu57xqWYSpUyAiahMWhFZkzZo1AICwsDDJ+pwxY4bBdmBgIDZt2gSVSiVZDOp8AlVuBtujVAoTZUJERFJgQWglCgoK8MEHH2DIkCGIj483uj+1Wo3XX38d/v7+cHFxwfHjx/G3v/0Nn376KUJDQ3Ho0CE4ODhIkHnnYy63L+W8qravsk5/hRAASm4bTygHcznnRETmigWhFdi/fz+mT58OV1dXZGZmws7O+AH6iYmJBtt+fn5Yt24dmpqakJGRgbVr1+L55583Ok5nZC7jwg4sDJWl34QNpViUexgiRIxSKVBSWYvk3CMIUimw6okAWWKayzknIjJX/Ja1cGVlZQgLC4MgCCgsLISPj4+s8eLj45GRkYHvvvvOYgtCcyHXk7mrY0bgxY1lBmMGx3n3xIro4VA4yvOUMRERyYsFoQUrLS3Fo48+ips3b+Kbb75BYGCg7DF79uwJAGhoaJA9FpmGwtEW6+ODOA8hEZEFYUFooUpLSxEaGoqmpiYUFhYiKCioQ+Lu3bsXAKBUKjskHpmOd29nFoJERBaCS9dZoOZisLGxEfn5+W2aCkapVEIQBFRVVd312JMnT+LXX39t0X706FEsWLAAAMxifWUiIiK6hVcILUxtbS1CQ0NRV1eHSZMmYevWrdi6davBMT169GjxUIhOpwPQtjkKd+/ejWeeeQbjx4/HwIED4ezsjOPHj2PLli1obGxEcnIyRo8eLdnPRERERPJiQWhhLl68iLq6OgC3ppopKChocUz//v0NCsL6+npUV1cjODgYXl5ed40REBCAGTNmYP/+/SgpKcGVK1fg7u6OyZMn4/nnn5d0nkMiIiKSHwtCC6NUKiGK97Z+WHFxMXQ6XatL3rVm2LBh+PTTT9uTHhEREXVCHENoxuLi4iAIgtHj9YqKiuDn54fw8HCJMmspJiYGgiAgLi5OthhERETUPoJ4r5eTyOTKy8uRk5Oj3/b19UVUVJTpEmqDrKwsVFRU6LfVajX8/f1NlxARERHpsSAkIiIisnK8ZUxERERk5VgQmiGtVgtBEPQvc5jzr3kMYfNLq9WaOiUiIiL6DxaEZiwyMhIpKSkG4wdFUUR+fj7mzJmDYcOGwdXVFQ4ODvDz88PSpUtx7do1SXPYuXMnIiMj0atXL9jZ2eH+++/H1KlTcfDgQYPjNBoNUlJSEBkZKWl8IiIiMh6nnTFjarUasbGxBm3Xr19HeHg47OzsEBISgokTJ+LatWsoLCzEggULkJOTg127dqF79+5Gx09PT8fChQvRt29fqNVq9OzZE2fPnsV3332HQ4cOwc/PT3+sRqOBRqOBVqtFbm6u0bGJOgrXbCYia8CC0MJ07doV6enpSEhIQI8ePfTtjY2N+NOf/oS8vDy89957eOWVV4yK89VXX2HhwoVQq9XIyMhoUWA2NTUZ1T9RW5y/fF22vusabmBB9mHsrazVtwWpFEif6gs3B1tZYro72cnSLxHR3bAgtDA2NjatTjBtY2ODN954A3l5edi1a5fRBeHrr78OZ2dnaLXaVq82tmUJPOocGm6Yb/E+Ysk22fruIgDO9jZ4f2YAAlVu2FdZh6TsQwh7Zzd0Ms3N8OPiifJ0LDMHW/59JzJ3/FtsRWxsbAAYX6xVVFTg6NGj0Gg0cHJyQn5+PioqKuDg4ICHHnrI4FYxdX4PJheaOoVOSScCS6cORcQwDwBAxDAPiBAxN6NMtpjm+llULYswdQpEZCQWhFZkzZo1AGD0WsP79+8HALi7u2Ps2LH44YcfDPY/8cQTWLNmDWxt5bmtRtRRAlVuBtujVAoTZUJEJC8WhFaioKAAH3zwAYYMGYL4+Hij+jp37hyAWwWmSqXCt99+i8DAQBw/fhzPP/88NmzYAE9PT7z11ltSpE4yM9fblID8V9T2VdbprxACQMlt4wnlYM6fBRGZNxaEVmD//v2YPn06XF1dkZmZCTs74wau63Q6/X8///xzDB8+HAAwfPhw5OTkYNCgQXjvvfewePFio2OR/Mx5/NeBhaGy9Z2woRSLcg9DhIhRKgVKKmuRnHsEQSoFVj0RIEtMc/4siMi88dvHwpWVlSEsLAyCIKCwsBA+Pj5G9+nq6goA8PLy0heDzXr16oWgoCBs27YNR48e5XrFJCs5n8pdHTMCL24sMxgzOM67J1ZED4fCkcMhiMiysCC0YKWlpXj00Udx8+ZNfPPNNwgMDJSk38GDBwOAwbQ2t2tuv3r1qiTxiExB4WiL9fFBnIeQiKwCC0ILVVpaitDQUDQ1NaGwsBBBQUGS9T169Gh0794dJ0+exLVr12Bvb2+w/+jRowAApVIpWUwiU/Hu7cxCkIgsHpeus0DNxWBjYyPy8/MxZsyYu75HqVRCEARUVVXd9VgnJyf8+c9/xpUrV7BkyRKDfevXr8eRI0cwduxYeHh43KEHIiIi6kx4hdDC1NbWIjQ0FHV1dZg0aRK2bt2KrVu3GhzTo0cPJCYmGrQ1PyjS1jkKly5dip07dyI9PR3FxcUYOXIkjh8/jry8PLi5ueGDDz6Q5OchIiIi+bEgtDAXL15EXV0dgFtTzRQUFLQ4pn///gYFYX19PaqrqxEcHAwvL682xXF3d8eePXuQlpaG7OxsfP/991AoFIiJiUFqaioGDBggyc9DRERE8mNBaGGUSiVE8d7W1SouLoZOp2t1ybvfo1AosGLFCqxYseKe3kdERESdC8cQmrG4uDgIgoDo6Gij+ikqKoKfnx/Cw8MlyqylmJgYCIKAuLg42WIQERFR+wjivV5OIpMrLy9HTk6OftvX1xdRUVGmS6gNsrKyUFFRod9Wq9Wco5CIiKiTYEFIREREZOV4y5iIiIjIyrEgNENarRaCIOhfxo4h7AjNYwibX1qt1tQpERER0X+wIDRjkZGRSElJaTF+MCsrC1FRUfD29oaLiwucnJzg4+ODxMREVFdXSxI7Ly8PL7zwAoKDg+Ho6AhBEJCamnrH4zUaDVJSUhAZGSlJfCIiIpIOp50xY2q1GrGxsS3as7OzcfDgQQQGBupXCykvL8fKlSuxbt06FBcXw8fHx6jYy5cvx65du+Di4oK+ffvixIkTv3u8RqOBRqOBVqtFbm6uUbGJOgrXMSYia8GC0AJ99NFHLdYXBoBPPvkEs2bNQmpqKjIzM42K8Ze//AV9+vTBoEGDsGnTJsyYMcOo/oju1fnL12Xru67hBhZkH8beylp9W5BKgfSpvnBzsJUtrruTnWx9ExH9HhaEFqi1YhAAHn/8ccyaNeuuV/PaYty4cUb3QZ1Dw40mU6fQLiOWbJOt7y4C4Gxvg/dnBiBQ5YZ9lXVIyj6EsHd2QyfjvAw/Lp4oX+cycrDlrxIic8e/xVZky5YtAG7NW0jU7MHkQlOn0OnoRGDp1KGIGHZryEXEMA+IEDE3o0zWuOb6WVQtizB1CkRkJBaEFiwnJwfl5eVoaGjAkSNHUFhYCJVKhcWLF5s6NaJOL1DlZrA9SqUwUSZERPJjQWjBcnJysG7dOv32yJEjsXHjRqhUKhNmRZ2Nud6mlPtq2r7KOv0VQgAouW08oVzM9bMgIvPHgtCCabVaaLVaXLhwAWVlZViwYAFGjBiBrKwsTJgwwdTpUSdhruO/DiwMla3vhA2lWJR7GCJEjFIpUFJZi+TcIwhSKbDqiQDZ4prrZ0FE5o/fPlbA1dUVISEhyM/Px+DBg/Hkk0+isrISNjY2pk6NqN3kfCJ3dcwIvLixzGDM4DjvnlgRPRwKR/meMiYiMhUWhFbExcUFo0ePRk5ODk6cOIEhQ4aYOiWiTknhaIv18UGch5CIrAYLQitz+vRpAEC3bvzoie7Gu7czC0Eisgpcus7CXL9+HT/88EOr+9auXYuSkhIMGjQI3t7eBvuUSiUEQUBVVVUHZElERESdCS8TWZirV69izJgx8PX1hb+/Pzw9PXHhwgWUlJSgtLQUTk5OWLt2bYv36XQ6AG2/cpiTk4OcnBwAQGVlpb6tuaAcO3YsZs2aZfwPRERERLJjQWhhHB0dkZaWhh07dmD79u2oqamBjY0NlEolEhMT8dJLL6Ffv34G76mvr0d1dTWCg4Ph5eXVpjjl5eUGU9oAwMGDB3Hw4EH9NgtCIiIi88CC0MLY2NggOTkZycnJbX5PcXExdDodkpKS2vye1NRUpKamtiNDIiIi6mw4htCMxcXFQRAEREdHG9VPUVER/Pz8EB4eLlFmLcXExEAQBMTFxckWg4iIiNpHEEVRxqXaSQ7l5eX68XvArbWJo6KiTJdQG2RlZaGiokK/rVar4e/vb7qEiIiISI8FIREREZGV4y1jIiIiIivHgtAMabVaCIKgfxk7hrAjNI8hbH5ptVpTp0RERET/wYLQjEVGRiIlJaXF+MGsrCxERUXB29sbLi4ucHJygo+PDxITE1FdXW103PPnz+PDDz/ElClTMGDAANjZ2aFnz56YPHkyCgsLW32PRqNBSkoKIiMjjY5PRERE0uK0M2ZMrVYjNja2RXt2djYOHjyIwMBAeHh4ALj1IMrKlSuxbt06FBcXw8fHp91xMzMzMWfOHHh6emLChAnw9PTEv/71L3z55ZcoKCjA22+/jfnz5xu8R6PRQKPRQKvVIjc3t92xiei/uNYyEUmFBaEF+uijj2Bvb9+i/ZNPPsGsWbOQmpqKzMzMdvf/wAMPYPPmzZg8eTK6dPnvReaFCxciKCgISUlJmDlzJvr27dvuGETm7vzl67L1XddwAwuyD2NvZa2+LUilQPpUX7g52Eoez93JTvI+iahzYUFogVorBgHg8ccfx6xZs3DixAmj+p8wYUKr7YMHD8b06dPx4Ycf4vvvv+/0U+GQeWm40WTqFO7JiCXbZOu7iwA429vg/ZkBCFS5YV9lHZKyDyHsnd3QyTBvxI+LJ0rfaQdysOWvOqK74d8SK7JlyxYAt+YtlIuNjQ2Atq+JTNRWDya3Pj7VGulEYOnUoYgYdmtISMQwD4gQMTejTJZ45n7uq5ZFmDoFok6Pv7UtWE5ODsrLy9HQ0IAjR46gsLAQKpUKixcvliXepUuX8MUXX8De3h7jxo2TJQYR3RKocjPYHqVSmCgTIrIELAgtWE5ODtatW6ffHjlyJDZu3AiVSiVLvNmzZ+Ps2bNYvHgx3N3dZYlB1svcblvKfVVtX2Wd/gohAJTcNp5QauZ27ono3rEgtGBarRZarRYXLlxAWVkZFixYgBEjRiArK+uO4wDbKykpCRkZGZg0aRKSkpIk7ZsIML9xYAcWhsrWd8KGUizKPQwRIkapFCiprEVy7hEEqRRY9USA5PHM7dwT0b3j33Ir4OrqipCQEOTn52Pw4MF48sknUVlZqR/vZ6y0tDS8+eabmDBhArKystC1a1dJ+iUyZ3I+mbs6ZgRe3FhmMGZwnHdPrIgeDoWj9E8ZE5HlY0FoRVxcXDB69Gjk5OTgxIkTGDJkiNF9pqWlITU1FSEhIcjLy0P37t0lyJSIfo/C0Rbr44M4DyERSYYFoZU5ffo0AGmeAk5NTUVaWhoefvhhbNmyBQ4ODkb3SURt593bmYUgEUmCS9dZmOvXr+OHH35odd/atWtRUlKCQYMGwdvb22CfUqmEIAioqqpqU5yUlBSkpaVh3LhxLAaJiIjMHK8QWpirV69izJgx8PX1hb+/Pzw9PXHhwgWUlJSgtLQUTk5OWLt2bYv36XQ6AG27cqjVarF48WJ069YNo0aNwttvv93imJCQEISEhBj98xAREZH8WBBaGEdHR6SlpWHHjh3Yvn07ampqYGNjA6VSicTERLz00kvo16+fwXvq6+tRXV2N4OBgeHl53TVG81XEpqYmLF++/I7HsSAkIiIyDywILYyNjQ2Sk5ORnJzc5vcUFxdDp9O1ebqY1NRUpKamtjNDIiIi6mw4htCMxcXFQRAEREdHG9VPUVER/Pz8EB4eLlFmLcXExEAQBMTFxckWg4iIiNpHEEVRhqXQSU7l5eXIycnRb/v6+iIqKsp0CbVBVlYWKioq9NtqtRr+/v6mS4iIiIj0WBASERERWTneMiYiIiKyciwIzZBWq4UgCPqXsWMIO0LzGMLml1arNXVKRERE9B98ytiMRUZGwt/fH76+vgbtWVlZyMjIwMGDB3H27FnodDr0798fjz76KF555RV4enoaFVcURbz66qvYt28ffvrpJ9TW1sLV1RUDBw5EfHw8nnzyyRbrJGs0GgwaNAjl5eXIzc01Kj4RERFJi2MIzZBWq0VcXBzWrl2L2NjYFvv//Oc/44cffkBgYCA8PDwA3HoQZceOHXB1dUVxcTF8fHzaHb+pqQlOTk4YOXIkhgwZgvvuuw91dXUoKChAVVUVJk6ciK+//hpdurS8AH233Ik6C64TTETWhFcILdBHH30Ee3v7Fu2ffPIJZs2ahdTUVGRmZra7/27duqG+vr5FjKamJoSFhaGwsBD5+fmIiIhodwyiuzl/+bos/dY13MCC7MPYW1mrbwtSKZA+1RduDrayxHR3spOlXyKitmJBaIFaKwYB4PHHH8esWbNw4sQJWWJ069YNarUaO3bskCQGya/hRpOpU2i3EUu2ydJvFwFwtrfB+zMDEKhyw77KOiRlH0LYO7uhk+l+yo+LJ8rTsYwcbPnrg8iS8G+0FdmyZQsAtBhzKBWdToeCggJZY5C0HkwuNHUKnY5OBJZOHYqIYbeGW0QM84AIEXMzymSLaY6fQ9Uy3gEgsiQsCC1YTk4OysvL0dDQgCNHjqCwsBAqlQqLFy+WLEbzEnY1NTXYvn07jh07htjYWDzyyCOSxSDqaIEqN4PtUSqFiTIhIuoYLAgtWE5ODtatW6ffHjlyJDZu3AiVSiVZjLS0NP3/C4KA+fPn480335Ssf5KXOd6qbCbnVbV9lXX6K4QAUHLbeEI5mPPnQESWgQWhBdNqtdBqtbhw4QLKysqwYMECjBgxAllZWZgwYYIkMURRhE6nw+nTp7F582a88cYb2LNnD77++mu4uLhIEoPkY87jwA4sDJWl34QNpViUexgiRIxSKVBSWYvk3CMIUimw6okAWWKa8+dARJaB30JWwNXVFSEhIcjPz8fgwYPx5JNPorKyssVcge3VpUsXeHl5Yfbs2XB3d8e0adOQnp6Ot956S5L+iVoj15O5q2NG4MWNZQZjBsd598SK6OFQOMrzlDERkamxILQiLi4uGD16NHJycnDixAkMGTJE8hhhYWEAgJ07d0reN1FHUDjaYn18EOchJCKrwoLQypw+fRrArSlizLF/oo7i3duZhSARWQ2uZWxhrl+/jh9++KHVfWvXrkVJSQkGDRoEb29vg31KpRKCIKCqququMY4dO4Zz5861aG9oaMDLL78MAJg8efK9J09EREQmwcs4Fubq1asYM2YMfH194e/vD09PT1y4cAElJSUoLS2Fk5MT1q5d2+J9Op0OQNuu7BUUFOC1115DSEgIBgwYAFdXV1RXVyM/Px/nz59HcHCwvjAkIiKizo8FoYVxdHREWloaduzYge3bt6OmpgY2NjZQKpVITEzESy+9hH79+hm8p76+HtXV1QgODoaXl9ddY4SGhiI+Ph7FxcXYt28fLl26BFdXV/j6+iI6OhqzZs3iLWMiIiIzwt/aFsbGxgbJyclITk5u83uKi4uh0+mQlJTUpuN9fX2xatWq9qZIREREnQzHEJqxuLg4CIKA6Ohoo/opKiqCn58fwsPDJcqspZiYGAiCgLi4ONliEBERUfsIoijKtFw7yaW8vBw5OTn6bV9fX0RFRZkuoTbIyspCRUWFflutVsPf3990CREREZEeC0IiIiIiK8dbxkRERERWjgWhGdJqtRAEQf8ydgxhR2geQ9j80mq1pk6JiIiI/oMFoRmLjIxESkpKi/GDWVlZiIqKgre3N1xcXODk5AQfHx8kJiaiurpallwyMzP1xd7GjRtb7NdoNEhJSUFkZKQs8YmIiKj9OO2MGVOr1YiNjW3Rnp2djYMHDyIwMBAeHh4Abj2IsnLlSqxbtw7FxcXw8fGRLI9z584hISEBjo6OuHLlSqvHaDQaaDQaaLVa5ObmShabLAPXDSYiMi0WhBboo48+gr29fYv2Tz75BLNmzUJqaioyMzMlizdnzhw4OjriqaeewvLlyyXrlzqP85evy9JvXcMNLMg+jL2Vtfq2IJUC6VN94eZgK0tMdyc7WfolIjJnLAgtUGvFIAA8/vjjmDVrFk6cOCFZrIyMDGRlZaGwsBDff/+9ZP1ai4YbTaZOoU1GLNkmS79dBMDZ3gbvzwxAoMoN+yrrkJR9CGHv7IZOpvkPflw8UZ6OO5iDLb++iUg6/EaxIlu2bAFwa95CKZw5cwYvvPACnn76aYSFhbEgbIcHkwtNnYJJ6URg6dShiBh2a2hDxDAPiBAxN6NMtpiWcs6rlkWYOgUisiAsCC1YTk4OysvL0dDQgCNHjqCwsBAqlQqLFy+WpP/nnnsO9vb2vE1MRglUuRlsj1IpTJQJEZH1YkFowXJycrBu3Tr99siRI7Fx40aoVCqj+/7000/x1VdfITc3Fz169DC6P2tlLrcv5byqtq+yTn+FEABKbhtPKAdzOedERB2JBaEF02q10Gq1uHDhAsrKyrBgwQKMGDECWVlZmDBhQrv7PX36NBITExEdHY0pU6ZImLH1MZdxYAcWhsrSb8KGUizKPQwRIkapFCiprEVy7hEEqRRY9USALDHN5ZwTEXUkfjNaAVdXV4SEhCA/Px+DBw/Gk08+icrKStjY2LSrv4SEBHTt2hXvvvuuxJlSZyXXk7mrY0bgxY1lBmMGx3n3xIro4VA4yvOUMRERtcSC0Iq4uLhg9OjRyMnJwYkTJzBkyJB29VNeXo6amhrcd999re6fMWMGZsyYgXfeeQeJiYlGZEyWTuFoi/XxQZyHkIjIxFgQWpnTp08DALp1a/9HHx0djZqamhbtpaWlKCsrw/jx4zFgwADJnmYmy+fd25mFIBGRCbEgtDDXr19HWVkZRo8e3WLf2rVrUVJSgkGDBsHb29tgn1KpxKlTp1BZWQmlUvm7MZYtW9Zqe2pqKsrKyvDss8+axfrKREREdAsLQgtz9epVjBkzBr6+vvD394enpycuXLiAkpISlJaWwsnJCWvXrm3xPp1OB8C4K4dERERknvjb38I4OjoiLS0NO3bswPbt21FTUwMbGxsolUokJibipZdeQr9+/QzeU19fj+rqagQHB8PLy8tEmRMREZGpsCC0MDY2NkhOTkZycnKb31NcXAydToekpCSjYqempiI1NdWoPoiIiKjjdTF1AtR+cXFxEATB6PF6RUVF8PPzQ3h4uESZtRQTEwNBEBAXFydbDCIiImofQRRFmZaQJ7mUl5cjJydHv+3r64uoqCjTJdQGWVlZqKio0G+r1Wr4+/ubLiEiIiLSY0FIREREZOV4y5iIiIjIyrEgNENarRaCIOhf5jDnX/MYwuaXVqs1dUpERET0HywIzVhkZCRSUlJajB/MyspCVFQUvL294eLiAicnJ/j4+CAxMRHV1dWSxA4JCTEo8G5/TZo0qcXxGo0GKSkpiIyMlCQ+ERERSYfTzpgxtVqN2NjYFu3Z2dk4ePAgAgMD4eHhAeDWgygrV67EunXrUFxcDB8fH0lySElJadE2aNCgFm0ajQYajQZarRa5ubmSxCaizo1rVBOZDxaEFuijjz6Cvb19i/ZPPvkEs2bNQmpqKjIzMyWJxXkHiczX+cvXZem3ruEGFmQfxt7KWn1bkEqB9Km+cHOwlTyeu5Od5H0SWRsWhBaotWIQAB5//HHMmjULJ06c6OCMiCxbw40mU6fQLiOWbJOl3y4C4Gxvg/dnBiBQ5YZ9lXVIyj6EsHd2QyfDvBY/Lp4ofacyc7Dlr1/qXPgn0ops2bIFwK15C6WyceNGVFZWwtHREYGBgRgzZoxkfROZiweTC02dQqeiE4GlU4ciYtitISsRwzwgQsTcjDJZ4pnj+a9aFmHqFIgMsCC0YDk5OSgvL0dDQwOOHDmCwsJCqFQqLF68WLIYM2bMMNgODAzEpk2boFKpJItBROYnUOVmsD1KpTBRJkTUFiwILVhOTg7WrVun3x45ciQ2btwoSbGmVqvx+uuvw9/fHy4uLjh+/Dj+9re/4dNPP0VoaCgOHToEBwcHo+MQmQNzvGUJyHtlbV9lnf4KIQCU3DaeUGrmev6JOhMWhBZMq9VCq9XiwoULKCsrw4IFCzBixAhkZWVhwoQJRvWdmJhosO3n54d169ahqakJGRkZWLt2LZ5//nmjYhCZC3MdD3ZgYags/SZsKMWi3MMQIWKUSoGSylok5x5BkEqBVU8ESB7PXM8/UWfCv0VWwNXVFSEhIcjPz8fgwYPx5JNPorKyEjY2NpLHio+PR0ZGBr777jsWhESdnFxP566OGYEXN5YZjBkc590TK6KHQ+Eo/VPGRGQ8FoRWxMXFBaNHj0ZOTg5OnDiBIUOGSB6jZ8+eAICGhgbJ+yYi86BwtMX6+CDOQ0hkRlgQWpnTp08DALp1k+ej37t3LwBAqVTK0j8RmQ/v3s4sBInMBJeuszDXr1/HDz/80Oq+tWvXoqSkBIMGDYK3t7fBPqVSCUEQUFVVddcYJ0+exK+//tqi/ejRo1iwYAEAmMX6ykRERHQLrxBamKtXr2LMmDHw9fWFv78/PD09ceHCBZSUlKC0tBROTk5Yu3Zti/fpdDoAbbtyuHv3bjzzzDMYP348Bg4cCGdnZxw/fhxbtmxBY2MjkpOTMXr0aMl/NiIiIpIHC0IL4+joiLS0NOzYsQPbt29HTU0NbGxsoFQqkZiYiJdeegn9+vUzeE99fT2qq6sRHBwMLy+vu8YICAjAjBkzsH//fpSUlODKlStwd3fH5MmT8fzzzyMsLEyuH4+IiIhkwILQwtjY2CA5ORnJycltfk9xcTF0Oh2SkpLadPywYcPw6aeftjdFIiIi6mQ4htCMxcXFQRAEo8frFRUVwc/PD+Hh4RJl1lJMTAwEQUBcXJxsMYiIiKh9BFEUZVhqnORUXl6OnJwc/bavry+ioqJMl1AbZGVloaKiQr+tVqvh7+9vuoSIiIhIjwUhERERkZXjLWMiIiIiK8eC0AxptVoIgqB/mcOcf81jCJtfWq3W1CkRERHRf7AgNGORkZFISUlpMX4wKysLUVFR8Pb2houLC5ycnODj44PExERUV1dLmsPOnTsRGRmJXr16wc7ODvfffz+mTp2KgwcPGhyn0WiQkpKCyMhISeMTERGR8TjtjBlTq9WIjY1t0Z6dnY2DBw8iMDAQHh4eAG49iLJy5UqsW7cOxcXF8PHxMTp+eno6Fi5ciL59+0KtVqNnz544e/YsvvvuOxw6dAh+fn76YzUaDTQaDbRaLXJzc42OTWTOuMYvEXU2LAgt0EcffQR7e/sW7Z988glmzZqF1NRUZGZmGhXjq6++wsKFC6FWq5GRkYHu3bsb7G9qajKqfyJTOn/5uiz91jXcwILsw9hbWatvC1IpkD7VF24OtrLEdHeyk6VfIrIsLAgtUGvFIAA8/vjjmDVrFk6cOGF0jNdffx3Ozs7QarUtikGgbUvgkWVruGG+/ygYsWSbLP12EQBnexu8PzMAgSo37KusQ1L2IYS9sxs6meZ7+HHxRHk67iAOtvwuIeoI/JtmRbZs2QLg1ryFxqioqMDRo0eh0Wjg5OSE/Px8VFRUwMHBAQ899JDBrWKyXg8mF5o6hU5HJwJLpw5FxLBbQzkihnlAhIi5GWWyxTT3z6FqWYSpUyCyCiwILVhOTg7Ky8vR0NCAI0eOoLCwECqVCosXLzaq3/379wMA3N3dMXbsWPzwww8G+5944gmsWbMGtrby3AIjMmeBKjeD7VEqhYkyISL6LxaEFiwnJwfr1q3Tb48cORIbN26ESqUyqt9z584BANasWQOVSoVvv/0WgYGBOH78OJ5//nls2LABnp6eeOutt4yKQ+bNnG9VynlVbV9lnf4KIQCU3DaeUA7m/DkQUcdhQWjBtFottFotLly4gLKyMixYsAAjRoxAVlYWJkyY0O5+dTqd/r+ff/45hg8fDgAYPnw4cnJyMGjQILz33ntYvHgx7Ow4oN1amfPYrwMLQ2XpN2FDKRblHoYIEaNUCpRU1iI59wiCVAqseiJAlpjm/DkQUcfhN4UVcHV1RUhICPLz8zF48GA8+eSTqKyshI2NTbv7AwAvLy99MdisV69eCAoKwrZt23D06FGuV0xmSa4nc1fHjMCLG8sMxgyO8+6JFdHDoXDkEAsiMh0WhFbExcUFo0ePRk5ODk6cOIEhQ4a0q5/BgwcDAHr06NHq/ub2q1evtqt/IkulcLTF+vggzkNIRJ0OC0Irc/r0aQDGTQszevRodO/eHSdPnsS1a9daTHNz9OhRAIBSqWx3DCJL5t3bmYUgEXUqXLrOwly/fr3FU7/N1q5di5KSEgwaNAje3t4G+5RKJQRBQFVV1V1jODk54c9//jOuXLmCJUuWGOxbv349jhw5grFjx+pXSSEiIqLOjVcILczVq1cxZswY+Pr6wt/fH56enrhw4QJKSkpQWloKJycnrF27tsX7mh8UaeuVw6VLl2Lnzp1IT09HcXExRo4ciePHjyMvLw9ubm744IMPJP25iIiISD4sCC2Mo6Mj0tLSsGPHDmzfvh01NTWwsbGBUqlEYmIiXnrpJfTr18/gPfX19aiurkZwcDC8vLzaFMfd3R179uxBWloasrOz8f3330OhUCAmJgapqakYMGCAHD8eERERyYAFoYWxsbFBcnIykpOT2/ye4uJi6HQ6JCUl3VMshUKBFStWYMWKFfeaJhEREXUiHENoxuLi4iAIAqKjo43qp6ioCH5+fggPD5cos5ZiYmIgCALi4uJki0FERETtI4iiKNOS6iSX8vJy5OTk6Ld9fX0RFRVluoTaICsrCxUVFfpttVrNOQqJiIg6CRaERERERFaOt4yJiIiIrBwLQjOk1WohCIL+ZewYwo7QPIaw+aXVak2dEhEREf0HC0IzFhkZiZSUlLuOH6yvr4enpycEQcCkSZMkiZ2Xl4cXXngBwcHBcHR0hCAISE1NvePxGo0GKSkpiIyMlCQ+ERERSYfTzpgxtVqN2NjYux43b948XLhwQdLYy5cvx65du+Di4oK+ffvixIkTv3u8RqOBRqOBVqtFbm6upLkQmRLXJSYiS8CC0MLl5eVh/fr1WLlyJebNmydZv3/5y1/Qp08fDBo0CJs2bcKMGTMk65tISucvX5el37qGG1iQfRh7K2v1bUEqBdKn+sLNwVbyeO5OdpL3SUTUjAWhBautrcWzzz6LmTNn4rHHHpO0IBw3bpxkfZF5abjRZOoU7smIJdtk6beLADjb2+D9mQEIVLlhX2UdkrIPIeyd3dDJMHfDj4snSt9pB3Cw5a8ZInPAv6kWbO7cubh58yZWrlyJS5cumTodshAPJheaOoVOQScCS6cORcQwDwBAxDAPiBAxN6NMlnjmet6rlkWYOgUiagMWhBYqOzsbn332GTZt2gR3d3cWhEQyCFS5GWyPUilMlAkRkXFYEFqgmpoazJ49G2q1GtOmTTN1OmRhzO3WpZxX1vZV1umvEAJAyW3jCaVmbuediMwLC0ILlJCQgMbGRqxevdrUqZAFMrcxYQcWhsrSb8KGUizKPQwRIkapFCiprEVy7hEEqRRY9USA5PHM7bwTkXnhN4yFyc3NRWZmJrRaLfr06WPqdIhMTq6nc1fHjMCLG8sMxgyO8+6JFdHDoXCU/iljIiI5sSC0MGVlt345xcbGtjpHYWFhIQRBgJ+fH8rLyzs2OSILonC0xfr4IM5DSEQWgQWhhQkICEB8fHyL9suXL2PTpk3w8vLCxIkT0a9fPxNkR2R5vHs7sxAkIrPHgtDCTJkyBVOmTGnRXlVVhU2bNsHHxwcff/xxi/1KpRKnTp1CZWUllEplB2RKREREnQULQgIA6HQ6AEC3bm37I5GTk4OcnBwAQGVlpb6tqqoKADB27FjMmjVL8jyJiIhIeiwICfX19aiurkZwcDC8vLza9J7y8nKsW7fOoO3gwYM4ePCgfpsFIRERkXlgQWgllEolRLH19bSKi4uh0+mQlJTU5v5SU1ORmpoqUXZERERkSl1MnQC1X1xcHARBQHR0tFH9FBUVwc/PD+Hh4RJl1lJMTAwEQUBcXJxsMYiIiKh9BPFOl42o0yovL9eP3wMAX19fREVFmS6hNsjKykJFRYV+W61Ww9/f33QJERERkR4LQiIiIiIrx1vGRERERFaOBaEZ0mq1EARB/zJ2DGFHaB5D2PzSarWmTomIiIj+gwWhGYuMjERKSspdxw/W19fD09MTgiBg0qRJRsc9f/48PvzwQ0yZMgUDBgyAnZ0devbsicmTJ6OwsLDV92g0GqSkpCAyMtLo+ERERCQtTjtjxtRqdavrFf/WvHnzcOHCBcniZmZmYs6cOfD09MSECRPg6emJf/3rX/jyyy9RUFCAt99+G/Pnzzd4j0ajgUajgVarRW5urmS5EFHH4JrNRJaNBaGFy8vLw/r167Fy5UrMmzdPkj4feOABbN68GZMnT0aXLv+9yLxw4UIEBQUhKSkJM2fORN++fSWJR0Rtc/7ydcn7rGu4gQXZh7G3slbfFqRSIH2qL9wcbCWPBwDuTnay9EtEd8aC0ILV1tbi2WefxcyZM/HYY49JVhBOmDCh1fbBgwdj+vTp+PDDD/H99993+qlwiH6r4UaTqVMwyogl2yTvs4sAONvb4P2ZAQhUuWFfZR2Ssg8h7J3d0Mk0R8WPiyfK07FMHGz5q5TMH/8UW7C5c+fi5s2bWLlyJS5dutQhMW1sbAC0fU1kos7kweTWx8BaM50ILJ06FBHDPAAAEcM8IELE3Iwy2WKa2+dQtSzC1CkQGY2/tS1UdnY2PvvsM2zatAnu7u4dUhBeunQJX3zxBezt7TFu3DjZ4xFRxwhUuRlsj1IpTJQJEcmFBaEFqqmpwezZs6FWqzFt2rQOizt79mycPXsWixcvhru7e4fFJZKKud2q/C25rqztq6zTXyEEgJLbxhPKwdw/ByJzxILQAiUkJKCxsRGrV6/usJhJSUnIyMjApEmTkJSU1GFxiaRk7mPBDiwMlbzPhA2lWJR7GCJEjFIpUFJZi+TcIwhSKbDqiQDJ4wHm/zkQmSP+rbMwubm5yMzMhFarRZ8+fTokZlpaGt58801MmDABWVlZ6Nq1a4fEJSJDcjyduzpmBF7cWGYwZnCcd0+siB4OhaM8TxkTUcdjQWhhyspufWnHxsa2OkdhYWEhBEGAn58fysvLjY6XlpaG1NRUhISEIC8vD927dze6TyLqPBSOtlgfH8R5CIksHAtCCxMQEID4+PgW7ZcvX8amTZvg5eWFiRMnol+/fkbHSk1NRVpaGh5++GFs2bIFDg4ORvdJRJ2Td29nFoJEFowFoYWZMmUKpkyZ0qK9qqoKmzZtgo+PDz7++OMW+5VKJU6dOoXKykoolcq7xklJScHixYsxbtw4FoNERERmjgUhAQB0Oh2Ats0fqNVqsXjxYnTr1g2jRo3C22+/3eKYkJAQhISESJ0mERERyYAFIaG+vh7V1dUIDg6Gl5fXXY+vqqoCADQ1NWH58uV3PI4FIRERkXlgQWgllEolRLH1daaKi4uh0+naPF1MamoqUlNTJcyOiIiITKmLqROg9ouLi4MgCIiOjjaqn6KiIvj5+SE8PFyizFqKiYmBIAiIi4uTLQYRERG1jyDe6bIRdVrl5eXIycnRb/v6+iIqKsp0CbVBVlYWKioq9NtqtRr+/v6mS4iIiIj0WBASERERWTneMiYiIiKyciwIzZBWq4UgCPqXsWMIO0LzGMLml1arNXVKRERE9B8sCM1YZGQkUlJS7jp+sL6+Hp6enhAEAZMmTTI6riiKeOWVVxASEoK+ffvC3t4evXv3xv/3//1/+OSTT9DY2NjiPRqNBikpKYiMjDQ6PhEREUmL086YMbVa3ep6xb81b948XLhwQbK4N2/exLvvvouRI0ciIiIC9913H+rq6lBQUIBZs2YhMzMTX3/9Nbp0+e+/NzQaDTQaDbRaLXJzcyXLhagZ19olImo/FoQWLi8vD+vXr8fKlSsxb948Sfrs1q0b6uvrYW9vb9De1NSEsLAwFBYWIj8/HxEREZLEI/N3/vJ12fqua7iBBdmHsbeyVt8WpFIgfaov3BxsJY3l7mQnaX9ERJ0FC0ILVltbi2effRYzZ87EY489JllBCKBFMQjcKhTVajV27NiBEydOSBaLgIYbTaZOwSgjlmyTre8uAuBsb4P3ZwYgUOWGfZV1SMo+hLB3dkMn8RwKPy6eKG2HHcDBll/zRHR3/KawYHPnzsXNmzexcuVKXLp0SfZ4Op0OBQUFAG7NjUjSeTC50NQpdFo6EVg6dSgihnkAACKGeUCEiLkZZZLHMsfPoWoZr9QT0d2xILRQ2dnZ+Oyzz7Bp0ya4u7vLVhA2L2FXU1OD7du349ixY4iNjcUjjzwiSzyi1gSq3Ay2R6kUJsqEiMg8sSC0QDU1NZg9ezbUajWmTZsma6y0tDT9/wuCgPnz5+PNN9+UNaY1MsdblbeT+8ravso6/RVCACi5bTyhlMz9cyAiuhMWhBYoISEBjY2NWL16teyxRFGETqfD6dOnsXnzZrzxxhvYs2cPvv76a7i4uMge31qY+ziwAwtDZes7YUMpFuUehggRo1QKlFTWIjn3CIJUCqx6IkDSWOb+ORAR3Qm/3SxMbm4uMjMzodVq0adPnw6J2aVLF3h5eWH27Nlwd3fHtGnTkJ6ejrfeeqtD4lPnJ+fTuatjRuDFjWUGYwbHeffEiujhUDhK+5QxEZGlYkFoYcrKbv1SjI2NbXWOwsLCQgiCAD8/P5SXl0sePywsDACwc+dOyfsmao3C0Rbr44M4DyERkRFYEFqYgIAAxMfHt2i/fPkyNm3aBC8vL0ycOBH9+vWTJf7p06cB3JqChqgjefd2ZiFIRNRO/K1tYaZMmYIpU6a0aK+qqsKmTZvg4+ODjz/+uMV+pVKJU6dOobKyEkql8ndjHDt2DAqFAr169TJob2howMsvvwwAmDx5cvt/CCIiIupQLAgJwK05BIG2XdkrKCjAa6+9hpCQEAwYMACurq6orq5Gfn4+zp8/j+DgYH1hSERERJ0fC0JCfX09qqurERwcDC8vr7seHxoaivj4eBQXF2Pfvn24dOkSXF1d4evri+joaMyaNYu3jImIiMwIf2tbCaVSCVFsfR2v4uJi6HQ6JCUltakvX19frFq1Ssr0iIiIyIS6mDoBar+4uDgIgoDo6Gij+ikqKoKfnx/Cw8MlyqylmJgYCIKAuLg42WIQERFR+wjinS4bUadVXl6OnJwc/bavry+ioqJMl1AbZGVloaKiQr+tVqvh7+9vuoSIiIhIjwUhERERkZXjLWMiIiIiK8eC0AxptVoIgqB/GTuGsCM0jyFsfmm1WlOnRERERP/BgtCMRUZGIiUl5a7jB+vr6+Hp6QlBEDBp0iRZcsnMzNQXexs3bmyxX6PRICUlBZGRkbLEJyIiovbjtDNmTK1Wt7pe8W/NmzcPFy5ckC2Pc+fOISEhAY6Ojrhy5Uqrx2g0Gmg0Gmi1WuTm5sqWC3Ucrh1MRGQ5WBBauLy8PKxfvx4rV67EvHnzZIkxZ84cODo64qmnnsLy5ctliUH35vzl67L1XddwAwuyD2NvZa2+LUilQPpUX7g52MoS093JTpZ+iYjoFhaEFqy2thbPPvssZs6ciccee0yWgjAjIwNZWVkoLCzE999/L3n/nVXDjSZTp/C7RizZJlvfXQTA2d4G788MQKDKDfsq65CUfQhh7+yGTqY5C35cPFGejk3IwZZfv0TUefAbyYLNnTsXN2/exMqVK3Hp0iXJ+z9z5gxeeOEFPP300wgLC7OqgvDB5EJTp2AyOhFYOnUoIoZ5AAAihnlAhIi5GWWyxbTE8121LMLUKRAR6bEgtFDZ2dn47LPPsGnTJri7u8tSED733HOwt7fnbWIrFKhyM9gepVKYKBMiIpICC0ILVFNTg9mzZ0OtVmPatGmyxPj000/x1VdfITc3Fz169JAlRmfW2W9hyn1FbV9lnf4KIQCU3DaeUA6d/XwTEZk7FoQWKCEhAY2NjVi9erUs/Z8+fRqJiYmIjo7GlClTZInR2XX28V8HFobK1nfChlIsyj0MESJGqRQoqaxFcu4RBKkUWPVEgCwxO/v5JiIyd/yWtTC5ubnIzMyEVqtFnz59ZImRkJCArl274t1335WlfzKenE/lro4ZgRc3lhmMGRzn3RMroodD4SjPU8ZERCQvFoQWpqzs1i/p2NjYVucoLCwshCAI8PPzQ3l5ebtilJeXo6amBvfdd1+r+2fMmIEZM2bgnXfeQWJiYrtiUOelcLTF+vggzkNIRGRBWBBamICAAMTHx7dov3z5MjZt2gQvLy9MnDgR/fr1a3eM6Oho1NTUtGgvLS1FWVkZxo8fjwEDBsDX17fdMajz8+7tzEKQiMhCsCC0MFOmTGl1XF9VVRU2bdoEHx8ffPzxxy32K5VKnDp1CpWVlVAqlb8bY9myZa22p6amoqysDM8++6xZrK9MREREt3AtYwIA6HQ6AEC3bvw3AhERkbVhQUior69HdXU1goOD4eXlZep0iIiIqIPxcpCVUCqVEMXW1xUrLi6GTqdDUlKSUTFSU1ORmppqVB9ERETU8XiF0IzFxcVBEASjx+sVFRXBz88P4eHhEmXWUkxMDARBQFxcnGwxiIiIqH0E8U6XjajTKi8vR05Ojn7b19cXUVFRpkuoDbKyslBRUaHfVqvV8Pf3N11CREREpMeCkIiIiMjK8ZYxERERkZVjQWiGtFotBEHQv8xhzr/mMYTNL61Wa+qUiIiI6D9YEJqxyMhIpKSk3HX8YH19PTw9PSEIAiZNmiRJ7JCQEIMC7/ZXazE0Gg1SUlIQGRkpSXwiIiKSDqedMWNqtbrV9Yp/a968ebhw4YIsOaSkpLRoGzRoUIs2jUYDjUYDrVaL3NxcWXIhsnRcP5qI5MKC0MLl5eVh/fr1WLlyJebNmyd5/5x3kOi/zl++Lku/dQ03sCD7MPZW1urbglQKpE/1hZuDrSwx3Z3sZOmXiDonFoQWrLa2Fs8++yxmzpyJxx57TJaCkEgODTeaTJ1Cu4xYsk2WfrsIgLO9Dd6fGYBAlRv2VdYhKfsQwt7ZDZ1M80T8uHiiPB3LxMGWv86IjMG/QRZs7ty5uHnzJlauXIlLly7JEmPjxo2orKyEo6MjAgMDMWbMGFnikHV5MLnQ1Cl0KjoRWDp1KCKGeQAAIoZ5QISIuRllssU0t8+galmEqVMgMmssCC1UdnY2PvvsM2zatAnu7u6yFYQzZsww2A4MDMSmTZugUqlkiUdkrQJVbgbbo1QKE2VCRJaIBaEFqqmpwezZs6FWqzFt2jRZYqjVarz++uvw9/eHi4sLjh8/jr/97W/49NNPERoaikOHDsHBwUGW2GT5zO12ZTM5r6rtq6zTXyEEgJLbxhPKwVw/AyJqHxaEFighIQGNjY1YvXq1bDESExMNtv38/LBu3To0NTUhIyMDa9euxfPPPy9bfLJs5joe7MDCUFn6TdhQikW5hyFCxCiVAiWVtUjOPYIglQKrngiQJaa5fgZE1D78G29hcnNzkZmZCa1Wiz59+nR4/Pj4eGRkZOC7775jQUhWR64nc1fHjMCLG8sMxgyO8+6JFdHDoXCU5yljIrIuLAgtTFnZrV8YsbGxrc5RWFhYCEEQ4Ofnh/Lycsnj9+zZEwDQ0NAged9E1krhaIv18UGch5CIZMOC0MIEBAQgPj6+Rfvly5exadMmeHl5YeLEiejXr58s8ffu3QsAUCqVsvRPZM28ezuzECQiWbAgtDBTpkzBlClTWrRXVVVh06ZN8PHxwccff9xiv1KpxKlTp1BZWXnXYu7kyZPo3r07PDw8DNqPHj2KBQsWAIBZrK9MREREt7AgJACATqcDAHTrdvc/Ert378YzzzyD8ePHY+DAgXB2dsbx48exZcsWNDY2Ijk5GaNHj5Y7ZSIiIpIIC0JCfX09qqurERwcDC8vr7seHxAQgBkzZmD//v0oKSnBlStX4O7ujsmTJ+P5559HWFhYB2RNREREUmFBaCWUSiVEsfU1roqLi6HT6ZCUlNSmvoYNG4ZPP/1UyvSIiIjIhLqYOgFqv7i4OAiCYPR4vaKiIvj5+SE8PFyizFqKiYmBIAiIi4uTLQYRERG1jyDe6bIRdVrl5eXIycnRb/v6+iIqKsp0CbVBVlYWKioq9NtqtRr+/v6mS4iIiIj0WBASERERWTneMiYiIiKyciwIzZBWq4UgCPqXOcz51zyGsPml1WpNnRIRERH9BwtCMxYZGYmUlJS7jh+sr6+Hp6cnBEHApEmTJM1h586diIyMRK9evWBnZ4f7778fU6dOxcGDBw2O02g0SElJQWRkpKTxiYiIyHicdsaMqdXqVtcr/q158+bhwoULksdPT0/HwoUL0bdvX6jVavTs2RNnz57Fd999h0OHDsHPz09/rEajgUajgVarRW5uruS5EFHH4ZrKRJaHBaGFy8vLw/r167Fy5UrMmzdPsn6/+uorLFy4EGq1GhkZGejevbvB/qamJsliEdG9OX/5uiz91jXcwILsw9hbWatvC1IpkD7VF24OtpLHc3eyk7xPImodC0ILVltbi2effRYzZ87EY489JmlB+Prrr8PZ2RlarbZFMQi0bQk8os6u4YZ5/sNmxJJtsvTbRQCc7W3w/swABKrcsK+yDknZhxD2zm7oZJiv4sfFE6XvtAM42PL7j8wP/9RasLlz5+LmzZtYuXIlLl26JFm/FRUVOHr0KDQaDZycnJCfn4+Kigo4ODjgoYceMrhVTGTOHkwuNHUKnYpOBJZOHYqIYR4AgIhhHhAhYm5GmSzxzPX8Vy2LMHUKRPeMBaGFys7OxmeffYZNmzbB3d1d0oJw//79AAB3d3eMHTsWP/zwg8H+J554AmvWrIGtrfS3kIjItAJVbgbbo1QKE2VCRFJiQWiBampqMHv2bKjVakybNk3y/s+dOwcAWLNmDVQqFb799lsEBgbi+PHjeP7557FhwwZ4enrirbfekjw2UUcy11uWcl5Z21dZp79CCAAlt40nlJq5nn8ic8SC0AIlJCSgsbERq1evlqV/nU6n/+/nn3+O4cOHAwCGDx+OnJwcDBo0CO+99x4WL14MOzsOCifzZa5jwQ4sDJWl34QNpViUexgiRIxSKVBSWYvk3CMIUimw6okAyeOZ6/knMkf822ZhcnNzkZmZCa1Wiz59+sgSw9XVFQDg5eWlLwab9erVC0FBQdi2bRuOHj3K9YqJTECup3NXx4zAixvLDMYMjvPuiRXRw6Fw5BARInPGgtDClJXd+qKOjY1tdY7CwsJCCIIAPz8/lJeXtyvG4MGDAQA9evRodX9z+9WrV9vVPxF1TgpHW6yPD+I8hEQWiAWhhQkICEB8fHyL9suXL2PTpk3w8vLCxIkT0a9fv3bHGD16NLp3746TJ0/i2rVrsLe3N9h/9OhRAIBSqWx3DCLqvLx7O7MQJLIwLAgtzJQpUzBlypQW7VVVVdi0aRN8fHzw8ccft9ivVCpx6tQpVFZW3rWQc3Jywp///Gd8+OGHWLJkCZYsWaLft379ehw5cgRjx46Fh4fH7/RCREREnQULQgLw3wdF2jqh9NKlS7Fz506kp6ejuLgYI0eOxPHjx5GXlwc3Nzd88MEHcqZLREREEupi6gTI9Orr61FdXY3g4GB4eXm16T3u7u7Ys2cP5s2bh5MnT2LlypXYu3cvYmJisH//fjz44IMyZ01ERERS4RVCK6FUKiGKra8tVVxcDJ1Oh6SkpHvqU6FQYMWKFVixYoUUKRIREZGJsCA0Y3FxcYiLi8P06dOxcePGdvdTVFQEPz8/hIeHS5idoZiYGGzYsOGe3yeKoqSrrBAREVkbZ2dnCILwu8cI4p0uG1GnVV5ejpycHP22r68voqKiTJdQG2RlZaGiokK/rVar2zRH4cWLF/XzHhIREdG9u3DhAlxcXH73GBaE1KlJeYXw4sWLuP/++/HLL7/c9S8GGeK5Mw7PX/vx3LUfz137Wdq5a8sVQt4ypk5NEATJ/zK6uLhYxF9wU+C5Mw7PX/vx3LUfz137WdO541PGRERERFaOBSERERGRlWNBSFbDzs4OKSkpsLOzM3UqZofnzjg8f+3Hc9d+PHftZ43njg+VEBEREVk5XiEkIiIisnIsCImIiIisHAtCIiIiIivHgpCIiIjIyrEgJCIiIrJyLAjJIp05cwazZs2Ch4cH7O3t8cADD2Dx4sW4ceOGUf0mJCRAEAQIgoAzZ85IlG3nI8X5O378OJYuXYqHHnoIffv2ha2tLe6//348+eSTOHbsmIzZy2/fvn0IDw+Hm5sbHB0dMWrUKGRkZNxTHzqdDu+99x6GDRuG7t2747777sO0adNw/PhxmbLuPIw9f8XFxfif//kfjBgxAu7u7rC3t8cf/vAHvPbaa6ivr5cv8U5Aij97t2tsbIS/vz8EQcAf/vAHCTPtfKQ6d5cuXUJKSgp8fX3h4OCAHj16ICAgAGlpaTJk3YFEIgvz66+/iv369RMFQRCnTp0qvvbaa2JwcLAIQJw0aZJ48+bNdvW7bds2URAE0dHRUQQg/vrrrxJn3jlIdf6mT58uAhB9fX3F2bNni6+++qo4efJkEYDYvXt3cffu3TL/JPLYsWOHaGtrKzo5OYmzZs0S/+d//kdUqVQiADE9Pb3N/TzzzDMiAPHBBx8UX3nlFfHJJ58U7ezsRFdXV/HIkSMy/gSmJcX56927t9i1a1fx4YcfFhMTE8WXXnpJHD58uAhAHDhwoHj27FmZfwrTkOrP3u0WLVqk/04bPHiwxBl3HlKdu1OnTokDBw4UBUEQH330UfHVV18VX3zxRTEiIkIcOnSojD+B/FgQksV58sknRQDiqlWr9G06nU586qmnRADimjVr7rnPixcviv379xc1Go348MMPW3RBKNX5W7t2rVheXt6i/bPPPtMXQuamsbFRHDhwoGhnZyeWlpbq2y9evCj6+PiI3bp1E3/66ae79vPtt9+KAMRx48aJ165d07c3/6PjoYcekiV/U5Pq/C1btkw8ffq0QZtOpxPnzJkjAhATEhIkz93UpDp3tztw4IDYrVs3ceXKlRZdEEp17pqamsTAwECxe/fu4rfffttqHHPGgpAsysWLF0U7OztxwIABok6nM9h3+vRpsUuXLuKYMWPuud9nnnlGVCgU4pkzZyy6IJTr/P3WAw88IAIQ//3vfxvdV0cqLCwUAYhxcXEt9m3cuFEEIL7xxht37WfGjBkiAHHXrl0t9k2aNEkEIP7jH/+QJOfORKrzdyenT58WAYg+Pj7GpNkpSX3url+/Lg4dOlQcO3asqNPpLLoglOrcNR+7aNEiOdI0uW7y3pAm6lh79uzB9evX8eijj0IQBIN9Hh4eGDp0KPbu3Ytr167B3t6+TX1+8803+Oijj/Dpp5+id+/ecqTdachx/lpjY2MDAOjWzby+gnbu3AkACAsLa7GvuW3Xrl1t6sfR0RHBwcEt9k2cOBEFBQXYtWsXHnjgAeMS7mSkOn93Yq5/rtpC6nOXmpqK48eP4+DBgy3+rlsaqc7dpk2bAACPP/44fvnlF2zZsgX19fUYOHAgJk+eDCcnJ+mSNgE+VEIWpXlAvre3d6v7vb29odPpcPLkyTb1d/HiRcyaNQvh4eH485//LFmenZXU5681JSUlOHLkCAIDA9GjR49292MKv3d+3Nzc0LNnz7s+FHLlyhX8+uuvUKlU6Nq1a4v9zX1b4sMlUpy/37NmzRoArf/iN3dSnrt9+/bhr3/9K9LS0izuHx2tkerc7d+/H8Cth5q8vb0xZ84cvPHGG5g2bRoGDBigLzzNFQtCsigXLlwAALi6ura638XFxeC4u0lMTMSFCxfwwQcfSJNgJyf1+Wut/6eeegpdunTBX//61/YlaUJtOT93Ozdyn+POTIrzdyfl5eVIS0tDr1698Oqrr7Y7x85KqnN3/fp1xMbGYvjw4fif//kfSXPsrKQ6d+fOnQMAvPDCC0hMTMQvv/yCf//731i5ciUuXLgAtVqNX3/9VbrEOxgLQuqUevbsqZ/epS0vOf5llp+fj7Vr1+Kvf/0rvLy8JO9fTp3h/P3WtWvXoNFocOzYMfzlL39BSEiI7DHJOlRWVuKPf/wjbt68iY0bN6Jnz56mTqnTWrRoEY4fP441a9a0eoWa7kyn0wEA/vjHP2LZsmXw8vJCz5498cILL+Cll17ChQsX8Mknn5g4y/azvIEWZBFmzJiBS5cutfn4Pn36APjvvwDv9K+9ixcvGhx3Jw0NDXjmmWcwfvx4PPvss23Oo7Mw9fn7revXr2Pq1Kn49ttv8cYbbyApKeme3t9ZtOX83O3cyHWOzYEU5++3Tp06hfHjx+Pf//43vvzyS4wfP97oPDsjKc5daWkp/va3v2HRokUYOnSo5Dl2VlL9uXN1dUVNTQ2mTJnSYt9jjz2Gt956S39b2RyxIKRO6d13323X++42/ur48ePo0qULBgwY8Lv9nDt3DtXV1aiurkaXLq1fSPfw8AAAlJWV4f9v7/5B2geiOIC/ghrFSoeAfwZpEdHBpVDBIkjVQR0ERZwcHEx1cRMRt1QQwdWxs4uIi1tR7OTgJChYXIyCxsGCU6ki4f0GibSmWn96aTX3/UCWu3C8exz6SJO7cDj8rXjdUu38FXp6eqLJyUlKpVK0srJCGxsb34rtNyjMTyQSKep7fHykbDZL/f39n47R2NhIbW1tZBgGWZbleEpT7j3Ov0xE/gpdX1/T0NAQmaZJu7u7ND4+LjTe30RE7s7OzsiyLEokEpRIJBz9l5eX5PP5KBAIeGqDb1Hrrru7m7LZbMl3n+22fD7/43irBT8Zg6dEo1FSFIUODg6ImYv67u/v6fz8nPr6+sp+IdvU1ESappW87KdpMzMzpGkaqarq2nwqTVT+bIXF4PLyMm1ubroRdsXEYjEiev3y/D27zb6n3Di5XI6Oj48dfalU6svj/DWi8kf0WgwODg7S3d0d7ezs0MTEhLhAfyERuevq6vrw7xrR6xMwTdNodnZWcPTVJWrdDQ8PExHRxcWFo89uC4VC3w2z+qq97w2AaP+7sXIul+NMJsM3NzdfGt/L+xAyi8tfPp/nkZERJiJeWlqqSOxue3l54Y6ODlYUhU9PT9/aCze4Ldw/8OHhgTOZjGO/xcKNqZ+fn9/aZdiYWkT+DMPgYDDINTU1vLe3V6nwq0pU7j5CHt6HUFTurq6uWFEUbm5u5tvb26JxwuEwExEfHh66Ph+3oCAEzzFNk9vb29nn8/HU1BSvrq6+Hb02OjrqOHotnU4zEXEsFvvS+F4vCEXlzy4gW1tbWdf1kpdhGJWbmCBHR0dcW1vLfr+f5+fni47AWl9fL7pX13UmItZ13TFOPB6X8ug6EfkLBoNMRByNRj9cW14kau2V4uWCkFlc7uxTXVRV5Xg8zouLixwKhZiIeGFhoUKzcQcKQvAk0zR5bm6OW1pauK6ujjs7O3ltba3omDAbCkInEfmz8/TZlU6nKzMhwU5OTnhsbIwDgQA3NDRwb28vb29vO+777B+LZVm8tbXFPT09rCgKq6rK09PTnjyh5L2f5q/cuvLyj18i1l4pXi8ImcXlbn9/nwcGBtjv93N9fT1HIhFOJpMuR+8+H/O7F4UAAAAAQCr4qAQAAABAcigIAQAAACSHghAAAABAcigIAQAAACSHghAAAABAcigIAQAAACSHghAAAABAcigIAQAAACSHghAAAABAcigIAQAAACSHghAAAABAcigIAQAAACT3D5h5tElnIaskAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x640 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_forest(trace, var_names=[\"alphabeta\"], combined=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32850285",
   "metadata": {},
   "source": [
    "## Model 2 with corner constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "80dc9379",
   "metadata": {
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [mu0, _alpha, _beta, _alphabeta, tau]\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:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\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='12000' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [12000/12000 00:10&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 10 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model(coords=coords) as m:\n",
    "    apa_data = pm.Data(\"apa_data\", data.apa, mutable=False)\n",
    "    time_idx_data = pm.Data(\"time_idx_data\", time_idx, dims=\"id\", mutable=False)\n",
    "    conc_idx_data = pm.Data(\"conc_idx_data\", conc_idx, dims=\"id\", mutable=False)\n",
    "\n",
    "    mu0 = pm.Normal(\"mu0\", 0, tau=0.0001)\n",
    "    _alpha = pm.Normal(\"_alpha\", 0, tau=0.0001, dims=\"conc\")\n",
    "    _beta = pm.Normal(\"_beta\", 0, tau=0.0001, dims=\"time\")\n",
    "    _alphabeta = pm.Normal(\"_alphabeta\", 0, tau=0.0001, dims=(\"conc\", \"time\"))\n",
    "    tau = pm.Gamma(\"tau\", 0.001, 0.001)\n",
    "    sigma = pm.Deterministic(\"sigma\", 1 / tau**0.5)\n",
    "\n",
    "    # corner constraints: sets the first element of a dimension to zero\n",
    "    alpha = pm.Deterministic(\"alpha\", st.set_subtensor(_alpha[0], 0), dims=\"conc\")\n",
    "    beta = pm.Deterministic(\"beta\", st.set_subtensor(_beta[0], 0), dims=\"time\")\n",
    "    _alphabeta = st.set_subtensor(_alphabeta[:, 0], 0)\n",
    "    alphabeta = pm.Deterministic(\n",
    "        \"alphabeta\", st.set_subtensor(_alphabeta[0, :], 0), dims=(\"conc\", \"time\")\n",
    "    )\n",
    "\n",
    "    mu = (\n",
    "        mu0\n",
    "        + alpha[conc_idx_data]\n",
    "        + beta[time_idx_data]\n",
    "        + alphabeta[conc_idx_data, time_idx_data]\n",
    "    )\n",
    "    pm.Normal(\"apa\", mu, tau=tau, observed=apa_data, dims=\"id\")\n",
    "\n",
    "    differences(alpha, coords[\"conc\"])\n",
    "    differences(beta, coords[\"time\"])\n",
    "\n",
    "    trace = pm.sample(2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c77d43a5",
   "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>mu0</th>\n",
       "      <td>0.247</td>\n",
       "      <td>0.124</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.469</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>tau</th>\n",
       "      <td>17.120</td>\n",
       "      <td>2.884</td>\n",
       "      <td>11.492</td>\n",
       "      <td>22.201</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>0.244</td>\n",
       "      <td>0.021</td>\n",
       "      <td>0.206</td>\n",
       "      <td>0.282</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[1]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[2]</th>\n",
       "      <td>0.257</td>\n",
       "      <td>0.176</td>\n",
       "      <td>-0.086</td>\n",
       "      <td>0.573</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[3]</th>\n",
       "      <td>-0.082</td>\n",
       "      <td>0.173</td>\n",
       "      <td>-0.400</td>\n",
       "      <td>0.243</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[4]</th>\n",
       "      <td>-0.025</td>\n",
       "      <td>0.176</td>\n",
       "      <td>-0.355</td>\n",
       "      <td>0.302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2]</th>\n",
       "      <td>-0.191</td>\n",
       "      <td>0.174</td>\n",
       "      <td>-0.511</td>\n",
       "      <td>0.140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3]</th>\n",
       "      <td>-0.186</td>\n",
       "      <td>0.174</td>\n",
       "      <td>-0.513</td>\n",
       "      <td>0.145</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[4]</th>\n",
       "      <td>0.771</td>\n",
       "      <td>0.175</td>\n",
       "      <td>0.445</td>\n",
       "      <td>1.100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[5]</th>\n",
       "      <td>-0.001</td>\n",
       "      <td>0.175</td>\n",
       "      <td>-0.347</td>\n",
       "      <td>0.309</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[6]</th>\n",
       "      <td>-0.146</td>\n",
       "      <td>0.177</td>\n",
       "      <td>-0.472</td>\n",
       "      <td>0.195</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 1]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 2]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 3]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 4]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 5]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[1, 6]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 1]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 2]</th>\n",
       "      <td>-0.256</td>\n",
       "      <td>0.246</td>\n",
       "      <td>-0.708</td>\n",
       "      <td>0.212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 3]</th>\n",
       "      <td>0.324</td>\n",
       "      <td>0.248</td>\n",
       "      <td>-0.137</td>\n",
       "      <td>0.787</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 4]</th>\n",
       "      <td>-0.912</td>\n",
       "      <td>0.246</td>\n",
       "      <td>-1.384</td>\n",
       "      <td>-0.458</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 5]</th>\n",
       "      <td>-0.406</td>\n",
       "      <td>0.251</td>\n",
       "      <td>-0.880</td>\n",
       "      <td>0.057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[2, 6]</th>\n",
       "      <td>-0.180</td>\n",
       "      <td>0.250</td>\n",
       "      <td>-0.656</td>\n",
       "      <td>0.287</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 1]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 2]</th>\n",
       "      <td>0.096</td>\n",
       "      <td>0.246</td>\n",
       "      <td>-0.370</td>\n",
       "      <td>0.552</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 3]</th>\n",
       "      <td>0.120</td>\n",
       "      <td>0.244</td>\n",
       "      <td>-0.334</td>\n",
       "      <td>0.583</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 4]</th>\n",
       "      <td>-0.682</td>\n",
       "      <td>0.248</td>\n",
       "      <td>-1.140</td>\n",
       "      <td>-0.199</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 5]</th>\n",
       "      <td>0.138</td>\n",
       "      <td>0.245</td>\n",
       "      <td>-0.335</td>\n",
       "      <td>0.586</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[3, 6]</th>\n",
       "      <td>0.069</td>\n",
       "      <td>0.246</td>\n",
       "      <td>-0.395</td>\n",
       "      <td>0.527</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 1]</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 2]</th>\n",
       "      <td>0.141</td>\n",
       "      <td>0.249</td>\n",
       "      <td>-0.324</td>\n",
       "      <td>0.603</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 3]</th>\n",
       "      <td>0.045</td>\n",
       "      <td>0.247</td>\n",
       "      <td>-0.428</td>\n",
       "      <td>0.502</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 4]</th>\n",
       "      <td>-0.753</td>\n",
       "      <td>0.250</td>\n",
       "      <td>-1.230</td>\n",
       "      <td>-0.296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 5]</th>\n",
       "      <td>-0.002</td>\n",
       "      <td>0.249</td>\n",
       "      <td>-0.473</td>\n",
       "      <td>0.449</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alphabeta[4, 6]</th>\n",
       "      <td>0.166</td>\n",
       "      <td>0.248</td>\n",
       "      <td>-0.283</td>\n",
       "      <td>0.639</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[1] - alpha[2]</th>\n",
       "      <td>-0.257</td>\n",
       "      <td>0.176</td>\n",
       "      <td>-0.573</td>\n",
       "      <td>0.086</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[1] - alpha[3]</th>\n",
       "      <td>0.082</td>\n",
       "      <td>0.173</td>\n",
       "      <td>-0.243</td>\n",
       "      <td>0.400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[1] - alpha[4]</th>\n",
       "      <td>0.025</td>\n",
       "      <td>0.176</td>\n",
       "      <td>-0.302</td>\n",
       "      <td>0.355</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[2] - alpha[3]</th>\n",
       "      <td>0.338</td>\n",
       "      <td>0.176</td>\n",
       "      <td>0.013</td>\n",
       "      <td>0.672</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[2] - alpha[4]</th>\n",
       "      <td>0.282</td>\n",
       "      <td>0.175</td>\n",
       "      <td>-0.041</td>\n",
       "      <td>0.619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha[3] - alpha[4]</th>\n",
       "      <td>-0.057</td>\n",
       "      <td>0.175</td>\n",
       "      <td>-0.378</td>\n",
       "      <td>0.276</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[2]</th>\n",
       "      <td>0.191</td>\n",
       "      <td>0.174</td>\n",
       "      <td>-0.140</td>\n",
       "      <td>0.511</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[3]</th>\n",
       "      <td>0.186</td>\n",
       "      <td>0.174</td>\n",
       "      <td>-0.145</td>\n",
       "      <td>0.513</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[4]</th>\n",
       "      <td>-0.771</td>\n",
       "      <td>0.175</td>\n",
       "      <td>-1.100</td>\n",
       "      <td>-0.445</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[5]</th>\n",
       "      <td>0.001</td>\n",
       "      <td>0.175</td>\n",
       "      <td>-0.309</td>\n",
       "      <td>0.347</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1] - beta[6]</th>\n",
       "      <td>0.146</td>\n",
       "      <td>0.177</td>\n",
       "      <td>-0.195</td>\n",
       "      <td>0.472</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2] - beta[3]</th>\n",
       "      <td>-0.005</td>\n",
       "      <td>0.176</td>\n",
       "      <td>-0.334</td>\n",
       "      <td>0.323</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2] - beta[4]</th>\n",
       "      <td>-0.962</td>\n",
       "      <td>0.177</td>\n",
       "      <td>-1.297</td>\n",
       "      <td>-0.637</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2] - beta[5]</th>\n",
       "      <td>-0.190</td>\n",
       "      <td>0.175</td>\n",
       "      <td>-0.526</td>\n",
       "      <td>0.127</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2] - beta[6]</th>\n",
       "      <td>-0.045</td>\n",
       "      <td>0.176</td>\n",
       "      <td>-0.360</td>\n",
       "      <td>0.308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3] - beta[4]</th>\n",
       "      <td>-0.956</td>\n",
       "      <td>0.173</td>\n",
       "      <td>-1.296</td>\n",
       "      <td>-0.645</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3] - beta[5]</th>\n",
       "      <td>-0.185</td>\n",
       "      <td>0.172</td>\n",
       "      <td>-0.515</td>\n",
       "      <td>0.126</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3] - beta[6]</th>\n",
       "      <td>-0.040</td>\n",
       "      <td>0.173</td>\n",
       "      <td>-0.363</td>\n",
       "      <td>0.285</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[4] - beta[5]</th>\n",
       "      <td>0.771</td>\n",
       "      <td>0.175</td>\n",
       "      <td>0.458</td>\n",
       "      <td>1.103</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[4] - beta[6]</th>\n",
       "      <td>0.916</td>\n",
       "      <td>0.177</td>\n",
       "      <td>0.600</td>\n",
       "      <td>1.264</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[5] - beta[6]</th>\n",
       "      <td>0.145</td>\n",
       "      <td>0.176</td>\n",
       "      <td>-0.184</td>\n",
       "      <td>0.474</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       mean     sd  hdi_3%  hdi_97%\n",
       "mu0                   0.247  0.124   0.007    0.469\n",
       "tau                  17.120  2.884  11.492   22.201\n",
       "sigma                 0.244  0.021   0.206    0.282\n",
       "alpha[1]              0.000  0.000   0.000    0.000\n",
       "alpha[2]              0.257  0.176  -0.086    0.573\n",
       "alpha[3]             -0.082  0.173  -0.400    0.243\n",
       "alpha[4]             -0.025  0.176  -0.355    0.302\n",
       "beta[1]               0.000  0.000   0.000    0.000\n",
       "beta[2]              -0.191  0.174  -0.511    0.140\n",
       "beta[3]              -0.186  0.174  -0.513    0.145\n",
       "beta[4]               0.771  0.175   0.445    1.100\n",
       "beta[5]              -0.001  0.175  -0.347    0.309\n",
       "beta[6]              -0.146  0.177  -0.472    0.195\n",
       "alphabeta[1, 1]       0.000  0.000   0.000    0.000\n",
       "alphabeta[1, 2]       0.000  0.000   0.000    0.000\n",
       "alphabeta[1, 3]       0.000  0.000   0.000    0.000\n",
       "alphabeta[1, 4]       0.000  0.000   0.000    0.000\n",
       "alphabeta[1, 5]       0.000  0.000   0.000    0.000\n",
       "alphabeta[1, 6]       0.000  0.000   0.000    0.000\n",
       "alphabeta[2, 1]       0.000  0.000   0.000    0.000\n",
       "alphabeta[2, 2]      -0.256  0.246  -0.708    0.212\n",
       "alphabeta[2, 3]       0.324  0.248  -0.137    0.787\n",
       "alphabeta[2, 4]      -0.912  0.246  -1.384   -0.458\n",
       "alphabeta[2, 5]      -0.406  0.251  -0.880    0.057\n",
       "alphabeta[2, 6]      -0.180  0.250  -0.656    0.287\n",
       "alphabeta[3, 1]       0.000  0.000   0.000    0.000\n",
       "alphabeta[3, 2]       0.096  0.246  -0.370    0.552\n",
       "alphabeta[3, 3]       0.120  0.244  -0.334    0.583\n",
       "alphabeta[3, 4]      -0.682  0.248  -1.140   -0.199\n",
       "alphabeta[3, 5]       0.138  0.245  -0.335    0.586\n",
       "alphabeta[3, 6]       0.069  0.246  -0.395    0.527\n",
       "alphabeta[4, 1]       0.000  0.000   0.000    0.000\n",
       "alphabeta[4, 2]       0.141  0.249  -0.324    0.603\n",
       "alphabeta[4, 3]       0.045  0.247  -0.428    0.502\n",
       "alphabeta[4, 4]      -0.753  0.250  -1.230   -0.296\n",
       "alphabeta[4, 5]      -0.002  0.249  -0.473    0.449\n",
       "alphabeta[4, 6]       0.166  0.248  -0.283    0.639\n",
       "alpha[1] - alpha[2]  -0.257  0.176  -0.573    0.086\n",
       "alpha[1] - alpha[3]   0.082  0.173  -0.243    0.400\n",
       "alpha[1] - alpha[4]   0.025  0.176  -0.302    0.355\n",
       "alpha[2] - alpha[3]   0.338  0.176   0.013    0.672\n",
       "alpha[2] - alpha[4]   0.282  0.175  -0.041    0.619\n",
       "alpha[3] - alpha[4]  -0.057  0.175  -0.378    0.276\n",
       "beta[1] - beta[2]     0.191  0.174  -0.140    0.511\n",
       "beta[1] - beta[3]     0.186  0.174  -0.145    0.513\n",
       "beta[1] - beta[4]    -0.771  0.175  -1.100   -0.445\n",
       "beta[1] - beta[5]     0.001  0.175  -0.309    0.347\n",
       "beta[1] - beta[6]     0.146  0.177  -0.195    0.472\n",
       "beta[2] - beta[3]    -0.005  0.176  -0.334    0.323\n",
       "beta[2] - beta[4]    -0.962  0.177  -1.297   -0.637\n",
       "beta[2] - beta[5]    -0.190  0.175  -0.526    0.127\n",
       "beta[2] - beta[6]    -0.045  0.176  -0.360    0.308\n",
       "beta[3] - beta[4]    -0.956  0.173  -1.296   -0.645\n",
       "beta[3] - beta[5]    -0.185  0.172  -0.515    0.126\n",
       "beta[3] - beta[6]    -0.040  0.173  -0.363    0.285\n",
       "beta[4] - beta[5]     0.771  0.175   0.458    1.103\n",
       "beta[4] - beta[6]     0.916  0.177   0.600    1.264\n",
       "beta[5] - beta[6]     0.145  0.176  -0.184    0.474"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace, var_names=\"~_\", filter_vars=\"like\", kind=\"stats\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8c1ff385",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: Wed Mar 22 2023\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.11.0\n",
      "IPython version      : 8.9.0\n",
      "\n",
      "aesara: 2.8.10\n",
      "aeppl : 0.1.1\n",
      "\n",
      "numpy   : 1.24.2\n",
      "pandas  : 1.5.3\n",
      "pytensor: 2.10.1\n",
      "arviz   : 0.14.0\n",
      "pymc    : 5.1.2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv -p aesara,aeppl"
   ]
  }
 ],
 "metadata": {
  "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.11.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "e197428f119775a30e0221ede525e07580bbbcb52f3c1ab01042e9594a2688a6"
   },
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}