Skip to content

Instantly share code, notes, and snippets.

@GStechschulte

GStechschulte/plot_comparisons_uncertainty_calculation.ipynb

Created June 9, 2023 15:07
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save GStechschulte/351e38d498d600c7f46423afb588634b to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save GStechschulte/351e38d498d600c7f46423afb588634b to your computer and use it in GitHub Desktop.
Download ZIP
plot_comparisons uncertainty in difference calculation
Display the source blob
Display the rendered blob
Raw
plot_comparisons_uncertainty_calculation.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The autoreload extension is already loaded. To reload it, use:\n",
" %reload_ext autoreload\n"
]
}
],
"source": [
"import arviz as az\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"import pandas as pd\n",
"import xarray as xr\n",
"\n",
"import bambi as bmb\n",
"from bambi.plots import plot_comparison\n",
"\n",
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### mtcars"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Load data\n",
"data = bmb.load_data('mtcars')\n",
"data[\"cyl\"] = data[\"cyl\"].replace({4: \"low\", 6: \"medium\", 8: \"high\"})\n",
"data[\"gear\"] = data[\"gear\"].replace({3: \"A\", 4: \"B\", 5: \"C\"})\n",
"data[\"cyl\"] = pd.Categorical(data[\"cyl\"], categories=[\"low\", \"medium\", \"high\"], ordered=True)\n",
"data[\"am\"] = pd.Categorical(data[\"am\"], categories=[0, 1], ordered=True)\n",
"#data[\"drat\"] = pd.Categorical(data[\"drat\"], ordered=True).codes"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Initializing NUTS using jitter+adapt_diag...\n",
"Multiprocess sampling (4 chains in 4 jobs)\n",
"NUTS: [mpg_sigma, Intercept, hp, drat, hp:drat, am, hp:am, drat:am, hp:drat:am]\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 23 seconds.\n"
]
}
],
"source": [
"# Define and fit the Bambi model\n",
"mt_model = bmb.Model(\"mpg ~ hp * drat * am\", data)\n",
"mt_idata = mt_model.fit(draws=1000, target_accept=0.95, random_seed=1234)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnoAAAFMCAYAAABcYREhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABh/0lEQVR4nO3deVxU9foH8M8Aww4ji2yCgqgsAu4LuKCJiIpLZWoaqTdNTXOrX2pmqLliWeauuZVrqZQmkmhCmSAoYqKmZiimIK6AKNvw/f3hZa7jDDAi6/B539e8bvOd55zznJkRHs53ORIhhAARERERaR2d6k6AiIiIiCoHCz0iIiIiLcVCj4iIiEhLsdAjIiIi0lIs9IiIiIi0FAs9IiIiIi3FQo+IiIhIS7HQIyIiItJSLPSIiIiItBQLPaoztmzZAolEongYGhrCzs4O3bt3x6JFi5CRkaGyzZw5cyCRSJTa8vPzMW7cONjb20NXVxctW7YEANy/fx9Dhw6FjY0NJBIJBg4cWAVnVXM5Oztj5MiR1Z1GtXB2dkZwcHCF7nPhwoX48ccfK3SfVD0kEgnmzJlTrm2jo6MhkUgQHR1doTlFRESUOyeq2fSqOwGiqrZ582a4u7ujoKAAGRkZOH78OJYsWYLPP/8cu3fvRkBAgCJ29OjRCAoKUtp+zZo1WLduHVasWIE2bdrA1NQUAPDZZ58hPDwcmzZtgqurKywtLav0vGqa8PBwmJubV3caWmPhwoUYNGhQnf8DgipHREQEVq1axWJPC7HQozrHy8sLbdu2VTx//fXXMXXqVHTu3BmvvfYarly5AltbWwCAo6MjHB0dlbZPTk6GkZERJk6cqNLu6uqK4cOHV1iuT548gZGRUYXtryq1atWqulMgqrMeP34MY2Pj6k6DagB23RIBaNiwIb744gtkZ2dj3bp1ivbnu24lEgm++eYbPHnyRNEFXNwlfOTIEVy8eFHRXty1kp+fj/nz58Pd3R0GBgaoX78+Ro0ahTt37ijlUNzdt2/fPrRq1QqGhoaYO3cuACA9PR1jx46Fo6Mj9PX14eLigrlz56KwsFCx/bVr1yCRSPD5559j2bJlcHFxgampKXx9fREXF6dyzidPnkS/fv1gZWUFQ0NDuLq6YsqUKUoxV65cwbBhw2BjYwMDAwN4eHhg1apVGr2nz3fdFnc57dy5E7NmzYKDgwPMzc0REBCAS5culbm/4s/izz//xBtvvAGZTAZLS0tMmzYNhYWFuHTpEoKCgmBmZgZnZ2eEhYUpbV98/G3btmHatGmws7ODkZER/P39cebMGZXjbdiwAc2aNYOBgQE8PT2xY8cOjBw5Es7OzhqdP/D0qqaPjw8MDQ3RuHFjfP311yoxWVlZ+PDDD+Hi4gJ9fX00aNAAU6ZMQU5OjiJGIpEgJycHW7duVXy/unXrhqysLOjp6WHp0qWK2Lt370JHRwcymUzp+zFp0iTUr18fQghF25EjR9CjRw+Ym5vD2NgYnTp1wtGjR1Vy1OR78LKf799//41Ro0ahadOmMDY2RoMGDdCvXz+cO3dO7XF27NiB6dOnw97eHqampujXrx9u376N7OxsvPvuu7C2toa1tTVGjRqFR48elXn8qKgoDBgwAI6OjjA0NESTJk0wduxY3L17Vymu+Ht4/vx5vPnmm5DJZLC1tcV//vMfZGZmKsVmZWVhzJgxsLKygqmpKYKCgnD58uUycyn2119/ISgoCMbGxrC2tsa4ceOQnZ2tEtetWzd4eXnht99+g5+fH4yNjfGf//wHALB7924EBgbC3t4eRkZG8PDwwIwZM5S+XyNHjlR8ns8Ob7l27ZrGuVINJojqiM2bNwsAIiEhQe3rjx49Erq6uqJHjx6KttDQUPHsP5PY2FjRp08fYWRkJGJjY0VsbKxIT08XsbGxolWrVqJx48aK9szMTCGXy0VQUJAwMTERc+fOFVFRUeKbb74RDRo0EJ6enuLx48eKfTdq1EjY29uLxo0bi02bNoljx46J+Ph4kZaWJpycnESjRo3EunXrxJEjR8Rnn30mDAwMxMiRIxXbp6SkCADC2dlZBAUFiR9//FH8+OOPwtvbW1hYWIiHDx8qYiMjI4VUKhU+Pj5iy5Yt4tdffxWbNm0SQ4cOVcScP39eyGQy4e3tLb799ltx+PBh8cEHHwgdHR0xZ86cMt/vRo0aiREjRiieHzt2TJHf8OHDxcGDB8XOnTtFw4YNRdOmTUVhYWGp+yv+LNzc3MRnn30moqKixEcffSQAiIkTJwp3d3fx9ddfi6ioKDFq1CgBQOzdu1fl+E5OTmLAgAHiwIEDYtu2baJJkybC3NxcXL16VRG7bt06AUC8/vrr4ueffxbbt28XzZo1E40aNRKNGjXS6NwbNGggGjZsKDZt2iQiIiLE8OHDBQCxdOlSRVxOTo5o2bKlsLa2FsuWLRNHjhwRy5cvFzKZTLzyyiuiqKhICPH0e2dkZCT69Omj+H6dP39eCCFEx44dRWBgoGKfu3btEoaGhkIikYg//vhD0e7h4SEGDx6seP7dd98JiUQiBg4cKPbt2ycOHDgggoODha6urjhy5IgiTtPvwct+vjExMeKDDz4Qe/bsETExMSI8PFwMHDhQGBkZib/++kvlOI0aNRIjR44UkZGRYu3atcLU1FR0795d9OzZU3z44Yfi8OHDYsmSJUJXV1e8//77ZX5ma9asEYsWLRL79+8XMTExYuvWraJFixbCzc1N5OfnK+Ke/R5++umnIioqSixbtkwYGBiIUaNGKeKKiopE9+7dhYGBgViwYIE4fPiwCA0NFY0bNxYARGhoaKn5pKenCxsbG9GgQQOxefNmxXeoYcOGAoA4duyYItbf319YWloKJycnsWLFCnHs2DERExMjhBDis88+E19++aU4ePCgiI6OFmvXrhUuLi6ie/fuiu3//vtvMWjQIAFA8f2KjY0Vubm5Zb5vVPOx0KM6o6xCTwghbG1thYeHh+L584WeEEKMGDFCmJiYqGzr7+8vmjdvrtS2c+dOlYJDCCESEhIEALF69WpFW6NGjYSurq64dOmSUuzYsWOFqampuH79ulL7559/LgAofuEXF3re3t5Kv1Tj4+MFALFz505Fm6urq3B1dRVPnjwp8b3o1auXcHR0FJmZmUrtEydOFIaGhuL+/fslblt8PuoKvT59+ijFff/994pfMKUp/iy++OILpfaWLVsKAGLfvn2KtoKCAlG/fn3x2muvqRy/devWigJKCCGuXbsmpFKpGD16tBBCCLlcLuzs7ESHDh2UjnP9+nUhlUo1LvQkEolISkpSau/Zs6cwNzcXOTk5QgghFi1aJHR0dFS+k3v27BEAREREhKLNxMRE6f0s9sknnwgjIyPFL+XRo0eLoKAg4ePjI+bOnSuEEOLmzZsCgFi/fr0Q4mmBaWlpKfr166e0L7lcLlq0aCHat2+vaNP0e/Cyn+/zCgsLRX5+vmjatKmYOnWqor34OM/nPmXKFAFATJo0Sal94MCBwtLS8oWOXVRUJAoKCsT169cFAPHTTz8pXiv+HoaFhSlt89577wlDQ0PFd+vQoUMCgFi+fLlS3IIFCzQq9KZPn17id0hdoQdAHD16VKPziomJEQDE2bNnFa9NmDBB5WcdaQd23RI9QzzTrVURfv75Z9SrVw/9+vVDYWGh4tGyZUvY2dmpzJzz8fFBs2bNVPbRvXt3ODg4KO2jd+/eAICYmBil+L59+0JXV1dpnwBw/fp1AMDly5dx9epVvPPOOzA0NFSbd25uLo4ePYpXX30VxsbGSsft06cPcnNz1XYHa6J///4q5/xsfmV5fjarh4cHJBKJ4v0AAD09PTRp0kTtPocNG6bUHd+oUSP4+fnh2LFjAIBLly4hPT0dgwcPVtquYcOG6NSpk0Y5AkDz5s3RokULlWNnZWUhMTERwNPP1svLCy1btlR6j3v16qXxzMoePXrgyZMnOHHiBICn3bE9e/ZEQEAAoqKiFG0AFBONTpw4gfv372PEiBFKxy0qKkJQUBASEhKQk5NTru9BeT/fwsJCLFy4EJ6entDX14eenh709fVx5coVXLx4USVe3fcAePr9f779/v37ZXbfZmRkYNy4cXBycoKenh6kUikaNWoEAGqPr+48c3NzFbP3i79Pz4/ZHTZsWKl5FDt27FiJ3yF1LCws8Morr6i0//PPPxg2bBjs7Oygq6sLqVQKf3//Es+LtA8nYxD9V05ODu7duwdvb+8K2+ft27fx8OFD6Ovrq339+fE/9vb2avdx4MABSKVSjfZhZWWl9NzAwADA04kdABRjA5+fZPKse/fuobCwECtWrMCKFSs0Oq6mysqvLM/PZtbX14exsbFK0aqvr4+srCyV7e3s7NS2nT17FsDTcwegmJDzLFtbW6SkpGiUZ0nHefYYt2/fxt9//63xZ6tO8ZisI0eOwMnJCdeuXUPPnj3x77//YsWKFXj06BGOHDmCxo0bw8XFRXFcABg0aFCJ+71//z50dHRe+HtQ3s932rRpWLVqFaZPnw5/f39YWFhAR0cHo0ePVrutuu9Bae25ubmKGfLPKyoqQmBgIG7duoXZs2fD29sbJiYmKCoqQseOHdUev6zzvHfvHvT09FTi1H0v1Ll3757i89Jke3U/Ox49eoQuXbrA0NAQ8+fPR7NmzWBsbIwbN27gtdde0/jfHNVuLPSI/uvgwYOQy+Xo1q1bhe3T2toaVlZWiIyMVPu6mZmZ0vPn1+wr3oePjw8WLFigdh8ODg4vlFP9+vUBAP/++2+JMRYWFtDV1UVISAgmTJigNkbdL6HaID09XW1b8S/k4v8vLobK2vZFj/PsMaytrWFkZIRNmzap3Ye1tXWZx9HX10fnzp1x5MgRODo6ws7ODt7e3mjcuDGAp5MXjh49qnQFrHi/K1asQMeOHdXu19bWFoWFhVX2Pdi2bRvefvttLFy4UKn97t27qFevXoUcoyTJyck4e/YstmzZghEjRija//7773Lv08rKCoWFhbh3755Ssafpd8jKyqrU79Dz1P3s+PXXX3Hr1i1ER0crruIBwMOHDzXKgbQDCz0iAKmpqfjwww8hk8kwduzYCttvcHAwdu3aBblcjg4dOpR7HxEREXB1dYWFhcVL59SsWTO4urpi06ZNmDZtmuJKxLOMjY3RvXt3nDlzBj4+PiVekayNdu7ciWnTpil+MV6/fh0nTpzA22+/DQBwc3ODnZ0dvv/+e0ybNk2xXWpqKk6cOKFxYX3+/HmcPXtWqettx44dMDMzQ+vWrQE8/WwXLlwIKyurMgsmAwODEq/ABAQEYObMmTAzM1N0z5qYmKBjx45YsWIFbt26pbQ+ZKdOnVCvXj1cuHBBZZmgZ+nr61fZ90Aikah8Fw8ePIibN2+iSZMmlXbc4mMDUDn+szPwX1T37t0RFhaG7du3Y9KkSYr2HTt2vND26r5DmnqR83r2imRtXdKJ1GOhR3VOcnKyYpxRRkYGfv/9d2zevBm6uroIDw9XXPGqCEOHDsX27dvRp08fTJ48Ge3bt4dUKsW///6LY8eOYcCAAXj11VdL3ce8efMQFRUFPz8/TJo0CW5ubsjNzcW1a9cQERGBtWvXltoNq86qVavQr18/dOzYEVOnTkXDhg2RmpqKX375Bdu3bwcALF++HJ07d0aXLl0wfvx4ODs7Izs7G3///TcOHDiAX3/9tdzvS3XKyMjAq6++ijFjxiAzMxOhoaEwNDTEzJkzAQA6OjqYO3cuxo4di0GDBuE///kPHj58iLlz58Le3h46OpoNbXZwcED//v0xZ84c2NvbY9u2bYiKisKSJUsU65tNmTIFe/fuRdeuXTF16lT4+PigqKgIqampOHz4MD744APFHwje3t6Ijo7GgQMHYG9vDzMzM7i5uQF4Ok5PLpfj6NGj2Lp1qyKHgIAAhIaGQiKRKI3fMjU1xYoVKzBixAjcv38fgwYNgo2NDe7cuYOzZ8/izp07WLNmDYCq+x4EBwdjy5YtcHd3h4+PD06fPo2lS5e+8He7PNzd3eHq6ooZM2ZACAFLS0scOHBAMcaxPAIDA9G1a1d89NFHyMnJQdu2bfHHH3/gu+++02j7KVOmYNOmTejbty/mz58PW1tbbN++HX/99ZfGOfj5+cHCwgLjxo1DaGgopFIptm/frhim8KziIStLlixB7969oaurq3V/5NVVLPSozhk1ahSAp1cr6tWrBw8PD0yfPh2jR4+u0CIPAHR1dbF//34sX74c3333HRYtWgQ9PT04OjrC399fo/GA9vb2OHXqFD777DMsXboU//77L8zMzODi4oKgoKByXeXr1asXfvvtN8ybNw+TJk1Cbm4uHB0dlQaYe3p6IjExEZ999hk++eQTZGRkoF69emjatCn69OnzwsesKRYuXIiEhASMGjUKWVlZaN++PXbt2gVXV1dFzLvvvguJRIKwsDC8+uqrcHZ2xowZM/DTTz8hNTVVo+O0bNkSo0aNQmhoKK5cuQIHBwcsW7YMU6dOVcSYmJjg999/x+LFi7F+/XqkpKTAyMgIDRs2REBAgNKafcuXL8eECRMwdOhQPH78GP7+/orJGq1atYK1tTXu3r2rdOWuuNBr1aqVylixt956Cw0bNkRYWBjGjh2L7Oxs2NjYoGXLlkrrH1bV92D58uWQSqVYtGgRHj16hNatW2Pfvn345JNPKuwYJZFKpThw4AAmT56MsWPHQk9PDwEBAThy5AgaNmxYrn3q6Ohg//79mDZtGsLCwpCfn49OnTohIiIC7u7uZW5vZ2eHmJgYTJ48GePHj4exsTFeffVVrFy5EgMGDNAoBysrKxw8eBAffPAB3nrrLZiYmGDAgAHYvXu34qpysWHDhuGPP/7A6tWrMW/ePAghkJKS8kLrRlLNJBEVPc2QiKgGio6ORvfu3fHDDz+UOgmhJA8fPkSzZs0wcOBArF+/vhIyJCKqeLyiR0T0nPT0dCxYsADdu3eHlZUVrl+/ji+//BLZ2dmYPHlydadHRKQxFnpERM8xMDDAtWvX8N577+H+/fswNjZGx44dsXbtWjRv3ry60yMi0hi7bomIiIi0FO+MQURERKSlWOgRaandu3ejefPmMDIygkQiQVJSUqUda8uWLZBIJDh16lSlHaM63bp1C3PmzKnU97AsO3bswFdffaVxvEQiUTw+//xzpdc++eQTBAcHo0GDBpBIJEqzbJ91/vx5vPfee/D19YWJiUmZt2XbtWsXWrZsCUNDQzg4OGDKlCll3nqsNMePH8fo0aPRpk0bGBgYQCKR4Nq1aypxxd+/kh6LFy/WKPb5xYhbtmypeO35W64R1RYs9Ii00J07dxASEgJXV1dERkYiNjZW5R66pLlbt25h7ty5tarQA4B33nkHsbGxKvdb/fLLL3Hv3j3079+/1HXSTp06hR9//BGWlpbo0aNHqcfavn073nzzTbRr1w6HDh1CaGgotmzZgtdee+2Fcn7W0aNHFUuc+Pn5lRjXt29fxMbGqjx69uwJAGrXqty8ebNK/PNL0Hz33XeIjY3V+LZlRDURJ2MQaaHLly+joKAAb731ltKtj17G48ePFQv9Uulqynvl6Oio9hZn2dnZioWfS1vANyQkRHFLsD179uDAgQNq4+RyOf7v//4PgYGB2LBhA4Cnd3YwMzPD8OHDcejQIfTu3fuF8589ezZCQ0MBAJ9//nmJVxPr16+vsgZmTk4OYmNj0blzZ8XC0s/y8vJC27ZtSz1+8TqX6u4eQ1Rb8IoekZYZOXIkOnfuDAAYMmQIJBKJ0v179+/fD19fXxgbG8PMzAw9e/ZEbGys0j7mzJkDiUSCxMREDBo0CBYWFkoLCpckOzsb48ePV9zj97XXXsOtW7eUYpydnREcHIzw8HD4+PjA0NAQjRs3xtdff63R+RUVFWHFihVo2bIljIyMUK9ePXTs2BH79+9XigkLC4O7uzsMDAxgY2ODt99+W+X+vt26dYOXlxcSEhLQpUsXGBsbo3Hjxli8eDGKiooAPF1/r127dgCeLrZd3JU3Z84cxfttamqKc+fOITAwEGZmZoqrX1FRURgwYAAcHR1haGiIJk2aYOzYsbh7965SHnfu3MG7774LJycnGBgYoH79+ujUqROOHDmiyPPgwYO4fv26UldjeWl6dw9N4+Li4pCWlqZYjLzYG2+8AVNTU4SHh79wji9yfHV2796NR48eYfTo0eXeB5E2YKFHpGVmz56NVatWAXh6F4jY2FisXr0awNPuvwEDBsDc3Bw7d+7Exo0b8eDBA3Tr1g3Hjx9X2ddrr72GJk2a4IcffsDatWvLPPbo0aMhlUqxY8cOhIWFITo6Gm+99ZZKXFJSEqZMmYKpU6ciPDwcfn5+mDx5sspYMnVGjhyJyZMno127dti9ezd27dqF/v37K43dGj9+PKZPn46ePXti//79+OyzzxAZGQk/Pz+VIis9PR3Dhw/HW2+9hf3796N3796YOXMmtm3bBgBo3bo1Nm/eDODp2Lbibr5nC4j8/Hz0798fr7zyCn766SfMnTsXAHD16lX4+vpizZo1OHz4MD799FOcPHkSnTt3RkFBgWL7kJAQ/Pjjj/j0009x+PBhfPPNNwgICMC9e/cAAKtXr0anTp1gZ2en1NVYUyQnJwMAfHx8lNqlUinc3d0Vr1eljRs3wtzcHG+88Yba14ODg6GrqwtLS0u89tpr1ZIjUZUQRKR1jh07JgCIH374QdEml8uFg4OD8Pb2FnK5XNGenZ0tbGxshJ+fn6ItNDRUABCffvqpRsfbvHmzACDee+89pfawsDABQKSlpSnaGjVqJCQSiUhKSlKK7dmzpzA3Nxc5OTklHue3334TAMSsWbNKjLl48aLaXE6ePCkAiI8//ljR5u/vLwCIkydPKsV6enqKXr16KZ4nJCQIAGLz5s0qxxsxYoQAIDZt2lRiTkIIUVRUJAoKCsT169cFAPHTTz8pXjM1NRVTpkwpdfu+ffuKRo0alRrzLAAiNDS0zDgTExMxYsSIMuN++OEHAUAcO3ZM5bUFCxaofM7FAgMDRbNmzTTIuHRLly4VAERKSkqZscXfgbFjx6q8dujQITFr1ixx4MABERMTI1auXCkcHR2FiYmJyneyWKNGjUTfvn1f9hSIqgWv6BHVEZcuXcKtW7cQEhKi1CVmamqK119/HXFxcXj8+LHSNq+//voLHePZe+UC/7vCc/36daX25s2bo0WLFkptw4YNQ1ZWFhITE0vc/6FDhwAAEyZMKDHm2LFjAKAyk7R9+/bw8PDA0aNHldrt7OzQvn17lbyfz7ks6t6rjIwMjBs3Dk5OTtDT04NUKkWjRo0AABcvXlTKbcuWLZg/fz7i4uKUrvbVJiV1J79MN3N5bNy4EQDUdtsGBQVh/vz5CA4ORteuXTFhwgT8/vvvkEgk+PTTT6s0T6KqwEKPqI4o7ga0t7dXec3BwQFFRUV48OCBUru62NI8P2uxeBD7kydPlNrVzWIsbivOU507d+5AV1e31FmQZZ3n8/t/PufivJ/PuTTGxsYwNzdXaisqKkJgYCD27duHjz76CEePHkV8fDzi4uIAKL8nu3fvxogRI/DNN9/A19cXlpaWePvtt1WW+6ipit9DdZ/d/fv3YWlpWWW5FBQU4Ntvv0WLFi3KnGxRzNnZGZ07d1Z8NkTahIUeUR1R/Ms4LS1N5bVbt25BR0cHFhYWSu2VdSVGXQFT3Kau8CpWv359yOXyUgugss7T2tr6RdMtk7r3KTk5GWfPnsXSpUvx/vvvo1u3bmjXrp3a87O2tsZXX32Fa9eu4fr161i0aBH27dtX4vp2NU3x7NRz584ptRcWFuKvv/6Cl5dXleXy888/IyMj44UnYQghXmryB1FNxW81UR3h5uaGBg0aYMeOHRDP3PkwJycHe/fuVczErQrnz5/H2bNnldp27NgBMzMztG7dusTtipfoWLNmTYkxr7zyCgAoJlMUS0hIwMWLF8tcD06dkq5Mlqa4+Ht+aY5169aVul3Dhg0xceJE9OzZU6kb+0WvMlalDh06wN7eHlu2bFFq37NnDx49evRSa+m9qI0bN8LQ0FBl7cDSpKSk4I8//lC7FA1Rbcd19IjqCB0dHYSFhWH48OEIDg7G2LFjkZeXh6VLl+Lhw4dKdw+obA4ODujfvz/mzJkDe3t7bNu2DVFRUViyZEmpxWaXLl0QEhKC+fPn4/bt2wgODoaBgQHOnDkDY2NjvP/++3Bzc8O7776LFStWQEdHB71798a1a9cwe/ZsODk5YerUqS+cr6urK4yMjLB9+3Z4eHjA1NQUDg4OcHBwKHEbd3d3uLq6YsaMGRBCwNLSEgcOHEBUVJRSXGZmJrp3745hw4bB3d0dZmZmSEhIQGRkpFKB5O3tjX379mHNmjVo06YNdHR0NO6afF5MTAzu3LkD4OkaeNevX8eePXsAAP7+/oo16R4/foyIiAgAUHRrxsTE4O7duzAxMVEU3rq6uggLC0NISAjGjh2LN998E1euXMFHH32Enj17IigoSOn4EokE/v7+pd5lA3jaVR8TEwPgf1cLDx06pFg37/k1Im/duoXIyEgMGTJE5ep0sYCAAHTt2hU+Pj4wNzfHuXPnEBYWBolEgs8++0yj94+oVqnmySBEVAnUzbot9uOPP4oOHToIQ0NDYWJiInr06CH++OMPpZjiWbd37tzR6HjFs24TEhLU5vHsTM3iGYx79uwRzZs3F/r6+sLZ2VksW7ZMo2PJ5XLx5ZdfCi8vL6Gvry9kMpnw9fUVBw4cUIpZsmSJaNasmZBKpcLa2lq89dZb4saNG0r78vf3F82bN1c5xogRI1RmuO7cuVO4u7sLqVSqNKN1xIgRwsTERG2uFy5cED179hRmZmbCwsJCvPHGGyI1NVVp+9zcXDFu3Djh4+MjzM3NhZGRkXBzcxOhoaFKM5Dv378vBg0aJOrVqyckEoko68c3Spl1WzzbWN3j2c8qJSWlxDh1M4B37NghfHx8hL6+vrCzsxOTJk0S2dnZSjHZ2dkCgBg6dGip+Qvxv++Puoe/v79KfPHs319//bXEfU6ZMkV4enoKMzMzoaenJxwcHMRbb70lLl26VOI2nHVLtZlEiGf6cIiIKpmzszO8vLzw888/V3cqWk0ikWD27Nn49NNPoaurW+UzX0sSERGB4OBgnD17VjG2r6aSy+UQQqBJkyb8zlKtxTF6RERa6rPPPoNUKsUXX3xR3akoHDt2DEOHDq3xRR4AtGnTBlKp9IWX2iGqSThGj4hICyUkJCj+28nJqRozUbZ06dLqTkFjO3bsUKwtWa9evepNhqic2HVLREREpKXYdUtERESkpVjoEREREWkpFnpEREREWoqFHhEREZGW4qzbClBUVIRbt27BzMysxqxVRURERNpJCIHs7Gw4ODiUeY9mFnoV4NatWzVq+QIiIiLSfjdu3ICjo2OpMSz0KoCZmRmAp2+4ubl5NWdDRERE2iwrKwtOTk6K+qM0LPQqQHF3rbm5OQs9IiIiqhKaDBfjZAwiIiIiLcVCj4iIiEhLsdAjIiIi0lIs9IiIiIi0FAs9IiIiIi3FQo+IiIhIS7HQIyIiItJSLPSIiEhjj/ML4TzjIJxnHMTj/MLqToeIysBCj4iIiEhLsdAjIiIi0lIs9IiIiIi0FAs9IiIiIi3FQo+IiIhIS7HQIyIiItJSLPSIiIiItBQLPSIiIiItxUKPiIiISEux0CMiIiLSUiz0iIiIiLRUrSv0Vq9eDRcXFxgaGqJNmzb4/fffS42PiYlBmzZtYGhoiMaNG2Pt2rUqMXv37oWnpycMDAzg6emJ8PDwykqfiIiIqMrUqkJv9+7dmDJlCmbNmoUzZ86gS5cu6N27N1JTU9XGp6SkoE+fPujSpQvOnDmDjz/+GJMmTcLevXsVMbGxsRgyZAhCQkJw9uxZhISEYPDgwTh58mRVnRYRERFRpZAIIUR1J6GpDh06oHXr1lizZo2izcPDAwMHDsSiRYtU4qdPn479+/fj4sWLirZx48bh7NmziI2NBQAMGTIEWVlZOHTokCImKCgIFhYW2Llzp0Z5ZWVlQSaTITMzE+bm5uU9PSKiGu9xfiE8P/0FAHBhXi8Y6+tVc0ZEdc+L1B215opefn4+Tp8+jcDAQKX2wMBAnDhxQu02sbGxKvG9evXCqVOnUFBQUGpMSfsEgLy8PGRlZSk9iIiIiGqaWlPo3b17F3K5HLa2tkrttra2SE9PV7tNenq62vjCwkLcvXu31JiS9gkAixYtgkwmUzycnJzKc0pERERElarWFHrFJBKJ0nMhhEpbWfHPt7/oPmfOnInMzEzF48aNGxrnT0RERFRVas3gCmtra+jq6qpcacvIyFC5IlfMzs5Obbyenh6srKxKjSlpnwBgYGAAAwOD8pwGEVGtZqyvh2uL+1Z3GkSkoVpzRU9fXx9t2rRBVFSUUntUVBT8/PzUbuPr66sSf/jwYbRt2xZSqbTUmJL2SURERFRb1JoregAwbdo0hISEoG3btvD19cX69euRmpqKcePGAXjapXrz5k18++23AJ7OsF25ciWmTZuGMWPGIDY2Fhs3blSaTTt58mR07doVS5YswYABA/DTTz/hyJEjOH78eLWcozqc5UZERETlUasqhiFDhuDevXuYN28e0tLS4OXlhYiICDRq1AgAkJaWprSmnouLCyIiIjB16lSsWrUKDg4O+Prrr/H6668rYvz8/LBr1y588sknmD17NlxdXbF792506NChys+PiIiIqCLVqnX0aqrKXkePV/SIiIiomFauo0dEREREL4aFHhEREZGWYqFHREREpKVY6BERERFpKRZ6RERERFqKhR4RERGRlmKhR0RERKSlWOgREVHdk58DzJE9feTnVHc2RJWGhR4RERGRlmKhR0RERKSlWOgRERERaSkWekRERERaioUeERERkZbSq+4EiIiIqpy+CTAns7qzIKp0vKJHREREpKV4RY+IiDQmLxKIT7mPjOxc2JgZor2LJXR1JNWdFhGVgIUeERFpJDI5DXMPXEBaZq6izV5miNB+ngjysq/GzIioJOy6JSKiMkUmp2H8tkSlIg8A0jNzMX5bIiKT06opMyIqDa/o1QLG+nq4trhvdadBRHWUvEhg7oELEGpeEwAkAOYeuICennbsxiWqYXhFj4iIShWfcl/lSt6zBIC0zFzEp9yvuqSISCMs9IiIqFQZ2SUXeeWJI6Kqw0KPiIhKZWNmWKFxRFR1WOgREVGp2rtYwl5miJJG30nwdPZtexfLqkyLiDTAQo/oBT0ueAzvrd7w3uqNxwWPqzsdokqnqyNBaD9PAFAp9oqfh/bzrF0TMYrkQMrvwLk9T/+/SF7dGRFVCs66JSKiMgV52WPNW61V1tGzq43r6F3YD0ROB7Ju/a/N3AEIWgJ49q++vIgqAQs9IiLSSJCXPXp62tXuO2Nc2A98/zbw/GIxWWlP2wd/y2KPtAoLPSIi0piujgS+rlbVnUb5FMmfXskrbUXAyBmAe19AR7eKkyOqHByjR0REdcP1E8rdtSoEkHXzaRyRlmChR1UnPweYI3v6yM+p7myIqK55dLti44hqARZ6RERUN5jaVmwcUS3AQo+IiOqGRn5PZ9eWtiKgeYOncbUEl3uisrDQIyKiukFH9+kSKgBKXBEwaDEnYpBWYaFXC8iLBGKv3sNPSTcRe/Ue5EXqZowREVGZPPs/XULF/Ll1/8wduLQKaSUur1LDRSanqSxQal8bFyglIqopPPs/XULl+omnEy9MbZ921/JKHmkhFno1WGRyGsZvS1RZ8Sk9MxfjtyVizVutWewREZWHji7g0qW6syCqdOy6raHkRQJzD1wocVlPAJh74ELt6sZ99l6S107w3pJERESVrNYUeg8ePEBISAhkMhlkMhlCQkLw8OHDUrcRQmDOnDlwcHCAkZERunXrhvPnzyvFdOvWDRKJROkxdOjQSjwTzcSn3Ffqrn2eAJCWmYv4lPtVl9TLuLAfWNXuf893DAK+8nraTkRE5WIsNca5EedwbsQ5GEuNqzsdqoFqTaE3bNgwJCUlITIyEpGRkUhKSkJISEip24SFhWHZsmVYuXIlEhISYGdnh549eyI7O1spbsyYMUhLS1M81q1bV5mnopGM7JKLvPLEVavie0tmpyu3F99bspYVe/zBSkQ1hbxIjoT0BET8E4GE9ATI2VNCz6kVY/QuXryIyMhIxMXFoUOHDgCADRs2wNfXF5cuXYKbm5vKNkIIfPXVV5g1axZee+01AMDWrVtha2uLHTt2YOzYsYpYY2Nj2NnZVc3JaMjGzLBC46oN7y1JRFQpjlw/gsXxi3H78f/u5GFrbIsZ7WcgoFFANWZGNUmtuKIXGxsLmUymKPIAoGPHjpDJZDhxQv09CVNSUpCeno7AwEBFm4GBAfz9/VW22b59O6ytrdG8eXN8+OGHKlf8npeXl4esrCylR0Vr72IJe5lhact6wl5miPYulhV+7ArFe0sSEVW4I9ePYFr0NKUiDwAyHmdgWvQ0HLl+pJoyo5qmVhR66enpsLGxUWm3sbFBenq6mi2gaLe1Vb6Vja2trdI2w4cPx86dOxEdHY3Zs2dj7969iiuAJVm0aJFirKBMJoOTk9OLnlKZdHUkCO3nCaDEZT0R2s8TujollYI1BO8tSURUoeRFciyOXwyhpqekuG1J/BJ24xKAai705syZozIR4vnHqVOnAAASiWpBI4RQ2/6s519/fpsxY8YgICAAXl5eGDp0KPbs2YMjR44gMTGxxH3OnDkTmZmZiseNGzde5LQ1FuRljzVvtYadTLl71k5mWHuWVuG9JYmIKlRiRqLKlbxnCQikP05HYkbJv8eo7qjWMXoTJ04sc4ars7Mz/vzzT9y+rfqlvnPnjsoVu2LFY+7S09Nhb/+/gigjI6PEbQCgdevWkEqluHLlClq3bq02xsDAAAYGBqXmXVGCvOzR09MO8Sn3kZGdCxuzp921Nf5KXrHie0tmpUH9OD3J09dr0b0l5UVyJGYk4s7jO6hvXB+tbVpDl+MLiaiK3Hl8p0LjagIhl+PxqdMovHMHevXrw7htG0h0+XO1IlRroWdtbQ1ra+sy43x9fZGZmYn4+Hi0b98eAHDy5ElkZmbCz099geDi4gI7OztERUWhVatWAID8/HzExMRgyZIlarcBgPPnz6OgoECpOKxuujoS+LpaVXca5VN8b8nv31bzYu27tyQHPxNRdatvXL9C46pb1uHDuL1wEQqfGValZ2cH249nwvyZcfZUPrVijJ6HhweCgoIwZswYxMXFIS4uDmPGjEFwcLDSjFt3d3eEh4cDeNplO2XKFCxcuBDh4eFITk7GyJEjYWxsjGHDhgEArl69innz5uHUqVO4du0aIiIi8MYbb6BVq1bo1KlTtZyrViq+t6TZczOba9m9JbVt8LOQy5FzMh6ZPx9Ezsl4CDnH8xDVBq1tWsPW2BaSEqbrSSCBnbEdWtuo75WqSbIOH8bNyVOUijwAKLx9GzcnT0HW4cPVlJn2qBXLqwBPZ8ZOmjRJMYu2f//+WLlypVLMpUuXkJmZqXj+0Ucf4cmTJ3jvvffw4MEDdOjQAYcPH4aZmRkAQF9fH0ePHsXy5cvx6NEjODk5oW/fvggNDYUuLxlXLM/+QONuwOL/TlwZtgdo8kqtuZJX1uBnCSRYEr8E3Z2614puXP4FTVR76eroYkb7GZgWPQ0SSJR+LhUXf9PbT6/xP4uEXI7bCxcBQs2wHiEAiQS3Fy6CWY8e7MZ9CRIh1L3D9CKysrIgk8mQmZkJc3Pz6k6HKkFCegL+88t/yozb1GsT2tm1KzOuOhX/Ba3yw/W/k5QaLP+KxR5RLaBuKImdsR2mt59eK4aS5JyMR+qIEWXGNdy6FSYd2ldBRrXHi9QdteaKHlF10pbBz/wLmkh7BDQKQHen7rV2cljhHc1+XmoaR+qx0CPSgLYMfn586rTKWBglQqAwPR2PT53mX9BEtYCujm6N70UoiV59zX5eahpH6tWKyRhE1U1bBj/zL2giqimM27aBnp2dYtiICokEenZ2MG7bpmoT0zIs9Ig0UDz4GYBKsVebBj/zL2giqikkurqw/Xjmf588V+z997ntxzM5jOQlsdAj0lBAowAs67YMNsbKt+OzNbbFsm7LasXgZ/4FTUQ1iXlgIBos/wp6z93IQM/WlhPDKghn3VYAzrqtW2r7nTEUs24B5UkZnHVLRNWEd8Z4MS9Sd7DQqwAs9Ki24Tp6RES1F5dXIaJSmQcGwqxHD/4FTUSk5VjoEdVREl1drVlCpahIIO3KQ+Rk5cHE3AD2TetBR6eEcYhERHUICz0iqtWunsnA77uvIOdhnqLNpJ4BugxpCtdWNqVsSUSk/TjrlohqratnMhC5LlmpyAOAnId5iFyXjKtnMqopMyKimoGFHhHVSkVFAr/vvlJqzPHvr6CoiPPNiKjuYqFHRLVS2pWHKlfynvfoQR7SrjysmoSIiGogjtGrQnK5HAUFBdWdBlGNJZVKoavhzN+crNKLvBeNIyLSRiz0qoAQAunp6Xj48GF1p0JU49WrVw92dnaQlHT3jv8yMTfQaH+axhERVYSCPDnWT44BALy73B9Sg+pdtoqFXhUoLvJsbGxgbGxc5i8worpICIHHjx8jI+PpBAp7e/tS4+2b1oNJPYNSu29NLZ4utUJEVFex0KtkcrlcUeRZWVlVdzpENZqRkREAICMjAzY2NqV24+roSNBlSFNErksuMabz4KZcT4+I6jROxqhkxWPyjI2NqzkTotqh+N+KJuNZXVvZIGisF0zqKXfPmloYIGisF9fRI6I6j1f0qgi7a4k086L/Vlxb2cClRX3eGYOISA0WekRU6+noSNDAzaK60yAiUlq789aVB3DytKrWPzzZdUtarVu3bpgyZUp1p0FERHXA1TMZ2DEnTvH855V/4tuPT1TrXXpY6BGVYc6cOWjZsmV1p0FERDVY8S0ZH2fmK7VX9y0ZWehRnZWfn192EBERURlq8i0ZWeiR1sjJycHbb78NU1NT2Nvb44svvlB63dnZGfPnz8fIkSMhk8kwZswYAMD06dPRrFkzGBsbo3Hjxpg9e7ZixueWLVswd+5cnD17FhKJBBKJBFu2bKnqUyMiohqsJt+SsVyTMSwsLNTOjJNIJDA0NESTJk0wcuRIjBo16qUTJNLU//3f/+HYsWMIDw+HnZ0dPv74Y5w+fVqp23Xp0qWYPXs2PvnkE0WbmZkZtmzZAgcHB5w7dw5jxoyBmZkZPvroIwwZMgTJycmIjIzEkSNHAAAymayqT42IiGqwmnxLxnIVep9++ikWLFiA3r17o3379hBCICEhAZGRkZgwYQJSUlIwfvx4FBYWKq6aEFWmR48eYePGjfj222/Rs2dPAMDWrVvh6OioFPfKK6/gww8/VGp7tuhzdnbGBx98gN27d+Ojjz6CkZERTE1NoaenBzs7u8o/ESIiqnVq8i0Zy1XoHT9+HPPnz8e4ceOU2tetW4fDhw9j79698PHxwddff81Cj6rE1atXkZ+fD19fX0WbpaUl3NzclOLatm2rsu2ePXvw1Vdf4e+//8ajR49QWFgIc3PzSs+ZiIi0Q02+JWO5xuj98ssvCAgIUGnv0aMHfvnlFwBAnz598M8//7xcdkQaEkKzAa4mJiZKz+Pi4jB06FD07t0bP//8M86cOYNZs2ZxokYtU5Anx6pxv2LVuF9RkCev7nSIqI4pviVjaarrlozlKvQsLS1x4MABlfYDBw7A0tISwNOB8WZmZi+XHZGGmjRpAqlUiri4/61f9ODBA1y+fLnU7f744w80atQIs2bNQtu2bdG0aVNcv35dKUZfXx9yOYsHIiIqWfEtGY1l+krt1X1LxnJ13c6ePRvjx4/HsWPH0L59e0gkEsTHxyMiIgJr164FAERFRcHf379CkyUqiampKd555x383//9H6ysrGBra4tZs2ZBR6f0v2WaNGmC1NRU7Nq1C+3atcPBgwcRHh6uFOPs7IyUlBQkJSXB0dERZmZmMDCo+nEWRERUs7m2soGjuyW+mfobACB4ok/tvDPGmDFjEBMTAxMTE+zbtw979uyBsbExYmJi8M477wCAYkA7UVVZunQpunbtiv79+yMgIACdO3dGmzZtSt1mwIABmDp1KiZOnIiWLVvixIkTmD17tlLM66+/jqCgIHTv3h3169fHzp07K/M0iIioFnu2qHNoalHt992WCE0HN1GJsrKyIJPJkJmZqTKIPzc3FykpKXBxcYGhoWE1ZUhUe5Tn30xBnhzrJ8cAAN5d7g+pgW5lpkhEVKKq+HlUWt3xvHJ13QKAXC5HeHg4Ll68CIlEAg8PDwwYMAB6euXeJRERERFVoHJVZcnJyRgwYADS09MVy1dcvnwZ9evXx/79++Ht7V2hSRIRERHRiyvXGL3Ro0ejefPm+Pfff5GYmIjExETcuHEDPj4+ePfddys6RyIiIiIqh3IVemfPnsWiRYtgYWGhaLOwsMCCBQuQlJRUUbkpefDgAUJCQiCTySCTyRASEoKHDx+Wus2+ffvQq1cvWFtbQyKRqM0tLy8P77//PqytrWFiYoL+/fvj33//rZRzICIiIqpK5Sr03NzccPv2bZX2jIwMNGnS5KWTUmfYsGFISkpCZGQkIiMjkZSUhJCQkFK3ycnJQadOnbB48eISY6ZMmYLw8HDs2rULx48fx6NHjxAcHMx104iIiKjWK9cYvYULF2LSpEmYM2cOOnbsCODpHQbmzZuHJUuWICsrSxFbEbeSunjxIiIjIxEXF4cOHToAADZs2ABfX19cunRJ5TZXxYoLwWvXrql9PTMzExs3bsR3332nuNPHtm3b4OTkhCNHjqBXr14vnTsRERFRdSlXoRccHAwAGDx4MCSSp+vDFK/S0q9fP8VziURSIVfGYmNjIZPJFEUeAHTs2BEymQwnTpwosdAry+nTp1FQUIDAwEBFm4ODA7y8vHDixAkWekRERFSrlavQO3bsWEXnUar09HTY2KjeOsTGxgbp6ekvtV99fX2lsYYAYGtrW+p+8/LykJf3vxsXP3sFk4iIiKimKNcYPX9/f3To0AFGRkbIzs5GZmam0sPf31/xKM2cOXMgkUhKfZw6dQoAFFcOn1V81bCilbXfRYsWKSaFyGQyODk5VXgORERERC+rXFf0IiMj8fbbb+Pu3bsqr71Id+3EiRMxdOjQUmOcnZ3x559/qp38cefOHdja2mqWtBp2dnbIz8/HgwcPlK7qZWRkwM/Pr8TtZs6ciWnTpimeZ2VlaV2xt2jRIuzbtw9//fUXjIyM4OfnhyVLlpTaTT5y5Ehs3bpVpd3T0xPnz59XPH/48CFmzZqFffv24cGDB3BxccEXX3yBPn36VMq5EBER1VXlKvQmTpyIN954A59++ulLFVrW1tawtrYuM87X1xeZmZmIj49H+/btAQAnT55EZmZmqQVZWdq0aQOpVIqoqCgMHjwYAJCWlobk5GSEhYWVuJ2BgUGV39ReXiQQn3IfGdm5sDEzRHsXS+hW4v3zYmJiMGHCBLRr1w6FhYWYNWsWAgMDceHCBZiYmKjdZvny5UoznAsLC9GiRQu88cYbirb8/Hz07NkTNjY22LNnDxwdHXHjxg2YmZlV2rkQERHVVeUq9DIyMjBt2rSXKvJehIeHB4KCgjBmzBisW7cOAPDuu+8iODhY6QqTu7s7Fi1ahFdffRUAcP/+faSmpuLWrVsAgEuXLgF4eiXPzs4OMpkM77zzDj744ANYWVnB0tISH374Iby9vRWzcGuCyOQ0zD1wAWmZuYo2e5khQvt5IsjLvnKOGRmp9Hzz5s2wsbHB6dOn0bVrV7XbFHdlF/vxxx/x4MEDjBo1StG2adMm3L9/HydOnIBUKgUANGrUqBLOgIiIiMo1Rm/QoEGIjo6u4FRKt337dnh7eyMwMBCBgYHw8fHBd999pxRz6dIlZGZmKp7v378frVq1Qt++fQEAQ4cORatWrbB27VpFzJdffomBAwdi8ODB6NSpE4yNjXHgwAHo6taMm6JHJqdh/LZEpSIPANIzczF+WyIik9OqJI/i99XS0lLjbTZu3IiAgAClQm7//v3w9fXFhAkTYGtrCy8vLyxcuJDrFlYDIZfjSXIyniQnQ/D9JyLSSuW6ordy5Uq88cYb+P333+Ht7a24MlNs0qRJFZLcsywtLbFt27ZSY4qXeCk2cuRIjBw5stRtDA0NsWLFCqxYseJlU9SIEAJPCjT7pSovEgjdfx5CzWsCgATAnP0X0KmJtUbduEZS3XJNXhFCYNq0aejcuTO8vLw02iYtLQ2HDh3Cjh07lNr/+ecf/Prrrxg+fDgiIiJw5coVTJgwAYWFhfj0009fODciIiIqWbkKvR07duCXX36BkZERoqOjlYoHiURSKYWetnhSIIfnp79UyL4EgPSsXHjPOaxR/IV5vWCs/+If+cSJE/Hnn3/i+PHjGm+zZcsW1KtXDwMHDlRqLyoqgo2NDdavXw9dXV20adMGt27dwtKlS1noERERVbByFXqffPIJ5s2bhxkzZkBHp1y9v1RLvP/++9i/fz9+++03ODo6arSNEAKbNm1CSEgI9PX1lV6zt7eHVCpV6hr38PBAeno68vPzVeKJiIio/MpV6OXn52PIkCEs8srBSKqLC/M0u+NGfMp9jNycUGbcllHt0N6l7LFzRlLNxx0KIfD+++8jPDwc0dHRcHFx0XjbmJgY/P3333jnnXdUXuvUqRN27NiBoqIixffn8uXLsLe3Z5FHRERUwcpVqY0YMQK7d++u6FzqBIlEAmN9PY0eXZrWh73MECWNqpPg6ezbLk3ra7S/FxmfN2HCBGzbtg07duyAmZkZ0tPTkZ6ejidPnihiZs6cibfffltl240bN6JDhw5qx/ONHz8e9+7dw+TJk3H58mUcPHgQCxcuxIQJEzTOjYiIiDRTrit6crkcYWFh+OWXX+Dj46MyGWPZsmUVklxdp6sjQWg/T4zflggJoDQpo7hkC+3nWSnr6a1ZswYA0K1bN6X2zZs3Kya4pKWlITU1Ven1zMxM7N27F8uXL1e7XycnJxw+fBhTp06Fj48PGjRogMmTJ2P69OkVfg5ERER1XbkKvXPnzqFVq1YAgOTkZKXXKuOWZHVZkJc91rzVWmUdPbtKXkfv+RnM6mzZskWlTSaT4fHjx6Vu5+vri7i4uPKmRkRERBoqV6F37Nixis6DShHkZY+ennZVemcMIiIiqv3KVehR1dPVkcDX1aq60yAiIqJahNNmiYiIiLQUCz0iIiIiLcVCj4iIiEhLsdAjIiIi0lKcjEFUR0l0dWGkZlFrIiLSHryiR0RERKSlWOgRERERaSkWekRERERaioUeqbVmzRr4+PjA3Nwc5ubm8PX1xaFDh0rdJiYmBm3atIGhoSEaN26MtWvXqsR89dVXcHNzg5GREZycnDB16lTk5uaq2RsRERG9LE7GqC2K5MD1E8Cj24CpLdDID9DRrbTDOTo6YvHixWjSpAkAYOvWrRgwYADOnDmD5s2bq8SnpKSgT58+GDNmDLZt24Y//vgD7733HurXr4/XX38dALB9+3bMmDEDmzZtgp+fHy5fvoyRI0cCAL788stKOxciIqK6ioVebXBhPxA5Hci69b82cwcgaAng2b9SDtmvXz+l5wsWLMCaNWsQFxenttBbu3YtGjZsiK+++goA4OHhgVOnTuHzzz9XFHqxsbHo1KkThg0bBgBwdnbGm2++ifj4+Eo5ByIiorqOXbc13YX9wPdvKxd5AJCV9rT9wv5KT0Eul2PXrl3IycmBr6+v2pjY2FgEBgYqtfXq1QunTp1CQUEBAKBz5844ffq0orD7559/EBERgb59+1buCRAREdVRvKJX1YQACh5rFlskBw59BECo2xEAydMrfY27adaNKzUGJBKNUz137hx8fX2Rm5sLU1NThIeHw9PTU21seno6bG1tldpsbW1RWFiIu3fvwt7eHkOHDsWdO3fQuXNnCCFQWFiI8ePHY8aMGRrnRERERJpjoVfVCh4DCx0qaGfi6ZW+xU6ahX98C9A30Xjvbm5uSEpKwsOHD7F3716MGDECMTExJRZ7kueKSCGEUnt0dDQWLFiA1atXo0OHDvj7778xefJk2NvbY/bs2RrnRURERJphoUcl0tfXV0zGaNu2LRISErB8+XKsW7dOJdbOzg7p6elKbRkZGdDT04OVlRUAYPbs2QgJCcHo0aMBAN7e3sjJycG7776LWbNmQUeHIwmIiIgqEgu9qiY1fnplTRPXTwDbB5UdN3zP01m4mhz7JQghkJeXp/Y1X19fHDhwQKnt8OHDaNu2LaRSKQDg8ePHKsWcrq4uhBCKq39ERERUcVjoVTWJRPPuU9dXns6uzUqD+nF6kqevu75S4UutfPzxx+jduzecnJyQnZ2NXbt2ITo6GpGRkQCAmTNn4ubNm/j2228BAOPGjcPKlSsxbdo0jBkzBrGxsdi4cSN27typ2Ge/fv2wbNkytGrVStF1O3v2bPTv3x+6upW3VAwREVFdxUKvJtPRfbqEyvdvA5BAudj773i4oMWVsp7e7du3ERISgrS0NMhkMvj4+CAyMhI9e/YEAKSlpSE1NVUR7+LigoiICEydOhWrVq2Cg4MDvv76a8XSKgDwySefQCKR4JNPPsHNmzdRv3599OvXDwsWLKjw/ImIiIiFXs3n2R8Y/G0J6+gtrrR19DZu3Fjq61u2bFFp8/f3R2JiYonb6OnpITQ0FKGhoS+bHhEREWmAhV5t4NkfcO9bpXfGICIiotqPhV5toaMLuHSp7iyIiIioFuF6FkRERERaioUeERERkZZioUdERESkpVjoEREREWkpFnpEREREWoqzbomo1pMa6GLC2leqOw0iohqHV/SIiIiItFStKfQePHiAkJAQyGQyyGQyhISE4OHDh6Vus2/fPvTq1QvW1taQSCRISkpSienWrRskEonSY+jQoZVzEkRERERVqNYUesOGDUNSUhIiIyMRGRmJpKQkhISElLpNTk4OOnXqhMWLF5caN2bMGKSlpSke69atq8jUK4S8SI6E9ARE/BOBhPQEyIvklXq8NWvWwMfHB+bm5jA3N4evry8OHTpU6jZ5eXmYNWsWGjVqBAMDA7i6umLTpk2K1wsKCjBv3jy4urrC0NAQLVq0QGRkpNI+5syZo1J429nZKcUIITBnzhw4ODjAyMgI3bp1w/nz51Vyef/992FtbQ0TExP0798f//77r1KMJn88pKamol+/fjAxMYG1tTUmTZqE/Px8pZhz587B398fRkZGaNCgAebNmwchhFJMTEwM2rRpA0NDQzRu3Bhr165Vef/27t0LT09PGBgYwNPTE+Hh4Soxq1evhouLCwwNDdGmTRv8/vvvNfa9uXTpErp37w5bW1vFeX/yyScoKChQOS8iIqokoha4cOGCACDi4uIUbbGxsQKA+Ouvv8rcPiUlRQAQZ86cUXnN399fTJ48+aXyy8zMFABEZmamymtPnjwRFy5cEE+ePCn3/qOuRYke3/cQXlu8FI8e3/cQUdeiXibtUu3fv18cPHhQXLp0SVy6dEl8/PHHQiqViuTk5BK36d+/v+jQoYOIiooSKSkp4uTJk+KPP/5QvP7RRx8JBwcHcfDgQXH16lWxevVqYWhoKBITExUxoaGhonnz5iItLU3xyMjIUDrO4sWLhZmZmdi7d684d+6cGDJkiLC3txdZWVmKmHHjxokGDRqIqKgokZiYKLp37y5atGghCgsLFTFBQUHCy8tLnDhxQpw4cUJ4eXmJ4OBgxeuFhYXCy8tLdO/eXSQmJoqoqCjh4OAgJk6cqIjJzMwUtra2YujQoeLcuXNi7969wszMTHz++eeKmH/++UcYGxuLyZMniwsXLogNGzYIqVQq9uzZo4g5ceKE0NXVFQsXLhQXL14UCxcuFHp6ekrf+V27dgmpVCo2bNggLly4ICZPnixMTEzE9evXa+R7c/XqVbFp0yaRlJQkrl27Jn766SdhY2MjZs6cqfb7U6wi/s0QEVWX/NxCsXLsUbFy7FGRn1tY9gblUFrd8bxaUeht3LhRyGQylXaZTCY2bdpU5vZlFXrW1tbCyspKeHp6ig8++EDpl6ImKrPQi7oWJby3eCsVeV5bvIT3Fm/hvcW7Uou951lYWIhvvvlG7WuHDh0SMplM3Lt3r8Tt7e3txcqVK5XaBgwYIIYPH654HhoaKlq0aFHiPoqKioSdnZ1YvHixoi03N1fIZDKxdu1aIYQQDx8+FFKpVOzatUsRc/PmTaGjoyMiIyOFEJr98RARESF0dHTEzZs3FTE7d+4UBgYGis969erVQiaTidzcXEXMokWLhIODgygqKhJCPC1w3d3dlc5j7NixomPHjorngwcPFkFBQUoxvXr1EkOHDlU8b9++vRg3bpxSjLu7u5gxY0aNfG/UmTp1qujcuXOJrwvBQo+IareaVujViq7b9PR02NjYqLTb2NggPT39pfY9fPhw7Ny5E9HR0Zg9ezb27t2L1157rdRt8vLykJWVpfTQlBACjwsea/TIzsvGovhFEBCq+/nv/xbHL0Z2XrZG+xNCdT+akMvl2LVrF3JycuDr66s2Zv/+/Wjbti3CwsLQoEEDNGvWDB9++CGePHmiiMnLy4OhoaHSdkZGRjh+/LhS25UrV+Dg4AAXFxcMHToU//zzj+K1lJQUpKenIzAwUNFmYGAAf39/nDhxAgBw+vRpFBQUKMU4ODjAy8tLERMbGwuZTIYOHTooYjp27AiZTKYU4+XlBQcHB0VMr169kJeXh9OnTyti/P39YWBgoBRz69YtXLt2TRHzbC7FMadOnVJ0Y5YUU5xLfn4+Tp8+rRITGBioiKlp783z/v77b0RGRsLf31/t60RE2qB4FYAJa1+B1EC3utOp3uVV5syZg7lz55Yak5CQAACQSCQqrwkh1La/iDFjxij+28vLC02bNkXbtm2RmJiI1q1bq91m0aJFZeZdkieFT9BhR4eyAzV0+/Ft+O3y0yj25LCTMJYaa7zvc+fOwdfXF7m5uTA1NUV4eDg8PT3Vxv7zzz84fvw4DA0NER4ejrt37+K9997D/fv3FeP0evXqhWXLlqFr165wdXXF0aNH8dNPP0Eu/994ww4dOuDbb79Fs2bNcPv2bcyfPx9+fn44f/48rKysFIW9ra2t0vFtbW1x/fp1AE//MNDX14eFhYVKTPH2mvzxkJ6ernIcCwsL6OvrK8U4OzurHKf4NRcXF7X7sbW1RWFhIe7evQt7e/sSY4qPc/fuXcjl8lJjatp7U8zPzw+JiYnIy8vDu+++i3nz5qnsm4iIKke1XtGbOHEiLl68WOrDy8sLdnZ2uH37tsr2d+7cUfll87Jat24NqVSKK1eulBgzc+ZMZGZmKh43btyo0BxqCjc3NyQlJSEuLg7jx4/HiBEjcOHCBbWxRUVFkEgk2L59O9q3b48+ffpg2bJl2LJli+Kq3vLly9G0aVO4u7tDX18fEydOxKhRo6Cr+7+/eHr37o3XX38d3t7eCAgIwMGDBwEAW7duVTre8wW+JkX/8zGa/PFQnpjiK6cVEfN8W0XFPK+y3hsA2L17NxITE7Fjxw4cPHgQn3/+eam5EBFRxanWK3rW1tawtrYuM87X1xeZmZmIj49H+/btAQAnT55EZmYm/Pw0u5qlqfPnz6OgoAD29vYlxhgYGCh11b0IIz0jnBx2UqPY07dP472j75UZt7rHarSxbaPRsV+Evr4+mjRpAgBo27YtEhISsHz5crWzku3t7dGgQQPIZDJFm4eHB4QQ+Pfff9G0aVPUr18fP/74I3Jzc3Hv3j04ODhgxowZcHFxKTEHExMTeHt7Kwrv4hm46enpSp9RRkaGoui3s7NDfn4+Hjx4oHTlKiMjQ/F90eSPBzs7O5w8qfxZPXjwAAUFBUoxz1/BysjIAIAyY/T09GBlZVVqTPE+rK2toaurW2pMTXtvijk5OQEAPD09IZfL8e677+KDDz5QKvCJiKhy1Ioxeh4eHggKCsKYMWMQFxeHuLg4jBkzBsHBwXBzc1PEubu7Ky1Jcf/+fSQlJSmuQl26dAlJSUmKX5ZXr17FvHnzcOrUKVy7dg0RERF444030KpVK3Tq1KlSzkUikcBYaqzRw8/BD7bGtpBA/dUYCSSwM7aDn4OfRvt72W5uIQTy8vLUvtapUyfcunULjx49UrRdvnwZOjo6cHR0VIo1NDREgwYNUFhYiL1792LAgAElHjMvLw8XL15UFC4uLi6ws7NDVFSUIiY/Px8xMTGKQqVNmzaQSqVKMWlpaUhOTlbEPPvHQ7Hn/3jw9fVFcnIy0tLSFDGHDx+GgYEB2rRpo4j57bfflJYVOXz4MBwcHBRdur6+vkq5FMe0bdsWUqm01JjiXPT19dGmTRuVmKioKEVMTXtv1BFCoKCgoNzjRYmI6AVV9EyQynLv3j0xfPhwYWZmJszMzMTw4cPFgwcPlGIAiM2bNyueb968WQBQeYSGhgohhEhNTRVdu3YVlpaWQl9fX7i6uopJkyaVOnNUnaqYdfv8zNvKnnU7c+ZM8dtvv4mUlBTx559/io8//ljo6OiIw4cPCyGEmDFjhggJCVHEZ2dnC0dHRzFo0CBx/vx5ERMTI5o2bSpGjx6tiImLixN79+4VV69eFb/99pt45ZVXhIuLi9Ln+MEHH4jo6Gjxzz//iLi4OBEcHCzMzMzEtWvXFDGLFy8WMplM7Nu3T5w7d068+eabapcQcXR0FEeOHBGJiYnilVdeUbuEiI+Pj4iNjRWxsbHC29tb7RIiPXr0EImJieLIkSPC0dFRaQmRhw8fCltbW/Hmm2+Kc+fOiX379glzc3O1y6tMnTpVXLhwQWzcuFFleZU//vhD6OrqisWLF4uLFy+KxYsXl7i8ysaNG8WFCxfElClThImJSY19b7Zt2yZ2794tLly4IK5evSq+//570aBBA6VZ1upw1i0RUem0bnmVmq461tEL+D6gUpdW+c9//iMaNWok9PX1Rf369UWPHj0URZ4QQowYMUL4+/srbXPx4kUREBAgjIyMhKOjo5g2bZp4/Pix4vXo6Gjh4eEhDAwMhJWVlQgJCVFankMIoVj3TSqVCgcHB/Haa6+J8+fPK8UUFRWJ0NBQYWdnJwwMDETXrl3FuXPnlGKePHkiJk6cKCwtLYWRkZEIDg4WqampSjGa/PFw/fp10bdvX2FkZCQsLS3FxIkTlZZSEUKIP//8U3Tp0kUYGBgIOzs7MWfOHMXSKs+ee6tWrYS+vr5wdnYWa9asUXnPf/jhB+Hm5iakUqlwd3cXe/fuVYlZtWqV4nNp3bq1iImJqbHvza5du0Tr1q2FqampMDExEZ6enmLhwoVl/ltgoUdEVLoXKfQkQrAP5WVlZWVBJpMhMzMT5ubmSq/l5uYiJSVFcTeD8pIXyZGYkYg7j++gvnF9tLZpDV0djnEi7VNR/2aIiLRVaXXH86p1MgZpTldHF+3s2lV3GkRERFSL1IrJGERERET04ljoEREREWkpFnpEREREWoqFHhEREZGWYqFHREREpKVY6BERERFpKRZ6RERERFqKhR4RERGRlmKhR0RERKSlWOjVEkIuR87JeGT+fBA5J+Mh5PJKPd6aNWvg4+MDc3NzmJubw9fXF4cOHSp1m1WrVsHDwwNGRkZwc3PDt99+q/T6hg0b0KVLF1hYWMDCwgIBAQGIj4+vzNMgIiKq03gLtFog6/Bh3F64CIXp6Yo2PTs72H48E+aBgZVyTEdHRyxevBhNmjQBAGzduhUDBgzAmTNn0Lx5c5X4NWvWYObMmdiwYQPatWuH+Ph4jBkzBhYWFujXrx8AIDo6Gm+++Sb8/PxgaGiIsLAwBAYG4vz582jQoEGlnAcREVFdJhFCiOpOorYr7ebCL3uD9qzDh3Fz8hTg+Y9JIgEANFj+VaUVe8+ztLTE0qVL8c4776i85ufnh06dOmHp0qWKtilTpuDUqVM4fvy42v3J5XJYWFhg5cqVePvttystb6pdXvbfDBGRtiut7nger+hVMSEExJMnmsXK5bg9f4Fqkfd0R4AEuL1gIUx8fSHR1S1zfxIjI0j+WyC+CLlcjh9++AE5OTnw9fVVG5OXl6fyS9nIyAjx8fEoKCiAVCpV2ebx48coKCiApaXlC+dEREREZWOhV8XEkye41LpNBe0MKLx9G5fbtdco3C3xNCTGxhrv/ty5c/D19UVubi5MTU0RHh4OT09PtbG9evXCN998g4EDB6J169Y4ffo0Nm3ahIKCAty9exf29vYq28yYMQMNGjRAQECAxjkRERGR5ljoUYnc3NyQlJSEhw8fYu/evRgxYgRiYmLUFnuzZ89Geno6OnbsCCEEbG1tMXLkSISFhUFXzdXGsLAw7Ny5E9HR0eyeIyIiqiQco1cBXmSM3ot03T4+dQo33h1bZpzT+nUwbtu2zLjydt0WCwgIgKurK9atW1diTEFBAW7fvg17e3usX78e06dPx8OHD6Gj878J3p9//jnmz5+PI0eOoK0GeVPdwjF6RESl4xi9GkwikWjcfWrSqRP07OxQePu2+nF6Egn0bG1h0qmTRmP0XpYQAnl5eaXGSKVSODo6AgB27dqF4OBgpSJv6dKlmD9/Pn755RcWeURERJWMhV4NJtHVhe3HM5/OupVIlIu9/16Zs/14ZqUUeR9//DF69+4NJycnZGdnY9euXYiOjkZkZCQAYObMmbh586ZirbzLly8jPj4eHTp0wIMHD7Bs2TIkJydj69atin2GhYVh9uzZ2LFjB5ydnZH+3+ViTE1NYWpqWuHnQEREVNdxweQazjwwEA2WfwU9W1uldj1b20pdWuX27dsICQmBm5sbevTogZMnTyIyMhI9e/YEAKSlpSE1NVURL5fL8cUXX6BFixbo2bMncnNzceLECTg7OytiVq9ejfz8fAwaNAj29vaKx+eff14p50BERFTXcYxeBajMdfSKCbkcj0+dRuGdO9CrXx/GbdtUSXctUVXjGD0iotJxjJ4WkujqwqSDZsuoEBEREQHsuiUiIiLSWiz0iIiIiLQUCz0iIiIiLcVCj4iIiEhLsdAjIiIi0lIs9IiIiIi0FAs9IiIiIi3FQo+IiIhIS7HQIyIiItJSLPRqiaIigZuXHuByQjpuXnqAoqKqvXPdokWLIJFIMGXKlFLjVq1aBQ8PDxgZGcHNzQ3ffvut0uv79u1D27ZtUa9ePZiYmKBly5b47rvvKjFzIiKiuou3QKsFrp7JwO+7ryDnYZ6izaSeAboMaQrXVjaVfvyEhASsX78ePj4+pcatWbMGM2fOxIYNG9CuXTvEx8djzJgxsLCwQL9+/QAAlpaWmDVrFtzd3aGvr4+ff/4Zo0aNgo2NDXr16lXp50JERFSX8IpeDXf1TAYi1yUrFXkAkPMwD5HrknH1TEalHv/Ro0cYPnw4NmzYAAsLi1Jjv/vuO4wdOxZDhgxB48aNMXToULzzzjtYsmSJIqZbt2549dVX4eHhAVdXV0yePBk+Pj44fvx4pZ4HERFRXcRCr4oJIVCQJ9fokfekEL/vvlzq/n7ffQV5Two12p8QL97dO2HCBPTt2xcBAQFlxubl5cHQ0FCpzcjICPHx8SgoKFCJF0Lg6NGjuHTpErp27frCuREREVHpak3X7YMHDzBp0iTs378fANC/f3+sWLEC9erVUxtfUFCATz75BBEREfjnn38gk8kQEBCAxYsXw8HBQRGXl5eHDz/8EDt37sSTJ0/Qo0cPrF69Go6OjpVyHoX5RVg/OabC9pfzMA/fTP1No9h3l/tDaqCr8b537dqFxMREJCQkaBTfq1cvfPPNNxg4cCBat26N06dPY9OmTSgoKMDdu3dhb28PAMjMzESDBg2Ql5cHXV1drF69Gj179tQ4LyIiItJMrbmiN2zYMCQlJSEyMhKRkZFISkpCSEhIifGPHz9GYmIiZs+ejcTEROzbtw+XL19G//79leKmTJmC8PBw7Nq1C8ePH8ejR48QHBwMuVxe2adUo924cQOTJ0/Gtm3bVK7SlWT27Nno3bs3OnbsCKlUigEDBmDkyJEAAF3d/xWYZmZmSEpKQkJCAhYsWIBp06YhOjq6Es6CiIiobpOI8vTnVbGLFy/C09MTcXFx6NChAwAgLi4Ovr6++Ouvv+Dm5qbRfhISEtC+fXtcv34dDRs2RGZmJurXr4/vvvsOQ4YMAQDcunULTk5OiIiI0HhyQFZWFmQyGTIzM2Fubq70Wm5uLlJSUuDi4gJDQ0MIIVCYX6TRfm9deYifV54tMy54Ygs4NK1XZpyevg4kEolGx/7xxx/x6quvKhVocrkcEokEOjo6iqtx6hQUFOD27duwt7fH+vXrMX36dDx8+BA6Our/rhg9ejRu3LiBX375RaPcSLs9/2+GiIiUlVZ3PK9WdN3GxsZCJpMpijwA6NixI2QyGU6cOKFxoZeZmQmJRKLo7j19+jQKCgoQGBioiHFwcICXlxdOnDhRYqGXl5eHvLz/TY7IysrS+FwkEonG3adOnpYwqWegMhHjWaYWBnDytISOjmYFnKZ69OiBc+fOKbWNGjUK7u7umD59eolFHgBIpVJF1/euXbsQHBxcYpEHPB2r9+z7SUREdUvR48e41LoNAMAt8TR0jI2rOSPtUSsKvfT0dNjYqC4jYmNjg/T0dI32kZubixkzZmDYsGGK6jc9PR36+voqs0ltbW1L3e+iRYswd+7cFziD8tHRkaDLkKaIXJdcYkznwU0rvMgDnnavenl5KbWZmJjAyspK0T5z5kzcvHlTsVbe5cuXER8fjw4dOuDBgwdYtmwZkpOTsXXrVsU+Fi1ahLZt28LV1RX5+fmIiIjAt99+izVr1lT4ORAREdV11TpGb86cOZBIJKU+Tp06BQBquxyFEBp1RRYUFGDo0KEoKirC6tWry4wva78zZ85EZmam4nHjxo0y91lerq1sEDTWCyb1DJTaTS0MEDTWq0rW0StJWloaUlNTFc/lcjm++OILtGjRAj179kRubi5OnDgBZ2dnRUxOTg7ee+89NG/eHH5+ftizZw+2bduG0aNHV8MZEBERabdqvaI3ceJEDB06tNQYZ2dn/Pnnn7h9+7bKa3fu3IGtrW2p2xcUFGDw4MFISUnBr7/+qtSXbWdnh/z8fDx48EDpql5GRgb8/PxK3KeBgQEMDAxKfL2iubaygUuL+ki78hA5WXkwMTeAfdN6lXIlrzTPT5jYsmWL0nMPDw+cOXOm1H3Mnz8f8+fPr+DMiIiISJ1qLfSsra1hbW1dZpyvry8yMzMRHx+P9u3bAwBOnjyJzMzMUguy4iLvypUrOHbsGKysrJReb9OmDaRSKaKiojB48GAAT69SJScnIyws7CXOrOLp6EjQwK30BYuJiIiInlUrllfx8PBAUFAQxowZg7i4OMTFxWHMmDEIDg5Wmojh7u6O8PBwAEBhYSEGDRqEU6dOYfv27ZDL5UhPT0d6ejry8/MBADKZDO+88w4++OADHD16FGfOnMFbb70Fb29vjRYIJiIiIqrJasVkDADYvn07Jk2apJgh279/f6xcuVIp5tKlS8jMzAQA/Pvvv4rFlVu2bKkUd+zYMXTr1g0A8OWXX0JPTw+DBw9WLJi8ZcuWUmeVEhEREdUGtabQs7S0xLZt20qNeXZJQGdnZ41u+WVoaIgVK1ZgxYoVL50jERERUU1Sawo9IiIi0k46xsbw+OtidaehlWrFGD1tUAtuQEJUI/DfChFRxWGhV8mkUimAp/feJaKyFf9bKf63Q0RE5ceu20qmq6uLevXqISMjAwBgbGys8f1mieoSIQQeP36MjIwM1KtXjxOiiIgqAAu9KmBnZwcAimKPiEpWr149xb8ZIiJ6OSz0qoBEIoG9vT1sbGxQUFBQ3ekQ1VhSqZRX8oiIKhALvSqkq6vLX2JERERUZTgZg4iIiEhLsdAjIiIi0lIs9IiIiIi0FMfoVYDiBV6zsrKqORMiIiLSdsX1hiYLzLPQqwDZ2dkAACcnp2rOhIiIiOqK7OxsyGSyUmMkgvcbemlFRUW4desWzMzMuBhyHZGVlQUnJyfcuHED5ubm1Z0OEdVh/HlU9wghkJ2dDQcHB+jolD4Kj1f0KoCOjg4cHR2rOw2qBubm5vzBSkQ1An8e1S1lXckrxskYRERERFqKhR4RERGRlmKhR1QOBgYGCA0NhYGBQXWnQkR1HH8eUWk4GYOIiIhIS/GKHhEREZGWYqFHREREpKVY6BERERFpKRZ6ROWwevVquLi4wNDQEG3atMHvv/9e3SkRUR3z22+/oV+/fnBwcIBEIsGPP/5Y3SlRDcRCj+gF7d69G1OmTMGsWbNw5swZdOnSBb1790Zqamp1p0ZEdUhOTg5atGiBlStXVncqVINx1i3RC+rQoQNat26NNWvWKNo8PDwwcOBALFq0qBozI6K6SiKRIDw8HAMHDqzuVKiG4RU9oheQn5+P06dPIzAwUKk9MDAQJ06cqKasiIiI1GOhR/QC7t69C7lcDltbW6V2W1tbpKenV1NWRERE6rHQIyoHiUSi9FwIodJGRERU3VjoEb0Aa2tr6Orqqly9y8jIULnKR0REVN1Y6BG9AH19fbRp0wZRUVFK7VFRUfDz86umrIiIiNTTq+4EiGqbadOmISQkBG3btoWvry/Wr1+P1NRUjBs3rrpTI6I65NGjR/j7778Vz1NSUpCUlARLS0s0bNiwGjOjmoTLqxCVw+rVqxEWFoa0tDR4eXnhyy+/RNeuXas7LSKqQ6Kjo9G9e3eV9hEjRmDLli1VnxDVSCz0iIiIiLQUx+gRERERaSkWekRERERaioUeERERkZZioUdERESkpVjoEREREWkpFnpEREREWoqFHhEREZGWYqFHREREpKVY6BERERFpKRZ6RERERFqKhR4RERGRlmKhR0RUSSIjI9G5c2fUq1cPVlZWCA4OxtWrVwEA165dg0Qiwffff48uXbrAyMgI7dq1w+XLl5GQkIC2bdvC1NQUQUFBuHPnTjWfCRHVVhIhhKjuJIiItNHevXshkUjg7e2NnJwcfPrpp7h27RqSkpKQmpoKFxcXuLu746uvvkLDhg3xn//8B/n5+TA3N8f8+fNhbGyMwYMHIyAgAGvWrKnu0yGiWoiFHhFRFblz5w5sbGxw7tw5mJqawsXFBd988w3eeecdAMCuXbvw5ptv4ujRo3jllVcAAIsXL8aWLVvw119/VWfqRFRLseuWiKiSXL16FcOGDUPjxo1hbm4OFxcXAEBqaqoixsfHR/Hftra2AABvb2+ltoyMjCrKmIi0jV51J0BEpK369esHJycnbNiwAQ4ODigqKoKXlxfy8/MVMVKpVPHfEolEbVtRUVHVJU1EWoWFHhFRJbh37x4uXryIdevWoUuXLgCA48ePV3NWRFTXsNAjIqoEFhYWsLKywvr162Fvb4/U1FTMmDGjutMiojqGY/SIiCqBjo4Odu3ahdOnT8PLywtTp07F0qVLqzstIqpjOOuWiIiISEvxih4RERGRlmKhR0RERKSlWOgRERERaSkWekRERERaioUeERERkZZioUdERESkpVjoEREREWkpFnpEREREWoqFHhEREZGWYqFHREREpKVY6BERERFpKRZ6RERERFrq/wF8QwAtjkiQagAAAABJRU5ErkJggg==",
"text/plain": [
"<Figure size 700x300 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# HDI does not include the mean contrast comparison...\n",
"fig, ax, comparisons_df, contrast_df, idata_new = plot_comparison(\n",
" mt_model,\n",
" mt_idata,\n",
" contrast_predictor=\"hp\",\n",
" conditional=[\"am\", \"drat\"]\n",
")\n",
"fig.set_size_inches(7, 3)\n",
"plt.title(\"Difference in mpg between am and drat\\nfor hp contrast [110, 175]\");"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>estimate</th>\n",
" <th>hdi_3%</th>\n",
" <th>hdi_97%</th>\n",
" <th>term</th>\n",
" <th>contrast</th>\n",
" <th>am</th>\n",
" <th>drat</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>-0.046206</td>\n",
" <td>-0.003756</td>\n",
" <td>-0.026338</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>2.760</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.049976</td>\n",
" <td>-0.037245</td>\n",
" <td>-0.056671</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.080</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>-0.057222</td>\n",
" <td>-0.047034</td>\n",
" <td>-0.031519</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.695</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.059873</td>\n",
" <td>-0.119706</td>\n",
" <td>-0.104732</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.920</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>-0.071772</td>\n",
" <td>-0.142725</td>\n",
" <td>-0.086004</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>4.930</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>-0.028507</td>\n",
" <td>-0.001350</td>\n",
" <td>0.017412</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>2.760</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>-0.036580</td>\n",
" <td>-0.014789</td>\n",
" <td>-0.020630</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.080</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>-0.052097</td>\n",
" <td>-0.038705</td>\n",
" <td>-0.045089</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.695</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>-0.057774</td>\n",
" <td>-0.207495</td>\n",
" <td>-0.198384</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.920</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>-0.083257</td>\n",
" <td>-0.148901</td>\n",
" <td>-0.072400</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>4.930</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" estimate hdi_3% hdi_97% term contrast am drat\n",
"0 -0.046206 -0.003756 -0.026338 hp [146, 147] 0 2.760\n",
"1 -0.049976 -0.037245 -0.056671 hp [146, 147] 0 3.080\n",
"2 -0.057222 -0.047034 -0.031519 hp [146, 147] 0 3.695\n",
"3 -0.059873 -0.119706 -0.104732 hp [146, 147] 0 3.920\n",
"4 -0.071772 -0.142725 -0.086004 hp [146, 147] 0 4.930\n",
"5 -0.028507 -0.001350 0.017412 hp [146, 147] 1 2.760\n",
"6 -0.036580 -0.014789 -0.020630 hp [146, 147] 1 3.080\n",
"7 -0.052097 -0.038705 -0.045089 hp [146, 147] 1 3.695\n",
"8 -0.057774 -0.207495 -0.198384 hp [146, 147] 1 3.920\n",
"9 -0.083257 -0.148901 -0.072400 hp [146, 147] 1 4.930"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"contrast_df"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using xarray\n",
"\n",
"For computing comparisons of posterior samples"
]
},
{
"cell_type": "code",
"execution_count": 260,
"metadata": {},
"outputs": [],
"source": [
"compare_xr = comparisons_df.iloc[:, :3].to_xarray().to_array().T"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"If we select only one chain and draw, then we can create an xarray Dataset with the posterior samples and data used for performing predictions."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"ename": "ValueError",
"evalue": "Could not convert tuple of form (dims, data[, attrs, encoding]): (['preds'], array([[[19.59478712, 19.54856056, 19.23140441, ..., 22.49626788,\n 22.51163257, 22.44590673],\n [18.27013994, 18.22625734, 18.58204345, ..., 23.3491659 ,\n 19.44539658, 19.35357745],\n [18.85953525, 18.81616193, 18.31847112, ..., 23.20774708,\n 22.53151746, 22.45055377],\n ...,\n [16.95947443, 16.91816056, 17.61175804, ..., 22.68251523,\n 22.94911059, 22.84984865],\n [16.75874405, 16.71122187, 17.75769402, ..., 23.31604767,\n 24.31327371, 24.22783439],\n [15.98276926, 15.93614123, 16.72108478, ..., 21.95560998,\n 17.56669211, 17.48011585]],\n\n [[15.13806506, 15.07475069, 16.57548916, ..., 22.93274638,\n 21.45973418, 21.40293702],\n [14.03187854, 13.96100312, 15.45925062, ..., 23.96888318,\n 18.56936995, 18.5055643 ],\n [15.43032887, 15.38152436, 15.72765387, ..., 24.00778199,\n 31.43343708, 31.39218212],\n ...,\n [15.46060449, 15.44254019, 16.22415641, ..., 22.5952011 ,\n 22.00937089, 21.92098704],\n [15.93206177, 15.91515816, 16.90845999, ..., 22.31407781,\n 25.07906851, 24.98810065],\n [15.49355372, 15.45764709, 17.06462284, ..., 23.43403333,\n 20.42975062, 20.34075804]],\n\n [[16.32471939, 16.30094022, 17.43097525, ..., 23.00163648,\n 26.80332891, 26.72146453],\n [17.33521783, 17.3134523 , 17.95562079, ..., 23.03445699,\n 21.59493384, 21.51720511],\n [17.54063129, 17.49845413, 17.65990214, ..., 22.36994524,\n 27.4844157 , 27.42609916],\n ...,\n [16.01830855, 15.99339059, 16.59490184, ..., 22.48199671,\n 27.81130108, 27.72439551],\n [17.4080136 , 17.34345667, 17.8184839 , ..., 22.84697905,\n 23.71841269, 23.66731893],\n [16.42590678, 16.38630959, 17.46862252, ..., 22.37950074,\n 22.93086901, 22.83010099]],\n\n [[14.73567007, 14.71271276, 16.69253713, ..., 23.71015323,\n 27.15856241, 27.08063007],\n [12.98946737, 12.96790069, 14.86294855, ..., 23.11624176,\n 28.56426266, 28.53871872],\n [18.6172853 , 18.59698247, 17.67934574, ..., 22.3743203 ,\n 21.99053511, 21.91634019],\n ...,\n [17.82254531, 17.77186059, 17.60227204, ..., 22.91981166,\n 26.42983377, 26.39747279],\n [16.01207202, 15.96297043, 17.00412761, ..., 24.1316353 ,\n 16.16548566, 16.05931968],\n [16.35202156, 16.29547816, 17.84546165, ..., 23.57569998,\n 25.46476917, 25.40356825]]])) to Variable.",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"File \u001b[0;32m~/miniforge3/envs/bambinos/lib/python3.11/site-packages/xarray/core/variable.py:126\u001b[0m, in \u001b[0;36mas_variable\u001b[0;34m(obj, name)\u001b[0m\n\u001b[1;32m 125\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m--> 126\u001b[0m obj \u001b[39m=\u001b[39m Variable(\u001b[39m*\u001b[39;49mobj)\n\u001b[1;32m 127\u001b[0m \u001b[39mexcept\u001b[39;00m (\u001b[39mTypeError\u001b[39;00m, \u001b[39mValueError\u001b[39;00m) \u001b[39mas\u001b[39;00m error:\n\u001b[1;32m 128\u001b[0m \u001b[39m# use .format() instead of % because it handles tuples consistently\u001b[39;00m\n",
"File \u001b[0;32m~/miniforge3/envs/bambinos/lib/python3.11/site-packages/xarray/core/variable.py:362\u001b[0m, in \u001b[0;36mVariable.__init__\u001b[0;34m(self, dims, data, attrs, encoding, fastpath)\u001b[0m\n\u001b[1;32m 361\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_data \u001b[39m=\u001b[39m as_compatible_data(data, fastpath\u001b[39m=\u001b[39mfastpath)\n\u001b[0;32m--> 362\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_dims \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_parse_dimensions(dims)\n\u001b[1;32m 363\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_attrs \u001b[39m=\u001b[39m \u001b[39mNone\u001b[39;00m\n",
"File \u001b[0;32m~/miniforge3/envs/bambinos/lib/python3.11/site-packages/xarray/core/variable.py:663\u001b[0m, in \u001b[0;36mVariable._parse_dimensions\u001b[0;34m(self, dims)\u001b[0m\n\u001b[1;32m 662\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mlen\u001b[39m(dims) \u001b[39m!=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mndim:\n\u001b[0;32m--> 663\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m 664\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mdimensions \u001b[39m\u001b[39m{\u001b[39;00mdims\u001b[39m}\u001b[39;00m\u001b[39m must have the same length as the \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 665\u001b[0m \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mnumber of data dimensions, ndim=\u001b[39m\u001b[39m{\u001b[39;00m\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mndim\u001b[39m}\u001b[39;00m\u001b[39m\"\u001b[39m\n\u001b[1;32m 666\u001b[0m )\n\u001b[1;32m 667\u001b[0m \u001b[39mreturn\u001b[39;00m dims\n",
"\u001b[0;31mValueError\u001b[0m: dimensions ('preds',) must have the same length as the number of data dimensions, ndim=3",
"\nDuring handling of the above exception, another exception occurred:\n",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[19], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m compare_DS \u001b[39m=\u001b[39m xr\u001b[39m.\u001b[39;49mDataset(\n\u001b[1;32m 2\u001b[0m \u001b[39m# can't do this because shape is not consistent with coords\u001b[39;49;00m\n\u001b[1;32m 3\u001b[0m \u001b[39m# {\u001b[39;49;00m\n\u001b[1;32m 4\u001b[0m \u001b[39m# \"mpg_preds\": ([\"preds\"], idata_new.posterior.mpg_mean.sel(chain=0, draw=0).values)\u001b[39;49;00m\n\u001b[1;32m 5\u001b[0m \u001b[39m# },\u001b[39;49;00m\n\u001b[1;32m 6\u001b[0m {\n\u001b[1;32m 7\u001b[0m \u001b[39m\"\u001b[39;49m\u001b[39mmpg_preds\u001b[39;49m\u001b[39m\"\u001b[39;49m: ([\u001b[39m\"\u001b[39;49m\u001b[39mpreds\u001b[39;49m\u001b[39m\"\u001b[39;49m], idata_new\u001b[39m.\u001b[39;49mposterior\u001b[39m.\u001b[39;49mmpg_mean\u001b[39m.\u001b[39;49mvalues)\n\u001b[1;32m 8\u001b[0m },\n\u001b[1;32m 9\u001b[0m coords\u001b[39m=\u001b[39;49m{\n\u001b[1;32m 10\u001b[0m \u001b[39m\"\u001b[39;49m\u001b[39mdrat\u001b[39;49m\u001b[39m\"\u001b[39;49m: ([\u001b[39m\"\u001b[39;49m\u001b[39mdrat\u001b[39;49m\u001b[39m\"\u001b[39;49m], drat),\n\u001b[1;32m 11\u001b[0m \u001b[39m\"\u001b[39;49m\u001b[39mam\u001b[39;49m\u001b[39m\"\u001b[39;49m: ([\u001b[39m\"\u001b[39;49m\u001b[39mam\u001b[39;49m\u001b[39m\"\u001b[39;49m], am),\n\u001b[1;32m 12\u001b[0m \u001b[39m\"\u001b[39;49m\u001b[39mhp\u001b[39;49m\u001b[39m\"\u001b[39;49m: ([\u001b[39m\"\u001b[39;49m\u001b[39mhp\u001b[39;49m\u001b[39m\"\u001b[39;49m], hp),\n\u001b[1;32m 13\u001b[0m }\n\u001b[1;32m 14\u001b[0m \u001b[39m# cant do this because the group by would not take into account the group levels of each covariate?\u001b[39;49;00m\n\u001b[1;32m 15\u001b[0m \u001b[39m# coords={\u001b[39;49;00m\n\u001b[1;32m 16\u001b[0m \u001b[39m# \"covariates\": ([\"drat\", \"am\", \"hp\"], comparisons_df.iloc[:, :3].values)\u001b[39;49;00m\n\u001b[1;32m 17\u001b[0m \u001b[39m# }\u001b[39;49;00m\n\u001b[1;32m 18\u001b[0m )\n",
"File \u001b[0;32m~/miniforge3/envs/bambinos/lib/python3.11/site-packages/xarray/core/dataset.py:612\u001b[0m, in \u001b[0;36mDataset.__init__\u001b[0;34m(self, data_vars, coords, attrs)\u001b[0m\n\u001b[1;32m 609\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(coords, Dataset):\n\u001b[1;32m 610\u001b[0m coords \u001b[39m=\u001b[39m coords\u001b[39m.\u001b[39mvariables\n\u001b[0;32m--> 612\u001b[0m variables, coord_names, dims, indexes, _ \u001b[39m=\u001b[39m merge_data_and_coords(\n\u001b[1;32m 613\u001b[0m data_vars, coords, compat\u001b[39m=\u001b[39;49m\u001b[39m\"\u001b[39;49m\u001b[39mbroadcast_equals\u001b[39;49m\u001b[39m\"\u001b[39;49m\n\u001b[1;32m 614\u001b[0m )\n\u001b[1;32m 616\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_attrs \u001b[39m=\u001b[39m \u001b[39mdict\u001b[39m(attrs) \u001b[39mif\u001b[39;00m attrs \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39melse\u001b[39;00m \u001b[39mNone\u001b[39;00m\n\u001b[1;32m 617\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_close \u001b[39m=\u001b[39m \u001b[39mNone\u001b[39;00m\n",
"File \u001b[0;32m~/miniforge3/envs/bambinos/lib/python3.11/site-packages/xarray/core/merge.py:564\u001b[0m, in \u001b[0;36mmerge_data_and_coords\u001b[0;34m(data_vars, coords, compat, join)\u001b[0m\n\u001b[1;32m 562\u001b[0m objects \u001b[39m=\u001b[39m [data_vars, coords]\n\u001b[1;32m 563\u001b[0m explicit_coords \u001b[39m=\u001b[39m coords\u001b[39m.\u001b[39mkeys()\n\u001b[0;32m--> 564\u001b[0m \u001b[39mreturn\u001b[39;00m merge_core(\n\u001b[1;32m 565\u001b[0m objects,\n\u001b[1;32m 566\u001b[0m compat,\n\u001b[1;32m 567\u001b[0m join,\n\u001b[1;32m 568\u001b[0m explicit_coords\u001b[39m=\u001b[39;49mexplicit_coords,\n\u001b[1;32m 569\u001b[0m indexes\u001b[39m=\u001b[39;49mIndexes(indexes, coords),\n\u001b[1;32m 570\u001b[0m )\n",
"File \u001b[0;32m~/miniforge3/envs/bambinos/lib/python3.11/site-packages/xarray/core/merge.py:744\u001b[0m, in \u001b[0;36mmerge_core\u001b[0;34m(objects, compat, join, combine_attrs, priority_arg, explicit_coords, indexes, fill_value)\u001b[0m\n\u001b[1;32m 740\u001b[0m coerced \u001b[39m=\u001b[39m coerce_pandas_values(objects)\n\u001b[1;32m 741\u001b[0m aligned \u001b[39m=\u001b[39m deep_align(\n\u001b[1;32m 742\u001b[0m coerced, join\u001b[39m=\u001b[39mjoin, copy\u001b[39m=\u001b[39m\u001b[39mFalse\u001b[39;00m, indexes\u001b[39m=\u001b[39mindexes, fill_value\u001b[39m=\u001b[39mfill_value\n\u001b[1;32m 743\u001b[0m )\n\u001b[0;32m--> 744\u001b[0m collected \u001b[39m=\u001b[39m collect_variables_and_indexes(aligned, indexes\u001b[39m=\u001b[39;49mindexes)\n\u001b[1;32m 745\u001b[0m prioritized \u001b[39m=\u001b[39m _get_priority_vars_and_indexes(aligned, priority_arg, compat\u001b[39m=\u001b[39mcompat)\n\u001b[1;32m 746\u001b[0m variables, out_indexes \u001b[39m=\u001b[39m merge_collected(\n\u001b[1;32m 747\u001b[0m collected, prioritized, compat\u001b[39m=\u001b[39mcompat, combine_attrs\u001b[39m=\u001b[39mcombine_attrs\n\u001b[1;32m 748\u001b[0m )\n",
"File \u001b[0;32m~/miniforge3/envs/bambinos/lib/python3.11/site-packages/xarray/core/merge.py:354\u001b[0m, in \u001b[0;36mcollect_variables_and_indexes\u001b[0;34m(list_of_mappings, indexes)\u001b[0m\n\u001b[1;32m 351\u001b[0m indexes_\u001b[39m.\u001b[39mpop(name, \u001b[39mNone\u001b[39;00m)\n\u001b[1;32m 352\u001b[0m append_all(coords_, indexes_)\n\u001b[0;32m--> 354\u001b[0m variable \u001b[39m=\u001b[39m as_variable(variable, name\u001b[39m=\u001b[39;49mname)\n\u001b[1;32m 355\u001b[0m \u001b[39mif\u001b[39;00m name \u001b[39min\u001b[39;00m indexes:\n\u001b[1;32m 356\u001b[0m append(name, variable, indexes[name])\n",
"File \u001b[0;32m~/miniforge3/envs/bambinos/lib/python3.11/site-packages/xarray/core/variable.py:129\u001b[0m, in \u001b[0;36mas_variable\u001b[0;34m(obj, name)\u001b[0m\n\u001b[1;32m 126\u001b[0m obj \u001b[39m=\u001b[39m Variable(\u001b[39m*\u001b[39mobj)\n\u001b[1;32m 127\u001b[0m \u001b[39mexcept\u001b[39;00m (\u001b[39mTypeError\u001b[39;00m, \u001b[39mValueError\u001b[39;00m) \u001b[39mas\u001b[39;00m error:\n\u001b[1;32m 128\u001b[0m \u001b[39m# use .format() instead of % because it handles tuples consistently\u001b[39;00m\n\u001b[0;32m--> 129\u001b[0m \u001b[39mraise\u001b[39;00m error\u001b[39m.\u001b[39m\u001b[39m__class__\u001b[39m(\n\u001b[1;32m 130\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mCould not convert tuple of form \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 131\u001b[0m \u001b[39m\"\u001b[39m\u001b[39m(dims, data[, attrs, encoding]): \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m 132\u001b[0m \u001b[39m\"\u001b[39m\u001b[39m{}\u001b[39;00m\u001b[39m to Variable.\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m.\u001b[39mformat(obj)\n\u001b[1;32m 133\u001b[0m )\n\u001b[1;32m 134\u001b[0m \u001b[39melif\u001b[39;00m utils\u001b[39m.\u001b[39mis_scalar(obj):\n\u001b[1;32m 135\u001b[0m obj \u001b[39m=\u001b[39m Variable([], obj)\n",
"\u001b[0;31mValueError\u001b[0m: Could not convert tuple of form (dims, data[, attrs, encoding]): (['preds'], array([[[19.59478712, 19.54856056, 19.23140441, ..., 22.49626788,\n 22.51163257, 22.44590673],\n [18.27013994, 18.22625734, 18.58204345, ..., 23.3491659 ,\n 19.44539658, 19.35357745],\n [18.85953525, 18.81616193, 18.31847112, ..., 23.20774708,\n 22.53151746, 22.45055377],\n ...,\n [16.95947443, 16.91816056, 17.61175804, ..., 22.68251523,\n 22.94911059, 22.84984865],\n [16.75874405, 16.71122187, 17.75769402, ..., 23.31604767,\n 24.31327371, 24.22783439],\n [15.98276926, 15.93614123, 16.72108478, ..., 21.95560998,\n 17.56669211, 17.48011585]],\n\n [[15.13806506, 15.07475069, 16.57548916, ..., 22.93274638,\n 21.45973418, 21.40293702],\n [14.03187854, 13.96100312, 15.45925062, ..., 23.96888318,\n 18.56936995, 18.5055643 ],\n [15.43032887, 15.38152436, 15.72765387, ..., 24.00778199,\n 31.43343708, 31.39218212],\n ...,\n [15.46060449, 15.44254019, 16.22415641, ..., 22.5952011 ,\n 22.00937089, 21.92098704],\n [15.93206177, 15.91515816, 16.90845999, ..., 22.31407781,\n 25.07906851, 24.98810065],\n [15.49355372, 15.45764709, 17.06462284, ..., 23.43403333,\n 20.42975062, 20.34075804]],\n\n [[16.32471939, 16.30094022, 17.43097525, ..., 23.00163648,\n 26.80332891, 26.72146453],\n [17.33521783, 17.3134523 , 17.95562079, ..., 23.03445699,\n 21.59493384, 21.51720511],\n [17.54063129, 17.49845413, 17.65990214, ..., 22.36994524,\n 27.4844157 , 27.42609916],\n ...,\n [16.01830855, 15.99339059, 16.59490184, ..., 22.48199671,\n 27.81130108, 27.72439551],\n [17.4080136 , 17.34345667, 17.8184839 , ..., 22.84697905,\n 23.71841269, 23.66731893],\n [16.42590678, 16.38630959, 17.46862252, ..., 22.37950074,\n 22.93086901, 22.83010099]],\n\n [[14.73567007, 14.71271276, 16.69253713, ..., 23.71015323,\n 27.15856241, 27.08063007],\n [12.98946737, 12.96790069, 14.86294855, ..., 23.11624176,\n 28.56426266, 28.53871872],\n [18.6172853 , 18.59698247, 17.67934574, ..., 22.3743203 ,\n 21.99053511, 21.91634019],\n ...,\n [17.82254531, 17.77186059, 17.60227204, ..., 22.91981166,\n 26.42983377, 26.39747279],\n [16.01207202, 15.96297043, 17.00412761, ..., 24.1316353 ,\n 16.16548566, 16.05931968],\n [16.35202156, 16.29547816, 17.84546165, ..., 23.57569998,\n 25.46476917, 25.40356825]]])) to Variable."
]
}
],
"source": [
"compare_DS = xr.Dataset(\n",
" # can't do this because shape is not consistent with coords\n",
" # {\n",
" # \"mpg_preds\": ([\"preds\"], idata_new.posterior.mpg_mean.sel(chain=0, draw=0).values)\n",
" # },\n",
" {\n",
" \"mpg_preds\": ([\"preds\"], idata_new.posterior.mpg_mean.values)\n",
" },\n",
" coords={\n",
" \"drat\": ([\"drat\"], drat),\n",
" \"am\": ([\"am\"], am),\n",
" \"hp\": ([\"hp\"], hp),\n",
" }\n",
"# cant do this because the group by would not take into account the group levels of each covariate?\n",
"# coords={\n",
"# \"covariates\": ([\"drat\", \"am\", \"hp\"], comparisons_df.iloc[:, :3].values)\n",
"# }\n",
")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using Pandas\n",
"\n",
"Since posterior is shape (4, 1000, 20) where 20 is the row number of the data used for prediction, we can create a DataFrame `posterior_with_obs` with the posterior samples and data used for performing predictions."
]
},
{
"cell_type": "code",
"execution_count": 31,
"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>chain</th>\n",
" <th>draw</th>\n",
" <th>mpg_obs</th>\n",
" <th>mpg_mean</th>\n",
" <th>am</th>\n",
" <th>drat</th>\n",
" <th>hp</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>19.594787</td>\n",
" <td>0</td>\n",
" <td>2.760</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>19.548561</td>\n",
" <td>0</td>\n",
" <td>2.760</td>\n",
" <td>147</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>2</td>\n",
" <td>19.231404</td>\n",
" <td>0</td>\n",
" <td>3.080</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>3</td>\n",
" <td>19.175185</td>\n",
" <td>0</td>\n",
" <td>3.080</td>\n",
" <td>147</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>4</td>\n",
" <td>18.533028</td>\n",
" <td>0</td>\n",
" <td>3.695</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79995</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>15</td>\n",
" <td>23.168502</td>\n",
" <td>1</td>\n",
" <td>3.695</td>\n",
" <td>147</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79996</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>16</td>\n",
" <td>23.633115</td>\n",
" <td>1</td>\n",
" <td>3.920</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79997</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>17</td>\n",
" <td>23.575700</td>\n",
" <td>1</td>\n",
" <td>3.920</td>\n",
" <td>147</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79998</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>18</td>\n",
" <td>25.464769</td>\n",
" <td>1</td>\n",
" <td>4.930</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79999</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>19</td>\n",
" <td>25.403568</td>\n",
" <td>1</td>\n",
" <td>4.930</td>\n",
" <td>147</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>80000 rows × 7 columns</p>\n",
"</div>"
],
"text/plain": [
" chain draw mpg_obs mpg_mean am drat hp\n",
"0 0 0 0 19.594787 0 2.760 146\n",
"1 0 0 1 19.548561 0 2.760 147\n",
"2 0 0 2 19.231404 0 3.080 146\n",
"3 0 0 3 19.175185 0 3.080 147\n",
"4 0 0 4 18.533028 0 3.695 146\n",
"... ... ... ... ... .. ... ...\n",
"79995 3 999 15 23.168502 1 3.695 147\n",
"79996 3 999 16 23.633115 1 3.920 146\n",
"79997 3 999 17 23.575700 1 3.920 147\n",
"79998 3 999 18 25.464769 1 4.930 146\n",
"79999 3 999 19 25.403568 1 4.930 147\n",
"\n",
"[80000 rows x 7 columns]"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"comparisons_df.index.name = \"mpg_obs\"\n",
"comparisons_df = comparisons_df.iloc[:, :3].reset_index()\n",
"\n",
"idata_df = idata_new.posterior.mpg_mean.to_dataframe().reset_index()\n",
"\n",
"posterior_with_obs = pd.merge(\n",
" idata_df,\n",
" comparisons_df,\n",
" how=\"left\",\n",
" on=\"mpg_obs\"\n",
")\n",
"\n",
"posterior_with_obs"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-2.90745548, 2.81403525])"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAoEAAAEXCAYAAADShJ57AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAn3ElEQVR4nO3daXBU15338V8LSQ0i0jXQlloyYgkPD4JBjAkOAhxGeMyiJDJ4qcKAkRcosB8CWBgKQzwTBKlHYvGAa4YQXMRl4gwTTcwyY3syCjIEFUQCBB4VAhlSTrCHrRFgqSWMaAl0nhc83KHRTtQguN9PVb/oc//nnnPPAfHj9iKXMcYIAAAAjhJ2rycAAACAu48QCAAA4ECEQAAAAAciBAIAADgQIRAAAMCBCIEAAAAORAgEAABwIEIgAACAAxECAQAAHIgQCAAA4EAhDYEVFRXKyMiQZVmyLEsZGRmqrKxsto8xRllZWUpISFCXLl00ZswYHTt2LKgmEAho3rx58ng86tq1qyZOnKjTp0+3aexLly4pLS1NCQkJcrvdSkxM1Ny5c1VVVRV0nt/85jd69NFHFRUVpd69e2vNmjV/0ZoAAAB0BCENgdOmTVNJSYny8vKUl5enkpISZWRkNNtn9erVWrt2rdavX6/i4mJ5vV6NGzdO1dXVdk1mZqZ27Nih3Nxc7du3T5cvX1Z6erquX7/e6rHDwsI0adIkffTRR/rjH/+ozZs369NPP9Vrr71m1/znf/6nXnjhBb322ms6evSoNmzYYM8NAADgfuYyxphQnPjzzz/XoEGDtH//fqWkpEiS9u/fr5EjR+r48eMaMGBAgz7GGCUkJCgzM1NvvvmmpBt3/eLi4rRq1Sq9+uqr8vv9evjhh/WrX/1Kzz//vCTp7NmzSkxM1G9/+1tNmDDhjsaWpH/8x3/UmjVrdOrUKUk3gmRdXZ0+/PBDu+add97RP/zDP+i///u/5XK5WlyH+vp6nT17VtHR0a2qBwAAuFPGGFVXVyshIUFhYc3f6wsP1SSKiopkWZYdwiRpxIgRsixLhYWFjQaxkydPyufzafz48Xab2+1WamqqCgsL9eqrr+rw4cOqq6sLqklISNDgwYNVWFioCRMm3NHYZ8+e1fbt25Wammq3BQIBRUVFBdV16dJFp0+f1ldffaU+ffo0OE8gEFAgELCfnzlzRoMGDWphtQAAANrPqVOn1LNnz2ZrQhYCfT6fYmNjG7THxsbK5/M12UeS4uLigtrj4uL01Vdf2TWRkZHq1q1bg5qb/dsy9tSpU/Xv//7vqqmp0VNPPaVf/OIX9rEJEyZowYIFevnll/XEE0/oiy++0DvvvCNJOnfuXKMhMCcnR8uXL2/QfurUKcXExDR63QBwqyu11zT8/+6SJB1860lFRYbsRzWAB0xVVZUSExMVHR3dYm2bf7JkZWU1GnJuVVxcLEmNvvxpjGnxZdHbj7emz+01rR173bp1WrZsmU6cOKEf//jHeuONN7RhwwZJ0qxZs/SnP/1J6enpqqurU0xMjF5//XVlZWWpU6dOjc5j6dKleuONN+znNzcjJiaGEAigVcJrrynMfeNViJiYGEIggDZrzVvQ2vyTZe7cuZoyZUqzNX369NGRI0d0/vz5BscuXLjQ4E7fTV6vV9KNO3nx8fF2e3l5ud3H6/WqtrZWFRUVQXcDy8vLNWrUKLumtWN7vV55vV4lJSWpR48eGj16tP7+7/9e8fHxcrlcWrVqlbKzs+Xz+fTwww9r165d9jU2xu12y+12N7U0AAAAHUKbPx3s8XiUlJTU7KNz584aOXKk/H6/Dh48aPc9cOCA/H6/HdZu17dvX3m9XuXn59tttbW1KigosPsMGzZMERERQTXnzp3T0aNH7Zo7GVu6cadQUtB7+iSpU6dOeuSRRxQZGalf//rXGjlyZKMvNwMAANwvQvYaw8CBA5WWlqZZs2bp3XfflSTNnj1b6enpQR/MSEpKUk5Ojp555hm5XC5lZmYqOztb/fv3V//+/ZWdna2oqChNmzZNkmRZlmbOnKmFCxeqR48e6t69uxYtWqTk5GSNHTu21WP/9re/1fnz5/Xd735X3/rWt1RWVqbFixfr8ccft+/yXbx4UVu3btWYMWN09epVvf/++/rwww9VUFAQqmUDAAC4K0L6RpMtW7Zo/vz59id5J06c2OA79k6cOCG/328/X7x4sWpqajRnzhxVVFQoJSVFO3fuDHqD47p16xQeHq7JkyerpqZGTz75pDZv3hz0Pr2Wxu7SpYs2bdqkBQsWKBAIKDExUc8++6yWLFkSNL9f/vKXWrRokYwxGjlypPbs2aPhw4e33yIBAADcAyH7nkDcUFVVJcuy5Pf7+WAIgFa5UntNg37yO0lS2YoJfDAEQKu1JXfwu4MBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADhTQEVlRUKCMjQ5ZlybIsZWRkqLKystk+xhhlZWUpISFBXbp00ZgxY3Ts2LGgmkAgoHnz5snj8ahr166aOHGiTp8+fcdjX7p0ST179pTL5WpQU1paqtTUVHXp0kWPPPKIVqxYIWNMW5cCAACgQwlpCJw2bZpKSkqUl5envLw8lZSUKCMjo9k+q1ev1tq1a7V+/XoVFxfL6/Vq3Lhxqq6utmsyMzO1Y8cO5ebmat++fbp8+bLS09N1/fr1Oxp75syZGjJkSIP2qqoqjRs3TgkJCSouLtY//dM/6e2339batWvvcEUAAAA6CBMiZWVlRpLZv3+/3VZUVGQkmePHjzfap76+3ni9XrNy5Uq77erVq8ayLLNx40ZjjDGVlZUmIiLC5Obm2jVnzpwxYWFhJi8vr81jb9iwwaSmpppdu3YZSaaioiLomGVZ5urVq3ZbTk6OSUhIMPX19a1aB7/fbyQZv9/fqnoA+CZQZ3q/+Ynp/eYn5ptA3b2eDoD7SFtyR8juBBYVFcmyLKWkpNhtI0aMkGVZKiwsbLTPyZMn5fP5NH78eLvN7XYrNTXV7nP48GHV1dUF1SQkJGjw4MF2TWvHLisr04oVK/TBBx8oLKzhUhQVFSk1NVVut9tumzBhgs6ePasvv/yy0WsIBAKqqqoKegAAAHQ0IQuBPp9PsbGxDdpjY2Pl8/ma7CNJcXFxQe1xcXH2MZ/Pp8jISHXr1q3ZmpbGDgQCmjp1qtasWaNevXo1OZ/G5nLrXG+Xk5Njvw/RsiwlJiY2WgcAAHAvtTkEZmVlyeVyNfs4dOiQJMnlcjXob4xptP1Wtx9vTZ/ba1oae+nSpRo4cKCmT5/e5rk0df6b5/X7/fbj1KlTzZ4fAADgXghva4e5c+dqypQpzdb06dNHR44c0fnz5xscu3DhQoO7azd5vV5JN+6yxcfH2+3l5eV2H6/Xq9raWlVUVATdDSwvL9eoUaPsmpbG3r17t0pLS7V161ZJ/xPuPB6P3nrrLS1fvlxer7fBHb/y8nJJDe9W3uR2u4NePgYAAOiI2hwCPR6PPB5Pi3UjR46U3+/XwYMHNXz4cEnSgQMH5Pf77bB2u759+8rr9So/P19Dhw6VJNXW1qqgoECrVq2SJA0bNkwRERHKz8/X5MmTJUnnzp3T0aNHtXr16laPvW3bNtXU1NhjFxcXa8aMGdq7d6/69etnn+fHP/6xamtrFRkZKUnauXOnEhIS1KdPnzatGwAAQIcSyk+opKWlmSFDhpiioiJTVFRkkpOTTXp6elDNgAEDzPbt2+3nK1euNJZlme3bt5vS0lIzdepUEx8fb6qqquya1157zfTs2dN8+umn5rPPPjN/+7d/a/76r//aXLt2rU1j3+r3v/99g08HV1ZWmri4ODN16lRTWlpqtm/fbmJiYszbb7/d6jXg08EA2opPBwO4U23JHW2+E9gWW7Zs0fz58+1P8k6cOFHr168Pqjlx4oT8fr/9fPHixaqpqdGcOXNUUVGhlJQU7dy5U9HR0XbNunXrFB4ersmTJ6umpkZPPvmkNm/erE6dOrVp7JZYlqX8/Hz96Ec/0mOPPaZu3brpjTfe0BtvvNHmtQAAAOhIXMbw6y9CqaqqSpZlye/3KyYm5l5PB8B94ErtNQ36ye8kSWUrJigqMqT/XwfwAGlL7uB3BwMAADgQIRAAAMCBCIEAAAAORAgEAABwIEIgAACAAxECAQAAHIgQCAAA4ECEQAAAAAciBAIAADgQIRAAAMCBCIEAAAAORAgEAABwIEIgAACAAxECAQAAHIgQCAAA4ECEQAAAAAciBAIAADgQIRAAAMCBCIEAAAAORAgEAABwIEIgAACAAxECAQAAHIgQCAAA4ECEQAAAAAciBAIAADgQIRAAAMCBCIEAAAAORAgEAABwIEIgAACAAxECAQAAHIgQCAAA4ECEQAAAAAciBAIAADgQIRAAAMCBCIEAAAAORAgEAABwIEIgAACAAxECAQAAHIgQCAAA4ECEQAAAAAcKaQisqKhQRkaGLMuSZVnKyMhQZWVls32MMcrKylJCQoK6dOmiMWPG6NixY0E1gUBA8+bNk8fjUdeuXTVx4kSdPn36jse+dOmSevbsKZfLFVRz9epVvfzyy0pOTlZ4eLiefvrpO1gFAACAjiekIXDatGkqKSlRXl6e8vLyVFJSooyMjGb7rF69WmvXrtX69etVXFwsr9ercePGqbq62q7JzMzUjh07lJubq3379uny5ctKT0/X9evX72jsmTNnasiQIQ3ar1+/ri5dumj+/PkaO3bsHa4CAABAB2RCpKyszEgy+/fvt9uKioqMJHP8+PFG+9TX1xuv12tWrlxpt129etVYlmU2btxojDGmsrLSREREmNzcXLvmzJkzJiwszOTl5bV57A0bNpjU1FSza9cuI8lUVFQ0OreXXnrJTJo0qU1rYIwxfr/fSDJ+v7/NfQE40zeBOtP7zU9M7zc/Md8E6u71dADcR9qSO0J2J7CoqEiWZSklJcVuGzFihCzLUmFhYaN9Tp48KZ/Pp/Hjx9ttbrdbqampdp/Dhw+rrq4uqCYhIUGDBw+2a1o7dllZmVasWKEPPvhAYWHtsxSBQEBVVVVBDwAAgI4mZCHQ5/MpNja2QXtsbKx8Pl+TfSQpLi4uqD0uLs4+5vP5FBkZqW7dujVb09LYgUBAU6dO1Zo1a9SrV682Xl3TcnJy7PchWpalxMTEdjs3AABAe2lzCMzKypLL5Wr2cejQIUmSy+Vq0N8Y02j7rW4/3po+t9e0NPbSpUs1cOBATZ8+vdnzttXSpUvl9/vtx6lTp9r1/AAAAO0hvK0d5s6dqylTpjRb06dPHx05ckTnz59vcOzChQsN7vTd5PV6Jd24kxcfH2+3l5eX2328Xq9qa2tVUVERdDewvLxco0aNsmtaGnv37t0qLS3V1q1bJd0IiJLk8Xj01ltvafny5c1eY1Pcbrfcbvcd9QUAALhb2hwCPR6PPB5Pi3UjR46U3+/XwYMHNXz4cEnSgQMH5Pf77bB2u759+8rr9So/P19Dhw6VJNXW1qqgoECrVq2SJA0bNkwRERHKz8/X5MmTJUnnzp3T0aNHtXr16laPvW3bNtXU1NhjFxcXa8aMGdq7d6/69evX1mUBAAC4r7Q5BLbWwIEDlZaWplmzZundd9+VJM2ePVvp6ekaMGCAXZeUlKScnBw988wzcrlcyszMVHZ2tvr376/+/fsrOztbUVFRmjZtmiTJsizNnDlTCxcuVI8ePdS9e3ctWrRIycnJ9te4tGbs24PexYsX7b4PPfSQ3V5WVqba2lp9/fXXqq6uVklJiSTp0Ucfbfc1AwAAuFtCFgIlacuWLZo/f779Sd6JEydq/fr1QTUnTpyQ3++3ny9evFg1NTWaM2eOKioqlJKSop07dyo6OtquWbduncLDwzV58mTV1NToySef1ObNm9WpU6c2jd0aP/jBD/TVV1/Zz2/eobz58jEAAMD9yGVIMyFVVVUly7Lk9/sVExNzr6cD4D5wpfaaBv3kd5KkshUTFBUZ0v+vA3iAtCV38LuDAQAAHIgQCAAdzPX6/3mB5sCfvw56DgDthRAIAB1I3tFzGru2wH7+yuZifW/VbuUdPXcPZwXgQUQIBIAOIu/oOf2ff/5M56sCQe0+/1X9n3/+jCAIoF0RAgGgA7heb7T84zI19sLvzbblH5fx0jCAdkMIBIAO4ODJr3XOf7XJ40bSOf9VHTz59d2bFIAHGiEQADqA8uqmA+Cd1AFASwiBANABxEZ3btc6AGgJIRAAOoDhfbsr3uosVxPHXZLirc4a3rf73ZwWgAcYIRAAOoBOYS4te2pQo8duBsNlTw1Sp7CmYiIAtA0hEAA6iLTB8fr59O8oLsYd1O61Ouvn07+jtMHx92hmAB5E/EJKAOhA0gbH6/H/5VFy1k5J0vsvf1d/878f5g4ggHbHnUAA6GBuDXwp3+5OAAQQEoRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIFCGgIrKiqUkZEhy7JkWZYyMjJUWVnZbB9jjLKyspSQkKAuXbpozJgxOnbsWFBNIBDQvHnz5PF41LVrV02cOFGnT5++47EvXbqknj17yuVyBdXs2bNHkyZNUnx8vLp27apHH31UW7ZsuZOlAAAA6FBCGgKnTZumkpIS5eXlKS8vTyUlJcrIyGi2z+rVq7V27VqtX79excXF8nq9GjdunKqrq+2azMxM7dixQ7m5udq3b58uX76s9PR0Xb9+/Y7GnjlzpoYMGdKgvbCwUEOGDNG2bdt05MgRzZgxQy+++KI+/vjjO1wRAACADsKESFlZmZFk9u/fb7cVFRUZSeb48eON9qmvrzder9esXLnSbrt69aqxLMts3LjRGGNMZWWliYiIMLm5uXbNmTNnTFhYmMnLy2vz2Bs2bDCpqalm165dRpKpqKho9rp+8IMfmFdeeaV1i2CM8fv9RpLx+/2t7gPA2b4J1Jneb35ier/5ifkmUHevpwPgPtKW3BGyO4FFRUWyLEspKSl224gRI2RZlgoLCxvtc/LkSfl8Po0fP95uc7vdSk1NtfscPnxYdXV1QTUJCQkaPHiwXdPascvKyrRixQp98MEHCgtr3VL4/X517969yeOBQEBVVVVBDwAAgI4mZCHQ5/MpNja2QXtsbKx8Pl+TfSQpLi4uqD0uLs4+5vP5FBkZqW7dujVb09LYgUBAU6dO1Zo1a9SrV69WXdPWrVtVXFysV155pcmanJwc+32IlmUpMTGxVecGAAC4m9ocArOysuRyuZp9HDp0SJLkcrka9DfGNNp+q9uPt6bP7TUtjb106VINHDhQ06dPb/a8N+3Zs0cvv/yyNm3apL/6q79qsm7p0qXy+/3249SpU606PwAAwN0U3tYOc+fO1ZQpU5qt6dOnj44cOaLz5883OHbhwoUGd/pu8nq9km7cyYuPj7fby8vL7T5er1e1tbWqqKgIuhtYXl6uUaNG2TUtjb17926VlpZq69atkm4EREnyeDx66623tHz5crtfQUGBnnrqKa1du1Yvvvhis9fudrvldrubrQEAALjX2hwCPR6PPB5Pi3UjR46U3+/XwYMHNXz4cEnSgQMH5Pf77bB2u759+8rr9So/P19Dhw6VJNXW1qqgoECrVq2SJA0bNkwRERHKz8/X5MmTJUnnzp3T0aNHtXr16laPvW3bNtXU1NhjFxcXa8aMGdq7d6/69etnt+/Zs0fp6elatWqVZs+e3aa1AgAA6KjaHAJba+DAgUpLS9OsWbP07rvvSpJmz56t9PR0DRgwwK5LSkpSTk6OnnnmGblcLmVmZio7O1v9+/dX//79lZ2draioKE2bNk2SZFmWZs6cqYULF6pHjx7q3r27Fi1apOTkZI0dO7bVY98a9CTp4sWLdt+HHnpI0o0A+MMf/lCvv/66nnvuOfv9hJGRkc1+OAQAAKCjC+n3BG7ZskXJyckaP368xo8fryFDhuhXv/pVUM2JEyfk9/vt54sXL1ZmZqbmzJmjxx57TGfOnNHOnTsVHR1t16xbt05PP/20Jk+erMcff1xRUVH6+OOP1alTpzaN3ZLNmzfrypUrysnJUXx8vP149tln73BFAAAAOgaXuflmOIREVVWVLMuS3+9XTEzMvZ4OgPvAldprGvST30mSylZMUFRkyF60AfCAaUvu4HcHAwAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAABwppCKyoqFBGRoYsy5JlWcrIyFBlZWWzfYwxysrKUkJCgrp06aIxY8bo2LFjQTWBQEDz5s2Tx+NR165dNXHiRJ0+ffqOx7506ZJ69uwpl8sVVHPixAk98cQTiouLU+fOnfXtb39bf/d3f6e6uro7WQ4AAIAOI6QhcNq0aSopKVFeXp7y8vJUUlKijIyMZvusXr1aa9eu1fr161VcXCyv16tx48apurrarsnMzNSOHTuUm5urffv26fLly0pPT9f169fvaOyZM2dqyJAhDdojIiL04osvaufOnTpx4oTeeecdbdq0ScuWLbvDFQEAAOggTIiUlZUZSWb//v12W1FRkZFkjh8/3mif+vp64/V6zcqVK+22q1evGsuyzMaNG40xxlRWVpqIiAiTm5tr15w5c8aEhYWZvLy8No+9YcMGk5qaanbt2mUkmYqKimava8GCBeZ73/te6xbBGOP3+40k4/f7W90HgLN9E6gzvd/8xPR+8xPzTaDuXk8HwH2kLbkjZHcCi4qKZFmWUlJS7LYRI0bIsiwVFhY22ufkyZPy+XwaP3683eZ2u5Wammr3OXz4sOrq6oJqEhISNHjwYLumtWOXlZVpxYoV+uCDDxQW1vJSfPHFF8rLy1NqamqTNYFAQFVVVUEPAACAjiZkIdDn8yk2NrZBe2xsrHw+X5N9JCkuLi6oPS4uzj7m8/kUGRmpbt26NVvT0tiBQEBTp07VmjVr1KtXr2avZdSoUercubP69++v0aNHa8WKFU3W5uTk2O9DtCxLiYmJzZ4bAADgXmhzCMzKypLL5Wr2cejQIUmSy+Vq0N8Y02j7rW4/3po+t9e0NPbSpUs1cOBATZ8+vdnzStK//uu/6rPPPtO//Mu/6D/+4z/09ttvN1m7dOlS+f1++3Hq1KkWzw8AAHC3hbe1w9y5czVlypRma/r06aMjR47o/PnzDY5duHChwZ2+m7xer6Qbd/Li4+Pt9vLycruP1+tVbW2tKioqgu4GlpeXa9SoUXZNS2Pv3r1bpaWl2rp1q6QbAVGSPB6P3nrrLS1fvtzud/Nu3qBBg3T9+nXNnj1bCxcuVKdOnRqM4Xa75Xa7m1oaAACADqHNIdDj8cjj8bRYN3LkSPn9fh08eFDDhw+XJB04cEB+v98Oa7fr27evvF6v8vPzNXToUElSbW2tCgoKtGrVKknSsGHDFBERofz8fE2ePFmSdO7cOR09elSrV69u9djbtm1TTU2NPXZxcbFmzJihvXv3ql+/fk1elzFGdXV1dmgEAAC4H7U5BLbWwIEDlZaWplmzZundd9+VJM2ePVvp6ekaMGCAXZeUlKScnBw988wzcrlcyszMVHZ2tvr376/+/fsrOztbUVFRmjZtmiTJsizNnDlTCxcuVI8ePdS9e3ctWrRIycnJGjt2bKvHvj3oXbx40e770EMPSZK2bNmiiIgIJScny+126/Dhw1q6dKmef/55hYeHbOkAAABCLqRJZsuWLZo/f779Sd6JEydq/fr1QTUnTpyQ3++3ny9evFg1NTWaM2eOKioqlJKSop07dyo6OtquWbduncLDwzV58mTV1NToySef1ObNm4Nenm3N2C0JDw/XqlWr9Mc//lHGGPXu3Vs/+tGPtGDBgjavBQAAQEfiMryuGVJVVVWyLEt+v18xMTH3ejoA7gNXaq9p0E9+J0kqWzFBUZG88gCgddqSO/jdwQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHAgQiAAAIADEQIBAAAciBAIAADgQIRAAAAAByIEAgAAOBAhEAAAwIEIgQAAAA5ECAQAAHCg8Hs9gQedMUaSVFVVdY9nAuB+caX2muoDVyTd+NlxLZIf1QBa52beuJk/muMyranCHTt9+rQSExPv9TQAAICDnDp1Sj179my2hhAYYvX19Tp79qyio6Plcrnu9XQ6rKqqKiUmJurUqVOKiYm519NxNPaiY2AfOg72omNgH1rHGKPq6molJCQoLKz5d/3xGkOIhYWFtZjE8T9iYmL4y91BsBcdA/vQcbAXHQP70DLLslpVxwdDAAAAHIgQCAAA4ECEQHQIbrdby5Ytk9vtvtdTcTz2omNgHzoO9qJjYB/aHx8MAQAAcCDuBAIAADgQIRAAAMCBCIEAAAAORAgEAABwIEIg7oqKigplZGTIsixZlqWMjAxVVlY228cYo6ysLCUkJKhLly4aM2aMjh071mTt97//fblcLv3bv/1b+1/AAyQUe/H1119r3rx5GjBggKKiotSrVy/Nnz9ffr8/xFdzf9mwYYP69u2rzp07a9iwYdq7d2+z9QUFBRo2bJg6d+6sb3/729q4cWODmm3btmnQoEFyu90aNGiQduzYEarpPzDaex82bdqk0aNHq1u3burWrZvGjh2rgwcPhvISHhih+DtxU25urlwul55++ul2nvUDxAB3QVpamhk8eLApLCw0hYWFZvDgwSY9Pb3ZPitXrjTR0dFm27ZtprS01Dz//PMmPj7eVFVVNahdu3at+f73v28kmR07doToKh4ModiL0tJS8+yzz5qPPvrIfPHFF2bXrl2mf//+5rnnnrsbl3RfyM3NNREREWbTpk2mrKzMvP7666Zr167mq6++arT+z3/+s4mKijKvv/66KSsrM5s2bTIRERFm69atdk1hYaHp1KmTyc7ONp9//rnJzs424eHhZv/+/Xfrsu47odiHadOmmZ/97Gfmv/7rv8znn39uXnnlFWNZljl9+vTduqz7Uij24qYvv/zSPPLII2b06NFm0qRJIb6S+xchECFXVlZmJAX9w1RUVGQkmePHjzfap76+3ni9XrNy5Uq77erVq8ayLLNx48ag2pKSEtOzZ09z7tw5QmALQr0Xt/rNb35jIiMjTV1dXftdwH1s+PDh5rXXXgtqS0pKMkuWLGm0fvHixSYpKSmo7dVXXzUjRoywn0+ePNmkpaUF1UyYMMFMmTKlnWb94AnFPtzu2rVrJjo62vzyl7/8yyf8AAvVXly7ds08/vjj5he/+IV56aWXCIHN4OVghFxRUZEsy1JKSordNmLECFmWpcLCwkb7nDx5Uj6fT+PHj7fb3G63UlNTg/pcuXJFU6dO1fr16+X1ekN3EQ+IUO7F7fx+v2JiYhQezq8or62t1eHDh4PWUJLGjx/f5BoWFRU1qJ8wYYIOHTqkurq6Zmua2xcnC9U+3O7KlSuqq6tT9+7d22fiD6BQ7sWKFSv08MMPa+bMme0/8QcMIRAh5/P5FBsb26A9NjZWPp+vyT6SFBcXF9QeFxcX1GfBggUaNWqUJk2a1I4zfnCFci9udenSJf30pz/Vq6+++hfO+MFw8eJFXb9+vU1r6PP5Gq2/du2aLl682GxNU+d0ulDtw+2WLFmiRx55RGPHjm2fiT+AQrUXf/jDH/Tee+9p06ZNoZn4A4YQiDuWlZUll8vV7OPQoUOSJJfL1aC/MabR9lvdfvzWPh999JF2796td955p30u6D52r/fiVlVVVfrhD3+oQYMGadmyZX/BVT14WruGzdXf3t7WcyI0+3DT6tWr9etf/1rbt29X586d22G2D7b23Ivq6mpNnz5dmzZtksfjaf/JPoB4nQZ3bO7cuZoyZUqzNX369NGRI0d0/vz5BscuXLjQ4H91N918adfn8yk+Pt5uLy8vt/vs3r1bf/rTn/TQQw8F9X3uuec0evRo7dmzpw1Xc3+713txU3V1tdLS0vStb31LO3bsUERERFsv5YHk8XjUqVOnBnc4GlvDm7xeb6P14eHh6tGjR7M1TZ3T6UK1Dze9/fbbys7O1qeffqohQ4a07+QfMKHYi2PHjunLL7/UU089ZR+vr6+XJIWHh+vEiRPq169fO1/J/Y07gbhjHo9HSUlJzT46d+6skSNHyu/3B31lwoEDB+T3+zVq1KhGz923b195vV7l5+fbbbW1tSooKLD7LFmyREeOHFFJSYn9kKR169bp/fffD92Fd0D3ei+kG3cAx48fr8jISH300UfcBblFZGSkhg0bFrSGkpSfn9/kuo8cObJB/c6dO/XYY4/Z4bqpmqbO6XSh2gdJWrNmjX76058qLy9Pjz32WPtP/gETir1ISkpSaWlp0L8JEydO1BNPPKGSkhIlJiaG7HruW/foAylwmLS0NDNkyBBTVFRkioqKTHJycoOvJRkwYIDZvn27/XzlypXGsiyzfft2U1paaqZOndrkV8TcJD4d3KJQ7EVVVZVJSUkxycnJ5osvvjDnzp2zH9euXbur19dR3fw6jPfee8+UlZWZzMxM07VrV/Pll18aY4xZsmSJycjIsOtvfh3GggULTFlZmXnvvfcafB3GH/7wB9OpUyezcuVK8/nnn5uVK1fyFTEtCMU+rFq1ykRGRpqtW7cG/dmvrq6+69d3PwnFXtyOTwc3jxCIu+LSpUvmhRdeMNHR0SY6Otq88MILpqKiIqhGknn//fft5/X19WbZsmXG6/Uat9tt/uZv/saUlpY2Ow4hsGWh2Ivf//73RlKjj5MnT96dC7sP/OxnPzO9e/c2kZGR5jvf+Y4pKCiwj7300ksmNTU1qH7Pnj1m6NChJjIy0vTp08f8/Oc/b3DODz/80AwYMMBERESYpKQks23btlBfxn2vvfehd+/ejf7ZX7Zs2V24mvtbKP5O3IoQ2DyXMf//XZUAAABwDN4TCAAA4ECEQAAAAAciBAIAADgQIRAAAMCBCIEAAAAORAgEAABwIEIgAACAAxECAQAAHIgQCAAA4ECEQAAAAAciBAIAADgQIRAAAMCB/h9KGb3q3qLzUwAAAABJRU5ErkJggg==",
"text/plain": [
"<Figure size 700x300 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# only look at am=0 and drat=2.76\n",
"# not showing entire vline for some reason\n",
"diffs = posterior_with_obs.query(\"am == 0 and drat == 2.760\")[\"mpg_mean\"].diff().dropna().reset_index(drop=True)\n",
"mean_diff = diffs.mean()\n",
"hdi = az.hdi(diffs.to_numpy(), hdi_prob=0.94)\n",
"\n",
"plt.figure(figsize=(7, 3))\n",
"plt.scatter(0, mean_diff)\n",
"plt.axvline(0, hdi[0], hdi[1])\n",
"hdi"
]
},
{
"cell_type": "code",
"execution_count": 27,
"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>chain</th>\n",
" <th>draw</th>\n",
" <th>mpg_obs</th>\n",
" <th>mpg_mean</th>\n",
" <th>am</th>\n",
" <th>drat</th>\n",
" <th>hp</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>19.594787</td>\n",
" <td>0</td>\n",
" <td>2.760</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>19.548561</td>\n",
" <td>0</td>\n",
" <td>2.760</td>\n",
" <td>147</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>2</td>\n",
" <td>19.231404</td>\n",
" <td>0</td>\n",
" <td>3.080</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>3</td>\n",
" <td>19.175185</td>\n",
" <td>0</td>\n",
" <td>3.080</td>\n",
" <td>147</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>4</td>\n",
" <td>18.533028</td>\n",
" <td>0</td>\n",
" <td>3.695</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79995</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>15</td>\n",
" <td>23.168502</td>\n",
" <td>1</td>\n",
" <td>3.695</td>\n",
" <td>147</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79996</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>16</td>\n",
" <td>23.633115</td>\n",
" <td>1</td>\n",
" <td>3.920</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79997</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>17</td>\n",
" <td>23.575700</td>\n",
" <td>1</td>\n",
" <td>3.920</td>\n",
" <td>147</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79998</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>18</td>\n",
" <td>25.464769</td>\n",
" <td>1</td>\n",
" <td>4.930</td>\n",
" <td>146</td>\n",
" </tr>\n",
" <tr>\n",
" <th>79999</th>\n",
" <td>3</td>\n",
" <td>999</td>\n",
" <td>19</td>\n",
" <td>25.403568</td>\n",
" <td>1</td>\n",
" <td>4.930</td>\n",
" <td>147</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>80000 rows × 7 columns</p>\n",
"</div>"
],
"text/plain": [
" chain draw mpg_obs mpg_mean am drat hp\n",
"0 0 0 0 19.594787 0 2.760 146\n",
"1 0 0 1 19.548561 0 2.760 147\n",
"2 0 0 2 19.231404 0 3.080 146\n",
"3 0 0 3 19.175185 0 3.080 147\n",
"4 0 0 4 18.533028 0 3.695 146\n",
"... ... ... ... ... .. ... ...\n",
"79995 3 999 15 23.168502 1 3.695 147\n",
"79996 3 999 16 23.633115 1 3.920 146\n",
"79997 3 999 17 23.575700 1 3.920 147\n",
"79998 3 999 18 25.464769 1 4.930 146\n",
"79999 3 999 19 25.403568 1 4.930 147\n",
"\n",
"[80000 rows x 7 columns]"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_with_obs"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"posterior_diff = pd.DataFrame((posterior_with_obs\n",
" .groupby([\"am\", \"drat\", \"chain\", \"draw\"])[\"mpg_mean\"]\n",
" .diff()\n",
" #.dropna()\n",
" .reset_index(drop=True)\n",
" )).rename(columns={\"mpg_mean\": \"mpg_diff\"})\n",
"\n",
"posterior_diff_obs = pd.concat([posterior_diff, posterior_with_obs], axis=1).dropna(axis=0, how=\"any\")\n",
"\n",
"posterior_diff_mean = pd.DataFrame((posterior_diff_obs[[\"mpg_diff\", \"chain\", \"draw\", \"mpg_obs\"]]\n",
" .groupby([\"mpg_obs\"])[\"mpg_diff\"]\n",
" .mean()\n",
" .reset_index(drop=True)\n",
"))\n",
"\n",
"posterior_diff_std = pd.DataFrame((posterior_diff_obs[[\"mpg_diff\", \"chain\", \"draw\", \"mpg_obs\"]]\n",
" .groupby([\"mpg_obs\"])[\"mpg_diff\"]\n",
" .std()\n",
" .reset_index(drop=True)\n",
"))"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"# posterior_diff = pd.DataFrame((posterior_with_obs\n",
"# .groupby([\"am\", \"drat\"])[\"mpg_mean\"]\n",
"# .diff()\n",
"# #.dropna()\n",
"# .reset_index(drop=True)\n",
"# )).rename(columns={\"mpg_mean\": \"mpg_diff\"})\n",
"\n",
"# posterior_diff_obs = pd.concat([posterior_diff, posterior_with_obs], axis=1).dropna(axis=0, how=\"any\")\n",
"\n",
"# posterior_diff_mean = pd.DataFrame((posterior_diff_obs[[\"mpg_diff\", \"chain\", \"draw\", \"mpg_obs\"]]\n",
"# .groupby([\"mpg_obs\"])[\"mpg_diff\"]\n",
"# .mean()\n",
"# .reset_index(drop=True)\n",
"# ))\n",
"\n",
"# posterior_diff_std = pd.DataFrame((posterior_diff_obs[[\"mpg_diff\", \"chain\", \"draw\", \"mpg_obs\"]]\n",
"# .groupby([\"mpg_obs\"])[\"mpg_diff\"]\n",
"# .std()\n",
"# .reset_index(drop=True)\n",
"# ))"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [],
"source": [
"contrast_df[\"estimate_full\"] = posterior_diff_mean[\"mpg_diff\"]\n",
"contrast_df[\"std_full_lower\"] = posterior_diff_mean[\"mpg_diff\"] - posterior_diff_std[\"mpg_diff\"]\n",
"contrast_df[\"std_full_higher\"] = posterior_diff_mean[\"mpg_diff\"] + posterior_diff_std[\"mpg_diff\"]"
]
},
{
"cell_type": "code",
"execution_count": 48,
"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>estimate</th>\n",
" <th>hdi_3%</th>\n",
" <th>hdi_97%</th>\n",
" <th>term</th>\n",
" <th>contrast</th>\n",
" <th>am</th>\n",
" <th>drat</th>\n",
" <th>estimate_full</th>\n",
" <th>std_full_lower</th>\n",
" <th>std_full_higher</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>-0.046206</td>\n",
" <td>-0.003756</td>\n",
" <td>-0.026338</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>2.760</td>\n",
" <td>-0.046206</td>\n",
" <td>-0.066607</td>\n",
" <td>-0.025805</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.049976</td>\n",
" <td>-0.037245</td>\n",
" <td>-0.056671</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.080</td>\n",
" <td>-0.049976</td>\n",
" <td>-0.065686</td>\n",
" <td>-0.034266</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>-0.057222</td>\n",
" <td>-0.047034</td>\n",
" <td>-0.031519</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.695</td>\n",
" <td>-0.057222</td>\n",
" <td>-0.073846</td>\n",
" <td>-0.040598</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.059873</td>\n",
" <td>-0.119706</td>\n",
" <td>-0.104732</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.920</td>\n",
" <td>-0.059873</td>\n",
" <td>-0.079970</td>\n",
" <td>-0.039775</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>-0.071772</td>\n",
" <td>-0.142725</td>\n",
" <td>-0.086004</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>4.930</td>\n",
" <td>-0.071772</td>\n",
" <td>-0.113567</td>\n",
" <td>-0.029977</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>-0.028507</td>\n",
" <td>-0.001350</td>\n",
" <td>0.017412</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>2.760</td>\n",
" <td>-0.028507</td>\n",
" <td>-0.057614</td>\n",
" <td>0.000601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>-0.036580</td>\n",
" <td>-0.014789</td>\n",
" <td>-0.020630</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.080</td>\n",
" <td>-0.036580</td>\n",
" <td>-0.058501</td>\n",
" <td>-0.014660</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>-0.052097</td>\n",
" <td>-0.038705</td>\n",
" <td>-0.045089</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.695</td>\n",
" <td>-0.052097</td>\n",
" <td>-0.064127</td>\n",
" <td>-0.040067</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>-0.057774</td>\n",
" <td>-0.207495</td>\n",
" <td>-0.198384</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.920</td>\n",
" <td>-0.057774</td>\n",
" <td>-0.069740</td>\n",
" <td>-0.045809</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>-0.083257</td>\n",
" <td>-0.148901</td>\n",
" <td>-0.072400</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>4.930</td>\n",
" <td>-0.083257</td>\n",
" <td>-0.113862</td>\n",
" <td>-0.052652</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" estimate hdi_3% hdi_97% term contrast am drat estimate_full \n",
"0 -0.046206 -0.003756 -0.026338 hp [146, 147] 0 2.760 -0.046206 \\\n",
"1 -0.049976 -0.037245 -0.056671 hp [146, 147] 0 3.080 -0.049976 \n",
"2 -0.057222 -0.047034 -0.031519 hp [146, 147] 0 3.695 -0.057222 \n",
"3 -0.059873 -0.119706 -0.104732 hp [146, 147] 0 3.920 -0.059873 \n",
"4 -0.071772 -0.142725 -0.086004 hp [146, 147] 0 4.930 -0.071772 \n",
"5 -0.028507 -0.001350 0.017412 hp [146, 147] 1 2.760 -0.028507 \n",
"6 -0.036580 -0.014789 -0.020630 hp [146, 147] 1 3.080 -0.036580 \n",
"7 -0.052097 -0.038705 -0.045089 hp [146, 147] 1 3.695 -0.052097 \n",
"8 -0.057774 -0.207495 -0.198384 hp [146, 147] 1 3.920 -0.057774 \n",
"9 -0.083257 -0.148901 -0.072400 hp [146, 147] 1 4.930 -0.083257 \n",
"\n",
" std_full_lower std_full_higher \n",
"0 -0.066607 -0.025805 \n",
"1 -0.065686 -0.034266 \n",
"2 -0.073846 -0.040598 \n",
"3 -0.079970 -0.039775 \n",
"4 -0.113567 -0.029977 \n",
"5 -0.057614 0.000601 \n",
"6 -0.058501 -0.014660 \n",
"7 -0.064127 -0.040067 \n",
"8 -0.069740 -0.045809 \n",
"9 -0.113862 -0.052652 "
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"contrast_df"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [],
"source": [
"drat_unique = contrast_df[\"drat\"].unique()\n",
"main_levels = len(drat_unique)\n",
"idx_main = np.arange(main_levels)"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [],
"source": [
"vals = contrast_df[[\"std_full_lower\", \"std_full_higher\"]].values"
]
},
{
"cell_type": "code",
"execution_count": 76,
"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>estimate</th>\n",
" <th>hdi_3%</th>\n",
" <th>hdi_97%</th>\n",
" <th>term</th>\n",
" <th>contrast</th>\n",
" <th>am</th>\n",
" <th>drat</th>\n",
" <th>estimate_full</th>\n",
" <th>std_full_lower</th>\n",
" <th>std_full_higher</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>-0.046206</td>\n",
" <td>-0.003756</td>\n",
" <td>-0.026338</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>2.760</td>\n",
" <td>-0.046206</td>\n",
" <td>-0.066607</td>\n",
" <td>-0.025805</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>-0.049976</td>\n",
" <td>-0.037245</td>\n",
" <td>-0.056671</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.080</td>\n",
" <td>-0.049976</td>\n",
" <td>-0.065686</td>\n",
" <td>-0.034266</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>-0.057222</td>\n",
" <td>-0.047034</td>\n",
" <td>-0.031519</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.695</td>\n",
" <td>-0.057222</td>\n",
" <td>-0.073846</td>\n",
" <td>-0.040598</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-0.059873</td>\n",
" <td>-0.119706</td>\n",
" <td>-0.104732</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>3.920</td>\n",
" <td>-0.059873</td>\n",
" <td>-0.079970</td>\n",
" <td>-0.039775</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>-0.071772</td>\n",
" <td>-0.142725</td>\n",
" <td>-0.086004</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>0</td>\n",
" <td>4.930</td>\n",
" <td>-0.071772</td>\n",
" <td>-0.113567</td>\n",
" <td>-0.029977</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>-0.028507</td>\n",
" <td>-0.001350</td>\n",
" <td>0.017412</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>2.760</td>\n",
" <td>-0.028507</td>\n",
" <td>-0.057614</td>\n",
" <td>0.000601</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>-0.036580</td>\n",
" <td>-0.014789</td>\n",
" <td>-0.020630</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.080</td>\n",
" <td>-0.036580</td>\n",
" <td>-0.058501</td>\n",
" <td>-0.014660</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>-0.052097</td>\n",
" <td>-0.038705</td>\n",
" <td>-0.045089</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.695</td>\n",
" <td>-0.052097</td>\n",
" <td>-0.064127</td>\n",
" <td>-0.040067</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>-0.057774</td>\n",
" <td>-0.207495</td>\n",
" <td>-0.198384</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>3.920</td>\n",
" <td>-0.057774</td>\n",
" <td>-0.069740</td>\n",
" <td>-0.045809</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>-0.083257</td>\n",
" <td>-0.148901</td>\n",
" <td>-0.072400</td>\n",
" <td>hp</td>\n",
" <td>[146, 147]</td>\n",
" <td>1</td>\n",
" <td>4.930</td>\n",
" <td>-0.083257</td>\n",
" <td>-0.113862</td>\n",
" <td>-0.052652</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" estimate hdi_3% hdi_97% term contrast am drat estimate_full \n",
"0 -0.046206 -0.003756 -0.026338 hp [146, 147] 0 2.760 -0.046206 \\\n",
"1 -0.049976 -0.037245 -0.056671 hp [146, 147] 0 3.080 -0.049976 \n",
"2 -0.057222 -0.047034 -0.031519 hp [146, 147] 0 3.695 -0.057222 \n",
"3 -0.059873 -0.119706 -0.104732 hp [146, 147] 0 3.920 -0.059873 \n",
"4 -0.071772 -0.142725 -0.086004 hp [146, 147] 0 4.930 -0.071772 \n",
"5 -0.028507 -0.001350 0.017412 hp [146, 147] 1 2.760 -0.028507 \n",
"6 -0.036580 -0.014789 -0.020630 hp [146, 147] 1 3.080 -0.036580 \n",
"7 -0.052097 -0.038705 -0.045089 hp [146, 147] 1 3.695 -0.052097 \n",
"8 -0.057774 -0.207495 -0.198384 hp [146, 147] 1 3.920 -0.057774 \n",
"9 -0.083257 -0.148901 -0.072400 hp [146, 147] 1 4.930 -0.083257 \n",
"\n",
" std_full_lower std_full_higher \n",
"0 -0.066607 -0.025805 \n",
"1 -0.065686 -0.034266 \n",
"2 -0.073846 -0.040598 \n",
"3 -0.079970 -0.039775 \n",
"4 -0.113567 -0.029977 \n",
"5 -0.057614 0.000601 \n",
"6 -0.058501 -0.014660 \n",
"7 -0.064127 -0.040067 \n",
"8 -0.069740 -0.045809 \n",
"9 -0.113862 -0.052652 "
]
},
"execution_count": 76,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# interval contains mean estimate\n",
"contrast_df"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import seaborn.objects as so\n",
"\n",
"(\n",
" so.Plot(contrast_df, x=\"am\", y=\"estimate_full\", color=\"drat\")\n",
" .add(so.Range(), so.Dodge())\n",
" .add(so.Dot(), so.Agg(), so.Dodge())\n",
")\n",
"\n",
"#plt.axvline(0, vals[:, 0], vals[:, 1])\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Why are the vlines showing like this? Looking at the df above, the estimate_full is always within the lower and higher interval....."
]
},
{
"cell_type": "code",
"execution_count": 114,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABFcAAAHACAYAAABnBYehAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABW6ElEQVR4nO3deXhU5f3//9dkTyCZREI2EkgAyw4igZgIAlY2ZdOqIH5SRaWiUATajx8QW1JbAVERWxAqgmhdghVRWjGFWkERQthSIdAoGAjbsJlMgGwkOb8/+DK/jllIMpNMhjwf15Xr8tznfc+8TzSReXGf+5gMwzAEAAAAAACAevFwdQMAAAAAAADujHAFAAAAAADAAYQrAAAAAAAADiBcAQAAAAAAcADhCgAAAAAAgAMIVwAAAAAAABxAuAIAAAAAAOAAwhUAAAAAAAAHeLm6getBRUWFTp48qcDAQJlMJle3AwAAAMBFDMPQhQsXFBUVJQ+Ppvt32eXl5bp8+bKr2wCaLG9vb3l6eta6nnDFCU6ePKmYmBhXtwEAAACgiTh27Jiio6Nd3UYlhmHIYrEoPz/f1a0ATV5wcLAiIiJqtYiCcMUJAgMDJV35BRoUFOTibgAAAAC4SkFBgWJiYmyfEZqaq8FKWFiYAgICWHkPVMEwDBUWFurMmTOSpMjIyGvOIVxxgqu/kIKCgghXAAAAADTJ0KK8vNwWrLRq1crV7QBNmr+/vyTpzJkzCgsLu+YtQk33JkAAAAAAgNNc3WMlICDAxZ0A7uHqz0pt9iciXAEAAACAZqQprqoBmqK6/KwQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAgFoxDEPpR/M0e8NBTf7wG83ecFDpR/NkGEaDvef8+fPVt29fBQYGKiwsTGPHjlV2dnaNcx5++GGZTKZKX926dbOry8/P15QpUxQZGSk/Pz916dJFGzZsaLBrwfWLpwUBAAAAAK4py3JBE1P3atdxq934C/86pPhos94c31vdIpz/COotW7ZoypQp6tu3r8rKyjRnzhwNHTpUBw4cUIsWLaqc8+qrr2rBggW247KyMvXq1Uv33Xefbay0tFRDhgxRWFiYPvzwQ0VHR+vYsWNN9jHaaNoIVwAAAAAANcqyXNCApV8rv6jqp6bsOm7VgKVf66sptzo9YElLS7M7fvPNNxUWFqbdu3frtttuq3KO2WyW2Wy2HX/88cfKy8vTxIkTbWOrVq3SDz/8oG3btsnb21uS1K5dO6f2juaD24IAAAAAANUyDEMTU/dWG6xclV90WY+syWzQW4QkyWq9snLmhhtuqPWclStX6o477rALT9avX6/ExERNmTJF4eHh6t69u+bNm6fy8nKn94zrH+EKAAAAAKBaO3LzK90KVJ2dx/KVkZvfYL0YhqGZM2eqf//+6t69e63mnDp1Sp999pkee+wxu/Hvv/9eH374ocrLy7VhwwY9++yzevnll/X88883ROu4znFbEAAAAACgWp9kWepU/3GWRQntQhqkl6lTp+qbb77R1q1baz1n9erVCg4O1tixY+3GKyoqFBYWptdff12enp7q06ePTp48qRdffFG//e1vndw5rneEKwAAAACAauUV1nw7kKP1tfXLX/5S69ev15dffqno6OhazTEMQ6tWrVJycrJ8fHzszkVGRsrb21uenp62sS5dushisai0tLRSPVATbgsCAAAAAFQrJMC7QeuvxTAMTZ06VR999JH+9a9/KS4urtZzt2zZokOHDunRRx+tdO7WW2/VoUOHVFFRYRv79ttvFRkZSbCCOnO7cOW1115TXFyc/Pz81KdPH3311Vc11m/ZskV9+vSRn5+f2rdvr+XLl1eqWbt2rbp27SpfX1917dpV69ata6j2AQAAAMCtjOkWUaf6sXWsv5YpU6bonXfe0XvvvafAwEBZLBZZLBYVFRXZambPnq2f//znleauXLlSCQkJVe7P8sQTT+j8+fN66qmn9O233+rTTz/VvHnzNGXKFKf2j+bBrcKVNWvWaPr06ZozZ4727t2rAQMGaMSIEcrNza2yPicnR3feeacGDBigvXv36plnntG0adO0du1aW8327ds1btw4JScn69///reSk5N1//33a8eOHY11WQAAAADQZCW0DVZ8tPnahZL6xgSrX9tgp77/smXLZLVaNWjQIEVGRtq+1qxZY6s5depUpc+FVqtVa9eurXLViiTFxMRo48aN2rlzp3r27Klp06bpqaee0qxZs5zaP5oHk9HQz8lyooSEBN18881atmyZbaxLly4aO3as5s+fX6n+//7v/7R+/XodPHjQNjZ58mT9+9//1vbt2yVJ48aNU0FBgT777DNbzfDhwxUSEqL333+/Vn0VFBTIbDbLarUqKCiovpcHAAAANEmlZRVa+sUhSdKUwR3l4+VWf0fbqJryZ4Pi4mLl5OTY7gSoiyzLBQ1Y+nWNj2MO9vfWV1NuVbeIQEdbBZqEuvzMuM1vxdLSUu3evVtDhw61Gx86dKi2bdtW5Zzt27dXqh82bJh27dqly5cv11hT3WtKUklJiQoKCuy+AAAAAOB61S0iUF9NubXaFSx9Y4IJVtCsuc3Tgs6dO6fy8nKFh4fbjYeHh8tiqfrRYBaLpcr6srIynTt3TpGRkdXWVPeakjR//nz97ne/q+eVAAAAAID76RYRqB1PDVBGbr4+zrIor/CyQgK8NbZbhPq1DZbJZHJ1i4DLuE24ctWPf2ANw6jxh7iq+h+P1/U1Z8+erZkzZ9qOCwoKFBMTc+3mAQAAAMCNmUwmJbQLUUK7EFe3AjQpbhOuhIaGytPTs9KKkjNnzlRaeXJVRERElfVeXl5q1apVjTXVvaYk+fr6ytfXtz6XAQAAAAAArjNus+eKj4+P+vTpo02bNtmNb9q0SUlJSVXOSUxMrFS/ceNGxcfHy9vbu8aa6l4TAAAAAADgv7nNyhVJmjlzppKTkxUfH6/ExES9/vrrys3N1eTJkyVduV3nxIkTevvttyVdeTLQkiVLNHPmTE2aNEnbt2/XypUr7Z4C9NRTT+m2227TCy+8oDFjxuiTTz7RP//5T23dutUl1wgAAAAAANyLW4Ur48aN0/nz5/Xcc8/p1KlT6t69uzZs2KB27dpJqvxs87i4OG3YsEEzZszQ0qVLFRUVpT/+8Y/62c9+ZqtJSkpSamqqnn32Wf3mN79Rhw4dtGbNGiUkJDT69QEAAAAAAPdjMq7u8Ip6a8rPsgcAAAAcVVpWoaVfHJIkTRncUT5ebrO7QKNryp8NiouLlZOTo7i4OPn5+bm6HaDJq8vPDL8VAQAAAAAAHEC4AgAAAAAA4ADCFQAAAABArRiGoeJTO/TD1jk69/kU/bB1jopP7VBD7jaxbNky9ezZU0FBQQoKClJiYqI+++yzGuds2bJFffr0kZ+fn9q3b6/ly5dXqlm8eLE6deokf39/xcTEaMaMGSouLm6oy8B1zq02tAUAAAAAuEbp+Syd3fiYSk/vthu37npRPuF91HroG/Jp1c3p7xsdHa0FCxaoY8eOkqS33npLY8aM0d69e9WtW+X3y8nJ0Z133qlJkybpnXfe0ddff60nn3xSrVu3tj3c5N1339WsWbO0atUqJSUl6dtvv9XDDz8sSXrllVecfg24/hGuAAAAAABqVHo+S6c+GKyKkvyqz5/erVMfDFbk/V84PWAZNWqU3fHzzz+vZcuWKT09vcpwZfny5Wrbtq0WL14sSerSpYt27dqll156yRaubN++XbfeeqsmTJggSYqNjdUDDzygjIwMp/aO5oPbggAAAAAA1TIMQ2c3PlZtsHJVRUm+zm6c1KC3CJWXlys1NVWXLl1SYmJilTXbt2/X0KFD7caGDRumXbt26fLly5Kk/v37a/fu3bYw5fvvv9eGDRt01113NVjvuL6xcgUAAAAAUK0SS0alW4GqU3p6l0osO+UX2c+pPezbt0+JiYkqLi5Wy5YttW7dOnXt2rXKWovFovDwcLux8PBwlZWV6dy5c4qMjNT48eN19uxZ9e/fX4ZhqKysTE888YRmzZrl1L7RfLByBQAAAABQrcLD6+tY/4nTe+jUqZMyMzOVnp6uJ554Qg899JAOHDhQbb3JZLI7vrqa5ur45s2b9fzzz+u1117Tnj179NFHH+nvf/+7fv/73zu9dzQPrFwBAAAAAFTrWrcDOVpfGz4+PrYNbePj47Vz5069+uqr+vOf/1ypNiIiQhaLxW7szJkz8vLyUqtWrSRJv/nNb5ScnKzHHntMktSjRw9dunRJv/jFLzRnzhx5eLAOAXXDfzEAAAAAgGp5+AY3aH19GIahkpKSKs8lJiZq06ZNdmMbN25UfHy8vL29JUmFhYWVAhRPT08ZhtGge8bg+kW4AgAAAACoVkCH0XWsH+PU93/mmWf01Vdf6ciRI9q3b5/mzJmjzZs368EHH5QkzZ49Wz//+c9t9ZMnT9bRo0c1c+ZMHTx4UKtWrdLKlSv161//2lYzatQoLVu2TKmpqcrJydGmTZv0m9/8RqNHj5anp6dT+0fzwG1BAAAAAIBq+Ub0k094n1ptausTHi/fiL5Off/Tp08rOTlZp06dktlsVs+ePZWWlqYhQ4ZIkk6dOqXc3FxbfVxcnDZs2KAZM2Zo6dKlioqK0h//+EfbY5gl6dlnn5XJZNKzzz6rEydOqHXr1ho1apSef/55p/aO5sNksObJYQUFBTKbzbJarQoKCnJ1OwAAAIBTlZZVaOkXhyRJUwZ3lI8XC+Cr05Q/GxQXFysnJ0dxcXHy8/Or09zS81k69cHgGvdT8fANVuT9X8inVTcHOwWahrr8zPBbEQAAAABQI59W3a4EJ+F9qj4fHk+wgmaN24IAAAAAANfk06qbosZvU4llpwoPf6KKknx5+AYroMMY+Ub0rfT4Y6A5IVwBAAAAANSKyWSSX2Q/+UX2c3UrQJPCbUEAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAmqxly5apZ8+eCgoKUlBQkBITE/XZZ5/VOKekpERz5sxRu3bt5Ovrqw4dOmjVqlW285cvX9Zzzz2nDh06yM/PT7169VJaWprda6SkpMhkMtl9RURE2NUYhqGUlBRFRUXJ399fgwYNUlZWVqVefvnLXyo0NFQtWrTQ6NGjdfz4cbuavLw8JScny2w2y2w2Kzk5Wfn5+XY1ubm5GjVqlFq0aKHQ0FBNmzZNpaWldjX79u3TwIED5e/vrzZt2ui5556TYRh2NVu2bFGfPn3k5+en9u3ba/ny5ZW+f2vXrlXXrl3l6+urrl27at26dZVqXnvtNcXFxcnPz099+vTRV1991WS/N9nZ2Ro8eLDCw8Nt1/3ss8/q8uXLla6rvghXAAAAAAC1YhiG0s8c1TO7NuiJbR/qmV0blH7maKUP8M4UHR2tBQsWaNeuXdq1a5duv/12jRkzptIH9f92//336/PPP9fKlSuVnZ2t999/X507d7adf/bZZ/XnP/9Zf/rTn3TgwAFNnjxZd999t/bu3Wv3Ot26ddOpU6dsX/v27bM7v3DhQi1atEhLlizRzp07FRERoSFDhujChQu2munTp2vdunVKTU3V1q1bdfHiRY0cOVLl5eW2mgkTJigzM1NpaWlKS0tTZmamkpOTbefLy8t111136dKlS9q6datSU1O1du1a/epXv7LVFBQUaMiQIYqKitLOnTv1pz/9SS+99JIWLVpkq8nJydGdd96pAQMGaO/evXrmmWc0bdo0rV271lazfft2jRs3TsnJyfr3v/+t5ORk3X///dqxY4etZs2aNZo+fbrmzJmjvXv3asCAARoxYoRyc3Ob5PfG29tbP//5z7Vx40ZlZ2dr8eLFWrFihebOnVvpv516M+Awq9VqSDKsVqurWwEAAACcruRyubFoY7axaGO2UXK53NXtNGlN+bNBUVGRceDAAaOoqKhe8/f/cMrot/4Vw2PVryp99Vv/irH/h1NO7rh6ISEhxhtvvFHluc8++8wwm83G+fPnq50fGRlpLFmyxG5szJgxxoMPPmg7njt3rtGrV69qX6OiosKIiIgwFixYYBsrLi42zGazsXz5csMwDCM/P9/w9vY2UlNTbTUnTpwwPDw8jLS0NMMwDOPAgQOGJCM9Pd1Ws337dkOS8Z///McwDMPYsGGD4eHhYZw4ccJW8/777xu+vr62/9Zee+01w2w2G8XFxbaa+fPnG1FRUUZFRYVhGIbx9NNPG507d7a7jscff9y45ZZbbMf333+/MXz4cLuaYcOGGePHj7cd9+vXz5g8ebJdTefOnY1Zs2Y1ye9NVWbMmGH079+/2vOGUbefGVauAAAAAABqlJVn0W0blmrXueNVnt917rhu27BUWXmWBu2jvLxcqampunTpkhITE6usWb9+veLj47Vw4UK1adNGP/nJT/TrX/9aRUVFtpqSkhL5+fnZzfP399fWrVvtxr777jtFRUUpLi5O48eP1/fff287l5OTI4vFoqFDh9rGfH19NXDgQG3btk2StHv3bl2+fNmuJioqSt27d7fVbN++XWazWQkJCbaaW265RWaz2a6me/fuioqKstUMGzZMJSUl2r17t61m4MCB8vX1tas5efKkjhw5Yqv5716u1uzatct2i0x1NVd7KS0t1e7duyvVDB061FbT1L43P3bo0CGlpaVp4MCBVZ6vD8IVAAAAAEC1DMPQI1tTlV9aVGNdfmmRHt26pkFuEdq3b59atmwpX19fTZ48WevWrVPXrl2rrP3++++1detW7d+/X+vWrdPixYv14YcfasqUKbaaYcOGadGiRfruu+9UUVGhTZs26ZNPPtGpU6dsNQkJCXr77bf1j3/8QytWrJDFYlFSUpLOnz8vSbJYrgRJ4eHhdu8fHh5uO2exWOTj46OQkJAaa8LCwipdR1hYmF3Nj98nJCREPj4+NdZcPb5WTVlZmc6dO1djzdXXOHfunMrLy6953U3pe3NVUlKS/Pz8dOONN2rAgAF67rnnKr12fRGuAAAAAACqteNsbrUrVn5s57ljyjiXe+3COurUqZMyMzOVnp6uJ554Qg899JAOHDhQZW1FRYVMJpPeffdd9evXT3feeacWLVqk1atX21avvPrqq7rxxhvVuXNn+fj4aOrUqZo4caI8PT1trzNixAj97Gc/U48ePXTHHXfo008/lSS99dZbdu9nMpnsjg3DqDT2Yz+uqareGTVXgy5n1Px4zFk1P9ZQ3xvpyl4xe/bs0XvvvadPP/1UL730Uo291AXhCgAAAACgWutzq984tiqfHK1bfW34+PioY8eOio+P1/z589WrVy+9+uqrVdZGRkaqTZs2MpvNtrEuXbrIMAzbk2hat26tjz/+WJcuXdLRo0f1n//8Ry1btlRcXFy1PbRo0UI9evTQd999J0m2Jwf9eHXEmTNnbCspIiIiVFpaqry8vBprTp8+Xen9zp49a1fz4/fJy8vT5cuXa6w5c+aMJF2zxsvLS61ataqx5uprhIaGytPT85rX3ZS+N1fFxMSoa9eueuCBB7RgwQKlpKTYbZ7rCMIVAAAAAEC18koLG7S+PgzDUElJSZXnbr31Vp08eVIXL160jX377bfy8PBQdHS0Xa2fn5/atGmjsrIyrV27VmPGjKn2PUtKSnTw4EFFRkZKkuLi4hQREaFNmzbZakpLS7VlyxYlJSVJkvr06SNvb2+7mlOnTmn//v22msTERFmtVmVkZNhqduzYIavValezf/9+u9uWNm7cKF9fX/Xp08dW8+WXX9o9gnjjxo2KiopSbGysrea/e7laEx8fL29v7xprrvbi4+OjPn36VKrZtGmTraapfW+qYhiGLl++7Lzb2K655S2uqSnvCA4AAAA4iqcF1V5T/mxQ36cFzd75aZVPCKrua/bOT53a9+zZs40vv/zSyMnJMb755hvjmWeeMTw8PIyNGzcahmEYs2bNMpKTk231Fy5cMKKjo417773XyMrKMrZs2WLceOONxmOPPWarSU9PN9auXWscPnzY+PLLL43bb7/diIuLM/Ly8mw1v/rVr4zNmzcb33//vZGenm6MHDnSCAwMNI4cOWKrWbBggWE2m42PPvrI2Ldvn/HAAw8YkZGRRkFBga1m8uTJRnR0tPHPf/7T2LNnj3H77bcbvXr1MsrKymw1w4cPN3r27Gls377d2L59u9GjRw9j5MiRtvNlZWVG9+7djZ/+9KfGnj17jH/+859GdHS0MXXqVFtNfn6+ER4ebjzwwAPGvn37jI8++sgICgoyXnrpJVvN999/bwQEBBgzZswwDhw4YKxcudLw9vY2PvzwQ1vN119/bXh6ehoLFiwwDh48aCxYsMDw8vKye2JPamqq4e3tbaxcudI4cOCAMX36dKNFixZN9nvzzjvvGGvWrDEOHDhgHD582Pjggw+MNm3a2D0dqip1+ZkhXHGCpvwLFAAAAHAU4UrtNeXPBvUNV7afPlKncCX9zJFrv2gdPPLII0a7du0MHx8fo3Xr1sZPf/pTW7BiGIbx0EMPGQMHDrSbc/DgQeOOO+4w/P39jejoaGPmzJlGYWGh7fzmzZuNLl26GL6+vkarVq2M5ORku0f5GoZhjBs3zoiMjDS8vb2NqKgo45577jGysrLsaioqKoy5c+caERERhq+vr3HbbbcZ+/bts6spKioypk6datxwww2Gv7+/MXLkSCM3N9eu5vz588aDDz5oBAYGGoGBgcaDDz5oF/QYhmEcPXrUuOuuuwx/f3/jhhtuMKZOnWr32GXDMIxvvvnGGDBggOHr62tEREQYKSkptscw//e19+7d2/Dx8TFiY2ONZcuWVfqe//WvfzU6depkeHt7G507dzbWrl1bqWbp0qW2fy8333yzsWXLlib7vUlNTTVuvvlmo2XLlkaLFi2Mrl27GvPmzbvmz0JdfmZMhtEAWzk3MwUFBTKbzbJarQoKCnJ1OwAAAIBTlZZVaOkXhyRJUwZ3lI8XuwtUpyl/NiguLlZOTo7i4uIqPYa4JoZh6Ja/v1qrTW37hsZo+8hp19y0FHAHdfmZ4bciAAAAAKBaJpNJq/qPV7CPf411wT7+Wtl/HMEKmiXCFQAAAABAjbqFROjLO6coPjS6yvN9Q2P05Z1T1C0kopE7A5oGL1c3AAAAAABo+rqFRCh95FPKOJerT45mKa+0UCE+ARrTrpv6hbZlxQqaNcIVAAAAAECtmEwmJbRup4TW7VzdCtCkcFsQAAAAAACAA1i5ArfEjvUAAAAAgKaCT6QAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAACoFcMwVHgoXac/mK1Tqyfr9AezVXgoXYZhNNh7Llu2TD179lRQUJCCgoKUmJiozz77rMY5S5cuVZcuXeTv769OnTrp7bfftju/YsUKDRgwQCEhIQoJCdEdd9yhjIyMBrsGXP94WhAAAAAA4JqKj2fp5BsPqzhnl934+U8XyC8uXlGPrZZfdDenv290dLQWLFigjh07SpLeeustjRkzRnv37lW3bpXfb9myZZo9e7ZWrFihvn37KiMjQ5MmTVJISIhGjRolSdq8ebMeeOABJSUlyc/PTwsXLtTQoUOVlZWlNm3aOP0acP0zGQ0ZMTYTBQUFMpvNslqtCgoKcnU7zQKPYgYAAGg8/Nmr9pryZ4Pi4mLl5OQoLi5Ofn5+dZt7PEtHnu+visL8ams8AoIVO2drgwQsP3bDDTfoxRdf1KOPPlrpXFJSkm699Va9+OKLtrHp06dr165d2rp1a5WvV15erpCQEC1ZskQ///nPG6xvuJe6/MzwWxEAAAAAUC3DMHTyjYdrDFYkqaIwXyffmNigtwiVl5crNTVVly5dUmJiYpU1JSUllT4I+/v7KyMjQ5cvX65yTmFhoS5fvqwbbrjB6T2jeSBcAQAAAABUq+jwjkq3AlWnOGenir53/t4l+/btU8uWLeXr66vJkydr3bp16tq1a5W1w4YN0xtvvKHdu3fLMAzt2rVLq1at0uXLl3Xu3Lkq58yaNUtt2rTRHXfc4fTe0TwQrgAAAAAAqnVhzyd1q9/9sdN76NSpkzIzM5Wenq4nnnhCDz30kA4cOFBl7W9+8xuNGDFCt9xyi7y9vTVmzBg9/PDDkiRPT89K9QsXLtT777+vjz76qM63SwFXEa4AAAAAAKpVUZjXoPW14ePjo44dOyo+Pl7z589Xr1699Oqrr1ZZ6+/vr1WrVqmwsFBHjhxRbm6uYmNjFRgYqNDQULval156SfPmzdPGjRvVs2dPp/eN5sNtwpW8vDwlJyfLbDbLbDYrOTlZ+fn5Nc4xDEMpKSmKioqSv7+/Bg0apKysLNv5H374Qb/85S/VqVMnBQQEqG3btpo2bZqsVmsDXw0AAAAAuAePgJAGra8PwzBUUlJSY423t7eio6Pl6emp1NRUjRw5Uh4e//9H4BdffFG///3vlZaWpvj4+IZuGdc5twlXJkyYoMzMTKWlpSktLU2ZmZlKTk6ucc7ChQu1aNEiLVmyRDt37lRERISGDBmiCxcuSJJOnjypkydP6qWXXtK+ffu0evVqpaWlVbnjNAAAAAA0R4E3j6lbfZ+xTn3/Z555Rl999ZWOHDmiffv2ac6cOdq8ebMefPBBSdLs2bPtnvDz7bff6p133tF3332njIwMjR8/Xvv379e8efNsNQsXLtSzzz6rVatWKTY2VhaLRRaLRRcvXnRq72g+vFzdQG0cPHhQaWlpSk9PV0JCgiRpxYoVSkxMVHZ2tjp16lRpjmEYWrx4sebMmaN77rlH0pXnoYeHh+u9997T448/ru7du2vt2rW2OR06dNDzzz+v//mf/1FZWZm8vNzi2wMAAAAADca/Q4L84uJrtamtX1xf+bfv59T3P336tJKTk3Xq1CmZzWb17NlTaWlpGjJkiCTp1KlTys3NtdWXl5fr5ZdfVnZ2try9vTV48GBt27ZNsbGxtprXXntNpaWluvfee+3ea+7cuUpJSXFq/2ge3CI92L59u8xmsy1YkaRbbrlFZrNZ27ZtqzJcycnJkcVi0dChQ21jvr6+GjhwoLZt26bHH3+8yve6+jz6moKVkpISuyVoBQUF9bksAAAAAGjyTCaToh5brSPP96/xccweAcGKeuxNmUwmp77/ypUrazy/evVqu+MuXbpo7969Nc45cuSIg10B9tzitiCLxaKwsLBK42FhYbJYLNXOkaTw8HC78fDw8GrnnD9/Xr///e+rDV6umj9/vm3vF7PZrJiYmNpcBgAAAAC4Jb/oboqds1V+cVXvTeIX1/fK+ehujdwZ0DS4dOVKSkqKfve739VYs3PnTkmqMv00DOOaqeiPz1c3p6CgQHfddZe6du2quXPn1vias2fP1syZM+3mErAAAAAAuJ75RXdT3NwMFX2foQu7P1ZFYZ48AkIU2Ges/Nv3c/qKFcCduDRcmTp1qsaPH19jTWxsrL755hudPn260rmzZ89WWplyVUREhKQrK1giIyNt42fOnKk058KFCxo+fLhatmypdevWydvbu8aefH195evrW2MNAAAAAFxvTCaTAjokKKBDwrWLgWbEpeFKaGhopeeMVyUxMVFWq1UZGRnq1+/K5kg7duyQ1WpVUlJSlXPi4uIUERGhTZs2qXfv3pKk0tJSbdmyRS+88IKtrqCgQMOGDZOvr6/Wr18vPz8/J1wZAAAAAABoLtxiz5UuXbpo+PDhmjRpktLT05Wenq5JkyZp5MiRdpvZdu7cWevWrZN0JVGdPn265s2bp3Xr1mn//v16+OGHFRAQoAkTJki6smJl6NChunTpklauXKmCggLbI7jKy8tdcq0AAAAAAMC9uMXTgiTp3Xff1bRp02xP/xk9erSWLFliV5OdnS2r1Wo7fvrpp1VUVKQnn3xSeXl5SkhI0MaNGxUYGChJ2r17t3bs2CFJ6tixo91r5eTk2D2qCwAAAAAAoCpuE67ccMMNeuedd2qsMQzD7thkMiklJaXa55QPGjSo0hwAAAAAAIC6cIvbggAAAAAAAJoqwhUAAAAAAAAHEK4AAAAAAAA4gHAFAAAAAFArhmHoYq5VJ9IO6ei6/+hE2iFdzLU26l6W8+fPtz0dtiZLly5Vly5d5O/vr06dOuntt9+2O//RRx8pPj5ewcHBatGihW666Sb95S9/acDOcT1zmw1tAQAAAACuU3T6oo789YAKT1ywG7dsOaqANoGKva+r/MNbNmgPO3fu1Ouvv66ePXvWWLds2TLNnj1bK1asUN++fZWRkaFJkyYpJCREo0aNknTloSlz5sxR586d5ePjo7///e+aOHGiwsLCNGzYsAa9Dlx/WLkCAAAAAKhR0emLyl6+u1KwclXhiQvKXr5bRacvNlgPFy9e1IMPPqgVK1YoJCSkxtq//OUvevzxxzVu3Di1b99e48eP16OPPqoXXnjBVjNo0CDdfffd6tKlizp06KCnnnpKPXv21NatWxvsGnD9IlwBAAAAAFTLMAwd+esBlReX1VhXXlymIx8ebLBbhKZMmaK77rpLd9xxxzVrS0pK5OfnZzfm7++vjIwMXb58uVK9YRj6/PPPlZ2drdtuu81pPaP54LYgAAAAAEC1Lh0rqHbFyo8VHi9Q4bECtWhrdmoPqamp2rNnj3bu3Fmr+mHDhumNN97Q2LFjdfPNN2v37t1atWqVLl++rHPnzikyMlKSZLVa1aZNG5WUlMjT01OvvfaahgwZ4tTe0TwQrgAAAAAAqmU9cLZO9fkHzjo1XDl27Jieeuopbdy4sdJqlOr85je/kcVi0S233CLDMBQeHq6HH35YCxculKenp60uMDBQmZmZunjxoj7//HPNnDlT7du316BBg5zWP5oHbgsCAAAAAFSrrKjm24Ecrb+W3bt368yZM+rTp4+8vLzk5eWlLVu26I9//KO8vLxUXl5eaY6/v79WrVqlwsJCHTlyRLm5uYqNjVVgYKBCQ0NtdR4eHurYsaNuuukm/epXv9K9996r+fPnO7V/NA+sXAEAAAAAVMvLv24fG+tafy0//elPtW/fPruxiRMnqnPnzvq///s/u5UoP+bt7a3o6GhJV24tGjlypDw8ql9jYBiGSkpKnNM4mhXCletIaVmFln5xSJI0ZXBH+XixMAkAAACAY8xdW8uy5Wit64O7tnbq+wcGBqp79+52Yy1atFCrVq1s47Nnz9aJEyf09ttvS5K+/fZbZWRkKCEhQXl5eVq0aJH279+vt956y/Ya8+fPV3x8vDp06KDS0lJt2LBBb7/9tpYtW+bU/tE8EK4AAAAAAKrVIiZIAW0Ca7WpbUB0kAJighqhK3unTp1Sbm6u7bi8vFwvv/yysrOz5e3trcGDB2vbtm2KjY211Vy6dElPPvmkjh8/Ln9/f3Xu3FnvvPOOxo0b1+j9w/0RrgAAAAAuZhiGduTm65Msi/IKLyskwFtjukUooW2wTCaTq9tDM2cymRR7X1dlL99d4+OYPf28FHtvl0b5b3bz5s12x6tXr7Y77tKli/bu3Vvja/zhD3/QH/7wByd3huaKcAUAAABwoSzLBU1M3atdx6124y/865Dio816c3xvdYsIdFF3wBX+4S3VaXIfHfnrgSpXsAREByn23i7yD2/pgu4A1yNcAQAAAFwky3JBA5Z+rfyiy1We33XcqgFLv9ZXU24lYIHL+Ye3VOcpfVV4rED5B86qrKhMXv5eCu7aWgExQayyQrNGuAIAAAC4gGEYmpi6t9pg5ar8ost6ZE2m0qf158MrXM5kMqlFW7NatDW7uhWgSeFxMgAAAIAL7MjNr3QrUHV2HstXRm5+wzYEAKg3whUAAADABT7JstSp/uM61gMAGg/hCgAAAOACeYU13w7kaD1QHcMwXN0C4Bbq8rNCuAIAAAC4QEiAd4PWAz/m7X3lv6HCwkIXdwK4h6s/K1d/dmrChrYAAACAC4zpFqEX/nWo1vVju0U0YDdoDjw9PRUcHKwzZ85IkgICAtgkGaiCYRgqLCzUmTNnFBwcLE9Pz2vOIVwBAAAAXCChbbDio8212tS2b0yw+rUNbvimcN2LiLgS0l0NWABULzg42PYzcy2EKwAAAIALmEwmvTm+twYs/brGxzEH+3tr1bibWGEApzCZTIqMjFRYWJguX2YfH6A63t7etVqxchXhChpFaVmFln5xZdnrlMEd5ePFdj8AAADdIgL11ZRbNTF1b5UrWPrGBGvVuJvULSLQBd3heubp6VmnD44Aaka4AgAAALhQt4hA7XhqgDJy8/VxlkV5hZcVEuCtsd0i1K9tMCtWAMANEK4AAAAALmYymZTQLkQJ7UJc3QoAoB64NwMAAAAAAMABhCsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAAAAO8HJ1AwAAAACk0rIKLf3ikCRpyuCO8vHi70EBwF0QrgAAAADXEcMwVGLJUOHh9aooyZeHb7ACOoyWb0Q/mUwmV7cHANclwhUAAADgOlF6PktnNz6m0tO77catu16UT3gftR76hnxadXNRdwBw/WKtIQAAAHAdKD2fpVMfDK4UrNjOn9595fz5rEbuDACuf4QrAAAAgJszDENnNz6mipL8GusqSvJ1duMkGYbROI0BQDNBuAIAAAC4uRJLRrUrVn6s9PQulVh2NnBHANC8EK4AAAAAbq7w8Po61n/SQJ0AQPNEuAIAAAC4uWvdDuRoPQCgZoQrAAAAgJvz8A1u0HoAQM0IVwAAAAA3F9BhdB3rxzRQJwDQPBGuAAAAAG7ON6KffML71KrWJzxevhF9G7gjAGheCFcAAAAAN2cymdR66BvXvN3HwzdYrYeukMlkapzGAKCZIFwBAAAArgM+rbop8v4vql3B4hMef+V8q26N3BkAXP/cJlzJy8tTcnKyzGazzGazkpOTlZ+fX+McwzCUkpKiqKgo+fv7a9CgQcrKyqq2dsSIETKZTPr444+dfwEAAABAA/Np1U1R47cpctxWmeP/V4E9Jskc/7+KHLdVUeO/JlgBgAbi5eoGamvChAk6fvy40tLSJEm/+MUvlJycrL/97W/Vzlm4cKEWLVqk1atX6yc/+Yn+8Ic/aMiQIcrOzlZgYKBd7eLFi1keCQAAALdnMpnkF9lPfpH9XN0KADQbbhGuHDx4UGlpaUpPT1dCQoIkacWKFUpMTFR2drY6depUaY5hGFq8eLHmzJmje+65R5L01ltvKTw8XO+9954ef/xxW+2///1vLVq0SDt37lRkZGTjXBQAAAAAALguuMVtQdu3b5fZbLYFK5J0yy23yGw2a9u2bVXOycnJkcVi0dChQ21jvr6+GjhwoN2cwsJCPfDAA1qyZIkiIiJq1U9JSYkKCgrsvgAAAAAAQPPkFuGKxWJRWFhYpfGwsDBZLJZq50hSeHi43Xh4eLjdnBkzZigpKUljxoypdT/z58+37f1iNpsVExNT67kAAABAQystq9Arm77VK5u+VWlZhavbAYDrnkvDlZSUFJlMphq/du3aJUlV7odiGMY190n58fn/nrN+/Xr961//0uLFi+vU9+zZs2W1Wm1fx44dq9N8AAAAoLki+AFwPXLpnitTp07V+PHja6yJjY3VN998o9OnT1c6d/bs2UorU666eouPxWKx20flzJkztjn/+te/dPjwYQUHB9vN/dnPfqYBAwZo8+bNVb62r6+vfH19a+z7elBaVqGlXxySJE0Z3FE+Xm6x0AkAAAAAgEbl0nAlNDRUoaGh16xLTEyU1WpVRkaG+vW7suv5jh07ZLValZSUVOWcuLg4RUREaNOmTerdu7ckqbS0VFu2bNELL7wgSZo1a5Yee+wxu3k9evTQK6+8olGjRjlyaQAAAAAAoJlwi6cFdenSRcOHD9ekSZP05z//WdKVRzGPHDnS7klBnTt31vz583X33XfLZDJp+vTpmjdvnm688UbdeOONmjdvngICAjRhwgRJV1a3VLWJbdu2bRUXF9c4FwcAAAAAANyaW4QrkvTuu+9q2rRptqf/jB49WkuWLLGryc7OltVqtR0//fTTKioq0pNPPqm8vDwlJCRo48aNCgwMbNTeAQAAAADA9cttwpUbbrhB77zzTo01hmHYHZtMJqWkpCglJaXW7/Pj1wAAAAAAAKgJO5QCAAAAAAA4gHAFAAAAAADAAbW+LWj9+vW1ftHRo0fXqxkAAAAAAAB3U+twZezYsbWqM5lMKi8vr28/AAAAAAAAbqXW4UpFRUVD9gEAAAAAAOCW2HMFAAAAAADAAbVeufLHP/6x1i86bdq0ejWD+jMMQxm5efr8u7MqulyhUyVluqdHpBLaBstkMrm6PQAAAAAArlu1DldeeeWVWtWZTCbClUaWZbmgial7teu4VWaPK0HKF8fy9dLmw4qPNuvN8b3VLSLQxV0CAAAAAHB9qnW4kpOT05B9oJ6yLBc0YOnXyi+6XOX5XcetGrD0a3015VYCFgAAAAAAGgB7rrgxwzA0MXVvtcHKVflFl/XImkwZhtFInQEAAAAA0HzUeuXKf3vkkUdqPL9q1ap6NYO62ZGbr13HrbWq3XksXxm5+UpoF9LAXQEAAAAA0LzUK1zJy8uzO758+bL279+v/Px83X777U5pDNf2SZalTvUfZ1kIV6pRWlahpV8ckiRNGdxRPl4s6gIAAAAA1E69wpV169ZVGquoqNCTTz6p9u3bO9wUaievsObbgRytBwAAAAAA1+a0v5738PDQjBkzav1UITguJMC7QesBAAAAAMC1OfXeh8OHD6usrMyZL4kajOkWUaf6sXWsBwAAAAAA11av24Jmzpxpd2wYhk6dOqVPP/1UDz30kFMaw7UltA1WfLS5Vpva9o0JVr+2wQ3fFAAAAAAAzUy9wpU9e/bIZDLZjj08PNS6dWu9/PLL13ySEJzHZDLpzfG9NWDp1zU+jjnY31urxt1k9+8MAAAAAAA4R63DlfXr12vEiBHy9vbW5s2bG7Al1EW3iEB9NeVWTUzdW+UKlr4xwVo17iZ1iwh0QXcAAAAAAFz/ah2u3H333bJYLGrdurU8PT116tQphYWFNWRvqKVuEYHa8dQAbcv5QfPT/qOiyxXq06GVftYjUv3aBtdrxYphGCq2ZKgoN1NGebHyfM/JfOMo+Ub0YwUMAAAAAAD/pdbhSuvWrZWenq5Ro0bJMAw+YDcxJpNJfduG6Kc3tpYkTRncUT5e9duvuPR8ls5ufEyXLP9WycWfSZKsl9aqcM9C+YT3Ueuhb8inVTen9Q4AAAAAgDur9afvyZMna8yYMfL09JTJZFJERIQ8PT2r/IL7Kj2fpVMfDFbp6d1Vnz+9+8r581mN3BkAAAAAAE1TrVeupKSkaPz48Tp06JBGjx6tN998U8HBwQ3YGhqbYRg6u/ExVZTk11hXUZKvsxsnKWr816xgAgAAAAA0e3V6WlDnzp3VuXNnzZ07V/fdd58CAgJqrP/6668VHx8vX19fh5pE4yixZFS7YuXHSk/vUollp/wi+zVwVwAAAAAANG312pRj7ty51wxWJGnEiBE6ceJEfd4CLlB4eH0d6z9poE4AAAAAAHAfdVq5UleGYTTky8PJrnU7kKP1AAAAcD+GYWjH2aPaYM1QYUWJLHu+092x3ZXQui23iAPA/9Og4Qrci4dvcIPWAwAAwL1k5Vn0yNZU7Tp33DaWnnVQL2V9ofjQaK3qP17dQiJc2CEANA31e1YvrksBHUbXsX5MA3UCAAAAV8vKs+i2DUvtgpX/tuvccd22Yamy8iyN3BkAND2EK7Dxjegnn/A+tar1CY+Xb0TfBu4IAAAArmAYhh7Zmqr80qIa6/JLi/To1jVsBwCg2WvQcIV7MN2LyWRS66FvXPN2Hw/fYLUeuoJ/vwAAANepHWdzq12x8mM7zx1TxrncBu4IAJq2Bg1XSLDdj0+rboq8/4tqV7D4hMdfOd+qWyN3BgAAgMayPjerTvWfHK1bPQBcb+q9oW1ZWZk2b96sw4cPa8KECQoMDNTJkycVFBSkli1bSpIuXLjgtEbReHxadVPU+G26cCJDvv/YK6O8WObO7WS+cbR8I/qyYgUAAOA6l1da2KD1AHC9qVe4cvToUQ0fPly5ubkqKSnRkCFDFBgYqIULF6q4uFjLly93dp9oZCaTSX4RfeXfNkSSFJLUUT5ebNEDAADQHIT4BDRoPQBcb+r1afmpp55SfHy88vLy5O/vbxu/++679fnnnzutOQAAAACNb3Tbut0CPqad624ZLy2r0CubvtUrm75VaVmFy/oA0LzVa+XK1q1b9fXXX8vHx8duvF27djpx4oRTGgMAAADgGgmt2yo+NLpWm9r2DY1Rv9C2jdAVADRd9Vq5UlFRofLy8krjx48fV2BgoMNNAQAAAHAdk8mkVf3HK9jHv8a6YB9/rew/jj35ADR79QpXhgwZosWLF9uOTSaTLl68qLlz5+rOO+90Vm8AAAAAXKRbSIS+vHOK4kOjqzzfNzRGX945Rd1CIhq5MwBoeup1W9Arr7yiwYMHq2vXriouLtaECRP03XffKTQ0VO+//76zewQAAADgAt1CIpQ+8il9bTmq3329VYUVJeofG6F74rqrX2hbVqwAwP9Tr3AlKipKmZmZSk1N1e7du1VRUaFHH31UDz74oN0GtwAAAADcm8lkUr/WbXWnuZ8kacrNPEUSAH6sXuHKl19+qaSkJE2cOFETJ060jZeVlenLL7/Ubbfd5rQGAQAAAAAAmrJ6Rc6DBw/WDz/8UGncarVq8ODBDjcFAAAAAADgLuq1csUwjCrvrzx//rxatGjhcFO4vhiGoR1nj2qDNUOFFSWy7PlOd8d2V0Jr7tMFAAAAALi/OoUr99xzj6Qr910+/PDD8vX1tZ0rLy/XN998o6SkJOd2CLeWlWfRI1tTtevccdtYetZBvZT1heJDo7Wq/3h2mAcAAAAAuLU6hStms1nSlZUIgYGBdpvX+vj46JZbbtGkSZOc2yHcVlaeRbdtWKr80qIqz+86d1y3bVjKI/wAAACaCcMwVHw4XTfufUvepVad+yFWwX3Gyr9DAiuaAbi1OoUrb775piQpNjZWv/71r7kFCNUyDEOPbE2tNli5Kr+0SI9uXaPtI6fxP1QAAIDrWPHxLJ1842EV5+xS+/83lv+dlL/hBfnFxSvqsdXyi+7m0h4BoL7qtaHt3LlzCVZQox1nc+1uBarJznPHlHEut4E7AgAAgKsUH8/Skef7qzhnV9Xnc3ZdOX88q5E7AwDnqNeGtpL04Ycf6oMPPlBubq5KS0vtzu3Zs8fhxuDe1ufW7X+MnxzNUkLrdg3UDQAAAFzFMAydfONhVRTm11hXUZivk29MVNzcHaxoBuB26rVy5Y9//KMmTpyosLAw7d27V/369VOrVq30/fffa8SIEc7uEW4or7SwQesBAADgHooO76h2xcqPFefsVNH3GQ3cEQA4X73Clddee02vv/66lixZIh8fHz399NPatGmTpk2bJqvV6uwe4YZCfAIatB4AAACNx8fLQzOG/EQzhvxEPl51+whxYc8ndavf/XGd6gGgKahXuJKbm2t75LK/v78uXLggSUpOTtb777/vvO7gtka3rdtmZGPa1b7+/99l/mV13fFbnVv7jAoPpcswjLq2CQAAgAZWUZjXoPUA0BTUa8+ViIgInT9/Xu3atVO7du2Unp6uXr16KScnhw+4kCQltG6r+NDoWm1q2zc0Rv1C29bqddllHgAAwL14BIQ0aD0ANAX1Wrly++23629/+5sk6dFHH9WMGTM0ZMgQjRs3TnfffbdTG7wqLy9PycnJMpvNMpvNSk5OVn5+fo1zDMNQSkqKoqKi5O/vr0GDBikrq/JGq9u3b9ftt9+uFi1aKDg4WIMGDVJRUc2PEEbNTCaTVvUfr2Af/xrrgn38tbL/uFptWsYu8wAAAO4n8OYxdavvM7ZhGgGABlSvcOX111/XnDlzJEmTJ0/W6tWr1aVLF/3ud7/TsmXLnNrgVRMmTFBmZqbS0tKUlpamzMxMJScn1zhn4cKFWrRokZYsWaKdO3cqIiJCQ4YMsd3GJF0JVoYPH66hQ4cqIyNDO3fu1NSpU+XhUa9vDf5Lt5AIfXnnFMWHRld5vm9ojL68c4q6hURc87Xquss8K6gAAACaBv8OCfKLi69VrV9cX/m379fAHQGA89XrtiAPDw+78OH+++/X/fff77SmfuzgwYNKS0tTenq6EhISJEkrVqxQYmKisrOz1alTp0pzDMPQ4sWLNWfOHN1zzz2SpLfeekvh4eF677339Pjjj0uSZsyYoWnTpmnWrFm2uTfeeGODXUtz0y0kQukjn9LXlqP63ddbVVhRov6xEbonrrv6hbat9WP26rPLfECHBEdaBwAAgBOYTCZFPbZaR57vX+NflHkEBCvqsTd5DDMAt1Tv5RnFxcXKyMjQ3//+d61fv97uy9m2b98us9lsC1Yk6ZZbbpHZbNa2bduqnJOTkyOLxaKhQ4faxnx9fTVw4EDbnDNnzmjHjh0KCwtTUlKSwsPDNXDgQG3dutXp19CcmUwm9WvdVnea++nekAH6/c0jlNC6XZ3+x8ku8wAAAO7LL7qbYudsrXYFi19c3yvn2TsPgJuq18qVtLQ0/fznP9e5c+cqnTOZTCovL3e4sf9msVgUFhZWaTwsLEwWi6XaOZIUHh5uNx4eHq6jR49Kkr7//ntJUkpKil566SXddNNNevvtt/XTn/5U+/fvr3YFS0lJiUpKSmzHBQUFdb8o1Am7zAMAALg3v+huipuboQvfpevLj1fLu9SqHjfGKjj+bvm378eKFQBurV4rV6ZOnar77rtPp06dUkVFhd1XXYKVlJQUmUymGr927bpyK0hVv2wNw7jmL+Efn//vORUVFZKkxx9/XBMnTlTv3r31yiuvqFOnTlq1alW1rzl//nzbxrpms1kxMTG1vmbUD7vMAwAAuD+TySS/9gn6rvevdCDhOYX+bJ4COiQQrABwe/VauXLmzBnNnDmz0qqQupo6darGjx9fY01sbKy++eYbnT59utK5s2fPVttDRMSVTVItFosiIyNt42fOnLHNuTretWtXu7ldunRRbm5utT3Nnj1bM2fOtB0XFBQQsDSwwJvH6PynC2pfzy7zAAAAAIBGUq9w5d5779XmzZvVoUMHh948NDRUoaGh16xLTEyU1WpVRkaG+vW7snv4jh07ZLValZSUVOWcuLg4RUREaNOmTerdu7ckqbS0VFu2bNELL7wg6UpwExUVpezsbLu53377rUaMGFFtP76+vvL19a3VNcI5ru4yX5tNbdllHgAAAADQmOoVrixZskT33XefvvrqK/Xo0UPe3t5256dNm+aU5q7q0qWLhg8frkmTJunPf/6zJOkXv/iFRo4cafekoM6dO2v+/Pm6++67ZTKZNH36dM2bN0833nijbrzxRs2bN08BAQGaMGGCpCvLEv/3f/9Xc+fOVa9evXTTTTfprbfe0n/+8x99+OGHTr0GOIZd5gEAAAAATVW9wpX33ntP//jHP+Tv76/NmzfbfZA1mUxOD1ck6d1339W0adNsT/8ZPXq0lixZYleTnZ0tq9VqO3766adVVFSkJ598Unl5eUpISNDGjRsVGBhoq5k+fbqKi4s1Y8YM/fDDD+rVq5c2bdrk8KocON/VXeZPvvFwlStY/OL6KuqxN9llHgAAAADQqOoVrjz77LN67rnnNGvWLHl41PtpznVyww036J133qmxxjAMu2OTyaSUlBSlpKTUOG/WrFmaNWuWoy2iEVzdZb7o+wxd2P2xKgrz5BEQosA+Y9llHgAAAADgEvUKV0pLSzVu3LhGC1aA/2YymRTQIUEBHRJc3QoAAAAAAPV7FPNDDz2kNWvWOLsXAAAAAAAAt1OvlSvl5eVauHCh/vGPf6hnz56VNrRdtGiRU5oDAAAAAABo6uoVruzbt8/2eOP9+/fbnWPPCwAAAAAA0JzUK1z54osvnN0HAAAAAACAW2JHWgAAAAAAAAfUeuXKPffco9WrVysoKEj33HNPjbUfffSRw40BAAAAuP74eHloxpCfuLoNAHCqWocrZrPZtp9KUFAQe6sAAAAAAACoDuHKm2++afvn1atXN0QvAAAAAAAAbqdee67cfvvtys/PrzReUFCg22+/3dGeAAAAAAAA3Ea9wpXNmzertLS00nhxcbG++uorh5sCAAAAAABwF3V6FPM333xj++cDBw7IYrHYjsvLy5WWlqY2bdo4rzsAAACgmWCjVwBwX3UKV2666SaZTCaZTKYqb//x9/fXn/70J6c1BwAAAAAA0NTVKVzJycmRYRhq3769MjIy1Lp1a9s5Hx8fhYWFydPT0+lNAgAAAAAANFV1ClfatWsnSaqoqGiQZuAYlpICAAAAAND46rWh7VtvvaVPP/3Udvz0008rODhYSUlJOnr0qNOaAwAAAAAAaOrqFa7MmzdP/v7+kqTt27dryZIlWrhwoUJDQzVjxgynNggAAAAAANCU1em2oKuOHTumjh07SpI+/vhj3XvvvfrFL36hW2+9VYMGDXJmfwAAAAAAAE1avVautGzZUufPn5ckbdy4UXfccYckyc/PT0VFRc7rDgAAAAAAoImr18qVIUOG6LHHHlPv3r317bff6q677pIkZWVl2Ta9BQAAAAAAaA7qtXJl6dKlSkpK0rlz5/TRRx+pVatWkqTdu3drwoQJTm0QAAAAAACgKatXuBIcHKz77rtPLVq0UEpKik6cOCFJ6tChgwYOHOjUBgEAAAAAAJqyeoUra9eu1fDhwxUQEKC9e/eqpKREknTx4kXNmzfPqQ0CAAAAAAA0ZfUKV/7whz9o+fLlWrFihby9vW3jSUlJ2rNnj9OaAxqLj5eHZgz5iWYM+Yl8vOr1YwEAAAAAaKbqtaFtdna2brvttkrjQUFBys/Pd7QnoFEZhqFLxwpkPXBWZUVl8vL3krlra7WICZLJZHJ1ewAAAACAJq5e4UpkZKQOHTqk2NhYu/GtW7eqffv2zugLaBRFpy/qyF8PqPDEBbtxy5ajCmgTqNj7uso/vKWLugMAAAAAuIN63f/w+OOP66mnntKOHTtkMpl08uRJvfvuu/r1r3+tJ5980tk9Ag2i6PRFZS/fXSlYuarwxAVlL9+totMXG7kzAAAAAIA7qdfKlaefflpWq1WDBw9WcXGxbrvtNvn6+urXv/61pk6d6uweAaczDENH/npA5cVlNdaVF5fpyIcH1fnJeG4RAgAAAABUqV7hiiQ9//zzmjNnjg4cOKCKigp17dpVLVty+wTcw6VjBdWuWPmxwuMFKjxWoBZtzQ3cFQAAAADAHdU7XJGkgIAAxcfHO6sXoNFYD5ytU33+gbOEKwAAAACAKvHMWTRLZUU13w7kaD0AAAAAoPkgXEGz5OVft0Vbda0HAAAAADQfhCtolsxdW9epPriO9QAAAACA5oNwBc1Si5ggBbQJrFVtQHSQAmKCGrgjAAAAAIC7IlxBs2QymRR7X1d5+tV8u4+nn5di7+3CY5gBAAAAANUiXEGz5R/eUp0m96l2BUtAdJA6Te4j/3AeMQ4AAAAAqB67dKJaPl4emjHkJ65uo0H5h7dU5yl9VXisQPkHzqqsqExe/l4K7tpaATFBrFgBAAAAAFwT4QqaPZPJpBZtzWrR1uzqVgAAAAAAbojbggAAAAAAABxAuAIAAAAAAOAAwhUAAAAAAAAHEK4AAAAAAAA4gHAFAAAAAADAAYQrAAAAAAAADiBcAQAAAAAAcADhCgAAAAAAgAMIVwAAAAAAABxAuAIAAAAAAOAAwhUAAAAAAAAHEK4AAAAAAAA4wG3Clby8PCUnJ8tsNstsNis5OVn5+fk1zjEMQykpKYqKipK/v78GDRqkrKwsuxqLxaLk5GRFRESoRYsWuvnmm/Xhhx824JUAAAAAAIDriduEKxMmTFBmZqbS0tKUlpamzMxMJScn1zhn4cKFWrRokZYsWaKdO3cqIiJCQ4YM0YULF2w1ycnJys7O1vr167Vv3z7dc889GjdunPbu3dvQlwQAAAAAAK4DbhGuHDx4UGlpaXrjjTeUmJioxMRErVixQn//+9+VnZ1d5RzDMLR48WLNmTNH99xzj7p376633npLhYWFeu+992x127dv1y9/+Uv169dP7du317PPPqvg4GDt2bOnsS4PAAAAAAC4MbcIV7Zv3y6z2ayEhATb2C233CKz2axt27ZVOScnJ0cWi0VDhw61jfn6+mrgwIF2c/r37681a9bohx9+UEVFhVJTU1VSUqJBgwY12PUAAAAAAIDrh5erG6gNi8WisLCwSuNhYWGyWCzVzpGk8PBwu/Hw8HAdPXrUdrxmzRqNGzdOrVq1kpeXlwICArRu3Tp16NCh2n5KSkpUUlJiOy4oKKjT9QAAAAAAgOuHS1eupKSkyGQy1fi1a9cuSZLJZKo03zCMKsf/24/P/3jOs88+q7y8PP3zn//Url27NHPmTN13333at29fta85f/5828a6ZrNZMTExdblsAAAAoEH5eHloxpCfaMaQn8jHyy0WqwOAW3PpypWpU6dq/PjxNdbExsbqm2++0enTpyudO3v2bKWVKVdFRERIurKCJTIy0jZ+5swZ25zDhw9ryZIl2r9/v7p16yZJ6tWrl7766istXbpUy5cvr/K1Z8+erZkzZ9qOCwoKCFgAAAAAAGimXBquhIaGKjQ09Jp1iYmJslqtysjIUL9+/SRJO3bskNVqVVJSUpVz4uLiFBERoU2bNql3796SpNLSUm3ZskUvvPCCJKmwsFCS5OFhn+Z7enqqoqKi2n58fX3l6+t77QsEAAAAAADXPbdYI9ilSxcNHz5ckyZNUnp6utLT0zVp0iSNHDlSnTp1stV17txZ69atk3TldqDp06dr3rx5Wrdunfbv36+HH35YAQEBmjBhgq2+Y8eOevzxx5WRkaHDhw/r5Zdf1qZNmzR27FhXXCoAAAAAAHAzbrGhrSS9++67mjZtmu3pP6NHj9aSJUvsarKzs2W1Wm3HTz/9tIqKivTkk08qLy9PCQkJ2rhxowIDAyVJ3t7e2rBhg2bNmqVRo0bp4sWL6tixo9566y3deeedjXdxAAAAAOrMMAwVHrMq9nC+vC5XyHL5sG7oHqYWMUHX3JsRAJzJZBiG4eom3F1BQYHMZrOsVquCgoJc3U6TVFpWoaVfHJIkTRnckY3VAAAA4JCi0xd15K8HVHjiQqVzAW0CFXtfV/mHt2z0vvhsADRPfMIFAAAA4FaKTl9U9vLdVQYrklR44oKyl+9W0emLjdwZgOaKcAUAAACA2zAMQ0f+ekDlxWU11pUXl+nIhwfFQn0AjYFwBQAAAIDbuHSsoNoVKz9WeLxAhccKGrgjACBcAQAAAOBGrAfO1qk+v471AFAfhCsAAAAA3EZZUc23AzlaDwD1QbgCAAAAwG14+Xs1aD0A1AfhCgAAAAC3Ye7auk71wXWsB4D6IFwBAAAA4DZaxAQpoE1grWoDooMUEBPUwB0BAOEKAAAAADdiMpkUe19XefrVfLuPp5+XYu/tIpPJ1EidAWjOCFcAAAAAuBX/8JbqNLlPtStYAqKD1GlyH/mHt2zkzgA0V+zuBAAAAMDt+Ie3VOcpfWU9kq8v/vGdvC5XqFeHVmrVPUwBMUGsWAHQqAhXAAAAALglk8mkgBizjnQIliQNG9xBPl4szgfQ+PjNAwAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAA9wmXMnLy1NycrLMZrPMZrOSk5OVn59f45yPPvpIw4YNU2hoqEwmkzIzMyvVlJSU6Je//KVCQ0PVokULjR49WsePH2+YiwAAAAAAANcdtwlXJkyYoMzMTKWlpSktLU2ZmZlKTk6ucc6lS5d06623asGCBdXWTJ8+XevWrVNqaqq2bt2qixcvauTIkSovL3f2JQAAAAAAgOuQl6sbqI2DBw8qLS1N6enpSkhIkCStWLFCiYmJys7OVqdOnaqcdzV8OXLkSJXnrVarVq5cqb/85S+64447JEnvvPOOYmJi9M9//lPDhg1z/sUAAAAAAIDrilusXNm+fbvMZrMtWJGkW265RWazWdu2bav36+7evVuXL1/W0KFDbWNRUVHq3r17ja9bUlKigoICuy8AAAAAANA8uUW4YrFYFBYWVmk8LCxMFovFodf18fFRSEiI3Xh4eHiNrzt//nzb3i9ms1kxMTH17gEAAAAAALg3l4YrKSkpMplMNX7t2rVLkmQymSrNNwyjynFHXet1Z8+eLavVavs6duyY03sAAAAAAADuwaV7rkydOlXjx4+vsSY2NlbffPONTp8+Xenc2bNnFR4eXu/3j4iIUGlpqfLy8uxWr5w5c0ZJSUnVzvP19ZWvr2+93xcAAAAAAFw/XBquhIaGKjQ09Jp1iYmJslqtysjIUL9+/SRJO3bskNVqrTEEuZY+ffrI29tbmzZt0v333y9JOnXqlPbv36+FCxfW+3UBAAAAAEDz4RZ7rnTp0kXDhw/XpEmTlJ6ervT0dE2aNEkjR460e1JQ586dtW7dOtvxDz/8oMzMTB04cECSlJ2drczMTNt+KmazWY8++qh+9atf6fPPP9fevXv1P//zP+rRo4ft6UEAAAAAAAA1cYtwRZLeffdd9ejRQ0OHDtXQoUPVs2dP/eUvf7Gryc7OltVqtR2vX79evXv31l133SVJGj9+vHr37q3ly5fbal555RWNHTtW999/v2699VYFBATob3/7mzw9PRvnwgAAAAAAgFszGYZhuLoJd1dQUCCz2Syr1aqgoCBXt9MklZZVaOkXhyRJUwZ3lI+X2+R6AAAAaMKa2p8z+WwANE98wgUAAAAAAHAA4QoAAAAAAIADCFcAAAAAAAAcQLgCAAAAAADgAMIVAAAAAAAABxCuAAAAAAAAOIBwBQAAAAAAwAGEKwAAAAAAAA4gXAEAAAAAAHAA4QoAAAAAAIADCFcAAAAAAAAcQLgCAAAAAADgAMIVAAAAAAAABxCuAAAAAAAAOIBwBQAAAAAAwAGEKwAAAAAAAA4gXAEAAAAAAHAA4QoAAAAAAIADCFcAAAAAAAAcQLgCAAAAAADgAMIVAAAAAAAABxCuAAAAAAAAOIBwBQAAAAAAwAGEKwAAAAAAAA4gXAEAAAAAAHAA4QoAAAAAAIADCFcAAAAAAAAcQLgCAAAAAADgAMIVAAAAAAAABxCuAAAAAAAAOIBwBQAAAAAAwAGEKwAAAAAAAA4gXAEAAAAAAHAA4QoAAAAAAIADvFzdAJoHHy8PzRjyE1e3AQAAAACA0xGuAAAAAHBb/CUegKaA24IAAAAAAAAcQLgCAAAAAADgAMIVAAAAAAAABxCuAAAAAAAAOIBwBQAAAAAAwAGEKwAAAAAAAA4gXAEAAAAAAHAA4QoAAAAAAIADCFcAAAAAAAAcQLgCAAAAAADgAMIVAAAAAAAABxCuAAAAAAAAOIBwBQAAAAAAwAGEKwAAAAAAAA4gXAEAAAAAAHAA4QoAAAAAAIADvFzdwPXAMAxJUkFBgYs7AQAAAOBKVz8TXP2MAKB5IFxxggsXLkiSYmJiXNwJAAAAgKbgwoULMpvNrm4DQCMxGUSqDquoqNDJkycVGBgok8nk6naarIKCAsXExOjYsWMKCgpydTsAAACoA/4sVzuGYejChQuKioqShwe7MADNBStXnMDDw0PR0dGubsNtBAUF8T9kAAAAN8Wf5a6NFStA80OUCgAAAAAA4ADCFQAAAAAAAAcQrqDR+Pr6au7cufL19XV1KwAAAKgj/iwHANVjQ1sAAAAAAAAHsHIFAAAAAADAAYQrAAAAAAAADiBcAQAAAAAAcADhCgAAAAAAgAMIV9BoXnvtNcXFxcnPz099+vTRV1995eqWAAAAcA1ffvmlRo0apaioKJlMJn388ceubgkAmhzCFTSKNWvWaPr06ZozZ4727t2rAQMGaMSIEcrNzXV1awAAAKjBpUuX1KtXLy1ZssTVrQBAk8WjmNEoEhISdPPNN2vZsmW2sS5dumjs2LGaP3++CzsDAABAbZlMJq1bt05jx451dSsA0KSwcgUNrrS0VLt379bQoUPtxocOHapt27a5qCsAAAAAAJyDcAUN7ty5cyovL1d4eLjdeHh4uCwWi4u6AgAAAADAOQhX0GhMJpPdsWEYlcYAAAAAAHA3hCtocKGhofL09Ky0SuXMmTOVVrMAAAAAAOBuCFfQ4Hx8fNSnTx9t2rTJbnzTpk1KSkpyUVcAAAAAADiHl6sbQPMwc+ZMJScnKz4+XomJiXr99deVm5uryZMnu7o1AAAA1ODixYs6dOiQ7TgnJ0eZmZm64YYb1LZtWxd2BgBNB49iRqN57bXXtHDhQp06dUrdu3fXK6+8ottuu83VbQEAAKAGmzdv1uDBgyuNP/TQQ1q9enXjNwQATRDhCgAAAAAAgAPYcwUAAAAAAMABhCsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAAAAOIFwBAKAW0tLS1L9/fwUHB6tVq1YaOXKkDh8+LEk6cuSITCaTPvjgAw0YMED+/v7q27evvv32W+3cuVPx8fFq2bKlhg8frrNnz7r4SgAAAOBsJsMwDFc3AQBAU7d27VqZTCb16NFDly5d0m9/+1sdOXJEmZmZys3NVVxcnDp37qzFixerbdu2euSRR1RaWqqgoCD94Q9/UEBAgO6//37dcccdWrZsmasvBwAAAE5EuAIAQD2cPXtWYWFh2rdvn1q2bKm4uDi98cYbevTRRyVJqampeuCBB/T555/r9ttvlyQtWLBAq1ev1n/+8x9Xtg4AAAAn47YgAABq4fDhw5owYYLat2+voKAgxcXFSZJyc3NtNT179rT9c3h4uCSpR48edmNnzpxppI4BAADQWLxc3QAAAO5g1KhRiomJ0YoVKxQVFaWKigp1795dpaWlthpvb2/bP5tMpirHKioqGq9pAAAANArCFQAAruH8+fM6ePCg/vznP2vAgAGSpK1bt7q4KwAAADQVhCsAAFxDSEiIWrVqpddff12RkZHKzc3VrFmzXN0WAAAAmgj2XAEA4Bo8PDyUmpqq3bt3q3v37poxY4ZefPFFV7cFAACAJoKnBQEAAAAAADiAlSsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHEC4AgAAAAAA4ADCFQAAAAAAAAcQrgAAAAAAADiAcAUAAAAAAMABhCsAAAAAAAAOIFwBAAAAAABwAOEKAAAAAACAAwhXAAAAAAAAHPD/AXXJDiTAaDVlAAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 1000x500 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(10, 5))\n",
"ax = sns.pointplot(\n",
" data=contrast_df,\n",
" x=\"am\",\n",
" y=\"estimate_full\",\n",
" hue=\"drat\",\n",
" dodge=True,\n",
" join=False,\n",
" palette=\"colorblind\"\n",
")\n",
"\n",
"x_coords = []\n",
"y_coords = []\n",
"for point_pair in ax.collections:\n",
" for x, y in point_pair.get_offsets():\n",
" x_coords.append(x)\n",
" y_coords.append(y)\n",
"\n",
"plt.vlines(x_coords, contrast_df[\"std_full_lower\"], contrast_df[\"std_full_higher\"], alpha=0.5)\n",
"plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "bambinos",
"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"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment